PWG2/FEMTOSCOPY/AliFemto/Infrastructure/AliFemtoParticle.cxx

   1 /***************************************************************************
   2  *
   3  * $Id$
   4  *
   5  * Author: Mike Lisa, Ohio State, lisa@mps.ohio-state.edu
   6  ***************************************************************************
   7  *
   8  * Description: part of STAR HBT Framework: AliFemtoMaker package
   9  *   Particle objects are part of the PicoEvent, which is what is
  10  *   stored in the EventMixingBuffers
  11  *   A Track object gets converted to a Particle object if it
  12  *   passes the ParticleCut of an Analysis
  13  *
  14  ***************************************************************************
  15  *
  16  * $Log$
  17  * Revision 1.1.1.1  2007/04/25 15:38:41  panos
  18  * Importing the HBT code dir
  19  *
  20  * Revision 1.1.1.1  2007/03/07 10:14:49  mchojnacki
  21  * First version on CVS
  22  *
  23  * Revision 1.22  2003/09/02 17:58:32  perev
  24  * gcc 3.2 updates + WarnOff
  25  *
  26  * Revision 1.21  2003/05/07 15:30:43  magestro
  27  * Fixed bug related to finding merged hits (commit for Fabrice)
  28  *
  29  * Revision 1.20  2003/01/14 09:41:26  renault
  30  * changes on average separation calculation, hit shared finder and memory optimisation
  31  * for Z,U and Sectors variables.
  32  *
  33  * Revision 1.19  2002/12/12 17:01:49  kisiel
  34  * Hidden Information handling and purity calculation
  35  *
  36  * Revision 1.18  2002/11/19 23:36:00  renault
  37  * Enable calculation of exit/entrance separation for V0 daughters
  38  *
  39  * Revision 1.17  2001/12/14 23:11:30  fretiere
  40  * Add class HitMergingCut. Add class fabricesPairCut = HitMerginCut + pair purity cuts. Add TpcLocalTransform function which convert to local tpc coord (not pretty). Modify AliFemtoTrack, AliFemtoParticle, AliFemtoHiddenInfo, AliFemtoPair to handle the hit information and cope with my code
  41  *
  42  * Revision 1.16  2001/05/25 23:23:59  lisa
  43  * Added in AliFemtoKink stuff
  44  *
  45  * Revision 1.15  2001/04/03 21:04:36  kisiel
  46  *
  47  *
  48  *   Changes needed to make the Theoretical code
  49  *   work. The main code is the ThCorrFctn directory.
  50  *   The most visible change is the addition of the
  51  *   HiddenInfo to AliFemtoPair.
  52  *
  53  * Revision 1.14  2000/10/05 23:09:05  lisa
  54  * Added kT-dependent radii to mixed-event simulator AND implemented AverageSeparation Cut and CorrFctn
  55  *
  56  * Revision 1.13  2000/08/31 22:31:31  laue
  57  * AliFemtoAnalysis: output changed (a little bit less)
  58  * AliFemtoEvent: new version, members for reference mult added
  59  * AliFemtoIOBinary: new IO for new AliFemtoEvent version
  60  * AliFemtoTypes: TTree typedef to AliFemtoTTree added
  61  * AliFemtoVertexAnalysis: overflow and underflow added
  62  *
  63  * Revision 1.12  2000/07/23 13:52:56  laue
  64  * NominalExitPoint set to (-9999.,-9999.-9999.) if helix.at()
  65  * returns nan (not a number).
  66  *
  67  * Revision 1.11  2000/07/19 17:18:48  laue
  68  * Calculation of Entrance and Exit point added in AliFemtoParticle constructor
  69  *
  70  * Revision 1.10  2000/07/17 20:03:17  lisa
  71  * Implemented tools for addressing and assessing trackmerging
  72  *
  73  * Revision 1.9  2000/07/16 21:38:23  laue
  74  * AliFemtoCoulomb.cxx AliFemtoSectoredAnalysis.cxx : updated for standalone version
  75  * AliFemtoV0.cc AliFemtoV0.h : some cast to prevent compiling warnings
  76  * AliFemtoParticle.cc AliFemtoParticle.h : pointers fTrack,fV0 initialized to 0
  77  * AliFemtoIOBinary.cc : some printouts in #ifdef STHBTDEBUG
  78  * AliFemtoEvent.cc : B-Field set to 0.25Tesla, we have to think about a better
  79  *                 solution
  80  *
  81  * Revision 1.8  2000/05/03 17:44:43  laue
  82  * AliFemtoEvent, AliFemtoTrack & AliFemtoV0 declared friend to AliFemtoIOBinary
  83  * AliFemtoParticle updated for V0 pos,neg track Id
  84  *
  85  * Revision 1.7  2000/04/03 16:21:51  laue
  86  * some include files changed
  87  * Multi track cut added
  88  *
  89  * Revision 1.6  1999/12/11 15:58:29  lisa
  90  * Add vertex decay position datum and accessor to AliFemtoParticle to allow pairwise cuts on seperation of V0s
  91  *
  92  * Revision 1.5  1999/09/17 22:38:02  lisa
  93  * first full integration of V0s into AliFemto framework
  94  *
  95  * Revision 1.4  1999/09/01 19:04:53  lisa
  96  * update Particle class AND add parity cf and Randys Coulomb correction
  97  *
  98  * Revision 1.3  1999/07/06 22:33:23  lisa
  99  * Adjusted all to work in pro and new - dev itself is broken
 100  *
 101  * Revision 1.2  1999/06/29 17:50:27  fisyak
 102  * formal changes to account new StEvent, does not complie yet
 103  *
 104  * Revision 1.1.1.1  1999/06/29 16:02:57  lisa
 105  * Installation of AliFemtoMaker
 106  *
 107  **************************************************************************/
 108
 109 #include "Infrastructure/AliFemtoParticle.h"
 110 //#include "math_constants.h"
 111 #ifdef __CC5__
 112   #include <math.h>
 113 #else
 114   #include <cmath>
 115 #endif
 116
 117
 118 #include "TMath.h"
 119 using namespace TMath;
 120
 121 double AliFemtoParticle::fPrimPimPar0= 9.05632e-01;
 122 double AliFemtoParticle::fPrimPimPar1= -2.26737e-01;
 123 double AliFemtoParticle::fPrimPimPar2= -1.03922e-01;
 124 double AliFemtoParticle::fPrimPipPar0= 9.09616e-01;
 125 double AliFemtoParticle::fPrimPipPar1= -9.00511e-02;
 126 double AliFemtoParticle::fPrimPipPar2= -6.02940e-02;
 127 double AliFemtoParticle::fPrimPmPar0= 0.;
 128 double AliFemtoParticle::fPrimPmPar1= 0.;
 129 double AliFemtoParticle::fPrimPmPar2= 0.;
 130 double AliFemtoParticle::fPrimPpPar0= 0.;
 131 double AliFemtoParticle::fPrimPpPar1= 0.;
 132 double AliFemtoParticle::fPrimPpPar2= 0.;
 133
 134 int TpcLocalTransform(AliFmThreeVectorD& xgl,
 135                       int& iSector,
 136                       int& iPadrow,
 137                       float& xlocal,
 138                       double& ttPhi);
 139
 140
 141 //_____________________
 142 AliFemtoParticle::AliFemtoParticle() :
 143   fTpcV0NegPosSample(0),
 144   fV0NegZ(0),
 145   fV0NegU(0),
 146   fV0NegSect(0),
 147   fTrack(0), fV0(0), fKink(0), fXi(0),
 148   fFourMomentum(0),
 149   fHelix(),
 150   fNominalTpcExitPoint(0),
 151   fNominalTpcEntrancePoint(0),
 152   fHiddenInfo(0),
 153   fPrimaryVertex(0),
 154   fSecondaryVertex(0),
 155   fHelixV0Pos(),
 156   fTpcV0PosEntrancePoint(0),
 157   fTpcV0PosExitPoint(0),
 158   fHelixV0Neg(),
 159   fTpcV0NegEntrancePoint(0),
 160   fTpcV0NegExitPoint(0)
 161 {
 162   /* no-op for default */
 163   //  cout << "Created particle " << this << endl;
 164 }
 165 //_____________________
 166 AliFemtoParticle::AliFemtoParticle(const AliFemtoParticle& aParticle):
 167   fTpcV0NegPosSample(0),
 168   fV0NegZ(0),
 169   fV0NegU(0),
 170   fV0NegSect(0),
 171   fTrack(0), fV0(0), fKink(0), fXi(0),
 172   fFourMomentum(0),
 173   fHelix(),
 174   fNominalTpcExitPoint(0),
 175   fNominalTpcEntrancePoint(0),
 176   fHiddenInfo(0),
 177   fPrimaryVertex(0),
 178   fSecondaryVertex(0),
 179   fHelixV0Pos(),
 180   fTpcV0PosEntrancePoint(0),
 181   fTpcV0PosExitPoint(0),
 182   fHelixV0Neg(),
 183   fTpcV0NegEntrancePoint(0),
 184   fTpcV0NegExitPoint(0)
 185 {
 186   if (aParticle.fTrack)
 187     fTrack = new AliFemtoTrack(*aParticle.fTrack);
 188   if (aParticle.fV0)
 189     fV0    = new AliFemtoV0(*aParticle.fV0);
 190   if (aParticle.fKink)
 191     fKink  = new AliFemtoKink(*aParticle.fKink);
 192   if (aParticle.fXi)
 193     fXi    = new AliFemtoXi(*aParticle.fXi);
 194   fFourMomentum = aParticle.fFourMomentum;
 195   fHelix = aParticle.fHelix;
 196
 197   for (int iter=0; iter<11; iter++)
 198     fNominalPosSample[iter] = aParticle.fNominalPosSample[iter];
 199
 200   if (aParticle.fTpcV0NegPosSample) {
 201     fTpcV0NegPosSample = (AliFemtoThreeVector *) malloc(sizeof(AliFemtoThreeVector) * 11);
 202     for (int iter=0; iter<11; iter++)
 203       fTpcV0NegPosSample[iter] = aParticle.fTpcV0NegPosSample[iter];
 204   }
 205
 206   if (aParticle.fV0NegZ) {
 207     fV0NegZ = (float *) malloc(sizeof(float) * 45);
 208     for (int iter=0; iter<11; iter++)
 209       fV0NegZ[iter] = aParticle.fV0NegZ[iter];
 210   }
 211   if (aParticle.fV0NegU) {
 212     fV0NegU = (float *) malloc(sizeof(float) * 45);
 213     for (int iter=0; iter<11; iter++)
 214       fV0NegU[iter] = aParticle.fV0NegU[iter];
 215   }
 216   if (aParticle.fV0NegSect) {
 217     fV0NegSect = (int *) malloc(sizeof(int) * 45);
 218     for (int iter=0; iter<11; iter++)
 219       fV0NegSect[iter] = aParticle.fV0NegSect[iter];
 220   }
 221
 222   fPrimaryVertex = aParticle.fPrimaryVertex;
 223   fSecondaryVertex = aParticle.fSecondaryVertex;
 224   CalculatePurity();
 225   if(aParticle.fHiddenInfo){
 226     fHiddenInfo= aParticle.HiddenInfo()->clone();
 227   }
 228
 229   fNominalTpcEntrancePoint = aParticle.fNominalTpcEntrancePoint;
 230   fNominalTpcExitPoint     = aParticle.fNominalTpcExitPoint;
 231
 232   for (int iter=0; iter<6; iter++)
 233     fPurity[iter] = aParticle.fPurity[iter];
 234
 235   fHelixV0Pos = aParticle.fHelixV0Pos;
 236   fTpcV0PosEntrancePoint = aParticle.fTpcV0PosEntrancePoint;
 237   fTpcV0PosExitPoint     = aParticle.fTpcV0PosExitPoint;
 238   fHelixV0Neg = aParticle.fHelixV0Neg;
 239   fTpcV0NegEntrancePoint = aParticle.fTpcV0NegEntrancePoint;
 240   fTpcV0NegExitPoint     = aParticle.fTpcV0NegExitPoint;
 241 }
 242 //_____________________
 243 AliFemtoParticle::~AliFemtoParticle(){
 244   //  cout << "Issuing delete for AliFemtoParticle." << endl;
 245
 246   if (fTrack) delete fTrack;
 247   if (fV0) {
 248     delete[] fTpcV0NegPosSample;
 249     delete[] fV0NegZ;
 250     delete[] fV0NegU;
 251     delete[] fV0NegSect;
 252     delete fV0;
 253   }
 254   if (fKink) delete fKink;
 255   //  cout << "Trying to delete HiddenInfo: " << fHiddenInfo << endl;
 256   if (fHiddenInfo)
 257     {
 258       //      cout << "Deleting HiddenInfo." << endl;
 259       delete fHiddenInfo;
 260     }
 261   //  cout << "Deleted particle " << this << endl;
 262 }
 263 //_____________________
 264 AliFemtoParticle::AliFemtoParticle(const AliFemtoTrack* const hbtTrack,const double& mass) :
 265   fTpcV0NegPosSample(0),
 266   fV0NegZ(0),
 267   fV0NegU(0),
 268   fV0NegSect(0),
 269   fTrack(0), fV0(0), fKink(0), fXi(0),
 270   fFourMomentum(0),
 271   fHelix(),
 272   fNominalTpcExitPoint(0),
 273   fNominalTpcEntrancePoint(0),
 274   fHiddenInfo(0),
 275   fPrimaryVertex(0),
 276   fSecondaryVertex(0),
 277   fHelixV0Pos(),
 278   fTpcV0PosEntrancePoint(0),
 279   fTpcV0PosExitPoint(0),
 280   fHelixV0Neg(),
 281   fTpcV0NegEntrancePoint(0),
 282   fTpcV0NegExitPoint(0)
 283 {
 284
 285
 286   // I know there is a better way to do this...
 287   fTrack = new AliFemtoTrack(*hbtTrack);
 288   AliFemtoThreeVector temp = hbtTrack->P();
 289   fFourMomentum.setVect(temp);
 290   double ener = ::sqrt(temp.mag2()+mass*mass);
 291   fFourMomentum.setE(ener);
 292 //  fMap[0] = hbtTrack->TopologyMap(0);
 293  // fMap[1] = hbtTrack->TopologyMap(1);
 294  // fNhits = hbtTrack->NHits();
 295   fHelix = hbtTrack->Helix();
 296   //CalculateNominalTpcExitAndEntrancePoints();
 297
 298
 299   fPrimaryVertex.setX(0.);
 300   fPrimaryVertex.setY(0.);
 301   fPrimaryVertex.setZ(0.);
 302   fSecondaryVertex.setX(0.);
 303   fSecondaryVertex.setY(0.);
 304   fSecondaryVertex.setZ(0.);
 305   /* TO JA ODZNACZYLEM NIE WIEM DLACZEGO
 306   CalculateTpcExitAndEntrancePoints(&fHelix,&fPrimaryVertex,
 307                                     &fSecondaryVertex,
 308                                     &fNominalTpcEntrancePoint,
 309                                     &fNominalTpcExitPoint,
 310                                     &mNominalPosSample[0],
 311                                     &fZ[0],
 312                                     &fU[0],
 313                                     &fSect[0]);
 314   */
 315   CalculatePurity();
 316   // ***
 317   fHiddenInfo= 0;
 318   if(hbtTrack->ValidHiddenInfo()){
 319     fHiddenInfo= hbtTrack->getHiddenInfo()->getParticleHiddenInfo()->clone();
 320   }
 321   // ***
 322   //  cout << "Created particle " << this << endl;
 323
 324 }
 325 //_____________________
 326 AliFemtoParticle::AliFemtoParticle(const AliFemtoV0* const hbtV0,const double& mass) :
 327   fTpcV0NegPosSample(0),
 328   fV0NegZ(0),
 329   fV0NegU(0),
 330   fV0NegSect(0),
 331   fTrack(0), fV0(0), fKink(0),  fXi(0),
 332   fFourMomentum(0),
 333   fHelix(),
 334   fNominalTpcExitPoint(0),
 335   fNominalTpcEntrancePoint(0),
 336   fHiddenInfo(0),
 337   fPrimaryVertex(0),
 338   fSecondaryVertex(0),
 339   fHelixV0Pos(),
 340   fTpcV0PosEntrancePoint(0),
 341   fTpcV0PosExitPoint(0),
 342   fHelixV0Neg(),
 343   fTpcV0NegEntrancePoint(0),
 344   fTpcV0NegExitPoint(0)
 345 {
 346   fV0 = new AliFemtoV0(*hbtV0);
 347  //fMap[0]= 0;
 348   //fMap[1]= 0;
 349   // I know there is a better way to do this...
 350   AliFemtoThreeVector temp = hbtV0->momV0();
 351   fFourMomentum.setVect(temp);
 352   double ener = ::sqrt(temp.mag2()+mass*mass);
 353   fFourMomentum.setE(ener);
 354   // Calculating TpcEntrancePoint for Positive V0 daugther
 355   fPrimaryVertex = hbtV0->primaryVertex();
 356   fSecondaryVertex = hbtV0->decayVertexV0();
 357   fHelixV0Pos = hbtV0->HelixPos();
 358
 359   fTpcV0NegPosSample = new AliFemtoThreeVector[45];//for V0Neg
 360   fV0NegZ = new float[45];//for V0Neg
 361   fV0NegU = new float[45];//for V0Neg
 362   fV0NegSect = new int[45];//for V0Neg
 363   CalculateTpcExitAndEntrancePoints(&fHelixV0Pos,&fPrimaryVertex,
 364                                     &fSecondaryVertex,
 365                                     &fTpcV0PosEntrancePoint,
 366                                     &fTpcV0PosExitPoint,
 367                                     &fNominalPosSample[0],
 368                                     &fZ[0],
 369                                     &fU[0],&fSect[0]);
 370   fHelixV0Neg = hbtV0->HelixNeg();
 371
 372   CalculateTpcExitAndEntrancePoints(&fHelixV0Neg,
 373                                     &fPrimaryVertex,
 374                                     &fSecondaryVertex,
 375                                     &fTpcV0NegEntrancePoint,
 376                                     &fTpcV0NegExitPoint,
 377                                     &fTpcV0NegPosSample[0],
 378                                     &fV0NegZ[0],
 379                                     &fV0NegU[0],&fV0NegSect[0]);
 380
 381   // ***
 382   fHiddenInfo= 0;
 383   if(hbtV0->ValidHiddenInfo()){
 384     fHiddenInfo= hbtV0->getHiddenInfo()->clone();
 385   }
 386   // ***
 387 }
 388 //_____________________
 389 AliFemtoParticle::AliFemtoParticle(const AliFemtoKink* const hbtKink,const double& mass) :
 390   fTpcV0NegPosSample(0),
 391   fV0NegZ(0),
 392   fV0NegU(0),
 393   fV0NegSect(0),
 394   fTrack(0), fV0(0), fKink(0), fXi(0),
 395   fFourMomentum(0),
 396   fHelix(),
 397   fNominalTpcExitPoint(0),
 398   fNominalTpcEntrancePoint(0),
 399   fHiddenInfo(0),
 400   fPrimaryVertex(0),
 401   fSecondaryVertex(0),
 402   fHelixV0Pos(),
 403   fTpcV0PosEntrancePoint(0),
 404   fTpcV0PosExitPoint(0),
 405   fHelixV0Neg(),
 406   fTpcV0NegEntrancePoint(0),
 407   fTpcV0NegExitPoint(0)
 408 {
 409   fKink = new AliFemtoKink(*hbtKink);
 410  // fMap[0]= 0;
 411   //fMap[1]= 0;
 412   // I know there is a better way to do this...
 413   AliFemtoThreeVector temp = hbtKink->Parent().P();
 414   fFourMomentum.setVect(temp);
 415   double ener = ::sqrt(temp.mag2()+mass*mass);
 416   fFourMomentum.setE(ener);
 417 }
 418
 419 //_____________________
 420 AliFemtoParticle::AliFemtoParticle(const AliFemtoXi* const hbtXi, const double& mass) :
 421   fTpcV0NegPosSample(0),
 422   fV0NegZ(0),
 423   fV0NegU(0),
 424   fV0NegSect(0),
 425   fTrack(0), fV0(0), fKink(0), fXi(0),
 426   fFourMomentum(0),
 427   fHelix(),
 428   fNominalTpcExitPoint(0),
 429   fNominalTpcEntrancePoint(0),
 430   fHiddenInfo(0),
 431   fPrimaryVertex(0),
 432   fSecondaryVertex(0),
 433   fHelixV0Pos(),
 434   fTpcV0PosEntrancePoint(0),
 435   fTpcV0PosExitPoint(0),
 436   fHelixV0Neg(),
 437   fTpcV0NegEntrancePoint(0),
 438   fTpcV0NegExitPoint(0)
 439 {
 440   fXi = new AliFemtoXi(*hbtXi);
 441  // fMap[0]= 0;
 442   //fMap[1]= 0;
 443   AliFemtoThreeVector temp;// = hbtXi->mMofXi;
 444   fFourMomentum.setVect(temp);
 445   double ener = ::sqrt(temp.mag2()+mass*mass);
 446   fFourMomentum.setE(ener);
 447   fHiddenInfo = 0;
 448 }
 449 //_____________________
 450 AliFemtoParticle& AliFemtoParticle::operator=(const AliFemtoParticle& aParticle)
 451 {
 452   if (this == &aParticle)
 453     return *this;
 454
 455   if (aParticle.fTrack)
 456     fTrack = new AliFemtoTrack(*aParticle.fTrack);
 457   if (aParticle.fV0)
 458     fV0    = new AliFemtoV0(*aParticle.fV0);
 459   if (aParticle.fKink)
 460     fKink  = new AliFemtoKink(*aParticle.fKink);
 461   if (aParticle.fXi)
 462     fXi    = new AliFemtoXi(*aParticle.fXi);
 463   fFourMomentum = aParticle.fFourMomentum;
 464   fHelix = aParticle.fHelix;
 465
 466   for (int iter=0; iter<11; iter++)
 467     fNominalPosSample[iter] = aParticle.fNominalPosSample[iter];
 468
 469   if (fTpcV0NegPosSample) delete fTpcV0NegPosSample;
 470   if (aParticle.fTpcV0NegPosSample) {
 471     fTpcV0NegPosSample = (AliFemtoThreeVector *) malloc(sizeof(AliFemtoThreeVector) * 11);
 472     for (int iter=0; iter<11; iter++)
 473       fTpcV0NegPosSample[iter] = aParticle.fTpcV0NegPosSample[iter];
 474   }
 475
 476   if (fV0NegZ) delete fV0NegZ;
 477   if (aParticle.fV0NegZ) {
 478     fV0NegZ = (float *) malloc(sizeof(float) * 45);
 479     for (int iter=0; iter<11; iter++)
 480       fV0NegZ[iter] = aParticle.fV0NegZ[iter];
 481   }
 482   if (fV0NegU) delete fV0NegU;
 483   if (aParticle.fV0NegU) {
 484     fV0NegU = (float *) malloc(sizeof(float) * 45);
 485     for (int iter=0; iter<11; iter++)
 486       fV0NegU[iter] = aParticle.fV0NegU[iter];
 487   }
 488   if (fV0NegSect) delete fV0NegSect;
 489   if (aParticle.fV0NegSect) {
 490     fV0NegSect = (int *) malloc(sizeof(int) * 45);
 491     for (int iter=0; iter<11; iter++)
 492       fV0NegSect[iter] = aParticle.fV0NegSect[iter];
 493   }
 494
 495   fPrimaryVertex = aParticle.fPrimaryVertex;
 496   fSecondaryVertex = aParticle.fSecondaryVertex;
 497   CalculatePurity();
 498   if(aParticle.fHiddenInfo){
 499     fHiddenInfo= aParticle.fHiddenInfo->clone();
 500   }
 501
 502   fNominalTpcEntrancePoint = aParticle.fNominalTpcEntrancePoint;
 503   fNominalTpcExitPoint     = aParticle.fNominalTpcExitPoint;
 504
 505   if (fHiddenInfo) delete fHiddenInfo;
 506   if (aParticle.fHiddenInfo)
 507     fHiddenInfo = aParticle.HiddenInfo()->clone();
 508
 509   for (int iter=0; iter<6; iter++)
 510     fPurity[iter] = aParticle.fPurity[iter];
 511
 512   fHelixV0Pos = aParticle.fHelixV0Pos;
 513   fTpcV0PosEntrancePoint = aParticle.fTpcV0PosEntrancePoint;
 514   fTpcV0PosExitPoint     = aParticle.fTpcV0PosExitPoint;
 515   fHelixV0Neg = aParticle.fHelixV0Neg;
 516   fTpcV0NegEntrancePoint = aParticle.fTpcV0NegEntrancePoint;
 517   fTpcV0NegExitPoint     = aParticle.fTpcV0NegExitPoint;
 518
 519   return *this;
 520 }
 521 //_____________________
 522 const AliFemtoThreeVector& AliFemtoParticle::NominalTpcExitPoint() const{
 523   // in future, may want to calculate this "on demand" only, sot this routine may get more sophisticated
 524   // for now, we calculate Exit and Entrance points upon instantiation
 525   return fNominalTpcExitPoint;
 526 }
 527 //_____________________
 528 const AliFemtoThreeVector& AliFemtoParticle::NominalTpcEntrancePoint() const{
 529   // in future, may want to calculate this "on demand" only, sot this routine may get more sophisticated
 530   // for now, we calculate Exit and Entrance points upon instantiation
 531   return fNominalTpcEntrancePoint;
 532 }
 533 //_____________________
 534 void AliFemtoParticle::CalculatePurity(){
 535   double tPt = fFourMomentum.perp();
 536   // pi -
 537   fPurity[0] = fPrimPimPar0*(1.-exp((tPt-fPrimPimPar1)/fPrimPimPar2));
 538   fPurity[0] *= fTrack->PidProbPion();
 539   // pi+
 540   fPurity[1] = fPrimPipPar0*(1.-exp((tPt-fPrimPipPar1)/fPrimPipPar2));
 541   fPurity[1] *= fTrack->PidProbPion();
 542   // K-
 543   fPurity[2] = fTrack->PidProbKaon();
 544   // K+
 545   fPurity[3] = fTrack->PidProbKaon();
 546   // pbar
 547   fPurity[4] = fTrack->PidProbProton();
 548   // p
 549   fPurity[5] = fTrack->PidProbProton();
 550 }
 551
 552 double AliFemtoParticle::GetPionPurity()
 553 {
 554   if (fTrack->Charge()>0)
 555     return fPurity[1];
 556   else
 557     return fPurity[0];
 558 }
 559 double AliFemtoParticle::GetKaonPurity()
 560 {
 561   if (fTrack->Charge()>0)
 562     return fPurity[3];
 563   else
 564     return fPurity[2];
 565 }
 566 double AliFemtoParticle::GetProtonPurity()
 567 {
 568   if (fTrack->Charge()>0)
 569     return fPurity[5];
 570   else
 571     return fPurity[4];
 572 }
 573
 574 void AliFemtoParticle::CalculateTpcExitAndEntrancePoints(AliFmPhysicalHelixD* tHelix,
 575                                                        AliFemtoThreeVector*  PrimVert,
 576                                                        AliFemtoThreeVector*  SecVert,
 577                                                        AliFemtoThreeVector* tmpTpcEntrancePoint,
 578                                                        AliFemtoThreeVector* tmpTpcExitPoint,
 579                                                        AliFemtoThreeVector* tmpPosSample,
 580                                                        float* tmpZ,
 581                                                        float* tmpU,
 582                                                        int* tmpSect){
 583   // this calculates the exit point of a secondary track,
 584   // either through the endcap or through the Outer Field Cage
 585   // We assume the track to start at tHelix.origin-PrimaryVertex
 586   // it also calculates the entrance point of the secondary track,
 587   // which is the point at which it crosses the
 588   // inner field cage
 589   //  static AliFemtoThreeVector ZeroVec(0.,0.,0.);
 590   AliFemtoThreeVector ZeroVec(0.,0.,0.);
 591 //   ZeroVec.setX(tHelix->origin().x()-PrimVert->x());
 592 //   ZeroVec.setY(tHelix->origin().y()-PrimVert->y());
 593 //   ZeroVec.setZ(tHelix->origin().z()-PrimVert->z());
 594   ZeroVec.setX(SecVert->x()-PrimVert->x());
 595   ZeroVec.setY(SecVert->y()-PrimVert->y());
 596   ZeroVec.setZ(SecVert->z()-PrimVert->z());
 597   double dip, curv, phase;
 598   int h;
 599   curv = tHelix->curvature();
 600   dip  = tHelix->dipAngle();
 601   phase= tHelix->phase();
 602   h    = tHelix->h();
 603
 604   AliFmHelixD hel(curv,dip,phase,ZeroVec,h);
 605
 606   pairD candidates;
 607   double sideLength;  // this is how much length to go to leave through sides of TPC
 608   double endLength;  // this is how much length to go to leave through endcap of TPC
 609   // figure out how far to go to leave through side...
 610   candidates = hel.pathLength(200.0);  // bugfix MAL jul00 - 200cm NOT 2cm
 611   sideLength = (candidates.first > 0) ? candidates.first : candidates.second;
 612
 613   static AliFemtoThreeVector WestEnd(0.,0.,200.);  // bugfix MAL jul00 - 200cm NOT 2cm
 614   static AliFemtoThreeVector EastEnd(0.,0.,-200.); // bugfix MAL jul00 - 200cm NOT 2cm
 615   static AliFemtoThreeVector EndCapNormal(0.,0.,1.0);
 616
 617   endLength = hel.pathLength(WestEnd,EndCapNormal);
 618   if (endLength < 0.0) endLength = hel.pathLength(EastEnd,EndCapNormal);
 619
 620   if (endLength < 0.0) cout <<
 621                          "AliFemtoParticle::CalculateTpcExitAndEntrancePoints(): "
 622                             << "Hey -- I cannot find an exit point out endcaps" << endl;
 623   // OK, firstExitLength will be the shortest way out of the detector...
 624   double firstExitLength = (endLength < sideLength) ? endLength : sideLength;
 625   // now then, let's return the POSITION at which particle leaves TPC...
 626   *tmpTpcExitPoint = hel.at(firstExitLength);
 627   // Finally, calculate the position at which the track crosses the inner field cage
 628   candidates = hel.pathLength(50.0);  // bugfix MAL jul00 - 200cm NOT 2cm
 629
 630   sideLength = (candidates.first > 0) ? candidates.first : candidates.second;
 631   //  cout << "sideLength 2 ="<<sideLength << endl;
 632   *tmpTpcEntrancePoint = hel.at(sideLength);
 633   // This is the secure way !
 634   if (IsNaN(tmpTpcEntrancePoint->x()) ||
 635       IsNaN(tmpTpcEntrancePoint->y()) ||
 636       IsNaN(tmpTpcEntrancePoint->z()) ){
 637     cout << "tmpTpcEntrancePoint NAN"<< endl;
 638     cout << "tmpNominalTpcEntrancePoint = " <<tmpTpcEntrancePoint<< endl;
 639     tmpTpcEntrancePoint->setX(-9999.);
 640     tmpTpcEntrancePoint->setY(-9999.);
 641     tmpTpcEntrancePoint->setZ(-9999.);
 642   }
 643
 644   if (IsNaN(tmpTpcExitPoint->x()) ||
 645       IsNaN(tmpTpcExitPoint->y()) ||
 646       IsNaN(tmpTpcExitPoint->z()) ) {
 647 //     cout << "tmpTpcExitPoint NAN set at (-9999,-9999,-9999)"<< endl;
 648 //     cout << "tmpTpcExitPoint X= " <<tmpTpcExitPoint->x()<< endl;
 649 //     cout << "tmpTpcExitPoint Y= " <<tmpTpcExitPoint->y()<< endl;
 650 //     cout << "tmpTpcExitPoint Z= " <<tmpTpcExitPoint->z()<< endl;
 651     tmpTpcExitPoint->setX(-9999.);
 652     tmpTpcExitPoint->setY(-9999.);
 653     tmpTpcExitPoint->setZ(-9999.);
 654   }
 655
 656
 657 //   if (IsNaN(tmpTpcExitPoint->x())) *tmpTpcExitPoint = AliFemtoThreeVector(-9999.,-9999.,-9999);
 658 //   if (IsNaN(tmpTpcEntrancetPoint->x())) *tmpTpcEntrancePoint = AliFemtoThreeVector(-9999.,-9999.,-9999);
 659   //  cout << "tmpTpcEntrancePoint"<<*tmpTpcEntrancePoint << endl;
 660
 661   // 03Oct00 - mal.  OK, let's try something a little more
 662   // along the lines of NA49 and E895 strategy.
 663   // calculate the "nominal" position at N radii (say N=11)
 664   // within the TPC, and for a pair cut
 665   // use the average separation of these N
 666   int irad = 0;
 667   candidates = hel.pathLength(50.0);
 668   sideLength = (candidates.first > 0) ? candidates.first : candidates.second;
 669   while (irad<11 && !IsNaN(sideLength)){
 670     float radius = 50.0 + irad*15.0;
 671     candidates = hel.pathLength(radius);
 672     sideLength = (candidates.first > 0) ? candidates.first : candidates.second;
 673     tmpPosSample[irad] = hel.at(sideLength);
 674     if(IsNaN(tmpPosSample[irad].x()) ||
 675        IsNaN(tmpPosSample[irad].y()) ||
 676        IsNaN(tmpPosSample[irad].z())
 677        ){
 678       cout << "tmpPosSample for radius=" << radius << " NAN"<< endl;
 679       cout << "tmpPosSample=(" <<tmpPosSample[irad]<<")"<< endl;
 680       tmpPosSample[irad] =  AliFemtoThreeVector(-9999.,-9999.,-9999);
 681     }
 682     irad++;
 683     if (irad<11){
 684       float radius = 50.0 + irad*15.0;
 685       candidates = hel.pathLength(radius);
 686       sideLength = (candidates.first > 0) ? candidates.first : candidates.second;
 687     }
 688    }
 689    for (int i = irad; i<11; i++)
 690      {
 691        tmpPosSample[i] =  AliFemtoThreeVector(-9999.,-9999.,-9999);
 692      }
 693
 694   static float tRowRadius[45] = {60,64.8,69.6,74.4,79.2,84,88.8,93.6,98.8,
 695                                  104,109.2,114.4,119.6,127.195,129.195,131.195,
 696                                  133.195,135.195,137.195,139.195,141.195,
 697                                  143.195,145.195,147.195,149.195,151.195,
 698                                  153.195,155.195,157.195,159.195,161.195,
 699                                  163.195,165.195,167.195,169.195,171.195,
 700                                  173.195,175.195,177.195,179.195,181.195,
 701                                  183.195,185.195,187.195,189.195};
 702   int tRow,tSect,tOutOfBound;
 703   double tLength,tPhi;
 704   float tU;
 705   AliFemtoThreeVector tPoint;
 706   AliFmThreeVectorD tn(0,0,0);
 707   AliFmThreeVectorD tr(0,0,0);
 708   int ti =0;
 709   // test to enter the loop
 710   candidates =  hel.pathLength(tRowRadius[ti]);
 711   tLength = (candidates.first > 0) ? candidates.first : candidates.second;
 712   if (IsNaN(tLength)){
 713     cout <<"tLength Init tmp NAN" << endl;
 714     cout <<"padrow number= "<<ti << "not reached" << endl;
 715     cout << "*** DO NOT ENTER THE LOOP***" << endl;
 716     tmpSect[ti]=-1;//sector
 717   }
 718   // end test
 719   while(ti<45 && !IsNaN(tLength)){
 720     candidates =  hel.pathLength(tRowRadius[ti]);
 721     tLength = (candidates.first > 0) ? candidates.first : candidates.second;
 722     if (IsNaN(tLength)){
 723       cout <<"tLength loop 1st NAN" << endl;
 724       cout <<"padrow number=  " << ti << " not reached" << endl;
 725       cout << "*** THIS IS AN ERROR SHOULDN'T  LOOP ***" << endl;
 726       tmpSect[ti]=-1;//sector
 727     }
 728     tPoint = hel.at(tLength);
 729     // Find which sector it is on
 730     TpcLocalTransform(tPoint,tmpSect[ti],tRow,tU,tPhi);
 731     if (IsNaN(tmpSect[ti])){
 732       cout <<"***ERROR tmpSect"<< endl;
 733     }
 734     if (IsNaN(tRow)){
 735       cout <<"***ERROR tRow"<< endl;
 736     }
 737     if (IsNaN(tU)){
 738       cout <<"***ERROR tU"<< endl;
 739     }
 740     if (IsNaN(tPhi)){
 741       cout <<"***ERROR tPhi"<< endl;
 742     }
 743     // calculate crossing plane
 744     tn.setX(cos(tPhi));
 745     tn.setY(sin(tPhi));
 746     tr.setX(tRowRadius[ti]*cos(tPhi));
 747     tr.setY(tRowRadius[ti]*sin(tPhi));
 748     // find crossing point
 749     tLength = hel.pathLength(tr,tn);
 750     if (IsNaN(tLength)){
 751       cout <<"tLength loop 2nd  NAN" << endl;
 752       cout <<"padrow number=  " << ti << " not reached" << endl;
 753       tmpSect[ti]=-2;//sector
 754     }
 755     tPoint = hel.at(tLength);
 756     tmpZ[ti] = tPoint.z();
 757     tOutOfBound = TpcLocalTransform(tPoint,tSect,tRow,tmpU[ti],tPhi);
 758     if (IsNaN(tSect)){
 759       cout <<"***ERROR tSect 2"<< endl;
 760     }
 761     if (IsNaN(tRow)){
 762       cout <<"***ERROR tRow 2"<< endl;
 763     }
 764     if (IsNaN(tmpU[ti])){
 765       cout <<"***ERROR tmpU[ti] 2"<< endl;
 766     }
 767     if (IsNaN(tPhi)){
 768       cout <<"***ERROR tPhi 2 "<< endl;
 769     }
 770     if(tOutOfBound || (tmpSect[ti] == tSect && tRow!=(ti+1))){
 771       tmpSect[ti]=-2;
 772       //          cout << "missed once"<< endl;
 773     }
 774     else{
 775       if(tmpSect[ti] != tSect){
 776         // Try again on the other sector
 777         tn.setX(cos(tPhi));
 778         tn.setY(sin(tPhi));
 779         tr.setX(tRowRadius[ti]*cos(tPhi));
 780         tr.setY(tRowRadius[ti]*sin(tPhi));
 781         // find crossing point
 782         tLength = hel.pathLength(tr,tn);
 783         tPoint = hel.at(tLength);
 784         if (IsNaN(tLength)){
 785           cout <<"tLength loop 3rd NAN" << endl;
 786           cout <<"padrow number=  "<< ti << " not reached" << endl;
 787           tmpSect[ti]=-1;//sector
 788         }
 789         tmpZ[ti] = tPoint.z();
 790         tmpSect[ti] = tSect;
 791         tOutOfBound = TpcLocalTransform(tPoint,tSect,tRow,tmpU[ti],tPhi);
 792         if (IsNaN(tSect)){
 793           cout <<"***ERROR tSect 3"<< endl;
 794         }
 795         if (IsNaN(tRow)){
 796           cout <<"***ERROR tRow 3"<< endl;
 797         }
 798         if (IsNaN(tmpU[ti])){
 799           cout <<"***ERROR tmpU[ti] 3"<< endl;
 800         }
 801         if (IsNaN(tPhi)){
 802           cout <<"***ERROR tPhi 3 "<< endl;
 803         }
 804         if(tOutOfBound || tSect!= tmpSect[ti] || tRow!=(ti+1)){
 805           tmpSect[ti]=-1;
 806         }
 807       }
 808     }
 809     if (IsNaN(tmpSect[ti])){
 810       cout << "*******************ERROR***************************" << endl;
 811       cout <<"AliFemtoParticle--Fctn tmpSect=" << tmpSect[ti] << endl;
 812       cout << "*******************ERROR***************************" << endl;
 813     }
 814     if (IsNaN(tmpU[ti])){
 815       cout << "*******************ERROR***************************" << endl;
 816       cout <<"AliFemtoParticle--Fctn tmpU=" << tmpU[ti] << endl;
 817       cout << "*******************ERROR***************************" << endl;
 818     }
 819     if (IsNaN(tmpZ[ti])){
 820       cout << "*******************ERROR***************************" << endl;
 821       cout <<"AliFemtoParticle--Fctn tmpZ=" << tmpZ[ti] << endl;
 822       cout << "*******************ERROR***************************" << endl;
 823     }
 824     // If padrow ti not reached all other beyond are not reached
 825     // in this case set sector to -1
 826     if (tmpSect[ti]==-1){
 827       for (int tj=ti; tj<45;tj++){
 828         tmpSect[tj] = -1;
 829         ti=45;
 830       }
 831     }
 832     ti++;
 833     if (ti<45){
 834       candidates =  hel.pathLength(tRowRadius[ti]);
 835       tLength = (candidates.first > 0) ? candidates.first : candidates.second;}
 836   }
 837 }
 838 //_____________________
 839 const AliFemtoThreeVector& AliFemtoParticle::TpcV0PosExitPoint() const{
 840   return fTpcV0PosExitPoint;
 841 }
 842 //_____________________
 843 const AliFemtoThreeVector& AliFemtoParticle::TpcV0PosEntrancePoint() const{
 844   return fTpcV0PosEntrancePoint;
 845 }
 846 //______________________
 847 const AliFemtoThreeVector& AliFemtoParticle::TpcV0NegExitPoint() const{
 848   return fTpcV0NegExitPoint;
 849 }
 850 //_____________________
 851 const AliFemtoThreeVector& AliFemtoParticle::TpcV0NegEntrancePoint() const{
 852   return fTpcV0NegEntrancePoint;
 853 }
 854 //______________________