Bring AliFemto up to date with latest code developements
[u/mrichter/AliRoot.git] / PWG2 / FEMTOSCOPY / AliFemto / AliFemtoParticle.cxx
1 ///////////////////////////////////////////////////////////////////////////
2 //                                                                       //
3 // AliFemtoParticle: main class halding all the necessary information    //
4 // about particle that is required during femtoscopic analysis           //
5 // This includes all the information about the quality of the track,     //
6 // its identification as well as track chracteristics with connection    //
7 // to the detector parts, e.g. entrance and exit points.                 //
8 //                                                                       //
9 ///////////////////////////////////////////////////////////////////////////
10 #include "AliFemtoParticle.h"
11 //#include "math_constants.h"
12 #ifdef __CC5__
13   #include <math.h>
14 #else
15   #include <cmath>
16 #endif
17
18
19 #include "TMath.h"
20 using namespace TMath;
21
22 double AliFemtoParticle::fgPrimPimPar0= 9.05632e-01;
23 double AliFemtoParticle::fgPrimPimPar1= -2.26737e-01;
24 double AliFemtoParticle::fgPrimPimPar2= -1.03922e-01;
25 double AliFemtoParticle::fgPrimPipPar0= 9.09616e-01;
26 double AliFemtoParticle::fgPrimPipPar1= -9.00511e-02;
27 double AliFemtoParticle::fgPrimPipPar2= -6.02940e-02;
28 double AliFemtoParticle::fgPrimPmPar0= 0.;
29 double AliFemtoParticle::fgPrimPmPar1= 0.;
30 double AliFemtoParticle::fgPrimPmPar2= 0.;
31 double AliFemtoParticle::fgPrimPpPar0= 0.;
32 double AliFemtoParticle::fgPrimPpPar1= 0.;
33 double AliFemtoParticle::fgPrimPpPar2= 0.;
34
35 int TpcLocalTransform(AliFmThreeVectorD& xgl, 
36                       int& iSector, 
37                       int& iPadrow, 
38                       float& xlocal,
39                       double& ttPhi);
40
41
42 //_____________________
43 AliFemtoParticle::AliFemtoParticle() : 
44 //   fTpcV0NegPosSample(0),
45 //   fV0NegZ(0),
46 //   fV0NegU(0),
47 //   fV0NegSect(0),
48   fTrack(0), fV0(0), fKink(0), fXi(0), 
49   fFourMomentum(0),
50   fHelix(),
51 //   fNominalTpcExitPoint(0),
52 //   fNominalTpcEntrancePoint(0),
53   fHiddenInfo(0),
54   fPrimaryVertex(0),
55   fSecondaryVertex(0),
56   fHelixV0Pos(),
57   fTpcV0PosEntrancePoint(0),
58   fTpcV0PosExitPoint(0),
59   fHelixV0Neg(),
60   fTpcV0NegEntrancePoint(0),
61   fTpcV0NegExitPoint(0)
62 {
63   // Default constructor
64   /* no-op for default */
65   //  cout << "Created particle " << this << endl;
66 }
67 //_____________________
68 AliFemtoParticle::AliFemtoParticle(const AliFemtoParticle& aParticle):
69 //   fTpcV0NegPosSample(0),
70 //   fV0NegZ(0),
71 //   fV0NegU(0),
72 //   fV0NegSect(0),
73   fTrack(0), fV0(0), fKink(0), fXi(0),
74   fFourMomentum(0),
75   fHelix(),
76 //   fNominalTpcExitPoint(0),
77 //   fNominalTpcEntrancePoint(0),
78   fHiddenInfo(0), 
79   fPrimaryVertex(0),
80   fSecondaryVertex(0),
81   fHelixV0Pos(),
82   fTpcV0PosEntrancePoint(0),
83   fTpcV0PosExitPoint(0),
84   fHelixV0Neg(),
85   fTpcV0NegEntrancePoint(0),
86   fTpcV0NegExitPoint(0)
87 {
88   // Copy constructor
89   if (aParticle.fTrack)
90     fTrack = new AliFemtoTrack(*aParticle.fTrack);
91   if (aParticle.fV0)
92     fV0    = new AliFemtoV0(*aParticle.fV0);
93   if (aParticle.fKink)
94     fKink  = new AliFemtoKink(*aParticle.fKink);
95   if (aParticle.fXi)
96     fXi    = new AliFemtoXi(*aParticle.fXi);
97   fFourMomentum = aParticle.fFourMomentum;
98   fHelix = aParticle.fHelix;
99
100 //   for (int iter=0; iter<11; iter++)
101 //     fNominalPosSample[iter] = aParticle.fNominalPosSample[iter];
102
103 //   if (aParticle.fTpcV0NegPosSample) {
104 //     fTpcV0NegPosSample = (AliFemtoThreeVector *) malloc(sizeof(AliFemtoThreeVector) * 11);
105 //     for (int iter=0; iter<11; iter++)
106 //       fTpcV0NegPosSample[iter] = aParticle.fTpcV0NegPosSample[iter];
107 //   }
108
109 //   if (aParticle.fV0NegZ) {
110 //     fV0NegZ = (float *) malloc(sizeof(float) * 45);
111 //     for (int iter=0; iter<11; iter++)
112 //       fV0NegZ[iter] = aParticle.fV0NegZ[iter];
113 //   }
114 //   if (aParticle.fV0NegU) {
115 //     fV0NegU = (float *) malloc(sizeof(float) * 45);
116 //     for (int iter=0; iter<11; iter++)
117 //       fV0NegU[iter] = aParticle.fV0NegU[iter];
118 //   }
119 //   if (aParticle.fV0NegSect) {
120 //     fV0NegSect = (int *) malloc(sizeof(int) * 45);
121 //     for (int iter=0; iter<11; iter++)
122 //       fV0NegSect[iter] = aParticle.fV0NegSect[iter];
123 //   }
124
125   fPrimaryVertex = aParticle.fPrimaryVertex;
126   fSecondaryVertex = aParticle.fSecondaryVertex;
127   CalculatePurity();
128   if(aParticle.fHiddenInfo){
129     fHiddenInfo= aParticle.HiddenInfo()->Clone();
130   }
131   
132 //   fNominalTpcEntrancePoint = aParticle.fNominalTpcEntrancePoint;
133 //   fNominalTpcExitPoint     = aParticle.fNominalTpcExitPoint;
134   
135   for (int iter=0; iter<6; iter++)
136     fPurity[iter] = aParticle.fPurity[iter];
137   
138   fHelixV0Pos = aParticle.fHelixV0Pos;
139   fTpcV0PosEntrancePoint = aParticle.fTpcV0PosEntrancePoint;
140   fTpcV0PosExitPoint     = aParticle.fTpcV0PosExitPoint;
141   fHelixV0Neg = aParticle.fHelixV0Neg;
142   fTpcV0NegEntrancePoint = aParticle.fTpcV0NegEntrancePoint;
143   fTpcV0NegExitPoint     = aParticle.fTpcV0NegExitPoint;
144 }
145 //_____________________
146 AliFemtoParticle::~AliFemtoParticle(){
147   //  cout << "Issuing delete for AliFemtoParticle." << endl;
148
149   if (fTrack) delete fTrack;
150   if (fV0) {
151 //     delete[] fTpcV0NegPosSample;
152 //     delete[] fV0NegZ;
153 //     delete[] fV0NegU;
154 //     delete[] fV0NegSect;
155     delete fV0;
156   }
157   if (fKink) delete fKink;
158   //  cout << "Trying to delete HiddenInfo: " << fHiddenInfo << endl;
159   if (fHiddenInfo) 
160     {
161       //      cout << "Deleting HiddenInfo." << endl;
162       delete fHiddenInfo;
163     }
164   //  cout << "Deleted particle " << this << endl;
165 }
166 //_____________________
167 AliFemtoParticle::AliFemtoParticle(const AliFemtoTrack* const hbtTrack,const double& mass) : 
168 //   fTpcV0NegPosSample(0),
169 //   fV0NegZ(0),
170 //   fV0NegU(0),
171 //   fV0NegSect(0),
172   fTrack(0), fV0(0), fKink(0), fXi(0),
173   fFourMomentum(0),
174   fHelix(),
175 //   fNominalTpcExitPoint(0),
176 //   fNominalTpcEntrancePoint(0),
177   fHiddenInfo(0), 
178   fPrimaryVertex(0),
179   fSecondaryVertex(0),
180   fHelixV0Pos(),
181   fTpcV0PosEntrancePoint(0),
182   fTpcV0PosExitPoint(0),
183   fHelixV0Neg(),
184   fTpcV0NegEntrancePoint(0),
185   fTpcV0NegExitPoint(0)
186 {
187   // Constructor from normal track
188   
189   // I know there is a better way to do this...
190   fTrack = new AliFemtoTrack(*hbtTrack);
191   AliFemtoThreeVector temp = hbtTrack->P();
192   fFourMomentum.setVect(temp);
193   double ener = ::sqrt(temp.mag2()+mass*mass);
194   fFourMomentum.setE(ener);
195 //  fMap[0] = hbtTrack->TopologyMap(0);
196  // fMap[1] = hbtTrack->TopologyMap(1);
197  // fNhits = hbtTrack->NHits();
198   fHelix = hbtTrack->Helix();
199   //CalculateNominalTpcExitAndEntrancePoints();
200
201  
202   fPrimaryVertex.setX(0.);
203   fPrimaryVertex.setY(0.);
204   fPrimaryVertex.setZ(0.);
205   fSecondaryVertex.setX(0.);
206   fSecondaryVertex.setY(0.);
207   fSecondaryVertex.setZ(0.);
208   /* TO JA ODZNACZYLEM NIE WIEM DLACZEGO
209   CalculateTpcExitAndEntrancePoints(&fHelix,&fPrimaryVertex,
210                                     &fSecondaryVertex,
211                                     &fNominalTpcEntrancePoint,
212                                     &fNominalTpcExitPoint,
213                                     &mNominalPosSample[0],
214                                     &fZ[0],
215                                     &fU[0],
216                                     &fSect[0]);
217   */
218   CalculatePurity();
219   // ***
220   fHiddenInfo= 0;
221   if(hbtTrack->ValidHiddenInfo()){
222     fHiddenInfo= hbtTrack->GetHiddenInfo()->Clone();
223   }
224   // ***
225   //  cout << "Created particle " << this << endl;
226
227 }
228 //_____________________
229 AliFemtoParticle::AliFemtoParticle(const AliFemtoV0* const hbtV0,const double& mass) : 
230 //   fTpcV0NegPosSample(0),
231 //   fV0NegZ(0),
232 //   fV0NegU(0),
233 //   fV0NegSect(0),
234   fTrack(0), fV0(0), fKink(0),  fXi(0),
235   fFourMomentum(0),
236   fHelix(),
237 //   fNominalTpcExitPoint(0),
238 //   fNominalTpcEntrancePoint(0),
239   fHiddenInfo(0),
240   fPrimaryVertex(0),
241   fSecondaryVertex(0),
242   fHelixV0Pos(),
243   fTpcV0PosEntrancePoint(0),
244   fTpcV0PosExitPoint(0),
245   fHelixV0Neg(),
246   fTpcV0NegEntrancePoint(0),
247   fTpcV0NegExitPoint(0)
248 {
249   // Constructor from V0
250   fV0 = new AliFemtoV0(*hbtV0);
251  //fMap[0]= 0;
252   //fMap[1]= 0;
253   // I know there is a better way to do this...
254   AliFemtoThreeVector temp = hbtV0->MomV0();
255   fFourMomentum.setVect(temp);
256   double ener = ::sqrt(temp.mag2()+mass*mass);
257   fFourMomentum.setE(ener);
258   // Calculating TpcEntrancePoint for Positive V0 daugther
259   fPrimaryVertex = hbtV0->PrimaryVertex();
260   fSecondaryVertex = hbtV0->DecayVertexV0();
261   fHelixV0Pos = hbtV0->HelixPos();
262
263 //   fTpcV0NegPosSample = new AliFemtoThreeVector[45];//for V0Neg
264 //   fV0NegZ = new float[45];//for V0Neg
265 //   fV0NegU = new float[45];//for V0Neg
266 //   fV0NegSect = new int[45];//for V0Neg
267 //   CalculateTpcExitAndEntrancePoints(&fHelixV0Pos,&fPrimaryVertex,
268 //                                  &fSecondaryVertex,
269 //                                  &fTpcV0PosEntrancePoint,
270 //                                  &fTpcV0PosExitPoint,
271 //                                  &fNominalPosSample[0],
272 //                                  &fZ[0],
273 //                                  &fU[0],&fSect[0]);
274   fHelixV0Neg = hbtV0->HelixNeg();
275
276 //   CalculateTpcExitAndEntrancePoints(&fHelixV0Neg,
277 //                                  &fPrimaryVertex,
278 //                                  &fSecondaryVertex,
279 //                                  &fTpcV0NegEntrancePoint,
280 //                                  &fTpcV0NegExitPoint,
281 //                                  &fTpcV0NegPosSample[0],
282 //                                  &fV0NegZ[0],
283 //                                  &fV0NegU[0],&fV0NegSect[0]);
284
285   // ***
286   fHiddenInfo= 0;
287   if(hbtV0->ValidHiddenInfo()){
288     fHiddenInfo= hbtV0->GetHiddenInfo()->Clone();
289   }
290   // ***
291 }
292 //_____________________
293 AliFemtoParticle::AliFemtoParticle(const AliFemtoKink* const hbtKink,const double& mass) : 
294 //   fTpcV0NegPosSample(0),
295 //   fV0NegZ(0),
296 //   fV0NegU(0),
297 //   fV0NegSect(0),
298   fTrack(0), fV0(0), fKink(0), fXi(0),
299   fFourMomentum(0),
300   fHelix(),
301 //   fNominalTpcExitPoint(0),
302 //   fNominalTpcEntrancePoint(0),
303   fHiddenInfo(0),
304   fPrimaryVertex(0),
305   fSecondaryVertex(0),
306   fHelixV0Pos(),
307   fTpcV0PosEntrancePoint(0),
308   fTpcV0PosExitPoint(0),
309   fHelixV0Neg(),
310   fTpcV0NegEntrancePoint(0),
311   fTpcV0NegExitPoint(0)
312 {
313   // Constructor from Kink
314   fKink = new AliFemtoKink(*hbtKink);
315  // fMap[0]= 0;
316   //fMap[1]= 0;
317   // I know there is a better way to do this...
318   AliFemtoThreeVector temp = hbtKink->Parent().P();
319   fFourMomentum.setVect(temp);
320   double ener = ::sqrt(temp.mag2()+mass*mass);
321   fFourMomentum.setE(ener);
322 }
323
324 //_____________________
325 AliFemtoParticle::AliFemtoParticle(const AliFemtoXi* const hbtXi, const double& mass) :
326 //   fTpcV0NegPosSample(0),
327 //   fV0NegZ(0),
328 //   fV0NegU(0),
329 //   fV0NegSect(0),
330   fTrack(0), fV0(0), fKink(0), fXi(0),
331   fFourMomentum(0),
332   fHelix(),
333 //   fNominalTpcExitPoint(0),
334 //   fNominalTpcEntrancePoint(0),
335   fHiddenInfo(0), 
336   fPrimaryVertex(0),
337   fSecondaryVertex(0),
338   fHelixV0Pos(),
339   fTpcV0PosEntrancePoint(0),
340   fTpcV0PosExitPoint(0),
341   fHelixV0Neg(),
342   fTpcV0NegEntrancePoint(0),
343   fTpcV0NegExitPoint(0)
344 {
345   // Constructor from Xi
346   fXi = new AliFemtoXi(*hbtXi);
347  // fMap[0]= 0;
348   //fMap[1]= 0;
349   AliFemtoThreeVector temp;// = hbtXi->mMofXi;
350   fFourMomentum.setVect(temp);
351   double ener = ::sqrt(temp.mag2()+mass*mass);
352   fFourMomentum.setE(ener);
353   fHiddenInfo = 0;
354 }
355 //_____________________
356 AliFemtoParticle& AliFemtoParticle::operator=(const AliFemtoParticle& aParticle)
357 {
358   // assignment operator
359   if (this == &aParticle)
360     return *this;
361
362   if (aParticle.fTrack)
363     fTrack = new AliFemtoTrack(*aParticle.fTrack);
364   if (aParticle.fV0)
365     fV0    = new AliFemtoV0(*aParticle.fV0);
366   if (aParticle.fKink)
367     fKink  = new AliFemtoKink(*aParticle.fKink);
368   if (aParticle.fXi)
369     fXi    = new AliFemtoXi(*aParticle.fXi);
370   fFourMomentum = aParticle.fFourMomentum;
371   fHelix = aParticle.fHelix;
372
373 //   for (int iter=0; iter<11; iter++)
374 //     fNominalPosSample[iter] = aParticle.fNominalPosSample[iter];
375
376 //   if (fTpcV0NegPosSample) delete fTpcV0NegPosSample;
377 //   if (aParticle.fTpcV0NegPosSample) {
378 //     fTpcV0NegPosSample = (AliFemtoThreeVector *) malloc(sizeof(AliFemtoThreeVector) * 11);
379 //     for (int iter=0; iter<11; iter++)
380 //       fTpcV0NegPosSample[iter] = aParticle.fTpcV0NegPosSample[iter];
381 //   }
382
383 //   if (fV0NegZ) delete fV0NegZ;
384 //   if (aParticle.fV0NegZ) {
385 //     fV0NegZ = (float *) malloc(sizeof(float) * 45);
386 //     for (int iter=0; iter<11; iter++)
387 //       fV0NegZ[iter] = aParticle.fV0NegZ[iter];
388 //   }
389 //   if (fV0NegU) delete fV0NegU;
390 //   if (aParticle.fV0NegU) {
391 //     fV0NegU = (float *) malloc(sizeof(float) * 45);
392 //     for (int iter=0; iter<11; iter++)
393 //       fV0NegU[iter] = aParticle.fV0NegU[iter];
394 //   }
395 //   if (fV0NegSect) delete fV0NegSect;
396 //   if (aParticle.fV0NegSect) {
397 //     fV0NegSect = (int *) malloc(sizeof(int) * 45);
398 //     for (int iter=0; iter<11; iter++)
399 //       fV0NegSect[iter] = aParticle.fV0NegSect[iter];
400 //   }
401
402   fPrimaryVertex = aParticle.fPrimaryVertex;
403   fSecondaryVertex = aParticle.fSecondaryVertex;
404   CalculatePurity();
405   if(aParticle.fHiddenInfo){
406     fHiddenInfo= aParticle.fHiddenInfo->Clone();
407   }
408   
409 //   fNominalTpcEntrancePoint = aParticle.fNominalTpcEntrancePoint;
410 //   fNominalTpcExitPoint     = aParticle.fNominalTpcExitPoint;
411  
412   if (fHiddenInfo) delete fHiddenInfo;
413   if (aParticle.fHiddenInfo) 
414     fHiddenInfo = aParticle.HiddenInfo()->Clone();
415   
416   for (int iter=0; iter<6; iter++)
417     fPurity[iter] = aParticle.fPurity[iter];
418   
419   fHelixV0Pos = aParticle.fHelixV0Pos;
420   fTpcV0PosEntrancePoint = aParticle.fTpcV0PosEntrancePoint;
421   fTpcV0PosExitPoint     = aParticle.fTpcV0PosExitPoint;
422   fHelixV0Neg = aParticle.fHelixV0Neg;
423   fTpcV0NegEntrancePoint = aParticle.fTpcV0NegEntrancePoint;
424   fTpcV0NegExitPoint     = aParticle.fTpcV0NegExitPoint;
425
426   return *this;
427 }
428 // //_____________________
429 // const AliFemtoThreeVector& AliFemtoParticle::NominalTpcExitPoint() const{
430 //   // in future, may want to calculate this "on demand" only, sot this routine may get more sophisticated
431 //   // for now, we calculate Exit and Entrance points upon instantiation
432 //   return fNominalTpcExitPoint;
433 // }
434 // //_____________________
435 // const AliFemtoThreeVector& AliFemtoParticle::NominalTpcEntrancePoint() const{
436 //   // in future, may want to calculate this "on demand" only, sot this routine may get more sophisticated
437 //   // for now, we calculate Exit and Entrance points upon instantiation
438 //   return fNominalTpcEntrancePoint;
439 // }
440 //_____________________
441 void AliFemtoParticle::CalculatePurity(){
442   // Calculate additional parameterized purity
443
444   double tPt = fFourMomentum.perp();
445   // pi -
446   fPurity[0] = fgPrimPimPar0*(1.-exp((tPt-fgPrimPimPar1)/fgPrimPimPar2));
447   fPurity[0] *= fTrack->PidProbPion();
448   // pi+
449   fPurity[1] = fgPrimPipPar0*(1.-exp((tPt-fgPrimPipPar1)/fgPrimPipPar2));
450   fPurity[1] *= fTrack->PidProbPion();
451   // K-
452   fPurity[2] = fTrack->PidProbKaon();
453   // K+
454   fPurity[3] = fTrack->PidProbKaon();
455   // pbar
456   fPurity[4] = fTrack->PidProbProton();
457   // p
458   fPurity[5] = fTrack->PidProbProton();
459 }
460
461 double AliFemtoParticle::GetPionPurity()
462 {
463   // Get full pion purity
464   if (fTrack->Charge()>0)
465     return fPurity[1];
466   else
467     return fPurity[0];
468 }
469 double AliFemtoParticle::GetKaonPurity()
470 {
471   // Get full kaon purity
472   if (fTrack->Charge()>0)
473     return fPurity[3];
474   else
475     return fPurity[2];
476 }
477 double AliFemtoParticle::GetProtonPurity()
478 {
479   // Get full proton purity
480   if (fTrack->Charge()>0)
481     return fPurity[5];
482   else
483     return fPurity[4];
484 }
485
486 // void AliFemtoParticle::CalculateTpcExitAndEntrancePoints(AliFmPhysicalHelixD* tHelix,
487 //                                                     AliFemtoThreeVector*  PrimVert,
488 //                                                     AliFemtoThreeVector*  SecVert,
489 //                                                     AliFemtoThreeVector* tmpTpcEntrancePoint,
490 //                                                     AliFemtoThreeVector* tmpTpcExitPoint,
491 //                                                     AliFemtoThreeVector* tmpPosSample,
492 //                                                     float* tmpZ,
493 //                                                     float* tmpU,
494 //                                                     int* tmpSect){
495 //   // this calculates the exit point of a secondary track, 
496 //   // either through the endcap or through the Outer Field Cage
497 //   // We assume the track to start at tHelix.origin-PrimaryVertex
498 //   // it also calculates the entrance point of the secondary track, 
499 //   // which is the point at which it crosses the
500 //   // inner field cage
501 //   //  static AliFemtoThreeVector ZeroVec(0.,0.,0.);
502 //   AliFemtoThreeVector tZeroVec(0.,0.,0.);
503 // //   tZeroVec.setX(tHelix->origin().x()-PrimVert->x());
504 // //   tZeroVec.setY(tHelix->origin().y()-PrimVert->y());
505 // //   tZeroVec.setZ(tHelix->origin().z()-PrimVert->z());
506 //   tZeroVec.setX(SecVert->x()-PrimVert->x());
507 //   tZeroVec.setY(SecVert->y()-PrimVert->y());
508 //   tZeroVec.setZ(SecVert->z()-PrimVert->z());
509 //   double dip, curv, phase;
510 //   int h;
511 //   curv = tHelix->Curvature();
512 //   dip  = tHelix->DipAngle();
513 //   phase= tHelix->Phase();
514 //   h    = tHelix->H();
515   
516 //   AliFmHelixD hel(curv,dip,phase,tZeroVec,h);
517
518 //   pairD candidates;
519 //   double sideLength;  // this is how much length to go to leave through sides of TPC
520 //   double endLength;  // this is how much length to go to leave through endcap of TPC
521 //   // figure out how far to go to leave through side...
522 //   candidates = hel.PathLength(200.0);  // bugfix MAL jul00 - 200cm NOT 2cm
523 //   sideLength = (candidates.first > 0) ? candidates.first : candidates.second;
524
525 //   static AliFemtoThreeVector tWestEnd(0.,0.,200.);  // bugfix MAL jul00 - 200cm NOT 2cm
526 //   static AliFemtoThreeVector tEastEnd(0.,0.,-200.); // bugfix MAL jul00 - 200cm NOT 2cm
527 //   static AliFemtoThreeVector tEndCapNormal(0.,0.,1.0);
528
529 //   endLength = hel.PathLength(tWestEnd,tEndCapNormal);
530 //   if (endLength < 0.0) endLength = hel.PathLength(tEastEnd,tEndCapNormal);
531
532 //   if (endLength < 0.0) cout << 
533 //                       "AliFemtoParticle::CalculateTpcExitAndEntrancePoints(): "
534 //                             << "Hey -- I cannot find an exit point out endcaps" << endl;
535 //   // OK, firstExitLength will be the shortest way out of the detector...
536 //   double firstExitLength = (endLength < sideLength) ? endLength : sideLength;
537 //   // now then, let's return the POSITION at which particle leaves TPC...
538 //   *tmpTpcExitPoint = hel.At(firstExitLength);
539 //   // Finally, calculate the position at which the track crosses the inner field cage
540 //   candidates = hel.PathLength(50.0);  // bugfix MAL jul00 - 200cm NOT 2cm
541
542 //   sideLength = (candidates.first > 0) ? candidates.first : candidates.second;
543 //   //  cout << "sideLength 2 ="<<sideLength << endl;
544 //   *tmpTpcEntrancePoint = hel.At(sideLength);
545 //   // This is the secure way !  
546 //   if (IsNaN(tmpTpcEntrancePoint->x()) || 
547 //       IsNaN(tmpTpcEntrancePoint->y()) || 
548 //       IsNaN(tmpTpcEntrancePoint->z()) ){ 
549 //     cout << "tmpTpcEntrancePoint NAN"<< endl; 
550 //     cout << "tmpNominalTpcEntrancePoint = " <<tmpTpcEntrancePoint<< endl;
551 //     tmpTpcEntrancePoint->setX(-9999.);
552 //     tmpTpcEntrancePoint->setY(-9999.);
553 //     tmpTpcEntrancePoint->setZ(-9999.);
554 //   } 
555     
556 //   if (IsNaN(tmpTpcExitPoint->x()) || 
557 //       IsNaN(tmpTpcExitPoint->y()) || 
558 //       IsNaN(tmpTpcExitPoint->z()) ) {
559 // //     cout << "tmpTpcExitPoint NAN set at (-9999,-9999,-9999)"<< endl; 
560 // //     cout << "tmpTpcExitPoint X= " <<tmpTpcExitPoint->x()<< endl;
561 // //     cout << "tmpTpcExitPoint Y= " <<tmpTpcExitPoint->y()<< endl;
562 // //     cout << "tmpTpcExitPoint Z= " <<tmpTpcExitPoint->z()<< endl;
563 //     tmpTpcExitPoint->setX(-9999.);
564 //     tmpTpcExitPoint->setY(-9999.);
565 //     tmpTpcExitPoint->setZ(-9999.);
566 //   }
567
568
569 // //   if (IsNaN(tmpTpcExitPoint->x())) *tmpTpcExitPoint = AliFemtoThreeVector(-9999.,-9999.,-9999); 
570 // //   if (IsNaN(tmpTpcEntrancetPoint->x())) *tmpTpcEntrancePoint = AliFemtoThreeVector(-9999.,-9999.,-9999); 
571 //   //  cout << "tmpTpcEntrancePoint"<<*tmpTpcEntrancePoint << endl;
572
573 //   // 03Oct00 - mal.  OK, let's try something a little more 
574 //   // along the lines of NA49 and E895 strategy.
575 //   // calculate the "nominal" position at N radii (say N=11) 
576 //   // within the TPC, and for a pair cut
577 //   // use the average separation of these N
578 //   int irad = 0;
579 //   candidates = hel.PathLength(50.0);
580 //   sideLength = (candidates.first > 0) ? candidates.first : candidates.second;
581 //   while (irad<11 && !IsNaN(sideLength)){ 
582 //     float radius = 50.0 + irad*15.0;
583 //     candidates = hel.PathLength(radius);
584 //     sideLength = (candidates.first > 0) ? candidates.first : candidates.second;
585 //     tmpPosSample[irad] = hel.At(sideLength);
586 //     if(IsNaN(tmpPosSample[irad].x()) ||
587 //        IsNaN(tmpPosSample[irad].y()) ||
588 //        IsNaN(tmpPosSample[irad].z()) 
589 //        ){
590 //       cout << "tmpPosSample for radius=" << radius << " NAN"<< endl; 
591 //       cout << "tmpPosSample=(" <<tmpPosSample[irad]<<")"<< endl;
592 //       tmpPosSample[irad] =  AliFemtoThreeVector(-9999.,-9999.,-9999);
593 //     }
594 //     irad++;
595 //     if (irad<11){
596 //       float radius = 50.0 + irad*15.0;
597 //       candidates = hel.PathLength(radius);
598 //       sideLength = (candidates.first > 0) ? candidates.first : candidates.second;
599 //     }
600 //    }
601 //    for (int i = irad; i<11; i++)
602 //      {
603 //        tmpPosSample[i] =  AliFemtoThreeVector(-9999.,-9999.,-9999);   
604 //      }
605
606 //   static float tRowRadius[45] = {60,64.8,69.6,74.4,79.2,84,88.8,93.6,98.8, 
607 //                               104,109.2,114.4,119.6,127.195,129.195,131.195,
608 //                               133.195,135.195,137.195,139.195,141.195,
609 //                               143.195,145.195,147.195,149.195,151.195,
610 //                               153.195,155.195,157.195,159.195,161.195,
611 //                               163.195,165.195,167.195,169.195,171.195,
612 //                               173.195,175.195,177.195,179.195,181.195,
613 //                               183.195,185.195,187.195,189.195};
614 //   int tRow,tSect,tOutOfBound;
615 //   double tLength,tPhi;
616 //   float tU;
617 //   AliFemtoThreeVector tPoint;
618 //   AliFmThreeVectorD tn(0,0,0);
619 //   AliFmThreeVectorD tr(0,0,0);
620 //   int ti =0;
621 //   // test to enter the loop
622 //   candidates =  hel.PathLength(tRowRadius[ti]);
623 //   tLength = (candidates.first > 0) ? candidates.first : candidates.second;
624 //   if (IsNaN(tLength)){
625 //     cout <<"tLength Init tmp NAN" << endl;
626 //     cout <<"padrow number= "<<ti << "not reached" << endl;
627 //     cout << "*** DO NOT ENTER THE LOOP***" << endl;
628 //     tmpSect[ti]=-1;//sector
629 //   }
630 //   // end test
631 //   while(ti<45 && !IsNaN(tLength)){
632 //     candidates =  hel.PathLength(tRowRadius[ti]);
633 //     tLength = (candidates.first > 0) ? candidates.first : candidates.second;
634 //     if (IsNaN(tLength)){
635 //       cout <<"tLength loop 1st NAN" << endl;
636 //       cout <<"padrow number=  " << ti << " not reached" << endl;
637 //       cout << "*** THIS IS AN ERROR SHOULDN'T  LOOP ***" << endl;
638 //       tmpSect[ti]=-1;//sector
639 //     }
640 //     tPoint = hel.At(tLength);
641 //     // Find which sector it is on
642 //     TpcLocalTransform(tPoint,tmpSect[ti],tRow,tU,tPhi);
643 //     if (IsNaN(tmpSect[ti])){
644 //       cout <<"***ERROR tmpSect"<< endl; 
645 //     }
646 //     if (IsNaN(tRow)){
647 //       cout <<"***ERROR tRow"<< endl;
648 //     }
649 //     if (IsNaN(tU)){
650 //       cout <<"***ERROR tU"<< endl;
651 //     }
652 //     if (IsNaN(tPhi)){
653 //       cout <<"***ERROR tPhi"<< endl;
654 //     }  
655 //     // calculate crossing plane
656 //     tn.setX(cos(tPhi));
657 //     tn.setY(sin(tPhi));       
658 //     tr.setX(tRowRadius[ti]*cos(tPhi));
659 //     tr.setY(tRowRadius[ti]*sin(tPhi));
660 //     // find crossing point
661 //     tLength = hel.PathLength(tr,tn); 
662 //     if (IsNaN(tLength)){
663 //       cout <<"tLength loop 2nd  NAN" << endl;
664 //       cout <<"padrow number=  " << ti << " not reached" << endl;
665 //       tmpSect[ti]=-2;//sector
666 //     }
667 //     tPoint = hel.At(tLength);
668 //     tmpZ[ti] = tPoint.z();
669 //     tOutOfBound = TpcLocalTransform(tPoint,tSect,tRow,tmpU[ti],tPhi);
670 //     if (IsNaN(tSect)){
671 //       cout <<"***ERROR tSect 2"<< endl; 
672 //     }
673 //     if (IsNaN(tRow)){
674 //       cout <<"***ERROR tRow 2"<< endl;
675 //     }
676 //     if (IsNaN(tmpU[ti])){
677 //       cout <<"***ERROR tmpU[ti] 2"<< endl;
678 //     }
679 //     if (IsNaN(tPhi)){
680 //       cout <<"***ERROR tPhi 2 "<< endl;
681 //     }  
682 //     if(tOutOfBound || (tmpSect[ti] == tSect && tRow!=(ti+1))){
683 //       tmpSect[ti]=-2;
684 //       //       cout << "missed once"<< endl;
685 //     }
686 //     else{
687 //       if(tmpSect[ti] != tSect){
688 //      // Try again on the other sector
689 //      tn.setX(cos(tPhi));
690 //      tn.setY(sin(tPhi));       
691 //      tr.setX(tRowRadius[ti]*cos(tPhi));
692 //      tr.setY(tRowRadius[ti]*sin(tPhi));
693 //      // find crossing point
694 //      tLength = hel.PathLength(tr,tn);
695 //      tPoint = hel.At(tLength);
696 //      if (IsNaN(tLength)){
697 //        cout <<"tLength loop 3rd NAN" << endl;
698 //        cout <<"padrow number=  "<< ti << " not reached" << endl;
699 //        tmpSect[ti]=-1;//sector
700 //      }
701 //      tmpZ[ti] = tPoint.z();
702 //      tmpSect[ti] = tSect;
703 //      tOutOfBound = TpcLocalTransform(tPoint,tSect,tRow,tmpU[ti],tPhi);
704 //      if (IsNaN(tSect)){
705 //        cout <<"***ERROR tSect 3"<< endl; 
706 //      }
707 //      if (IsNaN(tRow)){
708 //        cout <<"***ERROR tRow 3"<< endl;
709 //      }
710 //      if (IsNaN(tmpU[ti])){
711 //        cout <<"***ERROR tmpU[ti] 3"<< endl;
712 //      }
713 //      if (IsNaN(tPhi)){
714 //        cout <<"***ERROR tPhi 3 "<< endl;
715 //      }  
716 //      if(tOutOfBound || tSect!= tmpSect[ti] || tRow!=(ti+1)){
717 //        tmpSect[ti]=-1;
718 //      }
719 //       }
720 //     }
721 //     if (IsNaN(tmpSect[ti])){
722 //       cout << "*******************ERROR***************************" << endl;
723 //       cout <<"AliFemtoParticle--Fctn tmpSect=" << tmpSect[ti] << endl;
724 //       cout << "*******************ERROR***************************" << endl;
725 //     }
726 //     if (IsNaN(tmpU[ti])){
727 //       cout << "*******************ERROR***************************" << endl;
728 //       cout <<"AliFemtoParticle--Fctn tmpU=" << tmpU[ti] << endl;
729 //       cout << "*******************ERROR***************************" << endl;
730 //     }
731 //     if (IsNaN(tmpZ[ti])){
732 //       cout << "*******************ERROR***************************" << endl;
733 //       cout <<"AliFemtoParticle--Fctn tmpZ=" << tmpZ[ti] << endl;
734 //       cout << "*******************ERROR***************************" << endl;
735 //     }
736 //     // If padrow ti not reached all other beyond are not reached
737 //     // in this case set sector to -1
738 //     if (tmpSect[ti]==-1){
739 //       for (int tj=ti; tj<45;tj++){
740 //      tmpSect[tj] = -1;
741 //      ti=45;
742 //       }
743 //     }
744 //     ti++;
745 //     if (ti<45){
746 //       candidates =  hel.PathLength(tRowRadius[ti]);
747 //       tLength = (candidates.first > 0) ? candidates.first : candidates.second;}
748 //   }
749 // }
750 //_____________________
751 const AliFemtoThreeVector& AliFemtoParticle::TpcV0PosExitPoint() const{
752   return fTpcV0PosExitPoint;
753 }
754 //_____________________
755 const AliFemtoThreeVector& AliFemtoParticle::TpcV0PosEntrancePoint() const{
756   return fTpcV0PosEntrancePoint;
757 }
758 //______________________
759 const AliFemtoThreeVector& AliFemtoParticle::TpcV0NegExitPoint() const{
760   return fTpcV0NegExitPoint;
761 }
762 //_____________________
763 const AliFemtoThreeVector& AliFemtoParticle::TpcV0NegEntrancePoint() const{
764   return fTpcV0NegEntrancePoint;
765 }
766 //______________________