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