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