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