EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
FairGeaneUtil.cxx
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file FairGeaneUtil.cxx
1 // ------------------------------------------------------------------
2 // Version of June 10th
3 // modified to work properly in q/p variables instead of 1/p
4 // ------------------------------------------------------------------
5 
6 #include <iostream>
7 #include "FairGeaneUtil.h"
8 #include "TGeoTorus.h"
9 #include "TMath.h"
10 
11 using namespace std;
12 
13 FairGeaneUtil::FairGeaneUtil() : TObject() { }
14 
15 FairGeaneUtil::~FairGeaneUtil() { }
16 
17 void FairGeaneUtil::FromPtToSC(Double_t PC[3], Double_t RC[15],
18  // output
19  Double_t* PD, Double_t* RD, Int_t& IERR)
20 {
21 
22 // -------------------------------------------------------------
23 //
24 // TRANSFORM ERROR MATRIX
25 //
26 // FROM SC VARIABLES (q/Pt,LAMBDA,PHI,YT,ZT)
27 // FROM SC VARIABLES (q/P,LAMBDA,PHI,YT,ZT)
28 //
29 // INPUT
30 // CH Charge of the paticle
31 // PC[3] (q/Pt, Lambda, Phi, Yt, Zt)
32 // RC[15] error matrix in upper triangular form
33 //
34 // OUTPUT
35 // PD[3] (q/P, Lambda, Phi, Yt, Zt)
36 // RD[15] error matrix in upper triangular form
37 // IERR =1 when track is perp to X-axis of Master
38 //
39 // Author EMC Collaboration
40 // Translated in C/Rot by A. Rotondi (June 2007)
41 //
42 // -------------------------------------------------------------------
43 
44  Double_t A[5][5], S[15], COSL, SINL;
45  Double_t Vec[25];
46 
47  IERR = 0;
48  memset(RD,0,sizeof(*RD));
49 
50  COSL = cos(PC[1]);
51 
52  if(TMath::Abs(COSL) < 1.e-7) {
53  IERR =1;
54  return;
55  }
56 
57  SINL = sin(PC[1]);
58 
59  PD[0] = PC[0]*COSL;
60  PD[1] = PC[1];
61  PD[2] = PC[2];
62 
63 
64  for(Int_t I=0; I<5; I++) {
65  for(Int_t K=0; K<5; K++) {
66  A[I][K]=0.;
67  }
68  }
69 
70  // copy input in an internal vector S
71  for(Int_t I=0; I<15; I++) { S[I]=RC[I]; }
72 
73  A[0][0] = COSL;
74  A[1][1] = 1.0;
75  A[2][2] = 1.0;
76  A[3][3] = 1.0;
77  A[4][4] = 1.0;
78 
79  A[0][1] = -PC[0]*SINL;
80 
81  // transformation
82  FromMatToVec(A,Vec);
83  SymmProd(Vec,S,S);
84 
85  // copy the result in S in the output vector
86  for(Int_t I=0; I<15; I++) { RD[I]=S[I]; }
87 }
88 
89 void FairGeaneUtil::FromPtToSD(Double_t PD[3], Double_t RD[15], Double_t H[3], Int_t CH,
90  Double_t SPU, Double_t DJ[2], Double_t DK[2],
91  // output
92  Int_t& IERR, Double_t* PC, Double_t* RC)
93 {
94 
95 //
96 // **************************************************************************
97 //
98 // *** TRANSFORMS ERROR MATRIX
99 // FROM VARIABLES (q/Pt,V',W',V,W)
100 // FROM VARIABLES (q/P, V',W',V,W)
101 //
102 // input
103 // PD[3] (q/P, Lambda, Phi, Yt, Zt)
104 // RD[15] error matrix in upper triangular form
105 // H[3] magnetic field
106 // CH CHARGE OF PARTICLE
107 // CHARGE AND MAGNETIC FIELD ARE NEEDED
108 // FOR CORRELATION TERMS (V',YT),(V',ZT),(W',YT),(W',ZT)
109 // THESE CORRELATION TERMS APPEAR BECAUSE RC IS ASSUMED
110 // TO BE THE ERROR MATRIX FOR FIXED S (PATH LENGTH)
111 // AND RD FOR FIXED U
112 // SPU SIGN OF U-COMPONENT OF PARTICLE MOMENTUM
113 // spu = sign[p·(DJ x DK)]
114 // DJ[2] UNIT VECTOR IN V-DIRECTION
115 // DK[2] UNIT VECTOR IN W-DIRECTION OF DETECTOR SYSTEM
116 //
117 // output
118 // PC[3] (q/Pt, Lambda, Phi, Yt, Zt)
119 // RC[15] error matrix in upper triangular form
120 // IERR =1 when track is perp to X-axis of Master
121 //
122 // Author EMC Collaboration
123 // Translated in C/Rot by A. Rotondi (June 2007)
124 //
125 // -------------------------------------------------------------------
126 //
127 
128 
129 
130  Double_t A[5][5], S[15], TN[3], COSL, SINL, COSL1;
131 
132  Double_t PM, HM, UN[3], VN[3], DI[3], TVW[3];
133  Double_t UJ, UK, VJ, VK, HA, HAM, Q, SINZ, COSZ;
134 
135  Double_t CFACT8 = 2.997925e-04;
136 
137  Double_t Vec[25];
138 
139 
140  // ------------------------------------------------------------------
141 
142  IERR=0;
143  TVW[0]=1./TMath::Sqrt(1.+PD[1]*PD[1]+PD[2]*PD[2]);
144  if(SPU < 0.) { TVW[0]=-TVW[0]; }
145  TVW[1]=PD[1]*TVW[0];
146  TVW[2]=PD[2]*TVW[0];
147 
148  DI[0]=DJ[1]*DK[2]-DJ[2]*DK[1];
149  DI[1]=DJ[2]*DK[0]-DJ[0]*DK[2];
150  DI[2]=DJ[0]*DK[1]-DJ[1]*DK[0];
151 
152  for(Int_t I=0; I<3; I++) {
153  TN[I]=TVW[0]*DI[I]+TVW[1]*DJ[I]+TVW[2]*DK[I];
154  }
155 
156  COSL=TMath::Sqrt(TMath::Abs(1.-TN[2]*TN[2]));
157  if(COSL < 1.e-30) { COSL = 1.e-30; }
158  COSL1=1./COSL;
159  SINL = TN[2];
160 
161  PC[0]=PD[0]*COSL;
162  PC[1]=PD[1];
163  PC[2]=PD[2];
164  PM=PC[0];
165 
166  if(TMath::Abs (TN[0]) < 1.E-30) { TN[0] = 1.e-30; }
167 
168  UN[0]=-TN[1]*COSL1;
169  UN[1]=TN[0]*COSL1;
170  UN[2]=0.;
171 
172  VN[0]=-TN[2]*UN[1];
173  VN[1]=TN[2]*UN[0];
174  VN[2]=COSL;
175 
176  UJ=UN[0]*DJ[0]+UN[1]*DJ[1]+UN[2]*DJ[2];
177  UK=UN[0]*DK[0]+UN[1]*DK[1]+UN[2]*DK[2];
178  VJ=VN[0]*DJ[0]+VN[1]*DJ[1]+VN[2]*DJ[2];
179  VK=VN[0]*DK[0]+VN[1]*DK[1]+VN[2]*DK[2];
180 
181 
182  for(Int_t I=0; I<5; I++) {
183  for(Int_t K=0; K<5; K++) {
184  A[I][K]=0.;
185  }
186  }
187 
188  // copy input in an internal vector S
189  for(Int_t I=0; I<15; I++) { S[I]=RD[I]; }
190 
191  if(CH != 0.) {
192  HA=TMath::Sqrt(H[0]*H[0]+H[1]*H[1]+H[2]*H[2]);
193  HAM=HA*PM;
194  if(HAM != 0.) {
195  HM=CH/HA;
196 
197  Q=-HAM*CFACT8;
198 
199  SINZ=-(H[0]*UN[0]+H[1]*UN[1]+H[2]*UN[2])*HM;
200  COSZ= (H[0]*VN[0]+H[1]*VN[1]+H[2]*VN[2])*HM;
201  A[0][3] = Q*TVW[1]*SINZ*(SINL*PD[0]);
202  A[0][4] = Q*TVW[2]*SINZ*(SINL*PD[0]);
203  }
204  }
205 
206  A[0][0] = COSL;
207  A[1][1] = 1.;
208  A[2][2] = 1.;
209  A[3][3] = 1.;
210  A[4][4] = 1.;
211 
212  A[0][1] = -TVW[0]*VJ*(SINL*PD[0]);
213  A[0][2] = -TVW[0]*VK*(SINL*PD[0]);
214 
215 
216  // transformation
217  FromMatToVec(A,Vec);
218  SymmProd(Vec,S,S);
219 
220  // copy the result in S in the output vector
221  for(Int_t I=0; I<15; I++) { RC[I]=S[I]; }
222 }
223 
224 
225 void FairGeaneUtil::FromSCToPt(Double_t PC[3], Double_t RC[15],
226  // output
227  Double_t* PD, Double_t* RD, Int_t& IERR)
228 {
229 
230 // -------------------------------------------------------------
231 //
232 // TRANSFORM ERROR MATRIX
233 //
234 // FROM SC VARIABLES (q/P,LAMBDA,PHI,YT,ZT)
235 // FROM SC VARIABLES (q/Pt,LAMBDA,PHI,YT,ZT)
236 //
237 // INPUT
238 // PC[3] (q/P, Lambda, Phi, Yt, Zt)
239 // RC[15] error matrix in upper triangular form
240 //
241 // OUTPUT
242 // PD[3] (q/Pt, Lambda, Phi, Yt, Zt)
243 // RD[15] error matrix in upper triangular form
244 // IERR =1 when track is perp to X-axis of Master
245 //
246 // Author EMC Collaboration
247 // Translated in C/Rot by A. Rotondi (June 2007)
248 //
249 // -------------------------------------------------------------------
250 
251 
252  Double_t A[5][5], S[15], COSL, COSL1;
253  Double_t TANL;
254  Double_t Vec[25];
255 
256  IERR = 0;
257  memset(RD,0,sizeof(*RD));
258 
259  COSL = cos(PC[1]);
260  if(TMath::Abs(COSL) < 1.e-7) {
261  IERR = 1;
262  return;
263  }
264  COSL1 = 1./COSL;
265  TANL = tan(PC[1]);
266 
267  PD[0] = PC[0]*COSL1;
268  PD[1] = PC[1];
269  PD[2] = PC[2];
270 
271  for(Int_t I=0; I<5; I++) {
272  for(Int_t K=0; K<5; K++) {
273  A[I][K]=0.;
274  }
275  }
276 
277  // copy input in an internal vector S
278  for(Int_t I=0; I<15; I++) { S[I]=RC[I]; }
279 
280  A[0][0] = COSL1;
281  A[1][1] = 1.0;
282  A[2][2] = 1.0;
283  A[3][3] = 1.0;
284  A[4][4] = 1.0;
285 
286  A[0][1] = PD[0]*TANL;
287 
288  // transformation
289  FromMatToVec(A,Vec);
290  SymmProd(Vec,S,S);
291 
292  // copy the result in S in the output vector
293  for(Int_t I=0; I<15; I++) { RD[I]=S[I]; }
294 
295 }
296 
297 void FairGeaneUtil::FromSCToSD(Double_t PC[3], Double_t RC[15], Double_t H[3], Int_t CH,
298  Double_t DJ[3], Double_t DK[3],
299  // output
300  Int_t& IERR, Double_t& SPU,
301  Double_t* PD, Double_t* RD)
302 {
303 
304 // ----------------------------------------------------------------------
305 // Transform Error Matrix
306 // FROM SC (transverse system) VARIABLES (q/P,LAMBDA,PHI,YT,ZT)
307 // TO SD (detector system) VARIABLES (q/P,V',W',V,W)
308 
309 // Authors: A. Haas and W. Wittek
310 // Translated in CRoot by A. Rotondi and A. Fontana June 2007
311 
312 // INPUT
313 // PC(3) q/P,LAMBDA,PHI
314 // H(3) MAGNETIC FIELD
315 // RC(15) ERROR MATRIX IN SC VARIABLES (TRIANGLE)
316 // CH CHARGE OF PARTICLE
317 // CHARGE AND MAGNETIC FIELD ARE NEEDED
318 // FOR CORRELATION TERMS (V',YT),(V',ZT),(W',YT),(W',ZT)
319 // THESE CORRELATION TERMS APPEAR BECAUSE RC IS ASSUMED
320 // TO BE THE ERROR MATRIX FOR FIXED S (PATH LENGTH)
321 // AND RD FOR FIXED U
322 // DJ(3) UNIT VECTOR IN V-DIRECTION
323 // DK(3) UNIT VECTOR IN W-DIRECTION OF DETECTOR SYSTEM
324 
325 // OUTPUT
326 // RD(15) ERROR MATRIX IN q/P,V',W',V,W (TRIANGLE)
327 // PD(3) q/P,V',W'
328 // IERR = 1 PARTICLE MOVES PERPENDICULAR TO U-AXIS
329 // ( V',W' ARE NOT DEFINED )
330 // SPU SIGN OF U-COMPONENT OF PARTICLE MOMENTUM
331 // spu = sign[p·(DJ x DK)]
332 //
333 // ------------------------------------------------------------------------
334 
335  Double_t A[5][5], S[15], TN[3], COSL, COSP;
336  Double_t UN[3], VN[3], DI[3], TVW[3];
337 
338  Double_t Vec[25];
339  Double_t CFACT8= 2.997925e-04;
340  Double_t T1R, T2R, T3R, SINP, SINZ, COSZ, HA, HM, HAM;
341  Double_t Q, UI, VI, UJ, UK, VJ, VK;
342 
343  IERR=0;
344  memset(RD,0,sizeof(*RD));
345  memset(PD,0,sizeof(*PD));
346 
347  COSL=TMath::Cos(PC[1]);
348  SINP=TMath::Sin(PC[2]);
349  COSP=TMath::Cos(PC[2]);
350 
351  TN[0]=COSL*COSP;
352  TN[1]=COSL*SINP;
353  TN[2]=TMath::Sin(PC[1]);
354 
355  // DI = DJ x DK
356  DI[0]=DJ[1]*DK[2]-DJ[2]*DK[1];
357  DI[1]=DJ[2]*DK[0]-DJ[0]*DK[2];
358  DI[2]=DJ[0]*DK[1]-DJ[1]*DK[0];
359 
360  TVW[0]=TN[0]*DI[0]+TN[1]*DI[1]+TN[2]*DI[2];
361  SPU=1.;
362  if(TVW[0] < 0.) { SPU=-1.; }
363  TVW[1]=TN[0]*DJ[0]+TN[1]*DJ[1]+TN[2]*DJ[2];
364  TVW[2]=TN[0]*DK[0]+TN[1]*DK[1]+TN[2]*DK[2];
365 
366  // track lies in the detector plane: stop calculations
367  if(TMath::Abs(TVW[0]) < 1.e-7) {
368  IERR = 1;
369  return;
370  }
371 
372  T1R=1./TVW[0];
373  PD[0] = PC[0];
374  PD[1] = TVW[1]*T1R;
375  PD[2] = TVW[2]*T1R;
376 
377  UN[0] = -SINP;
378  UN[1] = COSP;
379  UN[2] = 0.;
380 
381  VN[0] =-TN[2]*UN[1];
382  VN[1] = TN[2]*UN[0];
383  VN[2] = COSL;
384 
385  UJ=UN[0]*DJ[0]+UN[1]*DJ[1]+UN[2]*DJ[2];
386  UK=UN[0]*DK[0]+UN[1]*DK[1]+UN[2]*DK[2];
387  VJ=VN[0]*DJ[0]+VN[1]*DJ[1]+VN[2]*DJ[2];
388  VK=VN[0]*DK[0]+VN[1]*DK[1]+VN[2]*DK[2];
389 
390  for(Int_t I=0; I<5; I++) {
391  for(Int_t K=0; K<5; K++) {
392  A[I][K]=0.;
393  }
394  }
395 
396  // copy input in an internal vector S
397  for(Int_t I=0; I<15; I++) { S[I]=RC[I]; }
398 
399  if(CH != 0.) {
400  //charged particles
401  HA=TMath::Sqrt(H[0]*H[0]+H[1]*H[1]+H[2]*H[2]);
402  HAM=HA*PC[0];
403  if(HAM != 0.) {
404  // ... in a magnetic field
405  HM=CH/HA;
406  Q=-HAM*CFACT8;
407  SINZ=-(H[0]*UN[0]+H[1]*UN[1]+H[2]*UN[2])*HM;
408  COSZ= (H[0]*VN[0]+H[1]*VN[1]+H[2]*VN[2])*HM;
409  T3R=Q*pow(T1R,3);
410  UI=UN[0]*DI[0]+UN[1]*DI[1]+UN[2]*DI[2];
411  VI=VN[0]*DI[0]+VN[1]*DI[1]+VN[2]*DI[2];
412  A[1][3] = -UI*(VK*COSZ-UK*SINZ)*T3R;
413  A[1][4] = -VI*(VK*COSZ-UK*SINZ)*T3R;
414  A[2][3] = UI*(VJ*COSZ-UJ*SINZ)*T3R;
415  A[2][4] = VI*(VJ*COSZ-UJ*SINZ)*T3R;
416  }
417  }
418 
419  T2R=T1R*T1R;
420 
421  // Transformation matrix from SC to SD
422 
423  A[0][0] = 1.;
424  A[1][1] = -UK*T2R;
425  A[1][2] = VK*COSL*T2R;
426  A[2][1] = UJ*T2R;
427  A[2][2] = -VJ*COSL*T2R;
428  A[3][3] = VK*T1R;
429  A[3][4] = -UK*T1R;
430  A[4][3] = -VJ*T1R;
431  A[4][4] = UJ*T1R;
432 
433  // transformation
434  FromMatToVec(A,Vec);
435  SymmProd(Vec,S,S);
436 
437  // copy the result in S in the output vector
438  for(Int_t I=0; I<15; I++) { RD[I]=S[I]; }
439 
440 }
441 
442 
443 void FairGeaneUtil::FromSD1ToSD2(Double_t PD1[2], Double_t RD1[15],Double_t H[2],
444  Int_t CH, Double_t SP1,
445  Double_t DJ1[2], Double_t DK1[2],
446  Double_t DJ2[2], Double_t DK2[2],
447  // output
448  Int_t& IERR, Double_t& SP2,
449  Double_t* PD2, Double_t* RD2)
450 {
451 
452 // -------------------------------------------------------------------------
453 //
454 // TRANSFORMS ERROR MATRIX
455 // FROM VARIABLES (q/P,V1',W1',V1,W1)
456 // TO VARIABLES (q/P,V2',W2',V2,W2)
457 //
458 // Authors: A. Haas and W. Wittek
459 // Translated in C/Root by A. Fontana and A. Rotondi (June 2007)
460 //
461 //
462 // INPUT
463 // PD1[2] q/P,V1',W1'
464 // H[2] MAGNETIC FIELD
465 // RD1(15) ERROR MATRIX IN 1/P,V1',W1',V1,W1 (Triangular)
466 // CH CHARGE OF PARTICLE
467 // CHARGE AND MAGNETIC FIELD ARE NEEDED
468 // FOR CORRELATION TERMS (V2',V1),(V2',W1),(W2',V1),(W2',W1)
469 // THESE CORRELATION TERMS APPEAR BECAUSE RD1 IS ASSUMED
470 // TO BE THE ERROR MATRIX FOR FIXED U1
471 // AND RD2 FOR FIXED U2
472 // SP1 SIGN OF U1-COMPONENT OF PARTICLE MOMENTUM INPUT
473 // DJ1[2] UNIT VECTOR IN V1-DIRECTION
474 // DK1[2] UNIT VECTOR IN W1-DIRECTION OF SYSTEM 1
475 // DJ2[2] UNIT VECTOR IN V2-DIRECTION
476 // DK2[2] UNIT VECTOR IN W2-DIRECTION OF SYSTEM 2
477 //
478 //
479 // OUTPUT
480 // PD2[3] q/P,V2',W2'
481 // RD2[15] ERROR MATRIX IN 1/P,V2',W2',V2,W2 (Triangular)
482 // SP2 SIGN OF U2-COMPONENT OF PARTICLE MOMENTUM OUTPUT
483 // IERR = 0 TRANSFORMATION OK
484 // = 1 MOMENTUM PERPENDICULAR TO U2-DIRECTION
485 // (V2',W2' NOT DEFINed)
486 // = 2 MOMENTUM PERPENDICULAR TO X-AXIS
487 //
488 // ----------------------------------------------------------------------
489 
490  Double_t A[5][5], S[15], TN[3], COSL, COSL1;
491  Double_t SINZ, COSZ, UN[3], VN[3];
492 
493  Double_t Vec[25];
494  Double_t PM, TR, TS, TT, HA, HM, HAM, Q;
495  Double_t UJ1, UK1, UJ2, UK2, VJ1, VJ2, VK1, VK2;
496  Double_t SJ1I2, SK1I2, SK2U, SK2V, SJ2U, SJ2V;
497  Double_t DI1[3], DI2[3], TVW1[3], TVW2[3];
498  Double_t CFACT8= 2.997925e-04;
499 
500  IERR=0;
501  memset(PD2,0,sizeof(*PD2));
502  memset(RD2,0,sizeof(*RD2));
503 
504  PM=PD1[0];
505  TVW1[0]=1./sqrt(1.+PD1[1]*PD1[1]+PD1[2]*PD1[2]);
506  if(SP1 < 0.) { TVW1[0]=-TVW1[0]; }
507  TVW1[1]=PD1[1]*TVW1[0];
508  TVW1[2]=PD1[2]*TVW1[0];
509 
510  DI1[0]=DJ1[1]*DK1[2]-DJ1[2]*DK1[1];
511  DI1[1]=DJ1[2]*DK1[0]-DJ1[0]*DK1[2];
512  DI1[2]=DJ1[0]*DK1[1]-DJ1[1]*DK1[0];
513 
514  for(Int_t I=0; I<3; I++) {
515  TN[I]=TVW1[0]*DI1[I]+TVW1[1]*DJ1[I]+TVW1[2]*DK1[I];
516  }
517 
518  DI2[0]=DJ2[1]*DK2[2]-DJ2[2]*DK2[1];
519  DI2[1]=DJ2[2]*DK2[0]-DJ2[0]*DK2[2];
520  DI2[2]=DJ2[0]*DK2[1]-DJ2[1]*DK2[0];
521 
522  TVW2[0]=TN[0]*DI2[0]+TN[1]*DI2[1]+TN[2]*DI2[2];
523  TVW2[1]=TN[0]*DJ2[0]+TN[1]*DJ2[1]+TN[2]*DJ2[2];
524  TVW2[2]=TN[0]*DK2[0]+TN[1]*DK2[1]+TN[2]*DK2[2];
525 
526  if(TMath::Abs(TVW2[0]) < 1.e-7) {
527  // track lies in the v-w plane: stop calculations
528  IERR = 1;
529  return;
530  }
531  TR=1./TVW2[0];
532  SP2=1;
533  if(TVW2[0] < 0.) { SP2=-1; }
534  PD2[0]=PD1[0];
535  PD2[1]=TVW2[1]*TR;
536  PD2[2]=TVW2[2]*TR;
537 
538  COSL=sqrt(TMath::Abs(1.-TN[2]*TN[2]));
539  if(TMath::Abs(COSL) < 1.e-7) {
540  // track perp to X-axis of Master: stop calculations
541  IERR=2;
542  return;
543  }
544  COSL1=1./COSL;
545  UN[0]=-TN[1]*COSL1;
546  UN[1]=TN[0]*COSL1;
547  UN[2]=0.;
548 
549  VN[0]=-TN[2]*UN[1];
550  VN[1]=TN[2]*UN[0];
551  VN[2]=COSL;
552 
553  UJ1=UN[0]*DJ1[0]+UN[1]*DJ1[1]+UN[2]*DJ1[2];
554  UK1=UN[0]*DK1[0]+UN[1]*DK1[1]+UN[2]*DK1[2];
555  VJ1=VN[0]*DJ1[0]+VN[1]*DJ1[1]+VN[2]*DJ1[2];
556  VK1=VN[0]*DK1[0]+VN[1]*DK1[1]+VN[2]*DK1[2];
557 
558  UJ2=UN[0]*DJ2[0]+UN[1]*DJ2[1]+UN[2]*DJ2[2];
559  UK2=UN[0]*DK2[0]+UN[1]*DK2[1]+UN[2]*DK2[2];
560  VJ2=VN[0]*DJ2[0]+VN[1]*DJ2[1]+VN[2]*DJ2[2];
561  VK2=VN[0]*DK2[0]+VN[1]*DK2[1]+VN[2]*DK2[2];
562 
563  // reset working vectors and matrices
564  for(Int_t I=0; I<5; I++) {
565  for(Int_t K=0; K<5; K++) {
566  A[I][K]=0.;
567  }
568  }
569  for(Int_t J=0; J<15; J++) { S[J]=RD1[J]; }
570 
571  if(CH != 0.) {
572  // a charged particle
573  HA=sqrt(H[0]*H[0]+H[1]*H[1]+H[2]*H[2]);
574  HAM=HA*PM;
575  if(HAM != 0.) {
576  // ...in a magnetic field
577  HM=CH/HA;
578  Q=-HAM*CFACT8;
579  TT=-Q*pow(TR,3);
580  SJ1I2=DJ1[0]*DI2[0]+DJ1[1]*DI2[1]+DJ1[2]*DI2[2];
581  SK1I2=DK1[0]*DI2[0]+DK1[1]*DI2[1]+DK1[2]*DI2[2];
582  SK2U=DK2[0]*UN[0]+DK2[1]*UN[1]+DK2[2]*UN[2];
583  SK2V=DK2[0]*VN[0]+DK2[1]*VN[1]+DK2[2]*VN[2];
584  SJ2U=DJ2[0]*UN[0]+DJ2[1]*UN[1]+DJ2[2]*UN[2];
585  SJ2V=DJ2[0]*VN[0]+DJ2[1]*VN[1]+DJ2[2]*VN[2];
586 
587  SINZ=-(H[0]*UN[0]+H[1]*UN[1]+H[2]*UN[2])*HM;
588  COSZ= (H[0]*VN[0]+H[1]*VN[1]+H[2]*VN[2])*HM;
589  A[1][3] = -TT*SJ1I2*(SK2U*SINZ-SK2V*COSZ);
590  A[1][4] = -TT*SK1I2*(SK2U*SINZ-SK2V*COSZ);
591  A[2][3] = TT*SJ1I2*(SJ2U*SINZ-SJ2V*COSZ);
592  A[2][4] = TT*SK1I2*(SJ2U*SINZ-SJ2V*COSZ);
593 
594  }
595  }
596  A[0][0] = 1.;
597  A[3][3] = TR*(UJ1*VK2-VJ1*UK2);
598  A[3][4] = TR*(UK1*VK2-VK1*UK2);
599  A[4][3] = TR*(VJ1*UJ2-UJ1*VJ2);
600  A[4][4] = TR*(VK1*UJ2-UK1*VJ2);
601 
602  TS=TR*TVW1[0];
603  A[1][1] = A[3][3]*TS;
604  A[1][2] = A[3][4]*TS;
605  A[2][1] = A[4][3]*TS;
606  A[2][2] = A[4][4]*TS;
607 
608  // transformation A*SA
609  FromMatToVec(A,Vec);
610  SymmProd(Vec,S,S);
611 
612  // final error (covariance) matrix in upper-triangular form
613 
614  // std::cout<<"SD1toSD2: "<<PD2[0]<<","<<PD2[1]<<","<<PD2[2]<<","<<std::endl;
615 
616  for(Int_t J=0; J<15; J++) {
617  RD2[J]=S[J];
618  // std::cout<<S[J]<<" ";
619  }
620 
621 }
622 
623 void FairGeaneUtil::FromSDToPt(Double_t PD[3], Double_t RD[15], Double_t H[3],
624  Int_t CH, Double_t SPU, Double_t DJ[3], Double_t DK[3],
625  // output
626  Int_t& IERR, Double_t* PC, Double_t* RC)
627 {
628 //
629 // -------------------------------------------------------------
630 //
631 // TRANSFORM ERROR MATRIX
632 // FROM VARIABLES (q/P,V',W',V,W)
633 // FROM VARIABLES (q/Pt,V',W',V,W)
634 //
635 // input
636 // PD[3] (q/P, Lambda, Phi, Yt, Zt)
637 // RD[15] error matrix in upper triangular form
638 // H[3] magnetic field
639 // CH CHARGE OF PARTICLE
640 // CHARGE AND MAGNETIC FIELD ARE NEEDED
641 // FOR CORRELATION TERMS (V',YT),(V',ZT),(W',YT),(W',ZT)
642 // THESE CORRELATION TERMS APPEAR BECAUSE RC IS ASSUMED
643 // TO BE THE ERROR MATRIX FOR FIXED S (PATH LENGTH)
644 // AND RD FOR FIXED U
645 // SPU SIGN OF U-COMPONENT OF PARTICLE MOMENTUM
646 // spu = sign[p·(DJ x DK)]
647 // DJ(3) UNIT VECTOR IN V-DIRECTION
648 // DK(3) UNIT VECTOR IN W-DIRECTION OF DETECTOR SYSTEM
649 //
650 // output
651 // PC[3] (q/Pt, Lambda, Phi, Yt, Zt)
652 // RC[15] error matrix in upper triangular form
653 // IERR =1 when track is perp to X-axis of Master
654 //
655 // Author EMC Collaboration
656 // Translated in C/Rot by A. Rotondi (June 2007)
657 //
658 // -------------------------------------------------------------------
659 //
660 
661 //
662 //
663 
664 
665  Double_t A[5][5], S[15], TN[3], TANL, COSL, COSL1;
666 
667  Double_t PM, HM, UN[3], VN[3], DI[3], TVW[3];
668 
669  Double_t CFACT8 = 2.997925e-04;
670 
671  Double_t UJ, UK, VJ, VK, HA, HAM, Q, SINZ, COSZ;
672  Double_t Vec[25];
673 
674  // ------------------------------------------------------------------
675 
676  IERR=0;
677 
678  PM=PD[0];
679  TVW[0]=1./TMath::Sqrt(1.+PD[1]*PD[1]+PD[2]*PD[2]);
680  if(SPU<0.) { TVW[0]=-TVW[0]; }
681  TVW[1]=PD[1]*TVW[0];
682  TVW[2]=PD[2]*TVW[0];
683 
684  DI[0]=DJ[1]*DK[2]-DJ[2]*DK[1];
685  DI[1]=DJ[2]*DK[0]-DJ[0]*DK[2];
686  DI[2]=DJ[0]*DK[1]-DJ[1]*DK[0];
687 
688  for(Int_t I=0; I<3; I++) {
689  TN[I]=TVW[0]*DI[I]+TVW[1]*DJ[I]+TVW[2]*DK[I];
690  }
691 
692  COSL=TMath::Sqrt(TMath::Abs(1.-TN[2]*TN[2]));
693  if(COSL < 1.e-30) { COSL = 1.e-30; }
694  COSL1=1./COSL;
695  TANL = TN[2]*COSL1;
696 
697  PC[0]=PD[0]*COSL1;
698  PC[1]=PD[1];
699  PC[2]=PD[2];
700 
701  if(TMath::Abs (TN[0])< 1.e-30) { TN[0] = 1.e-30; }
702 
703  UN[0]=-TN[1]*COSL1;
704  UN[1]=TN[0]*COSL1;
705  UN[2]=0.;
706 
707  VN[0]=-TN[2]*UN[1];
708  VN[1]=TN[2]*UN[0];
709  VN[2]=COSL;
710 
711  UJ=UN[0]*DJ[0]+UN[1]*DJ[1]+UN[2]*DJ[2];
712  UK=UN[0]*DK[0]+UN[1]*DK[1]+UN[2]*DK[2];
713  VJ=VN[0]*DJ[0]+VN[1]*DJ[1]+VN[2]*DJ[2];
714  VK=VN[0]*DK[0]+VN[1]*DK[1]+VN[2]*DK[2];
715 
716  for(Int_t I=0; I<5; I++) {
717  for(Int_t K=0; K<5; K++) {
718  A[I][K]=0.;
719  }
720  }
721 
722  // copy input in an internal vector S
723  for(Int_t I=0; I<15; I++) { S[I]=RD[I]; }
724 
725 
726  if(CH != 0.) {
727  HA=TMath::Sqrt(H[0]*H[0]+H[1]*H[1]+H[2]*H[2]);
728  HAM=HA*PM;
729  if(HAM != 0.) {
730  HM=CH/HA;
731 
732  Q=-HAM*CFACT8;
733 
734  SINZ=-(H[0]*UN[0]+H[1]*UN[1]+H[2]*UN[2])*HM;
735  COSZ= (H[0]*VN[0]+H[1]*VN[1]+H[2]*VN[2])*HM;
736  A[0][3] = -Q*TVW[1]*SINZ*(TANL*PC[0]);
737  A[0][4] = -Q*TVW[2]*SINZ*(TANL*PC[0]);
738  }
739  }
740 
741  A[0][0] = COSL1;
742  A[1][1] = 1.;
743  A[2][2] = 1.;
744  A[3][3] = 1.;
745  A[4][4] = 1.;
746 
747  A[0][1] = TVW[0]*VJ*(TANL*PC[0]);
748  A[0][2] = TVW[0]*VK*(TANL*PC[0]);
749 
750  // transformation
751  FromMatToVec(A,Vec);
752  SymmProd(Vec,S,S);
753 
754  // copy the result in S in the output vector
755  for(Int_t I=0; I<15; I++) { RC[I]=S[I]; }
756 }
757 
758 
759 
760 void FairGeaneUtil::FromSDToSC(Double_t PD[3], Double_t RD[15], Double_t H[3], Int_t CH,
761  Double_t SPU, Double_t DJ[3], Double_t DK[3],
762  // output
763  Int_t& IERR, Double_t* PC, Double_t* RC)
764 {
765 
766 // ----------------------------------------------------------------------
767 //
768 // Tranform error matrix
769 // FROM SD (detector plane) VARIABLES (q/P,V',W',V,W)
770 // TO SC (transverse system) VARIABLES (q/P,LAMBDA,PHI,YT,ZT)
771 //
772 // Authors: A. Haas and W. Wittek
773 // Translated in C/Root by A. Rotondi and A. Fontana (June 2007)
774 //
775 //
776 // INPUT
777 // PD(3) q/P,V',W'
778 // H(3) MAGNETIC FIELD
779 // RD(15) ERROR MATRIX IN 1/P,V',W',V,W (Triangular)
780 // CH CHARGE OF PARTICLE
781 // CHARGE AND MAGNETIC FIELD ARE NEEDED
782 // FOR CORRELATION TERMS (LAMBDA,V),(LAMBDA,W),(PHI,V),(PHI,W)
783 // THESE CORRELATION TERMS APPEAR BECAUSE RC IS ASSUMED
784 // TO BE THE ERROR MATRIX FOR FIXED S (PATH LENGTH)
785 // AND RD FOR FIXED U
786 // SPU SIGN OF U-COMPONENT OF PARTICLE MOMENTUM
787 // spu = sign[p·(DJ x DK)]
788 // DJ(3) UNIT VECTOR IN V-DIRECTION
789 // DK(3) UNIT VECTOR IN W-DIRECTION OF DETECTOR SYSTEM
790 //
791 
792 // OUTPUT
793 // PC(3) q/P,LAMBDA,PHI
794 // RC(15) ERROR MATRIX IN SC VARIABLES
795 // IERR NOT USED
796 //
797 // ---------------------------------------------------------------------
798 
799 
800  Double_t A[5][5], S[15], TN[3];
801  Double_t SINZ, COSZ, COSL, COSL1;
802  Double_t UN[3], VN[3], DI[3], TVW[3];
803 
804  Double_t Vec[25];
805  Double_t CFACT8= 2.997925e-04;
806  Double_t HA, HM, HAM, Q, UJ, UK, VJ, VK;
807 
808  Double_t PM;
809 
810  IERR=0;
811  memset(RC,0,sizeof(*RC));
812  memset(PC,0,sizeof(*PC));
813 
814  PM=PD[0];
815  TVW[0]=1./TMath::Sqrt(1.+PD[1]*PD[1]+PD[2]*PD[2]);
816 
817  if(SPU < 0.) { TVW[0]=-TVW[0]; }
818  TVW[1]=PD[1]*TVW[0];
819  TVW[2]=PD[2]*TVW[0];
820 
821  DI[0]=DJ[1]*DK[2]-DJ[2]*DK[1];
822  DI[1]=DJ[2]*DK[0]-DJ[0]*DK[2];
823  DI[2]=DJ[0]*DK[1]-DJ[1]*DK[0];
824 
825  for(Int_t I=0; I<3; I++) {
826  TN[I]=TVW[0]*DI[I]+TVW[1]*DJ[I]+TVW[2]*DK[I];
827  }
828 
829  PC[0] = PD[0];
830  PC[1] = TMath::ASin(TN[2]);
831  if(TMath::Abs(TN[0]) < 1.e-30) { TN[0] = 1.e-30; }
832 
833  PC[2] = TMath::ATan2(TN[1],TN[0]);
834 
835  COSL=TMath::Sqrt(TMath::Abs(1.-TN[2]*TN[2]));
836  if(COSL < 1.e-30) { COSL = 1.e-30; }
837  COSL1 = 1./COSL;
838  UN[0] = -TN[1]*COSL1;
839  UN[1] = TN[0]*COSL1;
840  UN[2] = 0.;
841 
842  VN[0] = -TN[2]*UN[1];
843  VN[1] = TN[2]*UN[0];
844  VN[2] = COSL;
845 
846  UJ = UN[0]*DJ[0]+UN[1]*DJ[1]+UN[2]*DJ[2];
847  UK = UN[0]*DK[0]+UN[1]*DK[1]+UN[2]*DK[2];
848  VJ = VN[0]*DJ[0]+VN[1]*DJ[1]+VN[2]*DJ[2];
849  VK = VN[0]*DK[0]+VN[1]*DK[1]+VN[2]*DK[2];
850 
851  // prepare matrices and vectors
852 
853  for(Int_t I=0; I<5; I++) {
854  for(Int_t K=0; K<5; K++) {
855  A[I][K]=0.;
856  }
857  }
858  for(Int_t J=0; J<15; J++) { S[J]=RD[J]; }
859 
860  if(CH != 0.) {
861  // charged particle
862  HA=TMath::Sqrt(H[0]*H[0]+H[1]*H[1]+H[2]*H[2]);
863  HAM=HA*PM;
864  if(HAM != 0.) {
865  // .... in a magnetic field
866  HM=CH/HA;
867  Q=-HAM*CFACT8;
868  SINZ=-(H[0]*UN[0]+H[1]*UN[1]+H[2]*UN[2])*HM;
869  COSZ= (H[0]*VN[0]+H[1]*VN[1]+H[2]*VN[2])*HM;
870  A[1][3] = -Q*TVW[1]*SINZ;
871  A[1][4] = -Q*TVW[2]*SINZ;
872  A[2][3] = -Q*TVW[1]*COSZ*COSL1;
873  A[2][4] = -Q*TVW[2]*COSZ*COSL1;
874  }
875  }
876  A[0][0] = 1.;
877  A[1][1] = TVW[0]*VJ;
878  A[1][2] = TVW[0]*VK;
879  A[2][1] = TVW[0]*UJ*COSL1;
880  A[2][2] = TVW[0]*UK*COSL1;
881  A[3][3] = UJ;
882  A[3][4] = UK;
883  A[4][3] = VJ;
884  A[4][4] = VK;
885 
886  // transformation matrix
887  FromMatToVec(A, Vec);
888  SymmProd(Vec,S,S);
889 
890  // copy the result in S in the output vector
891  for(Int_t I=0; I<15; I++) { RC[I]=S[I]; }
892 
893 }
894 
895 void FairGeaneUtil::FromVec15ToMat25(Double_t V[15], fiveMat& A)
896 {
897  //
898  // ------------------------------------------------------
899  // Passage from a 15-dim vector to a symmetric 5x5 matrix
900  // following the upper triangular convention
901  //
902  // Author A. Rotondi June 2007
903  //
904  // ------------------------------------------------------
905 
906  A[0][0] = V[0];
907  A[0][1] = V[1];
908  A[0][2] = V[2];
909  A[0][3] = V[3];
910  A[0][4] = V[4];
911 
912  A[1][1] = V[5];
913  A[1][2] = V[6];
914  A[1][3] = V[7];
915  A[1][4] = V[8];
916 
917  A[2][2] = V[9];
918  A[2][3] = V[10];
919  A[2][4] = V[11];
920 
921  A[3][3] = V[12];
922  A[3][4] = V[13];
923 
924  A[4][4] = V[14];
925 
926  for(Int_t I=0; I<5; I++) {
927  for(Int_t k=I; k<5; k++) {
928  A[k][I] = A[I][k];
929  }
930  }
931 }
932 void FairGeaneUtil::FromVecToMat(fiveMat& A, Double_t V[25])
933 {
934  //
935  // ------------------------------------------------------
936  // Passage from 25-dim vector to a symmetric 5x5 matrix
937  // (FoRTRAN column convention)
938  //
939  // Author A. Rotondi June 2007
940  //
941  // ------------------------------------------------------
942 
943  A[0][0] = V[0];
944  A[1][0] = V[1];
945  A[2][0] = V[2];
946  A[3][0] = V[3];
947  A[4][0] = V[4];
948 
949  A[0][1] = V[5];;
950  A[1][1] = V[6];;
951  A[2][1] = V[7];;
952  A[3][1] = V[8];;
953  A[4][1] = V[9];;
954 
955  A[0][2] = V[10] ;
956  A[1][2] = V[11] ;
957  A[2][2] = V[12] ;
958  A[3][2] = V[13] ;
959  A[4][2] = V[14] ;
960 
961  A[0][3] = V[15] ;
962  A[1][3] = V[16] ;
963  A[2][3] = V[17] ;
964  A[3][3] = V[18] ;
965  A[4][3] = V[19] ;
966 
967  A[0][4] = V[20] ;
968  A[1][4] = V[21] ;
969  A[2][4] = V[22] ;
970  A[3][4] = V[23] ;
971  A[4][4] = V[24] ;
972 
973 }
974 
975 void FairGeaneUtil::FromMarsToSC(Double_t PD[3], Double_t RD[6][6], Double_t H[3], Int_t CH,
976  // output
977  Double_t* PC, Double_t* RC)
978 {
979 // ----------------------------------------------------------------------
980 //
981 // Tranform error matrix
982 // FROM MASTER VARIABLES (px, py,pz, x, y, z)
983 // TO SC (transverse system) VARIABLES (q/p, lambda, phi, yt, zt)
984 // momentum along x axis (GEANE convention)
985 //
986 // Method: the matrix in MARS is transformed in a detctor SD system
987 // with the detector plane coincident with the SC transverse plane.
988 // Then, the SD to SC routine is used.
989 // In this way the track length s is not modified by the transformation.
990 //
991 // Authors: A. Rotondi and A. Fontana (June 2007)
992 // rewritten by A. Rotondi and Lia Lavezzi (may 2008)
993 //
994 // INPUT
995 // PD(3) px, py, pz
996 // RD(6,6) ERROR MATRIX from MASTER Reference System
997 // covariances of (px, py, pz, x, y, z)
998 //
999 // CH CHARGE OF PARTICLE
1000 // H(3) MAGNETIC FIELD components
1001 //
1002 // OUTPUT
1003 // PC(3) q/p, lambda, phi
1004 // RC(15) ERROR MATRIX IN SC VARIABLES
1005 //
1006 //
1007 // ---------------------------------------------------------------------
1008 
1009 
1010 // Double_t PDD[3], RDD[15];
1011 
1012  Double_t SPU, DJ[3], DK[3], PM, PM3, PT;
1013  Int_t IERR;
1014  Double_t clam, slam, cphi, sphi, PC1[3], RC1[15];
1015  // ------------------------------------------------------------------
1016 
1017  // reset
1018 
1019  IERR=0;
1020  memset(RC,0,sizeof(*RC));
1021  memset(PC,0,sizeof(*PC));
1022 
1023  PM = TMath::Sqrt(PD[0]*PD[0]+PD[1]*PD[1]+PD[2]*PD[2]);
1024  PM3 = TMath::Power(PM,3);
1025  PT = TMath::Sqrt(PD[0]*PD[0]+PD[1]*PD[1]);
1026 
1027  // prepare the director cosines of a virtual dtector system
1028  // lying on the SC transverse plane
1029 
1030  slam = PD[2]/PM;
1031  clam = TMath::Sqrt(1.-slam*slam);
1032 
1033 
1034  if(TMath::Abs(clam)<1.e-10) {
1035  clam = 0.;
1036  slam = 1.;
1037  cphi = 0.;
1038  sphi =-1.;
1039  } else {
1040  cphi = PD[0]/(clam*PM);
1041  sphi = PD[1]/(clam*PM);
1042  }
1043 
1044  DJ[0] = -sphi;
1045  DJ[1] = cphi;
1046  DJ[2] = 0.;
1047 
1048  DK[0] = -slam*cphi;
1049  DK[1] = -slam*sphi;
1050  DK[2] = clam;
1051 
1052  FromMarsToSD(PD, RD, H, CH, DJ, DK,
1053  IERR, SPU, PC1, RC1) ;
1054 
1055 
1056  // from SD to SC
1057 
1058  FromSDToSC(PC1, RC1, H, CH, SPU, DJ, DK, IERR, PC, RC);
1059 }
1060 
1061 
1062 
1063 
1064 void FairGeaneUtil::FromSCToMars(Double_t PC[3], Double_t RC[15], Double_t H[3], Int_t CH,
1065  // output
1066  Double_t* PD, sixMat& RD)
1067 {
1068 // ------------------------------------------------------------------------
1069 //
1070 // Transform error matrix
1071 //
1072 // FROM SC (transverse system) (q/P,LAMBDA,PHI,YT,ZT)
1073 // TO MASTER Reference System (MARS) (px, py, pz, z, y, z)
1074 //
1075 // Method: the SC system is transformed in a SD system with the
1076 // detector plane coincident with the transverse one.
1077 // Then, the SD to MARS routine is used.
1078 //
1079 // Authors: A. Rotondi and A. Fontana (June 2007)
1080 // rewritten by A. Rotondi and Lia Lavezzi (may 2008)
1081 // In this way the track length s is not modified by the transformation.
1082 //
1083 // INPUT
1084 // PC(3) q/p lambda phi
1085 // RC(15) ERROR MATRIX IN SC VARIABLES
1086 // covariances of (px, py, pz, x, y, z)
1087 // H(3) Magnetic field components
1088 // CH CHARGE OF PARTICLE
1089 //
1090 // OUTPUT
1091 // PD(3) px, py, pz
1092 // RD(6,6) ERROR MATRIX in MASTER Reference System
1093 //
1094 //
1095 // ---------------------------------------------------------------------------
1096 
1097 // Double_t M65[6][5], M65T[5][6], AJ[5][6];
1098  //Double_t DJ[3], DK[3], RCM[5][5];
1099  Double_t DJ[3], DK[3];
1100  Double_t PDD[3], RDD[15];
1101 
1102  Int_t IERR;
1103 // Double_t SPU, PM, PM2, PM3, PVW, PVW3;
1104  Double_t SPU;
1105  Double_t clam, slam, cphi, sphi;
1106  // -------------------------------------------------------------------------
1107 
1108  IERR=0;
1109  memset(PD,0,sizeof(*PD));
1110 
1111  // reset matrices
1112  for(Int_t I=0; I<6; I++) {
1113  for(Int_t K=0; K<6; K++) {
1114  RD[I][K] = 0.;
1115  }
1116  }
1117 
1118  // go from SC to SD with the same plane
1119  // prepare the director cosines of a virtual dtector system
1120  // lying on the SC transverse plane
1121 
1122  clam = TMath::Cos(PC[1]);
1123  slam = TMath::Sin(PC[1]);
1124  cphi = TMath::Cos(PC[2]);
1125  sphi = TMath::Sin(PC[2]);
1126 
1127 
1128  // momentum is perpendicular to the x-y plane
1129  if(TMath::Abs(clam)<1.e-15) {
1130  clam = 0.;
1131  slam = 1.;
1132  cphi = 0.;
1133  sphi =-1.;
1134  }
1135 
1136  // GEANE routines to go on SD
1137 
1138  DJ[0] = -sphi;
1139  DJ[1] = cphi;
1140  DJ[2] = 0.;
1141 
1142  DK[0] = -slam*cphi;
1143  DK[1] = -slam*sphi;
1144  DK[2] = clam;
1145 
1146  FromSCToSD(PC, RC, H, CH, DJ, DK, IERR, SPU, PDD, RDD);
1147 
1148  if(IERR !=1) {
1149 
1150  FromSDToMars( PDD, RDD, H, CH, SPU, DJ, DK,
1151  PD, RD);
1152  }
1153 
1154 }
1155 
1156 void FairGeaneUtil::FromMarsToSD(Double_t PD[3], Double_t RD[6][6],
1157  Double_t H[3], Int_t CH,
1158  Double_t DJ1[3], Double_t DK1[3],
1159  // output
1160  Int_t& IERR, Double_t& SP1,
1161  Double_t* PC, Double_t* RC)
1162 {
1163 
1164 // ----------------------------------------------------------------------
1165 //
1166 // Tranform error matrix
1167 // FROM MASTER VARIABLES (px, py,pz, x, y, z)
1168 // TO SD (transverse or local system)
1169 // VARIABLES (q/p, v', w', v, w)
1170 //
1171 // Method: the MARS system is rotated to a local cartesia system
1172 // with the x-y plane on the v-w one of SD. Hence eq (79) of the
1173 // report CMS 2006/001 is used to go from canonical to SD variables.
1174 // In this way the track length variation and the magnetic field
1175 // effects are correctly taken into account.
1176 //
1177 // Authors: A. Rotondi and A. Fontana (July 2007)
1178 // Rewritten by A. Rotondi and Lia Lavezzi (may 2008)
1179 // corrected for the q/p variable by A. Rotondi (10 june 2007)
1180 //
1181 // INPUT
1182 // PD(3) px, py, pz in MARS
1183 // RD(6,6) ERROR MATRIX from MASTER Reference System
1184 // covariances of (px, py, pz, x, y, z)
1185 //
1186 // CH CHARGE OF PARTICLE
1187 // H(3) MAGNETIC FIELD components
1188 // DJ1(3) Director cosines of axis v in MARS
1189 // DK1(3) Director cosines of axis w in MARS
1190 //
1191 // OUTPUt
1192 // IERR = 0 TRANSFORMATION OK
1193 // = 1 MOMENTUM LIES IN THE DETECTOR PLANE
1194 //
1195 // PC(3) q/p, v', w'
1196 // RC(15) ERROR MATRIX IN SD VARIABLES
1197 // SP1 SIGN OF U-COMPONENT OF PARTICLE MOMENTUM
1198 // SP1 = sign[p.(DJ x DK)]
1199 //
1200 //
1201 // ---------------------------------------------------------------------
1202 
1203 
1204 // Double_t PDD[3], RDD[15];
1205  Double_t PDD[3];
1206  Double_t M56[5][6], M56T[6][5], AJ[5][5], AJT[6][5];
1207  Double_t R6[6][6], RLC[6][6];
1208 // Double_t SPU, PM, PM3, PT;
1209  Double_t PM, PM3, PT;
1210  Double_t Rot[3][3], Rmat[6][6], Rtra[6][6];
1211  // ------------------------------------------------------------------
1212 
1213  // reset
1214 
1215  IERR=0;
1216  memset(RC,0,sizeof(*RC));
1217  memset(PC,0,sizeof(*PC));
1218 
1219  for(Int_t I=0; I<5; I++) {
1220  for(Int_t K=0; K<6; K++) {
1221  AJT[K][I]=0.;
1222  M56[I][K]=0.;
1223  M56T[K][I]= 0.;
1224  if(K != 5) { AJ[I][K]=0.; }
1225  }
1226  }
1227 
1228  for(Int_t i=0; i<6; i++) {
1229  for(Int_t k=0; k<6; k++) {
1230  R6[i][k] = 0.;
1231  RLC[i][k] = 0;
1232  Rmat[i][k] = 0.;
1233  Rtra[i][k] = 0.;
1234  }
1235  }
1236 
1237  // to local cartesian frame
1238 
1239  TVector3 MI(1.,0.,0.);
1240  TVector3 MJ(0.,1.,0.);
1241  TVector3 MK(0.,0.,1.);
1242  //local cartesian
1243  TVector3 DI(DJ1[0],DJ1[1],DJ1[2]);
1244  TVector3 DJ(DK1[0],DK1[1],DK1[2]);
1245  TVector3 DK= DI.Cross(DJ);
1246 
1247  // rotation: D.. final versors, M.. initial
1248  // vectors are column matrices
1249 
1250  Rot[0][0] = DI.Dot(MI);
1251  Rot[0][1] = DI.Dot(MJ);
1252  Rot[0][2] = DI.Dot(MK);
1253  Rot[1][0] = DJ.Dot(MI);
1254  Rot[1][1] = DJ.Dot(MJ);
1255  Rot[1][2] = DJ.Dot(MK);
1256  Rot[2][0] = DK.Dot(MI);
1257  Rot[2][1] = DK.Dot(MJ);
1258  Rot[2][2] = DK.Dot(MK);
1259 
1260  TMatrixT<double> PD1(3,1);
1261  PD1[0][0] = PD[0];
1262  PD1[1][0] = PD[1];
1263  PD1[2][0] = PD[2];
1264 
1265  TMatrixT<double> Rot1(3,3);
1266  for(int i = 0; i < 3; i++) for(int j = 0; j < 3; j++) {
1267  Rot1[i][j] = Rot[i][j];
1268  }
1269 
1270  // momentum components in the lcl cartesian frame
1271  // with the x-y plane on the detectorr one
1272 
1273  TMatrixT<double> Rot1xPD1(Rot1,TMatrixT<double>::kMult,PD1);
1274  PDD[0] = Rot1xPD1[0][0];
1275  PDD[1] = Rot1xPD1[1][0];
1276  PDD[2] = Rot1xPD1[2][0];
1277 
1278  // jacobian for error matrix in MARS
1279 
1280  Rmat[0][0] = Rot[0][0];
1281  Rmat[0][1] = Rot[0][1];
1282  Rmat[0][2] = Rot[0][2];
1283  Rmat[1][0] = Rot[1][0];
1284  Rmat[1][1] = Rot[1][1];
1285  Rmat[1][2] = Rot[1][2];
1286  Rmat[2][0] = Rot[2][0];
1287  Rmat[2][1] = Rot[2][1];
1288  Rmat[2][2] = Rot[2][2];
1289 
1290  Rmat[3][3] = Rot[0][0];
1291  Rmat[3][4] = Rot[0][1];
1292  Rmat[3][5] = Rot[0][2];
1293  Rmat[4][3] = Rot[1][0];
1294  Rmat[4][4] = Rot[1][1];
1295  Rmat[4][5] = Rot[1][2];
1296  Rmat[5][3] = Rot[2][0];
1297  Rmat[5][4] = Rot[2][1];
1298  Rmat[5][5] = Rot[2][2];
1299 
1300 // transposed of the jacobian
1301 
1302  for(Int_t I=0; I<6; I++) {
1303  for(Int_t K=0; K<6; K++) {
1304  Rtra[K][I]=Rmat[I][K];
1305  }
1306  }
1307 
1308  // product (J)(RD)(J+)
1309 
1310  for(Int_t i=0; i<6; i++) {
1311  for(Int_t k=0; k<6; k++) {
1312  for(Int_t l=0; l<6; l++) {
1313  R6[i][k] += RD[i][l]*Rtra[l][k];
1314  }
1315  }
1316  }
1317 
1318  for(Int_t i=0; i<6; i++) {
1319  for(Int_t k=0; k<6; k++) {
1320  for(Int_t l=0; l<6; l++) {
1321  RLC[i][k] += Rmat[i][l]*R6[l][k];
1322  }
1323  }
1324  }
1325  // go from local cartesian to the detector system SD
1326  // x-y of local are v and w of SD.
1327  // use of eq. (79) of the CMS report 2006/001 (Strandlie and Wittek)
1328  //
1329 
1330  PM = TMath::Sqrt(PDD[0]*PDD[0]+PDD[1]*PDD[1]+PDD[2]*PDD[2]);
1331  PM3 = TMath::Power(PM,3);
1332  PT = TMath::Sqrt(PDD[0]*PDD[0]+PDD[1]*PDD[1]);
1333 
1334  //
1335  // check if track lies in the x-y plane of the local cartesian frame
1336  // if so, IERR=1 and exit
1337  //
1338 
1339  if(TMath::Abs(PDD[2]) < 1.e-08) {
1340  IERR = 1;
1341  }
1342 
1343  else {
1344 
1345  // output q/p, v' and w' in SD
1346 
1347  PC[0] = CH/PM;
1348  PC[1] = PDD[0]/PDD[2];
1349  PC[2] = PDD[1]/PDD[2];;
1350 
1351  // Jacobian Mars --> SD parallel to Mars x-y plane
1352  // eq (79) of CMS note 2006/001 (Strandle and Wittek)
1353 
1354  M56[0][0] = - CH*PDD[0]/PM3;
1355  M56[0][1] = - CH*PDD[1]/PM3;
1356  M56[0][2] = - CH*PDD[2]/PM3;
1357 
1358  M56[1][0] = 1./PDD[2];
1359  M56[1][1] = 0.;
1360  M56[1][2] = - PDD[0]/(PDD[2]*PDD[2]);
1361 
1362  M56[2][0] = 0.;
1363  M56[2][1] = 1./PDD[2];
1364  M56[2][2] = - PDD[1]/(PDD[2]*PDD[2]);
1365 
1366  M56[3][3] = 1.;
1367  M56[4][4] = 1.;
1368 
1369 
1370  // matrix multiplication with the Jacobian
1371 
1372  for(Int_t k=0; k<5; k++) {
1373  for(Int_t l=0; l<6; l++) {
1374  M56T[l][k]= M56[k][l];
1375  }
1376  }
1377 
1378  for(Int_t i=0; i<6; i++) {
1379  for(Int_t k=0; k<5; k++) {
1380  for(Int_t l=0; l<6; l++) {
1381  AJT[i][k] += RLC[i][l]*M56T[l][k];
1382  }
1383  }
1384  }
1385 
1386  for(Int_t i=0; i<5; i++) {
1387  for(Int_t k=0; k<5; k++) {
1388  for(Int_t l=0; l<6; l++) {
1389  AJ[i][k] += M56[i][l]*AJT[l][k];
1390  }
1391  }
1392  }
1393 
1394  // SD format output erro matrix RC(15)
1395  FromMat25ToVec15(AJ,RC);
1396 
1397  SP1 = TMath::Sign(1., PD[0]*(DJ1[1]*DK1[2]-DJ1[2]*DK1[1])+
1398  PD[1]*(DJ1[2]*DK1[0]-DJ1[0]*DK1[2])+
1399  PD[2]*(DJ1[0]*DK1[1]-DJ1[1]*DK1[0]) );
1400  }
1401 }
1402 
1403 
1404 
1405 void FairGeaneUtil::FromSDToMars(Double_t PC[3], Double_t RC[15],
1406  Double_t H[3], Int_t CH,
1407  Double_t SP1, Double_t DJ1[3], Double_t DK1[3],
1408  // output
1409  Double_t* PD, sixMat& RD)
1410 {
1411 // ------------------------------------------------------------------------
1412 //
1413 // Transform error matrix
1414 //
1415 // FROM SD (detector system system) (q/P,v'.w'.v.w)
1416 // TO MASTER Reference System (MARS) (px, py, pz, z, y, z)
1417 //
1418 // Method: eq (80) of the report CMS 2006/001 is used
1419 // to go from SD to the local cartesian frame.
1420 // Then the error matrix in MARS is obtained through the jacobian
1421 // between the local cartesian frame and MARS.
1422 //
1423 // Authors: A. Rotondi and A. Fontana (June 2007)
1424 // Rewritten by A. Rotondi and Lia Lavezzi (may 2008)
1425 // corrected for the q/p variable by A. Rotondi (10 june 2007)
1426 //
1427 // INPUT
1428 // PC(3) q/p v' w'
1429 // RC(15) ERROR MATRIX IN SD VARIABLES
1430 // covariances of (px, py, pz, x, y, z)
1431 // H(3) Magnetic field components
1432 // CH CHARGE OF PARTICLE
1433 // DJ1(3) Director cosines of axis v in MARS
1434 // DK1(3) Director cosines of axis w in MARS
1435 // SP1 SIGN OF U-COMPONENT OF PARTICLE MOMENTUM
1436 // spu = sign[p·(DJ1 x DK1)]
1437 //
1438 // OUTPUT
1439 // PD(3) px, py, pz
1440 // RD(6,6) ERROR (Covariance) MATRIX in MASTER Reference System
1441 //
1442 //
1443 // ---------------------------------------------------------------------------
1444 
1445  Double_t M65[6][5], M65T[5][6], AJ[5][6];
1446  Double_t RCM[5][5];
1447  Double_t PDD[3];
1448  Double_t Rot[3][3];
1449 
1450  // TVector3 PD1, PD2;
1451  TVector3 PD2;
1452 
1453  Double_t RD1[6][6], Rmat[6][6], AJJ[6][6], Rtra[6][6];
1454 
1455  Double_t SPU, PM, PM2, PVW, PVW3;
1456 
1457  // -------------------------------------------------------------------------
1458 
1459  memset(PD,0,sizeof(*PD));
1460 
1461  // reset matrices
1462  for(Int_t I=0; I<5; I++) {
1463  for(Int_t K=0; K<6; K++) {
1464  AJ[I][K]=0.;
1465  M65[K][I]=0.;
1466  M65T[I][K]=0.;
1467  RD[I][K] = 0.;
1468  RD1[I][K] = 0.;
1469  Rmat[I][K] = 0.;
1470  Rtra[I][K] = 0.;
1471  AJJ[I][K] = 0.;
1472  }
1473  }
1474  for(Int_t I=0; I<6; I++) {
1475  RD[5][I] = 0.;
1476  RD1[5][I] = 0.;
1477  Rmat[5][I] = 0.;
1478  Rtra[5][I] = 0.;
1479  AJJ[5][I] = 0.;
1480  }
1481 
1482  SPU = SP1;
1483 
1484  // jacobian from SD to local cartesian
1485 
1486  PM = 1.e+30;
1487  if(PC[0] != 0.) { PM =CH/PC[0]; }
1488  PM2 = PM*PM;
1489  PVW = TMath::Sqrt(1.+PC[1]*PC[1]+PC[2]*PC[2]);
1490  PVW3 = TMath::Power(PVW,3);
1491 
1492  // output px, py, pz (cartesian)
1493 
1494  PDD[0] = SPU*PM*PC[1]/PVW;
1495  PDD[1] = SPU*PM*PC[2]/PVW;
1496  PDD[2] = SPU*PM/PVW ;
1497 
1498  // Jacobian SD --> Mars
1499  // eq (80) of CMS note 2006/001 (Strandlie and Wittek)
1500 
1501  M65[0][0] = - SPU*PM2*PC[1]/(CH*PVW);
1502  M65[1][0] = - SPU*PM2*PC[2]/(CH*PVW);
1503  M65[2][0] = - SPU*PM2/(CH*PVW);
1504 
1505  M65[0][1] = SPU*PM*(1.+PC[2]*PC[2])/PVW3;
1506  M65[1][1] = - SPU*PM*PC[1]*PC[2]/PVW3;
1507  M65[2][1] = - SPU*PM*PC[1]/PVW3;
1508 
1509  M65[0][2] = - SPU*PM*PC[1]*PC[2]/PVW3;
1510  M65[1][2] = SPU*PM*(1.+PC[1]*PC[1])/PVW3;
1511  M65[2][2] = - SPU*PM*PC[2]/PVW3;
1512 
1513  M65[3][3] = 1.;
1514  M65[4][4] = 1.;
1515 
1516  FromVec15ToMat25(RC,RCM);
1517 
1518 
1519  // transposed of the jacobian
1520 
1521  for(Int_t I=0; I<6; I++) {
1522  for(Int_t K=0; K<5; K++) {
1523  M65T[K][I]=M65[I][K];
1524  }
1525  }
1526 
1527  // product (J)(RCM)(J+)
1528 
1529  for(Int_t i=0; i<5; i++) {
1530  for(Int_t k=0; k<6; k++) {
1531  for(Int_t l=0; l<5; l++) {
1532  AJ[i][k] += RCM[i][l]*M65T[l][k];
1533  }
1534  }
1535  }
1536 
1537  for(Int_t i=0; i<6; i++) {
1538  for(Int_t k=0; k<6; k++) {
1539  for(Int_t l=0; l<5; l++) {
1540  RD1[i][k] += M65[i][l]*AJ[l][k];
1541  }
1542  }
1543  }
1544 
1545 
1546  // now go from the local cartesian to MARS
1547 
1548  // MARS
1549  TVector3 MI(1.,0.,0.);
1550  TVector3 MJ(0.,1.,0.);
1551  TVector3 MK(0.,0.,1.);
1552  //local cartesian
1553  TVector3 DI(DJ1[0],DJ1[1],DJ1[2]);
1554  TVector3 DJ(DK1[0],DK1[1],DK1[2]);
1555  TVector3 DK= DI.Cross(DJ);
1556 
1557  // rotation: M.. final versors, D.. initial
1558  // vectors are column matrices
1559 
1560  Rot[0][0] = MI.Dot(DI);
1561  Rot[0][1] = MI.Dot(DJ);
1562  Rot[0][2] = MI.Dot(DK);
1563  Rot[1][0] = MJ.Dot(DI);
1564  Rot[1][1] = MJ.Dot(DJ);
1565  Rot[1][2] = MJ.Dot(DK);
1566  Rot[2][0] = MK.Dot(DI);
1567  Rot[2][1] = MK.Dot(DJ);
1568  Rot[2][2] = MK.Dot(DK);
1569 
1570  TMatrixT<double> PD1(3,1);
1571  PD1[0][0] = PDD[0];
1572  PD1[1][0] = PDD[1];
1573  PD1[2][0] = PDD[2];
1574 
1575  TMatrixT<double> Rot1(3,3);
1576  for(int i = 0; i < 3; i++) for(int j = 0; j < 3; j++) {
1577  Rot1[i][j] = Rot[i][j];
1578  }
1579 
1580  TMatrixT<double> Rot1xPD1(Rot1,TMatrixT<double>::kMult,PD1);
1581  PD[0] = Rot1xPD1[0][0];
1582  PD[1] = Rot1xPD1[1][0];
1583  PD[2] = Rot1xPD1[2][0];
1584 
1585  // jacobian for error matrix in MARS
1586 
1587  Rmat[0][0] = Rot[0][0];
1588  Rmat[0][1] = Rot[0][1];
1589  Rmat[0][2] = Rot[0][2];
1590  Rmat[1][0] = Rot[1][0];
1591  Rmat[1][1] = Rot[1][1];
1592  Rmat[1][2] = Rot[1][2];
1593  Rmat[2][0] = Rot[2][0];
1594  Rmat[2][1] = Rot[2][1];
1595  Rmat[2][2] = Rot[2][2];
1596 
1597  Rmat[3][3] = Rot[0][0];
1598  Rmat[3][4] = Rot[0][1];
1599  Rmat[3][5] = Rot[0][2];
1600  Rmat[4][3] = Rot[1][0];
1601  Rmat[4][4] = Rot[1][1];
1602  Rmat[4][5] = Rot[1][2];
1603  Rmat[5][3] = Rot[2][0];
1604  Rmat[5][4] = Rot[2][1];
1605  Rmat[5][5] = Rot[2][2];
1606 
1607  // transposed of the jacobian
1608 
1609  for(Int_t I=0; I<6; I++) {
1610  for(Int_t K=0; K<6; K++) {
1611  Rtra[K][I]=Rmat[I][K];
1612  }
1613  }
1614 
1615  for(Int_t i=0; i<6; i++) {
1616  for(Int_t k=0; k<6; k++) {
1617  for(Int_t l=0; l<6; l++) {
1618  AJJ[i][k] += RD1[i][l]*Rtra[l][k];
1619  }
1620  }
1621  }
1622 
1623  for(Int_t i=0; i<6; i++) {
1624  for(Int_t k=0; k<6; k++) {
1625  for(Int_t l=0; l<6; l++) {
1626  RD[i][k] += Rmat[i][l]*AJJ[l][k];
1627  }
1628  }
1629  }
1630 
1631 }
1632 
1633 
1634 
1635 void FairGeaneUtil::FromMat25ToVec15(Double_t A[5][5], Double_t* V)
1636 {
1637  //
1638  // ------------------------------------------------------
1639  // Passage from a symmetric 5x5 matrix to a 15-dim vector
1640  // following the upper triangular convention
1641  //
1642  // Author A. Rotondi June 2007
1643  //
1644  // ------------------------------------------------------
1645 
1646  V[0] = A[0][0];
1647  V[1] = A[0][1];
1648  V[2] = A[0][2];
1649  V[3] = A[0][3];
1650  V[4] = A[0][4];
1651 
1652  V[5] = A[1][1];
1653  V[6] = A[1][2];
1654  V[7] = A[1][3];
1655  V[8] = A[1][4];
1656 
1657  V[ 9] = A[2][2];
1658  V[10] = A[2][3];
1659  V[11] = A[2][4];
1660 
1661  V[12] = A[3][3];
1662  V[13] = A[3][4];
1663 
1664  V[14] = A[4][4];
1665 
1666 
1667 }
1668 
1669 void FairGeaneUtil::FromMatToVec(Double_t A[5][5], Double_t* V)
1670 {
1671  //
1672  // ------------------------------------------------------
1673  // Passage from a 5x5 matrix to a 25-dim vector
1674  //
1675  // Author A. Rotondi June 2007
1676  //
1677  // ------------------------------------------------------
1678 
1679  V[0] = A[0][0];
1680  V[1] = A[1][0];
1681  V[2] = A[2][0];
1682  V[3] = A[3][0];
1683  V[4] = A[4][0];
1684 
1685  V[5] = A[0][1];
1686  V[6] = A[1][1];
1687  V[7] = A[2][1];
1688  V[8] = A[3][1];
1689  V[9] = A[4][1];
1690 
1691  V[10] = A[0][2];
1692  V[11] = A[1][2];
1693  V[12] = A[2][2];
1694  V[13] = A[3][2];
1695  V[14] = A[4][2];
1696 
1697 
1698  V[15] = A[0][3];
1699  V[16] = A[1][3];
1700  V[17] = A[2][3];
1701  V[18] = A[3][3];
1702  V[19] = A[4][3];
1703 
1704  V[20] = A[0][4];
1705  V[21] = A[1][4];
1706  V[22] = A[2][4];
1707  V[23] = A[3][4];
1708  V[24] = A[4][4];
1709 
1710 }
1711 
1712 void FairGeaneUtil::SymmProd(Double_t A[25], Double_t S[15], Double_t* R)
1713 {
1714  //
1715  // ---------------------------------------------------------
1716  // TRANSFORMATION OF SYMMETRIC 5X5 MATRIX S:
1717  // A*S*AT -> R.
1718  // A is a 25-dim vector corresonding to the 5x5 transformation
1719  // matrix
1720  // S and R ARE SYMMETRIC ERROR (COVARIANCE) MATRICES STORED
1721  // IN TRIANGULAR FORM.
1722  //
1723  // INPUT
1724  // A[25] transformation matrix 5x5
1725  // S[15] error matrix in triangular form
1726  //
1727  // OUTPUT
1728  // R[15] error matrix in triangular form
1729  //
1730  // NB: S AND R MAY WELL BE THE SAME MATRIX.
1731  //
1732  // Author: A. Haas (Freiburg University) 5/7/81
1733  // transported with modifications
1734  // in C/Root by A. Rotondi (June 2007)
1735  //
1736  // * ------------------------------------------------------
1737 
1738  Double_t Q[15],T1,T2,T3,T4,T5;
1739 
1740  for(Int_t i=0; i<15; i++) { Q[i]=S[i]; }
1741 
1742  Int_t K =0;
1743  for(Int_t J=0; J<5; J++) {
1744  T1=A[J ];
1745  T2=A[J+ 5];
1746  T3=A[J+10];
1747  T4=A[J+15];
1748  T5=A[J+20];
1749  for(Int_t I=J; I<5; I++) {
1750  R[K]=A[I ]*(Q[0]*T1+Q[1]*T2+Q[ 2]*T3+Q[ 3]*T4+Q[ 4]*T5)
1751  +A[I+ 5]*(Q[1]*T1+Q[5]*T2+Q[ 6]*T3+Q[ 7]*T4+Q[ 8]*T5)
1752  +A[I+10]*(Q[2]*T1+Q[6]*T2+Q[ 9]*T3+Q[10]*T4+Q[11]*T5)
1753  +A[I+15]*(Q[3]*T1+Q[7]*T2+Q[10]*T3+Q[12]*T4+Q[13]*T5)
1754  +A[I+20]*(Q[4]*T1+Q[8]*T2+Q[11]*T3+Q[13]*T4+Q[14]*T5);
1755  K++;
1756  }
1757  }
1758 }
1759 
1760 // ------------------------- modifiche 27 jul 2007 --------------------------
1761 
1762 TVector3 FairGeaneUtil::FromMARSToSDCoord(TVector3 xyz, TVector3 o, TVector3 di, TVector3 dj, TVector3 dk)
1763 {
1764 
1765  TMatrixT<double> matrix(3,3);
1766  // rotation matrix (u,v,w) = (fmatrix) * (x,y,z)
1767  matrix[0][0] = di[0];
1768  matrix[0][1] = di[1];
1769  matrix[0][2] = di[2];
1770  matrix[1][0] = dj[0];
1771  matrix[1][1] = dj[1];
1772  matrix[1][2] = dj[2];
1773  matrix[2][0] = dk[0];
1774  matrix[2][1] = dk[1];
1775  matrix[2][2] = dk[2];
1776 
1777 
1778 
1779  TMatrixT<double> xyzvec(3,1);
1780  xyzvec[0][0] = xyz.X();
1781  xyzvec[1][0] = xyz.Y();
1782  xyzvec[2][0] = xyz.Z();
1783 
1784 
1785  TMatrixT<double> origin(3,1);
1786  origin[0][0] = o.X();
1787  origin[1][0] = o.Y();
1788  origin[2][0] = o.Z();
1789 
1790  xyzvec -= origin;
1791 
1792  TMatrixT<double> uvwvec(matrix, TMatrixT<double>::kMult, xyzvec);
1793 
1794  TVector3 uvw = TVector3(uvwvec[0][0], uvwvec[1][0], uvwvec[2][0]);
1795  return uvw;
1796 }
1797 
1798 TVector3 FairGeaneUtil::FromSDToMARSCoord(TVector3 uvw, TVector3 o, TVector3 di, TVector3 dj, TVector3 dk)
1799 {
1800  TMatrixT<double> matrix(3,3);
1801  matrix[0][0] = di[0];
1802  matrix[0][1] = di[1];
1803  matrix[0][2] = di[2];
1804  matrix[1][0] = dj[0];
1805  matrix[1][1] = dj[1];
1806  matrix[1][2] = dj[2];
1807  matrix[2][0] = dk[0];
1808  matrix[2][1] = dk[1];
1809  matrix[2][2] = dk[2];
1810 
1811 
1812  TMatrixT<double> uvwvec(3,1);
1813  uvwvec[0][0] = uvw.X();
1814  uvwvec[1][0] = uvw.Y();
1815  uvwvec[2][0] = uvw.Z();
1816 
1817  TMatrixT<double> uvwrot(matrix, TMatrixT<double>::kTransposeMult, uvwvec);
1818 
1819  TMatrixT<double> origin(3,1);
1820  origin[0][0] = o.X();
1821  origin[1][0] = o.Y();
1822  origin[2][0] = o.Z();
1823 
1824  TMatrixT<double> xyzvec(3,1);
1825  xyzvec = uvwrot + origin;
1826 
1827 
1828  TVector3 xyz = TVector3(xyzvec[0][0], xyzvec[1][0], xyzvec[2][0]);
1829 
1830  return xyz;
1831 }
1832 
1833 
1834 ClassImp(FairGeaneUtil)
1835