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   motherids = new Int_t[fStack->GetNtrack()];
 232   for (int ip=0; ip<fStack->GetNtrack(); ip++) motherids[ip] = 0;
 233
 234   // Read in mother ids
 235   TParticle *motherpart;
 236   for (int ip=0; ip<fStack->GetNtrack(); ip++) {
 237     motherpart = fStack->Particle(ip);
 238     if (motherpart->GetDaughter(0) > 0)
 239       motherids[motherpart->GetDaughter(0)] = ip;
 240     if (motherpart->GetDaughter(1) > 0)
 241       motherids[motherpart->GetDaughter(1)] = ip;
 242
 243 //     if (motherpart->GetPdgCode() == 211) {
 244 //       cout << "Mother " << ip << " has daughters "
 245 //         << motherpart->GetDaughter(0) << " "
 246 //         << motherpart->GetDaughter(1) << " "
 247 //         << motherpart->Vx() << " "
 248 //         << motherpart->Vy() << " "
 249 //         << motherpart->Vz() << " "
 250 //         << endl;
 251
 252 //     }
 253   }
 254
 255   for (int i=0;i<nofTracks;i++)
 256     {
 257       //      cout << "Reading track " << i << endl;
 258       bool  tGoodMomentum=true; //flaga to chcek if we can read momentum of this track
 259
 260       AliFemtoTrack* trackCopy = new AliFemtoTrack();
 261       const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
 262       //      const AliESDfriendTrack *tESDfriendTrack = esdtrack->GetFriendTrack();
 263
 264       trackCopy->SetCharge((short)esdtrack->GetSign());
 265
 266       //in aliroot we have AliPID
 267       //0-electron 1-muon 2-pion 3-kaon 4-proton 5-photon 6-pi0 7-neutron 8-kaon0 9-eleCon
 268       //we use only 5 first
 269       double esdpid[5];
 270       esdtrack->GetESDpid(esdpid);
 271       trackCopy->SetPidProbElectron(esdpid[0]);
 272       trackCopy->SetPidProbMuon(esdpid[1]);
 273       trackCopy->SetPidProbPion(esdpid[2]);
 274       trackCopy->SetPidProbKaon(esdpid[3]);
 275       trackCopy->SetPidProbProton(esdpid[4]);
 276
 277       double pxyz[3];
 278       double rxyz[3];
 279       double impact[2];
 280       double covimpact[3];
 281
 282       if (fUseTPCOnly) {
 283         if (!esdtrack->GetTPCInnerParam()) {
 284           cout << "No TPC inner param !" << endl;
 285           delete trackCopy;
 286           continue;
 287         }
 288
 289         AliExternalTrackParam *param = new AliExternalTrackParam(*esdtrack->GetTPCInnerParam());
 290         param->GetXYZ(rxyz);
 291         param->PropagateToDCA(fEvent->GetPrimaryVertexTPC(), (fEvent->GetMagneticField()), 10000, impact, covimpact);
 292         param->GetPxPyPz(pxyz);//reading noconstarined momentum
 293
 294         if (fRotateToEventPlane) {
 295           double tPhi = TMath::ATan2(pxyz[1], pxyz[0]);
 296           double tRad = TMath::Hypot(pxyz[0], pxyz[1]);
 297
 298           pxyz[0] = tRad*TMath::Cos(tPhi - tReactionPlane);
 299           pxyz[1] = tRad*TMath::Sin(tPhi - tReactionPlane);
 300         }
 301
 302         AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
 303         if (v.mag() < 0.0001) {
 304           //      cout << "Found 0 momentum ???? " << pxyz[0] << " " << pxyz[1] << " " << pxyz[2] << endl;
 305           delete trackCopy;
 306           continue;
 307         }
 308
 309         trackCopy->SetP(v);//setting momentum
 310         trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
 311
 312         const AliFmThreeVectorD kP(pxyz[0],pxyz[1],pxyz[2]);
 313         const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
 314         //setting helix I do not if it is ok
 315         AliFmPhysicalHelixD helix(kP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge()));
 316         trackCopy->SetHelix(helix);
 317
 318         //some stuff which could be useful
 319         trackCopy->SetImpactD(impact[0]);
 320         trackCopy->SetImpactZ(impact[1]);
 321         trackCopy->SetCdd(covimpact[0]);
 322         trackCopy->SetCdz(covimpact[1]);
 323         trackCopy->SetCzz(covimpact[2]);
 324         trackCopy->SetSigmaToVertex(GetSigmaToVertex(impact, covimpact));
 325
 326         delete param;
 327       }
 328       else {
 329         if (fConstrained==true)
 330           tGoodMomentum=esdtrack->GetConstrainedPxPyPz(pxyz); //reading constrained momentum
 331         else
 332           tGoodMomentum=esdtrack->GetPxPyPz(pxyz);//reading noconstarined momentum
 333
 334
 335         if (fRotateToEventPlane) {
 336           double tPhi = TMath::ATan2(pxyz[1], pxyz[0]);
 337           double tRad = TMath::Hypot(pxyz[0], pxyz[1]);
 338
 339           pxyz[0] = tRad*TMath::Cos(tPhi - tReactionPlane);
 340           pxyz[1] = tRad*TMath::Sin(tPhi - tReactionPlane);
 341         }
 342
 343         AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
 344         if (v.mag() < 0.0001) {
 345           //      cout << "Found 0 momentum ???? "  << pxyz[0] << " " << pxyz[1] << " " << pxyz[2] << endl;
 346           delete trackCopy;
 347           continue;
 348         }
 349
 350         trackCopy->SetP(v);//setting momentum
 351         trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
 352         const AliFmThreeVectorD kP(pxyz[0],pxyz[1],pxyz[2]);
 353         const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
 354         //setting helix I do not if it is ok
 355         AliFmPhysicalHelixD helix(kP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge()));
 356         trackCopy->SetHelix(helix);
 357
 358         //some stuff which could be useful
 359         float imp[2];
 360         float cim[3];
 361         esdtrack->GetImpactParameters(imp,cim);
 362
 363         impact[0] = imp[0];
 364         impact[1] = imp[1];
 365         covimpact[0] = cim[0];
 366         covimpact[1] = cim[1];
 367         covimpact[2] = cim[2];
 368
 369         trackCopy->SetImpactD(impact[0]);
 370         trackCopy->SetImpactZ(impact[1]);
 371         trackCopy->SetCdd(covimpact[0]);
 372         trackCopy->SetCdz(covimpact[1]);
 373         trackCopy->SetCzz(covimpact[2]);
 374         trackCopy->SetSigmaToVertex(GetSigmaToVertex(impact,covimpact));
 375       }
 376
 377       trackCopy->SetTrackId(esdtrack->GetID());
 378       trackCopy->SetFlags(esdtrack->GetStatus());
 379       trackCopy->SetLabel(esdtrack->GetLabel());
 380
 381       trackCopy->SetITSchi2(esdtrack->GetITSchi2());
 382       trackCopy->SetITSncls(esdtrack->GetNcls(0));
 383       trackCopy->SetTPCchi2(esdtrack->GetTPCchi2());
 384       trackCopy->SetTPCncls(esdtrack->GetTPCNcls());
 385       trackCopy->SetTPCnclsF(esdtrack->GetTPCNclsF());
 386       trackCopy->SetTPCsignalN((short)esdtrack->GetTPCsignalN()); //due to bug in aliesdtrack class
 387       trackCopy->SetTPCsignalS(esdtrack->GetTPCsignalSigma());
 388
 389
 390       trackCopy->SetTPCClusterMap(esdtrack->GetTPCClusterMap());
 391       trackCopy->SetTPCSharedMap(esdtrack->GetTPCSharedMap());
 392
 393       double xtpc[3];
 394       esdtrack->GetInnerXYZ(xtpc);
 395       xtpc[2] -= fV1[2];
 396       trackCopy->SetNominalTPCEntrancePoint(xtpc);
 397
 398       esdtrack->GetOuterXYZ(xtpc);
 399       xtpc[2] -= fV1[2];
 400       trackCopy->SetNominalTPCExitPoint(xtpc);
 401
 402       int indexes[3];
 403       for (int ik=0; ik<3; ik++) {
 404         indexes[ik] = esdtrack->GetKinkIndex(ik);
 405       }
 406       trackCopy->SetKinkIndexes(indexes);
 407
 408       // Fill the hidden information with the simulated data
 409       if (TMath::Abs(esdtrack->GetLabel()) < fStack->GetNtrack()) {
 410         TParticle *tPart = fStack->Particle(TMath::Abs(esdtrack->GetLabel()));
 411
 412         // Check the mother information
 413
 414         // Using the new way of storing the freeze-out information
 415         // Final state particle is stored twice on the stack
 416         // one copy (mother) is stored with original freeze-out information
 417         //   and is not tracked
 418         // the other one (daughter) is stored with primary vertex position
 419         //   and is tracked
 420
 421         // Freeze-out coordinates
 422         double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0;
 423         fpx = tPart->Vx() - fV1[0];
 424         fpy = tPart->Vy() - fV1[1];
 425         fpz = tPart->Vz() - fV1[2];
 426         fpt = tPart->T();
 427
 428         AliFemtoModelGlobalHiddenInfo *tInfo = new AliFemtoModelGlobalHiddenInfo();
 429         tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz);
 430
 431         fpx *= 1e13;
 432         fpy *= 1e13;
 433         fpz *= 1e13;
 434         fpt *= 1e13;
 435
 436         //      cout << "Looking for mother ids " << endl;
 437         if (motherids[TMath::Abs(esdtrack->GetLabel())]>0) {
 438           //    cout << "Got mother id" << endl;
 439           TParticle *mother = fStack->Particle(motherids[TMath::Abs(esdtrack->GetLabel())]);
 440           // Check if this is the same particle stored twice on the stack
 441           if ((mother->GetPdgCode() == tPart->GetPdgCode() || (mother->Px() == tPart->Px()))) {
 442             // It is the same particle
 443             // Read in the original freeze-out information
 444             // and convert it from to [fm]
 445
 446             // EPOS style
 447             fpx = mother->Vx()*1e13*0.197327;
 448             fpy = mother->Vy()*1e13*0.197327;
 449             fpz = mother->Vz()*1e13*0.197327;
 450             fpt = mother->T() *1e13*0.197327;
 451
 452
 453             // Therminator style
 454 //          fpx = mother->Vx()*1e13;
 455 //          fpy = mother->Vy()*1e13;
 456 //          fpz = mother->Vz()*1e13;
 457 //          fpt = mother->T() *1e13*3e10;
 458
 459           }
 460         }
 461
 462         if (fRotateToEventPlane) {
 463           double tPhi = TMath::ATan2(fpy, fpx);
 464           double tRad = TMath::Hypot(fpx, fpy);
 465
 466           fpx = tRad*TMath::Cos(tPhi - tReactionPlane);
 467           fpy = tRad*TMath::Sin(tPhi - tReactionPlane);
 468         }
 469
 470         tInfo->SetPDGPid(tPart->GetPdgCode());
 471
 472         if (fRotateToEventPlane) {
 473           double tPhi = TMath::ATan2(tPart->Py(), tPart->Px());
 474           double tRad = TMath::Hypot(tPart->Px(), tPart->Py());
 475
 476           tInfo->SetTrueMomentum(tRad*TMath::Cos(tPhi - tReactionPlane),
 477                                  tRad*TMath::Sin(tPhi - tReactionPlane),
 478                                  tPart->Pz());
 479         }
 480         else
 481           tInfo->SetTrueMomentum(tPart->Px(), tPart->Py(), tPart->Pz());
 482         Double_t mass2 = (tPart->Energy() *tPart->Energy() -
 483                           tPart->Px()*tPart->Px() -
 484                           tPart->Py()*tPart->Py() -
 485                           tPart->Pz()*tPart->Pz());
 486         if (mass2>0.0)
 487           tInfo->SetMass(TMath::Sqrt(mass2));
 488         else
 489           tInfo->SetMass(0.0);
 490
 491         tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt);
 492         trackCopy->SetHiddenInfo(tInfo);
 493       }
 494       else {
 495         AliFemtoModelGlobalHiddenInfo *tInfo = new AliFemtoModelGlobalHiddenInfo();
 496         tInfo->SetMass(0.0);
 497         double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0;
 498         fpx = fV1[0]*1e13;
 499         fpy = fV1[1]*1e13;
 500         fpz = fV1[2]*1e13;
 501         fpt = 0.0;
 502         tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt);
 503
 504         tInfo->SetTrueMomentum(pxyz[0],pxyz[1],pxyz[2]);
 505
 506         trackCopy->SetHiddenInfo(tInfo);
 507       }
 508       //      cout << "Got freeze-out " << fpx << " " << fpy << " " << fpz << " " << fpt << " " <<  mass2 << " " << tPart->GetPdgCode() << endl;
 509
 510       //decision if we want this track
 511       //if we using diffrent labels we want that this label was use for first time
 512       //if we use hidden info we want to have match between sim data and ESD
 513       if (tGoodMomentum==true)
 514         {
 515           hbtEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
 516           realnofTracks++;//real number of tracks
 517           //      delete trackCopy;
 518         }
 519       else
 520         {
 521           delete  trackCopy;
 522         }
 523
 524     }
 525
 526   delete motherids;
 527
 528   hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event
 529   fCurEvent++;
 530   //  cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
 531   return hbtEvent;
 532 }
 533 //___________________
 534 void AliFemtoEventReaderESDChainKine::SetESDSource(AliESDEvent *aESD)
 535 {
 536   // The chain loads the ESD for us
 537   // You must provide the address where it can be found
 538   fEvent = aESD;
 539 }
 540 //___________________
 541 void AliFemtoEventReaderESDChainKine::SetStackSource(AliStack *aStack)
 542 {
 543   // The chain loads the stack for us
 544   // You must provide the address where it can be found
 545   fStack = aStack;
 546 }
 547 //___________________
 548 void AliFemtoEventReaderESDChainKine::SetGenEventHeader(AliGenEventHeader *aGenHeader)
 549 {
 550   // The chain loads the generator event header for us
 551   // You must provide the address where it can be found
 552   fGenHeader = aGenHeader;
 553 }
 554
 555 //__________________
 556 void AliFemtoEventReaderESDChainKine::SetUseTPCOnly(const bool usetpconly)
 557 {
 558   fUseTPCOnly=usetpconly;
 559 }
 560 void AliFemtoEventReaderESDChainKine::SetRotateToEventPlane(short dorotate)
 561 {
 562   fRotateToEventPlane=dorotate;
 563 }
 564
 565 bool AliFemtoEventReaderESDChainKine::GetUseTPCOnly() const
 566 {
 567   return fUseTPCOnly;
 568 }
 569 Float_t AliFemtoEventReaderESDChainKine::GetSigmaToVertex(double *impact, double *covar)
 570 {
 571   // Calculates the number of sigma to the vertex.
 572
 573   Float_t b[2];
 574   Float_t bRes[2];
 575   Float_t bCov[3];
 576
 577   b[0] = impact[0];
 578   b[1] = impact[1];
 579   bCov[0] = covar[0];
 580   bCov[1] = covar[1];
 581   bCov[2] = covar[2];
 582
 583   bRes[0] = TMath::Sqrt(bCov[0]);
 584   bRes[1] = TMath::Sqrt(bCov[2]);
 585
 586   // -----------------------------------
 587   // How to get to a n-sigma cut?
 588   //
 589   // The accumulated statistics from 0 to d is
 590   //
 591   // ->  Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
 592   // ->  1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
 593   //
 594   // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-x**2)/2)
 595   // Can this be expressed in a different way?
 596
 597   if (bRes[0] == 0 || bRes[1] ==0)
 598     return -1;
 599
 600   Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
 601
 602   // stupid rounding problem screws up everything:
 603   // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
 604   if (TMath::Exp(-d * d / 2) < 1e-10)
 605     return 1000;
 606
 607   d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
 608   return d;
 609 }