]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGCF/FEMTOSCOPY/AliFemto/AliFemtoEventReaderAOD.cxx
Merge branch 'master_patch'
[u/mrichter/AliRoot.git] / PWGCF / FEMTOSCOPY / AliFemto / AliFemtoEventReaderAOD.cxx
index 6d79692a11197c9a4f9190cc46d6b90e6d3e2e25..6a86baf443918e933b9d75308a112fb97da20e98 100644 (file)
@@ -29,6 +29,9 @@
 #include "AliPID.h"
 
 #include "AliAODpidUtil.h"
+#include "AliAnalysisUtils.h"
+#include "assert.h"
+#include "AliGenHijingEventHeader.h"
 
 ClassImp(AliFemtoEventReaderAOD)
 
@@ -63,6 +66,11 @@ AliFemtoEventReaderAOD::AliFemtoEventReaderAOD():
   fMagFieldSign(1),
   fisEPVZ(kTRUE),
   fpA2013(kFALSE),
+  fisPileUp(kFALSE),
+  fMVPlp(kFALSE),
+  fMinVtxContr(0),
+  fMinPlpContribMV(0),
+  fMinPlpContribSPD(0),
   fDCAglobalTrack(kFALSE),
   fFlatCent(kFALSE)
 {
@@ -95,6 +103,11 @@ AliFemtoEventReaderAOD::AliFemtoEventReaderAOD(const AliFemtoEventReaderAOD &aRe
   fMagFieldSign(1),
   fisEPVZ(kTRUE),
   fpA2013(kFALSE),
+  fisPileUp(kFALSE),
+  fMVPlp(kFALSE),
+  fMinVtxContr(0),
+  fMinPlpContribMV(0),
+  fMinPlpContribSPD(0),
   fDCAglobalTrack(kFALSE),
   fFlatCent(kFALSE)
 {
@@ -117,6 +130,11 @@ AliFemtoEventReaderAOD::AliFemtoEventReaderAOD(const AliFemtoEventReaderAOD &aRe
   fEstEventMult = aReader.fEstEventMult;
   fUsePreCent = aReader.fUsePreCent;
   fpA2013 = aReader.fpA2013;
+  fisPileUp = aReader.fisPileUp;
+  fMVPlp = aReader.fMVPlp;
+  fMinVtxContr = aReader.fMinVtxContr;
+  fMinPlpContribMV = aReader.fMinPlpContribMV;
+  fMinPlpContribSPD = aReader.fMinPlpContribSPD;
   fDCAglobalTrack = aReader.fDCAglobalTrack;
 
 }
@@ -161,6 +179,11 @@ AliFemtoEventReaderAOD& AliFemtoEventReaderAOD::operator=(const AliFemtoEventRea
   fUsePreCent = aReader.fUsePreCent;
   fEstEventMult = aReader.fEstEventMult;
   fpA2013 = aReader.fpA2013;
+  fisPileUp = aReader.fisPileUp;
+  fMVPlp = aReader.fMVPlp;
+  fMinVtxContr = aReader.fMinVtxContr;
+  fMinPlpContribMV = aReader.fMinPlpContribMV;
+  fMinPlpContribSPD = aReader.fMinPlpContribSPD;
   fDCAglobalTrack = aReader.fDCAglobalTrack;
   fFlatCent= aReader.fFlatCent;
 
@@ -246,22 +269,24 @@ AliFemtoEvent* AliFemtoEventReaderAOD::ReturnHbtEvent()
   fTree->GetEvent(fCurEvent);//getting next event
   //  cout << "Read event " << fEvent << " from file " << fTree << endl;
 
-  hbtEvent = new AliFemtoEvent;
+  //hbtEvent = new AliFemtoEvent;
 
-  CopyAODtoFemtoEvent(hbtEvent);
+  hbtEvent = CopyAODtoFemtoEvent();
   fCurEvent++;
 
 
   return hbtEvent;
 }
 
-void AliFemtoEventReaderAOD::CopyAODtoFemtoEvent(AliFemtoEvent *tEvent)
+AliFemtoEvent* AliFemtoEventReaderAOD::CopyAODtoFemtoEvent()
 {
 
   // A function that reads in the AOD event
   // and transfers the neccessary information into
   // the internal AliFemtoEvent
 
+  AliFemtoEvent *tEvent = new AliFemtoEvent();
+
   // setting global event characteristics
   tEvent->SetRunNumber(fEvent->GetRunNumber());
   tEvent->SetMagneticField(fEvent->GetMagneticField()*kilogauss);//to check if here is ok
@@ -289,56 +314,59 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoEvent(AliFemtoEvent *tEvent)
     }
   }
 
-  tEvent->SetReactionPlaneAngle(fEvent->GetHeader()->GetQTheta(0)/2.0);
-
-  Int_t *motherids=0;
-  if (mcP) {
-    motherids = new Int_t[((AliAODMCParticle *) mcP->At(mcP->GetEntries()-1))->GetLabel()];
-    for (int ip=0; ip<mcP->GetEntries(); ip++) motherids[ip] = 0;
-
-    // Read in mother ids
-    AliAODMCParticle *motherpart;
-    for (int ip=0; ip<mcP->GetEntries(); ip++) {
-      motherpart = (AliAODMCParticle *) mcP->At(ip);
-      if (motherpart->GetDaughter(0) > 0)
-        motherids[motherpart->GetDaughter(0)] = ip;
-      if (motherpart->GetDaughter(1) > 0)
-        motherids[motherpart->GetDaughter(1)] = ip;
+  AliAODHeader * header = dynamic_cast<AliAODHeader*>(fEvent->GetHeader());
+  assert(header&&"Not a standard AOD");
+
+  tEvent->SetReactionPlaneAngle(header->GetQTheta(0)/2.0);
+  // Int_t *motherids=0;
+  // if (mcP) {
+  //   const int motherTabSize = ((AliAODMCParticle *) mcP->At(mcP->GetEntries()-1))->GetLabel();
+  //   motherids = new int[motherTabSize+1];
+  //   for (int ip=0; ip<motherTabSize+1; ip++) motherids[ip] = 0;
+
+  //   // Read in mother ids
+  //   AliAODMCParticle *motherpart;
+  //   for (int ip=0; ip<mcP->GetEntries(); ip++) {
+  //     motherpart = (AliAODMCParticle *) mcP->At(ip);
+  //     if (motherpart->GetDaughter(0) > 0)
+  //       motherids[motherpart->GetDaughter(0)] = ip;
+  //     if (motherpart->GetDaughter(1) > 0)
+  //       motherids[motherpart->GetDaughter(1)] = ip;
+  //   }
+  // }
+
+  //AliAnalysisUtils
+  if(fisPileUp||fpA2013)
+    {
+      AliAnalysisUtils *anaUtil=new AliAnalysisUtils();
+      if(fMinVtxContr)
+       anaUtil->SetMinVtxContr(fMinVtxContr);
+      if(fpA2013)
+       if(anaUtil->IsVertexSelected2013pA(fEvent)==kFALSE)
+         {
+           delete tEvent;
+           return NULL; //Vertex rejection for pA analysis.
+         }
+      if(fMVPlp) anaUtil->SetUseMVPlpSelection(kTRUE);
+      else anaUtil->SetUseMVPlpSelection(kFALSE);
+      if(fMinPlpContribMV) anaUtil->SetMinPlpContribMV(fMinPlpContribMV);
+      if(fMinPlpContribSPD) anaUtil->SetMinPlpContribSPD(fMinPlpContribSPD);
+      if(fisPileUp)
+       if(anaUtil->IsPileUpEvent(fEvent)) { delete tEvent;return NULL;} //Pile-up rejection.
+      delete anaUtil;
     }
-  }
 
   // Primary Vertex position
-  //  double fV1[3];
-  if(fpA2013)
-  {
-    const AliAODVertex* trkVtx = fEvent->GetPrimaryVertex();
-    if (!trkVtx || trkVtx->GetNContributors()<=0)  return;
-    TString vtxTtl = trkVtx->GetTitle();
-    if (!vtxTtl.Contains("VertexerTracks")) return;
-    Float_t zvtx = trkVtx->GetZ();
-    const AliAODVertex* spdVtx = fEvent->GetPrimaryVertexSPD();
-    if (spdVtx->GetNContributors()<=0)  return;
-    TString vtxTyp = spdVtx->GetTitle();
-    Double_t cov[6]={0};
-    spdVtx->GetCovarianceMatrix(cov);
-    Double_t zRes = TMath::Sqrt(cov[5]);
-    if (vtxTyp.Contains("vertexer:Z") && (zRes>0.25))  return;
-    if (TMath::Abs(spdVtx->GetZ() - trkVtx->GetZ())>0.5)  return;
-
-    if (TMath::Abs(zvtx) > 10)  return;
-  }
-  fEvent->GetPrimaryVertex()->GetPosition(fV1);
+  const AliAODVertex* aodvertex = (AliAODVertex*) fEvent->GetPrimaryVertex();
+  if(!aodvertex || aodvertex->GetNContributors() < 1) { delete tEvent;return NULL;} //Bad vertex, skip event.
 
+  aodvertex->GetPosition(fV1);
   AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]);
   tEvent->SetPrimVertPos(vertex);
 
   //starting to reading tracks
   int nofTracks=0;  //number of reconstructed tracks in event
 
-  // Check to see whether the additional info exists
-  //   if (fPWG2AODTracks)
-  //     nofTracks=fPWG2AODTracks->GetEntries();
-  //   else
   nofTracks=fEvent->GetNumberOfTracks();
 
   AliEventplane *ep = fEvent->GetEventplane();
@@ -357,8 +385,8 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoEvent(AliFemtoEvent *tEvent)
         (cent->GetCentralityPercentile("V0M")*10 > fCentRange[1]))
     {
       // cout << "Centrality " << cent->GetCentralityPercentile("V0M") << " outside of preselection range " << fCentRange[0] << " - " << fCentRange[1] << endl;
-
-      return;
+      delete tEvent;
+      return NULL;
     }
   }
 
@@ -366,7 +394,7 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoEvent(AliFemtoEvent *tEvent)
 //flatten centrality dist.
   if(percent < 9){
     if(fFlatCent){
-      if(RejectEventCentFlat(fEvent->GetMagneticField(),percent)) return;
+      if(RejectEventCentFlat(fEvent->GetMagneticField(),percent)) { delete tEvent; return NULL;}
     }
   }
 
@@ -378,7 +406,8 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoEvent(AliFemtoEvent *tEvent)
 
   // looking for global tracks and saving their numbers to copy from them PID information to TPC-only tracks in the main loop over tracks
   for (int i=0;i<nofTracks;i++) {
-    const AliAODTrack *aodtrack=fEvent->GetTrack(i);
+    const AliAODTrack *aodtrack=dynamic_cast<const AliAODTrack*>(fEvent->GetTrack(i));
+    assert(aodtrack&&"Not a standard AOD");
     if (!aodtrack->TestFilterBit(fFilterBit)) {
       if(aodtrack->GetID() < 0) continue;
       labels[aodtrack->GetID()] = i;
@@ -388,7 +417,7 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoEvent(AliFemtoEvent *tEvent)
   int tNormMult = 0;
   for (int i=0;i<nofTracks;i++)
   {
-    AliFemtoTrack* trackCopy = new AliFemtoTrack();
+    AliFemtoTrack* trackCopy;// = new AliFemtoTrack();
 
 //       if (fPWG2AODTracks) {
 //     // Read tracks from the additional pwg2 specific AOD part
@@ -515,20 +544,21 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoEvent(AliFemtoEvent *tEvent)
     // No additional information exists
     // Read in the normal AliAODTracks
 
-    // const AliAODTrack *aodtrack=fEvent->GetTrack(i); // getting the AODtrack directly
-    AliAODTrack *aodtrack=fEvent->GetTrack(i); // getting the AODtrack directly
+    // const AliAODTrack *aodtrack=dynamic_cast<AliAODTrack*>(fEvent->GetTrack(i));
+    AliAODTrack *aodtrack=dynamic_cast<AliAODTrack*>(fEvent->GetTrack(i));
+    assert(aodtrack&&"Not a standard AOD"); // getting the AODtrack directly
 
 
 
     if (aodtrack->IsPrimaryCandidate()) tracksPrim++;
 
     if (fFilterBit && !aodtrack->TestFilterBit(fFilterBit)) {
-      delete trackCopy;
+      //delete trackCopy;
       continue;
     }
 
     if (fFilterMask && !aodtrack->TestFilterBit(fFilterMask)) {
-      delete trackCopy;
+      //delete trackCopy;
       continue;
     }
 
@@ -543,39 +573,44 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoEvent(AliFemtoEvent *tEvent)
 
     //counting particles to set multiplicity
     if(fEstEventMult==kGlobalCount){
-      double impact[2];
-      double covimpact[3];
-      AliAODTrack* trk_clone = (AliAODTrack*)aodtrack->Clone("trk_clone");
-      if (trk_clone->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),10000,impact,covimpact)) {
-        if(impact[0]<0.2 && TMath::Abs(impact[1]+fV1[2])<2.0)
-          //if (aodtrack->IsPrimaryCandidate()) //? instead of kinks?
-          if (aodtrack->Chi2perNDF() < 4.0)
-            if (aodtrack->Pt() > 0.15 && aodtrack->Pt() < 20)
-              if (aodtrack->GetTPCNcls() > 70)
-                if (aodtrack->Eta() < 0.8)
-                  tNormMult++;
-      }
+      AliAODTrack* trk_clone = (AliAODTrack*)aodtrack->Clone("trk_clone"); //no DCA cut for global count
+      //if (aodtrack->IsPrimaryCandidate()) //? instead of kinks?
+      if (aodtrack->Chi2perNDF() < 4.0)
+       if (aodtrack->Pt() > 0.15 && aodtrack->Pt() < 20)
+         if (aodtrack->GetTPCNcls() > 70)
+           if (aodtrack->Eta() < 0.8)
+             tNormMult++;
       delete trk_clone;
     }
 
-    CopyAODtoFemtoTrack(aodtrack, trackCopy);
+    trackCopy = CopyAODtoFemtoTrack(aodtrack);
 
     // copying PID information from the correspondent track
     // const AliAODTrack *aodtrackpid = fEvent->GetTrack(labels[-1-fEvent->GetTrack(i)->GetID()]);
 
 
     AliAODTrack *aodtrackpid;
-    if((fFilterBit ==  (1 << (7))) || fFilterMask==128) //for TPC Only tracks we have to copy PID information from corresponding global tracks
-      aodtrackpid = fEvent->GetTrack(labels[-1-fEvent->GetTrack(i)->GetID()]);
-    else
-      aodtrackpid = fEvent->GetTrack(i);
+    if((fFilterBit ==  (1 << (7))) || fFilterMask == 128) {//for TPC Only tracks we have to copy PID information from corresponding global tracks
+      aodtrackpid = dynamic_cast<AliAODTrack*>(fEvent->GetTrack(labels[-1-fEvent->GetTrack(i)->GetID()]));
+    }
+    else {
+      aodtrackpid = dynamic_cast<AliAODTrack*>(fEvent->GetTrack(i));
+    }
+    assert(aodtrackpid&&"Not a standard AOD");
+
     CopyPIDtoFemtoTrack(aodtrackpid, trackCopy);
 
     if (mcP) {
       // Fill the hidden information with the simulated data
       //         Int_t pLabel = aodtrack->GetLabel();
-      AliAODMCParticle *tPart = GetParticleWithLabel(mcP, (TMath::Abs(aodtrack->GetLabel())));
-
+      //      AliAODMCParticle *tPart = GetParticleWithLabel(mcP, (TMath::Abs(aodtrack->GetLabel())));
+      AliAODMCParticle *tPart;
+      if(aodtrack->GetLabel() > -1 ) {
+       tPart = (AliAODMCParticle*)mcP->At(aodtrack->GetLabel());
+      }
+      else {
+       tPart = NULL;
+      }
       AliFemtoModelGlobalHiddenInfo *tInfo = new AliFemtoModelGlobalHiddenInfo();
       double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0;
       if (!tPart) {
@@ -613,10 +648,13 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoEvent(AliFemtoEvent *tEvent)
         //       fpt *= 1e13;
 
         //      cout << "Looking for mother ids " << endl;
-        if (motherids[TMath::Abs(aodtrack->GetLabel())]>0) {
-          //   cout << "Got mother id" << endl;
-          AliAODMCParticle *mother = GetParticleWithLabel(mcP, motherids[TMath::Abs(aodtrack->GetLabel())]);
-          // Check if this is the same particle stored twice on the stack
+
+        //if (motherids[TMath::Abs(aodtrack->GetLabel())]>0) {
+        if(tPart->GetMother() > -1) { //MC particle has a mother
+       //      cout << "Got mother id" << endl;
+         //          AliAODMCParticle *mother = GetParticleWithLabel(mcP, motherids[TMath::Abs(aodtrack->GetLabel())]);
+          AliAODMCParticle *mother = (AliAODMCParticle*)mcP->At(tPart->GetMother());
+         // Check if this is the same particle stored twice on the stack
           if (mother) {
             if ((mother->GetPdgCode() == tPart->GetPdgCode() || (mother->Px() == tPart->Px()))) {
               // It is the same particle
@@ -637,6 +675,9 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoEvent(AliFemtoEvent *tEvent)
               //             fpt = mother->T() *1e13*3e10;
 
             }
+           else { //particle's mother exists and the information about it can be added to hiddeninfo:
+             tInfo->SetMotherPdgCode(mother->GetPdgCode());
+           }
           }
         }
 
@@ -791,7 +832,7 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoEvent(AliFemtoEvent *tEvent)
     tEvent->SetNormalizedMult(multV0);
   }
 
-  if (mcP) delete [] motherids;
+  // if (mcP) delete [] motherids;
 
   // cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
 
@@ -806,28 +847,61 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoEvent(AliFemtoEvent *tEvent)
       if(aodv0->GetCharge()!=0) continue;
       if(aodv0->ChargeProng(0)==aodv0->ChargeProng(1)) continue;
       if(aodv0->CosPointingAngle(fV1)<0.998) continue;
-      AliFemtoV0* trackCopyV0 = new AliFemtoV0();
-      count_pass++;
-      CopyAODtoFemtoV0(aodv0, trackCopyV0);
+
+      AliAODTrack* daughterTrackPos = (AliAODTrack*)aodv0->GetDaughter(0); //getting positive daughter track
+      AliAODTrack* daughterTrackNeg = (AliAODTrack*)aodv0->GetDaughter(1); //getting negative daughter track
+      if(!daughterTrackPos) continue; //daughter tracks must exist
+      if(!daughterTrackNeg) continue;
+      if(daughterTrackNeg->Charge() == daughterTrackPos->Charge() ) continue; //and have different charge
+
+      AliFemtoV0* trackCopyV0 = CopyAODtoFemtoV0(aodv0);
+      if(mcP) {
+       daughterTrackPos->SetAODEvent(fEvent);
+       daughterTrackNeg->SetAODEvent(fEvent);
+       if(daughterTrackPos->GetLabel() > 0 && daughterTrackNeg->GetLabel() > 0 ) {
+         AliAODMCParticle* mcParticlePos = (AliAODMCParticle*)mcP->At(daughterTrackPos->GetLabel());
+         AliAODMCParticle* mcParticleNeg = (AliAODMCParticle*)mcP->At(daughterTrackNeg->GetLabel() );
+         if((mcParticlePos!=NULL) && (mcParticleNeg!=NULL)){
+           int motherOfPosID = mcParticlePos->GetMother();
+           int motherOfNegID = mcParticleNeg->GetMother();
+           if ((motherOfPosID > -1) && (motherOfPosID == motherOfNegID)){
+             AliFemtoModelHiddenInfo *tInfo = new AliFemtoModelHiddenInfo();
+             // Both daughter tracks refer to the same mother, we can continue
+             AliAODMCParticle *v0 = (AliAODMCParticle*)mcP->At(motherOfPosID); //our V0 particle
+
+             tInfo->SetPDGPid(v0->GetPdgCode());
+             int v0MotherId = v0->GetMother();
+             if(v0MotherId>-1) { //V0 particle has a mother
+               AliAODMCParticle* motherOfV0 = (AliAODMCParticle*)mcP->At(v0MotherId);
+               tInfo->SetMotherPdgCode(motherOfV0->GetPdgCode());
+             }
+             trackCopyV0->SetHiddenInfo(tInfo);
+           }
+         }
+       }
+      }
       tEvent->V0Collection()->push_back(trackCopyV0);
+      count_pass++;
       //cout<<"Pushback v0 to v0collection"<<endl;
     }
   }
 
+  return tEvent;
 }
 
-void AliFemtoEventReaderAOD::CopyAODtoFemtoTrack(AliAODTrack *tAodTrack,
-                                                 AliFemtoTrack *tFemtoTrack
+AliFemtoTrack* AliFemtoEventReaderAOD::CopyAODtoFemtoTrack(AliAODTrack *tAodTrack
                                                  //                                             AliPWG2AODTrack *tPWG2AODTrack
   )
 {
   // Copy the track information from the AOD into the internal AliFemtoTrack
   // If it exists, use the additional information from the PWG2 AOD
+  AliFemtoTrack *tFemtoTrack = new AliFemtoTrack();
 
   // Primary Vertex position
 
   fEvent->GetPrimaryVertex()->GetPosition(fV1);
   //  fEvent->GetPrimaryVertex()->GetXYZ(fV1);
+  tFemtoTrack->SetPrimaryVertex(fV1);
 
   tFemtoTrack->SetCharge(tAodTrack->Charge());
 
@@ -853,23 +927,28 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoTrack(AliAODTrack *tAodTrack,
   tAodTrack->GetCovMatrix(covmat);
 
   if (!fDCAglobalTrack) {
-    double impact[2];
-    double covimpact[3];
 
-    AliAODTrack* trk_clone = (AliAODTrack*)tAodTrack->Clone("trk_clone");
-    // if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
-    // if (!tAodTrack->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),10000,impact,covimpact)) {
-    if (!trk_clone->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),10000,impact,covimpact)) {
-      //cout << "sth went wrong with dca propagation" << endl;
-      tFemtoTrack->SetImpactD(-1000.0);
-      tFemtoTrack->SetImpactZ(-1000.0);
+    tFemtoTrack->SetImpactD(tAodTrack->DCA());
+    tFemtoTrack->SetImpactZ(tAodTrack->ZAtDCA());
 
-    }
-    else {
-      tFemtoTrack->SetImpactD(impact[0]);
-      tFemtoTrack->SetImpactZ(impact[1]);
-    }
-    delete trk_clone;
+
+    // double impact[2];
+    // double covimpact[3];
+
+    // AliAODTrack* trk_clone = (AliAODTrack*)tAodTrack->Clone("trk_clone");
+    // // if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
+    // // if (!tAodTrack->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),10000,impact,covimpact)) {
+    // if (!trk_clone->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),10000,impact,covimpact)) {
+    //   //cout << "sth went wrong with dca propagation" << endl;
+    //   tFemtoTrack->SetImpactD(-1000.0);
+    //   tFemtoTrack->SetImpactZ(-1000.0);
+
+    // }
+    // else {
+    //   tFemtoTrack->SetImpactD(impact[0]);
+    //   tFemtoTrack->SetImpactZ(impact[1]);
+    // }
+    // delete trk_clone;
 
   }
 
@@ -967,11 +1046,13 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoTrack(AliAODTrack *tAodTrack,
     tFemtoTrack->SetITSHitOnLayer(ii,tAodTrack->HasPointOnITSLayer(ii));
   }
 
-
+  return tFemtoTrack;
 }
 
-void AliFemtoEventReaderAOD::CopyAODtoFemtoV0(AliAODv0 *tAODv0, AliFemtoV0 *tFemtoV0)
+AliFemtoV0* AliFemtoEventReaderAOD::CopyAODtoFemtoV0(AliAODv0 *tAODv0 )
 {
+  AliFemtoV0 *tFemtoV0 = new AliFemtoV0();
+
   tFemtoV0->SetdecayLengthV0(tAODv0->DecayLength(fV1));
   tFemtoV0->SetdecayVertexV0X(tAODv0->DecayVertexV0X());
   tFemtoV0->SetdecayVertexV0Y(tAODv0->DecayVertexV0Y());
@@ -1000,6 +1081,7 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoV0(AliAODv0 *tAODv0, AliFemtoV0 *tFem
   //void SetTrackTopologyMapPos(unsigned int word, const unsigned long& m);
   //void SetTrackTopologyMapNeg(unsigned int word, const unsigned long& m);
 
+
   tFemtoV0->SetmomV0X(tAODv0->MomV0X());
   tFemtoV0->SetmomV0Y(tAODv0->MomV0Y());
   tFemtoV0->SetmomV0Z(tAODv0->MomV0Z());
@@ -1038,6 +1120,7 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoV0(AliAODv0 *tAODv0, AliFemtoV0 *tFem
   tFemtoV0->SetCosPointingAngle(tAODv0->CosPointingAngle(fV1));
   //tFemtoV0->SetYV0(tAODv0->Y());
 
+
   //void SetdedxNeg(float x);
   //void SeterrdedxNeg(float x);//Gael 04Fev2002
   //void SetlendedxNeg(float x);//Gael 04Fev2002
@@ -1051,6 +1134,7 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoV0(AliAODv0 *tAODv0, AliFemtoV0 *tFem
   AliAODTrack *trackpos = (AliAODTrack*)tAODv0->GetDaughter(0);
   AliAODTrack *trackneg = (AliAODTrack*)tAODv0->GetDaughter(1);
 
+
   if(trackpos && trackneg)
   {
     tFemtoV0->SetEtaPos(trackpos->Eta());
@@ -1105,6 +1189,7 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoV0(AliAODv0 *tAODv0, AliFemtoV0 *tFem
     tmpVec.SetX(tpcExitNeg[0]); tmpVec.SetY(tpcExitNeg[1]); tmpVec.SetZ(tpcExitNeg[2]);
     tFemtoV0->SetNominalTpcExitPointNeg(tmpVec);
 
+
     AliFemtoThreeVector vecTpcPos[9];
     AliFemtoThreeVector vecTpcNeg[9];
     for(int i=0;i<9;i++)
@@ -1121,9 +1206,22 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoV0(AliAODv0 *tAODv0, AliFemtoV0 *tFem
     tFemtoV0->SetdedxPos(trackpos->GetTPCsignal());
     tFemtoV0->SetdedxNeg(trackneg->GetTPCsignal());
 
-    if((tFemtoV0->StatusPos()&AliESDtrack::kTOFpid)==0 || (tFemtoV0->StatusPos()&AliESDtrack::kTIME)==0 || (tFemtoV0->StatusPos()&AliESDtrack::kTOFout)==0 )
+
+    Float_t probMisPos = 1.0;
+    Float_t probMisNeg = 1.0;
+
+    if (tFemtoV0->StatusPos() & AliESDtrack::kTOFout & AliESDtrack::kTIME) {  //AliESDtrack::kTOFpid=0x8000
+      probMisPos = fAODpidUtil->GetTOFMismatchProbability(trackpos);
+    }
+    if (tFemtoV0->StatusNeg() & AliESDtrack::kTOFout & AliESDtrack::kTIME) {  //AliESDtrack::kTOFpid=0x8000
+      probMisNeg = fAODpidUtil->GetTOFMismatchProbability(trackneg);
+    }
+
+    if(// (tFemtoV0->StatusPos()& AliESDtrack::kTOFpid)==0 ||
+       (tFemtoV0->StatusPos()&AliESDtrack::kTIME)==0 || (tFemtoV0->StatusPos()&AliESDtrack::kTOFout)==0 || probMisPos > 0.01)
     {
-      if((tFemtoV0->StatusNeg()&AliESDtrack::kTOFpid)==0 || (tFemtoV0->StatusNeg()&AliESDtrack::kTIME)==0 || (tFemtoV0->StatusNeg()&AliESDtrack::kTOFout)==0 )
+      if(// (tFemtoV0->StatusNeg()&AliESDtrack::kTOFpid)==0 ||
+         (tFemtoV0->StatusNeg()&AliESDtrack::kTIME)==0 || (tFemtoV0->StatusNeg()&AliESDtrack::kTOFout)==0 || probMisNeg > 0.01)
            {
              tFemtoV0->SetPosNSigmaTOFK(-1000);
              tFemtoV0->SetNegNSigmaTOFK(-1000);
@@ -1142,15 +1240,20 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoV0(AliAODv0 *tAODv0, AliFemtoV0 *tFem
     }
     else
     {
-      tFemtoV0->SetPosNSigmaTOFK(fAODpidUtil->NumberOfSigmasTOF(trackpos,AliPID::kKaon));
-      tFemtoV0->SetNegNSigmaTOFK(fAODpidUtil->NumberOfSigmasTOF(trackneg,AliPID::kKaon));
-      tFemtoV0->SetPosNSigmaTOFP(fAODpidUtil->NumberOfSigmasTOF(trackpos,AliPID::kProton));
-      tFemtoV0->SetNegNSigmaTOFP(fAODpidUtil->NumberOfSigmasTOF(trackneg,AliPID::kProton));
-      tFemtoV0->SetPosNSigmaTOFPi(fAODpidUtil->NumberOfSigmasTOF(trackpos,AliPID::kPion));
-      tFemtoV0->SetNegNSigmaTOFPi(fAODpidUtil->NumberOfSigmasTOF(trackneg,AliPID::kPion));
-
+      if(trackpos->IsOn(AliESDtrack::kTOFout & AliESDtrack::kTIME)) {
+        tFemtoV0->SetPosNSigmaTOFK(fAODpidUtil->NumberOfSigmasTOF(trackpos,AliPID::kKaon));
+        tFemtoV0->SetPosNSigmaTOFP(fAODpidUtil->NumberOfSigmasTOF(trackpos,AliPID::kProton));
+        tFemtoV0->SetPosNSigmaTOFPi(fAODpidUtil->NumberOfSigmasTOF(trackpos,AliPID::kPion));
+      }
+      if(trackneg->IsOn(AliESDtrack::kTOFout & AliESDtrack::kTIME)) {
+        tFemtoV0->SetNegNSigmaTOFK(fAODpidUtil->NumberOfSigmasTOF(trackneg,AliPID::kKaon));
+        tFemtoV0->SetNegNSigmaTOFP(fAODpidUtil->NumberOfSigmasTOF(trackneg,AliPID::kProton));
+        tFemtoV0->SetNegNSigmaTOFPi(fAODpidUtil->NumberOfSigmasTOF(trackneg,AliPID::kPion));
+      }
       double TOFSignalPos = trackpos->GetTOFsignal();
       double TOFSignalNeg = trackneg->GetTOFsignal();
+      TOFSignalPos -= fAODpidUtil->GetTOFResponse().GetStartTime(trackpos->P());
+      TOFSignalNeg -= fAODpidUtil->GetTOFResponse().GetStartTime(trackneg->P());
       double pidPos[5];
       double pidNeg[5];
       trackpos->GetIntegratedTimes(pidPos);
@@ -1163,13 +1266,17 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoV0(AliAODv0 *tAODv0, AliFemtoV0 *tFem
       tFemtoV0->SetTOFKaonTimeNeg(TOFSignalNeg-pidNeg[3]);
       tFemtoV0->SetTOFProtonTimeNeg(TOFSignalNeg-pidNeg[4]);
     }
+
   }
   else
   {
+
     tFemtoV0->SetStatusPos(999);
     tFemtoV0->SetStatusNeg(999);
   }
+
   tFemtoV0->SetOnFlyStatusV0(tAODv0->GetOnFlyStatus());
+  return tFemtoV0;
 }
 
 void AliFemtoEventReaderAOD::SetFilterBit(UInt_t ibit)
@@ -1235,27 +1342,60 @@ void AliFemtoEventReaderAOD::CopyPIDtoFemtoTrack(AliAODTrack *tAodTrack,
 {
 
   if (fDCAglobalTrack) {
-    // // copying DCA information (taking it from global tracks gives better resolution than from TPC-only)
 
-    double impact[2];
-    double covimpact[3];
+    // code from Michael and Prabhat from AliAnalysisTaskDptDptCorrelations
+    const AliAODVertex* vertex = (AliAODVertex*) fEvent->GetPrimaryVertex();
+    float vertexX  = -999.;
+    float vertexY  = -999.;
+    float vertexZ  = -999.;
+
+    if(vertex) {
+      Double32_t fCov[6];
+      vertex->GetCovarianceMatrix(fCov);
+      if(vertex->GetNContributors() > 0) {
+        if(fCov[5] != 0) {
+          vertexX = vertex->GetX();
+          vertexY = vertex->GetY();
+          vertexZ = vertex->GetZ();
 
-    AliAODTrack* trk_clone = (AliAODTrack*)tAodTrack->Clone("trk_clone");
-    // if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
-    // if (!tAodTrack->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),10000,impact,covimpact)) {
-    if (!trk_clone->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),10000,impact,covimpact)) {
+        }
+      }
+    }
 
-    // if (!tAodTrack->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),10000,impact,covimpact)) {
-      //cout << "sth went wrong with dca propagation" << endl;
-      tFemtoTrack->SetImpactD(-1000.0);
-      tFemtoTrack->SetImpactZ(-1000.0);
+    Double_t pos[3];
+    tAodTrack->GetXYZ(pos);
 
-    }
-    else {
-      tFemtoTrack->SetImpactD(impact[0]);
-      tFemtoTrack->SetImpactZ(impact[1]);
-    }
-    delete trk_clone;
+    Double_t DCAX = pos[0] - vertexX;
+    Double_t DCAY = pos[1] - vertexY;
+    Double_t DCAZ = pos[2] - vertexZ;
+
+    Double_t DCAXY = TMath::Sqrt((DCAX*DCAX) + (DCAY*DCAY));
+
+    tFemtoTrack->SetImpactD(DCAXY);
+    tFemtoTrack->SetImpactZ(DCAZ);
+
+
+    // // // copying DCA information (taking it from global tracks gives better resolution than from TPC-only)
+
+    // double impact[2];
+    // double covimpact[3];
+
+    // AliAODTrack* trk_clone = (AliAODTrack*)tAodTrack->Clone("trk_clone");
+    // // if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
+    // // if (!tAodTrack->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),10000,impact,covimpact)) {
+    // if (!trk_clone->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),10000,impact,covimpact)) {
+
+    // // if (!tAodTrack->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),10000,impact,covimpact)) {
+    //   //cout << "sth went wrong with dca propagation" << endl;
+    //   tFemtoTrack->SetImpactD(-1000.0);
+    //   tFemtoTrack->SetImpactZ(-1000.0);
+
+    // }
+    // else {
+    //   tFemtoTrack->SetImpactD(impact[0]);
+    //   tFemtoTrack->SetImpactZ(impact[1]);
+    // }
+    // delete trk_clone;
   }
 
   double aodpid[10];
@@ -1276,7 +1416,7 @@ void AliFemtoEventReaderAOD::CopyPIDtoFemtoTrack(AliAODTrack *tAodTrack,
   Float_t probMis = 1.0;
 
   //what is that code? for what do we need it? nsigma values are not enough?
-  if (tAodTrack->GetStatus() & AliESDtrack::kTOFpid) {  //AliESDtrack::kTOFpid=0x8000
+  if (tAodTrack->GetStatus() & AliESDtrack::kTOFout & AliESDtrack::kTIME) {  //AliESDtrack::kTOFpid=0x8000
     tTOF = tAodTrack->GetTOFsignal();
     tAodTrack->GetIntegratedTimes(aodpid);
 
@@ -1294,6 +1434,7 @@ void AliFemtoEventReaderAOD::CopyPIDtoFemtoTrack(AliAODTrack *tAodTrack,
   float nsigmaTPCK=-1000.;
   float nsigmaTPCPi=-1000.;
   float nsigmaTPCP=-1000.;
+  float nsigmaTPCE=-1000.;
 
   //   cout<<"in reader fESDpid"<<fESDpid<<endl;
 
@@ -1301,11 +1442,13 @@ void AliFemtoEventReaderAOD::CopyPIDtoFemtoTrack(AliAODTrack *tAodTrack,
     nsigmaTPCK = fAODpidUtil->NumberOfSigmasTPC(tAodTrack,AliPID::kKaon);
     nsigmaTPCPi = fAODpidUtil->NumberOfSigmasTPC(tAodTrack,AliPID::kPion);
     nsigmaTPCP = fAODpidUtil->NumberOfSigmasTPC(tAodTrack,AliPID::kProton);
+    nsigmaTPCE = fAODpidUtil->NumberOfSigmasTPC(tAodTrack,AliPID::kElectron);
   }
 
   tFemtoTrack->SetNSigmaTPCPi(nsigmaTPCPi);
   tFemtoTrack->SetNSigmaTPCK(nsigmaTPCK);
   tFemtoTrack->SetNSigmaTPCP(nsigmaTPCP);
+  tFemtoTrack->SetNSigmaTPCE(nsigmaTPCE);
 
   tFemtoTrack->SetTPCchi2(tAodTrack->Chi2perNDF());
   tFemtoTrack->SetTPCncls(tAodTrack->GetTPCNcls());
@@ -1321,18 +1464,21 @@ void AliFemtoEventReaderAOD::CopyPIDtoFemtoTrack(AliAODTrack *tAodTrack,
   float nsigmaTOFPi=-1000.;
   float nsigmaTOFK=-1000.;
   float nsigmaTOFP=-1000.;
+  float nsigmaTOFE=-1000.;
 
-  if ((tAodTrack->GetStatus() & AliESDtrack::kTOFpid) && //AliESDtrack::kTOFpid=0x8000
+  if (// (tAodTrack->GetStatus() & AliESDtrack::kTOFpid) &&
+      //AliESDtrack::kTOFpid=0x8000
       (tAodTrack->GetStatus() & AliESDtrack::kTOFout) && //AliESDtrack::kTOFout=0x2000
       (tAodTrack->GetStatus() & AliESDtrack::kTIME) //AliESDtrack::kTIME=0x80000000
       && probMis < 0.01) // TOF mismatch probaility
   {
-    if(tAodTrack->IsOn(AliESDtrack::kTOFpid)) //AliESDtrack::kTOFpid=0x8000
+    if(tAodTrack->IsOn(AliESDtrack::kTOFout & AliESDtrack::kTIME)) //AliESDtrack::kTOFpid=0x8000
          {
 
            nsigmaTOFPi = fAODpidUtil->NumberOfSigmasTOF(tAodTrack,AliPID::kPion);
            nsigmaTOFK = fAODpidUtil->NumberOfSigmasTOF(tAodTrack,AliPID::kKaon);
            nsigmaTOFP = fAODpidUtil->NumberOfSigmasTOF(tAodTrack,AliPID::kProton);
+           nsigmaTOFE = fAODpidUtil->NumberOfSigmasTOF(tAodTrack,AliPID::kElectron);
 
            Double_t len=200;// esdtrack->GetIntegratedLength(); !!!!!
            Double_t tof=tAodTrack->GetTOFsignal();
@@ -1343,6 +1489,7 @@ void AliFemtoEventReaderAOD::CopyPIDtoFemtoTrack(AliAODTrack *tAodTrack,
   tFemtoTrack->SetNSigmaTOFPi(nsigmaTOFPi);
   tFemtoTrack->SetNSigmaTOFK(nsigmaTOFK);
   tFemtoTrack->SetNSigmaTOFP(nsigmaTOFP);
+  tFemtoTrack->SetNSigmaTOFE(nsigmaTOFE);
 
 
   //////////////////////////////////////
@@ -1462,6 +1609,18 @@ void AliFemtoEventReaderAOD::SetpA2013(Bool_t pa2013)
   fpA2013 = pa2013;
 }
 
+void AliFemtoEventReaderAOD::SetUseMVPlpSelection(Bool_t mvplp)
+{
+  fMVPlp = mvplp;
+}
+
+void AliFemtoEventReaderAOD::SetIsPileUpEvent(Bool_t ispileup)
+{
+  fisPileUp = ispileup;
+}
+
+
+
 void AliFemtoEventReaderAOD::SetDCAglobalTrack(Bool_t dcagt)
 {
   fDCAglobalTrack = dcagt;
@@ -1481,7 +1640,7 @@ bool AliFemtoEventReaderAOD::RejectEventCentFlat(float MagField, float CentPerce
                              {.828,.793,.776,.772,.775,.796,.788,.804,.839}};
   int weightBinCent = (int) CentPercent;
   if(fRandomNumber->Rndm() > kCentWeight[weightBinSign][weightBinCent]) RejectEvent = kTRUE;
-
+  delete fRandomNumber;
   return RejectEvent;
 }