LArRotationalTransformationPlugin.cc
Go to the documentation of this file.
1 /**
2  * @file larpandoracontent/LArPlugins/LArRotationalTransformationPlugin.cc
3  *
4  * @brief Implementation of the rotational transformation plugin class.
5  *
6  * $Log: $
7  */
8 
9 #include "Geometry/LArTPC.h"
10 
11 #include "Helpers/XmlHelper.h"
12 
13 #include "Managers/GeometryManager.h"
14 
15 #include "Objects/CartesianVector.h"
16 
17 #include "Pandora/Pandora.h"
18 
20 
21 #include <cmath>
22 
23 namespace lar_content
24 {
25 
26 using namespace pandora;
27 
29  m_thetaU(0.),
30  m_thetaV(0.),
31  m_thetaW(0.),
32  m_sinU(0.),
33  m_sinV(0.),
34  m_sinW(0.),
35  m_cosU(0.),
36  m_cosV(0.),
37  m_cosW(0.),
38  m_sinVminusU(0.),
39  m_sinWminusV(0.),
40  m_sinUminusW(0.),
41  m_maxAngularDiscrepancyU(0.03),
42  m_maxAngularDiscrepancyV(0.03),
43  m_maxAngularDiscrepancyW(0.03),
44  m_maxSigmaDiscrepancy(0.01)
45 {
46 }
47 
48 //------------------------------------------------------------------------------------------------------------------------------------------
49 
50 double LArRotationalTransformationPlugin::UVtoW(const double u, const double v) const
51 {
52  return (-1. * (u * m_sinWminusV + v * m_sinUminusW) / m_sinVminusU);
53 }
54 
55 //------------------------------------------------------------------------------------------------------------------------------------------
56 
57 double LArRotationalTransformationPlugin::VWtoU(const double v, const double w) const
58 {
59  return (-1. * (v * m_sinUminusW + w * m_sinVminusU) / m_sinWminusV);
60 }
61 
62 //------------------------------------------------------------------------------------------------------------------------------------------
63 
64 double LArRotationalTransformationPlugin::WUtoV(const double w, const double u) const
65 {
66  return (-1. * (u * m_sinWminusV + w * m_sinVminusU) / m_sinUminusW);
67 }
68 
69 //------------------------------------------------------------------------------------------------------------------------------------------
70 
71 double LArRotationalTransformationPlugin::UVtoY(const double u, const double v) const
72 {
73  return ((u * m_cosV - v * m_cosU) / m_sinVminusU);
74 }
75 
76 //------------------------------------------------------------------------------------------------------------------------------------------
77 
78 double LArRotationalTransformationPlugin::UVtoZ(const double u, const double v) const
79 {
80  return ((u * m_sinV - v * m_sinU) / m_sinVminusU);
81 }
82 
83 //------------------------------------------------------------------------------------------------------------------------------------------
84 
85 double LArRotationalTransformationPlugin::UWtoY(const double u, const double w) const
86 {
87  return ((w * m_cosU - u * m_cosW) / m_sinUminusW);
88 }
89 
90 //------------------------------------------------------------------------------------------------------------------------------------------
91 
92 double LArRotationalTransformationPlugin::UWtoZ(const double u, const double w) const
93 {
94  return ((w * m_sinU - u * m_sinW) / m_sinUminusW);
95 }
96 
97 //------------------------------------------------------------------------------------------------------------------------------------------
98 
99 double LArRotationalTransformationPlugin::VWtoY(const double v, const double w) const
100 {
101  return ((v * m_cosW - w * m_cosV) / m_sinWminusV);
102 }
103 
104 //------------------------------------------------------------------------------------------------------------------------------------------
105 
106 double LArRotationalTransformationPlugin::VWtoZ(const double v, const double w) const
107 {
108  return ((v * m_sinW - w * m_sinV) / m_sinWminusV);
109 }
110 
111 //------------------------------------------------------------------------------------------------------------------------------------------
112 
113 double LArRotationalTransformationPlugin::YZtoU(const double y, const double z) const
114 {
115  return (z * m_cosU - y * m_sinU);
116 }
117 
118 //------------------------------------------------------------------------------------------------------------------------------------------
119 
120 double LArRotationalTransformationPlugin::YZtoV(const double y, const double z) const
121 {
122  return (z * m_cosV - y * m_sinV);
123 }
124 
125 //------------------------------------------------------------------------------------------------------------------------------------------
126 
127 double LArRotationalTransformationPlugin::YZtoW(const double y, const double z) const
128 {
129  return (z * m_cosW - y * m_sinW);
130 }
131 
132 //------------------------------------------------------------------------------------------------------------------------------------------
133 
134 void LArRotationalTransformationPlugin::GetMinChiSquaredYZ(const double u, const double v, const double w, const double sigmaU,
135  const double sigmaV, const double sigmaW, double &y, double &z, double &chiSquared) const
136 {
137  const double sigmaU2(sigmaU * sigmaU), sigmaV2(sigmaV * sigmaV), sigmaW2(sigmaW * sigmaW);
138 
139  // Obtain expression for chi2, differentiate wrt y and z, set both results to zero and solve simultaneously. Here just paste-in result.
140  y = (sigmaW2 * v * m_cosU * m_cosV * m_sinU - sigmaW2 * u * m_cosV * m_cosV * m_sinU + sigmaV2 * w * m_cosU * m_cosW * m_sinU -
141  sigmaV2 * u * m_cosW * m_cosW * m_sinU - sigmaW2 * v * m_cosU * m_cosU * m_sinV + sigmaW2 * u * m_cosU * m_cosV * m_sinV +
142  sigmaU2 * w * m_cosV * m_cosW * m_sinV - sigmaU2 * v * m_cosW * m_cosW * m_sinV - sigmaV2 * w * m_cosU * m_cosU * m_sinW -
143  sigmaU2 * w * m_cosV * m_cosV * m_sinW + sigmaV2 * u * m_cosU * m_cosW * m_sinW + sigmaU2 * v * m_cosV * m_cosW * m_sinW) /
144  (sigmaW2 * m_cosV * m_cosV * m_sinU * m_sinU + sigmaV2 * m_cosW * m_cosW * m_sinU * m_sinU - 2. * sigmaW2 * m_cosU * m_cosV * m_sinU * m_sinV +
145  sigmaW2 * m_cosU * m_cosU * m_sinV * m_sinV + sigmaU2 * m_cosW * m_cosW * m_sinV * m_sinV -
146  2. * sigmaV2 * m_cosU * m_cosW * m_sinU * m_sinW - 2. * sigmaU2 * m_cosV * m_cosW * m_sinV * m_sinW +
147  sigmaV2 * m_cosU * m_cosU * m_sinW * m_sinW + sigmaU2 * m_cosV * m_cosV * m_sinW * m_sinW);
148 
149  z = (sigmaW2 * v * m_cosV * m_sinU * m_sinU + sigmaV2 * w * m_cosW * m_sinU * m_sinU - sigmaW2 * v * m_cosU * m_sinU * m_sinV -
150  sigmaW2 * u * m_cosV * m_sinU * m_sinV + sigmaW2 * u * m_cosU * m_sinV * m_sinV + sigmaU2 * w * m_cosW * m_sinV * m_sinV -
151  sigmaV2 * w * m_cosU * m_sinU * m_sinW - sigmaV2 * u * m_cosW * m_sinU * m_sinW - sigmaU2 * w * m_cosV * m_sinV * m_sinW -
152  sigmaU2 * v * m_cosW * m_sinV * m_sinW + sigmaV2 * u * m_cosU * m_sinW * m_sinW + sigmaU2 * v * m_cosV * m_sinW * m_sinW) /
153  (sigmaW2 * m_cosV * m_cosV * m_sinU * m_sinU + sigmaV2 * m_cosW * m_cosW * m_sinU * m_sinU - 2. * sigmaW2 * m_cosU * m_cosV * m_sinU * m_sinV +
154  sigmaW2 * m_cosU * m_cosU * m_sinV * m_sinV + sigmaU2 * m_cosW * m_cosW * m_sinV * m_sinV -
155  2. * sigmaV2 * m_cosU * m_cosW * m_sinU * m_sinW - 2. * sigmaU2 * m_cosV * m_cosW * m_sinV * m_sinW +
156  sigmaV2 * m_cosU * m_cosU * m_sinW * m_sinW + sigmaU2 * m_cosV * m_cosV * m_sinW * m_sinW);
157 
158  const double deltaU(u - LArRotationalTransformationPlugin::YZtoU(y, z));
159  const double deltaV(v - LArRotationalTransformationPlugin::YZtoV(y, z));
160  const double deltaW(w - LArRotationalTransformationPlugin::YZtoW(y, z));
161  chiSquared = ((deltaU * deltaU) / sigmaU2) + ((deltaV * deltaV) / sigmaV2) + ((deltaW * deltaW) / sigmaW2);
162 }
163 
164 //------------------------------------------------------------------------------------------------------------------------------------------
165 
166 void LArRotationalTransformationPlugin::GetMinChiSquaredYZ(const double u, const double v, const double w, const double sigmaU, const double sigmaV,
167  const double sigmaW, const double uFit, const double vFit, const double wFit, const double sigmaFit, double &y, double &z, double &chiSquared) const
168 {
169  const double sigmaU2(sigmaU * sigmaU), sigmaV2(sigmaV * sigmaV), sigmaW2(sigmaW * sigmaW), sigmaFit2(sigmaFit * sigmaFit);
170 
171  // Obtain expression for chi2, differentiate wrt y and z, set both results to zero and solve simultaneously. Here just paste-in result.
172  y = (vFit * sigmaU2 * sigmaV2 * sigmaW2 * m_cosU * m_cosV * m_sinU + vFit * sigmaV2 * sigmaW2 * sigmaFit2 * m_cosU * m_cosV * m_sinU +
173  sigmaU2 * sigmaW2 * sigmaFit2 * v * m_cosU * m_cosV * m_sinU + sigmaW2 * sigmaFit2 * sigmaFit2 * v * m_cosU * m_cosV * m_sinU -
174  uFit * sigmaU2 * sigmaV2 * sigmaW2 * m_cosV * m_cosV * m_sinU - uFit * sigmaU2 * sigmaW2 * sigmaFit2 * m_cosV * m_cosV * m_sinU -
175  sigmaV2 * sigmaW2 * sigmaFit2 * u * m_cosV * m_cosV * m_sinU - sigmaW2 * sigmaFit2 * sigmaFit2 * u * m_cosV * m_cosV * m_sinU +
176  wFit * sigmaU2 * sigmaV2 * sigmaW2 * m_cosU * m_cosW * m_sinU + wFit * sigmaV2 * sigmaW2 * sigmaFit2 * m_cosU * m_cosW * m_sinU +
177  sigmaU2 * sigmaV2 * sigmaFit2 * w * m_cosU * m_cosW * m_sinU + sigmaV2 * sigmaFit2 * sigmaFit2 * w * m_cosU * m_cosW * m_sinU -
178  uFit * sigmaU2 * sigmaV2 * sigmaW2 * m_cosW * m_cosW * m_sinU - uFit * sigmaU2 * sigmaV2 * sigmaFit2 * m_cosW * m_cosW * m_sinU -
179  sigmaV2 * sigmaW2 * sigmaFit2 * u * m_cosW * m_cosW * m_sinU - sigmaV2 * sigmaFit2 * sigmaFit2 * u * m_cosW * m_cosW * m_sinU -
180  vFit * sigmaU2 * sigmaV2 * sigmaW2 * m_cosU * m_cosU * m_sinV - vFit * sigmaV2 * sigmaW2 * sigmaFit2 * m_cosU * m_cosU * m_sinV -
181  sigmaU2 * sigmaW2 * sigmaFit2 * v * m_cosU * m_cosU * m_sinV - sigmaW2 * sigmaFit2 * sigmaFit2 * v * m_cosU * m_cosU * m_sinV +
182  uFit * sigmaU2 * sigmaV2 * sigmaW2 * m_cosU * m_cosV * m_sinV + uFit * sigmaU2 * sigmaW2 * sigmaFit2 * m_cosU * m_cosV * m_sinV +
183  sigmaV2 * sigmaW2 * sigmaFit2 * u * m_cosU * m_cosV * m_sinV + sigmaW2 * sigmaFit2 * sigmaFit2 * u * m_cosU * m_cosV * m_sinV +
184  wFit * sigmaU2 * sigmaV2 * sigmaW2 * m_cosV * m_cosW * m_sinV + wFit * sigmaU2 * sigmaW2 * sigmaFit2 * m_cosV * m_cosW * m_sinV +
185  sigmaU2 * sigmaV2 * sigmaFit2 * w * m_cosV * m_cosW * m_sinV + sigmaU2 * sigmaFit2 * sigmaFit2 * w * m_cosV * m_cosW * m_sinV -
186  vFit * sigmaU2 * sigmaV2 * sigmaW2 * m_cosW * m_cosW * m_sinV - vFit * sigmaU2 * sigmaV2 * sigmaFit2 * m_cosW * m_cosW * m_sinV -
187  sigmaU2 * sigmaW2 * sigmaFit2 * v * m_cosW * m_cosW * m_sinV - sigmaU2 * sigmaFit2 * sigmaFit2 * v * m_cosW * m_cosW * m_sinV -
188  wFit * sigmaU2 * sigmaV2 * sigmaW2 * m_cosU * m_cosU * m_sinW - wFit * sigmaV2 * sigmaW2 * sigmaFit2 * m_cosU * m_cosU * m_sinW -
189  sigmaU2 * sigmaV2 * sigmaFit2 * w * m_cosU * m_cosU * m_sinW - sigmaV2 * sigmaFit2 * sigmaFit2 * w * m_cosU * m_cosU * m_sinW -
190  wFit * sigmaU2 * sigmaV2 * sigmaW2 * m_cosV * m_cosV * m_sinW - wFit * sigmaU2 * sigmaW2 * sigmaFit2 * m_cosV * m_cosV * m_sinW -
191  sigmaU2 * sigmaV2 * sigmaFit2 * w * m_cosV * m_cosV * m_sinW - sigmaU2 * sigmaFit2 * sigmaFit2 * w * m_cosV * m_cosV * m_sinW +
192  uFit * sigmaU2 * sigmaV2 * sigmaW2 * m_cosU * m_cosW * m_sinW + uFit * sigmaU2 * sigmaV2 * sigmaFit2 * m_cosU * m_cosW * m_sinW +
193  sigmaV2 * sigmaW2 * sigmaFit2 * u * m_cosU * m_cosW * m_sinW + sigmaV2 * sigmaFit2 * sigmaFit2 * u * m_cosU * m_cosW * m_sinW +
194  vFit * sigmaU2 * sigmaV2 * sigmaW2 * m_cosV * m_cosW * m_sinW + vFit * sigmaU2 * sigmaV2 * sigmaFit2 * m_cosV * m_cosW * m_sinW +
195  sigmaU2 * sigmaW2 * sigmaFit2 * v * m_cosV * m_cosW * m_sinW + sigmaU2 * sigmaFit2 * sigmaFit2 * v * m_cosV * m_cosW * m_sinW) /
196  (sigmaU2 * sigmaV2 * sigmaW2 * m_cosV * m_cosV * m_sinU * m_sinU + sigmaU2 * sigmaW2 * sigmaFit2 * m_cosV * m_cosV * m_sinU * m_sinU +
197  sigmaV2 * sigmaW2 * sigmaFit2 * m_cosV * m_cosV * m_sinU * m_sinU + sigmaW2 * sigmaFit2 * sigmaFit2 * m_cosV * m_cosV * m_sinU * m_sinU +
198  sigmaU2 * sigmaV2 * sigmaW2 * m_cosW * m_cosW * m_sinU * m_sinU + sigmaU2 * sigmaV2 * sigmaFit2 * m_cosW * m_cosW * m_sinU * m_sinU +
199  sigmaV2 * sigmaW2 * sigmaFit2 * m_cosW * m_cosW * m_sinU * m_sinU + sigmaV2 * sigmaFit2 * sigmaFit2 * m_cosW * m_cosW * m_sinU * m_sinU -
200  2. * sigmaU2 * sigmaV2 * sigmaW2 * m_cosU * m_cosV * m_sinU * m_sinV - 2. * sigmaU2 * sigmaW2 * sigmaFit2 * m_cosU * m_cosV * m_sinU * m_sinV -
201  2. * sigmaV2 * sigmaW2 * sigmaFit2 * m_cosU * m_cosV * m_sinU * m_sinV -
202  2. * sigmaW2 * sigmaFit2 * sigmaFit2 * m_cosU * m_cosV * m_sinU * m_sinV + sigmaU2 * sigmaV2 * sigmaW2 * m_cosU * m_cosU * m_sinV * m_sinV +
203  sigmaU2 * sigmaW2 * sigmaFit2 * m_cosU * m_cosU * m_sinV * m_sinV + sigmaV2 * sigmaW2 * sigmaFit2 * m_cosU * m_cosU * m_sinV * m_sinV +
204  sigmaW2 * sigmaFit2 * sigmaFit2 * m_cosU * m_cosU * m_sinV * m_sinV + sigmaU2 * sigmaV2 * sigmaW2 * m_cosW * m_cosW * m_sinV * m_sinV +
205  sigmaU2 * sigmaV2 * sigmaFit2 * m_cosW * m_cosW * m_sinV * m_sinV + sigmaU2 * sigmaW2 * sigmaFit2 * m_cosW * m_cosW * m_sinV * m_sinV +
206  sigmaU2 * sigmaFit2 * sigmaFit2 * m_cosW * m_cosW * m_sinV * m_sinV - 2. * sigmaU2 * sigmaV2 * sigmaW2 * m_cosU * m_cosW * m_sinU * m_sinW -
207  2. * sigmaU2 * sigmaV2 * sigmaFit2 * m_cosU * m_cosW * m_sinU * m_sinW - 2. * sigmaV2 * sigmaW2 * sigmaFit2 * m_cosU * m_cosW * m_sinU * m_sinW -
208  2. * sigmaV2 * sigmaFit2 * sigmaFit2 * m_cosU * m_cosW * m_sinU * m_sinW -
209  2. * sigmaU2 * sigmaV2 * sigmaW2 * m_cosV * m_cosW * m_sinV * m_sinW - 2. * sigmaU2 * sigmaV2 * sigmaFit2 * m_cosV * m_cosW * m_sinV * m_sinW -
210  2. * sigmaU2 * sigmaW2 * sigmaFit2 * m_cosV * m_cosW * m_sinV * m_sinW -
211  2. * sigmaU2 * sigmaFit2 * sigmaFit2 * m_cosV * m_cosW * m_sinV * m_sinW +
212  sigmaU2 * sigmaV2 * sigmaW2 * m_cosU * m_cosU * m_sinW * m_sinW + sigmaU2 * sigmaV2 * sigmaFit2 * m_cosU * m_cosU * m_sinW * m_sinW +
213  sigmaV2 * sigmaW2 * sigmaFit2 * m_cosU * m_cosU * m_sinW * m_sinW + sigmaV2 * sigmaFit2 * sigmaFit2 * m_cosU * m_cosU * m_sinW * m_sinW +
214  sigmaU2 * sigmaV2 * sigmaW2 * m_cosV * m_cosV * m_sinW * m_sinW + sigmaU2 * sigmaV2 * sigmaFit2 * m_cosV * m_cosV * m_sinW * m_sinW +
215  sigmaU2 * sigmaW2 * sigmaFit2 * m_cosV * m_cosV * m_sinW * m_sinW + sigmaU2 * sigmaFit2 * sigmaFit2 * m_cosV * m_cosV * m_sinW * m_sinW);
216 
217  z = (vFit * sigmaU2 * sigmaV2 * sigmaW2 * m_cosV * m_sinU * m_sinU + vFit * sigmaV2 * sigmaW2 * sigmaFit2 * m_cosV * m_sinU * m_sinU +
218  sigmaU2 * sigmaW2 * sigmaFit2 * v * m_cosV * m_sinU * m_sinU + sigmaW2 * sigmaFit2 * sigmaFit2 * v * m_cosV * m_sinU * m_sinU +
219  wFit * sigmaU2 * sigmaV2 * sigmaW2 * m_cosW * m_sinU * m_sinU + wFit * sigmaV2 * sigmaW2 * sigmaFit2 * m_cosW * m_sinU * m_sinU +
220  sigmaU2 * sigmaV2 * sigmaFit2 * w * m_cosW * m_sinU * m_sinU + sigmaV2 * sigmaFit2 * sigmaFit2 * w * m_cosW * m_sinU * m_sinU -
221  vFit * sigmaU2 * sigmaV2 * sigmaW2 * m_cosU * m_sinU * m_sinV - vFit * sigmaV2 * sigmaW2 * sigmaFit2 * m_cosU * m_sinU * m_sinV -
222  sigmaU2 * sigmaW2 * sigmaFit2 * v * m_cosU * m_sinU * m_sinV - sigmaW2 * sigmaFit2 * sigmaFit2 * v * m_cosU * m_sinU * m_sinV -
223  uFit * sigmaU2 * sigmaV2 * sigmaW2 * m_cosV * m_sinU * m_sinV - uFit * sigmaU2 * sigmaW2 * sigmaFit2 * m_cosV * m_sinU * m_sinV -
224  sigmaV2 * sigmaW2 * sigmaFit2 * u * m_cosV * m_sinU * m_sinV - sigmaW2 * sigmaFit2 * sigmaFit2 * u * m_cosV * m_sinU * m_sinV +
225  uFit * sigmaU2 * sigmaV2 * sigmaW2 * m_cosU * m_sinV * m_sinV + uFit * sigmaU2 * sigmaW2 * sigmaFit2 * m_cosU * m_sinV * m_sinV +
226  sigmaV2 * sigmaW2 * sigmaFit2 * u * m_cosU * m_sinV * m_sinV + sigmaW2 * sigmaFit2 * sigmaFit2 * u * m_cosU * m_sinV * m_sinV +
227  wFit * sigmaU2 * sigmaV2 * sigmaW2 * m_cosW * m_sinV * m_sinV + wFit * sigmaU2 * sigmaW2 * sigmaFit2 * m_cosW * m_sinV * m_sinV +
228  sigmaU2 * sigmaV2 * sigmaFit2 * w * m_cosW * m_sinV * m_sinV + sigmaU2 * sigmaFit2 * sigmaFit2 * w * m_cosW * m_sinV * m_sinV -
229  wFit * sigmaU2 * sigmaV2 * sigmaW2 * m_cosU * m_sinU * m_sinW - wFit * sigmaV2 * sigmaW2 * sigmaFit2 * m_cosU * m_sinU * m_sinW -
230  sigmaU2 * sigmaV2 * sigmaFit2 * w * m_cosU * m_sinU * m_sinW - sigmaV2 * sigmaFit2 * sigmaFit2 * w * m_cosU * m_sinU * m_sinW -
231  uFit * sigmaU2 * sigmaV2 * sigmaW2 * m_cosW * m_sinU * m_sinW - uFit * sigmaU2 * sigmaV2 * sigmaFit2 * m_cosW * m_sinU * m_sinW -
232  sigmaV2 * sigmaW2 * sigmaFit2 * u * m_cosW * m_sinU * m_sinW - sigmaV2 * sigmaFit2 * sigmaFit2 * u * m_cosW * m_sinU * m_sinW -
233  wFit * sigmaU2 * sigmaV2 * sigmaW2 * m_cosV * m_sinV * m_sinW - wFit * sigmaU2 * sigmaW2 * sigmaFit2 * m_cosV * m_sinV * m_sinW -
234  sigmaU2 * sigmaV2 * sigmaFit2 * w * m_cosV * m_sinV * m_sinW - sigmaU2 * sigmaFit2 * sigmaFit2 * w * m_cosV * m_sinV * m_sinW -
235  vFit * sigmaU2 * sigmaV2 * sigmaW2 * m_cosW * m_sinV * m_sinW - vFit * sigmaU2 * sigmaV2 * sigmaFit2 * m_cosW * m_sinV * m_sinW -
236  sigmaU2 * sigmaW2 * sigmaFit2 * v * m_cosW * m_sinV * m_sinW - sigmaU2 * sigmaFit2 * sigmaFit2 * v * m_cosW * m_sinV * m_sinW +
237  uFit * sigmaU2 * sigmaV2 * sigmaW2 * m_cosU * m_sinW * m_sinW + uFit * sigmaU2 * sigmaV2 * sigmaFit2 * m_cosU * m_sinW * m_sinW +
238  sigmaV2 * sigmaW2 * sigmaFit2 * u * m_cosU * m_sinW * m_sinW + sigmaV2 * sigmaFit2 * sigmaFit2 * u * m_cosU * m_sinW * m_sinW +
239  vFit * sigmaU2 * sigmaV2 * sigmaW2 * m_cosV * m_sinW * m_sinW + vFit * sigmaU2 * sigmaV2 * sigmaFit2 * m_cosV * m_sinW * m_sinW +
240  sigmaU2 * sigmaW2 * sigmaFit2 * v * m_cosV * m_sinW * m_sinW + sigmaU2 * sigmaFit2 * sigmaFit2 * v * m_cosV * m_sinW * m_sinW) /
241  (sigmaU2 * sigmaV2 * sigmaW2 * m_cosV * m_cosV * m_sinU * m_sinU + sigmaU2 * sigmaW2 * sigmaFit2 * m_cosV * m_cosV * m_sinU * m_sinU +
242  sigmaV2 * sigmaW2 * sigmaFit2 * m_cosV * m_cosV * m_sinU * m_sinU + sigmaW2 * sigmaFit2 * sigmaFit2 * m_cosV * m_cosV * m_sinU * m_sinU +
243  sigmaU2 * sigmaV2 * sigmaW2 * m_cosW * m_cosW * m_sinU * m_sinU + sigmaU2 * sigmaV2 * sigmaFit2 * m_cosW * m_cosW * m_sinU * m_sinU +
244  sigmaV2 * sigmaW2 * sigmaFit2 * m_cosW * m_cosW * m_sinU * m_sinU + sigmaV2 * sigmaFit2 * sigmaFit2 * m_cosW * m_cosW * m_sinU * m_sinU -
245  2. * sigmaU2 * sigmaV2 * sigmaW2 * m_cosU * m_cosV * m_sinU * m_sinV - 2. * sigmaU2 * sigmaW2 * sigmaFit2 * m_cosU * m_cosV * m_sinU * m_sinV -
246  2. * sigmaV2 * sigmaW2 * sigmaFit2 * m_cosU * m_cosV * m_sinU * m_sinV -
247  2. * sigmaW2 * sigmaFit2 * sigmaFit2 * m_cosU * m_cosV * m_sinU * m_sinV + sigmaU2 * sigmaV2 * sigmaW2 * m_cosU * m_cosU * m_sinV * m_sinV +
248  sigmaU2 * sigmaW2 * sigmaFit2 * m_cosU * m_cosU * m_sinV * m_sinV + sigmaV2 * sigmaW2 * sigmaFit2 * m_cosU * m_cosU * m_sinV * m_sinV +
249  sigmaW2 * sigmaFit2 * sigmaFit2 * m_cosU * m_cosU * m_sinV * m_sinV + sigmaU2 * sigmaV2 * sigmaW2 * m_cosW * m_cosW * m_sinV * m_sinV +
250  sigmaU2 * sigmaV2 * sigmaFit2 * m_cosW * m_cosW * m_sinV * m_sinV + sigmaU2 * sigmaW2 * sigmaFit2 * m_cosW * m_cosW * m_sinV * m_sinV +
251  sigmaU2 * sigmaFit2 * sigmaFit2 * m_cosW * m_cosW * m_sinV * m_sinV - 2. * sigmaU2 * sigmaV2 * sigmaW2 * m_cosU * m_cosW * m_sinU * m_sinW -
252  2. * sigmaU2 * sigmaV2 * sigmaFit2 * m_cosU * m_cosW * m_sinU * m_sinW - 2. * sigmaV2 * sigmaW2 * sigmaFit2 * m_cosU * m_cosW * m_sinU * m_sinW -
253  2. * sigmaV2 * sigmaFit2 * sigmaFit2 * m_cosU * m_cosW * m_sinU * m_sinW -
254  2. * sigmaU2 * sigmaV2 * sigmaW2 * m_cosV * m_cosW * m_sinV * m_sinW - 2. * sigmaU2 * sigmaV2 * sigmaFit2 * m_cosV * m_cosW * m_sinV * m_sinW -
255  2. * sigmaU2 * sigmaW2 * sigmaFit2 * m_cosV * m_cosW * m_sinV * m_sinW -
256  2. * sigmaU2 * sigmaFit2 * sigmaFit2 * m_cosV * m_cosW * m_sinV * m_sinW +
257  sigmaU2 * sigmaV2 * sigmaW2 * m_cosU * m_cosU * m_sinW * m_sinW + sigmaU2 * sigmaV2 * sigmaFit2 * m_cosU * m_cosU * m_sinW * m_sinW +
258  sigmaV2 * sigmaW2 * sigmaFit2 * m_cosU * m_cosU * m_sinW * m_sinW + sigmaV2 * sigmaFit2 * sigmaFit2 * m_cosU * m_cosU * m_sinW * m_sinW +
259  sigmaU2 * sigmaV2 * sigmaW2 * m_cosV * m_cosV * m_sinW * m_sinW + sigmaU2 * sigmaV2 * sigmaFit2 * m_cosV * m_cosV * m_sinW * m_sinW +
260  sigmaU2 * sigmaW2 * sigmaFit2 * m_cosV * m_cosV * m_sinW * m_sinW + sigmaU2 * sigmaFit2 * sigmaFit2 * m_cosV * m_cosV * m_sinW * m_sinW);
261 
262  const double outputU(LArRotationalTransformationPlugin::YZtoU(y, z));
263  const double outputV(LArRotationalTransformationPlugin::YZtoV(y, z));
264  const double outputW(LArRotationalTransformationPlugin::YZtoW(y, z));
265 
266  const double deltaU(u - outputU), deltaV(v - outputV), deltaW(w - outputW);
267  const double deltaUFit(uFit - outputU), deltaVFit(vFit - outputV), deltaWFit(wFit - outputW);
268 
269  chiSquared = ((deltaU * deltaU) / sigmaU2) + ((deltaV * deltaV) / sigmaV2) + ((deltaW * deltaW) / sigmaW2) +
270  ((deltaUFit * deltaUFit) / sigmaFit2) + ((deltaVFit * deltaVFit) / sigmaFit2) + ((deltaWFit * deltaWFit) / sigmaFit2);
271 }
272 
273 //------------------------------------------------------------------------------------------------------------------------------------------
274 
276 {
277  const LArTPCMap &larTPCMap(this->GetPandora().GetGeometry()->GetLArTPCMap());
278 
279  if (larTPCMap.empty())
280  {
281  std::cout << "LArRotationalTransformationPlugin::Initialize - LArTPC description not registered with Pandora as required " << std::endl;
282  return STATUS_CODE_NOT_INITIALIZED;
283  }
284 
285  const LArTPC *const pFirstLArTPC(larTPCMap.begin()->second);
286  m_thetaU = pFirstLArTPC->GetWireAngleU();
287  m_thetaV = pFirstLArTPC->GetWireAngleV();
288  m_thetaW = pFirstLArTPC->GetWireAngleW();
289  const double sigmaUVW(pFirstLArTPC->GetSigmaUVW());
290 
291  m_sinU = std::sin(m_thetaU);
292  m_sinV = std::sin(m_thetaV);
293  m_sinW = std::sin(m_thetaW);
294  m_cosU = std::cos(m_thetaU);
295  m_cosV = std::cos(m_thetaV);
296  m_cosW = std::cos(m_thetaW);
297 
298  m_sinVminusU = std::sin(m_thetaV - m_thetaU);
299  m_sinWminusV = std::sin(m_thetaW - m_thetaV);
300  m_sinUminusW = std::sin(m_thetaU - m_thetaW);
301 
302  if ((std::fabs(m_sinVminusU) < std::numeric_limits<double>::epsilon()) ||
303  (std::fabs(m_sinWminusV) < std::numeric_limits<double>::epsilon()) || (std::fabs(m_sinUminusW) < std::numeric_limits<double>::epsilon()))
304  {
305  std::cout << "LArRotationalTransformationPlugin::Initialize - Equal wire angles; Plugin does not support provided LArTPC configurations. "
306  << std::endl;
307  return STATUS_CODE_INVALID_PARAMETER;
308  }
309 
310  for (const LArTPCMap::value_type &mapEntry : larTPCMap)
311  {
312  const LArTPC *const pLArTPC(mapEntry.second);
313 
314  if ((std::fabs(m_thetaU - pLArTPC->GetWireAngleU()) > m_maxAngularDiscrepancyU) ||
315  (std::fabs(m_thetaV - pLArTPC->GetWireAngleV()) > m_maxAngularDiscrepancyV) ||
316  (std::fabs(m_thetaW - pLArTPC->GetWireAngleW()) > m_maxAngularDiscrepancyW) ||
317  (std::fabs(sigmaUVW - pLArTPC->GetSigmaUVW()) > m_maxSigmaDiscrepancy))
318  {
319  std::cout << "LArRotationalTransformationPlugin::Initialize - Dissimilar drift volumes; Plugin does not support provided LArTPC configurations. "
320  << std::endl;
321  return STATUS_CODE_INVALID_PARAMETER;
322  }
323  }
324 
325  return STATUS_CODE_SUCCESS;
326 }
327 
328 //------------------------------------------------------------------------------------------------------------------------------------------
329 
330 pandora::StatusCode LArRotationalTransformationPlugin::ReadSettings(const pandora::TiXmlHandle xmlHandle)
331 {
332  PANDORA_RETURN_RESULT_IF_AND_IF(
333  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MaxAngularDiscrepancyU", m_maxAngularDiscrepancyU));
334 
335  PANDORA_RETURN_RESULT_IF_AND_IF(
336  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MaxAngularDiscrepancyV", m_maxAngularDiscrepancyV));
337 
338  PANDORA_RETURN_RESULT_IF_AND_IF(
339  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MaxAngularDiscrepancyW", m_maxAngularDiscrepancyW));
340 
341  PANDORA_RETURN_RESULT_IF_AND_IF(
342  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MaxSigmaDiscrepancy", m_maxSigmaDiscrepancy));
343 
344  return STATUS_CODE_SUCCESS;
345 }
346 
347 } // namespace lar_content
virtual void GetMinChiSquaredYZ(const double u, const double v, const double w, const double sigmaU, const double sigmaV, const double sigmaW, double &y, double &z, double &chiSquared) const
virtual double UVtoY(const double u, const double v) const
virtual double UVtoW(const double u, const double v) const
virtual double YZtoW(const double y, const double z) const
virtual double UVtoZ(const double u, const double v) const
virtual double VWtoU(const double v, const double w) const
double m_maxAngularDiscrepancyU
Maximum allowed difference between u wire angles between LArTPCs.
virtual double UWtoZ(const double u, const double w) const
Header file for the rotational transformation plugin class.
virtual double YZtoV(const double y, const double z) const
double m_maxAngularDiscrepancyV
Maximum allowed difference between v wire angles between LArTPCs.
double m_maxAngularDiscrepancyW
Maximum allowed difference between w wire angles between LArTPCs.
pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
virtual double VWtoZ(const double v, const double w) const
virtual double WUtoV(const double w, const double u) const
virtual double YZtoU(const double y, const double z) const
double m_maxSigmaDiscrepancy
Maximum allowed difference between like wire sigma values between LArTPCs.
virtual double VWtoY(const double v, const double w) const
QTextStream & endl(QTextStream &s)
GeomAnalyzerI * GetGeometry(void)
Definition: gAtmoEvGen.cxx:433
virtual double UWtoY(const double u, const double w) const