PWG2/FEMTOSCOPY/AliFemto/AliFemtoEventReaderESDChainKine.cxx

   1 ////////////////////////////////////////////////////////////////////////////////
   2 //                                                                            //
   3 // AliFemtoEventReaderESDChainKine - the reader class for the Alice ESD and   //
   4 // the model Kinematics information tailored for the Task framework and the   //
   5 // Reads in AliESDfriend to create shared hit/quality information             //
   6 // Authors: Adam Kisiel kisiel@mps.ohio-state.edu                             //
   7 //                                                                            //
   8 ////////////////////////////////////////////////////////////////////////////////
   9 #include "AliFemtoEventReaderESDChainKine.h"
  10
  11 #include "TFile.h"
  12 #include "TTree.h"
  13 #include "TList.h"
  14 #include "AliESDEvent.h"
  15 #include "AliESDtrack.h"
  16 #include "AliESDVertex.h"
  17
  18 #include "AliFmPhysicalHelixD.h"
  19 #include "AliFmThreeVectorF.h"
  20
  21 #include "SystemOfUnits.h"
  22
  23 #include "AliFemtoEvent.h"
  24
  25 #include "TParticle.h"
  26 #include "AliFemtoModelHiddenInfo.h"
  27 #include "AliFemtoModelGlobalHiddenInfo.h"
  28 #include "AliGenHijingEventHeader.h"
  29 #include "AliGenCocktailEventHeader.h"
  30
  31 #include "AliVertexerTracks.h"
  32
  33 ClassImp(AliFemtoEventReaderESDChainKine)
  34
  35 #if !(ST_NO_NAMESPACES)
  36   using namespace units;
  37 #endif
  38
  39 using namespace std;
  40 //____________________________
  41 AliFemtoEventReaderESDChainKine::AliFemtoEventReaderESDChainKine():
  42   fFileName(" "),
  43   fConstrained(true),
  44   fUseTPCOnly(false),
  45   fNumberofEvent(0),
  46   fCurEvent(0),
  47   fCurFile(0),
  48   fEvent(0x0),
  49   fStack(0x0),
  50   fGenHeader(0x0),
  51   fRotateToEventPlane(0)
  52 {
  53   //constructor with 0 parameters , look at default settings
  54 }
  55
  56 //__________________
  57 AliFemtoEventReaderESDChainKine::AliFemtoEventReaderESDChainKine(const AliFemtoEventReaderESDChainKine& aReader):
  58   AliFemtoEventReader(aReader),
  59   fFileName(" "),
  60   fConstrained(true),
  61   fUseTPCOnly(false),
  62   fNumberofEvent(0),
  63   fCurEvent(0),
  64   fCurFile(0),
  65   fEvent(0x0),
  66   fStack(0x0),
  67   fGenHeader(0x0),
  68   fRotateToEventPlane(0)
  69 {
  70   // Copy constructor
  71   fConstrained = aReader.fConstrained;
  72   fUseTPCOnly = aReader.fUseTPCOnly;
  73   fNumberofEvent = aReader.fNumberofEvent;
  74   fCurEvent = aReader.fCurEvent;
  75   fCurFile = aReader.fCurFile;
  76   fEvent = new AliESDEvent();
  77   fStack = aReader.fStack;
  78   fRotateToEventPlane = aReader.fRotateToEventPlane;
  79 }
  80 //__________________
  81 AliFemtoEventReaderESDChainKine::~AliFemtoEventReaderESDChainKine()
  82 {
  83   //Destructor
  84   delete fEvent;
  85 }
  86
  87 //__________________
  88 AliFemtoEventReaderESDChainKine& AliFemtoEventReaderESDChainKine::operator=(const AliFemtoEventReaderESDChainKine& aReader)
  89 {
  90   // Assignment operator
  91   if (this == &aReader)
  92     return *this;
  93
  94   fConstrained = aReader.fConstrained;
  95   fUseTPCOnly = aReader.fUseTPCOnly;
  96   fNumberofEvent = aReader.fNumberofEvent;
  97   fCurEvent = aReader.fCurEvent;
  98   fCurFile = aReader.fCurFile;
  99   if (fEvent) delete fEvent;
 100   fEvent = new AliESDEvent();
 101   fStack = aReader.fStack;
 102   fGenHeader = aReader.fGenHeader;
 103   fRotateToEventPlane = aReader.fRotateToEventPlane;
 104   return *this;
 105 }
 106 //__________________
 107 // Simple report
 108 AliFemtoString AliFemtoEventReaderESDChainKine::Report()
 109 {
 110   AliFemtoString temp = "\n This is the AliFemtoEventReaderESDChainKine\n";
 111   return temp;
 112 }
 113
 114 //__________________
 115 void AliFemtoEventReaderESDChainKine::SetConstrained(const bool constrained)
 116 {
 117   // Select whether to read constrained or not constrained momentum
 118   fConstrained=constrained;
 119 }
 120 //__________________
 121 bool AliFemtoEventReaderESDChainKine::GetConstrained() const
 122 {
 123   // Check whether we read constrained or not constrained momentum
 124   return fConstrained;
 125 }
 126 //__________________
 127 AliFemtoEvent* AliFemtoEventReaderESDChainKine::ReturnHbtEvent()
 128 {
 129   // Get the event, read all the relevant information
 130   // and fill the AliFemtoEvent class
 131   // Returns a valid AliFemtoEvent
 132   AliFemtoEvent *hbtEvent = 0;
 133   string tFriendFileName;
 134
 135   // Get the friend information
 136   cout << "AliFemtoEventReaderESDChainKine::Starting to read event: "<<fCurEvent<<endl;
 137   //  fEvent->SetESDfriend(fEventFriend);
 138
 139   hbtEvent = new AliFemtoEvent;
 140   //setting basic things
 141   //  hbtEvent->SetEventNumber(fEvent->GetEventNumber());
 142   hbtEvent->SetRunNumber(fEvent->GetRunNumber());
 143   //hbtEvent->SetNumberOfTracks(fEvent->GetNumberOfTracks());
 144   hbtEvent->SetMagneticField(fEvent->GetMagneticField()*kilogauss);//to check if here is ok
 145   hbtEvent->SetZDCN1Energy(fEvent->GetZDCN1Energy());
 146   hbtEvent->SetZDCP1Energy(fEvent->GetZDCP1Energy());
 147   hbtEvent->SetZDCN2Energy(fEvent->GetZDCN2Energy());
 148   hbtEvent->SetZDCP2Energy(fEvent->GetZDCP2Energy());
 149   hbtEvent->SetZDCEMEnergy(fEvent->GetZDCEMEnergy());
 150   hbtEvent->SetZDCParticipants(fEvent->GetZDCParticipants());
 151   hbtEvent->SetTriggerMask(fEvent->GetTriggerMask());
 152   hbtEvent->SetTriggerCluster(fEvent->GetTriggerCluster());
 153
 154   //Vertex
 155   double fV1[3];
 156   double fVCov[6];
 157   if (fUseTPCOnly) {
 158     fEvent->GetPrimaryVertexTPC()->GetXYZ(fV1);
 159     fEvent->GetPrimaryVertexTPC()->GetCovMatrix(fVCov);
 160     if (!fEvent->GetPrimaryVertexTPC()->GetStatus())
 161       fVCov[4] = -1001.0;
 162   }
 163   else {
 164     if (fEvent->GetPrimaryVertex()) {
 165       fEvent->GetPrimaryVertex()->GetXYZ(fV1);
 166       fEvent->GetPrimaryVertex()->GetCovMatrix(fVCov);
 167
 168       if (!fEvent->GetPrimaryVertex()->GetStatus()) {
 169         // Get the vertex from SPD
 170         fEvent->GetPrimaryVertexSPD()->GetXYZ(fV1);
 171         fEvent->GetPrimaryVertexSPD()->GetCovMatrix(fVCov);
 172
 173
 174         if (!fEvent->GetPrimaryVertexSPD()->GetStatus())
 175           fVCov[4] = -1001.0;
 176         else {
 177           fEvent->GetPrimaryVertexSPD()->GetXYZ(fV1);
 178           fEvent->GetPrimaryVertexSPD()->GetCovMatrix(fVCov);
 179         }
 180       }
 181     }
 182     else {
 183       if (fEvent->GetPrimaryVertexSPD()) {
 184         fEvent->GetPrimaryVertexSPD()->GetXYZ(fV1);
 185         fEvent->GetPrimaryVertexSPD()->GetCovMatrix(fVCov);
 186       }
 187     }
 188     if ((!fEvent->GetPrimaryVertex()) && (!fEvent->GetPrimaryVertexSPD()))
 189       {
 190         cout << "No vertex found !!!" << endl;
 191         fV1[0] = 10000.0;
 192         fV1[1] = 10000.0;
 193         fV1[2] = 10000.0;
 194         fVCov[4] = -1001.0;
 195       }
 196   }
 197
 198   AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]);
 199
 200   hbtEvent->SetPrimVertPos(vertex);
 201   hbtEvent->SetPrimVertCov(fVCov);
 202
 203   Double_t tReactionPlane = 0;
 204
 205   AliGenHijingEventHeader *hdh = dynamic_cast<AliGenHijingEventHeader *> (fGenHeader);
 206   if (!hdh) {
 207     AliGenCocktailEventHeader *cdh = dynamic_cast<AliGenCocktailEventHeader *> (fGenHeader);
 208     if (cdh) {
 209       TList *tGenHeaders = cdh->GetHeaders();
 210       for (int ihead = 0; ihead<tGenHeaders->GetEntries(); ihead++) {
 211         hdh = dynamic_cast<AliGenHijingEventHeader *> (fGenHeader);
 212         if (hdh) break;
 213       }
 214     }
 215   }
 216
 217   if (hdh)
 218     {
 219       tReactionPlane = hdh->ReactionPlaneAngle();
 220       cout << "Got reaction plane " << tReactionPlane << endl;
 221     }
 222
 223   hbtEvent->SetReactionPlaneAngle(tReactionPlane);
 224
 225   //starting to reading tracks
 226   int nofTracks=0;  //number of reconstructed tracks in event
 227   nofTracks=fEvent->GetNumberOfTracks();
 228   int realnofTracks=0;//number of track which we use ina analysis
 229
 230   Int_t *motherids;
 231   if (fStack) {
 232     motherids = new Int_t[fStack->GetNtrack()];
 233     for (int ip=0; ip<fStack->GetNtrack(); ip++) motherids[ip] = 0;
 234
 235     // Read in mother ids
 236     TParticle *motherpart;
 237     for (int ip=0; ip<fStack->GetNtrack(); ip++) {
 238       motherpart = fStack->Particle(ip);
 239       if (motherpart->GetDaughter(0) > 0)
 240         motherids[motherpart->GetDaughter(0)] = ip;
 241       if (motherpart->GetDaughter(1) > 0)
 242         motherids[motherpart->GetDaughter(1)] = ip;
 243
 244 //     if (motherpart->GetPdgCode() == 211) {
 245 //       cout << "Mother " << ip << " has daughters "
 246 //         << motherpart->GetDaughter(0) << " "
 247 //         << motherpart->GetDaughter(1) << " "
 248 //         << motherpart->Vx() << " "
 249 //         << motherpart->Vy() << " "
 250 //         << motherpart->Vz() << " "
 251 //         << endl;
 252
 253 //     }
 254     }
 255   }
 256   else {
 257     printf ("No Stack ???\n");
 258     return 0;
 259   }
 260
 261   for (int i=0;i<nofTracks;i++)
 262     {
 263       //      cout << "Reading track " << i << endl;
 264       bool  tGoodMomentum=true; //flaga to chcek if we can read momentum of this track
 265
 266       AliFemtoTrack* trackCopy = new AliFemtoTrack();
 267       const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
 268       //      const AliESDfriendTrack *tESDfriendTrack = esdtrack->GetFriendTrack();
 269
 270       trackCopy->SetCharge((short)esdtrack->GetSign());
 271
 272       //in aliroot we have AliPID
 273       //0-electron 1-muon 2-pion 3-kaon 4-proton 5-photon 6-pi0 7-neutron 8-kaon0 9-eleCon
 274       //we use only 5 first
 275       double esdpid[5];
 276       esdtrack->GetESDpid(esdpid);
 277       trackCopy->SetPidProbElectron(esdpid[0]);
 278       trackCopy->SetPidProbMuon(esdpid[1]);
 279       trackCopy->SetPidProbPion(esdpid[2]);
 280       trackCopy->SetPidProbKaon(esdpid[3]);
 281       trackCopy->SetPidProbProton(esdpid[4]);
 282
 283       double pxyz[3];
 284       double rxyz[3];
 285       double impact[2];
 286       double covimpact[3];
 287
 288       if (fUseTPCOnly) {
 289         if (!esdtrack->GetTPCInnerParam()) {
 290           cout << "No TPC inner param !" << endl;
 291           delete trackCopy;
 292           continue;
 293         }
 294
 295         AliExternalTrackParam *param = new AliExternalTrackParam(*esdtrack->GetTPCInnerParam());
 296         param->GetXYZ(rxyz);
 297         param->PropagateToDCA(fEvent->GetPrimaryVertexTPC(), (fEvent->GetMagneticField()), 10000, impact, covimpact);
 298         param->GetPxPyPz(pxyz);//reading noconstarined momentum
 299
 300         if (fRotateToEventPlane) {
 301           double tPhi = TMath::ATan2(pxyz[1], pxyz[0]);
 302           double tRad = TMath::Hypot(pxyz[0], pxyz[1]);
 303
 304           pxyz[0] = tRad*TMath::Cos(tPhi - tReactionPlane);
 305           pxyz[1] = tRad*TMath::Sin(tPhi - tReactionPlane);
 306         }
 307
 308         AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
 309         if (v.Mag() < 0.0001) {
 310           //      cout << "Found 0 momentum ???? " << pxyz[0] << " " << pxyz[1] << " " << pxyz[2] << endl;
 311           delete trackCopy;
 312           continue;
 313         }
 314
 315         trackCopy->SetP(v);//setting momentum
 316         trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
 317
 318         const AliFmThreeVectorD kP(pxyz[0],pxyz[1],pxyz[2]);
 319         const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
 320         //setting helix I do not if it is ok
 321         AliFmPhysicalHelixD helix(kP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge()));
 322         trackCopy->SetHelix(helix);
 323
 324         //some stuff which could be useful
 325         trackCopy->SetImpactD(impact[0]);
 326         trackCopy->SetImpactZ(impact[1]);
 327         trackCopy->SetCdd(covimpact[0]);
 328         trackCopy->SetCdz(covimpact[1]);
 329         trackCopy->SetCzz(covimpact[2]);
 330         trackCopy->SetSigmaToVertex(GetSigmaToVertex(impact, covimpact));
 331
 332         delete param;
 333       }
 334       else {
 335         if (fConstrained==true)
 336           tGoodMomentum=esdtrack->GetConstrainedPxPyPz(pxyz); //reading constrained momentum
 337         else
 338           tGoodMomentum=esdtrack->GetPxPyPz(pxyz);//reading noconstarined momentum
 339
 340
 341         if (fRotateToEventPlane) {
 342           double tPhi = TMath::ATan2(pxyz[1], pxyz[0]);
 343           double tRad = TMath::Hypot(pxyz[0], pxyz[1]);
 344
 345           pxyz[0] = tRad*TMath::Cos(tPhi - tReactionPlane);
 346           pxyz[1] = tRad*TMath::Sin(tPhi - tReactionPlane);
 347         }
 348
 349         AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
 350         if (v.Mag() < 0.0001) {
 351           //      cout << "Found 0 momentum ???? "  << pxyz[0] << " " << pxyz[1] << " " << pxyz[2] << endl;
 352           delete trackCopy;
 353           continue;
 354         }
 355
 356         trackCopy->SetP(v);//setting momentum
 357         trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
 358         const AliFmThreeVectorD kP(pxyz[0],pxyz[1],pxyz[2]);
 359         const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
 360         //setting helix I do not if it is ok
 361         AliFmPhysicalHelixD helix(kP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge()));
 362         trackCopy->SetHelix(helix);
 363
 364         //some stuff which could be useful
 365         float imp[2];
 366         float cim[3];
 367         esdtrack->GetImpactParameters(imp,cim);
 368
 369         impact[0] = imp[0];
 370         impact[1] = imp[1];
 371         covimpact[0] = cim[0];
 372         covimpact[1] = cim[1];
 373         covimpact[2] = cim[2];
 374
 375         trackCopy->SetImpactD(impact[0]);
 376         trackCopy->SetImpactZ(impact[1]);
 377         trackCopy->SetCdd(covimpact[0]);
 378         trackCopy->SetCdz(covimpact[1]);
 379         trackCopy->SetCzz(covimpact[2]);
 380         trackCopy->SetSigmaToVertex(GetSigmaToVertex(impact,covimpact));
 381       }
 382
 383       trackCopy->SetTrackId(esdtrack->GetID());
 384       trackCopy->SetFlags(esdtrack->GetStatus());
 385       trackCopy->SetLabel(esdtrack->GetLabel());
 386
 387       trackCopy->SetITSchi2(esdtrack->GetITSchi2());
 388       trackCopy->SetITSncls(esdtrack->GetNcls(0));
 389       trackCopy->SetTPCchi2(esdtrack->GetTPCchi2());
 390       trackCopy->SetTPCncls(esdtrack->GetTPCNcls());
 391       trackCopy->SetTPCnclsF(esdtrack->GetTPCNclsF());
 392       trackCopy->SetTPCsignalN((short)esdtrack->GetTPCsignalN()); //due to bug in aliesdtrack class
 393       trackCopy->SetTPCsignalS(esdtrack->GetTPCsignalSigma());
 394
 395
 396       trackCopy->SetTPCClusterMap(esdtrack->GetTPCClusterMap());
 397       trackCopy->SetTPCSharedMap(esdtrack->GetTPCSharedMap());
 398
 399       double xtpc[3];
 400       esdtrack->GetInnerXYZ(xtpc);
 401       xtpc[2] -= fV1[2];
 402       trackCopy->SetNominalTPCEntrancePoint(xtpc);
 403
 404       esdtrack->GetOuterXYZ(xtpc);
 405       xtpc[2] -= fV1[2];
 406       trackCopy->SetNominalTPCExitPoint(xtpc);
 407
 408       int indexes[3];
 409       for (int ik=0; ik<3; ik++) {
 410         indexes[ik] = esdtrack->GetKinkIndex(ik);
 411       }
 412       trackCopy->SetKinkIndexes(indexes);
 413
 414       // Fill the hidden information with the simulated data
 415       if (TMath::Abs(esdtrack->GetLabel()) < fStack->GetNtrack()) {
 416         TParticle *tPart = fStack->Particle(TMath::Abs(esdtrack->GetLabel()));
 417
 418         // Check the mother information
 419
 420         // Using the new way of storing the freeze-out information
 421         // Final state particle is stored twice on the stack
 422         // one copy (mother) is stored with original freeze-out information
 423         //   and is not tracked
 424         // the other one (daughter) is stored with primary vertex position
 425         //   and is tracked
 426
 427         // Freeze-out coordinates
 428         double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0;
 429         fpx = tPart->Vx() - fV1[0];
 430         fpy = tPart->Vy() - fV1[1];
 431         fpz = tPart->Vz() - fV1[2];
 432         fpt = tPart->T();
 433
 434         AliFemtoModelGlobalHiddenInfo *tInfo = new AliFemtoModelGlobalHiddenInfo();
 435         tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz);
 436
 437         fpx *= 1e13;
 438         fpy *= 1e13;
 439         fpz *= 1e13;
 440         fpt *= 1e13;
 441
 442         //      cout << "Looking for mother ids " << endl;
 443         if (motherids[TMath::Abs(esdtrack->GetLabel())]>0) {
 444           //    cout << "Got mother id" << endl;
 445           TParticle *mother = fStack->Particle(motherids[TMath::Abs(esdtrack->GetLabel())]);
 446           // Check if this is the same particle stored twice on the stack
 447           if ((mother->GetPdgCode() == tPart->GetPdgCode() || (mother->Px() == tPart->Px()))) {
 448             // It is the same particle
 449             // Read in the original freeze-out information
 450             // and convert it from to [fm]
 451
 452             // EPOS style
 453             fpx = mother->Vx()*1e13*0.197327;
 454             fpy = mother->Vy()*1e13*0.197327;
 455             fpz = mother->Vz()*1e13*0.197327;
 456             fpt = mother->T() *1e13*0.197327;
 457
 458
 459             // Therminator style
 460 //          fpx = mother->Vx()*1e13;
 461 //          fpy = mother->Vy()*1e13;
 462 //          fpz = mother->Vz()*1e13;
 463 //          fpt = mother->T() *1e13*3e10;
 464
 465           }
 466         }
 467
 468         if (fRotateToEventPlane) {
 469           double tPhi = TMath::ATan2(fpy, fpx);
 470           double tRad = TMath::Hypot(fpx, fpy);
 471
 472           fpx = tRad*TMath::Cos(tPhi - tReactionPlane);
 473           fpy = tRad*TMath::Sin(tPhi - tReactionPlane);
 474         }
 475
 476         tInfo->SetPDGPid(tPart->GetPdgCode());
 477
 478         if (fRotateToEventPlane) {
 479           double tPhi = TMath::ATan2(tPart->Py(), tPart->Px());
 480           double tRad = TMath::Hypot(tPart->Px(), tPart->Py());
 481
 482           tInfo->SetTrueMomentum(tRad*TMath::Cos(tPhi - tReactionPlane),
 483                                  tRad*TMath::Sin(tPhi - tReactionPlane),
 484                                  tPart->Pz());
 485         }
 486         else
 487           tInfo->SetTrueMomentum(tPart->Px(), tPart->Py(), tPart->Pz());
 488         Double_t mass2 = (tPart->Energy() *tPart->Energy() -
 489                           tPart->Px()*tPart->Px() -
 490                           tPart->Py()*tPart->Py() -
 491                           tPart->Pz()*tPart->Pz());
 492         if (mass2>0.0)
 493           tInfo->SetMass(TMath::Sqrt(mass2));
 494         else
 495           tInfo->SetMass(0.0);
 496
 497         tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt);
 498         trackCopy->SetHiddenInfo(tInfo);
 499       }
 500       else {
 501         AliFemtoModelGlobalHiddenInfo *tInfo = new AliFemtoModelGlobalHiddenInfo();
 502         tInfo->SetMass(0.0);
 503         double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0;
 504         fpx = fV1[0]*1e13;
 505         fpy = fV1[1]*1e13;
 506         fpz = fV1[2]*1e13;
 507         fpt = 0.0;
 508         tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt);
 509
 510         tInfo->SetTrueMomentum(pxyz[0],pxyz[1],pxyz[2]);
 511
 512         trackCopy->SetHiddenInfo(tInfo);
 513       }
 514       //      cout << "Got freeze-out " << fpx << " " << fpy << " " << fpz << " " << fpt << " " <<  mass2 << " " << tPart->GetPdgCode() << endl;
 515
 516       //decision if we want this track
 517       //if we using diffrent labels we want that this label was use for first time
 518       //if we use hidden info we want to have match between sim data and ESD
 519       if (tGoodMomentum==true)
 520         {
 521           hbtEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
 522           realnofTracks++;//real number of tracks
 523           //      delete trackCopy;
 524         }
 525       else
 526         {
 527           delete  trackCopy;
 528         }
 529
 530     }
 531
 532   delete [] motherids;
 533
 534   hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event
 535   fCurEvent++;
 536   //  cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
 537   return hbtEvent;
 538 }
 539 //___________________
 540 void AliFemtoEventReaderESDChainKine::SetESDSource(AliESDEvent *aESD)
 541 {
 542   // The chain loads the ESD for us
 543   // You must provide the address where it can be found
 544   fEvent = aESD;
 545 }
 546 //___________________
 547 void AliFemtoEventReaderESDChainKine::SetStackSource(AliStack *aStack)
 548 {
 549   // The chain loads the stack for us
 550   // You must provide the address where it can be found
 551   fStack = aStack;
 552 }
 553 //___________________
 554 void AliFemtoEventReaderESDChainKine::SetGenEventHeader(AliGenEventHeader *aGenHeader)
 555 {
 556   // The chain loads the generator event header for us
 557   // You must provide the address where it can be found
 558   fGenHeader = aGenHeader;
 559 }
 560
 561 //__________________
 562 void AliFemtoEventReaderESDChainKine::SetUseTPCOnly(const bool usetpconly)
 563 {
 564   fUseTPCOnly=usetpconly;
 565 }
 566 void AliFemtoEventReaderESDChainKine::SetRotateToEventPlane(short dorotate)
 567 {
 568   fRotateToEventPlane=dorotate;
 569 }
 570
 571 bool AliFemtoEventReaderESDChainKine::GetUseTPCOnly() const
 572 {
 573   return fUseTPCOnly;
 574 }
 575 Float_t AliFemtoEventReaderESDChainKine::GetSigmaToVertex(double *impact, double *covar)
 576 {
 577   // Calculates the number of sigma to the vertex.
 578
 579   Float_t b[2];
 580   Float_t bRes[2];
 581   Float_t bCov[3];
 582
 583   b[0] = impact[0];
 584   b[1] = impact[1];
 585   bCov[0] = covar[0];
 586   bCov[1] = covar[1];
 587   bCov[2] = covar[2];
 588
 589   bRes[0] = TMath::Sqrt(bCov[0]);
 590   bRes[1] = TMath::Sqrt(bCov[2]);
 591
 592   // -----------------------------------
 593   // How to get to a n-sigma cut?
 594   //
 595   // The accumulated statistics from 0 to d is
 596   //
 597   // ->  Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
 598   // ->  1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
 599   //
 600   // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-x**2)/2)
 601   // Can this be expressed in a different way?
 602
 603   if (bRes[0] == 0 || bRes[1] ==0)
 604     return -1;
 605
 606   Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
 607
 608   // stupid rounding problem screws up everything:
 609   // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
 610   if (TMath::Exp(-d * d / 2) < 1e-10)
 611     return 1000;
 612
 613   d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
 614   return d;
 615 }