Adding Tab values for QpPb analysis

[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / macros / AddTaskDvsMultiplicity.C

diff --git a/PWGHF/vertexingHF/macros/AddTaskDvsMultiplicity.C b/PWGHF/vertexingHF/macros/AddTaskDvsMultiplicity.C
old mode 100644 (file)
new mode 100755 (executable)
index 27fb991..b9f2eb5
--- a/PWGHF/vertexingHF/macros/AddTaskDvsMultiplicity.C
+++ b/PWGHF/vertexingHF/macros/AddTaskDvsMultiplicity.C
@@ -5,7 +5,13 @@ AliAnalysisTaskSEDvsMultiplicity *AddTaskDvsMultiplicity(Int_t system=0,
                                                         TString finDirname="Loose",
                                                         TString filename="",
                                                         TString finAnObjname="AnalysisCuts", 
-                                                        TString estimatorFilename="")
+                                                        TString estimatorFilename="",
+                                                        Double_t refMult=9.26,
+                                                        Bool_t subtractDau=kFALSE,
+                                                        Bool_t NchWeight=kFALSE,
+                                                        Int_t recoEstimator = AliAnalysisTaskSEDvsMultiplicity::kNtrk10,
+                                                        Int_t MCEstimator = AliAnalysisTaskSEDvsMultiplicity::kEta10,
+                                                        Bool_t isPPbData=kFALSE)
 {
   //
   // Test macro for the AliAnalysisTaskSE for D+ candidates
@@ -62,35 +68,75 @@ AliAnalysisTaskSEDvsMultiplicity *AddTaskDvsMultiplicity(Int_t system=0,
     Name="DStar";
   }
 
-
-  AliAnalysisTaskSEDvsMultiplicity *dMultTask = new AliAnalysisTaskSEDvsMultiplicity("dMultAnalysis",pdgMeson,analysiscuts);
+  AliAnalysisTaskSEDvsMultiplicity *dMultTask = new AliAnalysisTaskSEDvsMultiplicity("dMultAnalysis",pdgMeson,analysiscuts,isPPbData);
   dMultTask->SetReadMC(readMC);  
   dMultTask->SetDebugLevel(0);
   dMultTask->SetUseBit(kTRUE);
   dMultTask->SetDoImpactParameterHistos(kFALSE);
+  dMultTask->SetSubtractTrackletsFromDaughters(subtractDau);
+  dMultTask->SetMultiplicityEstimator(recoEstimator);
+  dMultTask->SetMCPrimariesEstimator(MCEstimator);
+  dMultTask->SetMCOption(MCOption);
+  if(isPPbData) dMultTask->SetIsPPbData();
+
+  if(NchWeight){
+    TH1F *hNchPrimaries = NULL; 
+    if(isPPbData) hNchPrimaries = (TH1F*)filecuts->Get("hNtrUnCorrEvWithDWeight");
+    else hNchPrimaries = (TH1F*)filecuts->Get("hGenPrimaryParticlesInelGt0");
+    if(hNchPrimaries) {
+      dMultTask->UseMCNchWeight(true);
+      dMultTask->SetHistoNchWeight(hNchPrimaries);
+    } else {
+      AliFatal("Histogram for multiplicity weights not found");
+      return 0x0;
+    }
+  }
 
   if(pdgMeson==421) { 
     dMultTask->SetMassLimits(1.5648,2.1648);
-    dMultTask->SetMassBins(300);
+    dMultTask->SetNMassBins(200);
   }else if(pdgMeson==411)dMultTask->SetMassLimits(pdgMeson,0.2);
   
   if(estimatorFilename.EqualTo("") ) {
     printf("Estimator file not provided, multiplcity corrected histograms will not be filled\n");
   } else{
-    const Char_t* periodNames[4] = {"LHC10b", "LHC10c", "LHC10d", "LHC10e"};
-    TProfile* multEstimatorAvg[4];                       
+                     
     TFile* fileEstimator=TFile::Open(estimatorFilename.Data());
     if(!fileEstimator)  {
-      AliFatal("File with multiplicity estimator not found\n"); 
+      AliFatal("File with multiplicity estimator not found\n");
       return;
     }
-    for(Int_t ip=0; ip<4; ip++) {
-      multEstimatorAvg[ip] = (TProfile*)(fileEstimator->Get(Form("SPDmult10_%s",periodNames[ip]))->Clone(Form("SPDmult10_%s_clone",periodNames[ip])));  
+
+    dMultTask->SetReferenceMultiplcity(refMult);
+
+    if (isPPbData) {    //Only use two profiles if pPb
+      const Char_t* periodNames[2] = {"LHC13b", "LHC13c"};
+      TProfile* multEstimatorAvg[2];
+      for(Int_t ip=0; ip<2; ip++) {
+       multEstimatorAvg[ip] = (TProfile*)(fileEstimator->Get(Form("SPDmult10_%s",periodNames[ip]))->Clone(Form("SPDmult10_%s_clone",periodNames[ip])));
+       if (!multEstimatorAvg[ip]) {
+         AliFatal(Form("Multiplicity estimator for %s not found! Please check your estimator file",periodNames[ip]));
+         return;
+       }
+      }
+      dMultTask->SetMultiplVsZProfileLHC13b(multEstimatorAvg[0]);
+      dMultTask->SetMultiplVsZProfileLHC13c(multEstimatorAvg[1]);
+    }
+    else {
+      const Char_t* periodNames[4] = {"LHC10b", "LHC10c", "LHC10d", "LHC10e"};
+      TProfile* multEstimatorAvg[4];                       
+      for(Int_t ip=0; ip<4; ip++) {
+       multEstimatorAvg[ip] = (TProfile*)(fileEstimator->Get(Form("SPDmult10_%s",periodNames[ip]))->Clone(Form("SPDmult10_%s_clone",periodNames[ip])));
+       if (!multEstimatorAvg[ip]) {
+         AliFatal(Form("Multiplicity estimator for %s not found! Please check your estimator file",periodNames[ip]));
+         return;
+       }
+      }
+      dMultTask->SetMultiplVsZProfileLHC10b(multEstimatorAvg[0]);
+      dMultTask->SetMultiplVsZProfileLHC10c(multEstimatorAvg[1]);
+      dMultTask->SetMultiplVsZProfileLHC10d(multEstimatorAvg[2]);
+      dMultTask->SetMultiplVsZProfileLHC10e(multEstimatorAvg[3]);
     }
-    dMultTask->SetMultiplVsZProfileLHC10b(multEstimatorAvg[0]);
-    dMultTask->SetMultiplVsZProfileLHC10c(multEstimatorAvg[1]);
-    dMultTask->SetMultiplVsZProfileLHC10d(multEstimatorAvg[2]);
-    dMultTask->SetMultiplVsZProfileLHC10e(multEstimatorAvg[3]);
   }
   mgr->AddTask(dMultTask);
   
@@ -100,16 +146,18 @@ AliAnalysisTaskSEDvsMultiplicity *AddTaskDvsMultiplicity(Int_t system=0,
   TString outname = "coutput";
   TString cutsname = "coutputCuts";
   TString normname = "coutputNorm";
+  TString profname = "coutputProf";
   
   inname += Name.Data();
   outname += Name.Data();
   cutsname += Name.Data();
   normname += Name.Data();
+  profname += Name.Data();
   inname += finDirname.Data();
   outname += finDirname.Data();
   cutsname += finDirname.Data();
   normname += finDirname.Data();
-
+  profname += finDirname.Data();
 
   AliAnalysisDataContainer *cinput = mgr->CreateContainer(inname,TChain::Class(),AliAnalysisManager::kInputContainer);
 
@@ -121,6 +169,7 @@ AliAnalysisTaskSEDvsMultiplicity *AddTaskDvsMultiplicity(Int_t system=0,
   AliAnalysisDataContainer *coutputCuts = mgr->CreateContainer(cutsname,TList::Class(),AliAnalysisManager::kOutputContainer,outputfile.Data());
   AliAnalysisDataContainer *coutput = mgr->CreateContainer(outname,TList::Class(),AliAnalysisManager::kOutputContainer,outputfile.Data());
   AliAnalysisDataContainer *coutputNorm = mgr->CreateContainer(normname,TList::Class(),AliAnalysisManager::kOutputContainer,outputfile.Data());
+  AliAnalysisDataContainer *coutputProf = mgr->CreateContainer(profname,TList::Class(),AliAnalysisManager::kOutputContainer,outputfile.Data());
   
   mgr->ConnectInput(dMultTask,0,mgr->GetCommonInputContainer());
   
@@ -129,6 +178,8 @@ AliAnalysisTaskSEDvsMultiplicity *AddTaskDvsMultiplicity(Int_t system=0,
   mgr->ConnectOutput(dMultTask,2,coutputCuts);
 
   mgr->ConnectOutput(dMultTask,3,coutputNorm);  
- 
+
+  mgr->ConnectOutput(dMultTask,4,coutputProf);
+
   return dMultTask;
 }