AliPhysics  v5-07-15-01 (b3d7633)
 All Classes Namespaces Files Functions Variables Enumerations Enumerator Macros
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 AliAnalysisTaskSEDvsEventShape
20 // AliAnalysisTaskSE for the D meson vs. Event shape analysis in different mutiplicity window
21 // Authors: Renu Bala, Manoj Bhanudas Jadhav
23 
24 #include <TClonesArray.h>
25 #include <TCanvas.h>
26 #include <TList.h>
27 #include <TString.h>
28 #include <TDatabasePDG.h>
29 #include <TH1F.h>
30 #include <TH2F.h>
31 #include <TH3F.h>
32 #include <THnSparse.h>
33 #include <TProfile.h>
34 #include "AliAnalysisManager.h"
35 #include "AliRDHFCuts.h"
38 #include "AliRDHFCutsD0toKpi.h"
39 #include "AliAODHandler.h"
40 #include "AliAODEvent.h"
41 #include "AliAODVertex.h"
42 #include "AliAODTrack.h"
43 #include "AliAODRecoDecayHF.h"
44 #include "AliAODRecoCascadeHF.h"
45 #include "AliAnalysisVertexingHF.h"
46 #include "AliAnalysisTaskSE.h"
49 #include "AliVertexingHFUtils.h"
50 #include "AliAODVZERO.h"
51 #include "AliESDUtils.h"
53 
54 //________________________________________________________________________
56 AliAnalysisTaskSE(),
57 fOutput(0),
58 fListCuts(0),
59 fOutputCounters(0),
60 fListProfiles(0),
61 fHistNEvents(0),
62 fHistNtrEta16vsNtrEta1EvSel(0),
63 fHistNtrEta05vsNtrEta1EvSel(0),
64 fHistNtrEta03vsNtrEta1EvSel(0),
65 fHistNtrEtaV0AvsNtrEta1EvSel(0),
66 fHistNtrEtaV0MvsNtrEta1EvSel(0),
67 fHistNtrEtaV0AvsV0AEqEvSel(0),
68 fHistNtrEtaV0MvsV0MEqEvSel(0),
69 fHistNtrCorrEta1vsNtrRawEta1EvSel(0),
70 fHistMultCorrvsMultRawEvSel(0),
71 fHistNtrEta16vsNtrEta1EvWithCand(0),
72 fHistNtrEta05vsNtrEta1EvWithCand(0),
73 fHistNtrEta03vsNtrEta1EvWithCand(0),
74 fHistNtrEtaV0AvsNtrEta1EvWithCand(0),
75 fHistNtrEtaV0MvsNtrEta1EvWithCand(0),
76 fHistNtrEtaV0AvsV0AEqEvWithCand(0),
77 fHistNtrEtaV0MvsV0MEqEvWithCand(0),
78 fHistNtrCorrEta1vsNtrRawEta1EvWithCand(0),
79 fHistMultCorrvsMultRawEvWithCand(0),
80 fHistNtrEta16vsNtrEta1EvWithD(0),
81 fHistNtrEta05vsNtrEta1EvWithD(0),
82 fHistNtrEta03vsNtrEta1EvWithD(0),
83 fHistNtrEtaV0AvsNtrEta1EvWithD(0),
84 fHistNtrEtaV0MvsNtrEta1EvWithD(0),
85 fHistNtrEtaV0AvsV0AEqEvWithD(0),
86 fHistNtrEtaV0MvsV0MEqEvWithD(0),
87 fHistNtrCorrEta1vsNtrRawEta1EvWithD(0),
88 fHistMultCorrvsMultRawEvWithD(0),
89 fHistNtrVsZvtx(0),
90 fHistNtrCorrVsZvtx(0),
91 fHistNtrVsSo(0),
92 fHistNtrCorrVsSo(0),
93 fHistNtrVsNchMC(0),
94 fHistNtrCorrVsNchMC(0),
95 fHistNtrVsNchMCPrimary(0),
96 fHistNtrCorrVsNchMCPrimary(0),
97 fHistNtrVsNchMCPhysicalPrimary(0),
98 fHistNtrCorrVsNchMCPhysicalPrimary(0),
99 fHistGenPrimaryParticlesInelGt0(0),
100 fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary(0),
101 fHistNtrUnCorrPSSel(0),
102 fHistNtrUnCorrPSTrigSel(0),
103 fHistNtrUnCorrPSTrigPileUpSel(0),
104 fHistNtrUnCorrPSTrigPileUpVtxSel(0),
105 fHistNtrUnCorrPSTrigPileUpVtxContSel(0),
106 fHistNtrUnCorrPSTrigPileUpVtxRangeSel(0),
107 fHistNtrUnCorrPSTrigPileUpVtxRangeCentrSel(0),
108 fHistNtrUnCorrEvSel(0),
109 fHistNtrUnCorrEvWithCand(0),
110 fHistNtrUnCorrEvWithD(0),
111 fHistNtrCorrPSSel(0),
112 fHistNtrCorrEvSel(0),
113 fHistNtrCorrEvWithCand(0),
114 fHistNtrCorrEvWithD(0),
115 fSparseSpherocity(0),
116 fSparseSpherocitywithNoPid(0),
117 fPtVsMassVsMultMC(0),
118 fUpmasslimit(1.965),
119 fLowmasslimit(1.765),
120 fNMassBins(200),
121 fRDCutsAnalysis(0),
122 fCounterC(0),
123 fCounterU(0),
124 fCounterCandidates(0),
125 fDoImpPar(kFALSE),
126 fNImpParBins(400),
127 fLowerImpPar(-2000.),
128 fHigherImpPar(2000.),
129 fReadMC(kFALSE),
130 fMCOption(0),
131 fisPPbData(kFALSE),
132 fUseBit(kTRUE),
133 fSubtractTrackletsFromDau(kFALSE),
134 fKeepCorrPlots(kFALSE),
135 fUseNchWeight(0),
136 fHistoMCNch(0),
137 fHistoMeasNch(0),
138 fRefMult(9.26),
139 fPdgMeson(411),
140 fMultiplicityEstimator(kNtrk10),
141 fMCPrimariesEstimator(kEta10),
142 fFillSoSparseChecks(0),
143 fptMin(0.15),
144 fptMax(10.),
145 fetaMin(-0.8),
146 fetaMax(0.8),
147 ffiltbit1(256),
148 ffiltbit2(512),
149 fminMult(3),
150 fphiStepSizeDeg(0.1),
151 fCalculateSphericity(kFALSE),
152 fDoVZER0ParamVertexCorr(1)
153 {
154  // Default constructor
155  for(Int_t i=0; i<5; i++) fHistMassPtImpPar[i]=0;
156  for(Int_t i=0; i<4; i++) fMultEstimatorAvg[i]=0;
157 }
158 
159 //________________________________________________________________________
160 AliAnalysisTaskSEDvsEventShapes::AliAnalysisTaskSEDvsEventShapes(const char *name, Int_t pdgMeson,AliRDHFCuts *cuts, Bool_t switchPPb):
161 AliAnalysisTaskSE(name),
162 fOutput(0),
163 fListCuts(0),
164 fOutputCounters(0),
165 fListProfiles(0),
166 fHistNEvents(0),
167 fHistNtrEta16vsNtrEta1EvSel(0),
168 fHistNtrEta05vsNtrEta1EvSel(0),
169 fHistNtrEta03vsNtrEta1EvSel(0),
170 fHistNtrEtaV0AvsNtrEta1EvSel(0),
171 fHistNtrEtaV0MvsNtrEta1EvSel(0),
172 fHistNtrEtaV0AvsV0AEqEvSel(0),
173 fHistNtrEtaV0MvsV0MEqEvSel(0),
174 fHistNtrCorrEta1vsNtrRawEta1EvSel(0),
175 fHistMultCorrvsMultRawEvSel(0),
176 fHistNtrEta16vsNtrEta1EvWithCand(0),
177 fHistNtrEta05vsNtrEta1EvWithCand(0),
178 fHistNtrEta03vsNtrEta1EvWithCand(0),
179 fHistNtrEtaV0AvsNtrEta1EvWithCand(0),
180 fHistNtrEtaV0MvsNtrEta1EvWithCand(0),
181 fHistNtrEtaV0AvsV0AEqEvWithCand(0),
182 fHistNtrEtaV0MvsV0MEqEvWithCand(0),
183 fHistNtrCorrEta1vsNtrRawEta1EvWithCand(0),
184 fHistMultCorrvsMultRawEvWithCand(0),
185 fHistNtrEta16vsNtrEta1EvWithD(0),
186 fHistNtrEta05vsNtrEta1EvWithD(0),
187 fHistNtrEta03vsNtrEta1EvWithD(0),
188 fHistNtrEtaV0AvsNtrEta1EvWithD(0),
189 fHistNtrEtaV0MvsNtrEta1EvWithD(0),
190 fHistNtrEtaV0AvsV0AEqEvWithD(0),
191 fHistNtrEtaV0MvsV0MEqEvWithD(0),
192 fHistNtrCorrEta1vsNtrRawEta1EvWithD(0),
193 fHistMultCorrvsMultRawEvWithD(0),
194 fHistNtrVsZvtx(0),
195 fHistNtrCorrVsZvtx(0),
196 fHistNtrVsSo(0),
197 fHistNtrCorrVsSo(0),
198 fHistNtrVsNchMC(0),
199 fHistNtrCorrVsNchMC(0),
200 fHistNtrVsNchMCPrimary(0),
201 fHistNtrCorrVsNchMCPrimary(0),
202 fHistNtrVsNchMCPhysicalPrimary(0),
203 fHistNtrCorrVsNchMCPhysicalPrimary(0),
204 fHistGenPrimaryParticlesInelGt0(0),
205 fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary(0),
206 fHistNtrUnCorrPSSel(0),
207 fHistNtrUnCorrPSTrigSel(0),
208 fHistNtrUnCorrPSTrigPileUpSel(0),
209 fHistNtrUnCorrPSTrigPileUpVtxSel(0),
210 fHistNtrUnCorrPSTrigPileUpVtxContSel(0),
211 fHistNtrUnCorrPSTrigPileUpVtxRangeSel(0),
212 fHistNtrUnCorrPSTrigPileUpVtxRangeCentrSel(0),
213 fHistNtrUnCorrEvSel(0),
214 fHistNtrUnCorrEvWithCand(0),
215 fHistNtrUnCorrEvWithD(0),
216 fHistNtrCorrPSSel(0),
217 fHistNtrCorrEvSel(0),
218 fHistNtrCorrEvWithCand(0),
219 fHistNtrCorrEvWithD(0),
220 fSparseSpherocity(0),
221 fSparseSpherocitywithNoPid(0),
222 fPtVsMassVsMultMC(0),
223 fUpmasslimit(1.965),
224 fLowmasslimit(1.765),
225 fNMassBins(200),
226 fRDCutsAnalysis(cuts),
227 fCounterC(0),
228 fCounterU(0),
229 fCounterCandidates(0),
230 fDoImpPar(kFALSE),
231 fNImpParBins(400),
232 fLowerImpPar(-2000.),
233 fHigherImpPar(2000.),
234 fReadMC(kFALSE),
235 fMCOption(0),
236 fisPPbData(switchPPb),
237 fUseBit(kTRUE),
238 fSubtractTrackletsFromDau(kFALSE),
239 fKeepCorrPlots(kFALSE),
240 fUseNchWeight(0),
241 fHistoMCNch(0),
242 fHistoMeasNch(0),
243 fRefMult(9.26),
244 fPdgMeson(pdgMeson),
245 fMultiplicityEstimator(kNtrk10),
246 fMCPrimariesEstimator(kEta10),
247 fFillSoSparseChecks(0),
248 fptMin(0.15),
249 fptMax(10.),
250 fetaMin(-0.8),
251 fetaMax(0.8),
252 ffiltbit1(256),
253 ffiltbit2(512),
254 fminMult(3),
255 fphiStepSizeDeg(0.1),
256 fCalculateSphericity(kFALSE),
257 fDoVZER0ParamVertexCorr(1)
258 {
259  //
260  // Standard constructor
261  //
262 
263  for(Int_t i=0; i<5; i++) fHistMassPtImpPar[i]=0;
264  for(Int_t i=0; i<4; i++) fMultEstimatorAvg[i]=0;
265  if(fPdgMeson==413){
266  fNMassBins=200;
267  SetMassLimits(0.12,0.2);
268  }else{
269  fNMassBins=200;
271  }
272  // Default constructor
273  // Otput slot #1 writes into a TList container
274  DefineOutput(1,TList::Class()); //My private output
275  // Output slot #2 writes cut to private output
276  DefineOutput(2,TList::Class());
277  // Output slot #3 writes cut to private output
278  DefineOutput(3,TList::Class());
279  // Output slot #4 writes cut to private output
280  DefineOutput(4,TList::Class());
281 }
282 //________________________________________________________________________
284 {
285  //
286  // Destructor
287  //
288  delete fOutput;
289  delete fHistNEvents;
290  delete fListCuts;
291  delete fListProfiles;
292  delete fRDCutsAnalysis;
293  delete fCounterC;
294  delete fCounterU;
295  delete fCounterCandidates;
296 
297  for(Int_t i=0; i<4; i++) {
298  if (fMultEstimatorAvg[i]) delete fMultEstimatorAvg[i];
299  }
300 
301  for(Int_t i=0; i<5; i++){
302  delete fHistMassPtImpPar[i];
303  }
304  if(fHistoMCNch) delete fHistoMCNch;
305  if(fHistoMeasNch) delete fHistoMeasNch;
306 }
307 
308 //_________________________________________________________________
309 void AliAnalysisTaskSEDvsEventShapes::SetMassLimits(Double_t lowlimit, Double_t uplimit){
310  // set invariant mass limits
311  if(uplimit>lowlimit){
312  fLowmasslimit = lowlimit;
313  fUpmasslimit = uplimit;
314  }else{
315  AliError("Wrong mass limits: upper value should be larger than lower one");
316  }
317 }
318 //_________________________________________________________________
320  // set invariant mass limits
321  Double_t mass=TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdg))->Mass();
322  SetMassLimits(mass-range,mass+range);
323 }
324 //________________________________________________________________________
326  //
327  // Initialization
328  //
329  printf("AnalysisTaskSEDvsMultiplicity_0::Init() \n");
330 
331  if(fUseNchWeight && !fReadMC){ AliFatal("Nch weights can only be used in MC mode"); return; }
332  if(fUseNchWeight && !fHistoMCNch){ AliFatal("Nch weights can only be used without histogram"); return; }
333  if(fUseNchWeight==1 && !fHistoMeasNch) {//Nch weights
334  if(fisPPbData){ AliFatal("Nch weights can only be used with MC and data histogram in pPb"); return; }
335  else CreateMeasuredNchHisto();
336  }
337  if(fUseNchWeight==2 && !fHistoMeasNch){ AliFatal("Ntrk weights can only be used with MC and data histogram"); return; } //for pp, pPb Ntrk weights
338 
339  fListCuts=new TList();
340  fListCuts->SetOwner();
341  fListCuts->SetName("CutsList");
342 
343  if(fPdgMeson==411){
344  AliRDHFCutsDplustoKpipi* copycut=new AliRDHFCutsDplustoKpipi(*(static_cast<AliRDHFCutsDplustoKpipi*>(fRDCutsAnalysis)));
345  copycut->SetName("AnalysisCutsDplus");
346  fListCuts->Add(copycut);
347  }else if(fPdgMeson==421){
348  AliRDHFCutsD0toKpi* copycut=new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fRDCutsAnalysis)));
349  copycut->SetName("AnalysisCutsDzero");
350  fListCuts->Add(copycut);
351  }else if(fPdgMeson==413){
352  AliRDHFCutsDStartoKpipi* copycut=new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fRDCutsAnalysis)));
353  copycut->SetName("AnalysisCutsDStar");
354  fListCuts->Add(copycut);
355  }
358 
359  PostData(2,fListCuts);
360 
361  fListProfiles = new TList();
362  fListProfiles->SetOwner();
363  TString period[4];
364  Int_t nProfiles=4;
365  if (fisPPbData) {period[0]="LHC13b"; period[1]="LHC13c"; nProfiles = 2;}
366  else {period[0]="LHC10b"; period[1]="LHC10c"; period[2]="LHC10d"; period[3]="LHC10e"; nProfiles = 4;}
367 
368  for(Int_t i=0; i<nProfiles; i++){
369  if(fMultEstimatorAvg[i]){
370  TProfile* hprof=new TProfile(*fMultEstimatorAvg[i]);
371  hprof->SetName(Form("ProfileTrkVsZvtx%s\n",period[i].Data()));
372  fListProfiles->Add(hprof);
373  }
374  }
375 
376  PostData(4,fListProfiles);
377 
378  return;
379 }
380 
381 //________________________________________________________________________
383 {
384  // Create the output container
385  //
386  if(fDebug > 1) printf("AnalysisTaskSEDvsMultiplicity::UserCreateOutputObjects() \n");
387 
388  // Several histograms are more conveniently managed in a TList
389  fOutput = new TList();
390  fOutput->SetOwner();
391  fOutput->SetName("OutputHistos");
392 
393  Int_t nMultBins = 200;
394  Float_t firstMultBin = -0.5;
395  Float_t lastMultBin = 199.5;
396  Int_t nMultBinsNtrk = nMultBins;
397  Float_t lastMultBinNtrk = lastMultBin;
398  Int_t nMultBinsV0 = 400;
399  Float_t lastMultBinV0 = 799.5;
400  const char *estimatorName="tracklets";
401  if(fisPPbData) {
402  nMultBinsNtrk = 375;
403  lastMultBinNtrk = 374.5;
404  nMultBins = nMultBinsNtrk;
405  lastMultBin = lastMultBinNtrk;
406  }
408  nMultBins = nMultBinsV0;
409  lastMultBin = lastMultBinV0;
410  estimatorName = "vzero";
411  }
412 
413  fHistNtrUnCorrPSSel = new TH1F("hNtrUnCorrPSSel",Form("Uncorrected %s multiplicity for PS selected events; %s ; Entries",estimatorName,estimatorName),nMultBins,firstMultBin,lastMultBin);
414  fHistNtrUnCorrPSTrigSel = new TH1F("hNtrUnCorrPSTrigSel",Form("Uncorrected %s multiplicity for PS + trigger name selected events; %s ; Entries",estimatorName,estimatorName),nMultBins,firstMultBin,lastMultBin);
415  fHistNtrUnCorrPSTrigPileUpSel = new TH1F("hNtrUnCorrPSTrigPileUpSel",Form("Uncorrected %s multiplicity for PS + trigger name + pileup selected events; %s ; Entries",estimatorName,estimatorName),nMultBins,firstMultBin,lastMultBin);
416  fHistNtrUnCorrPSTrigPileUpVtxSel = new TH1F("hNtrUnCorrPSTrigPileUpVtxSel",Form("Uncorrected %s multiplicity for PS + trigger name + pileup + with-vertex selected events; %s ; Entries",estimatorName,estimatorName),nMultBins,firstMultBin,lastMultBin);
417  fHistNtrUnCorrPSTrigPileUpVtxContSel = new TH1F("hNtrUnCorrPSTrigPileUpVtxContSel",Form("Uncorrected %s multiplicity for PS + trigger name + pileup + with-vertex-contrib selected events; %s ; Entries",estimatorName,estimatorName),nMultBins,firstMultBin,lastMultBin);
418  fHistNtrUnCorrPSTrigPileUpVtxRangeSel = new TH1F("hNtrUnCorrPSTrigPileUpVtxRangeSel",Form("Uncorrected %s multiplicity for PS + trigger name + pileup + with-vertex-contrib-range selected events; %s ; Entries",estimatorName,estimatorName),nMultBins,firstMultBin,lastMultBin);
419  fHistNtrUnCorrPSTrigPileUpVtxRangeCentrSel = new TH1F("hNtrUnCorrPSTrigPileUpVtxRangeCentrSel",Form("Uncorrected %s multiplicity for PS + trigger name + pileup + with-vertex-contrib-range + centrality selected events; %s ; Entries",estimatorName,estimatorName),nMultBins,firstMultBin,lastMultBin);
420  fHistNtrUnCorrEvSel = new TH1F("hNtrUnCorrEvSel",Form("Uncorrected %s multiplicity for selected events; %s ; Entries",estimatorName,estimatorName),nMultBins,firstMultBin,lastMultBin);
421  fHistNtrUnCorrEvWithCand = new TH1F("hNtrUnCorrEvWithCand",Form("Uncorrected %s multiplicity for events with D candidates; %s ; Entries",estimatorName,estimatorName),nMultBins,firstMultBin,lastMultBin);// Total multiplicity
422  fHistNtrUnCorrEvWithD = new TH1F("hNtrUnCorrEvWithD",Form("Uncorrected %s multiplicity for events with D in mass region ; %s ; Entries",estimatorName,estimatorName),nMultBins,firstMultBin,lastMultBin); //
423  fHistNtrCorrPSSel = new TH1F("hNtrCorrPSSel",Form("Corrected %s multiplicity for PS selected events; %s ; Entries",estimatorName,estimatorName),nMultBins,firstMultBin,lastMultBin);
424  fHistNtrCorrEvSel = new TH1F("hNtrCorrEvSel",Form("Corrected %s multiplicity for selected events; %s ; Entries",estimatorName,estimatorName),nMultBins,firstMultBin,lastMultBin);
425  fHistNtrCorrEvWithCand = new TH1F("hNtrCorrEvWithCand", Form("%s multiplicity for events with D candidates; %s ; Entries",estimatorName,estimatorName),nMultBins,firstMultBin,lastMultBin);// Total multiplicity
426  fHistNtrCorrEvWithD = new TH1F("hNtrCorrEvWithD", Form("%s multiplicity for events with D in mass region ; %s ; Entries",estimatorName,estimatorName),nMultBins,firstMultBin,lastMultBin); //
427 
428  if(fKeepCorrPlots){
429  fHistNtrEta16vsNtrEta1EvSel = new TH2F("hNtrEta16vsNtrEta1EvSel","Uncorrected Eta1.6 vs Eta1.0 (events selected); Ntracklets #eta<1.0; Ntracklets #eta<1.6",nMultBinsNtrk,firstMultBin,lastMultBinNtrk,nMultBinsNtrk,firstMultBin,lastMultBinNtrk); //eta 1.6 vs eta 1.0 histogram
430  fHistNtrEta05vsNtrEta1EvSel = new TH2F("hNtrEta05vsNtrEta1EvSel","Uncorrected Eta0.5 vs Eta1.0 (events selected); Ntracklets #eta<1.0; Ntracklets #eta<0.5",nMultBinsNtrk,firstMultBin,lastMultBinNtrk,nMultBinsNtrk,firstMultBin,lastMultBinNtrk); //eta 0.5 vs eta 1.0 histogram
431  fHistNtrEta03vsNtrEta1EvSel = new TH2F("hNtrEta03vsNtrEta1EvSel","Uncorrected Eta0.3 vs Eta1.0 (events selected); Ntracklets #eta<1.0; Ntracklets #eta<0.3",nMultBinsNtrk,firstMultBin,lastMultBinNtrk,nMultBinsNtrk,firstMultBin,lastMultBinNtrk); //eta 0.3 vs eta 1.0 histogram
432  fHistNtrEtaV0AvsNtrEta1EvSel = new TH2F("hNtrEtaV0AvsNtrEta1EvSel","Uncorrected Eta-V0A vs Eta1.0 (events selected); Ntracklets #eta<1.0; Multiplicity V0A",nMultBinsNtrk,firstMultBin,lastMultBinNtrk,nMultBinsV0,firstMultBin,lastMultBinV0); //eta V0A vs eta 1.0 histogram
433  fHistNtrEtaV0MvsNtrEta1EvSel = new TH2F("hNtrEtaV0MvsNtrEta1EvSel","Uncorrected Eta-V0M vs Eta1.0 (events selected); Ntracklets #eta<1.0; Multiplicity V0A+V0C",nMultBinsNtrk,firstMultBin,lastMultBinNtrk,nMultBinsV0,firstMultBin,lastMultBinV0); //eta V0M vs eta 1.0 histogram
434  fHistNtrEtaV0AvsV0AEqEvSel = new TH2F("hNtrEtaV0AvsV0AEqEvSel","Corrected V0A vs corrected V0A-Equalized (events selected); Vzero-A; Vzero-A Equalized",nMultBinsV0,firstMultBin,lastMultBinV0,nMultBinsV0,firstMultBin,lastMultBinV0); // comparison V0A - V0Aeq
435  fHistNtrEtaV0MvsV0MEqEvSel = new TH2F("hNtrEtaV0MvsV0MEqEvSel","Corrected V0M vs corrected V0M-Equalized (events selected); Vzero-M; Vzero-M Equalized",nMultBinsV0,firstMultBin,lastMultBinV0,nMultBinsV0,firstMultBin,lastMultBinV0); // comparison V0M - V0Meq
436  fHistNtrCorrEta1vsNtrRawEta1EvSel = new TH2F("hNtrCorrEta1vsNtrRawEta1EvSel","Corrected Eta1 vs Eta1.0 (events selected); Ntracklets #eta<1.0 corrected; Ntracklets #eta<1",nMultBinsNtrk,firstMultBin,lastMultBinNtrk,nMultBinsNtrk,firstMultBin,lastMultBinNtrk); //eta 1.6 vs eta 1.0 histogram
437  fHistMultCorrvsMultRawEvSel = new TH2F("hMultCorrvsMultRawEvSel",Form("Corrected multiplicity vs uncorrected multiplicity (events selected); %s corrected; %s",estimatorName,estimatorName),nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); // corrected vs uncorrected multiplicity
438 
439  fHistNtrEta16vsNtrEta1EvWithCand = new TH2F("hNtrEta16vsNtrEta1EvWithCand","Uncorrected Eta1.6 vs Eta1.0 (events selected with a D candidate); Ntracklets #eta<1.0; Ntracklets #eta<1.6",nMultBinsNtrk,firstMultBin,lastMultBinNtrk,nMultBinsNtrk,firstMultBin,lastMultBinNtrk); //eta 1.6 vs eta 1.0 histogram
440  fHistNtrEta05vsNtrEta1EvWithCand = new TH2F("hNtrEta05vsNtrEta1EvWithCand","Uncorrected Eta0.5 vs Eta1.0 (events selected with a D candidate); Ntracklets #eta<1.0; Ntracklets #eta<0.5",nMultBinsNtrk,firstMultBin,lastMultBinNtrk,nMultBinsNtrk,firstMultBin,lastMultBinNtrk); //eta 0.5 vs eta 1.0 histogram
441  fHistNtrEta03vsNtrEta1EvWithCand = new TH2F("hNtrEta03vsNtrEta1EvWithCand","Uncorrected Eta0.3 vs Eta1.0 (events selected with a D candidate); Ntracklets #eta<1.0; Ntracklets #eta<0.3",nMultBinsNtrk,firstMultBin,lastMultBinNtrk,nMultBinsNtrk,firstMultBin,lastMultBinNtrk); //eta 0.3 vs eta 1.0 histogram
442  fHistNtrEtaV0AvsNtrEta1EvWithCand = new TH2F("hNtrEtaV0AvsNtrEta1EvWithCand","Uncorrected Eta-V0A vs Eta1.0 (events selected with a D candidate); Ntracklets #eta<1.0; Multiplicity V0A",nMultBinsNtrk,firstMultBin,lastMultBinNtrk,nMultBinsV0,firstMultBin,lastMultBinV0); //eta V0A vs eta 1.0 histogram
443  fHistNtrEtaV0MvsNtrEta1EvWithCand = new TH2F("hNtrEtaV0MvsNtrEta1EvWithCand","Uncorrected Eta-V0M vs Eta1.0 (events selected with a D candidate); Ntracklets #eta<1.0; Multiplicity V0A+V0C",nMultBinsNtrk,firstMultBin,lastMultBinNtrk,nMultBinsV0,firstMultBin,lastMultBinV0); //eta V0M vs eta 1.0 histogram
444  fHistNtrEtaV0AvsV0AEqEvWithCand = new TH2F("hNtrEtaV0AvsV0AEqEvWithCand","Corrected V0A vs corrected V0A-Equalized (events selected with a D candidate); Vzero-A; Vzero-A Equalized",nMultBinsV0,firstMultBin,lastMultBinV0,nMultBinsV0,firstMultBin,lastMultBinV0); // comparison V0A - V0Aeq
445  fHistNtrEtaV0MvsV0MEqEvWithCand = new TH2F("hNtrEtaV0MvsV0MEqEvWithCand","Corrected V0M vs corrected V0M-Equalized (events selected with a D candidate); Vzero-M; Vzero-M Equalized",nMultBinsV0,firstMultBin,lastMultBinV0,nMultBinsV0,firstMultBin,lastMultBinV0); // comparison V0M - V0Meq
446  fHistNtrCorrEta1vsNtrRawEta1EvWithCand = new TH2F("hNtrCorrEta1vsNtrRawEta1EvWithCand","Corrected Eta1 vs Eta1.0 (events selected with a D candidate); Ntracklets #eta<1.0 corrected; Ntracklets #eta<1",nMultBinsNtrk,firstMultBin,lastMultBinNtrk,nMultBinsNtrk,firstMultBin,lastMultBinNtrk); //eta 1.6 vs eta 1.0 histogram
447  fHistMultCorrvsMultRawEvWithCand = new TH2F("hMultCorrvsMultRawEvWithCand",Form("Corrected multiplicity vs uncorrected multiplicity (events selected) with a D candidate; %s corrected; %s",estimatorName,estimatorName),nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); // corrected vs uncorrected multiplicity
448 
449 
450  fHistNtrEta16vsNtrEta1EvWithD = new TH2F("hNtrEta16vsNtrEta1EvWithD","Uncorrected Eta1.6 vs Eta1.0 (events selected with D in mass range); Ntracklets #eta<1.0; Ntracklets #eta<1.6",nMultBinsNtrk,firstMultBin,lastMultBinNtrk,nMultBinsNtrk,firstMultBin,lastMultBinNtrk); //eta 1.6 vs eta 1.0 histogram
451  fHistNtrEta05vsNtrEta1EvWithD = new TH2F("hNtrEta05vsNtrEta1EvWithD","Uncorrected Eta0.5 vs Eta1.0 (events selected with D in mass range); Ntracklets #eta<1.0; Ntracklets #eta<0.5",nMultBinsNtrk,firstMultBin,lastMultBinNtrk,nMultBinsNtrk,firstMultBin,lastMultBinNtrk); //eta 0.5 vs eta 1.0 histogram
452  fHistNtrEta03vsNtrEta1EvWithD = new TH2F("hNtrEta03vsNtrEta1EvWithD","Uncorrected Eta0.3 vs Eta1.0 (events selected with D in mass range); Ntracklets #eta<1.0; Ntracklets #eta<0.3",nMultBinsNtrk,firstMultBin,lastMultBinNtrk,nMultBinsNtrk,firstMultBin,lastMultBinNtrk); //eta 0.3 vs eta 1.0 histogram
453  fHistNtrEtaV0AvsNtrEta1EvWithD = new TH2F("hNtrEtaV0AvsNtrEta1EvWithD","Uncorrected Eta-V0A vs Eta1.0 (events selected with D in mass range); Ntracklets #eta<1.0; Multiplicity V0A",nMultBinsNtrk,firstMultBin,lastMultBinNtrk,nMultBinsV0,firstMultBin,lastMultBinV0); //eta V0A vs eta 1.0 histogram
454  fHistNtrEtaV0MvsNtrEta1EvWithD = new TH2F("hNtrEtaV0MvsNtrEta1EvWithD","Uncorrected Eta-V0M vs Eta1.0 (events selected with D in mass range); Ntracklets #eta<1.0; Multiplicity V0A+V0C",nMultBinsNtrk,firstMultBin,lastMultBinNtrk,nMultBinsV0,firstMultBin,lastMultBinV0); //eta V0M vs eta 1.0 histogram
455  fHistNtrEtaV0AvsV0AEqEvWithD = new TH2F("hNtrEtaV0AvsV0AEqEvWithD","Corrected V0A vs corrected V0A-Equalized (events selected with D in mass range); Vzero-A; Vzero-A Equalized",nMultBinsV0,firstMultBin,lastMultBinV0,nMultBinsV0,firstMultBin,lastMultBinV0); // comparison V0A - V0Aeq
456  fHistNtrEtaV0MvsV0MEqEvWithD = new TH2F("hNtrEtaV0MvsV0MEqEvWithD","Corrected V0M vs corrected V0M-Equalized (events selected with D in mass range); Vzero-M; Vzero-M Equalized",nMultBinsV0,firstMultBin,lastMultBinV0,nMultBinsV0,firstMultBin,lastMultBinV0); // comparison V0M - V0Meq
457  fHistNtrCorrEta1vsNtrRawEta1EvWithD = new TH2F("hNtrCorrEta1vsNtrRawEta1EvWithD","Corrected Eta1 vs Eta1.0 (events selected with D in mass range); Ntracklets #eta<1.0 corrected; Ntracklets #eta<1",nMultBinsNtrk,firstMultBin,lastMultBinNtrk,nMultBinsNtrk,firstMultBin,lastMultBinNtrk); //eta 1.6 vs eta 1.0 histogram
458  fHistMultCorrvsMultRawEvWithD = new TH2F("hMultCorrvsMultRawEvWithD",Form("Corrected multiplicity vs uncorrected multiplicity (events selected with D in mass range); %s corrected; %s",estimatorName,estimatorName),nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); // corrected vs uncorrected multiplicity
459 
460  }
461  fHistNtrVsZvtx = new TH2F("hNtrVsZvtx",Form("N%s vs VtxZ; VtxZ;N_{%s};",estimatorName,estimatorName),300,-15,15,nMultBins,firstMultBin,lastMultBin); //
462  fHistNtrCorrVsZvtx = new TH2F("hNtrCorrVsZvtx",Form("N%s vs VtxZ; VtxZ;N_{%s};",estimatorName,estimatorName),300,-15,15,nMultBins,firstMultBin,lastMultBin); //
463 
464  TString HistoNtrName;
465  TString HistoNtrCorrName;
466  TString ParNameNtr;
468  HistoNtrName = "hNtrVsSpheri";
469  HistoNtrCorrName = "hNtrCorrVsSpheri";
470  ParNameNtr = "Spheri";
471  }
472  else{
473  HistoNtrName = "hNtrVsSphero";
474  HistoNtrCorrName = "hNtrCorrVsSphero";
475  ParNameNtr = "Sphero";
476  }
477 
478  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); //
479  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); //
480 
481  fHistNtrVsNchMC = new TH2F("hNtrVsNchMC",Form("N%s vs NchMC; Nch;N_{%s};",estimatorName,estimatorName),nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); //
482  fHistNtrCorrVsNchMC = new TH2F("hNtrCorrVsNchMC",Form("N%s vs Nch; Nch;N_{%s};",estimatorName,estimatorName),nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); //
483 
484  fHistNtrVsNchMCPrimary = new TH2F("hNtrVsNchMCPrimary",Form("N%s vs Nch (Primary); Nch (Primary);N_{%s};",estimatorName,estimatorName),nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); //
485  fHistNtrCorrVsNchMCPrimary = new TH2F("hNtrCorrVsNchMCPrimary",Form("N%s vs Nch (Primary); Nch(Primary) ;N_{%s};",estimatorName,estimatorName),nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); //
486 
487  fHistNtrVsNchMCPhysicalPrimary = new TH2F("hNtrVsNchMCPhysicalPrimary",Form("N%s vs Nch (Physical Primary); Nch (Physical Primary);N_{%s};",estimatorName,estimatorName),nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); //
488  fHistNtrCorrVsNchMCPhysicalPrimary = new TH2F("hNtrCorrVsMCPhysicalPrimary",Form("N%s vs Nch (Physical Primary); Nch (Physical Primary);N_{%s};",estimatorName,estimatorName),nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); //
489 
490  fHistGenPrimaryParticlesInelGt0 = new TH1F("hGenPrimaryParticlesInelGt0","Multiplcity of generated charged particles ; Nparticles ; Entries",nMultBins,firstMultBin,lastMultBin);
491 
492  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);
493 
494  fHistNtrUnCorrPSSel->Sumw2();
495  fHistNtrUnCorrPSTrigSel->Sumw2();
501  fHistNtrUnCorrEvSel->Sumw2();
502  fHistNtrUnCorrEvWithCand->Sumw2();
503  fHistNtrUnCorrEvWithD->Sumw2();
504  fHistNtrCorrPSSel->Sumw2();
505  fHistNtrCorrEvSel->Sumw2();
506  fHistNtrCorrEvWithCand->Sumw2();
507  fHistNtrCorrEvWithD->Sumw2();
523  if(fKeepCorrPlots){
551  }
552  fOutput->Add(fHistNtrVsZvtx);
554  fOutput->Add(fHistNtrVsSo);
556 
557  fOutput->Add(fHistNtrVsNchMC);
565 
566  fHistNEvents = new TH1F("fHistNEvents", "number of events ",11,-0.5,10.5);
567  fHistNEvents->GetXaxis()->SetBinLabel(1,"nEvents total");
568  fHistNEvents->GetXaxis()->SetBinLabel(2,"nEvents with Z vertex");
569  fHistNEvents->GetXaxis()->SetBinLabel(3,"nEvents selected");
570  fHistNEvents->GetXaxis()->SetBinLabel(4,"Rejected due to trigger");
571  fHistNEvents->GetXaxis()->SetBinLabel(5,"Rejected due to phys sel");
572  fHistNEvents->GetXaxis()->SetBinLabel(6,"Rejected due to vertex cuts");
573  fHistNEvents->GetXaxis()->SetBinLabel(7,"Rejected due to pileup");
574  fHistNEvents->GetXaxis()->SetBinLabel(8,"Total no. of candidate");
575  fHistNEvents->GetXaxis()->SetBinLabel(9,"no. of cand wo bitmask");
576  fHistNEvents->GetXaxis()->SetBinLabel(10,"D after cuts (No PID)");
577  fHistNEvents->GetXaxis()->SetBinLabel(11,"D after cuts + PID)");
578  fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);
579  fHistNEvents->Sumw2();
580  fHistNEvents->SetMinimum(0);
581  fOutput->Add(fHistNEvents);
582 
583  fPtVsMassVsMultMC=new TH3F("hPtVsMassvsMultMC", "D true candidates: p_{t} vs mass vs tracklets multiplicity; Tracklets; Mass M [GeV/c^{2}]; p_{t} [GeV/c]",nMultBins,firstMultBin,lastMultBin,fNMassBins,fLowmasslimit,fUpmasslimit,48,0.,24.);
584 
585  // With flag fFillSoSparseChecks to fill THnSparse with MultUncorr and NoPid cases ( 0 = only Mult, 1 = Mult and multUncorr, 2 = NoPid and 3 is All)
586  Int_t nbinsSo[4]={48, fNMassBins, 20, nMultBins};
587  Double_t xminSo[4]={0., fLowmasslimit,0., firstMultBin};
588  Double_t xmaxSo[4]={24., fUpmasslimit, 1., lastMultBin};
589 
590  TString HistoName;
591  TString HistoNameNoPid;
592  TString ParName;
593 
595  HistoName = "hSparseSphericity";
596  HistoNameNoPid = "hSparseSphericitywithNoPid";
597  ParName = "Sphericity";
598  }
599  else{
600  HistoName = "hSparseSpherocity";
601  HistoNameNoPid = "hSparseSpherocitywithNoPid";
602  ParName = "Spherocity";
603  }
604  if(fFillSoSparseChecks == 1 || fFillSoSparseChecks == 3){
605  Int_t nbinsSowithMultUncorr[5]={48, fNMassBins, 20, nMultBins, nMultBins};
606  Double_t xminSowithMultUncorr[5]={0., fLowmasslimit,0., firstMultBin, firstMultBin};
607  Double_t xmaxSowithMultUncorr[5]={24., fUpmasslimit, 1., lastMultBin, lastMultBin};
608  fSparseSpherocity = new THnSparseD(HistoName.Data(), Form("D candidates:; p_{T} [GeV/c]; InvMass [GeV/c^{2}]; %s; Multipicity; MultipicityUncorr;", ParName.Data()), 5 , nbinsSowithMultUncorr, xminSowithMultUncorr, xmaxSowithMultUncorr);
609  }
610  else{
611  fSparseSpherocity = new THnSparseD(HistoName.Data(), Form("D candidates:; p_{T} [GeV/c]; InvMass [GeV/c^{2}]; %s; Multipicity;", ParName.Data()), 4 , nbinsSo, xminSo, xmaxSo);
612  }
613 
614  if(fFillSoSparseChecks == 2|| fFillSoSparseChecks == 3) fSparseSpherocitywithNoPid = new THnSparseD(HistoNameNoPid.Data(), Form("D candidates with NoPID:; p_{T} [GeV/c]; InvMass [GeV/c^{2}]; %s; Multipicity;", ParName.Data()), 4 , nbinsSo, xminSo, xmaxSo);
615 
619 
621 
622  fCounterC = new AliNormalizationCounter("NormCounterCorrMult");
623  fCounterC->SetStudyMultiplicity(kTRUE,1.);
624  fCounterC->SetStudySpherocity(kTRUE,10.);
625  fCounterC->Init();
626 
627  fCounterU = new AliNormalizationCounter("NormCounterUnCorrMult");
628  fCounterU->SetStudyMultiplicity(kTRUE,1.);
629  fCounterU->SetStudySpherocity(kTRUE,10.);
630  fCounterU->Init();
631 
632  fCounterCandidates = new AliNormalizationCounter("NormCounterCorrMultCandidates");
636 
637 
638  fOutputCounters = new TList();
639  fOutputCounters->SetOwner();
640  fOutputCounters->SetName("OutputCounters");
644 
645 
646  PostData(1,fOutput);
647  PostData(2,fListCuts);
648  PostData(3,fOutputCounters);
649  PostData(4,fListProfiles);
650 
651  return;
652 }
653 
654 //________________________________________________________________________
656 {
657  // Execute analysis for current event:
658  // heavy flavor candidates association to MC truth
659 
660  AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
661 
662  TClonesArray *arrayCand = 0;
663  TString arrayName="";
664  UInt_t pdgDau[3];
665  Int_t nDau=0;
666  Int_t selbit=0;
667  if(fPdgMeson==411){
668  arrayName="Charm3Prong";
669  pdgDau[0]=211; pdgDau[1]=321; pdgDau[2]=211;
670  nDau=3;
672  }else if(fPdgMeson==421){
673  arrayName="D0toKpi";
674  pdgDau[0]=211; pdgDau[1]=321; pdgDau[2]=0;
675  nDau=2;
677  }else if(fPdgMeson==413){
678  arrayName="Dstar";
679  pdgDau[0]=321; pdgDau[1]=211; pdgDau[2]=0; // Quoting here D0 daughters (D* ones on another variable later)
680  nDau=2;
682  }
683 
684  if(!aod && AODEvent() && IsStandardAOD()) {
685  // In case there is an AOD handler writing a standard AOD, use the AOD
686  // event in memory rather than the input (ESD) event.
687  aod = dynamic_cast<AliAODEvent*> (AODEvent());
688  // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
689  // have to taken from the AOD event hold by the AliAODExtension
690  AliAODHandler* aodHandler = (AliAODHandler*)
691  ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
692  if(aodHandler->GetExtensions()) {
693  AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
694  AliAODEvent *aodFromExt = ext->GetAOD();
695  arrayCand=(TClonesArray*)aodFromExt->GetList()->FindObject(arrayName.Data());
696  }
697  } else if(aod) {
698  arrayCand=(TClonesArray*)aod->GetList()->FindObject(arrayName.Data());
699  }
700 
701  if(!aod || !arrayCand) {
702  printf("AliAnalysisTaskSEDvsEventShapes::UserExec: Charm3Prong branch not found!\n");
703  return;
704  }
705 
706  if(fisPPbData && fReadMC){
707  Int_t runnumber = aod->GetRunNumber();
708  if(aod->GetTriggerMask()==0 &&
709  (runnumber>=195344 && runnumber<=195677)){
710  AliDebug(3,"Event rejected because of null trigger mask");
711  return;
712  }
713  }
714 
715  // fix for temporary bug in ESDfilter
716  // the AODs with null vertex pointer didn't pass the PhysSel
717  if(!aod->GetPrimaryVertex()||TMath::Abs(aod->GetMagneticField())<0.001) return;
718 
719  Int_t countTreta1=0, countTreta03=0, countTreta05=0, countTreta16=0;
720  AliAODTracklets* tracklets=aod->GetTracklets();
721  Int_t nTr=tracklets->GetNumberOfTracklets();
722  for(Int_t iTr=0; iTr<nTr; iTr++){
723  Double_t theta=tracklets->GetTheta(iTr);
724  Double_t eta=-TMath::Log(TMath::Tan(theta/2.));
725  if(eta>-0.3 && eta<0.3) countTreta03++;
726  if(eta>-0.5 && eta<0.5) countTreta05++;
727  if(eta>-1.0 && eta<1.0) countTreta1++;
728  if(eta>-1.6 && eta<1.6) countTreta16++;
729  }
730 
731 
732  Int_t vzeroMult=0, vzeroMultA=0, vzeroMultC=0;
733  Int_t vzeroMultEq=0, vzeroMultAEq=0, vzeroMultCEq=0;
734  AliAODVZERO *vzeroAOD = (AliAODVZERO*)aod->GetVZEROData();
735  if(vzeroAOD) {
736  vzeroMultA = static_cast<Int_t>(vzeroAOD->GetMTotV0A());
737  vzeroMultC = static_cast<Int_t>(vzeroAOD->GetMTotV0C());
738  vzeroMult = vzeroMultA + vzeroMultC;
739  vzeroMultAEq = static_cast<Int_t>(AliVertexingHFUtils::GetVZEROAEqualizedMultiplicity(aod));
740  vzeroMultCEq = static_cast<Int_t>(AliVertexingHFUtils::GetVZEROCEqualizedMultiplicity(aod));
741  vzeroMultEq = vzeroMultAEq + vzeroMultCEq;
742  }
743 
744  Int_t countMult = countTreta1;
745  if(fMultiplicityEstimator==kNtrk03) { countMult = countTreta03; }
746  else if(fMultiplicityEstimator==kNtrk05) { countMult = countTreta05; }
747  else if(fMultiplicityEstimator==kNtrk10to16) { countMult = countTreta16 - countTreta1; }
748  else if(fMultiplicityEstimator==kVZERO) { countMult = vzeroMult; }
749  else if(fMultiplicityEstimator==kVZEROA) { countMult = vzeroMultA; }
750  else if(fMultiplicityEstimator==kVZEROEq) { countMult = vzeroMultEq; }
751  else if(fMultiplicityEstimator==kVZEROAEq) { countMult = vzeroMultAEq; }
752 
753  Double_t spherocity;
754  if(fCalculateSphericity){ //When kTRUE, it calculates Sphericity and THnSparse filled for sphericity
756  }
757  else{
759  }
760  Double_t St=1;
761  // printf("hello \n");
762  fCounterU->StoreEvent(aod,fRDCutsAnalysis,fReadMC,countMult,spherocity);
763  fHistNEvents->Fill(0); // count event
764 
765  Double_t countTreta1corr=countTreta1;
766  Double_t countCorr=countMult;
767  AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
768  // In case of VZERO multiplicity, consider the zvtx correction flag
769  // fDoVZER0ParamVertexCorr: 0= none, 1= usual d2h, 2=AliESDUtils
770  Bool_t isDataDrivenZvtxCorr=kTRUE;
771  Bool_t isVtxOk=kFALSE;
772  Int_t vzeroMultACorr=vzeroMultA, vzeroMultCCorr=vzeroMultC, vzeroMultCorr=vzeroMult;
773  Int_t vzeroMultAEqCorr=vzeroMultAEq, vzeroMultCEqCorr=vzeroMultCEq, vzeroMultEqCorr=vzeroMultEq;
774  if(vtx1){
775  if(vtx1->GetNContributors()>0){
776  fHistNEvents->Fill(1);
777  isVtxOk=kTRUE;
778  }
779  }
780  if(isVtxOk){
784  // do not correct
785  isDataDrivenZvtxCorr=kFALSE;
786  } else if (fDoVZER0ParamVertexCorr==2){
787  // use AliESDUtils correction
788  Float_t zvtx = vtx1->GetZ();
789  isDataDrivenZvtxCorr=kFALSE;
790  vzeroMultACorr = static_cast<Int_t>(AliESDUtils::GetCorrV0A(vzeroMultA,zvtx));
791  vzeroMultCCorr = static_cast<Int_t>(AliESDUtils::GetCorrV0C(vzeroMultC,zvtx));
792  vzeroMultCorr = vzeroMultACorr + vzeroMultCCorr;
793  vzeroMultAEqCorr = static_cast<Int_t>(AliESDUtils::GetCorrV0A(vzeroMultAEq,zvtx));
794  vzeroMultCEqCorr =static_cast<Int_t>( AliESDUtils::GetCorrV0C(vzeroMultCEq,zvtx));
795  vzeroMultEqCorr = vzeroMultAEqCorr + vzeroMultCEqCorr;
796  if(fMultiplicityEstimator==kVZERO) { countCorr = vzeroMultCorr; }
797  else if(fMultiplicityEstimator==kVZEROA) { countCorr = vzeroMultACorr; }
798  else if(fMultiplicityEstimator==kVZEROEq) { countCorr = vzeroMultEqCorr; }
799  else if(fMultiplicityEstimator==kVZEROAEq) { countCorr = vzeroMultAEqCorr; }
800  }
801  }
802  }
803  // Data driven multiplicity z-vertex correction
804  if(isVtxOk && isDataDrivenZvtxCorr){
805  TProfile* estimatorAvg = GetEstimatorHistogram(aod);
806  if(estimatorAvg){
807  countTreta1corr=static_cast<Int_t>(AliVertexingHFUtils::GetCorrectedNtracklets(estimatorAvg,countTreta1,vtx1->GetZ(),fRefMult));
808  countCorr=static_cast<Int_t>(AliVertexingHFUtils::GetCorrectedNtracklets(estimatorAvg,countMult,vtx1->GetZ(),fRefMult));
809  }
810  }
811 
812  fCounterC->StoreEvent(aod,fRDCutsAnalysis,fReadMC,countCorr,spherocity);
813 
814  Bool_t isEvSel=fRDCutsAnalysis->IsEventSelected(aod);
815 
816  if(fRDCutsAnalysis->GetWhyRejection()==5) fHistNEvents->Fill(3);
817  if(fRDCutsAnalysis->GetWhyRejection()==7) fHistNEvents->Fill(4);
818  if(fRDCutsAnalysis->GetWhyRejection()==6) fHistNEvents->Fill(5);
819  if(fRDCutsAnalysis->GetWhyRejection()==1) fHistNEvents->Fill(6);
820 
821  Bool_t isEvPSRejected = fRDCutsAnalysis->IsEventRejectedDuePhysicsSelection();
822  Bool_t isEvTrigNameRejected = fRDCutsAnalysis->IsEventRejectedDueToTrigger();
823  Bool_t isEvPileUpRejected = fRDCutsAnalysis->IsEventRejectedDueToPileup();
824  Bool_t isEvNoVtxRejected = fRDCutsAnalysis->IsEventRejectedDueToNotRecoVertex();
825  Bool_t isEvVtxContribRejected = fRDCutsAnalysis->IsEventRejectedDueToVertexContributors();
827  Bool_t isEvCentralityRejected = fRDCutsAnalysis->IsEventRejectedDueToCentrality();
828  if(!isEvPSRejected){
829  fHistNtrUnCorrPSSel->Fill(countMult);
830  fHistNtrCorrPSSel->Fill(countCorr);
831  if(!isEvTrigNameRejected){
832  fHistNtrUnCorrPSTrigSel->Fill(countMult);
833  if(!isEvPileUpRejected){
834  fHistNtrUnCorrPSTrigPileUpSel->Fill(countMult);
835  if(!isEvNoVtxRejected){
836  fHistNtrUnCorrPSTrigPileUpVtxSel->Fill(countMult);
837  if(!isEvVtxContribRejected){
838  fHistNtrUnCorrPSTrigPileUpVtxContSel->Fill(countMult);
839  if(!isEvVtxRangeRejected){
840  fHistNtrUnCorrPSTrigPileUpVtxRangeSel->Fill(countMult);
841  if(!isEvCentralityRejected){
843  }
844  }
845  }
846  }
847  }
848  }
849  }
850 
851  if(!isEvSel)return;
852  if(fKeepCorrPlots){
853  fHistNtrEta16vsNtrEta1EvSel->Fill(countTreta1,countTreta16);
854  fHistNtrEta05vsNtrEta1EvSel->Fill(countTreta1,countTreta05);
855  fHistNtrEta03vsNtrEta1EvSel->Fill(countTreta1,countTreta03);
856  fHistNtrEtaV0AvsNtrEta1EvSel->Fill(countTreta1,vzeroMultA);
857  fHistNtrEtaV0MvsNtrEta1EvSel->Fill(countTreta1,vzeroMult);
858  fHistNtrEtaV0AvsV0AEqEvSel->Fill(vzeroMultA,vzeroMultAEq);
859  fHistNtrEtaV0MvsV0MEqEvSel->Fill(vzeroMult,vzeroMultEq);
860  fHistNtrCorrEta1vsNtrRawEta1EvSel->Fill(countTreta1,countTreta1corr);
861  fHistMultCorrvsMultRawEvSel->Fill(countMult,countCorr);
862  }
863  if(vtx1){
864  fHistNtrVsZvtx->Fill(vtx1->GetZ(),countMult);
865  fHistNtrCorrVsZvtx->Fill(vtx1->GetZ(),countCorr);
866  fHistNtrVsSo->Fill(spherocity,countMult);
867  fHistNtrCorrVsSo->Fill(spherocity,countCorr);
868  }
869 
870  TClonesArray *arrayMC=0;
871  AliAODMCHeader *mcHeader=0;
872 
873  Double_t nchWeight=1.0;
874 
875  // load MC particles
876  if(fReadMC){
877 
878  arrayMC = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
879  if(!arrayMC) {
880  printf("AliAnalysisTaskSEDvsEventShapes::UserExec: MC particles branch not found!\n");
881  return;
882  }
883  // load MC header
884  mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
885  if(!mcHeader) {
886  printf("AliAnalysisTaskSEDvsEventShapes::UserExec: MC header branch not found!\n");
887  return;
888  }
889 
890 
891  Int_t nChargedMCEta10=0, nChargedMCEta03=0, nChargedMCEta05=0, nChargedMCEta16=0, nChargedMCEtam37tm17=0, nChargedMCEta28t51=0;
892  Int_t nChargedMCPrimaryEta10=0, nChargedMCPrimaryEta03=0, nChargedMCPrimaryEta05=0, nChargedMCPrimaryEta16=0, nChargedMCPrimaryEtam37tm17=0, nChargedMCPrimaryEta28t51=0;
893  Int_t nChargedMCPhysicalPrimaryEta10=0, nChargedMCPhysicalPrimaryEta03=0, nChargedMCPhysicalPrimaryEta05=0, nChargedMCPhysicalPrimaryEta16=0, nChargedMCPhysicalPrimaryEtam37tm17=0, nChargedMCPhysicalPrimaryEta28t51=0;
894  for(Int_t i=0; i<arrayMC->GetEntriesFast(); i++){
895  AliAODMCParticle *part=(AliAODMCParticle*)arrayMC->UncheckedAt(i);
896  Int_t charge = part->Charge();
897  Double_t eta = part->Eta();
898  Bool_t isPrim = part->IsPrimary();
899  Bool_t isPhysPrim = part->IsPhysicalPrimary();
900  if(charge!=0) {
901  if(eta>-0.3 && eta< 0.3) {
902  nChargedMCEta03++;
903  if(isPrim) nChargedMCPrimaryEta03++;
904  if(isPhysPrim) nChargedMCPhysicalPrimaryEta03++;
905  }
906  if(eta>-0.5 && eta< 0.5) {
907  nChargedMCEta05++;
908  if(isPrim) nChargedMCPrimaryEta05++;
909  if(isPhysPrim) nChargedMCPhysicalPrimaryEta05++;
910  }
911  if(eta>-1.0 && eta< 1.0) {
912  nChargedMCEta10++;
913  if(isPrim) nChargedMCPrimaryEta10++;
914  if(isPhysPrim) nChargedMCPhysicalPrimaryEta10++;
915  }
916  if(eta>-1.6 && eta< 1.6) {
917  nChargedMCEta16++;
918  if(isPrim) nChargedMCPrimaryEta16++;
919  if(isPhysPrim) nChargedMCPhysicalPrimaryEta16++;
920  }
921  if(eta>-3.7 && eta<-1.7) {
922  nChargedMCEtam37tm17++;
923  if(isPrim) nChargedMCPrimaryEtam37tm17++;
924  if(isPhysPrim) nChargedMCPhysicalPrimaryEtam37tm17++;
925  }
926  if(eta> 2.8 && eta< 5.1) {
927  nChargedMCEta28t51++;
928  if(isPrim) nChargedMCPrimaryEta28t51++;
929  if(isPhysPrim) nChargedMCPhysicalPrimaryEta28t51++;
930  }
931  }
932  }
933  Int_t nChargedMC=nChargedMCEta10;
934  Int_t nChargedMCPrimary=nChargedMCPrimaryEta10;
935  Int_t nChargedMCPhysicalPrimary=nChargedMCPhysicalPrimaryEta10;
936 
937  // Compute the Nch weights (reference is Ntracklets within |eta|<1.0)
938  if(fUseNchWeight>0){
939 
940  Double_t tmpweight = 1.0;
941  Double_t tmpXweight=nChargedMCPhysicalPrimary; // Nch weights
942  if(fUseNchWeight==2) tmpXweight=countMult; // Ntrk weights
943 
944  if(tmpXweight<=0) tmpweight = 0.0;
945  else{
946  Double_t pMeas = fHistoMeasNch->GetBinContent(fHistoMeasNch->FindBin(tmpXweight));
947  //printf(" pMeas=%2.2f and histo MCNch %s \n",pMeas,fHistoMCNch);
948  Double_t pMC = fHistoMCNch->GetBinContent(fHistoMCNch->FindBin(tmpXweight));
949  tmpweight = pMC>0 ? pMeas/pMC : 0.;
950  }
951  nchWeight *= tmpweight;
952  AliDebug(2,Form("Using Nch weights, Mult=%f Weight=%f\n",tmpXweight,nchWeight));
953 
954  }
955 
956  // Now recompute the variables in case another MC estimator is considered
958  nChargedMC = nChargedMCEta16 - nChargedMCEta10;
959  nChargedMCPrimary = nChargedMCPrimaryEta16 - nChargedMCPrimaryEta10;
960  nChargedMCPhysicalPrimary = nChargedMCPhysicalPrimaryEta16 - nChargedMCPhysicalPrimaryEta10;
961  } else if(fMCPrimariesEstimator==kEta05){
962  nChargedMC = nChargedMCEta05;
963  nChargedMCPrimary = nChargedMCPrimaryEta05;
964  nChargedMCPhysicalPrimary = nChargedMCPhysicalPrimaryEta05;
965  } else if(fMCPrimariesEstimator==kEta03){
966  nChargedMC = nChargedMCEta03;
967  nChargedMCPrimary = nChargedMCPrimaryEta03;
968  nChargedMCPhysicalPrimary = nChargedMCPhysicalPrimaryEta03;
969  } else if(fMCPrimariesEstimator==kEtaVZERO){
970  nChargedMC = nChargedMCEtam37tm17 + nChargedMCEta28t51;
971  nChargedMCPrimary = nChargedMCPrimaryEtam37tm17 + nChargedMCPrimaryEta28t51;
972  nChargedMCPhysicalPrimary = nChargedMCPhysicalPrimaryEtam37tm17 + nChargedMCPhysicalPrimaryEta28t51;
973  } else if(fMCPrimariesEstimator==kEtaVZEROA){
974  nChargedMC = nChargedMCEta28t51;
975  nChargedMCPrimary = nChargedMCPrimaryEta28t51;
976  nChargedMCPhysicalPrimary = nChargedMCPhysicalPrimaryEta28t51;
977  }
978 
979  // Here fill the MC correlation plots
980  if(nChargedMCPhysicalPrimary>0){ // INEL>0 for |eta|<1
981  fHistGenPrimaryParticlesInelGt0->Fill(nChargedMCPhysicalPrimary,nchWeight);
982  }
983 
984  fHistNtrVsNchMC->Fill(nChargedMC,countMult,nchWeight);
985  fHistNtrCorrVsNchMC->Fill(nChargedMC,countCorr,nchWeight);
986 
987  fHistNtrVsNchMCPrimary->Fill(nChargedMCPrimary,countMult,nchWeight);
988  fHistNtrCorrVsNchMCPrimary->Fill(nChargedMCPrimary,countCorr,nchWeight);
989 
990  fHistNtrVsNchMCPhysicalPrimary->Fill(nChargedMCPhysicalPrimary,countMult,nchWeight);
991  fHistNtrCorrVsNchMCPhysicalPrimary->Fill(nChargedMCPhysicalPrimary,countCorr,nchWeight);
992 
993  fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary->Fill(nChargedMC,nChargedMCPrimary,nChargedMCPhysicalPrimary,nchWeight);
994  }
995 
996  Int_t nCand = arrayCand->GetEntriesFast();
997  Int_t nSelectedNoPID=0,nSelectedPID=0,nSelectedInMassPeak=0;
998  Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();
999  Double_t mDplusPDG = TDatabasePDG::Instance()->GetParticle(411)->Mass();
1000  Double_t mDstarPDG = TDatabasePDG::Instance()->GetParticle(413)->Mass();
1001 
1002  // pdg of daughters needed for D* too
1003  UInt_t pdgDgDStartoD0pi[2]={421,211};
1004 
1005  Double_t aveMult=0.;
1006  Double_t nSelCand=0.;
1007  for (Int_t iCand = 0; iCand < nCand; iCand++) {
1008  AliAODRecoDecayHF *d = (AliAODRecoDecayHF*)arrayCand->UncheckedAt(iCand);
1009  AliAODRecoCascadeHF *dCascade = NULL;
1010  if(fPdgMeson==413) dCascade = (AliAODRecoCascadeHF*)d;
1011 
1012  fHistNEvents->Fill(7);
1013  if(fUseBit && !d->HasSelectionBit(selbit)){
1014  fHistNEvents->Fill(8);
1015  continue;
1016  }
1017 
1018  Double_t ptCand = d->Pt();
1019  Double_t rapid=d->Y(fPdgMeson);
1020  Bool_t isFidAcc=fRDCutsAnalysis->IsInFiducialAcceptance(ptCand,rapid);
1021  if(!isFidAcc) continue;
1022 
1023  Int_t labD=-1;
1024  if(fReadMC) {
1025  if(fPdgMeson==413){
1026  labD = dCascade->MatchToMC(fPdgMeson,421,(Int_t*)pdgDgDStartoD0pi,(Int_t*)pdgDau,arrayMC);
1027  } else {
1028  labD = d->MatchToMC(fPdgMeson,arrayMC,nDau,(Int_t*)pdgDau);
1029  }
1030  FillMCMassHistos(arrayMC,labD, countMult,nchWeight);
1031  }
1032 
1033  Int_t passAllCuts=fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kAll,aod);
1034  Int_t passTopolCuts=fRDCutsAnalysis->GetIsSelectedCuts();
1035  if(passTopolCuts==0) continue;
1036  nSelectedNoPID++;
1037  fHistNEvents->Fill(9);
1038  if(passAllCuts){
1039  nSelectedPID++;
1040  fHistNEvents->Fill(10);
1041  }
1042  Double_t multForCand = countCorr;
1043 
1045  // For the D* case, subtract only the D0 daughter tracks <=== FIXME !!
1046  AliAODRecoDecayHF2Prong* d0fromDstar = NULL;
1047  if(fPdgMeson==413) d0fromDstar = (AliAODRecoDecayHF2Prong*)dCascade->Get2Prong();
1048 
1049  for(Int_t iDau=0; iDau<nDau; iDau++){
1050  AliAODTrack *t = NULL;
1051  if(fPdgMeson==413){ t = (AliAODTrack*)d0fromDstar->GetDaughter(iDau); }
1052  else{ t = (AliAODTrack*)d->GetDaughter(iDau); }
1053  if(!t) continue;
1054  if(t->HasPointOnITSLayer(0) && t->HasPointOnITSLayer(1)){
1055  if(multForCand>0) multForCand-=1;
1056  }
1057  }
1058  }
1059  Bool_t isPrimary=kTRUE;
1060  Double_t trueImpParXY=9999.;
1061  Double_t impparXY=d->ImpParXY()*10000.;
1062  Double_t dlen=0.1; //FIXME
1063  Double_t mass[2];
1064  if(fPdgMeson==411){
1065  mass[0]=d->InvMass(nDau,pdgDau);
1066  mass[1]=-1.;
1067  if(TMath::Abs(mass[0]-mDplusPDG)<0.02) nSelectedInMassPeak++; //20 MeV for now... FIXME
1068  }else if(fPdgMeson==421){
1069  UInt_t pdgdaughtersD0[2]={211,321};//pi,K
1070  UInt_t pdgdaughtersD0bar[2]={321,211};//K,pi
1071  mass[0]=d->InvMass(2,pdgdaughtersD0);
1072  mass[1]=d->InvMass(2,pdgdaughtersD0bar);
1073  if(TMath::Abs(mass[0]-mD0PDG)<0.02 || TMath::Abs(mass[1]-mD0PDG)<0.02 ) nSelectedInMassPeak++; //20 MeV for now... FIXME
1074  }else if(fPdgMeson==413){
1075  // FIXME
1076  mass[0]=dCascade->DeltaInvMass();
1077  mass[1]=-1.;
1078  if(TMath::Abs(mass[0]-(mDstarPDG-mD0PDG))<0.0015) nSelectedInMassPeak++; //1 MeV for now... FIXME
1079  }
1080  for(Int_t iHyp=0; iHyp<2; iHyp++){
1081  if(mass[iHyp]<0.) continue; // for D+ and D* we have 1 mass hypothesis
1082  Double_t invMass=mass[iHyp];
1083  Double_t arrayForSparse[5]={invMass,ptCand,impparXY,dlen,multForCand};
1084 
1085  if(fReadMC){
1086 
1087  if(fPdgMeson==413){
1088  labD = dCascade->MatchToMC(fPdgMeson,421,(Int_t*)pdgDgDStartoD0pi,(Int_t*)pdgDau,arrayMC);
1089  } else {
1090  labD = d->MatchToMC(fPdgMeson,arrayMC,nDau,(Int_t*)pdgDau);
1091  }
1092 
1093  Bool_t fillHisto=fDoImpPar;
1094  if(labD>=0){
1095  AliAODMCParticle *partD = (AliAODMCParticle*)arrayMC->At(labD);
1096  Int_t code=partD->GetPdgCode();
1097  if(CheckOrigin(arrayMC,partD)==5) isPrimary=kFALSE;
1098  if(code<0 && iHyp==0) fillHisto=kFALSE;
1099  if(code>0 && iHyp==1) fillHisto=kFALSE;
1100  if(!isPrimary){
1101  if(fPdgMeson==411){
1102  trueImpParXY=AliVertexingHFUtils::GetTrueImpactParameterDplus(mcHeader,arrayMC,partD)*10000.;
1103  }else if(fPdgMeson==421){
1104  trueImpParXY=AliVertexingHFUtils::GetTrueImpactParameterDzero(mcHeader,arrayMC,partD)*10000.;
1105  }else if(fPdgMeson==413){
1106  trueImpParXY=0.;
1107  }
1108  Double_t arrayForSparseTrue[5]={invMass,ptCand,trueImpParXY,dlen,multForCand};
1109  if(fillHisto && passAllCuts){
1110  fHistMassPtImpPar[2]->Fill(arrayForSparse);
1111  fHistMassPtImpPar[3]->Fill(arrayForSparseTrue);
1112  }
1113  }else{
1114  if(fillHisto && passAllCuts) fHistMassPtImpPar[1]->Fill(arrayForSparse);
1115  }
1116  }else{
1117  if(fillHisto && passAllCuts)fHistMassPtImpPar[4]->Fill(arrayForSparse);
1118  }
1119 
1120  if(TMath::Abs(labD)==fPdgMeson && fMCOption==2) continue;
1121  if(TMath::Abs(labD)!=fPdgMeson && fMCOption==1) continue;
1122 
1123  }
1124 
1125  if(fPdgMeson==421){
1126  if(iHyp==0 && !(passTopolCuts&1)) continue; // candidate not passing as D0
1127  if(iHyp==1 && !(passTopolCuts&2)) continue; // candidate not passing as D0bar
1128  }
1129 
1130  if(fFillSoSparseChecks == 2 || fFillSoSparseChecks == 3){ //Filling THnSparse for Spherocity without PID
1131  fSparseSpherocitywithNoPid->Sumw2();
1132  Double_t arrayForSparseSoNoPid[4]={ptCand, invMass, spherocity, multForCand};
1133  fSparseSpherocitywithNoPid->Fill(arrayForSparseSoNoPid);
1134  }
1135  if(fPdgMeson==421){
1136  if(iHyp==0 && !(passAllCuts&1)) continue; // candidate not passing as D0
1137  if(iHyp==1 && !(passAllCuts&2)) continue; // candidate not passing as D0bar
1138  }
1139  if(passAllCuts){
1140  aveMult+=multForCand;
1141  nSelCand+=1.;
1142 
1143  if(fFillSoSparseChecks == 1 || fFillSoSparseChecks == 3){
1144  fSparseSpherocity->Sumw2();
1145  Double_t arrayForSparseSowithMultUnncorr[5]={ptCand, invMass, spherocity, multForCand, (Double_t)countTreta1};
1146  fSparseSpherocity->Fill(arrayForSparseSowithMultUnncorr);
1147  }
1148  else{
1149  fSparseSpherocity->Sumw2();
1150  Double_t arrayForSparseSo[4]={ptCand, invMass, spherocity, multForCand};
1151  fSparseSpherocity->Fill(arrayForSparseSo);
1152  }
1153  if(fDoImpPar){
1154  fHistMassPtImpPar[0]->Fill(arrayForSparse);
1155  }
1156 
1157  }
1158 
1159  }
1160  }
1161  if(fSubtractTrackletsFromDau && nSelCand>0){
1162  aveMult/=nSelCand;
1163  fCounterCandidates->StoreEvent(aod,fRDCutsAnalysis,fReadMC,(Int_t)(aveMult+0.5001),spherocity);
1164  }else{
1165  fCounterCandidates->StoreEvent(aod,fRDCutsAnalysis,fReadMC,(Int_t)countCorr,spherocity);
1166  }
1167 
1168 
1169  fCounterCandidates->StoreCandidates(aod,nSelectedNoPID,kTRUE);
1170  fCounterCandidates->StoreCandidates(aod,nSelectedPID,kFALSE);
1171  fHistNtrUnCorrEvSel->Fill(countMult,nchWeight);
1172  fHistNtrCorrEvSel->Fill(countCorr,nchWeight);
1173  if(nSelectedPID>0) {
1174  fHistNtrUnCorrEvWithCand->Fill(countMult,nchWeight);
1175  fHistNtrCorrEvWithCand->Fill(countCorr,nchWeight);
1176  if(fKeepCorrPlots){
1177  fHistNtrEta16vsNtrEta1EvWithCand->Fill(countTreta1,countTreta16);
1178  fHistNtrEta05vsNtrEta1EvWithCand->Fill(countTreta1,countTreta05);
1179  fHistNtrEta03vsNtrEta1EvWithCand->Fill(countTreta1,countTreta03);
1180  fHistNtrEtaV0AvsNtrEta1EvWithCand->Fill(countTreta1,vzeroMultA);
1181  fHistNtrEtaV0MvsNtrEta1EvWithCand->Fill(countTreta1,vzeroMult);
1182  fHistNtrEtaV0AvsV0AEqEvWithCand->Fill(vzeroMultA,vzeroMultAEq);
1183  fHistNtrEtaV0MvsV0MEqEvWithCand->Fill(vzeroMult,vzeroMultEq);
1184  fHistNtrCorrEta1vsNtrRawEta1EvWithCand->Fill(countTreta1,countTreta1corr);
1185  fHistMultCorrvsMultRawEvWithCand->Fill(countMult,countCorr);
1186  }
1187  }
1188  if(nSelectedInMassPeak>0) {
1189  fHistNtrUnCorrEvWithD->Fill(countMult,nchWeight);
1190  fHistNtrCorrEvWithD->Fill(countCorr,nchWeight);
1191  if(fKeepCorrPlots){
1192  fHistNtrEta16vsNtrEta1EvWithD->Fill(countTreta1,countTreta16);
1193  fHistNtrEta05vsNtrEta1EvWithD->Fill(countTreta1,countTreta05);
1194  fHistNtrEta03vsNtrEta1EvWithD->Fill(countTreta1,countTreta03);
1195  fHistNtrEtaV0AvsNtrEta1EvWithD->Fill(countTreta1,vzeroMultA);
1196  fHistNtrEtaV0MvsNtrEta1EvWithD->Fill(countTreta1,vzeroMult);
1197  fHistNtrEtaV0AvsV0AEqEvWithD->Fill(vzeroMultA,vzeroMultAEq);
1198  fHistNtrEtaV0MvsV0MEqEvWithD->Fill(vzeroMult,vzeroMultEq);
1199  fHistNtrCorrEta1vsNtrRawEta1EvWithD->Fill(countTreta1,countTreta1corr);
1200  fHistMultCorrvsMultRawEvWithD->Fill(countMult,countCorr);
1201  }
1202  }
1203 
1204  PostData(1,fOutput);
1205  PostData(2,fListCuts);
1206  PostData(3,fOutputCounters);
1207 
1208  return;
1209 }
1210 
1211 //________________________________________________________________________
1213  // Histos for impact paramter study
1214  // mass . pt , impact parameter , decay length , multiplicity
1215 
1216  Int_t nbins[5]={fNMassBins,200,fNImpParBins,50,100};
1217  Double_t xmin[5]={fLowmasslimit,0.,fLowerImpPar,0.,0.};
1218  Double_t xmax[5]={fUpmasslimit,20.,fHigherImpPar,1.,100.};
1219 
1220  fHistMassPtImpPar[0]=new THnSparseF("hMassPtImpParAll",
1221  "Mass vs. pt vs.imppar - All",
1222  5,nbins,xmin,xmax);
1223  fHistMassPtImpPar[1]=new THnSparseF("hMassPtImpParPrompt",
1224  "Mass vs. pt vs.imppar - promptD",
1225  5,nbins,xmin,xmax);
1226  fHistMassPtImpPar[2]=new THnSparseF("hMassPtImpParBfeed",
1227  "Mass vs. pt vs.imppar - DfromB",
1228  5,nbins,xmin,xmax);
1229  fHistMassPtImpPar[3]=new THnSparseF("hMassPtImpParTrueBfeed",
1230  "Mass vs. pt vs.true imppar -DfromB",
1231  5,nbins,xmin,xmax);
1232  fHistMassPtImpPar[4]=new THnSparseF("hMassPtImpParBkg",
1233  "Mass vs. pt vs.imppar - backgr.",
1234  5,nbins,xmin,xmax);
1235  for(Int_t i=0; i<5;i++){
1236  fOutput->Add(fHistMassPtImpPar[i]);
1237  }
1238 }
1239 
1240 //________________________________________________________________________
1242 {
1243  // Terminate analysis
1244  //
1245  if(fDebug > 1) printf("AnalysisTaskSEDvsMultiplicity: Terminate() \n");
1246 
1247  fOutput = dynamic_cast<TList*> (GetOutputData(1));
1248  if (!fOutput) {
1249  printf("ERROR: fOutput not available\n");
1250  return;
1251  }
1252 
1253  fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
1254  if(!fHistNEvents){
1255  printf("ERROR: fHistNEvents not available\n");
1256  return;
1257  }
1258  printf("Number of Analyzed Events = %d\n",(Int_t)fHistNEvents->GetBinContent(3));
1259 
1260  return;
1261 }
1262 //_________________________________________________________________________________________________
1263 Int_t AliAnalysisTaskSEDvsEventShapes::CheckOrigin(TClonesArray* arrayMC, AliAODMCParticle *mcPartCandidate) const {
1264  //
1265  // checking whether the mother of the particles come from a charm or a bottom quark
1266  //
1267 
1268  Int_t pdgGranma = 0;
1269  Int_t mother = 0;
1270  mother = mcPartCandidate->GetMother();
1271  Int_t istep = 0;
1272  Int_t abspdgGranma =0;
1273  Bool_t isFromB=kFALSE;
1274  // Bool_t isQuarkFound=kFALSE;
1275  while (mother >0 ){
1276  istep++;
1277  AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
1278  if (mcGranma){
1279  pdgGranma = mcGranma->GetPdgCode();
1280  abspdgGranma = TMath::Abs(pdgGranma);
1281  if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
1282  isFromB=kTRUE;
1283  }
1284  // if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
1285  mother = mcGranma->GetMother();
1286  }else{
1287  AliError("Failed casting the mother particle!");
1288  break;
1289  }
1290  }
1291 
1292  if(isFromB) return 5;
1293  else return 4;
1294 }
1295 
1296 
1297 
1298 //____________________________________________________________________________
1300  // Get Estimator Histogram from period event->GetRunNumber();
1301  //
1302  // If you select SPD tracklets in |eta|<1 you should use type == 1
1303  //
1304 
1305  Int_t runNo = event->GetRunNumber();
1306  Int_t period = -1; // pp: 0-LHC10b, 1-LHC10c, 2-LHC10d, 3-LHC10e
1307  // pPb: 0-LHC13b, 1-LHC13c
1308  if (fisPPbData) {
1309  if (runNo>195343 && runNo<195484) period = 0;
1310  if (runNo>195528 && runNo<195678) period = 1;
1311  if (period < 0 || period > 1) return 0;
1312  }
1313  else {
1314  if(runNo>114930 && runNo<117223) period = 0;
1315  if(runNo>119158 && runNo<120830) period = 1;
1316  if(runNo>122373 && runNo<126438) period = 2;
1317  if(runNo>127711 && runNo<130841) period = 3;
1318  if(period<0 || period>3) return 0;
1319 
1320  }
1321 
1322  return fMultEstimatorAvg[period];
1323 }
1324 
1325 //__________________________________________________________________________________________________
1327  // creates historgam with measured multiplcity distribution in pp 7 TeV collisions (from Eur. Phys. J. C (2010) 68: 345–354)
1328  //
1329  // for Nch > 70 the points were obtainedwith a double NBD distribution
1330  // 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
1331  // fit1->SetParameter(1,1.63); // k parameter
1332  // fit1->SetParameter(2,12.8); // mean multiplicity
1333  Double_t nchbins[82]={0.50,1.50,2.50,3.50,4.50,5.50,6.50,7.50,8.50,9.50,
1334  10.50,11.50,12.50,13.50,14.50,15.50,16.50,17.50,18.50,19.50,
1335  20.50,21.50,22.50,23.50,24.50,25.50,26.50,27.50,28.50,29.50,
1336  30.50,31.50,32.50,33.50,34.50,35.50,36.50,37.50,38.50,39.50,
1337  40.50,41.50,42.50,43.50,44.50,45.50,46.50,47.50,48.50,49.50,
1338  50.50,51.50,52.50,53.50,54.50,55.50,56.50,57.50,58.50,59.50,
1339  60.50,62.50,64.50,66.50,68.50,70.50,72.50,74.50,76.50,78.50,
1340  80.50,82.50,84.50,86.50,88.50,90.50,92.50,94.50,96.50,98.50,
1341  100.50,102.50};
1342  Double_t pch[81]={0.062011,0.072943,0.070771,0.067245,0.062834,0.057383,0.051499,0.04591,0.041109,0.036954,
1343  0.03359,0.030729,0.028539,0.026575,0.024653,0.0229,0.021325,0.019768,0.018561,0.017187,
1344  0.01604,0.014836,0.013726,0.012576,0.011481,0.010393,0.009502,0.008776,0.008024,0.007452,
1345  0.006851,0.006428,0.00594,0.005515,0.005102,0.00469,0.004162,0.003811,0.003389,0.003071,
1346  0.002708,0.002422,0.002184,0.001968,0.00186,0.00165,0.001577,0.001387,0.001254,0.001118,
1347  0.001037,0.000942,0.000823,0.000736,0.000654,0.000579,0.000512,0.00049,0.00045,0.000355,
1348  0.000296,0.000265,0.000193,0.00016,0.000126,0.0000851, 0.0000676,0.0000537,0.0000426, 0.0000338,
1349  0.0000268,0.0000213,0.0000166,0.0000133,0.0000106,0.00000837,0.00000662, 0.00000524,0.00000414, 0.00000327,
1350  0.00000258};
1351 
1352  if(fHistoMeasNch) delete fHistoMeasNch;
1353  fHistoMeasNch=new TH1F("hMeaseNch","",81,nchbins);
1354  for(Int_t i=0; i<81; i++){
1355  fHistoMeasNch->SetBinContent(i+1,pch[i]);
1356  fHistoMeasNch->SetBinError(i+1,0.);
1357  }
1358 }
1359 
1360 //__________________________________________________________________________________________________
1361 void AliAnalysisTaskSEDvsEventShapes::FillMCMassHistos(TClonesArray *arrayMC, Int_t labD, Int_t countMult,Double_t nchWeight)
1362 {
1363  //
1364  // Function to fill the true MC signal
1365  //
1366 
1367  if(labD>=0){
1368  AliAODMCParticle *partD = (AliAODMCParticle*)arrayMC->At(labD);
1369  Double_t mass = partD->M();
1370  Double_t pt = partD->Pt();
1371  fPtVsMassVsMultMC->Fill(countMult,mass,pt,nchWeight);
1372  }
1373 
1374 }
Int_t charge
Bool_t IsEventRejectedDueToCentrality() const
Definition: AliRDHFCuts.h:310
Int_t pdg
TH2F * fHistMultCorrvsMultRawEvWithD
hist. for Ntracklets in eta<1 with and w/o corrections for events with a candidate in D mass peak ...
Bool_t IsEventRejectedDueToZVertexOutsideFiducialRegion() const
Definition: AliRDHFCuts.h:304
TH2F * fHistNtrEta05vsNtrEta1EvWithD
hist. for Ntracklets in eta<1.6 vs. eta<1. for events with a candidate in D mass peak ...
TH2F * fHistNtrEta05vsNtrEta1EvSel
hist. for Ntracklets in eta<1.6 vs. eta<1.
TH2F * fHistNtrCorrVsSo
hist of ntracklets vs So
TH2F * fHistNtrEtaV0AvsNtrEta1EvWithCand
hist. for Ntracklets in eta<0.3 vs. eta<1. for events with a candidate
Bool_t IsEventRejectedDueToNotRecoVertex() const
Definition: AliRDHFCuts.h:298
TH1F * fHistNtrUnCorrPSTrigPileUpVtxContSel
hist. of ntracklets for physics selection + trigger name + pileup + with-vertex selected events ...
void StoreCandidates(AliVEvent *, Int_t nCand=0, Bool_t flagFilter=kTRUE)
TH2F * fHistNtrEta03vsNtrEta1EvWithD
hist. for Ntracklets in eta<0.5 vs. eta<1. for events with a candidate in D mass peak ...
Double_t DeltaInvMass() const
TH2F * fHistNtrCorrEta1vsNtrRawEta1EvWithCand
hist. for V0M raw mult vs V0M equalized multiplicity for events with a candidate
AliNormalizationCounter * fCounterCandidates
Counter for normalization, uncorrected multiplicity.
TH2F * fHistNtrCorrEta1vsNtrRawEta1EvWithD
hist. for V0M raw mult vs V0M equalized multiplicity with a candidate in D mass peak ...
Int_t GetIsSelectedCuts() const
Definition: AliRDHFCuts.h:333
TH1F * fHistNtrUnCorrPSSel
hist of Nch (generated) vs Nch (Primary) vs Nch (Physical Primary)
TH2F * fHistNtrEta03vsNtrEta1EvSel
hist. for Ntracklets in eta<0.5 vs. eta<1.
Bool_t HasSelectionBit(Int_t i) const
TH2F * fHistNtrEta16vsNtrEta1EvWithCand
hist. for multiplicity with and w/o corrections
Int_t MatchToMC(Int_t pdgabs, Int_t pdgabs2prong, Int_t *pdgDg, Int_t *pdgDg2prong, TClonesArray *mcArray, Bool_t isV0=kFALSE) const
TH1F * fHistNtrUnCorrPSTrigPileUpVtxRangeSel
hist. of ntracklets for physics selection + trigger name + pileup + with-vertex-contrib selected even...
TH2F * fHistNtrEtaV0AvsV0AEqEvWithCand
hist. for Ntracklets in eta-V0M vs. eta<1. for events with a candidate
Double_t mass
TH2F * fHistMultCorrvsMultRawEvWithCand
hist. for Ntracklets in eta<1 with and w/o corrections for events with a candidate ...
TH2F * fHistNtrEtaV0AvsV0AEqEvWithD
hist. for Ntracklets in eta-V0M vs. eta<1. for events with a candidate in D mass peak ...
THnSparseD * fSparseSpherocity
hist. of Pt vs Mult vs. mass (MC true candidates before reconstruction)
Double_t ImpParXY() const
TH2F * fHistNtrVsNchMCPhysicalPrimary
hist of ntracklets vs Nch (Primary)
TH2F * fHistNtrCorrVsNchMCPhysicalPrimary
hist of ntracklets vs Nch (Physical Primary)
Bool_t IsEventRejectedDueToVertexContributors() const
Definition: AliRDHFCuts.h:301
TH3F * fPtVsMassVsMultMC
hist. of ntracklets for evnts with a candidate in D mass peak
Int_t CheckOrigin(TClonesArray *arrayMC, AliAODMCParticle *mcPartCandidate) const
Int_t GetWhyRejection() const
Definition: AliRDHFCuts.h:293
TH2F * fHistNtrVsNchMC
hist of ntracklets vs So
TH1F * fHistNtrCorrPSSel
hist. of ntracklets for evnts with a candidate in D mass peak
TH2F * fHistNtrEtaV0AvsNtrEta1EvWithD
hist. for Ntracklets in eta<0.3 vs. eta<1. for events with a candidate in D mass peak ...
TH2F * fHistNtrCorrVsZvtx
hist of ntracklets vs Zvertex
TH2F * fHistNtrEta16vsNtrEta1EvWithD
hist. for multiplicity with and w/o corrections for events with a candidate
THnSparseF * fHistMassPtImpPar[5]
THnSparse histograms for Spherocity studies.
TH2F * fHistNtrVsZvtx
hist. for multiplicity with and w/o corrections for events with a candidate in D mass peak ...
TH1F * fHistNtrUnCorrPSTrigPileUpVtxRangeCentrSel
hist. of ntracklets for physics selection + trigger name + pileup + with-vertex-contrib-range selecte...
Bool_t fDoImpPar
Counter for normalization, corrected multiplicity for candidates.
TList * fListProfiles
list send on output slot 3
static Double_t GetVZEROCEqualizedMultiplicity(AliAODEvent *ev)
Class for cuts on AOD reconstructed D+->Kpipi.
void SetStudyMultiplicity(Bool_t flag, Float_t etaRange)
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)
Functions for event shape variables.
TH2F * fHistNtrCorrVsNchMC
hist of ntracklets vs Nch (Generated)
TH2F * fHistNtrEtaV0MvsNtrEta1EvWithD
hist. for Ntracklets in eta-V0A vs. eta<1. for events with a candidate in D mass peak ...
TH2F * fHistNtrEtaV0MvsV0MEqEvSel
hist. for V0A raw mult vs V0A equalized multiplicity
THnSparseD * fSparseSpherocitywithNoPid
THnSparse histograms for Spherocity studies.
TH2F * fHistNtrEtaV0MvsV0MEqEvWithD
hist. for V0A raw mult vs V0A equalized multiplicity with a candidate in D mass peak ...
TH2F * fHistNtrVsNchMCPrimary
hist of ntracklets vs Nch (Generated)
TH3F * fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary
hist. of geenrated multiplcity
AliNormalizationCounter * fCounterU
Counter for normalization, corrected multiplicity.
TH2F * fHistNtrEtaV0AvsV0AEqEvSel
hist. for Ntracklets in eta-V0M vs. eta<1.
void SetStudySpherocity(Bool_t flag, Double_t nsteps=100.)
TH2F * fHistNtrEta05vsNtrEta1EvWithCand
hist. for Ntracklets in eta<1.6 vs. eta<1. for events with a candidate
Bool_t IsEventRejectedDueToPileup() const
Definition: AliRDHFCuts.h:307
TH1F * fHistNtrCorrEvSel
hist. of ntracklets for physics selection only selected events
void SetMassLimits(Double_t lowlimit, Double_t uplimit)
TH2F * fHistNtrEtaV0MvsV0MEqEvWithCand
hist. for V0A raw mult vs V0A equalized multiplicity for events with a candidate
TH2F * fHistNtrEta16vsNtrEta1EvSel
hist. for No. of events
TH2F * fHistNtrVsSo
hist of ntracklets vs Zvertex
TH1F * fHistNtrUnCorrEvWithD
hist. of ntracklets for evnts with a candidate
TH2F * fHistNtrEta03vsNtrEta1EvWithCand
hist. for Ntracklets in eta<0.5 vs. eta<1. for events with a candidate
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)
TH1F * fHistNtrCorrEvWithD
hist. of ntracklets for evnts with a candidate
Bool_t IsEventRejectedDuePhysicsSelection() const
Definition: AliRDHFCuts.h:319
void FillMCMassHistos(TClonesArray *arrayMC, Int_t labD, Int_t countMult, Double_t nchWeight)
TProfile * GetEstimatorHistogram(const AliVEvent *event)
TH1F * fHistNtrCorrEvWithCand
hist. of ntracklets for selected events
ClassImp(AliAnalysisTaskSEDvsEventShapes) AliAnalysisTaskSEDvsEventShapes
static Double_t GetCorrectedNtracklets(TProfile *estimatorAvg, Double_t uncorrectedNacc, Double_t vtxZ, Double_t refMult)
Bool_t IsEventSelected(AliVEvent *event)
void StoreEvent(AliVEvent *, AliRDHFCuts *, Bool_t mc=kFALSE, Int_t multiplicity=-9999, Double_t spherocity=-99.)
Double_t fUpmasslimit
histograms for impact paramter studies
Bool_t IsSelected(TObject *obj)
Definition: AliRDHFCuts.h:273
TH1F * fHistNtrUnCorrPSTrigSel
hist. of ntracklets for physics selection only selected events
TH1F * fHistNtrUnCorrEvSel
hist. of ntracklets for physics selection + trigger name + pileup + with-vertex-contrib-range + centr...
TList * fListCuts
list send on output slot 1
TH1F * fHistGenPrimaryParticlesInelGt0
hist of ntracklets vs Nch (Physical Primary)
static Double_t GetVZEROAEqualizedMultiplicity(AliAODEvent *ev)
Utilities for V0 multiplicity checks.
Bool_t IsEventRejectedDueToTrigger() const
Definition: AliRDHFCuts.h:295
TH1F * fHistNtrUnCorrPSTrigPileUpSel
hist. of ntracklets for physics selection + trigger name selected events
static Double_t GetTrueImpactParameterDzero(AliAODMCHeader *mcHeader, TClonesArray *arrayMC, AliAODMCParticle *partDp)
Functions for computing true impact parameter of D meson.
const Int_t nbins
TH2F * fHistMultCorrvsMultRawEvSel
hist. for Ntracklets in eta<1 with and w/o corrections
TH2F * fHistNtrCorrVsNchMCPrimary
hist of ntracklets vs Nch (Primary)
TH1F * fHistNtrUnCorrPSTrigPileUpVtxSel
hist. of ntracklets for physics selection + trigger name + pileup selected events ...
AliAODRecoDecayHF2Prong * Get2Prong() const
static Double_t GetTrueImpactParameterDplus(AliAODMCHeader *mcHeader, TClonesArray *arrayMC, AliAODMCParticle *partDp)
TH2F * fHistNtrEtaV0MvsNtrEta1EvSel
hist. for Ntracklets in eta-V0A vs. eta<1.
virtual Bool_t IsInFiducialAcceptance(Double_t, Double_t) const
Definition: AliRDHFCuts.h:290
TH2F * fHistNtrEtaV0MvsNtrEta1EvWithCand
hist. for Ntracklets in eta-V0A vs. eta<1. for events with a candidate
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)
TH2F * fHistNtrEtaV0AvsNtrEta1EvSel
hist. for Ntracklets in eta<0.3 vs. eta<1.
TH1F * fHistNtrUnCorrEvWithCand
hist. of ntracklets for selected events
TH2F * fHistNtrCorrEta1vsNtrRawEta1EvSel
hist. for V0M raw mult vs V0M equalized multiplicity