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