41b425f7d62ab462183368def7722f0fc102e2d2
[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 //______________________