PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoEventReaderESDKine.cxx

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