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