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 //____________________________
  36 AliFemtoEventReaderESDChain::AliFemtoEventReaderESDChain():
  37   fFileName(" "),
  38   fConstrained(true),
  39   fReadInner(false),
  40   fUseTPCOnly(false),
  41   fNumberofEvent(0),
  42   fCurEvent(0),
  43   fCurFile(0),
  44   fEvent(0x0),
  45   fTrackType(kGlobal),
  46   fEstEventMult(kReferenceITSTPC),
  47   fEventTrig(AliVEvent::kMB), //trigger
  48   fESDpid(0),
  49   fIsPidOwner(0),
  50   fReadV0(0),
  51   fMagFieldSign(0),
  52   fpA2013(kFALSE)
  53 {
  54   //constructor with 0 parameters , look at default settings
  55   //   fClusterPerPadrow = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
  56   //   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
  57   //     fClusterPerPadrow[tPad] = new list<Int_t>();
  58   //   }
  59   //   fSharedList = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
  60   //   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
  61   //     fSharedList[tPad] = new list<Int_t>();
  62   //   }
  63 }
  64
  65 //__________________
  66 AliFemtoEventReaderESDChain::AliFemtoEventReaderESDChain(const AliFemtoEventReaderESDChain& aReader):
  67   AliFemtoEventReader(aReader),
  68   fFileName(" "),
  69   fConstrained(true),
  70   fReadInner(false),
  71   fUseTPCOnly(false),
  72   fNumberofEvent(0),
  73   fCurEvent(0),
  74   fCurFile(0),
  75   fEvent(0x0),
  76   fTrackType(kGlobal),
  77   fEstEventMult(kReferenceITSTPC),
  78   fEventTrig(AliVEvent::kMB), //trigger
  79   fESDpid(0),
  80   fIsPidOwner(0),
  81   fReadV0(0),
  82   fMagFieldSign(0),
  83   fpA2013(kFALSE)
  84 {
  85   // Copy constructor
  86   fConstrained = aReader.fConstrained;
  87   fReadInner = aReader.fReadInner;
  88   fUseTPCOnly = aReader.fUseTPCOnly;
  89   fNumberofEvent = aReader.fNumberofEvent;
  90   fCurEvent = aReader.fCurEvent;
  91   fCurFile = aReader.fCurFile;
  92   //  fEvent = new AliESD(*aReader.fEvent);
  93   fEvent = new AliESDEvent();
  94   fTrackType = aReader.fTrackType;
  95   fEstEventMult = aReader.fEstEventMult;
  96   fEventTrig = aReader.fEventTrig; //trigger
  97   fReadV0 = aReader.fReadV0;
  98   fMagFieldSign = aReader.fMagFieldSign;
  99   fpA2013 = aReader.fpA2013;
 100
 101   //   fEventFriend = aReader.fEventFriend;
 102   //   fClusterPerPadrow = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
 103   //   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
 104   //     fClusterPerPadrow[tPad] = new list<Int_t>();
 105   //     list<Int_t>::iterator iter;
 106   //     for (iter=aReader.fClusterPerPadrow[tPad]->begin(); iter!=aReader.fClusterPerPadrow[tPad]->end(); iter++) {
 107   //       fClusterPerPadrow[tPad]->push_back(*iter);
 108   //     }
 109   //   }
 110   //   fSharedList = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
 111   //   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
 112   //     fSharedList[tPad] = new list<Int_t>();
 113   //     list<Int_t>::iterator iter;
 114   //     for (iter=aReader.fSharedList[tPad]->begin(); iter!=aReader.fSharedList[tPad]->end(); iter++) {
 115   //       fSharedList[tPad]->push_back(*iter);
 116   //     }
 117   //   }
 118 }
 119 //__________________
 120 AliFemtoEventReaderESDChain::~AliFemtoEventReaderESDChain()
 121 {
 122   //Destructor
 123   delete fEvent;
 124
 125   //   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
 126   //     fClusterPerPadrow[tPad]->clear();
 127   //     delete fClusterPerPadrow[tPad];
 128   //   }
 129   //   delete [] fClusterPerPadrow;
 130   //   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
 131   //     fSharedList[tPad]->clear();
 132   //     delete fSharedList[tPad];
 133   //   }
 134   //   delete [] fSharedList;
 135
 136 }
 137
 138 //__________________
 139 AliFemtoEventReaderESDChain& AliFemtoEventReaderESDChain::operator=(const AliFemtoEventReaderESDChain& aReader)
 140 {
 141   // Assignment operator
 142   if (this == &aReader)
 143     return *this;
 144
 145   fConstrained = aReader.fConstrained;
 146   fReadInner = aReader.fReadInner;
 147   fUseTPCOnly = aReader.fUseTPCOnly;
 148   fNumberofEvent = aReader.fNumberofEvent;
 149   fCurEvent = aReader.fCurEvent;
 150   fCurFile = aReader.fCurFile;
 151   if (fEvent) delete fEvent;
 152   fEvent = new AliESDEvent();
 153   fTrackType = aReader.fTrackType;
 154   fEstEventMult = aReader.fEstEventMult;
 155
 156   fReadV0 = aReader.fReadV0;
 157   fMagFieldSign = aReader.fMagFieldSign;
 158   fpA2013 = aReader.fpA2013;
 159   //  fEventFriend = aReader.fEventFriend;
 160
 161   //   if (fClusterPerPadrow) {
 162   //     for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
 163   //       fClusterPerPadrow[tPad]->clear();
 164   //       delete fClusterPerPadrow[tPad];
 165   //     }
 166   //     delete [] fClusterPerPadrow;
 167   //   }
 168
 169   //   if (fSharedList) {
 170   //     for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
 171   //       fSharedList[tPad]->clear();
 172   //       delete fSharedList[tPad];
 173   //     }
 174   //     delete [] fSharedList;
 175   //   }
 176
 177   //   fClusterPerPadrow = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
 178   //   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
 179   //     fClusterPerPadrow[tPad] = new list<Int_t>();
 180   //     list<Int_t>::iterator iter;
 181   //     for (iter=aReader.fClusterPerPadrow[tPad]->begin(); iter!=aReader.fClusterPerPadrow[tPad]->end(); iter++) {
 182   //       fClusterPerPadrow[tPad]->push_back(*iter);
 183   //     }
 184   //   }
 185   //   fSharedList = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
 186   //   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
 187   //     fSharedList[tPad] = new list<Int_t>();
 188   //     list<Int_t>::iterator iter;
 189   //     for (iter=aReader.fSharedList[tPad]->begin(); iter!=aReader.fSharedList[tPad]->end(); iter++) {
 190   //       fSharedList[tPad]->push_back(*iter);
 191   //     }
 192   //   }
 193
 194   return *this;
 195 }
 196 //__________________
 197 // Simple report
 198 AliFemtoString AliFemtoEventReaderESDChain::Report()
 199 {
 200   AliFemtoString temp = "\n This is the AliFemtoEventReaderESDChain\n";
 201   return temp;
 202 }
 203
 204 //__________________
 205 void AliFemtoEventReaderESDChain::SetConstrained(const bool constrained)
 206 {
 207   // Select whether to read constrained or not constrained momentum
 208   fConstrained=constrained;
 209 }
 210 //__________________
 211 bool AliFemtoEventReaderESDChain::GetConstrained() const
 212 {
 213   // Check whether we read constrained or not constrained momentum
 214   return fConstrained;
 215 }
 216 //__________________
 217 void AliFemtoEventReaderESDChain::SetReadTPCInner(const bool readinner)
 218 {
 219   fReadInner=readinner;
 220 }
 221
 222 bool AliFemtoEventReaderESDChain::GetReadTPCInner() const
 223 {
 224   return fReadInner;
 225 }
 226
 227 //__________________
 228 void AliFemtoEventReaderESDChain::SetUseTPCOnly(const bool usetpconly)
 229 {
 230   fUseTPCOnly=usetpconly;
 231 }
 232
 233 bool AliFemtoEventReaderESDChain::GetUseTPCOnly() const
 234 {
 235   return fUseTPCOnly;
 236 }
 237
 238 void AliFemtoEventReaderESDChain::SetUseMultiplicity(EstEventMult aType)
 239 {
 240   fEstEventMult = aType;
 241 }
 242
 243 AliFemtoEvent* AliFemtoEventReaderESDChain::ReturnHbtEvent()
 244 {
 245   // Get the event, read all the relevant information
 246   // and fill the AliFemtoEvent class
 247   // Returns a valid AliFemtoEvent
 248   AliFemtoEvent *hbtEvent = 0;
 249
 250
 251   hbtEvent = new AliFemtoEvent;
 252
 253   CopyESDtoFemtoEvent(hbtEvent);
 254
 255   fCurEvent++;
 256
 257
 258   return hbtEvent;
 259 }
 260
 261 void AliFemtoEventReaderESDChain::CopyESDtoFemtoEvent(AliFemtoEvent *hbtEvent)
 262 {
 263
 264
 265  //string tFriendFileName;
 266   // Get the friend information
 267   if (Debug()>1) cout<<"starting to read event "<<fCurEvent<<endl;
 268   //  fEvent->SetESDfriend(fEventFriend);
 269   if(fEvent->GetAliESDOld())fEvent->CopyFromOldESD();
 270
 271
 272
 273    //setting basic things
 274   //  hbtEvent->SetEventNumber(fEvent->GetEventNumber());
 275   hbtEvent->SetRunNumber(fEvent->GetRunNumber());
 276   //hbtEvent->SetNumberOfTracks(fEvent->GetNumberOfTracks());
 277   hbtEvent->SetMagneticField(fEvent->GetMagneticField()*kilogauss);//to check if here is ok
 278   hbtEvent->SetZDCN1Energy(fEvent->GetZDCN1Energy());
 279   hbtEvent->SetZDCP1Energy(fEvent->GetZDCP1Energy());
 280   hbtEvent->SetZDCN2Energy(fEvent->GetZDCN2Energy());
 281   hbtEvent->SetZDCP2Energy(fEvent->GetZDCP2Energy());
 282   hbtEvent->SetZDCEMEnergy(fEvent->GetZDCEMEnergy());
 283   hbtEvent->SetZDCParticipants(fEvent->GetZDCParticipants());
 284   hbtEvent->SetTriggerMask(fEvent->GetTriggerMask());
 285   //  hbtEvent->SetTriggerCluster(fEvent->GetTriggerCluster());
 286
 287   if ((fEvent->IsTriggerClassFired("CINT1WU-B-NOPF-ALL")) ||
 288       (fEvent->IsTriggerClassFired("CINT1B-ABCE-NOPF-ALL")) ||
 289       (fEvent->IsTriggerClassFired("CINT1-B-NOPF-ALLNOTRD")) ||
 290       (fEvent->IsTriggerClassFired("CINT1-B-NOPF-FASTNOTRD")))
 291     hbtEvent->SetTriggerCluster(1);
 292   else if ((fEvent->IsTriggerClassFired("CSH1WU-B-NOPF-ALL")) ||
 293            (fEvent->IsTriggerClassFired("CSH1-B-NOPF-ALLNOTRD")))
 294     hbtEvent->SetTriggerCluster(2);
 295   else
 296     hbtEvent->SetTriggerCluster(0);
 297
 298   //Vertex
 299   double fV1[3];
 300   double fVCov[6];
 301
 302   if(fpA2013)
 303     {
 304       const AliESDVertex* trkVtx = fEvent->GetPrimaryVertex();
 305       if (!trkVtx || trkVtx->GetNContributors()<=0) return;
 306       TString vtxTtl = trkVtx->GetTitle();
 307       if (!vtxTtl.Contains("VertexerTracks")) return;
 308       Float_t zvtx = trkVtx->GetZ();
 309       const AliESDVertex* spdVtx = fEvent->GetVertex();
 310       if (spdVtx->GetNContributors()<=0) return;
 311       TString vtxTyp = spdVtx->GetTitle();
 312       Double_t cov[6]={0};
 313       spdVtx->GetCovarianceMatrix(cov);
 314       Double_t zRes = TMath::Sqrt(cov[5]);
 315       if (vtxTyp.Contains("vertexer:Z") && (zRes>0.25)) return;
 316       if (TMath::Abs(spdVtx->GetZ() - trkVtx->GetZ())>0.5) return;
 317
 318       if (TMath::Abs(zvtx) > 10) return;
 319     }
 320   if (fUseTPCOnly) {
 321     fEvent->GetPrimaryVertexTPC()->GetXYZ(fV1);
 322     fEvent->GetPrimaryVertexTPC()->GetCovMatrix(fVCov);
 323     if (!fEvent->GetPrimaryVertexTPC()->GetStatus())
 324       fVCov[4] = -1001.0;
 325   }
 326   else {
 327     fEvent->GetPrimaryVertex()->GetXYZ(fV1);
 328     fEvent->GetPrimaryVertex()->GetCovMatrix(fVCov);
 329     if (!fEvent->GetPrimaryVertex()->GetStatus())
 330       fVCov[4] = -1001.0;
 331   }
 332
 333   AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]);
 334   hbtEvent->SetPrimVertPos(vertex);
 335   hbtEvent->SetPrimVertCov(fVCov);
 336
 337   Int_t spdetaonecount = 0;
 338
 339   for (int iter=0; iter<fEvent->GetMultiplicity()->GetNumberOfTracklets(); iter++)
 340     if (fabs(fEvent->GetMultiplicity()->GetEta(iter)) < 1.0)
 341       spdetaonecount++;
 342
 343   //  hbtEvent->SetSPDMult(fEvent->GetMultiplicity()->GetNumberOfTracklets());
 344   hbtEvent->SetSPDMult(spdetaonecount);
 345
 346   //starting to reading tracks
 347   int nofTracks=0;  //number of reconstructed tracks in event
 348   nofTracks=fEvent->GetNumberOfTracks();
 349   int realnofTracks=0;//number of track which we use ina analysis
 350
 351   //   // Clear the shared cluster list
 352   //   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
 353   //     fClusterPerPadrow[tPad]->clear();
 354   //   }
 355   //   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
 356   //     fSharedList[tPad]->clear();
 357   //   }
 358
 359
 360   //   for (int i=0;i<nofTracks;i++) {
 361   //     const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
 362
 363   //     list<Int_t>::iterator tClustIter;
 364
 365   //     Int_t tTrackIndices[AliESDfriendTrack::kMaxTPCcluster];
 366   //     Int_t tNClusters = esdtrack->GetTPCclusters(tTrackIndices);
 367   //     for (int tNcl=0; tNcl<AliESDfriendTrack::kMaxTPCcluster; tNcl++) {
 368   //       if (tTrackIndices[tNcl] >= 0) {
 369   //    tClustIter = find(fClusterPerPadrow[tNcl]->begin(), fClusterPerPadrow[tNcl]->end(), tTrackIndices[tNcl]);
 370   //    if (tClustIter == fClusterPerPadrow[tNcl]->end()) {
 371   //      fClusterPerPadrow[tNcl]->push_back(tTrackIndices[tNcl]);
 372   //    }
 373   //    else {
 374   //      fSharedList[tNcl]->push_back(tTrackIndices[tNcl]);
 375   //    }
 376   //       }
 377   //     }
 378
 379   //   }
 380
 381   int tNormMult = 0;
 382   int tNormMultPos = 0;
 383   int tNormMultNeg = 0;
 384
 385   Float_t tTotalPt = 0.0;
 386
 387   Float_t b[2];
 388   Float_t bCov[3];
 389
 390   Int_t tTracklet=0, tITSTPC=0;
 391
 392   //W-AliESDEvent::EstimateMultiplicity: This obsolete method will be eliminated soon. Use AliESDtrackCuts::GetReferenceMultiplicity
 393   //fEvent->EstimateMultiplicity(tTracklet, tITSTPC, tITSPure, 1.2);
 394
 395   hbtEvent->SetMultiplicityEstimateITSTPC(tITSTPC);
 396   hbtEvent->SetMultiplicityEstimateTracklets(tTracklet);
 397   //  hbtEvent->SetMultiplicityEstimateITSPure(tITSPure);
 398   hbtEvent->SetMultiplicityEstimateITSPure(fEvent->GetMultiplicity()->GetNumberOfITSClusters(1));
 399
 400   for (int i=0;i<nofTracks;i++)
 401     {
 402       bool  tGoodMomentum=true; //flaga to chcek if we can read momentum of this track
 403
 404       const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
 405       //      const AliESDfriendTrack *tESDfriendTrack = esdtrack->GetFriendTrack();
 406
 407       if ((esdtrack->GetStatus() & AliESDtrack::kTPCrefit) &&
 408           (esdtrack->GetStatus() & AliESDtrack::kITSrefit)) {
 409         if (esdtrack->GetTPCNcls() > 70)
 410           if (esdtrack->GetTPCchi2()/esdtrack->GetTPCNcls() < 4.0) {
 411             if (esdtrack->Pt() > 0.15 && esdtrack->Pt() < 20)
 412               if (TMath::Abs(esdtrack->Eta()) < 0.8) {
 413                 esdtrack->GetImpactParameters(b,bCov);
 414                 if ((b[0]<0.2) && (b[1] < 2.0)) {
 415                   tNormMult++;
 416                   tTotalPt += esdtrack->Pt();
 417                 }
 418               }
 419           }
 420       }
 421
 422       hbtEvent->SetZDCEMEnergy(tTotalPt);
 423       //       if (esdtrack->GetStatus() & AliESDtrack::kTPCrefit)
 424       //        if (esdtrack->GetTPCNcls() > 80)
 425       //          if (esdtrack->GetTPCchi2()/esdtrack->GetTPCNcls() < 6.0)
 426       //            if (esdtrack->GetConstrainedParam())
 427       //              if (fabs(esdtrack->GetConstrainedParam()->Eta()) < 0.5)
 428       //                if (esdtrack->GetConstrainedParam()->Pt() < 1.0) {
 429       //                  if (esdtrack->GetSign() > 0)
 430       //                    tNormMultPos++;
 431       //                  else if (esdtrack->GetSign() < 0)
 432       //                    tNormMultNeg--;
 433       //                }
 434
 435       // If reading ITS-only tracks, reject all with TPC
 436       if (fTrackType == kITSOnly) {
 437         if (esdtrack->GetStatus() & AliESDtrack::kTPCrefit) continue;
 438         if (!(esdtrack->GetStatus() & AliESDtrack::kITSrefit)) continue;
 439         if (esdtrack->GetStatus() & AliESDtrack::kTPCin) continue;
 440         UChar_t iclm = esdtrack->GetITSClusterMap();
 441         Int_t incls = 0;
 442         for (int iter=0; iter<6; iter++) if (iclm&(1<<iter)) incls++;
 443         if (incls<=3) {
 444           if (Debug()>1) cout << "Rejecting track with " << incls << " clusters" << endl;
 445           continue;
 446         }
 447       }
 448
 449       AliFemtoTrack* trackCopy = new AliFemtoTrack();
 450       trackCopy->SetCharge((short)esdtrack->GetSign());
 451
 452       //in aliroot we have AliPID
 453       //0-electron 1-muon 2-pion 3-kaon 4-proton 5-photon 6-pi0 7-neutron 8-kaon0 9-eleCon
 454       //we use only 5 first
 455       double esdpid[5];
 456       //       esdtrack->GetESDpid(esdpid);
 457       esdtrack->GetTPCpid(esdpid);
 458       trackCopy->SetPidProbElectron(esdpid[0]);
 459       trackCopy->SetPidProbMuon(esdpid[1]);
 460       trackCopy->SetPidProbPion(esdpid[2]);
 461       trackCopy->SetPidProbKaon(esdpid[3]);
 462       trackCopy->SetPidProbProton(esdpid[4]);
 463
 464       esdpid[0] = -100000.0;
 465       esdpid[1] = -100000.0;
 466       esdpid[2] = -100000.0;
 467       esdpid[3] = -100000.0;
 468       esdpid[4] = -100000.0;
 469
 470       double tTOF = 0.0;
 471
 472       if (esdtrack->GetStatus()&AliESDtrack::kTOFpid) {
 473         tTOF = esdtrack->GetTOFsignal();
 474         esdtrack->GetIntegratedTimes(esdpid);
 475       }
 476
 477       trackCopy->SetTofExpectedTimes(tTOF-esdpid[2], tTOF-esdpid[3], tTOF-esdpid[4]);
 478
 479       //////  TPC ////////////////////////////////////////////
 480
 481       float nsigmaTPCK=-1000.;
 482       float nsigmaTPCPi=-1000.;
 483       float nsigmaTPCP=-1000.;
 484
 485
 486       if ((fESDpid) && (esdtrack->IsOn(AliESDtrack::kTPCpid))){
 487         nsigmaTPCK = fESDpid->NumberOfSigmasTPC(esdtrack,AliPID::kKaon);
 488         nsigmaTPCPi = fESDpid->NumberOfSigmasTPC(esdtrack,AliPID::kPion);
 489         nsigmaTPCP = fESDpid->NumberOfSigmasTPC(esdtrack,AliPID::kProton);
 490
 491       }
 492       trackCopy->SetNSigmaTPCPi(nsigmaTPCPi);
 493       trackCopy->SetNSigmaTPCK(nsigmaTPCK);
 494       trackCopy->SetNSigmaTPCP(nsigmaTPCP);
 495
 496       ///// TOF ///////////////////////////////////////////////
 497
 498         float vp=-1000.;
 499         float nsigmaTOFPi=-1000.;
 500         float nsigmaTOFK=-1000.;
 501         float nsigmaTOFP=-1000.;
 502
 503         if ((esdtrack->GetStatus()&AliESDtrack::kTOFpid) &&
 504             (esdtrack->GetStatus()&AliESDtrack::kTOFout) &&
 505             (esdtrack->GetStatus()&AliESDtrack::kTIME))
 506           {
 507
 508             //if ((esdtrack->GetStatus()&AliESDtrack::kTOFpid) &&
 509             //(esdtrack->GetStatus()&AliESDtrack::kTOFout) &&
 510             //(esdtrack->GetStatus()&AliESDtrack::kTIME)){
 511             // collect info from ESDpid class
 512
 513             if ((fESDpid) && (esdtrack->IsOn(AliESDtrack::kTOFpid))) {
 514
 515
 516               double tZero = fESDpid->GetTOFResponse().GetStartTime(esdtrack->P());
 517
 518               nsigmaTOFPi = fESDpid->NumberOfSigmasTOF(esdtrack,AliPID::kPion,tZero);
 519               nsigmaTOFK = fESDpid->NumberOfSigmasTOF(esdtrack,AliPID::kKaon,tZero);
 520               nsigmaTOFP = fESDpid->NumberOfSigmasTOF(esdtrack,AliPID::kProton,tZero);
 521
 522               Double_t len=esdtrack->GetIntegratedLength();
 523               Double_t tof=esdtrack->GetTOFsignal();
 524               if(tof > 0.) vp=len/tof/0.03;
 525             }
 526           }
 527
 528         trackCopy->SetVTOF(vp);
 529         trackCopy->SetNSigmaTOFPi(nsigmaTOFPi);
 530         trackCopy->SetNSigmaTOFK(nsigmaTOFK);
 531         trackCopy->SetNSigmaTOFP(nsigmaTOFP);
 532
 533         double pxyz[3];
 534         double rxyz[3];
 535         double impact[2];
 536         double covimpact[3];
 537
 538         if (fUseTPCOnly) {
 539           if (!esdtrack->GetTPCInnerParam()) {
 540             delete trackCopy;
 541             continue;
 542           }
 543
 544
 545           AliExternalTrackParam *param = new AliExternalTrackParam(*esdtrack->GetTPCInnerParam());
 546           param->GetXYZ(rxyz);
 547           param->PropagateToDCA(fEvent->GetPrimaryVertexTPC(), (fEvent->GetMagneticField()), 10000, impact, covimpact);
 548           param->GetPxPyPz(pxyz);//reading noconstarined momentum
 549
 550           if (fReadInner == true) {
 551             AliFemtoModelHiddenInfo *tInfo = new AliFemtoModelHiddenInfo();
 552             tInfo->SetPDGPid(211);
 553             tInfo->SetTrueMomentum(pxyz[0], pxyz[1], pxyz[2]);
 554             tInfo->SetMass(0.13957);
 555             //    tInfo->SetEmissionPoint(rxyz[0], rxyz[1], rxyz[2], 0.0);
 556             //    tInfo->SetEmissionPoint(fV1[0], fV1[1], fV1[2], 0.0);
 557             tInfo->SetEmissionPoint(rxyz[0]-fV1[0], rxyz[1]-fV1[1], rxyz[2]-fV1[2], 0.0);
 558             trackCopy->SetHiddenInfo(tInfo);
 559           }
 560
 561           AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
 562           if (v.Mag() < 0.0001) {
 563             //  cout << "Found 0 momentum ???? " <<endl;
 564             delete trackCopy;
 565             continue;
 566           }
 567           trackCopy->SetP(v);//setting momentum
 568           trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
 569
 570           const AliFmThreeVectorD kP(pxyz[0],pxyz[1],pxyz[2]);
 571           const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
 572           //setting helix I do not if it is ok
 573           AliFmPhysicalHelixD helix(kP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge()));
 574           trackCopy->SetHelix(helix);
 575
 576           //some stuff which could be useful
 577           trackCopy->SetImpactD(impact[0]);
 578           trackCopy->SetImpactZ(impact[1]);
 579           trackCopy->SetCdd(covimpact[0]);
 580           trackCopy->SetCdz(covimpact[1]);
 581           trackCopy->SetCzz(covimpact[2]);
 582           trackCopy->SetSigmaToVertex(GetSigmaToVertex(impact, covimpact));
 583
 584           delete param;
 585         }
 586         else {
 587           if (fReadInner == true) {
 588
 589             if (esdtrack->GetTPCInnerParam()) {
 590               AliExternalTrackParam *param = new AliExternalTrackParam(*esdtrack->GetInnerParam());
 591               //trackCopy->SetInnerMomentum(param->P());
 592               trackCopy->SetInnerMomentum(esdtrack->GetTPCmomentum());
 593               param->GetXYZ(rxyz);
 594               //            param->PropagateToDCA(fEvent->GetPrimaryVertex(), (fEvent->GetMagneticField()), 10000);
 595               param->GetPxPyPz(pxyz);//reading noconstarined momentum
 596               delete param;
 597
 598               AliFemtoModelHiddenInfo *tInfo = new AliFemtoModelHiddenInfo();
 599               tInfo->SetPDGPid(211);
 600               tInfo->SetTrueMomentum(pxyz[0], pxyz[1], pxyz[2]);
 601               tInfo->SetMass(0.13957);
 602               //            tInfo->SetEmissionPoint(rxyz[0], rxyz[1], rxyz[2], 0.0);
 603               //tInfo->SetEmissionPoint(fV1[0], fV1[1], fV1[2], 0.0);
 604               tInfo->SetEmissionPoint(rxyz[0]-fV1[0], rxyz[1]-fV1[1], rxyz[2]-fV1[2], 0.0);
 605               trackCopy->SetHiddenInfo(tInfo);
 606             }
 607           }
 608
 609           if (fTrackType == kGlobal) {
 610             if (fConstrained==true)
 611               tGoodMomentum=esdtrack->GetConstrainedPxPyPz(pxyz); //reading constrained momentum
 612             else
 613               tGoodMomentum=esdtrack->GetPxPyPz(pxyz);//reading noconstarined momentum
 614           }
 615           else if (fTrackType == kTPCOnly) {
 616             if (esdtrack->GetTPCInnerParam())
 617               esdtrack->GetTPCInnerParam()->GetPxPyPz(pxyz);
 618             else {
 619               delete trackCopy;
 620               continue;
 621             }
 622           }
 623           else if (fTrackType == kITSOnly) {
 624             if (fConstrained==true)
 625               tGoodMomentum=esdtrack->GetConstrainedPxPyPz(pxyz); //reading constrained momentum
 626             else
 627               tGoodMomentum=esdtrack->GetPxPyPz(pxyz);//reading noconstarined momentum
 628           }
 629
 630
 631           AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
 632           if (v.Mag() < 0.0001) {
 633             //  cout << "Found 0 momentum ???? " <<endl;
 634             delete trackCopy;
 635             continue;
 636           }
 637           trackCopy->SetP(v);//setting momentum
 638           trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
 639           const AliFmThreeVectorD kP(pxyz[0],pxyz[1],pxyz[2]);
 640           const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
 641           //setting helix I do not if it is ok
 642           AliFmPhysicalHelixD helix(kP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge()));
 643           trackCopy->SetHelix(helix);
 644
 645           //some stuff which could be useful
 646           float imp[2];
 647           float cim[3];
 648           // if (fTrackType == kTPCOnly) {
 649           //   esdtrack->GetTPCInnerParam()->GetImpactParameters(imp,cim);
 650           // }
 651           // else {
 652           esdtrack->GetImpactParameters(imp,cim);
 653           // }
 654
 655           impact[0] = imp[0];
 656           impact[1] = imp[1];
 657           covimpact[0] = cim[0];
 658           covimpact[1] = cim[1];
 659           covimpact[2] = cim[2];
 660
 661           trackCopy->SetImpactD(impact[0]);
 662           trackCopy->SetImpactZ(impact[1]);
 663           trackCopy->SetCdd(covimpact[0]);
 664           trackCopy->SetCdz(covimpact[1]);
 665           trackCopy->SetCzz(covimpact[2]);
 666           trackCopy->SetSigmaToVertex(GetSigmaToVertex(impact,covimpact));
 667         }
 668
 669         trackCopy->SetTrackId(esdtrack->GetID());
 670         trackCopy->SetFlags(esdtrack->GetStatus());
 671         trackCopy->SetLabel(esdtrack->GetLabel());
 672
 673         trackCopy->SetITSchi2(esdtrack->GetITSchi2());
 674         if (esdtrack->GetITSFakeFlag())
 675           trackCopy->SetITSncls(-esdtrack->GetNcls(0));
 676         else
 677           trackCopy->SetITSncls(esdtrack->GetNcls(0));
 678         trackCopy->SetTPCchi2(esdtrack->GetTPCchi2());
 679         trackCopy->SetTPCncls(esdtrack->GetTPCNcls());
 680         trackCopy->SetTPCnclsF(esdtrack->GetTPCNclsF());
 681         trackCopy->SetTPCsignal(esdtrack->GetTPCsignal());
 682         trackCopy->SetTPCsignalN((short)esdtrack->GetTPCsignalN()); //due to bug in aliesdtrack class
 683         trackCopy->SetTPCsignalS(esdtrack->GetTPCsignalSigma());
 684
 685         trackCopy->SetTPCClusterMap(esdtrack->GetTPCClusterMap());
 686         trackCopy->SetTPCSharedMap(esdtrack->GetTPCSharedMap());
 687
 688         double xtpc[3];
 689         esdtrack->GetInnerXYZ(xtpc);
 690         xtpc[2] -= fV1[2];
 691         trackCopy->SetNominalTPCEntrancePoint(xtpc);
 692
 693         esdtrack->GetOuterXYZ(xtpc);
 694         xtpc[2] -= fV1[2];
 695         trackCopy->SetNominalTPCExitPoint(xtpc);
 696
 697         int indexes[3];
 698         for (int ik=0; ik<3; ik++) {
 699           indexes[ik] = esdtrack->GetKinkIndex(ik);
 700         }
 701         trackCopy->SetKinkIndexes(indexes);
 702
 703         for (int ii=0; ii<6; ii++){
 704           trackCopy->SetITSHitOnLayer(ii,esdtrack->HasPointOnITSLayer(ii));
 705         }
 706
 707         //decision if we want this track
 708         //if we using diffrent labels we want that this label was use for first time
 709         //if we use hidden info we want to have match between sim data and ESD
 710         if (tGoodMomentum==true)
 711           {
 712             hbtEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
 713             realnofTracks++;//real number of tracks
 714             //    delete trackCopy;
 715           }
 716         else
 717           {
 718             delete  trackCopy;
 719           }
 720
 721     }
 722
 723   hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event
 724
 725   AliCentrality *cent = fEvent->GetCentrality();
 726   if (cent) {
 727     hbtEvent->SetCentralityV0(cent->GetCentralityPercentile("V0M"));
 728     hbtEvent->SetCentralityV0A(cent->GetCentralityPercentile("V0A"));
 729     hbtEvent->SetCentralityV0C(cent->GetCentralityPercentile("V0C"));
 730     hbtEvent->SetCentralityZNA(cent->GetCentralityPercentile("ZNA"));
 731     hbtEvent->SetCentralityZNC(cent->GetCentralityPercentile("ZNC"));
 732     hbtEvent->SetCentralityCL1(cent->GetCentralityPercentile("CL1"));
 733     hbtEvent->SetCentralityCL0(cent->GetCentralityPercentile("CL0"));
 734     hbtEvent->SetCentralityTKL(cent->GetCentralityPercentile("TKL"));
 735     hbtEvent->SetCentralityFMD(cent->GetCentralityPercentile("FMD"));
 736     hbtEvent->SetCentralityFMD(cent->GetCentralityPercentile("NPA"));
 737     //    hbtEvent->SetCentralitySPD1(cent->GetCentralityPercentile("CL1"));
 738     hbtEvent->SetCentralityTrk(cent->GetCentralityPercentile("TRK"));
 739     hbtEvent->SetCentralityCND(cent->GetCentralityPercentile("CND"));
 740
 741     if (Debug()>1) printf("  FemtoReader Got Event with %f %f %f %f\n", cent->GetCentralityPercentile("V0M"), 0.0, cent->GetCentralityPercentile("CL1"), 0.0);
 742   }
 743
 744   if (fEstEventMult == kGlobalCount)
 745     hbtEvent->SetNormalizedMult(tNormMult);
 746   else if (fEstEventMult == kReferenceITSTPC)
 747     hbtEvent->SetNormalizedMult(AliESDtrackCuts::GetReferenceMultiplicity(fEvent,AliESDtrackCuts::kTrackletsITSTPC,1.0));
 748   else if(fEstEventMult == kReferenceITSSA)
 749     hbtEvent->SetNormalizedMult(AliESDtrackCuts::GetReferenceMultiplicity(fEvent,AliESDtrackCuts::kTrackletsITSSA,1.0));
 750   else if(fEstEventMult == kReferenceTracklets)
 751     hbtEvent->SetNormalizedMult(AliESDtrackCuts::GetReferenceMultiplicity(fEvent,AliESDtrackCuts::kTracklets,1.0));
 752   else if (fEstEventMult == kSPDLayer1)
 753     hbtEvent->SetNormalizedMult(fEvent->GetMultiplicity()->GetNumberOfITSClusters(1));
 754   else if (fEstEventMult == kVZERO)
 755     {
 756       Float_t multV0 = 0;
 757       for (Int_t i=0; i<64; i++)
 758         multV0 += fEvent->GetVZEROData()->GetMultiplicity(i);
 759       hbtEvent->SetNormalizedMult(multV0);
 760     }
 761   else if (fEstEventMult == kCentrality) {
 762     // centrality between 0 (central) and 1 (very peripheral)
 763
 764     if (cent) {
 765       if (cent->GetCentralityPercentile("V0M") < 0.00001)
 766         hbtEvent->SetNormalizedMult(-1);
 767       else
 768         hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("V0M")));
 769       if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(),
 770                              10.0*cent->GetCentralityPercentile("V0M"), lrint(10.0*cent->GetCentralityPercentile("V0M")));
 771     }
 772   }
 773   else if (fEstEventMult == kCentralityV0A) {
 774     // centrality between 0 (central) and 1 (very peripheral)
 775
 776     if (cent) {
 777       if (cent->GetCentralityPercentile("V0A") < 0.00001)
 778         hbtEvent->SetNormalizedMult(-1);
 779       else
 780         hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("V0A")));
 781       if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(),
 782                              10.0*cent->GetCentralityPercentile("V0A"), lrint(10.0*cent->GetCentralityPercentile("V0A")));
 783     }
 784   }
 785   else if (fEstEventMult == kCentralityV0C) {
 786     // centrality between 0 (central) and 1 (very peripheral)
 787
 788     if (cent) {
 789       if (cent->GetCentralityPercentile("V0C") < 0.00001)
 790         hbtEvent->SetNormalizedMult(-1);
 791       else
 792         hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("V0C")));
 793       if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(),
 794                              10.0*cent->GetCentralityPercentile("V0C"), lrint(10.0*cent->GetCentralityPercentile("V0C")));
 795     }
 796   }
 797   else if (fEstEventMult == kCentralityZNA) {
 798     // centrality between 0 (central) and 1 (very peripheral)
 799
 800     if (cent) {
 801       if (cent->GetCentralityPercentile("ZNA") < 0.00001)
 802         hbtEvent->SetNormalizedMult(-1);
 803       else
 804         hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("ZNA")));
 805       if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(),
 806                              10.0*cent->GetCentralityPercentile("ZNA"), lrint(10.0*cent->GetCentralityPercentile("ZNA")));
 807     }
 808   }
 809   else if (fEstEventMult == kCentralityZNC) {
 810     // centrality between 0 (central) and 1 (very peripheral)
 811
 812     if (cent) {
 813       if (cent->GetCentralityPercentile("ZNC") < 0.00001)
 814         hbtEvent->SetNormalizedMult(-1);
 815       else
 816         hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("ZNC")));
 817       if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(),
 818                              10.0*cent->GetCentralityPercentile("ZNC"), lrint(10.0*cent->GetCentralityPercentile("ZNC")));
 819     }
 820   }
 821   else if (fEstEventMult == kCentralityCL1) {
 822     // centrality between 0 (central) and 1 (very peripheral)
 823
 824     if (cent) {
 825       if (cent->GetCentralityPercentile("CL1") < 0.00001)
 826         hbtEvent->SetNormalizedMult(-1);
 827       else
 828         hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("CL1")));
 829       if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(),
 830                              10.0*cent->GetCentralityPercentile("CL1"), lrint(10.0*cent->GetCentralityPercentile("CL1")));
 831     }
 832   }
 833   else if (fEstEventMult == kCentralityCL0) {
 834     // centrality between 0 (central) and 1 (very peripheral)
 835
 836     if (cent) {
 837       if (cent->GetCentralityPercentile("CL0") < 0.00001)
 838         hbtEvent->SetNormalizedMult(-1);
 839       else
 840         hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("CL0")));
 841       if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(),
 842                              10.0*cent->GetCentralityPercentile("CL0"), lrint(10.0*cent->GetCentralityPercentile("CL0")));
 843     }
 844   }
 845   else if (fEstEventMult == kCentralityTRK) {
 846     // centrality between 0 (central) and 1 (very peripheral)
 847
 848     if (cent) {
 849       if (cent->GetCentralityPercentile("TRK") < 0.00001)
 850         hbtEvent->SetNormalizedMult(-1);
 851       else
 852         hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("TRK")));
 853       if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(),
 854                              10.0*cent->GetCentralityPercentile("TRK"), lrint(10.0*cent->GetCentralityPercentile("TRK")));
 855     }
 856   }
 857   else if (fEstEventMult == kCentralityTKL) {
 858     // centrality between 0 (central) and 1 (very peripheral)
 859
 860     if (cent) {
 861       if (cent->GetCentralityPercentile("TKL") < 0.00001)
 862         hbtEvent->SetNormalizedMult(-1);
 863       else
 864         hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("TKL")));
 865       if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(),
 866                              10.0*cent->GetCentralityPercentile("TKL"), lrint(10.0*cent->GetCentralityPercentile("TKL")));
 867     }
 868   }
 869   else if (fEstEventMult == kCentralityNPA) {
 870     // centrality between 0 (central) and 1 (very peripheral)
 871
 872     if (cent) {
 873       if (cent->GetCentralityPercentile("NPA") < 0.00001)
 874         hbtEvent->SetNormalizedMult(-1);
 875       else
 876         hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("NPA")));
 877       if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(),
 878                              10.0*cent->GetCentralityPercentile("NPA"), lrint(10.0*cent->GetCentralityPercentile("NPA")));
 879     }
 880   }
 881   else if (fEstEventMult == kCentralityCND) {
 882     // centrality between 0 (central) and 1 (very peripheral)
 883
 884     if (cent) {
 885       if (cent->GetCentralityPercentile("CND") < 0.00001)
 886         hbtEvent->SetNormalizedMult(-1);
 887       else
 888         hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("CND")));
 889       if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(),
 890                              10.0*cent->GetCentralityPercentile("CND"), lrint(10.0*cent->GetCentralityPercentile("CND")));
 891     }
 892   }
 893   else if (fEstEventMult == kCentralityFMD) {
 894     // centrality between 0 (central) and 1 (very peripheral)
 895
 896     if (cent) {
 897       if (cent->GetCentralityPercentile("FMD") < 0.00001)
 898         hbtEvent->SetNormalizedMult(-1);
 899       else
 900         hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("FMD")));
 901       if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(),
 902                              10.0*cent->GetCentralityPercentile("FMD"), lrint(10.0*cent->GetCentralityPercentile("FMD")));
 903     }
 904   }
 905
 906
 907
 908   if (tNormMultPos > tNormMultNeg)
 909     hbtEvent->SetZDCParticipants(tNormMultPos);
 910   else
 911     hbtEvent->SetZDCParticipants(tNormMultNeg);
 912
 913   AliEventplane* ep = fEvent->GetEventplane();
 914   if (ep) {
 915     hbtEvent->SetEP(ep);
 916     hbtEvent->SetReactionPlaneAngle(ep->GetEventplane("Q"));
 917   }
 918
 919   //V0
 920   if(fReadV0)
 921     {
 922       for (Int_t i = 0; i < fEvent->GetNumberOfV0s(); i++) {
 923         AliESDv0* esdv0 = fEvent->GetV0(i);
 924         if (!esdv0) continue;
 925         //if(esdv0->GetNDaughters()>2) continue;
 926         //if(esdv0->GetNProngs()>2) continue;
 927         if(esdv0->Charge()!=0) continue;
 928         AliESDtrack *trackPos = fEvent->GetTrack(esdv0->GetPindex());
 929         if(!trackPos) continue;
 930         AliESDtrack *trackNeg = fEvent->GetTrack(esdv0->GetNindex());
 931         if(!trackNeg) continue;
 932         if(trackPos->Charge()==trackNeg->Charge()) continue;
 933         AliFemtoV0* trackCopyV0 = new AliFemtoV0();
 934         CopyESDtoFemtoV0(esdv0, trackCopyV0, fEvent);
 935         hbtEvent->V0Collection()->push_back(trackCopyV0);
 936         //cout<<"Pushback v0 to v0collection"<<endl;
 937       }
 938     }
 939
 940
 941   if (Debug()>1) cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
 942 }
 943
 944 //___________________
 945 void AliFemtoEventReaderESDChain::SetESDSource(AliESDEvent *aESD)
 946 {
 947   // The chain loads the ESD for us
 948   // You must provide the address where it can be found
 949   fEvent = aESD;
 950 }
 951 //___________________
 952 // void AliFemtoEventReaderESDChain::SetESDfriendSource(AliESDfriend *aFriend)
 953 // {
 954 //   // We need the ESD tree to obtain
 955 //   // information about the friend file location
 956 //   fEventFriend = aFriend;
 957 // }
 958
 959 //____________________________________________________________________
 960 Float_t AliFemtoEventReaderESDChain::GetSigmaToVertex(double *impact, double *covar)
 961 {
 962   // Calculates the number of sigma to the vertex.
 963
 964   Float_t b[2];
 965   Float_t bRes[2];
 966   Float_t bCov[3];
 967
 968   b[0] = impact[0];
 969   b[1] = impact[1];
 970   bCov[0] = covar[0];
 971   bCov[1] = covar[1];
 972   bCov[2] = covar[2];
 973
 974   bRes[0] = TMath::Sqrt(bCov[0]);
 975   bRes[1] = TMath::Sqrt(bCov[2]);
 976
 977   // -----------------------------------
 978   // How to get to a n-sigma cut?
 979   //
 980   // The accumulated statistics from 0 to d is
 981   //
 982   // ->  Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
 983   // ->  1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
 984   //
 985   // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-x**2)/2)
 986   // Can this be expressed in a different way?
 987
 988   if (bRes[0] == 0 || bRes[1] ==0)
 989     return -1;
 990
 991   Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
 992
 993   // stupid rounding problem screws up everything:
 994   // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
 995   if (TMath::Exp(-d * d / 2) < 1e-10)
 996     return 1000;
 997
 998   d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
 999   return d;
1000 }
1001
1002 void AliFemtoEventReaderESDChain::CopyESDtoFemtoV0(AliESDv0 *tESDv0, AliFemtoV0 *tFemtoV0, AliESDEvent *tESDevent)
1003 {
1004   double fPrimaryVtxPosition[3];
1005   tESDevent->GetPrimaryVertex()->GetXYZ(fPrimaryVtxPosition);
1006   tFemtoV0->SetdcaV0ToPrimVertex(tESDv0->GetD(fPrimaryVtxPosition[0],fPrimaryVtxPosition[1],fPrimaryVtxPosition[2]));
1007
1008   //tFemtoV0->SetdecayLengthV0(tESDv0->DecayLengthV0()); //wrocic do tego
1009   //tFemtoV0->SetdecayVertexV0X(tESDv0->DecayVertexV0X());
1010   //tFemtoV0->SetdecayVertexV0Y(tESDv0->DecayVertexV0Y());
1011   //tFemtoV0->SetdecayVertexV0Z(tESDv0->DecayVertexV0Z()); //nie ma w AliESDv0
1012   //AliFemtoThreeVector decayvertex(tESDv0->DecayVertexV0X(),tESDv0->DecayVertexV0Y(),tESDv0->DecayVertexV0Z());
1013   //tFemtoV0->SetdecayVertexV0(decayvertex);
1014   tFemtoV0->SetdcaV0Daughters(tESDv0->GetDcaV0Daughters());
1015   tFemtoV0->SetmomV0X(tESDv0->Px());
1016   tFemtoV0->SetmomV0Y(tESDv0->Py());
1017   tFemtoV0->SetmomV0Z(tESDv0->Pz());
1018   AliFemtoThreeVector momv0(tESDv0->Px(),tESDv0->Py(),tESDv0->Pz());
1019   tFemtoV0->SetmomV0(momv0);
1020   tFemtoV0->SetalphaV0(tESDv0->AlphaV0());
1021   tFemtoV0->SetptArmV0(tESDv0->PtArmV0());
1022   //tFemtoV0->SeteLambda(tESDv0->ELambda());
1023   //tFemtoV0->SeteK0Short(tESDv0->EK0Short());
1024   //tFemtoV0->SetePosProton(tESDv0->EPosProton());
1025   //tFemtoV0->SeteNegProton(tESDv0->ENegProton());
1026   tFemtoV0->SetmassLambda(tESDv0->GetEffMass(4,2));
1027   tFemtoV0->SetmassAntiLambda(tESDv0->GetEffMass(2,4));
1028   tFemtoV0->SetmassK0Short(tESDv0->GetEffMass(2,2));
1029   //tFemtoV0->SetrapLambda(tESDv0->RapLambda());
1030   //tFemtoV0->SetrapK0Short(tESDv0->RapK0Short());
1031   tFemtoV0->SetptV0(tESDv0->Pt());
1032   tFemtoV0->SetptotV0(tESDv0->P());
1033   tFemtoV0->SetEtaV0(tESDv0->Eta());
1034   tFemtoV0->SetPhiV0(tESDv0->Phi());
1035   tFemtoV0->SetCosPointingAngle(tESDv0->GetV0CosineOfPointingAngle(fPrimaryVtxPosition[0],fPrimaryVtxPosition[1], fPrimaryVtxPosition[2]));
1036   tFemtoV0->SetYV0(tESDv0->Y());
1037
1038   AliESDtrack *trackpos = tESDevent->GetTrack(tESDv0->GetPindex()); //AliAODTrack *trackpos = (AliAODTrack*)tESDv0->GetDaughter(0);
1039   AliESDtrack *trackneg = tESDevent->GetTrack(tESDv0->GetNindex()); //AliAODTrack *trackneg = (AliAODTrack*)tESDv0->GetDaughter(1);
1040
1041   if(trackpos && trackneg)
1042     {
1043       tFemtoV0->SetdcaPosToPrimVertex(TMath::Abs(trackpos->GetD(fPrimaryVtxPosition[0],fPrimaryVtxPosition[1],tESDevent->GetMagneticField())));
1044       tFemtoV0->SetdcaNegToPrimVertex(TMath::Abs(trackneg->GetD(fPrimaryVtxPosition[0],fPrimaryVtxPosition[1],tESDevent->GetMagneticField())));
1045       double MomPos[3];
1046       trackpos->PxPyPz(MomPos);
1047       tFemtoV0->SetmomPosX(MomPos[0]);
1048       tFemtoV0->SetmomPosY(MomPos[1]);
1049       tFemtoV0->SetmomPosZ(MomPos[2]);
1050       AliFemtoThreeVector mompos(MomPos[0],MomPos[1],MomPos[2]);
1051       tFemtoV0->SetmomPos(mompos);
1052
1053       double MomNeg[3];
1054       trackneg->PxPyPz(MomNeg);
1055       tFemtoV0->SetmomNegX(MomNeg[0]);
1056       tFemtoV0->SetmomNegY(MomNeg[1]);
1057       tFemtoV0->SetmomNegZ(MomNeg[2]);
1058       AliFemtoThreeVector momneg(MomNeg[0],MomNeg[1],MomNeg[2]);
1059       tFemtoV0->SetmomNeg(momneg);
1060
1061       tFemtoV0->SetptPos(trackpos->Pt());
1062       tFemtoV0->SetptotPos(trackpos->P());
1063       tFemtoV0->SetptNeg(trackneg->Pt());
1064       tFemtoV0->SetptotNeg(trackneg->P());
1065
1066       tFemtoV0->SetidNeg(trackneg->GetID());
1067       //cout<<"tESDv0->GetNegID(): "<<tESDv0->GetNegID()<<endl;
1068       //cout<<"tFemtoV0->IdNeg(): "<<tFemtoV0->IdNeg()<<endl;
1069       tFemtoV0->SetidPos(trackpos->GetID());
1070
1071       tFemtoV0->SetEtaPos(trackpos->Eta());
1072       tFemtoV0->SetEtaNeg(trackneg->Eta());
1073
1074       tFemtoV0->SetEtaPos(trackpos->Eta()); //tESDv0->PseudoRapPos()
1075       tFemtoV0->SetEtaNeg(trackneg->Eta()); //tESDv0->PseudoRapNeg()
1076       tFemtoV0->SetTPCNclsPos(trackpos->GetTPCNcls());
1077       tFemtoV0->SetTPCNclsNeg(trackneg->GetTPCNcls());
1078       tFemtoV0->SetTPCclustersPos(trackpos->GetTPCClusterMap());
1079       tFemtoV0->SetTPCclustersNeg(trackneg->GetTPCClusterMap());
1080       tFemtoV0->SetTPCsharingPos(trackpos->GetTPCSharedMap());
1081       tFemtoV0->SetTPCsharingNeg(trackneg->GetTPCSharedMap());
1082       tFemtoV0->SetNdofPos(trackpos->GetTPCchi2()/trackpos->GetTPCNcls());
1083       tFemtoV0->SetNdofNeg(trackneg->GetTPCchi2()/trackneg->GetTPCNcls());
1084       tFemtoV0->SetStatusPos(trackpos->GetStatus());
1085       tFemtoV0->SetStatusNeg(trackneg->GetStatus());
1086
1087       float bfield = 5*fMagFieldSign;
1088       float globalPositionsAtRadiiPos[9][3];
1089       GetGlobalPositionAtGlobalRadiiThroughTPC(trackpos,bfield,globalPositionsAtRadiiPos);
1090       double tpcEntrancePos[3]={globalPositionsAtRadiiPos[0][0],globalPositionsAtRadiiPos[0][1],globalPositionsAtRadiiPos[0][2]};
1091       double tpcExitPos[3]={globalPositionsAtRadiiPos[7][0],globalPositionsAtRadiiPos[7][1],globalPositionsAtRadiiPos[7][2]};
1092
1093       float globalPositionsAtRadiiNeg[9][3];
1094       GetGlobalPositionAtGlobalRadiiThroughTPC(trackneg,bfield,globalPositionsAtRadiiNeg);
1095       double tpcEntranceNeg[3]={globalPositionsAtRadiiNeg[0][0],globalPositionsAtRadiiNeg[0][1],globalPositionsAtRadiiNeg[0][2]};
1096       double tpcExitNeg[3]={globalPositionsAtRadiiNeg[7][0],globalPositionsAtRadiiNeg[7][1],globalPositionsAtRadiiNeg[7][2]};
1097
1098       AliFemtoThreeVector tmpVec;
1099       tmpVec.SetX(tpcEntrancePos[0]); tmpVec.SetX(tpcEntrancePos[1]); tmpVec.SetX(tpcEntrancePos[2]);
1100       tFemtoV0->SetNominalTpcEntrancePointPos(tmpVec);
1101
1102       tmpVec.SetX(tpcExitPos[0]); tmpVec.SetX(tpcExitPos[1]); tmpVec.SetX(tpcExitPos[2]);
1103       tFemtoV0->SetNominalTpcExitPointPos(tmpVec);
1104
1105       tmpVec.SetX(tpcEntranceNeg[0]); tmpVec.SetX(tpcEntranceNeg[1]); tmpVec.SetX(tpcEntranceNeg[2]);
1106       tFemtoV0->SetNominalTpcEntrancePointNeg(tmpVec);
1107
1108       tmpVec.SetX(tpcExitNeg[0]); tmpVec.SetX(tpcExitNeg[1]); tmpVec.SetX(tpcExitNeg[2]);
1109       tFemtoV0->SetNominalTpcExitPointNeg(tmpVec);
1110
1111       AliFemtoThreeVector vecTpcPos[9];
1112       AliFemtoThreeVector vecTpcNeg[9];
1113       for(int i=0;i<9;i++)
1114         {
1115           vecTpcPos[i].SetX(globalPositionsAtRadiiPos[i][0]); vecTpcPos[i].SetY(globalPositionsAtRadiiPos[i][1]); vecTpcPos[i].SetZ(globalPositionsAtRadiiPos[i][2]);
1116           vecTpcNeg[i].SetX(globalPositionsAtRadiiNeg[i][0]); vecTpcNeg[i].SetY(globalPositionsAtRadiiNeg[i][1]); vecTpcNeg[i].SetZ(globalPositionsAtRadiiNeg[i][2]);
1117         }
1118       tFemtoV0->SetNominalTpcPointPos(vecTpcPos);
1119       tFemtoV0->SetNominalTpcPointNeg(vecTpcNeg);
1120
1121       tFemtoV0->SetTPCMomentumPos(trackpos->GetTPCInnerParam()->P()); //trackpos->GetTPCmomentum();
1122       tFemtoV0->SetTPCMomentumNeg(trackneg->GetTPCInnerParam()->P()); //trackneg->GetTPCmomentum();
1123
1124       tFemtoV0->SetdedxPos(trackpos->GetTPCsignal());
1125       tFemtoV0->SetdedxNeg(trackneg->GetTPCsignal());
1126
1127
1128       if (fESDpid) {
1129         tFemtoV0->SetPosNSigmaTPCK(fESDpid->NumberOfSigmasTPC(trackpos,AliPID::kKaon));
1130         tFemtoV0->SetNegNSigmaTPCK(fESDpid->NumberOfSigmasTPC(trackneg,AliPID::kKaon));
1131         tFemtoV0->SetPosNSigmaTPCP(fESDpid->NumberOfSigmasTPC(trackpos,AliPID::kProton));
1132         tFemtoV0->SetNegNSigmaTPCP(fESDpid->NumberOfSigmasTPC(trackneg,AliPID::kProton));
1133         tFemtoV0->SetPosNSigmaTPCPi(fESDpid->NumberOfSigmasTPC(trackpos,AliPID::kPion));
1134         tFemtoV0->SetNegNSigmaTPCPi(fESDpid->NumberOfSigmasTPC(trackneg,AliPID::kPion));
1135       }
1136       else {
1137         tFemtoV0->SetPosNSigmaTPCK(-1000);
1138         tFemtoV0->SetNegNSigmaTPCK(-1000);
1139         tFemtoV0->SetPosNSigmaTPCP(-1000);
1140         tFemtoV0->SetNegNSigmaTPCP(-1000);
1141         tFemtoV0->SetPosNSigmaTPCPi(-1000);
1142         tFemtoV0->SetNegNSigmaTPCPi(-1000);
1143       }
1144
1145       if((tFemtoV0->StatusPos()&AliESDtrack::kTOFpid)==0 || (tFemtoV0->StatusPos()&AliESDtrack::kTIME)==0 || (tFemtoV0->StatusPos()&AliESDtrack::kTOFout)==0)
1146         {
1147           if((tFemtoV0->StatusNeg()&AliESDtrack::kTOFpid)==0 || (tFemtoV0->StatusNeg()&AliESDtrack::kTIME)==0 || (tFemtoV0->StatusNeg()&AliESDtrack::kTOFout)==0)
1148             {
1149               tFemtoV0->SetPosNSigmaTOFK(-1000);
1150               tFemtoV0->SetNegNSigmaTOFK(-1000);
1151               tFemtoV0->SetPosNSigmaTOFP(-1000);
1152               tFemtoV0->SetNegNSigmaTOFP(-1000);
1153               tFemtoV0->SetPosNSigmaTOFPi(-1000);
1154               tFemtoV0->SetNegNSigmaTOFPi(-1000);
1155             }
1156         }
1157       else
1158         {
1159           if (fESDpid) {
1160             tFemtoV0->SetPosNSigmaTOFK(fESDpid->NumberOfSigmasTOF(trackpos,AliPID::kKaon));
1161             tFemtoV0->SetNegNSigmaTOFK(fESDpid->NumberOfSigmasTOF(trackneg,AliPID::kKaon));
1162             tFemtoV0->SetPosNSigmaTOFP(fESDpid->NumberOfSigmasTOF(trackpos,AliPID::kProton));
1163             tFemtoV0->SetNegNSigmaTOFP(fESDpid->NumberOfSigmasTOF(trackneg,AliPID::kProton));
1164             tFemtoV0->SetPosNSigmaTOFPi(fESDpid->NumberOfSigmasTOF(trackpos,AliPID::kPion));
1165             tFemtoV0->SetNegNSigmaTOFPi(fESDpid->NumberOfSigmasTOF(trackneg,AliPID::kPion));
1166           }
1167           else {
1168             tFemtoV0->SetPosNSigmaTOFK(-1000);
1169             tFemtoV0->SetNegNSigmaTOFK(-1000);
1170             tFemtoV0->SetPosNSigmaTOFP(-1000);
1171             tFemtoV0->SetNegNSigmaTOFP(-1000);
1172             tFemtoV0->SetPosNSigmaTOFPi(-1000);
1173             tFemtoV0->SetNegNSigmaTOFPi(-1000);
1174           }
1175         }
1176     }
1177   else
1178     {
1179       tFemtoV0->SetStatusPos(999);
1180       tFemtoV0->SetStatusNeg(999);
1181     }
1182   tFemtoV0->SetOnFlyStatusV0(tESDv0->GetOnFlyStatus());
1183 }
1184
1185 void AliFemtoEventReaderESDChain::SetReadTrackType(ReadTrackType aType)
1186 {
1187   fTrackType = aType;
1188 }
1189
1190 //trigger
1191 void AliFemtoEventReaderESDChain::SetEventTrigger(UInt_t eventtrig)
1192 {
1193   fEventTrig = eventtrig;
1194 }
1195
1196 //V0 reading
1197 void AliFemtoEventReaderESDChain::SetReadV0(bool a)
1198 {
1199   fReadV0 = a;
1200 }
1201
1202 void AliFemtoEventReaderESDChain::SetMagneticFieldSign(int s)
1203 {
1204   if(s>0)
1205     fMagFieldSign = 1;
1206   else if(s<0)
1207     fMagFieldSign = -1;
1208   else
1209     fMagFieldSign = 0;
1210 }
1211
1212 void AliFemtoEventReaderESDChain::GetGlobalPositionAtGlobalRadiiThroughTPC(AliESDtrack *track, Float_t bfield, Float_t globalPositionsAtRadii[9][3])
1213 {
1214   // Gets the global position of the track at nine different radii in the TPC
1215   // track is the track you want to propagate
1216   // bfield is the magnetic field of your event
1217   // globalPositionsAtRadii is the array of global positions in the radii and xyz
1218
1219   // Initialize the array to something indicating there was no propagation
1220   for(Int_t i=0;i<9;i++){
1221     for(Int_t j=0;j<3;j++){
1222       globalPositionsAtRadii[i][j]=-9999.;
1223     }
1224   }
1225
1226   // Make a copy of the track to not change parameters of the track
1227   AliExternalTrackParam etp; etp.CopyFromVTrack(track);
1228   //printf("\nAfter CopyFromVTrack\n");
1229   //etp.Print();
1230
1231   // The global position of the the track
1232   Double_t xyz[3]={-9999.,-9999.,-9999.};
1233
1234   // Counter for which radius we want
1235   Int_t iR=0;
1236   // The radii at which we get the global positions
1237   // IROC (OROC) from 84.1 cm to 132.1 cm (134.6 cm to 246.6 cm)
1238   Float_t Rwanted[9]={85.,105.,125.,145.,165.,185.,205.,225.,245.};
1239   // The global radius we are at
1240   Float_t globalRadius=0;
1241
1242   // Propagation is done in local x of the track
1243   for (Float_t x = etp.GetX();x<247.;x+=1.){ // GetX returns local coordinates
1244     // Starts at the tracks fX and goes outwards. x = 245 is the outer radial limit
1245     // of the TPC when the track is straight, i.e. has inifinite pt and doesn't get bent.
1246     // If the track's momentum is smaller than infinite, it will develop a y-component, which
1247     // adds to the global radius
1248
1249     // Stop if the propagation was not succesful. This can happen for low pt tracks
1250     // that don't reach outer radii
1251     if(!etp.PropagateTo(x,bfield))break;
1252     etp.GetXYZ(xyz); // GetXYZ returns global coordinates
1253     globalRadius = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]); //Idea to speed up: compare squared radii
1254
1255     // Roughly reached the radius we want
1256     if(globalRadius > Rwanted[iR]){
1257
1258       // Bigger loop has bad precision, we're nearly one centimeter too far, go back in small steps.
1259       while (globalRadius>Rwanted[iR]){
1260         x-=.1;
1261         //      printf("propagating to x %5.2f\n",x);
1262         if(!etp.PropagateTo(x,bfield))break;
1263         etp.GetXYZ(xyz); // GetXYZ returns global coordinates
1264         globalRadius = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]); //Idea to speed up: compare squared radii
1265       }
1266       //printf("At Radius:%05.2f (local x %5.2f). Setting position to x %4.1f y %4.1f z %4.1f\n",globalRadius,x,xyz[0],xyz[1],xyz[2]);
1267       globalPositionsAtRadii[iR][0]=xyz[0];
1268       globalPositionsAtRadii[iR][1]=xyz[1];
1269       globalPositionsAtRadii[iR][2]=xyz[2];
1270       // Indicate we want the next radius
1271       iR+=1;
1272     }
1273     if(iR>=8){
1274       // TPC edge reached
1275       return;
1276     }
1277   }
1278 }
1279
1280
1281 void AliFemtoEventReaderESDChain::SetpA2013(Bool_t pA2013)
1282 {
1283   fpA2013 = pA2013;
1284 }