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