AliPhysics  58f3d52 (58f3d52)
 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  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  }
77 
78  AliAnalysisTaskSEDvsMultiplicity *dMultTask = new AliAnalysisTaskSEDvsMultiplicity("dMultAnalysis",pdgMeson,analysiscuts,isPPbData);
79  dMultTask->SetReadMC(readMC);
80  dMultTask->SetDebugLevel(0);
81  dMultTask->SetUseBit(kTRUE);
82  dMultTask->SetDoImpactParameterHistos(kFALSE);
83  dMultTask->SetSubtractTrackletsFromDaughters(subtractDau);
84  dMultTask->SetMultiplicityEstimator(recoEstimator);
85  dMultTask->SetMCPrimariesEstimator(MCEstimator);
86  dMultTask->SetMCOption(MCOption);
87  if(isPPbData) dMultTask->SetIsPPbData();
88 
89  if(NchWeight){
90  TH1F *hNchPrimaries = NULL;
91  TH1F *hMeasNchPrimaries = NULL;
92  if(NchWeight==1){
93  if(isPPbData) {
94  hNchPrimaries = (TH1F*)filecuts->Get("hNtrUnCorrEvWithDWeight"); // MC distribution
95  }
96  else hNchPrimaries = (TH1F*)filecuts->Get("hGenPrimaryParticlesInelGt0");
97  if(hNchPrimaries) {
98  dMultTask->UseMCNchWeight(NchWeight);
99  dMultTask->SetHistoNchWeight(hNchPrimaries);
100  } else {
101  AliFatal("Histogram for Nch multiplicity weights not found");
102  return 0x0;
103  }
104  hMeasNchPrimaries = (TH1F*)filecuts->Get("hMeasNtrUnCorrEvWithD"); // data distribution
105  if(hMeasNchPrimaries) {
106  dMultTask->SetMeasuredNchHisto(hMeasNchPrimaries);
107  }
108  }
109  else if(NchWeight==2){
110  hNchPrimaries = (TH1F*)filecuts->Get("hNtrUnCorrEvWithDWeight"); // MC distribution
111  hMeasNchPrimaries = (TH1F*)filecuts->Get("hMeasNtrUnCorrEvWithD"); // data distribution
112  if(hNchPrimaries && hMeasNchPrimaries) {
113  dMultTask->UseMCNchWeight(NchWeight);
114  dMultTask->SetHistoNchWeight(hNchPrimaries);
115  dMultTask->SetMeasuredNchHisto(hMeasNchPrimaries);
116  } else {
117  AliFatal("Histogram for Ntrk multiplicity weights not found");
118  return 0x0;
119  }
120  }
121  }
122 
123 
124  if(pdgMeson==421) {
125  dMultTask->SetMassLimits(1.5648,2.1648);
126  dMultTask->SetNMassBins(200);
127  }else if(pdgMeson==411)dMultTask->SetMassLimits(pdgMeson,0.2);
128 
129  if(estimatorFilename.EqualTo("") ) {
130  printf("Estimator file not provided, multiplcity corrected histograms will not be filled\n");
131  } else{
132 
133  TFile* fileEstimator=TFile::Open(estimatorFilename.Data());
134  if(!fileEstimator) {
135  AliFatal("File with multiplicity estimator not found\n");
136  return;
137  }
138 
139  dMultTask->SetReferenceMultiplcity(refMult);
140 
141  const Char_t* profilebasename="SPDmult10";
142  if(recoEstimator==AliAnalysisTaskSEDvsMultiplicity::kVZEROA || recoEstimator==AliAnalysisTaskSEDvsMultiplicity::kVZEROAEq) profilebasename="VZEROAmult";
143  else if(recoEstimator==AliAnalysisTaskSEDvsMultiplicity::kVZERO || recoEstimator==AliAnalysisTaskSEDvsMultiplicity::kVZEROEq) profilebasename="VZEROMmult";
144  cout<<endl<<endl<<" profilebasename="<<profilebasename<<endl<<endl;
145 
146  if (isPPbData) {
147  if(year == 16) {
148  const Char_t* periodNames[4] = {"LHC16q_265499to265525_265309to265387", "LHC16q_265435","LHC16q_265388to265427","LHC16t_267163to267166"};
149  TProfile* multEstimatorAvg[4];
150  for(Int_t ip=0; ip<4; ip++) {
151  cout<< " Trying to get "<<Form("%s_%s",profilebasename,periodNames[ip])<<endl;
152  multEstimatorAvg[ip] = (TProfile*)(fileEstimator->Get(Form("%s_%s",profilebasename,periodNames[ip]))->Clone(Form("%s_%s_clone",profilebasename,periodNames[ip])));
153  if (!multEstimatorAvg[ip]) {
154  AliFatal(Form("Multiplicity estimator for %s not found! Please check your estimator file",periodNames[ip]));
155  return;
156  }
157  }
158  dMultTask->SetMultiplVsZProfileLHC16qt1stBunch(multEstimatorAvg[0]);
159  dMultTask->SetMultiplVsZProfileLHC16qt2ndBunch(multEstimatorAvg[1]);
160  dMultTask->SetMultiplVsZProfileLHC16qt3rdBunch(multEstimatorAvg[2]);
161  dMultTask->SetMultiplVsZProfileLHC16qt4thBunch(multEstimatorAvg[3]);
162  }
163  else {
164  //Only use two profiles if pPb 2013
165  const Char_t* periodNames[2] = {"LHC13b", "LHC13c"};
166  TProfile* multEstimatorAvg[2];
167  for(Int_t ip=0; ip<2; ip++) {
168  cout<< " Trying to get "<<Form("%s_%s",profilebasename,periodNames[ip])<<endl;
169  multEstimatorAvg[ip] = (TProfile*)(fileEstimator->Get(Form("%s_%s",profilebasename,periodNames[ip]))->Clone(Form("%s_%s_clone",profilebasename,periodNames[ip])));
170  if (!multEstimatorAvg[ip]) {
171  AliFatal(Form("Multiplicity estimator for %s not found! Please check your estimator file",periodNames[ip]));
172  return;
173  }
174  }
175  dMultTask->SetMultiplVsZProfileLHC13b(multEstimatorAvg[0]);
176  dMultTask->SetMultiplVsZProfileLHC13c(multEstimatorAvg[1]);
177  }
178  }
179  else {
180  const Char_t* periodNames[4] = {"LHC10b", "LHC10c", "LHC10d", "LHC10e"};
181  TProfile* multEstimatorAvg[4];
182  for(Int_t ip=0; ip<4; ip++) {
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  AliFatal(Form("Multiplicity estimator for %s not found! Please check your estimator file",periodNames[ip]));
186  return;
187  }
188  }
189  dMultTask->SetMultiplVsZProfileLHC10b(multEstimatorAvg[0]);
190  dMultTask->SetMultiplVsZProfileLHC10c(multEstimatorAvg[1]);
191  dMultTask->SetMultiplVsZProfileLHC10d(multEstimatorAvg[2]);
192  dMultTask->SetMultiplVsZProfileLHC10e(multEstimatorAvg[3]);
193  }
194  }
195  mgr->AddTask(dMultTask);
196 
197  // Create containers for input/output
198 
199  TString inname = "cinput";
200  TString outname = "coutput";
201  TString cutsname = "coutputCuts";
202  TString normname = "coutputNorm";
203  TString profname = "coutputProf";
204 
205  inname += Name.Data();
206  outname += Name.Data();
207  cutsname += Name.Data();
208  normname += Name.Data();
209  profname += Name.Data();
210  inname += finDirname.Data();
211  outname += finDirname.Data();
212  cutsname += finDirname.Data();
213  normname += finDirname.Data();
214  profname += finDirname.Data();
215 
216  AliAnalysisDataContainer *cinput = mgr->CreateContainer(inname,TChain::Class(),AliAnalysisManager::kInputContainer);
217 
218  TString outputfile = AliAnalysisManager::GetCommonFileName();
219  outputfile += ":PWG3_D2H_DMult_";
220  outputfile += Name.Data();
221  outputfile += finDirname.Data();
222 
223  AliAnalysisDataContainer *coutputCuts = mgr->CreateContainer(cutsname,TList::Class(),AliAnalysisManager::kOutputContainer,outputfile.Data());
224  AliAnalysisDataContainer *coutput = mgr->CreateContainer(outname,TList::Class(),AliAnalysisManager::kOutputContainer,outputfile.Data());
225  AliAnalysisDataContainer *coutputNorm = mgr->CreateContainer(normname,TList::Class(),AliAnalysisManager::kOutputContainer,outputfile.Data());
226  AliAnalysisDataContainer *coutputProf = mgr->CreateContainer(profname,TList::Class(),AliAnalysisManager::kOutputContainer,outputfile.Data());
227 
228  mgr->ConnectInput(dMultTask,0,mgr->GetCommonInputContainer());
229 
230  mgr->ConnectOutput(dMultTask,1,coutput);
231 
232  mgr->ConnectOutput(dMultTask,2,coutputCuts);
233 
234  mgr->ConnectOutput(dMultTask,3,coutputNorm);
235 
236  mgr->ConnectOutput(dMultTask,4,coutputProf);
237 
238  return dMultTask;
239 }
const char * filename
Definition: TestFCM.C:1
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, Int_t year=16)
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