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   for (int ip=0; ip<6; ip++) fPurity[ip] = 0.0;
  58   //  cout << "Created particle " << this << endl;
  59 }
  60 //_____________________
  61 AliFemtoParticle::AliFemtoParticle(const AliFemtoParticle& aParticle):
  62 //   fTpcV0NegPosSample(0),
  63 //   fV0NegZ(0),
  64 //   fV0NegU(0),
  65 //   fV0NegSect(0),
  66   fTrack(0), fV0(0), fKink(0), fXi(0),
  67   fFourMomentum(0),
  68   fHelix(),
  69 //   fNominalTpcExitPoint(0),
  70 //   fNominalTpcEntrancePoint(0),
  71   fHiddenInfo(0),
  72   fPrimaryVertex(0),
  73   fSecondaryVertex(0),
  74   fHelixV0Pos(),
  75   fTpcV0PosEntrancePoint(0),
  76   fTpcV0PosExitPoint(0),
  77   fHelixV0Neg(),
  78   fTpcV0NegEntrancePoint(0),
  79   fTpcV0NegExitPoint(0)
  80 {
  81   // Copy constructor
  82   if (aParticle.fTrack)
  83     fTrack = new AliFemtoTrack(*aParticle.fTrack);
  84   if (aParticle.fV0)
  85     fV0    = new AliFemtoV0(*aParticle.fV0);
  86   if (aParticle.fKink)
  87     fKink  = new AliFemtoKink(*aParticle.fKink);
  88   if (aParticle.fXi)
  89     fXi    = new AliFemtoXi(*aParticle.fXi);
  90   fFourMomentum = aParticle.fFourMomentum;
  91   fHelix = aParticle.fHelix;
  92   for (int ip=0; ip<6; ip++) fPurity[ip] = 0.0;
  93
  94 //   for (int iter=0; iter<11; iter++)
  95 //     fNominalPosSample[iter] = aParticle.fNominalPosSample[iter];
  96
  97 //   if (aParticle.fTpcV0NegPosSample) {
  98 //     fTpcV0NegPosSample = (AliFemtoThreeVector *) malloc(sizeof(AliFemtoThreeVector) * 11);
  99 //     for (int iter=0; iter<11; iter++)
 100 //       fTpcV0NegPosSample[iter] = aParticle.fTpcV0NegPosSample[iter];
 101 //   }
 102
 103 //   if (aParticle.fV0NegZ) {
 104 //     fV0NegZ = (float *) malloc(sizeof(float) * 45);
 105 //     for (int iter=0; iter<11; iter++)
 106 //       fV0NegZ[iter] = aParticle.fV0NegZ[iter];
 107 //   }
 108 //   if (aParticle.fV0NegU) {
 109 //     fV0NegU = (float *) malloc(sizeof(float) * 45);
 110 //     for (int iter=0; iter<11; iter++)
 111 //       fV0NegU[iter] = aParticle.fV0NegU[iter];
 112 //   }
 113 //   if (aParticle.fV0NegSect) {
 114 //     fV0NegSect = (int *) malloc(sizeof(int) * 45);
 115 //     for (int iter=0; iter<11; iter++)
 116 //       fV0NegSect[iter] = aParticle.fV0NegSect[iter];
 117 //   }
 118
 119   fPrimaryVertex = aParticle.fPrimaryVertex;
 120   fSecondaryVertex = aParticle.fSecondaryVertex;
 121   CalculatePurity();
 122   if(aParticle.fHiddenInfo){
 123     fHiddenInfo= aParticle.HiddenInfo()->Clone();
 124   }
 125
 126 //   fNominalTpcEntrancePoint = aParticle.fNominalTpcEntrancePoint;
 127 //   fNominalTpcExitPoint     = aParticle.fNominalTpcExitPoint;
 128
 129   for (int iter=0; iter<6; iter++)
 130     fPurity[iter] = aParticle.fPurity[iter];
 131
 132   fHelixV0Pos = aParticle.fHelixV0Pos;
 133   fTpcV0PosEntrancePoint = aParticle.fTpcV0PosEntrancePoint;
 134   fTpcV0PosExitPoint     = aParticle.fTpcV0PosExitPoint;
 135   fHelixV0Neg = aParticle.fHelixV0Neg;
 136   fTpcV0NegEntrancePoint = aParticle.fTpcV0NegEntrancePoint;
 137   fTpcV0NegExitPoint     = aParticle.fTpcV0NegExitPoint;
 138 }
 139 //_____________________
 140 AliFemtoParticle::~AliFemtoParticle(){
 141   //  cout << "Issuing delete for AliFemtoParticle." << endl;
 142
 143   if (fTrack) delete fTrack;
 144   if (fV0) {
 145 //     delete[] fTpcV0NegPosSample;
 146 //     delete[] fV0NegZ;
 147 //     delete[] fV0NegU;
 148 //     delete[] fV0NegSect;
 149     delete fV0;
 150   }
 151   if (fKink) delete fKink;
 152   if (fXi) delete fXi;
 153   //  cout << "Trying to delete HiddenInfo: " << fHiddenInfo << endl;
 154   if (fHiddenInfo)
 155     {
 156       //      cout << "Deleting HiddenInfo." << endl;
 157       delete fHiddenInfo;
 158     }
 159   //  cout << "Deleted particle " << this << endl;
 160 }
 161 //_____________________
 162 AliFemtoParticle::AliFemtoParticle(const AliFemtoTrack* const hbtTrack,const double& mass) :
 163 //   fTpcV0NegPosSample(0),
 164 //   fV0NegZ(0),
 165 //   fV0NegU(0),
 166 //   fV0NegSect(0),
 167   fTrack(0), fV0(0), fKink(0), fXi(0),
 168   fFourMomentum(0),
 169   fHelix(),
 170 //   fNominalTpcExitPoint(0),
 171 //   fNominalTpcEntrancePoint(0),
 172   fHiddenInfo(0),
 173   fPrimaryVertex(0),
 174   fSecondaryVertex(0),
 175   fHelixV0Pos(),
 176   fTpcV0PosEntrancePoint(0),
 177   fTpcV0PosExitPoint(0),
 178   fHelixV0Neg(),
 179   fTpcV0NegEntrancePoint(0),
 180   fTpcV0NegExitPoint(0)
 181 {
 182   // Constructor from normal track
 183
 184   // I know there is a better way to do this...
 185   fTrack = new AliFemtoTrack(*hbtTrack);
 186   AliFemtoThreeVector temp = hbtTrack->P();
 187   fFourMomentum.SetVect(temp);
 188   double ener = ::sqrt(temp.Mag2()+mass*mass);
 189   fFourMomentum.SetE(ener);
 190 //  fMap[0] = hbtTrack->TopologyMap(0);
 191  // fMap[1] = hbtTrack->TopologyMap(1);
 192  // fNhits = hbtTrack->NHits();
 193   fHelix = hbtTrack->Helix();
 194   //CalculateNominalTpcExitAndEntrancePoints();
 195
 196
 197   fPrimaryVertex.SetX(0.);
 198   fPrimaryVertex.SetY(0.);
 199   fPrimaryVertex.SetZ(0.);
 200   fSecondaryVertex.SetX(0.);
 201   fSecondaryVertex.SetY(0.);
 202   fSecondaryVertex.SetZ(0.);
 203   /* TO JA ODZNACZYLEM NIE WIEM DLACZEGO
 204   CalculateTpcExitAndEntrancePoints(&fHelix,&fPrimaryVertex,
 205                                     &fSecondaryVertex,
 206                                     &fNominalTpcEntrancePoint,
 207                                     &fNominalTpcExitPoint,
 208                                     &mNominalPosSample[0],
 209                                     &fZ[0],
 210                                     &fU[0],
 211                                     &fSect[0]);
 212   */
 213   CalculatePurity();
 214   // ***
 215   fHiddenInfo= 0;
 216   if(hbtTrack->ValidHiddenInfo()){
 217     fHiddenInfo= hbtTrack->GetHiddenInfo()->Clone();
 218   }
 219   // ***
 220   //  cout << "Created particle " << this << endl;
 221
 222 }
 223 //_____________________
 224 AliFemtoParticle::AliFemtoParticle(const AliFemtoV0* const hbtV0,const double& mass) :
 225 //   fTpcV0NegPosSample(0),
 226 //   fV0NegZ(0),
 227 //   fV0NegU(0),
 228 //   fV0NegSect(0),
 229   fTrack(0), fV0(0), fKink(0),  fXi(0),
 230   fFourMomentum(0),
 231   fHelix(),
 232 //   fNominalTpcExitPoint(0),
 233 //   fNominalTpcEntrancePoint(0),
 234   fHiddenInfo(0),
 235   fPrimaryVertex(0),
 236   fSecondaryVertex(0),
 237   fHelixV0Pos(),
 238   fTpcV0PosEntrancePoint(0),
 239   fTpcV0PosExitPoint(0),
 240   fHelixV0Neg(),
 241   fTpcV0NegEntrancePoint(0),
 242   fTpcV0NegExitPoint(0)
 243 {
 244   // Constructor from V0
 245   fV0 = new AliFemtoV0(*hbtV0);
 246  //fMap[0]= 0;
 247   //fMap[1]= 0;
 248   // I know there is a better way to do this...
 249   AliFemtoThreeVector temp = hbtV0->MomV0();
 250   fFourMomentum.SetVect(temp);
 251   double ener = ::sqrt(temp.Mag2()+mass*mass);
 252   fFourMomentum.SetE(ener);
 253   // Calculating TpcEntrancePoint for Positive V0 daugther
 254   fPrimaryVertex = hbtV0->PrimaryVertex();
 255   fSecondaryVertex = hbtV0->DecayVertexV0();
 256   fHelixV0Pos = hbtV0->HelixPos();
 257
 258 //   fTpcV0NegPosSample = new AliFemtoThreeVector[45];//for V0Neg
 259 //   fV0NegZ = new float[45];//for V0Neg
 260 //   fV0NegU = new float[45];//for V0Neg
 261 //   fV0NegSect = new int[45];//for V0Neg
 262 //   CalculateTpcExitAndEntrancePoints(&fHelixV0Pos,&fPrimaryVertex,
 263 //                                  &fSecondaryVertex,
 264 //                                  &fTpcV0PosEntrancePoint,
 265 //                                  &fTpcV0PosExitPoint,
 266 //                                  &fNominalPosSample[0],
 267 //                                  &fZ[0],
 268 //                                  &fU[0],&fSect[0]);
 269   fHelixV0Neg = hbtV0->HelixNeg();
 270
 271 //   CalculateTpcExitAndEntrancePoints(&fHelixV0Neg,
 272 //                                  &fPrimaryVertex,
 273 //                                  &fSecondaryVertex,
 274 //                                  &fTpcV0NegEntrancePoint,
 275 //                                  &fTpcV0NegExitPoint,
 276 //                                  &fTpcV0NegPosSample[0],
 277 //                                  &fV0NegZ[0],
 278 //                                  &fV0NegU[0],&fV0NegSect[0]);
 279
 280   // ***
 281   fHiddenInfo= 0;
 282   if(hbtV0->ValidHiddenInfo()){
 283     fHiddenInfo= hbtV0->GetHiddenInfo()->Clone();
 284   }
 285   for (int ip=0; ip<6; ip++) fPurity[ip] = 0.0;
 286   // ***
 287 }
 288 //_____________________
 289 AliFemtoParticle::AliFemtoParticle(const AliFemtoKink* const hbtKink,const double& mass) :
 290 //   fTpcV0NegPosSample(0),
 291 //   fV0NegZ(0),
 292 //   fV0NegU(0),
 293 //   fV0NegSect(0),
 294   fTrack(0), fV0(0), fKink(0), fXi(0),
 295   fFourMomentum(0),
 296   fHelix(),
 297 //   fNominalTpcExitPoint(0),
 298 //   fNominalTpcEntrancePoint(0),
 299   fHiddenInfo(0),
 300   fPrimaryVertex(0),
 301   fSecondaryVertex(0),
 302   fHelixV0Pos(),
 303   fTpcV0PosEntrancePoint(0),
 304   fTpcV0PosExitPoint(0),
 305   fHelixV0Neg(),
 306   fTpcV0NegEntrancePoint(0),
 307   fTpcV0NegExitPoint(0)
 308 {
 309   // Constructor from Kink
 310   fKink = new AliFemtoKink(*hbtKink);
 311  // fMap[0]= 0;
 312   //fMap[1]= 0;
 313   // I know there is a better way to do this...
 314   AliFemtoThreeVector temp = hbtKink->Parent().P();
 315   fFourMomentum.SetVect(temp);
 316   double ener = ::sqrt(temp.Mag2()+mass*mass);
 317   fFourMomentum.SetE(ener);
 318   for (int ip=0; ip<6; ip++) fPurity[ip] = 0.0;
 319 }
 320
 321 //_____________________
 322 AliFemtoParticle::AliFemtoParticle(const AliFemtoXi* const hbtXi, const double& mass) :
 323 //   fTpcV0NegPosSample(0),
 324 //   fV0NegZ(0),
 325 //   fV0NegU(0),
 326 //   fV0NegSect(0),
 327   fTrack(0), fV0(0), fKink(0), fXi(0),
 328   fFourMomentum(0),
 329   fHelix(),
 330 //   fNominalTpcExitPoint(0),
 331 //   fNominalTpcEntrancePoint(0),
 332   fHiddenInfo(0),
 333   fPrimaryVertex(0),
 334   fSecondaryVertex(0),
 335   fHelixV0Pos(),
 336   fTpcV0PosEntrancePoint(0),
 337   fTpcV0PosExitPoint(0),
 338   fHelixV0Neg(),
 339   fTpcV0NegEntrancePoint(0),
 340   fTpcV0NegExitPoint(0)
 341 {
 342   // Constructor from Xi
 343   fXi = new AliFemtoXi(*hbtXi);
 344  // fMap[0]= 0;
 345   //fMap[1]= 0;
 346   AliFemtoThreeVector temp;// = hbtXi->mMofXi;
 347   fFourMomentum.SetVect(temp);
 348   double ener = ::sqrt(temp.Mag2()+mass*mass);
 349   fFourMomentum.SetE(ener);
 350   fHiddenInfo = 0;
 351   for (int ip=0; ip<6; ip++) fPurity[ip] = 0.0;
 352
 353 }
 354 //_____________________
 355 AliFemtoParticle& AliFemtoParticle::operator=(const AliFemtoParticle& aParticle)
 356 {
 357   // assignment operator
 358   if (this != &aParticle) {
 359
 360     if (fTrack) delete fTrack;
 361     if (aParticle.fTrack) {
 362       fTrack = new AliFemtoTrack(*aParticle.fTrack);
 363       CalculatePurity();
 364     }
 365     if (fV0) delete fV0;
 366     if (aParticle.fV0)
 367       fV0    = new AliFemtoV0(*aParticle.fV0);
 368     if (fKink) delete fKink;
 369     if (aParticle.fKink)
 370       fKink  = new AliFemtoKink(*aParticle.fKink);
 371     if (fXi) delete fXi;
 372     if (aParticle.fXi)
 373       fXi    = new AliFemtoXi(*aParticle.fXi);
 374     fFourMomentum = aParticle.fFourMomentum;
 375     fHelix = aParticle.fHelix;
 376
 377     //   for (int iter=0; iter<11; iter++)
 378     //     fNominalPosSample[iter] = aParticle.fNominalPosSample[iter];
 379
 380     //   if (fTpcV0NegPosSample) delete fTpcV0NegPosSample;
 381     //   if (aParticle.fTpcV0NegPosSample) {
 382     //     fTpcV0NegPosSample = (AliFemtoThreeVector *) malloc(sizeof(AliFemtoThreeVector) * 11);
 383     //     for (int iter=0; iter<11; iter++)
 384     //       fTpcV0NegPosSample[iter] = aParticle.fTpcV0NegPosSample[iter];
 385     //   }
 386
 387     //   if (fV0NegZ) delete fV0NegZ;
 388     //   if (aParticle.fV0NegZ) {
 389     //     fV0NegZ = (float *) malloc(sizeof(float) * 45);
 390     //     for (int iter=0; iter<11; iter++)
 391     //       fV0NegZ[iter] = aParticle.fV0NegZ[iter];
 392     //   }
 393     //   if (fV0NegU) delete fV0NegU;
 394     //   if (aParticle.fV0NegU) {
 395     //     fV0NegU = (float *) malloc(sizeof(float) * 45);
 396     //     for (int iter=0; iter<11; iter++)
 397     //       fV0NegU[iter] = aParticle.fV0NegU[iter];
 398     //   }
 399     //   if (fV0NegSect) delete fV0NegSect;
 400     //   if (aParticle.fV0NegSect) {
 401     //     fV0NegSect = (int *) malloc(sizeof(int) * 45);
 402     //     for (int iter=0; iter<11; iter++)
 403     //       fV0NegSect[iter] = aParticle.fV0NegSect[iter];
 404     //   }
 405
 406     fPrimaryVertex = aParticle.fPrimaryVertex;
 407     fSecondaryVertex = aParticle.fSecondaryVertex;
 408 //     if (fHiddenInfo) delete fHiddenInfo;
 409 //     if(aParticle.fHiddenInfo){
 410 //       fHiddenInfo= aParticle.fHiddenInfo->Clone();
 411 //     }
 412
 413     //   fNominalTpcEntrancePoint = aParticle.fNominalTpcEntrancePoint;
 414     //   fNominalTpcExitPoint     = aParticle.fNominalTpcExitPoint;
 415
 416     if (fHiddenInfo) delete fHiddenInfo;
 417     if (aParticle.fHiddenInfo)
 418       fHiddenInfo = aParticle.HiddenInfo()->Clone();
 419
 420     for (int iter=0; iter<6; iter++)
 421       fPurity[iter] = aParticle.fPurity[iter];
 422
 423     fHelixV0Pos = aParticle.fHelixV0Pos;
 424     fTpcV0PosEntrancePoint = aParticle.fTpcV0PosEntrancePoint;
 425     fTpcV0PosExitPoint     = aParticle.fTpcV0PosExitPoint;
 426     fHelixV0Neg = aParticle.fHelixV0Neg;
 427     fTpcV0NegEntrancePoint = aParticle.fTpcV0NegEntrancePoint;
 428     fTpcV0NegExitPoint     = aParticle.fTpcV0NegExitPoint;
 429
 430   }
 431
 432   return *this;
 433 }
 434 // //_____________________
 435 // const AliFemtoThreeVector& AliFemtoParticle::NominalTpcExitPoint() 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 fNominalTpcExitPoint;
 439 // }
 440 // //_____________________
 441 // const AliFemtoThreeVector& AliFemtoParticle::NominalTpcEntrancePoint() const{
 442 //   // in future, may want to calculate this "on demand" only, sot this routine may get more sophisticated
 443 //   // for now, we calculate Exit and Entrance points upon instantiation
 444 //   return fNominalTpcEntrancePoint;
 445 // }
 446 //_____________________
 447 void AliFemtoParticle::CalculatePurity(){
 448   // Calculate additional parameterized purity
 449
 450   double tPt = fFourMomentum.Perp();
 451   // pi -
 452   fPurity[0] = fgPrimPimPar0*(1.-exp((tPt-fgPrimPimPar1)/fgPrimPimPar2));
 453   fPurity[0] *= fTrack->PidProbPion();
 454   // pi+
 455   fPurity[1] = fgPrimPipPar0*(1.-exp((tPt-fgPrimPipPar1)/fgPrimPipPar2));
 456   fPurity[1] *= fTrack->PidProbPion();
 457   // K-
 458   fPurity[2] = fTrack->PidProbKaon();
 459   // K+
 460   fPurity[3] = fTrack->PidProbKaon();
 461   // pbar
 462   fPurity[4] = fTrack->PidProbProton();
 463   // p
 464   fPurity[5] = fTrack->PidProbProton();
 465 }
 466
 467 double AliFemtoParticle::GetPionPurity()
 468 {
 469   // Get full pion purity
 470   if (fTrack->Charge()>0)
 471     return fPurity[1];
 472   else
 473     return fPurity[0];
 474 }
 475 double AliFemtoParticle::GetKaonPurity()
 476 {
 477   // Get full kaon purity
 478   if (fTrack->Charge()>0)
 479     return fPurity[3];
 480   else
 481     return fPurity[2];
 482 }
 483 double AliFemtoParticle::GetProtonPurity()
 484 {
 485   // Get full proton purity
 486   if (fTrack->Charge()>0)
 487     return fPurity[5];
 488   else
 489     return fPurity[4];
 490 }
 491
 492 // void AliFemtoParticle::CalculateTpcExitAndEntrancePoints(AliFmPhysicalHelixD* tHelix,
 493 //                                                     AliFemtoThreeVector*  PrimVert,
 494 //                                                     AliFemtoThreeVector*  SecVert,
 495 //                                                     AliFemtoThreeVector* tmpTpcEntrancePoint,
 496 //                                                     AliFemtoThreeVector* tmpTpcExitPoint,
 497 //                                                     AliFemtoThreeVector* tmpPosSample,
 498 //                                                     float* tmpZ,
 499 //                                                     float* tmpU,
 500 //                                                     int* tmpSect){
 501 //   // this calculates the exit point of a secondary track,
 502 //   // either through the endcap or through the Outer Field Cage
 503 //   // We assume the track to start at tHelix.origin-PrimaryVertex
 504 //   // it also calculates the entrance point of the secondary track,
 505 //   // which is the point at which it crosses the
 506 //   // inner field cage
 507 //   //  static AliFemtoThreeVector ZeroVec(0.,0.,0.);
 508 //   AliFemtoThreeVector tZeroVec(0.,0.,0.);
 509 // //   tZeroVec.SetX(tHelix->origin().x()-PrimVert->x());
 510 // //   tZeroVec.SetY(tHelix->origin().y()-PrimVert->y());
 511 // //   tZeroVec.SetZ(tHelix->origin().z()-PrimVert->z());
 512 //   tZeroVec.SetX(SecVert->x()-PrimVert->x());
 513 //   tZeroVec.SetY(SecVert->y()-PrimVert->y());
 514 //   tZeroVec.SetZ(SecVert->z()-PrimVert->z());
 515 //   double dip, curv, phase;
 516 //   int h;
 517 //   curv = tHelix->Curvature();
 518 //   dip  = tHelix->DipAngle();
 519 //   phase= tHelix->Phase();
 520 //   h    = tHelix->H();
 521
 522 //   AliFmHelixD hel(curv,dip,phase,tZeroVec,h);
 523
 524 //   pairD candidates;
 525 //   double sideLength;  // this is how much length to go to leave through sides of TPC
 526 //   double endLength;  // this is how much length to go to leave through endcap of TPC
 527 //   // figure out how far to go to leave through side...
 528 //   candidates = hel.PathLength(200.0);  // bugfix MAL jul00 - 200cm NOT 2cm
 529 //   sideLength = (candidates.first > 0) ? candidates.first : candidates.second;
 530
 531 //   static AliFemtoThreeVector tWestEnd(0.,0.,200.);  // bugfix MAL jul00 - 200cm NOT 2cm
 532 //   static AliFemtoThreeVector tEastEnd(0.,0.,-200.); // bugfix MAL jul00 - 200cm NOT 2cm
 533 //   static AliFemtoThreeVector tEndCapNormal(0.,0.,1.0);
 534
 535 //   endLength = hel.PathLength(tWestEnd,tEndCapNormal);
 536 //   if (endLength < 0.0) endLength = hel.PathLength(tEastEnd,tEndCapNormal);
 537
 538 //   if (endLength < 0.0) cout <<
 539 //                       "AliFemtoParticle::CalculateTpcExitAndEntrancePoints(): "
 540 //                             << "Hey -- I cannot find an exit point out endcaps" << endl;
 541 //   // OK, firstExitLength will be the shortest way out of the detector...
 542 //   double firstExitLength = (endLength < sideLength) ? endLength : sideLength;
 543 //   // now then, let's return the POSITION at which particle leaves TPC...
 544 //   *tmpTpcExitPoint = hel.At(firstExitLength);
 545 //   // Finally, calculate the position at which the track crosses the inner field cage
 546 //   candidates = hel.PathLength(50.0);  // bugfix MAL jul00 - 200cm NOT 2cm
 547
 548 //   sideLength = (candidates.first > 0) ? candidates.first : candidates.second;
 549 //   //  cout << "sideLength 2 ="<<sideLength << endl;
 550 //   *tmpTpcEntrancePoint = hel.At(sideLength);
 551 //   // This is the secure way !
 552 //   if (IsNaN(tmpTpcEntrancePoint->x()) ||
 553 //       IsNaN(tmpTpcEntrancePoint->y()) ||
 554 //       IsNaN(tmpTpcEntrancePoint->z()) ){
 555 //     cout << "tmpTpcEntrancePoint NAN"<< endl;
 556 //     cout << "tmpNominalTpcEntrancePoint = " <<tmpTpcEntrancePoint<< endl;
 557 //     tmpTpcEntrancePoint->SetX(-9999.);
 558 //     tmpTpcEntrancePoint->SetY(-9999.);
 559 //     tmpTpcEntrancePoint->SetZ(-9999.);
 560 //   }
 561
 562 //   if (IsNaN(tmpTpcExitPoint->x()) ||
 563 //       IsNaN(tmpTpcExitPoint->y()) ||
 564 //       IsNaN(tmpTpcExitPoint->z()) ) {
 565 // //     cout << "tmpTpcExitPoint NAN Set at (-9999,-9999,-9999)"<< endl;
 566 // //     cout << "tmpTpcExitPoint X= " <<tmpTpcExitPoint->x()<< endl;
 567 // //     cout << "tmpTpcExitPoint Y= " <<tmpTpcExitPoint->y()<< endl;
 568 // //     cout << "tmpTpcExitPoint Z= " <<tmpTpcExitPoint->z()<< endl;
 569 //     tmpTpcExitPoint->SetX(-9999.);
 570 //     tmpTpcExitPoint->SetY(-9999.);
 571 //     tmpTpcExitPoint->SetZ(-9999.);
 572 //   }
 573
 574
 575 // //   if (IsNaN(tmpTpcExitPoint->x())) *tmpTpcExitPoint = AliFemtoThreeVector(-9999.,-9999.,-9999);
 576 // //   if (IsNaN(tmpTpcEntrancetPoint->x())) *tmpTpcEntrancePoint = AliFemtoThreeVector(-9999.,-9999.,-9999);
 577 //   //  cout << "tmpTpcEntrancePoint"<<*tmpTpcEntrancePoint << endl;
 578
 579 //   // 03Oct00 - mal.  OK, let's try something a little more
 580 //   // along the lines of NA49 and E895 strategy.
 581 //   // calculate the "nominal" position at N radii (say N=11)
 582 //   // within the TPC, and for a pair cut
 583 //   // use the average separation of these N
 584 //   int irad = 0;
 585 //   candidates = hel.PathLength(50.0);
 586 //   sideLength = (candidates.first > 0) ? candidates.first : candidates.second;
 587 //   while (irad<11 && !IsNaN(sideLength)){
 588 //     float radius = 50.0 + irad*15.0;
 589 //     candidates = hel.PathLength(radius);
 590 //     sideLength = (candidates.first > 0) ? candidates.first : candidates.second;
 591 //     tmpPosSample[irad] = hel.At(sideLength);
 592 //     if(IsNaN(tmpPosSample[irad].x()) ||
 593 //        IsNaN(tmpPosSample[irad].y()) ||
 594 //        IsNaN(tmpPosSample[irad].z())
 595 //        ){
 596 //       cout << "tmpPosSample for radius=" << radius << " NAN"<< endl;
 597 //       cout << "tmpPosSample=(" <<tmpPosSample[irad]<<")"<< endl;
 598 //       tmpPosSample[irad] =  AliFemtoThreeVector(-9999.,-9999.,-9999);
 599 //     }
 600 //     irad++;
 601 //     if (irad<11){
 602 //       float radius = 50.0 + irad*15.0;
 603 //       candidates = hel.PathLength(radius);
 604 //       sideLength = (candidates.first > 0) ? candidates.first : candidates.second;
 605 //     }
 606 //    }
 607 //    for (int i = irad; i<11; i++)
 608 //      {
 609 //        tmpPosSample[i] =  AliFemtoThreeVector(-9999.,-9999.,-9999);
 610 //      }
 611
 612 //   static float tRowRadius[45] = {60,64.8,69.6,74.4,79.2,84,88.8,93.6,98.8,
 613 //                               104,109.2,114.4,119.6,127.195,129.195,131.195,
 614 //                               133.195,135.195,137.195,139.195,141.195,
 615 //                               143.195,145.195,147.195,149.195,151.195,
 616 //                               153.195,155.195,157.195,159.195,161.195,
 617 //                               163.195,165.195,167.195,169.195,171.195,
 618 //                               173.195,175.195,177.195,179.195,181.195,
 619 //                               183.195,185.195,187.195,189.195};
 620 //   int tRow,tSect,tOutOfBound;
 621 //   double tLength,tPhi;
 622 //   float tU;
 623 //   AliFemtoThreeVector tPoint;
 624 //   AliFmThreeVectorD tn(0,0,0);
 625 //   AliFmThreeVectorD tr(0,0,0);
 626 //   int ti =0;
 627 //   // test to enter the loop
 628 //   candidates =  hel.PathLength(tRowRadius[ti]);
 629 //   tLength = (candidates.first > 0) ? candidates.first : candidates.second;
 630 //   if (IsNaN(tLength)){
 631 //     cout <<"tLength Init tmp NAN" << endl;
 632 //     cout <<"padrow number= "<<ti << "not reached" << endl;
 633 //     cout << "*** DO NOT ENTER THE LOOP***" << endl;
 634 //     tmpSect[ti]=-1;//sector
 635 //   }
 636 //   // end test
 637 //   while(ti<45 && !IsNaN(tLength)){
 638 //     candidates =  hel.PathLength(tRowRadius[ti]);
 639 //     tLength = (candidates.first > 0) ? candidates.first : candidates.second;
 640 //     if (IsNaN(tLength)){
 641 //       cout <<"tLength loop 1st NAN" << endl;
 642 //       cout <<"padrow number=  " << ti << " not reached" << endl;
 643 //       cout << "*** THIS IS AN ERROR SHOULDN'T  LOOP ***" << endl;
 644 //       tmpSect[ti]=-1;//sector
 645 //     }
 646 //     tPoint = hel.At(tLength);
 647 //     // Find which sector it is on
 648 //     TpcLocalTransform(tPoint,tmpSect[ti],tRow,tU,tPhi);
 649 //     if (IsNaN(tmpSect[ti])){
 650 //       cout <<"***ERROR tmpSect"<< endl;
 651 //     }
 652 //     if (IsNaN(tRow)){
 653 //       cout <<"***ERROR tRow"<< endl;
 654 //     }
 655 //     if (IsNaN(tU)){
 656 //       cout <<"***ERROR tU"<< endl;
 657 //     }
 658 //     if (IsNaN(tPhi)){
 659 //       cout <<"***ERROR tPhi"<< endl;
 660 //     }
 661 //     // calculate crossing plane
 662 //     tn.SetX(cos(tPhi));
 663 //     tn.SetY(sin(tPhi));
 664 //     tr.SetX(tRowRadius[ti]*cos(tPhi));
 665 //     tr.SetY(tRowRadius[ti]*sin(tPhi));
 666 //     // find crossing point
 667 //     tLength = hel.PathLength(tr,tn);
 668 //     if (IsNaN(tLength)){
 669 //       cout <<"tLength loop 2nd  NAN" << endl;
 670 //       cout <<"padrow number=  " << ti << " not reached" << endl;
 671 //       tmpSect[ti]=-2;//sector
 672 //     }
 673 //     tPoint = hel.At(tLength);
 674 //     tmpZ[ti] = tPoint.z();
 675 //     tOutOfBound = TpcLocalTransform(tPoint,tSect,tRow,tmpU[ti],tPhi);
 676 //     if (IsNaN(tSect)){
 677 //       cout <<"***ERROR tSect 2"<< endl;
 678 //     }
 679 //     if (IsNaN(tRow)){
 680 //       cout <<"***ERROR tRow 2"<< endl;
 681 //     }
 682 //     if (IsNaN(tmpU[ti])){
 683 //       cout <<"***ERROR tmpU[ti] 2"<< endl;
 684 //     }
 685 //     if (IsNaN(tPhi)){
 686 //       cout <<"***ERROR tPhi 2 "<< endl;
 687 //     }
 688 //     if(tOutOfBound || (tmpSect[ti] == tSect && tRow!=(ti+1))){
 689 //       tmpSect[ti]=-2;
 690 //       //       cout << "missed once"<< endl;
 691 //     }
 692 //     else{
 693 //       if(tmpSect[ti] != tSect){
 694 //      // Try again on the other sector
 695 //      tn.SetX(cos(tPhi));
 696 //      tn.SetY(sin(tPhi));
 697 //      tr.SetX(tRowRadius[ti]*cos(tPhi));
 698 //      tr.SetY(tRowRadius[ti]*sin(tPhi));
 699 //      // find crossing point
 700 //      tLength = hel.PathLength(tr,tn);
 701 //      tPoint = hel.At(tLength);
 702 //      if (IsNaN(tLength)){
 703 //        cout <<"tLength loop 3rd NAN" << endl;
 704 //        cout <<"padrow number=  "<< ti << " not reached" << endl;
 705 //        tmpSect[ti]=-1;//sector
 706 //      }
 707 //      tmpZ[ti] = tPoint.z();
 708 //      tmpSect[ti] = tSect;
 709 //      tOutOfBound = TpcLocalTransform(tPoint,tSect,tRow,tmpU[ti],tPhi);
 710 //      if (IsNaN(tSect)){
 711 //        cout <<"***ERROR tSect 3"<< endl;
 712 //      }
 713 //      if (IsNaN(tRow)){
 714 //        cout <<"***ERROR tRow 3"<< endl;
 715 //      }
 716 //      if (IsNaN(tmpU[ti])){
 717 //        cout <<"***ERROR tmpU[ti] 3"<< endl;
 718 //      }
 719 //      if (IsNaN(tPhi)){
 720 //        cout <<"***ERROR tPhi 3 "<< endl;
 721 //      }
 722 //      if(tOutOfBound || tSect!= tmpSect[ti] || tRow!=(ti+1)){
 723 //        tmpSect[ti]=-1;
 724 //      }
 725 //       }
 726 //     }
 727 //     if (IsNaN(tmpSect[ti])){
 728 //       cout << "*******************ERROR***************************" << endl;
 729 //       cout <<"AliFemtoParticle--Fctn tmpSect=" << tmpSect[ti] << endl;
 730 //       cout << "*******************ERROR***************************" << endl;
 731 //     }
 732 //     if (IsNaN(tmpU[ti])){
 733 //       cout << "*******************ERROR***************************" << endl;
 734 //       cout <<"AliFemtoParticle--Fctn tmpU=" << tmpU[ti] << endl;
 735 //       cout << "*******************ERROR***************************" << endl;
 736 //     }
 737 //     if (IsNaN(tmpZ[ti])){
 738 //       cout << "*******************ERROR***************************" << endl;
 739 //       cout <<"AliFemtoParticle--Fctn tmpZ=" << tmpZ[ti] << endl;
 740 //       cout << "*******************ERROR***************************" << endl;
 741 //     }
 742 //     // If padrow ti not reached all other beyond are not reached
 743 //     // in this case Set sector to -1
 744 //     if (tmpSect[ti]==-1){
 745 //       for (int tj=ti; tj<45;tj++){
 746 //      tmpSect[tj] = -1;
 747 //      ti=45;
 748 //       }
 749 //     }
 750 //     ti++;
 751 //     if (ti<45){
 752 //       candidates =  hel.PathLength(tRowRadius[ti]);
 753 //       tLength = (candidates.first > 0) ? candidates.first : candidates.second;}
 754 //   }
 755 // }
 756 //_____________________
 757 const AliFemtoThreeVector& AliFemtoParticle::TpcV0PosExitPoint() const{
 758   return fTpcV0PosExitPoint;
 759 }
 760 //_____________________
 761 const AliFemtoThreeVector& AliFemtoParticle::TpcV0PosEntrancePoint() const{
 762   return fTpcV0PosEntrancePoint;
 763 }
 764 //______________________
 765 const AliFemtoThreeVector& AliFemtoParticle::TpcV0NegExitPoint() const{
 766   return fTpcV0NegExitPoint;
 767 }
 768 //_____________________
 769 const AliFemtoThreeVector& AliFemtoParticle::TpcV0NegEntrancePoint() const{
 770   return fTpcV0NegEntrancePoint;
 771 }
 772 //______________________