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