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