1 /**************************************************************************
2 * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 //*************************************************************************
19 // Class AliAnalysisTaskSEDvsMultiplicity
20 // AliAnalysisTaskSE for the D meson vs. multiplcity analysis
21 // Authors: Renu Bala, Zaida Conesa del Valle, Francesco Prino
22 /////////////////////////////////////////////////////////////
24 #include <TClonesArray.h>
28 #include <TDatabasePDG.h>
32 #include <THnSparse.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 ClassImp(AliAnalysisTaskSEDvsMultiplicity)
54 //________________________________________________________________________
55 AliAnalysisTaskSEDvsMultiplicity::AliAnalysisTaskSEDvsMultiplicity():
62 fHistNtrEta16vsNtrEta1EvSel(0),
63 fHistNtrEta05vsNtrEta1EvSel(0),
64 fHistNtrEta03vsNtrEta1EvSel(0),
65 fHistNtrEtaV0AvsNtrEta1EvSel(0),
66 fHistNtrEtaV0MvsNtrEta1EvSel(0),
67 fHistNtrCorrEta1vsNtrRawEta1EvSel(0),
68 fHistNtrEta16vsNtrEta1EvWithCand(0),
69 fHistNtrEta05vsNtrEta1EvWithCand(0),
70 fHistNtrEta03vsNtrEta1EvWithCand(0),
71 fHistNtrEtaV0AvsNtrEta1EvWithCand(0),
72 fHistNtrEtaV0MvsNtrEta1EvWithCand(0),
73 fHistNtrCorrEta1vsNtrRawEta1EvWithCand(0),
74 fHistNtrEta16vsNtrEta1EvWithD(0),
75 fHistNtrEta05vsNtrEta1EvWithD(0),
76 fHistNtrEta03vsNtrEta1EvWithD(0),
77 fHistNtrEtaV0AvsNtrEta1EvWithD(0),
78 fHistNtrEtaV0MvsNtrEta1EvWithD(0),
79 fHistNtrCorrEta1vsNtrRawEta1EvWithD(0),
81 fHistNtrCorrVsZvtx(0),
83 fHistNtrCorrVsNchMC(0),
84 fHistNtrVsNchMCPrimary(0),
85 fHistNtrCorrVsNchMCPrimary(0),
86 fHistNtrVsNchMCPhysicalPrimary(0),
87 fHistNtrCorrVsNchMCPhysicalPrimary(0),
88 fHistGenPrimaryParticlesInelGt0(0),
89 fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary(0),
90 fHistNtrUnCorrPSSel(0),
91 fHistNtrUnCorrEvSel(0),
92 fHistNtrUnCorrEvWithCand(0),
93 fHistNtrUnCorrEvWithD(0),
96 fHistNtrCorrEvWithCand(0),
97 fHistNtrCorrEvWithD(0),
99 fPtVsMassVsMultNoPid(0),
100 fPtVsMassVsMultUncorr(0),
101 fPtVsMassVsMultPart(0),
102 fPtVsMassVsMultAntiPart(0),
103 fPtVsMassVsMultMC(0),
105 fLowmasslimit(1.765),
112 fLowerImpPar(-2000.),
113 fHigherImpPar(2000.),
118 fSubtractTrackletsFromDau(kFALSE),
119 fKeepCorrPlots(kFALSE),
120 fUseNchWeight(kFALSE),
125 fMultiplicityEstimator(kNtrk10),
126 fMCPrimariesEstimator(kEta10)
128 // Default constructor
129 for(Int_t i=0; i<5; i++) fHistMassPtImpPar[i]=0;
130 for(Int_t i=0; i<4; i++) fMultEstimatorAvg[i]=0;
133 //________________________________________________________________________
134 AliAnalysisTaskSEDvsMultiplicity::AliAnalysisTaskSEDvsMultiplicity(const char *name, Int_t pdgMeson,AliRDHFCuts *cuts, Bool_t switchPPb):
135 AliAnalysisTaskSE(name),
141 fHistNtrEta16vsNtrEta1EvSel(0),
142 fHistNtrEta05vsNtrEta1EvSel(0),
143 fHistNtrEta03vsNtrEta1EvSel(0),
144 fHistNtrEtaV0AvsNtrEta1EvSel(0),
145 fHistNtrEtaV0MvsNtrEta1EvSel(0),
146 fHistNtrCorrEta1vsNtrRawEta1EvSel(0),
147 fHistNtrEta16vsNtrEta1EvWithCand(0),
148 fHistNtrEta05vsNtrEta1EvWithCand(0),
149 fHistNtrEta03vsNtrEta1EvWithCand(0),
150 fHistNtrEtaV0AvsNtrEta1EvWithCand(0),
151 fHistNtrEtaV0MvsNtrEta1EvWithCand(0),
152 fHistNtrCorrEta1vsNtrRawEta1EvWithCand(0),
153 fHistNtrEta16vsNtrEta1EvWithD(0),
154 fHistNtrEta05vsNtrEta1EvWithD(0),
155 fHistNtrEta03vsNtrEta1EvWithD(0),
156 fHistNtrEtaV0AvsNtrEta1EvWithD(0),
157 fHistNtrEtaV0MvsNtrEta1EvWithD(0),
158 fHistNtrCorrEta1vsNtrRawEta1EvWithD(0),
160 fHistNtrCorrVsZvtx(0),
162 fHistNtrCorrVsNchMC(0),
163 fHistNtrVsNchMCPrimary(0),
164 fHistNtrCorrVsNchMCPrimary(0),
165 fHistNtrVsNchMCPhysicalPrimary(0),
166 fHistNtrCorrVsNchMCPhysicalPrimary(0),
167 fHistGenPrimaryParticlesInelGt0(0),
168 fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary(0),
169 fHistNtrUnCorrPSSel(0),
170 fHistNtrUnCorrEvSel(0),
171 fHistNtrUnCorrEvWithCand(0),
172 fHistNtrUnCorrEvWithD(0),
173 fHistNtrCorrPSSel(0),
174 fHistNtrCorrEvSel(0),
175 fHistNtrCorrEvWithCand(0),
176 fHistNtrCorrEvWithD(0),
178 fPtVsMassVsMultNoPid(0),
179 fPtVsMassVsMultUncorr(0),
180 fPtVsMassVsMultPart(0),
181 fPtVsMassVsMultAntiPart(0),
182 fPtVsMassVsMultMC(0),
184 fLowmasslimit(1.765),
186 fRDCutsAnalysis(cuts),
191 fLowerImpPar(-2000.),
192 fHigherImpPar(2000.),
195 fisPPbData(switchPPb),
197 fSubtractTrackletsFromDau(kFALSE),
198 fKeepCorrPlots(kFALSE),
199 fUseNchWeight(kFALSE),
204 fMultiplicityEstimator(kNtrk10),
205 fMCPrimariesEstimator(kEta10)
208 // Standard constructor
211 for(Int_t i=0; i<5; i++) fHistMassPtImpPar[i]=0;
212 for(Int_t i=0; i<4; i++) fMultEstimatorAvg[i]=0;
215 SetMassLimits(0.12,0.2);
218 SetMassLimits(fPdgMeson,0.1);
220 // Default constructor
221 // Otput slot #1 writes into a TList container
222 DefineOutput(1,TList::Class()); //My private output
223 // Output slot #2 writes cut to private output
224 DefineOutput(2,TList::Class());
225 // Output slot #3 writes cut to private output
226 DefineOutput(3,TList::Class());
227 // Output slot #4 writes cut to private output
228 DefineOutput(4,TList::Class());
230 //________________________________________________________________________
231 AliAnalysisTaskSEDvsMultiplicity::~AliAnalysisTaskSEDvsMultiplicity()
239 delete fListProfiles;
240 delete fRDCutsAnalysis;
243 for(Int_t i=0; i<4; i++) {
244 if (fMultEstimatorAvg[i]) delete fMultEstimatorAvg[i];
247 for(Int_t i=0; i<5; i++){
248 delete fHistMassPtImpPar[i];
250 if(fHistoMCNch) delete fHistoMCNch;
251 if(fHistoMeasNch) delete fHistoMeasNch;
254 //_________________________________________________________________
255 void AliAnalysisTaskSEDvsMultiplicity::SetMassLimits(Double_t lowlimit, Double_t uplimit){
256 // set invariant mass limits
257 if(uplimit>lowlimit){
258 fLowmasslimit = lowlimit;
259 fUpmasslimit = uplimit;
261 AliError("Wrong mass limits: upper value should be larger than lower one");
264 //_________________________________________________________________
265 void AliAnalysisTaskSEDvsMultiplicity::SetMassLimits(Int_t pdg, Double_t range){
266 // set invariant mass limits
267 Double_t mass=TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdg))->Mass();
268 SetMassLimits(mass-range,mass+range);
270 //________________________________________________________________________
271 void AliAnalysisTaskSEDvsMultiplicity::Init(){
275 printf("AnalysisTaskSEDvsMultiplicity::Init() \n");
277 if(fUseNchWeight && !fReadMC){ AliFatal("Nch weights can only be used in MC mode"); return; }
278 if(fUseNchWeight && !fHistoMCNch){ AliFatal("Nch weights can only be used without histogram"); return; }
280 fListCuts=new TList();
281 fListCuts->SetOwner();
282 fListCuts->SetName("CutsList");
286 AliRDHFCutsDplustoKpipi* copycut=new AliRDHFCutsDplustoKpipi(*(static_cast<AliRDHFCutsDplustoKpipi*>(fRDCutsAnalysis)));
287 copycut->SetName("AnalysisCutsDplus");
288 fListCuts->Add(copycut);
289 }else if(fPdgMeson==421){
290 AliRDHFCutsD0toKpi* copycut=new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fRDCutsAnalysis)));
291 copycut->SetName("AnalysisCutsDzero");
292 fListCuts->Add(copycut);
293 }else if(fPdgMeson==413){
294 AliRDHFCutsDStartoKpipi* copycut=new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fRDCutsAnalysis)));
295 copycut->SetName("AnalysisCutsDStar");
296 fListCuts->Add(copycut);
298 PostData(2,fListCuts);
300 fListProfiles = new TList();
301 fListProfiles->SetOwner();
304 if (fisPPbData) {period[0]="LHC13b"; period[1]="LHC13c"; nProfiles = 2;}
305 else {period[0]="LHC10b"; period[1]="LHC10c"; period[2]="LHC10d"; period[3]="LHC10e"; nProfiles = 4;}
307 for(Int_t i=0; i<nProfiles; i++){
308 if(fMultEstimatorAvg[i]){
309 TProfile* hprof=new TProfile(*fMultEstimatorAvg[i]);
310 hprof->SetName(Form("ProfileTrkVsZvtx%s\n",period[i].Data()));
311 fListProfiles->Add(hprof);
314 PostData(4,fListProfiles);
319 //________________________________________________________________________
320 void AliAnalysisTaskSEDvsMultiplicity::UserCreateOutputObjects()
322 // Create the output container
324 if(fDebug > 1) printf("AnalysisTaskSEDvsMultiplicity::UserCreateOutputObjects() \n");
326 // Several histograms are more conveniently managed in a TList
327 fOutput = new TList();
329 fOutput->SetName("OutputHistos");
331 Int_t nMultBins = 200;
332 Float_t firstMultBin = -0.5;
333 Float_t lastMultBin = 199.5;
334 Int_t nMultBinsNtrk = nMultBins;
335 Float_t lastMultBinNtrk = lastMultBin;
336 Int_t nMultBinsV0 = 400;
337 Float_t lastMultBinV0 = 799.5;
340 lastMultBinNtrk = 374.5;
341 nMultBins = nMultBinsNtrk;
342 lastMultBin = lastMultBinNtrk;
344 if(fMultiplicityEstimator==kVZERO || fMultiplicityEstimator==kVZEROA) {
345 nMultBins = nMultBinsV0;
346 lastMultBin = lastMultBinV0;
349 fHistNtrUnCorrPSSel = new TH1F("hNtrUnCorrPSSel","Uncorrected tracklets multiplicity for PS selected events; Tracklets ; Entries",nMultBins,firstMultBin,lastMultBin);
350 fHistNtrUnCorrEvSel = new TH1F("hNtrUnCorrEvSel","Uncorrected tracklets multiplicity for selected events; Tracklets ; Entries",nMultBins,firstMultBin,lastMultBin);
351 fHistNtrUnCorrEvWithCand = new TH1F("hNtrUnCorrEvWithCand", "Uncorrected Tracklets multiplicity for events with D candidates; Tracklets ; Entries",nMultBins,firstMultBin,lastMultBin);// Total multiplicity
352 fHistNtrUnCorrEvWithD = new TH1F("hNtrUnCorrEvWithD","Uncorrected Tracklets multiplicity for events with D in mass region ; Tracklets ; Entries",nMultBins,firstMultBin,lastMultBin); //
353 fHistNtrCorrPSSel = new TH1F("hNtrCorrPSSel","Corrected tracklets multiplicity for PS selected events; Tracklets ; Entries",nMultBins,firstMultBin,lastMultBin);
354 fHistNtrCorrEvSel = new TH1F("hNtrCorrEvSel","Corrected tracklets multiplicity for selected events; Tracklets ; Entries",nMultBins,firstMultBin,lastMultBin);
355 fHistNtrCorrEvWithCand = new TH1F("hNtrCorrEvWithCand", "Tracklets multiplicity for events with D candidates; Tracklets ; Entries",nMultBins,firstMultBin,lastMultBin);// Total multiplicity
356 fHistNtrCorrEvWithD = new TH1F("hNtrCorrEvWithD", "Tracklets multiplicity for events with D in mass region ; Tracklets ; Entries",nMultBins,firstMultBin,lastMultBin); //
359 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
360 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
361 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
362 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
363 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
364 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
366 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
367 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
368 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
369 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
370 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
371 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
373 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
374 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
375 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
376 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
377 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
378 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
380 fHistNtrVsZvtx = new TH2F("hNtrVsZvtx","Ntracklet vs VtxZ; VtxZ;N_{tracklet};",300,-15,15,nMultBins,firstMultBin,lastMultBin); //
381 fHistNtrCorrVsZvtx = new TH2F("hNtrCorrVsZvtx","Ntracklet vs VtxZ; VtxZ;N_{tracklet};",300,-15,15,nMultBins,firstMultBin,lastMultBin); //
383 fHistNtrVsNchMC = new TH2F("hNtrVsNchMC","Ntracklet vs NchMC; Nch;N_{tracklet};",nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); //
384 fHistNtrCorrVsNchMC = new TH2F("hNtrCorrVsNchMC","Ntracklet vs Nch; Nch;N_{tracklet};",nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); //
386 fHistNtrVsNchMCPrimary = new TH2F("hNtrVsNchMCPrimary","Ntracklet vs Nch (Primary); Nch (Primary);N_{tracklet};",nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); //
387 fHistNtrCorrVsNchMCPrimary = new TH2F("hNtrCorrVsNchMCPrimary","Ntracklet vs Nch (Primary); Nch(Primary) ;N_{tracklet};",nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); //
389 fHistNtrVsNchMCPhysicalPrimary = new TH2F("hNtrVsNchMCPhysicalPrimary","Ntracklet vs Nch (Physical Primary); Nch (Physical Primary);N_{tracklet};",nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); //
390 fHistNtrCorrVsNchMCPhysicalPrimary = new TH2F("hNtrCorrVsMCPhysicalPrimary","Ntracklet vs Nch (Physical Primary); Nch (Physical Primary);N_{tracklet};",nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); //
392 fHistGenPrimaryParticlesInelGt0 = new TH1F("hGenPrimaryParticlesInelGt0","Multiplcity of generated charged particles ; Nparticles ; Entries",nMultBins,firstMultBin,lastMultBin);
394 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);
396 fHistNtrUnCorrPSSel->Sumw2();
397 fHistNtrUnCorrEvSel->Sumw2();
398 fHistNtrUnCorrEvWithCand->Sumw2();
399 fHistNtrUnCorrEvWithD->Sumw2();
400 fHistNtrCorrPSSel->Sumw2();
401 fHistNtrCorrEvSel->Sumw2();
402 fHistNtrCorrEvWithCand->Sumw2();
403 fHistNtrCorrEvWithD->Sumw2();
404 fHistGenPrimaryParticlesInelGt0->Sumw2();
405 fOutput->Add(fHistNtrUnCorrPSSel);
406 fOutput->Add(fHistNtrUnCorrEvSel);
407 fOutput->Add(fHistNtrUnCorrEvWithCand);
408 fOutput->Add(fHistNtrUnCorrEvWithD);
409 fOutput->Add(fHistNtrCorrPSSel);
410 fOutput->Add(fHistNtrCorrEvSel);
411 fOutput->Add(fHistNtrCorrEvWithCand);
412 fOutput->Add(fHistNtrCorrEvWithD);
414 fOutput->Add(fHistNtrEta16vsNtrEta1EvSel);
415 fOutput->Add(fHistNtrEta05vsNtrEta1EvSel);
416 fOutput->Add(fHistNtrEta03vsNtrEta1EvSel);
417 fOutput->Add(fHistNtrEtaV0AvsNtrEta1EvSel);
418 fOutput->Add(fHistNtrEtaV0MvsNtrEta1EvSel);
419 fOutput->Add(fHistNtrCorrEta1vsNtrRawEta1EvSel);
420 fOutput->Add(fHistNtrEta16vsNtrEta1EvWithCand);
421 fOutput->Add(fHistNtrEta05vsNtrEta1EvWithCand);
422 fOutput->Add(fHistNtrEta03vsNtrEta1EvWithCand);
423 fOutput->Add(fHistNtrEtaV0AvsNtrEta1EvWithCand);
424 fOutput->Add(fHistNtrEtaV0MvsNtrEta1EvWithCand);
425 fOutput->Add(fHistNtrCorrEta1vsNtrRawEta1EvWithCand);
426 fOutput->Add(fHistNtrEta16vsNtrEta1EvWithD);
427 fOutput->Add(fHistNtrEta05vsNtrEta1EvWithD);
428 fOutput->Add(fHistNtrEta03vsNtrEta1EvWithD);
429 fOutput->Add(fHistNtrEtaV0AvsNtrEta1EvWithD);
430 fOutput->Add(fHistNtrEtaV0MvsNtrEta1EvWithD);
431 fOutput->Add(fHistNtrCorrEta1vsNtrRawEta1EvWithD);
433 fOutput->Add(fHistNtrVsZvtx);
434 fOutput->Add(fHistNtrCorrVsZvtx);
436 fOutput->Add(fHistNtrVsNchMC);
437 fOutput->Add(fHistNtrCorrVsNchMC);
438 fOutput->Add(fHistNtrVsNchMCPrimary);
439 fOutput->Add(fHistNtrCorrVsNchMCPrimary);
440 fOutput->Add(fHistNtrVsNchMCPhysicalPrimary);
441 fOutput->Add(fHistNtrCorrVsNchMCPhysicalPrimary);
442 fOutput->Add(fHistGenPrimaryParticlesInelGt0);
443 fOutput->Add(fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary);
446 fHistNEvents = new TH1F("fHistNEvents", "number of events ",11,-0.5,10.5);
447 fHistNEvents->GetXaxis()->SetBinLabel(1,"nEvents total");
448 fHistNEvents->GetXaxis()->SetBinLabel(2,"nEvents with Z vertex");
449 fHistNEvents->GetXaxis()->SetBinLabel(3,"nEvents selected");
450 fHistNEvents->GetXaxis()->SetBinLabel(4,"Rejected due to trigger");
451 fHistNEvents->GetXaxis()->SetBinLabel(5,"Rejected due to phys sel");
452 fHistNEvents->GetXaxis()->SetBinLabel(6,"Rejected due to vertex cuts");
453 fHistNEvents->GetXaxis()->SetBinLabel(7,"Rejected due to pileup");
454 fHistNEvents->GetXaxis()->SetBinLabel(8,"Total no. of candidate");
455 fHistNEvents->GetXaxis()->SetBinLabel(9,"no. of cand wo bitmask");
456 fHistNEvents->GetXaxis()->SetBinLabel(10,"D after cuts (No PID)");
457 fHistNEvents->GetXaxis()->SetBinLabel(11,"D after cuts + PID)");
458 fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);
459 fHistNEvents->Sumw2();
460 fHistNEvents->SetMinimum(0);
461 fOutput->Add(fHistNEvents);
463 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.);
465 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.);
467 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.);
469 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.);
471 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.);
473 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.);
475 fOutput->Add(fPtVsMassVsMult);
476 fOutput->Add(fPtVsMassVsMultUncorr);
477 fOutput->Add(fPtVsMassVsMultNoPid);
478 fOutput->Add(fPtVsMassVsMultPart);
479 fOutput->Add(fPtVsMassVsMultAntiPart);
480 fOutput->Add(fPtVsMassVsMultMC);
482 if(fDoImpPar) CreateImpactParameterHistos();
484 fCounter = new AliNormalizationCounter("NormCounterCorrMult");
485 fCounter->SetStudyMultiplicity(kTRUE,1.);
488 fCounterU = new AliNormalizationCounter("NormCounterUnCorrMult");
489 fCounterU->SetStudyMultiplicity(kTRUE,1.);
492 fOutputCounters = new TList();
493 fOutputCounters->SetOwner();
494 fOutputCounters->SetName("OutputCounters");
495 fOutputCounters->Add(fCounter);
496 fOutputCounters->Add(fCounterU);
499 PostData(2,fListCuts);
500 PostData(3,fOutputCounters);
501 PostData(4,fListProfiles);
503 if(fUseNchWeight) CreateMeasuredNchHisto();
510 //________________________________________________________________________
511 void AliAnalysisTaskSEDvsMultiplicity::UserExec(Option_t */*option*/)
513 // Execute analysis for current event:
514 // heavy flavor candidates association to MC truth
516 AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
518 // AliAODTracklets* tracklets = aod->GetTracklets();
519 //Int_t ntracklets = tracklets->GetNumberOfTracklets();
522 TClonesArray *arrayCand = 0;
523 TString arrayName="";
528 arrayName="Charm3Prong";
529 pdgDau[0]=211; pdgDau[1]=321; pdgDau[2]=211;
531 selbit=AliRDHFCuts::kDplusCuts;
532 }else if(fPdgMeson==421){
534 pdgDau[0]=211; pdgDau[1]=321; pdgDau[2]=0;
536 selbit=AliRDHFCuts::kD0toKpiCuts;
537 }else if(fPdgMeson==413){
539 pdgDau[0]=321; pdgDau[1]=211; pdgDau[2]=0; // Quoting here D0 daughters (D* ones on another variable later)
541 selbit=AliRDHFCuts::kDstarCuts;
544 if(!aod && AODEvent() && IsStandardAOD()) {
545 // In case there is an AOD handler writing a standard AOD, use the AOD
546 // event in memory rather than the input (ESD) event.
547 aod = dynamic_cast<AliAODEvent*> (AODEvent());
548 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
549 // have to taken from the AOD event hold by the AliAODExtension
550 AliAODHandler* aodHandler = (AliAODHandler*)
551 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
552 if(aodHandler->GetExtensions()) {
553 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
554 AliAODEvent *aodFromExt = ext->GetAOD();
555 arrayCand=(TClonesArray*)aodFromExt->GetList()->FindObject(arrayName.Data());
558 arrayCand=(TClonesArray*)aod->GetList()->FindObject(arrayName.Data());
561 if(!aod || !arrayCand) {
562 printf("AliAnalysisTaskSEDvsMultiplicity::UserExec: Charm3Prong branch not found!\n");
566 if(fisPPbData && fReadMC){
567 Int_t runnumber = aod->GetRunNumber();
568 if(aod->GetTriggerMask()==0 &&
569 (runnumber>=195344 && runnumber<=195677)){
570 AliDebug(3,"Event rejected because of null trigger mask");
576 // fix for temporary bug in ESDfilter
577 // the AODs with null vertex pointer didn't pass the PhysSel
578 if(!aod->GetPrimaryVertex()||TMath::Abs(aod->GetMagneticField())<0.001) return;
580 Int_t countTreta1=AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aod,-1.,1.);
581 Int_t countTreta03=AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aod,-0.3,0.3);
582 Int_t countTreta05=AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aod,-0.5,0.5);
583 Int_t countTreta16=AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aod,-1.6,1.6);
587 AliAODVZERO *vzeroAOD = (AliAODVZERO*)aod->GetVZEROData();
589 vzeroMult = vzeroAOD->GetMTotV0A() + vzeroAOD->GetMTotV0C();
590 vzeroMultA = vzeroAOD->GetMTotV0A();
593 Int_t countMult = countTreta1;
594 if(fMultiplicityEstimator==kNtrk03) { countMult = countTreta03; }
595 if(fMultiplicityEstimator==kNtrk05) { countMult = countTreta05; }
596 if(fMultiplicityEstimator==kNtrk10to16) { countMult = countTreta16 - countTreta1; }
597 if(fMultiplicityEstimator==kVZERO) { countMult = vzeroMult; }
598 if(fMultiplicityEstimator==kVZEROA) { countMult = vzeroMultA; }
601 fCounterU->StoreEvent(aod,fRDCutsAnalysis,fReadMC,countMult);
602 fHistNEvents->Fill(0); // count event
604 Double_t countTreta1corr=countTreta1;
605 Double_t countCorr=countMult;
606 AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
608 // FIX ME: No correction to the VZERO !!
609 if(vtx1 && (fMultiplicityEstimator!=kVZERO) && (fMultiplicityEstimator!=kVZEROA)){
610 if(vtx1->GetNContributors()>0){
611 fHistNEvents->Fill(1);
612 TProfile* estimatorAvg = GetEstimatorHistogram(aod);
614 countTreta1corr=AliVertexingHFUtils::GetCorrectedNtracklets(estimatorAvg,countTreta1,vtx1->GetZ(),fRefMult);
615 countCorr=AliVertexingHFUtils::GetCorrectedNtracklets(estimatorAvg,countMult,vtx1->GetZ(),fRefMult);
621 Bool_t isEvSel=fRDCutsAnalysis->IsEventSelected(aod);
623 if(fRDCutsAnalysis->GetWhyRejection()==5) fHistNEvents->Fill(3);
624 if(fRDCutsAnalysis->GetWhyRejection()==7) fHistNEvents->Fill(4);
625 if(fRDCutsAnalysis->GetWhyRejection()==6) fHistNEvents->Fill(5);
626 if(fRDCutsAnalysis->GetWhyRejection()==1) fHistNEvents->Fill(6);
628 Bool_t isEvPSRejected = fRDCutsAnalysis->IsEventRejectedDuePhysicsSelection();
630 fHistNtrUnCorrPSSel->Fill(countMult);
631 fHistNtrCorrPSSel->Fill(countCorr);
636 fHistNtrEta16vsNtrEta1EvSel->Fill(countTreta1,countTreta16);
637 fHistNtrEta05vsNtrEta1EvSel->Fill(countTreta1,countTreta05);
638 fHistNtrEta03vsNtrEta1EvSel->Fill(countTreta1,countTreta03);
639 fHistNtrEtaV0AvsNtrEta1EvSel->Fill(countTreta1,vzeroMultA);
640 fHistNtrEtaV0MvsNtrEta1EvSel->Fill(countTreta1,vzeroMult);
641 fHistNtrCorrEta1vsNtrRawEta1EvSel->Fill(countTreta1,countTreta1corr);
644 fHistNtrVsZvtx->Fill(vtx1->GetZ(),countMult);
645 fHistNtrCorrVsZvtx->Fill(vtx1->GetZ(),countCorr);
648 TClonesArray *arrayMC=0;
649 AliAODMCHeader *mcHeader=0;
651 Double_t nchWeight=1.0;
656 arrayMC = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
658 printf("AliAnalysisTaskSEDvsMultiplicity::UserExec: MC particles branch not found!\n");
662 mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
664 printf("AliAnalysisTaskSEDvsMultiplicity::UserExec: MC header branch not found!\n");
669 Int_t nChargedMC=AliVertexingHFUtils::GetGeneratedMultiplicityInEtaRange(arrayMC,-1.0,1.0);
670 Int_t nChargedMCPrimary=AliVertexingHFUtils::GetGeneratedPrimariesInEtaRange(arrayMC,-1.0,1.0);
671 Int_t nChargedMCPhysicalPrimary=AliVertexingHFUtils::GetGeneratedPhysicalPrimariesInEtaRange(arrayMC,-1.0,1.0);
673 // Compute the Nch weights (reference is Ntracklets within |eta|<1.0)
675 Double_t tmpweight = 1.0;
676 if(nChargedMCPhysicalPrimary<=0) tmpweight = 0.0;
678 Double_t pMeas = fHistoMeasNch->GetBinContent(fHistoMeasNch->FindBin(nChargedMCPhysicalPrimary));
679 // printf(" pMeas=%2.2f and histo MCNch %s \n",pMeas,fHistoMCNch);
680 Double_t pMC = fHistoMCNch->GetBinContent(fHistoMCNch->FindBin(nChargedMCPhysicalPrimary));
681 tmpweight = pMC>0 ? pMeas/pMC : 0.;
683 nchWeight *= tmpweight;
684 AliDebug(2,Form("Using Nch weights, Mult=%d Weight=%f\n",nChargedMCPhysicalPrimary,nchWeight));
687 // Now recompute the variables in case another MC estimator is considered
688 Int_t nChargedMCEta10 = nChargedMC;
689 Int_t nChargedMCEta16 = AliVertexingHFUtils::GetGeneratedMultiplicityInEtaRange(arrayMC,-1.6,1.6);
690 Int_t nChargedMCEta05 = AliVertexingHFUtils::GetGeneratedMultiplicityInEtaRange(arrayMC,-0.5,0.5);
691 Int_t nChargedMCEta03 = AliVertexingHFUtils::GetGeneratedMultiplicityInEtaRange(arrayMC,-0.3,0.3);
692 Int_t nChargedMCEtam37tm17 = AliVertexingHFUtils::GetGeneratedMultiplicityInEtaRange(arrayMC,-3.7,-1.7);
693 Int_t nChargedMCEta28t51 = AliVertexingHFUtils::GetGeneratedMultiplicityInEtaRange(arrayMC,2.8,5.1);
694 Int_t nChargedMCPrimaryEta10 = nChargedMCPrimary;
695 Int_t nChargedMCPrimaryEta16 = AliVertexingHFUtils::GetGeneratedPrimariesInEtaRange(arrayMC,-1.6,1.6);
696 Int_t nChargedMCPrimaryEta05 = AliVertexingHFUtils::GetGeneratedPrimariesInEtaRange(arrayMC,-0.5,0.5);
697 Int_t nChargedMCPrimaryEta03 = AliVertexingHFUtils::GetGeneratedPrimariesInEtaRange(arrayMC,-0.3,0.3);
698 Int_t nChargedMCPrimaryEtam37tm17 = AliVertexingHFUtils::GetGeneratedPrimariesInEtaRange(arrayMC,-3.7,-1.7);
699 Int_t nChargedMCPrimaryEta28t51 = AliVertexingHFUtils::GetGeneratedPrimariesInEtaRange(arrayMC,2.8,5.1);
700 Int_t nChargedMCPhysicalPrimaryEta10 = nChargedMCPhysicalPrimary;
701 Int_t nChargedMCPhysicalPrimaryEta16 = AliVertexingHFUtils::GetGeneratedPhysicalPrimariesInEtaRange(arrayMC,-1.6,1.6);
702 Int_t nChargedMCPhysicalPrimaryEta05 = AliVertexingHFUtils::GetGeneratedPhysicalPrimariesInEtaRange(arrayMC,-0.5,0.5);
703 Int_t nChargedMCPhysicalPrimaryEta03 = AliVertexingHFUtils::GetGeneratedPhysicalPrimariesInEtaRange(arrayMC,-0.3,0.3);
704 Int_t nChargedMCPhysicalPrimaryEtam37tm17 = AliVertexingHFUtils::GetGeneratedPhysicalPrimariesInEtaRange(arrayMC,-3.7,-1.7);
705 Int_t nChargedMCPhysicalPrimaryEta28t51 = AliVertexingHFUtils::GetGeneratedPhysicalPrimariesInEtaRange(arrayMC,2.8,5.1);
706 if(fMCPrimariesEstimator==kEta10to16){
707 nChargedMC = nChargedMCEta16 - nChargedMCEta10;
708 nChargedMCPrimary = nChargedMCPrimaryEta16 - nChargedMCPrimaryEta10;
709 nChargedMCPhysicalPrimary = nChargedMCPhysicalPrimaryEta16 - nChargedMCPhysicalPrimaryEta10;
710 } else if(fMCPrimariesEstimator==kEta05){
711 nChargedMC = nChargedMCEta05;
712 nChargedMCPrimary = nChargedMCPrimaryEta05;
713 nChargedMCPhysicalPrimary = nChargedMCPhysicalPrimaryEta05;
714 } else if(fMCPrimariesEstimator==kEta03){
715 nChargedMC = nChargedMCEta03;
716 nChargedMCPrimary = nChargedMCPrimaryEta03;
717 nChargedMCPhysicalPrimary = nChargedMCPhysicalPrimaryEta03;
718 } else if(fMCPrimariesEstimator==kEtaVZERO){
719 nChargedMC = nChargedMCEtam37tm17 + nChargedMCEta28t51;
720 nChargedMCPrimary = nChargedMCPrimaryEtam37tm17 + nChargedMCPrimaryEta28t51;
721 nChargedMCPhysicalPrimary = nChargedMCPhysicalPrimaryEtam37tm17 + nChargedMCPhysicalPrimaryEta28t51;
722 } else if(fMCPrimariesEstimator==kEtaVZEROA){
723 nChargedMC = nChargedMCEta28t51;
724 nChargedMCPrimary = nChargedMCPrimaryEta28t51;
725 nChargedMCPhysicalPrimary = nChargedMCPhysicalPrimaryEta28t51;
728 // Here fill the MC correlation plots
729 if(nChargedMCPhysicalPrimary>0){ // INEL>0 for |eta|<1
730 fHistGenPrimaryParticlesInelGt0->Fill(nChargedMCPhysicalPrimary,nchWeight);
733 fHistNtrVsNchMC->Fill(nChargedMC,countMult,nchWeight);
734 fHistNtrCorrVsNchMC->Fill(nChargedMC,countCorr,nchWeight);
736 fHistNtrVsNchMCPrimary->Fill(nChargedMCPrimary,countMult,nchWeight);
737 fHistNtrCorrVsNchMCPrimary->Fill(nChargedMCPrimary,countCorr,nchWeight);
739 fHistNtrVsNchMCPhysicalPrimary->Fill(nChargedMCPhysicalPrimary,countMult,nchWeight);
740 fHistNtrCorrVsNchMCPhysicalPrimary->Fill(nChargedMCPhysicalPrimary,countCorr,nchWeight);
742 fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary->Fill(nChargedMC,nChargedMCPrimary,nChargedMCPhysicalPrimary,nchWeight);
745 Int_t nCand = arrayCand->GetEntriesFast();
746 Int_t nSelectedNoPID=0,nSelectedPID=0,nSelectedInMassPeak=0;
747 Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();
748 Double_t mDplusPDG = TDatabasePDG::Instance()->GetParticle(411)->Mass();
749 Double_t mDstarPDG = TDatabasePDG::Instance()->GetParticle(413)->Mass();
751 // pdg of daughters needed for D* too
752 UInt_t pdgDgDStartoD0pi[2]={421,211};
755 Double_t nSelCand=0.;
756 for (Int_t iCand = 0; iCand < nCand; iCand++) {
757 AliAODRecoDecayHF *d = (AliAODRecoDecayHF*)arrayCand->UncheckedAt(iCand);
758 AliAODRecoCascadeHF *dCascade = NULL;
759 if(fPdgMeson==413) dCascade = (AliAODRecoCascadeHF*)d;
761 fHistNEvents->Fill(7);
762 if(fUseBit && !d->HasSelectionBit(selbit)){
763 fHistNEvents->Fill(8);
767 Double_t ptCand = d->Pt();
768 Double_t rapid=d->Y(fPdgMeson);
769 Bool_t isFidAcc=fRDCutsAnalysis->IsInFiducialAcceptance(ptCand,rapid);
770 if(!isFidAcc) continue;
775 labD = dCascade->MatchToMC(fPdgMeson,421,(Int_t*)pdgDgDStartoD0pi,(Int_t*)pdgDau,arrayMC);
777 labD = d->MatchToMC(fPdgMeson,arrayMC,nDau,(Int_t*)pdgDau);
779 FillMCMassHistos(arrayMC,labD, countMult,nchWeight);
782 Int_t passAllCuts=fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kAll,aod);
783 Int_t passTopolCuts=fRDCutsAnalysis->GetIsSelectedCuts();
784 if(passTopolCuts==0) continue;
786 fHistNEvents->Fill(9);
789 fHistNEvents->Fill(10);
791 Double_t multForCand = countCorr;
793 if(fSubtractTrackletsFromDau){
794 // For the D* case, subtract only the D0 daughter tracks <=== FIXME !!
795 AliAODRecoDecayHF2Prong* d0fromDstar = NULL;
796 if(fPdgMeson==413) d0fromDstar = (AliAODRecoDecayHF2Prong*)dCascade->Get2Prong();
798 for(Int_t iDau=0; iDau<nDau; iDau++){
799 AliAODTrack *t = NULL;
800 if(fPdgMeson==413){ t = (AliAODTrack*)d0fromDstar->GetDaughter(iDau); }
801 else{ t = (AliAODTrack*)d->GetDaughter(iDau); }
803 if(t->HasPointOnITSLayer(0) && t->HasPointOnITSLayer(1)){
804 if(multForCand>0) multForCand-=1;
808 Bool_t isPrimary=kTRUE;
809 Double_t trueImpParXY=9999.;
810 Double_t impparXY=d->ImpParXY()*10000.;
811 Double_t dlen=0.1; //FIXME
814 mass[0]=d->InvMass(nDau,pdgDau);
816 if(TMath::Abs(mass[0]-mDplusPDG)<0.02) nSelectedInMassPeak++; //20 MeV for now... FIXME
817 }else if(fPdgMeson==421){
818 UInt_t pdgdaughtersD0[2]={211,321};//pi,K
819 UInt_t pdgdaughtersD0bar[2]={321,211};//K,pi
820 mass[0]=d->InvMass(2,pdgdaughtersD0);
821 mass[1]=d->InvMass(2,pdgdaughtersD0bar);
822 if(TMath::Abs(mass[0]-mD0PDG)<0.02 || TMath::Abs(mass[1]-mD0PDG)<0.02 ) nSelectedInMassPeak++; //20 MeV for now... FIXME
823 }else if(fPdgMeson==413){
825 mass[0]=dCascade->DeltaInvMass();
827 if(TMath::Abs(mass[0]-(mDstarPDG-mD0PDG))<0.0015) nSelectedInMassPeak++; //1 MeV for now... FIXME
829 for(Int_t iHyp=0; iHyp<2; iHyp++){
830 if(mass[iHyp]<0.) continue; // for D+ and D* we have 1 mass hypothesis
831 Double_t invMass=mass[iHyp];
832 Double_t arrayForSparse[5]={invMass,ptCand,impparXY,dlen,multForCand};
837 labD = dCascade->MatchToMC(fPdgMeson,421,(Int_t*)pdgDgDStartoD0pi,(Int_t*)pdgDau,arrayMC);
839 labD = d->MatchToMC(fPdgMeson,arrayMC,nDau,(Int_t*)pdgDau);
842 Bool_t fillHisto=fDoImpPar;
844 AliAODMCParticle *partD = (AliAODMCParticle*)arrayMC->At(labD);
845 Int_t code=partD->GetPdgCode();
846 if(CheckOrigin(arrayMC,partD)==5) isPrimary=kFALSE;
847 if(code<0 && iHyp==0) fillHisto=kFALSE;
848 if(code>0 && iHyp==1) fillHisto=kFALSE;
851 trueImpParXY=AliVertexingHFUtils::GetTrueImpactParameterDplus(mcHeader,arrayMC,partD)*10000.;
852 }else if(fPdgMeson==421){
853 trueImpParXY=AliVertexingHFUtils::GetTrueImpactParameterDzero(mcHeader,arrayMC,partD)*10000.;
854 }else if(fPdgMeson==413){
855 trueImpParXY=0.; /// FIXME
857 Double_t arrayForSparseTrue[5]={invMass,ptCand,trueImpParXY,dlen,multForCand};
858 if(fillHisto && passAllCuts){
859 fHistMassPtImpPar[2]->Fill(arrayForSparse);
860 fHistMassPtImpPar[3]->Fill(arrayForSparseTrue);
863 if(fillHisto && passAllCuts) fHistMassPtImpPar[1]->Fill(arrayForSparse);
866 if(fillHisto && passAllCuts)fHistMassPtImpPar[4]->Fill(arrayForSparse);
869 if(TMath::Abs(labD)==fPdgMeson && fMCOption==2) continue;
870 if(TMath::Abs(labD)!=fPdgMeson && fMCOption==1) continue;
875 if(iHyp==0 && !(passTopolCuts&1)) continue; // candidate not passing as D0
876 if(iHyp==1 && !(passTopolCuts&2)) continue; // candidate not passing as D0bar
879 fPtVsMassVsMultNoPid->Fill(multForCand,invMass,ptCand);
882 if(iHyp==0 && !(passAllCuts&1)) continue; // candidate not passing as D0
883 if(iHyp==1 && !(passAllCuts&2)) continue; // candidate not passing as D0bar
886 aveMult+=multForCand;
888 fPtVsMassVsMult->Fill(multForCand,invMass,ptCand,nchWeight);
889 fPtVsMassVsMultUncorr->Fill(countTreta1,invMass,ptCand,nchWeight);
890 // Add separation between part antipart
892 if(d->GetCharge()>0) fPtVsMassVsMultPart->Fill(multForCand,invMass,ptCand,nchWeight);
893 else fPtVsMassVsMultAntiPart->Fill(multForCand,invMass,ptCand,nchWeight);
894 }else if(fPdgMeson==421){
895 if(passAllCuts&1) fPtVsMassVsMultPart->Fill(multForCand,invMass,ptCand,nchWeight);
896 if(passAllCuts&2) fPtVsMassVsMultAntiPart->Fill(multForCand,invMass,ptCand,nchWeight);
897 }else if(fPdgMeson==413){
898 if(d->GetCharge()>0) fPtVsMassVsMultPart->Fill(multForCand,invMass,ptCand,nchWeight);
899 else fPtVsMassVsMultAntiPart->Fill(multForCand,invMass,ptCand,nchWeight);
903 fHistMassPtImpPar[0]->Fill(arrayForSparse);
910 if(fSubtractTrackletsFromDau && nSelCand>0){
912 fCounter->StoreEvent(aod,fRDCutsAnalysis,fReadMC,(Int_t)(aveMult+0.5001));
914 fCounter->StoreEvent(aod,fRDCutsAnalysis,fReadMC,(Int_t)countCorr);
918 fCounter->StoreCandidates(aod,nSelectedNoPID,kTRUE);
919 fCounter->StoreCandidates(aod,nSelectedPID,kFALSE);
920 fHistNtrUnCorrEvSel->Fill(countMult,nchWeight);
921 fHistNtrCorrEvSel->Fill(countCorr,nchWeight);
923 fHistNtrUnCorrEvWithCand->Fill(countMult,nchWeight);
924 fHistNtrCorrEvWithCand->Fill(countCorr,nchWeight);
926 fHistNtrEta16vsNtrEta1EvWithCand->Fill(countTreta1,countTreta16);
927 fHistNtrEta05vsNtrEta1EvWithCand->Fill(countTreta1,countTreta05);
928 fHistNtrEta03vsNtrEta1EvWithCand->Fill(countTreta1,countTreta03);
929 fHistNtrEtaV0AvsNtrEta1EvWithCand->Fill(countTreta1,vzeroMultA);
930 fHistNtrEtaV0MvsNtrEta1EvWithCand->Fill(countTreta1,vzeroMult);
931 fHistNtrCorrEta1vsNtrRawEta1EvWithCand->Fill(countTreta1,countTreta1corr);
934 if(nSelectedInMassPeak>0) {
935 fHistNtrUnCorrEvWithD->Fill(countMult,nchWeight);
936 fHistNtrCorrEvWithD->Fill(countCorr,nchWeight);
938 fHistNtrEta16vsNtrEta1EvWithD->Fill(countTreta1,countTreta16);
939 fHistNtrEta05vsNtrEta1EvWithD->Fill(countTreta1,countTreta05);
940 fHistNtrEta03vsNtrEta1EvWithD->Fill(countTreta1,countTreta03);
941 fHistNtrEtaV0AvsNtrEta1EvWithD->Fill(countTreta1,vzeroMultA);
942 fHistNtrEtaV0MvsNtrEta1EvWithD->Fill(countTreta1,vzeroMult);
943 fHistNtrCorrEta1vsNtrRawEta1EvWithD->Fill(countTreta1,countTreta1corr);
948 PostData(2,fListCuts);
949 PostData(3,fOutputCounters);
954 //________________________________________________________________________
955 void AliAnalysisTaskSEDvsMultiplicity::CreateImpactParameterHistos(){
956 // Histos for impact paramter study
957 // mass . pt , impact parameter , decay length , multiplicity
959 Int_t nbins[5]={fNMassBins,200,fNImpParBins,50,100};
960 Double_t xmin[5]={fLowmasslimit,0.,fLowerImpPar,0.,0.};
961 Double_t xmax[5]={fUpmasslimit,20.,fHigherImpPar,1.,100.};
963 fHistMassPtImpPar[0]=new THnSparseF("hMassPtImpParAll",
964 "Mass vs. pt vs.imppar - All",
966 fHistMassPtImpPar[1]=new THnSparseF("hMassPtImpParPrompt",
967 "Mass vs. pt vs.imppar - promptD",
969 fHistMassPtImpPar[2]=new THnSparseF("hMassPtImpParBfeed",
970 "Mass vs. pt vs.imppar - DfromB",
972 fHistMassPtImpPar[3]=new THnSparseF("hMassPtImpParTrueBfeed",
973 "Mass vs. pt vs.true imppar -DfromB",
975 fHistMassPtImpPar[4]=new THnSparseF("hMassPtImpParBkg",
976 "Mass vs. pt vs.imppar - backgr.",
978 for(Int_t i=0; i<5;i++){
979 fOutput->Add(fHistMassPtImpPar[i]);
983 //________________________________________________________________________
984 void AliAnalysisTaskSEDvsMultiplicity::Terminate(Option_t */*option*/)
986 // Terminate analysis
988 if(fDebug > 1) printf("AnalysisTaskSEDvsMultiplicity: Terminate() \n");
990 fOutput = dynamic_cast<TList*> (GetOutputData(1));
992 printf("ERROR: fOutput not available\n");
996 fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
998 printf("ERROR: fHistNEvents not available\n");
1001 printf("Number of Analyzed Events = %d\n",(Int_t)fHistNEvents->GetBinContent(3));
1005 //_________________________________________________________________________________________________
1006 Int_t AliAnalysisTaskSEDvsMultiplicity::CheckOrigin(TClonesArray* arrayMC, AliAODMCParticle *mcPartCandidate) const {
1008 // checking whether the mother of the particles come from a charm or a bottom quark
1011 Int_t pdgGranma = 0;
1013 mother = mcPartCandidate->GetMother();
1015 Int_t abspdgGranma =0;
1016 Bool_t isFromB=kFALSE;
1017 // Bool_t isQuarkFound=kFALSE;
1020 AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
1022 pdgGranma = mcGranma->GetPdgCode();
1023 abspdgGranma = TMath::Abs(pdgGranma);
1024 if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
1027 // if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
1028 mother = mcGranma->GetMother();
1030 AliError("Failed casting the mother particle!");
1035 if(isFromB) return 5;
1041 //____________________________________________________________________________
1042 TProfile* AliAnalysisTaskSEDvsMultiplicity::GetEstimatorHistogram(const AliVEvent* event){
1043 // Get Estimator Histogram from period event->GetRunNumber();
1045 // If you select SPD tracklets in |eta|<1 you should use type == 1
1048 Int_t runNo = event->GetRunNumber();
1049 Int_t period = -1; // pp: 0-LHC10b, 1-LHC10c, 2-LHC10d, 3-LHC10e
1050 // pPb: 0-LHC13b, 1-LHC13c
1052 if (runNo>195343 && runNo<195484) period = 0;
1053 if (runNo>195528 && runNo<195678) period = 1;
1054 if (period < 0 || period > 1) return 0;
1057 if(runNo>114930 && runNo<117223) period = 0;
1058 if(runNo>119158 && runNo<120830) period = 1;
1059 if(runNo>122373 && runNo<126438) period = 2;
1060 if(runNo>127711 && runNo<130841) period = 3;
1061 if(period<0 || period>3) return 0;
1065 return fMultEstimatorAvg[period];
1068 //__________________________________________________________________________________________________
1069 void AliAnalysisTaskSEDvsMultiplicity::CreateMeasuredNchHisto(){
1070 // creates historgam with measured multiplcity distribution in pp 7 TeV collisions (from Eur. Phys. J. C (2010) 68: 345–354)
1072 // for Nch > 70 the points were obtainedwith a double NBD distribution
1073 // 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
1074 // fit1->SetParameter(1,1.63); // k parameter
1075 // fit1->SetParameter(2,12.8); // mean multiplicity
1076 Double_t nchbins[82]={0.50,1.50,2.50,3.50,4.50,5.50,6.50,7.50,8.50,9.50,
1077 10.50,11.50,12.50,13.50,14.50,15.50,16.50,17.50,18.50,19.50,
1078 20.50,21.50,22.50,23.50,24.50,25.50,26.50,27.50,28.50,29.50,
1079 30.50,31.50,32.50,33.50,34.50,35.50,36.50,37.50,38.50,39.50,
1080 40.50,41.50,42.50,43.50,44.50,45.50,46.50,47.50,48.50,49.50,
1081 50.50,51.50,52.50,53.50,54.50,55.50,56.50,57.50,58.50,59.50,
1082 60.50,62.50,64.50,66.50,68.50,70.50,72.50,74.50,76.50,78.50,
1083 80.50,82.50,84.50,86.50,88.50,90.50,92.50,94.50,96.50,98.50,
1085 Double_t pch[81]={0.062011,0.072943,0.070771,0.067245,0.062834,0.057383,0.051499,0.04591,0.041109,0.036954,
1086 0.03359,0.030729,0.028539,0.026575,0.024653,0.0229,0.021325,0.019768,0.018561,0.017187,
1087 0.01604,0.014836,0.013726,0.012576,0.011481,0.010393,0.009502,0.008776,0.008024,0.007452,
1088 0.006851,0.006428,0.00594,0.005515,0.005102,0.00469,0.004162,0.003811,0.003389,0.003071,
1089 0.002708,0.002422,0.002184,0.001968,0.00186,0.00165,0.001577,0.001387,0.001254,0.001118,
1090 0.001037,0.000942,0.000823,0.000736,0.000654,0.000579,0.000512,0.00049,0.00045,0.000355,
1091 0.000296,0.000265,0.000193,0.00016,0.000126,0.0000851, 0.0000676,0.0000537,0.0000426, 0.0000338,
1092 0.0000268,0.0000213,0.0000166,0.0000133,0.0000106,0.00000837,0.00000662, 0.00000524,0.00000414, 0.00000327,
1095 if(fHistoMeasNch) delete fHistoMeasNch;
1096 fHistoMeasNch=new TH1F("hMeaseNch","",81,nchbins);
1097 for(Int_t i=0; i<81; i++){
1098 fHistoMeasNch->SetBinContent(i+1,pch[i]);
1099 fHistoMeasNch->SetBinError(i+1,0.);
1103 //__________________________________________________________________________________________________
1104 void AliAnalysisTaskSEDvsMultiplicity::FillMCMassHistos(TClonesArray *arrayMC, Int_t labD, Int_t countMult,Double_t nchWeight)
1107 // Function to fill the true MC signal
1111 AliAODMCParticle *partD = (AliAODMCParticle*)arrayMC->At(labD);
1112 Double_t mass = partD->M();
1113 Double_t pt = partD->Pt();
1114 fPtVsMassVsMultMC->Fill(countMult,mass,pt,nchWeight);