Speedup of v2 task (Giacomo)

diff --git a/PWGHF/vertexingHF/charmFlow/AliAnalysisTaskSEHFv2.cxx b/PWGHF/vertexingHF/charmFlow/AliAnalysisTaskSEHFv2.cxx
index 93b81af68c34a7031eb905dbf63e0a9bfad5ae3c..3f6aa8d98bd0b4bf4e42af710d5e1e67431440c4 100644 (file)
--- a/PWGHF/vertexingHF/charmFlow/AliAnalysisTaskSEHFv2.cxx
+++ b/PWGHF/vertexingHF/charmFlow/AliAnalysisTaskSEHFv2.cxx
@@ -83,18 +83,16 @@ AliAnalysisTaskSE(),
   fhEventsInfo(0),
   fOutput(0),
   fRDCuts(0),
   fhEventsInfo(0),
   fOutput(0),
   fRDCuts(0),

-  fParHist(0),
-  fLowmasslimit(1.765),
-  fUpmasslimit(1.965),
+  fLowmasslimit(1.669),
+  fUpmasslimit(2.069),

   fNPtBins(1),
   fNPtBins(1),

-  fNPhiBinLims(2),
-  fPhiBins(0),

   fNMassBins(200),
   fReadMC(kFALSE),    
   fUseAfterBurner(kFALSE),
   fDecChannel(0),
   fAfterBurner(0),
   fNMassBins(200),
   fReadMC(kFALSE),    
   fUseAfterBurner(kFALSE),
   fDecChannel(0),
   fAfterBurner(0),

-  fUseV0EP(kFALSE),
+  fEventPlaneMeth(kTPCVZERO),
+  fEventPlanesComp(10),

   fV0EPorder(2),
   fMinCentr(20),
   fMaxCentr(80),
   fV0EPorder(2),
   fMinCentr(20),
   fMaxCentr(80),


@@ -104,23 +102,21 @@ AliAnalysisTaskSE(),
 }
 
 //________________________________________________________________________
 }
 
 //________________________________________________________________________

-AliAnalysisTaskSEHFv2::AliAnalysisTaskSEHFv2(const char *name,AliRDHFCuts *rdCuts,Int_t decaychannel, Int_t nlimsphibin, Float_t *phibinlimits):
+AliAnalysisTaskSEHFv2::AliAnalysisTaskSEHFv2(const char *name,AliRDHFCuts *rdCuts,Int_t decaychannel):

   AliAnalysisTaskSE(name),
   fhEventsInfo(0),
   fOutput(0),
   fRDCuts(rdCuts),
   AliAnalysisTaskSE(name),
   fhEventsInfo(0),
   fOutput(0),
   fRDCuts(rdCuts),

-  fParHist(0),

   fLowmasslimit(0),
   fUpmasslimit(0),
   fNPtBins(1),
   fLowmasslimit(0),
   fUpmasslimit(0),
   fNPtBins(1),

-  fNPhiBinLims(2),
-  fPhiBins(0),

   fNMassBins(200),
   fReadMC(kFALSE),
   fUseAfterBurner(kFALSE),
   fDecChannel(decaychannel),
   fAfterBurner(0),
   fNMassBins(200),
   fReadMC(kFALSE),
   fUseAfterBurner(kFALSE),
   fDecChannel(decaychannel),
   fAfterBurner(0),

-  fUseV0EP(kFALSE),
+  fEventPlaneMeth(kTPCVZERO),
+  fEventPlanesComp(10),

   fV0EPorder(2),
   fMinCentr(20),
   fMaxCentr(80),
   fV0EPorder(2),
   fMinCentr(20),
   fMaxCentr(80),


@@ -152,11 +148,6 @@ AliAnalysisTaskSEHFv2::AliAnalysisTaskSEHFv2(const char *name,AliRDHFCuts *rdCut
   if(pdg==413) SetMassLimits((Float_t)0.135,(Float_t)0.165);
   else SetMassLimits((Float_t)0.2,pdg); //check range
   fNPtBins=fRDCuts->GetNPtBins();
   if(pdg==413) SetMassLimits((Float_t)0.135,(Float_t)0.165);
   else SetMassLimits((Float_t)0.2,pdg); //check range
   fNPtBins=fRDCuts->GetNPtBins();

-  if(nlimsphibin>2) fNPhiBinLims=nlimsphibin;
-  else AliInfo("At least 2 limits in Delta phi needed");
-
-  fPhiBins=new Float_t[fNPhiBinLims];
-  for(Int_t i=0;i<fNPhiBinLims;i++) fPhiBins[i]=phibinlimits[i];

 
   if(fDebug>1)fRDCuts->PrintAll();
   // Output slot #1 writes into a TH1F container
 
   if(fDebug>1)fRDCuts->PrintAll();
   // Output slot #1 writes into a TH1F container


@@ -186,29 +177,9 @@ AliAnalysisTaskSEHFv2::~AliAnalysisTaskSEHFv2()
   delete fOutput;
   delete fhEventsInfo;
   delete fRDCuts;
   delete fOutput;
   delete fhEventsInfo;
   delete fRDCuts;

-  delete fParHist;
-
-  for(Int_t i=0;i<6;i++){
-    delete fHistvzero[i];
-  }
-

   delete fAfterBurner;
 }
 //_________________________________________________________________
   delete fAfterBurner;
 }
 //_________________________________________________________________

-void  AliAnalysisTaskSEHFv2::SetVZEROParHist(TH2D** histPar){
-  // Set the histograms for VZERO EP parameters
-  for(Int_t i=0;i<6;i++)fHistvzero[i]=(TH2D*)histPar[i]->Clone();
-  for(Int_t i=0;i<6;i++){
-    if(!fHistvzero[i]){
-      printf("No VZERO histograms!\n");
-      fUseV0EP=kFALSE;
-      return;
-    }
-  }
-  DefineOutput(4,TList::Class());
-  fUseV0EP=kTRUE;
-}
-//_________________________________________________________________

 void  AliAnalysisTaskSEHFv2::SetMassLimits(Float_t range, Int_t pdg){
   // Set limits for mass spectra plots
   Float_t mass=0;
 void  AliAnalysisTaskSEHFv2::SetMassLimits(Float_t range, Int_t pdg){
   // Set limits for mass spectra plots
   Float_t mass=0;


@@ -272,15 +243,13 @@ void AliAnalysisTaskSEHFv2::UserCreateOutputObjects()
  
   if(fDebug > 1) printf("AnalysisTaskSEHFv2::UserCreateOutputObjects() \n");
 
  
   if(fDebug > 1) printf("AnalysisTaskSEHFv2::UserCreateOutputObjects() \n");
 

-  fhEventsInfo=new TH1F(GetOutputSlot(1)->GetContainer()->GetName(), "Number of AODs scanned",8,-0.5,7.5);
+  fhEventsInfo=new TH1F(GetOutputSlot(1)->GetContainer()->GetName(), "Number of AODs scanned",6,-0.5,5.5);

   fhEventsInfo->GetXaxis()->SetBinLabel(1,"nEventsAnal");
   fhEventsInfo->GetXaxis()->SetBinLabel(2,"nEvSelected");
   fhEventsInfo->GetXaxis()->SetBinLabel(3,"nCandidatesSelected");
   fhEventsInfo->GetXaxis()->SetBinLabel(4,"out of pt bounds");
   fhEventsInfo->GetXaxis()->SetBinLabel(1,"nEventsAnal");
   fhEventsInfo->GetXaxis()->SetBinLabel(2,"nEvSelected");
   fhEventsInfo->GetXaxis()->SetBinLabel(3,"nCandidatesSelected");
   fhEventsInfo->GetXaxis()->SetBinLabel(4,"out of pt bounds");

-  fhEventsInfo->GetXaxis()->SetBinLabel(5,"Pile-up Rej");
-  fhEventsInfo->GetXaxis()->SetBinLabel(6,"N. of 0SMH");
-  fhEventsInfo->GetXaxis()->SetBinLabel(7,Form("Ev Sel in Centr %.0f-%.0f%s",fRDCuts->GetMinCentrality(),fRDCuts->GetMaxCentrality(),"%"));
-  fhEventsInfo->GetXaxis()->SetBinLabel(8,"mismatch lab");
+  fhEventsInfo->GetXaxis()->SetBinLabel(5,Form("Ev Sel in Centr %.0f-%.0f%s",fRDCuts->GetMinCentrality(),fRDCuts->GetMaxCentrality(),"%"));
+  fhEventsInfo->GetXaxis()->SetBinLabel(6,"mismatch lab");

   fhEventsInfo->GetXaxis()->SetNdivisions(1,kFALSE);
 
 
   fhEventsInfo->GetXaxis()->SetNdivisions(1,kFALSE);
 
 


@@ -291,53 +260,19 @@ void AliAnalysisTaskSEHFv2::UserCreateOutputObjects()
   
   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+5;icentr<=fMaxCentr;icentr=icentr+5){
     TString centrname;centrname.Form("centr%d_%d",icentr-5,icentr);

-    for(Int_t i=0;i<fNPtBins;i++){
-      
-      TString hname;
-      TString title;
-      hname.Form("hPhi_pt%d",i);hname.Append(centrname.Data());
-      title.Form("Phi distribution (Pt bin %d %s);#phi;Entries",i,centrname.Data());
-     
-      TH1F* hPhi=new TH1F(hname.Data(),title.Data(),96,0.,2*TMath::Pi());
-      hPhi->Sumw2();
-      fOutput->Add(hPhi);
-     
-      for(Int_t j=0;j<fNPhiBinLims-1;j++){
-       
-       hname.Form("hMass_pt%dphi%d",i,j);hname.Append(centrname.Data());
-       title.Form("Invariant mass (Pt bin %d, Phi bin %d, %s);M(GeV/c^{2});Entries",i,j,centrname.Data());
-       
-       TH1F* hMass=new TH1F(hname.Data(),title.Data(),fNMassBins,fLowmasslimit,fUpmasslimit);
-       hMass->Sumw2();
-       fOutput->Add(hMass);
-       
-       
-       if(fReadMC){
-         hname.Form("hSgn_pt%dphi%d",i,j);hname.Append(centrname.Data());
-         title.Form("Invariant mass S (Pt bin %d, Phi bin %d, %s);M(GeV/c^{2});Entries",i,j,centrname.Data());
-         TH1F* hSgn=new TH1F(hname.Data(),title.Data(),fNMassBins,fLowmasslimit,fUpmasslimit);
-         hSgn->Sumw2();
-         fOutput->Add(hSgn);
-        
-         hname.Form("hBkg_pt%dphi%d",i,j);hname.Append(centrname.Data());
-         title.Form("Invariant mass B (Pt bin %d, Phi bin %d, %s);M(GeV/c^{2});Entries",i,j,centrname.Data());
-         TH1F* hBkg=new TH1F(hname.Data(),title.Data(),fNMassBins,fLowmasslimit,fUpmasslimit);
-         hBkg->Sumw2();
-         fOutput->Add(hBkg);
-        
-         if((fDecChannel != AliAnalysisTaskSEHFv2::kDplustoKpipi) && (fDecChannel != AliAnalysisTaskSEHFv2::kDstartoKpipi)){
-           hname.Form("hRfl_pt%dphi%d",i,j);hname.Append(centrname.Data());
-           title.Form("Invariant mass Reflections (Pt bin %d, Phi bin %d %s);M(GeV/c^{2});Entries",i,j,centrname.Data());
-           TH1F* hRfl=new TH1F(hname.Data(),title.Data(),fNMassBins,fLowmasslimit,fUpmasslimit);
-           hRfl->Sumw2();
-           fOutput->Add(hRfl);
-         }
-       }
-      }

 
 

-      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);
+      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
+      

       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);
       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);


@@ -355,22 +290,18 @@ void AliAnalysisTaskSEHFv2::UserCreateOutputObjects()
        }
       }
     }
        }
       }
     }

-  
-    TH2F* hMphi=new TH2F(Form("hMphi%s",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);

 
 

-    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);
+
+    //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* hEvPlaneneg=new TH1F(Form("hEvPlaneneg%s",centrname.Data()),Form("Event plane angle %s;#phi Ev Plane;Entries",centrname.Data()),200,0.,TMath::Pi());
     fOutput->Add(hEvPlaneneg);
        
 
     TH1F* hEvPlaneneg=new TH1F(Form("hEvPlaneneg%s",centrname.Data()),Form("Event plane angle %s;#phi Ev Plane;Entries",centrname.Data()),200,0.,TMath::Pi());
     fOutput->Add(hEvPlaneneg);
        

-       TH1F* hEvPlanepos=new TH1F(Form("hEvPlanepos%s",centrname.Data()),Form("Event plane angle %s;#phi Ev Plane;Entries",centrname.Data()),200,0.,TMath::Pi());
+    TH1F* hEvPlanepos=new TH1F(Form("hEvPlanepos%s",centrname.Data()),Form("Event plane angle %s;#phi Ev Plane;Entries",centrname.Data()),200,0.,TMath::Pi());

     fOutput->Add(hEvPlanepos);
 
     fOutput->Add(hEvPlanepos);
 

-    //TH1F* hEvPlaneCheck=new TH1F(Form("hEvPlaneCheck%s",centrname.Data()),Form("Event plane angle - Event plane angle per candidate %s;(#phi(Ev Plane) - #phi(Ev Plane Candidate))/#phi(EvPlane);Entries",centrname.Data()),200,-0.2,0.2);
-    //fOutput->Add(hEvPlaneCheck);
-

     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);
 
     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);
 


@@ -378,20 +309,9 @@ void AliAnalysisTaskSEHFv2::UserCreateOutputObjects()
     fOutput->Add(hEvPlaneReso);
   }
   
     fOutput->Add(hEvPlaneReso);
   }
   

-  TH1F* hPhiBins=new TH1F("hPhiBins","Bins in #Delta#phi used in this analysis;#phi bin;n jobs",fNPhiBinLims-1,fPhiBins);
-  fOutput->Add(hPhiBins);
-  for(Int_t k=0;k<fNPhiBinLims-1;k++)hPhiBins->SetBinContent(k+1,1);

   PostData(1,fhEventsInfo);
   PostData(2,fOutput);
   PostData(1,fhEventsInfo);
   PostData(2,fOutput);

-  if(fUseV0EP){
-    fParHist = new TList();
-    fParHist->SetOwner();
-    fParHist->SetName("VZEROcorr");
-    for(Int_t i=0;i<6;i++){
-      fParHist->Add((TH2D*)fHistvzero[i]);
-    }
-    PostData(4,fParHist);
-  }
+

   return;
 }
 
   return;
 }
 


@@ -486,25 +406,30 @@ void AliAnalysisTaskSEHFv2::UserExec(Option_t */*option*/)
   Int_t nCand = arrayProng->GetEntriesFast();
   if(fDebug>2) printf("Number of D2H: %d\n",nCand);
 
   Int_t nCand = arrayProng->GetEntriesFast();
   if(fDebug>2) printf("Number of D2H: %d\n",nCand);
 

-  // trigger class for PbPb C0SMH-B-NOPF-ALLNOTRD, C0SMH-B-NOPF-ALL
-  TString trigclass=aod->GetFiredTriggerClasses();
-  if(trigclass.Contains("C0SMH-B-NOPF-ALLNOTRD") || trigclass.Contains("C0SMH-B-NOPF-ALL")) fhEventsInfo->Fill(5);
+  Bool_t isEvSel=fRDCuts->IsEventSelected(aod);
+  if(!isEvSel){
+    if(!fRDCuts->IsEventRejectedDueToCentrality())fhEventsInfo->Fill(4);
+    return;
+  }

 
 

-  if(fRDCuts->IsEventSelectedInCentrality(aod)>0) return;
-  else fhEventsInfo->Fill(6);
-  if(fRDCuts->IsEventSelected(aod)) fhEventsInfo->Fill(1);
-  else{
-    if(fRDCuts->GetWhyRejection()==1) // rejected for pileup
-      fhEventsInfo->Fill(4);
+  fhEventsInfo->Fill(1);
+  AliEventplane *pl=aod->GetEventplane();
+  if(!pl){
+    AliError("AliAnalysisTaskSEHFv2::UserExec:no eventplane! v2 analysis without eventplane not possible!\n");

     return;
   }
     return;
   }

- 
-  AliEventplane *pl=0x0;
-  TVector2* q=0x0;
-  Double_t rpangleevent=0;
+  
+  //Event plane
+  Double_t eventplane=0;
+  Double_t rpangleTPC=0;
+  Double_t rpangleVZERO=0;
+  Double_t planereso=0;
+  Double_t deltaPsi=0;

   Double_t rpangleeventneg=0;
   Double_t rpangleeventpos=0;
   Double_t rpangleeventneg=0;
   Double_t rpangleeventpos=0;

-  Double_t eventplane=0;
+
+  //For candidate removal from TPC EP
+  TVector2* q=0x0;

   TVector2 *qsub1=0x0;
   TVector2 *qsub2=0x0;
   
   TVector2 *qsub1=0x0;
   TVector2 *qsub2=0x0;
   


@@ -518,100 +443,107 @@ void AliAnalysisTaskSEHFv2::UserExec(Option_t */*option*/)
     }
   }
   TString centrbinname=Form("centr%d_%d",icentr-5,icentr);
     }
   }
   TString centrbinname=Form("centr%d_%d",icentr-5,icentr);

+

   if(fReadMC){
   if(fReadMC){

-    fUseV0EP=kTRUE;
+    fEventPlaneMeth=kVZERO;

     TRandom3 *g = new TRandom3(0);
     TRandom3 *g = new TRandom3(0);

-    rpangleevent=g->Rndm()*TMath::Pi();
+    rpangleVZERO=g->Rndm()*TMath::Pi();

     delete g;g=0x0;
     delete g;g=0x0;

-    eventplane=rpangleevent;
-       ((TH1F*)fOutput->FindObject(Form("hEvPlanepos%s",centrbinname.Data())))->Fill(rpangleevent);
-    if(fUseAfterBurner)fAfterBurner->SetEventPlane((Double_t)rpangleevent);
+    eventplane=rpangleVZERO;
+    if(fUseAfterBurner)fAfterBurner->SetEventPlane((Double_t)eventplane);

   }else{
   }else{

-    if(fUseV0EP){
-      rpangleevent=GetEventPlaneFromV0(aod);
-      eventplane=rpangleevent;
-         ((TH1F*)fOutput->FindObject(Form("hEvPlanepos%s",centrbinname.Data())))->Fill(rpangleevent);
-    }else{
-      // event plane and resolution 
-      //--------------------------------------------------------------------------
-      // extracting Q vectors and calculating v2 + resolution
-      pl = aod->GetHeader()->GetEventplaneP();
-      if(!pl){
-       AliError("AliAnalysisTaskSEHFv2::UserExec:no eventplane! v2 analysis without eventplane not possible!\n");
-       return;
+    if(fEventPlaneMeth!=kTPC){
+      //VZERO EP and resolution
+      rpangleVZERO=pl->GetEventplane("V0",aod,fV0EPorder);
+      if(fEventPlaneMeth!=kTPCVZERO){
+       //Using V0A/C to determine resolution
+       rpangleeventneg= pl->GetEventplane("V0A",aod,fV0EPorder);
+       rpangleeventpos= pl->GetEventplane("V0C",aod,fV0EPorder);
+       deltaPsi =rpangleeventneg-rpangleeventpos;
+       if(TMath::Abs(deltaPsi)>TMath::Pi()/2.){// difference of subevents reaction plane angle cannot be bigger than phi/2
+         if(deltaPsi>0.) deltaPsi-=TMath::Pi();
+         else deltaPsi +=TMath::Pi();
+       } 
+       eventplane=rpangleVZERO;

       }
       }

+    }
+    if(fEventPlaneMeth!=kVZERO){
+      // TPC event plane and resolution 
+
+      // extracting Q vectors and calculating v2 + resolution

       if(fEtaGap){
         qsub1 = pl->GetQsub1();
        qsub2 = pl->GetQsub2();
        rpangleeventpos = qsub1->Phi()/2.;
        rpangleeventneg = qsub2->Phi()/2.;
       if(fEtaGap){
         qsub1 = pl->GetQsub1();
        qsub2 = pl->GetQsub2();
        rpangleeventpos = qsub1->Phi()/2.;
        rpangleeventneg = qsub2->Phi()/2.;

-       ((TH1F*)fOutput->FindObject(Form("hEvPlanepos%s",centrbinname.Data())))->Fill(rpangleeventpos);
-       ((TH1F*)fOutput->FindObject(Form("hEvPlaneneg%s",centrbinname.Data())))->Fill(rpangleeventneg);
-       }
+      }

       else if(!fEtaGap){
        q = pl->GetQVector();
       else if(!fEtaGap){
        q = pl->GetQVector();

-        rpangleevent = pl->GetEventplane("Q");
-       ((TH1F*)fOutput->FindObject(Form("hEvPlanepos%s",centrbinname.Data())))->Fill(rpangleevent); // reaction plane angle without autocorrelations removal
-        } 
-       
-      Double_t deltaPsi = pl->GetQsubRes();
-      if(TMath::Abs(deltaPsi)>TMath::Pi()/2.){
-       if(deltaPsi>0.) deltaPsi-=TMath::Pi();
-       else deltaPsi +=TMath::Pi();
-      } // difference of subevents reaction plane angle cannot be bigger than phi/2
-      Double_t planereso = TMath::Cos(2.*deltaPsi); // reaction plane resolution
-      //--------------------------------------------------------------------------
-      ((TH1F*)fOutput->FindObject(Form("hEvPlaneReso%s",centrbinname.Data())))->Fill(planereso);    
+        rpangleTPC = pl->GetEventplane("Q");
+      } 
+      if(fEventPlaneMeth!=kVZEROTPC){
+       deltaPsi = pl->GetQsubRes();
+       if(TMath::Abs(deltaPsi)>TMath::Pi()/2.){
+         if(deltaPsi>0.) deltaPsi-=TMath::Pi();
+         else deltaPsi +=TMath::Pi();
+       } // difference of subevents reaction plane angle cannot be bigger than phi/2
+       planereso = TMath::Cos(2.*deltaPsi); // reaction plane resolution
+      }

     }
   }
     }
   }

-  
+  if(TMath::Abs(rpangleTPC-rpangleVZERO)>fEventPlanesComp)return;
+  //Filling EP-related histograms
+  planereso = TMath::Cos(2.*deltaPsi); // reaction plane resolution
+  ((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("hEvPlanepos%s",centrbinname.Data())))->Fill(rpangleeventpos); //Angle of first subevent
+  ((TH1F*)fOutput->FindObject(Form("hEvPlaneneg%s",centrbinname.Data())))->Fill(rpangleeventneg); //Angle of second subevent

 
 

+  //Loop on D candidates

   for (Int_t iCand = 0; iCand < nCand; iCand++) {
   for (Int_t iCand = 0; iCand < nCand; iCand++) {

-    

     d=(AliAODRecoDecayHF*)arrayProng->UncheckedAt(iCand);
     d=(AliAODRecoDecayHF*)arrayProng->UncheckedAt(iCand);

-    
-    Bool_t isFidAcc = fRDCuts->IsInFiducialAcceptance(d->Pt(),d->Y(absPdgMom));
-    Int_t isSelected= fRDCuts->IsSelected(d,AliRDHFCuts::kCandidate,aod);

     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);
     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(isSelected && isFidAcc && isSelBit) {
-      fhEventsInfo->Fill(2); // candidate selected
-      if(fDebug>3) printf("+++++++Is Selected\n");
+    if(!isSelBit)continue;
+    Int_t ptbin=fRDCuts->PtBin(d->Pt());
+    if(ptbin<0) {
+      fhEventsInfo->Fill(3);
+      continue;
+    }
+    Bool_t isFidAcc = fRDCuts->IsInFiducialAcceptance(d->Pt(),d->Y(absPdgMom));
+    if(!isFidAcc)continue;    
+    Int_t isSelected= fRDCuts->IsSelected(d,AliRDHFCuts::kCandidate,aod);
+    if(!isSelected)continue;
+
+    fhEventsInfo->Fill(2); // candidate selected
+    if(fDebug>3) printf("+++++++Is Selected\n");

       
       

-      Float_t* invMass=0x0;
-      Int_t nmasses;
-      CalculateInvMasses(d,invMass,nmasses);
-
-      Int_t ptbin=fRDCuts->PtBin(d->Pt());
-      if(ptbin==-1) {
-       fhEventsInfo->Fill(3);
-       continue;
-      }
+    Float_t* invMass=0x0;
+    Int_t nmasses;
+    CalculateInvMasses(d,invMass,nmasses);

 
 

-      if(!fUseV0EP) {
-       eventplane = GetEventPlaneForCandidate(d,q,pl,qsub1,qsub2); // remove autocorrelations
-       ((TH1F*)fOutput->FindObject(Form("hEvPlaneCand%s",centrbinname.Data())))->Fill(rpangleevent-eventplane);
-       //((TH1F*)fOutput->FindObject(Form("hEvPlaneCheck%s",centrbinname.Data())))->Fill((rpangleevent-eventplane)/100.*rpangleevent);
-      }
+    if(fEventPlaneMeth%2==0){
+      eventplane = GetEventPlaneForCandidate(d,q,pl,qsub1,qsub2); // remove autocorrelations
+      ((TH1F*)fOutput->FindObject(Form("hEvPlaneCand%s",centrbinname.Data())))->Fill(rpangleTPC-eventplane);
+    }

 
 

-      Float_t phi=d->Phi();
-      ((TH1F*)fOutput->FindObject(Form("hPhi_pt%d%s",ptbin,centrbinname.Data())))->Fill(phi);
-      
-      if(fReadMC&&fUseAfterBurner)phi=fAfterBurner->GetNewAngle(d,arrayMC);
-      Float_t deltaphi=GetPhi0Pi(phi-eventplane);
-
-      //fill the histograms with the appropriate method
-      if(fDecChannel==0)FillDplus(d,arrayMC,ptbin,deltaphi,invMass,isSelected,icentr);
-      if(fDecChannel==1)FillD02p(d,arrayMC,ptbin,deltaphi,invMass,isSelected,icentr);
-      if(fDecChannel==2)FillDstar(d,arrayMC,ptbin,deltaphi,invMass,isSelected,icentr);
-    }// end if selected
+    Double_t phi=d->Phi(); 
+    if(fReadMC&&fUseAfterBurner)phi=fAfterBurner->GetNewAngle(d,arrayMC);
+    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);
+    
+    delete invMass;

   }
   
   PostData(1,fhEventsInfo);
   PostData(2,fOutput);
   }
   
   PostData(1,fhEventsInfo);
   PostData(2,fOutput);

-  //PostData(4,flowEvent);
+

   return;
 }
 
   return;
 }
 


@@ -650,9 +582,10 @@ void AliAnalysisTaskSEHFv2::CalculateInvMasses(AliAODRecoDecayHF* d,Float_t*& ma
 
 //******************************************************************************
 
 
 //******************************************************************************
 

-//Methods to fill the istograms with MC information, one for each candidate
+//Methods to fill the histograms, one for each channel

 //NB: the implementation for each candidate is responsibility of the corresponding developer
 
 //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){
   //D+ channel
   if(!isSel){
 void AliAnalysisTaskSEHFv2::FillDplus(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t ptbin,Float_t deltaphi, const Float_t* masses,Int_t isSel,Int_t icentr){
   //D+ channel
   if(!isSel){


@@ -664,11 +597,8 @@ void AliAnalysisTaskSEHFv2::FillDplus(AliAODRecoDecayHF* d,TClonesArray *arrayMC
     return;
   }
  
     return;
   }
  

-  Int_t phibin=GetPhiBin(deltaphi);
-
-  ((TH1F*)fOutput->FindObject(Form("hMass_pt%dphi%dcentr%d_%d",ptbin,phibin,icentr-5,icentr)))->Fill(masses[0]);

   ((TH2F*)fOutput->FindObject(Form("hMc2phi_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
   ((TH2F*)fOutput->FindObject(Form("hMc2phi_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);

-  ((TH2F*)fOutput->FindObject(Form("hMphicentr%d_%d",icentr-5,icentr)))->Fill(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 pdgdaughters[3] = {211,321,211};
 
   ((TH2F*)fOutput->FindObject(Form("hMPtCandcentr%d_%d",icentr-5,icentr)))->Fill(d->Pt(),masses[0]);
   Int_t pdgdaughters[3] = {211,321,211};
 


@@ -681,18 +611,16 @@ void AliAnalysisTaskSEHFv2::FillDplus(AliAODRecoDecayHF* d,TClonesArray *arrayMC
       lab = d->MatchToMC(411,arrayMC,3,pdgdaughters);
     }
     if(lab>=0){ //signal
       lab = d->MatchToMC(411,arrayMC,3,pdgdaughters);
     }
     if(lab>=0){ //signal

-      ((TH1F*)fOutput->FindObject(Form("hSgn_pt%dphi%dcentr%d_%d",ptbin,phibin,icentr-5,icentr)))->Fill(masses[0]);

       ((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]);
     } else{ //background
       ((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]);
     } else{ //background

-      ((TH1F*)fOutput->FindObject(Form("hBkg_pt%dphi%dcentr%d_%d",ptbin,phibin,icentr-5,icentr)))->Fill(masses[0]);

       ((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("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]);
     } 
   }   
 }
 

-
+//******************************************************************************

 void AliAnalysisTaskSEHFv2::FillD02p(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t ptbin,Float_t deltaphi, const Float_t* masses,Int_t isSel,Int_t icentr){
 
   //D0->Kpi channel
 void AliAnalysisTaskSEHFv2::FillD02p(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t ptbin,Float_t deltaphi, const Float_t* masses,Int_t isSel,Int_t icentr){
 
   //D0->Kpi channel


@@ -702,21 +630,17 @@ void AliAnalysisTaskSEHFv2::FillD02p(AliAODRecoDecayHF* d,TClonesArray *arrayMC,
     if(fDebug>3)AliWarning("Masses not calculated\n");
     return;
   }
     if(fDebug>3)AliWarning("Masses not calculated\n");
     return;
   }

-  Int_t phibin=GetPhiBin(deltaphi);

   if(isSel==1 || isSel==3) {
   if(isSel==1 || isSel==3) {

-    ((TH1F*)fOutput->FindObject(Form("hMass_pt%dphi%dcentr%d_%d",ptbin,phibin,icentr-5,icentr)))->Fill(masses[0]);

     ((TH2F*)fOutput->FindObject(Form("hMc2phi_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
     ((TH2F*)fOutput->FindObject(Form("hMc2phi_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);

-    ((TH2F*)fOutput->FindObject(Form("hMphicentr%d_%d",icentr-5,icentr)))->Fill(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]);
   }
   if(isSel>=2) {
     ((TH2F*)fOutput->FindObject(Form("hMPtCandcentr%d_%d",icentr-5,icentr)))->Fill(d->Pt(),masses[0]);
   }
   if(isSel>=2) {

-    ((TH1F*)fOutput->FindObject(Form("hMass_pt%dphi%dcentr%d_%d",ptbin,phibin,icentr-5,icentr)))->Fill(masses[1]);

     ((TH2F*)fOutput->FindObject(Form("hMc2phi_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[1]);
     ((TH2F*)fOutput->FindObject(Form("hMc2phi_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[1]);

-    ((TH2F*)fOutput->FindObject(Form("hMphicentr%d_%d",icentr-5,icentr)))->Fill(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("hMPtCandcentr%d_%d",icentr-5,icentr)))->Fill(d->Pt(),masses[1]);
   }
 

-

   //MC histograms
   if(fReadMC){
 
   //MC histograms
   if(fReadMC){
 


@@ -738,19 +662,16 @@ void AliAnalysisTaskSEHFv2::FillD02p(AliAODRecoDecayHF* d,TClonesArray *arrayMC,
        Int_t pdgMC = dMC->GetPdgCode();
        
        if(pdgMC==prongPdgPlus) {
        Int_t pdgMC = dMC->GetPdgCode();
        
        if(pdgMC==prongPdgPlus) {

-         ((TH1F*)fOutput->FindObject(Form("hSgn_pt%dphi%dcentr%d_%d",ptbin,phibin,icentr-5,icentr)))->Fill(masses[0]);

          ((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]);
        }
        else {
          ((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]);
        }
        else {

-         ((TH1F*)fOutput->FindObject(Form("hRfl_pt%dphi%dcentr%d_%d",ptbin,phibin,icentr-5,icentr)))->Fill(masses[0]);

          ((TH2F*)fOutput->FindObject(Form("hMc2phiR_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
          ((TH2F*)fOutput->FindObject(Form("hMc2phiR_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);

-         ((TH2F*)fOutput->FindObject(Form("hMphiRcentr%d_%d",icentr-5,icentr)))->Fill(deltaphi,masses[0]);
+         ((TH2F*)fOutput->FindObject(Form("hMphiR_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(deltaphi,masses[0]);

        }
       } else {
        }
       } else {

-       ((TH1F*)fOutput->FindObject(Form("hBkg_pt%dphi%dcentr%d_%d",ptbin,phibin,icentr-5,icentr)))->Fill(masses[0]);

        ((TH2F*)fOutput->FindObject(Form("hMc2phiB_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
        ((TH2F*)fOutput->FindObject(Form("hMc2phiB_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);

-       ((TH2F*)fOutput->FindObject(Form("hMphicentr%d_%d",icentr-5,icentr)))->Fill(deltaphi,masses[0]);
+       ((TH2F*)fOutput->FindObject(Form("hMphi_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(deltaphi,masses[0]);

       }
     }
     if(isSel>=2){ //D0bar
       }
     }
     if(isSel>=2){ //D0bar


@@ -759,24 +680,21 @@ void AliAnalysisTaskSEHFv2::FillD02p(AliAODRecoDecayHF* d,TClonesArray *arrayMC,
        Int_t pdgMC = dMC->GetPdgCode();
        
        if(pdgMC==prongPdgMinus) {
        Int_t pdgMC = dMC->GetPdgCode();
        
        if(pdgMC==prongPdgMinus) {

-         ((TH1F*)fOutput->FindObject(Form("hSgn_pt%dphi%dcentr%d_%d",ptbin,phibin,icentr-5,icentr)))->Fill(masses[1]);

          ((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]);
        }
        else {
          ((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]);
        }
        else {

-         ((TH1F*)fOutput->FindObject(Form("hRfl_pt%dphi%dcentr%d_%d",ptbin,phibin,icentr-5,icentr)))->Fill(masses[1]);

          ((TH2F*)fOutput->FindObject(Form("hMc2phiR_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[1]);
          ((TH2F*)fOutput->FindObject(Form("hMc2phiR_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[1]);

-         ((TH2F*)fOutput->FindObject(Form("hMphicentr%d_%d",icentr-5,icentr)))->Fill(deltaphi,masses[1]);
+         ((TH2F*)fOutput->FindObject(Form("hMphi_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(deltaphi,masses[1]);

        }
       } else {
        }
       } else {

-       ((TH1F*)fOutput->FindObject(Form("hBkg_pt%dphi%dcentr%d_%d",ptbin,phibin,icentr-5,icentr)))->Fill(masses[1]);

        ((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("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]);
       }
     }
   }
 }

-
+//******************************************************************************

 void AliAnalysisTaskSEHFv2::FillDstar(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t ptbin,Float_t deltaphi, const Float_t* masses,Int_t isSel,Int_t icentr){
   //D* channel
   if(!isSel){
 void AliAnalysisTaskSEHFv2::FillDstar(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t ptbin,Float_t deltaphi, const Float_t* masses,Int_t isSel,Int_t icentr){
   //D* channel
   if(!isSel){


@@ -787,11 +705,9 @@ void AliAnalysisTaskSEHFv2::FillDstar(AliAODRecoDecayHF* d,TClonesArray *arrayMC
     if(fDebug>3)AliWarning("Masses not calculated\n");
     return;
   }
     if(fDebug>3)AliWarning("Masses not calculated\n");
     return;
   }

-  Int_t phibin=GetPhiBin(deltaphi);
-  
-  ((TH1F*)fOutput->FindObject(Form("hMass_pt%dphi%dcentr%d_%d",ptbin,phibin,icentr-5,icentr)))->Fill(masses[0]);
+

   ((TH2F*)fOutput->FindObject(Form("hMc2phi_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
   ((TH2F*)fOutput->FindObject(Form("hMc2phi_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);

-  ((TH2F*)fOutput->FindObject(Form("hMphicentr%d_%d",icentr-5,icentr)))->Fill(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 pdgDgDStartoD0pi[2]={421,211};
   Int_t pdgDgD0toKpi[2]={321,211};
   ((TH2F*)fOutput->FindObject(Form("hMPtCandcentr%d_%d",icentr-5,icentr)))->Fill(d->Pt(),masses[0]);
   Int_t pdgDgDStartoD0pi[2]={421,211};
   Int_t pdgDgD0toKpi[2]={321,211};


@@ -800,44 +716,24 @@ void AliAnalysisTaskSEHFv2::FillDstar(AliAODRecoDecayHF* d,TClonesArray *arrayMC
     Int_t lab=-1;
     lab = ((AliAODRecoCascadeHF*)d)->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,arrayMC);
     if(lab>=0){ //signal
     Int_t lab=-1;
     lab = ((AliAODRecoCascadeHF*)d)->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,arrayMC);
     if(lab>=0){ //signal

-      ((TH1F*)fOutput->FindObject(Form("hSgn_pt%dphi%dcentr%d_%d",ptbin,phibin,icentr-5,icentr)))->Fill(masses[0]);

       ((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]);
     } else{ //background
       ((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]);
     } else{ //background

-      ((TH1F*)fOutput->FindObject(Form("hBkg_pt%dphi%dcentr%d_%d",ptbin,phibin,icentr-5,icentr)))->Fill(masses[0]);

       ((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("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]);
     } 
   }

-  
-

 }
 
 }
 

-

 //________________________________________________________________________
 //________________________________________________________________________

-Int_t AliAnalysisTaskSEHFv2::GetPhiBin(Float_t deltaphi) const{
-
-  //give the bin corresponding to the value of deltaphi according to the binning requested in the constructor
-  Int_t phibin=0;
-  for(Int_t i=0;i<fNPhiBinLims-1;i++) {
-    if(deltaphi>=fPhiBins[i] && deltaphi<fPhiBins[i+1]) {
-      phibin=i;
-      break;
-    }
+void AliAnalysisTaskSEHFv2::SetEventPlaneMethod(Int_t method){
+  if(method>3||method<0){
+    AliWarning("No EP method associated to the selection, setting to TPC EP\n");
+    method=kTPCVZERO;

   }
   }

-  return phibin;
+  fEventPlaneMeth=method;

 }
 }

-//________________________________________________________________________
-// Float_t AliAnalysisTaskSEHFv2::GetPhi02Pi(Float_t phi){
-//   Float_t result=phi;
-//   while(result<0){
-//     result=result+2*TMath::Pi();
-//   }
-//   while(result>TMath::Pi()*2){
-//     result=result-2*TMath::Pi();
-//   }
-//   return result;
-// }
+

 //________________________________________________________________________
 Float_t AliAnalysisTaskSEHFv2::GetPhi0Pi(Float_t phi){
   // Sets the phi angle in the range 0-pi
 //________________________________________________________________________
 Float_t AliAnalysisTaskSEHFv2::GetPhi0Pi(Float_t phi){
   // Sets the phi angle in the range 0-pi


@@ -848,7 +744,7 @@ Float_t AliAnalysisTaskSEHFv2::GetPhi0Pi(Float_t phi){
   while(result>TMath::Pi()){
     result=result-TMath::Pi();
   }
   while(result>TMath::Pi()){
     result=result-TMath::Pi();
   }

-  return result;
+   return result;

 }
 
 //________________________________________________________________________
 }
 
 //________________________________________________________________________


@@ -876,8 +772,6 @@ Float_t AliAnalysisTaskSEHFv2::GetEventPlaneForCandidate(AliAODRecoDecayHF* d, c
     }
   }
   
     }
   }
   

-  
- 

   if(fDecChannel==2){
     //D* -- Yifei, Alessandro,Robert
     AliAODRecoDecayHF2Prong* theD0particle = ((AliAODRecoCascadeHF*)d)->Get2Prong();
   if(fDecChannel==2){
     //D* -- Yifei, Alessandro,Robert
     AliAODRecoDecayHF2Prong* theD0particle = ((AliAODRecoCascadeHF*)d)->Get2Prong();


@@ -945,204 +839,64 @@ Float_t AliAnalysisTaskSEHFv2::GetEventPlaneForCandidate(AliAODRecoDecayHF* d, c
   return qcopy.Phi()/2.;
 
 }
   return qcopy.Phi()/2.;
 
 }

-//________________________________________________________________________
-Float_t AliAnalysisTaskSEHFv2::GetEventPlaneFromV0(AliAODEvent *aodEvent){
-  // Compute event plane for VZERO
-
-  Int_t centr=fRDCuts->GetCentrality(aodEvent);
-  centr=centr-centr%10;
-  //temporary fix
-  if(centr<20)centr=20;
-  if(centr>70)centr=70;
-  //end temporary fix
-  Int_t binx=0;
-  Int_t iParHist=(centr-20)/10;
-
-  TString name;name.Form("parhist%d_%d",centr,centr+10);
-
-  if(fDebug>15)printf("EPfromV0 centr %d, iparhist %d (%p-%p)\n",centr,iParHist,fParHist->FindObject(name.Data()),fParHist->At(iParHist));
-
-  Int_t runnumber=aodEvent->GetRunNumber();
-  if(fParHist->At(iParHist)){
-    for(Int_t i=1;i<=((TH2D*)fParHist->At(iParHist))->GetNbinsX()&&binx<=0;i++){
-      Int_t run=atoi(((TH2D*)fParHist->At(iParHist))->GetXaxis()->GetBinLabel(i));
-      if(run>=runnumber)binx=i;
-    }
-  }else{
-    fhEventsInfo->Fill(7);
-  }
+// //________________________________________________________________________
+// Float_t AliAnalysisTaskSEHFv2::GetEventPlaneFromV0(AliAODEvent *aodEvent){
+//   // Compute event plane for VZERO - Obsolete, used for 2010 data
+
+//   Int_t centr=fRDCuts->GetCentrality(aodEvent);
+//   centr=centr-centr%10;
+//   //temporary fix
+//   if(centr<20)centr=20;
+//   if(centr>70)centr=70;
+//   //end temporary fix
+//   Int_t binx=0;
+//   Int_t iParHist=(centr-20)/10;
+
+//   TString name;name.Form("parhist%d_%d",centr,centr+10);
+
+//   if(fDebug>15)printf("EPfromV0 centr %d, iparhist %d (%p-%p)\n",centr,iParHist,fParHist->FindObject(name.Data()),fParHist->At(iParHist));
+
+//   Int_t runnumber=aodEvent->GetRunNumber();
+//   if(fParHist->At(iParHist)){
+//     for(Int_t i=1;i<=((TH2D*)fParHist->At(iParHist))->GetNbinsX()&&binx<=0;i++){
+//       Int_t run=atoi(((TH2D*)fParHist->At(iParHist))->GetXaxis()->GetBinLabel(i));
+//       if(run>=runnumber)binx=i;
+//     }
+//   }else{
+//     fhEventsInfo->Fill(7);
+//   }

 
 

-  AliFlowTrackCuts* cutsRP = AliFlowTrackCuts::GetStandardVZEROOnlyTrackCuts();
-  cutsRP->SetEvent(aodEvent, MCEvent());//, 0x0);
-  cutsRP->SetName( Form("rp_cuts") );
-  AliFlowTrackCuts* dummy = new AliFlowTrackCuts("null_cuts");
-  dummy->SetParamType(AliFlowTrackCuts::kGlobal);
-  dummy->SetPtRange(+1,-1); // select nothing QUICK
-  dummy->SetEtaRange(+1,-1); // select nothing VZERO
-  dummy->SetEvent(aodEvent,MCEvent());
-
-  //////////////// construct the flow event container ////////////
-  AliFlowEvent flowEvent(cutsRP,dummy);
-  flowEvent.SetReferenceMultiplicity( 64 );
-  for(Int_t i=0;i<64&&binx>0;i++){
-    AliFlowTrack *flowTrack=flowEvent.GetTrack(i);
-    Double_t inte=((TH2D*)fParHist->At(iParHist))->Integral(binx,binx,i+1,i+1);
-    if(inte>0)flowTrack->SetWeight(flowTrack->Weight()/inte);
-  }
-  if(fDebug>15)printf("EPfromV0 flow tracks weights done\n");
+//   AliFlowTrackCuts* cutsRP = AliFlowTrackCuts::GetStandardVZEROOnlyTrackCuts();
+//   cutsRP->SetEvent(aodEvent, MCEvent());//, 0x0);
+//   cutsRP->SetName( Form("rp_cuts") );
+//   AliFlowTrackCuts* dummy = new AliFlowTrackCuts("null_cuts");
+//   dummy->SetParamType(AliFlowTrackCuts::kGlobal);
+//   dummy->SetPtRange(+1,-1); // select nothing QUICK
+//   dummy->SetEtaRange(+1,-1); // select nothing VZERO
+//   dummy->SetEvent(aodEvent,MCEvent());
+
+//   //////////////// construct the flow event container ////////////
+//   AliFlowEvent flowEvent(cutsRP,dummy);
+//   flowEvent.SetReferenceMultiplicity( 64 );
+//   for(Int_t i=0;i<64&&binx>0;i++){
+//     AliFlowTrack *flowTrack=flowEvent.GetTrack(i);
+//     Double_t inte=((TH2D*)fParHist->At(iParHist))->Integral(binx,binx,i+1,i+1);
+//     if(inte>0)flowTrack->SetWeight(flowTrack->Weight()/inte);
+//   }
+//   if(fDebug>15)printf("EPfromV0 flow tracks weights done\n");

 
 

-  AliFlowVector qvec=flowEvent.GetQ(fV0EPorder);
-  Double_t angleEP=(1./(Double_t)fV0EPorder)*qvec.Phi();
-  if(fDebug>15)printf("EPfromV0 phi %f\n",angleEP);
-  return angleEP;
-}
+//   AliFlowVector qvec=flowEvent.GetQ(fV0EPorder);
+//   Double_t angleEP=(1./(Double_t)fV0EPorder)*qvec.Phi();
+//   if(fDebug>15)printf("EPfromV0 phi %f\n",angleEP);
+//   return angleEP;
+// }

 //________________________________________________________________________
 void AliAnalysisTaskSEHFv2::Terminate(Option_t */*option*/)
 {
 //________________________________________________________________________
 void AliAnalysisTaskSEHFv2::Terminate(Option_t */*option*/)
 {

-  

   // Terminate analysis
   //
   if(fDebug > 1) printf("AnalysisTaskSEHFv2: Terminate() \n");
   // Terminate analysis
   //
   if(fDebug > 1) printf("AnalysisTaskSEHFv2: Terminate() \n");

-  /*
-    fhEventsInfo = dynamic_cast<TH1F*> (GetOutputData(1));
-    if(!fhEventsInfo){
-    printf("Error: hEventsInfo not available\n");
-    return;
-    }
-
-    fOutput = dynamic_cast<TList*> (GetOutputData(2));
-    if (!fOutput) {     
-    printf("ERROR: fOutput not available\n");
-    return;
-    }
-    switch(fDecChannel){
-    case 0:
-    fRDCuts = dynamic_cast<AliRDHFCutsDplustoKpipi*> (GetOutputData(3));
-    break;
-    case 1:
-    fRDCuts = dynamic_cast<AliRDHFCutsD0toKpi*> (GetOutputData(3));
-    break;
-    case 2:
-    fRDCuts = dynamic_cast<AliRDHFCutsDStartoKpipi*> (GetOutputData(3));
-    break;
-    default:
-    break;
-    }
-    if (!fRDCuts) {     
-    printf("ERROR: fRDCuts not available\n");
-    return;
-    }
-   
-    // TCanvas* cvex=new TCanvas("example","Example Mass");
-    // cvex->cd();
-    // ((TH1F*)fOutput->FindObject("hMass_pt0phi0"))->Draw();
-    // TCanvas* cv2d=new TCanvas("2d","mass-cos2phi");
-    // cv2d->cd();
-    // ((TH2F*)fOutput->FindObject("hMc2phi"))->Draw("colz");
-    // TCanvas *cstat=new TCanvas("cstat","Stat");
-    // cstat->SetGridy();
-    // cstat->cd();
-    // fhEventsInfo->Draw("htext0");
-   
-    //  TMultiGraph *multig = new TMultiGraph();
-    if(fDecChannel==2)return;
-    TGraphErrors *g[fNPtBins];
-    TH1F *h = new TH1F("h","h",100,0.,1.);
-    TString hname;
-    TString gname;
-    for(Int_t ipt = 0;ipt<fNPtBins;ipt++){
-    g[ipt] = new TGraphErrors(fNPhiBinLims);
-    gname.Form("hMass_pt%d",ipt);
-    g[ipt]->SetTitle(gname.Data());
-    for(Int_t iphi = 0;iphi<fNPhiBinLims;iphi++){  
-    hname.Form("hMass_pt%dphi%d",ipt,iphi);
-    h = (TH1F*)fOutput->FindObject("hMass_pt0phi0");
-    AliHFMassFitter fitter(h,fLowmasslimit,fUpmasslimit,2,0);
-    Int_t pdg=0;
-    switch(fDecChannel){
-    case 0:
-    pdg=411;
-    break;
-    case 1:
-    pdg=421;
-    break;
-    case 2:
-    pdg=413;
-    break;
-    default:
-    break;
-    }
-    fitter.SetInitialGaussianMean(TDatabasePDG::Instance()->GetParticle(pdg)->Mass());
-    fitter.SetInitialGaussianSigma(0.012);
-    fitter.InitNtuParam(Form("ntuPtbin%d",ipt));
-    Double_t signal=0, errSignal=0;
-    if(fitter.MassFitter(kFALSE)){
-    fitter.Signal(3,signal,errSignal);
-    }
-    g[ipt]->SetPoint(iphi,fPhiBins[iphi],signal);
-    g[ipt]->SetPointError(iphi,fPhiBins[iphi],errSignal);
-    }//end loop over phi
-     //     multig->Add(g[ipt],"ap");
-     }//end loop on pt
-     TCanvas *cdndphi = new TCanvas("dN/d#phi","dN/d#phi");
-     cdndphi->Divide(1,fNPtBins); 
-     for(Int_t ipt = 0;ipt<fNPtBins;ipt++){
-     cdndphi->cd(ipt+1);
-     g[ipt]->Draw("AP");
-     }
-  */
+ 

   return;
 }
 
   return;
 }
 

-//-------------------------------------------
-/*
-  Float_t GetEventPlaneFromVZERO(){
-  AliAODHFUtil *info = (AliAODHFUtil*) aodEvent->GetList()->FindObject("fHFUtilInfo");
-  if (!info) return -999.;
-  //============= FIX ONLY FOR AOD033
-  Double_t *par0;
-  Double_t par0_137161[64] = { 6.71e-02 , 6.86e-02 , 7.06e-02 , 6.32e-02 , 
-  5.91e-02 , 6.07e-02 , 5.78e-02 , 5.73e-02 , 5.91e-02 , 6.22e-02 , 
-  5.90e-02 , 6.11e-02 , 5.55e-02 , 5.29e-02 , 5.19e-02 , 5.56e-02 , 
-  6.25e-02 , 7.03e-02 , 5.64e-02 , 5.81e-02 , 4.57e-02 , 5.30e-02 , 
-  5.13e-02 , 6.43e-02 , 6.27e-02 , 6.48e-02 , 6.07e-02 , 1.01e-01 , 
-  6.68e-02 , 7.16e-02 , 6.36e-02 , 5.95e-02 , 2.52e-02 , 2.82e-02 , 
-  2.56e-02 , 2.86e-02 , 2.82e-02 , 2.10e-02 , 2.13e-02 , 2.32e-02 , 
-  2.75e-02 , 4.34e-02 , 3.78e-02 , 4.52e-02 , 4.11e-02 , 3.89e-02 , 
-  4.10e-02 , 3.73e-02 , 4.51e-02 , 5.07e-02 , 5.42e-02 , 4.74e-02 , 
-  4.33e-02 , 4.44e-02 , 4.64e-02 , 3.01e-02 , 6.38e-02 , 5.26e-02 , 
-  4.99e-02 , 5.26e-02 , 5.47e-02 , 3.84e-02 , 5.00e-02 , 5.20e-02 };
-  Double_t par0_137366[64] = { 7.12e-02 , 7.34e-02 , 7.39e-02 , 6.54e-02 , 6.11e-02 , 6.31e-02 , 6.15e-02 , 
-  6.00e-02 , 6.10e-02 , 6.49e-02 , 6.17e-02 , 6.33e-02 , 6.00e-02 , 5.48e-02 , 
-  5.44e-02 , 5.81e-02 , 6.49e-02 , 7.07e-02 , 5.91e-02 , 6.18e-02 , 4.82e-02 , 
-  5.67e-02 , 5.36e-02 , 6.60e-02 , 6.37e-02 , 6.78e-02 , 6.31e-02 , 1.04e-01 , 
-  6.91e-02 , 7.32e-02 , 6.61e-02 , 6.16e-02 , 2.64e-02 , 2.81e-02 , 2.64e-02 , 
-  2.85e-02 , 2.87e-02 , 2.18e-02 , 2.19e-02 , 2.43e-02 , 2.81e-02 , 4.37e-02 , 
-  3.90e-02 , 4.66e-02 , 4.24e-02 , 4.09e-02 , 4.21e-02 , 3.88e-02 , 4.83e-02 , 
-  5.23e-02 , 5.44e-02 , 4.85e-02 , 4.42e-02 , 4.58e-02 , 4.74e-02 , 3.14e-02 , 
-  6.31e-02 , 5.30e-02 , 5.01e-02 , 5.33e-02 , 5.70e-02 , 3.95e-02 , 4.98e-02 , 5.31e-02 };
-  if (aodEvent->GetRunNumber() <= 137165)
-  par0=par0_137161;
-  else
-  par0=par0_137366;
-  Float_t vChCorr[64];
-  for(int i=0; i!=64; ++i)
-  vChCorr[i] = (info->GetVZEROChannel(i))/par0[i]/64.;
-  //============= END OF FIX AOD033
-  Float_t multR[8];
-  for(int i=0; i!=8; ++i) multR[i] = 0;
-  for(int i=0; i!=64; ++i)
-  multR[i/8] += vChCorr[i];
-  for(int i=0; i!=8; ++i) 
-  if(multR[i]) {
-  double Qx=0, Qy=0;
-  for(int j=0; j!=8; ++j) {
-  double phi = TMath::PiOver4()*(0.5+j);
-  Qx += TMath::Cos(2*phi)*vChCorr[i*8+j]/multR[i];
-  Qy += TMath::Sin(2*phi)*vChCorr[i*8+j]/multR[i];
-  }
-  }
-  return 0.5*TMath::ATan2(Qy,Qx)+TMath::PiOver2();
-  }
-
-*/

diff --git a/PWGHF/vertexingHF/charmFlow/AliAnalysisTaskSEHFv2.h b/PWGHF/vertexingHF/charmFlow/AliAnalysisTaskSEHFv2.h
index a1d5ff3ab028bc270f813102ac0135cd39356b59..7c7839a8d3c7b6c9b672495afb5825595e3391f1 100644 (file)
--- a/PWGHF/vertexingHF/charmFlow/AliAnalysisTaskSEHFv2.h
+++ b/PWGHF/vertexingHF/charmFlow/AliAnalysisTaskSEHFv2.h
@@ -34,9 +34,10 @@ class AliAnalysisTaskSEHFv2 : public AliAnalysisTaskSE
  public:
 
   enum DecChannel{kDplustoKpipi,kD0toKpi,kDstartoKpipi}; //more particles can be added
  public:
 
   enum DecChannel{kDplustoKpipi,kD0toKpi,kDstartoKpipi}; //more particles can be added

+  enum EventPlaneMeth{kTPC,kVZERO,kTPCVZERO,kVZEROTPC}; //Event plane to be calculated in the task

 
   AliAnalysisTaskSEHFv2();
 
   AliAnalysisTaskSEHFv2();

-  AliAnalysisTaskSEHFv2(const char *name, AliRDHFCuts *rdCuts, Int_t decaychannel,Int_t nbinsphi, Float_t *phibinlimits);
+  AliAnalysisTaskSEHFv2(const char *name, AliRDHFCuts *rdCuts, Int_t decaychannel);

  
   virtual ~AliAnalysisTaskSEHFv2();
 
  
   virtual ~AliAnalysisTaskSEHFv2();
 


@@ -44,22 +45,28 @@ class AliAnalysisTaskSEHFv2 : public AliAnalysisTaskSE
   void SetMassLimits(Float_t range,Int_t pdg);
   void SetMassLimits(Float_t lowlimit, Float_t uplimit);
   void SetNMassBins(Int_t nbins){fNMassBins=nbins;}
   void SetMassLimits(Float_t range,Int_t pdg);
   void SetMassLimits(Float_t lowlimit, Float_t uplimit);
   void SetNMassBins(Int_t nbins){fNMassBins=nbins;}

-  //void SetUseV0EP(Bool_t flagV0EP){fUseV0EP=flagV0EP;}

   void SetV0EventPlaneOrder(Int_t n){fV0EPorder=n;}
   void SetMinCentrality(Int_t mincentr){fMinCentr=mincentr;}
   void SetMaxCentrality(Int_t maxcentr){fMaxCentr=maxcentr;}
   void SetUseAfterBurner(Bool_t ab){fUseAfterBurner=ab;}
   void SetAfterBurner(AliHFAfterBurner *ab){fAfterBurner=ab;}
   void SetEtaGapFeatureForEventplaneFromTracks (Bool_t etaGap) {fEtaGap = etaGap;}
   void SetV0EventPlaneOrder(Int_t n){fV0EPorder=n;}
   void SetMinCentrality(Int_t mincentr){fMinCentr=mincentr;}
   void SetMaxCentrality(Int_t maxcentr){fMaxCentr=maxcentr;}
   void SetUseAfterBurner(Bool_t ab){fUseAfterBurner=ab;}
   void SetAfterBurner(AliHFAfterBurner *ab){fAfterBurner=ab;}
   void SetEtaGapFeatureForEventplaneFromTracks (Bool_t etaGap) {fEtaGap = etaGap;}

-  void SetVZEROParHist(TH2D** h);
-
+  void SetEventPlaneMethod(Int_t epmethod);
+  void SetTPCEPOnly(){SetEventPlaneMethod(kTPC);}
+  void SetVZEROEPOnly(){SetEventPlaneMethod(kVZERO);}
+  void SetTPCEP(){SetEventPlaneMethod(kTPCVZERO);}
+  void SetVZEROEP(){SetEventPlaneMethod(kVZEROTPC);}
+  void SetEventPlanesCompatibility(Float_t comp) {fEventPlanesComp=comp;}
+
+  Int_t GetEventPlaneMethod()const {return fEventPlaneMeth;}
+  Float_t GetEventPlanesCompatibility()const {return fEventPlanesComp;}

   Float_t GetUpperMassLimit()const {return fUpmasslimit;}
   Float_t GetLowerMassLimit()const {return fLowmasslimit;}
   Int_t GetNMassBins()const {return fNMassBins;}
   Float_t GetUpperMassLimit()const {return fUpmasslimit;}
   Float_t GetLowerMassLimit()const {return fLowmasslimit;}
   Int_t GetNMassBins()const {return fNMassBins;}

-  Int_t GetPhiBin(Float_t deltaphi) const;

   //Float_t GetPhi02Pi(Float_t phi);
   Float_t GetPhi0Pi(Float_t phi);
   AliHFAfterBurner *GetAfterBurner()const {return fAfterBurner;}
   //Float_t GetPhi02Pi(Float_t phi);
   Float_t GetPhi0Pi(Float_t phi);
   AliHFAfterBurner *GetAfterBurner()const {return fAfterBurner;}

+

   // Implementation of interface methods
   virtual void UserCreateOutputObjects();
   virtual void LocalInit();// {Init();}
   // Implementation of interface methods
   virtual void UserCreateOutputObjects();
   virtual void LocalInit();// {Init();}


@@ -77,25 +84,22 @@ class AliAnalysisTaskSEHFv2 : public AliAnalysisTaskSE
   void FillD02p(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t ptbin, Float_t dphi, const Float_t* masses, Int_t isSel,Int_t icentr);
   void FillDstar(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t ptbin, Float_t dphi, const Float_t* masses,Int_t isSel,Int_t icentr);
   Float_t GetEventPlaneForCandidate(AliAODRecoDecayHF* d, const TVector2* q,AliEventplane *pl,const TVector2* qsub1,const TVector2* qsub2);
   void FillD02p(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t ptbin, Float_t dphi, const Float_t* masses, Int_t isSel,Int_t icentr);
   void FillDstar(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t ptbin, Float_t dphi, const Float_t* masses,Int_t isSel,Int_t icentr);
   Float_t GetEventPlaneForCandidate(AliAODRecoDecayHF* d, const TVector2* q,AliEventplane *pl,const TVector2* qsub1,const TVector2* qsub2);

-  Float_t GetEventPlaneFromV0(AliAODEvent *aodEvent);
+  //  Float_t GetEventPlaneFromV0(AliAODEvent *aodEvent);

 
 
   TH1F* fhEventsInfo;           //! histogram send on output slot 1
   TList   *fOutput;             //! list send on output slot 2
   AliRDHFCuts *fRDCuts;         //cut values (saved in slot 3)
 
 
   TH1F* fhEventsInfo;           //! histogram send on output slot 1
   TList   *fOutput;             //! list send on output slot 2
   AliRDHFCuts *fRDCuts;         //cut values (saved in slot 3)

-  TList *fParHist;               //list for VZERO EP parameters (slot 4)
-  TH2D *fHistvzero[6];            //histograms for VZERO EP parameters

   Float_t fLowmasslimit;        //lower inv mass limit for histos
   Float_t fUpmasslimit;         //upper inv mass limit for histos
   Int_t fNPtBins;               //number of pt bins
   Float_t fLowmasslimit;        //lower inv mass limit for histos
   Float_t fUpmasslimit;         //upper inv mass limit for histos
   Int_t fNPtBins;               //number of pt bins

-  Int_t fNPhiBinLims;           //number of delta phi bins limits (= number of bins +1)
-  Float_t *fPhiBins;            //[fNPhiBinLims] limits of each phi bin

   Int_t fNMassBins;             //number of bins in the mass histograms
   Bool_t fReadMC;               //flag for access to MC
   Bool_t fUseAfterBurner;      //enable afterburning
   Int_t fDecChannel;            //decay channel identifier
   AliHFAfterBurner *fAfterBurner;//Afterburner options
   Int_t fNMassBins;             //number of bins in the mass histograms
   Bool_t fReadMC;               //flag for access to MC
   Bool_t fUseAfterBurner;      //enable afterburning
   Int_t fDecChannel;            //decay channel identifier
   AliHFAfterBurner *fAfterBurner;//Afterburner options

-  Bool_t fUseV0EP;              //flag to select EP method
+  Int_t fEventPlaneMeth;         //flag to select EP method
+  Float_t fEventPlanesComp;     // Maximum distance between TPC/VZERO event planes

   Int_t  fV0EPorder;            //harmonic for VZERO event plane
   Int_t fMinCentr;              //minimum centrality
   Int_t fMaxCentr;              //maximum centrality
   Int_t  fV0EPorder;            //harmonic for VZERO event plane
   Int_t fMinCentr;              //minimum centrality
   Int_t fMaxCentr;              //maximum centrality