1 /***************************************************************************
5 * Author: Mike Lisa, Ohio State, lisa@mps.ohio-state.edu
6 ***************************************************************************
8 * Description: part of STAR HBT Framework: AliFemtoMaker package
9 * Particle objects are part of the PicoEvent, which is what is
10 * stored in the EventMixingBuffers
11 * A Track object gets converted to a Particle object if it
12 * passes the ParticleCut of an Analysis
14 ***************************************************************************
17 * Revision 1.1.1.1 2007/04/25 15:38:41 panos
18 * Importing the HBT code dir
20 * Revision 1.1.1.1 2007/03/07 10:14:49 mchojnacki
21 * First version on CVS
23 * Revision 1.22 2003/09/02 17:58:32 perev
24 * gcc 3.2 updates + WarnOff
26 * Revision 1.21 2003/05/07 15:30:43 magestro
27 * Fixed bug related to finding merged hits (commit for Fabrice)
29 * Revision 1.20 2003/01/14 09:41:26 renault
30 * changes on average separation calculation, hit shared finder and memory optimisation
31 * for Z,U and Sectors variables.
33 * Revision 1.19 2002/12/12 17:01:49 kisiel
34 * Hidden Information handling and purity calculation
36 * Revision 1.18 2002/11/19 23:36:00 renault
37 * Enable calculation of exit/entrance separation for V0 daughters
39 * Revision 1.17 2001/12/14 23:11:30 fretiere
40 * Add class HitMergingCut. Add class fabricesPairCut = HitMerginCut + pair purity cuts. Add TpcLocalTransform function which convert to local tpc coord (not pretty). Modify AliFemtoTrack, AliFemtoParticle, AliFemtoHiddenInfo, AliFemtoPair to handle the hit information and cope with my code
42 * Revision 1.16 2001/05/25 23:23:59 lisa
43 * Added in AliFemtoKink stuff
45 * Revision 1.15 2001/04/03 21:04:36 kisiel
48 * Changes needed to make the Theoretical code
49 * work. The main code is the ThCorrFctn directory.
50 * The most visible change is the addition of the
51 * HiddenInfo to AliFemtoPair.
53 * Revision 1.14 2000/10/05 23:09:05 lisa
54 * Added kT-dependent radii to mixed-event simulator AND implemented AverageSeparation Cut and CorrFctn
56 * Revision 1.13 2000/08/31 22:31:31 laue
57 * AliFemtoAnalysis: output changed (a little bit less)
58 * AliFemtoEvent: new version, members for reference mult added
59 * AliFemtoIOBinary: new IO for new AliFemtoEvent version
60 * AliFemtoTypes: TTree typedef to AliFemtoTTree added
61 * AliFemtoVertexAnalysis: overflow and underflow added
63 * Revision 1.12 2000/07/23 13:52:56 laue
64 * NominalExitPoint set to (-9999.,-9999.-9999.) if helix.at()
65 * returns nan (not a number).
67 * Revision 1.11 2000/07/19 17:18:48 laue
68 * Calculation of Entrance and Exit point added in AliFemtoParticle constructor
70 * Revision 1.10 2000/07/17 20:03:17 lisa
71 * Implemented tools for addressing and assessing trackmerging
73 * Revision 1.9 2000/07/16 21:38:23 laue
74 * AliFemtoCoulomb.cxx AliFemtoSectoredAnalysis.cxx : updated for standalone version
75 * AliFemtoV0.cc AliFemtoV0.h : some cast to prevent compiling warnings
76 * AliFemtoParticle.cc AliFemtoParticle.h : pointers fTrack,fV0 initialized to 0
77 * AliFemtoIOBinary.cc : some printouts in #ifdef STHBTDEBUG
78 * AliFemtoEvent.cc : B-Field set to 0.25Tesla, we have to think about a better
81 * Revision 1.8 2000/05/03 17:44:43 laue
82 * AliFemtoEvent, AliFemtoTrack & AliFemtoV0 declared friend to AliFemtoIOBinary
83 * AliFemtoParticle updated for V0 pos,neg track Id
85 * Revision 1.7 2000/04/03 16:21:51 laue
86 * some include files changed
87 * Multi track cut added
89 * Revision 1.6 1999/12/11 15:58:29 lisa
90 * Add vertex decay position datum and accessor to AliFemtoParticle to allow pairwise cuts on seperation of V0s
92 * Revision 1.5 1999/09/17 22:38:02 lisa
93 * first full integration of V0s into AliFemto framework
95 * Revision 1.4 1999/09/01 19:04:53 lisa
96 * update Particle class AND add parity cf and Randys Coulomb correction
98 * Revision 1.3 1999/07/06 22:33:23 lisa
99 * Adjusted all to work in pro and new - dev itself is broken
101 * Revision 1.2 1999/06/29 17:50:27 fisyak
102 * formal changes to account new StEvent, does not complie yet
104 * Revision 1.1.1.1 1999/06/29 16:02:57 lisa
105 * Installation of AliFemtoMaker
107 **************************************************************************/
109 #include "Infrastructure/AliFemtoParticle.h"
110 //#include "math_constants.h"
119 using namespace TMath;
121 double AliFemtoParticle::fPrimPimPar0= 9.05632e-01;
122 double AliFemtoParticle::fPrimPimPar1= -2.26737e-01;
123 double AliFemtoParticle::fPrimPimPar2= -1.03922e-01;
124 double AliFemtoParticle::fPrimPipPar0= 9.09616e-01;
125 double AliFemtoParticle::fPrimPipPar1= -9.00511e-02;
126 double AliFemtoParticle::fPrimPipPar2= -6.02940e-02;
127 double AliFemtoParticle::fPrimPmPar0= 0.;
128 double AliFemtoParticle::fPrimPmPar1= 0.;
129 double AliFemtoParticle::fPrimPmPar2= 0.;
130 double AliFemtoParticle::fPrimPpPar0= 0.;
131 double AliFemtoParticle::fPrimPpPar1= 0.;
132 double AliFemtoParticle::fPrimPpPar2= 0.;
134 int TpcLocalTransform(AliFmThreeVectorD& xgl,
141 //_____________________
142 AliFemtoParticle::AliFemtoParticle() :
143 fTpcV0NegPosSample(0),
147 fTrack(0), fV0(0), fKink(0), fXi(0),
150 fNominalTpcExitPoint(0),
151 fNominalTpcEntrancePoint(0),
156 fTpcV0PosEntrancePoint(0),
157 fTpcV0PosExitPoint(0),
159 fTpcV0NegEntrancePoint(0),
160 fTpcV0NegExitPoint(0)
162 /* no-op for default */
163 // cout << "Created particle " << this << endl;
165 //_____________________
166 AliFemtoParticle::AliFemtoParticle(const AliFemtoParticle& aParticle):
167 fTpcV0NegPosSample(0),
171 fTrack(0), fV0(0), fKink(0), fXi(0),
174 fNominalTpcExitPoint(0),
175 fNominalTpcEntrancePoint(0),
180 fTpcV0PosEntrancePoint(0),
181 fTpcV0PosExitPoint(0),
183 fTpcV0NegEntrancePoint(0),
184 fTpcV0NegExitPoint(0)
186 if (aParticle.fTrack)
187 fTrack = new AliFemtoTrack(*aParticle.fTrack);
189 fV0 = new AliFemtoV0(*aParticle.fV0);
191 fKink = new AliFemtoKink(*aParticle.fKink);
193 fXi = new AliFemtoXi(*aParticle.fXi);
194 fFourMomentum = aParticle.fFourMomentum;
195 fHelix = aParticle.fHelix;
197 for (int iter=0; iter<11; iter++)
198 fNominalPosSample[iter] = aParticle.fNominalPosSample[iter];
200 if (aParticle.fTpcV0NegPosSample) {
201 fTpcV0NegPosSample = (AliFemtoThreeVector *) malloc(sizeof(AliFemtoThreeVector) * 11);
202 for (int iter=0; iter<11; iter++)
203 fTpcV0NegPosSample[iter] = aParticle.fTpcV0NegPosSample[iter];
206 if (aParticle.fV0NegZ) {
207 fV0NegZ = (float *) malloc(sizeof(float) * 45);
208 for (int iter=0; iter<11; iter++)
209 fV0NegZ[iter] = aParticle.fV0NegZ[iter];
211 if (aParticle.fV0NegU) {
212 fV0NegU = (float *) malloc(sizeof(float) * 45);
213 for (int iter=0; iter<11; iter++)
214 fV0NegU[iter] = aParticle.fV0NegU[iter];
216 if (aParticle.fV0NegSect) {
217 fV0NegSect = (int *) malloc(sizeof(int) * 45);
218 for (int iter=0; iter<11; iter++)
219 fV0NegSect[iter] = aParticle.fV0NegSect[iter];
222 fPrimaryVertex = aParticle.fPrimaryVertex;
223 fSecondaryVertex = aParticle.fSecondaryVertex;
225 if(aParticle.fHiddenInfo){
226 fHiddenInfo= aParticle.HiddenInfo()->clone();
229 fNominalTpcEntrancePoint = aParticle.fNominalTpcEntrancePoint;
230 fNominalTpcExitPoint = aParticle.fNominalTpcExitPoint;
232 for (int iter=0; iter<6; iter++)
233 fPurity[iter] = aParticle.fPurity[iter];
235 fHelixV0Pos = aParticle.fHelixV0Pos;
236 fTpcV0PosEntrancePoint = aParticle.fTpcV0PosEntrancePoint;
237 fTpcV0PosExitPoint = aParticle.fTpcV0PosExitPoint;
238 fHelixV0Neg = aParticle.fHelixV0Neg;
239 fTpcV0NegEntrancePoint = aParticle.fTpcV0NegEntrancePoint;
240 fTpcV0NegExitPoint = aParticle.fTpcV0NegExitPoint;
242 //_____________________
243 AliFemtoParticle::~AliFemtoParticle(){
244 // cout << "Issuing delete for AliFemtoParticle." << endl;
246 if (fTrack) delete fTrack;
248 delete[] fTpcV0NegPosSample;
254 if (fKink) delete fKink;
255 // cout << "Trying to delete HiddenInfo: " << fHiddenInfo << endl;
258 // cout << "Deleting HiddenInfo." << endl;
261 // cout << "Deleted particle " << this << endl;
263 //_____________________
264 AliFemtoParticle::AliFemtoParticle(const AliFemtoTrack* const hbtTrack,const double& mass) :
265 fTpcV0NegPosSample(0),
269 fTrack(0), fV0(0), fKink(0), fXi(0),
272 fNominalTpcExitPoint(0),
273 fNominalTpcEntrancePoint(0),
278 fTpcV0PosEntrancePoint(0),
279 fTpcV0PosExitPoint(0),
281 fTpcV0NegEntrancePoint(0),
282 fTpcV0NegExitPoint(0)
286 // I know there is a better way to do this...
287 fTrack = new AliFemtoTrack(*hbtTrack);
288 AliFemtoThreeVector temp = hbtTrack->P();
289 fFourMomentum.setVect(temp);
290 double ener = ::sqrt(temp.mag2()+mass*mass);
291 fFourMomentum.setE(ener);
292 // fMap[0] = hbtTrack->TopologyMap(0);
293 // fMap[1] = hbtTrack->TopologyMap(1);
294 // fNhits = hbtTrack->NHits();
295 fHelix = hbtTrack->Helix();
296 //CalculateNominalTpcExitAndEntrancePoints();
299 fPrimaryVertex.setX(0.);
300 fPrimaryVertex.setY(0.);
301 fPrimaryVertex.setZ(0.);
302 fSecondaryVertex.setX(0.);
303 fSecondaryVertex.setY(0.);
304 fSecondaryVertex.setZ(0.);
305 /* TO JA ODZNACZYLEM NIE WIEM DLACZEGO
306 CalculateTpcExitAndEntrancePoints(&fHelix,&fPrimaryVertex,
308 &fNominalTpcEntrancePoint,
309 &fNominalTpcExitPoint,
310 &mNominalPosSample[0],
318 if(hbtTrack->ValidHiddenInfo()){
319 fHiddenInfo= hbtTrack->getHiddenInfo()->getParticleHiddenInfo()->clone();
322 // cout << "Created particle " << this << endl;
325 //_____________________
326 AliFemtoParticle::AliFemtoParticle(const AliFemtoV0* const hbtV0,const double& mass) :
327 fTpcV0NegPosSample(0),
331 fTrack(0), fV0(0), fKink(0), fXi(0),
334 fNominalTpcExitPoint(0),
335 fNominalTpcEntrancePoint(0),
340 fTpcV0PosEntrancePoint(0),
341 fTpcV0PosExitPoint(0),
343 fTpcV0NegEntrancePoint(0),
344 fTpcV0NegExitPoint(0)
346 fV0 = new AliFemtoV0(*hbtV0);
349 // I know there is a better way to do this...
350 AliFemtoThreeVector temp = hbtV0->momV0();
351 fFourMomentum.setVect(temp);
352 double ener = ::sqrt(temp.mag2()+mass*mass);
353 fFourMomentum.setE(ener);
354 // Calculating TpcEntrancePoint for Positive V0 daugther
355 fPrimaryVertex = hbtV0->primaryVertex();
356 fSecondaryVertex = hbtV0->decayVertexV0();
357 fHelixV0Pos = hbtV0->HelixPos();
359 fTpcV0NegPosSample = new AliFemtoThreeVector[45];//for V0Neg
360 fV0NegZ = new float[45];//for V0Neg
361 fV0NegU = new float[45];//for V0Neg
362 fV0NegSect = new int[45];//for V0Neg
363 CalculateTpcExitAndEntrancePoints(&fHelixV0Pos,&fPrimaryVertex,
365 &fTpcV0PosEntrancePoint,
367 &fNominalPosSample[0],
370 fHelixV0Neg = hbtV0->HelixNeg();
372 CalculateTpcExitAndEntrancePoints(&fHelixV0Neg,
375 &fTpcV0NegEntrancePoint,
377 &fTpcV0NegPosSample[0],
379 &fV0NegU[0],&fV0NegSect[0]);
383 if(hbtV0->ValidHiddenInfo()){
384 fHiddenInfo= hbtV0->getHiddenInfo()->clone();
388 //_____________________
389 AliFemtoParticle::AliFemtoParticle(const AliFemtoKink* const hbtKink,const double& mass) :
390 fTpcV0NegPosSample(0),
394 fTrack(0), fV0(0), fKink(0), fXi(0),
397 fNominalTpcExitPoint(0),
398 fNominalTpcEntrancePoint(0),
403 fTpcV0PosEntrancePoint(0),
404 fTpcV0PosExitPoint(0),
406 fTpcV0NegEntrancePoint(0),
407 fTpcV0NegExitPoint(0)
409 fKink = new AliFemtoKink(*hbtKink);
412 // I know there is a better way to do this...
413 AliFemtoThreeVector temp = hbtKink->Parent().P();
414 fFourMomentum.setVect(temp);
415 double ener = ::sqrt(temp.mag2()+mass*mass);
416 fFourMomentum.setE(ener);
419 //_____________________
420 AliFemtoParticle::AliFemtoParticle(const AliFemtoXi* const hbtXi, const double& mass) :
421 fTpcV0NegPosSample(0),
425 fTrack(0), fV0(0), fKink(0), fXi(0),
428 fNominalTpcExitPoint(0),
429 fNominalTpcEntrancePoint(0),
434 fTpcV0PosEntrancePoint(0),
435 fTpcV0PosExitPoint(0),
437 fTpcV0NegEntrancePoint(0),
438 fTpcV0NegExitPoint(0)
440 fXi = new AliFemtoXi(*hbtXi);
443 AliFemtoThreeVector temp;// = hbtXi->mMofXi;
444 fFourMomentum.setVect(temp);
445 double ener = ::sqrt(temp.mag2()+mass*mass);
446 fFourMomentum.setE(ener);
449 //_____________________
450 AliFemtoParticle& AliFemtoParticle::operator=(const AliFemtoParticle& aParticle)
452 if (this == &aParticle)
455 if (aParticle.fTrack)
456 fTrack = new AliFemtoTrack(*aParticle.fTrack);
458 fV0 = new AliFemtoV0(*aParticle.fV0);
460 fKink = new AliFemtoKink(*aParticle.fKink);
462 fXi = new AliFemtoXi(*aParticle.fXi);
463 fFourMomentum = aParticle.fFourMomentum;
464 fHelix = aParticle.fHelix;
466 for (int iter=0; iter<11; iter++)
467 fNominalPosSample[iter] = aParticle.fNominalPosSample[iter];
469 if (fTpcV0NegPosSample) delete fTpcV0NegPosSample;
470 if (aParticle.fTpcV0NegPosSample) {
471 fTpcV0NegPosSample = (AliFemtoThreeVector *) malloc(sizeof(AliFemtoThreeVector) * 11);
472 for (int iter=0; iter<11; iter++)
473 fTpcV0NegPosSample[iter] = aParticle.fTpcV0NegPosSample[iter];
476 if (fV0NegZ) delete fV0NegZ;
477 if (aParticle.fV0NegZ) {
478 fV0NegZ = (float *) malloc(sizeof(float) * 45);
479 for (int iter=0; iter<11; iter++)
480 fV0NegZ[iter] = aParticle.fV0NegZ[iter];
482 if (fV0NegU) delete fV0NegU;
483 if (aParticle.fV0NegU) {
484 fV0NegU = (float *) malloc(sizeof(float) * 45);
485 for (int iter=0; iter<11; iter++)
486 fV0NegU[iter] = aParticle.fV0NegU[iter];
488 if (fV0NegSect) delete fV0NegSect;
489 if (aParticle.fV0NegSect) {
490 fV0NegSect = (int *) malloc(sizeof(int) * 45);
491 for (int iter=0; iter<11; iter++)
492 fV0NegSect[iter] = aParticle.fV0NegSect[iter];
495 fPrimaryVertex = aParticle.fPrimaryVertex;
496 fSecondaryVertex = aParticle.fSecondaryVertex;
498 if(aParticle.fHiddenInfo){
499 fHiddenInfo= aParticle.fHiddenInfo->clone();
502 fNominalTpcEntrancePoint = aParticle.fNominalTpcEntrancePoint;
503 fNominalTpcExitPoint = aParticle.fNominalTpcExitPoint;
505 if (fHiddenInfo) delete fHiddenInfo;
506 if (aParticle.fHiddenInfo)
507 fHiddenInfo = aParticle.HiddenInfo()->clone();
509 for (int iter=0; iter<6; iter++)
510 fPurity[iter] = aParticle.fPurity[iter];
512 fHelixV0Pos = aParticle.fHelixV0Pos;
513 fTpcV0PosEntrancePoint = aParticle.fTpcV0PosEntrancePoint;
514 fTpcV0PosExitPoint = aParticle.fTpcV0PosExitPoint;
515 fHelixV0Neg = aParticle.fHelixV0Neg;
516 fTpcV0NegEntrancePoint = aParticle.fTpcV0NegEntrancePoint;
517 fTpcV0NegExitPoint = aParticle.fTpcV0NegExitPoint;
521 //_____________________
522 const AliFemtoThreeVector& AliFemtoParticle::NominalTpcExitPoint() const{
523 // in future, may want to calculate this "on demand" only, sot this routine may get more sophisticated
524 // for now, we calculate Exit and Entrance points upon instantiation
525 return fNominalTpcExitPoint;
527 //_____________________
528 const AliFemtoThreeVector& AliFemtoParticle::NominalTpcEntrancePoint() const{
529 // in future, may want to calculate this "on demand" only, sot this routine may get more sophisticated
530 // for now, we calculate Exit and Entrance points upon instantiation
531 return fNominalTpcEntrancePoint;
533 //_____________________
534 void AliFemtoParticle::CalculatePurity(){
535 double tPt = fFourMomentum.perp();
537 fPurity[0] = fPrimPimPar0*(1.-exp((tPt-fPrimPimPar1)/fPrimPimPar2));
538 fPurity[0] *= fTrack->PidProbPion();
540 fPurity[1] = fPrimPipPar0*(1.-exp((tPt-fPrimPipPar1)/fPrimPipPar2));
541 fPurity[1] *= fTrack->PidProbPion();
543 fPurity[2] = fTrack->PidProbKaon();
545 fPurity[3] = fTrack->PidProbKaon();
547 fPurity[4] = fTrack->PidProbProton();
549 fPurity[5] = fTrack->PidProbProton();
552 double AliFemtoParticle::GetPionPurity()
554 if (fTrack->Charge()>0)
559 double AliFemtoParticle::GetKaonPurity()
561 if (fTrack->Charge()>0)
566 double AliFemtoParticle::GetProtonPurity()
568 if (fTrack->Charge()>0)
574 void AliFemtoParticle::CalculateTpcExitAndEntrancePoints(AliFmPhysicalHelixD* tHelix,
575 AliFemtoThreeVector* PrimVert,
576 AliFemtoThreeVector* SecVert,
577 AliFemtoThreeVector* tmpTpcEntrancePoint,
578 AliFemtoThreeVector* tmpTpcExitPoint,
579 AliFemtoThreeVector* tmpPosSample,
583 // this calculates the exit point of a secondary track,
584 // either through the endcap or through the Outer Field Cage
585 // We assume the track to start at tHelix.origin-PrimaryVertex
586 // it also calculates the entrance point of the secondary track,
587 // which is the point at which it crosses the
589 // static AliFemtoThreeVector ZeroVec(0.,0.,0.);
590 AliFemtoThreeVector ZeroVec(0.,0.,0.);
591 // ZeroVec.setX(tHelix->origin().x()-PrimVert->x());
592 // ZeroVec.setY(tHelix->origin().y()-PrimVert->y());
593 // ZeroVec.setZ(tHelix->origin().z()-PrimVert->z());
594 ZeroVec.setX(SecVert->x()-PrimVert->x());
595 ZeroVec.setY(SecVert->y()-PrimVert->y());
596 ZeroVec.setZ(SecVert->z()-PrimVert->z());
597 double dip, curv, phase;
599 curv = tHelix->curvature();
600 dip = tHelix->dipAngle();
601 phase= tHelix->phase();
604 AliFmHelixD hel(curv,dip,phase,ZeroVec,h);
607 double sideLength; // this is how much length to go to leave through sides of TPC
608 double endLength; // this is how much length to go to leave through endcap of TPC
609 // figure out how far to go to leave through side...
610 candidates = hel.pathLength(200.0); // bugfix MAL jul00 - 200cm NOT 2cm
611 sideLength = (candidates.first > 0) ? candidates.first : candidates.second;
613 static AliFemtoThreeVector WestEnd(0.,0.,200.); // bugfix MAL jul00 - 200cm NOT 2cm
614 static AliFemtoThreeVector EastEnd(0.,0.,-200.); // bugfix MAL jul00 - 200cm NOT 2cm
615 static AliFemtoThreeVector EndCapNormal(0.,0.,1.0);
617 endLength = hel.pathLength(WestEnd,EndCapNormal);
618 if (endLength < 0.0) endLength = hel.pathLength(EastEnd,EndCapNormal);
620 if (endLength < 0.0) cout <<
621 "AliFemtoParticle::CalculateTpcExitAndEntrancePoints(): "
622 << "Hey -- I cannot find an exit point out endcaps" << endl;
623 // OK, firstExitLength will be the shortest way out of the detector...
624 double firstExitLength = (endLength < sideLength) ? endLength : sideLength;
625 // now then, let's return the POSITION at which particle leaves TPC...
626 *tmpTpcExitPoint = hel.at(firstExitLength);
627 // Finally, calculate the position at which the track crosses the inner field cage
628 candidates = hel.pathLength(50.0); // bugfix MAL jul00 - 200cm NOT 2cm
630 sideLength = (candidates.first > 0) ? candidates.first : candidates.second;
631 // cout << "sideLength 2 ="<<sideLength << endl;
632 *tmpTpcEntrancePoint = hel.at(sideLength);
633 // This is the secure way !
634 if (IsNaN(tmpTpcEntrancePoint->x()) ||
635 IsNaN(tmpTpcEntrancePoint->y()) ||
636 IsNaN(tmpTpcEntrancePoint->z()) ){
637 cout << "tmpTpcEntrancePoint NAN"<< endl;
638 cout << "tmpNominalTpcEntrancePoint = " <<tmpTpcEntrancePoint<< endl;
639 tmpTpcEntrancePoint->setX(-9999.);
640 tmpTpcEntrancePoint->setY(-9999.);
641 tmpTpcEntrancePoint->setZ(-9999.);
644 if (IsNaN(tmpTpcExitPoint->x()) ||
645 IsNaN(tmpTpcExitPoint->y()) ||
646 IsNaN(tmpTpcExitPoint->z()) ) {
647 // cout << "tmpTpcExitPoint NAN set at (-9999,-9999,-9999)"<< endl;
648 // cout << "tmpTpcExitPoint X= " <<tmpTpcExitPoint->x()<< endl;
649 // cout << "tmpTpcExitPoint Y= " <<tmpTpcExitPoint->y()<< endl;
650 // cout << "tmpTpcExitPoint Z= " <<tmpTpcExitPoint->z()<< endl;
651 tmpTpcExitPoint->setX(-9999.);
652 tmpTpcExitPoint->setY(-9999.);
653 tmpTpcExitPoint->setZ(-9999.);
657 // if (IsNaN(tmpTpcExitPoint->x())) *tmpTpcExitPoint = AliFemtoThreeVector(-9999.,-9999.,-9999);
658 // if (IsNaN(tmpTpcEntrancetPoint->x())) *tmpTpcEntrancePoint = AliFemtoThreeVector(-9999.,-9999.,-9999);
659 // cout << "tmpTpcEntrancePoint"<<*tmpTpcEntrancePoint << endl;
661 // 03Oct00 - mal. OK, let's try something a little more
662 // along the lines of NA49 and E895 strategy.
663 // calculate the "nominal" position at N radii (say N=11)
664 // within the TPC, and for a pair cut
665 // use the average separation of these N
667 candidates = hel.pathLength(50.0);
668 sideLength = (candidates.first > 0) ? candidates.first : candidates.second;
669 while (irad<11 && !IsNaN(sideLength)){
670 float radius = 50.0 + irad*15.0;
671 candidates = hel.pathLength(radius);
672 sideLength = (candidates.first > 0) ? candidates.first : candidates.second;
673 tmpPosSample[irad] = hel.at(sideLength);
674 if(IsNaN(tmpPosSample[irad].x()) ||
675 IsNaN(tmpPosSample[irad].y()) ||
676 IsNaN(tmpPosSample[irad].z())
678 cout << "tmpPosSample for radius=" << radius << " NAN"<< endl;
679 cout << "tmpPosSample=(" <<tmpPosSample[irad]<<")"<< endl;
680 tmpPosSample[irad] = AliFemtoThreeVector(-9999.,-9999.,-9999);
684 float radius = 50.0 + irad*15.0;
685 candidates = hel.pathLength(radius);
686 sideLength = (candidates.first > 0) ? candidates.first : candidates.second;
689 for (int i = irad; i<11; i++)
691 tmpPosSample[i] = AliFemtoThreeVector(-9999.,-9999.,-9999);
694 static float tRowRadius[45] = {60,64.8,69.6,74.4,79.2,84,88.8,93.6,98.8,
695 104,109.2,114.4,119.6,127.195,129.195,131.195,
696 133.195,135.195,137.195,139.195,141.195,
697 143.195,145.195,147.195,149.195,151.195,
698 153.195,155.195,157.195,159.195,161.195,
699 163.195,165.195,167.195,169.195,171.195,
700 173.195,175.195,177.195,179.195,181.195,
701 183.195,185.195,187.195,189.195};
702 int tRow,tSect,tOutOfBound;
705 AliFemtoThreeVector tPoint;
706 AliFmThreeVectorD tn(0,0,0);
707 AliFmThreeVectorD tr(0,0,0);
709 // test to enter the loop
710 candidates = hel.pathLength(tRowRadius[ti]);
711 tLength = (candidates.first > 0) ? candidates.first : candidates.second;
713 cout <<"tLength Init tmp NAN" << endl;
714 cout <<"padrow number= "<<ti << "not reached" << endl;
715 cout << "*** DO NOT ENTER THE LOOP***" << endl;
716 tmpSect[ti]=-1;//sector
719 while(ti<45 && !IsNaN(tLength)){
720 candidates = hel.pathLength(tRowRadius[ti]);
721 tLength = (candidates.first > 0) ? candidates.first : candidates.second;
723 cout <<"tLength loop 1st NAN" << endl;
724 cout <<"padrow number= " << ti << " not reached" << endl;
725 cout << "*** THIS IS AN ERROR SHOULDN'T LOOP ***" << endl;
726 tmpSect[ti]=-1;//sector
728 tPoint = hel.at(tLength);
729 // Find which sector it is on
730 TpcLocalTransform(tPoint,tmpSect[ti],tRow,tU,tPhi);
731 if (IsNaN(tmpSect[ti])){
732 cout <<"***ERROR tmpSect"<< endl;
735 cout <<"***ERROR tRow"<< endl;
738 cout <<"***ERROR tU"<< endl;
741 cout <<"***ERROR tPhi"<< endl;
743 // calculate crossing plane
746 tr.setX(tRowRadius[ti]*cos(tPhi));
747 tr.setY(tRowRadius[ti]*sin(tPhi));
748 // find crossing point
749 tLength = hel.pathLength(tr,tn);
751 cout <<"tLength loop 2nd NAN" << endl;
752 cout <<"padrow number= " << ti << " not reached" << endl;
753 tmpSect[ti]=-2;//sector
755 tPoint = hel.at(tLength);
756 tmpZ[ti] = tPoint.z();
757 tOutOfBound = TpcLocalTransform(tPoint,tSect,tRow,tmpU[ti],tPhi);
759 cout <<"***ERROR tSect 2"<< endl;
762 cout <<"***ERROR tRow 2"<< endl;
764 if (IsNaN(tmpU[ti])){
765 cout <<"***ERROR tmpU[ti] 2"<< endl;
768 cout <<"***ERROR tPhi 2 "<< endl;
770 if(tOutOfBound || (tmpSect[ti] == tSect && tRow!=(ti+1))){
772 // cout << "missed once"<< endl;
775 if(tmpSect[ti] != tSect){
776 // Try again on the other sector
779 tr.setX(tRowRadius[ti]*cos(tPhi));
780 tr.setY(tRowRadius[ti]*sin(tPhi));
781 // find crossing point
782 tLength = hel.pathLength(tr,tn);
783 tPoint = hel.at(tLength);
785 cout <<"tLength loop 3rd NAN" << endl;
786 cout <<"padrow number= "<< ti << " not reached" << endl;
787 tmpSect[ti]=-1;//sector
789 tmpZ[ti] = tPoint.z();
791 tOutOfBound = TpcLocalTransform(tPoint,tSect,tRow,tmpU[ti],tPhi);
793 cout <<"***ERROR tSect 3"<< endl;
796 cout <<"***ERROR tRow 3"<< endl;
798 if (IsNaN(tmpU[ti])){
799 cout <<"***ERROR tmpU[ti] 3"<< endl;
802 cout <<"***ERROR tPhi 3 "<< endl;
804 if(tOutOfBound || tSect!= tmpSect[ti] || tRow!=(ti+1)){
809 if (IsNaN(tmpSect[ti])){
810 cout << "*******************ERROR***************************" << endl;
811 cout <<"AliFemtoParticle--Fctn tmpSect=" << tmpSect[ti] << endl;
812 cout << "*******************ERROR***************************" << endl;
814 if (IsNaN(tmpU[ti])){
815 cout << "*******************ERROR***************************" << endl;
816 cout <<"AliFemtoParticle--Fctn tmpU=" << tmpU[ti] << endl;
817 cout << "*******************ERROR***************************" << endl;
819 if (IsNaN(tmpZ[ti])){
820 cout << "*******************ERROR***************************" << endl;
821 cout <<"AliFemtoParticle--Fctn tmpZ=" << tmpZ[ti] << endl;
822 cout << "*******************ERROR***************************" << endl;
824 // If padrow ti not reached all other beyond are not reached
825 // in this case set sector to -1
826 if (tmpSect[ti]==-1){
827 for (int tj=ti; tj<45;tj++){
834 candidates = hel.pathLength(tRowRadius[ti]);
835 tLength = (candidates.first > 0) ? candidates.first : candidates.second;}
838 //_____________________
839 const AliFemtoThreeVector& AliFemtoParticle::TpcV0PosExitPoint() const{
840 return fTpcV0PosExitPoint;
842 //_____________________
843 const AliFemtoThreeVector& AliFemtoParticle::TpcV0PosEntrancePoint() const{
844 return fTpcV0PosEntrancePoint;
846 //______________________
847 const AliFemtoThreeVector& AliFemtoParticle::TpcV0NegExitPoint() const{
848 return fTpcV0NegExitPoint;
850 //_____________________
851 const AliFemtoThreeVector& AliFemtoParticle::TpcV0NegEntrancePoint() const{
852 return fTpcV0NegEntrancePoint;
854 //______________________