PWGCF/FEMTOSCOPY/AliFemto/AliFemtoEventReaderESDChain.cxx

   1 ////////////////////////////////////////////////////////////////////////////////
   2 ///                                                                          ///
   3 /// AliFemtoEventReaderESDChain - the reader class for the Alice ESD         ///
   4 /// tailored for the Task framework                                 ///
   5 /// Reads in AliESDfriend to create shared hit/quality information           ///
   6 /// Authors: Adam Kisiel kisiel@mps.ohio-state.edu                           ///
   7 ///                                                                          ///
   8 ////////////////////////////////////////////////////////////////////////////////
   9 #include "AliFemtoEventReaderESDChain.h"
  10
  11 #include "TFile.h"
  12 #include "TTree.h"
  13 #include "AliESDEvent.h"
  14 #include "AliESDtrack.h"
  15 #include "AliESDVertex.h"
  16 #include "AliMultiplicity.h"
  17 #include "AliCentrality.h"
  18 #include "AliEventplane.h"
  19 #include "AliESDVZERO.h"
  20 #include "AliFmPhysicalHelixD.h"
  21 #include "AliFmThreeVectorF.h"
  22 #include "SystemOfUnits.h"
  23 #include "AliFemtoEvent.h"
  24 #include "AliFemtoModelHiddenInfo.h"
  25 #include "AliPID.h"
  26
  27 ClassImp(AliFemtoEventReaderESDChain)
  28
  29 #if !(ST_NO_NAMESPACES)
  30   using namespace units;
  31 #endif
  32
  33 using namespace std;
  34 //____________________________
  35 AliFemtoEventReaderESDChain::AliFemtoEventReaderESDChain():
  36   fFileName(" "),
  37   fConstrained(true),
  38   fReadInner(false),
  39   fUseTPCOnly(false),
  40   fNumberofEvent(0),
  41   fCurEvent(0),
  42   fCurFile(0),
  43   fEvent(0x0),
  44   fUsePhysicsSel(kFALSE),
  45   fSelect(0x0),
  46   fTrackType(kGlobal),
  47   fEstEventMult(kITSTPC),
  48   fEventTrig(AliVEvent::kMB), //trigger
  49   fESDpid(0),
  50   fIsPidOwner(0)
  51 {
  52   //constructor with 0 parameters , look at default settings
  53   //   fClusterPerPadrow = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
  54   //   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
  55   //     fClusterPerPadrow[tPad] = new list<Int_t>();
  56   //   }
  57   //   fSharedList = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
  58   //   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
  59   //     fSharedList[tPad] = new list<Int_t>();
  60   //   }
  61 }
  62
  63 //__________________
  64 AliFemtoEventReaderESDChain::AliFemtoEventReaderESDChain(const AliFemtoEventReaderESDChain& aReader):
  65   AliFemtoEventReader(aReader),
  66   fFileName(" "),
  67   fConstrained(true),
  68   fReadInner(false),
  69   fUseTPCOnly(false),
  70   fNumberofEvent(0),
  71   fCurEvent(0),
  72   fCurFile(0),
  73   fEvent(0x0),
  74   fUsePhysicsSel(kFALSE),
  75   fSelect(0x0),
  76   fTrackType(kGlobal),
  77   fEstEventMult(kITSTPC),
  78   fEventTrig(AliVEvent::kMB), //trigger
  79   fESDpid(0),
  80   fIsPidOwner(0)
  81 {
  82   // Copy constructor
  83   fConstrained = aReader.fConstrained;
  84   fReadInner = aReader.fReadInner;
  85   fUseTPCOnly = aReader.fUseTPCOnly;
  86   fNumberofEvent = aReader.fNumberofEvent;
  87   fCurEvent = aReader.fCurEvent;
  88   fCurFile = aReader.fCurFile;
  89   //  fEvent = new AliESD(*aReader.fEvent);
  90   fEvent = new AliESDEvent();
  91   fUsePhysicsSel = aReader.fUsePhysicsSel;
  92   if (aReader.fUsePhysicsSel)
  93     fSelect = new AliPhysicsSelection();
  94   fTrackType = aReader.fTrackType;
  95   fEstEventMult = aReader.fEstEventMult;
  96   fEventTrig = aReader.fEventTrig; //trigger
  97
  98   //   fEventFriend = aReader.fEventFriend;
  99   //   fClusterPerPadrow = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
 100   //   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
 101   //     fClusterPerPadrow[tPad] = new list<Int_t>();
 102   //     list<Int_t>::iterator iter;
 103   //     for (iter=aReader.fClusterPerPadrow[tPad]->begin(); iter!=aReader.fClusterPerPadrow[tPad]->end(); iter++) {
 104   //       fClusterPerPadrow[tPad]->push_back(*iter);
 105   //     }
 106   //   }
 107   //   fSharedList = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
 108   //   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
 109   //     fSharedList[tPad] = new list<Int_t>();
 110   //     list<Int_t>::iterator iter;
 111   //     for (iter=aReader.fSharedList[tPad]->begin(); iter!=aReader.fSharedList[tPad]->end(); iter++) {
 112   //       fSharedList[tPad]->push_back(*iter);
 113   //     }
 114   //   }
 115 }
 116 //__________________
 117 AliFemtoEventReaderESDChain::~AliFemtoEventReaderESDChain()
 118 {
 119   //Destructor
 120   delete fEvent;
 121
 122   //   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
 123   //     fClusterPerPadrow[tPad]->clear();
 124   //     delete fClusterPerPadrow[tPad];
 125   //   }
 126   //   delete [] fClusterPerPadrow;
 127   //   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
 128   //     fSharedList[tPad]->clear();
 129   //     delete fSharedList[tPad];
 130   //   }
 131   //   delete [] fSharedList;
 132   if (fSelect) delete fSelect;
 133 }
 134
 135 //__________________
 136 AliFemtoEventReaderESDChain& AliFemtoEventReaderESDChain::operator=(const AliFemtoEventReaderESDChain& aReader)
 137 {
 138   // Assignment operator
 139   if (this == &aReader)
 140     return *this;
 141
 142   fConstrained = aReader.fConstrained;
 143   fReadInner = aReader.fReadInner;
 144   fUseTPCOnly = aReader.fUseTPCOnly;
 145   fNumberofEvent = aReader.fNumberofEvent;
 146   fCurEvent = aReader.fCurEvent;
 147   fCurFile = aReader.fCurFile;
 148   if (fEvent) delete fEvent;
 149   fEvent = new AliESDEvent();
 150   fTrackType = aReader.fTrackType;
 151   fEstEventMult = aReader.fEstEventMult;
 152
 153   fUsePhysicsSel = aReader.fUsePhysicsSel;
 154   if (aReader.fUsePhysicsSel)
 155     fSelect = new AliPhysicsSelection();
 156   //  fEventFriend = aReader.fEventFriend;
 157
 158   //   if (fClusterPerPadrow) {
 159   //     for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
 160   //       fClusterPerPadrow[tPad]->clear();
 161   //       delete fClusterPerPadrow[tPad];
 162   //     }
 163   //     delete [] fClusterPerPadrow;
 164   //   }
 165
 166   //   if (fSharedList) {
 167   //     for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
 168   //       fSharedList[tPad]->clear();
 169   //       delete fSharedList[tPad];
 170   //     }
 171   //     delete [] fSharedList;
 172   //   }
 173
 174   //   fClusterPerPadrow = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
 175   //   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
 176   //     fClusterPerPadrow[tPad] = new list<Int_t>();
 177   //     list<Int_t>::iterator iter;
 178   //     for (iter=aReader.fClusterPerPadrow[tPad]->begin(); iter!=aReader.fClusterPerPadrow[tPad]->end(); iter++) {
 179   //       fClusterPerPadrow[tPad]->push_back(*iter);
 180   //     }
 181   //   }
 182   //   fSharedList = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
 183   //   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
 184   //     fSharedList[tPad] = new list<Int_t>();
 185   //     list<Int_t>::iterator iter;
 186   //     for (iter=aReader.fSharedList[tPad]->begin(); iter!=aReader.fSharedList[tPad]->end(); iter++) {
 187   //       fSharedList[tPad]->push_back(*iter);
 188   //     }
 189   //   }
 190
 191   return *this;
 192 }
 193 //__________________
 194 // Simple report
 195 AliFemtoString AliFemtoEventReaderESDChain::Report()
 196 {
 197   AliFemtoString temp = "\n This is the AliFemtoEventReaderESDChain\n";
 198   return temp;
 199 }
 200
 201 //__________________
 202 void AliFemtoEventReaderESDChain::SetConstrained(const bool constrained)
 203 {
 204   // Select whether to read constrained or not constrained momentum
 205   fConstrained=constrained;
 206 }
 207 //__________________
 208 bool AliFemtoEventReaderESDChain::GetConstrained() const
 209 {
 210   // Check whether we read constrained or not constrained momentum
 211   return fConstrained;
 212 }
 213 //__________________
 214 void AliFemtoEventReaderESDChain::SetReadTPCInner(const bool readinner)
 215 {
 216   fReadInner=readinner;
 217 }
 218
 219 bool AliFemtoEventReaderESDChain::GetReadTPCInner() const
 220 {
 221   return fReadInner;
 222 }
 223
 224 //__________________
 225 void AliFemtoEventReaderESDChain::SetUseTPCOnly(const bool usetpconly)
 226 {
 227   fUseTPCOnly=usetpconly;
 228 }
 229
 230 bool AliFemtoEventReaderESDChain::GetUseTPCOnly() const
 231 {
 232   return fUseTPCOnly;
 233 }
 234
 235 void AliFemtoEventReaderESDChain::SetUsePhysicsSelection(const bool usephysics)
 236 {
 237   fUsePhysicsSel = usephysics;
 238   if (!fSelect) fSelect = new AliPhysicsSelection();
 239 }
 240
 241 void AliFemtoEventReaderESDChain::SetUseMultiplicity(EstEventMult aType)
 242 {
 243   fEstEventMult = aType;
 244 }
 245
 246 AliFemtoEvent* AliFemtoEventReaderESDChain::ReturnHbtEvent()
 247 {
 248   // Get the event, read all the relevant information
 249   // and fill the AliFemtoEvent class
 250   // Returns a valid AliFemtoEvent
 251   AliFemtoEvent *hbtEvent = 0;
 252   string tFriendFileName;
 253
 254   // Get the friend information
 255   if (Debug()>1) cout<<"starting to read event "<<fCurEvent<<endl;
 256   //  fEvent->SetESDfriend(fEventFriend);
 257   if(fEvent->GetAliESDOld())fEvent->CopyFromOldESD();
 258
 259   hbtEvent = new AliFemtoEvent;
 260
 261   if (fUsePhysicsSel) {
 262     hbtEvent->SetIsCollisionCandidate(fSelect->IsCollisionCandidate(fEvent));
 263     if (!(fSelect->IsCollisionCandidate(fEvent)))
 264       printf("Event not a collision candidate\n");
 265   }
 266   else
 267     hbtEvent->SetIsCollisionCandidate(kTRUE);
 268
 269   //setting basic things
 270   //  hbtEvent->SetEventNumber(fEvent->GetEventNumber());
 271   hbtEvent->SetRunNumber(fEvent->GetRunNumber());
 272   //hbtEvent->SetNumberOfTracks(fEvent->GetNumberOfTracks());
 273   hbtEvent->SetMagneticField(fEvent->GetMagneticField()*kilogauss);//to check if here is ok
 274   hbtEvent->SetZDCN1Energy(fEvent->GetZDCN1Energy());
 275   hbtEvent->SetZDCP1Energy(fEvent->GetZDCP1Energy());
 276   hbtEvent->SetZDCN2Energy(fEvent->GetZDCN2Energy());
 277   hbtEvent->SetZDCP2Energy(fEvent->GetZDCP2Energy());
 278   hbtEvent->SetZDCEMEnergy(fEvent->GetZDCEMEnergy());
 279   hbtEvent->SetZDCParticipants(fEvent->GetZDCParticipants());
 280   hbtEvent->SetTriggerMask(fEvent->GetTriggerMask());
 281   //  hbtEvent->SetTriggerCluster(fEvent->GetTriggerCluster());
 282
 283   if ((fEvent->IsTriggerClassFired("CINT1WU-B-NOPF-ALL")) ||
 284       (fEvent->IsTriggerClassFired("CINT1B-ABCE-NOPF-ALL")) ||
 285       (fEvent->IsTriggerClassFired("CINT1-B-NOPF-ALLNOTRD")) ||
 286       (fEvent->IsTriggerClassFired("CINT1-B-NOPF-FASTNOTRD")))
 287     hbtEvent->SetTriggerCluster(1);
 288   else if ((fEvent->IsTriggerClassFired("CSH1WU-B-NOPF-ALL")) ||
 289            (fEvent->IsTriggerClassFired("CSH1-B-NOPF-ALLNOTRD")))
 290     hbtEvent->SetTriggerCluster(2);
 291   else
 292     hbtEvent->SetTriggerCluster(0);
 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   Int_t spdetaonecount = 0;
 315
 316   for (int iter=0; iter<fEvent->GetMultiplicity()->GetNumberOfTracklets(); iter++)
 317     if (fabs(fEvent->GetMultiplicity()->GetEta(iter)) < 1.0)
 318       spdetaonecount++;
 319
 320   //  hbtEvent->SetSPDMult(fEvent->GetMultiplicity()->GetNumberOfTracklets());
 321   hbtEvent->SetSPDMult(spdetaonecount);
 322
 323   //starting to reading tracks
 324   int nofTracks=0;  //number of reconstructed tracks in event
 325   nofTracks=fEvent->GetNumberOfTracks();
 326   int realnofTracks=0;//number of track which we use ina analysis
 327
 328   //   // Clear the shared cluster list
 329   //   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
 330   //     fClusterPerPadrow[tPad]->clear();
 331   //   }
 332   //   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
 333   //     fSharedList[tPad]->clear();
 334   //   }
 335
 336
 337   //   for (int i=0;i<nofTracks;i++) {
 338   //     const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
 339
 340   //     list<Int_t>::iterator tClustIter;
 341
 342   //     Int_t tTrackIndices[AliESDfriendTrack::kMaxTPCcluster];
 343   //     Int_t tNClusters = esdtrack->GetTPCclusters(tTrackIndices);
 344   //     for (int tNcl=0; tNcl<AliESDfriendTrack::kMaxTPCcluster; tNcl++) {
 345   //       if (tTrackIndices[tNcl] >= 0) {
 346   //    tClustIter = find(fClusterPerPadrow[tNcl]->begin(), fClusterPerPadrow[tNcl]->end(), tTrackIndices[tNcl]);
 347   //    if (tClustIter == fClusterPerPadrow[tNcl]->end()) {
 348   //      fClusterPerPadrow[tNcl]->push_back(tTrackIndices[tNcl]);
 349   //    }
 350   //    else {
 351   //      fSharedList[tNcl]->push_back(tTrackIndices[tNcl]);
 352   //    }
 353   //       }
 354   //     }
 355
 356   //   }
 357
 358   int tNormMult = 0;
 359   int tNormMultPos = 0;
 360   int tNormMultNeg = 0;
 361
 362   Float_t tTotalPt = 0.0;
 363
 364   Float_t b[2];
 365   Float_t bCov[3];
 366
 367   Int_t tTracklet=0, tITSTPC=0, tITSPure=0;
 368
 369   fEvent->EstimateMultiplicity(tTracklet, tITSTPC, tITSPure, 1.2);
 370
 371   hbtEvent->SetMultiplicityEstimateITSTPC(tITSTPC);
 372   hbtEvent->SetMultiplicityEstimateTracklets(tTracklet);
 373   //  hbtEvent->SetMultiplicityEstimateITSPure(tITSPure);
 374   hbtEvent->SetMultiplicityEstimateITSPure(fEvent->GetMultiplicity()->GetNumberOfITSClusters(1));
 375
 376   for (int i=0;i<nofTracks;i++)
 377     {
 378       bool  tGoodMomentum=true; //flaga to chcek if we can read momentum of this track
 379
 380       const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
 381       //      const AliESDfriendTrack *tESDfriendTrack = esdtrack->GetFriendTrack();
 382
 383       if ((esdtrack->GetStatus() & AliESDtrack::kTPCrefit) &&
 384           (esdtrack->GetStatus() & AliESDtrack::kITSrefit)) {
 385         if (esdtrack->GetTPCNcls() > 70)
 386           if (esdtrack->GetTPCchi2()/esdtrack->GetTPCNcls() < 4.0) {
 387             if (TMath::Abs(esdtrack->Eta()) < 1.2) {
 388               esdtrack->GetImpactParameters(b,bCov);
 389               if ((b[0]<0.2) && (b[1] < 0.25)) {
 390                 tNormMult++;
 391                 tTotalPt += esdtrack->Pt();
 392               }
 393             }
 394           }
 395       }
 396       else if (esdtrack->GetStatus() & AliESDtrack::kTPCrefit) {
 397         if (esdtrack->GetTPCNcls() > 100)
 398           if (esdtrack->GetTPCchi2()/esdtrack->GetTPCNcls() < 4.0) {
 399             if (TMath::Abs(esdtrack->Eta()) < 1.2) {
 400               esdtrack->GetImpactParameters(b,bCov);
 401               if ((b[0]<2.4) && (b[1] < 3.2)) {
 402                 tNormMult++;
 403                 tTotalPt += esdtrack->Pt();
 404               }
 405             }
 406           }
 407       }
 408
 409       hbtEvent->SetZDCEMEnergy(tTotalPt);
 410       //       if (esdtrack->GetStatus() & AliESDtrack::kTPCrefit)
 411       //        if (esdtrack->GetTPCNcls() > 80)
 412       //          if (esdtrack->GetTPCchi2()/esdtrack->GetTPCNcls() < 6.0)
 413       //            if (esdtrack->GetConstrainedParam())
 414       //              if (fabs(esdtrack->GetConstrainedParam()->Eta()) < 0.5)
 415       //                if (esdtrack->GetConstrainedParam()->Pt() < 1.0) {
 416       //                  if (esdtrack->GetSign() > 0)
 417       //                    tNormMultPos++;
 418       //                  else if (esdtrack->GetSign() < 0)
 419       //                    tNormMultNeg--;
 420       //                }
 421
 422       // If reading ITS-only tracks, reject all with TPC
 423       if (fTrackType == kITSOnly) {
 424         if (esdtrack->GetStatus() & AliESDtrack::kTPCrefit) continue;
 425         if (!(esdtrack->GetStatus() & AliESDtrack::kITSrefit)) continue;
 426         if (esdtrack->GetStatus() & AliESDtrack::kTPCin) continue;
 427         UChar_t iclm = esdtrack->GetITSClusterMap();
 428         Int_t incls = 0;
 429         for (int iter=0; iter<6; iter++) if (iclm&(1<<iter)) incls++;
 430         if (incls<=3) {
 431           if (Debug()>1) cout << "Rejecting track with " << incls << " clusters" << endl;
 432           continue;
 433         }
 434       }
 435
 436       AliFemtoTrack* trackCopy = new AliFemtoTrack();
 437       trackCopy->SetCharge((short)esdtrack->GetSign());
 438
 439       //in aliroot we have AliPID
 440       //0-electron 1-muon 2-pion 3-kaon 4-proton 5-photon 6-pi0 7-neutron 8-kaon0 9-eleCon
 441       //we use only 5 first
 442       double esdpid[5];
 443       //       esdtrack->GetESDpid(esdpid);
 444       esdtrack->GetTPCpid(esdpid);
 445       trackCopy->SetPidProbElectron(esdpid[0]);
 446       trackCopy->SetPidProbMuon(esdpid[1]);
 447       trackCopy->SetPidProbPion(esdpid[2]);
 448       trackCopy->SetPidProbKaon(esdpid[3]);
 449       trackCopy->SetPidProbProton(esdpid[4]);
 450
 451       esdpid[0] = -100000.0;
 452       esdpid[1] = -100000.0;
 453       esdpid[2] = -100000.0;
 454       esdpid[3] = -100000.0;
 455       esdpid[4] = -100000.0;
 456
 457       double tTOF = 0.0;
 458
 459       if (esdtrack->GetStatus()&AliESDtrack::kTOFpid) {
 460         tTOF = esdtrack->GetTOFsignal();
 461         esdtrack->GetIntegratedTimes(esdpid);
 462       }
 463
 464       trackCopy->SetTofExpectedTimes(tTOF-esdpid[2], tTOF-esdpid[3], tTOF-esdpid[4]);
 465
 466       //////  TPC ////////////////////////////////////////////
 467
 468       float nsigmaTPCK=-1000.;
 469       float nsigmaTPCPi=-1000.;
 470       float nsigmaTPCP=-1000.;
 471
 472
 473       if ((fESDpid) && (esdtrack->IsOn(AliESDtrack::kTPCpid))){
 474         nsigmaTPCK = fESDpid->NumberOfSigmasTPC(esdtrack,AliPID::kKaon);
 475         nsigmaTPCPi = fESDpid->NumberOfSigmasTPC(esdtrack,AliPID::kPion);
 476         nsigmaTPCP = fESDpid->NumberOfSigmasTPC(esdtrack,AliPID::kProton);
 477
 478       }
 479       trackCopy->SetNSigmaTPCPi(nsigmaTPCPi);
 480       trackCopy->SetNSigmaTPCK(nsigmaTPCK);
 481       trackCopy->SetNSigmaTPCP(nsigmaTPCP);
 482
 483       ///// TOF ///////////////////////////////////////////////
 484
 485         float vp=-1000.;
 486         float nsigmaTOFPi=-1000.;
 487         float nsigmaTOFK=-1000.;
 488         float nsigmaTOFP=-1000.;
 489
 490         if ((esdtrack->GetStatus()&AliESDtrack::kTOFpid) &&
 491             (esdtrack->GetStatus()&AliESDtrack::kTOFout) &&
 492             (esdtrack->GetStatus()&AliESDtrack::kTIME) &&
 493             !(esdtrack->GetStatus()&AliESDtrack::kTOFmismatch))
 494           {
 495
 496             //if ((esdtrack->GetStatus()&AliESDtrack::kTOFpid) &&
 497             //(esdtrack->GetStatus()&AliESDtrack::kTOFout) &&
 498             //(esdtrack->GetStatus()&AliESDtrack::kTIME)){
 499             // collect info from ESDpid class
 500
 501             if ((fESDpid) && (esdtrack->IsOn(AliESDtrack::kTOFpid))) {
 502
 503
 504               double tZero = fESDpid->GetTOFResponse().GetStartTime(esdtrack->P());
 505
 506               nsigmaTOFPi = fESDpid->NumberOfSigmasTOF(esdtrack,AliPID::kPion,tZero);
 507               nsigmaTOFK = fESDpid->NumberOfSigmasTOF(esdtrack,AliPID::kKaon,tZero);
 508               nsigmaTOFP = fESDpid->NumberOfSigmasTOF(esdtrack,AliPID::kProton,tZero);
 509
 510               Double_t len=esdtrack->GetIntegratedLength();
 511               Double_t tof=esdtrack->GetTOFsignal();
 512               if(tof > 0.) vp=len/tof/0.03;
 513             }
 514           }
 515
 516         trackCopy->SetVTOF(vp);
 517         trackCopy->SetNSigmaTOFPi(nsigmaTOFPi);
 518         trackCopy->SetNSigmaTOFK(nsigmaTOFK);
 519         trackCopy->SetNSigmaTOFP(nsigmaTOFP);
 520
 521         double pxyz[3];
 522         double rxyz[3];
 523         double impact[2];
 524         double covimpact[3];
 525
 526         if (fUseTPCOnly) {
 527           if (!esdtrack->GetTPCInnerParam()) {
 528             delete trackCopy;
 529             continue;
 530           }
 531
 532
 533           AliExternalTrackParam *param = new AliExternalTrackParam(*esdtrack->GetTPCInnerParam());
 534           param->GetXYZ(rxyz);
 535           param->PropagateToDCA(fEvent->GetPrimaryVertexTPC(), (fEvent->GetMagneticField()), 10000, impact, covimpact);
 536           param->GetPxPyPz(pxyz);//reading noconstarined momentum
 537
 538           if (fReadInner == true) {
 539             AliFemtoModelHiddenInfo *tInfo = new AliFemtoModelHiddenInfo();
 540             tInfo->SetPDGPid(211);
 541             tInfo->SetTrueMomentum(pxyz[0], pxyz[1], pxyz[2]);
 542             tInfo->SetMass(0.13957);
 543             //    tInfo->SetEmissionPoint(rxyz[0], rxyz[1], rxyz[2], 0.0);
 544             //    tInfo->SetEmissionPoint(fV1[0], fV1[1], fV1[2], 0.0);
 545             tInfo->SetEmissionPoint(rxyz[0]-fV1[0], rxyz[1]-fV1[1], rxyz[2]-fV1[2], 0.0);
 546             trackCopy->SetHiddenInfo(tInfo);
 547           }
 548
 549           AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
 550           if (v.Mag() < 0.0001) {
 551             //  cout << "Found 0 momentum ???? " <<endl;
 552             delete trackCopy;
 553             continue;
 554           }
 555           trackCopy->SetP(v);//setting momentum
 556           trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
 557
 558           const AliFmThreeVectorD kP(pxyz[0],pxyz[1],pxyz[2]);
 559           const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
 560           //setting helix I do not if it is ok
 561           AliFmPhysicalHelixD helix(kP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge()));
 562           trackCopy->SetHelix(helix);
 563
 564           //some stuff which could be useful
 565           trackCopy->SetImpactD(impact[0]);
 566           trackCopy->SetImpactZ(impact[1]);
 567           trackCopy->SetCdd(covimpact[0]);
 568           trackCopy->SetCdz(covimpact[1]);
 569           trackCopy->SetCzz(covimpact[2]);
 570           trackCopy->SetSigmaToVertex(GetSigmaToVertex(impact, covimpact));
 571
 572           delete param;
 573         }
 574         else {
 575           if (fReadInner == true) {
 576
 577             if (esdtrack->GetTPCInnerParam()) {
 578               AliExternalTrackParam *param = new AliExternalTrackParam(*esdtrack->GetTPCInnerParam());
 579               param->GetXYZ(rxyz);
 580               //            param->PropagateToDCA(fEvent->GetPrimaryVertex(), (fEvent->GetMagneticField()), 10000);
 581               param->GetPxPyPz(pxyz);//reading noconstarined momentum
 582               delete param;
 583
 584               AliFemtoModelHiddenInfo *tInfo = new AliFemtoModelHiddenInfo();
 585               tInfo->SetPDGPid(211);
 586               tInfo->SetTrueMomentum(pxyz[0], pxyz[1], pxyz[2]);
 587               tInfo->SetMass(0.13957);
 588               //            tInfo->SetEmissionPoint(rxyz[0], rxyz[1], rxyz[2], 0.0);
 589               //tInfo->SetEmissionPoint(fV1[0], fV1[1], fV1[2], 0.0);
 590               tInfo->SetEmissionPoint(rxyz[0]-fV1[0], rxyz[1]-fV1[1], rxyz[2]-fV1[2], 0.0);
 591               trackCopy->SetHiddenInfo(tInfo);
 592             }
 593           }
 594
 595           if (fTrackType == kGlobal) {
 596             if (fConstrained==true)
 597               tGoodMomentum=esdtrack->GetConstrainedPxPyPz(pxyz); //reading constrained momentum
 598             else
 599               tGoodMomentum=esdtrack->GetPxPyPz(pxyz);//reading noconstarined momentum
 600           }
 601           else if (fTrackType == kTPCOnly) {
 602             if (esdtrack->GetTPCInnerParam())
 603               esdtrack->GetTPCInnerParam()->GetPxPyPz(pxyz);
 604             else {
 605               delete trackCopy;
 606               continue;
 607             }
 608           }
 609           else if (fTrackType == kITSOnly) {
 610             if (fConstrained==true)
 611               tGoodMomentum=esdtrack->GetConstrainedPxPyPz(pxyz); //reading constrained momentum
 612             else
 613               tGoodMomentum=esdtrack->GetPxPyPz(pxyz);//reading noconstarined momentum
 614           }
 615
 616
 617           AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
 618           if (v.Mag() < 0.0001) {
 619             //  cout << "Found 0 momentum ???? " <<endl;
 620             delete trackCopy;
 621             continue;
 622           }
 623           trackCopy->SetP(v);//setting momentum
 624           trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
 625           const AliFmThreeVectorD kP(pxyz[0],pxyz[1],pxyz[2]);
 626           const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
 627           //setting helix I do not if it is ok
 628           AliFmPhysicalHelixD helix(kP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge()));
 629           trackCopy->SetHelix(helix);
 630
 631           //some stuff which could be useful
 632           float imp[2];
 633           float cim[3];
 634           // if (fTrackType == kTPCOnly) {
 635           //   esdtrack->GetTPCInnerParam()->GetImpactParameters(imp,cim);
 636           // }
 637           // else {
 638           esdtrack->GetImpactParameters(imp,cim);
 639           // }
 640
 641           impact[0] = imp[0];
 642           impact[1] = imp[1];
 643           covimpact[0] = cim[0];
 644           covimpact[1] = cim[1];
 645           covimpact[2] = cim[2];
 646
 647           trackCopy->SetImpactD(impact[0]);
 648           trackCopy->SetImpactZ(impact[1]);
 649           trackCopy->SetCdd(covimpact[0]);
 650           trackCopy->SetCdz(covimpact[1]);
 651           trackCopy->SetCzz(covimpact[2]);
 652           trackCopy->SetSigmaToVertex(GetSigmaToVertex(impact,covimpact));
 653         }
 654
 655         trackCopy->SetTrackId(esdtrack->GetID());
 656         trackCopy->SetFlags(esdtrack->GetStatus());
 657         trackCopy->SetLabel(esdtrack->GetLabel());
 658
 659         trackCopy->SetITSchi2(esdtrack->GetITSchi2());
 660         if (esdtrack->GetITSFakeFlag())
 661           trackCopy->SetITSncls(-esdtrack->GetNcls(0));
 662         else
 663           trackCopy->SetITSncls(esdtrack->GetNcls(0));
 664         trackCopy->SetTPCchi2(esdtrack->GetTPCchi2());
 665         trackCopy->SetTPCncls(esdtrack->GetTPCNcls());
 666         trackCopy->SetTPCnclsF(esdtrack->GetTPCNclsF());
 667         trackCopy->SetTPCsignal(esdtrack->GetTPCsignal());
 668         trackCopy->SetTPCsignalN((short)esdtrack->GetTPCsignalN()); //due to bug in aliesdtrack class
 669         trackCopy->SetTPCsignalS(esdtrack->GetTPCsignalSigma());
 670
 671         trackCopy->SetTPCClusterMap(esdtrack->GetTPCClusterMap());
 672         trackCopy->SetTPCSharedMap(esdtrack->GetTPCSharedMap());
 673
 674         double xtpc[3];
 675         esdtrack->GetInnerXYZ(xtpc);
 676         xtpc[2] -= fV1[2];
 677         trackCopy->SetNominalTPCEntrancePoint(xtpc);
 678
 679         esdtrack->GetOuterXYZ(xtpc);
 680         xtpc[2] -= fV1[2];
 681         trackCopy->SetNominalTPCExitPoint(xtpc);
 682
 683         int indexes[3];
 684         for (int ik=0; ik<3; ik++) {
 685           indexes[ik] = esdtrack->GetKinkIndex(ik);
 686         }
 687         trackCopy->SetKinkIndexes(indexes);
 688         //decision if we want this track
 689         //if we using diffrent labels we want that this label was use for first time
 690         //if we use hidden info we want to have match between sim data and ESD
 691         if (tGoodMomentum==true)
 692           {
 693             hbtEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
 694             realnofTracks++;//real number of tracks
 695             //    delete trackCopy;
 696           }
 697         else
 698           {
 699             delete  trackCopy;
 700           }
 701
 702     }
 703
 704   hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event
 705
 706   AliCentrality *cent = fEvent->GetCentrality();
 707   if (cent) {
 708     hbtEvent->SetCentralityV0(cent->GetCentralityPercentile("V0M"));
 709     //    hbtEvent->SetCentralityFMD(cent->GetCentralityPercentile("FMD"));
 710     hbtEvent->SetCentralitySPD1(cent->GetCentralityPercentile("CL1"));
 711     //    hbtEvent->SetCentralityTrk(cent->GetCentralityPercentile("TRK"));
 712
 713     if (Debug()>1) printf("  FemtoReader Got Event with %f %f %f %f\n", cent->GetCentralityPercentile("V0M"), 0.0, cent->GetCentralityPercentile("CL1"), 0.0);
 714   }
 715
 716   if (fEstEventMult == kGlobalCount)
 717     hbtEvent->SetNormalizedMult(tNormMult);
 718   else if (fEstEventMult == kTracklet)
 719     hbtEvent->SetNormalizedMult(tTracklet);
 720   else if (fEstEventMult == kITSTPC)
 721     hbtEvent->SetNormalizedMult(tITSTPC);
 722   else if (fEstEventMult == kITSPure)
 723     hbtEvent->SetNormalizedMult(tITSPure);
 724   else if (fEstEventMult == kSPDLayer1)
 725     hbtEvent->SetNormalizedMult(fEvent->GetMultiplicity()->GetNumberOfITSClusters(1));
 726   else if (fEstEventMult == kV0Centrality) {
 727     // centrality between 0 (central) and 1 (very peripheral)
 728
 729     if (cent) {
 730       if (cent->GetCentralityPercentile("V0M") < 0.00001)
 731         hbtEvent->SetNormalizedMult(-1);
 732       else
 733         hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("V0M")));
 734       if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(),
 735                              10.0*cent->GetCentralityPercentile("V0M"), lrint(10.0*cent->GetCentralityPercentile("V0M")));
 736     }
 737   }
 738
 739   if (tNormMultPos > tNormMultNeg)
 740     hbtEvent->SetZDCParticipants(tNormMultPos);
 741   else
 742     hbtEvent->SetZDCParticipants(tNormMultNeg);
 743
 744   AliEventplane* ep = fEvent->GetEventplane();
 745   if (ep) {
 746     hbtEvent->SetEP(ep);
 747     hbtEvent->SetReactionPlaneAngle(ep->GetEventplane("Q"));
 748   }
 749
 750   fCurEvent++;
 751   if (Debug()>1) cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
 752
 753   return hbtEvent;
 754 }
 755 //___________________
 756 void AliFemtoEventReaderESDChain::SetESDSource(AliESDEvent *aESD)
 757 {
 758   // The chain loads the ESD for us
 759   // You must provide the address where it can be found
 760   fEvent = aESD;
 761 }
 762 //___________________
 763 // void AliFemtoEventReaderESDChain::SetESDfriendSource(AliESDfriend *aFriend)
 764 // {
 765 //   // We need the ESD tree to obtain
 766 //   // information about the friend file location
 767 //   fEventFriend = aFriend;
 768 // }
 769
 770 //____________________________________________________________________
 771 Float_t AliFemtoEventReaderESDChain::GetSigmaToVertex(double *impact, double *covar)
 772 {
 773   // Calculates the number of sigma to the vertex.
 774
 775   Float_t b[2];
 776   Float_t bRes[2];
 777   Float_t bCov[3];
 778
 779   b[0] = impact[0];
 780   b[1] = impact[1];
 781   bCov[0] = covar[0];
 782   bCov[1] = covar[1];
 783   bCov[2] = covar[2];
 784
 785   bRes[0] = TMath::Sqrt(bCov[0]);
 786   bRes[1] = TMath::Sqrt(bCov[2]);
 787
 788   // -----------------------------------
 789   // How to get to a n-sigma cut?
 790   //
 791   // The accumulated statistics from 0 to d is
 792   //
 793   // ->  Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
 794   // ->  1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
 795   //
 796   // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-x**2)/2)
 797   // Can this be expressed in a different way?
 798
 799   if (bRes[0] == 0 || bRes[1] ==0)
 800     return -1;
 801
 802   Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
 803
 804   // stupid rounding problem screws up everything:
 805   // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
 806   if (TMath::Exp(-d * d / 2) < 1e-10)
 807     return 1000;
 808
 809   d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
 810   return d;
 811 }
 812
 813 void AliFemtoEventReaderESDChain::SetReadTrackType(ReadTrackType aType)
 814 {
 815   fTrackType = aType;
 816 }
 817
 818 //trigger
 819 void AliFemtoEventReaderESDChain::SetEventTrigger(UInt_t eventtrig)
 820 {
 821   fEventTrig = eventtrig;
 822 }
 823
 824
 825
 826
 827
 828