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