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