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