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 //______________________