Added classes and tasks for D meson v2 analysis. Just for reference, not inclueded...
authorprino <prino@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 16 Aug 2011 23:09:57 +0000 (23:09 +0000)
committerprino <prino@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 16 Aug 2011 23:09:57 +0000 (23:09 +0000)
PWG3/vertexingHF/charmFlow/AddTaskFlowD2H.C [new file with mode: 0644]
PWG3/vertexingHF/charmFlow/AddTaskHFv2.C [new file with mode: 0644]
PWG3/vertexingHF/charmFlow/AliAnalysisTaskFlowD2H.cxx [new file with mode: 0644]
PWG3/vertexingHF/charmFlow/AliAnalysisTaskFlowD2H.h [new file with mode: 0644]
PWG3/vertexingHF/charmFlow/AliAnalysisTaskSEHFv2.cxx [new file with mode: 0644]
PWG3/vertexingHF/charmFlow/AliAnalysisTaskSEHFv2.h [new file with mode: 0644]
PWG3/vertexingHF/charmFlow/AliEventPlaneResolution.cxx [new file with mode: 0644]
PWG3/vertexingHF/charmFlow/AliEventPlaneResolution.h [new file with mode: 0644]
PWG3/vertexingHF/charmFlow/DmesonsFlowAnalysis.C [new file with mode: 0644]
PWG3/vertexingHF/charmFlow/Extractv2from2Dhistos.C [new file with mode: 0644]

diff --git a/PWG3/vertexingHF/charmFlow/AddTaskFlowD2H.C b/PWG3/vertexingHF/charmFlow/AddTaskFlowD2H.C
new file mode 100644 (file)
index 0000000..b498906
--- /dev/null
@@ -0,0 +1,167 @@
+void AddTaskFlowD2H( int nDmeson, TString fileName, int RPsource,
+                     int EtaGapSP=3, double aMax=-0.2, double bMin=+0.2,
+                     int minCent=30, int maxCent=50, TString suffixName="C" )
+{
+  suffixName += Form("%d%d_%s",minCent,maxCent,fileName.Data());
+  //-R-P---c-u-t-s--------------------------------------------------------------
+  AliFlowTrackCuts* cutsRP;
+  double aMin, aMax, bMin, bMax;    // For SP method
+  switch (RPsource) {
+    case (0):
+      cutsRP = (AliFlowTrackCuts*) 
+               AliFlowTrackCuts::GetStandardTPCStandaloneTrackCuts2010();
+      cutsRP->SetEtaRange( -0.8, +0.8 );
+      cutsRP->SetPtRange( 0.2, 5.0 );
+      cutsRP->SetAODfilterBit(1);
+      cutsRP->SetName( Form("rp_cuts_%s",suffixName.Data()) );
+      aMin=-0.8, bMax=0.8;
+      break;
+    case (1):
+      cutsRP = (AliFlowTrackCuts*) 
+               AliFlowTrackCuts::GetStandardVZEROOnlyTrackCuts();
+      aMin = -3.7; aMax = -2.2 ; bMin = 3.4; bMax = +5.1;
+      break;
+  }
+  fileName.Append(".root");
+  //-P-O-I---c-u-t-s------------------------------------------------------------
+  double flowBands[5][2];
+  switch (nDmeson) { //thanks to ROOT for allowing this
+    case (AliRDHFCuts::kD0toKpiCuts):
+      AliRDHFCutsD0toKpi *cutsPOI = new AliRDHFCutsD0toKpi( Form("D0toKpi_%s",suffixName.Data()) );
+      flowBands[0][0]=1.700; flowBands[0][1]=1.750;
+      flowBands[1][0]=1.750; flowBands[1][1]=1.800;
+      flowBands[2][0]=1.830; flowBands[2][1]=1.900;
+      flowBands[3][0]=1.930; flowBands[3][1]=1.970;
+      flowBands[4][0]=1.970; flowBands[4][1]=2.030;
+      break;
+    case (AliRDHFCuts::kDstarCuts):
+      AliRDHFCutsDStartoKpipi *cutsPOI = new AliRDHFCutsDStartoKpipi( Form("DStartoKpipi_%s",suffixName.Data()) );
+      flowBands[0][0]=0.140; flowBands[0][1]=0.142;
+      flowBands[1][0]=0.144; flowBands[1][1]=0.147;
+      flowBands[2][0]=0.150; flowBands[2][1]=0.155;
+      flowBands[3][0]=0.155; flowBands[3][1]=0.160;
+      flowBands[4][0]=0.160; flowBands[4][1]=0.165;
+      break;
+  }
+  cutsPOI->SetStandardCutsPbPb2010();
+  cutsPOI->SetUseCentrality(AliRDHFCuts::kCentV0M);
+  cutsPOI->SetMinCentrality(minCent); cutsPOI->SetMaxCentrality(maxCent);
+  cutsPOI->SetUseAOD049();
+
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+  AliAnalysisDataContainer *cinput1 = mgr->GetCommonInputContainer();
+
+  //-----------------FLOWD2H TASK----------------------------
+  AliAnalysisTaskFlowD2H *taskSel = new AliAnalysisTaskFlowD2H(
+                                    Form("FlowD2H_%s",suffixName.Data()),
+                                    cutsRP, cutsPOI, nDmeson);
+  taskSel->SelectCollisionCandidates(AliVEvent::kMB);
+  taskSel->SetDebug();
+  taskSel->SetFlowEtaRangeAB( aMin, aMax, bMin, bMax );
+  TString Qvector;
+  switch (EtaGapSP) {
+    case 0:
+      taskSel->SetPOIEtaRange( -1.5, aMax );
+      Qvector = "Qb";
+      break;
+    case 1:
+      taskSel->SetPOIEtaRange( bMin, +1.5 );
+      Qvector = "Qa";
+      break;
+    case 3:
+      taskSel->SetPOIEtaRange( -1.5, 1.5 );
+      Qvector = "QaQb"; // dummy
+      break;
+  }
+  taskSel->SetFlowPtRange( 2, 16 );
+  for(int i=0; i!=5; ++i)
+    taskSel->SetFlowBandRange( i, flowBands[i][0], flowBands[i][1] );
+  mgr->AddTask(taskSel);
+  AliAnalysisDataContainer *coutputCandidatesQA = mgr->CreateContainer(
+     Form("QA_%s",suffixName.Data()),TList::Class(),
+     AliAnalysisManager::kOutputContainer, fileName+":QACandidates");
+  mgr->ConnectOutput(taskSel,1,coutputCandidatesQA);
+  mgr->ConnectInput (taskSel,0,cinput1);
+  AliAnalysisDataContainer *coutputCandidates[5];
+  for(int r=0; r!=5; ++r) {
+    coutputCandidates[r] = mgr->CreateContainer(
+       Form("Flow_MassBand%d_%s",r,suffixName.Data()),
+       AliFlowEventSimple::Class(), AliAnalysisManager::kExchangeContainer);
+    mgr->ConnectOutput(taskSel,2+r,coutputCandidates[r]);
+  }
+  // Scalar Product
+  AliAnalysisTaskScalarProduct *taskSP[5];
+  AliAnalysisDataContainer *coutputSP[5];
+  for(int r=0; r!=5; ++r) {
+    taskSP[r] = new AliAnalysisTaskScalarProduct(
+             Form("SP_MassBand%d_%s",r,suffixName.Data()),kFALSE);
+    taskSP[r]->SetRelDiffMsub(1.0);
+    taskSP[r]->SetApplyCorrectionForNUA(kTRUE);
+    taskSP[r]->SetTotalQvector( Qvector.Data() );
+    mgr->AddTask( taskSP[r] );
+    TString outputSP = fileName + Form(":outputSP_MassBand%d",r);
+    coutputSP[r] = mgr->CreateContainer(
+         Form("cobjSP_MassBand%d_%s",r,suffixName.Data()),TList::Class(),
+         AliAnalysisManager::kOutputContainer,outputSP.Data());
+    mgr->ConnectInput(taskSP[r],0,coutputCandidates[r]);
+    mgr->ConnectOutput(taskSP[r],1,coutputSP[r]);
+  }
+  if(RPsource!=0)
+    return;
+  if(EtaGapSP!=3)
+    return;
+  cout << "HERE!\n";
+  // Event Plane
+  AliAnalysisTaskEventPlane *taskEP[5];
+  AliAnalysisDataContainer *coutputEP[5];
+  for(int r=0; r!=5; ++r) {
+    taskEP[r] = new AliAnalysisTaskEventPlane(
+             Form("EP_MassBand%d_%s",r,suffixName.Data()),kFALSE);
+    taskEP[r]->SetApplyCorrectionForNUA(kTRUE);
+    taskEP[r]->SetNormalizationType( 0 );
+    mgr->AddTask( taskEP[r] );
+    TString outputEP = fileName + Form(":outputEP_MassBand%d",r);
+    coutputEP[r] = mgr->CreateContainer(
+         Form("cobjEP_MassBand%d_%s",r,suffixName.Data()),TList::Class(),
+         AliAnalysisManager::kOutputContainer,outputEP.Data());
+    mgr->ConnectInput(taskEP[r],0,coutputCandidates[r]);
+    mgr->ConnectOutput(taskEP[r],1,coutputEP[r]);
+  }
+  // Event Plane 2
+  AliAnalysisTaskEventPlane *taskEP2[5];
+  AliAnalysisDataContainer *coutputEP2[5];
+  for(int r=0; r!=5; ++r) {
+    taskEP2[r] = new AliAnalysisTaskEventPlane(
+             Form("EP2_MassBand%d_%s",r,suffixName.Data()),kFALSE);
+    taskEP2[r]->SetApplyCorrectionForNUA(kTRUE);
+    taskEP2[r]->SetNormalizationType( 0 );
+    mgr->AddTask( taskEP2[r] );
+    TString outputEP2 = fileName + Form(":outputEP2_MassBand%d",r);
+    coutputEP2[r] = mgr->CreateContainer(
+         Form("cobjEP2_MassBand%d_%s",r,suffixName.Data()),TList::Class(),
+         AliAnalysisManager::kOutputContainer,outputEP2.Data());
+    mgr->ConnectInput(taskEP2[r],0,coutputCandidates[r]);
+    mgr->ConnectOutput(taskEP2[r],1,coutputEP2[r]);
+  }
+  // Q-Cumulants
+  AliAnalysisTaskQCumulants *taskQC[5];
+  AliAnalysisDataContainer *coutputQC[5];
+  for(int r=0; r!=5; ++r) {
+    taskQC[r] = new AliAnalysisTaskQCumulants(
+             Form("QC_MassBand%d_%s",r,suffixName.Data()),kFALSE);
+    taskQC[r]->SetCalculateCumulantsVsM(kFALSE);
+    taskQC[r]->SetnBinsMult(10000);
+    taskQC[r]->SetMinMult(0.);
+    taskQC[r]->SetMaxMult(10000.);
+    taskQC[r]->SetApplyCorrectionForNUA(kTRUE);
+    taskQC[r]->SetFillMultipleControlHistograms(kFALSE);
+    mgr->AddTask( taskQC[r] );
+    TString outputQC = fileName + Form(":outputQC_MassBand%d",r);
+    coutputQC[r] = mgr->CreateContainer(
+         Form("cobjQC_MassBand%d_%s",r,suffixName.Data()),TList::Class(),
+         AliAnalysisManager::kOutputContainer,outputQC.Data());
+    mgr->ConnectInput(taskQC[r],0,coutputCandidates[r]);
+    mgr->ConnectOutput(taskQC[r],1,coutputQC[r]);
+  }
+}
+
diff --git a/PWG3/vertexingHF/charmFlow/AddTaskHFv2.C b/PWG3/vertexingHF/charmFlow/AddTaskHFv2.C
new file mode 100644 (file)
index 0000000..215701b
--- /dev/null
@@ -0,0 +1,121 @@
+AliAnalysisTaskSEHFv2 *AddTaskHFv2(TString filename="DplustoKpipiCuts.root", AliAnalysisTaskSEHFv2::DecChannel decCh=AliAnalysisTaskSEHFv2::kD0toKpi,Bool_t readMC=kFALSE,TString name="")
+{
+  //
+  // Test macro for the AliAnalysisTaskSE for  D 
+  // mesons v2 analysis with event plane method
+  // Authors: Chiara Bianchin, cbianchi@pd.infn.it, 
+  //          Robert Grajcarek, grajcarek@physi.uni-heidelberg.de
+  //          Giacomo Ortona, ortona@to.infn.it,
+  //          Carlos Perez Lara, carlos.eugenio.perez.lara@cern.ch
+  //          Francesco Prino, prino@to.infn.it
+  // Get the pointer to the existing analysis manager via the static access method.
+  //============================================================================
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+  if (!mgr) {
+    ::Error("AddTaskHFv2", "No analysis manager to connect to.");
+    return NULL;
+  }
+  Bool_t stdcuts=kFALSE;
+  TFile* filecuts=new TFile(filename.Data());
+  if(!filecuts->IsOpen()){
+    cout<<"Input file not found:  using std cut object"<<endl;
+    stdcuts=kTRUE;
+  }
+  
+  AliRDHFCuts *analysiscuts=0x0;
+  TString suffix="";
+
+  TString cutsobjname="loosercuts";
+  //Analysis cuts
+  switch (decCh){
+  case 0:
+    cutsobjname="AnalysisCuts"; 
+    if(stdcuts){
+      analysiscuts = new AliRDHFCutsDplustoKpipi();
+      analysiscuts->SetStandardCutsPbPb2010();
+    }
+    else analysiscuts = (AliRDHFCutsDplustoKpipi*)filecuts->Get(cutsobjname);
+    suffix="Dplus";
+    break;
+  case 1:
+    cutsobjname="D0toKpiCuts";
+    if(stdcuts) {
+      analysiscuts = new AliRDHFCutsD0toKpi();
+      analysiscuts->SetStandardCutsPbPb2010();
+    }
+    else analysiscuts = (AliRDHFCutsD0toKpi*)filecuts->Get(cutsobjname);
+    suffix="D0";
+    break;
+  case 2:
+    cutsobjname="DStartoKpipiCuts";
+    if(stdcuts) {
+      analysiscuts = new AliRDHFCutsDStartoKpipi();
+      analysiscuts->SetStandardCutsPbPb2010();
+    }
+    else analysiscuts = (AliRDHFCutsDStartoKpipi*)filecuts->Get(cutsobjname);
+    suffix="Dstar";
+    break;
+  default:
+    cout<<"Not available"<<endl;
+    break;
+  }
+
+  if(!analysiscuts){
+    cout<<"Specific AliRDHFCuts not found"<<endl;
+    return;
+  }
+
+  suffix+=name;
+  const Int_t nphibins=4;
+  const Int_t nphibinlimits=nphibins+1;
+  Float_t pi=TMath::Pi();
+  Float_t philimits[nphibinlimits]={0., pi/4.,pi/2., 3./4.*pi, pi};
+
+  //histogram for V0
+  TFile *fpar = TFile::Open("VZEROParHist.root");
+  TH2D *hh[6];
+  for(Int_t i=0;i<6;i++){
+    TString name;name.Form("parhist%d_%d",(i+2)*10,(i+3)*10);
+    hh[i]=(TH2D*)fpar->Get(name.Data());
+  }
+
+  // Analysis task    
+  AliAnalysisTaskSEHFv2 *v2Task = new AliAnalysisTaskSEHFv2("HFv2Analysis",analysiscuts,decCh,nphibinlimits,philimits,hh);
+  v2Task->SetReadMC(readMC);
+
+  v2Task->SetDebugLevel(0);
+  
+  mgr->AddTask(v2Task);
+
+  // Create containers for input/output
+
+  TString contname=Form("cinputv2%s",suffix.Data());
+  AliAnalysisDataContainer *cinputv2 = mgr->CreateContainer(contname.Data(),TChain::Class(),AliAnalysisManager::kInputContainer);
+
+  TString outputfile = AliAnalysisManager::GetCommonFileName();
+  TString outputhistos = outputfile += ":PWG3_D2H_HFv2"; 
+
+  contname=Form("hEventsInfo%s",suffix.Data());
+  AliAnalysisDataContainer *coutputstat = mgr->CreateContainer(contname.Data(),TH1F::Class(),AliAnalysisManager::kOutputContainer,outputfile.Data());
+
+  contname=Form("coutputv2%s",suffix.Data());
+  AliAnalysisDataContainer *coutputv2 = mgr->CreateContainer(contname.Data(),TList::Class(),AliAnalysisManager::kOutputContainer,outputfile.Data());
+
+  contname=Form("cutobj%s",suffix.Data());
+  AliAnalysisDataContainer *cutobj = mgr->CreateContainer(contname.Data(),AliRDHFCuts::Class(),AliAnalysisManager::kOutputContainer,outputfile.Data());
+  
+  contname=Form("coutputVZEROpar%s",suffix.Data());
+  AliAnalysisDataContainer *coutputpar = mgr->CreateContainer(contname.Data(),TH2D::Class(),AliAnalysisManager::kOutputContainer,outputfile.Data());
+
+  mgr->ConnectInput(v2Task,0,mgr->GetCommonInputContainer());
+  
+  mgr->ConnectOutput(v2Task,1,coutputstat);
+  
+  mgr->ConnectOutput(v2Task,2,coutputv2);
+
+  mgr->ConnectOutput(v2Task,3,cutobj);
+  mgr->ConnectOutput(v2Task,4,coutputpar);
+
+  return v2Task;
+}
diff --git a/PWG3/vertexingHF/charmFlow/AliAnalysisTaskFlowD2H.cxx b/PWG3/vertexingHF/charmFlow/AliAnalysisTaskFlowD2H.cxx
new file mode 100644 (file)
index 0000000..d65121f
--- /dev/null
@@ -0,0 +1,420 @@
+/**************************************************************************
+* Copyright(c) 1998-2008,ALICE Experiment at CERN,All rights reserved.    *
+*                                                                         *
+* Author: The ALICE Off-line Project.                                     *
+* Contributors are mentioned in the code where appropriate.               *
+*                                                                         *
+* Permission to use,copy,modify and distribute this software and its      *
+* documentation strictly for non-commercial purposes is hereby granted    *
+* without fee,provided that the above copyright notice appears in all     *
+* copies and that both the copyright notice and this permission notice    *
+* appear in the supporting documentation. The authors make no claims      *
+* about the suitability of this software for any purpose. It is           *
+* provided "as is" without express or implied warranty.                   *
+**************************************************************************/
+
+//==============================================================================
+// FlowD2H main task:
+// >> Make flowEvent with RPcuts and POIcuts given in constructor and passes it
+//    to the daughter tasks.
+// >> The POIcuts are polymorphic based on the AliRDHFCuts class allowing the 
+//    use of all charmed candidates reconstructed in the central barrel.
+// Author: Carlos Perez (cperez@cern.ch)
+//==============================================================================
+
+#include "TChain.h"
+#include "TList.h"
+#include "TH2D.h"
+#include "TProfile.h"
+
+#include "AliAnalysisTaskSE.h"
+#include "AliAnalysisManager.h"
+
+#include "AliCentrality.h"
+
+#include "AliAODEvent.h"
+#include "AliAODInputHandler.h"
+#include "AliAODTrack.h"
+
+#include "TMath.h"
+#include "TObjArray.h"
+
+#include "AliFlowCandidateTrack.h"
+#include "AliFlowTrackCuts.h"
+#include "AliFlowEvent.h"
+#include "AliFlowCommonConstants.h"
+
+#include "AliRDHFCuts.h"
+#include "AliRDHFCutsD0toKpi.h"
+#include "AliRDHFCutsDStartoKpipi.h"
+#include "AliRDHFCutsDplustoKpipi.h"
+
+#include "AliAODRecoDecayHF2Prong.h"
+#include "AliAODRecoDecayHF3Prong.h"
+#include "AliAODRecoCascadeHF.h"
+
+#include "AliAnalysisTaskFlowD2H.h"
+
+ClassImp(AliAnalysisTaskFlowD2H)
+
+//_____________________________________________________________________________
+AliAnalysisTaskFlowD2H::AliAnalysisTaskFlowD2H() :
+  AliAnalysisTaskSE(), fCutsRP(NULL), fCutsPOI(NULL), fSource(0),
+  fDebug(kFALSE), fHList(NULL), fAnaCuts(NULL)
+{
+// Default constructor
+  for(int i=0; i!=2; ++i)
+    fEvent[i] = NULL;
+  for(int i=0; i!=2; ++i)
+    fMass[i] = NULL;
+  for(int i=0; i!=2; ++i)
+    fPOIEta[i] = 0;
+  for(int i=0; i!=4; ++i)
+    fFlowEta[i] = 0;
+  for(int x=0; x!=2; ++x)
+    fFlowPts[x] = 0;
+  for(int x=0; x!=2; ++x)
+    for(int m=0; m!=5; ++m)
+      fFlowBands[x][m] = 0;
+}
+
+//_____________________________________________________________________________
+AliAnalysisTaskFlowD2H::AliAnalysisTaskFlowD2H(const char *name,
+    AliFlowTrackCuts *cutsRPs, AliRDHFCuts *cutsPOIs, Int_t specie) :
+  AliAnalysisTaskSE(name), fCutsRP(cutsRPs), 
+  fCutsPOI(cutsPOIs), fSource(specie),
+  fDebug(kFALSE), fHList(NULL), fAnaCuts(NULL)
+{
+// Standard constructor
+  for(int i=0; i!=2; ++i)
+    fEvent[i] = NULL;
+  for(int i=0; i!=2; ++i)
+    fMass[i] = NULL;
+  for(int i=0; i!=2; ++i)
+    fPOIEta[i] = 0;
+  for(int i=0; i!=4; ++i)
+    fFlowEta[i] = 0;
+  for(int x=0; x!=2; ++x)
+    fFlowPts[x] = 0;
+  for(int x=0; x!=2; ++x)
+    for(int m=0; m!=5; ++m)
+      fFlowBands[x][m] = 0;
+
+  DefineInput( 0,TChain::Class());
+  DefineOutput(1,TList::Class());
+  DefineOutput(2,AliFlowEventSimple::Class()); // first band
+  DefineOutput(3,AliFlowEventSimple::Class()); // second band
+  DefineOutput(4,AliFlowEventSimple::Class()); // third band
+  DefineOutput(5,AliFlowEventSimple::Class()); // fourth band
+  DefineOutput(6,AliFlowEventSimple::Class()); // fifth band
+}
+
+//_____________________________________________________________________________
+AliAnalysisTaskFlowD2H::~AliAnalysisTaskFlowD2H()
+{
+  if(fHList)
+    delete fHList;
+}
+
+//_____________________________________________________________________________
+void AliAnalysisTaskFlowD2H::UserCreateOutputObjects()
+{
+  if(fDebug)
+    printf("DEBUGGER: Creating output\n");
+  fHList = new TList();
+  fHList->SetOwner();
+  AddHistograms();
+  PostData(1,fHList);
+
+  AliFlowCommonConstants* cc = AliFlowCommonConstants::GetMaster();
+  cc->SetNbinsMult(10000);
+  cc->SetMultMin(0);
+  cc->SetMultMax(10000);
+
+  cc->SetNbinsPt(fFlowPts[1]-fFlowPts[0]);
+  cc->SetPtMin(fFlowPts[0]);
+  cc->SetPtMax(fFlowPts[1]);
+
+  cc->SetNbinsPhi(180);
+  cc->SetPhiMin(0.0);
+  cc->SetPhiMax(TMath::TwoPi());
+
+  cc->SetNbinsEta(200);
+  cc->SetEtaMin(-5.0);
+  cc->SetEtaMax(+5.0);
+
+  cc->SetNbinsQ(500);
+  cc->SetQMin(0.0);
+  cc->SetQMax(3.0);
+}
+//_____________________________________________________________________________
+void AliAnalysisTaskFlowD2H::AddHistograms()
+{
+  TList *tEvents = new TList();
+    tEvents->SetName("Events");
+    tEvents->SetOwner();
+    for(int i=0; i!=2; ++i) {
+      fEvent[i] = new TH2D(Form("Event%d",i),"Events;V0M;Arb",10,0,100,10,0,12000);
+      tEvents->Add(fEvent[i]);
+    }
+    fAnaCuts=new TProfile("Cuts","Analysis Cuts", 10,0,10);
+      fAnaCuts->Fill(0.5,fPOIEta[0],1);fAnaCuts->GetXaxis()->SetBinLabel(1,"ETAm");
+      fAnaCuts->Fill(1.5,fPOIEta[1],1);fAnaCuts->GetXaxis()->SetBinLabel(2,"ETAM");
+    tEvents->Add(fAnaCuts);
+    fHList->Add(tEvents);
+  TList *tCandidates;
+    tCandidates = new TList();
+    tCandidates->SetOwner();
+    tCandidates->SetName(Form("Candidates%d",fSource));
+    Double_t dBinMin=0, dBinMax=3;
+    Int_t nBins=3;
+    switch (fSource) {
+      case (AliRDHFCuts::kD0toKpiCuts): // 360/72 (bw=5MeV)
+        dBinMin = 1.695; dBinMax = 2.055; nBins=72; break;
+      case (AliRDHFCuts::kDstarCuts): // 36/72 (bw=0.5MeV)
+        dBinMin = 0.137; dBinMax = 0.173; nBins=72; break;
+      case (AliRDHFCuts::kDplusCuts): // 360/72 (bw=5MeV)
+        dBinMin = 1.695; dBinMax = 2.055; nBins=72; break;
+    }
+    for(int i=0; i!=2; ++i) {
+      fMass[i] = new TH2D( Form("Mass%d",i),
+                           Form("Mass%d;Mass [GeV];Pt [GeV]",i),
+                           nBins, dBinMin, dBinMax,
+                           fFlowPts[1]-fFlowPts[0], fFlowPts[0], fFlowPts[1] );
+      tCandidates->Add(fMass[i]);
+    }
+    fHList->Add(tCandidates);
+  if (fDebug) printf("DEBUGGER: Created histos for DMeson %d\n", fSource );
+}
+//_____________________________________________________________________________
+void AliAnalysisTaskFlowD2H::NotifyRun()
+{
+}
+//_____________________________________________________________________________
+void AliAnalysisTaskFlowD2H::UserExec(Option_t *)
+{
+  AliAODEvent *fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
+
+  if(!fAOD)
+    return;
+  Int_t iMulti = fAOD->GetNTracks();  /// TO DO
+  fEvent[0]->Fill( fCutsPOI->GetCentrality(fAOD), iMulti );
+  if(!fCutsPOI->IsEventSelected(fAOD)) return;
+  if(fCutsPOI->IsEventSelectedInCentrality(fAOD)>0) return;
+  fEvent[1]->Fill( fCutsPOI->GetCentrality(fAOD), iMulti );
+
+  if (fDebug) printf("Event selected\n");
+  fCutsRP->SetEvent(fAOD,MCEvent());
+  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(fAOD,MCEvent());
+  AliFlowEvent *flowEvent[5];
+  for(int r=0; r!=5; ++r) {
+    flowEvent[r] = new AliFlowEvent(fCutsRP,dummy);
+    flowEvent[r]->SetReferenceMultiplicity( iMulti );
+    flowEvent[r]->DefineDeadZone(0,0,0,0);
+    flowEvent[r]->TagSubeventsInEta( fFlowEta[0], fFlowEta[1],
+                                     fFlowEta[2], fFlowEta[3] );
+  }
+  delete dummy;
+  if (fDebug) printf(" ᶫFlowEvent has %d RPs\n", flowEvent[0]->NumberOfTracks() );
+
+  switch (fSource) {
+    case (AliRDHFCuts::kD0toKpiCuts):
+      FillD0toKpi(fAOD,flowEvent); break;
+    case (AliRDHFCuts::kDstarCuts):
+      FillDStartoKpipi(fAOD,flowEvent); break;
+    case (AliRDHFCuts::kDplusCuts):
+      FillDplustoKpipi(fAOD,flowEvent); break;
+  }
+
+  for(int m=0; m!=5; ++m)
+    PostData(2+m,flowEvent[m]);
+  PostData(1,fHList);
+
+}
+//______________________________________________________________________________
+void AliAnalysisTaskFlowD2H::FillD0toKpi(AliAODEvent *theAOD, 
+                                         AliFlowEvent *theMB[5] )
+{
+  TList *listHF = (TList*) theAOD->GetList();
+  if(!listHF) return;
+  TClonesArray *listDzero = (TClonesArray*) listHF->FindObject("D0toKpi");
+  if(!listDzero) return;
+  int nEntries = listDzero->GetEntriesFast();
+  if( !nEntries ) return;
+  AliRDHFCutsD0toKpi *fCutsD0toKpi = (AliRDHFCutsD0toKpi*) fCutsPOI;
+  if (fDebug) printf("  ᶫ%d candidates found. Looping...\n", nEntries);
+  const Int_t ndght=2;
+  Int_t nIDs[ndght];
+  for( int iEntry=0; iEntry!=nEntries; ++iEntry ) {
+    AliAODRecoDecayHF2Prong *D0 = (AliAODRecoDecayHF2Prong*) listDzero->UncheckedAt( iEntry );
+    if( !D0 ) continue;
+    if( !D0->HasSelectionBit(AliRDHFCuts::kD0toKpiCuts) ) continue;
+    if( !fCutsD0toKpi->IsInFiducialAcceptance(D0->Pt(),D0->Y(421)) )continue;
+    int topCut  = fCutsD0toKpi->IsSelected( D0, AliRDHFCuts::kAll, theAOD );
+    int nLevel=topCut>0?1:0;
+    if( (D0->Eta()<fPOIEta[0])||(D0->Eta()>fPOIEta[1]) )
+      nLevel=0;
+    Double_t MassD0=topCut>1?D0->InvMassD0bar():D0->InvMassD0();
+    for(int h=0; h!=nLevel+1; ++h)
+      fMass[h]->Fill(MassD0,D0->Pt());
+    if( (D0->Pt()<fFlowPts[0]) || (D0->Pt()>fFlowPts[1]) ) continue;
+    AliAODTrack* iT;
+    for(Int_t i=0; i!=ndght; ++i) {
+      iT = (AliAODTrack*)D0->GetDaughter(i);
+      nIDs[i] = iT->GetID();
+    }
+    // Candidates Insertion (done in filling method: faster)
+    if(nLevel)
+      for(Int_t r=0; r!=5; ++r)
+        if( (MassD0>=fFlowBands[0][r]) && (MassD0<fFlowBands[1][r]) ) {
+          AliFlowCandidateTrack *sTrack = (AliFlowCandidateTrack*)
+                                 MakeTrack(MassD0, D0->Pt(), D0->Phi(),
+                                           D0->Eta(), ndght, nIDs);
+          if(fDebug) printf("   ᶫInjecting D0 candidate on band %d \n", r);
+          for(Int_t iDau=0; iDau!=ndght; ++iDau)
+            for(Int_t iRPs=0; iRPs!=theMB[r]->NumberOfTracks(); ++iRPs) {
+              AliFlowTrack *iRP = (AliFlowTrack*) (theMB[r]->GetTrack(iRPs));
+              if(!iRP->InRPSelection()) continue;
+              if( fabs(sTrack->GetIDDaughter(iDau)) == fabs(iRP->GetID()) ) {
+                sTrack->SetDaughter(iDau,iRP);
+                iRP->SetForRPSelection(kFALSE);
+                if(fDebug) printf("    ᶫdaughter%d with fID %d was removed from this RP set\n", iDau, sTrack->GetIDDaughter(iDau));
+              }
+            }
+          theMB[r]->AddTrack(sTrack);
+        }
+    // END of Candidates Insertion
+  }
+}
+//______________________________________________________________________________
+void AliAnalysisTaskFlowD2H::FillDStartoKpipi(AliAODEvent *theAOD, 
+                                              AliFlowEvent *theMB[5] )
+{
+  TList *listHF = (TList*) theAOD->GetList();
+  if(!listHF) return;
+  TClonesArray *listDstar = (TClonesArray*) listHF->FindObject("Dstar");
+  if(!listDstar) return;
+  int nEntries = listDstar->GetEntriesFast();
+  if( !nEntries ) return;
+  AliRDHFCutsDStartoKpipi *fCutsDStartoKpipi = (AliRDHFCutsDStartoKpipi*) fCutsPOI;
+  if (fDebug) printf("  ᶫ%d candidates found. Looping...\n", nEntries);
+  const Int_t ndght=3;
+  Int_t nIDs[ndght];
+  for( int iEntry=0; iEntry!=nEntries; ++iEntry ) {
+    AliAODRecoCascadeHF *DS = (AliAODRecoCascadeHF*) listDstar->UncheckedAt( iEntry );
+    if( !DS ) continue;
+    if( !DS->HasSelectionBit(AliRDHFCuts::kDstarCuts) ) continue;
+    if( !fCutsDStartoKpipi->IsInFiducialAcceptance(DS->Pt(),DS->YDstar()) )continue;
+    int topCut = fCutsDStartoKpipi->IsSelected( DS, AliRDHFCuts::kAll );
+    int nLevel=topCut>0?1:0;
+    if( (DS->Eta()<fPOIEta[0])||(DS->Eta()>fPOIEta[1]) )
+      nLevel=0;
+    Double_t MassDS=DS->DeltaInvMass();
+    for(int h=0; h!=nLevel+1; ++h)
+      fMass[h]->Fill(MassDS,DS->Pt());
+    if( (DS->Pt()<fFlowPts[0]) || (DS->Pt()>fFlowPts[1]) ) continue;
+    AliAODRecoDecayHF2Prong *D0 = (AliAODRecoDecayHF2Prong*)DS->Get2Prong();
+    nIDs[0] = ((AliAODTrack*)D0->GetDaughter(0))->GetID();
+    nIDs[1] = ((AliAODTrack*)D0->GetDaughter(1))->GetID();
+    nIDs[2] = ((AliAODTrack*)DS->GetBachelor() )->GetID();
+    // Candidates Insertion (done in filling method: faster)
+    if(nLevel)
+      for(Int_t r=0; r!=5; ++r)
+        if( (MassDS>=fFlowBands[0][r]) && (MassDS<fFlowBands[1][r]) ) {
+          AliFlowCandidateTrack *sTrack = (AliFlowCandidateTrack*)
+                                 MakeTrack(MassDS, DS->Pt(), DS->Phi(),
+                                           DS->Eta(), ndght, nIDs);
+          if(fDebug) printf("   ᶫInjecting DStar candidate on band %d \n", r);
+          for(Int_t iDau=0; iDau!=ndght; ++iDau)
+            for(Int_t iRPs=0; iRPs!=theMB[r]->NumberOfTracks(); ++iRPs) {
+              AliFlowTrack *iRP = (AliFlowTrack*) (theMB[r]->GetTrack(iRPs));
+              if(!iRP->InRPSelection()) continue;
+              if( fabs(sTrack->GetIDDaughter(iDau)) == fabs(iRP->GetID()) ) {
+                sTrack->SetDaughter(iDau,iRP);
+                iRP->SetForRPSelection(kFALSE);
+                if(fDebug) printf("    ᶫdaughter%d with fID %d was removed from this RP set\n", iDau, sTrack->GetIDDaughter(iDau));
+              }
+            }
+          theMB[r]->AddTrack(sTrack);
+        }
+    // END of Candidates Insertion
+  }
+}
+//______________________________________________________________________________
+void AliAnalysisTaskFlowD2H::FillDplustoKpipi(AliAODEvent *theAOD, 
+                                              AliFlowEvent *theMB[5] )
+{
+  TList *listHF = (TList*) theAOD->GetList();
+  if(!listHF) return;
+  TClonesArray *listDplus = (TClonesArray*) listHF->FindObject("Charm3Prong");
+  if(!listDplus) return;
+  int nEntries = listDplus->GetEntriesFast();
+  if( !nEntries ) return;
+  AliRDHFCutsDplustoKpipi *fCutsDStartoKpipi = (AliRDHFCutsDplustoKpipi*) fCutsPOI;
+  if (fDebug) printf("  ᶫ%d candidates found. Looping...\n", nEntries);
+  const Int_t ndght=3;
+  Int_t nIDs[ndght];
+  for( int iEntry=0; iEntry!=nEntries; ++iEntry ) {
+    AliAODRecoDecayHF3Prong *Dp = (AliAODRecoDecayHF3Prong*) listDplus->UncheckedAt( iEntry );
+    if( !Dp ) continue;
+    if( !Dp->HasSelectionBit(AliRDHFCuts::kDplusCuts) ) continue;
+    if( !fCutsDStartoKpipi->IsInFiducialAcceptance(Dp->Pt(),Dp->YDplus()) )continue;
+    int topCut = fCutsDStartoKpipi->IsSelected( Dp, AliRDHFCuts::kAll );
+    int nLevel=topCut>0?1:0;
+    if( (Dp->Eta()<fPOIEta[0])||(Dp->Eta()>fPOIEta[1]) )
+      nLevel=0;
+    Double_t MassDp=Dp->InvMassDplus();
+    for(int h=0; h!=nLevel+1; ++h)
+      fMass[h]->Fill(MassDp,Dp->Pt());
+    if( (Dp->Pt()<fFlowPts[0]) || (Dp->Pt()>fFlowPts[1]) ) continue;
+    nIDs[0] = ((AliAODTrack*)Dp->GetDaughter(0))->GetID();
+    nIDs[1] = ((AliAODTrack*)Dp->GetDaughter(1))->GetID();
+    nIDs[2] = ((AliAODTrack*)Dp->GetDaughter(2))->GetID();
+    // Candidates Insertion (done in filling method: faster)
+    if(nLevel)
+      for(Int_t r=0; r!=5; ++r)
+        if( (MassDp>=fFlowBands[0][r]) && (MassDp<fFlowBands[1][r]) ) {
+          AliFlowCandidateTrack *sTrack = (AliFlowCandidateTrack*)
+                                 MakeTrack(MassDp, Dp->Pt(), Dp->Phi(),
+                                           Dp->Eta(), ndght, nIDs);
+          if(fDebug) printf("   ᶫInjecting DStar candidate on band %d \n", r);
+          for(Int_t iDau=0; iDau!=ndght; ++iDau)
+            for(Int_t iRPs=0; iRPs!=theMB[r]->NumberOfTracks(); ++iRPs) {
+              AliFlowTrack *iRP = (AliFlowTrack*) (theMB[r]->GetTrack(iRPs));
+              if(!iRP->InRPSelection()) continue;
+              if( fabs(sTrack->GetIDDaughter(iDau)) == fabs(iRP->GetID()) ) {
+                sTrack->SetDaughter(iDau,iRP);
+                iRP->SetForRPSelection(kFALSE);
+                if(fDebug) printf("    ᶫdaughter%d with fID %d was removed from this RP set\n", iDau, sTrack->GetIDDaughter(iDau));
+              }
+            }
+          theMB[r]->AddTrack(sTrack);
+        }
+    // END of Candidates Insertion
+  }
+}
+//______________________________________________________________________________
+AliFlowCandidateTrack* AliAnalysisTaskFlowD2H::MakeTrack( Double_t mass, 
+                          Double_t pt, Double_t phi, Double_t eta, 
+                          Int_t nDau, Int_t *iID ) {
+  AliFlowCandidateTrack *sTrack = new AliFlowCandidateTrack();
+  sTrack->SetMass(mass);
+  sTrack->SetPt(pt);
+  sTrack->SetPhi(phi);
+  sTrack->SetEta(eta);
+  for(Int_t iDau=0; iDau!=nDau; ++iDau)
+    sTrack->AddDaughter(iID[iDau]);
+  sTrack->SetForPOISelection(kTRUE);
+  sTrack->SetForRPSelection(kFALSE);
+  return sTrack;
+}
+//_____________________________________________________________________________
+void AliAnalysisTaskFlowD2H::Terminate(Option_t *)
+{
+
+}
+
diff --git a/PWG3/vertexingHF/charmFlow/AliAnalysisTaskFlowD2H.h b/PWG3/vertexingHF/charmFlow/AliAnalysisTaskFlowD2H.h
new file mode 100644 (file)
index 0000000..1f5384c
--- /dev/null
@@ -0,0 +1,85 @@
+/* Copyright(c) 1998-1999, ALICExperiment at CERN, All rights reserved. *
+* See cxx source for full Copyright notice */
+/* $Id: $ */
+
+#ifndef AliAnalysisTaskFlowD2H_H
+#define AliAnalysisTaskFlowD2H_H
+
+#include "AliAnalysisTaskSE.h"
+
+//==============================================================================
+// FlowD2H main task:
+// >> Make flowEvent with RPcuts and POIcuts given in constructor and passes it
+//    to the daughter tasks.
+// >> The POIcuts are polymorphic based on the AliRDHFCuts class allowing the 
+//    use of all charmed candidates reconstructed in the central barrel.
+// Author: Carlos Perez (cperez@cern.ch)
+//==============================================================================
+
+class TList;
+class TH2D;
+class TProfile;
+
+class AliAODEvent;
+
+class AliFlowEvent;
+class AliFlowCandidateTrack;
+class AliFlowTrackCuts;
+
+class AliRDHFCuts;
+class AliRDHFCutsD0toKpi;
+
+class AliAnalysisTaskFlowD2H : public AliAnalysisTaskSE {
+  public:
+    AliAnalysisTaskFlowD2H();
+    AliAnalysisTaskFlowD2H( const Char_t *name, AliFlowTrackCuts *cutsRPs, 
+                            AliRDHFCuts *cutsPOIs, Int_t specie );
+    void SetDebug() {fDebug = true;}
+    void SetPOIEtaRange( Double_t minEta, Double_t maxEta )
+                          { fPOIEta[0] = minEta; fPOIEta[1] = maxEta; }
+    void SetFlowEtaRangeAB( Double_t minA, Double_t maxA, Double_t minB, Double_t maxB )
+                          { fFlowEta[0] = minA; fFlowEta[1] = maxA;
+                            fFlowEta[2] = minB; fFlowEta[3] = maxB; }
+    void SetFlowPtRange( Int_t minPt, Int_t maxPt )
+                       { fFlowPts[0] = minPt; fFlowPts[1] = maxPt; }
+    void SetFlowBandRange( Int_t band, Double_t minMass, Double_t maxMass )
+                         { fFlowBands[0][band] = minMass;
+                           fFlowBands[1][band] = maxMass; } //TODO
+    virtual ~AliAnalysisTaskFlowD2H();
+    virtual void UserCreateOutputObjects();
+    virtual void UserExec(Option_t *);
+    virtual void Terminate(Option_t *);
+    virtual void NotifyRun();
+
+  private:
+    AliAnalysisTaskFlowD2H(const AliAnalysisTaskFlowD2H& analysisTask);
+    AliAnalysisTaskFlowD2H& operator=(const AliAnalysisTaskFlowD2H& analysisTask);
+    void AddHistograms();
+    void FillD0toKpi(      AliAODEvent *aod, AliFlowEvent *mb[5]);
+    void FillD0toKpipipi(  AliAODEvent *aod, AliFlowEvent *mb[5]);
+    void FillDStartoKpipi( AliAODEvent *aod, AliFlowEvent *mb[5]);
+    void FillDplustoKpipi( AliAODEvent *aod, AliFlowEvent *mb[5]);
+    void FillDstoKKpi(     AliAODEvent *aod, AliFlowEvent *mb[5]);
+    void FillJpsitoee(     AliAODEvent *aod, AliFlowEvent *mb[5]);
+    void FillLctoV0(       AliAODEvent *aod, AliFlowEvent *mb[5]);
+    void FillLctopKpi(     AliAODEvent *aod, AliFlowEvent *mb[5]);
+    AliFlowCandidateTrack* MakeTrack( Double_t mass, Double_t pt,
+                                      Double_t phi, Double_t eta,
+                                      Int_t nDaughters, Int_t *iID );
+    AliFlowTrackCuts *fCutsRP;  // cuts for RPs
+    AliRDHFCuts      *fCutsPOI; // cuts for POIs
+    Int_t  fSource;             // AliRDHFCuts::ESele
+    Bool_t fDebug;              // fully talkative task
+    TList *fHList;    // List for histos
+    TProfile *fAnaCuts; // store analysis related cuts
+    TH2D  *fEvent[2]; // Events histogram
+    TH2D  *fMass[2];  // Mass spectra
+    Double_t fPOIEta[2];        // Eta cut for POI
+    Double_t fFlowEta[4];       // SP subEvents
+    Int_t fFlowPts[2];          // Pt range
+    Double_t fFlowBands[2][5];  // Mass bands TODO
+
+  ClassDef(AliAnalysisTaskFlowD2H, 1);
+};
+
+#endif
diff --git a/PWG3/vertexingHF/charmFlow/AliAnalysisTaskSEHFv2.cxx b/PWG3/vertexingHF/charmFlow/AliAnalysisTaskSEHFv2.cxx
new file mode 100644 (file)
index 0000000..e6c697e
--- /dev/null
@@ -0,0 +1,1078 @@
+/**************************************************************************
+ * Copyright(c) 1998-2010, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/////////////////////////////////////////////////////////////
+//
+// AliAnalysisTaskSEHFv2 gives the needed tools for the D 
+// mesons v2 analysis with event plane method
+// Authors: Chiara Bianchin, cbianchi@pd.infn.it, 
+// Robert Grajcarek, grajcarek@physi.uni-heidelberg.de
+// Giacomo Ortona, ortona@to.infn.it
+// Carlos Eugenio Perez Lara, carlos.eugenio.perez.lara@cern.ch
+// 
+/////////////////////////////////////////////////////////////
+
+#include <Riostream.h>
+#include <TClonesArray.h>
+#include <TCanvas.h>
+#include <TList.h>
+#include <TFile.h>
+#include <TString.h>
+#include <TH1F.h>
+#include <TH2F.h>
+#include <TGraphErrors.h>
+#include <TGraph.h>
+#include <TDatabasePDG.h>
+#include <TRandom3.h>
+#include <TVector2.h>
+#include <TArrayF.h>
+
+#include <AliLog.h>
+#include <AliAnalysisDataSlot.h>
+#include <AliAnalysisDataContainer.h>
+#include "AliAnalysisManager.h"
+#include "AliAODHandler.h"
+#include "AliAODEvent.h"
+#include "AliAODVertex.h"
+#include "AliAODTrack.h"
+#include "AliAODMCHeader.h"
+#include "AliAODMCParticle.h"
+#include "AliAODRecoDecayHF3Prong.h"
+#include "AliAODRecoDecayHF.h"
+#include "AliAODRecoDecayHF2Prong.h"
+#include "AliAODRecoDecayHF4Prong.h"
+#include "AliAODRecoCascadeHF.h"
+#include "AliAODHFUtil.h"
+
+#include "AliAnalysisTaskSE.h"
+#include "AliRDHFCutsDplustoKpipi.h"
+#include "AliRDHFCutsD0toKpipipi.h"
+#include "AliRDHFCutsDstoKKpi.h"
+#include "AliRDHFCutsDStartoKpipi.h"
+#include "AliRDHFCutsD0toKpi.h"
+#include "AliRDHFCutsLctopKpi.h"
+
+#include "AliHFMassFitter.h"
+#include "AliEventplane.h"
+#include "AliFlowTrack.h"
+#include "AliFlowVector.h"
+#include "AliFlowTrackCuts.h"
+#include "AliFlowEvent.h"
+
+#include "AliAnalysisTaskSEHFv2.h"
+
+ClassImp(AliAnalysisTaskSEHFv2)
+
+
+//________________________________________________________________________
+AliAnalysisTaskSEHFv2::AliAnalysisTaskSEHFv2():
+AliAnalysisTaskSE(),
+  fhEventsInfo(0),
+  fOutput(0),
+  fRDCuts(0),
+  fParHist(0),
+  fLowmasslimit(1.765),
+  fUpmasslimit(1.965),
+  fNPtBins(1),
+  fNPhiBinLims(2),
+  fPhiBins(0),
+  fCentLowLimit(0),
+  fCentUpLimit(100),
+  fNMassBins(200),
+  fReadMC(kFALSE),
+  fDecChannel(0),
+  fUseV0EP(kFALSE),
+  fV0EPorder(2)
+{
+  // Default constructor
+}
+
+//________________________________________________________________________
+AliAnalysisTaskSEHFv2::AliAnalysisTaskSEHFv2(const char *name,AliRDHFCuts *rdCuts,Int_t decaychannel, Int_t nlimsphibin, Float_t *phibinlimits,TH2D **histPar):
+  AliAnalysisTaskSE(name),
+  fhEventsInfo(0),
+  fOutput(0),
+  fRDCuts(rdCuts),
+  fParHist(0),
+  fLowmasslimit(0),
+  fUpmasslimit(0),
+  fNPtBins(1),
+  fNPhiBinLims(2),
+  fPhiBins(0),
+  fCentLowLimit(0),
+  fCentUpLimit(100),
+  fNMassBins(200),
+  fReadMC(kFALSE),
+  fDecChannel(decaychannel),
+  fUseV0EP(kFALSE),
+  fV0EPorder(2)
+{
+
+  Int_t pdg=421;
+  switch(fDecChannel){
+  case 0:
+    pdg=411;
+    break;
+  case 1:
+    pdg=421;
+    break;
+  case 2:
+    pdg=413;
+    break;
+  case 3:
+    pdg=431;
+    break;
+  case 4:
+    pdg=421;
+    break;
+  case 5:
+    pdg=4122;
+    break;
+  }
+  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");
+  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
+  DefineOutput(1,TH1F::Class());   //Info on the number of events etc.
+  // Output slot #2 writes into a TList container
+  DefineOutput(2,TList::Class());  //Main output
+  // Output slot #3 writes into a AliRDHFCuts container (cuts)
+  switch(fDecChannel){
+  case 0:
+    DefineOutput(3,AliRDHFCutsDplustoKpipi::Class());  //Cut object for Dplus
+    break;
+  case 1:
+    DefineOutput(3,AliRDHFCutsD0toKpi::Class());  //Cut object for D0
+    break;
+  case 2:
+    DefineOutput(3,AliRDHFCutsDStartoKpipi::Class());  //Cut object for D*
+    break;
+  }
+  //DefineOutput(4,AliFlowEventSimple::Class());
+  DefineOutput(4,TList::Class());
+}
+
+//________________________________________________________________________
+AliAnalysisTaskSEHFv2::~AliAnalysisTaskSEHFv2()
+{
+  // Destructor
+  if (fOutput) {
+    delete fOutput;
+    fOutput = 0;
+  }
+
+  if(fhEventsInfo){
+    delete fhEventsInfo;
+    fhEventsInfo=0;
+  }
+
+  if(fRDCuts) {
+    delete fRDCuts;
+    fRDCuts= 0;
+  } 
+
+  if(fParHist) {
+    delete fParHist;
+    fParHist= 0;
+  } 
+  for(Int_t i=0;i<6;i++){
+    if(fHistvzero[i]) {
+      delete fHistvzero[i];
+      fHistvzero[i]=0x0;
+    }
+  } 
+  
+}
+
+//_________________________________________________________________
+void  AliAnalysisTaskSEHFv2::SetMassLimits(Float_t range, Int_t pdg){
+  Float_t mass=0;
+  Int_t abspdg=TMath::Abs(pdg);
+  mass=TDatabasePDG::Instance()->GetParticle(abspdg)->Mass();
+  fUpmasslimit = mass+range;
+  fLowmasslimit = mass-range;
+}
+//_________________________________________________________________
+void  AliAnalysisTaskSEHFv2::SetMassLimits(Float_t lowlimit, Float_t uplimit){
+  if(uplimit>lowlimit)
+    {
+      fUpmasslimit = uplimit;
+      fLowmasslimit = lowlimit;
+    }
+}
+
+
+//________________________________________________________________________
+void AliAnalysisTaskSEHFv2::LocalInit()
+{
+  // Initialization
+
+  if(fDebug > 1) printf("AnalysisTaskSEHFv2::Init() \n");
+
+  fRDCuts->SetMinCentrality(10);
+  fRDCuts->SetMaxCentrality(80);
+
+  switch(fDecChannel){
+  case 0:
+    {
+      AliRDHFCutsDplustoKpipi* copycut=new AliRDHFCutsDplustoKpipi(*(static_cast<AliRDHFCutsDplustoKpipi*>(fRDCuts)));
+      // Post the data
+      PostData(3,copycut);
+    }
+    break;
+  case 1:
+    {
+      AliRDHFCutsD0toKpi* copycut=new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fRDCuts)));
+      // Post the data
+      PostData(3,copycut);
+    }
+    break;
+  case 2:
+    {
+      AliRDHFCutsDStartoKpipi* copycut=new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fRDCuts)));
+      // Post the data
+      PostData(3,copycut);
+    }
+    break;
+  default:
+    return;
+  }
+  return;
+}
+//________________________________________________________________________
+void AliAnalysisTaskSEHFv2::UserCreateOutputObjects()
+{
+  // Create the output container
+  if(fDebug > 1) printf("AnalysisTaskSEHFv2::UserCreateOutputObjects() \n");
+
+  fhEventsInfo=new TH1F(GetOutputSlot(1)->GetContainer()->GetName(), "Number of AODs scanned",7,-0.5,6.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(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()->SetNdivisions(1,kFALSE);
+
+
+  // Several histograms are more conveniently managed in a TList
+  fOutput = new TList();
+  fOutput->SetOwner();
+  fOutput->SetName("MainOutput");
+  
+  for(Int_t icentr=10;icentr<=80;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);
+
+      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);
+       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* 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);
+
+    TH1F* hEvPlane=new TH1F(Form("hEvPlane%s",centrname.Data()),Form("Event plane angle %s;#phi Ev Plane;Entries",centrname.Data()),200,0.,TMath::Pi());
+    fOutput->Add(hEvPlane);
+
+    //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* 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);
+  }
+  
+  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);
+  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;
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskSEHFv2::UserExec(Option_t */*option*/)
+{
+  // Execute analysis for current event:
+  // heavy flavor candidates association to MC truth
+   
+  AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
+  if(fDebug>2) printf("Analysing decay %d\n",fDecChannel);
+  // Post the data already here
+  PostData(1,fhEventsInfo);
+  PostData(2,fOutput);
+
+  TClonesArray *arrayProng =0;
+  Int_t absPdgMom=0;
+  if(!aod && AODEvent() && IsStandardAOD()) { 
+    // In case there is an AOD handler writing a standard AOD, use the AOD 
+    // event in memory rather than the input (ESD) event.    
+    aod = dynamic_cast<AliAODEvent*> (AODEvent());
+    // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
+    // have to taken from the AOD event hold by the AliAODExtension
+    AliAODHandler* aodHandler = (AliAODHandler*) 
+      ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
+    if(aodHandler->GetExtensions()) {
+      
+      AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
+      AliAODEvent *aodFromExt = ext->GetAOD();
+   
+      
+      if(fDecChannel==0){
+       absPdgMom=411;
+       arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
+      }
+      if(fDecChannel==1){
+       absPdgMom=421;
+       arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
+      }
+      if(fDecChannel==2){
+       absPdgMom=413;
+       arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
+      }
+    }
+  } else {
+    if(fDecChannel==0){
+      absPdgMom=411;
+      arrayProng=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
+    }
+    if(fDecChannel==1){
+      absPdgMom=421;
+      arrayProng=(TClonesArray*)aod->GetList()->FindObject("D0toKpi");
+    }
+    if(fDecChannel==2){
+      absPdgMom=413;
+      arrayProng=(TClonesArray*)aod->GetList()->FindObject("Dstar");
+    }
+  }
+
+  if(!arrayProng) {
+    AliError("AliAnalysisTaskSEHFv2::UserExec:Branch not found!\n");
+    return;
+  }
+  
+  // fix for temporary bug in ESDfilter 
+  // the AODs with null vertex pointer didn't pass the PhysSel
+  if(!aod->GetPrimaryVertex() || TMath::Abs(aod->GetMagneticField())<0.001) return;
+
+  TClonesArray *arrayMC=0;
+  AliAODMCHeader *mcHeader=0;
+  
+  // load MC particles
+  if(fReadMC){
+    
+    arrayMC =  (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
+    if(!arrayMC) {
+      AliWarning("AliAnalysisTaskSEHFv2::UserExec:MC particles branch not found!\n");
+      //    return;
+    }
+    
+    // load MC header
+    mcHeader =  (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
+    if(!mcHeader) {
+      AliError("AliAnalysisTaskSEHFv2::UserExec:MC header branch not found!\n");
+      return;
+    }
+  }
+
+  fhEventsInfo->Fill(0); // count event
+
+  AliAODRecoDecayHF *d=0;
+
+  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);
+
+  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);
+    return;
+  }
+
+  AliEventplane *pl=0x0;
+  TVector2* q=0x0;
+  Double_t rpangleevent=0;
+  Double_t eventplane=0;
+  //determine centrality bin
+  Float_t centr=fRDCuts->GetCentrality(aod);
+  Int_t icentr=0;
+  for(Int_t ic=5;ic<=80;ic=ic+5){
+    if(ic>centr){
+      icentr=ic;
+      break;
+    }
+  }
+  TString centrbinname=Form("centr%d_%d",icentr-5,icentr);
+
+  if(fUseV0EP){
+    rpangleevent=GetEventPlaneFromV0(aod);
+    eventplane=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;
+    }
+    q = pl->GetQVector();
+    rpangleevent = pl->GetEventplane("Q"); // 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);
+  }
+  ((TH1F*)fOutput->FindObject(Form("hEvPlane%s",centrbinname.Data())))->Fill(rpangleevent);
+
+  for (Int_t iCand = 0; iCand < nCand; 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);
+    if(isSelected && isFidAcc && isSelBit) {
+      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;
+      }
+
+      if(!fUseV0EP) {
+       eventplane = GetEventPlaneForCandidate(d,q,pl); // 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);
+      }
+
+      Float_t phi=d->Phi();
+      ((TH1F*)fOutput->FindObject(Form("hPhi_pt%d%s",ptbin,centrbinname.Data())))->Fill(phi);
+      
+      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
+  }
+  
+  PostData(1,fhEventsInfo);
+  PostData(2,fOutput);
+  //PostData(4,flowEvent);
+  return;
+}
+
+//***************************************************************************
+
+// Methods used in the UserExec
+
+void AliAnalysisTaskSEHFv2::CalculateInvMasses(AliAODRecoDecayHF* d,Float_t*& masses,Int_t& nmasses){
+  //Calculates all the possible invariant masses for each candidate
+  //NB: the implementation for each candidate is responsibility of the corresponding developer
+
+  if(fDecChannel==0){
+    //D+ -- Giacomo
+    nmasses=1;
+    masses=new Float_t[nmasses];
+    Int_t pdgdaughters[3] = {211,321,211};
+    masses[0]=d->InvMass(3,(UInt_t*)pdgdaughters);
+  }
+  if(fDecChannel==1){
+    //D0 (Kpi)  -- Chiara
+    const Int_t ndght=2;
+    nmasses=2;
+    masses=new Float_t[nmasses];
+    Int_t pdgdaughtersD0[ndght]={211,321};//pi,K 
+    masses[0]=d->InvMass(ndght,(UInt_t*)pdgdaughtersD0); //D0
+    Int_t pdgdaughtersD0bar[ndght]={321,211};//K,pi 
+    masses[1]=d->InvMass(ndght,(UInt_t*)pdgdaughtersD0bar); //D0bar
+  }
+  if(fDecChannel==2){
+    //D* -- Robert,Yifei, Alessandro
+    nmasses=1;
+    masses=new Float_t[nmasses];
+    masses[0]=((AliAODRecoCascadeHF*)d)->DeltaInvMass();
+  } 
+}
+
+//******************************************************************************
+
+//Methods to fill the istograms with MC information, one for each candidate
+//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,Float_t* masses,Int_t isSel,Int_t icentr){
+  //D+ channel
+  if(!isSel){
+    if(fDebug>3)AliWarning("Candidate not selected\n");
+    return;
+  }
+  if(!masses){
+    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("hMphicentr%d_%d",icentr-5,icentr)))->Fill(deltaphi,masses[0]);
+  Int_t pdgdaughters[3] = {211,321,211};
+
+  if(fReadMC){
+    Int_t lab=-1;
+    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
+      ((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]);
+    } 
+  }   
+}
+
+
+void AliAnalysisTaskSEHFv2::FillD02p(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t ptbin,Float_t deltaphi,Float_t* masses,Int_t isSel,Int_t icentr){
+
+  //D0->Kpi channel
+
+  //mass histograms
+  if(!masses){
+    if(fDebug>3)AliWarning("Masses not calculated\n");
+    return;
+  }
+  Int_t phibin=GetPhiBin(deltaphi);
+  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("hMphicentr%d_%d",icentr-5,icentr)))->Fill(deltaphi,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("hMphicentr%d_%d",icentr-5,icentr)))->Fill(deltaphi,masses[1]);
+  }
+
+
+  //MC histograms
+  if(fReadMC){
+
+    Int_t matchtoMC=-1;
+
+    //D0
+    Int_t pdgdaughters[2];
+    pdgdaughters[0]=211;//pi 
+    pdgdaughters[1]=321;//K
+    Int_t nprongs=2;
+    Int_t absPdgMom=421;
+
+    matchtoMC = d->MatchToMC(absPdgMom,arrayMC,nprongs,pdgdaughters);
+
+    Int_t prongPdgPlus=421,prongPdgMinus=(-1)*421;
+    if((isSel==1 || isSel==3)){ //D0
+      if(matchtoMC>=0){
+       AliAODMCParticle *dMC = (AliAODMCParticle*)arrayMC->At(matchtoMC);
+       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 {
+         ((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("hMphiRcentr%d_%d",icentr-5,icentr)))->Fill(deltaphi,masses[0]);
+       }
+      } 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("hMphicentr%d_%d",icentr-5,icentr)))->Fill(deltaphi,masses[0]);
+      }
+    }
+    if(isSel>=2){ //D0bar
+      if(matchtoMC>=0){
+       AliAODMCParticle *dMC = (AliAODMCParticle*)arrayMC->At(matchtoMC);
+       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 {
+         ((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("hMphicentr%d_%d",icentr-5,icentr)))->Fill(deltaphi,masses[1]);
+       }
+      } 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]);
+      }
+    }
+  }
+}
+
+void AliAnalysisTaskSEHFv2::FillDstar(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t ptbin,Float_t deltaphi, Float_t* masses,Int_t isSel,Int_t icentr){
+  //D* channel
+  if(!isSel){
+    if(fDebug>3)AliWarning("Candidate not selected\n");
+    return;
+  }
+  if(!masses){
+    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("hMphicentr%d_%d",icentr-5,icentr)))->Fill(deltaphi,masses[0]);
+  Int_t pdgDgDStartoD0pi[2]={421,211};
+  Int_t pdgDgD0toKpi[2]={321,211};
+  
+  if(fReadMC){
+    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
+      ((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]);
+    } 
+  }
+  
+
+}
+
+
+//________________________________________________________________________
+Int_t AliAnalysisTaskSEHFv2::GetPhiBin(Float_t deltaphi){
+
+  //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;
+    }
+  }
+  return phibin;
+}
+//________________________________________________________________________
+// 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){
+  Float_t result=phi;
+  while(result<0){
+    result=result+TMath::Pi();
+  }
+  while(result>TMath::Pi()){
+    result=result-TMath::Pi();
+  }
+  return result;
+}
+
+//________________________________________________________________________
+Float_t AliAnalysisTaskSEHFv2::GetEventPlaneForCandidate(AliAODRecoDecayHF* d, TVector2* q,AliEventplane *pl){
+  // remove autocorrelations 
+  TArrayF* qx = pl->GetQContributionXArray();
+  TArrayF* qy = pl->GetQContributionYArray();
+  TVector2 qcopy = *q;
+  
+  if(fDecChannel==2){
+    //D* -- Yifei, Alessandro,Robert
+    AliAODRecoDecayHF2Prong* theD0particle = ((AliAODRecoCascadeHF*)d)->Get2Prong();
+    AliAODTrack *track0 = (AliAODTrack*)theD0particle->GetDaughter(0);
+    AliAODTrack *track1 = (AliAODTrack*)theD0particle->GetDaughter(1);  
+    AliAODTrack *track2 = ((AliAODRecoCascadeHF*)d)->GetBachelor(); 
+    // reduce global q vector
+
+    TVector2 q0;
+    if((track0->GetID()) < qx->fN){
+      q0.Set(qx->At(track0->GetID()),qy->At(track0->GetID()));}
+       
+    TVector2 q1;
+    if((track1->GetID()) < qx->fN){
+      q1.Set(qx->At(track1->GetID()),qy->At(track1->GetID()));}
+       
+    TVector2 q2;
+    if((track2->GetID()) < qx->fN){
+      q2.Set(qx->At(track2->GetID()),qy->At(track2->GetID()));}
+      
+    qcopy = qcopy -(q0+q1+q2);
+  
+  }
+  
+  // reduce Q vector for D+ and D0
+  
+  if(fDecChannel==1){    
+    //D0 -- Chiara
+    AliAODTrack *track0 = (AliAODTrack*)d->GetDaughter(0);
+    AliAODTrack *track1 = (AliAODTrack*)d->GetDaughter(1);  
+
+    TVector2 q0;
+    if((track0->GetID()) < qx->fN){
+      q0.Set(qx->At(track0->GetID()),qy->At(track0->GetID()));}
+       
+    TVector2 q1;
+    if((track1->GetID()) < qx->fN){
+      q1.Set(qx->At(track1->GetID()),qy->At(track1->GetID()));}
+       
+    qcopy = qcopy -(q0+q1);
+  }
+
+  if(fDecChannel==0){
+    //D+ -- Giacomo
+    AliAODTrack *track0 = (AliAODTrack*)d->GetDaughter(0);
+    AliAODTrack *track1 = (AliAODTrack*)d->GetDaughter(1);  
+    AliAODTrack *track2 = (AliAODTrack*)d->GetDaughter(2);  
+    
+    TVector2 q0;
+    if((track0->GetID()) < qx->fN){
+      q0.Set(qx->At(track0->GetID()),qy->At(track0->GetID()));}
+       
+    TVector2 q1;
+    if((track1->GetID()) < qx->fN){
+      q1.Set(qx->At(track1->GetID()),qy->At(track1->GetID()));}
+       
+    TVector2 q2;
+    if((track2->GetID()) < qx->fN){
+      q2.Set(qx->At(track2->GetID()),qy->At(track2->GetID()));}
+      
+    qcopy = qcopy -(q0+q1+q2);
+       
+  }
+
+  return qcopy.Phi()/2.;
+
+}
+//________________________________________________________________________
+Float_t AliAnalysisTaskSEHFv2::GetEventPlaneFromV0(AliAODEvent *aodEvent){
+
+  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");
+  
+  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*/)
+{
+  
+  // 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;
+}
+
+//-------------------------------------------
+/*
+  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/PWG3/vertexingHF/charmFlow/AliAnalysisTaskSEHFv2.h b/PWG3/vertexingHF/charmFlow/AliAnalysisTaskSEHFv2.h
new file mode 100644 (file)
index 0000000..58566cc
--- /dev/null
@@ -0,0 +1,97 @@
+#ifndef ALIANALYSISTASKSEHFV2_H
+#define ALIANALYSISTASKSEHFV2_H
+
+/* Copyright(c) 1998-2010, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+//*************************************************************************
+// AliAnalysisTaskSEHFv2 gives the needed tools for the D 
+// mesons v2 analysis 
+// Authors: Chiara Bianchin, cbianchi@pd.infn.it, 
+//          Robert Grajcarek, grajcarek@physi.uni-heidelberg.de
+//          Giacomo Ortona, ortona@to.infn.it,
+//          Carlos Perez Lara, carlos.eugenio.perez.lara@cern.ch
+//          Francesco Prino, prino@to.infn.it
+//
+//*************************************************************************
+
+#include "AliAnalysisTaskSE.h"
+#include "AliAnalysisVertexingHF.h"
+
+class TH1F;
+class TH2D;
+class AliMultiDimVector;
+class AliRDHFCuts;
+class TVector2;
+
+class AliAnalysisTaskSEHFv2 : public AliAnalysisTaskSE
+{
+
+ public:
+
+  enum DecChannel{kDplustoKpipi,kD0toKpi,kDstartoKpipi}; //more particles can be added
+
+  AliAnalysisTaskSEHFv2();
+  AliAnalysisTaskSEHFv2(const char *name, AliRDHFCuts *rdCuts, Int_t decaychannel,Int_t nbinsphi, Float_t *phibinlimits,TH2D** histPar);
+  virtual ~AliAnalysisTaskSEHFv2();
+
+  void SetReadMC(Bool_t readMC=kTRUE){fReadMC=readMC;}
+  void SetMassLimits(Float_t range,Int_t pdg);
+  void SetMassLimits(Float_t lowlimit, Float_t uplimit);
+  void SetNMassBins(Int_t nbins){fNMassBins=nbins;}
+  void SetUpperCentLimit(Float_t lim){fCentUpLimit = lim;}
+  void SetLowerCentLimit(Float_t lim){fCentLowLimit = lim;}
+  void SetUseV0EP(Bool_t flagV0EP){fUseV0EP=flagV0EP;}
+  void SetV0EventPlaneOrder(Int_t n){fV0EPorder=n;}
+
+  Float_t GetUpperMassLimit()const {return fUpmasslimit;}
+  Float_t GetLowerMassLimit()const {return fLowmasslimit;}
+  Int_t GetNMassBins()const {return fNMassBins;}
+  Int_t GetPhiBin(Float_t deltaphi);
+  //Float_t GetPhi02Pi(Float_t phi);
+  Float_t GetPhi0Pi(Float_t phi);
+  Float_t GetLowerCentLimit()const {return fCentLowLimit;}
+  Float_t GetUpperCentLimit()const {return fCentUpLimit;}
+  // Implementation of interface methods
+  virtual void UserCreateOutputObjects();
+  virtual void LocalInit();// {Init();}
+  virtual void UserExec(Option_t *option);
+  virtual void Terminate(Option_t *option);
+    
+ private:
+
+  AliAnalysisTaskSEHFv2(const AliAnalysisTaskSEHFv2 &source);
+  AliAnalysisTaskSEHFv2& operator=(const AliAnalysisTaskSEHFv2& source); 
+
+  void CalculateInvMasses(AliAODRecoDecayHF* d,Float_t* &masses,Int_t& nmasses);
+
+  void FillDplus(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t ptbin, Float_t dphi,Float_t* masses,Int_t isSel,Int_t icentr);
+  void FillD02p(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t ptbin, Float_t dphi,Float_t* masses, Int_t isSel,Int_t icentr);
+  void FillDstar(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t ptbin, Float_t dphi,Float_t* masses,Int_t isSel,Int_t icentr);
+  Float_t GetEventPlaneForCandidate(AliAODRecoDecayHF* d, TVector2* q,AliEventplane *pl);
+  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)
+  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
+  Int_t fNPhiBinLims;           //number of delta phi bins limits (= number of bins +1)
+  Float_t *fPhiBins;            //[fNPhiBinLims] limits of each phi bin
+  Float_t fCentLowLimit;        //lower centrality limit
+  Float_t fCentUpLimit;         //upper centrality limit
+  Int_t fNMassBins;             //number of bins in the mass histograms
+  Bool_t fReadMC;               //flag for access to MC
+  Int_t fDecChannel;            //decay channel identifier
+  Bool_t fUseV0EP;              //flag to select EP method
+  Int_t  fV0EPorder;            //harmonic for VZERO event plane
+
+  ClassDef(AliAnalysisTaskSEHFv2,1); // AliAnalysisTaskSE for the HF v2 analysis
+};
+
+#endif
diff --git a/PWG3/vertexingHF/charmFlow/AliEventPlaneResolution.cxx b/PWG3/vertexingHF/charmFlow/AliEventPlaneResolution.cxx
new file mode 100644 (file)
index 0000000..ab9a039
--- /dev/null
@@ -0,0 +1,132 @@
+/**************************************************************************
+ * Copyright(c) 2007-2009, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+#include <TMath.h>
+#include "AliEventPlaneResolution.h"
+
+/* $Id: AliITSOnlineSDDInjectors.cxx 48265 2011-03-09 23:36:06Z prino $ */
+
+///////////////////////////////////////////////////////////////////
+//                                                               //
+// Implementation of the class to compute event plane resolution //
+// for flow analyses                                             //
+// Origin: F.Prino, Torino, prino@to.infn.it                     //
+//                                                               //
+///////////////////////////////////////////////////////////////////
+
+ClassImp(AliEventPlaneResolution)
+
+//______________________________________________________________________
+AliEventPlaneResolution::AliEventPlaneResolution():TObject(),
+  fK(1),
+  fSubRes(1.)
+{
+  // Default contructor
+}
+
+
+//______________________________________________________________________
+AliEventPlaneResolution::AliEventPlaneResolution(Int_t k):
+  TObject(),
+  fK(k),
+  fSubRes(1.)
+{
+  // Standard constructor
+}
+
+
+//______________________________________________________________________
+Double_t AliEventPlaneResolution::Pol(Double_t x, Int_t k){
+  // compute chi from polynomial approximation
+  Double_t c[5];
+  if(k==1){ 
+    c[0]=0.626657; c[1]=0.; c[2]=-0.09694; c[3]=0.02754; c[4]=-0.002283;
+  }
+  else if(k==2){
+    c[0]=0.; c[1]=0.25; c[2]=-0.011414; c[3]=-0.034726; c[4]=0.006815;
+  }
+  else return -1;
+  return c[0]*x+c[1]*x*x+c[2]*x*x*x+c[3]*x*x*x*x+c[4]*x*x*x*x*x;
+}
+
+//______________________________________________________________________
+Double_t AliEventPlaneResolution:: ResolK1(Double_t x){
+  return TMath::Sqrt(TMath::Pi()/8)*x*TMath::Exp(-x*x/4)*(TMath::BesselI0(x*x/4)+TMath::BesselI1(x*x/4));
+}
+
+
+//______________________________________________________________________
+Double_t AliEventPlaneResolution::FindChi(Double_t res,  Int_t k){
+  // compute chi variable (=v2*sqrt(N)) from external values
+
+  Double_t x1=0;
+  Double_t x2=15;
+  Double_t y1,y2;
+  if(k==1){
+    y1=ResolK1(x1)-res;
+    y2=ResolK1(x2)-res;
+  }
+  else if(k==2){
+    y1=Pol(x1,2)-res;
+    y2=Pol(x2,2)-res;
+  }
+  else return -1;
+
+  if(y1*y2>0) return -1;
+  if(y1==0) return y1;
+  if(y2==0) return y2;
+  Double_t xmed,ymed;
+  Int_t jiter=0;
+  while((x2-x1)>0.0001){
+    xmed=0.5*(x1+x2);
+    if(k==1){
+      y1=ResolK1(x1)-res;
+      y2=ResolK1(x2)-res;
+      ymed=ResolK1(xmed)-res;
+    }
+    else if(k==2){
+      y1=Pol(x1,2)-res;
+      y2=Pol(x2,2)-res;
+      ymed=Pol(xmed,2)-res;
+    }
+    else return -1;
+    if((y1*ymed)<0) x2=xmed;
+    if((y2*ymed)<0)x1=xmed;
+    if(ymed==0) return xmed;
+    jiter++;
+  }
+  return 0.5*(x1+x2);
+}
+
+//______________________________________________________________________
+Double_t AliEventPlaneResolution::GetFullEvResol(Double_t resSub, Int_t k){
+  // computes event plane resolution starting from sub event resolution
+  Double_t chisub=FindChi(resSub,k);
+  Double_t chifull=chisub*TMath::Sqrt(2);
+  if(k==1) return ResolK1(chifull);
+  else if(k==2) return Pol(chifull,2);
+  else{
+    printf("k should be <=2\n");
+    return 1.;
+  }
+}
+
+//______________________________________________________________________
+Double_t AliEventPlaneResolution::GetFullEvResol(const TH1F* hSubEvCorr, Int_t k){
+  // computes event plane resolution starting from sub event correlation histogram
+  if(!hSubEvCorr) return 1.;
+  Double_t resSub=GetSubEvResol(hSubEvCorr);
+  return GetFullEvResol(resSub,k);
+}
diff --git a/PWG3/vertexingHF/charmFlow/AliEventPlaneResolution.h b/PWG3/vertexingHF/charmFlow/AliEventPlaneResolution.h
new file mode 100644 (file)
index 0000000..0dab03b
--- /dev/null
@@ -0,0 +1,53 @@
+#ifndef ALIEVENTPLANERESOLUTION_H
+#define ALIEVENTPLANERESOLUTION_H
+
+
+/* $Id: $ */
+
+///////////////////////////////////////////////////////////////////
+//                                                               //
+// Class to compute event plane resolution for flow analyses     //
+// Origin: F.Prino, Torino, prino@to.infn.it                     //
+//                                                               //
+///////////////////////////////////////////////////////////////////
+
+#include "TObject.h"
+#include "TH1F.h"
+
+class AliEventPlaneResolution : public TObject{
+ public:
+  AliEventPlaneResolution();
+  AliEventPlaneResolution(Int_t k);
+  virtual ~AliEventPlaneResolution() {};
+
+  void SetK(Int_t k){fK=k;}
+  void SetSubEvResol(Double_t res){fSubRes=res;}
+  void SetSubEventHisto(const TH1F* hSub){
+    fSubRes=GetSubEvResol(hSub);
+  }
+
+  Int_t GetK() const  {return fK;}
+  Double_t GetSubEvResol() const  {return fSubRes;}
+
+  Double_t Pol(Double_t x) const {return Pol(x,fK);}
+  Double_t FindChi() const {return FindChi(fSubRes,fK);}
+  Double_t GetFullEvResol() const {return GetFullEvResol(fSubRes,fK);}
+
+  static Double_t FindChi(Double_t res,  Int_t k=1);
+  static Double_t Pol(Double_t x, Int_t k);
+  static Double_t ResolK1(Double_t x);
+  static Double_t GetSubEvResol(const TH1F* hSubEvCorr){
+    if(hSubEvCorr) return TMath::Sqrt(hSubEvCorr->GetMean());
+    else return 1.;
+  }
+  static Double_t GetFullEvResol(Double_t resSub, Int_t k=1);
+  static Double_t GetFullEvResol(const TH1F* hSubEvCorr, Int_t k=1);
+
+ private:
+
+  Int_t fK;             // ratio of measured harmonic to event plane harmonic
+  Double_t fSubRes;     // sub-event resolution = sqrt(<cos[n(phiA-phiB)] >)
+
+  ClassDef(AliEventPlaneResolution,0) 
+};
+#endif
diff --git a/PWG3/vertexingHF/charmFlow/DmesonsFlowAnalysis.C b/PWG3/vertexingHF/charmFlow/DmesonsFlowAnalysis.C
new file mode 100644 (file)
index 0000000..11ca1b0
--- /dev/null
@@ -0,0 +1,1643 @@
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include <Riostream.h>
+#include <TFile.h>
+#include <TString.h>
+#include <TH2F.h>
+#include <TH1F.h>
+#include <TF1.h>
+#include <TGraph.h>
+#include <TGraphErrors.h>
+#include <TMultiGraph.h>
+#include <TDirectoryFile.h>
+#include <TList.h>
+#include <TCanvas.h>
+#include <TLegend.h>
+#include <TLegendEntry.h>
+#include <TPaveText.h>
+#include <TStyle.h>
+#include <TASImage.h>
+#include <TPad.h>
+#include <TROOT.h>
+#include <TDatabasePDG.h>
+#include <TParameter.h>
+#include <TLatex.h>
+
+#include <AliHFMassFitter.h>
+#include "AliRDHFCutsD0toKpi.h"
+#include "AliRDHFCutsDplustoKpipi.h"
+#include "AliRDHFCutsDStartoKpipi.h"
+
+#endif
+
+//methods for the analysis of AliAnalysisTaskSEHFv2 output
+//Authors: Chiara Bianchin cbianchi@pd.infn.it
+//         Giacomo Ortona  ortona@to.infn.it
+//         Francesco Prino prino@to.infn.it
+
+//global variables to be set
+// const Int_t nptbinsnew=6;
+// Float_t ptbinsnew[nptbinsnew+1]={2,3,4,5,8,12,999.};
+const Int_t nptbinsnew=3;
+Float_t ptbinsnew[nptbinsnew+1]={2,5,8,16.};
+//const Int_t nptbinsnew=4;
+//Float_t ptbinsnew[nptbinsnew+1]={2,4,5,8,999.};
+//Float_t ptbinsnew[nptbinsnew+1]={2,3,5,8,999.};
+// const Int_t nptbinsnew=5;
+// Float_t ptbinsnew[nptbinsnew+1]={2,3,4,5,8,999.};
+// const Int_t nptbinsnew=2;
+// Float_t ptbinsnew[nptbinsnew+1]={2,8,999.};
+Int_t fittype=0;
+//Int_t rebin[nptbinsnew]={4,4,4,5};
+Int_t rebin[nptbinsnew]={4,4,4};
+Double_t nsigma=3;
+
+//methods
+Bool_t ReadFile(TList* &list,TH1F* &hstat,AliRDHFCuts* &cutobj,TString listname,TString partname,TString path="./",TString filename="AnalysisResults.root");
+void WriteCanvas(TCanvas* cv,TString text);
+Int_t FindPtBin(Int_t nbins, Float_t* array,Float_t value);
+void InOutPic(TVirtualPad *c,Int_t inout=0,TString where="tr");//inout: 0=IN, 1=OUT 
+void PhiBinPic(TVirtualPad *c,Int_t angle,TString where);
+Double_t FindChi(Double_t res, Int_t k);
+Double_t Pol(Double_t x, Int_t k);
+Double_t ResolK1(Double_t x);
+Double_t ComputeResol(Double_t resSub, Int_t k);
+
+
+//methods implementation
+
+Bool_t ReadFile(TList* &list,TH1F* &hstat,AliRDHFCuts* &cutsobj,TString listname,TString partname,TString path,TString filename){
+
+  TString hstatname="hEventsInfo",dirname="PWG3_D2H_HFv2";
+  filename.Prepend(path);
+  listname+=partname;
+  hstatname+=partname;
+  // TString tmpsuff="NoCos"; //"Nod0d0"
+  // listname+=tmpsuff;
+  // hstatname+=tmpsuff;
+
+  TFile* f=new TFile(filename.Data());
+  if(!f){
+    cout<<filename.Data()<<" not found"<<endl;
+    return kFALSE;
+  }
+  TDirectoryFile* dir=(TDirectoryFile*)f->Get(dirname);
+  if(!f){
+    cout<<dirname.Data()<<" not found  in "<<filename.Data()<<endl;
+    return kFALSE;
+  }
+
+  list=(TList*)dir->Get(listname);
+  if(!list){
+    cout<<"List "<<listname.Data()<<" not found"<<endl;
+    dir->ls();
+    return kFALSE;
+  }
+
+  hstat=(TH1F*)dir->Get(hstatname);
+  if(!hstat){
+    cout<<hstatname.Data()<<" not found"<<endl;
+    return kFALSE;
+  }
+  
+  if(partname.Contains("D0")) cutsobj=((AliRDHFCutsD0toKpi*)dir->Get(dir->GetListOfKeys()->At(2)->GetName()));
+  if(partname.Contains("Dplus")) cutsobj=((AliRDHFCutsDplustoKpipi*)dir->Get(dir->GetListOfKeys()->At(2)->GetName()));
+  if(partname.Contains("Dstar")) cutsobj=((AliRDHFCutsDStartoKpipi*)dir->Get(dir->GetListOfKeys()->At(2)->GetName()));
+               
+  if(!cutsobj){
+    cout<<"Cut object not found! check position of the key!"<<endl;
+    return kFALSE;
+  }
+  
+  return kTRUE;
+}
+
+Int_t FindPtBin(Int_t nbins, Float_t* array,Float_t value){
+  for (Int_t i=0;i<nbins;i++){
+    //cout<<value<<" from "<<array[i]<<" to "<<array[i+1]<<"?"<<endl;
+    if(value>=array[i] && value<array[i+1]){
+      return i;
+    }
+  }
+  cout<<value<< " out of range "<<array[0]<<", "<<array[nbins]<<endl;
+  return -1;
+
+}
+
+void WriteCanvas(TCanvas* cv,TString text){
+  cv->SaveAs(Form("%s%s.png",cv->GetName(),text.Data()));
+}
+
+void InOutPic(TVirtualPad *c,Int_t inout,TString where){
+  TPad *myPadLogo = 0x0;
+  if(where=="tr"){
+    myPadLogo = new TPad("myPadPic", "Pad for In/Out pic",0.66,0.68,0.86,0.88);
+  }
+
+  myPadLogo->SetFillColor(0); 
+  myPadLogo->SetBorderMode(0);
+  myPadLogo->SetBorderSize(2);
+  myPadLogo->SetFrameBorderMode(0);
+  myPadLogo->SetLeftMargin(0.0);
+  myPadLogo->SetTopMargin(0.0);
+  myPadLogo->SetBottomMargin(0.0);
+  myPadLogo->SetRightMargin(0.0);
+  myPadLogo->Draw();
+  myPadLogo->cd();
+  TASImage *myAliceLogo = new TASImage(Form("/home/cbianchi/macros/sketch%s.png",inout==0 ? "IN" : "OUT"));
+  myAliceLogo->Draw();
+  c->cd();
+  TPaveText* t1=0x0;//=new TPaveText(0.62,0.63,0.89,0.71,"NDC");
+  if(where=="tr"){
+    t1=new TPaveText(0.65,0.58,0.85,0.67,"NDC");
+  }
+  t1->SetFillStyle(0);
+  t1->AddText(Form("%s-PLANE",inout==0 ? "IN" : "OUT-OF"));
+  if(inout==0)t1->SetTextColor(kBlue);
+  else t1->SetTextColor(kRed);
+  t1->Draw();
+}
+
+void PhiBinPic(TVirtualPad *c,Int_t angle,TString where){
+
+  TString picname="angles";
+  if(angle==0) picname.Append("0pi4");
+  if(angle==1) picname.Append("pi4pi2");
+  if(angle==2) picname.Append("pi232pi");
+  if(angle==3) picname.Append("32pipi");
+  picname.Append(".png");
+
+  TPad *myPadLogo = 0x0;
+  if(where=="tr"){
+    myPadLogo = new TPad("myPadPic", "Pad for bin of phi pic",0.5,0.75,0.8,0.8);
+  }
+
+  myPadLogo->SetFillColor(0); 
+  myPadLogo->SetBorderMode(0);
+  myPadLogo->SetBorderSize(2);
+  myPadLogo->SetFrameBorderMode(0);
+  myPadLogo->SetLeftMargin(0.0);
+  myPadLogo->SetTopMargin(0.0);
+  myPadLogo->SetBottomMargin(0.0);
+  myPadLogo->SetRightMargin(0.0);
+  myPadLogo->Draw();
+  myPadLogo->cd();
+  TASImage *myAliceLogo = new TASImage(Form("/home/cbianchi/macros/%s",picname.Data()));
+  c->cd();
+  myAliceLogo->Draw();
+  
+
+  /*
+  TPaveText* t1=0x0;//=new TPaveText(0.62,0.63,0.89,0.71,"NDC");
+  if(where=="tr"){
+    t1=new TPaveText(0.65,0.58,0.85,0.67,"NDC");
+  }
+  t1->SetFillStyle(0);
+  t1->AddText(Form("%s-PLANE",inout==0 ? "IN" : "OUT-OF"));
+  if(inout==0)t1->SetTextColor(kBlue);
+  else t1->SetTextColor(kRed);
+  t1->Draw();
+  */
+}
+
+
+void DmesonsSignalExtraction(TString partname="D0",Int_t mincentr=30,Int_t maxcentr=70,TString textleg="",TString path="./",Bool_t official=kFALSE){
+
+  //read the output of HFv2task and extract the signal in pt bins in plane and out of plane and in bins of phi pt integrated
+
+  gStyle->SetCanvasColor(0);
+  gStyle->SetTitleFillColor(0);
+  gStyle->SetStatColor(0);
+  gStyle->SetPalette(1);
+  gStyle->SetOptStat(0);
+
+  Int_t info=1;
+  if (official) info=0;
+
+  if(textleg==""){
+    textleg=Form("centr%d-%d",mincentr,maxcentr);
+  }
+  //  Double_t pi=TMath::Pi();
+  TString listname="coutputv2";
+
+  TList* list;
+  TH1F * hstat;
+  AliRDHFCuts* cutobj;
+
+  Bool_t isRead=ReadFile(list,hstat,cutobj,listname,partname,path);
+  if(!isRead) return;
+  if(!list || !hstat || !cutobj){
+    cout<<":-( null pointers..."<<endl;
+    return;
+  }
+  //uncomment when problems with the cut obj
+  /* 
+  const Int_t nptbins=12;
+  Float_t ptbins[nptbins+1]={0,1,2,3,4,5,6,8,12,16,999};
+  */
+
+  Float_t* ptbins=cutobj->GetPtBinLimits();
+  const Int_t nptbins=cutobj->GetNPtBins();
+
+  TH1F* hnphibins=(TH1F*)list->FindObject("hPhiBins");
+  Int_t nphibins=hnphibins->GetNbinsX();
+  Double_t phibins[nphibins],dx[nphibins];
+  cout<<"(";
+  for(Int_t i=0;i<nphibins;i++){
+    Double_t width=hnphibins->GetBinWidth(i+1);
+    phibins[i]=hnphibins->GetBinLowEdge(i+1)+width/2.;
+    dx[i]=width/2.;
+    cout<<phibins[i]<<",";
+  }
+  cout<<")"<<endl;
+  cout<<"Number of phi bins = "<<nphibins<<endl;
+
+  cout<<"Number of pt bins = "<<cutobj->GetNPtBins()<<endl;
+  Double_t deltapt[nptbinsnew],pt[nptbinsnew],ptledges[nptbinsnew+1];
+  for(Int_t j=0;j<nptbinsnew;j++){
+    deltapt[j]=(ptbinsnew[j+1]-ptbinsnew[j])/2.;
+    pt[j]=ptbinsnew[j]+deltapt[j];
+    ptledges[j]=ptbinsnew[j];
+
+    if(j==nptbinsnew-1 && deltapt[j]>10) {
+      pt[j]=ptbinsnew[j]+10;
+      deltapt[j]=10;
+    }
+  }
+  ptledges[nptbinsnew]=ptledges[nptbinsnew-1]+deltapt[nptbinsnew-1];
+
+  //invariant mass pt integrated, in bins of phi
+  TH1F** hmassptint=new TH1F*[nphibins];
+  //invariant mass in bins of pt and in/out plane
+  TH1F*** hmassnonintphi=new TH1F**[2];
+
+  //init
+  for(Int_t j=0;j<2;j++){
+    hmassnonintphi[j]=new TH1F*[nptbinsnew];
+    for(Int_t i=0;i<nptbinsnew;i++){
+      hmassnonintphi[j][i]=0x0;
+      cout<<"Init "<<hmassnonintphi[j][i]<<endl;
+    }
+  }
+  for(Int_t j=0;j<nphibins;j++){
+    hmassptint[j]=0x0;
+  }  
+
+
+  //output file
+  TFile* fout=new TFile(Form("HistoInputv2Calc%s.root",textleg.Data()),"recreate");
+  TClonesArray ptblim("TParameter<float>",nptbinsnew+2); //the first is the number of bins
+  new(ptblim[0])TParameter<float>("nptbins",nptbinsnew);
+  for(Int_t ival=0;ival<nptbinsnew+1;ival++){
+    TString name=Form("ptbin%d",ival);
+    new(ptblim[ival+1])TParameter<float>(name.Data(),ptbinsnew[ival]);
+  }
+  fout->cd();
+  ptblim.Write();
+  cutobj->Write();
+
+
+  //reading the hostograms and filling hmassnonintphi and hmassptint
+
+  Bool_t isInPlane=kFALSE, isOutOfPlane=kFALSE;
+  Int_t indx=-1;
+
+  for(Int_t i=0;i<nphibins;i++){ //loop on phi
+    cout<<" -- phi bin "<<i<<"/"<<nphibins-1<<endl;
+
+    for(Int_t j=0;j<nptbins;j++){ //loop on pt
+      cout<<"   ** pt bin "<<j<<"/"<<nptbins-1<<endl;
+  
+      isInPlane=kFALSE, isOutOfPlane=kFALSE;
+
+      TH1F* h=(TH1F*)list->FindObject(Form("hMass_pt%dphi%dcentr%d_%d",j,i,mincentr,mincentr+5));
+      if(!h){
+       cout<<"qui Histogram pt "<<j<<", phi "<<i<<", centr "<<mincentr<<"-"<<mincentr+5<<" not found"<<endl;
+       //list->ls();
+       continue;
+       //return;
+      }
+      for(Int_t icentr=mincentr+5;icentr<maxcentr;icentr=icentr+5){
+       TH1F* hcentr05=(TH1F*)list->FindObject(Form("hMass_pt%dphi%dcentr%d_%d",j,i,icentr,icentr+5));
+       if(!hcentr05){
+         cout<<"Histogram pt "<<j<<", phi "<<i<<", centr "<<icentr<<"-"<<icentr+5<<" not found"<<endl;
+         return;
+       }
+       h->Add(hcentr05);
+      }
+      Int_t kptnew=FindPtBin(nptbinsnew,ptbinsnew,ptbins[j]);
+      cout<<"----*** pt"<<kptnew<<endl;
+      if(kptnew==-1) continue;
+      //Filling in plane and out of plane
+      if((phibins[i]>0 && phibins[i]<= 0.79) || (phibins[i]> 2.36 && phibins[i]<3.14)) {isInPlane=kTRUE; indx=0;}
+      if(phibins[i]> 0.79 && phibins[i]<= 2.36) {isOutOfPlane=kTRUE; indx=1;cout<<" IS OUT of plane"<<endl;}
+      if(!hmassnonintphi[indx][kptnew]) {
+
+       cout<<"\tpt"<<j<<", phi"<<i<<" in phi "<<i<<" ptnew "<<kptnew<<endl;
+       hmassnonintphi[indx][kptnew]=(TH1F*)h->Clone(Form("hMass_ptnw%dphi%d",kptnew,isInPlane? 0 : 1));
+       hmassnonintphi[indx][kptnew]->SetTitle(Form("Mass ptnw%d phi%s",kptnew,isInPlane ? "In plane" : "Out of plane"));
+
+      }
+      else {
+       cout<<"\tadding pt"<<j<<", phi"<<i<<" in phi "<<i<<" ptnew "<<kptnew<<endl;
+       hmassnonintphi[indx][kptnew]->Add(h);
+      }
+
+
+      //pt integrated mass plots
+      if(!hmassptint[i]) {
+       cout<<"Filling histo ptint["<<i<<"] with hMass_pt"<<j<<"phi"<<i<<endl;
+       hmassptint[i]=(TH1F*)h->Clone(Form("hphi_%d",i));
+       hmassptint[i]->SetTitle(Form("Invariant Mass (Phi %d, Pt integr)",i));
+       //cout<<"Entries "<<hmassptint[i]->GetEntries()<<endl;
+      }
+      else {
+       cout<<"Adding histo hMass_pt"<<j<<"phi"<<i<<" to ptint["<<i<<"]"<<endl;
+       hmassptint[i]->Add(h);
+      }
+    } //end pt loop
+  } //end phi loop
+
+  //In-plane Out-of-plane anisotropy
+
+  //canvas
+  //mass in bins of pt
+  TCanvas* cvnewptbins=new TCanvas("cvnewptbins","Some bins in pt, in-plane/out-of-plane",1600,800);
+  cvnewptbins->Divide(nptbinsnew,2);
+  //chi square per ptbin
+  TCanvas* cvchinewptbins=new TCanvas("cvchinewptbins", "Chi vs pt");
+
+  //legend IN-PLANE OUT-OF-PLANE
+  TLegend* legsig=new TLegend(0.7,0.6,0.9,0.8,"");
+  legsig->SetBorderSize(0);
+  legsig->SetFillStyle(0);
+
+  TPaveText* txtoff=0x0;
+  if(official){
+    gStyle->SetOptTitle(0);
+    gStyle->SetFrameBorderMode(0);
+    txtoff=new TPaveText(0.1,0.53,0.55,0.83,"NDC");
+    txtoff->SetBorderSize(0);
+    txtoff->SetFillStyle(0);
+    txtoff->AddText(Form("D^{0}#rightarrow K^{-}#pi^{+} %.0f#times10^{6} events",hstat->Integral(0,1)/1e6));
+    txtoff->AddText("Pb-Pb @ #sqrt{s_{NN}} = 2.76TeV");
+    txtoff->AddText(Form("Centrality %d-%d%s",mincentr,maxcentr,"%"));
+    txtoff->SetTextColor(kCyan+3);
+  }
+
+  Double_t signal[2][nptbinsnew],background[2][nptbinsnew],significance[2][nptbinsnew];
+  Double_t errsgn[2][nptbinsnew],errbkg[2][nptbinsnew],errsignificance[2][nptbinsnew];
+  Int_t k=0;
+  //TString textinoutpl[2]={"0<#Delta#phi<#pi/4 && 3/4#pi<#Delta#phi<#pi","#pi/4<#Delta#phi<3/4#pi"};
+
+
+  //histograms
+  //chi square
+  TH1F** hchivsptnew=new TH1F*[2];
+
+
+  for(Int_t i=0;i<2;i++){ //in plane/out of plane
+    printf("Drawing %s\n",i==0 ? "in plane" : "out of plane");
+    hchivsptnew[i]=new TH1F(Form("hchivsptnew%d",i), "Chi square from fit;p_{t} (GeV/c);#chi^{2}",nptbinsnew,ptledges);
+    hchivsptnew[i]->SetMarkerColor(i+1);
+    hchivsptnew[i]->SetMarkerStyle(20);
+
+     for(Int_t j=0;j<nptbinsnew;j++){
+      TPaveText* pvbinphipt=new TPaveText(0.1,0.8,0.7,0.9,"NDC");
+      pvbinphipt->SetBorderSize(0);
+      pvbinphipt->SetFillStyle(0);
+      if(j!=nptbinsnew-1) pvbinphipt->AddText(Form("%.0f<p_{t}<%.0f GeV/c",ptbinsnew[j],ptbinsnew[j+1]));
+      else pvbinphipt->AddText(Form("p_{t}>%.0f GeV/c",ptbinsnew[j]));
+      k++;
+
+      if(!hmassnonintphi[i][j]) {
+       printf("hist %s pt %d not found\n",i==0 ? "In-plane" : "out-of-plane",j);
+       continue;
+      }
+
+      fout->cd();
+      hmassnonintphi[i][j]->Write();
+
+      cvnewptbins->cd(k);
+
+      //Fit mass histograms in bins of pt in-plane or out-of-plane
+      AliHFMassFitter fitter(hmassnonintphi[i][j],hmassnonintphi[i][j]->GetBinLowEdge(1),hmassnonintphi[i][j]->GetBinLowEdge(hmassnonintphi[i][j]->GetNbinsX()+1),rebin[j],fittype);
+
+      Bool_t ok=fitter.MassFitter(kFALSE);
+      if(ok) {
+       fitter.DrawHere(cvnewptbins->cd(k),nsigma,info);
+
+       fitter.Signal(nsigma,signal[i][j],errsgn[i][j]);
+       fitter.Background(nsigma,background[i][j],errbkg[i][j]);
+       fitter.Significance(nsigma,significance[i][j],errsignificance[i][j]);
+
+       //fill chi square
+       hchivsptnew[i]->SetBinContent(j,fitter.GetChiSquare());
+
+      }
+      else { //fit failed
+       hmassnonintphi[i][j]->Draw();
+       hchivsptnew[i]->SetBinContent(j,-1);
+      }
+
+      //draw mass
+      pvbinphipt->Draw();
+     }
+
+     //draw chi square
+     cvchinewptbins->cd();
+     if(i==0){
+       hchivsptnew[i]->Draw("ptext");
+       legsig->AddEntry(hchivsptnew[i],"In-plane","p");
+     }
+     else {
+       hchivsptnew[i]->Draw("psamestext");
+       legsig->AddEntry(hchivsptnew[i],"Out-on-plane","p");
+     }
+  }
+
+
+  //write png of mass as a function of pt in plane and out of plane and chi square
+  cvchinewptbins->cd();
+  legsig->Draw();
+  WriteCanvas(cvchinewptbins,textleg);
+  WriteCanvas(cvnewptbins,Form("inout%dptbinscentr%d_%d",nptbinsnew,mincentr,maxcentr));
+
+  //write canvas chi square
+  fout->cd();
+  cvchinewptbins->Write();
+
+  //graphs of yields
+  TGraphErrors** gr=new TGraphErrors*[2];
+  TGraphErrors** grb=new TGraphErrors*[2];
+  TGraphErrors** grsgf=new TGraphErrors*[2];
+  TGraph** grerr=new TGraph*[2];
+  TGraph** grberr=new TGraph*[2];
+  //canvas for these graphs
+  TCanvas* cvsig=new TCanvas("cvsig","Signal vs pt");
+  TCanvas* cvsgnf=new TCanvas("cvsgnf","Significance vs pt");
+  TCanvas* cvbkg=new TCanvas("cvbkg","Background vs pt");
+  TCanvas* cvsiger=new TCanvas("cvsiger","Error on Signal vs pt");
+  TCanvas* cvbkger=new TCanvas("cvbkger","Error on Background vs pt");
+
+  for(Int_t i=0;i<2;i++) { //in/out
+    gr[i]=new TGraphErrors(0);
+    grb[i]=new TGraphErrors(0);
+    grsgf[i]=new TGraphErrors(0);
+    grerr[i]=new TGraph(0);
+    grberr[i]=new TGraph(0);
+
+    for(Int_t j=0;j<nptbinsnew;j++){ //loop on ptbins
+
+      //signal
+      gr[i]->SetPoint(j,pt[j],signal[i][j]);
+      gr[i]->SetPointError(j,deltapt[j],errsgn[i][j]);
+      cout<<j<<" pt "<<pt[j]<<" pm "<<deltapt[j]<<", sgn "<<signal[i][j]<<" pm "<<errsgn[i][j]<<endl;
+
+      //background
+      grb[i]->SetPoint(j,pt[j],background[i][j]);
+      grb[i]->SetPointError(j,deltapt[j],errbkg[i][j]);
+      grsgf[i]->SetPoint(j,pt[j],significance[i][j]);
+      grsgf[i]->SetPointError(j,deltapt[j],errsignificance[i][j]);
+
+      //error on signal
+      grerr[i]->SetPoint(j,pt[j],errsgn[i][j]/signal[i][j]);
+      grberr[i]->SetPoint(j,pt[j],errbkg[i][j]/background[i][j]);
+    }
+
+    //titles etc.
+    gr[i]->SetTitle("Signal vs p_{t}");
+    gr[i]->SetName(Form("sig%s", i==0 ? "inplane" : "outofplane"));
+    gr[i]->GetXaxis()->SetTitle("p_{t}");
+    gr[i]->GetYaxis()->SetTitle("Signal");
+    gr[i]->SetMinimum(0);
+    gr[i]->SetMarkerStyle(20);
+    gr[i]->SetMarkerColor(i+1);
+    
+    grb[i]->SetTitle("Background vs p_{t}");
+    grb[i]->SetName(Form("bkg%s", i==0 ? "inplane" : "outofplane"));
+    grb[i]->GetXaxis()->SetTitle("p_{t} (GeV/c)");
+    grb[i]->GetYaxis()->SetTitle("Background");
+    grb[i]->SetMinimum(0);
+    grb[i]->SetMarkerStyle(20);
+    grb[i]->SetMarkerColor(i+1);
+
+    grsgf[i]->SetTitle("Significance vs p_{t};p_{t} (GeV/c);Significance");
+    grsgf[i]->SetName(Form("sgf%s", i==0 ? "inplane" : "outofplane"));
+    grsgf[i]->SetMinimum(0);
+    grsgf[i]->SetMarkerStyle(20);
+    grsgf[i]->SetMarkerColor(i+1);
+
+    grerr[i]->SetTitle("#epsilon_{r} Signal");
+    grerr[i]->SetName(Form("sigerr%s", i==0 ? "inplane" : "outofplane"));
+    grerr[i]->GetYaxis()->SetTitle("#epsilon_{r} Signal");
+    grerr[i]->GetXaxis()->SetTitle("p_{t}");
+    grerr[i]->SetMinimum(0);
+    grerr[i]->SetMarkerStyle(20);
+    grerr[i]->SetMarkerColor(i+1);
+
+    grberr[i]->SetTitle("#epsilon_{r} Background");
+    grberr[i]->SetName(Form("bkgerr%s", i==0 ? "inplane" : "outofplane"));
+    grberr[i]->GetYaxis()->SetTitle("#epsilon_{r} Background");
+    grberr[i]->GetXaxis()->SetTitle("p_{t}");
+    grberr[i]->SetMinimum(0);
+    grberr[i]->SetMarkerStyle(20);
+    grberr[i]->SetMarkerColor(i+1);
+
+    //draw signal, bkg, significance, errors
+   if(i==0) {
+      cvsig->cd();
+      gr[i]->Draw("AP");
+
+      cvbkg->cd();
+      grb[i]->Draw("AP");
+
+      cvsgnf->cd();
+      grsgf[i]->Draw("AP");
+
+      cvsiger->cd();
+      grerr[i]->Draw("AP");
+
+      cvbkger->cd();
+      grberr[i]->Draw("AP");
+    }
+    else {
+      cvsig->cd();
+      gr[i]->Draw("P");
+
+      cvbkg->cd();
+      grb[i]->Draw("P");
+
+      cvsgnf->cd();
+      grsgf[i]->Draw("P");
+
+      cvsiger->cd();
+      grerr[i]->Draw("P");
+
+      cvbkger->cd();
+      grberr[i]->Draw("P");
+   }
+
+   //write in output file
+   fout->cd();
+   gr[i]->Write();
+   grb[i]->Write();
+   grerr[i]->Write();
+   grsgf[i]->Write();
+  }
+
+  cvsig->cd();
+  legsig->Draw();
+  cvbkg->cd();
+  legsig->Draw();
+  cvsgnf->cd();
+  legsig->Draw();
+  cvsiger->cd();
+  legsig->Draw();
+  cvbkger->cd();
+  legsig->Draw();
+
+  //write png
+  WriteCanvas(cvsig,textleg);
+  WriteCanvas(cvbkg,textleg);
+  WriteCanvas(cvsiger,textleg);
+  WriteCanvas(cvbkger,textleg);
+  WriteCanvas(cvsgnf,textleg);
+
+  //fit again with sigma obtained from total inv mass (in-plane+out-of-plane)
+  Double_t sigmafullmass[nptbinsnew];
+
+  //canvas
+  TCanvas* cvmasstotnewptbins=new TCanvas("cvmasstotnewptbins", "Mass in some pt bins in+out",1600,400);
+  cvmasstotnewptbins->Divide(nptbinsnew,1);
+  TCanvas* cvnewptbins2=new TCanvas("cvnewptbins2", "Mass in some pt bins in/out of plane fixed sigma",1600,800);
+  cvnewptbins2->Divide(nptbinsnew,2);
+
+  for(Int_t i=0;i<nptbinsnew;i++){//pt bins
+
+    //pave text
+    TPaveText* pvbinpt=new TPaveText(0.1,0.8,0.7,0.9,"NDC");
+    pvbinpt->SetBorderSize(0);
+    pvbinpt->SetFillStyle(0);
+    if(i!= nptbinsnew-1) pvbinpt->AddText(Form("%.0f<p_{t}<%.0f GeV/c",ptbinsnew[i],ptbinsnew[i+1]));
+    else if(ptbinsnew[i+1]-ptbinsnew[i] > 10) pvbinpt->AddText(Form("p_{t}>%.0f GeV/c",ptbinsnew[i]));
+
+    //mass histogram integrated in phi
+    TH1F* hmassfull=(TH1F*)hmassnonintphi[0][i]->Clone();
+    hmassfull->Add(hmassnonintphi[1][i]);
+
+    AliHFMassFitter fitter(hmassfull,hmassfull->GetBinLowEdge(1),hmassfull->GetBinLowEdge(hmassfull->GetNbinsX()+1),rebin[i],fittype);
+    Bool_t ok=fitter.MassFitter(kFALSE);
+    if(ok){
+      fitter.DrawHere(cvmasstotnewptbins->cd(i+1),nsigma,info);
+      sigmafullmass[i]=fitter.GetSigma();
+    }else{
+      cvmasstotnewptbins->cd(i+1);
+      hmassfull->Draw();
+      sigmafullmass[i]=-2;
+    }
+    cvmasstotnewptbins->cd(i+1);
+    pvbinpt->Draw();
+
+
+    //in-plane
+    AliHFMassFitter fitter0(hmassnonintphi[0][i],hmassnonintphi[0][i]->GetBinLowEdge(1),hmassnonintphi[0][i]->GetBinLowEdge(hmassnonintphi[0][i]->GetNbinsX()+1),rebin[i],fittype);
+    if(sigmafullmass[i]!=-2) fitter0.SetFixGaussianSigma(sigmafullmass[i]);
+    ok=fitter0.MassFitter(kFALSE);
+    if(ok){
+      fitter0.DrawHere(cvnewptbins2->cd(i+1),nsigma,info);
+      fitter0.Signal(nsigma,signal[0][i],errsgn[0][i]);
+      fitter0.Background(nsigma,background[0][i],errbkg[0][i]);
+      fitter0.Significance(nsigma,significance[0][i],errsignificance[0][i]);
+
+    }else{
+      cvnewptbins2->cd(i+1);
+      hmassnonintphi[0][i]->Draw();
+    }
+    cvnewptbins2->cd(i+1);
+    pvbinpt->Draw();
+
+
+    if(official) {
+      cvnewptbins2->cd(i+1);
+      TPaveText* txtsignif=new TPaveText(0.1,0.1,0.6,0.2,"NDC");
+      txtsignif->SetBorderSize(0);
+      txtsignif->SetFillStyle(0);
+      txtsignif->SetTextColor(kGreen+3);
+      txtsignif->AddText(Form("Sgnf(%.0f#sigma) = %.1f#pm%.1f",nsigma,significance[0][i],errsignificance[0][i]));
+
+      txtsignif->Draw();
+      TH1F* htmp=(TH1F*)cvnewptbins2->cd(i+1)->FindObject("fhistoInvMass");
+      htmp->SetYTitle(Form("Entries/%.3f GeV/c^{2}",htmp->GetBinWidth(3)));
+      htmp->SetXTitle("D^{0} invariant mass (GeV/c^{2})");
+    }
+
+    if(official) InOutPic(cvnewptbins2->cd(i+1),i);
+
+    //out-of-plane
+    AliHFMassFitter fitter1(hmassnonintphi[1][i],hmassnonintphi[1][i]->GetBinLowEdge(1),hmassnonintphi[1][i]->GetBinLowEdge(hmassnonintphi[1][i]->GetNbinsX()+1),rebin[i],fittype);
+    if(sigmafullmass[i]!=-2) fitter1.SetFixGaussianSigma(sigmafullmass[i]);
+    ok=fitter1.MassFitter(kFALSE);
+    if(ok){
+      fitter1.DrawHere(cvnewptbins2->cd(nptbinsnew+i+1),nsigma,info);
+      fitter1.Signal(nsigma,signal[1][i],errsgn[1][i]);
+      fitter1.Background(nsigma,background[1][i],errbkg[1][i]);
+      fitter1.Significance(nsigma,significance[1][i],errsignificance[1][i]);
+    }else{
+      cvnewptbins2->cd(nptbinsnew+i+1);
+      hmassnonintphi[1][i]->Draw();
+    }
+    cvnewptbins2->cd(nptbinsnew+i+1);
+    pvbinpt->Draw();
+
+    if(official) {
+      cvnewptbins2->cd(nptbinsnew+i+1);
+      TPaveText* txtsignif=new TPaveText(0.1,0.1,0.6,0.2,"NDC");
+      txtsignif->SetBorderSize(0);
+      txtsignif->SetFillStyle(0);
+      txtsignif->SetTextColor(kGreen+3);
+      txtsignif->AddText(Form("Sgnf(%.0f#sigma) = %.1f#pm%.1f",nsigma,significance[1][i],errsignificance[1][i]));
+
+      txtsignif->Draw();
+      TH1F* htmp=(TH1F*)cvnewptbins2->cd(nptbinsnew+i+1)->FindObject("fhistoInvMass");
+      htmp->SetYTitle(Form("Entries/%.3f GeV/c^{2}",htmp->GetBinWidth(3)));
+      htmp->SetXTitle("D^{0} invariant mass (GeV/c^{2})");
+    }
+    if(official) InOutPic(cvnewptbins2->cd(nptbinsnew+i+1),i);
+
+  }
+
+  //other drawing details
+  if(official){
+    cvnewptbins2->cd(2);
+    txtoff->Draw();
+    gROOT->LoadMacro("/home/cbianchi/macros/ALICEPerformance.C");
+    gROOT->ProcessLine(Form("ALICEPerformance((TVirtualPad*)%p,\"today\",\"br\")",cvnewptbins2->cd(1)));
+  }
+
+
+  //png mass histograms
+  WriteCanvas(cvnewptbins2,Form("inout%dptbinscentr%d_%d",nptbinsnew,mincentr,maxcentr));
+  WriteCanvas(cvmasstotnewptbins,Form("totmass%dptbinscentr%d_%d",nptbinsnew,mincentr,maxcentr));
+  //graphs of yields
+  TGraphErrors** gr2=new TGraphErrors*[2];
+  TGraphErrors** grb2=new TGraphErrors*[2];
+  TGraphErrors** grsgf2=new TGraphErrors*[2];
+  TGraph** grerr2=new TGraph*[2];
+  TGraph** grberr2=new TGraph*[2];
+
+  //canvas for these graphs
+  TCanvas* cvsig2=new TCanvas("cvsig2","Signal vs pt (fixed sigma)");
+  TCanvas* cvbkg2=new TCanvas("cvbkg2","Background vs pt (fixed sigma)");
+  TCanvas* cvsgnf2=new TCanvas("cvsgnf2","Significance vs pt (fixed sigma)");
+  TCanvas* cvsiger2=new TCanvas("cvsiger2","Error on Signal vs pt (fixed sigma)");
+  TCanvas* cvbkger2=new TCanvas("cvbkger2","Error on Background vs pt (fixed sigma)");
+
+  for(Int_t i=0;i<2;i++) { //in/out
+
+    gr2[i]=new TGraphErrors(0);
+    grb2[i]=new TGraphErrors(0);
+    grsgf2[i]=new TGraphErrors(0);
+    grerr2[i]=new TGraph(0);
+    grberr2[i]=new TGraph(0);
+
+    for(Int_t j=0;j<nptbinsnew;j++){ //bins of pt
+
+      if(significance[i][j]>3.){
+       gr2[i]->SetPoint(j,pt[j],signal[i][j]);
+       gr2[i]->SetPointError(j,deltapt[j],errsgn[i][j]);
+       grsgf2[i]->SetPoint(j,pt[j],significance[i][j]);
+       grsgf2[i]->SetPointError(j,deltapt[j],errsignificance[i][j]);
+       grerr2[i]->SetPoint(j,pt[j],errsgn[i][j]/signal[i][j]);
+      }
+      cout<<j<<" pt "<<pt[j]<<" pm "<<deltapt[j]<<", sgn "<<signal[i][j]<<" pm "<<errsgn[i][j]<<endl;
+      grb2[i]->SetPoint(j,pt[j],background[i][j]);
+      grb2[i]->SetPointError(j,deltapt[j],errbkg[i][j]);
+      grberr2[i]->SetPoint(j,pt[j],errbkg[i][j]/background[i][j]);
+    }
+
+    gr2[i]->SetTitle("Signal vs p_{t} (#sigma fixed)");
+    gr2[i]->SetName(Form("sig%sfs", i==0 ? "inplane" : "outofplane"));
+    gr2[i]->GetXaxis()->SetTitle("p_{t} (GeV/c)");
+    gr2[i]->GetYaxis()->SetTitle("Signal");
+    gr2[i]->SetMinimum(0);
+    gr2[i]->SetMarkerStyle(20);
+    gr2[i]->SetMarkerColor(i+1);
+    
+    grb2[i]->SetTitle("Background vs p_{t} (#sigma fixed)");
+    grb2[i]->SetName(Form("bkg%sfs", i==0 ? "inplane" : "outofplane"));
+    grb2[i]->GetXaxis()->SetTitle("p_{t} (GeV/c)");
+    grb2[i]->GetYaxis()->SetTitle("Background");
+    grb2[i]->SetMinimum(0);
+    grb2[i]->SetMarkerStyle(20);
+    grb2[i]->SetMarkerColor(i+1);
+
+    grsgf2[i]->SetTitle("Significance vs p_{t} (#sigma fixed);p_{t} (GeV/c);Significance");
+    grsgf2[i]->SetName(Form("sgf%sfs", i==0 ? "inplane" : "outofplane"));
+    grsgf2[i]->SetMinimum(0);
+    grsgf2[i]->SetMarkerStyle(20);
+    grsgf2[i]->SetMarkerColor(i+1);
+
+    grerr2[i]->SetTitle("#epsilon_{r} Signal (#sigma fixed)");
+    grerr2[i]->SetName(Form("sgnerr%sfs", i==0 ? "inplane" : "outofplane"));
+    grerr2[i]->GetYaxis()->SetTitle("#epsilon_{r} Signal");
+    grerr2[i]->GetXaxis()->SetTitle("p_{t} (GeV/c)");
+    grerr2[i]->SetMinimum(0);
+    grerr2[i]->SetMarkerStyle(20);
+    grerr2[i]->SetMarkerColor(i+1);
+
+    grberr2[i]->SetTitle("#epsilon_{r} Background (#sigma fixed)");
+    grberr2[i]->SetName(Form("bkgerr%sfs", i==0 ? "inplane" : "outofplane"));
+    grberr2[i]->GetYaxis()->SetTitle("#epsilon_{r} Background");
+    grberr2[i]->GetXaxis()->SetTitle("p_{t} (GeV/c)");
+    grberr2[i]->SetMinimum(0);
+    grberr2[i]->SetMarkerStyle(20);
+    grberr2[i]->SetMarkerColor(i+1);
+
+   
+    if(i==0) {
+      cvsig2->cd();
+      gr2[i]->Draw("AP");
+
+      cvbkg2->cd();
+      grb2[i]->Draw("AP");
+
+      cvsgnf2->cd();
+      grsgf2[i]->Draw("AP");
+
+      cvsiger2->cd();
+      grerr2[i]->Draw("AP");
+
+      cvbkger2->cd();
+      grberr2[i]->Draw("AP");
+    }
+    else {
+      cvsig2->cd();
+      gr2[i]->Draw("P");
+
+      cvbkg2->cd();
+      grb2[i]->Draw("P");
+
+      cvsgnf2->cd();
+      grsgf2[i]->Draw("P");
+
+      cvsiger2->cd();
+      grerr2[i]->Draw("P");
+
+      cvbkger2->cd();
+      grberr2[i]->Draw("P");
+    }
+
+    //write in output file
+    fout->cd();
+    gr2[i]->Write();
+    grb2[i]->Write();
+    grerr2[i]->Write();
+    grsgf2[i]->Write();
+  }
+
+  cvsig2->cd();
+  legsig->Draw();
+  cvbkg2->cd();
+  legsig->Draw();
+  cvsgnf2->cd();
+  legsig->Draw();
+  cvsiger2->cd();
+  legsig->Draw();
+  cvbkger2->cd();
+  legsig->Draw();
+
+  //write png
+  WriteCanvas(cvsig2,textleg);
+  WriteCanvas(cvbkg2,textleg);
+  WriteCanvas(cvsgnf2,textleg);
+  WriteCanvas(cvsiger2,textleg);
+  WriteCanvas(cvbkger2,textleg);
+
+
+  //invariant mass in bins of phi, pt integrated
+
+  Double_t sigvsphiptint[nphibins],errsigvsphiptint[nphibins],sgnfvsphiptint[nphibins],errsgnfvsphiptint[nphibins];
+
+  //canvas
+  TCanvas* cvptint=new TCanvas("cvptint","Pt integrated mass plots",1600,600);
+  cvptint->Divide(nphibins,1);
+  TCanvas* cvphierr=new TCanvas("cvphierr","Error on Signal phi bins");
+  TCanvas* cvchiphi=new TCanvas("cvchiphi","Chi as a function of phi");
+  TCanvas* cvsgnvsphiptint=new TCanvas("cvsgnvsphiptint","Sigmal vs #Delta #phi");
+
+  //chi square and yields
+  TH1F* hchivsphi=(TH1F*)hnphibins->Clone("hchivsphi");
+  hchivsphi->SetTitle("#chi^{2} vs #phi;#phi (rad);#chi^{2}");
+  hchivsphi->SetMarkerStyle(20);
+  TGraphErrors* grphi=new TGraphErrors(0);
+  TGraph* grphierr=new TGraph(0);
+
+  for (Int_t i=0;i<nphibins;i++){ //bins of phi
+    cvptint->cd(i+1);
+
+    if(hmassptint[i]) {
+      AliHFMassFitter fitter(hmassptint[i],hmassptint[i]->GetBinLowEdge(1),hmassptint[i]->GetBinLowEdge(hmassptint[i]->GetNbinsX()+1),rebin[0],fittype);
+      Bool_t ok=fitter.MassFitter(kFALSE);
+
+      TPaveText* txtsignif=new TPaveText(0.1,0.1,0.6,0.2,"NDC");
+      txtsignif->SetBorderSize(0);
+      txtsignif->SetFillStyle(0);
+      txtsignif->SetTextColor(kGreen+3);
+
+      if(ok) {
+       fitter.DrawHere(cvptint->cd(i+1),nsigma,info);
+       fitter.Signal(nsigma,sigvsphiptint[i],errsigvsphiptint[i]);
+
+       fitter.Significance(nsigma,sgnfvsphiptint[i],errsgnfvsphiptint[i]);
+       txtsignif->AddText(Form("Sgnf(%.0f#sigma) = %.1f#pm%.1f",nsigma,sgnfvsphiptint[i],errsgnfvsphiptint[i]));
+
+       //fill chi square
+       hchivsphi->SetBinContent(i,fitter.GetChiSquare());
+       //fill yields vs phi
+       grphi->SetPoint(i,phibins[i],sigvsphiptint[i]);
+       grphi->SetPointError(i,dx[i],errsigvsphiptint[i]);
+       grphierr->SetPoint(i, phibins[i],errsigvsphiptint[i]/sigvsphiptint[i]);
+      }
+      else {
+       //retry fit if it fails
+       fitter.Reset();
+       fitter.SetHisto(hmassptint[i]);
+       fitter.SetRangeFit(hmassptint[i]->GetBinLowEdge(3),hmassptint[i]->GetBinLowEdge(hmassptint[i]->GetNbinsX()-2));
+       fitter.RebinMass(rebin[0]*2);
+       fitter.SetType(fittype,0);
+       ok=fitter.MassFitter(kFALSE);
+       if(ok){
+         fitter.DrawHere(cvptint->cd(i+1),nsigma,info);
+         TPaveText* pvnote=new TPaveText(0.1,0.7,0.8,0.9,"NDC");
+         pvnote->SetBorderSize(0);
+         pvnote->SetFillStyle(0);
+         pvnote->AddText("Refitted w/ different range and rebin");
+         cvptint->cd(i+1);
+         pvnote->Draw();
+
+         fitter.Signal(nsigma,sigvsphiptint[i],errsigvsphiptint[i]);
+
+         //fill chi square
+         hchivsphi->SetBinContent(i,fitter.GetChiSquare());
+         //fill yields vs phi
+         grphi->SetPoint(i,phibins[i],sigvsphiptint[i]);
+         grphi->SetPointError(i,dx[i],errsigvsphiptint[i]);
+         grphierr->SetPoint(i, phibins[i],errsigvsphiptint[i]/sigvsphiptint[i]);
+       } else hmassptint[i]->Draw(); //fit didn't work at all
+      }
+      cvptint->cd(i+1);
+
+      if(official){
+       txtsignif->Draw();
+       //PhiBinPic(cvptint->cd(i+1),i,"tr");
+      }
+
+    }
+
+  }
+
+  cvchiphi->cd();
+  hchivsphi->Draw("ptext");
+
+  grphi->SetTitle("Signal vs #phi");
+  grphi->SetName("gsigvsphi");
+  grphi->GetXaxis()->SetTitle("#phi (rad)");
+  grphi->GetYaxis()->SetTitle("Signal");
+  grphi->SetMinimum(0);
+  grphi->SetMarkerStyle(20);
+  cvsgnvsphiptint->cd();
+  grphi->Draw("AP");
+
+  grphierr->SetTitle("#epsilon_{r} Signal vs #phi");
+  grphierr->SetName("gsigerrvsphi");
+  grphierr->GetXaxis()->SetTitle("#phi (rad)");
+  grphierr->GetYaxis()->SetTitle("#epsilon_{r} Signal");
+  grphierr->SetMinimum(0);
+  grphierr->SetMarkerStyle(20);
+  cvphierr->cd();
+  grphierr->Draw("AP");
+
+  //fit dN/dDeltaphi
+  TF1 *flowFunc = new TF1("flow","[0]*(1.+2.*[1]*TMath::Cos(2.*x))");
+  grphi->Fit(flowFunc);
+
+  TPaveText* pvv2=new TPaveText(0.1,0.1,0.6,0.2,"NDC");
+  pvv2->SetBorderSize(0);
+  pvv2->SetFillStyle(0);
+  pvv2->AddText(Form("v'_{2} = %.3f#pm%.3f",flowFunc->GetParameter(1),flowFunc->GetParError(1)));
+
+  cvsgnvsphiptint->cd();
+  pvv2->Draw();
+
+  //write chi square, yields and errors in output file
+  fout->cd();
+  cvchiphi->Write();
+  grphierr->Write();
+  grphi->Write();
+
+  //other drawing details
+  if(official){
+    cvptint->cd(1);
+    txtoff->Draw();
+    gROOT->LoadMacro("/home/cbianchi/macros/ALICEPerformance.C");
+    gROOT->ProcessLine(Form("ALICEPerformance((TVirtualPad*)%p,\"today\",\"br\")",cvptint->cd(1)));
+  }
+
+  //write png mass histos yields and error
+  WriteCanvas(cvptint,textleg);
+  WriteCanvas(cvphierr,textleg);
+  WriteCanvas(cvsgnvsphiptint,textleg);
+
+  //statistics
+  TCanvas* cst=new TCanvas("cst","Stat");
+  cst->SetGridy();
+  cst->cd();
+  hstat->Draw("htext0");
+  WriteCanvas(cst,textleg);
+
+
+  //save event plane info in the output file
+  //fix to the latest version of the task!!!!!
+  TH1F* hresos=(TH1F*)list->FindObject(Form("hEvPlaneResocentr%d_%d",mincentr,mincentr+5));
+  /*
+  TH2F* hMphis=(TH2F*)list->FindObject(Form("hMphicentr%d_%d",mincentr,mincentr+5));
+  TH2F* hMc2phis=(TH2F*)list->FindObject(Form("hMc2phicentr%d_%d",mincentr,mincentr+5));
+  */
+  TString buildname=Form("centr%d_%d",mincentr,maxcentr);
+  hresos->SetName(Form("hEvPlaneReso%s",buildname.Data()));
+  /*
+  hMphis->SetName(Form("hMphi%s",buildname.Data()));
+  hMc2phis->SetName(Form("hMc2phi%s",buildname.Data()));
+  */
+  for(Int_t icentr=mincentr+5;icentr<=maxcentr;icentr=icentr+5){
+    TH1F* hreso=(TH1F*)list->FindObject(Form("hEvPlaneResocentr%d_%d",icentr,icentr+5));
+    hresos->Add(hreso);
+    // for(Int_t i=0;i<nptbins;i++){
+    // TH2F* hMphi=(TH2F*)list->FindObject(Form("hMphicentr_pt%d%d_%d",icentr,icentr+5));
+    // TH2F* hMc2phi=(TH2F*)list->FindObject(Form("hMc2phi_pt%dcentr%d_%d",icentr,icentr+5));
+    // hMphis->Add(hMphi);
+    // hMc2phis->Add(hMc2phi);
+    // }
+  }
+
+  fout->cd();
+  hresos->Write();
+  /*
+  hMphis->Write();
+  hMc2phis->Write();
+  */
+}
+
+void Dmesonsv2InOutAnisotropy(Int_t mincentr,Int_t maxcentr,Bool_t fixedsigma=kTRUE){
+
+  //Compute v2 with the method of the anisotropy between in-plane and out-of plane. Needs the output file from the previous method
+  //Author: Francesco Prino prino@to.infn.it
+
+  gStyle->SetCanvasColor(0);
+  gStyle->SetTitleFillColor(0);
+  gStyle->SetStatColor(0);
+  gStyle->SetOptTitle(0);
+
+  TString name=Form("centr%d-%d",mincentr,maxcentr),
+    nameh=Form("centr%d_%d",mincentr,maxcentr);
+
+  TFile* fil=new TFile(Form("HistoInputv2Calc%s.root",name.Data()));
+  if(!fil->IsOpen()){
+    cout<<"Input file not found"<<endl;
+    return;
+  }
+  
+  //histogram resolution with subevents
+  TH1F* hResolSubAB=(TH1F*)fil->Get(Form("hEvPlaneReso%s",nameh.Data()));
+  hResolSubAB->GetXaxis()->SetTitle("cos[2(#psi_{A}-#psi_{B})]");
+
+  //graphs yield
+  TGraphErrors* gsigin=(TGraphErrors*)fil->Get(Form("siginplane%s",fixedsigma ? "fs" : ""));
+  TGraphErrors* gsigout=(TGraphErrors*)fil->Get(Form("sigoutofplane%s",fixedsigma ? "fs" : ""));
+
+  TCanvas* c1=new TCanvas("c1","EvPlaneResol");
+  hResolSubAB->Draw();
+  Double_t resolSub=TMath::Sqrt(hResolSubAB->GetMean());
+  printf("Resolution on sub event  = %.4f\n",resolSub);
+  Double_t chisub=FindChi(resolSub,1);
+  printf("Chi Subevent             = %.4f\n",chisub);
+  Double_t chifull=chisub*TMath::Sqrt(2);
+  printf("Chi Full Event           = %.4f\n",chifull);
+  Double_t resolFull=ComputeResol(resolSub,1);
+  printf("Resolution on full event = %.4f\n",resolFull);
+  TLatex* tres=new TLatex(0.15,0.7,Form("Resolution on full event = %.4f\n",resolFull));
+  tres->SetNDC();
+  tres->Draw();
+  c1->SaveAs(Form("EvPlaneResol%s.png",nameh.Data()));
+
+  TCanvas* c2=new TCanvas("c2",Form("Yields %s",fixedsigma ? "fixed sigma" : ""));
+  c2->cd();
+  gsigin->Draw("AP");
+  gsigout->Draw("Psame");
+
+  Int_t nPtBins=gsigin->GetN();
+
+  //graphs results
+  TGraphErrors *gAnis=new TGraphErrors(0);
+  gAnis->SetName(Form("gAnisotropy%s",nameh.Data()));
+  gAnis->SetTitle("(Nin-Nout)/(Nin+Nout)");
+
+  TGraphErrors *gv2obs=new TGraphErrors(0);
+  gv2obs->SetName(Form("gv2obs%s",nameh.Data()));
+  gv2obs->SetTitle("v2obs");
+
+  TGraphErrors *gv2=new TGraphErrors(0);
+  gv2->SetName(Form("gv2%s",nameh.Data()));
+  gv2->SetTitle("v2");
+
+  TGraph *gv2err=new TGraph(0);
+  gv2err->SetName(Form("gv2err%s",nameh.Data()));
+
+  for(Int_t iBin=0; iBin<nPtBins; iBin++){
+    Double_t pt,nIn,ept,enIn;
+    Double_t nOut,enOut;
+    gsigin->GetPoint(iBin,pt,nIn);
+    ept=gsigin->GetErrorX(iBin);
+    enIn=gsigin->GetErrorY(iBin);
+    gsigout->GetPoint(iBin,pt,nOut);
+    enOut=gsigout->GetErrorY(iBin);
+    if(nIn<1e-6 || nOut<1e-6) continue;
+
+    Double_t anis=(nIn-nOut)/(nIn+nOut);
+    Double_t eAnis=TMath::Sqrt(enIn*enIn+enOut*enOut)/(nIn+nOut);
+    gAnis->SetPoint(iBin,pt,anis);
+    gAnis->SetPointError(iBin,ept,eAnis);
+
+    Double_t v2obs=anis*TMath::Pi()/4.;
+    Double_t ev2obs=eAnis*TMath::Pi()/4.;
+    gv2obs->SetPoint(iBin,pt,v2obs);
+    gv2obs->SetPointError(iBin,ept,ev2obs);
+
+    Double_t v2=v2obs/resolFull;
+    Double_t ev2=ev2obs/resolFull;
+    gv2->SetPoint(iBin,pt,v2);
+    gv2->SetPointError(iBin,ept,ev2);
+
+    gv2err->SetPoint(iBin,pt,ev2/TMath::Abs(v2));
+  }
+
+  TCanvas* c3=new TCanvas("c3","Anisotropy");
+  gAnis->SetMarkerStyle(20);
+  gAnis->GetXaxis()->SetTitle("p_{t} (GeV/c)");
+  gAnis->GetYaxis()->SetTitle("(N_{IN-PLANE}-N_{OUT-OF-PLANE})/(N_{IN-PLANE}+N_{OUT-OF-PLANE})");
+  gAnis->GetYaxis()->SetTitleOffset(1.1);
+  gAnis->Draw("AP");
+  gAnis->Draw("Psame");
+  c3->SaveAs(Form("Anisotropy%s%s.png",name.Data(),fixedsigma ? "fs" : ""));
+
+  TCanvas* c4=new TCanvas("c4","v2obs");
+  gv2obs->SetMarkerStyle(20);
+  gv2obs->GetXaxis()->SetTitle("p_{t} (GeV/c)");
+  gv2obs->GetYaxis()->SetTitle("v_{2}^{obs}");
+  gv2obs->GetYaxis()->SetTitleOffset(1.1);
+  gv2obs->Draw("AP");
+  gv2obs->Draw("Psame");
+  c4->SaveAs(Form("v2obsD%s%s.png",name.Data(),fixedsigma ? "fs" : ""));
+
+  TCanvas* c5=new TCanvas("c5","v2");
+  gv2->SetMarkerStyle(20);
+  gv2->GetXaxis()->SetTitle("p_{t} (GeV/c)");
+  gv2->GetYaxis()->SetTitle("v_{2}");
+  gv2->GetYaxis()->SetTitleOffset(1.1);
+  gv2->Draw("AP");
+  gv2->Draw("Psame");
+  c5->SaveAs(Form("v2D%s%s.png",name.Data(),fixedsigma ? "fs" : ""));
+
+  TCanvas* c6=new TCanvas("c6","error v2");
+  gv2err->SetTitle("Relative error v2;p_{t} (GeV/c);#epsilon_{r} v2");
+  gv2err->SetMarkerStyle(20);
+  gv2err->GetYaxis()->SetTitleOffset(1.1);
+  c6->cd();
+  gv2err->Draw("AP");
+  c6->SaveAs(Form("Errv2%s%s.png",name.Data(),fixedsigma ? "fs" : ""));
+
+  TFile* fout=new TFile(Form("Outputv2%s.root",name.Data()),"recreate");
+
+  fout->cd();
+  gAnis->Write();
+  gv2obs->Write();
+  gv2->Write();
+  gv2err->Write();
+  c1->Write(); //canvas resolution
+}
+
+//methods needed inside Dmesonsv2InOutAnisotropy
+//Author: Francesco Prino prino@to.infn.it
+
+Double_t ComputeResol(Double_t resSub, Int_t k){
+  Double_t chisub=FindChi(resSub,k);
+  Double_t chifull=chisub*TMath::Sqrt(2);
+  if(k==1) return ResolK1(chifull);
+  else if(k==2) return Pol(chifull,2);
+  else{
+    printf("k should be <=2\n");
+    return 1.;
+  }
+}
+
+Double_t FindChi(Double_t res, Int_t k){
+  Double_t x1=0;
+  Double_t x2=15;
+  Double_t y1,y2;
+  if(k==1){
+    y1=ResolK1(x1)-res;
+    y2=ResolK1(x2)-res;
+  }
+  else if(k==2){
+    y1=Pol(x1,2)-res;
+    y2=Pol(x2,2)-res;
+  }
+  else return -1;
+
+  if(y1*y2>0.) return -1;
+  if(y1==0.) return y1;
+  if(y2==0.) return y2;
+  Double_t xmed,ymed;
+  Int_t jiter=0;
+  while((x2-x1)>0.0001){
+    xmed=0.5*(x1+x2);
+    if(k==1){
+      y1=ResolK1(x1)-res;
+      y2=ResolK1(x2)-res;
+      ymed=ResolK1(xmed)-res;
+    }
+    else if(k==2){
+      y1=Pol(x1,2)-res;
+      y2=Pol(x2,2)-res;
+      ymed=Pol(xmed,2)-res;
+    }
+    else return -1;
+    if((y1*ymed)<0) x2=xmed;
+    if((y2*ymed)<0)x1=xmed;
+    if(ymed==0) return xmed;
+    jiter++;
+  }
+  return 0.5*(x1+x2);
+}
+
+Double_t Pol(Double_t x, Int_t k){
+  Double_t c[5];
+  if(k==1){ 
+    c[0]=0.626657; c[1]=0.; c[2]=-0.09694; c[3]=0.02754; c[4]=-0.002283;
+  }
+  else if(k==2){
+    c[0]=0.; c[1]=0.25; c[2]=-0.011414; c[3]=-0.034726; c[4]=0.006815;
+  }
+  else return -1;
+  return c[0]*x+c[1]*x*x+c[2]*x*x*x+c[3]*x*x*x*x+c[4]*x*x*x*x*x;
+}
+
+Double_t ResolK1(Double_t x){
+  return TMath::Sqrt(TMath::Pi()/8)*x*TMath::Exp(-x*x/4)*(TMath::BesselI0(x*x/4)+TMath::BesselI1(x*x/4));
+}
+
+
+void DrawEventPlane(Int_t mincentr=0,Int_t maxcentr=0,TString partname="D0",TString path="./"){
+
+ //draw the histograms correlated to the event plane
+
+  gStyle->SetCanvasColor(0);
+  gStyle->SetTitleFillColor(0);
+  gStyle->SetStatColor(0);
+  //gStyle->SetPalette(1);
+  gStyle->SetOptFit(1);
+
+  TString listname="coutputv2";
+  TList* list;
+  TH1F * hstat;
+  AliRDHFCuts* cutobj;
+
+  Bool_t isRead=ReadFile(list,hstat,cutobj,listname,partname,path);
+  if(!isRead) return;
+  if(!list || !hstat || !cutobj){
+    cout<<":-( null pointers..."<<endl;
+    return;
+  }
+  if(!(mincentr==0 && maxcentr==0)){ //draw the total in mincentr-maxcentr
+    TString suffixcentr=Form("centr%d_%d",mincentr,maxcentr);
+    TH1F* hevpls=(TH1F*)list->FindObject(Form("hEvPlanecentr%d_%d",mincentr,mincentr+5));
+    hevpls->SetName(Form("hEvPlane%s",suffixcentr.Data()));
+    hevpls->SetTitle(Form("Event Plane angle %s",suffixcentr.Data()));
+    TH1F* hevplresos=(TH1F*)list->FindObject(Form("hEvPlaneResocentr%d_%d",mincentr,mincentr+5));
+    hevplresos->SetName(Form("hEvPlaneReso%s",suffixcentr.Data()));
+    hevplresos->SetTitle(Form("Event Plane Resolution %s",suffixcentr.Data()));
+
+    for(Int_t icentr=mincentr+5;icentr<=maxcentr;icentr=icentr+5){
+      TH1F* h=(TH1F*)list->FindObject(Form("hEvPlanecentr%d_%d",icentr,icentr+5));
+      if(h)hevpls->Add(h);
+      else cout<<"skipping ev plane "<<icentr<<"_"<<icentr+5<<endl;
+      h=(TH1F*)list->FindObject(Form("hEvPlaneResocentr%d_%d",icentr,icentr+5));
+      if(h)hevplresos->Add(h);
+      else cout<<"skipping ev pl reso "<<icentr<<"_"<<icentr+5<<endl;
+    }
+
+    TCanvas* cvtotevpl=new TCanvas("cvtotevpl",Form("Ev plane %s",suffixcentr.Data()));
+    cvtotevpl->cd();
+    hevpls->Draw();
+    hevpls->Fit("pol0");
+
+    TCanvas* cvtotevplreso=new TCanvas("cvtotevplreso",Form("Ev plane Resolution %s",suffixcentr.Data()));
+    cvtotevplreso->cd();
+    hevplresos->Draw();
+    Double_t resolSub=TMath::Sqrt(hevplresos->GetMean());
+    Double_t resolFull=ComputeResol(resolSub,1);
+
+    TPaveText* pvreso=new TPaveText(0.1,0.1,0.6,0.2,"NDC");
+    pvreso->SetBorderSize(0);
+    pvreso->SetFillStyle(0);
+    pvreso->AddText(Form("Resolution on full event = %.4f\n",resolFull));
+    pvreso->Draw();
+
+    TFile* fout=new TFile(Form("EvPlanecentr%d-%d.root",mincentr,maxcentr),"recreate");
+    fout->cd();
+    hevpls->Write();
+    hevplresos->Write();
+  }
+  else{ //draw all in 5% centrality bins
+    for(Int_t i=0;i<list->GetEntries();i++){
+      TH1F* h=(TH1F*)list->At(i);
+      if(!h){
+       cout<<"Histogram "<<i<<" not found"<<endl;
+       continue;
+      }
+      TString hname=h->GetName();
+      if(hname.Contains("EvPlane")){
+       TCanvas* cv=new TCanvas(Form("cv%s",hname.Data()),hname.Data());
+       cv->cd();
+       h->Draw();
+
+       if(!hname.Contains("Reso")){
+         h->Fit("pol0");
+         
+       }else{
+         Double_t resolSub=TMath::Sqrt(h->GetMean());
+         Double_t resolFull=ComputeResol(resolSub,1);
+
+         TPaveText* pvreso=new TPaveText(0.1,0.1,0.6,0.2,"NDC");
+         pvreso->SetBorderSize(0);
+         pvreso->SetFillStyle(0);
+         pvreso->AddText(Form("Resolution on full event = %.4f\n",resolFull));
+         pvreso->Draw();
+
+       }
+      }
+    }
+  }
+
+}
+
+
+void DmesonsFlowAnalysisvsCentrality(Int_t mincentr, Int_t maxcentr, Int_t step=10){
+  //needs outputs from the other methods HistoInputv2Calccentr%d-%d.root,EvPlanecentr%d-%d.root, Outputv2centr%d-%d.root
+
+  //min(max)centr= min(max) value of centrality among those files - note that if there are "holes" they are skipped
+  //step= step in centrality between mincentr and maxcentr - default 10%
+
+  gStyle->SetCanvasColor(0);
+  gStyle->SetTitleFillColor(0);
+  gStyle->SetStatColor(0);
+  gStyle->SetOptStat(0);
+  gStyle->SetOptTitle(0);
+
+  //number of files expected
+  Int_t n=(maxcentr-mincentr)/step;
+  cout<<"Trying to read "<<n<<" files"<<endl;
+  
+  TString filename="",gname="";
+  Int_t icentr=0;
+
+  //v2
+  TCanvas* cvv2=new TCanvas("cvv2","v_{2}");
+  TCanvas* cvv2obs=new TCanvas("cvv2obs","v_{2}'");
+
+  //Event plane distribution and resolution
+  TCanvas* cvevplres=new TCanvas("cvevplres","Resolution");
+  cvevplres->SetLogy();
+  TH1F* hresolFull=new TH1F("hresolFull","Resolution vs centrality;centrality;resolution",n,mincentr,maxcentr);
+  TCanvas* cvres=new TCanvas("cvres", "Resolution vs centrality");
+  TCanvas* cvevplane=new TCanvas("cvevplane", "Event plane vs centrality");
+
+  //v2 in phi bins 
+  TCanvas* cvfitvsphi=new TCanvas("cvfitvsphi","Fit flow vs phi");
+  TH1F* hv2fcorr=new TH1F("hv2fcorr","v2 from fit corrected;centrality;v_{2}",n,mincentr,maxcentr);
+  hv2fcorr->SetMarkerStyle(30);
+  hv2fcorr->SetMarkerSize(2);
+  TH1F* hv2corr=new TH1F("hv2corr","v2 from asym corrected;centrality;v_{2}",n,mincentr,maxcentr);
+  hv2corr->SetMarkerStyle(29);
+  hv2corr->SetMarkerSize(1.7);
+  hv2corr->SetMarkerColor(2);
+  hv2corr->SetLineColor(2);
+  TCanvas* cvv2cor=new TCanvas("v2corr", "v2 from asymmetry corrected vs centrality");
+
+  TLegend* leg=new TLegend(0.7,0.6,0.9,0.8,"Centrality %");
+  leg->SetBorderSize(0);
+  leg->SetFillStyle(0);
+  TLegend* leg2=new TLegend(0.5,0.6,0.9,0.9,"");
+  leg2->SetBorderSize(0);
+  leg2->SetFillStyle(0);
+  TLegend* legresoval=new TLegend(0.2,0.6,0.5,0.9,"Full resolution");
+  legresoval->SetBorderSize(0);
+  legresoval->SetFillStyle(0);
+
+  Int_t nread=0;
+  for(Int_t i=0;i<n;i++){
+    icentr=mincentr+i*step;
+    filename=Form("Outputv2centr%d-%d.root",icentr,icentr+step);
+    TFile* f=new TFile(filename);
+    if(!f->IsOpen()){
+      cout<<filename.Data()<<" not found"<<endl;
+      continue;
+    }
+    nread++;
+    gname=Form("gv2centr%d_%d",icentr,icentr+step);
+    TGraphErrors* grv2=(TGraphErrors*)f->Get(gname);
+    if(!grv2){
+      cout<<gname<<" not found in "<<filename<<endl;
+      continue;
+    }
+    grv2->GetYaxis()->SetRangeUser(-0.2,0.4);
+    grv2->SetMarkerColor(i+1);
+    grv2->SetLineColor(i+1);
+    cvv2->cd();
+    if(nread==1)grv2->Draw("APL");
+    else grv2->Draw("PL");
+
+    leg->AddEntry(grv2,Form("%d-%d",icentr,icentr+step),"p");
+
+    gname=Form("gv2obscentr%d_%d",icentr,icentr+step);
+    TGraphErrors* grv2obs=(TGraphErrors*)f->Get(gname);
+    if(!grv2obs){
+      cout<<gname<<" not found in "<<filename<<endl;
+      continue;
+    }
+    grv2obs->GetYaxis()->SetRangeUser(-0.2,0.4);
+    grv2obs->SetMarkerColor(i+1);
+    grv2obs->SetLineColor(i+1);
+    cvv2obs->cd();
+    if(nread==1)grv2obs->Draw("APL");
+    else grv2obs->Draw("PL");
+
+    //draw ev plane reso
+    filename=Form("EvPlanecentr%d-%d.root",icentr,icentr+step);
+    TFile* f2=new TFile(filename);
+    if(!f2->IsOpen()){
+      cout<<filename.Data()<<" not found"<<endl;
+      continue;
+    }
+    TString hname=Form("hEvPlaneResocentr%d_%d",icentr,icentr+step);
+    TH1F* hevplres=(TH1F*)f2->Get(hname);
+    if(!hevplres){
+      cout<<hname<<" not found in "<<filename<<endl;
+      continue;
+    }
+    cvevplres->cd();
+    hevplres->SetLineColor(i+1);
+    if(nread==1) hevplres->Draw();
+    else hevplres->Draw("sames");
+
+    Double_t resolSub=TMath::Sqrt(hevplres->GetMean());
+    Double_t resolFull=ComputeResol(resolSub,1);
+    hresolFull->SetBinContent(i+1,resolFull);
+    TLegendEntry* lentry=legresoval->AddEntry(hevplres,Form("%.2f",resolFull),"");
+    lentry->SetTextColor(i+1);
+    cvevplres->cd();
+    legresoval->Draw();
+
+    hname=Form("hEvPlanecentr%d_%d",icentr,icentr+step);
+    TH1F* hevpl=(TH1F*)f2->Get(hname);
+    if(!hevpl){
+      cout<<hname<<" not found in "<<filename<<endl;
+      continue;
+    }
+    hevpl->Scale(1./hevpl->Integral());
+    cvevplane->cd();
+    hevpl->SetLineColor(i+1);
+    if(nread==1) hevpl->Draw();
+    else hevpl->Draw("sames");
+
+    filename=Form("HistoInputv2Calccentr%d-%d.root",icentr,icentr+step);
+    TFile* f3=new TFile(filename);
+    if(!f3->IsOpen()){
+      cout<<filename.Data()<<" not found"<<endl;
+      continue;
+    }
+    gname="gsigvsphi";
+    TGraphErrors* grdNdphi=(TGraphErrors*)f3->Get(gname);
+    if(!grdNdphi){
+      cout<<gname<<" not found in "<<filename<<endl;
+      continue;
+    }
+    TF1* ffit=(TF1*)grdNdphi->FindObject("flow");
+    ffit->SetName(Form("flow%d_%d",icentr,icentr+step));
+    ffit->SetLineColor(i+1);
+    // Double_t xmin,xmax;
+    // ffit->GetRange(xmin,xmax);
+    //ffit->Scale(1./ffit->Integral(xmin,xmax));
+    cvfitvsphi->cd();
+    if(nread==1){
+      ffit->SetMinimum(0);
+      ffit->Draw();
+    }
+    else ffit->Draw("sames");
+
+    Double_t v2obsf=ffit->GetParameter(1);
+    Double_t v2obsferr=ffit->GetParError(1);
+    Double_t v2fcorr=v2obsf/resolFull;
+    Double_t v2fcorrerr=v2obsferr/resolFull;
+    hv2fcorr->SetBinContent(i+1,v2fcorr);
+    hv2fcorr->SetBinError(i+1,v2fcorrerr);
+
+    gname="siginplane";
+    TGraphErrors* gsigin=(TGraphErrors*)f3->Get(gname);
+    gname="sigoutofplane";
+    TGraphErrors* gsigout=(TGraphErrors*)f3->Get(gname);
+    if(!gsigin || !gsigout){
+      cout<<"graphs in/out not found"<<endl;
+      continue;
+    }
+    Int_t nbins=((TParameter<float>*)f3->Get("nptbins"))->GetVal();
+    Double_t countin=0,countout=0,e2in=0,e2out=0;
+    for(Int_t k=0;k<nbins;k++){
+      Double_t x,y,ey;
+      gsigin->GetPoint(k,x,y);
+      ey=gsigin->GetErrorY(k);
+      e2in+=ey*ey;
+      countin+=y;
+      gsigout->GetPoint(k,x,y);
+      ey=gsigin->GetErrorY(k);
+      e2out+=ey*ey;
+      countout+=y;
+    }
+    if(countin<1e-6 || countout<1e-6) continue;
+    Double_t ein=TMath::Sqrt(e2in),eout=TMath::Sqrt(e2out);
+    Double_t anis=(countin-countout)/(countin+countout);
+    Double_t eAnis=TMath::Sqrt(ein*ein+eout*eout)/(countin+countout);
+    Double_t v2obs=anis*TMath::Pi()/4.;
+    Double_t ev2obs=eAnis*TMath::Pi()/4.;
+    Double_t v2corr=v2obs/resolFull;
+    Double_t v2correrr=ev2obs/resolFull;
+    hv2corr->SetBinContent(i+1,v2corr);
+    hv2corr->SetBinError(i+1,v2correrr);
+
+  }
+
+  cvv2->cd();
+  leg->Draw();
+  cvv2obs->cd();
+  leg->Draw();
+  cvevplres->cd();
+  leg->Draw();
+  cvevplane->cd();
+  leg->Draw();
+
+  cvres->cd();
+  hresolFull->Draw();
+
+  cvfitvsphi->cd();
+  leg->Draw();
+
+  cvv2cor->cd();
+  hv2corr->SetLineWidth(2);
+  hv2fcorr->SetLineWidth(2);
+  hv2corr->GetYaxis()->SetRangeUser(-0.2,0.4);
+  hv2corr->Draw();
+  hv2fcorr->Draw("sames");
+  leg2->AddEntry(hv2fcorr,"v_{2} from dN/d#phi","l");
+  leg2->AddEntry(hv2corr,"v_{2} from asymmetry","l");
+  leg2->Draw();
+
+  cvv2->SaveAs(Form("v2Dcentr%d-%d.png",mincentr, maxcentr));
+  cvv2obs->SaveAs(Form("v2obsDcentr%d-%d.png",mincentr, maxcentr));
+  cvevplres->SaveAs(Form("evplaneresocentr%d-%d.png",mincentr, maxcentr));
+  cvevplane->SaveAs(Form("evplanedistrcentr%d-%d.png",mincentr, maxcentr));
+  cvres->SaveAs(Form("evplresovscentr%d-%d.png",mincentr, maxcentr));
+  cvfitvsphi->SaveAs(Form("flowvsphicentr%d-%d.png",mincentr, maxcentr));
+  cvv2cor->SaveAs(Form("cmpv2centr%d-%d.png",mincentr, maxcentr));
+
+
+}
+
+void ReflectdNdphiPoints(Int_t mincentr,Int_t maxcentr){
+
+  gStyle->SetCanvasColor(0);
+  gStyle->SetTitleFillColor(0);
+  gStyle->SetOptTitle(0);
+  gStyle->SetStatColor(0);
+  gStyle->SetOptStat(0);
+
+  TString textleg=Form("centr%d-%d",mincentr,maxcentr);
+  TFile* fin=new TFile(Form("HistoInputv2Calc%s.root",textleg.Data()));
+  if(!fin->IsOpen()){
+    printf("HistoInputv2Calc%s.root not found",textleg.Data());
+    return;
+  }
+  TGraphErrors* grphi=(TGraphErrors*)fin->Get("gsigvsphi");
+  if(!grphi){
+    cout<<"gsigvsphi not found in "<<fin->GetName()<<endl;
+    return;
+  }
+
+  TGraphErrors* grphidouble=new TGraphErrors(0);
+  Int_t n=grphi->GetN();
+  //Int_t doublen=n*2;
+  Double_t xn=0,yn=0;
+  grphi->GetPoint(n-1,xn,yn);
+  xn+=grphi->GetErrorX(n-1);
+  cout<<"xn = "<<xn<<endl;
+  for(Int_t i=0;i<n;i++){
+    Double_t x,y;
+    grphi->GetPoint(i,x,y);
+    //x+=xn;
+    x=(-1)*x;
+    grphidouble->SetPoint(i,x,y);
+    grphidouble->SetPointError(i,grphi->GetErrorX(i),grphi->GetErrorY(i));
+  }
+  grphidouble->SetMarkerStyle(24);
+  grphidouble->SetName(Form("%sReflected",grphi->GetName()));
+  grphidouble->SetTitle(grphi->GetTitle());
+  grphidouble->GetXaxis()->SetRangeUser(0,6.5);
+
+  TMultiGraph* grtot=new TMultiGraph();
+  grtot->Add(grphi,"P");
+  grtot->Add(grphidouble,"P");
+  grtot->SetNameTitle("grphi2pi","dN/d#Delta#phi;#Delta#phi (rad);dN/d#Delta#phi");
+  TPaveText* pv=new TPaveText(0.3,0.7,0.7,0.9,"NDC");
+  pv->SetBorderSize(0);
+  pv->SetFillStyle(0);
+
+  TPaveText* pvcentr=new TPaveText(0.1,0.7,0.5,0.8,"NDC");
+  pvcentr->SetBorderSize(0);
+  pvcentr->SetFillStyle(0);
+  pvcentr->AddText(Form("Centrality %d-%d",mincentr,maxcentr));
+
+  TPaveText* pvwp=new TPaveText(0.5,0.1,0.9,0.3,"NDC");
+  pvwp->SetBorderSize(0);
+  pvwp->SetFillStyle(0);
+  pvwp->SetTextColor(kRed);
+  pvwp->SetTextFont(65);
+  pvwp->AddText("ALICE Work in progress");
+
+  pv->AddText("Open points are reflected");
+  TCanvas *cvphidouble=new TCanvas("cvphidouble","dN/dphi in 2*pi");
+  cvphidouble->cd();
+  grtot->Draw("AP");
+  pv->Draw();
+    
+  TCanvas *cvphi=new TCanvas("cvphi","dN/dphi");
+  cvphi->cd();
+  grphi->SetNameTitle(grphi->GetName(),"dN/d#Delta#phi;#Delta#phi (rad);dN/d#Delta#phi");
+  grphi->Draw("AP");
+  pvcentr->Draw();
+  pvwp->Draw();
+}
diff --git a/PWG3/vertexingHF/charmFlow/Extractv2from2Dhistos.C b/PWG3/vertexingHF/charmFlow/Extractv2from2Dhistos.C
new file mode 100644 (file)
index 0000000..2d6d212
--- /dev/null
@@ -0,0 +1,249 @@
+
+Double_t v2vsMass(Double_t *x, Double_t *par){
+  // Fit function for signal+background
+  // par[0] = S/B at the mass peak
+  // par[1] = mass
+  // par[2] = sigma
+  // par[3] = v2sig
+  // par[4] = v2back
+
+  Double_t fracsig=par[0]*TMath::Exp(-(x[0]-par[1])*(x[0]-par[1])/2./par[2]/par[2]);
+  Double_t fracbkg=1-fracsig;
+  return par[3]*fracsig+par[4]*fracbkg;
+}
+
+void Extractv2from2Dhistos(){
+
+  TFile* fil=new TFile("AnalysisResults.root");
+  TDirectoryFile* df=(TDirectoryFile*)fil->Get("PWG3_D2H_HFv2");
+  TList* lst=(TList*)df->Get("coutputv2D0");
+  Int_t minCent=30;
+  Int_t maxCent=50;
+  TH2F* hMassDphi=0x0;
+  for(Int_t iHisC=minCent; iHisC<=maxCent-5; iHisC+=5){    
+    TString hisname=Form("hMc2phicentr%d_%d",iHisC,iHisC+5);
+    TH2F* htmp=(TH2F*)lst->FindObject(hisname.Data());
+    if(iHisC==minCent) hMassDphi=(TH2F*)htmp->Clone("hMassCos2Dphi");
+    else hMassDphi->Add(htmp);
+    printf("Adding histogram %s\n",hisname.Data());
+  }
+
+  TH1F* hMass=(TH1F*)hMassDphi->ProjectionY();
+
+  Double_t sigmaRangeForSig=3.;
+  Double_t sigmaRangeForBkg=4.5;
+
+  gStyle->SetPalette(1);
+  gStyle->SetOptTitle(0);
+  TCanvas* c1=new TCanvas("c1");
+  c1->Divide(2,2);
+  c1->cd(1);
+  hMassDphi->Draw("colz");
+  c1->cd(2);
+
+  Int_t nMassBins=hMass->GetNbinsX();
+  Int_t hMinBin=3;
+  Int_t hMaxBin=nMassBins-2;
+  Double_t hmin=hMass->GetBinLowEdge(hMinBin);
+  Double_t hmax=hMass->GetBinLowEdge(hMaxBin)+hMass->GetBinWidth(hMaxBin);
+  Int_t factor4refl=0;
+  Float_t massD=TDatabasePDG::Instance()->GetParticle(421)->Mass();
+
+  AliHFMassFitter* fitter=new AliHFMassFitter(hMass,hmin,hmax,2,0,0);
+  fitter->SetReflectionSigmaFactor(factor4refl);
+  fitter->SetInitialGaussianMean(massD);
+  Bool_t out=fitter->MassFitter(0);
+  if(!out) return;
+  fitter->DrawHere(gPad);
+  Double_t sigfitter,esigfitter;
+  fitter->Signal(sigmaRangeForSig, sigfitter,esigfitter);
+  TH1F* hCos2PhiBkgLo=0x0;
+  TH1F* hCos2PhiBkgHi=0x0;
+  TH1F* hCos2PhiBkgLoScal=0x0;
+  TH1F* hCos2PhiBkgHiScal=0x0;
+  TH1F* hCos2PhiBkgAver=0x0;
+  TH1F* hCos2PhiSigReg=0x0;
+  TH1F* hCos2PhiSig=0x0;
+  Double_t mass=fitter->GetMean();
+  Double_t sigma=fitter->GetSigma();
+  TF1* fB1=fitter->GetBackgroundFullRangeFunc();
+  TF1* fB2=fitter->GetBackgroundRecalcFunc();
+  TF1* fSB=fitter->GetMassFunc();
+  Double_t minMassSig=mass-sigmaRangeForSig*sigma;
+  Double_t maxMassSig=mass+sigmaRangeForSig*sigma;
+  Int_t minBinSig=hMass->FindBin(minMassSig);
+  Int_t maxBinSig=hMass->FindBin(maxMassSig);
+  Double_t minMassSigBin=hMass->GetBinLowEdge(minBinSig);
+  Double_t maxMassSigBin=hMass->GetBinLowEdge(maxBinSig)+hMass->GetBinWidth(maxBinSig);
+  printf("Signal Fit Limits = %f %f\n",minMassSigBin,maxMassSigBin);
+  Double_t maxMassBkgLow=mass-sigmaRangeForBkg*sigma;
+  Int_t minBinBkgLow=hMinBin;
+  Int_t maxBinBkgLow=hMass->FindBin(maxMassBkgLow);
+  Double_t minMassBkgLowBin=hmin;
+  Double_t maxMassBkgLowBin=hMass->GetBinLowEdge(maxBinBkgLow)+hMass->GetBinWidth(maxBinBkgLow);
+  Double_t minMassBkgHi=mass+sigmaRangeForBkg*sigma;
+  Int_t minBinBkgHi=hMass->FindBin(minMassBkgHi);
+  Int_t maxBinBkgHi=hMaxBin;
+  Double_t minMassBkgHiBin=hMass->GetBinLowEdge(minBinBkgHi);
+  Double_t maxMassBkgHiBin=hmax;
+  printf("BKG Fit Limits = %f %f  && %f %f\n",minMassBkgLowBin,maxMassBkgLowBin,minMassBkgHiBin,maxMassBkgHiBin);
+  Double_t bkgSig=fB2->Integral(minMassSigBin,maxMassSigBin);
+  Double_t bkgLow=fB2->Integral(minMassBkgLowBin,maxMassBkgLowBin);
+  Double_t bkgHi=fB2->Integral(minMassBkgHiBin,maxMassBkgHiBin);
+  printf("Background integrals = %f %f %f\n",bkgLow,bkgSig,bkgHi);
+  TBox* bleft=new TBox(minMassBkgLowBin,0.,maxMassBkgLowBin,hMass->GetMaximum());
+  bleft->SetFillColor(kRed+1);
+  bleft->SetFillStyle(3002);
+  bleft->Draw();
+  TBox* bright=new TBox(minMassBkgHiBin,0.,maxMassBkgHiBin,hMass->GetMaximum());
+  bright->SetFillColor(kBlue+1);
+  bright->SetFillStyle(3002);
+  bright->Draw();
+  TBox* bsig=new TBox(minMassSigBin,0.,maxMassSigBin,hMass->GetMaximum()*2);
+  bsig->SetFillColor(1);
+  bsig->SetFillStyle(3002);
+  bsig->Draw();
+
+  TH1F* hCos2PhiBkgLo=(TH1F*)hMassDphi->ProjectionX("hCos2PhiBkgLo",minBinBkgLow,maxBinBkgLow);
+  TH1F* hCos2PhiBkgHi=(TH1F*)hMassDphi->ProjectionX("hCos2PhiBkgHi",minBinBkgHi,maxBinBkgHi);
+  TH1F* hCos2PhiSigReg=(TH1F*)hMassDphi->ProjectionX("hCos2PhiBkgSig",minBinSig,maxBinSig);
+  hCos2PhiBkgLo->Rebin(4);
+  hCos2PhiBkgHi->Rebin(4);
+  hCos2PhiSigReg->Rebin(4);
+  hCos2PhiSigReg->SetLineWidth(2);
+  hCos2PhiBkgLo->SetLineWidth(2);
+  hCos2PhiBkgHi->SetLineWidth(2);
+  hCos2PhiBkgLo->SetLineColor(kRed+1);
+  hCos2PhiBkgHi->SetLineColor(kBlue+1);
+  TH1F* hCos2PhiBkgLoScal=(TH1F*)hCos2PhiBkgLo->Clone("hCos2PhiBkgLoScal");
+  hCos2PhiBkgLoScal->Scale(bkgSig/bkgLow);
+  TH1F* hCos2PhiBkgHiScal=(TH1F*)hCos2PhiBkgHi->Clone("hCos2PhiBkgHiScal");
+  hCos2PhiBkgHiScal->Scale(bkgSig/bkgHi);
+  hCos2PhiBkgLoScal->SetLineWidth(2);
+  hCos2PhiBkgHiScal->SetLineWidth(2);
+  hCos2PhiBkgLoScal->SetLineColor(kRed+1);
+  hCos2PhiBkgHiScal->SetLineColor(kBlue+1);
+  TH1F* hCos2PhiBkgAver=(TH1F*)hCos2PhiBkgLoScal->Clone("hCos2PhiBkgAver");
+  hCos2PhiBkgAver->Add(hCos2PhiBkgHiScal);
+  hCos2PhiBkgAver->Scale(0.5);
+  hCos2PhiBkgAver->SetLineWidth(2);
+  hCos2PhiBkgAver->SetLineColor(kGreen+1);
+  TH1F* hCos2PhiSig=(TH1F*)hCos2PhiSigReg->Clone("hCos2PhiSig");
+  hCos2PhiSig->Add(hCos2PhiBkgAver,-1.);   
+  printf("v2=%f +- %f\n",hCos2PhiSig->GetMean(),hCos2PhiSig->GetMeanError());
+  c1->cd(3);
+  hCos2PhiSigReg->Draw();
+  hCos2PhiBkgLoScal->Draw("same");
+  hCos2PhiBkgHiScal->Draw("same");
+  hCos2PhiBkgAver->Draw("same");
+  TLegend* leg0=new TLegend(0.3,0.6,0.75,0.89);
+  leg0->SetFillColor(0);
+  TLegendEntry* ent=leg0->AddEntry(hCos2PhiBkgSig,"Signal region","L");
+  ent->SetTextColor(hCos2PhiBkgSig->GetLineColor());
+  ent=leg0->AddEntry(hCos2PhiBkgLoScal,"Left side band","L");
+  ent->SetTextColor(hCos2PhiBkgLoScal->GetLineColor());
+  ent=leg0->AddEntry(hCos2PhiBkgHiScal,"Right side band","L");
+  ent->SetTextColor(hCos2PhiBkgHiScal->GetLineColor());
+  ent=leg0->AddEntry(hCos2PhiBkgAver,"Average of side bands","L");
+  ent->SetTextColor(hCos2PhiBkgAver->GetLineColor());
+  leg0->Draw();
+  c1->cd(4);
+  hCos2PhiSig->Draw("EP");
+  TPaveText* t0= new TPaveText(0.15,0.70,0.45,0.89,"NDC");
+  t0->SetFillColor(0);
+  t0->AddText(Form("v2=%.3f+-%.3f\n",hCos2PhiSig->GetMean(),hCos2PhiSig->GetMeanError()));
+  t0->Draw();
+
+  printf("Signal from mass fitter = %f  Signal from subracted histo= %f\n",
+        sigfitter,hCos2PhiSig->Integral());
+
+  Int_t npars=fSB->GetNpar();
+  Double_t sigma=fSB->GetParameter(npars-1);
+  Double_t mass=fSB->GetParameter(npars-2);
+  Double_t integr=fSB->GetParameter(npars-3);
+  Double_t sOverAll=(fSB->Eval(mass)-fB2->Eval(mass))/fSB->Eval(mass);
+  printf("mass=%f  S+B=%f   bkg=%f S/(S+B)=%f\n",mass,fSB->Eval(mass),fB2->Eval(mass),sOverAll);
+  printf("Number of parameters: %d. Signal params: %f %f %f\n",npars,mass,sigma,integr);
+  TF1* fSig=new TF1("fSig","[0]/TMath::Sqrt(2.*TMath::Pi())/[2]*TMath::Exp(-(x-[1])*(x-[1])/2./[2]/[2])",hMassDphi->GetYaxis()->GetXmin(),hMassDphi->GetYaxis()->GetXmax());
+  fSig->SetParameters(integr,mass,sigma);
+  TH1F* hAverCos2Phi=new TH1F("hAverCos2Phi","",hMassDphi->GetNbinsY(),hMassDphi->GetYaxis()->GetXmin(),hMassDphi->GetYaxis()->GetXmax());
+  TH1F* hFractionSig=new TH1F("hFractionSig","",hMassDphi->GetNbinsY(),hMassDphi->GetYaxis()->GetXmin(),hMassDphi->GetYaxis()->GetXmax());
+  TH1F* hFractionBkg=new TH1F("hFractionBkg","",hMassDphi->GetNbinsY(),hMassDphi->GetYaxis()->GetXmin(),hMassDphi->GetYaxis()->GetXmax());
+
+  for(Int_t iBin=1; iBin<=hMassDphi->GetNbinsY(); iBin++){
+    TH1F* htemp=(TH1F*)hMassDphi->ProjectionX("htemp",iBin,iBin);
+    hAverCos2Phi->SetBinContent(iBin,htemp->GetMean());
+    hAverCos2Phi->SetBinError(iBin,htemp->GetMeanError());
+    Double_t sig=fSig->Eval(hFractionSig->GetBinCenter(iBin));
+    Double_t bkg=fB2->Eval(hFractionSig->GetBinCenter(iBin));
+    if(bkg<1 && sig<1){
+      hFractionSig->SetBinContent(iBin,0.);
+      hFractionSig->SetBinError(iBin,0.);
+      hFractionBkg->SetBinContent(iBin,1.);
+      hFractionBkg->SetBinError(iBin,0.);
+    }else{
+      Double_t fracs=sig/(sig+bkg);
+      Double_t fracb=bkg/(sig+bkg);
+      Double_t efracs=0.;//TMath::Sqrt(fracs*(1.-fracs)/(sig+bkg));
+      Double_t efracb=0.;//TMath::Sqrt(fracb*(1.-fracb)/(sig+bkg));
+       
+      hFractionSig->SetBinContent(iBin,fracs),
+      hFractionSig->SetBinError(iBin,efracs);
+      hFractionBkg->SetBinContent(iBin,fracb);      
+      hFractionBkg->SetBinError(iBin,efracb);
+    }
+    delete htemp;
+  }
+  
+  TF1* fv2=new TF1("fv2",v2vsMass,hMassDphi->GetYaxis()->GetXmin(),hMassDphi->GetYaxis()->GetXmax(),5);
+  fv2->SetParameter(0,sOverAll);
+  fv2->SetParameter(1,mass);
+  fv2->SetParameter(2,sigma);
+  fv2->SetParameter(3,0.2);
+  fv2->SetParameter(4,0.2);
+  fv2->FixParameter(0,sOverAll);
+  fv2->FixParameter(1,mass);
+  fv2->FixParameter(2,sigma);
+
+  hAverCos2Phi->Rebin(2);
+  hAverCos2Phi->Scale(0.5);
+
+  TCanvas* c2=new TCanvas("c2");
+  c2->Divide(2,2);
+  c2->cd(1);
+  hMassDphi->Draw("colz");
+  c2->cd(2);
+  hMass->Rebin(2);
+  hMass->SetMinimum(0.);
+  hMass->SetMarkerStyle(20);
+  hMass->Draw("E");
+  fSB->Draw("same");
+  fSig->Draw("same");
+  fB2->Draw("same");
+  c2->cd(3);
+  hFractionSig->SetMaximum(1.2);
+  hFractionSig->Draw();
+  hFractionSig->GetXaxis()->SetTitle("Mass (GeV/c^2)");
+  hFractionSig->GetYaxis()->SetTitle("Fraction");
+  hFractionBkg->SetLineColor(2);
+  hFractionBkg->Draw("same");
+  TLegend* leg1=new TLegend(0.15,0.15,0.35,0.35);
+  leg1->SetFillColor(0);
+  ent=leg1->AddEntry(hFractionSig,"S/(S+B)","L");
+  ent->SetTextColor(hFractionSig->GetLineColor());
+  ent=leg1->AddEntry(hFractionBkg,"B/(S+B)","L");
+  ent->SetTextColor(hFractionBkg->GetLineColor());
+  leg1->Draw();
+  c2->cd(4);
+  hAverCos2Phi->Fit(fv2);
+  hAverCos2Phi->GetXaxis()->SetTitle("Mass (GeV/c^2)");
+  hAverCos2Phi->GetYaxis()->SetTitle("v_2^{obs}");
+  TPaveText* t1= new TPaveText(0.55,0.70,0.89,0.89,"NDC");
+  t1->SetFillColor(0);
+  t1->AddText(Form("v2sig=%.3f+-%.3f\n",fv2->GetParameter(3),fv2->GetParError(3)));
+  t1->AddText(Form("v2bkg=%.3f+-%.3f\n",fv2->GetParameter(4),fv2->GetParError(4)));
+  t1->Draw();
+  printf("v2(from fit)=%f+-%f\n",fv2->GetParameter(3),fv2->GetParError(3));
+
+}