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