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