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