]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGHF/vertexingHF/charmFlow/AliAnalysisTaskSEHFv2.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff
? search:
summary | shortlog | log | commit | commitdiff | tree
raw | inline | side by side
Mem.leak fix: Get(...) only once per object
[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / charmFlow / AliAnalysisTaskSEHFv2.cxx
diff --git a/PWGHF/vertexingHF/charmFlow/AliAnalysisTaskSEHFv2.cxx b/PWGHF/vertexingHF/charmFlow/AliAnalysisTaskSEHFv2.cxx
index 86c786582dd38bf05f3a9870383f7d9522f13d25..dc2f701f92ad0e7225f912c7e30cd9ec58bed797 100644 (file)
--- a/PWGHF/vertexingHF/charmFlow/AliAnalysisTaskSEHFv2.cxx
+++ b/PWGHF/vertexingHF/charmFlow/AliAnalysisTaskSEHFv2.cxx
@@ -96,7 +96,9 @@ AliAnalysisTaskSE(),
   fV0EPorder(2),
   fMinCentr(20),
   fMaxCentr(80),
-  fEtaGap(kFALSE)
+  fEtaGap(kFALSE),
+  fSubEvents(2),
+  fCentBinSizePerMil(25)
 {
   // Default constructor
 }
@@ -120,7 +122,9 @@ AliAnalysisTaskSEHFv2::AliAnalysisTaskSEHFv2(const char *name,AliRDHFCuts *rdCut
   fV0EPorder(2),
   fMinCentr(20),
   fMaxCentr(80),
-  fEtaGap(kFALSE)
+  fEtaGap(kFALSE),
+  fSubEvents(2),
+  fCentBinSizePerMil(25)
 {
   // standard constructor
   Int_t pdg=421;
@@ -145,7 +149,7 @@ AliAnalysisTaskSEHFv2::AliAnalysisTaskSEHFv2(const char *name,AliRDHFCuts *rdCut
     break;
   }
   fAfterBurner = new AliHFAfterBurner(fDecChannel);
-  if(pdg==413) SetMassLimits((Float_t)0.135,(Float_t)0.165);
+  if(pdg==413) SetMassLimits((Float_t)0.1,(Float_t)0.2);
   else SetMassLimits((Float_t)0.2,pdg); //check range
   fNPtBins=fRDCuts->GetNPtBins();
 
@@ -259,62 +263,79 @@ void AliAnalysisTaskSEHFv2::UserCreateOutputObjects()
   fOutput->SetOwner();
   fOutput->SetName("MainOutput");
   
-  for(Int_t icentr=fMinCentr+5;icentr<=fMaxCentr;icentr=icentr+5){
-    TString centrname;centrname.Form("centr%d_%d",icentr-5,icentr);
+  for(Int_t icentr=fMinCentr*10+fCentBinSizePerMil;icentr<=fMaxCentr*10;icentr=icentr+fCentBinSizePerMil){
+    TString centrname;centrname.Form("centr%d_%d",icentr-fCentBinSizePerMil,icentr);
 
-      TH2F* hMPtCand=new TH2F(Form("hMPtCand%s",centrname.Data()),Form("Mass vs pt %s;pt (GeV);M (GeV/c^{2})",centrname.Data()),120,0,24.,fNMassBins,fLowmasslimit,fUpmasslimit);
+    TH2F* hMPtCand=new TH2F(Form("hMPtCand%s",centrname.Data()),Form("Mass vs pt %s;pt (GeV);M (GeV/c^{2})",centrname.Data()),120,0,24.,fNMassBins,fLowmasslimit,fUpmasslimit);
       fOutput->Add(hMPtCand);//For <pt> calculation
 
 
     //Candidate distributions
+
     for(Int_t i=0;i<fNPtBins;i++){ 
-      TH2F* hMc2phi=new TH2F(Form("hMc2phi_pt%d%s",i,centrname.Data()),Form("Mass vs cos2#Delta#phi (p_{t} bin %d %s);cos2#Delta#phi;M (GeV/c^{2})",i,centrname.Data()),100,-1.,1.,fNMassBins,fLowmasslimit,fUpmasslimit);
-      fOutput->Add(hMc2phi);//for 2D analysis
-      
-      TH2F* hMphi=new TH2F(Form("hMphi_pt%d%s",i,centrname.Data()),Form("Mass vs #Delta#phi %s;#Delta#phi;M (GeV/c^{2})",centrname.Data()),96,0,TMath::Pi(),fNMassBins,fLowmasslimit,fUpmasslimit);
-      fOutput->Add(hMphi);//for phi bins analysis
+
+      // Delta Phi histograms
+      Double_t maxphi = TMath::Pi();
+      if (fDecChannel == 2) maxphi = TMath::PiOver2();
+      TH2F* hMdeltaphi=new TH2F(Form("hMdeltaphi_pt%d%s",i,centrname.Data()),Form("Mass vs #Delta#phi %s;#Delta#phi;M (GeV/c^{2})",centrname.Data()),96,0,maxphi,fNMassBins,fLowmasslimit,fUpmasslimit);
+      fOutput->Add(hMdeltaphi);//for phi bins analysis
+      TH2F* hMc2deltaphi=new TH2F(Form("hMc2deltaphi_pt%d%s",i,centrname.Data()),Form("Mass vs cos2#Delta#phi (p_{t} bin %d %s);cos2#Delta#phi;M (GeV/c^{2})",i,centrname.Data()),100,-1.,1.,fNMassBins,fLowmasslimit,fUpmasslimit);
+      fOutput->Add(hMc2deltaphi);
+      TH2F* hMs2deltaphi=new TH2F(Form("hMs2deltaphi_pt%d%s",i,centrname.Data()),Form("Mass vs sin2#Delta#phi (p_{t} bin %d %s);sin2#Delta#phi;M (GeV/c^{2})",i,centrname.Data()),100,-1.,1.,fNMassBins,fLowmasslimit,fUpmasslimit);
+      fOutput->Add(hMs2deltaphi);
       
+      // phi histograms (for systematics)
+      TH2F* hCos2PhiMass=new TH2F(Form("hCos2phiMass_pt%d%s",i,centrname.Data()),Form("Mass vs cos(2#phi) %s;cos(2#phi);M (GeV/c^{2})",centrname.Data()),100,-1.,1.,fNMassBins,fLowmasslimit,fUpmasslimit);
+      fOutput->Add(hCos2PhiMass);
+      TH2F* hSin2PhiMass=new TH2F(Form("hSin2phiMass_pt%d%s",i,centrname.Data()),Form("Mass vs sin(2#phi) %s;sin(2#phi);M (GeV/c^{2})",centrname.Data()),100,-1.,1.,fNMassBins,fLowmasslimit,fUpmasslimit);
+      fOutput->Add(hSin2PhiMass);
+
+      // Histos using MC truth
       if (fReadMC){
-       TH2F* hMc2phiS=new TH2F(Form("hMc2phiS_pt%d%s",i,centrname.Data()),Form("Mass vs cos2#Delta#phi (p_{t} bin %d %s);cos2#Delta#phi;M (GeV/c^{2})",i,centrname.Data()),100,-1.,1.,fNMassBins,fLowmasslimit,fUpmasslimit);
-       fOutput->Add(hMc2phiS);
-       TH2F * hMphiS=new TH2F(Form("hMphiS_pt%d%s",i,centrname.Data()),Form("Mass vs #Delta#phi (p_{t} bin %d %s);#Delta#phi;M (GeV/c^{2})",i,centrname.Data()),96,0,2*TMath::Pi(),fNMassBins,fLowmasslimit,fUpmasslimit);
-       fOutput->Add(hMphiS);
-       TH2F* hMc2phiB=new TH2F(Form("hMc2phiB_pt%d%s",i,centrname.Data()),Form("Mass vs cos2#Delta#phi (p_{t} bin %d %s);cos2#Delta#phi;M (GeV/c^{2})",i,centrname.Data()),100,-1.,1.,fNMassBins,fLowmasslimit,fUpmasslimit);
-       fOutput->Add(hMc2phiB);
-       TH2F * hMphiB=new TH2F(Form("hMphiB_pt%d%s",i,centrname.Data()),Form("Mass vs #Delta#phi (p_{t} bin %d %s);#Delta#phi;M (GeV/c^{2})",i,centrname.Data()),96,0,2*TMath::Pi(),fNMassBins,fLowmasslimit,fUpmasslimit);
-       fOutput->Add(hMphiB);
+       TH2F* hMc2deltaphiS=new TH2F(Form("hMc2deltaphiS_pt%d%s",i,centrname.Data()),Form("Mass vs cos2#Delta#phi (p_{t} bin %d %s);cos2#Delta#phi;M (GeV/c^{2})",i,centrname.Data()),100,-1.,1.,fNMassBins,fLowmasslimit,fUpmasslimit);
+       fOutput->Add(hMc2deltaphiS);
+       TH2F * hMdeltaphiS=new TH2F(Form("hMdeltaphiS_pt%d%s",i,centrname.Data()),Form("Mass vs #Delta#phi (p_{t} bin %d %s);#Delta#phi;M (GeV/c^{2})",i,centrname.Data()),96,0,2*TMath::Pi(),fNMassBins,fLowmasslimit,fUpmasslimit);
+       fOutput->Add(hMdeltaphiS);
+       TH2F* hMc2deltaphiB=new TH2F(Form("hMc2deltaphiB_pt%d%s",i,centrname.Data()),Form("Mass vs cos2#Delta#phi (p_{t} bin %d %s);cos2#Delta#phi;M (GeV/c^{2})",i,centrname.Data()),100,-1.,1.,fNMassBins,fLowmasslimit,fUpmasslimit);
+       fOutput->Add(hMc2deltaphiB);
+       TH2F * hMdeltaphiB=new TH2F(Form("hMdeltaphiB_pt%d%s",i,centrname.Data()),Form("Mass vs #Delta#phi (p_{t} bin %d %s);#Delta#phi;M (GeV/c^{2})",i,centrname.Data()),96,0,2*TMath::Pi(),fNMassBins,fLowmasslimit,fUpmasslimit);
+       fOutput->Add(hMdeltaphiB);
        if((fDecChannel != AliAnalysisTaskSEHFv2::kDplustoKpipi) &&(fDecChannel != AliAnalysisTaskSEHFv2::kDstartoKpipi)){
-         TH2F* hMc2phiR=new TH2F(Form("hMc2phiR_pt%d%s",i,centrname.Data()),Form("Mass vs cos2#Delta#phi (p_{t} bin %d %s);cos2#Delta#phi;M (GeV/c^{2})",i,centrname.Data()),100,-1.,1.,fNMassBins,fLowmasslimit,fUpmasslimit);
-         fOutput->Add(hMc2phiR);
-         TH2F* hMphiR=new TH2F(Form("hMphiR_pt%d%s",i,centrname.Data()),Form("Mass vs #Delta#phi (p_{t} bin %d %s);#Delta#phi;M (GeV/c^{2})",i,centrname.Data()),96,0,2*TMath::Pi(),fNMassBins,fLowmasslimit,fUpmasslimit);
-         fOutput->Add(hMphiR);
+         TH2F* hMc2deltaphiR=new TH2F(Form("hMc2deltaphiR_pt%d%s",i,centrname.Data()),Form("Mass vs cos2#Delta#phi (p_{t} bin %d %s);cos2#Delta#phi;M (GeV/c^{2})",i,centrname.Data()),100,-1.,1.,fNMassBins,fLowmasslimit,fUpmasslimit);
+         fOutput->Add(hMc2deltaphiR);
+         TH2F* hMdeltaphiR=new TH2F(Form("hMdeltaphiR_pt%d%s",i,centrname.Data()),Form("Mass vs #Delta#phi (p_{t} bin %d %s);#Delta#phi;M (GeV/c^{2})",i,centrname.Data()),96,0,2*TMath::Pi(),fNMassBins,fLowmasslimit,fUpmasslimit);
+         fOutput->Add(hMdeltaphiR);
        }
       }
     }
 
 
-    //Event Plane
+    // Event Plane
     TH2F* hEvPlane=new TH2F(Form("hEvPlane%s",centrname.Data()),Form("VZERO/TPC Event plane angle %s;#phi Ev Plane (TPC);#phi Ev Plane (VZERO);Entries",centrname.Data()),200,0.,TMath::Pi(),200,0.,TMath::Pi());
     fOutput->Add(hEvPlane);
-
     TH1F* hEvPlaneA=new TH1F(Form("hEvPlaneA%s",centrname.Data()),Form("Event plane angle %s;#phi Ev Plane;Entries",centrname.Data()),200,0.,TMath::Pi());
     fOutput->Add(hEvPlaneA);
-
     TH1F* hEvPlaneB=new TH1F(Form("hEvPlaneB%s",centrname.Data()),Form("Event plane angle %s;#phi Ev Plane;Entries",centrname.Data()),200,0.,TMath::Pi());
     fOutput->Add(hEvPlaneB);
-
     TH1F* hEvPlaneCand=new TH1F(Form("hEvPlaneCand%s",centrname.Data()),Form("Event plane angle - Event plane angle per candidate %s;#phi(Ev Plane Candidate);Entries",centrname.Data()),200,-TMath::Pi(),TMath::Pi());
     fOutput->Add(hEvPlaneCand);
 
+    // histos for EP resolution
     TH1F* hEvPlaneReso=new TH1F(Form("hEvPlaneReso%s",centrname.Data()),Form("Event plane angle Resolution %s;cos2(#psi_{A}-#psi_{B});Entries",centrname.Data()),220,-1.1,1.1);
     fOutput->Add(hEvPlaneReso);
-    if(fEventPlaneMeth>kTPCVZERO){
-      TH1F* hEvPlaneReso2=new TH1F(Form("hEvPlaneReso2%s",centrname.Data()),Form("Event plane angle Resolution %s;cos2(#psi_{A}-#psi_{B});Entries",centrname.Data()),220,-1.1,1.1);
+    if(fEventPlaneMeth>=kTPCVZERO){
+      TH1F* hEvPlaneReso2=new TH1F(Form("hEvPlaneReso2%s",centrname.Data()),Form("Event plane angle Resolution %s;cos2(#psi_{A}-#psi_{C});Entries",centrname.Data()),220,-1.1,1.1);
       fOutput->Add(hEvPlaneReso2);
-      TH1F* hEvPlaneReso3=new TH1F(Form("hEvPlaneReso3%s",centrname.Data()),Form("Event plane angle Resolution %s;cos2(#psi_{A}-#psi_{B});Entries",centrname.Data()),220,-1.1,1.1);
+      TH1F* hEvPlaneReso3=new TH1F(Form("hEvPlaneReso3%s",centrname.Data()),Form("Event plane angle Resolution %s;cos2(#psi_{B}-#psi_{C});Entries",centrname.Data()),220,-1.1,1.1);
       fOutput->Add(hEvPlaneReso3);
     }
-  }
+
+    // histos for EPsystematics    
+    TH1F *hCos2EP=new TH1F(Form("hCos2EP%s",centrname.Data()),Form("cos(2PsiEP) %s;cos2(#psi_{EP});Entries",centrname.Data()),100,-1.,1.);
+    fOutput->Add(hCos2EP);
+    TH1F *hSin2EP=new TH1F(Form("hSin2EP%s",centrname.Data()),Form("sin(2PsiEP) %s;sin2(#psi_{EP});Entries",centrname.Data()),100,-1.,1.);
+    fOutput->Add(hSin2EP);
+   }
   
   PostData(1,fhEventsInfo);
   PostData(2,fOutput);
@@ -425,6 +446,7 @@ void AliAnalysisTaskSEHFv2::UserExec(Option_t */*option*/)
   AliEventplane *pl=aod->GetEventplane();
   if(!pl){
     AliError("AliAnalysisTaskSEHFv2::UserExec:no eventplane! v2 analysis without eventplane not possible!\n");
+    fhEventsInfo->Fill(6);
     return;
   }
   
@@ -434,6 +456,7 @@ void AliAnalysisTaskSEHFv2::UserExec(Option_t */*option*/)
   Double_t rpangleVZERO=0;
   Double_t planereso=0;
   Double_t deltaPsi=0;
+  Double_t rpangleeventC=0;
   Double_t rpangleeventB=0;
   Double_t rpangleeventA=0;
 
@@ -444,14 +467,17 @@ void AliAnalysisTaskSEHFv2::UserExec(Option_t */*option*/)
   
   //determine centrality bin
   Float_t centr=fRDCuts->GetCentrality(aod);
-  Int_t icentr=0;
-  for(Int_t ic=fMinCentr+5;ic<=fMaxCentr;ic=ic+5){
-    if(ic>centr){
+  Float_t centrPerMil=centr*10.;
+  Int_t icentr=-1;
+  for(Int_t ic=fMinCentr*10+fCentBinSizePerMil;ic<=fMaxCentr*10;ic=ic+fCentBinSizePerMil){
+    if(ic>centrPerMil){
       icentr=ic;
       break;
     }
   }
-  TString centrbinname=Form("centr%d_%d",icentr-5,icentr);
+  if(icentr==-1) return;
+
+  TString centrbinname=Form("centr%d_%d",icentr-fCentBinSizePerMil,icentr);
 
   if(fReadMC){
     TRandom3 *g = new TRandom3(0);
@@ -459,60 +485,44 @@ void AliAnalysisTaskSEHFv2::UserExec(Option_t */*option*/)
     delete g;g=0x0;
     if(fUseAfterBurner)fAfterBurner->SetEventPlane((Double_t)eventplane);
   }else{
+   // TPC event plane
+    rpangleTPC = pl->GetEventplane("Q");
+    if(rpangleTPC<0){
+      fhEventsInfo->Fill(6);
+      return;
+    }
+    rpangleeventA = rpangleTPC;
+    if(fSubEvents==2||fEventPlaneMeth==kVZERO){
+      qsub1 = pl->GetQsub1();
+      qsub2 = pl->GetQsub2();
+      if(!qsub1 || !qsub2){
+       fhEventsInfo->Fill(6);
+       return;
+      }
+      rpangleeventA = qsub1->Phi()/2.;
+      rpangleeventB = qsub2->Phi()/2.;
+    }
     if(fEventPlaneMeth!=kTPC){
       //VZERO EP and resolution
       rpangleVZERO=GetPhi0Pi(pl->GetEventplane("V0",aod,fV0EPorder));
-      if(fEventPlaneMeth>kTPCVZERO){
-       //Using V0A/C for VZERO resolution
-       rpangleeventA= GetPhi0Pi(pl->GetEventplane("V0A",aod,fV0EPorder));
-       rpangleeventB= GetPhi0Pi(pl->GetEventplane("V0C",aod,fV0EPorder));
-       if(fEventPlaneMeth==kVZEROA)rpangleVZERO=rpangleeventA;
-       else if(fEventPlaneMeth==kVZEROC)rpangleVZERO=rpangleeventB;
-       deltaPsi =rpangleeventA-rpangleeventB;
-       eventplane=rpangleVZERO;
+      rpangleeventC=rpangleVZERO;
+      if(fEventPlaneMeth==kVZEROA||fEventPlaneMeth==kVZEROC||(fEventPlaneMeth==kTPCVZERO&&fSubEvents==3)){
+       rpangleeventB=GetPhi0Pi(pl->GetEventplane("V0A",aod,fV0EPorder));
+       rpangleeventC=GetPhi0Pi(pl->GetEventplane("V0C",aod,fV0EPorder));
+       if(fEventPlaneMeth==kVZEROA)rpangleVZERO=rpangleeventB;
+       else if(fEventPlaneMeth==kVZEROC)rpangleVZERO=rpangleeventC;
       }
     }
-    // TPC event plane
-    rpangleTPC = pl->GetEventplane("Q");
-    if(fEventPlaneMeth<=kTPCVZERO){
-      // TPC resolution 
+
+
+    if(fEventPlaneMeth>kTPCVZERO)eventplane=rpangleVZERO;
+    else{
       q = pl->GetQVector();
-      if(fEtaGap){
-       qsub1 = pl->GetQsub1();
-       qsub2 = pl->GetQsub2();
-       if(!qsub1 || !qsub2){
-         fhEventsInfo->Fill(6);
-         return;
-       }
-       rpangleeventA = qsub1->Phi()/2.;
-       rpangleeventB = qsub2->Phi()/2.;
-      }
-      deltaPsi = pl->GetQsubRes();
-      //       planereso = TMath::Cos(2.*deltaPsi); // reaction plane resolution
+      eventplane=rpangleTPC;
     }
+    deltaPsi =rpangleeventA-rpangleeventB;
   }
 
-  //verify TPC EP
-  Double_t rpsubTPC=rpangleTPC;
-  // if(fEventPlaneMeth==kVZEROpos||fEventPlaneMeth==kVZEROneg){
-  //   TVector2 *qsub=0x0;
-  //   if(fEventPlaneMeth==kVZEROneg)qsub=pl->GetQsub2();
-  //   else qsub=pl->GetQsub1();
-  //   if(!qsub){
-  //     fhEventsInfo->Fill(6);
-  //     return;
-  //   }
-  //   if(qsub->X()==0 && qsub->Y()==0){
-  //     fhEventsInfo->Fill(6);
-  //     return;
-  //   }
-  //   rpsubTPC=qsub->Phi()/2.;
-  // }
-  if(rpsubTPC<0){
-    fhEventsInfo->Fill(6);
-    return;
-  }
-  
   if(TMath::Abs(deltaPsi)>TMath::Pi()/2.){
     if(deltaPsi>0.) deltaPsi-=TMath::Pi();
     else deltaPsi +=TMath::Pi();
@@ -524,23 +534,28 @@ void AliAnalysisTaskSEHFv2::UserExec(Option_t */*option*/)
   //Filling EP-related histograms
   ((TH2F*)fOutput->FindObject(Form("hEvPlane%s",centrbinname.Data())))->Fill(rpangleTPC,rpangleVZERO); // reaction plane angle without autocorrelations removal
   ((TH1F*)fOutput->FindObject(Form("hEvPlaneReso%s",centrbinname.Data())))->Fill(planereso); //RP resolution   
+  ((TH1F*)fOutput->FindObject(Form("hCos2EP%s",centrbinname.Data())))->Fill(TMath::Cos(2*eventplane)); // syst check
+  ((TH1F*)fOutput->FindObject(Form("hSin2EP%s",centrbinname.Data())))->Fill(TMath::Sin(2*eventplane)); // syst check
+    
   if(fEventPlaneMeth>kTPCVZERO||fEtaGap){
     ((TH1F*)fOutput->FindObject(Form("hEvPlaneA%s",centrbinname.Data())))->Fill(rpangleeventA); //Angle of first subevent
     ((TH1F*)fOutput->FindObject(Form("hEvPlaneB%s",centrbinname.Data())))->Fill(rpangleeventB); //Angle of second subevent
   }
-  if(fEventPlaneMeth>kTPCVZERO){
-    Double_t deltaSub=rpsubTPC-rpangleeventA;
+  if(fEventPlaneMeth>kTPCVZERO||fSubEvents==3){
+    Double_t deltaSub=rpangleeventA-rpangleeventC;
     if(TMath::Abs(deltaSub)>TMath::Pi()/2.){// difference of subevents reaction plane angle cannot be bigger than phi/2
       if(deltaSub>0.) deltaSub-=TMath::Pi();
       else deltaSub +=TMath::Pi();
     } 
-    ((TH1F*)fOutput->FindObject(Form("hEvPlaneReso2%s",centrbinname.Data())))->Fill(TMath::Cos(2.*deltaSub)); //RP resolution   
-    deltaSub =rpsubTPC-rpangleeventB;
+    TH1F* htmp1=(TH1F*)fOutput->FindObject(Form("hEvPlaneReso2%s",centrbinname.Data()));
+    if(htmp1) htmp1->Fill(TMath::Cos(2.*deltaSub)); //RP resolution   
+    deltaSub =rpangleeventB-rpangleeventC;
     if(TMath::Abs(deltaSub)>TMath::Pi()/2.){// difference of subevents reaction plane angle cannot be bigger than phi/2
       if(deltaSub>0.) deltaSub-=TMath::Pi();
       else deltaSub +=TMath::Pi();
     } 
-    ((TH1F*)fOutput->FindObject(Form("hEvPlaneReso3%s",centrbinname.Data())))->Fill(TMath::Cos(2.*deltaSub)); //RP resolution   
+    TH1F* htmp2=(TH1F*)fOutput->FindObject(Form("hEvPlaneReso3%s",centrbinname.Data()));
+    if(htmp2) htmp2->Fill(TMath::Cos(2.*deltaSub)); //RP resolution   
   }
 
   if(fDebug>2)printf("Loop on D candidates\n");
@@ -550,7 +565,7 @@ void AliAnalysisTaskSEHFv2::UserExec(Option_t */*option*/)
     Bool_t isSelBit=kTRUE;
     if(fDecChannel==0) isSelBit=d->HasSelectionBit(AliRDHFCuts::kDplusCuts);
     if(fDecChannel==1) isSelBit=d->HasSelectionBit(AliRDHFCuts::kD0toKpiCuts);
-    if(fDecChannel==2) isSelBit=d->HasSelectionBit(AliRDHFCuts::kD0fromDstarCuts);
+    if(fDecChannel==2) isSelBit=kTRUE;
     if(!isSelBit)continue;
     Int_t ptbin=fRDCuts->PtBin(d->Pt());
     if(ptbin<0) {
@@ -559,7 +574,7 @@ void AliAnalysisTaskSEHFv2::UserExec(Option_t */*option*/)
     }
     Bool_t isFidAcc = fRDCuts->IsInFiducialAcceptance(d->Pt(),d->Y(absPdgMom));
     if(!isFidAcc)continue;    
-    Int_t isSelected= fRDCuts->IsSelected(d,AliRDHFCuts::kCandidate,aod);
+    Int_t isSelected= fRDCuts->IsSelected(d,AliRDHFCuts::kAll,aod);
     if(!isSelected)continue;
 
     fhEventsInfo->Fill(2); // candidate selected
@@ -579,9 +594,9 @@ void AliAnalysisTaskSEHFv2::UserExec(Option_t */*option*/)
     Float_t deltaphi=GetPhi0Pi(phi-eventplane);
     
     //fill the histograms with the appropriate method
-    if(fDecChannel==0)FillDplus(d,arrayMC,ptbin,deltaphi,invMass,isSelected,icentr);
-    else if(fDecChannel==1)FillD02p(d,arrayMC,ptbin,deltaphi,invMass,isSelected,icentr);
-    else if(fDecChannel==2)FillDstar(d,arrayMC,ptbin,deltaphi,invMass,isSelected,icentr);
+    if(fDecChannel==0)FillDplus(d,arrayMC,ptbin,deltaphi,invMass,isSelected,icentr,phi);
+    else if(fDecChannel==1)FillD02p(d,arrayMC,ptbin,deltaphi,invMass,isSelected,icentr,phi);
+    else if(fDecChannel==2)FillDstar(d,arrayMC,ptbin,deltaphi,invMass,isSelected,icentr,phi);
     
     delete [] invMass;
   }
@@ -631,7 +646,7 @@ void AliAnalysisTaskSEHFv2::CalculateInvMasses(AliAODRecoDecayHF* d,Float_t*& ma
 //NB: the implementation for each candidate is responsibility of the corresponding developer
 
 //******************************************************************************
-void AliAnalysisTaskSEHFv2::FillDplus(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t ptbin,Float_t deltaphi, const Float_t* masses,Int_t isSel,Int_t icentr){
+void AliAnalysisTaskSEHFv2::FillDplus(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t ptbin,Float_t deltaphi, const Float_t* masses,Int_t isSel,Int_t icentr, Double_t phiD){
   //D+ channel
   if(!isSel){
     if(fDebug>3)AliWarning("Candidate not selected\n");
@@ -641,10 +656,15 @@ void AliAnalysisTaskSEHFv2::FillDplus(AliAODRecoDecayHF* d,TClonesArray *arrayMC
     if(fDebug>3)AliWarning("Masses not calculated\n");
     return;
   }
- 
-  ((TH2F*)fOutput->FindObject(Form("hMc2phi_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
-  ((TH2F*)fOutput->FindObject(Form("hMphi_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(deltaphi,masses[0]);
-  ((TH2F*)fOutput->FindObject(Form("hMPtCandcentr%d_%d",icentr-5,icentr)))->Fill(d->Pt(),masses[0]);
+  Int_t icentrmin=icentr-fCentBinSizePerMil;
+  ((TH2F*)fOutput->FindObject(Form("hMc2deltaphi_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
+  ((TH2F*)fOutput->FindObject(Form("hMs2deltaphi_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(TMath::Sin(2*deltaphi),masses[0]);
+  ((TH2F*)fOutput->FindObject(Form("hMdeltaphi_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(deltaphi,masses[0]);
+  ((TH2F*)fOutput->FindObject(Form("hMPtCandcentr%d_%d",icentrmin,icentr)))->Fill(d->Pt(),masses[0]);
+  ((TH2F*)fOutput->FindObject(Form("hCos2phiMass_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(TMath::Cos(2*phiD),masses[0]);
+  ((TH2F*)fOutput->FindObject(Form("hSin2phiMass_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(TMath::Sin(2*phiD),masses[0]);
+
+
   Int_t pdgdaughters[3] = {211,321,211};
 
   if(fReadMC){
@@ -656,17 +676,17 @@ void AliAnalysisTaskSEHFv2::FillDplus(AliAODRecoDecayHF* d,TClonesArray *arrayMC
       lab = d->MatchToMC(411,arrayMC,3,pdgdaughters);
     }
     if(lab>=0){ //signal
-      ((TH2F*)fOutput->FindObject(Form("hMc2phiS_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
-      ((TH2F*)fOutput->FindObject(Form("hMphiS_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(deltaphi,masses[0]);
+      ((TH2F*)fOutput->FindObject(Form("hMc2deltaphiS_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
+      ((TH2F*)fOutput->FindObject(Form("hMdeltaphiS_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(deltaphi,masses[0]);
     } else{ //background
-      ((TH2F*)fOutput->FindObject(Form("hMc2phiB_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
-      ((TH2F*)fOutput->FindObject(Form("hMphiB_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(deltaphi,masses[0]);
+      ((TH2F*)fOutput->FindObject(Form("hMc2deltaphiB_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
+      ((TH2F*)fOutput->FindObject(Form("hMdeltaphiB_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(deltaphi,masses[0]);
     } 
   }   
 }
 
 //******************************************************************************
-void AliAnalysisTaskSEHFv2::FillD02p(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t ptbin,Float_t deltaphi, const Float_t* masses,Int_t isSel,Int_t icentr){
+void AliAnalysisTaskSEHFv2::FillD02p(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t ptbin,Float_t deltaphi, const Float_t* masses,Int_t isSel,Int_t icentr,Double_t phiD){
 
   //D0->Kpi channel
 
@@ -675,15 +695,22 @@ void AliAnalysisTaskSEHFv2::FillD02p(AliAODRecoDecayHF* d,TClonesArray *arrayMC,
     if(fDebug>3)AliWarning("Masses not calculated\n");
     return;
   }
+  Int_t icentrmin=icentr-fCentBinSizePerMil;
   if(isSel==1 || isSel==3) {
-    ((TH2F*)fOutput->FindObject(Form("hMc2phi_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
-    ((TH2F*)fOutput->FindObject(Form("hMphi_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(deltaphi,masses[0]);
-    ((TH2F*)fOutput->FindObject(Form("hMPtCandcentr%d_%d",icentr-5,icentr)))->Fill(d->Pt(),masses[0]);
+    ((TH2F*)fOutput->FindObject(Form("hMc2deltaphi_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
+    ((TH2F*)fOutput->FindObject(Form("hMs2deltaphi_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(TMath::Sin(2*deltaphi),masses[0]);
+    ((TH2F*)fOutput->FindObject(Form("hMdeltaphi_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(deltaphi,masses[0]);
+    ((TH2F*)fOutput->FindObject(Form("hMPtCandcentr%d_%d",icentrmin,icentr)))->Fill(d->Pt(),masses[0]);
+    ((TH2F*)fOutput->FindObject(Form("hCos2phiMass_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(TMath::Cos(2*phiD),masses[0]);
+    ((TH2F*)fOutput->FindObject(Form("hSin2phiMass_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(TMath::Sin(2*phiD),masses[0]);
   }
   if(isSel>=2) {
-    ((TH2F*)fOutput->FindObject(Form("hMc2phi_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[1]);
-    ((TH2F*)fOutput->FindObject(Form("hMphi_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(deltaphi,masses[1]);
-    ((TH2F*)fOutput->FindObject(Form("hMPtCandcentr%d_%d",icentr-5,icentr)))->Fill(d->Pt(),masses[1]);
+    ((TH2F*)fOutput->FindObject(Form("hMc2deltaphi_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[1]);
+    ((TH2F*)fOutput->FindObject(Form("hMs2deltaphi_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(TMath::Sin(2*deltaphi),masses[1]);
+    ((TH2F*)fOutput->FindObject(Form("hMdeltaphi_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(deltaphi,masses[1]);
+    ((TH2F*)fOutput->FindObject(Form("hMPtCandcentr%d_%d",icentrmin,icentr)))->Fill(d->Pt(),masses[1]);
+    ((TH2F*)fOutput->FindObject(Form("hCos2phiMass_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(TMath::Cos(2*phiD),masses[1]);
+    ((TH2F*)fOutput->FindObject(Form("hSin2phiMass_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(TMath::Sin(2*phiD),masses[1]);
   }
 
   //MC histograms
@@ -707,16 +734,16 @@ void AliAnalysisTaskSEHFv2::FillD02p(AliAODRecoDecayHF* d,TClonesArray *arrayMC,
        Int_t pdgMC = dMC->GetPdgCode();
        
        if(pdgMC==prongPdgPlus) {
-         ((TH2F*)fOutput->FindObject(Form("hMc2phiS_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
-         ((TH2F*)fOutput->FindObject(Form("hMphiS_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(deltaphi,masses[0]);
+         ((TH2F*)fOutput->FindObject(Form("hMc2deltaphiS_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
+         ((TH2F*)fOutput->FindObject(Form("hMdeltaphiS_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(deltaphi,masses[0]);
        }
        else {
-         ((TH2F*)fOutput->FindObject(Form("hMc2phiR_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
-         ((TH2F*)fOutput->FindObject(Form("hMphiR_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(deltaphi,masses[0]);
+         ((TH2F*)fOutput->FindObject(Form("hMc2deltaphiR_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
+         ((TH2F*)fOutput->FindObject(Form("hMdeltaphiR_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(deltaphi,masses[0]);
        }
       } else {
-       ((TH2F*)fOutput->FindObject(Form("hMc2phiB_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
-       ((TH2F*)fOutput->FindObject(Form("hMphi_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(deltaphi,masses[0]);
+       ((TH2F*)fOutput->FindObject(Form("hMc2deltaphiB_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
+       ((TH2F*)fOutput->FindObject(Form("hMdeltaphi_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(deltaphi,masses[0]);
       }
     }
     if(isSel>=2){ //D0bar
@@ -725,22 +752,22 @@ void AliAnalysisTaskSEHFv2::FillD02p(AliAODRecoDecayHF* d,TClonesArray *arrayMC,
        Int_t pdgMC = dMC->GetPdgCode();
        
        if(pdgMC==prongPdgMinus) {
-         ((TH2F*)fOutput->FindObject(Form("hMc2phiS_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[1]);
-         ((TH2F*)fOutput->FindObject(Form("hMphiS_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(deltaphi,masses[1]);
+         ((TH2F*)fOutput->FindObject(Form("hMc2deltaphiS_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[1]);
+         ((TH2F*)fOutput->FindObject(Form("hMdeltaphiS_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(deltaphi,masses[1]);
        }
        else {
-         ((TH2F*)fOutput->FindObject(Form("hMc2phiR_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[1]);
-         ((TH2F*)fOutput->FindObject(Form("hMphi_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(deltaphi,masses[1]);
+         ((TH2F*)fOutput->FindObject(Form("hMc2deltaphiR_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[1]);
+         ((TH2F*)fOutput->FindObject(Form("hMdeltaphi_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(deltaphi,masses[1]);
        }
       } else {
-       ((TH2F*)fOutput->FindObject(Form("hMc2phiB_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[1]);
-       ((TH2F*)fOutput->FindObject(Form("hMphiB_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(deltaphi,masses[1]);
+       ((TH2F*)fOutput->FindObject(Form("hMc2deltaphiB_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[1]);
+       ((TH2F*)fOutput->FindObject(Form("hMdeltaphiB_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(deltaphi,masses[1]);
       }
     }
   }
 }
 //******************************************************************************
-void AliAnalysisTaskSEHFv2::FillDstar(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t ptbin,Float_t deltaphi, const Float_t* masses,Int_t isSel,Int_t icentr){
+void AliAnalysisTaskSEHFv2::FillDstar(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t ptbin,Float_t deltaphi, const Float_t* masses,Int_t isSel,Int_t icentr, Double_t phiD){
   //D* channel
   if(!isSel){
     if(fDebug>3)AliWarning("Candidate not selected\n");
@@ -751,9 +778,15 @@ void AliAnalysisTaskSEHFv2::FillDstar(AliAODRecoDecayHF* d,TClonesArray *arrayMC
     return;
   }
 
-  ((TH2F*)fOutput->FindObject(Form("hMc2phi_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
-  ((TH2F*)fOutput->FindObject(Form("hMphi_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(deltaphi,masses[0]);
-  ((TH2F*)fOutput->FindObject(Form("hMPtCandcentr%d_%d",icentr-5,icentr)))->Fill(d->Pt(),masses[0]);
+  Float_t deltaphi1 = deltaphi;
+  if(deltaphi1 > TMath::PiOver2()) deltaphi1 = TMath::Pi() - deltaphi1; 
+  Int_t icentrmin=icentr-fCentBinSizePerMil;
+  ((TH2F*)fOutput->FindObject(Form("hMc2deltaphi_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
+  ((TH2F*)fOutput->FindObject(Form("hMs2deltaphi_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(TMath::Sin(2*deltaphi),masses[0]);
+  ((TH2F*)fOutput->FindObject(Form("hMdeltaphi_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(deltaphi1,masses[0]);
+  ((TH2F*)fOutput->FindObject(Form("hMPtCandcentr%d_%d",icentrmin,icentr)))->Fill(d->Pt(),masses[0]); 
+  ((TH2F*)fOutput->FindObject(Form("hCos2phiMass_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(TMath::Cos(2*phiD),masses[0]);
+  ((TH2F*)fOutput->FindObject(Form("hSin2phiMass_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(TMath::Sin(2*phiD),masses[0]);
   Int_t pdgDgDStartoD0pi[2]={421,211};
   Int_t pdgDgD0toKpi[2]={321,211};
   
@@ -761,11 +794,11 @@ void AliAnalysisTaskSEHFv2::FillDstar(AliAODRecoDecayHF* d,TClonesArray *arrayMC
     Int_t lab=-1;
     lab = ((AliAODRecoCascadeHF*)d)->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,arrayMC);
     if(lab>=0){ //signal
-      ((TH2F*)fOutput->FindObject(Form("hMc2phiS_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
-      ((TH2F*)fOutput->FindObject(Form("hMphiS_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(deltaphi,masses[0]);
+      ((TH2F*)fOutput->FindObject(Form("hMc2deltaphiS_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
+      ((TH2F*)fOutput->FindObject(Form("hMdeltaphiS_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(deltaphi,masses[0]);
     } else{ //background
-      ((TH2F*)fOutput->FindObject(Form("hMc2phiB_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
-      ((TH2F*)fOutput->FindObject(Form("hMphiB_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(deltaphi,masses[0]);
+      ((TH2F*)fOutput->FindObject(Form("hMc2deltaphiB_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
+      ((TH2F*)fOutput->FindObject(Form("hMdeltaphiB_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(deltaphi,masses[0]);
     } 
   }
 }
@@ -778,7 +811,18 @@ void AliAnalysisTaskSEHFv2::SetEventPlaneMethod(Int_t method){
   }
   fEventPlaneMeth=method;
 }
-
+//________________________________________________________________________
+void AliAnalysisTaskSEHFv2::SetNTPCSubEvents(Int_t nsub){
+  if(nsub>3||nsub<2){
+    AliWarning("Only 2 or 3 subevents implemented. Setting 2 subevents\n");
+    nsub=2;
+  }
+  if(fEventPlaneMeth==kTPC&&nsub==3){
+    AliWarning("V0 must be enabled to run 3 TPC subevents. Enabling...\n");
+    fEventPlaneMeth=kTPCVZERO;
+  }
+  fSubEvents=nsub;
+}
 //________________________________________________________________________
 Float_t AliAnalysisTaskSEHFv2::GetPhi0Pi(Float_t phi){
   // Sets the phi angle in the range 0-pi
@@ -803,14 +847,12 @@ Float_t AliAnalysisTaskSEHFv2::GetEventPlaneForCandidate(AliAODRecoDecayHF* d, c
     qx = pl->GetQContributionXArray();
     qy = pl->GetQContributionYArray();
     qcopy = *q;
-    }
-  else {
-    if(d->Eta()>0.){
+  }else {
+    if(d->Eta()<0.){
       qx = pl->GetQContributionXArraysub1();
       qy = pl->GetQContributionYArraysub1();
       qcopy = *qsub1;
-    }
-    else{
+    }else{
       qx = pl->GetQContributionXArraysub2();
       qy = pl->GetQContributionYArraysub2();
       qcopy = *qsub2;