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