1 /**************************************************************************
2 * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
18 //*************************************************************************
19 // Class AliAnalysisTaskSEDvsMultiplicity
20 // AliAnalysisTaskSE for the D meson vs. multiplcity analysis
21 // Authors: Renu Bala, Zaida Conesa del Valle, Francesco Prino
22 /////////////////////////////////////////////////////////////
24 #include <TClonesArray.h>
28 #include <TDatabasePDG.h>
32 #include <THnSparse.h>
34 #include "AliAnalysisManager.h"
35 #include "AliRDHFCuts.h"
36 #include "AliRDHFCutsDplustoKpipi.h"
37 #include "AliRDHFCutsDStartoKpipi.h"
38 #include "AliRDHFCutsD0toKpi.h"
39 #include "AliAODHandler.h"
40 #include "AliAODEvent.h"
41 #include "AliAODVertex.h"
42 #include "AliAODTrack.h"
43 #include "AliAODRecoDecayHF.h"
44 #include "AliAODRecoCascadeHF.h"
45 #include "AliAnalysisVertexingHF.h"
46 #include "AliAnalysisTaskSE.h"
47 #include "AliAnalysisTaskSEDvsMultiplicity.h"
48 #include "AliNormalizationCounter.h"
49 #include "AliVertexingHFUtils.h"
50 #include "AliAODVZERO.h"
51 ClassImp(AliAnalysisTaskSEDvsMultiplicity)
54 //________________________________________________________________________
55 AliAnalysisTaskSEDvsMultiplicity::AliAnalysisTaskSEDvsMultiplicity():
62 fHistNtrEta16vsNtrEta1(0),
63 fHistNtrEta05vsNtrEta1(0),
64 fHistNtrEta03vsNtrEta1(0),
65 fHistNtrEtaV0AvsNtrEta1(0),
66 fHistNtrEtaV0MvsNtrEta1(0),
67 fHistNtrCorrEta1vsNtrRawEta1(0),
69 fHistNtrCorrVsZvtx(0),
71 fHistNtrCorrVsNchMC(0),
72 fHistNtrVsNchMCPrimary(0),
73 fHistNtrCorrVsNchMCPrimary(0),
74 fHistNtrVsNchMCPhysicalPrimary(0),
75 fHistNtrCorrVsNchMCPhysicalPrimary(0),
76 fHistGenPrimaryParticlesInelGt0(0),
77 fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary(0),
78 fHistNtrUnCorrEvSel(0),
79 fHistNtrUnCorrEvWithCand(0),
80 fHistNtrUnCorrEvWithD(0),
82 fHistNtrCorrEvWithCand(0),
83 fHistNtrCorrEvWithD(0),
85 fPtVsMassVsMultNoPid(0),
86 fPtVsMassVsMultUncorr(0),
87 fPtVsMassVsMultPart(0),
88 fPtVsMassVsMultAntiPart(0),
104 fSubtractTrackletsFromDau(kFALSE),
105 fKeepCorrPlots(kFALSE),
106 fUseNchWeight(kFALSE),
111 fMultiplicityEstimator(kNtrk10),
112 fMCPrimariesEstimator(kEta10)
114 // Default constructor
115 for(Int_t i=0; i<5; i++) fHistMassPtImpPar[i]=0;
116 for(Int_t i=0; i<4; i++) fMultEstimatorAvg[i]=0;
119 //________________________________________________________________________
120 AliAnalysisTaskSEDvsMultiplicity::AliAnalysisTaskSEDvsMultiplicity(const char *name, Int_t pdgMeson,AliRDHFCuts *cuts, Bool_t switchPPb):
121 AliAnalysisTaskSE(name),
127 fHistNtrEta16vsNtrEta1(0),
128 fHistNtrEta05vsNtrEta1(0),
129 fHistNtrEta03vsNtrEta1(0),
130 fHistNtrEtaV0AvsNtrEta1(0),
131 fHistNtrEtaV0MvsNtrEta1(0),
132 fHistNtrCorrEta1vsNtrRawEta1(0),
134 fHistNtrCorrVsZvtx(0),
136 fHistNtrCorrVsNchMC(0),
137 fHistNtrVsNchMCPrimary(0),
138 fHistNtrCorrVsNchMCPrimary(0),
139 fHistNtrVsNchMCPhysicalPrimary(0),
140 fHistNtrCorrVsNchMCPhysicalPrimary(0),
141 fHistGenPrimaryParticlesInelGt0(0),
142 fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary(0),
143 fHistNtrUnCorrEvSel(0),
144 fHistNtrUnCorrEvWithCand(0),
145 fHistNtrUnCorrEvWithD(0),
146 fHistNtrCorrEvSel(0),
147 fHistNtrCorrEvWithCand(0),
148 fHistNtrCorrEvWithD(0),
150 fPtVsMassVsMultNoPid(0),
151 fPtVsMassVsMultUncorr(0),
152 fPtVsMassVsMultPart(0),
153 fPtVsMassVsMultAntiPart(0),
154 fPtVsMassVsMultMC(0),
156 fLowmasslimit(1.765),
158 fRDCutsAnalysis(cuts),
163 fLowerImpPar(-2000.),
164 fHigherImpPar(2000.),
167 fisPPbData(switchPPb),
169 fSubtractTrackletsFromDau(kFALSE),
170 fKeepCorrPlots(kFALSE),
171 fUseNchWeight(kFALSE),
176 fMultiplicityEstimator(kNtrk10),
177 fMCPrimariesEstimator(kEta10)
180 // Standard constructor
183 for(Int_t i=0; i<5; i++) fHistMassPtImpPar[i]=0;
184 for(Int_t i=0; i<4; i++) fMultEstimatorAvg[i]=0;
187 SetMassLimits(0.12,0.2);
190 SetMassLimits(fPdgMeson,0.1);
192 // Default constructor
193 // Otput slot #1 writes into a TList container
194 DefineOutput(1,TList::Class()); //My private output
195 // Output slot #2 writes cut to private output
196 DefineOutput(2,TList::Class());
197 // Output slot #3 writes cut to private output
198 DefineOutput(3,TList::Class());
199 // Output slot #4 writes cut to private output
200 DefineOutput(4,TList::Class());
202 //________________________________________________________________________
203 AliAnalysisTaskSEDvsMultiplicity::~AliAnalysisTaskSEDvsMultiplicity()
211 delete fListProfiles;
212 delete fRDCutsAnalysis;
215 for(Int_t i=0; i<4; i++) {
216 if (fMultEstimatorAvg[i]) delete fMultEstimatorAvg[i];
219 for(Int_t i=0; i<5; i++){
220 delete fHistMassPtImpPar[i];
222 if(fHistoMCNch) delete fHistoMCNch;
223 if(fHistoMeasNch) delete fHistoMeasNch;
226 //_________________________________________________________________
227 void AliAnalysisTaskSEDvsMultiplicity::SetMassLimits(Double_t lowlimit, Double_t uplimit){
228 // set invariant mass limits
229 if(uplimit>lowlimit){
230 fLowmasslimit = lowlimit;
231 fUpmasslimit = uplimit;
233 AliError("Wrong mass limits: upper value should be larger than lower one");
236 //_________________________________________________________________
237 void AliAnalysisTaskSEDvsMultiplicity::SetMassLimits(Int_t pdg, Double_t range){
238 // set invariant mass limits
239 Double_t mass=TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdg))->Mass();
240 SetMassLimits(mass-range,mass+range);
242 //________________________________________________________________________
243 void AliAnalysisTaskSEDvsMultiplicity::Init(){
247 printf("AnalysisTaskSEDvsMultiplicity::Init() \n");
249 if(fUseNchWeight && !fReadMC){ AliFatal("Nch weights can only be used in MC mode"); return; }
250 if(fUseNchWeight && !fHistoMCNch){ AliFatal("Nch weights can only be used without histogram"); return; }
252 fListCuts=new TList();
253 fListCuts->SetOwner();
254 fListCuts->SetName("CutsList");
258 AliRDHFCutsDplustoKpipi* copycut=new AliRDHFCutsDplustoKpipi(*(static_cast<AliRDHFCutsDplustoKpipi*>(fRDCutsAnalysis)));
259 copycut->SetName("AnalysisCutsDplus");
260 fListCuts->Add(copycut);
261 }else if(fPdgMeson==421){
262 AliRDHFCutsD0toKpi* copycut=new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fRDCutsAnalysis)));
263 copycut->SetName("AnalysisCutsDzero");
264 fListCuts->Add(copycut);
265 }else if(fPdgMeson==413){
266 AliRDHFCutsDStartoKpipi* copycut=new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fRDCutsAnalysis)));
267 copycut->SetName("AnalysisCutsDStar");
268 fListCuts->Add(copycut);
270 PostData(2,fListCuts);
272 fListProfiles = new TList();
273 fListProfiles->SetOwner();
276 if (fisPPbData) {period[0]="LHC13b"; period[1]="LHC13c"; nProfiles = 2;}
277 else {period[0]="LHC10b"; period[1]="LHC10c"; period[2]="LHC10d"; period[3]="LHC10e"; nProfiles = 4;}
279 for(Int_t i=0; i<nProfiles; i++){
280 if(fMultEstimatorAvg[i]){
281 TProfile* hprof=new TProfile(*fMultEstimatorAvg[i]);
282 hprof->SetName(Form("ProfileTrkVsZvtx%s\n",period[i].Data()));
283 fListProfiles->Add(hprof);
286 PostData(4,fListProfiles);
291 //________________________________________________________________________
292 void AliAnalysisTaskSEDvsMultiplicity::UserCreateOutputObjects()
294 // Create the output container
296 if(fDebug > 1) printf("AnalysisTaskSEDvsMultiplicity::UserCreateOutputObjects() \n");
298 // Several histograms are more conveniently managed in a TList
299 fOutput = new TList();
301 fOutput->SetName("OutputHistos");
303 Int_t nMultBins = 200;
304 Float_t firstMultBin = -0.5;
305 Float_t lastMultBin = 199.5;
310 if(fMultiplicityEstimator==kVZERO || fMultiplicityEstimator==kVZEROA) {
315 fHistNtrUnCorrEvSel = new TH1F("hNtrUnCorrEvSel","Uncorrected tracklets multiplicity for selected events; Tracklets ; Entries",nMultBins,firstMultBin,lastMultBin);
316 fHistNtrUnCorrEvWithCand = new TH1F("hNtrUnCorrEvWithCand", "Uncorrected Tracklets multiplicity for events with D candidates; Tracklets ; Entries",nMultBins,firstMultBin,lastMultBin);// Total multiplicity
317 fHistNtrUnCorrEvWithD = new TH1F("hNtrUnCorrEvWithD","Uncorrected Tracklets multiplicity for events with D in mass region ; Tracklets ; Entries",nMultBins,firstMultBin,lastMultBin); //
318 fHistNtrCorrEvSel = new TH1F("hNtrCorrEvSel","Corrected tracklets multiplicity for selected events; Tracklets ; Entries",nMultBins,firstMultBin,lastMultBin);
319 fHistNtrCorrEvWithCand = new TH1F("hNtrCorrEvWithCand", "Tracklets multiplicity for events with D candidates; Tracklets ; Entries",nMultBins,firstMultBin,lastMultBin);// Total multiplicity
320 fHistNtrCorrEvWithD = new TH1F("hNtrCorrEvWithD", "Tracklets multiplicity for events with D in mass region ; Tracklets ; Entries",nMultBins,firstMultBin,lastMultBin); //
323 fHistNtrEta16vsNtrEta1 = new TH2F("hNtrEta16vsNtrEta1","Uncorrected Eta1.6 vs Eta1.0; Ntracklets #eta<1.0; Ntracklets #eta<1.6",nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); //eta 1.6 vs eta 1.0 histogram
324 fHistNtrEta05vsNtrEta1 = new TH2F("hNtrEta05vsNtrEta1","Uncorrected Eta0.5 vs Eta1.0; Ntracklets #eta<1.0; Ntracklets #eta<0.5",nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); //eta 0.5 vs eta 1.0 histogram
325 fHistNtrEta03vsNtrEta1 = new TH2F("hNtrEta03vsNtrEta1","Uncorrected Eta0.3 vs Eta1.0; Ntracklets #eta<1.0; Ntracklets #eta<0.3",nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); //eta 0.3 vs eta 1.0 histogram
326 fHistNtrEtaV0AvsNtrEta1 = new TH2F("hNtrEtaV0AvsNtrEta1","Uncorrected Eta-V0A vs Eta1.0; Ntracklets #eta<1.0; Ntracklets #eta<0.3",nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); //eta V0A vs eta 1.0 histogram
327 fHistNtrEtaV0MvsNtrEta1 = new TH2F("hNtrEtaV0MvsNtrEta1","Uncorrected Eta-V0M vs Eta1.0; Ntracklets #eta<1.0; Ntracklets #eta<0.3",nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); //eta V0M vs eta 1.0 histogram
328 fHistNtrCorrEta1vsNtrRawEta1 = new TH2F("hNtrCorrEta1vsNtrRawEta1","Corrected Eta1 vs Eta1.0; Ntracklets #eta<1.0 corrected; Ntracklets #eta<1",nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); //eta 1.6 vs eta 1.0 histogram
330 fHistNtrVsZvtx = new TH2F("hNtrVsZvtx","Ntracklet vs VtxZ; VtxZ;N_{tracklet};",300,-15,15,nMultBins,firstMultBin,lastMultBin); //
331 fHistNtrCorrVsZvtx = new TH2F("hNtrCorrVsZvtx","Ntracklet vs VtxZ; VtxZ;N_{tracklet};",300,-15,15,nMultBins,firstMultBin,lastMultBin); //
333 fHistNtrVsNchMC = new TH2F("hNtrVsNchMC","Ntracklet vs NchMC; Nch;N_{tracklet};",nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); //
334 fHistNtrCorrVsNchMC = new TH2F("hNtrCorrVsNchMC","Ntracklet vs Nch; Nch;N_{tracklet};",nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); //
336 fHistNtrVsNchMCPrimary = new TH2F("hNtrVsNchMCPrimary","Ntracklet vs Nch (Primary); Nch (Primary);N_{tracklet};",nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); //
337 fHistNtrCorrVsNchMCPrimary = new TH2F("hNtrCorrVsNchMCPrimary","Ntracklet vs Nch (Primary); Nch(Primary) ;N_{tracklet};",nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); //
339 fHistNtrVsNchMCPhysicalPrimary = new TH2F("hNtrVsNchMCPhysicalPrimary","Ntracklet vs Nch (Physical Primary); Nch (Physical Primary);N_{tracklet};",nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); //
340 fHistNtrCorrVsNchMCPhysicalPrimary = new TH2F("hNtrCorrVsMCPhysicalPrimary","Ntracklet vs Nch (Physical Primary); Nch (Physical Primary);N_{tracklet};",nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); //
342 fHistGenPrimaryParticlesInelGt0 = new TH1F("hGenPrimaryParticlesInelGt0","Multiplcity of generated charged particles ; Nparticles ; Entries",nMultBins,firstMultBin,lastMultBin);
344 fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary = new TH3F("fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary", "MC: Nch (Physical Primary) vs Nch (Primary) vs Nch (Generated); Nch (Generated); Nch (Primary); Nch (Physical Primary)",nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin);
346 fHistNtrUnCorrEvSel->Sumw2();
347 fHistNtrUnCorrEvWithCand->Sumw2();
348 fHistNtrUnCorrEvWithD->Sumw2();
349 fHistNtrCorrEvSel->Sumw2();
350 fHistNtrCorrEvWithCand->Sumw2();
351 fHistNtrCorrEvWithD->Sumw2();
352 fHistGenPrimaryParticlesInelGt0->Sumw2();
353 fOutput->Add(fHistNtrUnCorrEvSel);
354 fOutput->Add(fHistNtrUnCorrEvWithCand);
355 fOutput->Add(fHistNtrUnCorrEvWithD);
356 fOutput->Add(fHistNtrCorrEvSel);
357 fOutput->Add(fHistNtrCorrEvWithCand);
358 fOutput->Add(fHistNtrCorrEvWithD);
360 fOutput->Add(fHistNtrEta16vsNtrEta1);
361 fOutput->Add(fHistNtrEta05vsNtrEta1);
362 fOutput->Add(fHistNtrEta03vsNtrEta1);
363 fOutput->Add(fHistNtrEtaV0AvsNtrEta1);
364 fOutput->Add(fHistNtrEtaV0MvsNtrEta1);
365 fOutput->Add(fHistNtrCorrEta1vsNtrRawEta1);
367 fOutput->Add(fHistNtrVsZvtx);
368 fOutput->Add(fHistNtrCorrVsZvtx);
370 fOutput->Add(fHistNtrVsNchMC);
371 fOutput->Add(fHistNtrCorrVsNchMC);
372 fOutput->Add(fHistNtrVsNchMCPrimary);
373 fOutput->Add(fHistNtrCorrVsNchMCPrimary);
374 fOutput->Add(fHistNtrVsNchMCPhysicalPrimary);
375 fOutput->Add(fHistNtrCorrVsNchMCPhysicalPrimary);
376 fOutput->Add(fHistGenPrimaryParticlesInelGt0);
377 fOutput->Add(fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary);
380 fHistNEvents = new TH1F("fHistNEvents", "number of events ",11,-0.5,10.5);
381 fHistNEvents->GetXaxis()->SetBinLabel(1,"nEvents total");
382 fHistNEvents->GetXaxis()->SetBinLabel(2,"nEvents with Z vertex");
383 fHistNEvents->GetXaxis()->SetBinLabel(3,"nEvents selected");
384 fHistNEvents->GetXaxis()->SetBinLabel(4,"Rejected due to trigger");
385 fHistNEvents->GetXaxis()->SetBinLabel(5,"Rejected due to phys sel");
386 fHistNEvents->GetXaxis()->SetBinLabel(6,"Rejected due to vertex cuts");
387 fHistNEvents->GetXaxis()->SetBinLabel(7,"Rejected due to pileup");
388 fHistNEvents->GetXaxis()->SetBinLabel(8,"Total no. of candidate");
389 fHistNEvents->GetXaxis()->SetBinLabel(9,"no. of cand wo bitmask");
390 fHistNEvents->GetXaxis()->SetBinLabel(10,"D after cuts (No PID)");
391 fHistNEvents->GetXaxis()->SetBinLabel(11,"D after cuts + PID)");
392 fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);
393 fHistNEvents->Sumw2();
394 fHistNEvents->SetMinimum(0);
395 fOutput->Add(fHistNEvents);
397 fPtVsMassVsMult=new TH3F("hPtVsMassvsMult", "D candidates: p_{t} vs mass vs tracklets multiplicity; Tracklets; Mass M [GeV/c^{2}]; p_{t} [GeV/c]",nMultBins,firstMultBin,lastMultBin,fNMassBins,fLowmasslimit,fUpmasslimit,48,0.,24.);
399 fPtVsMassVsMultNoPid=new TH3F("hPtVsMassvsMultNoPid", "D candidates: p_{t} vs mass vs tracklets multiplicity; Tracklets; Mass M [GeV/c^{2}]; p_{t} [GeV/c]",nMultBins,firstMultBin,lastMultBin,fNMassBins,fLowmasslimit,fUpmasslimit,48,0.,24.);
401 fPtVsMassVsMultUncorr=new TH3F("hPtVsMassvsMultUncorr", "D candidates: p_{t} vs mass vs tracklets multiplicity; Tracklets; Mass M [GeV/c^{2}]; p_{t} [GeV/c]",nMultBins,firstMultBin,lastMultBin,fNMassBins,fLowmasslimit,fUpmasslimit,48,0.,24.);
403 fPtVsMassVsMultPart=new TH3F("hPtVsMassvsMultPart", "D candidates: p_{t} vs mass vs tracklets multiplicity; Tracklets; Mass M [GeV/c^{2}]; p_{t} [GeV/c]",nMultBins,firstMultBin,lastMultBin,fNMassBins,fLowmasslimit,fUpmasslimit,48,0.,24.);
405 fPtVsMassVsMultAntiPart=new TH3F("hPtVsMassvsMultAntiPart", "D candidates: p_{t} vs mass vs tracklets multiplicity; Tracklets; Mass M [GeV/c^{2}]; p_{t} [GeV/c]",nMultBins,firstMultBin,lastMultBin,fNMassBins,fLowmasslimit,fUpmasslimit,48,0.,24.);
407 fPtVsMassVsMultMC=new TH3F("hPtVsMassvsMultMC", "D true candidates: p_{t} vs mass vs tracklets multiplicity; Tracklets; Mass M [GeV/c^{2}]; p_{t} [GeV/c]",nMultBins,firstMultBin,lastMultBin,fNMassBins,fLowmasslimit,fUpmasslimit,48,0.,24.);
409 fOutput->Add(fPtVsMassVsMult);
410 fOutput->Add(fPtVsMassVsMultUncorr);
411 fOutput->Add(fPtVsMassVsMultNoPid);
412 fOutput->Add(fPtVsMassVsMultPart);
413 fOutput->Add(fPtVsMassVsMultAntiPart);
414 fOutput->Add(fPtVsMassVsMultMC);
416 if(fDoImpPar) CreateImpactParameterHistos();
418 fCounter = new AliNormalizationCounter("NormCounterCorrMult");
419 fCounter->SetStudyMultiplicity(kTRUE,1.);
422 fCounterU = new AliNormalizationCounter("NormCounterUnCorrMult");
423 fCounterU->SetStudyMultiplicity(kTRUE,1.);
426 fOutputCounters = new TList();
427 fOutputCounters->SetOwner();
428 fOutputCounters->SetName("OutputCounters");
429 fOutputCounters->Add(fCounter);
430 fOutputCounters->Add(fCounterU);
433 PostData(2,fListCuts);
434 PostData(3,fOutputCounters);
435 PostData(4,fListProfiles);
437 if(fUseNchWeight) CreateMeasuredNchHisto();
444 //________________________________________________________________________
445 void AliAnalysisTaskSEDvsMultiplicity::UserExec(Option_t */*option*/)
447 // Execute analysis for current event:
448 // heavy flavor candidates association to MC truth
450 AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
452 // AliAODTracklets* tracklets = aod->GetTracklets();
453 //Int_t ntracklets = tracklets->GetNumberOfTracklets();
456 TClonesArray *arrayCand = 0;
457 TString arrayName="";
462 arrayName="Charm3Prong";
463 pdgDau[0]=211; pdgDau[1]=321; pdgDau[2]=211;
465 selbit=AliRDHFCuts::kDplusCuts;
466 }else if(fPdgMeson==421){
468 pdgDau[0]=211; pdgDau[1]=321; pdgDau[2]=0;
470 selbit=AliRDHFCuts::kD0toKpiCuts;
471 }else if(fPdgMeson==413){
473 pdgDau[0]=321; pdgDau[1]=211; pdgDau[2]=0; // Quoting here D0 daughters (D* ones on another variable later)
475 selbit=AliRDHFCuts::kDstarCuts;
478 if(!aod && AODEvent() && IsStandardAOD()) {
479 // In case there is an AOD handler writing a standard AOD, use the AOD
480 // event in memory rather than the input (ESD) event.
481 aod = dynamic_cast<AliAODEvent*> (AODEvent());
482 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
483 // have to taken from the AOD event hold by the AliAODExtension
484 AliAODHandler* aodHandler = (AliAODHandler*)
485 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
486 if(aodHandler->GetExtensions()) {
487 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
488 AliAODEvent *aodFromExt = ext->GetAOD();
489 arrayCand=(TClonesArray*)aodFromExt->GetList()->FindObject(arrayName.Data());
492 arrayCand=(TClonesArray*)aod->GetList()->FindObject(arrayName.Data());
495 if(!aod || !arrayCand) {
496 printf("AliAnalysisTaskSEDvsMultiplicity::UserExec: Charm3Prong branch not found!\n");
500 if(fisPPbData && fReadMC){
501 Int_t runnumber = aod->GetRunNumber();
502 if(aod->GetTriggerMask()==0 &&
503 (runnumber>=195344 && runnumber<=195677)){
504 AliDebug(3,"Event rejected because of null trigger mask");
510 // fix for temporary bug in ESDfilter
511 // the AODs with null vertex pointer didn't pass the PhysSel
512 if(!aod->GetPrimaryVertex()||TMath::Abs(aod->GetMagneticField())<0.001) return;
514 Int_t countTreta1=AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aod,-1.,1.);
515 Int_t countTreta03=AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aod,-0.3,0.3);
516 Int_t countTreta05=AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aod,-0.5,0.5);
517 Int_t countTreta16=AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aod,-1.6,1.6);
521 AliAODVZERO *vzeroAOD = (AliAODVZERO*)aod->GetVZEROData();
523 vzeroMult = vzeroAOD->GetMTotV0A() + vzeroAOD->GetMTotV0C();
524 vzeroMultA = vzeroAOD->GetMTotV0A();
527 Int_t countMult = countTreta1;
528 if(fMultiplicityEstimator==kNtrk03) { countMult = countTreta03; }
529 if(fMultiplicityEstimator==kNtrk05) { countMult = countTreta05; }
530 if(fMultiplicityEstimator==kNtrk10to16) { countMult = countTreta16 - countTreta1; }
531 if(fMultiplicityEstimator==kVZERO) { countMult = vzeroMult; }
532 if(fMultiplicityEstimator==kVZEROA) { countMult = vzeroMultA; }
535 fCounterU->StoreEvent(aod,fRDCutsAnalysis,fReadMC,countMult);
536 fHistNEvents->Fill(0); // count event
538 Double_t countTreta1corr=countTreta1;
539 Double_t countCorr=countMult;
540 AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
542 // FIX ME: No correction to the VZERO !!
543 if(vtx1 && (fMultiplicityEstimator!=kVZERO) && (fMultiplicityEstimator!=kVZEROA)){
544 if(vtx1->GetNContributors()>0){
545 fHistNEvents->Fill(1);
546 TProfile* estimatorAvg = GetEstimatorHistogram(aod);
548 countTreta1corr=AliVertexingHFUtils::GetCorrectedNtracklets(estimatorAvg,countTreta1,vtx1->GetZ(),fRefMult);
549 countCorr=AliVertexingHFUtils::GetCorrectedNtracklets(estimatorAvg,countMult,vtx1->GetZ(),fRefMult);
555 Bool_t isEvSel=fRDCutsAnalysis->IsEventSelected(aod);
557 if(fRDCutsAnalysis->GetWhyRejection()==5) fHistNEvents->Fill(3);
558 if(fRDCutsAnalysis->GetWhyRejection()==7) fHistNEvents->Fill(4);
559 if(fRDCutsAnalysis->GetWhyRejection()==6) fHistNEvents->Fill(5);
560 if(fRDCutsAnalysis->GetWhyRejection()==1) fHistNEvents->Fill(6);
565 fHistNtrEta16vsNtrEta1->Fill(countTreta1,countTreta16);
566 fHistNtrEta05vsNtrEta1->Fill(countTreta1,countTreta05);
567 fHistNtrEta03vsNtrEta1->Fill(countTreta1,countTreta03);
568 fHistNtrEtaV0AvsNtrEta1->Fill(countTreta1,vzeroMultA);
569 fHistNtrEtaV0MvsNtrEta1->Fill(countTreta1,vzeroMult);
570 fHistNtrCorrEta1vsNtrRawEta1->Fill(countTreta1,countTreta1corr);
573 fHistNtrVsZvtx->Fill(vtx1->GetZ(),countMult);
574 fHistNtrCorrVsZvtx->Fill(vtx1->GetZ(),countCorr);
577 TClonesArray *arrayMC=0;
578 AliAODMCHeader *mcHeader=0;
580 Double_t nchWeight=1.0;
585 arrayMC = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
587 printf("AliAnalysisTaskSEDvsMultiplicity::UserExec: MC particles branch not found!\n");
591 mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
593 printf("AliAnalysisTaskSEDvsMultiplicity::UserExec: MC header branch not found!\n");
598 Int_t nChargedMC=AliVertexingHFUtils::GetGeneratedMultiplicityInEtaRange(arrayMC,-1.0,1.0);
599 Int_t nChargedMCPrimary=AliVertexingHFUtils::GetGeneratedPrimariesInEtaRange(arrayMC,-1.0,1.0);
600 Int_t nChargedMCPhysicalPrimary=AliVertexingHFUtils::GetGeneratedPhysicalPrimariesInEtaRange(arrayMC,-1.0,1.0);
602 // Compute the Nch weights (reference is Ntracklets within |eta|<1.0)
604 Double_t tmpweight = 1.0;
605 if(nChargedMCPhysicalPrimary<=0) tmpweight = 0.0;
607 Double_t pMeas = fHistoMeasNch->GetBinContent(fHistoMeasNch->FindBin(nChargedMCPhysicalPrimary));
608 // printf(" pMeas=%2.2f and histo MCNch %s \n",pMeas,fHistoMCNch);
609 Double_t pMC = fHistoMCNch->GetBinContent(fHistoMCNch->FindBin(nChargedMCPhysicalPrimary));
610 tmpweight = pMC>0 ? pMeas/pMC : 0.;
612 nchWeight *= tmpweight;
613 AliDebug(2,Form("Using Nch weights, Mult=%d Weight=%f\n",nChargedMCPhysicalPrimary,nchWeight));
616 // Now recompute the variables in case another MC estimator is considered
617 Int_t nChargedMCEta10 = nChargedMC;
618 Int_t nChargedMCEta16 = AliVertexingHFUtils::GetGeneratedMultiplicityInEtaRange(arrayMC,-1.6,1.6);
619 Int_t nChargedMCEta05 = AliVertexingHFUtils::GetGeneratedMultiplicityInEtaRange(arrayMC,-0.5,0.5);
620 Int_t nChargedMCEta03 = AliVertexingHFUtils::GetGeneratedMultiplicityInEtaRange(arrayMC,-0.3,0.3);
621 Int_t nChargedMCEtam37tm17 = AliVertexingHFUtils::GetGeneratedMultiplicityInEtaRange(arrayMC,-3.7,-1.7);
622 Int_t nChargedMCEta28t51 = AliVertexingHFUtils::GetGeneratedMultiplicityInEtaRange(arrayMC,2.8,5.1);
623 Int_t nChargedMCPrimaryEta10 = nChargedMCPrimary;
624 Int_t nChargedMCPrimaryEta16 = AliVertexingHFUtils::GetGeneratedPrimariesInEtaRange(arrayMC,-1.6,1.6);
625 Int_t nChargedMCPrimaryEta05 = AliVertexingHFUtils::GetGeneratedPrimariesInEtaRange(arrayMC,-0.5,0.5);
626 Int_t nChargedMCPrimaryEta03 = AliVertexingHFUtils::GetGeneratedPrimariesInEtaRange(arrayMC,-0.3,0.3);
627 Int_t nChargedMCPrimaryEtam37tm17 = AliVertexingHFUtils::GetGeneratedPrimariesInEtaRange(arrayMC,-3.7,-1.7);
628 Int_t nChargedMCPrimaryEta28t51 = AliVertexingHFUtils::GetGeneratedPrimariesInEtaRange(arrayMC,2.8,5.1);
629 Int_t nChargedMCPhysicalPrimaryEta10 = nChargedMCPhysicalPrimary;
630 Int_t nChargedMCPhysicalPrimaryEta16 = AliVertexingHFUtils::GetGeneratedPhysicalPrimariesInEtaRange(arrayMC,-1.6,1.6);
631 Int_t nChargedMCPhysicalPrimaryEta05 = AliVertexingHFUtils::GetGeneratedPhysicalPrimariesInEtaRange(arrayMC,-0.5,0.5);
632 Int_t nChargedMCPhysicalPrimaryEta03 = AliVertexingHFUtils::GetGeneratedPhysicalPrimariesInEtaRange(arrayMC,-0.3,0.3);
633 Int_t nChargedMCPhysicalPrimaryEtam37tm17 = AliVertexingHFUtils::GetGeneratedPhysicalPrimariesInEtaRange(arrayMC,-3.7,-1.7);
634 Int_t nChargedMCPhysicalPrimaryEta28t51 = AliVertexingHFUtils::GetGeneratedPhysicalPrimariesInEtaRange(arrayMC,2.8,5.1);
635 if(fMCPrimariesEstimator==kEta10to16){
636 nChargedMC = nChargedMCEta16 - nChargedMCEta10;
637 nChargedMCPrimary = nChargedMCPrimaryEta16 - nChargedMCPrimaryEta10;
638 nChargedMCPhysicalPrimary = nChargedMCPhysicalPrimaryEta16 - nChargedMCPhysicalPrimaryEta10;
639 } else if(fMCPrimariesEstimator==kEta05){
640 nChargedMC = nChargedMCEta05;
641 nChargedMCPrimary = nChargedMCPrimaryEta05;
642 nChargedMCPhysicalPrimary = nChargedMCPhysicalPrimaryEta05;
643 } else if(fMCPrimariesEstimator==kEta03){
644 nChargedMC = nChargedMCEta03;
645 nChargedMCPrimary = nChargedMCPrimaryEta03;
646 nChargedMCPhysicalPrimary = nChargedMCPhysicalPrimaryEta03;
647 } else if(fMCPrimariesEstimator==kEtaVZERO){
648 nChargedMC = nChargedMCEtam37tm17 + nChargedMCEta28t51;
649 nChargedMCPrimary = nChargedMCPrimaryEtam37tm17 + nChargedMCPrimaryEta28t51;
650 nChargedMCPhysicalPrimary = nChargedMCPhysicalPrimaryEtam37tm17 + nChargedMCPhysicalPrimaryEta28t51;
651 } else if(fMCPrimariesEstimator==kEtaVZEROA){
652 nChargedMC = nChargedMCEta28t51;
653 nChargedMCPrimary = nChargedMCPrimaryEta28t51;
654 nChargedMCPhysicalPrimary = nChargedMCPhysicalPrimaryEta28t51;
657 // Here fill the MC correlation plots
658 if(nChargedMCPhysicalPrimary>0){ // INEL>0 for |eta|<1
659 fHistGenPrimaryParticlesInelGt0->Fill(nChargedMCPhysicalPrimary,nchWeight);
662 fHistNtrVsNchMC->Fill(nChargedMC,countMult,nchWeight);
663 fHistNtrCorrVsNchMC->Fill(nChargedMC,countCorr,nchWeight);
665 fHistNtrVsNchMCPrimary->Fill(nChargedMCPrimary,countMult,nchWeight);
666 fHistNtrCorrVsNchMCPrimary->Fill(nChargedMCPrimary,countCorr,nchWeight);
668 fHistNtrVsNchMCPhysicalPrimary->Fill(nChargedMCPhysicalPrimary,countMult,nchWeight);
669 fHistNtrCorrVsNchMCPhysicalPrimary->Fill(nChargedMCPhysicalPrimary,countCorr,nchWeight);
671 fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary->Fill(nChargedMC,nChargedMCPrimary,nChargedMCPhysicalPrimary,nchWeight);
674 Int_t nCand = arrayCand->GetEntriesFast();
675 Int_t nSelectedNoPID=0,nSelectedPID=0,nSelectedInMassPeak=0;
676 Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();
677 Double_t mDplusPDG = TDatabasePDG::Instance()->GetParticle(411)->Mass();
678 Double_t mDstarPDG = TDatabasePDG::Instance()->GetParticle(413)->Mass();
680 // pdg of daughters needed for D* too
681 UInt_t pdgDgDStartoD0pi[2]={421,211};
684 Double_t nSelCand=0.;
685 for (Int_t iCand = 0; iCand < nCand; iCand++) {
686 AliAODRecoDecayHF *d = (AliAODRecoDecayHF*)arrayCand->UncheckedAt(iCand);
687 AliAODRecoCascadeHF *dCascade = NULL;
688 if(fPdgMeson==413) dCascade = (AliAODRecoCascadeHF*)d;
690 fHistNEvents->Fill(7);
691 if(fUseBit && !d->HasSelectionBit(selbit)){
692 fHistNEvents->Fill(8);
696 Double_t ptCand = d->Pt();
697 Double_t rapid=d->Y(fPdgMeson);
698 Bool_t isFidAcc=fRDCutsAnalysis->IsInFiducialAcceptance(ptCand,rapid);
699 if(!isFidAcc) continue;
704 labD = dCascade->MatchToMC(fPdgMeson,421,(Int_t*)pdgDgDStartoD0pi,(Int_t*)pdgDau,arrayMC);
706 labD = d->MatchToMC(fPdgMeson,arrayMC,nDau,(Int_t*)pdgDau);
708 FillMCMassHistos(arrayMC,labD, countMult,nchWeight);
711 Int_t passAllCuts=fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kAll,aod);
712 Int_t passTopolCuts=fRDCutsAnalysis->GetIsSelectedCuts();
713 if(passTopolCuts==0) continue;
715 fHistNEvents->Fill(9);
718 fHistNEvents->Fill(10);
720 Double_t multForCand = countCorr;
722 if(fSubtractTrackletsFromDau){
723 // For the D* case, subtract only the D0 daughter tracks <=== FIXME !!
724 AliAODRecoDecayHF2Prong* d0fromDstar = NULL;
725 if(fPdgMeson==413) d0fromDstar = (AliAODRecoDecayHF2Prong*)dCascade->Get2Prong();
727 for(Int_t iDau=0; iDau<nDau; iDau++){
728 AliAODTrack *t = NULL;
729 if(fPdgMeson==413){ t = (AliAODTrack*)d0fromDstar->GetDaughter(iDau); }
730 else{ t = (AliAODTrack*)d->GetDaughter(iDau); }
732 if(t->HasPointOnITSLayer(0) && t->HasPointOnITSLayer(1)){
733 if(multForCand>0) multForCand-=1;
737 Bool_t isPrimary=kTRUE;
738 Double_t trueImpParXY=9999.;
739 Double_t impparXY=d->ImpParXY()*10000.;
740 Double_t dlen=0.1; //FIXME
743 mass[0]=d->InvMass(nDau,pdgDau);
745 if(TMath::Abs(mass[0]-mDplusPDG)<0.02) nSelectedInMassPeak++; //20 MeV for now... FIXME
746 }else if(fPdgMeson==421){
747 UInt_t pdgdaughtersD0[2]={211,321};//pi,K
748 UInt_t pdgdaughtersD0bar[2]={321,211};//K,pi
749 mass[0]=d->InvMass(2,pdgdaughtersD0);
750 mass[1]=d->InvMass(2,pdgdaughtersD0bar);
751 if(TMath::Abs(mass[0]-mD0PDG)<0.02 || TMath::Abs(mass[1]-mD0PDG)<0.02 ) nSelectedInMassPeak++; //20 MeV for now... FIXME
752 }else if(fPdgMeson==413){
754 mass[0]=dCascade->DeltaInvMass();
756 if(TMath::Abs(mass[0]-(mDstarPDG-mD0PDG))<0.0015) nSelectedInMassPeak++; //1 MeV for now... FIXME
758 for(Int_t iHyp=0; iHyp<2; iHyp++){
759 if(mass[iHyp]<0.) continue; // for D+ and D* we have 1 mass hypothesis
760 Double_t invMass=mass[iHyp];
761 Double_t arrayForSparse[5]={invMass,ptCand,impparXY,dlen,multForCand};
766 labD = dCascade->MatchToMC(fPdgMeson,421,(Int_t*)pdgDgDStartoD0pi,(Int_t*)pdgDau,arrayMC);
768 labD = d->MatchToMC(fPdgMeson,arrayMC,nDau,(Int_t*)pdgDau);
771 Bool_t fillHisto=fDoImpPar;
773 AliAODMCParticle *partD = (AliAODMCParticle*)arrayMC->At(labD);
774 Int_t code=partD->GetPdgCode();
775 if(CheckOrigin(arrayMC,partD)==5) isPrimary=kFALSE;
776 if(code<0 && iHyp==0) fillHisto=kFALSE;
777 if(code>0 && iHyp==1) fillHisto=kFALSE;
780 trueImpParXY=AliVertexingHFUtils::GetTrueImpactParameterDplus(mcHeader,arrayMC,partD)*10000.;
781 }else if(fPdgMeson==421){
782 trueImpParXY=AliVertexingHFUtils::GetTrueImpactParameterDzero(mcHeader,arrayMC,partD)*10000.;
783 }else if(fPdgMeson==413){
784 trueImpParXY=0.; /// FIXME
786 Double_t arrayForSparseTrue[5]={invMass,ptCand,trueImpParXY,dlen,multForCand};
787 if(fillHisto && passAllCuts){
788 fHistMassPtImpPar[2]->Fill(arrayForSparse);
789 fHistMassPtImpPar[3]->Fill(arrayForSparseTrue);
792 if(fillHisto && passAllCuts) fHistMassPtImpPar[1]->Fill(arrayForSparse);
795 if(fillHisto && passAllCuts)fHistMassPtImpPar[4]->Fill(arrayForSparse);
798 if(TMath::Abs(labD)==fPdgMeson && fMCOption==2) continue;
799 if(TMath::Abs(labD)!=fPdgMeson && fMCOption==1) continue;
804 if(iHyp==0 && !(passTopolCuts&1)) continue; // candidate not passing as D0
805 if(iHyp==1 && !(passTopolCuts&2)) continue; // candidate not passing as D0bar
808 fPtVsMassVsMultNoPid->Fill(multForCand,invMass,ptCand);
811 if(iHyp==0 && !(passAllCuts&1)) continue; // candidate not passing as D0
812 if(iHyp==1 && !(passAllCuts&2)) continue; // candidate not passing as D0bar
815 aveMult+=multForCand;
817 fPtVsMassVsMult->Fill(multForCand,invMass,ptCand,nchWeight);
818 fPtVsMassVsMultUncorr->Fill(countTreta1,invMass,ptCand,nchWeight);
819 // Add separation between part antipart
821 if(d->GetCharge()>0) fPtVsMassVsMultPart->Fill(multForCand,invMass,ptCand,nchWeight);
822 else fPtVsMassVsMultAntiPart->Fill(multForCand,invMass,ptCand,nchWeight);
823 }else if(fPdgMeson==421){
824 if(passAllCuts&1) fPtVsMassVsMultPart->Fill(multForCand,invMass,ptCand,nchWeight);
825 if(passAllCuts&2) fPtVsMassVsMultAntiPart->Fill(multForCand,invMass,ptCand,nchWeight);
826 }else if(fPdgMeson==413){
827 if(d->GetCharge()>0) fPtVsMassVsMultPart->Fill(multForCand,invMass,ptCand,nchWeight);
828 else fPtVsMassVsMultAntiPart->Fill(multForCand,invMass,ptCand,nchWeight);
832 fHistMassPtImpPar[0]->Fill(arrayForSparse);
839 if(fSubtractTrackletsFromDau && nSelCand>0){
841 fCounter->StoreEvent(aod,fRDCutsAnalysis,fReadMC,(Int_t)(aveMult+0.5001));
843 fCounter->StoreEvent(aod,fRDCutsAnalysis,fReadMC,(Int_t)countCorr);
847 fCounter->StoreCandidates(aod,nSelectedNoPID,kTRUE);
848 fCounter->StoreCandidates(aod,nSelectedPID,kFALSE);
849 fHistNtrUnCorrEvSel->Fill(countMult,nchWeight);
850 fHistNtrCorrEvSel->Fill(countCorr,nchWeight);
852 fHistNtrUnCorrEvWithCand->Fill(countMult,nchWeight);
853 fHistNtrCorrEvWithCand->Fill(countCorr,nchWeight);
855 if(nSelectedInMassPeak>0) {
856 fHistNtrUnCorrEvWithD->Fill(countMult,nchWeight);
857 fHistNtrCorrEvWithD->Fill(countCorr,nchWeight);
861 PostData(2,fListCuts);
862 PostData(3,fOutputCounters);
867 //________________________________________________________________________
868 void AliAnalysisTaskSEDvsMultiplicity::CreateImpactParameterHistos(){
869 // Histos for impact paramter study
870 // mass . pt , impact parameter , decay length , multiplicity
872 Int_t nbins[5]={fNMassBins,200,fNImpParBins,50,100};
873 Double_t xmin[5]={fLowmasslimit,0.,fLowerImpPar,0.,0.};
874 Double_t xmax[5]={fUpmasslimit,20.,fHigherImpPar,1.,100.};
876 fHistMassPtImpPar[0]=new THnSparseF("hMassPtImpParAll",
877 "Mass vs. pt vs.imppar - All",
879 fHistMassPtImpPar[1]=new THnSparseF("hMassPtImpParPrompt",
880 "Mass vs. pt vs.imppar - promptD",
882 fHistMassPtImpPar[2]=new THnSparseF("hMassPtImpParBfeed",
883 "Mass vs. pt vs.imppar - DfromB",
885 fHistMassPtImpPar[3]=new THnSparseF("hMassPtImpParTrueBfeed",
886 "Mass vs. pt vs.true imppar -DfromB",
888 fHistMassPtImpPar[4]=new THnSparseF("hMassPtImpParBkg",
889 "Mass vs. pt vs.imppar - backgr.",
891 for(Int_t i=0; i<5;i++){
892 fOutput->Add(fHistMassPtImpPar[i]);
896 //________________________________________________________________________
897 void AliAnalysisTaskSEDvsMultiplicity::Terminate(Option_t */*option*/)
899 // Terminate analysis
901 if(fDebug > 1) printf("AnalysisTaskSEDvsMultiplicity: Terminate() \n");
903 fOutput = dynamic_cast<TList*> (GetOutputData(1));
905 printf("ERROR: fOutput not available\n");
909 fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
911 printf("ERROR: fHistNEvents not available\n");
914 printf("Number of Analyzed Events = %d\n",(Int_t)fHistNEvents->GetBinContent(3));
918 //_________________________________________________________________________________________________
919 Int_t AliAnalysisTaskSEDvsMultiplicity::CheckOrigin(TClonesArray* arrayMC, AliAODMCParticle *mcPartCandidate) const {
921 // checking whether the mother of the particles come from a charm or a bottom quark
926 mother = mcPartCandidate->GetMother();
928 Int_t abspdgGranma =0;
929 Bool_t isFromB=kFALSE;
930 Bool_t isQuarkFound=kFALSE;
933 AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
935 pdgGranma = mcGranma->GetPdgCode();
936 abspdgGranma = TMath::Abs(pdgGranma);
937 if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
940 if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
941 mother = mcGranma->GetMother();
943 AliError("Failed casting the mother particle!");
948 if(isFromB) return 5;
954 //____________________________________________________________________________
955 TProfile* AliAnalysisTaskSEDvsMultiplicity::GetEstimatorHistogram(const AliVEvent* event){
956 // Get Estimator Histogram from period event->GetRunNumber();
958 // If you select SPD tracklets in |eta|<1 you should use type == 1
961 Int_t runNo = event->GetRunNumber();
962 Int_t period = -1; // pp: 0-LHC10b, 1-LHC10c, 2-LHC10d, 3-LHC10e
963 // pPb: 0-LHC13b, 1-LHC13c
965 if (runNo>195343 && runNo<195484) period = 0;
966 if (runNo>195528 && runNo<195678) period = 1;
967 if (period < 0 || period > 1) return 0;
970 if(runNo>114930 && runNo<117223) period = 0;
971 if(runNo>119158 && runNo<120830) period = 1;
972 if(runNo>122373 && runNo<126438) period = 2;
973 if(runNo>127711 && runNo<130841) period = 3;
974 if(period<0 || period>3) return 0;
978 return fMultEstimatorAvg[period];
981 //__________________________________________________________________________________________________
982 void AliAnalysisTaskSEDvsMultiplicity::CreateMeasuredNchHisto(){
983 // creates historgam with measured multiplcity distribution in pp 7 TeV collisions (from Eur. Phys. J. C (2010) 68: 345–354)
985 // for Nch > 70 the points were obtainedwith a double NBD distribution
986 // TF1 *fit1 = new TF1("fit1","[0]*(TMath::Gamma(x+[1])/(TMath::Gamma(x+1)*TMath::Gamma([1])))*(TMath::Power(([2]/[1]),x))*(TMath::Power((1+([2]/[1])),-x-[1]))"); fit1->SetParameter(0,1.);// normalization constant
987 // fit1->SetParameter(1,1.63); // k parameter
988 // fit1->SetParameter(2,12.8); // mean multiplicity
989 Double_t nchbins[82]={0.50,1.50,2.50,3.50,4.50,5.50,6.50,7.50,8.50,9.50,
990 10.50,11.50,12.50,13.50,14.50,15.50,16.50,17.50,18.50,19.50,
991 20.50,21.50,22.50,23.50,24.50,25.50,26.50,27.50,28.50,29.50,
992 30.50,31.50,32.50,33.50,34.50,35.50,36.50,37.50,38.50,39.50,
993 40.50,41.50,42.50,43.50,44.50,45.50,46.50,47.50,48.50,49.50,
994 50.50,51.50,52.50,53.50,54.50,55.50,56.50,57.50,58.50,59.50,
995 60.50,62.50,64.50,66.50,68.50,70.50,72.50,74.50,76.50,78.50,
996 80.50,82.50,84.50,86.50,88.50,90.50,92.50,94.50,96.50,98.50,
998 Double_t pch[81]={0.062011,0.072943,0.070771,0.067245,0.062834,0.057383,0.051499,0.04591,0.041109,0.036954,
999 0.03359,0.030729,0.028539,0.026575,0.024653,0.0229,0.021325,0.019768,0.018561,0.017187,
1000 0.01604,0.014836,0.013726,0.012576,0.011481,0.010393,0.009502,0.008776,0.008024,0.007452,
1001 0.006851,0.006428,0.00594,0.005515,0.005102,0.00469,0.004162,0.003811,0.003389,0.003071,
1002 0.002708,0.002422,0.002184,0.001968,0.00186,0.00165,0.001577,0.001387,0.001254,0.001118,
1003 0.001037,0.000942,0.000823,0.000736,0.000654,0.000579,0.000512,0.00049,0.00045,0.000355,
1004 0.000296,0.000265,0.000193,0.00016,0.000126,0.0000851, 0.0000676,0.0000537,0.0000426, 0.0000338,
1005 0.0000268,0.0000213,0.0000166,0.0000133,0.0000106,0.00000837,0.00000662, 0.00000524,0.00000414, 0.00000327,
1008 if(fHistoMeasNch) delete fHistoMeasNch;
1009 fHistoMeasNch=new TH1F("hMeaseNch","",81,nchbins);
1010 for(Int_t i=0; i<81; i++){
1011 fHistoMeasNch->SetBinContent(i+1,pch[i]);
1012 fHistoMeasNch->SetBinError(i+1,0.);
1016 //__________________________________________________________________________________________________
1017 void AliAnalysisTaskSEDvsMultiplicity::FillMCMassHistos(TClonesArray *arrayMC, Int_t labD, Int_t countMult,Double_t nchWeight)
1020 // Function to fill the true MC signal
1024 AliAODMCParticle *partD = (AliAODMCParticle*)arrayMC->At(labD);
1025 Double_t mass = partD->M();
1026 Double_t pt = partD->Pt();
1027 fPtVsMassVsMultMC->Fill(countMult,mass,pt,nchWeight);