AliPhysics  32b88a8 (32b88a8)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AddTaskDvsMultiplicity.C
Go to the documentation of this file.
2  Bool_t readMC=kFALSE,
3  Int_t MCOption=0,
4  Int_t pdgMeson=411,
5  TString finDirname="Loose",
6  TString filename="",
7  TString finAnObjname="AnalysisCuts",
8  TString estimatorFilename="",
9  Double_t refMult=9.26,
10  Bool_t subtractDau=kFALSE,
11  Int_t NchWeight=0,
14  Bool_t isPPbData=kFALSE)
15 {
16  //
17  // Macro for the AliAnalysisTaskSE for D candidates vs Multiplicity
18  // Invariant mass histogram in pt and multiplicity bins in a 3D histogram
19  // different estimators implemented
20  //==============================================================================
21 
22  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
23  if (!mgr) {
24  ::Error("AddTaskDvsMultiplicity", "No analysis manager to connect to.");
25  }
26 
27  Bool_t stdcuts=kFALSE;
28  TFile* filecuts;
29  if( filename.EqualTo("") ) {
30  stdcuts=kTRUE;
31  } else {
32  filecuts=TFile::Open(filename.Data());
33  if(!filecuts ||(filecuts&& !filecuts->IsOpen())){
34  AliFatal("Input file not found : check your cut object");
35  }
36  }
37 
38 
39  //Analysis Task
40  AliRDHFCuts *analysiscuts=0x0;
41 
42  TString Name="";
43  if(pdgMeson==411){
44  if(stdcuts) {
45  analysiscuts = new AliRDHFCutsDplustoKpipi();
46  if (system == 0) analysiscuts->SetStandardCutsPP2010();
47  else analysiscuts->SetStandardCutsPbPb2011();
48  }
49  else analysiscuts = (AliRDHFCutsDplustoKpipi*)filecuts->Get(finAnObjname);
50  Name="Dplus";
51  }else if(pdgMeson==421){
52  if(stdcuts) {
53  analysiscuts = new AliRDHFCutsD0toKpi();
54  if (system == 0) analysiscuts->SetStandardCutsPP2010();
55  else analysiscuts->SetStandardCutsPbPb2011();
56  }
57  else analysiscuts = (AliRDHFCutsD0toKpi*)filecuts->Get(finAnObjname);
58  Name="D0";
59  }else if(pdgMeson==413){
60  if(stdcuts) {
61  analysiscuts = new AliRDHFCutsDStartoKpipi();
62  if (system == 0) analysiscuts->SetStandardCutsPP2010();
63  else analysiscuts->SetStandardCutsPbPb2011();
64  }
65  else analysiscuts = (AliRDHFCutsDStartoKpipi*)filecuts->Get(finAnObjname);
66  Name="DStar";
67  }
68 
69  AliAnalysisTaskSEDvsMultiplicity *dMultTask = new AliAnalysisTaskSEDvsMultiplicity("dMultAnalysis",pdgMeson,analysiscuts,isPPbData);
70  dMultTask->SetReadMC(readMC);
71  dMultTask->SetDebugLevel(0);
72  dMultTask->SetUseBit(kTRUE);
73  dMultTask->SetDoImpactParameterHistos(kFALSE);
74  dMultTask->SetSubtractTrackletsFromDaughters(subtractDau);
75  dMultTask->SetMultiplicityEstimator(recoEstimator);
76  dMultTask->SetMCPrimariesEstimator(MCEstimator);
77  dMultTask->SetMCOption(MCOption);
78  if(isPPbData) dMultTask->SetIsPPbData();
79 
80  if(NchWeight){
81  TH1F *hNchPrimaries = NULL;
82  TH1F *hMeasNchPrimaries = NULL;
83  if(NchWeight==1){
84  if(isPPbData) {
85  hNchPrimaries = (TH1F*)filecuts->Get("hNtrUnCorrEvWithDWeight"); // MC distribution
86  }
87  else hNchPrimaries = (TH1F*)filecuts->Get("hGenPrimaryParticlesInelGt0");
88  if(hNchPrimaries) {
89  dMultTask->UseMCNchWeight(NchWeight);
90  dMultTask->SetHistoNchWeight(hNchPrimaries);
91  } else {
92  AliFatal("Histogram for Nch multiplicity weights not found");
93  return 0x0;
94  }
95  hMeasNchPrimaries = (TH1F*)filecuts->Get("hMeasNtrUnCorrEvWithD"); // data distribution
96  if(hMeasNchPrimaries) {
97  dMultTask->SetMeasuredNchHisto(hMeasNchPrimaries);
98  }
99  }
100  else if(NchWeight==2){
101  hNchPrimaries = (TH1F*)filecuts->Get("hNtrUnCorrEvWithDWeight"); // MC distribution
102  hMeasNchPrimaries = (TH1F*)filecuts->Get("hMeasNtrUnCorrEvWithD"); // data distribution
103  if(hNchPrimaries && hMeasNchPrimaries) {
104  dMultTask->UseMCNchWeight(NchWeight);
105  dMultTask->SetHistoNchWeight(hNchPrimaries);
106  dMultTask->SetMeasuredNchHisto(hMeasNchPrimaries);
107  } else {
108  AliFatal("Histogram for Ntrk multiplicity weights not found");
109  return 0x0;
110  }
111  }
112  }
113 
114 
115  if(pdgMeson==421) {
116  dMultTask->SetMassLimits(1.5648,2.1648);
117  dMultTask->SetNMassBins(200);
118  }else if(pdgMeson==411)dMultTask->SetMassLimits(pdgMeson,0.2);
119 
120  if(estimatorFilename.EqualTo("") ) {
121  printf("Estimator file not provided, multiplcity corrected histograms will not be filled\n");
122  } else{
123 
124  TFile* fileEstimator=TFile::Open(estimatorFilename.Data());
125  if(!fileEstimator) {
126  AliFatal("File with multiplicity estimator not found\n");
127  return;
128  }
129 
130  dMultTask->SetReferenceMultiplcity(refMult);
131 
132  const Char_t* profilebasename="SPDmult10";
133  if(recoEstimator==AliAnalysisTaskSEDvsMultiplicity::kVZEROA || recoEstimator==AliAnalysisTaskSEDvsMultiplicity::kVZEROAEq) profilebasename="VZEROAmult";
134  else if(recoEstimator==AliAnalysisTaskSEDvsMultiplicity::kVZERO || recoEstimator==AliAnalysisTaskSEDvsMultiplicity::kVZEROEq) profilebasename="VZEROMmult";
135  cout<<endl<<endl<<" profilebasename="<<profilebasename<<endl<<endl;
136 
137  if (isPPbData) { //Only use two profiles if pPb
138  const Char_t* periodNames[2] = {"LHC13b", "LHC13c"};
139  TProfile* multEstimatorAvg[2];
140  for(Int_t ip=0; ip<2; ip++) {
141  cout<< " Trying to get "<<Form("%s_%s",profilebasename,periodNames[ip])<<endl;
142  multEstimatorAvg[ip] = (TProfile*)(fileEstimator->Get(Form("%s_%s",profilebasename,periodNames[ip]))->Clone(Form("%s_%s_clone",profilebasename,periodNames[ip])));
143  if (!multEstimatorAvg[ip]) {
144  AliFatal(Form("Multiplicity estimator for %s not found! Please check your estimator file",periodNames[ip]));
145  return;
146  }
147  }
148  dMultTask->SetMultiplVsZProfileLHC13b(multEstimatorAvg[0]);
149  dMultTask->SetMultiplVsZProfileLHC13c(multEstimatorAvg[1]);
150  }
151  else {
152  const Char_t* periodNames[4] = {"LHC10b", "LHC10c", "LHC10d", "LHC10e"};
153  TProfile* multEstimatorAvg[4];
154  for(Int_t ip=0; ip<4; ip++) {
155  multEstimatorAvg[ip] = (TProfile*)(fileEstimator->Get(Form("%s_%s",profilebasename,periodNames[ip]))->Clone(Form("%s_%s_clone",profilebasename,periodNames[ip])));
156  if (!multEstimatorAvg[ip]) {
157  AliFatal(Form("Multiplicity estimator for %s not found! Please check your estimator file",periodNames[ip]));
158  return;
159  }
160  }
161  dMultTask->SetMultiplVsZProfileLHC10b(multEstimatorAvg[0]);
162  dMultTask->SetMultiplVsZProfileLHC10c(multEstimatorAvg[1]);
163  dMultTask->SetMultiplVsZProfileLHC10d(multEstimatorAvg[2]);
164  dMultTask->SetMultiplVsZProfileLHC10e(multEstimatorAvg[3]);
165  }
166  }
167  mgr->AddTask(dMultTask);
168 
169  // Create containers for input/output
170 
171  TString inname = "cinput";
172  TString outname = "coutput";
173  TString cutsname = "coutputCuts";
174  TString normname = "coutputNorm";
175  TString profname = "coutputProf";
176 
177  inname += Name.Data();
178  outname += Name.Data();
179  cutsname += Name.Data();
180  normname += Name.Data();
181  profname += Name.Data();
182  inname += finDirname.Data();
183  outname += finDirname.Data();
184  cutsname += finDirname.Data();
185  normname += finDirname.Data();
186  profname += finDirname.Data();
187 
188  AliAnalysisDataContainer *cinput = mgr->CreateContainer(inname,TChain::Class(),AliAnalysisManager::kInputContainer);
189 
190  TString outputfile = AliAnalysisManager::GetCommonFileName();
191  outputfile += ":PWG3_D2H_DMult_";
192  outputfile += Name.Data();
193  outputfile += finDirname.Data();
194 
195  AliAnalysisDataContainer *coutputCuts = mgr->CreateContainer(cutsname,TList::Class(),AliAnalysisManager::kOutputContainer,outputfile.Data());
196  AliAnalysisDataContainer *coutput = mgr->CreateContainer(outname,TList::Class(),AliAnalysisManager::kOutputContainer,outputfile.Data());
197  AliAnalysisDataContainer *coutputNorm = mgr->CreateContainer(normname,TList::Class(),AliAnalysisManager::kOutputContainer,outputfile.Data());
198  AliAnalysisDataContainer *coutputProf = mgr->CreateContainer(profname,TList::Class(),AliAnalysisManager::kOutputContainer,outputfile.Data());
199 
200  mgr->ConnectInput(dMultTask,0,mgr->GetCommonInputContainer());
201 
202  mgr->ConnectOutput(dMultTask,1,coutput);
203 
204  mgr->ConnectOutput(dMultTask,2,coutputCuts);
205 
206  mgr->ConnectOutput(dMultTask,3,coutputNorm);
207 
208  mgr->ConnectOutput(dMultTask,4,coutputProf);
209 
210  return dMultTask;
211 }
const char * filename
Definition: TestFCM.C:1
double Double_t
Definition: External.C:58
char Char_t
Definition: External.C:18
void SetMassLimits(Double_t lowlimit, Double_t uplimit)
Class for cuts on AOD reconstructed D+->Kpipi.
int Int_t
Definition: External.C:63
virtual void SetStandardCutsPbPb2011()
Definition: AliRDHFCuts.h:48
void UseMCNchWeight(Int_t flag)
Nch Ntrk weights on MC.
virtual void SetStandardCutsPP2010()
Definition: AliRDHFCuts.h:46
bool Bool_t
Definition: External.C:53
AliAnalysisTaskSEDvsMultiplicity * AddTaskDvsMultiplicity(Int_t system=0, Bool_t readMC=kFALSE, Int_t MCOption=0, Int_t pdgMeson=411, TString finDirname="Loose", TString filename="", TString finAnObjname="AnalysisCuts", TString estimatorFilename="", Double_t refMult=9.26, Bool_t subtractDau=kFALSE, Int_t NchWeight=0, Int_t recoEstimator=AliAnalysisTaskSEDvsMultiplicity::kNtrk10, Int_t MCEstimator=AliAnalysisTaskSEDvsMultiplicity::kEta10, Bool_t isPPbData=kFALSE)