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