647980887ea316c2ed770b6b94b2aa70a73024ee
[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     return *this;
360
361   if (aParticle.fTrack)
362     fTrack = new AliFemtoTrack(*aParticle.fTrack);
363   if (aParticle.fV0)
364     fV0    = new AliFemtoV0(*aParticle.fV0);
365   if (aParticle.fKink)
366     fKink  = new AliFemtoKink(*aParticle.fKink);
367   if (aParticle.fXi)
368     fXi    = new AliFemtoXi(*aParticle.fXi);
369   fFourMomentum = aParticle.fFourMomentum;
370   fHelix = aParticle.fHelix;
371
372 //   for (int iter=0; iter<11; iter++)
373 //     fNominalPosSample[iter] = aParticle.fNominalPosSample[iter];
374
375 //   if (fTpcV0NegPosSample) delete fTpcV0NegPosSample;
376 //   if (aParticle.fTpcV0NegPosSample) {
377 //     fTpcV0NegPosSample = (AliFemtoThreeVector *) malloc(sizeof(AliFemtoThreeVector) * 11);
378 //     for (int iter=0; iter<11; iter++)
379 //       fTpcV0NegPosSample[iter] = aParticle.fTpcV0NegPosSample[iter];
380 //   }
381
382 //   if (fV0NegZ) delete fV0NegZ;
383 //   if (aParticle.fV0NegZ) {
384 //     fV0NegZ = (float *) malloc(sizeof(float) * 45);
385 //     for (int iter=0; iter<11; iter++)
386 //       fV0NegZ[iter] = aParticle.fV0NegZ[iter];
387 //   }
388 //   if (fV0NegU) delete fV0NegU;
389 //   if (aParticle.fV0NegU) {
390 //     fV0NegU = (float *) malloc(sizeof(float) * 45);
391 //     for (int iter=0; iter<11; iter++)
392 //       fV0NegU[iter] = aParticle.fV0NegU[iter];
393 //   }
394 //   if (fV0NegSect) delete fV0NegSect;
395 //   if (aParticle.fV0NegSect) {
396 //     fV0NegSect = (int *) malloc(sizeof(int) * 45);
397 //     for (int iter=0; iter<11; iter++)
398 //       fV0NegSect[iter] = aParticle.fV0NegSect[iter];
399 //   }
400
401   fPrimaryVertex = aParticle.fPrimaryVertex;
402   fSecondaryVertex = aParticle.fSecondaryVertex;
403   CalculatePurity();
404   if(aParticle.fHiddenInfo){
405     fHiddenInfo= aParticle.fHiddenInfo->Clone();
406   }
407   
408 //   fNominalTpcEntrancePoint = aParticle.fNominalTpcEntrancePoint;
409 //   fNominalTpcExitPoint     = aParticle.fNominalTpcExitPoint;
410  
411   if (fHiddenInfo) delete fHiddenInfo;
412   if (aParticle.fHiddenInfo) 
413     fHiddenInfo = aParticle.HiddenInfo()->Clone();
414   
415   for (int iter=0; iter<6; iter++)
416     fPurity[iter] = aParticle.fPurity[iter];
417   
418   fHelixV0Pos = aParticle.fHelixV0Pos;
419   fTpcV0PosEntrancePoint = aParticle.fTpcV0PosEntrancePoint;
420   fTpcV0PosExitPoint     = aParticle.fTpcV0PosExitPoint;
421   fHelixV0Neg = aParticle.fHelixV0Neg;
422   fTpcV0NegEntrancePoint = aParticle.fTpcV0NegEntrancePoint;
423   fTpcV0NegExitPoint     = aParticle.fTpcV0NegExitPoint;
424
425   return *this;
426 }
427 // //_____________________
428 // const AliFemtoThreeVector& AliFemtoParticle::NominalTpcExitPoint() const{
429 //   // in future, may want to calculate this "on demand" only, sot this routine may get more sophisticated
430 //   // for now, we calculate Exit and Entrance points upon instantiation
431 //   return fNominalTpcExitPoint;
432 // }
433 // //_____________________
434 // const AliFemtoThreeVector& AliFemtoParticle::NominalTpcEntrancePoint() const{
435 //   // in future, may want to calculate this "on demand" only, sot this routine may get more sophisticated
436 //   // for now, we calculate Exit and Entrance points upon instantiation
437 //   return fNominalTpcEntrancePoint;
438 // }
439 //_____________________
440 void AliFemtoParticle::CalculatePurity(){
441   // Calculate additional parameterized purity
442
443   double tPt = fFourMomentum.Perp();
444   // pi -
445   fPurity[0] = fgPrimPimPar0*(1.-exp((tPt-fgPrimPimPar1)/fgPrimPimPar2));
446   fPurity[0] *= fTrack->PidProbPion();
447   // pi+
448   fPurity[1] = fgPrimPipPar0*(1.-exp((tPt-fgPrimPipPar1)/fgPrimPipPar2));
449   fPurity[1] *= fTrack->PidProbPion();
450   // K-
451   fPurity[2] = fTrack->PidProbKaon();
452   // K+
453   fPurity[3] = fTrack->PidProbKaon();
454   // pbar
455   fPurity[4] = fTrack->PidProbProton();
456   // p
457   fPurity[5] = fTrack->PidProbProton();
458 }
459
460 double AliFemtoParticle::GetPionPurity()
461 {
462   // Get full pion purity
463   if (fTrack->Charge()>0)
464     return fPurity[1];
465   else
466     return fPurity[0];
467 }
468 double AliFemtoParticle::GetKaonPurity()
469 {
470   // Get full kaon purity
471   if (fTrack->Charge()>0)
472     return fPurity[3];
473   else
474     return fPurity[2];
475 }
476 double AliFemtoParticle::GetProtonPurity()
477 {
478   // Get full proton purity
479   if (fTrack->Charge()>0)
480     return fPurity[5];
481   else
482     return fPurity[4];
483 }
484
485 // void AliFemtoParticle::CalculateTpcExitAndEntrancePoints(AliFmPhysicalHelixD* tHelix,
486 //                                                     AliFemtoThreeVector*  PrimVert,
487 //                                                     AliFemtoThreeVector*  SecVert,
488 //                                                     AliFemtoThreeVector* tmpTpcEntrancePoint,
489 //                                                     AliFemtoThreeVector* tmpTpcExitPoint,
490 //                                                     AliFemtoThreeVector* tmpPosSample,
491 //                                                     float* tmpZ,
492 //                                                     float* tmpU,
493 //                                                     int* tmpSect){
494 //   // this calculates the exit point of a secondary track, 
495 //   // either through the endcap or through the Outer Field Cage
496 //   // We assume the track to start at tHelix.origin-PrimaryVertex
497 //   // it also calculates the entrance point of the secondary track, 
498 //   // which is the point at which it crosses the
499 //   // inner field cage
500 //   //  static AliFemtoThreeVector ZeroVec(0.,0.,0.);
501 //   AliFemtoThreeVector tZeroVec(0.,0.,0.);
502 // //   tZeroVec.SetX(tHelix->origin().x()-PrimVert->x());
503 // //   tZeroVec.SetY(tHelix->origin().y()-PrimVert->y());
504 // //   tZeroVec.SetZ(tHelix->origin().z()-PrimVert->z());
505 //   tZeroVec.SetX(SecVert->x()-PrimVert->x());
506 //   tZeroVec.SetY(SecVert->y()-PrimVert->y());
507 //   tZeroVec.SetZ(SecVert->z()-PrimVert->z());
508 //   double dip, curv, phase;
509 //   int h;
510 //   curv = tHelix->Curvature();
511 //   dip  = tHelix->DipAngle();
512 //   phase= tHelix->Phase();
513 //   h    = tHelix->H();
514   
515 //   AliFmHelixD hel(curv,dip,phase,tZeroVec,h);
516
517 //   pairD candidates;
518 //   double sideLength;  // this is how much length to go to leave through sides of TPC
519 //   double endLength;  // this is how much length to go to leave through endcap of TPC
520 //   // figure out how far to go to leave through side...
521 //   candidates = hel.PathLength(200.0);  // bugfix MAL jul00 - 200cm NOT 2cm
522 //   sideLength = (candidates.first > 0) ? candidates.first : candidates.second;
523
524 //   static AliFemtoThreeVector tWestEnd(0.,0.,200.);  // bugfix MAL jul00 - 200cm NOT 2cm
525 //   static AliFemtoThreeVector tEastEnd(0.,0.,-200.); // bugfix MAL jul00 - 200cm NOT 2cm
526 //   static AliFemtoThreeVector tEndCapNormal(0.,0.,1.0);
527
528 //   endLength = hel.PathLength(tWestEnd,tEndCapNormal);
529 //   if (endLength < 0.0) endLength = hel.PathLength(tEastEnd,tEndCapNormal);
530
531 //   if (endLength < 0.0) cout << 
532 //                       "AliFemtoParticle::CalculateTpcExitAndEntrancePoints(): "
533 //                             << "Hey -- I cannot find an exit point out endcaps" << endl;
534 //   // OK, firstExitLength will be the shortest way out of the detector...
535 //   double firstExitLength = (endLength < sideLength) ? endLength : sideLength;
536 //   // now then, let's return the POSITION at which particle leaves TPC...
537 //   *tmpTpcExitPoint = hel.At(firstExitLength);
538 //   // Finally, calculate the position at which the track crosses the inner field cage
539 //   candidates = hel.PathLength(50.0);  // bugfix MAL jul00 - 200cm NOT 2cm
540
541 //   sideLength = (candidates.first > 0) ? candidates.first : candidates.second;
542 //   //  cout << "sideLength 2 ="<<sideLength << endl;
543 //   *tmpTpcEntrancePoint = hel.At(sideLength);
544 //   // This is the secure way !  
545 //   if (IsNaN(tmpTpcEntrancePoint->x()) || 
546 //       IsNaN(tmpTpcEntrancePoint->y()) || 
547 //       IsNaN(tmpTpcEntrancePoint->z()) ){ 
548 //     cout << "tmpTpcEntrancePoint NAN"<< endl; 
549 //     cout << "tmpNominalTpcEntrancePoint = " <<tmpTpcEntrancePoint<< endl;
550 //     tmpTpcEntrancePoint->SetX(-9999.);
551 //     tmpTpcEntrancePoint->SetY(-9999.);
552 //     tmpTpcEntrancePoint->SetZ(-9999.);
553 //   } 
554     
555 //   if (IsNaN(tmpTpcExitPoint->x()) || 
556 //       IsNaN(tmpTpcExitPoint->y()) || 
557 //       IsNaN(tmpTpcExitPoint->z()) ) {
558 // //     cout << "tmpTpcExitPoint NAN Set at (-9999,-9999,-9999)"<< endl; 
559 // //     cout << "tmpTpcExitPoint X= " <<tmpTpcExitPoint->x()<< endl;
560 // //     cout << "tmpTpcExitPoint Y= " <<tmpTpcExitPoint->y()<< endl;
561 // //     cout << "tmpTpcExitPoint Z= " <<tmpTpcExitPoint->z()<< endl;
562 //     tmpTpcExitPoint->SetX(-9999.);
563 //     tmpTpcExitPoint->SetY(-9999.);
564 //     tmpTpcExitPoint->SetZ(-9999.);
565 //   }
566
567
568 // //   if (IsNaN(tmpTpcExitPoint->x())) *tmpTpcExitPoint = AliFemtoThreeVector(-9999.,-9999.,-9999); 
569 // //   if (IsNaN(tmpTpcEntrancetPoint->x())) *tmpTpcEntrancePoint = AliFemtoThreeVector(-9999.,-9999.,-9999); 
570 //   //  cout << "tmpTpcEntrancePoint"<<*tmpTpcEntrancePoint << endl;
571
572 //   // 03Oct00 - mal.  OK, let's try something a little more 
573 //   // along the lines of NA49 and E895 strategy.
574 //   // calculate the "nominal" position at N radii (say N=11) 
575 //   // within the TPC, and for a pair cut
576 //   // use the average separation of these N
577 //   int irad = 0;
578 //   candidates = hel.PathLength(50.0);
579 //   sideLength = (candidates.first > 0) ? candidates.first : candidates.second;
580 //   while (irad<11 && !IsNaN(sideLength)){ 
581 //     float radius = 50.0 + irad*15.0;
582 //     candidates = hel.PathLength(radius);
583 //     sideLength = (candidates.first > 0) ? candidates.first : candidates.second;
584 //     tmpPosSample[irad] = hel.At(sideLength);
585 //     if(IsNaN(tmpPosSample[irad].x()) ||
586 //        IsNaN(tmpPosSample[irad].y()) ||
587 //        IsNaN(tmpPosSample[irad].z()) 
588 //        ){
589 //       cout << "tmpPosSample for radius=" << radius << " NAN"<< endl; 
590 //       cout << "tmpPosSample=(" <<tmpPosSample[irad]<<")"<< endl;
591 //       tmpPosSample[irad] =  AliFemtoThreeVector(-9999.,-9999.,-9999);
592 //     }
593 //     irad++;
594 //     if (irad<11){
595 //       float radius = 50.0 + irad*15.0;
596 //       candidates = hel.PathLength(radius);
597 //       sideLength = (candidates.first > 0) ? candidates.first : candidates.second;
598 //     }
599 //    }
600 //    for (int i = irad; i<11; i++)
601 //      {
602 //        tmpPosSample[i] =  AliFemtoThreeVector(-9999.,-9999.,-9999);   
603 //      }
604
605 //   static float tRowRadius[45] = {60,64.8,69.6,74.4,79.2,84,88.8,93.6,98.8, 
606 //                               104,109.2,114.4,119.6,127.195,129.195,131.195,
607 //                               133.195,135.195,137.195,139.195,141.195,
608 //                               143.195,145.195,147.195,149.195,151.195,
609 //                               153.195,155.195,157.195,159.195,161.195,
610 //                               163.195,165.195,167.195,169.195,171.195,
611 //                               173.195,175.195,177.195,179.195,181.195,
612 //                               183.195,185.195,187.195,189.195};
613 //   int tRow,tSect,tOutOfBound;
614 //   double tLength,tPhi;
615 //   float tU;
616 //   AliFemtoThreeVector tPoint;
617 //   AliFmThreeVectorD tn(0,0,0);
618 //   AliFmThreeVectorD tr(0,0,0);
619 //   int ti =0;
620 //   // test to enter the loop
621 //   candidates =  hel.PathLength(tRowRadius[ti]);
622 //   tLength = (candidates.first > 0) ? candidates.first : candidates.second;
623 //   if (IsNaN(tLength)){
624 //     cout <<"tLength Init tmp NAN" << endl;
625 //     cout <<"padrow number= "<<ti << "not reached" << endl;
626 //     cout << "*** DO NOT ENTER THE LOOP***" << endl;
627 //     tmpSect[ti]=-1;//sector
628 //   }
629 //   // end test
630 //   while(ti<45 && !IsNaN(tLength)){
631 //     candidates =  hel.PathLength(tRowRadius[ti]);
632 //     tLength = (candidates.first > 0) ? candidates.first : candidates.second;
633 //     if (IsNaN(tLength)){
634 //       cout <<"tLength loop 1st NAN" << endl;
635 //       cout <<"padrow number=  " << ti << " not reached" << endl;
636 //       cout << "*** THIS IS AN ERROR SHOULDN'T  LOOP ***" << endl;
637 //       tmpSect[ti]=-1;//sector
638 //     }
639 //     tPoint = hel.At(tLength);
640 //     // Find which sector it is on
641 //     TpcLocalTransform(tPoint,tmpSect[ti],tRow,tU,tPhi);
642 //     if (IsNaN(tmpSect[ti])){
643 //       cout <<"***ERROR tmpSect"<< endl; 
644 //     }
645 //     if (IsNaN(tRow)){
646 //       cout <<"***ERROR tRow"<< endl;
647 //     }
648 //     if (IsNaN(tU)){
649 //       cout <<"***ERROR tU"<< endl;
650 //     }
651 //     if (IsNaN(tPhi)){
652 //       cout <<"***ERROR tPhi"<< endl;
653 //     }  
654 //     // calculate crossing plane
655 //     tn.SetX(cos(tPhi));
656 //     tn.SetY(sin(tPhi));       
657 //     tr.SetX(tRowRadius[ti]*cos(tPhi));
658 //     tr.SetY(tRowRadius[ti]*sin(tPhi));
659 //     // find crossing point
660 //     tLength = hel.PathLength(tr,tn); 
661 //     if (IsNaN(tLength)){
662 //       cout <<"tLength loop 2nd  NAN" << endl;
663 //       cout <<"padrow number=  " << ti << " not reached" << endl;
664 //       tmpSect[ti]=-2;//sector
665 //     }
666 //     tPoint = hel.At(tLength);
667 //     tmpZ[ti] = tPoint.z();
668 //     tOutOfBound = TpcLocalTransform(tPoint,tSect,tRow,tmpU[ti],tPhi);
669 //     if (IsNaN(tSect)){
670 //       cout <<"***ERROR tSect 2"<< endl; 
671 //     }
672 //     if (IsNaN(tRow)){
673 //       cout <<"***ERROR tRow 2"<< endl;
674 //     }
675 //     if (IsNaN(tmpU[ti])){
676 //       cout <<"***ERROR tmpU[ti] 2"<< endl;
677 //     }
678 //     if (IsNaN(tPhi)){
679 //       cout <<"***ERROR tPhi 2 "<< endl;
680 //     }  
681 //     if(tOutOfBound || (tmpSect[ti] == tSect && tRow!=(ti+1))){
682 //       tmpSect[ti]=-2;
683 //       //       cout << "missed once"<< endl;
684 //     }
685 //     else{
686 //       if(tmpSect[ti] != tSect){
687 //      // Try again on the other sector
688 //      tn.SetX(cos(tPhi));
689 //      tn.SetY(sin(tPhi));       
690 //      tr.SetX(tRowRadius[ti]*cos(tPhi));
691 //      tr.SetY(tRowRadius[ti]*sin(tPhi));
692 //      // find crossing point
693 //      tLength = hel.PathLength(tr,tn);
694 //      tPoint = hel.At(tLength);
695 //      if (IsNaN(tLength)){
696 //        cout <<"tLength loop 3rd NAN" << endl;
697 //        cout <<"padrow number=  "<< ti << " not reached" << endl;
698 //        tmpSect[ti]=-1;//sector
699 //      }
700 //      tmpZ[ti] = tPoint.z();
701 //      tmpSect[ti] = tSect;
702 //      tOutOfBound = TpcLocalTransform(tPoint,tSect,tRow,tmpU[ti],tPhi);
703 //      if (IsNaN(tSect)){
704 //        cout <<"***ERROR tSect 3"<< endl; 
705 //      }
706 //      if (IsNaN(tRow)){
707 //        cout <<"***ERROR tRow 3"<< endl;
708 //      }
709 //      if (IsNaN(tmpU[ti])){
710 //        cout <<"***ERROR tmpU[ti] 3"<< endl;
711 //      }
712 //      if (IsNaN(tPhi)){
713 //        cout <<"***ERROR tPhi 3 "<< endl;
714 //      }  
715 //      if(tOutOfBound || tSect!= tmpSect[ti] || tRow!=(ti+1)){
716 //        tmpSect[ti]=-1;
717 //      }
718 //       }
719 //     }
720 //     if (IsNaN(tmpSect[ti])){
721 //       cout << "*******************ERROR***************************" << endl;
722 //       cout <<"AliFemtoParticle--Fctn tmpSect=" << tmpSect[ti] << endl;
723 //       cout << "*******************ERROR***************************" << endl;
724 //     }
725 //     if (IsNaN(tmpU[ti])){
726 //       cout << "*******************ERROR***************************" << endl;
727 //       cout <<"AliFemtoParticle--Fctn tmpU=" << tmpU[ti] << endl;
728 //       cout << "*******************ERROR***************************" << endl;
729 //     }
730 //     if (IsNaN(tmpZ[ti])){
731 //       cout << "*******************ERROR***************************" << endl;
732 //       cout <<"AliFemtoParticle--Fctn tmpZ=" << tmpZ[ti] << endl;
733 //       cout << "*******************ERROR***************************" << endl;
734 //     }
735 //     // If padrow ti not reached all other beyond are not reached
736 //     // in this case Set sector to -1
737 //     if (tmpSect[ti]==-1){
738 //       for (int tj=ti; tj<45;tj++){
739 //      tmpSect[tj] = -1;
740 //      ti=45;
741 //       }
742 //     }
743 //     ti++;
744 //     if (ti<45){
745 //       candidates =  hel.PathLength(tRowRadius[ti]);
746 //       tLength = (candidates.first > 0) ? candidates.first : candidates.second;}
747 //   }
748 // }
749 //_____________________
750 const AliFemtoThreeVector& AliFemtoParticle::TpcV0PosExitPoint() const{
751   return fTpcV0PosExitPoint;
752 }
753 //_____________________
754 const AliFemtoThreeVector& AliFemtoParticle::TpcV0PosEntrancePoint() const{
755   return fTpcV0PosEntrancePoint;
756 }
757 //______________________
758 const AliFemtoThreeVector& AliFemtoParticle::TpcV0NegExitPoint() const{
759   return fTpcV0NegExitPoint;
760 }
761 //_____________________
762 const AliFemtoThreeVector& AliFemtoParticle::TpcV0NegEntrancePoint() const{
763   return fTpcV0NegEntrancePoint;
764 }
765 //______________________