AliPhysics  b0b429f (b0b429f)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskSEDvsEventShapes.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
3  * *
4  * Author: The ALICE Off-line Project. *
5  * Contributors are mentioned in the code where appropriate. *
6  * *
7  * Permission to use, copy, modify and distribute this software and its *
8  * documentation strictly for non-commercial purposes is hereby granted *
9  * without fee, provided that the above copyright notice appears in all *
10  * copies and that both the copyright notice and this permission notice *
11  * appear in the supporting documentation. The authors make no claims *
12  * about the suitability of this software for any purpose. It is *
13  * provided "as is" without express or implied warranty. *
14  **************************************************************************/
15 
16 /* $Id$ */
17 
18 //*************************************************************************
19 // Class AliAnalysisTaskSEDvsEventShapes
20 // AliAnalysisTaskSE for the D meson vs. Event shape analysis in different mutiplicity window
21 // Authors: Renu Bala, Manoj Bhanudas Jadhav
22 //*************************************************************************
23 
24 #include <TClonesArray.h>
25 #include <TCanvas.h>
26 #include <TList.h>
27 #include <TMath.h>
28 #include <TString.h>
29 #include <TDatabasePDG.h>
30 #include <TH1F.h>
31 #include <TH2F.h>
32 #include <TH3F.h>
33 #include <THnSparse.h>
34 #include <TProfile.h>
35 #include "AliAnalysisManager.h"
36 #include "AliRDHFCuts.h"
37 #include "AliRDHFCutsDplustoKpipi.h"
38 #include "AliRDHFCutsDStartoKpipi.h"
39 #include "AliRDHFCutsD0toKpi.h"
40 #include "AliAODHandler.h"
41 #include "AliAODEvent.h"
42 #include "AliAODVertex.h"
43 #include "AliAODTrack.h"
44 #include "AliAODRecoDecayHF.h"
45 #include "AliAODRecoCascadeHF.h"
46 #include "AliAnalysisVertexingHF.h"
47 #include "AliAnalysisTaskSE.h"
48 #include "AliAnalysisTaskSEDvsEventShapes.h"
49 #include "AliNormalizationCounter.h"
50 #include "AliVertexingHFUtils.h"
51 #include "AliAODVZERO.h"
52 #include "AliESDUtils.h"
53 ClassImp(AliAnalysisTaskSEDvsEventShapes)
54 
55 //________________________________________________________________________
56 AliAnalysisTaskSEDvsEventShapes::AliAnalysisTaskSEDvsEventShapes():
57 AliAnalysisTaskSE(),
58 fOutput(0),
59 fListCuts(0),
60 fOutputCounters(0),
61 fListProfiles(0),
62 fOutputEffCorr(0),
63 fHistNEvents(0),
64 fHistNtrVsZvtx(0),
65 fHistNtrCorrVsZvtx(0),
66 fHistNtrVsnTrackEvWithCand(0),
67 fHistNtrVsSo(0),
68 fHistNtrCorrVsSo(0),
69 fHistNtrVsSpheri(0),
70 fHistNtrCorrVsSpheri(0),
71 fHistNtrVsNchMC(0),
72 fHistNtrCorrVsNchMC(0),
73 fHistNtrVsNchMCPrimary(0),
74 fHistNtrCorrVsNchMCPrimary(0),
75 fHistNtrVsNchMCPhysicalPrimary(0),
76 fHistNtrCorrVsNchMCPhysicalPrimary(0),
77 fHistGenPrimaryParticlesInelGt0(0),
78 fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary(0),
79 fHistNtrCorrPSSel(0),
80 fHistNtrCorrEvSel(0),
81 fHistNtrCorrEvWithCand(0),
82 fHistNtrCorrEvWithD(0),
83 fHistnTrackvsEtavsPhi(0),
84 fHistnTrackvsEtavsPhiEvWithCand(0),
85 fSparseEvtShape(0),
86 fSparseEvtShapewithNoPid(0),
87 fSparseEvtShapePrompt(0),
88 fSparseEvtShapeFeeddown(0),
89 fSparseEvtShapeRecSphero(0),
90 fMCAccGenPrompt(0),
91 fMCAccGenFeeddown(0),
92 fMCRecoPrompt(0),
93 fMCRecoFeeddown(0),
94 fMCRecoBothPromptFD(0),
95 fMCAccGenPromptSpheri(0),
96 fMCAccGenFeeddownSpheri(0),
97 fMCRecoPromptSpheri(0),
98 fMCRecoFeeddownSpheri(0),
99 fMCRecoBothPromptFDSpheri(0),
100 fMCAccGenPromptEvSel(0),
101 fMCAccGenFeeddownEvSel(0),
102 fUpmasslimit(1.965),
103 fLowmasslimit(1.765),
104 fNMassBins(200),
105 fRDCutsAnalysis(0),
106 fCounterC(0),
107 fCounterU(0),
108 fCounterCandidates(0),
109 fDoImpPar(kFALSE),
110 fNImpParBins(400),
111 fLowerImpPar(-2000.),
112 fHigherImpPar(2000.),
113 fReadMC(kFALSE),
114 fMCOption(0),
115 fisPPbData(kFALSE),
116 fUseBit(kTRUE),
117 fSubtractTrackletsFromDau(kFALSE),
118 fCalculateSphericity(kFALSE),
119 fRecomputeSpherocity(kFALSE),
120 fRemoveD0fromDstar(kFALSE),
121 fUseNchWeight(0),
122 fHistoMCNch(0),
123 fHistoMeasNch(0),
124 fUsePtWeight(kFALSE),
125 fWeight(1.),
126 fRefMult(9.26),
127 fPdgMeson(411),
128 fMultiplicityEstimator(kNtrk10),
129 fMCPrimariesEstimator(kEta10),
130 fDoVZER0ParamVertexCorr(1),
131 fFillSoSparseChecks(0),
132 fUseQuarkTag(kTRUE),
133 fFillTrackHisto(kFALSE),
134 fEtaAccCut(0.9),
135 fPtAccCut(0.1),
136 fetaMin(-0.8),
137 fetaMax(0.8),
138 fptMin(0.15),
139 fptMax(10.),
140 fminMult(3),
141 ffiltbit1(256),
142 ffiltbit2(512),
143 fphiStepSizeDeg(0.1)
144 {
145  // Default constructor
146  for(Int_t i=0; i<5; i++) fHistMassPtImpPar[i]=0;
147  for(Int_t i=0; i<4; i++) fMultEstimatorAvg[i]=0;
148 }
149 
150 //________________________________________________________________________
151 AliAnalysisTaskSEDvsEventShapes::AliAnalysisTaskSEDvsEventShapes(const char *name, Int_t pdgMeson,AliRDHFCuts *cuts, Bool_t switchPPb):
152 AliAnalysisTaskSE(name),
153 fOutput(0),
154 fListCuts(0),
155 fOutputCounters(0),
156 fListProfiles(0),
157 fOutputEffCorr(0),
158 fHistNEvents(0),
159 fHistNtrVsZvtx(0),
160 fHistNtrCorrVsZvtx(0),
161 fHistNtrVsnTrackEvWithCand(0),
162 fHistNtrVsSo(0),
163 fHistNtrCorrVsSo(0),
164 fHistNtrVsSpheri(0),
165 fHistNtrCorrVsSpheri(0),
166 fHistNtrVsNchMC(0),
167 fHistNtrCorrVsNchMC(0),
168 fHistNtrVsNchMCPrimary(0),
169 fHistNtrCorrVsNchMCPrimary(0),
170 fHistNtrVsNchMCPhysicalPrimary(0),
171 fHistNtrCorrVsNchMCPhysicalPrimary(0),
172 fHistGenPrimaryParticlesInelGt0(0),
173 fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary(0),
174 fHistNtrCorrPSSel(0),
175 fHistNtrCorrEvSel(0),
176 fHistNtrCorrEvWithCand(0),
177 fHistNtrCorrEvWithD(0),
178 fHistnTrackvsEtavsPhi(0),
179 fHistnTrackvsEtavsPhiEvWithCand(0),
180 fSparseEvtShape(0),
181 fSparseEvtShapewithNoPid(0),
182 fSparseEvtShapePrompt(0),
183 fSparseEvtShapeFeeddown(0),
184 fSparseEvtShapeRecSphero(0),
185 fMCAccGenPrompt(0),
186 fMCAccGenFeeddown(0),
187 fMCRecoPrompt(0),
188 fMCRecoFeeddown(0),
189 fMCRecoBothPromptFD(0),
190 fMCAccGenPromptSpheri(0),
191 fMCAccGenFeeddownSpheri(0),
192 fMCRecoPromptSpheri(0),
193 fMCRecoFeeddownSpheri(0),
194 fMCRecoBothPromptFDSpheri(0),
195 fMCAccGenPromptEvSel(0),
196 fMCAccGenFeeddownEvSel(0),
197 fUpmasslimit(1.965),
198 fLowmasslimit(1.765),
199 fNMassBins(200),
200 fRDCutsAnalysis(cuts),
201 fCounterC(0),
202 fCounterU(0),
203 fCounterCandidates(0),
204 fDoImpPar(kFALSE),
205 fNImpParBins(400),
206 fLowerImpPar(-2000.),
207 fHigherImpPar(2000.),
208 fReadMC(kFALSE),
209 fMCOption(0),
210 fisPPbData(switchPPb),
211 fUseBit(kTRUE),
212 fSubtractTrackletsFromDau(kFALSE),
213 fCalculateSphericity(kFALSE),
214 fRecomputeSpherocity(kFALSE),
215 fRemoveD0fromDstar(kFALSE),
216 fUseNchWeight(0),
217 fHistoMCNch(0),
218 fHistoMeasNch(0),
219 fUsePtWeight(kFALSE),
220 fWeight(1.),
221 fRefMult(9.26),
222 fPdgMeson(pdgMeson),
223 fMultiplicityEstimator(kNtrk10),
224 fMCPrimariesEstimator(kEta10),
225 fDoVZER0ParamVertexCorr(1),
226 fFillSoSparseChecks(0),
227 fUseQuarkTag(kTRUE),
228 fFillTrackHisto(kFALSE),
229 fEtaAccCut(0.9),
230 fPtAccCut(0.1),
231 fetaMin(-0.8),
232 fetaMax(0.8),
233 fptMin(0.15),
234 fptMax(10.),
235 fminMult(3),
236 ffiltbit1(256),
237 ffiltbit2(512),
238 fphiStepSizeDeg(0.1)
239 {
240  //
241  // Standard constructor
242  //
243 
244  for(Int_t i=0; i<5; i++) fHistMassPtImpPar[i]=0;
245  for(Int_t i=0; i<4; i++) fMultEstimatorAvg[i]=0;
246  if(fPdgMeson==413){
247  fNMassBins=200;
248  SetMassLimits(0.12,0.2);
249  }else{
250  fNMassBins=200;
251  SetMassLimits(fPdgMeson,0.1);
252  }
253  // Default constructor
254  // Otput slot #1 writes into a TList container
255  DefineOutput(1,TList::Class()); //My private output
256  // Output slot #2 writes cut to private output
257  DefineOutput(2,TList::Class());
258  // Output slot #3 writes cut to private output
259  DefineOutput(3,TList::Class());
260  // Output slot #4 writes cut to private output
261  DefineOutput(4,TList::Class());
262  // Output slot #5 writes cut to private output
263  DefineOutput(5,TList::Class());
264 }
265 //________________________________________________________________________
266 AliAnalysisTaskSEDvsEventShapes::~AliAnalysisTaskSEDvsEventShapes()
267 {
268  //
269  // Destructor
270  //
271  delete fOutput;
272  delete fHistNEvents;
273  delete fListCuts;
274  delete fListProfiles;
275  delete fOutputEffCorr;
276  delete fRDCutsAnalysis;
277  delete fCounterC;
278  delete fCounterU;
279  delete fCounterCandidates;
280 
281  for(Int_t i=0; i<4; i++) {
282  if (fMultEstimatorAvg[i]) delete fMultEstimatorAvg[i];
283  }
284 
285  for(Int_t i=0; i<5; i++){
286  delete fHistMassPtImpPar[i];
287  }
288  if(fHistoMCNch) delete fHistoMCNch;
289  if(fHistoMeasNch) delete fHistoMeasNch;
290 }
291 
292 //_________________________________________________________________
293 void AliAnalysisTaskSEDvsEventShapes::SetMassLimits(Double_t lowlimit, Double_t uplimit){
294  // set invariant mass limits
295  if(uplimit>lowlimit){
296  fLowmasslimit = lowlimit;
297  fUpmasslimit = uplimit;
298  }else{
299  AliError("Wrong mass limits: upper value should be larger than lower one");
300  }
301 }
302 //_________________________________________________________________
303 void AliAnalysisTaskSEDvsEventShapes::SetMassLimits(Int_t pdg, Double_t range){
304  // set invariant mass limits
305  Double_t mass=TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdg))->Mass();
306  SetMassLimits(mass-range,mass+range);
307 }
308 //________________________________________________________________________
309 void AliAnalysisTaskSEDvsEventShapes::Init(){
310  //
311  // Initialization
312  //
313  printf("AnalysisTaskSEDvsMultiplicity_0::Init() \n");
314 
315  if(fUseNchWeight && !fReadMC){ AliFatal("Nch weights can only be used in MC mode"); return; }
316  if(fUsePtWeight && !fReadMC){ AliFatal("pT weights can only be used in MC mode"); return; }
317  if(fUsePtWeight && fUseNchWeight) { AliInfo("Beware, using at the same time pT and Nch weights, please check"); }
318  if(fUseNchWeight && !fHistoMCNch){ AliFatal("Nch weights can only be used without histogram"); return; }
319  if(fUseNchWeight==1 && !fHistoMeasNch) {//Nch weights
320  if(fisPPbData){ AliFatal("Nch weights can only be used with MC and data histogram in pPb"); return; }
321  else CreateMeasuredNchHisto();
322  }
323  if(fUseNchWeight==2 && !fHistoMeasNch){ AliFatal("Ntrk weights can only be used with MC and data histogram"); return; } //for pp, pPb Ntrk weights
324 
325  fListCuts=new TList();
326  fListCuts->SetOwner();
327  fListCuts->SetName("CutsList");
328 
329  if(fPdgMeson==411){
330  AliRDHFCutsDplustoKpipi* copycut=new AliRDHFCutsDplustoKpipi(*(static_cast<AliRDHFCutsDplustoKpipi*>(fRDCutsAnalysis)));
331  copycut->SetName("AnalysisCutsDplus");
332  fListCuts->Add(copycut);
333  }else if(fPdgMeson==421){
334  AliRDHFCutsD0toKpi* copycut=new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fRDCutsAnalysis)));
335  copycut->SetName("AnalysisCutsDzero");
336  fListCuts->Add(copycut);
337  }else if(fPdgMeson==413){
338  AliRDHFCutsDStartoKpipi* copycut=new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fRDCutsAnalysis)));
339  copycut->SetName("AnalysisCutsDStar");
340  fListCuts->Add(copycut);
341  }
342  if(fHistoMeasNch) fListCuts->Add(fHistoMeasNch);
343  if(fHistoMCNch) fListCuts->Add(fHistoMCNch);
344 
345  PostData(2,fListCuts);
346 
347  fListProfiles = new TList();
348  fListProfiles->SetOwner();
349  TString period[4];
350  Int_t nProfiles=4;
351  if (fisPPbData) {period[0]="LHC13b"; period[1]="LHC13c"; nProfiles = 2;}
352  else {period[0]="LHC10b"; period[1]="LHC10c"; period[2]="LHC10d"; period[3]="LHC10e"; nProfiles = 4;}
353 
354  for(Int_t i=0; i<nProfiles; i++){
355  if(fMultEstimatorAvg[i]){
356  TProfile* hprof=new TProfile(*fMultEstimatorAvg[i]);
357  hprof->SetName(Form("ProfileTrkVsZvtx%s\n",period[i].Data()));
358  fListProfiles->Add(hprof);
359  }
360  }
361 
362  PostData(4,fListProfiles);
363 
364  return;
365 }
366 
367 //________________________________________________________________________
368 void AliAnalysisTaskSEDvsEventShapes::UserCreateOutputObjects()
369 {
370  // Create the output container
371  //
372  if(fDebug > 1) printf("AnalysisTaskSEDvsMultiplicity::UserCreateOutputObjects() \n");
373 
374  // Several histograms are more conveniently managed in a TList
375  fOutput = new TList();
376  fOutput->SetOwner();
377  fOutput->SetName("OutputHistos");
378 
379  fOutputEffCorr = new TList();
380  fOutputEffCorr->SetOwner();
381  fOutputEffCorr->SetName("OutputEffCorrHistos");
382 
383  Int_t nMultBins = 200;
384  Float_t firstMultBin = -0.5;
385  Float_t lastMultBin = 199.5;
386  Int_t nMultBinsNtrk = nMultBins;
387  Float_t lastMultBinNtrk = lastMultBin;
388  Int_t nMultBinsV0 = 400;
389  Float_t lastMultBinV0 = 799.5;
390  const char *estimatorName="tracklets";
391  if(fisPPbData) {
392  nMultBinsNtrk = 375;
393  lastMultBinNtrk = 374.5;
394  nMultBins = nMultBinsNtrk;
395  lastMultBin = lastMultBinNtrk;
396  }
397  if(fMultiplicityEstimator==kVZERO || fMultiplicityEstimator==kVZEROA || fMultiplicityEstimator==kVZEROEq || fMultiplicityEstimator==kVZEROAEq) {
398  nMultBins = nMultBinsV0;
399  lastMultBin = lastMultBinV0;
400  estimatorName = "vzero";
401  }
402 
403  fHistNtrCorrPSSel = new TH1F("hNtrCorrPSSel",Form("Corrected %s multiplicity for PS selected events; %s ; Entries",estimatorName,estimatorName),nMultBins,firstMultBin,lastMultBin);
404  fHistNtrCorrEvSel = new TH1F("hNtrCorrEvSel",Form("Corrected %s multiplicity for selected events; %s ; Entries",estimatorName,estimatorName),nMultBins,firstMultBin,lastMultBin);
405  fHistNtrCorrEvWithCand = new TH1F("hNtrCorrEvWithCand", Form("%s multiplicity for events with D candidates; %s ; Entries",estimatorName,estimatorName),nMultBins,firstMultBin,lastMultBin);// Total multiplicity
406  fHistNtrCorrEvWithD = new TH1F("hNtrCorrEvWithD", Form("%s multiplicity for events with D in mass region ; %s ; Entries",estimatorName,estimatorName),nMultBins,firstMultBin,lastMultBin); //
407 
408  fHistNtrVsZvtx = new TH2F("hNtrVsZvtx",Form("N%s vs VtxZ; VtxZ;N_{%s};",estimatorName,estimatorName),300,-15,15,nMultBins,firstMultBin,lastMultBin); //
409  fHistNtrCorrVsZvtx = new TH2F("hNtrCorrVsZvtx",Form("N%s vs VtxZ; VtxZ;N_{%s};",estimatorName,estimatorName),300,-15,15,nMultBins,firstMultBin,lastMultBin); //
410  fHistNtrVsnTrackEvWithCand = new TH2F("hNtrVsnTrackEvWithCand",Form("N%s vs nTracks; nTracks; N_{%s};",estimatorName,estimatorName),nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); //
411 
412  TString histoNtrName;
413  TString histoNtrCorrName;
414  TString parNameNtr;
415 
416  histoNtrName = "hNtrVsSphero";
417  histoNtrCorrName = "hNtrCorrVsSphero";
418  parNameNtr = "Sphero";
419 
420  fHistNtrVsSo = new TH2F(histoNtrName.Data(),Form("N_{%s} vs %s; %s; N_{%s};",estimatorName,parNameNtr.Data(),parNameNtr.Data(),estimatorName), 20, 0., 1., nMultBins,firstMultBin,lastMultBin); //
421  fHistNtrCorrVsSo = new TH2F(histoNtrCorrName.Data(),Form("N_{%s} vs %s; %s; N_{%s};",estimatorName,parNameNtr.Data(),parNameNtr.Data(),estimatorName), 20, 0., 1., nMultBins, firstMultBin,lastMultBin); //
422 
423  fHistnTrackvsEtavsPhi = new TH3F("hnTrackvsEtavsPhi", "Eta vs Phi vs nTracks; #eta; #phi[rad]; nTracks;", 100, -1.5, 1.5, 100, 0., 6.28,nMultBins,firstMultBin,lastMultBin);
424  fHistnTrackvsEtavsPhiEvWithCand = new TH3F("hnTrackvsEtavsPhiEvWithCand", "Eta vs Phi vs nTracks; #eta; #phi[rad]; nTracks;", 100, -1.5, 1.5, 100, 0., 6.28,nMultBins,firstMultBin,lastMultBin);
425 
426  TString histoNtrSphriName;
427  TString histoNtrCorrSphriName;
428  TString parNameNtrSphri;
429 
430  histoNtrSphriName = "hNtrVsSpheri";
431  histoNtrCorrSphriName = "hNtrCorrVsSpheri";
432  parNameNtrSphri = "Spheri";
433 
434  if(fCalculateSphericity){
435  fHistNtrVsSpheri = new TH2F(histoNtrSphriName.Data(),Form("N_{%s} vs %s; %s; N_{%s};",estimatorName,parNameNtrSphri.Data(),parNameNtrSphri.Data(),estimatorName), 20, 0., 1., nMultBins,firstMultBin,lastMultBin); //
436  fHistNtrCorrVsSpheri = new TH2F(histoNtrCorrSphriName.Data(),Form("N_{%s} vs %s; %s; N_{%s};",estimatorName,parNameNtrSphri.Data(),parNameNtrSphri.Data(),estimatorName), 20, 0., 1., nMultBins, firstMultBin,lastMultBin); //
437  }
438 
439  fHistNtrVsNchMC = new TH2F("hNtrVsNchMC",Form("N%s vs NchMC; Nch;N_{%s};",estimatorName,estimatorName),nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); //
440  fHistNtrCorrVsNchMC = new TH2F("hNtrCorrVsNchMC",Form("N%s vs Nch; Nch;N_{%s};",estimatorName,estimatorName),nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); //
441 
442  fHistNtrVsNchMCPrimary = new TH2F("hNtrVsNchMCPrimary",Form("N%s vs Nch (Primary); Nch (Primary);N_{%s};",estimatorName,estimatorName),nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); //
443  fHistNtrCorrVsNchMCPrimary = new TH2F("hNtrCorrVsNchMCPrimary",Form("N%s vs Nch (Primary); Nch(Primary) ;N_{%s};",estimatorName,estimatorName),nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); //
444 
445  fHistNtrVsNchMCPhysicalPrimary = new TH2F("hNtrVsNchMCPhysicalPrimary",Form("N%s vs Nch (Physical Primary); Nch (Physical Primary);N_{%s};",estimatorName,estimatorName),nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); //
446  fHistNtrCorrVsNchMCPhysicalPrimary = new TH2F("hNtrCorrVsMCPhysicalPrimary",Form("N%s vs Nch (Physical Primary); Nch (Physical Primary);N_{%s};",estimatorName,estimatorName),nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); //
447 
448  fHistGenPrimaryParticlesInelGt0 = new TH1F("hGenPrimaryParticlesInelGt0","Multiplcity of generated charged particles ; Nparticles ; Entries",nMultBins,firstMultBin,lastMultBin);
449 
450  fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary = new TH3F("fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary", "MC: Nch (Physical Primary) vs Nch (Primary) vs Nch (Generated); Nch (Generated); Nch (Primary); Nch (Physical Primary)",nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin);
451 
452  fOutput->Add(fHistNtrCorrPSSel);
453  fOutput->Add(fHistNtrCorrEvSel);
454  fOutput->Add(fHistNtrCorrEvWithCand);
455  fOutput->Add(fHistNtrCorrEvWithD);
456 
457  fOutput->Add(fHistNtrVsZvtx);
458  fOutput->Add(fHistNtrCorrVsZvtx);
459  fOutput->Add(fHistNtrVsnTrackEvWithCand);
460  fOutput->Add(fHistnTrackvsEtavsPhi);
461  fOutput->Add(fHistnTrackvsEtavsPhiEvWithCand);
462 
463  fOutput->Add(fHistNtrVsSo);
464  fOutput->Add(fHistNtrCorrVsSo);
465  if(fCalculateSphericity){
466  fOutput->Add(fHistNtrVsSpheri);
467  fOutput->Add(fHistNtrCorrVsSpheri);
468  }
469  fOutput->Add(fHistNtrVsNchMC);
470  fOutput->Add(fHistNtrCorrVsNchMC);
471  fOutput->Add(fHistNtrVsNchMCPrimary);
472  fOutput->Add(fHistNtrCorrVsNchMCPrimary);
473  fOutput->Add(fHistNtrVsNchMCPhysicalPrimary);
474  fOutput->Add(fHistNtrCorrVsNchMCPhysicalPrimary);
475  fOutput->Add(fHistGenPrimaryParticlesInelGt0);
476  fOutput->Add(fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary);
477 
478  fHistNEvents = new TH1F("fHistNEvents", "number of events ",11,-0.5,10.5);
479  fHistNEvents->GetXaxis()->SetBinLabel(1,"nEvents total");
480  fHistNEvents->GetXaxis()->SetBinLabel(2,"nEvents with Z vertex");
481  fHistNEvents->GetXaxis()->SetBinLabel(3,"nEvents selected");
482  fHistNEvents->GetXaxis()->SetBinLabel(4,"Rejected due to trigger");
483  fHistNEvents->GetXaxis()->SetBinLabel(5,"Rejected due to phys sel");
484  fHistNEvents->GetXaxis()->SetBinLabel(6,"Rejected due to vertex cuts");
485  fHistNEvents->GetXaxis()->SetBinLabel(7,"Rejected due to pileup");
486  fHistNEvents->GetXaxis()->SetBinLabel(8,"Total no. of candidate");
487  fHistNEvents->GetXaxis()->SetBinLabel(9,"no. of cand wo bitmask");
488  fHistNEvents->GetXaxis()->SetBinLabel(10,"D after cuts (No PID)");
489  fHistNEvents->GetXaxis()->SetBinLabel(11,"D after cuts + PID)");
490  fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);
491  fHistNEvents->SetMinimum(0);
492  fOutput->Add(fHistNEvents);
493 
494  // With flag fFillSoSparseChecks to fill THnSparse with MultUncorr and NoPid cases ( 0 = only Mult, 1 = Mult and multUncorr, 2 = NoPid and 3 is All)
495  Int_t nbinsSo[4]={48, fNMassBins, 20, nMultBins};
496  Double_t xminSo[4]={0., fLowmasslimit,0., firstMultBin};
497  Double_t xmaxSo[4]={24., fUpmasslimit, 1., lastMultBin};
498 
499  Int_t nbinsSoSpheri[5]={48, fNMassBins, 20, nMultBins, 20};
500  Double_t xminSoSpheri[5]={0., fLowmasslimit,0., firstMultBin, 0.};
501  Double_t xmaxSoSpheri[5]={24., fUpmasslimit, 1., lastMultBin, 1.};
502 
503  Int_t nbinsSowithMultUncorr[5]={48, fNMassBins, 20, nMultBins, nMultBins};
504  Double_t xminSowithMultUncorr[5]={0., fLowmasslimit,0., firstMultBin, firstMultBin};
505  Double_t xmaxSowithMultUncorr[5]={24., fUpmasslimit, 1., lastMultBin, lastMultBin};
506 
507  Int_t nbinsSoSpheriwithMultUncorr[6]={48, fNMassBins, 20, nMultBins, nMultBins, 20};
508  Double_t xminSoSpheriwithMultUncorr[6]={0., fLowmasslimit,0., firstMultBin, firstMultBin, 0.};
509  Double_t xmaxSoSpheriwithMultUncorr[6]={24., fUpmasslimit, 1., lastMultBin, lastMultBin, 1.};
510 
511  TString histoName = "hSparseEvtShape";
512  TString histoNameNoPid = "hSparseEvtShapewithNoPid";
513  TString parNameSo = "Spherocity";
514  TString parNameSpheri = "Sphericity";
515  TString histoNamePrompt = "hSparseEvtShapePrompt";
516  TString histoNameFeeddown = "hSparseEvtShapeFeeddown";
517  TString histoNameRecSphero = "hSparseEvtShapeRecSphero";
518 
519  if(fFillSoSparseChecks == 1 || fFillSoSparseChecks == 3){
520  if(fCalculateSphericity) fSparseEvtShape = new THnSparseD(histoName.Data(), Form("D candidates:; p_{T} [GeV/c]; InvMass [GeV/c^{2}]; %s; Multipicity; MultipicityUncorr; %s;", parNameSo.Data(), parNameSpheri.Data()), 6 , nbinsSoSpheriwithMultUncorr, xminSoSpheriwithMultUncorr, xmaxSoSpheriwithMultUncorr);
521  else fSparseEvtShape = new THnSparseD(histoName.Data(), Form("D candidates:; p_{T} [GeV/c]; InvMass [GeV/c^{2}]; %s; Multipicity; MultipicityUncorr;", parNameSo.Data()), 5 , nbinsSowithMultUncorr, xminSowithMultUncorr, xmaxSowithMultUncorr);
522  }
523  else{
524  if(fCalculateSphericity) fSparseEvtShape = new THnSparseD(histoName.Data(), Form("D candidates:; p_{T} [GeV/c]; InvMass [GeV/c^{2}]; %s; Multipicity; %s;", parNameSo.Data(), parNameSpheri.Data()), 5 , nbinsSoSpheri, xminSoSpheri, xmaxSoSpheri);
525  else fSparseEvtShape = new THnSparseD(histoName.Data(), Form("D candidates:; p_{T} [GeV/c]; InvMass [GeV/c^{2}]; %s; Multipicity;", parNameSo.Data()), 4 , nbinsSo, xminSo, xmaxSo);
526  }
527 
528  if(fFillSoSparseChecks == 2|| fFillSoSparseChecks == 3) {
529  if(fCalculateSphericity) fSparseEvtShapewithNoPid = new THnSparseD(histoNameNoPid.Data(), Form("D candidates with NoPID:; p_{T} [GeV/c]; InvMass [GeV/c^{2}]; %s; Multipicity; %s;", parNameSo.Data(), parNameSpheri.Data()), 5 , nbinsSoSpheri, xminSoSpheri, xmaxSoSpheri);
530  else fSparseEvtShapewithNoPid = new THnSparseD(histoNameNoPid.Data(), Form("D candidates with NoPID:; p_{T} [GeV/c]; InvMass [GeV/c^{2}]; %s; Multipicity;", parNameSo.Data()), 4 , nbinsSo, xminSo, xmaxSo);
531  }
532 
533  if(fRecomputeSpherocity) fSparseEvtShapeRecSphero = new THnSparseD(histoNameRecSphero.Data(), Form("D candidates:; p_{T} [GeV/c]; InvMass [GeV/c^{2}]; %s; Multipicity; RecSpherocity;", parNameSo.Data()), 5 , nbinsSoSpheri, xminSoSpheri, xmaxSoSpheri);
534 
535  fOutput->Add(fSparseEvtShape);
536  if(fFillSoSparseChecks == 2 || fFillSoSparseChecks == 3) fOutput->Add(fSparseEvtShapewithNoPid);
537  if(fRecomputeSpherocity) fOutput->Add(fSparseEvtShapeRecSphero);
538 
539  Int_t nbinsPrompt[4]={48, nMultBins, 20, 100};
540  Int_t nbinsFeeddown[4]={48, nMultBins, 20, 100};
541  Double_t xminPrompt[4] = {0.,firstMultBin, 0., -1.};
542  Double_t xmaxPrompt[4] = {24.,lastMultBin, 1., 1.};
543  Double_t xminFeeddown[4] = {0.,firstMultBin, 0., -1.};
544  Double_t xmaxFeeddown[4] = {24.,lastMultBin, 1., 1.};
545 
546  Int_t nbinsRecSpheroPrompt[5]={48, nMultBins, 20, 100, 20};
547  Int_t nbinsRecSpheroFeeddown[5]={48, nMultBins, 20, 100, 20};
548  Double_t xminRecSpheroPrompt[5] = {0.,firstMultBin, 0., -1., 0.};
549  Double_t xmaxRecSpheroPrompt[5] = {24.,lastMultBin, 1., 1., 1.};
550  Double_t xminRecSpheroFeeddown[5] = {0.,firstMultBin, 0., -1., 0.};
551  Double_t xmaxRecSpheroFeeddown[5] = {24.,lastMultBin, 1., 1., 1.};
552 
553  if(fReadMC){
554  if(fRecomputeSpherocity){
555  fSparseEvtShapePrompt = new THnSparseD(histoNamePrompt.Data(), Form("D candidates:; p_{T} [GeV/c]; InvMass [GeV/c^{2}]; %s; Multipicity; RecSpherocity;", parNameSo.Data()), 5 , nbinsSoSpheri, xminSoSpheri, xmaxSoSpheri);
556  fSparseEvtShapeFeeddown = new THnSparseD(histoNameFeeddown.Data(), Form("D candidates:; p_{T} [GeV/c]; InvMass [GeV/c^{2}]; %s; Multipicity; RecSpherocity;", parNameSo.Data()), 5 , nbinsSoSpheri, xminSoSpheri, xmaxSoSpheri);
557  }
558  else{
559  fSparseEvtShapePrompt = new THnSparseD(histoNamePrompt.Data(), Form("D candidates:; p_{T} [GeV/c]; InvMass [GeV/c^{2}]; %s; Multipicity;", parNameSo.Data()), 4 , nbinsSo, xminSo, xmaxSo);
560  fSparseEvtShapeFeeddown = new THnSparseD(histoNameFeeddown.Data(), Form("D candidates:; p_{T} [GeV/c]; InvMass [GeV/c^{2}]; %s; Multipicity;", parNameSo.Data()), 4 , nbinsSo, xminSo, xmaxSo);
561  }
562 
563  fOutput->Add(fSparseEvtShapePrompt);
564  fOutput->Add(fSparseEvtShapeFeeddown);
565 
566  //Prompt
567  if(fRecomputeSpherocity){
568  fMCAccGenPrompt = new THnSparseD("hMCAccGenPrompt", "kStepMCAcceptance:; p_{T} [GeV/c]; Multipicity; Spherocity; y; RecSpherocity; - promptD",5,nbinsRecSpheroPrompt,xminRecSpheroPrompt,xmaxRecSpheroPrompt);
569  fMCRecoPrompt = new THnSparseD("hMCRecoPrompt","kStepRecoPID:; p_{T} [GeV/c]; Multipicity; Spherocity; y; RecSpherocity; - promptD",5,nbinsRecSpheroPrompt,xminRecSpheroPrompt,xmaxRecSpheroPrompt);
570  }
571  else{
572  fMCAccGenPrompt = new THnSparseD("hMCAccGenPrompt","kStepMCAcceptance:; p_{T} [GeV/c]; Multipicity; Spherocity; y; - promptD",4,nbinsPrompt,xminPrompt,xmaxPrompt);
573  fMCRecoPrompt = new THnSparseD("hMCRecoPrompt","kStepRecoPID:; p_{T} [GeV/c]; Multipicity; Spherocity; y; - promptD",4,nbinsPrompt,xminPrompt,xmaxPrompt);
574  }
575  if(fCalculateSphericity){
576  fMCAccGenPromptSpheri = new THnSparseD("hMCAccGenPromptSpheri","kStepMCAcceptance:; p_{T} [GeV/c]; Multipicity; Sphericity; y; - promptD",4,nbinsPrompt,xminPrompt,xmaxPrompt);
577  fMCRecoPromptSpheri = new THnSparseD("hMCRecoPromptSpheri","kStepRecoPID:; p_{T} [GeV/c]; Multipicity; Sphericity; y; - promptD",4,nbinsPrompt,xminPrompt,xmaxPrompt);
578  }
579  fMCAccGenPromptEvSel = new THnSparseD("hMCAccGenPromptEvSel","kStepMCAcceptanceEvSel:; p_{T} [GeV/c]; Multipicity; Spherocity; y; - promptD",4,nbinsPrompt,xminPrompt,xmaxPrompt);
580 
581  //Feeddown
582  if(fRecomputeSpherocity){
583  fMCAccGenFeeddown = new THnSparseD("hMCAccGenBFeeddown","kStepMCAcceptance:; p_{T} [GeV/c]; Multipicity; Spherocity; y; RecSpherocity; - DfromB",5,nbinsRecSpheroFeeddown,xminRecSpheroFeeddown,xmaxRecSpheroFeeddown);
584  fMCRecoFeeddown = new THnSparseD("hMCRecoFeeddown","kStepRecoPID:; p_{T} [GeV/c]; Multipicity; Spherocity; y; RecSpherocity; - DfromB",5,nbinsRecSpheroFeeddown,xminRecSpheroFeeddown,xmaxRecSpheroFeeddown);
585  }
586  else{
587  fMCAccGenFeeddown = new THnSparseD("hMCAccGenBFeeddown","kStepMCAcceptance:; p_{T} [GeV/c]; Multipicity; Spherocity; y; - DfromB",4,nbinsFeeddown,xminFeeddown,xmaxFeeddown);
588  fMCRecoFeeddown = new THnSparseD("hMCRecoFeeddown","kStepRecoPID:; p_{T} [GeV/c]; Multipicity; Spherocity; y; - DfromB",4,nbinsFeeddown,xminFeeddown,xmaxFeeddown);
589  }
590  if(fCalculateSphericity){
591  fMCAccGenFeeddownSpheri = new THnSparseD("hMCAccGenBFeeddownSpheri","kStepMCAcceptance:; p_{T} [GeV/c]; Multipicity; Sphericity; y; - DfromB",4,nbinsFeeddown,xminFeeddown,xmaxFeeddown);
592  fMCRecoFeeddownSpheri = new THnSparseD("hMCRecoFeeddownSpheri","kStepRecoPID:; p_{T} [GeV/c]; Multipicity; Sphericity; y; - DfromB",4,nbinsFeeddown,xminFeeddown,xmaxFeeddown);
593  }
594  fMCAccGenFeeddownEvSel = new THnSparseD("hMCAccGenBFeeddownEvSel","kStepMCAcceptance:; p_{T} [GeV/c]; Multipicity; Spherocity; y; - DfromB",4,nbinsFeeddown,xminFeeddown,xmaxFeeddown);
595 
596  //BothPromptFeeddown
597  fMCRecoBothPromptFD = new THnSparseD("hMCRecoBothPromptFD","kStepRecoPID:; p_{T} [GeV/c]; Multipicity; Spherocity; y; - BothPromptFD",4,nbinsPrompt,xminPrompt,xmaxPrompt);
598 
599  if(fCalculateSphericity) fMCRecoBothPromptFDSpheri = new THnSparseD("hMCRecoBothPromptFDSpheri","kStepRecoPID:; p_{T} [GeV/c]; Multipicity; Sphericity; y; - BothPromptFD",4,nbinsPrompt,xminPrompt,xmaxPrompt);
600 
601  fOutputEffCorr->Add(fMCAccGenPrompt);
602  fOutputEffCorr->Add(fMCAccGenFeeddown);
603  fOutputEffCorr->Add(fMCRecoPrompt);
604  fOutputEffCorr->Add(fMCRecoFeeddown);
605  fOutputEffCorr->Add(fMCRecoBothPromptFD);
606  if(fCalculateSphericity){
607  fOutputEffCorr->Add(fMCAccGenPromptSpheri);
608  fOutputEffCorr->Add(fMCAccGenFeeddownSpheri);
609  fOutputEffCorr->Add(fMCRecoPromptSpheri);
610  fOutputEffCorr->Add(fMCRecoFeeddownSpheri);
611  fOutputEffCorr->Add(fMCRecoBothPromptFDSpheri);
612  }
613  fOutputEffCorr->Add(fMCAccGenPromptEvSel);
614  fOutputEffCorr->Add(fMCAccGenFeeddownEvSel);
615  }
616 
617  if(fDoImpPar) CreateImpactParameterHistos();
618 
619  fCounterC = new AliNormalizationCounter("NormCounterCorrMult");
620  fCounterC->SetStudyMultiplicity(kTRUE,1.);
621  fCounterC->SetStudySpherocity(kTRUE,10.);
622  fCounterC->Init();
623 
624  fCounterU = new AliNormalizationCounter("NormCounterUnCorrMult");
625  fCounterU->SetStudyMultiplicity(kTRUE,1.);
626  fCounterU->SetStudySpherocity(kTRUE,10.);
627  fCounterU->Init();
628 
629  fCounterCandidates = new AliNormalizationCounter("NormCounterCorrMultCandidates");
630  fCounterCandidates->SetStudyMultiplicity(kTRUE,1.);
631  fCounterCandidates->SetStudySpherocity(kTRUE,10.);
632  fCounterCandidates->Init();
633 
634  fOutputCounters = new TList();
635  fOutputCounters->SetOwner();
636  fOutputCounters->SetName("OutputCounters");
637  fOutputCounters->Add(fCounterC);
638  fOutputCounters->Add(fCounterU);
639  fOutputCounters->Add(fCounterCandidates);
640 
641  PostData(1,fOutput);
642  PostData(2,fListCuts);
643  PostData(3,fOutputCounters);
644  PostData(4,fListProfiles);
645  if(fReadMC) PostData(5,fOutputEffCorr);
646 
647  return;
648 }
649 
650 //________________________________________________________________________
651 void AliAnalysisTaskSEDvsEventShapes::UserExec(Option_t */*option*/)
652 {
653  // Execute analysis for current event:
654  // heavy flavor candidates association to MC truth
655 
656  AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
657 
658  TClonesArray *arrayCand = 0;
659  TString arrayName="";
660  UInt_t pdgDau[3];
661  Int_t nDau=0;
662  Int_t selbit=0;
663  if(fPdgMeson==411){
664  arrayName="Charm3Prong";
665  pdgDau[0]=211; pdgDau[1]=321; pdgDau[2]=211;
666  nDau=3;
667  selbit=AliRDHFCuts::kDplusCuts;
668  }else if(fPdgMeson==421){
669  arrayName="D0toKpi";
670  pdgDau[0]=211; pdgDau[1]=321; pdgDau[2]=0;
671  nDau=2;
672  selbit=AliRDHFCuts::kD0toKpiCuts;
673  }else if(fPdgMeson==413){
674  arrayName="Dstar";
675  pdgDau[0]=321; pdgDau[1]=211; pdgDau[2]=0; // Quoting here D0 daughters (D* ones on another variable later)
676  nDau=2;
677  selbit=AliRDHFCuts::kDstarCuts;
678  }
679 
680  if(!aod && AODEvent() && IsStandardAOD()) {
681  // In case there is an AOD handler writing a standard AOD, use the AOD
682  // event in memory rather than the input (ESD) event.
683  aod = dynamic_cast<AliAODEvent*> (AODEvent());
684  // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
685  // have to taken from the AOD event hold by the AliAODExtension
686  AliAODHandler* aodHandler = (AliAODHandler*)
687  ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
688  if(aodHandler->GetExtensions()) {
689  AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
690  AliAODEvent *aodFromExt = ext->GetAOD();
691  arrayCand=(TClonesArray*)aodFromExt->GetList()->FindObject(arrayName.Data());
692  }
693  } else if(aod) {
694  arrayCand=(TClonesArray*)aod->GetList()->FindObject(arrayName.Data());
695  }
696 
697  if(!aod || !arrayCand) {
698  printf("AliAnalysisTaskSEDvsEventShapes::UserExec: Charm3Prong branch not found!\n");
699  return;
700  }
701 
702  if(fisPPbData && fReadMC){
703  Int_t runnumber = aod->GetRunNumber();
704  if(aod->GetTriggerMask()==0 &&
705  (runnumber>=195344 && runnumber<=195677)){
706  AliDebug(3,"Event rejected because of null trigger mask");
707  return;
708  }
709  }
710 
711  // fix for temporary bug in ESDfilter
712  // the AODs with null vertex pointer didn't pass the PhysSel
713  if(!aod->GetPrimaryVertex()||TMath::Abs(aod->GetMagneticField())<0.001) return;
714 
715  Int_t countTreta1=0, countTreta03=0, countTreta05=0, countTreta16=0;
716  AliAODTracklets* tracklets=aod->GetTracklets();
717  Int_t nTr=tracklets->GetNumberOfTracklets();
718  for(Int_t iTr=0; iTr<nTr; iTr++){
719  Double_t theta=tracklets->GetTheta(iTr);
720  Double_t eta=-TMath::Log(TMath::Tan(theta/2.));
721  if(eta>-0.3 && eta<0.3) countTreta03++;
722  if(eta>-0.5 && eta<0.5) countTreta05++;
723  if(eta>-1.0 && eta<1.0) countTreta1++;
724  if(eta>-1.6 && eta<1.6) countTreta16++;
725  }
726 
727  Int_t vzeroMult=0, vzeroMultA=0, vzeroMultC=0;
728  Int_t vzeroMultEq=0, vzeroMultAEq=0, vzeroMultCEq=0;
729  AliAODVZERO *vzeroAOD = (AliAODVZERO*)aod->GetVZEROData();
730  if(vzeroAOD) {
731  vzeroMultA = static_cast<Int_t>(vzeroAOD->GetMTotV0A());
732  vzeroMultC = static_cast<Int_t>(vzeroAOD->GetMTotV0C());
733  vzeroMult = vzeroMultA + vzeroMultC;
734  vzeroMultAEq = static_cast<Int_t>(AliVertexingHFUtils::GetVZEROAEqualizedMultiplicity(aod));
735  vzeroMultCEq = static_cast<Int_t>(AliVertexingHFUtils::GetVZEROCEqualizedMultiplicity(aod));
736  vzeroMultEq = vzeroMultAEq + vzeroMultCEq;
737  }
738 
739  Int_t countMult = countTreta1;
740  if(fMultiplicityEstimator==kNtrk03) { countMult = countTreta03; }
741  else if(fMultiplicityEstimator==kNtrk05) { countMult = countTreta05; }
742  else if(fMultiplicityEstimator==kNtrk10to16) { countMult = countTreta16 - countTreta1; }
743  else if(fMultiplicityEstimator==kVZERO) { countMult = vzeroMult; }
744  else if(fMultiplicityEstimator==kVZEROA) { countMult = vzeroMultA; }
745  else if(fMultiplicityEstimator==kVZEROEq) { countMult = vzeroMultEq; }
746  else if(fMultiplicityEstimator==kVZEROAEq) { countMult = vzeroMultAEq; }
747 
748  Double_t spherocity;
749  Double_t sphericity = -0.5;
750  if(fCalculateSphericity){ //When kTRUE, it calculates Sphericity and THnSparse filled for sphericity
751  sphericity=AliVertexingHFUtils::GetSphericity(aod, fetaMin, fetaMax, fptMin, fptMax, ffiltbit1, ffiltbit2, fminMult);
752  }
753  spherocity=AliVertexingHFUtils::GetSpherocity(aod, fetaMin, fetaMax, fptMin, fptMax, ffiltbit1, ffiltbit2, fminMult, fphiStepSizeDeg);
754 
755  Double_t St=1;
756  fCounterU->StoreEvent(aod,fRDCutsAnalysis,fReadMC,countMult,spherocity);
757  fHistNEvents->Fill(0); // count event
758 
759  Double_t countTreta1corr=countTreta1;
760  Double_t countCorr=countMult;
761  AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
762  // In case of VZERO multiplicity, consider the zvtx correction flag
763  // fDoVZER0ParamVertexCorr: 0= none, 1= usual d2h, 2=AliESDUtils
764  Bool_t isDataDrivenZvtxCorr=kTRUE;
765  Bool_t isVtxOk=kFALSE;
766  Int_t vzeroMultACorr=vzeroMultA, vzeroMultCCorr=vzeroMultC, vzeroMultCorr=vzeroMult;
767  Int_t vzeroMultAEqCorr=vzeroMultAEq, vzeroMultCEqCorr=vzeroMultCEq, vzeroMultEqCorr=vzeroMultEq;
768  if(vtx1){
769  if(vtx1->GetNContributors()>0){
770  fHistNEvents->Fill(1);
771  isVtxOk=kTRUE;
772  }
773  }
774  if(isVtxOk){
775  if( (fMultiplicityEstimator==kVZERO) || (fMultiplicityEstimator==kVZEROA) ||
776  (fMultiplicityEstimator==kVZEROEq) || (fMultiplicityEstimator==kVZEROAEq) ){
777  if(fDoVZER0ParamVertexCorr==0){
778  // do not correct
779  isDataDrivenZvtxCorr=kFALSE;
780  } else if (fDoVZER0ParamVertexCorr==2){
781  // use AliESDUtils correction
782  Float_t zvtx = vtx1->GetZ();
783  isDataDrivenZvtxCorr=kFALSE;
784  vzeroMultACorr = static_cast<Int_t>(AliESDUtils::GetCorrV0A(vzeroMultA,zvtx));
785  vzeroMultCCorr = static_cast<Int_t>(AliESDUtils::GetCorrV0C(vzeroMultC,zvtx));
786  vzeroMultCorr = vzeroMultACorr + vzeroMultCCorr;
787  vzeroMultAEqCorr = static_cast<Int_t>(AliESDUtils::GetCorrV0A(vzeroMultAEq,zvtx));
788  vzeroMultCEqCorr =static_cast<Int_t>( AliESDUtils::GetCorrV0C(vzeroMultCEq,zvtx));
789  vzeroMultEqCorr = vzeroMultAEqCorr + vzeroMultCEqCorr;
790  if(fMultiplicityEstimator==kVZERO) { countCorr = vzeroMultCorr; }
791  else if(fMultiplicityEstimator==kVZEROA) { countCorr = vzeroMultACorr; }
792  else if(fMultiplicityEstimator==kVZEROEq) { countCorr = vzeroMultEqCorr; }
793  else if(fMultiplicityEstimator==kVZEROAEq) { countCorr = vzeroMultAEqCorr; }
794  }
795  }
796  }
797  // Data driven multiplicity z-vertex correction
798  if(isVtxOk && isDataDrivenZvtxCorr){
799  TProfile* estimatorAvg = GetEstimatorHistogram(aod);
800  if(estimatorAvg){
801  countTreta1corr=static_cast<Int_t>(AliVertexingHFUtils::GetCorrectedNtracklets(estimatorAvg,countTreta1,vtx1->GetZ(),fRefMult));
802  countCorr=static_cast<Int_t>(AliVertexingHFUtils::GetCorrectedNtracklets(estimatorAvg,countMult,vtx1->GetZ(),fRefMult));
803  }
804  }
805 
806  fCounterC->StoreEvent(aod,fRDCutsAnalysis,fReadMC,countCorr,spherocity);
807 
808  Bool_t isEvSel=fRDCutsAnalysis->IsEventSelected(aod);
809 
810  if(fRDCutsAnalysis->GetWhyRejection()==5) fHistNEvents->Fill(3);
811  if(fRDCutsAnalysis->GetWhyRejection()==7) fHistNEvents->Fill(4);
812  if(fRDCutsAnalysis->GetWhyRejection()==6) fHistNEvents->Fill(5);
813  if(fRDCutsAnalysis->GetWhyRejection()==1) fHistNEvents->Fill(6);
814 
815  Bool_t isEvPSRejected = fRDCutsAnalysis->IsEventRejectedDuePhysicsSelection();
816 
817  if(!isEvPSRejected){
818  fHistNtrCorrPSSel->Fill(countCorr);
819  }
820 
821  TClonesArray *arrayMC=0;
822  AliAODMCHeader *mcHeader=0;
823 
824  fWeight=1.;
825  Double_t nchWeight=1.0;
826 
827  Int_t nChargedMCEta10=0, nChargedMCEta03=0, nChargedMCEta05=0, nChargedMCEta16=0, nChargedMCEtam37tm17=0, nChargedMCEta28t51=0;
828  Int_t nChargedMCPrimaryEta10=0, nChargedMCPrimaryEta03=0, nChargedMCPrimaryEta05=0, nChargedMCPrimaryEta16=0, nChargedMCPrimaryEtam37tm17=0, nChargedMCPrimaryEta28t51=0;
829  Int_t nChargedMCPhysicalPrimaryEta10=0, nChargedMCPhysicalPrimaryEta03=0, nChargedMCPhysicalPrimaryEta05=0, nChargedMCPhysicalPrimaryEta16=0, nChargedMCPhysicalPrimaryEtam37tm17=0, nChargedMCPhysicalPrimaryEta28t51=0;
830  Int_t nChargedMC=0, nChargedMCPrimary=0, nChargedMCPhysicalPrimary=0;
831 
832  // load MC particles and get weight on Nch
833  if(fReadMC){
834 
835  arrayMC = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
836  if(!arrayMC) {
837  printf("AliAnalysisTaskSEDvsEventShapes::UserExec: MC particles branch not found!\n");
838  return;
839  }
840  // load MC header
841  mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
842  if(!mcHeader) {
843  printf("AliAnalysisTaskSEDvsEventShapes::UserExec: MC header branch not found!\n");
844  return;
845  }
846 
847  for(Int_t i=0; i<arrayMC->GetEntriesFast(); i++){
848  AliAODMCParticle *part=(AliAODMCParticle*)arrayMC->UncheckedAt(i);
849  Int_t charge = part->Charge();
850  Double_t eta = part->Eta();
851  Bool_t isPrim = part->IsPrimary();
852  Bool_t isPhysPrim = part->IsPhysicalPrimary();
853  if(charge!=0) {
854  if(eta>-0.3 && eta< 0.3) {
855  nChargedMCEta03++;
856  if(isPrim) nChargedMCPrimaryEta03++;
857  if(isPhysPrim) nChargedMCPhysicalPrimaryEta03++;
858  }
859  if(eta>-0.5 && eta< 0.5) {
860  nChargedMCEta05++;
861  if(isPrim) nChargedMCPrimaryEta05++;
862  if(isPhysPrim) nChargedMCPhysicalPrimaryEta05++;
863  }
864  if(eta>-1.0 && eta< 1.0) {
865  nChargedMCEta10++;
866  if(isPrim) nChargedMCPrimaryEta10++;
867  if(isPhysPrim) nChargedMCPhysicalPrimaryEta10++;
868  }
869  if(eta>-1.6 && eta< 1.6) {
870  nChargedMCEta16++;
871  if(isPrim) nChargedMCPrimaryEta16++;
872  if(isPhysPrim) nChargedMCPhysicalPrimaryEta16++;
873  }
874  if(eta>-3.7 && eta<-1.7) {
875  nChargedMCEtam37tm17++;
876  if(isPrim) nChargedMCPrimaryEtam37tm17++;
877  if(isPhysPrim) nChargedMCPhysicalPrimaryEtam37tm17++;
878  }
879  if(eta> 2.8 && eta< 5.1) {
880  nChargedMCEta28t51++;
881  if(isPrim) nChargedMCPrimaryEta28t51++;
882  if(isPhysPrim) nChargedMCPhysicalPrimaryEta28t51++;
883  }
884  }
885  }
886 
887  nChargedMC=nChargedMCEta10;
888  nChargedMCPrimary=nChargedMCPrimaryEta10;
889  nChargedMCPhysicalPrimary=nChargedMCPhysicalPrimaryEta10;
890 
891  // Compute the Nch weights (reference is Ntracklets within |eta|<1.0)
892  if(fUseNchWeight>0){
893 
894  Double_t tmpweight = 1.0;
895  Double_t tmpXweight=nChargedMCPhysicalPrimary; // Nch weights
896  if(fUseNchWeight==2) tmpXweight=countMult; // Ntrk weights
897 
898  if(tmpXweight<=0) tmpweight = 0.0;
899  else{
900  Double_t pMeas = fHistoMeasNch->GetBinContent(fHistoMeasNch->FindBin(tmpXweight));
901  //printf(" pMeas=%2.2f and histo MCNch %s \n",pMeas,fHistoMCNch);
902  Double_t pMC = fHistoMCNch->GetBinContent(fHistoMCNch->FindBin(tmpXweight));
903  tmpweight = pMC>0 ? pMeas/pMC : 0.;
904  }
905  nchWeight *= tmpweight;
906  AliDebug(2,Form("Using Nch weights, Mult=%f Weight=%f\n",tmpXweight,nchWeight));
907  }
908 
909  FillMCGenAccHistos(aod, arrayMC, mcHeader, countCorr, spherocity, sphericity, isEvSel, nchWeight);//Fill 2 separate THnSparses, one for prompt andf one for feeddown
910  }
911 
912  if(!isEvSel)return;
913 
914  if(vtx1){
915  fHistNtrVsZvtx->Fill(vtx1->GetZ(),countMult);
916  fHistNtrCorrVsZvtx->Fill(vtx1->GetZ(),countCorr);
917  fHistNtrVsSo->Fill(spherocity,countMult);
918  fHistNtrCorrVsSo->Fill(spherocity,countCorr);
919  if(fCalculateSphericity){
920  fHistNtrVsSpheri->Fill(sphericity,countMult);
921  fHistNtrCorrVsSpheri->Fill(sphericity,countCorr);
922  }
923  }
924 
925  if(fReadMC){
926  // Now recompute the variables in case another MC estimator is considered
927  if(fMCPrimariesEstimator==kEta10to16){
928  nChargedMC = nChargedMCEta16 - nChargedMCEta10;
929  nChargedMCPrimary = nChargedMCPrimaryEta16 - nChargedMCPrimaryEta10;
930  nChargedMCPhysicalPrimary = nChargedMCPhysicalPrimaryEta16 - nChargedMCPhysicalPrimaryEta10;
931  } else if(fMCPrimariesEstimator==kEta05){
932  nChargedMC = nChargedMCEta05;
933  nChargedMCPrimary = nChargedMCPrimaryEta05;
934  nChargedMCPhysicalPrimary = nChargedMCPhysicalPrimaryEta05;
935  } else if(fMCPrimariesEstimator==kEta03){
936  nChargedMC = nChargedMCEta03;
937  nChargedMCPrimary = nChargedMCPrimaryEta03;
938  nChargedMCPhysicalPrimary = nChargedMCPhysicalPrimaryEta03;
939  } else if(fMCPrimariesEstimator==kEtaVZERO){
940  nChargedMC = nChargedMCEtam37tm17 + nChargedMCEta28t51;
941  nChargedMCPrimary = nChargedMCPrimaryEtam37tm17 + nChargedMCPrimaryEta28t51;
942  nChargedMCPhysicalPrimary = nChargedMCPhysicalPrimaryEtam37tm17 + nChargedMCPhysicalPrimaryEta28t51;
943  } else if(fMCPrimariesEstimator==kEtaVZEROA){
944  nChargedMC = nChargedMCEta28t51;
945  nChargedMCPrimary = nChargedMCPrimaryEta28t51;
946  nChargedMCPhysicalPrimary = nChargedMCPhysicalPrimaryEta28t51;
947  }
948  // Here fill the MC correlation plots
949  if(nChargedMCPhysicalPrimary>0){ // INEL>0 for |eta|<1
950  fHistGenPrimaryParticlesInelGt0->Fill(nChargedMCPhysicalPrimary,nchWeight);
951  }
952 
953  fHistNtrVsNchMC->Fill(nChargedMC,countMult,nchWeight);
954  fHistNtrCorrVsNchMC->Fill(nChargedMC,countCorr,nchWeight);
955 
956  fHistNtrVsNchMCPrimary->Fill(nChargedMCPrimary,countMult,nchWeight);
957  fHistNtrCorrVsNchMCPrimary->Fill(nChargedMCPrimary,countCorr,nchWeight);
958 
959  fHistNtrVsNchMCPhysicalPrimary->Fill(nChargedMCPhysicalPrimary,countMult,nchWeight);
960  fHistNtrCorrVsNchMCPhysicalPrimary->Fill(nChargedMCPhysicalPrimary,countCorr,nchWeight);
961 
962  fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary->Fill(nChargedMC,nChargedMCPrimary,nChargedMCPhysicalPrimary,nchWeight);
963  }
964 
965  Int_t nCand = arrayCand->GetEntriesFast();
966  Int_t nSelectedNoPID=0,nSelectedPID=0,nSelectedInMassPeak=0;
967  Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();
968  Double_t mDplusPDG = TDatabasePDG::Instance()->GetParticle(411)->Mass();
969  Double_t mDstarPDG = TDatabasePDG::Instance()->GetParticle(413)->Mass();
970 
971  // pdg of daughters needed for D* too
972  UInt_t pdgDgDStartoD0pi[2]={421,211};
973 
974  Double_t aveMult=0.;
975  Double_t nSelCand=0.;
976  for (Int_t iCand = 0; iCand < nCand; iCand++) {
977  AliAODRecoDecayHF *d = (AliAODRecoDecayHF*)arrayCand->UncheckedAt(iCand);
978  AliAODRecoCascadeHF *dCascade = NULL;
979  if(fPdgMeson==413) dCascade = (AliAODRecoCascadeHF*)d;
980 
981  fHistNEvents->Fill(7);
982  if(fUseBit && !d->HasSelectionBit(selbit)){
983  fHistNEvents->Fill(8);
984  continue;
985  }
986 
987  Double_t ptCand = d->Pt();
988  Double_t rapid=d->Y(fPdgMeson);
989  Bool_t isFidAcc=fRDCutsAnalysis->IsInFiducialAcceptance(ptCand,rapid);
990  if(!isFidAcc) continue;
991 
992  Int_t passAllCuts=fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kAll,aod);
993  Int_t passTopolCuts=fRDCutsAnalysis->GetIsSelectedCuts();
994  if(passTopolCuts==0) continue;
995  nSelectedNoPID++;
996  fHistNEvents->Fill(9);
997  if(passAllCuts){
998  nSelectedPID++;
999  fHistNEvents->Fill(10);
1000  }
1001  Double_t multForCand = countCorr;
1002 
1003  //Weight according to ptcand, using FONLL
1004  Double_t ptWeight = 1.;
1005  if(fUsePtWeight) ptWeight = GetPtWeight(ptCand);
1006  fWeight = nchWeight*ptWeight;
1007 
1008  if(fSubtractTrackletsFromDau){
1009  // For the D* case, subtract only the D0 daughter tracks <=== FIXME !!
1010  AliAODRecoDecayHF2Prong* d0fromDstar = NULL;
1011  if(fPdgMeson==413) d0fromDstar = (AliAODRecoDecayHF2Prong*)dCascade->Get2Prong();
1012 
1013  for(Int_t iDau=0; iDau<nDau; iDau++){
1014  AliAODTrack *t = NULL;
1015  if(fPdgMeson==413){ t = (AliAODTrack*)d0fromDstar->GetDaughter(iDau); }
1016  else{ t = (AliAODTrack*)d->GetDaughter(iDau); }
1017  if(!t) continue;
1018  if(t->HasPointOnITSLayer(0) && t->HasPointOnITSLayer(1)){
1019  if(multForCand>0) multForCand-=1;
1020  }
1021  }
1022  }
1023 
1024  Bool_t isPrimary=kTRUE;
1025  Double_t trueImpParXY=9999.;
1026  Double_t impparXY=d->ImpParXY()*10000.;
1027  Double_t dlen=0.1; //FIXME
1028  Double_t mass[2];
1029  if(fPdgMeson==411){
1030  mass[0]=d->InvMass(nDau,pdgDau);
1031  mass[1]=-1.;
1032  if(TMath::Abs(mass[0]-mDplusPDG)<0.02) nSelectedInMassPeak++; //20 MeV for now... FIXME
1033  }else if(fPdgMeson==421){
1034  UInt_t pdgdaughtersD0[2]={211,321};//pi,K
1035  UInt_t pdgdaughtersD0bar[2]={321,211};//K,pi
1036  mass[0]=d->InvMass(2,pdgdaughtersD0);
1037  mass[1]=d->InvMass(2,pdgdaughtersD0bar);
1038  if(TMath::Abs(mass[0]-mD0PDG)<0.02 || TMath::Abs(mass[1]-mD0PDG)<0.02 ) nSelectedInMassPeak++; //20 MeV for now... FIXME
1039  }else if(fPdgMeson==413){
1040  // FIXME
1041  mass[0]=dCascade->DeltaInvMass();
1042  mass[1]=-1.;
1043  if(TMath::Abs(mass[0]-(mDstarPDG-mD0PDG))<0.0015) nSelectedInMassPeak++; //1 MeV for now... FIXME
1044  }
1045 
1046  Int_t labD0=-1;
1047  Double_t recSpherocity = -0.5;
1048 
1049  // subtract D-meson daughters from spherocity calculation !!
1050  if(fRecomputeSpherocity){
1051  Int_t totTrkToSkip = d->GetNDaughters();
1052  const Int_t nTrkToSkip = totTrkToSkip;
1053  Int_t idToSkip[nTrkToSkip];
1054  for(Int_t i=0; i<nTrkToSkip; i++) idToSkip[i]=-1;
1055 
1056  for(Int_t iDau=0; iDau<nTrkToSkip; iDau++){
1057  AliAODTrack *t = NULL;
1058  t = dynamic_cast<AliAODTrack*>(d->GetDaughter(iDau));
1059  if(!t) continue;
1060  idToSkip[iDau] = t->GetID();
1061  }
1062  recSpherocity=AliVertexingHFUtils::GetSpherocity(aod, fetaMin, fetaMax, fptMin, fptMax, ffiltbit1, ffiltbit2, fminMult, fphiStepSizeDeg, nTrkToSkip, idToSkip);
1063  }
1064 
1065  // remove D0 from Dstar at reconstruction !!
1066  if(fReadMC && fRemoveD0fromDstar){
1067  if(fPdgMeson==421){
1068  labD0 = d->MatchToMC(fPdgMeson,arrayMC,nDau,(Int_t*)pdgDau);
1069  if(labD0>=0){
1070  Bool_t keep=kTRUE;
1071  AliAODMCParticle* mcMoth = dynamic_cast<AliAODMCParticle*>(arrayMC->At(labD0));
1072  Int_t motherD0 = mcMoth->GetMother();
1073  AliAODMCParticle* mcMothD0 = dynamic_cast<AliAODMCParticle*>(arrayMC->At(motherD0));
1074  if(!mcMothD0) continue;
1075  if(TMath::Abs(mcMothD0->GetPdgCode())==413) keep=kFALSE;
1076  if(!keep) continue;
1077  }
1078  }
1079  }
1080 
1081  Int_t labD=-1;
1082  Int_t Origin = 0;
1083 
1084  for(Int_t iHyp=0; iHyp<2; iHyp++){
1085  if(mass[iHyp]<0.) continue; // for D+ and D* we have 1 mass hypothesis
1086  Double_t invMass=mass[iHyp];
1087  Double_t arrayForSparse[5]={invMass,ptCand,impparXY,dlen,multForCand};
1088 
1089  if(fReadMC){
1090  if(fPdgMeson==413){
1091  labD = dCascade->MatchToMC(fPdgMeson,421,(Int_t*)pdgDgDStartoD0pi,(Int_t*)pdgDau,arrayMC);
1092  } else {
1093  labD = d->MatchToMC(fPdgMeson,arrayMC,nDau,(Int_t*)pdgDau);
1094  }
1095 
1096  Bool_t fillHisto=fDoImpPar;
1097  if(labD>=0){
1098  AliAODMCParticle *partD = (AliAODMCParticle*)arrayMC->At(labD);
1099  Int_t code=partD->GetPdgCode();
1100  Origin = AliVertexingHFUtils::CheckOrigin(arrayMC,partD,fUseQuarkTag);
1101  if(Origin==5) isPrimary=kFALSE;
1102  if(code<0 && iHyp==0) fillHisto=kFALSE;
1103  if(code>0 && iHyp==1) fillHisto=kFALSE;
1104  if(!isPrimary){
1105  if(fPdgMeson==411){
1106  trueImpParXY=AliVertexingHFUtils::GetTrueImpactParameterDplus(mcHeader,arrayMC,partD)*10000.;
1107  }else if(fPdgMeson==421){
1108  trueImpParXY=AliVertexingHFUtils::GetTrueImpactParameterDzero(mcHeader,arrayMC,partD)*10000.;
1109  }else if(fPdgMeson==413){
1110  trueImpParXY=0.;
1111  }
1112  Double_t arrayForSparseTrue[5]={invMass,ptCand,trueImpParXY,dlen,multForCand};
1113  if(fillHisto && passAllCuts){
1114  fHistMassPtImpPar[2]->Fill(arrayForSparse, fWeight);
1115  fHistMassPtImpPar[3]->Fill(arrayForSparseTrue, fWeight);
1116  }
1117  }else{
1118  if(fillHisto && passAllCuts) fHistMassPtImpPar[1]->Fill(arrayForSparse, fWeight);
1119  }
1120  }else{
1121  if(fillHisto && passAllCuts)fHistMassPtImpPar[4]->Fill(arrayForSparse, fWeight);
1122  }
1123  if(TMath::Abs(labD)==fPdgMeson && fMCOption==2) continue;
1124  if(TMath::Abs(labD)!=fPdgMeson && fMCOption==1) continue;
1125  }
1126  if(fPdgMeson==421){
1127  if(iHyp==0 && !(passTopolCuts&1)) continue; // candidate not passing as D0
1128  if(iHyp==1 && !(passTopolCuts&2)) continue; // candidate not passing as D0bar
1129  }
1130  if(fFillSoSparseChecks == 2 || fFillSoSparseChecks == 3){ //Filling THnSparse for Spherocity without PID
1131  if(fCalculateSphericity){
1132  Double_t arrayForSparseSoNoPid[5]={ptCand, invMass, spherocity, multForCand, sphericity};
1133  fSparseEvtShapewithNoPid->Fill(arrayForSparseSoNoPid, fWeight);
1134  }else{
1135  Double_t arrayForSparseSoNoPid[4]={ptCand, invMass, spherocity, multForCand};
1136  fSparseEvtShapewithNoPid->Fill(arrayForSparseSoNoPid, fWeight);
1137  }
1138  }
1139  if(fPdgMeson==421){
1140  if(iHyp==0 && !(passAllCuts&1)) continue; // candidate not passing as D0
1141  if(iHyp==1 && !(passAllCuts&2)) continue; // candidate not passing as D0bar
1142  }
1143  if(passAllCuts){
1144  aveMult+=multForCand;
1145  nSelCand+=1.;
1146 
1147  if(fFillSoSparseChecks == 1 || fFillSoSparseChecks == 3){
1148  if(fCalculateSphericity){
1149  Double_t arrayForSparseSowithMultUnncorr[6]={ptCand, invMass, spherocity, multForCand, (Double_t)countTreta1, sphericity};
1150  fSparseEvtShape->Fill(arrayForSparseSowithMultUnncorr, fWeight);
1151  }else{
1152  Double_t arrayForSparseSowithMultUnncorr[5]={ptCand, invMass, spherocity, multForCand, (Double_t)countTreta1};
1153  fSparseEvtShape->Fill(arrayForSparseSowithMultUnncorr, fWeight);
1154  }
1155  }
1156  else{
1157  if(fCalculateSphericity){
1158  Double_t arrayForSparseSo[5]={ptCand, invMass, spherocity, multForCand, sphericity};
1159  fSparseEvtShape->Fill(arrayForSparseSo, fWeight);
1160  }else{
1161  Double_t arrayForSparseSo[4]={ptCand, invMass, spherocity, multForCand};
1162  fSparseEvtShape->Fill(arrayForSparseSo, fWeight);
1163  }
1164  }
1165 
1166  if(fRecomputeSpherocity){
1167  Double_t arrayForSparseRecSphero[5]={ptCand, invMass, spherocity, multForCand, recSpherocity};
1168  fSparseEvtShapeRecSphero->Fill(arrayForSparseRecSphero, fWeight);
1169  }
1170 
1171  if(fReadMC){
1172  if(fRecomputeSpherocity){
1173  Double_t arrayForSparseSoPromptFD[5]={ptCand, invMass, spherocity, multForCand, recSpherocity};
1174 
1175  if(Origin==4) fSparseEvtShapePrompt->Fill(arrayForSparseSoPromptFD, fWeight);
1176  else if(Origin==5) fSparseEvtShapeFeeddown->Fill(arrayForSparseSoPromptFD, fWeight);
1177  }else{
1178  Double_t arrayForSparseSoPromptFD[4]={ptCand, invMass, spherocity, multForCand};
1179 
1180  if(Origin==4) fSparseEvtShapePrompt->Fill(arrayForSparseSoPromptFD, fWeight);
1181  else if(Origin==5) fSparseEvtShapeFeeddown->Fill(arrayForSparseSoPromptFD, fWeight);
1182  }
1183  }
1184  if(labD>=0){
1185  Bool_t keepCase=kTRUE;
1186  if(fPdgMeson==421){
1187  AliAODMCParticle *partD = (AliAODMCParticle*)arrayMC->At(labD);
1188  Int_t code=partD->GetPdgCode();
1189  if(code<0 && iHyp==0) keepCase=kFALSE;
1190  if(code>0 && iHyp==1) keepCase=kFALSE;
1191  }
1192  if(keepCase) FillMCMassHistos(arrayMC,labD, multForCand, spherocity, sphericity, recSpherocity, nchWeight);
1193  }
1194  if(fDoImpPar) fHistMassPtImpPar[0]->Fill(arrayForSparse, fWeight);
1195  }
1196  }
1197  }
1198  if(fSubtractTrackletsFromDau && nSelCand>0){
1199  aveMult/=nSelCand;
1200  fCounterCandidates->StoreEvent(aod,fRDCutsAnalysis,fReadMC,(Int_t)(aveMult+0.5001),spherocity);
1201  }else{
1202  fCounterCandidates->StoreEvent(aod,fRDCutsAnalysis,fReadMC,(Int_t)countCorr,spherocity);
1203  }
1204 
1205  fCounterCandidates->StoreCandidates(aod,nSelectedNoPID,kTRUE);
1206  fCounterCandidates->StoreCandidates(aod,nSelectedPID,kFALSE);
1207  fHistNtrCorrEvSel->Fill(countCorr,nchWeight);
1208  if(nSelectedPID>0) fHistNtrCorrEvWithCand->Fill(countCorr,nchWeight);
1209  if(nSelectedInMassPeak>0) fHistNtrCorrEvWithD->Fill(countCorr,nchWeight);
1210 
1211  if(fFillTrackHisto) FillTrackControlHisto(aod, countCorr, nSelectedInMassPeak);
1212 
1213  PostData(1,fOutput);
1214  PostData(2,fListCuts);
1215  PostData(3,fOutputCounters);
1216  if(fReadMC) PostData(5,fOutputEffCorr);
1217 
1218  return;
1219 }
1220 
1221 //________________________________________________________________________
1222 void AliAnalysisTaskSEDvsEventShapes::CreateImpactParameterHistos()
1223 {
1224  // Histos for impact paramter study
1225  // mass . pt , impact parameter , decay length , multiplicity
1226 
1227  Int_t nbins[5]={fNMassBins,200,fNImpParBins,50,100};
1228  Double_t xmin[5]={fLowmasslimit,0.,fLowerImpPar,0.,0.};
1229  Double_t xmax[5]={fUpmasslimit,20.,fHigherImpPar,1.,100.};
1230 
1231  fHistMassPtImpPar[0]=new THnSparseF("hMassPtImpParAll",
1232  "Mass vs. pt vs.imppar - All",
1233  5,nbins,xmin,xmax);
1234  fHistMassPtImpPar[1]=new THnSparseF("hMassPtImpParPrompt",
1235  "Mass vs. pt vs.imppar - promptD",
1236  5,nbins,xmin,xmax);
1237  fHistMassPtImpPar[2]=new THnSparseF("hMassPtImpParBfeed",
1238  "Mass vs. pt vs.imppar - DfromB",
1239  5,nbins,xmin,xmax);
1240  fHistMassPtImpPar[3]=new THnSparseF("hMassPtImpParTrueBfeed",
1241  "Mass vs. pt vs.true imppar -DfromB",
1242  5,nbins,xmin,xmax);
1243  fHistMassPtImpPar[4]=new THnSparseF("hMassPtImpParBkg",
1244  "Mass vs. pt vs.imppar - backgr.",
1245  5,nbins,xmin,xmax);
1246  for(Int_t i=0; i<5;i++){
1247  fOutput->Add(fHistMassPtImpPar[i]);
1248  }
1249 }
1250 
1251 //________________________________________________________________________
1252 void AliAnalysisTaskSEDvsEventShapes::Terminate(Option_t */*option*/)
1253 {
1254  // Terminate analysis
1255  //
1256  if(fDebug > 1) printf("AnalysisTaskSEDvsMultiplicity: Terminate() \n");
1257 
1258  fOutput = dynamic_cast<TList*> (GetOutputData(1));
1259  if (!fOutput) {
1260  printf("ERROR: fOutput not available\n");
1261  return;
1262  }
1263 
1264  fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
1265  if(!fHistNEvents){
1266  printf("ERROR: fHistNEvents not available\n");
1267  return;
1268  }
1269  printf("Number of Analyzed Events = %d\n",(Int_t)fHistNEvents->GetBinContent(3));
1270 
1271  return;
1272 }
1273 
1274 //____________________________________________________________________________
1275 TProfile* AliAnalysisTaskSEDvsEventShapes::GetEstimatorHistogram(const AliVEvent* event)
1276 {
1277  // Get Estimator Histogram from period event->GetRunNumber();
1278  //
1279  // If you select SPD tracklets in |eta|<1 you should use type == 1
1280  //
1281 
1282  Int_t runNo = event->GetRunNumber();
1283  Int_t period = -1; // pp: 0-LHC10b, 1-LHC10c, 2-LHC10d, 3-LHC10e
1284  // pPb: 0-LHC13b, 1-LHC13c
1285  if (fisPPbData) {
1286  if (runNo>195343 && runNo<195484) period = 0;
1287  if (runNo>195528 && runNo<195678) period = 1;
1288  if (period < 0 || period > 1) return 0;
1289  }
1290  else {
1291  if(runNo>114930 && runNo<117223) period = 0;
1292  if(runNo>119158 && runNo<120830) period = 1;
1293  if(runNo>122373 && runNo<126438) period = 2;
1294  if(runNo>127711 && runNo<130851) period = 3;
1295  if(period<0 || period>3) return 0;
1296  }
1297 
1298  return fMultEstimatorAvg[period];
1299 }
1300 //__________________________________________________________________________________________________
1301 
1302 Double_t AliAnalysisTaskSEDvsEventShapes::GetPtWeight(Float_t pt)
1303 {
1304  //
1305  // calculating the weight to fill the container
1306  //
1307 
1308  // FNOLL central:
1309  // p0 = 1.63297e-01 --> 0.322643
1310  // p1 = 2.96275e+00
1311  // p2 = 2.30301e+00
1312  // p3 = 2.50000e+00
1313 
1314  // PYTHIA
1315  // p0 = 1.85906e-01 --> 0.36609
1316  // p1 = 1.94635e+00
1317  // p2 = 1.40463e+00
1318  // p3 = 2.50000e+00
1319 
1320  Double_t func1[4] = {0.322643,2.96275,2.30301,2.5};
1321  Double_t func2[4] = {0.36609,1.94635,1.40463,2.5};
1322 
1323  Double_t dndpt_func1 = dNdptFit(pt,func1);
1324  //if(fUseFlatPtWeight) dndpt_func1 = 1./30.;
1325  Double_t dndpt_func2 = dNdptFit(pt,func2);
1326  AliDebug(2,Form("pt = %f, FONLL = %f, Pythia = %f, ratio = %f",pt,dndpt_func1,dndpt_func2,dndpt_func1/dndpt_func2));
1327  return dndpt_func1/dndpt_func2;
1328 }
1329 
1330 //__________________________________________________________________________________________________
1331 Double_t AliAnalysisTaskSEDvsEventShapes::dNdptFit(Float_t pt, Double_t* par)
1332 {
1333  //
1334  // calculating dNdpt
1335  //
1336 
1337  Double_t denom = TMath::Power((pt/par[1]), par[3] );
1338  Double_t dNdpt = par[0]*pt/TMath::Power(1.+denom, par[2]);
1339 
1340  return dNdpt;
1341 }
1342 
1343 
1344 //__________________________________________________________________________________________________
1345 void AliAnalysisTaskSEDvsEventShapes::CreateMeasuredNchHisto()
1346 {
1347  // creates historgam with measured multiplcity distribution in pp 7 TeV collisions (from Eur. Phys. J. C (2010) 68: 345–354)
1348  //
1349  // for Nch > 70 the points were obtainedwith a double NBD distribution
1350  // TF1 *fit1 = new TF1("fit1","[0]*(TMath::Gamma(x+[1])/(TMath::Gamma(x+1)*TMath::Gamma([1])))*(TMath::Power(([2]/[1]),x))*(TMath::Power((1+([2]/[1])),-x-[1]))"); fit1->SetParameter(0,1.);// normalization constant
1351  // fit1->SetParameter(1,1.63); // k parameter
1352  // fit1->SetParameter(2,12.8); // mean multiplicity
1353  Double_t nchbins[82]={0.50,1.50,2.50,3.50,4.50,5.50,6.50,7.50,8.50,9.50,
1354  10.50,11.50,12.50,13.50,14.50,15.50,16.50,17.50,18.50,19.50,
1355  20.50,21.50,22.50,23.50,24.50,25.50,26.50,27.50,28.50,29.50,
1356  30.50,31.50,32.50,33.50,34.50,35.50,36.50,37.50,38.50,39.50,
1357  40.50,41.50,42.50,43.50,44.50,45.50,46.50,47.50,48.50,49.50,
1358  50.50,51.50,52.50,53.50,54.50,55.50,56.50,57.50,58.50,59.50,
1359  60.50,62.50,64.50,66.50,68.50,70.50,72.50,74.50,76.50,78.50,
1360  80.50,82.50,84.50,86.50,88.50,90.50,92.50,94.50,96.50,98.50,
1361  100.50,102.50};
1362  Double_t pch[81]={0.062011,0.072943,0.070771,0.067245,0.062834,0.057383,0.051499,0.04591,0.041109,0.036954,
1363  0.03359,0.030729,0.028539,0.026575,0.024653,0.0229,0.021325,0.019768,0.018561,0.017187,
1364  0.01604,0.014836,0.013726,0.012576,0.011481,0.010393,0.009502,0.008776,0.008024,0.007452,
1365  0.006851,0.006428,0.00594,0.005515,0.005102,0.00469,0.004162,0.003811,0.003389,0.003071,
1366  0.002708,0.002422,0.002184,0.001968,0.00186,0.00165,0.001577,0.001387,0.001254,0.001118,
1367  0.001037,0.000942,0.000823,0.000736,0.000654,0.000579,0.000512,0.00049,0.00045,0.000355,
1368  0.000296,0.000265,0.000193,0.00016,0.000126,0.0000851, 0.0000676,0.0000537,0.0000426, 0.0000338,
1369  0.0000268,0.0000213,0.0000166,0.0000133,0.0000106,0.00000837,0.00000662, 0.00000524,0.00000414, 0.00000327,
1370  0.00000258};
1371 
1372  if(fHistoMeasNch) delete fHistoMeasNch;
1373  fHistoMeasNch=new TH1F("hMeaseNch","",81,nchbins);
1374  for(Int_t i=0; i<81; i++){
1375  fHistoMeasNch->SetBinContent(i+1,pch[i]);
1376  fHistoMeasNch->SetBinError(i+1,0.);
1377  }
1378 }
1379 
1380 //________________________________________________________________________
1381 Bool_t AliAnalysisTaskSEDvsEventShapes::FillTrackControlHisto(AliAODEvent* aod, Int_t nSelTrkCorr, Int_t nSelectedEvwithCand)
1382 {
1384 
1385  Int_t nTracks=aod->GetNumberOfTracks();
1386  Int_t nSelTracks=0;
1387 
1388  Double_t* etaArr=new Double_t[nTracks];
1389  Double_t* phiArr=new Double_t[nTracks];
1390  Double_t sumpt=0.;
1391 
1392  for(Int_t it=0; it<nTracks; it++) {
1393  AliAODTrack *tr=dynamic_cast<AliAODTrack*>(aod->GetTrack(it));
1394  if(!tr) continue;
1395  Float_t eta = tr->Eta();
1396  Float_t pt = tr->Pt();
1397  Float_t phi = tr->Phi();
1398  if(eta<fetaMin || eta>fetaMax) continue;
1399  if(pt<fptMin || pt>fptMax) continue;
1400  Bool_t fb1 = tr->TestFilterBit(ffiltbit1);
1401  Bool_t fb2 = tr->TestFilterBit(ffiltbit2);
1402  if( !(fb1 || fb2) ) continue;
1403  if(ffiltbit1==0 || ffiltbit2==0){
1404  Int_t status=tr->GetStatus();
1405  if(!(status & AliAODTrack::kTPCrefit)) continue;
1406  }
1407  etaArr[nSelTracks]=eta;
1408  phiArr[nSelTracks]=phi;
1409  nSelTracks++;
1410  sumpt+=pt;
1411  }
1412 
1413  if(nSelTracks<fminMult) return kFALSE;
1414 
1415  if(nSelectedEvwithCand>0) fHistNtrVsnTrackEvWithCand->Fill(nSelTracks, nSelTrkCorr);
1416 
1417  for(Int_t iselt=0; iselt<nSelTracks; iselt++) {
1418  fHistnTrackvsEtavsPhi->Fill(etaArr[iselt], phiArr[iselt], nSelTracks);
1419  if(nSelectedEvwithCand>0) fHistnTrackvsEtavsPhiEvWithCand->Fill(etaArr[iselt], phiArr[iselt], nSelTracks);
1420  }
1421 
1422  return kTRUE;
1423 
1424 }
1425 
1426 //__________________________________________________________________________________________________
1427 void AliAnalysisTaskSEDvsEventShapes::FillMCMassHistos(TClonesArray *arrayMC, Int_t labD, Double_t countMult, Double_t spherocity, Double_t sphericity, Double_t recSpherocity, Double_t nchWeight)
1428 {
1429  // Function to fill the true MC signal
1430 
1431  AliAODMCParticle *partD = (AliAODMCParticle*)arrayMC->At(labD);
1432 
1433  Double_t mass = partD->M();
1434  Double_t pt = partD->Pt();
1435  Double_t rapid = partD->Y();
1436 
1437  //Weight according to pt, using FONLL
1438  Double_t ptWeight = 1.;
1439  if(fUsePtWeight) ptWeight = GetPtWeight(pt);
1440  fWeight = nchWeight*ptWeight;
1441 
1442  Int_t orig=AliVertexingHFUtils::CheckOrigin(arrayMC,partD,fUseQuarkTag); // Prompt = 4, FeedDown = 5
1443 
1444  //for prompt
1445  if(orig == 4){
1446  //fill histo for prompt
1447  Double_t arrayMCRecoRecSpheroPrompt[5] = {pt, countMult, spherocity, rapid, recSpherocity};
1448  Double_t arrayMCRecoPrompt[4] = {pt, countMult, spherocity, rapid};
1449  if(fRecomputeSpherocity) fMCRecoPrompt->Fill(arrayMCRecoRecSpheroPrompt, fWeight);
1450  else fMCRecoPrompt->Fill(arrayMCRecoPrompt, fWeight);
1451  if(fCalculateSphericity) {
1452  Double_t arrayMCRecoPromptSpheri[4] = {pt, countMult, sphericity, rapid};
1453  fMCRecoPromptSpheri->Fill(arrayMCRecoPromptSpheri, fWeight);
1454  }
1455  }
1456  //for FD
1457  else if(orig == 5){
1458  //fill histo for FD
1459  Double_t arrayMCRecoRecSpheroFeeddown[5] = {pt, countMult, spherocity, rapid, recSpherocity};
1460  Double_t arrayMCRecoFeeddown[4] = {pt, countMult, spherocity, rapid};
1461  if(fRecomputeSpherocity) fMCRecoFeeddown->Fill(arrayMCRecoRecSpheroFeeddown, fWeight);
1462  else fMCRecoFeeddown->Fill(arrayMCRecoFeeddown, fWeight);
1463  if(fCalculateSphericity){
1464  Double_t arrayMCRecoFeeddownSpheri[4] = {pt, countMult, sphericity, rapid};
1465  fMCRecoFeeddownSpheri->Fill(arrayMCRecoFeeddownSpheri, fWeight);
1466  }
1467  }
1468 
1469  Double_t arrayMCReco[4] = {pt, countMult, spherocity, rapid};
1470  fMCRecoBothPromptFD->Fill(arrayMCReco, fWeight);
1471  if(fCalculateSphericity){
1472  Double_t arrayMCRecoSpheri[4] = {pt, countMult, sphericity, rapid};
1473  fMCRecoBothPromptFDSpheri->Fill(arrayMCRecoSpheri, fWeight);
1474  }
1475 
1476 }
1477 
1478 //__________________________________________________________________________________________________
1479 void AliAnalysisTaskSEDvsEventShapes::FillMCGenAccHistos(AliAODEvent* aod, TClonesArray *arrayMC, AliAODMCHeader *mcHeader, Double_t countMult, Double_t spherocity, Double_t sphericity, Bool_t isEvSel, Double_t nchWeight)
1480 {
1481 
1483 
1484  Int_t nProng=2;
1485  Int_t totPart = arrayMC->GetEntriesFast(); //number of particles
1486  Int_t totTracks = aod->GetNumberOfTracks(); // number of tracks
1487  Double_t recSpherocity = -0.5;
1488 
1489  const Int_t nPart = totPart;
1490  Int_t trkToSkip[nPart];
1491  for(Int_t i=0; i<nPart; i++) trkToSkip[i]=-1; //stores ID of tracks at i'th label of particle
1492 
1493  if(fRecomputeSpherocity){
1494  for(Int_t it=0; it<totTracks; it++){
1495  AliAODTrack *t=dynamic_cast<AliAODTrack*>(aod->GetTrack(it));
1496  if(!t) continue;
1497  if(!(t->TestFilterMask(BIT(4)))) continue;
1498  //if(!(t->HasPointOnITSLayer(0)) && !(t->HasPointOnITSLayer(1))) continue;
1499 
1500  Int_t lab=TMath::Abs(t->GetLabel());
1501  Int_t id=t->GetID();
1502  Double_t pt = t->Pt();
1503  if(id<0) continue;
1504  if(pt<0.3) continue;
1505  AliAODMCParticle* genDup = dynamic_cast<AliAODMCParticle*>(arrayMC->At(lab));
1506  Double_t xver = genDup->Xv();
1507  Double_t yver = genDup->Yv();
1508  Double_t rver = TMath::Sqrt(xver*xver + yver*yver);
1509  if(rver>3) continue;
1510  trkToSkip[lab] = id;
1511  }
1512  }
1513 
1514  if(fPdgMeson==421){
1515  nProng=2;
1516  }else if(fPdgMeson==411 || fPdgMeson==431){
1517  nProng=3;
1518  }
1519 
1520  Int_t nTrkToSkip = nProng;
1521  Int_t idToSkip[nTrkToSkip];
1522  for(Int_t i=0; i<nTrkToSkip; i++) idToSkip[i]=-1;
1523 
1524  Double_t zMCVertex = mcHeader->GetVtxZ(); //vertex MC
1525 
1526  for(Int_t iPart=0; iPart<totPart; iPart++){
1527  AliAODMCParticle* mcGenPart = dynamic_cast<AliAODMCParticle*>(arrayMC->At(iPart));
1528 
1529  if (TMath::Abs(mcGenPart->GetPdgCode()) == fPdgMeson){
1530  Int_t orig=AliVertexingHFUtils::CheckOrigin(arrayMC,mcGenPart,fUseQuarkTag);//Prompt = 4, FeedDown = 5
1531 
1532  Int_t deca = 0;
1533  Bool_t isGoodDecay=kFALSE;
1534  Int_t labDau[4]={-1,-1,-1,-1};
1535  Bool_t isInAcc = kFALSE;
1536  Bool_t isFidAcc = kFALSE;
1537 
1538  if(fPdgMeson==421){
1539  deca=AliVertexingHFUtils::CheckD0Decay(arrayMC,mcGenPart,labDau);
1540  if(mcGenPart->GetNDaughters()!=2) continue;
1541  if(deca==1) isGoodDecay=kTRUE;
1542  }else if(fPdgMeson==411){
1543  deca=AliVertexingHFUtils::CheckDplusDecay(arrayMC,mcGenPart,labDau);
1544  if(deca>0) isGoodDecay=kTRUE;
1545  }else if(fPdgMeson==431){
1546  deca=AliVertexingHFUtils::CheckDsDecay(arrayMC,mcGenPart,labDau);
1547  if(deca==1) isGoodDecay=kTRUE;
1548  }else if(fPdgMeson==413){
1549  deca=AliVertexingHFUtils::CheckDstarDecay(arrayMC,mcGenPart,labDau);
1550  if(deca==1) isGoodDecay=kTRUE;
1551  }
1552 
1553  if(labDau[0]==-1){
1554  continue; //protection against unfilled array of labels
1555  }
1556 
1557  if(fRecomputeSpherocity && isGoodDecay){
1558  for(Int_t iDau=0; iDau<nTrkToSkip; iDau++){
1559  Int_t indexDau = TMath::Abs(mcGenPart->GetDaughter(iDau)); //index of daughter i.e. label
1560  idToSkip[iDau] = trkToSkip[indexDau];
1561  }
1562  recSpherocity=AliVertexingHFUtils::GetSpherocity(aod, fetaMin, fetaMax, fptMin, fptMax, ffiltbit1, ffiltbit2, fminMult, fphiStepSizeDeg, nTrkToSkip, idToSkip);
1563  }
1564  if(fPdgMeson==421){
1565  //Removal of D0 from Dstar at Generation !!
1566  if(fRemoveD0fromDstar){
1567  Int_t mother = mcGenPart->GetMother();
1568  AliAODMCParticle* mcMoth = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
1569  if(!mcMoth) continue;
1570  if(TMath::Abs(mcMoth->GetPdgCode())==413) continue;
1571  }
1572  }
1573 
1574  Double_t pt = mcGenPart->Pt();
1575  Double_t rapid = mcGenPart->Y();
1576 
1577  //Weight according to pt, using FONLL
1578  Double_t ptWeight = 1.;
1579  if(fUsePtWeight) ptWeight = GetPtWeight(pt);
1580  fWeight = nchWeight*ptWeight;
1581 
1582  isFidAcc=fRDCutsAnalysis->IsInFiducialAcceptance(pt,rapid);
1583  isInAcc=CheckGenAcc(arrayMC,nProng,labDau);
1584 
1585  if(isGoodDecay && TMath::Abs(zMCVertex) < fRDCutsAnalysis->GetMaxVtxZ() && isFidAcc && isInAcc) {
1586  //for prompt
1587  if(orig == 4){
1588  //fill histo for prompt
1589  Double_t arrayMCGenRecSpheroPrompt[5] = {pt, countMult, spherocity, rapid, recSpherocity};
1590  Double_t arrayMCGenPrompt[4] = {pt, countMult, spherocity, rapid};
1591  if(fRecomputeSpherocity) fMCAccGenPrompt->Fill(arrayMCGenRecSpheroPrompt, fWeight);
1592  else fMCAccGenPrompt->Fill(arrayMCGenPrompt, fWeight);
1593  if(isEvSel) fMCAccGenPromptEvSel->Fill(arrayMCGenPrompt, fWeight);
1594  if(fCalculateSphericity){
1595  Double_t arrayMCGenPromptSpheri[4] = {pt, countMult, sphericity, rapid};
1596  fMCAccGenPromptSpheri->Fill(arrayMCGenPromptSpheri, fWeight);
1597  }
1598  }
1599  //for FD
1600  else if(orig == 5){
1601  //fill histo for FD
1602  Double_t arrayMCGenRecSpheroFeeddown[5] = {pt, countMult, spherocity, rapid, recSpherocity};
1603  Double_t arrayMCGenFeeddown[4] = {pt, countMult, spherocity, rapid};
1604  if(fRecomputeSpherocity) fMCAccGenFeeddown->Fill(arrayMCGenRecSpheroFeeddown, fWeight);
1605  else fMCAccGenFeeddown->Fill(arrayMCGenFeeddown, fWeight);
1606  if(isEvSel) fMCAccGenFeeddownEvSel->Fill(arrayMCGenFeeddown, fWeight);
1607  if(fCalculateSphericity){
1608  Double_t arrayMCGenFeeddownSpheri[4] = {pt, countMult, sphericity, rapid};
1609  fMCAccGenFeeddownSpheri->Fill(arrayMCGenFeeddownSpheri, fWeight);
1610  }
1611  }
1612  else
1613  continue;
1614  }
1615  }
1616  }
1617 }
1618 
1619 //_________________________________________________________________
1620 Bool_t AliAnalysisTaskSEDvsEventShapes::CheckGenAcc(TClonesArray* arrayMC, Int_t nProng, Int_t *labDau)
1621 {
1623  for (Int_t iProng = 0; iProng<nProng; iProng++){
1624  AliAODMCParticle* mcPartDaughter=dynamic_cast<AliAODMCParticle*>(arrayMC->At(labDau[iProng]));
1625  if(!mcPartDaughter) return kFALSE;
1626  Double_t eta = mcPartDaughter->Eta();
1627  Double_t pt = mcPartDaughter->Pt();
1628  if (TMath::Abs(eta) > fEtaAccCut || pt < fPtAccCut) return kFALSE;
1629  }
1630  return kTRUE;
1631 }
AliAnalysisTaskSEDvsEventShapes::fetaMin
Double_t fetaMin
pt limits for acceptance step
Definition: AliAnalysisTaskSEDvsEventShapes.h:254
charge
Int_t charge
Definition: MuonTrackingEfficiency.C:54
AliAnalysisTaskSEDvsEventShapes::CheckGenAcc
Bool_t CheckGenAcc(TClonesArray *arrayMC, Int_t nProng, Int_t *labDau)
Definition: AliAnalysisTaskSEDvsEventShapes.cxx:1620
pdg
Int_t pdg
Definition: charmCutsOptimization.C:34
AliAnalysisTaskSEDvsEventShapes::fHistNtrCorrVsSo
TH2F * fHistNtrCorrVsSo
hist of ntracklets vs So
Definition: AliAnalysisTaskSEDvsEventShapes.h:168
AliVertexingHFUtils::CheckDplusDecay
static Int_t CheckDplusDecay(AliStack *stack, Int_t label, Int_t *arrayDauLab)
Definition: AliVertexingHFUtils.cxx:914
AliRDHFCutsD0toKpi.h
AliAnalysisTaskSEDvsEventShapes::fMCRecoPromptSpheri
THnSparseD * fMCRecoPromptSpheri
histo for StepMCGenAcc for D meson feeddown for Sphericity
Definition: AliAnalysisTaskSEDvsEventShapes.h:201
Double_t
double Double_t
Definition: External.C:58
TH3F
Definition: External.C:260
AliNormalizationCounter::StoreCandidates
void StoreCandidates(AliVEvent *, Int_t nCand=0, Bool_t flagFilter=kTRUE)
Definition: AliNormalizationCounter.cxx:261
TH2F
Definition: External.C:236
AliRDHFCutsDplustoKpipi.h
AliAODRecoCascadeHF::DeltaInvMass
Double_t DeltaInvMass() const
Definition: AliAODRecoCascadeHF.h:86
AliAnalysisTaskSEDvsEventShapes::fHistNtrCorrVsSpheri
TH2F * fHistNtrCorrVsSpheri
hist of ntracklets vs Spheri
Definition: AliAnalysisTaskSEDvsEventShapes.h:170
AliAnalysisTaskSEDvsEventShapes::fCounterCandidates
AliNormalizationCounter * fCounterCandidates
Counter for normalization, uncorrected multiplicity.
Definition: AliAnalysisTaskSEDvsEventShapes.h:217
AliRDHFCuts::GetIsSelectedCuts
Int_t GetIsSelectedCuts() const
Definition: AliRDHFCuts.h:342
AliAODRecoCascadeHF.h
AliVertexingHFUtils::CheckDstarDecay
static Int_t CheckDstarDecay(AliStack *stack, Int_t label, Int_t *arrayDauLab)
Definition: AliVertexingHFUtils.cxx:1515
AliAODRecoDecayHF::HasSelectionBit
Bool_t HasSelectionBit(Int_t i) const
Definition: AliAODRecoDecayHF.h:124
AliAODRecoCascadeHF::MatchToMC
Int_t MatchToMC(Int_t pdgabs, Int_t pdgabs2prong, Int_t *pdgDg, Int_t *pdgDg2prong, TClonesArray *mcArray, Bool_t isV0=kFALSE) const
Definition: AliAODRecoCascadeHF.cxx:118
AliAnalysisVertexingHF.h
mass
Double_t mass
Definition: charmCutsOptimization.C:35
AliVertexingHFUtils.h
AliAnalysisTaskSEDvsEventShapes::fMCAccGenFeeddown
THnSparseD * fMCAccGenFeeddown
histo for StepMCGenAcc for D meson prompt
Definition: AliAnalysisTaskSEDvsEventShapes.h:195
AliAODRecoDecayHF::ImpParXY
Double_t ImpParXY() const
Definition: AliAODRecoDecayHF.h:84
AliAnalysisTaskSEDvsEventShapes::fHistnTrackvsEtavsPhi
TH3F * fHistnTrackvsEtavsPhi
hist. of ntracklets for evnts with a candidate in D mass peak
Definition: AliAnalysisTaskSEDvsEventShapes.h:186
AliAnalysisTaskSEDvsEventShapes::fMCRecoFeeddown
THnSparseD * fMCRecoFeeddown
histo for StepMCReco for D meson feeddown
Definition: AliAnalysisTaskSEDvsEventShapes.h:197
AliAnalysisTaskSEDvsEventShapes::fSparseEvtShapeFeeddown
THnSparseD * fSparseEvtShapeFeeddown
THnSparse histograms for Prompt D0 vs. Spherocity.
Definition: AliAnalysisTaskSEDvsEventShapes.h:192
AliAnalysisTaskSEDvsEventShapes::fHistNtrVsNchMCPhysicalPrimary
TH2F * fHistNtrVsNchMCPhysicalPrimary
! hist of ntracklets vs Nch (Physical Primary)
Definition: AliAnalysisTaskSEDvsEventShapes.h:176
AliAnalysisTaskSEDvsEventShapes::fHistNtrCorrVsNchMCPhysicalPrimary
TH2F * fHistNtrCorrVsNchMCPhysicalPrimary
! hist of ntracklets vs Nch (Physical Primary)
Definition: AliAnalysisTaskSEDvsEventShapes.h:177
AliAnalysisTaskSEDvsEventShapes.h
AliVertexingHFUtils::CheckDsDecay
static Int_t CheckDsDecay(AliStack *stack, Int_t label, Int_t *arrayDauLab)
Definition: AliVertexingHFUtils.cxx:1252
AliRDHFCuts::GetWhyRejection
Int_t GetWhyRejection() const
Definition: AliRDHFCuts.h:299
AliAnalysisTaskSEDvsEventShapes::fMCAccGenFeeddownEvSel
THnSparseD * fMCAccGenFeeddownEvSel
histo for StepMCGenAcc for D meson prompt with Vertex selection (IsEvSel = kTRUE) ...
Definition: AliAnalysisTaskSEDvsEventShapes.h:206
AliAnalysisTaskSEDvsEventShapes::fHistNtrVsNchMC
TH2F * fHistNtrVsNchMC
hist of ntracklets vs Spheri
Definition: AliAnalysisTaskSEDvsEventShapes.h:172
AliVertexingHFUtils::CheckOrigin
static Int_t CheckOrigin(TClonesArray *arrayMC, AliAODMCParticle *mcPart, Bool_t searchUpToQuark=kTRUE)
Definition: AliVertexingHFUtils.cxx:637
AliRDHFCuts::GetMaxVtxZ
Double_t GetMaxVtxZ() const
Definition: AliRDHFCuts.h:242
AliAnalysisTaskSEDvsEventShapes::fHistNtrCorrVsZvtx
TH2F * fHistNtrCorrVsZvtx
hist of ntracklets vs Zvertex
Definition: AliAnalysisTaskSEDvsEventShapes.h:165
AliAnalysisTaskSEDvsEventShapes::fHistMassPtImpPar
THnSparseF * fHistMassPtImpPar[5]
histo for StepMCGenAcc for D meson feeddown with Vertex selection (IsEvSel = kTRUE) ...
Definition: AliAnalysisTaskSEDvsEventShapes.h:208
AliAnalysisTaskSEDvsEventShapes::fHistNtrVsZvtx
TH2F * fHistNtrVsZvtx
hist. for No. of events
Definition: AliAnalysisTaskSEDvsEventShapes.h:164
AliAnalysisTaskSEDvsEventShapes::fDoImpPar
Bool_t fDoImpPar
Counter for normalization, corrected multiplicity for candidates.
Definition: AliAnalysisTaskSEDvsEventShapes.h:219
AliAnalysisTaskSEDvsEventShapes::fListProfiles
TList * fListProfiles
list send on output slot 3
Definition: AliAnalysisTaskSEDvsEventShapes.h:159
AliAnalysisTaskSEDvsEventShapes::fMCRecoPrompt
THnSparseD * fMCRecoPrompt
histo for StepMCGenAcc for D meson feeddown
Definition: AliAnalysisTaskSEDvsEventShapes.h:196
AliVertexingHFUtils::GetVZEROCEqualizedMultiplicity
static Double_t GetVZEROCEqualizedMultiplicity(AliAODEvent *ev)
Definition: AliVertexingHFUtils.cxx:263
AliAnalysisTaskSEDvsEventShapes::fSparseEvtShapewithNoPid
THnSparseD * fSparseEvtShapewithNoPid
THnSparse histograms for Spherocity.
Definition: AliAnalysisTaskSEDvsEventShapes.h:190
AliRDHFCutsDplustoKpipi
Class for cuts on AOD reconstructed D+->Kpipi.
Definition: AliRDHFCutsDplustoKpipi.h:22
AliAnalysisTaskSEDvsEventShapes::fHistNtrVsnTrackEvWithCand
TH2F * fHistNtrVsnTrackEvWithCand
hist of ntracklets vs Zvertex
Definition: AliAnalysisTaskSEDvsEventShapes.h:166
AliNormalizationCounter::SetStudyMultiplicity
void SetStudyMultiplicity(Bool_t flag, Float_t etaRange)
Definition: AliNormalizationCounter.h:61
AliAnalysisTaskSEDvsEventShapes::FillMCGenAccHistos
void FillMCGenAccHistos(AliAODEvent *aod, TClonesArray *arrayMC, AliAODMCHeader *mcHeader, Double_t countMult, Double_t spherocity, Double_t sphericity, Bool_t isEvSel, Double_t nchWeight)
Definition: AliAnalysisTaskSEDvsEventShapes.cxx:1479
AliAnalysisTaskSEDvsEventShapes::fHistNtrCorrVsNchMC
TH2F * fHistNtrCorrVsNchMC
! hist of ntracklets vs Nch (Generated)
Definition: AliAnalysisTaskSEDvsEventShapes.h:173
Int_t
int Int_t
Definition: External.C:63
AliAnalysisTaskSEDvsEventShapes::fMCRecoFeeddownSpheri
THnSparseD * fMCRecoFeeddownSpheri
histo for StepMCReco for D meson feeddown for Sphericity
Definition: AliAnalysisTaskSEDvsEventShapes.h:202
UInt_t
unsigned int UInt_t
Definition: External.C:33
Float_t
float Float_t
Definition: External.C:68
AliAnalysisTaskSEDvsEventShapes::fMCRecoBothPromptFDSpheri
THnSparseD * fMCRecoBothPromptFDSpheri
histo for StepMCReco for D meson feeddown for Sphericity
Definition: AliAnalysisTaskSEDvsEventShapes.h:203
AliAnalysisTaskSEDvsEventShapes::fHistNtrVsNchMCPrimary
TH2F * fHistNtrVsNchMCPrimary
! hist of ntracklets vs Nch (Primary)
Definition: AliAnalysisTaskSEDvsEventShapes.h:174
AliAnalysisTaskSEDvsEventShapes::fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary
TH3F * fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary
! hist of Nch (generated) vs Nch (Primary) vs Nch (Physical Primary)
Definition: AliAnalysisTaskSEDvsEventShapes.h:179
AliAnalysisTaskSEDvsEventShapes::fFillTrackHisto
Bool_t fFillTrackHisto
flag for quark/hadron level identification of prompt and feeddown
Definition: AliAnalysisTaskSEDvsEventShapes.h:250
AliRDHFCuts.h
AliAnalysisTaskSEDvsEventShapes::fCounterU
AliNormalizationCounter * fCounterU
Counter for normalization, corrected multiplicity.
Definition: AliAnalysisTaskSEDvsEventShapes.h:216
AliVertexingHFUtils::CheckD0Decay
static Int_t CheckD0Decay(AliStack *stack, Int_t label, Int_t *arrayDauLab)
Definition: AliVertexingHFUtils.cxx:697
AliNormalizationCounter::SetStudySpherocity
void SetStudySpherocity(Bool_t flag, Double_t nsteps=100.)
Definition: AliNormalizationCounter.h:62
AliAnalysisTaskSEDvsEventShapes::fSparseEvtShapePrompt
THnSparseD * fSparseEvtShapePrompt
THnSparse histograms for D0 vs. Spherocity.
Definition: AliAnalysisTaskSEDvsEventShapes.h:191
AliAnalysisTaskSEDvsEventShapes::fHistNtrCorrEvSel
TH1F * fHistNtrCorrEvSel
hist. of ntracklets for physics selection only selected events
Definition: AliAnalysisTaskSEDvsEventShapes.h:182
AliAnalysisTaskSEDvsEventShapes::SetMassLimits
void SetMassLimits(Double_t lowlimit, Double_t uplimit)
Definition: AliAnalysisTaskSEDvsEventShapes.cxx:293
Data
Bool_t Data(TH1F *h, Double_t *rangefit, Bool_t writefit, Double_t &sgn, Double_t &errsgn, Double_t &bkg, Double_t &errbkg, Double_t &sgnf, Double_t &errsgnf, Double_t &sigmafit, Int_t &status)
Definition: charmCutsOptimization.C:371
AliAnalysisTaskSEDvsEventShapes::fHistNtrCorrEvWithD
TH1F * fHistNtrCorrEvWithD
hist. of ntracklets for evnts with a candidate
Definition: AliAnalysisTaskSEDvsEventShapes.h:184
AliRDHFCuts::IsEventRejectedDuePhysicsSelection
Bool_t IsEventRejectedDuePhysicsSelection() const
Definition: AliRDHFCuts.h:328
AliAnalysisTaskSEDvsEventShapes::GetEstimatorHistogram
TProfile * GetEstimatorHistogram(const AliVEvent *event)
Definition: AliAnalysisTaskSEDvsEventShapes.cxx:1275
AliAnalysisTaskSEDvsEventShapes::fHistNtrCorrEvWithCand
TH1F * fHistNtrCorrEvWithCand
hist. of ntracklets for selected events
Definition: AliAnalysisTaskSEDvsEventShapes.h:183
AliAnalysisTaskSEDvsEventShapes::fMCAccGenFeeddownSpheri
THnSparseD * fMCAccGenFeeddownSpheri
histo for StepMCGenAcc for D meson prompt for Sphericity
Definition: AliAnalysisTaskSEDvsEventShapes.h:200
ClassImp
ClassImp(AliAnalysisTaskSEDvsEventShapes) AliAnalysisTaskSEDvsEventShapes
Definition: AliAnalysisTaskSEDvsEventShapes.cxx:53
AliVertexingHFUtils::GetCorrectedNtracklets
static Double_t GetCorrectedNtracklets(TProfile *estimatorAvg, Double_t uncorrectedNacc, Double_t vtxZ, Double_t refMult)
Definition: AliVertexingHFUtils.cxx:490
AliAnalysisTaskSEDvsEventShapes::fMCRecoBothPromptFD
THnSparseD * fMCRecoBothPromptFD
histo for StepMCReco for D meson feeddown
Definition: AliAnalysisTaskSEDvsEventShapes.h:198
AliRDHFCuts::IsEventSelected
Bool_t IsEventSelected(AliVEvent *event)
Definition: AliRDHFCuts.cxx:479
AliAnalysisTaskSEDvsEventShapes::fMCAccGenPromptSpheri
THnSparseD * fMCAccGenPromptSpheri
histo for StepMCReco for D meson Both Prompt Feeddown
Definition: AliAnalysisTaskSEDvsEventShapes.h:199
AliNormalizationCounter::StoreEvent
void StoreEvent(AliVEvent *, AliRDHFCuts *, Bool_t mc=kFALSE, Int_t multiplicity=-9999, Double_t spherocity=-99.)
Definition: AliNormalizationCounter.cxx:163
AliAnalysisTaskSEDvsEventShapes::fUpmasslimit
Double_t fUpmasslimit
histograms for impact paramter studies
Definition: AliAnalysisTaskSEDvsEventShapes.h:210
AliAnalysisTaskSEDvsEventShapes::dNdptFit
Double_t dNdptFit(Float_t pt, Double_t *par)
Definition: AliAnalysisTaskSEDvsEventShapes.cxx:1331
AliNormalizationCounter.h
AliAODRecoDecayHF.h
AliAnalysisTaskSEDvsEventShapes::fMCAccGenPromptEvSel
THnSparseD * fMCAccGenPromptEvSel
histo for StepMCReco for D meson Both Prompt Feeddown for Sphericity
Definition: AliAnalysisTaskSEDvsEventShapes.h:205
AliAnalysisTaskSEDvsEventShapes::FillMCMassHistos
void FillMCMassHistos(TClonesArray *arrayMC, Int_t labD, Double_t countMult, Double_t spherocity, Double_t sphericity, Double_t recSpherocity, Double_t nchWeight)
Definition: AliAnalysisTaskSEDvsEventShapes.cxx:1427
AliAnalysisTaskSEDvsEventShapes::fHistnTrackvsEtavsPhiEvWithCand
TH3F * fHistnTrackvsEtavsPhiEvWithCand
! hist. of number of tracks passing track selection for spherocity calculation vs eta vs...
Definition: AliAnalysisTaskSEDvsEventShapes.h:187
AliAnalysisTaskSEDvsEventShapes::fPtAccCut
Double_t fPtAccCut
eta limits for acceptance step
Definition: AliAnalysisTaskSEDvsEventShapes.h:252
AliRDHFCuts::IsSelected
Bool_t IsSelected(TObject *obj)
Definition: AliRDHFCuts.h:276
AliRDHFCutsDStartoKpipi.h
Option_t
const char Option_t
Definition: External.C:48
AliAnalysisTaskSEDvsEventShapes::fListCuts
TList * fListCuts
list send on output slot 1
Definition: AliAnalysisTaskSEDvsEventShapes.h:157
AliAnalysisTaskSEDvsEventShapes::fSparseEvtShapeRecSphero
THnSparseD * fSparseEvtShapeRecSphero
THnSparse histograms for feeddown D0 vs. Spherocity.
Definition: AliAnalysisTaskSEDvsEventShapes.h:193
AliAnalysisTaskSEDvsEventShapes::fHistGenPrimaryParticlesInelGt0
TH1F * fHistGenPrimaryParticlesInelGt0
!hist. of geenrated multiplcity
Definition: AliAnalysisTaskSEDvsEventShapes.h:178
AliVertexingHFUtils::GetVZEROAEqualizedMultiplicity
static Double_t GetVZEROAEqualizedMultiplicity(AliAODEvent *ev)
Utilities for V0 multiplicity checks.
Definition: AliVertexingHFUtils.cxx:248
AliAnalysisTaskSEDvsEventShapes::FillTrackControlHisto
Bool_t FillTrackControlHisto(AliAODEvent *aod, Int_t nSelTrkCorr, Int_t nSelectedEvwithCand)
Definition: AliAnalysisTaskSEDvsEventShapes.cxx:1381
AliVertexingHFUtils::GetTrueImpactParameterDzero
static Double_t GetTrueImpactParameterDzero(AliAODMCHeader *mcHeader, TClonesArray *arrayMC, AliAODMCParticle *partDp)
Functions for computing true impact parameter of D meson.
Definition: AliVertexingHFUtils.cxx:361
nbins
const Int_t nbins
Definition: AverageDmesonRaa.C:63
Bool_t
bool Bool_t
Definition: External.C:53
ptWeight
Int_t ptWeight
Definition: ComputeEfficiencyFromCombinHF.C:26
AliAnalysisTaskSEDvsEventShapes::fHistNtrCorrVsNchMCPrimary
TH2F * fHistNtrCorrVsNchMCPrimary
! hist of ntracklets vs Nch (Primary)
Definition: AliAnalysisTaskSEDvsEventShapes.h:175
AliAODRecoCascadeHF::Get2Prong
AliAODRecoDecayHF2Prong * Get2Prong() const
Definition: AliAODRecoCascadeHF.h:43
AliAnalysisTaskSEDvsEventShapes::fMCAccGenPrompt
THnSparseD * fMCAccGenPrompt
THnSparse histograms for Both Prompt and feeddown D0 vs. Spherocity.
Definition: AliAnalysisTaskSEDvsEventShapes.h:194
AliVertexingHFUtils::GetTrueImpactParameterDplus
static Double_t GetTrueImpactParameterDplus(AliAODMCHeader *mcHeader, TClonesArray *arrayMC, AliAODMCParticle *partDp)
Definition: AliVertexingHFUtils.cxx:404
AliAnalysisTaskSEDvsEventShapes::fEtaAccCut
Double_t fEtaAccCut
flag for filling track control histograms
Definition: AliAnalysisTaskSEDvsEventShapes.h:251
AliRDHFCuts::IsInFiducialAcceptance
virtual Bool_t IsInFiducialAcceptance(Double_t, Double_t) const
Definition: AliRDHFCuts.h:296
AliVertexingHFUtils::GetSphericity
static Double_t GetSphericity(AliAODEvent *aod, Double_t etaMin=-0.8, Double_t etaMax=0.8, Double_t ptMin=0.15, Double_t ptMax=10., Int_t filtbit1=256, Int_t filtbit2=512, Int_t minMult=3)
Definition: AliVertexingHFUtils.cxx:1917
AliAnalysisTaskSEDvsEventShapes::fHistNtrVsSpheri
TH2F * fHistNtrVsSpheri
hist of ntracklets vs So
Definition: AliAnalysisTaskSEDvsEventShapes.h:169
AliAnalysisTaskSEDvsEventShapes::fHistNEvents
TH1F * fHistNEvents
list send on output slot 5
Definition: AliAnalysisTaskSEDvsEventShapes.h:162
AliVertexingHFUtils::GetSpherocity
static Double_t GetSpherocity(AliAODEvent *aod, Double_t etaMin=-0.8, Double_t etaMax=0.8, Double_t ptMin=0.15, Double_t ptMax=10., Int_t filtbit1=256, Int_t filtbit2=512, Int_t minMult=3, Double_t phiStepSizeDeg=0.1, Int_t nTrksToSkip=0, Int_t *idToSkip=0x0)
Functions for event shape variables.
Definition: AliVertexingHFUtils.cxx:1973