EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
FairGeanePro.cxx
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file FairGeanePro.cxx
1 // Class for the interface to propagate track parameters with GEANE
2 //
3 // Authors: M. Al-Turany, A. Fontana, L. Lavezzi and A. Rotondi
4 //
5 
6 #include "FairGeanePro.h"
7 #include "FairTrackParP.h"
8 #include "FairTrackParH.h"
9 #include "FairRunAna.h"
10 #include "FairField.h"
11 #include "FairGeaneUtil.h"
12 
13 #include "TGeant3TGeo.h"
14 #include "TVector3.h"
15 #include "TArrayD.h"
16 #include "TDatabasePDG.h"
17 #include "TGeoTorus.h"
18 #include "TMatrixFSym.h"
19 #include "TVirtualMC.h"
20 #include "TGeant3.h"
21 #include "FairGeaneApplication.h"
22 #include <iostream>
23 #include <cmath>
24 
25 using std::cout;
26 using std::endl;
27 
28 // ----- Default constructor -------------------------------------------
29 FairGeanePro::FairGeanePro()
30  : TNamed("Geane", "Propagate Tracks"),
31  //gMC3((TGeant3*) gMC),
32  fPropOption(""),
33  nepred(1),
34  fdbPDG(TDatabasePDG::Instance()),
35  //afErtrio(gMC3->fErtrio),
36  GeantCode(0),
37  ProMode(0),
38  VName(""),
39  VCopyNo(0),
40  VEnter(kTRUE),
41  fpoint(TVector3(0., 0., 0.)),
42  fwire1(TVector3(0., 0., 0.)),
43  fwire2(TVector3(0., 0., 0.)),
44  fPCA(0),
45  fRad(0.),
46  fDi(0.),
47  fvpf(TVector3(0., 0., 0.)),
48  fvwi(TVector3(0., 0., 0.)),
49  ftrklength(0.),
50  flag(0),
51  fApp(FairGeaneApplication::Instance())
52 {
53  gMC3 = (TGeant3*)gMC;//),
54  if(gMC3==NULL) {
55  std::cerr<<"FairGeanePro::TGeant3 has not been initialized! ABORTING!"<<std::endl;
56  throw;
57  }
58  afErtrio = gMC3->fErtrio;
59 
60  // nepred=1;
61  xlf[0] = 0.;
62  // fdbPDG= TDatabasePDG::Instance();
63  // fErrorMat= new TArrayD(15);
64  // afErtrio=gMC3->fErtrio;
65  // Pos=TVector3(0, 0 , 0);
66  // PosErr = TVector3(0,0,0);
67  // Mom=TVector3(0,0,0);
68  // fTrkPar= new FairTrackPar();
69  // ProMode=0;
70  // FairRunAna *fRun= FairRunAna::Instance();
71  // fField= fRun->GetField();
72 
73  // PCA stuff
74  // fPCA = 0;
75  /*
76  fpoint = TVector3(0., 0., 0.);
77  fwire1 = TVector3(0., 0., 0.);
78  fwire2 = TVector3(0., 0., 0.);
79  */
80 
81  for(int i = 0; i < 5; i++) for(int j = 0; j < 5; j++) { trpmat[i][j] = 0.; }
82 
83  // fApp = FairGeaneApplication::Instance();
84 
85  // fPropOption = "";
86  for(int i = 0; i < 15; i++) {
87  ein[i] = 0.;
88  if(i < 3) {
89  x2[i] = 0;
90  p2[i] = 0;
91  x1[i] = 0;
92  p1[i] = 0;
93  }
94  }
95 
96  pli[0] = 1.;
97  pli[1] = 0.;
98  pli[2] = 0.;
99  pli[3] = 0.;
100  pli[4] = 1.;
101  pli[5] = 0.;
102 
103  plo[0] = 0.;
104  plo[1] = 0.;
105  plo[2] = 0.;
106  plo[3] = 1.;
107  plo[4] = 0.;
108  plo[5] = 0.;
109  plo[6] = 0.;
110  plo[7] = 1.;
111  plo[8] = 0.;
112  plo[9] = 0.;
113  plo[10] = 0.;
114  plo[11] = 1.;
115 
116  // GeantCode = 0;
117  /*
118  VName = "";
119  VCopyNo = 0;
120  VEnter = kTRUE;
121  */
122 
123 }
124 
125 // ----- Destructor ----------------------------------------------------
126 FairGeanePro::~FairGeanePro() { }
127 
128 
129 
130 Bool_t FairGeanePro::Propagate(FairTrackParH* TParam, FairTrackParH* TEnd, Int_t PDG)
131 {
132  // Propagate a helix track and return a helix (SC system)
133 
134  //Bool_t NeedSDSC=kFALSE;
135  Int_t ch=1; // CHARGE OF PARTICLE
136 
137  Double_t fCov[15], fCovOut[15];
138  TParam->GetCovQ(fCov);
139 
140  Init(TParam);
141  Double_t Q = TParam->GetQ();
142  if (fabs(Q)>1.E-8) { ch= int (Q/TMath::Abs(Q)); }
143  if (ProMode==1) { //Propagate to Volume
144  //***** We have the right representation go further
145  for(Int_t i=0; i<15; i++) {
146  ein[i]=fCov[i];
147 
148  }
149  if(fPropOption.Contains("V")) {
150  //cout << "Propagate Helix to Volume" << endl;
151  Int_t option;
152  if(VEnter) { option =1; }
153  else { option =2; }
154  gMC3->Eufilv(1, ein, (Char_t*)VName.Data(), &VCopyNo, &option);
155  } else if(fPropOption.Contains("L")) {
156  if(fPCA == 0) {
157  gMC3->Eufill(nepred, ein,xlf);
158  } else if(fPCA != 0) {
159 
160  // max length estimate:
161  // we calculate the geometrical distance of the start point
162  // from the point/wire extremity and multiply it * 2
163  TVector3 start = TVector3(TParam->GetX(), TParam->GetY(), TParam->GetZ());
164  Double_t maxdistance = 0;
165  if(fPCA == 1) { maxdistance = (fpoint - start).Mag(); }
166  else if(fPCA == 2) {
167  Double_t distance1, distance2;
168  distance1 = (fwire1 - start).Mag();
169  distance2 = (fwire2 - start).Mag();
170  if(distance1 < distance2) { maxdistance = distance2; }
171  else { maxdistance = distance1; }
172  }
173  maxdistance *= 2.;
174 
175  // output
176  Int_t findpca = FindPCA(fPCA, PDG, fpoint, fwire1, fwire2, maxdistance, fRad, fvpf, fvwi, fDi, ftrklength);
177  if(findpca != 0) { return kFALSE; }
178 
179  // reset parameters
180  Init(TParam);
181  gMC3->Eufill(nepred, ein, &ftrklength);
182  }
183  }
184  } else if(ProMode ==3) {
185  cout << "Propagate Helix parameter to Plane is not implimented yet" << endl;
186  return kFALSE;
187  }
188  //Propagate
189  if(Propagate(PDG)==kFALSE) { return kFALSE; }
190 
191  for(Int_t i=0; i<15; i++) {
192  fCovOut[i]=afErtrio->errout[i];
193  if(i == 0) { fCovOut[i] = fCovOut[i] * ch * ch; }
194  if(i > 0 && i < 5) { fCovOut[i] = fCovOut[i] * ch; }
195  }
196 
197  // do not remove (useful for debug)
198  if(fabs(p2[0]) < 1e-9 && fabs(p2[1]) < 1e-9 && fabs(p2[2]) < 1e-9) { return kFALSE; }
199  TEnd->SetTrackPar(x2[0], x2[1], x2[2],p2[0],p2[1],p2[2], ch ,fCovOut );
200  return kTRUE;
201 
202 }
203 Bool_t FairGeanePro::Propagate(FairTrackParP* TStart, FairTrackParH* TEnd, Int_t PDG)
204 {
205  // Propagate a parabola track (SD system) and return a helix (SC system) (not used nor implemented)
206  cout << "FairGeanePro::Propagate(FairTrackParP *TParam, FairTrackParH &TEnd, Int_t PDG) : (not used nor implemented)" << endl;
207  return kFALSE;
208 }
209 
210 Bool_t FairGeanePro::Propagate(FairTrackParP* TStart, FairTrackParP* TEnd, Int_t PDG)
211 {
212  // Propagate a parabola track (SD system) and return a parabola (SD system) (not used nor implemented)
213  Int_t ch=1; // CHARGE OF PARTICLE
214  Double_t fCov[15], fCovOut[15];
215  TStart->GetCovQ(fCov);
216 
217  //printf("FairGeanePro::Propagate() #1\n");
218 
219  Init(TStart);
220  Double_t Q = TStart->GetQ() ;
221  if (Q!=0) { ch= int (Q/TMath::Abs(Q)); }
222 
223  if (ProMode==1) { //Propagate to Volume
224  cout << "Propagate Parabola parameter to Volume is not implimented yet" << endl;
225  return kFALSE;
226  } else if(ProMode ==3) {
228  for(Int_t i=0; i<15; i++) {
229  ein[i]=fCov[i];
230 
231  }
232 
233  if(fPropOption.Contains("P")) { gMC3->Eufilp(nepred, ein, pli, plo); }
234  else if(fPropOption.Contains("L")) {
235  if(fPCA == 0) { cout << "Propagate Parabola to Parabola in Length not yet implemented" << endl; }
236  else if(fPCA != 0) {
237 
238  // max length estimate:
239  // we calculate the geometrical distance of the start point
240  // from the point/wire extremity and multiply it * 2
241  TVector3 start = TVector3(TStart->GetX(), TStart->GetY(), TStart->GetZ());
242  Double_t maxdistance = 0;
243  if(fPCA == 1) { maxdistance = (fpoint - start).Mag(); }
244  else if(fPCA == 2) {
245  Double_t distance1, distance2;
246  distance1 = (fwire1 - start).Mag();
247  distance2 = (fwire2 - start).Mag();
248  if(distance1 < distance2) { maxdistance = distance2; }
249  else { maxdistance = distance1; }
250  }
251  maxdistance *= 2.;
252 
253  // output
254  Int_t findpca = FindPCA(fPCA, PDG, fpoint, fwire1, fwire2, maxdistance, fRad, fvpf, fvwi, fDi, ftrklength);
255  //std::cout<<" FairGeanePro::ftrklength="<<ftrklength<< std::endl;
256  if(findpca != 0) { return kFALSE; }
257 
258  // reset parameters
259  Init(TStart);
260 
261  // find plane
262  // unitary vector along distance
263  // fvpf on track, fvwi on wire
264  TVector3 fromwiretoextr = fvpf - fvwi;
265  fromwiretoextr.SetMag(1.);
266  if(fabs(fromwiretoextr.Mag()-1)>1E-4) {
267  std::cerr<<"fromwire.Mag()!=1"<<std::endl;
268  return kFALSE;
269  }
270 
271  //for wires:
272  // unitary vector along the wire
273  TVector3 wiredirection = fwire2 - fwire1;
274  if(fPCA==1) { // point
275  TVector3 mom(TStart->GetPx(),TStart->GetPy(),TStart->GetPz());
276  wiredirection=mom.Cross(fromwiretoextr);
277  }
278  wiredirection.SetMag(1.);
279  // check orthogonality
280  if(fabs(fromwiretoextr * wiredirection) > 1e-3) {
281  return kFALSE; // throw away the event
282  // wiredirection = (fromwiretoextr.Cross(wiredirection)).Cross(fromwiretoextr);
283  // wiredirection.SetMag(1.);
284  }
285 
286  TVector3 jver = TStart->GetJVer();;
287  TVector3 kver = TStart->GetKVer();
288  Bool_t backtracking = kFALSE;
289  if(fPropOption.Contains("B")) { backtracking = kTRUE; }
290  PropagateFromPlane(jver, kver);
291  PropagateToPlane(fvwi, fromwiretoextr, wiredirection);
292  if(backtracking == kTRUE) { fPropOption = "BPE"; }
293 
294  gMC3->Eufilp(nepred, ein, pli, plo);
295  }
296  }
297  }
298  //printf("FairGeanePro::Propagate() #2\n");
299  //Propagate
300  if(Propagate(PDG)==kFALSE) { return kFALSE; }
301  //printf("FairGeanePro::Propagate() #3\n");
302 
303  for(Int_t i=0; i<15; i++) {
304  fCovOut[i]=afErtrio->errout[i];
305  if(i == 0) { fCovOut[i] = fCovOut[i] * ch * ch; }
306  if(i > 0 && i < 5) { fCovOut[i] = fCovOut[i] * ch; }
307  }
308 
309  // plane
310  TVector3 origin(plo[6], plo[7], plo[8]);
311  TVector3 dj(plo[0], plo[1], plo[2]);
312  TVector3 dk(plo[3], plo[4], plo[5]);
313  TVector3 di(plo[9], plo[10], plo[11]); // = dj.Cross(dk);
314 
315 
316  if(fabs(p2[0]) < 1e-9 && fabs(p2[1]) < 1e-9 && fabs(p2[2]) < 1e-9) { return kFALSE; }
317 
318  TEnd->SetTrackPar(x2[0], x2[1], x2[2],p2[0],p2[1],p2[2], ch ,fCovOut, origin, di, dj, dk);
319 
320  return kTRUE;
321 
322 }
323 
324 
325 Bool_t FairGeanePro::Propagate(FairTrackParH* TStart, FairTrackParP* TEnd, Int_t PDG)
326 {
327  // Propagate a helix track (SC system) and return a parabola (SD system) (not used nor implemented)
328  cout << "FairGeanePro::Propagate(FairTrackParH *TParam, FairTrackParP &TEnd, Int_t PDG) not implemented" << endl;
329  return kFALSE;
330 }
331 
332 Bool_t FairGeanePro::Propagate(Float_t* X1, Float_t* P1, Float_t* X2, Float_t* P2,Int_t PDG)
333 {
334 // fApp->GeanePreTrack(X1, P1, PDG);
335  GeantCode=fdbPDG->ConvertPdgToGeant3(PDG);
336  xlf[0]=1000
337  ;
338  gMC3->Eufill(1, ein,xlf);
339  gMC3->Ertrak(X1,P1,X2,P2,GeantCode, "L");
340  if(X2[0]<-1.E29) { return kFALSE; }
341  if(gMC3->IsTrackOut()) { return kFALSE; }
342  return kTRUE;
343 }
344 
345 Bool_t FairGeanePro::Propagate(Int_t PDG)
346 {
347  // main propagate call to fortran ERTRAK
348 
349  //printf("FairGeanePro::Propagate(PDG) #1\n");
350  GeantCode=fdbPDG->ConvertPdgToGeant3(PDG);
351  //cout << " FairGeanePro::Propagate ---------------------------"<< " " << x1[0]<< " "<< x1[1]<< " "<< x1[2] << endl;
352  //fApp->GeanePreTrack(x1, p1, PDG);
353  gMC3->Ertrak(x1,p1,x2,p2,GeantCode, fPropOption.Data());
354  if(x2[0]<-1.E29) { return kFALSE; }
355  //printf("FairGeanePro::Propagate(PDG) #2\n");
356  if(gMC3->IsTrackOut()) { return kFALSE; }
357 
358  //printf("FairGeanePro::Propagate(PDG) #3\n");
359  ftrklength=gMC3->TrackLength();
360 
361  Double_t trasp[25];
362  for(int i = 0; i < 25; i++) {
363  // trasp[i] = afErtrio->ertrsp[i]; // single precision tr. mat.
364  trasp[i] = afErtrio->erdtrp[i]; // double precision tr. mat.
365  }
366  FairGeaneUtil fUtil;
367  fUtil.FromVecToMat(trpmat, trasp);
368  return kTRUE;
369 }
370 
371 void FairGeanePro::Init(FairTrackPar* TParam)
372 {
373  // starting and ending point initialization
374  x1[0]=TParam->GetX();
375  x1[1]=TParam->GetY();
376  x1[2]=TParam->GetZ();
377  p1[0]=TParam->GetPx();
378  p1[1]=TParam->GetPy();
379  p1[2]=TParam->GetPz();
380 
381  x2[0]=0;
382  x2[1]=0;
383  x2[2]=0;
384  p2[0]=0;
385  p2[1]=0;
386  p2[2]=0;
387 }
388 
389 Bool_t FairGeanePro::PropagateFromPlane(TVector3& v1, TVector3& v2)
390 {
391  // define initial plane (option "P")
392  TVector3 v1u=v1.Unit();
393  TVector3 v2u=v2.Unit();
394  pli[0]=v1u.X();
395  pli[1]=v1u.Y();
396  pli[2]=v1u.Z();
397  pli[3]=v2u.X();
398  pli[4]=v2u.Y();
399  pli[5]=v2u.Z();
400  return kTRUE;
401 }
402 
403 Bool_t FairGeanePro::PropagateToPlane(TVector3& v0, TVector3& v1, TVector3& v2)
404 {
405  // define final plane (option "P")
406  // uncomment to set the initial error to zero
407  // for(Int_t i=0;i<15;i++) ein[i]=0.00;
408  TVector3 v1u=v1.Unit();
409  TVector3 v2u=v2.Unit();
410 
411 
412  plo[0]=v1u.X();
413  plo[1]=v1u.Y();
414  plo[2]=v1u.Z();
415  plo[3]=v2u.X();
416  plo[4]=v2u.Y();
417  plo[5]=v2u.Z();
418 
419  plo[6]=v0.X();
420  plo[7]=v0.Y();
421  plo[8]=v0.Z();
422 
423  TVector3 v3=v1u.Cross(v2u);
424 
425  plo[9]=v3(0);
426  plo[10]=v3(1);
427  plo[11]=v3(2);
428 
429  fPropOption="PE";
430  ProMode=3; //need errors in representation 3 (SD)(see Geane doc)
431  // gMC3->Eufilp(nepred, ein, pli, plo);
432  return kTRUE;
433 }
434 Bool_t FairGeanePro::PropagateToVolume(TString VolName, Int_t CopyNo , Int_t option)
435 {
436  // define final volume (option "V")
437  for(Int_t i=0; i<15; i++) { ein[i]=0.00; }
438  VName= VolName;
439  VCopyNo= CopyNo;
440  if(option==1) { VEnter=kTRUE; }
441  else { VEnter=kFALSE; }
442  fPropOption="VE";
443  ProMode=1; //need errors in representation 1 (SC) (see Geane doc)
444  return kTRUE;
445 }
446 
447 Bool_t FairGeanePro::PropagateToLength(Float_t length)
448 {
449  // define final length (option "L")
450  xlf[0]=length;
451  fPropOption="LE";
452  ProMode=1; //need errors in representation 1 (SC)(see Geane doc)
453  return kTRUE;
454 }
455 Bool_t FairGeanePro::PropagateOnlyParameters()
456 {
457  int index = fPropOption.Index("E");
458  if(index == -1) { fPropOption.Append("O"); }
459  else { fPropOption.Replace(index, 1, "O"); }
460  return kTRUE;
461 }
462 
463 Bool_t FairGeanePro::SetWire(TVector3 extremity1, TVector3 extremity2)
464 {
465  // define wires for PCA extrapolation in STT
466  fwire1 = extremity1;
467  fwire2 = extremity2;
468  return kTRUE;
469 }
470 
471 Bool_t FairGeanePro::SetPoint(TVector3 pnt)
472 {
473  // define point for PCA extrapolation in TPC
474  fpoint = pnt;
475  return kTRUE;
476 }
477 
478 Bool_t FairGeanePro::PropagateToPCA(Int_t pca)
479 {
480  // through track length
481  fPropOption="LE";
482  ProMode=1; //need errors in representation 1 (SC)(see Geane doc)
483  fPCA = pca;
484  // initialization
485  fRad = 0.;
486  fDi = 0.;
487  fvpf = TVector3(0.,0.,0.);
488  fvwi = TVector3(0.,0.,0.);
489  ftrklength = 0;
490  return kTRUE;
491 }
492 
493 Bool_t FairGeanePro::PropagateToPCA(Int_t pca, Int_t dir)
494 {
495  // through track length
496  if(dir > 0) { fPropOption="LE"; }
497  else if(dir < 0) { fPropOption="BLE"; }
498  else { cout << "FairGeanePro::PropagateToPCA(int, int) ERROR: no direction set" << endl; }
499  ProMode=1; //need errors in representation 1 (SC)(see Geane doc)
500  fPCA = pca;
501  // initialization
502  fRad = 0.;
503  fDi = 0.;
504  fvpf = TVector3(0.,0.,0.);
505  fvwi = TVector3(0.,0.,0.);
506  ftrklength = 0;
507  return kTRUE;
508 }
509 
510 Bool_t FairGeanePro::ActualFindPCA(Int_t pca, FairTrackParP* par, Int_t dir)
511 {
512  Init(par);
513  // through track length
514  if(dir > 0) { fPropOption="LE"; }
515  else if(dir < 0) { fPropOption="BLE"; }
516  else { cout << "FairGeanePro::ActualFindPCA ERROR: no direction set" << endl; }
517  ProMode=1; //need errors in representation 1 (SC)(see Geane doc)
518  fPCA = pca;
519  // initialization
520  fRad = 0.;
521  fDi = 0.;
522  fvpf = TVector3(0.,0.,0.);
523  fvwi = TVector3(0.,0.,0.);
524  ftrklength = 0;
525  for(Int_t i=0; i<15; i++) { ein[i]=0.00; }
526  return kTRUE;
527 }
528 
529 Bool_t FairGeanePro::BackTrackToVertex()
530 {
531  // through track length
532  fPropOption="BLE";
533  ProMode=1; //need errors in representation 1 (SC)(see Geane doc)
534  fPCA = 1; // to point
535  // initialization (forse non necessario) CHECK!!!!
536  fRad = 0.;
537  fDi = 0.;
538  fvpf = TVector3(0.,0.,0.);
539  fvwi = TVector3(0.,0.,0.);
540  ftrklength = 0;
541  return kTRUE;
542 }
543 
544 Bool_t FairGeanePro::PropagateToVirtualPlaneAtPCA(Int_t pca)
545 {
546  // through track length
547  fPropOption="LE";
548  ProMode=3; //need errors in representation 3 (SD)(see Geane doc)
549  fPCA = pca;
550  // initialization
551  fRad = 0.;
552  fDi = 0.;
553  fvpf = TVector3(0.,0.,0.);
554  fvwi = TVector3(0.,0.,0.);
555  ftrklength = 0;
556  return kTRUE;
557 }
558 
559 Bool_t FairGeanePro::BackTrackToVirtualPlaneAtPCA(Int_t pca)
560 {
561  // through track length
562  fPropOption="BLE";
563  ProMode=3; //need errors in representation 3 (SD)(see Geane doc)
564  fPCA = pca;
565  // initialization
566  fRad = 0.;
567  fDi = 0.;
568  fvpf = TVector3(0.,0.,0.);
569  fvwi = TVector3(0.,0.,0.);
570  ftrklength = 0;
571  return kTRUE;
572 }
573 
574 //=====================
575 
576 int FairGeanePro::FindPCA(Int_t pca, Int_t PDGCode, TVector3 point, TVector3 wire1, TVector3 wire2, Double_t maxdistance, Double_t& Rad, TVector3& vpf, TVector3& vwi, Double_t& Di, Float_t& trklength)
577 {
578  // find the point of closest approach of the track to a point(measured position) or to a line(wire)
579 
580  // INPUT ----------------------------------------
581  // .. pca = ic = 1 closest approach to point
582  // = 2 closest approach to wire
583  // = 0 no closest approach
584  // .. PDGCode = pdg code of the particle
585  // .. point point with respect to which calculate the closest approach
586  // .. wire, wire2 line with respect to which calculate the closest approach
587  // .. maxdistance = geometrical distance[start - point/wire extr] * 2
588  // OUTPUT ----------------------------------------
589  // .. Rad radius if the found circle
590  // .. vpf: point of closest approach on track
591  // .. vwi: point of closest approach on wire
592  // .. Di : distance between track and wire in the PCA
593  // .. trklength : track length to add to the GEANE one
594 
595  Float_t pf[3] = {point.X(), point.Y(), point.Z()};
596  Float_t w1[3] = {wire1.X(), wire1.Y(), wire1.Z()};
597  Float_t w2[3] = {wire2.X(), wire2.Y(), wire2.Z()};
598 
599  GeantCode=fdbPDG->ConvertPdgToGeant3(PDGCode);
600 
601  // flags Rotondi's function
602  int flg=0;
603 
604  // cl track length to the three last points of closest approach
605  // dst assigned distance between initial point in ERTRAK and PFINAL along straight line (currently noy used)
606  Float_t cl[3],dst;
607 
608  // GEANE filled points
609  Float_t po1[3],po2[3],po3[3];
610 
611  // cl track length to the three last points of closest approach
612  Float_t clen[3];
613 
614  // track length to add to GEANE computed one
615  Double_t Le=0.0;
616  Double_t dist1,dist2;
617 
618  // initialization of some variables
619  dst = 999.;
620  cl[0] = 0;
621  cl[1] = 0;
622  cl[2] = 0;
623 
624  // GEANE filled points
625  po1[0]=0;
626  po1[1]=0;
627  po1[2]=0;
628  po2[0]=0;
629  po2[1]=0;
630  po2[2]=0;
631  po3[0]=0;
632  po3[1]=0;
633  po3[2]=0;
634 
635  gMC3->SetClose(pca,pf,dst,w1,w2,po1,po2,po3,cl);
636 
637  // maximum distance calculated 2 * geometric distance
638  // start point - end point (the point to which we want
639  // to find the PCA)
640  Float_t stdlength[1] = {maxdistance};
641 
642  gMC3->Eufill(1, ein, stdlength);
643 
644  //check needed for low momentum tracks
645  gMC3->Ertrak(x1,p1,x2,p2,GeantCode, fPropOption.Data());
646  if(x2[0]<-1.E29) { return 1; }
647  if(gMC3->IsTrackOut()) { return 1; }
648  gMC3->GetClose(po1,po2,po3,clen);
649 
650  // check on cases when only two steps are performed!
651  // in these cases po1[i] = 0 ==> let' s copy po2 into
652  // po1 in order to use only the two actually extrapolated
653  // points po2 and po3 to complete the PCA calculation
654  if(clen[0] == 0 && clen[1] == 0) {
655  po1[0] = po2[0];
656  po1[1] = po2[1];
657  po1[2] = po2[2];
658  }
659 
660  if(pca == 1) {
661  if((po1[0] == po2[0] && po1[1] == po2[1] && po1[2] == po2[2])
662  || (po2[0] == po3[0] && po2[1] == po3[1] && po2[2] == po3[2])) {
663  Int_t quitFlag=0;
664  Track2ToPoint(TVector3(po1),TVector3(po3),TVector3(pf),vpf,Di,Le,quitFlag);
665  if(quitFlag!=0) {std::cerr<<"ABORT"<<std::endl; return 1;} //abort
666  } else {
667  Track3ToPoint(TVector3(po1),TVector3(po2),TVector3(po3),TVector3(pf),vpf,flg,Di,Le,Rad);
668  if(flg==1) {
669  Int_t quitFlag=0;
670  Track2ToPoint(TVector3(po1),TVector3(po3),TVector3(pf),vpf,Di,Le,quitFlag);
671  if(quitFlag!=0) {std::cerr<<"ABORT"<<std::endl; return 1;} //abort
672  } else if(flg==2) {std::cerr<<"ABORT"<<std::endl; return 1;} //abort
673  }
674  // if the propagation to closest approach to a POINT is performed
675  // vwi is the point itself (with respect to which the PCA is calculated)
676  vwi = point;
677  } else if(pca == 2) {
678  if((po1[0] == po2[0] && po1[1] == po2[1] && po1[2] == po2[2])
679  || (po2[0] == po3[0] && po2[1] == po3[1] && po2[2] == po3[2])) {
680  Track2ToLine(TVector3(po1),TVector3(po3),TVector3(w1),
681  TVector3(w2),vpf,vwi,flg,Di,Le);
682  if(flg==1) {
683  dist1 = (vwi-TVector3(w1)).Mag();
684  dist2 = (vwi-TVector3(w2)).Mag();
685  Int_t quitFlag=0;
686  dist1<dist2?Track2ToPoint(TVector3(po1),TVector3(po3),TVector3(w1),vpf,Di,Le,quitFlag):
687  Track2ToPoint(TVector3(po1),TVector3(po3),TVector3(w2),vpf,Di,Le,quitFlag);
688  if(quitFlag!=0) {std::cerr<<"ABORT"<<std::endl; return 1;} //abort
689  } else if(flg==2) {
690  std::cerr<<"ABORT"<<std::endl;
691  return 1;
692  }
693  } else {
694  Track3ToLine(TVector3(po1),TVector3(po2),TVector3(po3),
695  TVector3(w1),TVector3(w2),vpf,vwi,flg,Di,Le,Rad);
696  if(flg==1) {
697  Track2ToLine(TVector3(po1),TVector3(po3),TVector3(w1),
698  TVector3(w2),vpf,vwi,flg,Di,Le);
699  if(flg==1) {
700  dist1 = (vwi-TVector3(w1)).Mag();
701  dist2 = (vwi-TVector3(w2)).Mag();
702  Int_t quitFlag=0;
703  dist1<dist2?Track2ToPoint(TVector3(po1),TVector3(po3),TVector3(w1),vpf,Di,Le,quitFlag):
704  Track2ToPoint(TVector3(po1),TVector3(po3),TVector3(w2),vpf,Di,Le,quitFlag);
705  if(quitFlag!=0) {std::cerr<<"ABORT"<<std::endl; return 1;} //abort
706  } else if(flg==2) {
707  std::cerr<<"ABORT"<<std::endl;
708  return 1;
709  }
710  } else if(flg==2) {
711  dist1 = (vwi-TVector3(w1)).Mag();
712  dist2 = (vwi-TVector3(w2)).Mag();
713 
714  dist1<dist2?Track3ToPoint(TVector3(po1),TVector3(po2),TVector3(po3),TVector3(w1),vpf,flg,Di,Le,Rad):
715  Track3ToPoint(TVector3(po1),TVector3(po2),TVector3(po3),TVector3(w2),vpf,flg,Di,Le,Rad);
716  if(flg==2) {std::cerr<<"ABORT"<<std::endl; return 1;} //abort
717  } else if(flg==3) {
718  dist1 = (vwi-TVector3(w1)).Mag();
719  dist2 = (vwi-TVector3(w2)).Mag();
720  Int_t quitFlag=0;
721  dist1<dist2?Track2ToPoint(TVector3(po1),TVector3(po3),TVector3(w1),vpf,Di,Le,quitFlag):
722  Track2ToPoint(TVector3(po1),TVector3(po3),TVector3(w2),vpf,Di,Le,quitFlag);
723  if(quitFlag!=0) {std::cerr<<"ABORT"<<std::endl; return 1;} //abort
724  } else if(flg==4) {
725  Track2ToLine(TVector3(po1),TVector3(po3),TVector3(w1),
726  TVector3(w2),vpf,vwi,flg,Di,Le);
727  if(flg==2) {
728  std::cerr<<"ABORT"<<std::endl;
729  return 1;
730  }
731  }
732 
733  }
734  }
735 
736  // calculated track length corresponding
737  // to the point of closest approach
738  trklength = clen[0]+Le;
739 
740  // PCA before starting point
741  if(trklength<0) { return 1; }
742  flag = flg;
743 
744  return 0;
745 }
746 
747 
748 void FairGeanePro::Track2ToLine( TVector3 X1, TVector3 X2, TVector3 w1,
749  TVector3 w2, TVector3& Pfinal, TVector3& Pwire,
750  Int_t& Iflag, Double_t& Dist, Double_t& Length)
751 {
752 
753  // Closest approach to a line from 2 GEANE points
754  //
755  // METHOD: the nearest points on the two lines
756  // x1,x2 and w1,w2 is found.
757  // The method is described in:
758  // http://softsurfer.com/Archive/algorithm_0106/algorithm_0106.htm
759  // http://www.geometrictools.com/Documentation/DistanceLine3Line3.pdf.
760  //
761  // INPUT: x1, x2 closest appoach GEANE points
762  // w1, w2 points of the line (wire) to approach
763  //
764  // OUTPUT: Pfinal point of closest approach on the track
765  // Pwire point of closest approach on the wire
766  // Dist distance between Pfian and w1
767  // Length arc length to add to the GEANE length of x1
768  // Iflag =1 when Pwire is outside [w1,w2]
769  // In this case, when w1 and w2 are the extremes
770  // of the wire, the user could remake the procedure
771  // by calling Track3ToPoint or Track2ToPoint, where the
772  // Point is w1 or w2;
773  // = 2 when the two lines are parallel and the solution
774  // does not exists.
775  //
776  // Authors: Andrea Fontana and Alberto Rotondi 20 MAy 2007
777  //
778 
779 
780  TVector3 x21, x32, w21;
781  TVector3 xw1, xw2;
782 
783  Double_t a1, b1, c1, d1, e1, t1, s1;
784  Double_t Delta1;
785 
786  Double_t Eps = 1.E-08;
787 
788  Iflag =0;
789 
790  // line-line distance
791  x21 = X2-X1;
792  w21 = w2-w1;
793 
794  xw1 = X1-w1;
795  xw2 = X2-w1;
796 
797  a1 = x21.Mag2();
798  b1 = x21.Dot(w21);
799  c1 = w21.Mag2();
800  d1 = xw1.Dot(x21);
801  e1 = xw1.Dot(w21);
802 
803  Delta1 = a1*c1-b1*b1;
804 
805  if(Delta1 > Eps) {
806  t1 = (a1*e1-b1*d1)/Delta1;
807  s1 = (b1*e1-c1*d1)/Delta1;
808 
809  Pfinal = (X1 + x21*s1);
810  Pwire = (w1 + w21*t1);
811  Length = s1*x21.Mag();
812  Dist= (Pfinal-Pwire).Mag();
813 
814  } else {
815  // lines are paralllel, no solution does exist
816  Pfinal.SetXYZ(0.,0.,0.);
817  Pwire.SetXYZ(0.,0.,0.);
818  Dist=0.;
819  Length=0.;
820  Iflag = 2;
821  return;
822  }
823  // flag when the point on the wire is outside (w1,w2)
824  if((((Pwire[0]<w1[0] && Pwire[0]<w2[0]) || (w2[0]<Pwire[0] && w1[0]<Pwire[0]))
825  && (fabs(Pwire[0]-w1[0]) > 1e-11 && fabs(Pwire[0]- w2[0]) > 1e-11))
826  || (((Pwire[1]<w1[1] && Pwire[1]<w2[1]) || (w2[1]<Pwire[1] && w1[1]<Pwire[1]))
827  && (fabs(Pwire[1]-w1[1]) > 1e-11 && fabs(Pwire[1]- w2[1]) > 1e-11))
828  || (((Pwire[2]<w1[2] && Pwire[2]<w2[2]) || (w2[2]<Pwire[2] && w1[2]<Pwire[2]))
829  && (fabs(Pwire[2]-w1[2]) > 1e-11 && fabs(Pwire[2]- w2[2]) > 1e-11))) {
830  Iflag=1;
831  }
832 }
833 
834 
835 
836 void FairGeanePro::Track2ToPoint( TVector3 X1, TVector3 X2, TVector3 w1, TVector3& Pfinal,
837  Double_t& Dist, Double_t& Length, Int_t& quitFlag)
838 {
839 
840  //
841  // Closest approach to a point from 2 GEANE points
842  //
843  // METHOD: the nearest point to w1 on the line x1,x2
844  // is found. Elementary vector calculus is used!
845  //
846  // INPUT: x1, x2 closest appoach GEANE points
847  // w1 point to approach w1
848  //
849  // OUTPUT: Pfinal point of closest approach
850  // Dist distance between Pfinal and w1
851  // Length arc length to add to the GEANE length of x1
852  // quitFlag error flag which will be set to 1 if points dont form line
853  // Authors: Andrea Fontana and Alberto Rotondi May 2007
854  //
855 
856  quitFlag = 0;
857 
858  TVector3 u21;
859 
860  Double_t a, t1;
861 
862  // w-point - x-line distance
863  Double_t d=(X2-X1).Mag();
864  if(fabs(d)<1.E-8) {
865  quitFlag=1;
866  return;
867  }
868  a= 1./d;
869 
870  u21 = (X2-X1).Unit();
871 
872  // output
873  Dist= ((w1-X1).Cross(u21)).Mag();
874  t1 = a*(w1-X1).Dot(u21);
875  Pfinal = (X2-X1)*t1 + X1;
876  Length = (X2-X1).Mag()*t1;
877 }
878 
879 
880 void FairGeanePro::Track3ToLine(TVector3 X1, TVector3 X2, TVector3 X3,
881  TVector3 w1, TVector3 w2,
882  // output
883  TVector3& Pfinal, TVector3& Wire,
884  Int_t& Iflag, Double_t& Dist,
885  Double_t& Length, Double_t& Radius)
886 {
887 
888  // Find the closest approach points between a curve (helix)
889  // and a line (wire)
890  //
891  // METHOD:the classical Eberly method is used: see
892  // http://www.geometrictools.com/Documentation/DistanceLine3Circle.pdf.
893  // (see also www.geometrictools.com for the other formulae used
894  // in this interface).
895  // The 4-degree polynomial resulting for the line parameter t is
896  // solved with the efficient SolveQuartic Root routine.
897  // The minimal distance solution is found by using our Track3ToPoint
898  // routine
899  //
900  // INPUT: x1, x2, x3 closest appoach GEANE points
901  // w1, w2 points of the line (wire) to approach
902  //
903  // OUTPUT: Pfinal track point of closest approach
904  // Wire line point of closest approach
905  // Iflag = 1, the points are on a straight line
906  // within the precision of the method (20 micron).
907  // In this case the user should recall Track2ToPoint
908  // = 2, Pwire is outside [w1,w2]
909  // In this case, when w1 and w2 are the extremes
910  // of the wire, the user could remake the procedure
911  // by calling Track3ToPoint, where the Point is w1 or w2.
912  // = 3 both conditions 1 and 2 are encountered
913  // = 4, the method failed for mathematical reasons.
914  // In this case the user should use Track2ToLine where
915  // x1=x1 and x2=x3.
916  // Dist distance between Pfinal and Wire
917  // Length arc length to add to the GEANE length of x1
918  // Radius radius of the found circle
919  //
920  // Authors: Andrea Fontana and Alberto Rotondi 20 June 2007
921  //
922 
923  TVector3 xp1, xp2, xp3, xp32;
924  TVector3 x21, x31;
925  TVector3 e1, e2, e3, aperp;
926  TVector3 Ppfinal, Pwire;
927  TVector3 wp1, wp2, wpt, xR, xpR;
928  TVector3 xw1, xw2;
929  TVector3 px;
930 
931  Double_t T[3][3], TM1[3][3];
932 
933  TVector3 N, M, D, B, Pointw;
934  Double_t a0, a1, b0, b1, c0, c1, c2;
935  Double_t d0, d1, d2, d3, d4, sol4[4], dmin;
936  Double_t Angle;
937  Double_t dx;
938  Int_t it, imin;
939 
940  Iflag = 0;
941 
942  // go to the circle plane: matrix of director cosines
943 
944  x21 = X2-X1;
945  e1 = x21.Unit();
946  T[0][0] = e1.X();
947  T[0][1] = e1.Y();
948  T[0][2] = e1.Z();
949 
950  x31 = X3-X1;
951  e3 = e1.Cross(x31);
952  // if the points are on the same line
953  if(e3.Mag() < 1e-8) {
954  Iflag = 1;
955  return;
956  }
957  e3 = e3.Unit();
958 
959  T[2][0] = e3.X();
960  T[2][1] = e3.Y();
961  T[2][2] = e3.Z();
962 
963  e2 = e3.Cross(e1);
964  T[1][0] = e2.X();
965  T[1][1] = e2.Y();
966  T[1][2] = e2.Z();
967 
968  // new coordinates
969  for(Int_t i=0; i<3; i++) {
970  xp1[i]=0.;
971  xp2[i]=0.;
972  xp3[i]=0.;
973  wp1[i]=0.;
974  wp2[i]=0.;
975  }
976  for(Int_t i=0; i<3; i++) {
977  xp1[i] = 0.;
978  for(Int_t j=0; j<3; j++) {
979  TM1[i][j] = T[j][i];
980  xp2[i] += T[i][j] * (X2[j]-X1[j]);
981  xp3[i] += T[i][j] * (X3[j]-X1[j]);
982  wp1[i] += T[i][j] * (w1[j]-X1[j]);
983  wp2[i] += T[i][j] * (w2[j]-X1[j]);
984  }
985  }
986 
987  // radius and center
988 
989  xp32= xp3 - xp2;
990  xpR[0] = 0.5*xp2[0];
991  if(fabs(xp3[1])<1.E-8) {
992  Iflag = 4;
993  return;
994  }
995  xpR[1] = 0.5*(xp32[0]*xp3[0]/xp3[1]+ xp3[1]);
996  xpR[2] = 0.;
997  Radius = sqrt(pow(xpR[0]-xp1[0],2)+pow(xpR[1]-xp1[1],2));
998 
999  // Eberly's method
1000 
1001  B = wp1;
1002  M = wp2 - wp1;
1003  D = wp1-xpR;
1004  N.SetXYZ(0.,0.,1.);
1005 
1006  a0 = M.Dot(D);
1007  a1 = M.Dot(M);
1008  b0 = M.Dot(D) -(N.Dot(M))*(N.Dot(D));
1009  b1 = M.Dot(M) -(N.Dot(M))*(N.Dot(M));
1010  c0 = D.Dot(D) -(N.Dot(D))*(N.Dot(D));
1011  c1 = b0;
1012  c2 = b1;
1013 
1014  d0 = a0*a0*c0 -b0*b0*Radius*Radius;
1015  d1 = 2.*(a0*a1*c0+a0*a0*c1-b0*b1*Radius*Radius);
1016  d2 = a1*a1*c0+4.*a0*a1*c1+a0*a0*c2-b1*b1*Radius*Radius;
1017  d3 = 2.*(a1*a1*c1+a0*a1*c2);
1018  d4 = a1*a1*c2;
1019 
1020  // solve the quartic equation
1021  for(Int_t k=0; k<4; k++) {
1022  sol4[k] =0.;
1023  }
1024  if(fabs(d4) < 1.E-12) {
1025  Iflag = 4;
1026  return;
1027  }
1028 
1029  TGeoTorus t;
1030  it = t.SolveQuartic(d3/d4,d2/d4,d1/d4,d0/d4,sol4);
1031 
1032  if(it==0) {
1033  Iflag = 4;
1034  return;
1035  }
1036 
1037  // select the right solution
1038  dmin = 1.e+08;
1039  imin=-1;
1040  for(Int_t j=0; j<it; j++) {
1041  Pointw[0] = B[0] + sol4[j]*M[0];
1042  Pointw[1] = B[1] + sol4[j]*M[1];
1043  Pointw[2] = B[2] + sol4[j]*M[2];
1044  Track3ToPoint(xp1,xp2,xp3, Pointw, px, Iflag, dx, Length, Radius);
1045  if(Iflag==2) {
1046  Iflag = 4;
1047  return;
1048  }
1049  if(dx<dmin) {
1050  dmin=dx;
1051  imin=j;
1052  }
1053  }
1054 
1055  // final solution
1056  Pwire[0] = B[0] + sol4[imin]*M[0];
1057  Pwire[1] = B[1] + sol4[imin]*M[1];
1058  Pwire[2] = B[2] + sol4[imin]*M[2];
1059  Track3ToPoint(xp1,xp2,xp3, Pwire, px, Iflag, dx, Length, Radius);
1060  if(Iflag==2) {
1061  Iflag = 4;
1062  return;
1063  }
1064 
1065  // output: distance and points in the circle plane reference
1066 
1067  Dist = dx;
1068  Ppfinal = px;
1069 
1070  //
1071  // back to lab coordinates
1072  //
1073 
1074  xR[0]=0.;
1075  xR[1]=0.;
1076  xR[2]=0.;
1077  Pfinal[0]=0.;
1078  Pfinal[1]=0.;
1079  Pfinal[2]=0.;
1080  Wire[0]=0.;
1081  Wire[1]=0.;
1082  Wire[2]=0.;
1083 
1084  for(Int_t i=0; i<3; i++) {
1085  for(Int_t j=0; j<3; j++) {
1086  Pfinal[i] += TM1[i][j] * Ppfinal[j];
1087  Wire[i] += TM1[i][j] * Pwire[j];
1088  xR[i] += TM1[i][j] * xpR[j];
1089  }
1090  }
1091  Pfinal = Pfinal+X1;
1092  Wire = Wire + X1;
1093  xR = xR + X1;
1094 
1095  double dx1=(X1-xR).Mag();
1096  double dx2=(Pfinal-xR).Mag();
1097  double dx12=dx1*dx2;
1098  if(fabs(dx12)<1.E-8) {
1099  Iflag = 4;
1100  return;
1101  }
1102 
1103  // now find the length
1104  Angle = TMath::ACos((X1-xR).Dot(Pfinal-xR)/(dx12));
1105  Length = Radius*Angle;
1106  if((X2-X1).Dot(Pfinal-X1) < 0.) { Length = -Length; }
1107 
1108  // flag straight points within 20 microns
1109  Double_t epsi=0;
1110  if(Radius>1E-8) { epsi = Radius*(1.-TMath::Cos(0.5*(X3-X1).Mag()/Radius)); }
1111  if(epsi < 0.0020) { Iflag=1; }
1112 
1113  // flag when the point on the wire is outside (w1,w2)
1114  if((((Wire[0]<w1[0] && Wire[0]<w2[0]) || (w2[0]<Wire[0] && w1[0]<Wire[0]))
1115  && (fabs(Wire[0]-w1[0]) > 1e-11 && fabs(Wire[0]- w2[0]) > 1e-11))
1116  || (((Wire[1]<w1[1] && Wire[1]<w2[1]) || (w2[1]<Wire[1] && w1[1]<Wire[1]))
1117  && (fabs(Wire[1]-w1[1]) > 1e-11 && fabs(Wire[1]- w2[1]) > 1e-11))
1118  || (((Wire[2]<w1[2] && Wire[2]<w2[2]) || (w2[2]<Wire[2] && w1[2]<Wire[2]))
1119  && (fabs(Wire[2]-w1[2]) > 1e-11 && fabs(Wire[2]- w2[2]) > 1e-11))) {
1120  Iflag=2;
1121  }
1122 }
1123 
1124 void FairGeanePro::Track3ToPoint( TVector3 X1, TVector3 X2, TVector3 X3, TVector3 w1,
1125  // output
1126  TVector3& Pfinal, Int_t& Iflag,
1127  Double_t& Dist, Double_t& Length, Double_t& Radius)
1128 {
1129 
1130  // Closest approach to a point from 3 GEANE points
1131  //
1132  // METHOD: first, we go on the circle plane to have
1133  // x1=(0,0), x2=(x2,0), x3=(x3,y3).
1134  // Then, the point on the circle is found as the
1135  // intersection between the circle and the line
1136  // joining the circle center and the projection of
1137  // the point on the circle plane.
1138  // The 3D distance is found between w1 and this
1139  // point on the circle.
1140  //
1141  // INPUT: x1, x2, x3 closest appoach GEANE points
1142  // w1 point to approach
1143  //
1144  // OUTPUT: Pfinal point of closest approach on the track
1145  // Dist distance between Pfinal and w1
1146  // Length arc length to add to the GEANE length of x1
1147  // Iflag when =1 the points are on a straight line
1148  // within the precision of the method (20 micron).
1149  // In this case the user should recall Track2ToPoint
1150  // if =2 mathematical errors
1151  // Authors: Andrea Fontana and Alberto Rotondi 20 May 2007
1152  //
1153 
1154 
1155  TVector3 xp1, xp2, xp3, xp32;
1156  TVector3 x21, x31;
1157  TVector3 e1, e2, e3;
1158  TVector3 Ppfinal, x32;
1159  TVector3 wp1, wpt, xR, xpR;
1160  TVector3 xw1, xw2;
1161 
1162  TVector3 xc1, xc2, xc3, wc1;
1163 
1164  Double_t m1, m3, Rt;
1165  Double_t T[3][3], TM1[3][3];
1166 
1167  Double_t Angle;
1168 
1169 
1170  Iflag = 0;
1171 
1172  // go to the circle plane with origin in x1 prime (xp1):
1173  // matrix of director cosines
1174 
1175  x21 = X2-X1;
1176 
1177  Double_t x21mag=x21.Mag();
1178  if(x21mag<1.E-8) {
1179  Iflag=2;
1180  return;
1181  }
1182 
1183  m1 = 1./x21mag;
1184  e1 = m1*x21;
1185  T[0][0] = e1.X();
1186  T[0][1] = e1.Y();
1187  T[0][2] = e1.Z();
1188 
1189  x31 = X3-X1;
1190  e3 = e1.Cross(x31);
1191 
1192  // if the points are on the same line
1193  if(e3.Mag() < 1e-8) {
1194  Iflag = 1;
1195  return;
1196  }
1197 
1198  m3 = 1./e3.Mag();
1199  e3 = m3*e3;
1200  T[2][0] = e3.X();
1201  T[2][1] = e3.Y();
1202  T[2][2] = e3.Z();
1203 
1204  e2 = e3.Cross(e1);
1205  T[1][0] = e2.X();
1206  T[1][1] = e2.Y();
1207  T[1][2] = e2.Z();
1208 
1209  // new coordinates
1210  for(Int_t i=0; i<3; i++) {
1211  xp1[i]=0.;
1212  xp2[i]=0.;
1213  xp3[i]=0.;
1214  wp1[i]=0.;
1215  }
1216  for(Int_t i=0; i<3; i++) {
1217  for(Int_t j=0; j<3; j++) {
1218  TM1[i][j] = T[j][i];
1219  xp1[i] += 0.;
1220  xp2[i] += T[i][j] * (X2[j]-X1[j]);
1221  xp3[i] += T[i][j] * (X3[j]-X1[j]);
1222  wp1[i] += T[i][j] * (w1[j]-X1[j]);
1223  }
1224  }
1225 
1226  // radius Radius and center xpR
1227 
1228  xp32= xp3 - xp2;
1229  xpR[0] = 0.5*xp2[0];
1230  if(fabs(xp3[1])<1.E-8) {
1231  Iflag = 2;
1232  return;
1233  }
1234  xpR[1] = 0.5*(xp32[0]*xp3[0]/xp3[1]+ xp3[1]);
1235  xpR[2] = 0.;
1236 
1237  Radius = sqrt( pow(xpR[0]-xp1[0],2) + pow(xpR[1]-xp1[1],2) );
1238 
1239  // distance and points
1240  wpt = wp1;
1241  wpt[2] =0.; // point projection on the circle plane
1242 
1243  Double_t dwp=(wpt-xpR).Mag();
1244  if(fabs(dwp)<1.E-8) {
1245  Iflag = 2;
1246  return;
1247  }
1248  Rt = Radius/dwp;
1249  Ppfinal = (wpt-xpR)*Rt + xpR;
1250  Dist = (wp1-Ppfinal).Mag();
1251 
1252  // back to lab coordinates:
1253  //from Ppfinal to Pfinal and from xpR to xR
1254 
1255  xR[0]=0.;
1256  xR[1]=0.;
1257  xR[2]=0.;
1258  Pfinal[0]=0.;
1259  Pfinal[1]=0.;
1260  Pfinal[2]=0.;
1261  for(Int_t i=0; i<3; i++) {
1262  for(Int_t j=0; j<3; j++) {
1263  Pfinal[i] += TM1[i][j] * Ppfinal[j];
1264  xR[i] += TM1[i][j] * xpR[j];
1265  }
1266  }
1267  Pfinal = Pfinal+X1;
1268  xR = xR +X1;
1269 
1270  // now find the length
1271  double dx1=(X1-xR).Mag();
1272  double dx2=(Pfinal-xR).Mag();
1273  double dx12=dx1*dx2;
1274  if(fabs(dx12)<1.E-8) {
1275  Iflag = 2;
1276  return;
1277  }
1278  // now find the length
1279  Angle = TMath::ACos((X1-xR).Dot(Pfinal-xR)/(dx12));
1280  Length = Radius*Angle;
1281 
1282  // flag straight points within 20 microns
1283  Double_t epsi=0;
1284  if(Radius>1E-8) { epsi = Radius*(1.-TMath::Cos(0.5*(X3-X1).Mag()/Radius)); }
1285  if(epsi < 0.0020) { Iflag=1; }
1286 }
1287 
1288 void FairGeanePro::GetTransportMatrix(Double_t trm[5][5])
1289 {
1290  for(int i = 0; i < 5; i++) for(int j = 0; j < 5; j++) { trm[i][j] = trpmat[i][j]; }
1291 }
1292 
1293 ClassImp(FairGeanePro)
1294 
1295