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