]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGHF/vertexingHF/AliAnalysisTaskSEDvsMultiplicity.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Merge branch 'master' into TRDdev
[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / AliAnalysisTaskSEDvsMultiplicity.cxx
CommitLineData
43093cfa 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>
19aaafff 33#include <TProfile.h>
43093cfa 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"
3cca4c06 50#include "AliAODVZERO.h"
43093cfa 51ClassImp(AliAnalysisTaskSEDvsMultiplicity)
52
53
54//________________________________________________________________________
55AliAnalysisTaskSEDvsMultiplicity::AliAnalysisTaskSEDvsMultiplicity():
56AliAnalysisTaskSE(),
57 fOutput(0),
58 fListCuts(0),
59 fOutputCounters(0),
6c06d0fe 60 fListProfiles(0),
43093cfa 61 fHistNEvents(0),
c1a78ad6 62 fHistNtrEta16vsNtrEta1(0),
cd465657 63 fHistNtrEta05vsNtrEta1(0),
64 fHistNtrEta03vsNtrEta1(0),
0429d713 65 fHistNtrEtaV0AvsNtrEta1(0),
66 fHistNtrEtaV0MvsNtrEta1(0),
c1a78ad6 67 fHistNtrCorrEta1vsNtrRawEta1(0),
68 fHistNtrVsZvtx(0),
69 fHistNtrCorrVsZvtx(0),
11457f09 70 fHistNtrVsNchMC(0),
71 fHistNtrCorrVsNchMC(0),
72 fHistNtrVsNchMCPrimary(0),
73 fHistNtrCorrVsNchMCPrimary(0),
74 fHistNtrVsNchMCPhysicalPrimary(0),
75 fHistNtrCorrVsNchMCPhysicalPrimary(0),
0c514eec 76 fHistGenPrimaryParticlesInelGt0(0),
511c4e3e 77 fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary(0),
f73a90da 78 fHistNtrUnCorrEvSel(0),
bba1a026 79 fHistNtrUnCorrEvWithCand(0),
80 fHistNtrUnCorrEvWithD(0),
c1a78ad6 81 fHistNtrCorrEvSel(0),
82 fHistNtrCorrEvWithCand(0),
83 fHistNtrCorrEvWithD(0),
43093cfa 84 fPtVsMassVsMult(0),
85 fPtVsMassVsMultNoPid(0),
86 fPtVsMassVsMultUncorr(0),
87 fPtVsMassVsMultPart(0),
88 fPtVsMassVsMultAntiPart(0),
98fa4fcd 89 fPtVsMassVsMultMC(0),
43093cfa 90 fUpmasslimit(1.965),
91 fLowmasslimit(1.765),
7826c36d 92 fNMassBins(200),
43093cfa 93 fRDCutsAnalysis(0),
94 fCounter(0),
95 fCounterU(0),
96 fDoImpPar(kFALSE),
97 fNImpParBins(400),
98 fLowerImpPar(-2000.),
99 fHigherImpPar(2000.),
100 fReadMC(kFALSE),
101 fMCOption(0),
bba1a026 102 fisPPbData(kFALSE),
43093cfa 103 fUseBit(kTRUE),
0d1e32a6 104 fSubtractTrackletsFromDau(kFALSE),
0429d713 105 fKeepCorrPlots(kFALSE),
3cca4c06 106 fUseNchWeight(kFALSE),
107 fHistoMCNch(0),
108 fHistoMeasNch(0),
f73a90da 109 fRefMult(9.26),
3cca4c06 110 fPdgMeson(411),
20550d50 111 fMultiplicityEstimator(kNtrk10),
112 fMCPrimariesEstimator(kEta10)
43093cfa 113{
114 // Default constructor
115 for(Int_t i=0; i<5; i++) fHistMassPtImpPar[i]=0;
116 for(Int_t i=0; i<4; i++) fMultEstimatorAvg[i]=0;
117}
118
119//________________________________________________________________________
b97c5dc0 120AliAnalysisTaskSEDvsMultiplicity::AliAnalysisTaskSEDvsMultiplicity(const char *name, Int_t pdgMeson,AliRDHFCuts *cuts, Bool_t switchPPb):
43093cfa 121 AliAnalysisTaskSE(name),
122 fOutput(0),
123 fListCuts(0),
124 fOutputCounters(0),
6c06d0fe 125 fListProfiles(0),
43093cfa 126 fHistNEvents(0),
c1a78ad6 127 fHistNtrEta16vsNtrEta1(0),
cd465657 128 fHistNtrEta05vsNtrEta1(0),
129 fHistNtrEta03vsNtrEta1(0),
0429d713 130 fHistNtrEtaV0AvsNtrEta1(0),
131 fHistNtrEtaV0MvsNtrEta1(0),
c1a78ad6 132 fHistNtrCorrEta1vsNtrRawEta1(0),
133 fHistNtrVsZvtx(0),
134 fHistNtrCorrVsZvtx(0),
11457f09 135 fHistNtrVsNchMC(0),
136 fHistNtrCorrVsNchMC(0),
137 fHistNtrVsNchMCPrimary(0),
138 fHistNtrCorrVsNchMCPrimary(0),
139 fHistNtrVsNchMCPhysicalPrimary(0),
140 fHistNtrCorrVsNchMCPhysicalPrimary(0),
0c514eec 141 fHistGenPrimaryParticlesInelGt0(0),
511c4e3e 142 fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary(0),
f73a90da 143 fHistNtrUnCorrEvSel(0),
bba1a026 144 fHistNtrUnCorrEvWithCand(0),
145 fHistNtrUnCorrEvWithD(0),
c1a78ad6 146 fHistNtrCorrEvSel(0),
147 fHistNtrCorrEvWithCand(0),
148 fHistNtrCorrEvWithD(0),
43093cfa 149 fPtVsMassVsMult(0),
150 fPtVsMassVsMultNoPid(0),
151 fPtVsMassVsMultUncorr(0),
152 fPtVsMassVsMultPart(0),
153 fPtVsMassVsMultAntiPart(0),
98fa4fcd 154 fPtVsMassVsMultMC(0),
43093cfa 155 fUpmasslimit(1.965),
156 fLowmasslimit(1.765),
7826c36d 157 fNMassBins(200),
43093cfa 158 fRDCutsAnalysis(cuts),
159 fCounter(0),
160 fCounterU(0),
161 fDoImpPar(kFALSE),
162 fNImpParBins(400),
163 fLowerImpPar(-2000.),
164 fHigherImpPar(2000.),
165 fReadMC(kFALSE),
166 fMCOption(0),
b97c5dc0 167 fisPPbData(switchPPb),
43093cfa 168 fUseBit(kTRUE),
0d1e32a6 169 fSubtractTrackletsFromDau(kFALSE),
0429d713 170 fKeepCorrPlots(kFALSE),
3cca4c06 171 fUseNchWeight(kFALSE),
172 fHistoMCNch(0),
173 fHistoMeasNch(0),
f73a90da 174 fRefMult(9.26),
3cca4c06 175 fPdgMeson(pdgMeson),
20550d50 176 fMultiplicityEstimator(kNtrk10),
177 fMCPrimariesEstimator(kEta10)
43093cfa 178{
179 //
180 // Standard constructor
181 //
b97c5dc0 182
43093cfa 183 for(Int_t i=0; i<5; i++) fHistMassPtImpPar[i]=0;
184 for(Int_t i=0; i<4; i++) fMultEstimatorAvg[i]=0;
7826c36d 185 if(fPdgMeson==413){
8ce9400f 186 fNMassBins=200;
187 SetMassLimits(0.12,0.2);
188 }else{
189 fNMassBins=200;
7826c36d 190 SetMassLimits(fPdgMeson,0.1);
191 }
43093cfa 192 // Default constructor
3cca4c06 193 // Otput slot #1 writes into a TList container
43093cfa 194 DefineOutput(1,TList::Class()); //My private output
195 // Output slot #2 writes cut to private output
196 DefineOutput(2,TList::Class());
197 // Output slot #3 writes cut to private output
7826c36d 198 DefineOutput(3,TList::Class());
6c06d0fe 199 // Output slot #4 writes cut to private output
200 DefineOutput(4,TList::Class());
7826c36d 201}
43093cfa 202//________________________________________________________________________
203AliAnalysisTaskSEDvsMultiplicity::~AliAnalysisTaskSEDvsMultiplicity()
204{
205 //
206 // Destructor
207 //
208 delete fOutput;
209 delete fHistNEvents;
210 delete fListCuts;
6c06d0fe 211 delete fListProfiles;
43093cfa 212 delete fRDCutsAnalysis;
213 delete fCounter;
214 delete fCounterU;
b97c5dc0 215 for(Int_t i=0; i<4; i++) {
216 if (fMultEstimatorAvg[i]) delete fMultEstimatorAvg[i];
217 }
218
43093cfa 219 for(Int_t i=0; i<5; i++){
220 delete fHistMassPtImpPar[i];
221 }
3cca4c06 222 if(fHistoMCNch) delete fHistoMCNch;
223 if(fHistoMeasNch) delete fHistoMeasNch;
224}
43093cfa 225
226//_________________________________________________________________
7826c36d 227void AliAnalysisTaskSEDvsMultiplicity::SetMassLimits(Double_t lowlimit, Double_t uplimit){
43093cfa 228 // set invariant mass limits
229 if(uplimit>lowlimit){
a8f9e77f 230 fLowmasslimit = lowlimit;
231 fUpmasslimit = uplimit;
7826c36d 232 }else{
233 AliError("Wrong mass limits: upper value should be larger than lower one");
43093cfa 234 }
235}
7826c36d 236//_________________________________________________________________
237void AliAnalysisTaskSEDvsMultiplicity::SetMassLimits(Int_t pdg, Double_t range){
238 // set invariant mass limits
239 Double_t mass=TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdg))->Mass();
240 SetMassLimits(mass-range,mass+range);
43093cfa 241}
242//________________________________________________________________________
243void AliAnalysisTaskSEDvsMultiplicity::Init(){
244 //
245 // Initialization
246 //
247 printf("AnalysisTaskSEDvsMultiplicity::Init() \n");
3cca4c06 248
249 if(fUseNchWeight && !fReadMC){ AliFatal("Nch weights can only be used in MC mode"); return; }
250 if(fUseNchWeight && !fHistoMCNch){ AliFatal("Nch weights can only be used without histogram"); return; }
43093cfa 251
252 fListCuts=new TList();
bce4056e 253 fListCuts->SetOwner();
254 fListCuts->SetName("CutsList");
255
43093cfa 256
257 if(fPdgMeson==411){
258 AliRDHFCutsDplustoKpipi* copycut=new AliRDHFCutsDplustoKpipi(*(static_cast<AliRDHFCutsDplustoKpipi*>(fRDCutsAnalysis)));
259 copycut->SetName("AnalysisCutsDplus");
260 fListCuts->Add(copycut);
261 }else if(fPdgMeson==421){
262 AliRDHFCutsD0toKpi* copycut=new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fRDCutsAnalysis)));
263 copycut->SetName("AnalysisCutsDzero");
264 fListCuts->Add(copycut);
265 }else if(fPdgMeson==413){
266 AliRDHFCutsDStartoKpipi* copycut=new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fRDCutsAnalysis)));
267 copycut->SetName("AnalysisCutsDStar");
268 fListCuts->Add(copycut);
269 }
270 PostData(2,fListCuts);
271
6c06d0fe 272 fListProfiles = new TList();
273 fListProfiles->SetOwner();
b97c5dc0 274 TString period[4];
3d87cb41 275 Int_t nProfiles=4;
276 if (fisPPbData) {period[0]="LHC13b"; period[1]="LHC13c"; nProfiles = 2;}
277 else {period[0]="LHC10b"; period[1]="LHC10c"; period[2]="LHC10d"; period[3]="LHC10e"; nProfiles = 4;}
b97c5dc0 278
3d87cb41 279 for(Int_t i=0; i<nProfiles; i++){
6c06d0fe 280 if(fMultEstimatorAvg[i]){
281 TProfile* hprof=new TProfile(*fMultEstimatorAvg[i]);
282 hprof->SetName(Form("ProfileTrkVsZvtx%s\n",period[i].Data()));
283 fListProfiles->Add(hprof);
284 }
285 }
286 PostData(4,fListProfiles);
287
43093cfa 288 return;
289}
290
291//________________________________________________________________________
292void AliAnalysisTaskSEDvsMultiplicity::UserCreateOutputObjects()
293{
294 // Create the output container
295 //
296 if(fDebug > 1) printf("AnalysisTaskSEDvsMultiplicity::UserCreateOutputObjects() \n");
297
298 // Several histograms are more conveniently managed in a TList
299 fOutput = new TList();
300 fOutput->SetOwner();
301 fOutput->SetName("OutputHistos");
302
3cca4c06 303 Int_t nMultBins = 200;
304 Float_t firstMultBin = -0.5;
305 Float_t lastMultBin = 199.5;
bba1a026 306 if(fisPPbData) {
307 nMultBins = 375;
308 lastMultBin = 374.5;
309 }
5da5a290 310 if(fMultiplicityEstimator==kVZERO || fMultiplicityEstimator==kVZEROA) {
6807ad9a 311 nMultBins = 400;
312 lastMultBin = 799.5;
3cca4c06 313 }
314
315 fHistNtrUnCorrEvSel = new TH1F("hNtrUnCorrEvSel","Uncorrected tracklets multiplicity for selected events; Tracklets ; Entries",nMultBins,firstMultBin,lastMultBin);
bba1a026 316 fHistNtrUnCorrEvWithCand = new TH1F("hNtrUnCorrEvWithCand", "Uncorrected Tracklets multiplicity for events with D candidates; Tracklets ; Entries",nMultBins,firstMultBin,lastMultBin);// Total multiplicity
317 fHistNtrUnCorrEvWithD = new TH1F("hNtrUnCorrEvWithD","Uncorrected Tracklets multiplicity for events with D in mass region ; Tracklets ; Entries",nMultBins,firstMultBin,lastMultBin); //
3cca4c06 318 fHistNtrCorrEvSel = new TH1F("hNtrCorrEvSel","Corrected tracklets multiplicity for selected events; Tracklets ; Entries",nMultBins,firstMultBin,lastMultBin);
319 fHistNtrCorrEvWithCand = new TH1F("hNtrCorrEvWithCand", "Tracklets multiplicity for events with D candidates; Tracklets ; Entries",nMultBins,firstMultBin,lastMultBin);// Total multiplicity
320 fHistNtrCorrEvWithD = new TH1F("hNtrCorrEvWithD", "Tracklets multiplicity for events with D in mass region ; Tracklets ; Entries",nMultBins,firstMultBin,lastMultBin); //
0429d713 321
322 if(fKeepCorrPlots){
323 fHistNtrEta16vsNtrEta1 = new TH2F("hNtrEta16vsNtrEta1","Uncorrected Eta1.6 vs Eta1.0; Ntracklets #eta<1.0; Ntracklets #eta<1.6",nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); //eta 1.6 vs eta 1.0 histogram
324 fHistNtrEta05vsNtrEta1 = new TH2F("hNtrEta05vsNtrEta1","Uncorrected Eta0.5 vs Eta1.0; Ntracklets #eta<1.0; Ntracklets #eta<0.5",nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); //eta 0.5 vs eta 1.0 histogram
325 fHistNtrEta03vsNtrEta1 = new TH2F("hNtrEta03vsNtrEta1","Uncorrected Eta0.3 vs Eta1.0; Ntracklets #eta<1.0; Ntracklets #eta<0.3",nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); //eta 0.3 vs eta 1.0 histogram
81cf4d8c 326 fHistNtrEtaV0AvsNtrEta1 = new TH2F("hNtrEtaV0AvsNtrEta1","Uncorrected Eta-V0A vs Eta1.0; Ntracklets #eta<1.0; Multiplicity V0A",nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); //eta V0A vs eta 1.0 histogram
327 fHistNtrEtaV0MvsNtrEta1 = new TH2F("hNtrEtaV0MvsNtrEta1","Uncorrected Eta-V0M vs Eta1.0; Ntracklets #eta<1.0; Multiplicity V0A+V0C",nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); //eta V0M vs eta 1.0 histogram
0429d713 328 fHistNtrCorrEta1vsNtrRawEta1 = new TH2F("hNtrCorrEta1vsNtrRawEta1","Corrected Eta1 vs Eta1.0; Ntracklets #eta<1.0 corrected; Ntracklets #eta<1",nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); //eta 1.6 vs eta 1.0 histogram
329 }
3cca4c06 330 fHistNtrVsZvtx = new TH2F("hNtrVsZvtx","Ntracklet vs VtxZ; VtxZ;N_{tracklet};",300,-15,15,nMultBins,firstMultBin,lastMultBin); //
331 fHistNtrCorrVsZvtx = new TH2F("hNtrCorrVsZvtx","Ntracklet vs VtxZ; VtxZ;N_{tracklet};",300,-15,15,nMultBins,firstMultBin,lastMultBin); //
332
333 fHistNtrVsNchMC = new TH2F("hNtrVsNchMC","Ntracklet vs NchMC; Nch;N_{tracklet};",nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); //
334 fHistNtrCorrVsNchMC = new TH2F("hNtrCorrVsNchMC","Ntracklet vs Nch; Nch;N_{tracklet};",nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); //
11457f09 335
3cca4c06 336 fHistNtrVsNchMCPrimary = new TH2F("hNtrVsNchMCPrimary","Ntracklet vs Nch (Primary); Nch (Primary);N_{tracklet};",nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); //
337 fHistNtrCorrVsNchMCPrimary = new TH2F("hNtrCorrVsNchMCPrimary","Ntracklet vs Nch (Primary); Nch(Primary) ;N_{tracklet};",nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); //
11457f09 338
3cca4c06 339 fHistNtrVsNchMCPhysicalPrimary = new TH2F("hNtrVsNchMCPhysicalPrimary","Ntracklet vs Nch (Physical Primary); Nch (Physical Primary);N_{tracklet};",nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); //
340 fHistNtrCorrVsNchMCPhysicalPrimary = new TH2F("hNtrCorrVsMCPhysicalPrimary","Ntracklet vs Nch (Physical Primary); Nch (Physical Primary);N_{tracklet};",nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); //
43093cfa 341
3cca4c06 342 fHistGenPrimaryParticlesInelGt0 = new TH1F("hGenPrimaryParticlesInelGt0","Multiplcity of generated charged particles ; Nparticles ; Entries",nMultBins,firstMultBin,lastMultBin);
43093cfa 343
3cca4c06 344 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);
511c4e3e 345
2813d7bc 346 fHistNtrUnCorrEvSel->Sumw2();
bba1a026 347 fHistNtrUnCorrEvWithCand->Sumw2();
348 fHistNtrUnCorrEvWithD->Sumw2();
c1a78ad6 349 fHistNtrCorrEvSel->Sumw2();
350 fHistNtrCorrEvWithCand->Sumw2();
351 fHistNtrCorrEvWithD->Sumw2();
0c514eec 352 fHistGenPrimaryParticlesInelGt0->Sumw2();
2813d7bc 353 fOutput->Add(fHistNtrUnCorrEvSel);
bba1a026 354 fOutput->Add(fHistNtrUnCorrEvWithCand);
355 fOutput->Add(fHistNtrUnCorrEvWithD);
c1a78ad6 356 fOutput->Add(fHistNtrCorrEvSel);
357 fOutput->Add(fHistNtrCorrEvWithCand);
358 fOutput->Add(fHistNtrCorrEvWithD);
0429d713 359 if(fKeepCorrPlots){
360 fOutput->Add(fHistNtrEta16vsNtrEta1);
361 fOutput->Add(fHistNtrEta05vsNtrEta1);
362 fOutput->Add(fHistNtrEta03vsNtrEta1);
363 fOutput->Add(fHistNtrEtaV0AvsNtrEta1);
364 fOutput->Add(fHistNtrEtaV0MvsNtrEta1);
365 fOutput->Add(fHistNtrCorrEta1vsNtrRawEta1);
366 }
c1a78ad6 367 fOutput->Add(fHistNtrVsZvtx);
368 fOutput->Add(fHistNtrCorrVsZvtx);
43093cfa 369
11457f09 370 fOutput->Add(fHistNtrVsNchMC);
371 fOutput->Add(fHistNtrCorrVsNchMC);
372 fOutput->Add(fHistNtrVsNchMCPrimary);
373 fOutput->Add(fHistNtrCorrVsNchMCPrimary);
374 fOutput->Add(fHistNtrVsNchMCPhysicalPrimary);
375 fOutput->Add(fHistNtrCorrVsNchMCPhysicalPrimary);
0c514eec 376 fOutput->Add(fHistGenPrimaryParticlesInelGt0);
511c4e3e 377 fOutput->Add(fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary);
11457f09 378
43093cfa 379
380 fHistNEvents = new TH1F("fHistNEvents", "number of events ",11,-0.5,10.5);
381 fHistNEvents->GetXaxis()->SetBinLabel(1,"nEvents total");
382 fHistNEvents->GetXaxis()->SetBinLabel(2,"nEvents with Z vertex");
383 fHistNEvents->GetXaxis()->SetBinLabel(3,"nEvents selected");
384 fHistNEvents->GetXaxis()->SetBinLabel(4,"Rejected due to trigger");
385 fHistNEvents->GetXaxis()->SetBinLabel(5,"Rejected due to phys sel");
386 fHistNEvents->GetXaxis()->SetBinLabel(6,"Rejected due to vertex cuts");
387 fHistNEvents->GetXaxis()->SetBinLabel(7,"Rejected due to pileup");
388 fHistNEvents->GetXaxis()->SetBinLabel(8,"Total no. of candidate");
389 fHistNEvents->GetXaxis()->SetBinLabel(9,"no. of cand wo bitmask");
390 fHistNEvents->GetXaxis()->SetBinLabel(10,"D after cuts (No PID)");
391 fHistNEvents->GetXaxis()->SetBinLabel(11,"D after cuts + PID)");
392 fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);
393 fHistNEvents->Sumw2();
394 fHistNEvents->SetMinimum(0);
395 fOutput->Add(fHistNEvents);
396
3cca4c06 397 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.);
43093cfa 398
3cca4c06 399 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.);
43093cfa 400
3cca4c06 401 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.);
43093cfa 402
3cca4c06 403 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.);
43093cfa 404
3cca4c06 405 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.);
43093cfa 406
20550d50 407 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.);
408
43093cfa 409 fOutput->Add(fPtVsMassVsMult);
410 fOutput->Add(fPtVsMassVsMultUncorr);
411 fOutput->Add(fPtVsMassVsMultNoPid);
412 fOutput->Add(fPtVsMassVsMultPart);
413 fOutput->Add(fPtVsMassVsMultAntiPart);
20550d50 414 fOutput->Add(fPtVsMassVsMultMC);
43093cfa 415
416 if(fDoImpPar) CreateImpactParameterHistos();
417
418 fCounter = new AliNormalizationCounter("NormCounterCorrMult");
419 fCounter->SetStudyMultiplicity(kTRUE,1.);
420 fCounter->Init();
421
422 fCounterU = new AliNormalizationCounter("NormCounterUnCorrMult");
423 fCounterU->SetStudyMultiplicity(kTRUE,1.);
424 fCounterU->Init();
425
426 fOutputCounters = new TList();
427 fOutputCounters->SetOwner();
428 fOutputCounters->SetName("OutputCounters");
429 fOutputCounters->Add(fCounter);
430 fOutputCounters->Add(fCounterU);
431
432 PostData(1,fOutput);
433 PostData(2,fListCuts);
434 PostData(3,fOutputCounters);
6c06d0fe 435 PostData(4,fListProfiles);
436
3cca4c06 437 if(fUseNchWeight) CreateMeasuredNchHisto();
43093cfa 438
439 return;
440}
441
442
443
444//________________________________________________________________________
445void AliAnalysisTaskSEDvsMultiplicity::UserExec(Option_t */*option*/)
446{
447 // Execute analysis for current event:
448 // heavy flavor candidates association to MC truth
449
43093cfa 450 AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
451
452 // AliAODTracklets* tracklets = aod->GetTracklets();
453 //Int_t ntracklets = tracklets->GetNumberOfTracklets();
454
455
456 TClonesArray *arrayCand = 0;
457 TString arrayName="";
458 UInt_t pdgDau[3];
459 Int_t nDau=0;
460 Int_t selbit=0;
461 if(fPdgMeson==411){
462 arrayName="Charm3Prong";
463 pdgDau[0]=211; pdgDau[1]=321; pdgDau[2]=211;
464 nDau=3;
465 selbit=AliRDHFCuts::kDplusCuts;
466 }else if(fPdgMeson==421){
467 arrayName="D0toKpi";
468 pdgDau[0]=211; pdgDau[1]=321; pdgDau[2]=0;
469 nDau=2;
470 selbit=AliRDHFCuts::kD0toKpiCuts;
471 }else if(fPdgMeson==413){
472 arrayName="Dstar";
8ce9400f 473 pdgDau[0]=321; pdgDau[1]=211; pdgDau[2]=0; // Quoting here D0 daughters (D* ones on another variable later)
474 nDau=2;
43093cfa 475 selbit=AliRDHFCuts::kDstarCuts;
476 }
6c8b22ea 477
43093cfa 478 if(!aod && AODEvent() && IsStandardAOD()) {
479 // In case there is an AOD handler writing a standard AOD, use the AOD
480 // event in memory rather than the input (ESD) event.
481 aod = dynamic_cast<AliAODEvent*> (AODEvent());
482 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
483 // have to taken from the AOD event hold by the AliAODExtension
484 AliAODHandler* aodHandler = (AliAODHandler*)
485 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
486 if(aodHandler->GetExtensions()) {
487 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
488 AliAODEvent *aodFromExt = ext->GetAOD();
489 arrayCand=(TClonesArray*)aodFromExt->GetList()->FindObject(arrayName.Data());
490 }
491 } else if(aod) {
492 arrayCand=(TClonesArray*)aod->GetList()->FindObject(arrayName.Data());
493 }
43093cfa 494
495 if(!aod || !arrayCand) {
496 printf("AliAnalysisTaskSEDvsMultiplicity::UserExec: Charm3Prong branch not found!\n");
497 return;
498 }
499
bba1a026 500 if(fisPPbData && fReadMC){
501 Int_t runnumber = aod->GetRunNumber();
502 if(aod->GetTriggerMask()==0 &&
503 (runnumber>=195344 && runnumber<=195677)){
504 AliDebug(3,"Event rejected because of null trigger mask");
505 return;
506 }
507 }
43093cfa 508
509
510 // fix for temporary bug in ESDfilter
511 // the AODs with null vertex pointer didn't pass the PhysSel
512 if(!aod->GetPrimaryVertex()||TMath::Abs(aod->GetMagneticField())<0.001) return;
513
514 Int_t countTreta1=AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aod,-1.,1.);
cd465657 515 Int_t countTreta03=AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aod,-0.3,0.3);
516 Int_t countTreta05=AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aod,-0.5,0.5);
517 Int_t countTreta16=AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aod,-1.6,1.6);
43093cfa 518
3cca4c06 519 Int_t vzeroMult=0;
cd465657 520 Int_t vzeroMultA=0;
3cca4c06 521 AliAODVZERO *vzeroAOD = (AliAODVZERO*)aod->GetVZEROData();
cd465657 522 if(vzeroAOD) {
523 vzeroMult = vzeroAOD->GetMTotV0A() + vzeroAOD->GetMTotV0C();
524 vzeroMultA = vzeroAOD->GetMTotV0A();
525 }
3cca4c06 526
527 Int_t countMult = countTreta1;
cd465657 528 if(fMultiplicityEstimator==kNtrk03) { countMult = countTreta03; }
529 if(fMultiplicityEstimator==kNtrk05) { countMult = countTreta05; }
530 if(fMultiplicityEstimator==kNtrk10to16) { countMult = countTreta16 - countTreta1; }
3cca4c06 531 if(fMultiplicityEstimator==kVZERO) { countMult = vzeroMult; }
cd465657 532 if(fMultiplicityEstimator==kVZEROA) { countMult = vzeroMultA; }
3cca4c06 533
534
535 fCounterU->StoreEvent(aod,fRDCutsAnalysis,fReadMC,countMult);
43093cfa 536 fHistNEvents->Fill(0); // count event
537
43093cfa 538 Double_t countTreta1corr=countTreta1;
3cca4c06 539 Double_t countCorr=countMult;
f1c9cf1c 540 AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
3cca4c06 541 // if(vtx1){
542 // FIX ME: No correction to the VZERO !!
0429d713 543 if(vtx1 && (fMultiplicityEstimator!=kVZERO) && (fMultiplicityEstimator!=kVZEROA)){
544 if(vtx1->GetNContributors()>0){
545 fHistNEvents->Fill(1);
f1c9cf1c 546 TProfile* estimatorAvg = GetEstimatorHistogram(aod);
547 if(estimatorAvg){
0429d713 548 countTreta1corr=AliVertexingHFUtils::GetCorrectedNtracklets(estimatorAvg,countTreta1,vtx1->GetZ(),fRefMult);
3cca4c06 549 countCorr=AliVertexingHFUtils::GetCorrectedNtracklets(estimatorAvg,countMult,vtx1->GetZ(),fRefMult);
f1c9cf1c 550 }
43093cfa 551 }
552 }
0429d713 553
43093cfa 554
555 Bool_t isEvSel=fRDCutsAnalysis->IsEventSelected(aod);
556
557 if(fRDCutsAnalysis->GetWhyRejection()==5) fHistNEvents->Fill(3);
558 if(fRDCutsAnalysis->GetWhyRejection()==7) fHistNEvents->Fill(4);
c1a78ad6 559 if(fRDCutsAnalysis->GetWhyRejection()==6) fHistNEvents->Fill(5);
560 if(fRDCutsAnalysis->GetWhyRejection()==1) fHistNEvents->Fill(6);
43093cfa 561
43093cfa 562
c1a78ad6 563 if(!isEvSel)return;
0429d713 564 if(fKeepCorrPlots){
565 fHistNtrEta16vsNtrEta1->Fill(countTreta1,countTreta16);
566 fHistNtrEta05vsNtrEta1->Fill(countTreta1,countTreta05);
567 fHistNtrEta03vsNtrEta1->Fill(countTreta1,countTreta03);
568 fHistNtrEtaV0AvsNtrEta1->Fill(countTreta1,vzeroMultA);
569 fHistNtrEtaV0MvsNtrEta1->Fill(countTreta1,vzeroMult);
570 fHistNtrCorrEta1vsNtrRawEta1->Fill(countTreta1,countTreta1corr);
571 }
c1a78ad6 572 if(vtx1){
3cca4c06 573 fHistNtrVsZvtx->Fill(vtx1->GetZ(),countMult);
574 fHistNtrCorrVsZvtx->Fill(vtx1->GetZ(),countCorr);
c1a78ad6 575 }
576
43093cfa 577 TClonesArray *arrayMC=0;
578 AliAODMCHeader *mcHeader=0;
579
3cca4c06 580 Double_t nchWeight=1.0;
581
43093cfa 582 // load MC particles
583 if(fReadMC){
584
585 arrayMC = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
586 if(!arrayMC) {
587 printf("AliAnalysisTaskSEDvsMultiplicity::UserExec: MC particles branch not found!\n");
588 return;
589 }
590 // load MC header
591 mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
592 if(!mcHeader) {
593 printf("AliAnalysisTaskSEDvsMultiplicity::UserExec: MC header branch not found!\n");
594 return;
595 }
11457f09 596
597
598 Int_t nChargedMC=AliVertexingHFUtils::GetGeneratedMultiplicityInEtaRange(arrayMC,-1.0,1.0);
599 Int_t nChargedMCPrimary=AliVertexingHFUtils::GetGeneratedPrimariesInEtaRange(arrayMC,-1.0,1.0);
600 Int_t nChargedMCPhysicalPrimary=AliVertexingHFUtils::GetGeneratedPhysicalPrimariesInEtaRange(arrayMC,-1.0,1.0);
20550d50 601
602 // Compute the Nch weights (reference is Ntracklets within |eta|<1.0)
3cca4c06 603 if(fUseNchWeight){
604 Double_t tmpweight = 1.0;
605 if(nChargedMCPhysicalPrimary<=0) tmpweight = 0.0;
606 else{
bce4056e 607 Double_t pMeas = fHistoMeasNch->GetBinContent(fHistoMeasNch->FindBin(nChargedMCPhysicalPrimary));
3cca4c06 608 // printf(" pMeas=%2.2f and histo MCNch %s \n",pMeas,fHistoMCNch);
bce4056e 609 Double_t pMC = fHistoMCNch->GetBinContent(fHistoMCNch->FindBin(nChargedMCPhysicalPrimary));
bf27d64a 610 tmpweight = pMC>0 ? pMeas/pMC : 0.;
3cca4c06 611 }
612 nchWeight *= tmpweight;
613 AliDebug(2,Form("Using Nch weights, Mult=%d Weight=%f\n",nChargedMCPhysicalPrimary,nchWeight));
614 }
20550d50 615
616 // Now recompute the variables in case another MC estimator is considered
617 Int_t nChargedMCEta10 = nChargedMC;
618 Int_t nChargedMCEta16 = AliVertexingHFUtils::GetGeneratedMultiplicityInEtaRange(arrayMC,-1.6,1.6);
cd465657 619 Int_t nChargedMCEta05 = AliVertexingHFUtils::GetGeneratedMultiplicityInEtaRange(arrayMC,-0.5,0.5);
620 Int_t nChargedMCEta03 = AliVertexingHFUtils::GetGeneratedMultiplicityInEtaRange(arrayMC,-0.3,0.3);
20550d50 621 Int_t nChargedMCEtam37tm17 = AliVertexingHFUtils::GetGeneratedMultiplicityInEtaRange(arrayMC,-3.7,-1.7);
622 Int_t nChargedMCEta28t51 = AliVertexingHFUtils::GetGeneratedMultiplicityInEtaRange(arrayMC,2.8,5.1);
623 Int_t nChargedMCPrimaryEta10 = nChargedMCPrimary;
624 Int_t nChargedMCPrimaryEta16 = AliVertexingHFUtils::GetGeneratedPrimariesInEtaRange(arrayMC,-1.6,1.6);
cd465657 625 Int_t nChargedMCPrimaryEta05 = AliVertexingHFUtils::GetGeneratedPrimariesInEtaRange(arrayMC,-0.5,0.5);
626 Int_t nChargedMCPrimaryEta03 = AliVertexingHFUtils::GetGeneratedPrimariesInEtaRange(arrayMC,-0.3,0.3);
20550d50 627 Int_t nChargedMCPrimaryEtam37tm17 = AliVertexingHFUtils::GetGeneratedPrimariesInEtaRange(arrayMC,-3.7,-1.7);
628 Int_t nChargedMCPrimaryEta28t51 = AliVertexingHFUtils::GetGeneratedPrimariesInEtaRange(arrayMC,2.8,5.1);
629 Int_t nChargedMCPhysicalPrimaryEta10 = nChargedMCPhysicalPrimary;
630 Int_t nChargedMCPhysicalPrimaryEta16 = AliVertexingHFUtils::GetGeneratedPhysicalPrimariesInEtaRange(arrayMC,-1.6,1.6);
cd465657 631 Int_t nChargedMCPhysicalPrimaryEta05 = AliVertexingHFUtils::GetGeneratedPhysicalPrimariesInEtaRange(arrayMC,-0.5,0.5);
632 Int_t nChargedMCPhysicalPrimaryEta03 = AliVertexingHFUtils::GetGeneratedPhysicalPrimariesInEtaRange(arrayMC,-0.3,0.3);
20550d50 633 Int_t nChargedMCPhysicalPrimaryEtam37tm17 = AliVertexingHFUtils::GetGeneratedPhysicalPrimariesInEtaRange(arrayMC,-3.7,-1.7);
634 Int_t nChargedMCPhysicalPrimaryEta28t51 = AliVertexingHFUtils::GetGeneratedPhysicalPrimariesInEtaRange(arrayMC,2.8,5.1);
635 if(fMCPrimariesEstimator==kEta10to16){
636 nChargedMC = nChargedMCEta16 - nChargedMCEta10;
637 nChargedMCPrimary = nChargedMCPrimaryEta16 - nChargedMCPrimaryEta10;
638 nChargedMCPhysicalPrimary = nChargedMCPhysicalPrimaryEta16 - nChargedMCPhysicalPrimaryEta10;
cd465657 639 } else if(fMCPrimariesEstimator==kEta05){
640 nChargedMC = nChargedMCEta05;
641 nChargedMCPrimary = nChargedMCPrimaryEta05;
642 nChargedMCPhysicalPrimary = nChargedMCPhysicalPrimaryEta05;
643 } else if(fMCPrimariesEstimator==kEta03){
644 nChargedMC = nChargedMCEta03;
645 nChargedMCPrimary = nChargedMCPrimaryEta03;
646 nChargedMCPhysicalPrimary = nChargedMCPhysicalPrimaryEta03;
20550d50 647 } else if(fMCPrimariesEstimator==kEtaVZERO){
648 nChargedMC = nChargedMCEtam37tm17 + nChargedMCEta28t51;
649 nChargedMCPrimary = nChargedMCPrimaryEtam37tm17 + nChargedMCPrimaryEta28t51;
650 nChargedMCPhysicalPrimary = nChargedMCPhysicalPrimaryEtam37tm17 + nChargedMCPhysicalPrimaryEta28t51;
cd465657 651 } else if(fMCPrimariesEstimator==kEtaVZEROA){
652 nChargedMC = nChargedMCEta28t51;
653 nChargedMCPrimary = nChargedMCPrimaryEta28t51;
654 nChargedMCPhysicalPrimary = nChargedMCPhysicalPrimaryEta28t51;
20550d50 655 }
656
657 // Here fill the MC correlation plots
0c514eec 658 if(nChargedMCPhysicalPrimary>0){ // INEL>0 for |eta|<1
3cca4c06 659 fHistGenPrimaryParticlesInelGt0->Fill(nChargedMCPhysicalPrimary,nchWeight);
0c514eec 660 }
11457f09 661
3cca4c06 662 fHistNtrVsNchMC->Fill(nChargedMC,countMult,nchWeight);
663 fHistNtrCorrVsNchMC->Fill(nChargedMC,countCorr,nchWeight);
664
665 fHistNtrVsNchMCPrimary->Fill(nChargedMCPrimary,countMult,nchWeight);
666 fHistNtrCorrVsNchMCPrimary->Fill(nChargedMCPrimary,countCorr,nchWeight);
11457f09 667
3cca4c06 668 fHistNtrVsNchMCPhysicalPrimary->Fill(nChargedMCPhysicalPrimary,countMult,nchWeight);
669 fHistNtrCorrVsNchMCPhysicalPrimary->Fill(nChargedMCPhysicalPrimary,countCorr,nchWeight);
11457f09 670
3cca4c06 671 fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary->Fill(nChargedMC,nChargedMCPrimary,nChargedMCPhysicalPrimary,nchWeight);
43093cfa 672 }
673
674 Int_t nCand = arrayCand->GetEntriesFast();
d2cdcb07 675 Int_t nSelectedNoPID=0,nSelectedPID=0,nSelectedInMassPeak=0;
676 Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();
677 Double_t mDplusPDG = TDatabasePDG::Instance()->GetParticle(411)->Mass();
678 Double_t mDstarPDG = TDatabasePDG::Instance()->GetParticle(413)->Mass();
8ce9400f 679
680 // pdg of daughters needed for D* too
681 UInt_t pdgDgDStartoD0pi[2]={421,211};
682
3cca4c06 683 Double_t aveMult=0.;
684 Double_t nSelCand=0.;
43093cfa 685 for (Int_t iCand = 0; iCand < nCand; iCand++) {
686 AliAODRecoDecayHF *d = (AliAODRecoDecayHF*)arrayCand->UncheckedAt(iCand);
8ce9400f 687 AliAODRecoCascadeHF *dCascade = NULL;
688 if(fPdgMeson==413) dCascade = (AliAODRecoCascadeHF*)d;
689
43093cfa 690 fHistNEvents->Fill(7);
691 if(fUseBit && !d->HasSelectionBit(selbit)){
692 fHistNEvents->Fill(8);
693 continue;
694 }
695
696 Double_t ptCand = d->Pt();
697 Double_t rapid=d->Y(fPdgMeson);
698 Bool_t isFidAcc=fRDCutsAnalysis->IsInFiducialAcceptance(ptCand,rapid);
a8f9e77f 699 if(!isFidAcc) continue;
20550d50 700
701 Int_t labD=-1;
702 if(fReadMC) {
8ce9400f 703 if(fPdgMeson==413){
704 labD = dCascade->MatchToMC(fPdgMeson,421,(Int_t*)pdgDgDStartoD0pi,(Int_t*)pdgDau,arrayMC);
705 } else {
706 labD = d->MatchToMC(fPdgMeson,arrayMC,nDau,(Int_t*)pdgDau);
707 }
20550d50 708 FillMCMassHistos(arrayMC,labD, countMult,nchWeight);
709 }
710
588dc385 711 Int_t passAllCuts=fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kAll,aod);
43093cfa 712 Int_t passTopolCuts=fRDCutsAnalysis->GetIsSelectedCuts();
713 if(passTopolCuts==0) continue;
714 nSelectedNoPID++;
715 fHistNEvents->Fill(9);
716 if(passAllCuts){
717 nSelectedPID++;
718 fHistNEvents->Fill(10);
719 }
3cca4c06 720 Double_t multForCand = countCorr;
8ce9400f 721
0d1e32a6 722 if(fSubtractTrackletsFromDau){
8ce9400f 723 // For the D* case, subtract only the D0 daughter tracks <=== FIXME !!
724 AliAODRecoDecayHF2Prong* d0fromDstar = NULL;
725 if(fPdgMeson==413) d0fromDstar = (AliAODRecoDecayHF2Prong*)dCascade->Get2Prong();
726
0d1e32a6 727 for(Int_t iDau=0; iDau<nDau; iDau++){
8ce9400f 728 AliAODTrack *t = NULL;
729 if(fPdgMeson==413){ t = (AliAODTrack*)d0fromDstar->GetDaughter(iDau); }
730 else{ t = (AliAODTrack*)d->GetDaughter(iDau); }
0d1e32a6 731 if(!t) continue;
732 if(t->HasPointOnITSLayer(0) && t->HasPointOnITSLayer(1)){
733 if(multForCand>0) multForCand-=1;
734 }
735 }
736 }
43093cfa 737 Bool_t isPrimary=kTRUE;
43093cfa 738 Double_t trueImpParXY=9999.;
739 Double_t impparXY=d->ImpParXY()*10000.;
740 Double_t dlen=0.1; //FIXME
3cca4c06 741 Double_t mass[2];
43093cfa 742 if(fPdgMeson==411){
743 mass[0]=d->InvMass(nDau,pdgDau);
744 mass[1]=-1.;
d2cdcb07 745 if(TMath::Abs(mass[0]-mDplusPDG)<0.02) nSelectedInMassPeak++; //20 MeV for now... FIXME
43093cfa 746 }else if(fPdgMeson==421){
747 UInt_t pdgdaughtersD0[2]={211,321};//pi,K
748 UInt_t pdgdaughtersD0bar[2]={321,211};//K,pi
749 mass[0]=d->InvMass(2,pdgdaughtersD0);
750 mass[1]=d->InvMass(2,pdgdaughtersD0bar);
d2cdcb07 751 if(TMath::Abs(mass[0]-mD0PDG)<0.02 || TMath::Abs(mass[1]-mD0PDG)<0.02 ) nSelectedInMassPeak++; //20 MeV for now... FIXME
43093cfa 752 }else if(fPdgMeson==413){
753 // FIXME
8ce9400f 754 mass[0]=dCascade->DeltaInvMass();
43093cfa 755 mass[1]=-1.;
8ce9400f 756 if(TMath::Abs(mass[0]-(mDstarPDG-mD0PDG))<0.0015) nSelectedInMassPeak++; //1 MeV for now... FIXME
43093cfa 757 }
758 for(Int_t iHyp=0; iHyp<2; iHyp++){
759 if(mass[iHyp]<0.) continue; // for D+ and D* we have 1 mass hypothesis
760 Double_t invMass=mass[iHyp];
0d1e32a6 761 Double_t arrayForSparse[5]={invMass,ptCand,impparXY,dlen,multForCand};
43093cfa 762
763 if(fReadMC){
764
8ce9400f 765 if(fPdgMeson==413){
766 labD = dCascade->MatchToMC(fPdgMeson,421,(Int_t*)pdgDgDStartoD0pi,(Int_t*)pdgDau,arrayMC);
767 } else {
768 labD = d->MatchToMC(fPdgMeson,arrayMC,nDau,(Int_t*)pdgDau);
769 }
770
43093cfa 771 Bool_t fillHisto=fDoImpPar;
772 if(labD>=0){
773 AliAODMCParticle *partD = (AliAODMCParticle*)arrayMC->At(labD);
774 Int_t code=partD->GetPdgCode();
775 if(CheckOrigin(arrayMC,partD)==5) isPrimary=kFALSE;
776 if(code<0 && iHyp==0) fillHisto=kFALSE;
777 if(code>0 && iHyp==1) fillHisto=kFALSE;
778 if(!isPrimary){
779 if(fPdgMeson==411){
780 trueImpParXY=AliVertexingHFUtils::GetTrueImpactParameterDplus(mcHeader,arrayMC,partD)*10000.;
781 }else if(fPdgMeson==421){
7826c36d 782 trueImpParXY=AliVertexingHFUtils::GetTrueImpactParameterDzero(mcHeader,arrayMC,partD)*10000.;
43093cfa 783 }else if(fPdgMeson==413){
784 trueImpParXY=0.; /// FIXME
785 }
0d1e32a6 786 Double_t arrayForSparseTrue[5]={invMass,ptCand,trueImpParXY,dlen,multForCand};
a8f9e77f 787 if(fillHisto && passAllCuts){
43093cfa 788 fHistMassPtImpPar[2]->Fill(arrayForSparse);
789 fHistMassPtImpPar[3]->Fill(arrayForSparseTrue);
790 }
791 }else{
a8f9e77f 792 if(fillHisto && passAllCuts) fHistMassPtImpPar[1]->Fill(arrayForSparse);
43093cfa 793 }
794 }else{
a8f9e77f 795 if(fillHisto && passAllCuts)fHistMassPtImpPar[4]->Fill(arrayForSparse);
43093cfa 796 }
20550d50 797
798 if(TMath::Abs(labD)==fPdgMeson && fMCOption==2) continue;
799 if(TMath::Abs(labD)!=fPdgMeson && fMCOption==1) continue;
800
43093cfa 801 }
802
f3fddcd9 803 if(fPdgMeson==421){
804 if(iHyp==0 && !(passTopolCuts&1)) continue; // candidate not passing as D0
805 if(iHyp==1 && !(passTopolCuts&2)) continue; // candidate not passing as D0bar
806 }
43093cfa 807
0d1e32a6 808 fPtVsMassVsMultNoPid->Fill(multForCand,invMass,ptCand);
19aaafff 809
810 if(fPdgMeson==421){
811 if(iHyp==0 && !(passAllCuts&1)) continue; // candidate not passing as D0
812 if(iHyp==1 && !(passAllCuts&2)) continue; // candidate not passing as D0bar
813 }
a8f9e77f 814 if(passAllCuts){
3cca4c06 815 aveMult+=multForCand;
816 nSelCand+=1.;
817 fPtVsMassVsMult->Fill(multForCand,invMass,ptCand,nchWeight);
818 fPtVsMassVsMultUncorr->Fill(countTreta1,invMass,ptCand,nchWeight);
a8f9e77f 819 // Add separation between part antipart
820 if(fPdgMeson==411){
3cca4c06 821 if(d->GetCharge()>0) fPtVsMassVsMultPart->Fill(multForCand,invMass,ptCand,nchWeight);
822 else fPtVsMassVsMultAntiPart->Fill(multForCand,invMass,ptCand,nchWeight);
a8f9e77f 823 }else if(fPdgMeson==421){
3cca4c06 824 if(passAllCuts&1) fPtVsMassVsMultPart->Fill(multForCand,invMass,ptCand,nchWeight);
825 if(passAllCuts&2) fPtVsMassVsMultAntiPart->Fill(multForCand,invMass,ptCand,nchWeight);
a8f9e77f 826 }else if(fPdgMeson==413){
8ce9400f 827 if(d->GetCharge()>0) fPtVsMassVsMultPart->Fill(multForCand,invMass,ptCand,nchWeight);
828 else fPtVsMassVsMultAntiPart->Fill(multForCand,invMass,ptCand,nchWeight);
a8f9e77f 829 }
830
831 if(fDoImpPar){
832 fHistMassPtImpPar[0]->Fill(arrayForSparse);
43093cfa 833 }
a8f9e77f 834
43093cfa 835 }
836
837 }
838 }
3cca4c06 839 if(fSubtractTrackletsFromDau && nSelCand>0){
840 aveMult/=nSelCand;
841 fCounter->StoreEvent(aod,fRDCutsAnalysis,fReadMC,(Int_t)(aveMult+0.5001));
842 }else{
20550d50 843 fCounter->StoreEvent(aod,fRDCutsAnalysis,fReadMC,(Int_t)countCorr);
3cca4c06 844 }
845
846
43093cfa 847 fCounter->StoreCandidates(aod,nSelectedNoPID,kTRUE);
848 fCounter->StoreCandidates(aod,nSelectedPID,kFALSE);
bba1a026 849 fHistNtrUnCorrEvSel->Fill(countMult,nchWeight);
850 fHistNtrCorrEvSel->Fill(countCorr,nchWeight);
851 if(nSelectedPID>0) {
852 fHistNtrUnCorrEvWithCand->Fill(countMult,nchWeight);
853 fHistNtrCorrEvWithCand->Fill(countCorr,nchWeight);
854 }
855 if(nSelectedInMassPeak>0) {
856 fHistNtrUnCorrEvWithD->Fill(countMult,nchWeight);
857 fHistNtrCorrEvWithD->Fill(countCorr,nchWeight);
858 }
d2cdcb07 859
43093cfa 860 PostData(1,fOutput);
861 PostData(2,fListCuts);
f3fddcd9 862 PostData(3,fOutputCounters);
6c06d0fe 863
43093cfa 864 return;
865}
3cca4c06 866
43093cfa 867//________________________________________________________________________
868void AliAnalysisTaskSEDvsMultiplicity::CreateImpactParameterHistos(){
869 // Histos for impact paramter study
7826c36d 870 // mass . pt , impact parameter , decay length , multiplicity
43093cfa 871
7826c36d 872 Int_t nbins[5]={fNMassBins,200,fNImpParBins,50,100};
873 Double_t xmin[5]={fLowmasslimit,0.,fLowerImpPar,0.,0.};
874 Double_t xmax[5]={fUpmasslimit,20.,fHigherImpPar,1.,100.};
43093cfa 875
876 fHistMassPtImpPar[0]=new THnSparseF("hMassPtImpParAll",
877 "Mass vs. pt vs.imppar - All",
878 5,nbins,xmin,xmax);
879 fHistMassPtImpPar[1]=new THnSparseF("hMassPtImpParPrompt",
880 "Mass vs. pt vs.imppar - promptD",
881 5,nbins,xmin,xmax);
882 fHistMassPtImpPar[2]=new THnSparseF("hMassPtImpParBfeed",
883 "Mass vs. pt vs.imppar - DfromB",
884 5,nbins,xmin,xmax);
885 fHistMassPtImpPar[3]=new THnSparseF("hMassPtImpParTrueBfeed",
886 "Mass vs. pt vs.true imppar -DfromB",
887 5,nbins,xmin,xmax);
888 fHistMassPtImpPar[4]=new THnSparseF("hMassPtImpParBkg",
889 "Mass vs. pt vs.imppar - backgr.",
890 5,nbins,xmin,xmax);
891 for(Int_t i=0; i<5;i++){
892 fOutput->Add(fHistMassPtImpPar[i]);
893 }
894}
895
896//________________________________________________________________________
897void AliAnalysisTaskSEDvsMultiplicity::Terminate(Option_t */*option*/)
898{
899 // Terminate analysis
900 //
901 if(fDebug > 1) printf("AnalysisTaskSEDvsMultiplicity: Terminate() \n");
902
903 fOutput = dynamic_cast<TList*> (GetOutputData(1));
904 if (!fOutput) {
905 printf("ERROR: fOutput not available\n");
906 return;
907 }
2247bb0c 908
43093cfa 909 fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
2247bb0c 910 if(!fHistNEvents){
911 printf("ERROR: fHistNEvents not available\n");
912 return;
913 }
c1a78ad6 914 printf("Number of Analyzed Events = %d\n",(Int_t)fHistNEvents->GetBinContent(3));
43093cfa 915
916 return;
917}
918//_________________________________________________________________________________________________
919Int_t AliAnalysisTaskSEDvsMultiplicity::CheckOrigin(TClonesArray* arrayMC, AliAODMCParticle *mcPartCandidate) const {
920 //
921 // checking whether the mother of the particles come from a charm or a bottom quark
922 //
923
924 Int_t pdgGranma = 0;
925 Int_t mother = 0;
926 mother = mcPartCandidate->GetMother();
927 Int_t istep = 0;
928 Int_t abspdgGranma =0;
929 Bool_t isFromB=kFALSE;
930 Bool_t isQuarkFound=kFALSE;
931 while (mother >0 ){
932 istep++;
933 AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
934 if (mcGranma){
935 pdgGranma = mcGranma->GetPdgCode();
936 abspdgGranma = TMath::Abs(pdgGranma);
937 if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
938 isFromB=kTRUE;
939 }
940 if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
941 mother = mcGranma->GetMother();
942 }else{
943 AliError("Failed casting the mother particle!");
944 break;
945 }
946 }
947
948 if(isFromB) return 5;
949 else return 4;
950}
951
952
953
954//____________________________________________________________________________
955TProfile* AliAnalysisTaskSEDvsMultiplicity::GetEstimatorHistogram(const AliVEvent* event){
956 // Get Estimator Histogram from period event->GetRunNumber();
957 //
958 // If you select SPD tracklets in |eta|<1 you should use type == 1
959 //
960
961 Int_t runNo = event->GetRunNumber();
b97c5dc0 962 Int_t period = -1; // pp: 0-LHC10b, 1-LHC10c, 2-LHC10d, 3-LHC10e
963 // pPb: 0-LHC13b, 1-LHC13c
964 if (fisPPbData) {
965 if (runNo>195343 && runNo<195484) period = 0;
966 if (runNo>195528 && runNo<195678) period = 1;
967 if (period < 0 || period > 1) return 0;
968 }
969 else {
970 if(runNo>114930 && runNo<117223) period = 0;
971 if(runNo>119158 && runNo<120830) period = 1;
972 if(runNo>122373 && runNo<126438) period = 2;
973 if(runNo>127711 && runNo<130841) period = 3;
974 if(period<0 || period>3) return 0;
975
976}
43093cfa 977
978 return fMultEstimatorAvg[period];
979}
3cca4c06 980
981//__________________________________________________________________________________________________
982void AliAnalysisTaskSEDvsMultiplicity::CreateMeasuredNchHisto(){
983 // creates historgam with measured multiplcity distribution in pp 7 TeV collisions (from Eur. Phys. J. C (2010) 68: 345–354)
bf27d64a 984 //
985 // for Nch > 70 the points were obtainedwith a double NBD distribution
986 // 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
987 // fit1->SetParameter(1,1.63); // k parameter
988 // fit1->SetParameter(2,12.8); // mean multiplicity
989 Double_t nchbins[82]={0.50,1.50,2.50,3.50,4.50,5.50,6.50,7.50,8.50,9.50,
3cca4c06 990 10.50,11.50,12.50,13.50,14.50,15.50,16.50,17.50,18.50,19.50,
991 20.50,21.50,22.50,23.50,24.50,25.50,26.50,27.50,28.50,29.50,
992 30.50,31.50,32.50,33.50,34.50,35.50,36.50,37.50,38.50,39.50,
993 40.50,41.50,42.50,43.50,44.50,45.50,46.50,47.50,48.50,49.50,
994 50.50,51.50,52.50,53.50,54.50,55.50,56.50,57.50,58.50,59.50,
bf27d64a 995 60.50,62.50,64.50,66.50,68.50,70.50,72.50,74.50,76.50,78.50,
996 80.50,82.50,84.50,86.50,88.50,90.50,92.50,94.50,96.50,98.50,
997 100.50,102.50};
998 Double_t pch[81]={0.062011,0.072943,0.070771,0.067245,0.062834,0.057383,0.051499,0.04591,0.041109,0.036954,
3cca4c06 999 0.03359,0.030729,0.028539,0.026575,0.024653,0.0229,0.021325,0.019768,0.018561,0.017187,
1000 0.01604,0.014836,0.013726,0.012576,0.011481,0.010393,0.009502,0.008776,0.008024,0.007452,
1001 0.006851,0.006428,0.00594,0.005515,0.005102,0.00469,0.004162,0.003811,0.003389,0.003071,
1002 0.002708,0.002422,0.002184,0.001968,0.00186,0.00165,0.001577,0.001387,0.001254,0.001118,
1003 0.001037,0.000942,0.000823,0.000736,0.000654,0.000579,0.000512,0.00049,0.00045,0.000355,
bf27d64a 1004 0.000296,0.000265,0.000193,0.00016,0.000126,0.0000851, 0.0000676,0.0000537,0.0000426, 0.0000338,
1005 0.0000268,0.0000213,0.0000166,0.0000133,0.0000106,0.00000837,0.00000662, 0.00000524,0.00000414, 0.00000327,
1006 0.00000258};
3cca4c06 1007
1008 if(fHistoMeasNch) delete fHistoMeasNch;
bf27d64a 1009 fHistoMeasNch=new TH1F("hMeaseNch","",81,nchbins);
1010 for(Int_t i=0; i<81; i++){
3cca4c06 1011 fHistoMeasNch->SetBinContent(i+1,pch[i]);
1012 fHistoMeasNch->SetBinError(i+1,0.);
1013 }
1014}
20550d50 1015
1016//__________________________________________________________________________________________________
1017void AliAnalysisTaskSEDvsMultiplicity::FillMCMassHistos(TClonesArray *arrayMC, Int_t labD, Int_t countMult,Double_t nchWeight)
1018{
1019 //
1020 // Function to fill the true MC signal
1021 //
1022
1023 if(labD>=0){
1024 AliAODMCParticle *partD = (AliAODMCParticle*)arrayMC->At(labD);
1025 Double_t mass = partD->M();
1026 Double_t pt = partD->Pt();
1027 fPtVsMassVsMultMC->Fill(countMult,mass,pt,nchWeight);
1028 }
1029
1030}