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