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