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