]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGHF/hfe/AliHFEmcQA.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Changes to corrections methods in CorrFctnDEtaDPhiCorrections, correction factor...
[u/mrichter/AliRoot.git] / PWGHF / hfe / AliHFEmcQA.cxx
CommitLineData
259c3296 1/**************************************************************************
2 * Copyright(c) 1998-1999, 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 **************************************************************************/
50685501 15//
16// QA class of Heavy Flavor quark and fragmeted/decayed particles
17// -Check kinematics of Heavy Quarks/hadrons, and decayed leptons
18// pT, rapidity
19// decay lepton kinematics w/wo acceptance
20// heavy hadron decay length, electron pT fraction carried from decay
21// -Check yield of Heavy Quarks/hadrons
22// Number of produced heavy quark
23// Number of produced hadron of given pdg code
24//
25//
26// Authors:
27// MinJung Kweon <minjung@physi.uni-heidelberg.de>
28//
259c3296 29
259c3296 30#include <TH2F.h>
259c3296 31#include <TH1F.h>
32#include <TH2F.h>
70da6c5a 33#include <TList.h>
259c3296 34#include <TParticle.h>
a8ef1999 35#include "TTreeStream.h"
259c3296 36
37#include <AliLog.h>
faee3b18 38#include <AliMCEvent.h>
9bcfd1ab 39#include <AliGenEventHeader.h>
0792aa82 40#include <AliAODMCParticle.h>
c2690925 41#include <AliStack.h>
259c3296 42
43#include "AliHFEmcQA.h"
70da6c5a 44#include "AliHFEtools.h"
c2690925 45#include "AliHFEcollection.h"
259c3296 46
259c3296 47ClassImp(AliHFEmcQA)
48
49//_______________________________________________________________________________________________
76d0b522 50AliHFEmcQA::AliHFEmcQA() :
51fMCEvent(NULL)
52 ,fMCHeader(NULL)
53 ,fMCArray(NULL)
54 ,fQAhistos(NULL)
55 ,fMCQACollection(NULL)
56 ,fNparents(0)
57 ,fCentrality(0)
58 ,fPerCentrality(-1)
59 ,fIsPbPb(kFALSE)
60 ,fIsppMultiBin(kFALSE)
61 ,fContainerStep(0)
62 ,fIsDebugStreamerON(kFALSE)
63 ,fRecPt(-999)
64 ,fRecEta(-999)
65 ,fRecPhi(-999)
66 ,fLyrhit(0)
67 ,fLyrstat(0)
68 ,fHfeImpactR(-999)
69 ,fHfeImpactnsigmaR(-999)
70 ,fTreeStream(NULL)
259c3296 71{
76d0b522 72 // Default constructor
e088e9f0 73 for(Int_t mom = 0; mom < 9; mom++){
74 fhD[mom] = NULL;
e088e9f0 75 }
bf892a6a 76 for(Int_t mom = 0; mom < 50; mom++){
77 fHeavyQuark[mom] = NULL;
78 }
79 for(Int_t mom = 0; mom < 2; mom++){
80 fIsHeavy[mom] = 0;
81 }
8c1c76e9 82 memset(fElecBackgroundFactor, 0, sizeof(Double_t) * kElecBgSpecies * kBgPtBins * kCentBins * kBgLevels);
9250ffbf 83 memset(fBinLimit, 0, sizeof(Double_t) * (kBgPtBins+1));
259c3296 84}
85
86//_______________________________________________________________________________________________
87AliHFEmcQA::AliHFEmcQA(const AliHFEmcQA&p):
76d0b522 88 TObject(p)
89 ,fMCEvent(NULL)
90 ,fMCHeader(NULL)
91 ,fMCArray(NULL)
92 ,fQAhistos(p.fQAhistos)
93 ,fMCQACollection(p.fMCQACollection)
94 ,fNparents(p.fNparents)
95 ,fCentrality(0)
96 ,fPerCentrality(-1)
97 ,fIsPbPb(kFALSE)
98 ,fIsppMultiBin(kFALSE)
99 ,fContainerStep(0)
100 ,fIsDebugStreamerON(kFALSE)
101 ,fRecPt(-999)
102 ,fRecEta(-999)
103 ,fRecPhi(-999)
104 ,fLyrhit(0)
105 ,fLyrstat(0)
106 ,fHfeImpactR(0)
107 ,fHfeImpactnsigmaR(0)
108 ,fTreeStream(NULL)
259c3296 109{
76d0b522 110 // Copy constructor
e088e9f0 111 for(Int_t mom = 0; mom < 9; mom++){
112 fhD[mom] = NULL;
e088e9f0 113 }
bf892a6a 114 for(Int_t mom = 0; mom < 50; mom++){
115 fHeavyQuark[mom] = NULL;
116 }
117 for(Int_t mom = 0; mom < 2; mom++){
118 fIsHeavy[mom] = 0;
119 }
8c1c76e9 120 memset(fElecBackgroundFactor, 0, sizeof(Double_t) * kElecBgSpecies * kBgPtBins * kCentBins * kBgLevels);
9250ffbf 121 memset(fBinLimit, 0, sizeof(Double_t) * (kBgPtBins+1));
259c3296 122}
259c3296 123//_______________________________________________________________________________________________
124AliHFEmcQA&
125AliHFEmcQA::operator=(const AliHFEmcQA &)
126{
127 // Assignment operator
76d0b522 128
259c3296 129 AliInfo("Not yet implemented.");
130 return *this;
131}
132
133//_______________________________________________________________________________________________
134AliHFEmcQA::~AliHFEmcQA()
135{
76d0b522 136 // Destructor
137
138 if(fTreeStream && fIsDebugStreamerON) delete fTreeStream;
139 AliInfo("Analysis Done.");
259c3296 140}
259c3296 141//_______________________________________________________________________________________________
50685501 142void AliHFEmcQA::PostAnalyze() const
259c3296 143{
76d0b522 144 //
145 // Post analysis
146 //
259c3296 147}
9250ffbf 148//_______________________________________________________________________________________________
149void AliHFEmcQA::SetBackgroundWeightFactor(Double_t *elecBackgroundFactor, Double_t *binLimit)
150{
76d0b522 151 //
152 // copy background weighting factors into data member
153 //
8c1c76e9 154
155 memcpy(fElecBackgroundFactor,elecBackgroundFactor,sizeof(Double_t) * kElecBgSpecies * kBgPtBins * kCentBins * kBgLevels);
156 memcpy(fBinLimit,binLimit,sizeof(Double_t) * (kBgPtBins+1));
9250ffbf 157}
e3fc062d 158//__________________________________________
159void AliHFEmcQA::CreatDefaultHistograms(TList * const qaList)
160{
161 //
162 // make default histograms
163 //
164
165 if(!qaList) return;
76d0b522 166
e3fc062d 167 fQAhistos = qaList;
168 fQAhistos->SetName("MCqa");
76d0b522 169
a8ef1999 170 CreateHistograms(AliHFEmcQA::kCharm); // create histograms for charm
171 CreateHistograms(AliHFEmcQA::kBeauty); // create histograms for beauty
172 CreateHistograms(AliHFEmcQA::kOthers); // create histograms for beauty
76d0b522 173
174 // prepare 2D(pt vs Y) histogram for D spectra, we consider following 9 particles
e17c1f86 175 const Int_t nbspecies = 9;
176 TString kDspecies[nbspecies];
177 kDspecies[0]="411"; //D+
178 kDspecies[1]="421"; //D0
179 kDspecies[2]="431"; //Ds+
180 kDspecies[3]="4122"; //Lambdac+
181 kDspecies[4]="4132"; //Ksic0
182 kDspecies[5]="4232"; //Ksic+
183 kDspecies[6]="4332"; //OmegaC0
184 kDspecies[7]="413"; //D*(2010)+
185 kDspecies[8]="423"; //D*(2007)0
ccc37cdc 186
187 const Double_t kPtbound[2] = {0.1, 20.}; //bin taken for considering inclusive e analysis binning
c2690925 188 Int_t iBin[2];
e17c1f86 189 iBin[0] = 44; // bins in pt for log binning
190 iBin[1] = 23; // bins in pt for pi0 measurement binning
ccc37cdc 191 Double_t* binEdges[1];
192 binEdges[0] = AliHFEtools::MakeLogarithmicBinning(iBin[0], kPtbound[0], kPtbound[1]);
e97c2edf 193
194 // bin size is chosen to consider ALICE D measurement
e17c1f86 195 const Int_t nptbins = 15;
196 const Int_t nybins = 9;
197 Double_t xbins[nptbins+1]={0,0.5,1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,12,16,24,32,40,50}; //pt binning for the final 7 TeV D measurement
198 Double_t ybins[nybins+1]={-7.5,-1.0,-0.9,-0.8,-0.5,0.5,0.8,0.9,1.0,7.5}; // y binning
e97c2edf 199 TString hname;
e17c1f86 200 for (Int_t iDmeson=0; iDmeson<nbspecies; iDmeson++){
e97c2edf 201 hname = "Dmeson"+kDspecies[iDmeson];
e17c1f86 202 fhD[iDmeson] = new TH2F(hname,hname+";p_{T} (GeV/c)",nptbins,xbins,nybins,ybins);
e97c2edf 203 if(fQAhistos) fQAhistos->Add(fhD[iDmeson]);
204 }
205
c2690925 206 const Double_t kPtRange[24] = {0.,0.3,0.4,0.5,0.6,0.8,1.,1.2,1.4,1.6,1.8,2.,2.2,2.4,2.6,2.8,3.,3.5,4.,5.,6.,7.,20.,30.}; // to cope with Ana's bin
207
8c1c76e9 208 Int_t kNcent;
209 if(fIsPbPb) kNcent=11;
210 else
211 {
212 if(fIsppMultiBin) kNcent=8;
213 else kNcent = 1;
214 }
215
c2690925 216 fMCQACollection = new AliHFEcollection("TaskMCQA", "MC QA histos for meason pt spectra");
8c1c76e9 217
218 for(Int_t centbin=0; centbin<kNcent; centbin++)
219 {
220 fMCQACollection->CreateTH1Farray(Form("pionspectra_centrbin%i",centbin), "pion yields: MC p_{t} ", iBin[1],kPtRange);
221 fMCQACollection->CreateTH1Farray(Form("etaspectra_centrbin%i",centbin), "eta yields: MC p_{t} ", iBin[1],kPtRange);
222 fMCQACollection->CreateTH1Farray(Form("omegaspectra_centrbin%i",centbin), "omega yields: MC p_{t} ", iBin[1],kPtRange);
223 fMCQACollection->CreateTH1Farray(Form("phispectra_centrbin%i",centbin), "phi yields: MC p_{t} ", iBin[1],kPtRange);
224 fMCQACollection->CreateTH1Farray(Form("etapspectra_centrbin%i",centbin), "etap yields: MC p_{t} ", iBin[1],kPtRange);
225 fMCQACollection->CreateTH1Farray(Form("rhospectra_centrbin%i",centbin), "rho yields: MC p_{t} ", iBin[1],kPtRange);
226
afb48e1d 227 fMCQACollection->CreateTH1Farray(Form("pionspectrapri_centrbin%i",centbin), "pion yields: MC p_{t} ", iBin[1],kPtRange);
228 fMCQACollection->CreateTH1Farray(Form("etaspectrapri_centrbin%i",centbin), "eta yields: MC p_{t} ", iBin[1],kPtRange);
229 fMCQACollection->CreateTH1Farray(Form("pionspectrasec_centrbin%i",centbin), "pion yields: MC p_{t} ", iBin[1],kPtRange);
230 fMCQACollection->CreateTH1Farray(Form("etaspectrasec_centrbin%i",centbin), "eta yields: MC p_{t} ", iBin[1],kPtRange);
231 fMCQACollection->CreateTH1Farray(Form("pionspectraLogpri_centrbin%i",centbin), "pion yields: MC p_{t} ", iBin[1],kPtRange);
232 fMCQACollection->CreateTH1Farray(Form("etaspectraLogpri_centrbin%i",centbin), "eta yields: MC p_{t} ", iBin[1],kPtRange);
233 fMCQACollection->CreateTH1Farray(Form("pionspectraLogsec_centrbin%i",centbin), "pion yields: MC p_{t} ", iBin[1],kPtRange);
234 fMCQACollection->CreateTH1Farray(Form("etaspectraLogsec_centrbin%i",centbin), "eta yields: MC p_{t} ", iBin[1],kPtRange);
235
8c1c76e9 236 fMCQACollection->CreateTH1F(Form("pionspectraLog_centrbin%i",centbin), "pion yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
237 fMCQACollection->CreateTH1F(Form("etaspectraLog_centrbin%i",centbin), "eta yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
238 fMCQACollection->CreateTH1F(Form("omegaspectraLog_centrbin%i",centbin), "omega yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
239 fMCQACollection->CreateTH1F(Form("phispectraLog_centrbin%i",centbin), "phi yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
240 fMCQACollection->CreateTH1F(Form("etapspectraLog_centrbin%i",centbin), "etap yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
241 fMCQACollection->CreateTH1F(Form("rhospectraLog_centrbin%i",centbin), "rho yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
242 fMCQACollection->CreateTH1F(Form("kaonspectraLog_centrbin%i",centbin), "kaon yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
243 fMCQACollection->CreateTH1F(Form("k0LspectraLog_centrbin%i",centbin), "k0L yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
244
245 fMCQACollection->CreateTH1F(Form("piondaughters_centrbin%i",centbin), "pion yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
246 fMCQACollection->CreateTH1F(Form("etadaughters_centrbin%i",centbin), "eta yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
247 fMCQACollection->CreateTH1F(Form("omegadaughters_centrbin%i",centbin), "omega yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
248 fMCQACollection->CreateTH1F(Form("phidaughters_centrbin%i",centbin), "phi yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
249 fMCQACollection->CreateTH1F(Form("etapdaughters_centrbin%i",centbin), "etap yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
250 fMCQACollection->CreateTH1F(Form("rhodaughters_centrbin%i",centbin), "rho yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
251 }
c2690925 252
253 fQAhistos->Add(fMCQACollection->GetList());
254
a8ef1999 255 if(!fTreeStream && fIsDebugStreamerON){
256 fTreeStream = new TTreeSRedirector(Form("HFEmcqadebugTree%s.root", GetName()));
257 }
258
e3fc062d 259}
260
259c3296 261//__________________________________________
5a516e0a 262void AliHFEmcQA::CreateHistograms(Int_t kquark)
259c3296 263{
264 // create histograms
265
faee3b18 266 if (!(kquark == kCharm || kquark == kBeauty || kquark == kOthers)) {
259c3296 267 AliDebug(1, "This task is only for heavy quark QA, return\n");
268 return;
269 }
dbe3abbe 270 Int_t iq = kquark - kCharm;
259c3296 271
272 TString kqTypeLabel[fgkqType];
dbe3abbe 273 if (kquark == kCharm){
274 kqTypeLabel[kQuark]="c";
275 kqTypeLabel[kantiQuark]="cbar";
276 kqTypeLabel[kHadron]="cHadron";
277 kqTypeLabel[keHadron]="ceHadron";
278 kqTypeLabel[kDeHadron]="nullHadron";
279 kqTypeLabel[kElectron]="ce";
280 kqTypeLabel[kElectron2nd]="nulle";
281 } else if (kquark == kBeauty){
282 kqTypeLabel[kQuark]="b";
283 kqTypeLabel[kantiQuark]="bbar";
284 kqTypeLabel[kHadron]="bHadron";
285 kqTypeLabel[keHadron]="beHadron";
286 kqTypeLabel[kDeHadron]="bDeHadron";
287 kqTypeLabel[kElectron]="be";
288 kqTypeLabel[kElectron2nd]="bce";
faee3b18 289 } else if (kquark == kOthers){
290 kqTypeLabel[kGamma-4]="gammae";
291 kqTypeLabel[kPi0-4]="pi0e";
292 kqTypeLabel[kElse-4]="elsee";
293 kqTypeLabel[kMisID-4]="miside";
259c3296 294 }
3a72645a 295
a8ef1999 296 TString kqEtaRangeLabel[fgkEtaRanges];
297 kqEtaRangeLabel[0] = "mcqa_";
298 kqEtaRangeLabel[1] = "mcqa_barrel_";
299 kqEtaRangeLabel[2] = "mcqa_unitY_";
300
3a72645a 301 const Double_t kPtbound[2] = {0.1, 20.}; //bin taken for considering inclusive e analysis binning
e17c1f86 302 const Int_t nptbinning1 = 35;
303 Int_t iBin[2];
3a72645a 304 iBin[0] = 44; // bins in pt
e17c1f86 305 iBin[1] = nptbinning1; // bins in pt
faee3b18 306 Double_t* binEdges[1];
307 binEdges[0] = AliHFEtools::MakeLogarithmicBinning(iBin[0], kPtbound[0], kPtbound[1]);
c2690925 308
e17c1f86 309 // new binning for final electron analysis
310 const Double_t kPtbinning1[nptbinning1+1] = {0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4, 1.5, 1.75, 2., 2.25, 2.5, 2.75, 3., 3.5, 4., 4.5, 5., 5.5, 6., 7., 8., 10., 12., 14., 16., 18., 20.};
c2690925 311
e17c1f86 312 const Int_t ndptbins = 500;
313 Double_t xcorrbin[ndptbins+1];
314 for (int icorrbin = 0; icorrbin< ndptbins+1; icorrbin++){
c2690925 315 xcorrbin[icorrbin]=icorrbin*0.1;
316 }
317
11ff28c5 318 Int_t fCentrmax = 0;
319 if(!fIsPbPb&&!fIsppMultiBin) fCentrmax=0;
320 else fCentrmax = kCentBins;
259c3296 321 TString hname;
faee3b18 322 if(kquark == kOthers){
a8ef1999 323 for (Int_t icut = 0; icut < fgkEtaRanges; icut++ ){
11ff28c5 324 for (Int_t iqType = 0; iqType < 4; iqType++ ){
325 for(Int_t icentr = 0; icentr < (fCentrmax+1); icentr++)
326 {
327 hname = kqEtaRangeLabel[icut]+"Pt_"+kqTypeLabel[iqType];
328 fHist[iq][iqType][icut][icentr].fPt = new TH1F(Form("%sCentr_%i",hname.Data(),icentr),Form("%sCentr_%i;p_{T} (GeV/c)",hname.Data(),icentr),60,0.25,30.25);
329 hname = kqEtaRangeLabel[icut]+"Y_"+kqTypeLabel[iqType];
330 fHist[iq][iqType][icut][icentr].fY = new TH1F(Form("%sCentr_%i",hname.Data(),icentr),Form("%sCentr_%i;y",hname.Data(),icentr),150,-7.5,7.5);
331 hname = kqEtaRangeLabel[icut]+"Eta_"+kqTypeLabel[iqType];
332 fHist[iq][iqType][icut][icentr].fEta = new TH1F(Form("%sCentr_%i",hname.Data(),icentr),Form("%sCentr_%i;eta",hname.Data(),icentr),150,-7.5,7.5);
333 // Fill List
334 if(fQAhistos) fHist[iq][iqType][icut][icentr].FillList(fQAhistos);
335 }
336 }
a8ef1999 337 }
338 return;
faee3b18 339 }
a8ef1999 340 for (Int_t icut = 0; icut < fgkEtaRanges; icut++ ){
341 for (Int_t iqType = 0; iqType < fgkqType; iqType++ ){
11ff28c5 342 if (iqType < keHadron && icut > 0) continue; // don't duplicate histogram for quark and hadron
343 for(Int_t icentr = 0; icentr<(fCentrmax+1); icentr++)
344 {
345 hname = kqEtaRangeLabel[icut]+"PdgCode_"+kqTypeLabel[iqType];
346 fHist[iq][iqType][icut][icentr].fPdgCode = new TH1F(Form("%sCentr_%i",hname.Data(),icentr),Form("%sCentr_%i;PdgCode",hname.Data(),icentr),20001,-10000.5,10000.5);
347 hname = kqEtaRangeLabel[icut]+"Pt_"+kqTypeLabel[iqType];
348 fHist[iq][iqType][icut][icentr].fPt = new TH1F(Form("%sCentr_%i",hname.Data(),icentr),Form("%sCentr_%i;p_{T} (GeV/c)",hname.Data(),icentr),iBin[1],kPtbinning1); // new binning
349 hname = kqEtaRangeLabel[icut]+"Y_"+kqTypeLabel[iqType];
350 fHist[iq][iqType][icut][icentr].fY = new TH1F(Form("%sCentr_%i",hname.Data(),icentr),Form("%sCentr_%i;y",hname.Data(),icentr),150,-7.5,7.5);
351 hname = kqEtaRangeLabel[icut]+"Eta_"+kqTypeLabel[iqType];
352 fHist[iq][iqType][icut][icentr].fEta = new TH1F(Form("%sCentr_%i",hname.Data(),icentr),Form("%sCentr_%i;eta",hname.Data(),icentr),150,-7.5,7.5);
353 // Fill List
354 if(fQAhistos) fHist[iq][iqType][icut][icentr].FillList(fQAhistos);
355 }
a8ef1999 356 }
259c3296 357 }
358
a8ef1999 359 for (Int_t icut = 0; icut < fgkEtaRanges; icut++ ){
360 hname = kqEtaRangeLabel[icut]+"PtCorr_"+kqTypeLabel[kQuark];
361 fHistComm[iq][icut].fPtCorr = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1);
362 hname = kqEtaRangeLabel[icut]+"PtCorrDp_"+kqTypeLabel[kQuark];
363 fHistComm[iq][icut].fPtCorrDp= new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1);
364 hname = kqEtaRangeLabel[icut]+"PtCorrD0_"+kqTypeLabel[kQuark];
365 fHistComm[iq][icut].fPtCorrD0 = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1);
366 hname = kqEtaRangeLabel[icut]+"PtCorrDrest_"+kqTypeLabel[kQuark];
367 fHistComm[iq][icut].fPtCorrDrest = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1);
cedf0381 368
a8ef1999 369 hname = kqEtaRangeLabel[icut]+"ePtRatio_"+kqTypeLabel[kQuark];
370 fHistComm[iq][icut].fePtRatio = new TH2F(hname,hname+";p_{T} (GeV/c);momentum fraction",200,0,20,100,0,1);
371 hname = kqEtaRangeLabel[icut]+"DePtRatio_"+kqTypeLabel[kQuark];
372 fHistComm[iq][icut].fDePtRatio = new TH2F(hname,hname+";p_{T} (GeV/c);momentum fraction",100,0,20,100,0,1);
373 hname = kqEtaRangeLabel[icut]+"eDistance_"+kqTypeLabel[kQuark];
374 fHistComm[iq][icut].feDistance= new TH2F(hname,hname+";p_{T} (GeV/c);distance (cm)",100,0,20,200,0,2);
375 hname = kqEtaRangeLabel[icut]+"DeDistance_"+kqTypeLabel[kQuark];
376 fHistComm[iq][icut].fDeDistance= new TH2F(hname,hname+";p_{T} (GeV/c);distance (cm)",100,0,20,200,0,2);
cedf0381 377
378 if(icut <1){
379 hname = kqEtaRangeLabel[icut]+"PtCorrDinein_"+kqTypeLabel[kQuark];
380 fHistComm[iq][icut].fPtCorrDinein = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
381 hname = kqEtaRangeLabel[icut]+"PtCorrDineout_"+kqTypeLabel[kQuark];
382 fHistComm[iq][icut].fPtCorrDineout = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
383 hname = kqEtaRangeLabel[icut]+"PtCorrDoutein_"+kqTypeLabel[kQuark];
384 fHistComm[iq][icut].fPtCorrDoutein = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
385 hname = kqEtaRangeLabel[icut]+"PtCorrDouteout_"+kqTypeLabel[kQuark];
386 fHistComm[iq][icut].fPtCorrDouteout = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
387
388 hname = kqEtaRangeLabel[icut]+"PtCorrDpDinein_"+kqTypeLabel[kQuark];
389 fHistComm[iq][icut].fPtCorrDpDinein = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
390 hname = kqEtaRangeLabel[icut]+"PtCorrDpDineout_"+kqTypeLabel[kQuark];
391 fHistComm[iq][icut].fPtCorrDpDineout = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
392 hname = kqEtaRangeLabel[icut]+"PtCorrDpDoutein_"+kqTypeLabel[kQuark];
393 fHistComm[iq][icut].fPtCorrDpDoutein = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
394 hname = kqEtaRangeLabel[icut]+"PtCorrDpDouteout_"+kqTypeLabel[kQuark];
395 fHistComm[iq][icut].fPtCorrDpDouteout = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
396
397 hname = kqEtaRangeLabel[icut]+"PtCorrD0Dinein_"+kqTypeLabel[kQuark];
398 fHistComm[iq][icut].fPtCorrD0Dinein = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
399 hname = kqEtaRangeLabel[icut]+"PtCorrD0Dineout_"+kqTypeLabel[kQuark];
400 fHistComm[iq][icut].fPtCorrD0Dineout = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
401 hname = kqEtaRangeLabel[icut]+"PtCorrD0Doutein_"+kqTypeLabel[kQuark];
402 fHistComm[iq][icut].fPtCorrD0Doutein = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
403 hname = kqEtaRangeLabel[icut]+"PtCorrD0Douteout_"+kqTypeLabel[kQuark];
404 fHistComm[iq][icut].fPtCorrD0Douteout = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
405
406 hname = kqEtaRangeLabel[icut]+"PtCorrDrestDinein_"+kqTypeLabel[kQuark];
407 fHistComm[iq][icut].fPtCorrDrestDinein = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
408 hname = kqEtaRangeLabel[icut]+"PtCorrDrestDineout_"+kqTypeLabel[kQuark];
409 fHistComm[iq][icut].fPtCorrDrestDineout = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
410 hname = kqEtaRangeLabel[icut]+"PtCorrDrestDoutein_"+kqTypeLabel[kQuark];
411 fHistComm[iq][icut].fPtCorrDrestDoutein = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
412 hname = kqEtaRangeLabel[icut]+"PtCorrDrestDouteout_"+kqTypeLabel[kQuark];
413 fHistComm[iq][icut].fPtCorrDrestDouteout = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
414
415 hname = kqEtaRangeLabel[icut]+"fEtaCorrD_"+kqTypeLabel[kQuark];
416 fHistComm[iq][icut].fEtaCorrD = new TH2F(hname,hname+";D Y;e eta",200,-10,10,200,-10,10);
417 hname = kqEtaRangeLabel[icut]+"fEtaCorrDp_"+kqTypeLabel[kQuark];
418 fHistComm[iq][icut].fEtaCorrDp = new TH2F(hname,hname+";D Y;e eta",200,-10,10,200,-10,10);
419 hname = kqEtaRangeLabel[icut]+"fEtaCorrD0_"+kqTypeLabel[kQuark];
420 fHistComm[iq][icut].fEtaCorrD0 = new TH2F(hname,hname+";D Y;e eta",200,-10,10,200,-10,10);
421 hname = kqEtaRangeLabel[icut]+"fEtaCorrDrest_"+kqTypeLabel[kQuark];
422 fHistComm[iq][icut].fEtaCorrDrest = new TH2F(hname,hname+";D Y;e eta",200,-10,10,200,-10,10);
423
424 hname = kqEtaRangeLabel[icut]+"fEtaCorrGD_"+kqTypeLabel[kQuark];
425 fHistComm[iq][icut].fEtaCorrGD = new TH2F(hname,hname+";D Y;e eta",200,-10,10,200,-10,10);
426 hname = kqEtaRangeLabel[icut]+"fEtaCorrGDp_"+kqTypeLabel[kQuark];
427 fHistComm[iq][icut].fEtaCorrGDp = new TH2F(hname,hname+";D Y;e eta",200,-10,10,200,-10,10);
428 hname = kqEtaRangeLabel[icut]+"fEtaCorrGD0_"+kqTypeLabel[kQuark];
429 fHistComm[iq][icut].fEtaCorrGD0 = new TH2F(hname,hname+";D Y;e eta",200,-10,10,200,-10,10);
430 hname = kqEtaRangeLabel[icut]+"fEtaCorrGDrest_"+kqTypeLabel[kQuark];
431 fHistComm[iq][icut].fEtaCorrGDrest = new TH2F(hname,hname+";D Y;e eta",200,-10,10,200,-10,10);
432
433 hname = kqEtaRangeLabel[icut]+"fEtaCorrB_"+kqTypeLabel[kQuark];
434 fHistComm[iq][icut].fEtaCorrB = new TH2F(hname,hname+";B Y;e eta",200,-10,10,200,-10,10);
435 hname = kqEtaRangeLabel[icut]+"fEtaCorrGB_"+kqTypeLabel[kQuark];
436 fHistComm[iq][icut].fEtaCorrGB = new TH2F(hname,hname+";B Y;e eta",200,-10,10,200,-10,10);
437
438 hname = kqEtaRangeLabel[icut]+"PtCorrBinein_"+kqTypeLabel[kQuark];
439 fHistComm[iq][icut].fPtCorrBinein = new TH2F(hname,hname+";B p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
440 hname = kqEtaRangeLabel[icut]+"PtCorrBineout_"+kqTypeLabel[kQuark];
441 fHistComm[iq][icut].fPtCorrBineout = new TH2F(hname,hname+";B p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
442 hname = kqEtaRangeLabel[icut]+"PtCorrBoutein_"+kqTypeLabel[kQuark];
443 fHistComm[iq][icut].fPtCorrBoutein = new TH2F(hname,hname+";B p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
444 hname = kqEtaRangeLabel[icut]+"PtCorrBouteout_"+kqTypeLabel[kQuark];
445 fHistComm[iq][icut].fPtCorrBouteout = new TH2F(hname,hname+";B p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
446 }
a8ef1999 447 if(fQAhistos) fHistComm[iq][icut].FillList(fQAhistos);
dbe3abbe 448 }
a8ef1999 449
cedf0381 450
a8ef1999 451 hname = kqEtaRangeLabel[0]+"Nq_"+kqTypeLabel[kQuark];
452 fHistComm[iq][0].fNq = new TH1F(hname,hname,50,-0.5,49.5);
453 hname = kqEtaRangeLabel[0]+"ProcessID_"+kqTypeLabel[kQuark];
454 fHistComm[iq][0].fProcessID = new TH1F(hname,hname,21,-10.5,10.5);
e17c1f86 455
259c3296 456}
457
458//__________________________________________
459void AliHFEmcQA::Init()
460{
461 // called at begining every event
462
463 for (Int_t i=0; i<2; i++){
464 fIsHeavy[i] = 0;
465 }
466
467 fNparents = 7;
468
e17c1f86 469 fParentSelect[0][0] = 411; //D+
470 fParentSelect[0][1] = 421; //D0
471 fParentSelect[0][2] = 431; //Ds+
472 fParentSelect[0][3] = 4122; //Lambdac+
473 fParentSelect[0][4] = 4132; //Ksic0
474 fParentSelect[0][5] = 4232; //Ksic+
475 fParentSelect[0][6] = 4332; //OmegaC0
476
477 fParentSelect[1][0] = 511; //B0
478 fParentSelect[1][1] = 521; //B+
479 fParentSelect[1][2] = 531; //Bs0
480 fParentSelect[1][3] = 5122; //Lambdab0
481 fParentSelect[1][4] = 5132; //Ksib-
482 fParentSelect[1][5] = 5232; //Ksib0
483 fParentSelect[1][6] = 5332; //Omegab-
259c3296 484
11ff28c5 485
259c3296 486}
487
9250ffbf 488//__________________________________________
c2690925 489void AliHFEmcQA::GetMesonKine()
490{
491 //
492 // get meson pt spectra
493 //
494
495 AliVParticle *mctrack2 = NULL;
496 AliMCParticle *mctrack0 = NULL;
497 AliVParticle *mctrackdaugt= NULL;
498 AliMCParticle *mctrackd= NULL;
499 Int_t id1=0, id2=0;
500
8c1c76e9 501
afb48e1d 502 if(fCentrality>=11) {
503 AliWarning(Form("Centrality out of histogram array limits: %d", fCentrality));
504 return;
505 }
11ff28c5 506 if(!fIsPbPb&&!fIsppMultiBin) fCentrality=0;
8c1c76e9 507
c2690925 508 for(Int_t imc = 0; imc <fMCEvent->GetNumberOfPrimaries(); imc++){
509 if(!(mctrack2 = fMCEvent->GetTrack(imc))) continue;
510 TParticle* mcpart0 = fMCEvent->Stack()->Particle(imc);
511 if(!mcpart0) continue;
512 mctrack0 = dynamic_cast<AliMCParticle *>(mctrack2);
f9f097c0 513 if(!mctrack0) continue;
8c1c76e9 514
11ff28c5 515// if(!fIsPbPb&&!fIsppMultiBin) fCentrality=0;
8c1c76e9 516
b8bec44f 517 if(TMath::Abs(mctrack0->PdgCode()) == 111) // pi0
c2690925 518 {
519 if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
8c1c76e9 520 fMCQACollection->Fill(Form("pionspectra_centrbin%i",fCentrality),mctrack0->Pt());
521 fMCQACollection->Fill(Form("pionspectraLog_centrbin%i",fCentrality),mctrack0->Pt());
afb48e1d 522 if(imc>fMCEvent->GetNumberOfPrimaries()) {
523 fMCQACollection->Fill(Form("pionspectrasec_centrbin%i",fCentrality),mctrack0->Pt());
524 fMCQACollection->Fill(Form("pionspectraLogsec_centrbin%i",fCentrality),mctrack0->Pt());
525 }
526 else {
527 fMCQACollection->Fill(Form("pionspectrapri_centrbin%i",fCentrality),mctrack0->Pt());
528 fMCQACollection->Fill(Form("pionspectraLogpri_centrbin%i",fCentrality),mctrack0->Pt());
529 }
c2690925 530 }
531 id1=mctrack0->GetFirstDaughter();
532 id2=mctrack0->GetLastDaughter();
533 if(!((id2-id1)==2)) continue;
534 for(int idx=id1; idx<=id2; idx++){
535 if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
9250ffbf 536 if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
b8bec44f 537 if(TMath::Abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
8c1c76e9 538 fMCQACollection->Fill(Form("piondaughters_centrbin%i",fCentrality),mctrackd->Pt());
c2690925 539 }
540 }
b8bec44f 541 else if(TMath::Abs(mctrack0->PdgCode()) == 221) // eta
c2690925 542 {
543 if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
8c1c76e9 544 fMCQACollection->Fill(Form("etaspectra_centrbin%i",fCentrality),mctrack0->Pt());
545 fMCQACollection->Fill(Form("etaspectraLog_centrbin%i",fCentrality),mctrack0->Pt());
afb48e1d 546 if(imc>fMCEvent->GetNumberOfPrimaries()) {
547 fMCQACollection->Fill(Form("etaspectrasec_centrbin%i",fCentrality),mctrack0->Pt());
548 fMCQACollection->Fill(Form("etaspectraLogsec_centrbin%i",fCentrality),mctrack0->Pt());
549 }
550 else {
551 fMCQACollection->Fill(Form("etaspectrapri_centrbin%i",fCentrality),mctrack0->Pt());
552 fMCQACollection->Fill(Form("etaspectraLogpri_centrbin%i",fCentrality),mctrack0->Pt());
553 }
c2690925 554 }
555 id1=mctrack0->GetFirstDaughter();
556 id2=mctrack0->GetLastDaughter();
557 if(!((id2-id1)==2||(id2-id1)==3)) continue;
558 for(int idx=id1; idx<=id2; idx++){
559 if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
9250ffbf 560 if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
b8bec44f 561 if(TMath::Abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
8c1c76e9 562 fMCQACollection->Fill(Form("etadaughters_centrbin%i",fCentrality),mctrackd->Pt());
c2690925 563 }
564 }
b8bec44f 565 else if(TMath::Abs(mctrack0->PdgCode()) == 223) // omega
c2690925 566 {
567 if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
8c1c76e9 568 fMCQACollection->Fill(Form("omegaspectra_centrbin%i",fCentrality),mctrack0->Pt());
569 fMCQACollection->Fill(Form("omegaspectraLog_centrbin%i",fCentrality),mctrack0->Pt());
c2690925 570 }
571 id1=mctrack0->GetFirstDaughter();
572 id2=mctrack0->GetLastDaughter();
573 if(!((id2-id1)==1||(id2-id1)==2)) continue;
574 for(int idx=id1; idx<=id2; idx++){
575 if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
9250ffbf 576 if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
b8bec44f 577 if(TMath::Abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
8c1c76e9 578 fMCQACollection->Fill(Form("omegadaughters_centrbin%i",fCentrality),mctrackd->Pt());
c2690925 579 }
580 }
b8bec44f 581 else if(TMath::Abs(mctrack0->PdgCode()) == 333) // phi
c2690925 582 {
583 if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
8c1c76e9 584 fMCQACollection->Fill(Form("phispectra_centrbin%i",fCentrality),mctrack0->Pt());
585 fMCQACollection->Fill(Form("phispectraLog_centrbin%i",fCentrality),mctrack0->Pt());
c2690925 586 }
587 id1=mctrack0->GetFirstDaughter();
588 id2=mctrack0->GetLastDaughter();
589 if(!((id2-id1)==1)) continue;
590 for(int idx=id1; idx<=id2; idx++){
591 if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
9250ffbf 592 if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
b8bec44f 593 if(TMath::Abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
8c1c76e9 594 fMCQACollection->Fill(Form("phidaughters_centrbin%i",fCentrality),mctrackd->Pt());
c2690925 595 }
596 }
b8bec44f 597 else if(TMath::Abs(mctrack0->PdgCode()) == 331) // eta prime
c2690925 598 {
599 if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
8c1c76e9 600 fMCQACollection->Fill(Form("etapspectra_centrbin%i",fCentrality),mctrack0->Pt());
601 fMCQACollection->Fill(Form("etapspectraLog_centrbin%i",fCentrality),mctrack0->Pt());
c2690925 602 }
603 id1=mctrack0->GetFirstDaughter();
604 id2=mctrack0->GetLastDaughter();
605 if(!((id2-id1)==2||(id2-id1)==3)) continue;
606 for(int idx=id1; idx<=id2; idx++){
607 if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
9250ffbf 608 if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
b8bec44f 609 if(TMath::Abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
8c1c76e9 610 fMCQACollection->Fill(Form("etapdaughters_centrbin%i",fCentrality),mctrackd->Pt());
c2690925 611 }
612 }
b8bec44f 613 else if(TMath::Abs(mctrack0->PdgCode()) == 113) // rho
c2690925 614 {
615 if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
8c1c76e9 616 fMCQACollection->Fill(Form("rhospectra_centrbin%i",fCentrality),mctrack0->Pt());
617 fMCQACollection->Fill(Form("rhospectraLog_centrbin%i",fCentrality),mctrack0->Pt());
c2690925 618 }
619 id1=mctrack0->GetFirstDaughter();
620 id2=mctrack0->GetLastDaughter();
621 if(!((id2-id1)==1)) continue;
622 for(int idx=id1; idx<=id2; idx++){
623 if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
9250ffbf 624 if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
b8bec44f 625 if(TMath::Abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
8c1c76e9 626 fMCQACollection->Fill(Form("rhodaughters_centrbin%i",fCentrality),mctrackd->Pt());
627 }
628 }
b8bec44f 629 else if(TMath::Abs(mctrack0->PdgCode()) == 321) // kaon+-
8c1c76e9 630 {
631 if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
632 fMCQACollection->Fill(Form("kaonspectraLog_centrbin%i",fCentrality),mctrack0->Pt());
633 }
634 }
b8bec44f 635 else if(TMath::Abs(mctrack0->PdgCode()) == 130) // k0L
8c1c76e9 636 {
637 if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
638 fMCQACollection->Fill(Form("k0LspectraLog_centrbin%i",fCentrality),mctrack0->Pt());
c2690925 639 }
640 }
641 }
642
643}
259c3296 644//__________________________________________
5a516e0a 645void AliHFEmcQA::GetQuarkKine(TParticle *part, Int_t iTrack, Int_t kquark)
259c3296 646{
647 // get heavy quark kinematics
648
dbe3abbe 649 if (kquark != kCharm && kquark != kBeauty) {
259c3296 650 AliDebug(1, "This task is only for heavy quark QA, return\n");
651 return;
652 }
dbe3abbe 653 Int_t iq = kquark - kCharm;
259c3296 654
faee3b18 655 if (iTrack < 0 || !part) {
656 AliDebug(1, "Stack label is negative or no mcparticle, return\n");
259c3296 657 return;
658 }
659
11ff28c5 660 if(!fIsPbPb&&!fIsppMultiBin) fCentrality=0;
661
faee3b18 662 AliMCParticle *mctrack = NULL;
259c3296 663 Int_t partPdgcode = TMath::Abs(part->GetPdgCode());
664
665 // select heavy hadron or not fragmented heavy quark
666 if ( int(partPdgcode/100.)==kquark || int(partPdgcode/1000.)==kquark || (partPdgcode==kquark && (part->GetNDaughters()==0 && iTrack>5)) ){
667
668 TParticle *partMother;
669 Int_t iLabel;
670
671 if (partPdgcode == kquark){ // in case of not fragmented heavy quark
672 partMother = part;
673 iLabel = iTrack;
674 } else{ // in case of heavy hadron, start to search for mother heavy parton
675 iLabel = part->GetFirstMother();
faee3b18 676 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return;
677 if (iLabel>-1) { partMother = mctrack->Particle(); }
259c3296 678 else {
679 AliDebug(1, "Stack label is negative, return\n");
680 return;
681 }
682 }
683
684 // heavy parton selection as a mother of heavy hadron
685 // if the heavy particle comes from string which is denoted as particle status 12|12|12...12|11,[PYTHIA p.60]
686 // in this case, the mother of heavy particle can be one of the fragmented parton of the string
687 // should I make a condition that partMother should be quark or diquark? -> not necessary
b8bec44f 688 if ( TMath::Abs(partMother->GetPdgCode()) == kquark || (partMother->GetStatusCode() == 11) ){
689 //if ( TMath::Abs(partMother->GetPdgCode()) == kquark || (partMother->GetStatusCode() == 11 || partMother->GetStatusCode() == 12) ){
259c3296 690
b8bec44f 691 if ( TMath::Abs(partMother->GetPdgCode()) != kquark ){
259c3296 692 // search fragmented partons in the same string
dbe3abbe 693 Bool_t isSameString = kTRUE;
259c3296 694 for (Int_t i=1; i<fgkMaxIter; i++){
695 iLabel = iLabel - 1;
faee3b18 696 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return;
697 if (iLabel>-1) { partMother = mctrack->Particle(); }
259c3296 698 else {
699 AliDebug(1, "Stack label is negative, return\n");
700 return;
701 }
b8bec44f 702 if ( TMath::Abs(partMother->GetPdgCode()) == kquark ) break;
dbe3abbe 703 if ( partMother->GetStatusCode() != 12 ) isSameString = kFALSE;
704 if (!isSameString) return;
259c3296 705 }
706 }
707 AliDebug(1, "Can not find heavy parton of this heavy hadron in the string, return\n");
b8bec44f 708 if (TMath::Abs(partMother->GetPdgCode()) != kquark) return;
259c3296 709
710 if (fIsHeavy[iq] >= 50) return;
711 fHeavyQuark[fIsHeavy[iq]] = partMother;
712 fIsHeavy[iq]++;
713
714 // fill kinematics for heavy parton
715 if (partMother->GetPdgCode() > 0) { // quark
11ff28c5 716 fHist[iq][kQuark][0][fCentrality].fPdgCode->Fill(partMother->GetPdgCode());
717 fHist[iq][kQuark][0][fCentrality].fPt->Fill(partMother->Pt());
718 fHist[iq][kQuark][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMother));
719 fHist[iq][kQuark][0][fCentrality].fEta->Fill(partMother->Eta());
259c3296 720 } else{ // antiquark
11ff28c5 721 fHist[iq][kantiQuark][0][fCentrality].fPdgCode->Fill(partMother->GetPdgCode());
722 fHist[iq][kantiQuark][0][fCentrality].fPt->Fill(partMother->Pt());
723 fHist[iq][kantiQuark][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMother));
724 fHist[iq][kantiQuark][0][fCentrality].fEta->Fill(partMother->Eta());
259c3296 725 }
726
727 } // end of heavy parton slection loop
728
729 } // end of heavy hadron or quark selection
730
731}
732
733//__________________________________________
5a516e0a 734void AliHFEmcQA::EndOfEventAna(Int_t kquark)
259c3296 735{
736 // end of event analysis
737
dbe3abbe 738 if (kquark != kCharm && kquark != kBeauty) {
259c3296 739 AliDebug(1, "This task is only for heavy quark QA, return\n");
740 return;
741 }
dbe3abbe 742 Int_t iq = kquark - kCharm;
259c3296 743
744
745 // # of heavy quark per event
746 AliDebug(1,Form("Number of heavy quark in this event = %d \n",fIsHeavy[iq]));
dbe3abbe 747 fHistComm[iq][0].fNq->Fill(fIsHeavy[iq]);
259c3296 748
749 Int_t motherID[fgkMaxGener];
750 Int_t motherType[fgkMaxGener];
751 Int_t motherLabel[fgkMaxGener];
752 Int_t ancestorPdg[fgkMaxGener];
753 Int_t ancestorLabel[fgkMaxGener];
754
755 for (Int_t i = 0; i < fgkMaxGener; i++){ // initialization
756 motherID[i] = 0;
757 motherType[i] = 0;
758 motherLabel[i] = 0;
759 ancestorPdg[i] = 0;
760 ancestorLabel[i] = 0;
761 }
762
3a72645a 763
259c3296 764 // check history of found heavy quarks
765 for (Int_t i = 0; i < fIsHeavy[iq]; i++){
766
3a72645a 767 if(!fHeavyQuark[i]) return;
768
259c3296 769 ancestorLabel[0] = i;
770 ancestorPdg[0] = fHeavyQuark[i]->GetPdgCode();
771 ancestorLabel[1] = fHeavyQuark[i]->GetFirstMother();
772
773 AliDebug(1,Form("pdg code= %d\n",ancestorPdg[0]));
774 AliDebug(1,Form("ancestor label= %d\n",ancestorLabel[1]));
775
776 Int_t ig = 1;
777 while (ancestorLabel[ig] != -1){
778 // in case there is mother, get mother's pdg code and grandmother's label
779 IdentifyMother(ancestorLabel[ig], ancestorPdg[ig], ancestorLabel[ig+1]);
780 // if mother is still heavy, find again mother's ancestor
781 if (ancestorPdg[ig-1] == ancestorPdg[ig]) {
782 ig++;
783 continue; // if it is from same heavy
784 }
785 // if the heavy's mother is not heavy, check the mother's label to know if it comes from inital or final parton shower
786 if (IsFromInitialShower(ancestorLabel[ig],motherID[i],motherType[i],motherLabel[i])) break;
787 if (IsFromFinalParton(ancestorLabel[ig],motherID[i],motherType[i],motherLabel[i])) break;
788 // if it is not the above case, something is strange
dbe3abbe 789 ReportStrangeness(motherID[i],motherType[i],motherLabel[i]);
259c3296 790 break;
791 }
792 if (ancestorLabel[ig] == -1){ // from hard scattering
793 HardScattering(kquark, motherID[i],motherType[i], motherLabel[i]);
794 }
795
796 } // end of found heavy quark loop
797
798
799 // check process type
800 Int_t processID = 0;
801 for (Int_t i = 0; i < fIsHeavy[iq]; i++){
802 AliDebug(1,Form("Mother ID= %d type= %d label= %d\n",motherID[i],motherType[i],motherLabel[i]));
803 }
804
805
806 Int_t nheavypair = Int_t(fIsHeavy[iq]/2.);
807 for (Int_t ipair = 0; ipair < nheavypair; ipair++){
808
809 Int_t id1 = ipair*2;
810 Int_t id2 = ipair*2 + 1;
811
812 if (motherType[id1] == 2 && motherType[id2] == 2){
dbe3abbe 813 if (motherLabel[id1] == motherLabel[id2]) processID = kGluonSplitting; // gluon spliting
259c3296 814 else processID = -9;
815 }
816 else if (motherType[id1] == -1 && motherType[id2] == -1) {
817 if (motherLabel[id1] == -1 && motherLabel[id2] == -1) {
dbe3abbe 818 if (motherID[id1] == fgkGluon) processID = kPairCreationFromg; // gluon fusion
819 else processID = kPairCreationFromq; // q-qbar pair creation
259c3296 820 }
821 else processID = -8;
822 }
823 else if (motherType[id1] == -1 || motherType[id2] == -1) {
824 if ((motherLabel[id1] == -1 || motherLabel[id2] == -1) && (motherLabel[id1]*motherLabel[id2] == -2 || motherLabel[id1]*motherLabel[id2] == -3)) {
dbe3abbe 825 if(motherID[id1]*motherID[id2] == kquark*fgkGluon) processID = kFlavourExitation; // flavour exitation
826 else processID = kLightQuarkShower;
259c3296 827 }
828 else processID = -7;
829 }
830 else if (motherType[id1] == -2 || motherType[id2] == -2) {
dbe3abbe 831 if (motherLabel[id1] == motherLabel[id2]) processID = kInitialPartonShower; // initial parton shower
259c3296 832 else processID = -6;
833
834 }
835 else processID = -5;
836
837 if (nheavypair >1) AliDebug(1,Form("Multi pair found : process ID = %d\n",processID));
dbe3abbe 838 else fHistComm[iq][0].fProcessID->Fill(processID);
259c3296 839 AliDebug(1,Form("Process ID = %d\n",processID));
840 } // end of # heavy quark pair loop
841
842}
843
844//__________________________________________
5a516e0a 845void AliHFEmcQA::GetHadronKine(TParticle* mcpart, Int_t kquark)
dbe3abbe 846{
847 // decay electron kinematics
848
849 if (kquark != kCharm && kquark != kBeauty) {
850 AliDebug(1, "This task is only for heavy quark QA, return\n");
851 return;
852 }
853 Int_t iq = kquark - kCharm;
854
faee3b18 855 if(!mcpart){
856 AliDebug(1, "no mc particle, return\n");
857 return;
858 }
dbe3abbe 859
860 Int_t iLabel = mcpart->GetFirstMother();
861 if (iLabel<0){
862 AliDebug(1, "Stack label is negative, return\n");
863 return;
864 }
865
11ff28c5 866 if(!fIsPbPb&&!fIsppMultiBin) fCentrality=0;
867
dbe3abbe 868 TParticle *partCopy = mcpart;
869 Int_t pdgcode = mcpart->GetPdgCode();
870 Int_t pdgcodeCopy = pdgcode;
871
faee3b18 872 AliMCParticle *mctrack = NULL;
873
dbe3abbe 874 // if the mother is charmed hadron
75d81601 875 Bool_t isDirectCharm = kFALSE;
b8bec44f 876 if ( int(TMath::Abs(pdgcode)/100.) == kCharm || int(TMath::Abs(pdgcode)/1000.) == kCharm ) {
dbe3abbe 877
878 // iterate until you find B hadron as a mother or become top ancester
879 for (Int_t i=1; i<fgkMaxIter; i++){
880
881 Int_t jLabel = mcpart->GetFirstMother();
882 if (jLabel == -1){
75d81601 883 isDirectCharm = kTRUE;
dbe3abbe 884 break; // if there is no ancester
885 }
886 if (jLabel < 0){ // safety protection
887 AliDebug(1, "Stack label is negative, return\n");
888 return;
889 }
890 // if there is an ancester
faee3b18 891 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return;
892 TParticle* mother = mctrack->Particle();
dbe3abbe 893 Int_t motherPDG = mother->GetPdgCode();
894
895 for (Int_t j=0; j<fNparents; j++){
b8bec44f 896 if (TMath::Abs(motherPDG)==fParentSelect[1][j]) return; // return if this hadron is originated from b
dbe3abbe 897 }
898
899 mcpart = mother;
900 } // end of iteration
901 } // end of if
75d81601 902 if((isDirectCharm == kTRUE && kquark == kCharm) || kquark == kBeauty) {
dbe3abbe 903 for (Int_t i=0; i<fNparents; i++){
b8bec44f 904 if (TMath::Abs(pdgcodeCopy)==fParentSelect[iq][i]){
dbe3abbe 905
906 // fill hadron kinematics
11ff28c5 907 fHist[iq][kHadron][0][fCentrality].fPdgCode->Fill(pdgcodeCopy);
908 fHist[iq][kHadron][0][fCentrality].fPt->Fill(partCopy->Pt());
909 fHist[iq][kHadron][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partCopy));
910 fHist[iq][kHadron][0][fCentrality].fEta->Fill(partCopy->Eta());
dbe3abbe 911
ccc37cdc 912 if(iq==0) {
913 fhD[i]->Fill(partCopy->Pt(),AliHFEtools::GetRapidity(partCopy));
ccc37cdc 914 }
dbe3abbe 915 }
916 }
ccc37cdc 917 // I also want to store D* info to compare with D* measurement
b8bec44f 918 if (TMath::Abs(pdgcodeCopy)==413 && iq==0) { //D*+
ccc37cdc 919 fhD[7]->Fill(partCopy->Pt(),AliHFEtools::GetRapidity(partCopy));
ccc37cdc 920 }
b8bec44f 921 if (TMath::Abs(pdgcodeCopy)==423 && iq==0) { //D*0
ccc37cdc 922 fhD[8]->Fill(partCopy->Pt(),AliHFEtools::GetRapidity(partCopy));
ccc37cdc 923 }
dbe3abbe 924 } // end of if
925}
926
927//__________________________________________
5a516e0a 928void AliHFEmcQA::GetDecayedKine(TParticle* mcpart, Int_t kquark, Int_t kdecayed)
259c3296 929{
930 // decay electron kinematics
931
faee3b18 932 if (!(kquark == kCharm || kquark == kBeauty || kquark == kOthers)){
259c3296 933 AliDebug(1, "This task is only for heavy quark QA, return\n");
934 return;
935 }
dbe3abbe 936 Int_t iq = kquark - kCharm;
50685501 937 Bool_t isFinalOpenCharm = kFALSE;
259c3296 938
faee3b18 939 if(!mcpart){
940 AliDebug(1, "no mcparticle, return\n");
941 return;
259c3296 942 }
943
11ff28c5 944 if(!fIsPbPb&&!fIsppMultiBin) fCentrality=0;
945
a8ef1999 946 Double_t eabsEta = TMath::Abs(mcpart->Eta());
947 Double_t eabsY = TMath::Abs(AliHFEtools::GetRapidity(mcpart));
948
faee3b18 949 if(kquark==kOthers){
950 Int_t esource = -1;
b8bec44f 951 if ( TMath::Abs(mcpart->GetPdgCode()) != kdecayed ) esource = kMisID-4;
faee3b18 952 else esource =GetSource(mcpart)-4; // return for the cases kGamma=4, kPi0=5, kElse=6
953 if(esource==0|| esource==1 || esource==2 || esource==3){
11ff28c5 954 fHist[iq][esource][0][fCentrality].fPt->Fill(mcpart->Pt());
955 fHist[iq][esource][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
956 fHist[iq][esource][0][fCentrality].fEta->Fill(mcpart->Eta());
a8ef1999 957 if(eabsEta<0.9){
11ff28c5 958 fHist[iq][esource][1][fCentrality].fPt->Fill(mcpart->Pt());
959 fHist[iq][esource][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
960 fHist[iq][esource][1][fCentrality].fEta->Fill(mcpart->Eta());
a8ef1999 961 }
962 if(eabsY<0.5){
11ff28c5 963 fHist[iq][esource][2][fCentrality].fPt->Fill(mcpart->Pt());
964 fHist[iq][esource][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
965 fHist[iq][esource][2][fCentrality].fEta->Fill(mcpart->Eta());
a8ef1999 966 }
faee3b18 967 return;
968 }
969 else {
970 AliDebug(1, "e source is out of defined ranges, return\n");
971 return;
972 }
973 }
259c3296 974
b8bec44f 975 if ( TMath::Abs(mcpart->GetPdgCode()) != kdecayed ) return;
259c3296 976
977 Int_t iLabel = mcpart->GetFirstMother();
978 if (iLabel<0){
979 AliDebug(1, "Stack label is negative, return\n");
980 return;
981 }
982
faee3b18 983 AliMCParticle *mctrack = NULL;
984 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return;
985 TParticle *partMother = mctrack->Particle();
dbe3abbe 986 TParticle *partMotherCopy = partMother;
259c3296 987 Int_t maPdgcode = partMother->GetPdgCode();
dbe3abbe 988 Int_t maPdgcodeCopy = maPdgcode;
259c3296 989
9bcfd1ab 990 // get mc primary vertex
faee3b18 991 /*
9bcfd1ab 992 TArrayF mcPrimVtx;
993 if(fMCHeader) fMCHeader->PrimaryVertex(mcPrimVtx);
994
259c3296 995 // get electron production vertex
996 TLorentzVector ePoint;
997 mcpart->ProductionVertex(ePoint);
998
9bcfd1ab 999 // calculated production vertex to primary vertex (in xy plane)
1000 Float_t decayLxy = TMath::Sqrt((mcPrimVtx[0]-ePoint[0])*(mcPrimVtx[0]-ePoint[0])+(mcPrimVtx[1]-ePoint[1])*(mcPrimVtx[1]-ePoint[1]));
faee3b18 1001 */
1002 Float_t decayLxy = 0;
9bcfd1ab 1003
259c3296 1004 // if the mother is charmed hadron
dbe3abbe 1005 Bool_t isMotherDirectCharm = kFALSE;
b8bec44f 1006 if ( int(TMath::Abs(maPdgcode)/100.) == kCharm || int(TMath::Abs(maPdgcode)/1000.) == kCharm ) {
259c3296 1007
0792aa82 1008 for (Int_t i=0; i<fNparents; i++){
b8bec44f 1009 if (TMath::Abs(maPdgcode)==fParentSelect[0][i]){
50685501 1010 isFinalOpenCharm = kTRUE;
0792aa82 1011 }
1012 }
50685501 1013 if (!isFinalOpenCharm) return ;
0792aa82 1014
259c3296 1015 // iterate until you find B hadron as a mother or become top ancester
1016 for (Int_t i=1; i<fgkMaxIter; i++){
1017
1018 Int_t jLabel = partMother->GetFirstMother();
1019 if (jLabel == -1){
dbe3abbe 1020 isMotherDirectCharm = kTRUE;
259c3296 1021 break; // if there is no ancester
1022 }
1023 if (jLabel < 0){ // safety protection
1024 AliDebug(1, "Stack label is negative, return\n");
1025 return;
1026 }
1027
1028 // if there is an ancester
faee3b18 1029 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return;
1030 TParticle* grandMa = mctrack->Particle();
259c3296 1031 Int_t grandMaPDG = grandMa->GetPdgCode();
1032
fc0de2a0 1033 for (Int_t j=0; j<fNparents; j++){
b8bec44f 1034 if (TMath::Abs(grandMaPDG)==fParentSelect[1][j]){
259c3296 1035
dbe3abbe 1036 if (kquark == kCharm) return;
259c3296 1037 // fill electron kinematics
11ff28c5 1038 fHist[iq][kElectron2nd][0][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1039 fHist[iq][kElectron2nd][0][fCentrality].fPt->Fill(mcpart->Pt());
1040 fHist[iq][kElectron2nd][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1041 fHist[iq][kElectron2nd][0][fCentrality].fEta->Fill(mcpart->Eta());
259c3296 1042
1043 // fill mother hadron kinematics
11ff28c5 1044 fHist[iq][kDeHadron][0][fCentrality].fPdgCode->Fill(grandMaPDG);
1045 fHist[iq][kDeHadron][0][fCentrality].fPt->Fill(grandMa->Pt());
1046 fHist[iq][kDeHadron][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(grandMa));
1047 fHist[iq][kDeHadron][0][fCentrality].fEta->Fill(grandMa->Eta());
a8ef1999 1048
1049 if(eabsEta<0.9){
11ff28c5 1050 fHist[iq][kElectron2nd][1][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1051 fHist[iq][kElectron2nd][1][fCentrality].fPt->Fill(mcpart->Pt());
1052 fHist[iq][kElectron2nd][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1053 fHist[iq][kElectron2nd][1][fCentrality].fEta->Fill(mcpart->Eta());
a8ef1999 1054
1055 // fill mother hadron kinematics
11ff28c5 1056 fHist[iq][kDeHadron][1][fCentrality].fPdgCode->Fill(grandMaPDG);
1057 fHist[iq][kDeHadron][1][fCentrality].fPt->Fill(grandMa->Pt());
1058 fHist[iq][kDeHadron][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(grandMa));
1059 fHist[iq][kDeHadron][1][fCentrality].fEta->Fill(grandMa->Eta());
a8ef1999 1060 }
1061
1062 if(eabsY<0.5){
11ff28c5 1063 fHist[iq][kElectron2nd][2][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1064 fHist[iq][kElectron2nd][2][fCentrality].fPt->Fill(mcpart->Pt());
1065 fHist[iq][kElectron2nd][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1066 fHist[iq][kElectron2nd][2][fCentrality].fEta->Fill(mcpart->Eta());
a8ef1999 1067
1068 // fill mother hadron kinematics
11ff28c5 1069 fHist[iq][kDeHadron][2][fCentrality].fPdgCode->Fill(grandMaPDG);
1070 fHist[iq][kDeHadron][2][fCentrality].fPt->Fill(grandMa->Pt());
1071 fHist[iq][kDeHadron][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(grandMa));
1072 fHist[iq][kDeHadron][2][fCentrality].fEta->Fill(grandMa->Eta());
a8ef1999 1073 }
1074
cedf0381 1075 //mj: to calculate B to e eta correlation to calculate total heavy quark cross section
1076 Int_t kLabel0 = grandMa->GetFirstMother();
1077 Bool_t isGGrandmaYes = kFALSE;
1078 Double_t ggmrapidwstmp=0;
1079 if (!(kLabel0 < 0)){ // safety protection
1080 if((mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(kLabel0))))){
1081 TParticle* ggrandMatmp = mctrack->Particle();
1082 Int_t ggrandMaPDGtmp = ggrandMatmp->GetPdgCode();
b8bec44f 1083 if ( int(TMath::Abs(ggrandMaPDGtmp)/100.) == kBeauty || int(TMath::Abs(ggrandMaPDGtmp)/1000.) == kBeauty) isGGrandmaYes = kTRUE;
cedf0381 1084 ggmrapidwstmp = AliHFEtools::GetRapidity(ggrandMatmp);
1085 }
1086 }
1087
1088 Double_t gmrapidwstmp0 = AliHFEtools::GetRapidity(grandMa);
1089 Double_t eetawstmp0 = mcpart->Eta();
1090
1091 Double_t gmrapidtmp0 = TMath::Abs(gmrapidwstmp0);
1092 Double_t eetatmp0 = TMath::Abs(eetawstmp0);
1093
1094 fHistComm[iq][0].fEtaCorrB->Fill(gmrapidwstmp0,eetawstmp0);
1095 if(isGGrandmaYes) fHistComm[iq][0].fEtaCorrGB->Fill(ggmrapidwstmp,eetawstmp0);
1096 else fHistComm[iq][0].fEtaCorrGB->Fill(gmrapidwstmp0,eetawstmp0);
1097
1098 if(gmrapidtmp0<0.5 && eetatmp0<0.5 ) fHistComm[iq][0].fPtCorrBinein->Fill(grandMa->Pt(),mcpart->Pt());
1099 else if(gmrapidtmp0<0.5 && eetatmp0>0.5 ) fHistComm[iq][0].fPtCorrBineout->Fill(grandMa->Pt(),mcpart->Pt());
1100 else if(gmrapidtmp0>0.5 && eetatmp0<0.5 ) fHistComm[iq][0].fPtCorrBoutein->Fill(grandMa->Pt(),mcpart->Pt());
1101 else if(gmrapidtmp0>0.5 && eetatmp0>0.5 ) fHistComm[iq][0].fPtCorrBouteout->Fill(grandMa->Pt(),mcpart->Pt());
1102 //======================================================================================
1103
259c3296 1104 // ratio between pT of electron and pT of mother B hadron
a8ef1999 1105 if(grandMa->Pt()) {
1106 fHistComm[iq][0].fDePtRatio->Fill(grandMa->Pt(),mcpart->Pt()/grandMa->Pt());
1107 if(eabsEta<0.9){
1108 fHistComm[iq][1].fDePtRatio->Fill(grandMa->Pt(),mcpart->Pt()/grandMa->Pt());
1109 }
1110 if(eabsY<0.5){
1111 fHistComm[iq][2].fDePtRatio->Fill(grandMa->Pt(),mcpart->Pt()/grandMa->Pt());
1112 }
1113 }
259c3296 1114
9bcfd1ab 1115 // distance between electron production point and primary vertex
a8ef1999 1116 fHistComm[iq][0].fDeDistance->Fill(grandMa->Pt(),decayLxy);
1117 if(eabsEta<0.9){
1118 fHistComm[iq][1].fDeDistance->Fill(grandMa->Pt(),decayLxy);
1119 }
1120 if(eabsY<0.5){
1121 fHistComm[iq][2].fDeDistance->Fill(grandMa->Pt(),decayLxy);
1122 }
259c3296 1123 return;
1124 }
dc306130 1125 }
259c3296 1126
1127 partMother = grandMa;
1128 } // end of iteration
1129 } // end of if
dbe3abbe 1130 if((isMotherDirectCharm == kTRUE && kquark == kCharm) || kquark == kBeauty) {
259c3296 1131 for (Int_t i=0; i<fNparents; i++){
b8bec44f 1132 if (TMath::Abs(maPdgcodeCopy)==fParentSelect[iq][i]){
259c3296 1133
11ff28c5 1134 fHist[iq][kElectron][0][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1135 fHist[iq][kElectron][0][fCentrality].fPt->Fill(mcpart->Pt());
1136 fHist[iq][kElectron][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1137 fHist[iq][kElectron][0][fCentrality].fEta->Fill(mcpart->Eta());
a8ef1999 1138
259c3296 1139 // fill mother hadron kinematics
11ff28c5 1140 fHist[iq][keHadron][0][fCentrality].fPdgCode->Fill(maPdgcodeCopy);
1141 fHist[iq][keHadron][0][fCentrality].fPt->Fill(partMotherCopy->Pt());
1142 fHist[iq][keHadron][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
1143 fHist[iq][keHadron][0][fCentrality].fEta->Fill(partMotherCopy->Eta());
a8ef1999 1144
1145 if(eabsEta<0.9){
11ff28c5 1146 fHist[iq][kElectron][1][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1147 fHist[iq][kElectron][1][fCentrality].fPt->Fill(mcpart->Pt());
1148 fHist[iq][kElectron][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1149 fHist[iq][kElectron][1][fCentrality].fEta->Fill(mcpart->Eta());
a8ef1999 1150
1151 // fill mother hadron kinematics
11ff28c5 1152 fHist[iq][keHadron][1][fCentrality].fPdgCode->Fill(maPdgcodeCopy);
1153 fHist[iq][keHadron][1][fCentrality].fPt->Fill(partMotherCopy->Pt());
1154 fHist[iq][keHadron][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
1155 fHist[iq][keHadron][1][fCentrality].fEta->Fill(partMotherCopy->Eta());
a8ef1999 1156 }
1157
1158 if(eabsY<0.5){
11ff28c5 1159 fHist[iq][kElectron][2][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1160 fHist[iq][kElectron][2][fCentrality].fPt->Fill(mcpart->Pt());
1161 fHist[iq][kElectron][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1162 fHist[iq][kElectron][2][fCentrality].fEta->Fill(mcpart->Eta());
a8ef1999 1163
1164 // fill mother hadron kinematics
11ff28c5 1165 fHist[iq][keHadron][2][fCentrality].fPdgCode->Fill(maPdgcodeCopy);
1166 fHist[iq][keHadron][2][fCentrality].fPt->Fill(partMotherCopy->Pt());
1167 fHist[iq][keHadron][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
1168 fHist[iq][keHadron][2][fCentrality].fEta->Fill(partMotherCopy->Eta());
a8ef1999 1169 }
259c3296 1170
ccc37cdc 1171 // ratio between pT of electron and pT of mother B or direct D hadron
a8ef1999 1172 if(partMotherCopy->Pt()) {
1173 fHistComm[iq][0].fePtRatio->Fill(partMotherCopy->Pt(),mcpart->Pt()/partMotherCopy->Pt());
1174 if(eabsEta<0.9){
1175 fHistComm[iq][1].fePtRatio->Fill(partMotherCopy->Pt(),mcpart->Pt()/partMotherCopy->Pt());
1176 }
1177 if(eabsY<0.5){
1178 fHistComm[iq][2].fePtRatio->Fill(partMotherCopy->Pt(),mcpart->Pt()/partMotherCopy->Pt());
1179 }
1180 }
1181 fHistComm[iq][0].fPtCorr->Fill(partMotherCopy->Pt(),mcpart->Pt());
1182 if(eabsEta<0.9){
1183 fHistComm[iq][1].fPtCorr->Fill(partMotherCopy->Pt(),mcpart->Pt());
1184 }
1185 if(eabsY<0.5){
1186 fHistComm[iq][2].fPtCorr->Fill(partMotherCopy->Pt(),mcpart->Pt());
1187 }
1188 if(TMath::Abs(partMotherCopy->GetPdgCode())==411) {
1189 fHistComm[iq][0].fPtCorrDp->Fill(partMotherCopy->Pt(),mcpart->Pt());
1190 if(eabsEta<0.9){
1191 fHistComm[iq][1].fPtCorrDp->Fill(partMotherCopy->Pt(),mcpart->Pt());
1192 }
1193 if(eabsY<0.5){
1194 fHistComm[iq][2].fPtCorrDp->Fill(partMotherCopy->Pt(),mcpart->Pt());
1195 }
1196 }
1197 else if(TMath::Abs(partMotherCopy->GetPdgCode())==421) {
1198 fHistComm[iq][0].fPtCorrD0->Fill(partMotherCopy->Pt(),mcpart->Pt());
1199 if(eabsEta<0.9){
1200 fHistComm[iq][1].fPtCorrD0->Fill(partMotherCopy->Pt(),mcpart->Pt());
1201 }
1202 if(eabsY<0.5){
1203 fHistComm[iq][2].fPtCorrD0->Fill(partMotherCopy->Pt(),mcpart->Pt());
1204 }
1205 }
1206 else {
1207 fHistComm[iq][0].fPtCorrDrest->Fill(partMotherCopy->Pt(),mcpart->Pt());
1208 if(eabsEta<0.9){
1209 fHistComm[iq][1].fPtCorrDrest->Fill(partMotherCopy->Pt(),mcpart->Pt());
1210 }
1211 if(eabsY<0.5){
1212 fHistComm[iq][2].fPtCorrDrest->Fill(partMotherCopy->Pt(),mcpart->Pt());
1213 }
1214 }
259c3296 1215
cedf0381 1216 //mj: to calculate D to e eta correlation to calculate total heavy quark cross section
1217 Int_t kLabel = partMotherCopy->GetFirstMother();
1218 Bool_t isGrandmaYes = kFALSE;
1219 Double_t gmrapidwstmp =0;
1220 if (!(kLabel < 0)){ // safety protection
1221 if((mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(kLabel))))){
1222 TParticle* grandMatmp = mctrack->Particle();
1223 Int_t grandMaPDGtmp = grandMatmp->GetPdgCode();
b8bec44f 1224 if ( int(TMath::Abs(grandMaPDGtmp)/100.) == kCharm || int(TMath::Abs(grandMaPDGtmp)/1000.) == kCharm ) isGrandmaYes = kTRUE;
1225 if ( int(TMath::Abs(grandMaPDGtmp)/100.) == kBeauty || int(TMath::Abs(grandMaPDGtmp)/1000.) == kBeauty) isGrandmaYes = kTRUE;
cedf0381 1226 gmrapidwstmp = AliHFEtools::GetRapidity(grandMatmp);
1227 }
1228 }
1229
1230 Double_t mrapidwstmp = AliHFEtools::GetRapidity(partMotherCopy);
1231 Double_t eetawstmp = mcpart->Eta();
1232
1233 Double_t mrapidtmp = TMath::Abs(mrapidwstmp);
1234 Double_t eetatmp = TMath::Abs(eetawstmp);
1235
1236 fHistComm[iq][0].fEtaCorrD->Fill(mrapidwstmp,eetawstmp);
1237 if(isGrandmaYes) fHistComm[iq][0].fEtaCorrGD->Fill(gmrapidwstmp,eetawstmp);
1238 else fHistComm[iq][0].fEtaCorrGD->Fill(mrapidwstmp,eetawstmp);
1239
1240 if(mrapidtmp<0.5 && eetatmp<0.5 ) fHistComm[iq][0].fPtCorrDinein->Fill(partMotherCopy->Pt(),mcpart->Pt());
1241 else if(mrapidtmp<0.5 && eetatmp>0.5 ) fHistComm[iq][0].fPtCorrDineout->Fill(partMotherCopy->Pt(),mcpart->Pt());
1242 else if(mrapidtmp>0.5 && eetatmp<0.5 ) fHistComm[iq][0].fPtCorrDoutein->Fill(partMotherCopy->Pt(),mcpart->Pt());
1243 else if(mrapidtmp>0.5 && eetatmp>0.5 ) fHistComm[iq][0].fPtCorrDouteout->Fill(partMotherCopy->Pt(),mcpart->Pt());
1244 if(TMath::Abs(partMotherCopy->GetPdgCode())==411) {
1245 fHistComm[iq][0].fEtaCorrDp->Fill(mrapidwstmp,eetawstmp);
1246 if(isGrandmaYes) fHistComm[iq][0].fEtaCorrGDp->Fill(gmrapidwstmp,eetawstmp);
1247 else fHistComm[iq][0].fEtaCorrGDp->Fill(mrapidwstmp,eetawstmp);
1248 if(mrapidtmp<0.5 && eetatmp<0.5 ) fHistComm[iq][0].fPtCorrDpDinein->Fill(partMotherCopy->Pt(),mcpart->Pt());
1249 else if(mrapidtmp<0.5 && eetatmp>0.5 ) fHistComm[iq][0].fPtCorrDpDineout->Fill(partMotherCopy->Pt(),mcpart->Pt());
1250 else if(mrapidtmp>0.5 && eetatmp<0.5 ) fHistComm[iq][0].fPtCorrDpDoutein->Fill(partMotherCopy->Pt(),mcpart->Pt());
1251 else if(mrapidtmp>0.5 && eetatmp>0.5 ) fHistComm[iq][0].fPtCorrDpDouteout->Fill(partMotherCopy->Pt(),mcpart->Pt());
1252 }
1253 else if(TMath::Abs(partMotherCopy->GetPdgCode())==421) {
1254 fHistComm[iq][0].fEtaCorrD0->Fill(mrapidwstmp,eetawstmp);
1255 if(isGrandmaYes) fHistComm[iq][0].fEtaCorrGD0->Fill(gmrapidwstmp,eetawstmp);
1256 else fHistComm[iq][0].fEtaCorrGD0->Fill(mrapidwstmp,eetawstmp);
1257 if(mrapidtmp<0.5 && eetatmp<0.5 ) fHistComm[iq][0].fPtCorrD0Dinein->Fill(partMotherCopy->Pt(),mcpart->Pt());
1258 else if(mrapidtmp<0.5 && eetatmp>0.5 ) fHistComm[iq][0].fPtCorrD0Dineout->Fill(partMotherCopy->Pt(),mcpart->Pt());
1259 else if(mrapidtmp>0.5 && eetatmp<0.5 ) fHistComm[iq][0].fPtCorrD0Doutein->Fill(partMotherCopy->Pt(),mcpart->Pt());
1260 else if(mrapidtmp>0.5 && eetatmp>0.5 ) fHistComm[iq][0].fPtCorrD0Douteout->Fill(partMotherCopy->Pt(),mcpart->Pt());
1261 }
1262 else {
1263 fHistComm[iq][0].fEtaCorrDrest->Fill(mrapidwstmp,eetawstmp);
1264 if(isGrandmaYes) fHistComm[iq][0].fEtaCorrGDrest->Fill(gmrapidwstmp,eetawstmp);
1265 else fHistComm[iq][0].fEtaCorrGDrest->Fill(mrapidwstmp,eetawstmp);
1266 if(mrapidtmp<0.5 && eetatmp<0.5 ) fHistComm[iq][0].fPtCorrDrestDinein->Fill(partMotherCopy->Pt(),mcpart->Pt());
1267 else if(mrapidtmp<0.5 && eetatmp>0.5 ) fHistComm[iq][0].fPtCorrDrestDineout->Fill(partMotherCopy->Pt(),mcpart->Pt());
1268 else if(mrapidtmp>0.5 && eetatmp<0.5 ) fHistComm[iq][0].fPtCorrDrestDoutein->Fill(partMotherCopy->Pt(),mcpart->Pt());
1269 else if(mrapidtmp>0.5 && eetatmp>0.5 ) fHistComm[iq][0].fPtCorrDrestDouteout->Fill(partMotherCopy->Pt(),mcpart->Pt());
1270 }
1271
9bcfd1ab 1272 // distance between electron production point and primary vertex
a8ef1999 1273 fHistComm[iq][0].feDistance->Fill(partMotherCopy->Pt(),decayLxy);
1274 if(eabsEta<0.9){
1275 fHistComm[iq][1].feDistance->Fill(partMotherCopy->Pt(),decayLxy);
1276 }
1277 if(eabsY<0.5){
1278 fHistComm[iq][2].feDistance->Fill(partMotherCopy->Pt(),decayLxy);
1279 }
259c3296 1280 }
1281 }
1282 } // end of if
1283}
1284
0792aa82 1285//____________________________________________________________________
5a516e0a 1286void AliHFEmcQA::GetDecayedKine(AliAODMCParticle *mcpart, Int_t kquark, Int_t kdecayed)
0792aa82 1287{
1288 // decay electron kinematics
1289
1290 if (kquark != kCharm && kquark != kBeauty) {
1291 AliDebug(1, "This task is only for heavy quark QA, return\n");
1292 return;
1293 }
1294
1295 Int_t iq = kquark - kCharm;
50685501 1296 Bool_t isFinalOpenCharm = kFALSE;
0792aa82 1297
faee3b18 1298 if(!mcpart){
1299 AliDebug(1, "no mcparticle, return\n");
1300 return;
1301 }
1302
b8bec44f 1303 if ( TMath::Abs(mcpart->GetPdgCode()) != kdecayed ) return;
0792aa82 1304
a8ef1999 1305 Double_t eabsEta = TMath::Abs(mcpart->Eta());
1306 Double_t eabsY = TMath::Abs(AliHFEtools::GetRapidity(mcpart));
1307
0792aa82 1308 // mother
1309 Int_t iLabel = mcpart->GetMother();
1310 if (iLabel<0){
1311 AliDebug(1, "Stack label is negative, return\n");
1312 return;
1313 }
1314
11ff28c5 1315 if(!fIsPbPb&&!fIsppMultiBin) fCentrality=0;
1316
0792aa82 1317 AliAODMCParticle *partMother = (AliAODMCParticle*)fMCArray->At(iLabel);
1318 AliAODMCParticle *partMotherCopy = partMother;
1319 Int_t maPdgcode = partMother->GetPdgCode();
1320 Int_t maPdgcodeCopy = maPdgcode;
1321
1322 Bool_t isMotherDirectCharm = kFALSE;
b8bec44f 1323 if ( int(TMath::Abs(maPdgcode)/100.) == kCharm || int(TMath::Abs(maPdgcode)/1000.) == kCharm ) {
0792aa82 1324
1325 for (Int_t i=0; i<fNparents; i++){
b8bec44f 1326 if (TMath::Abs(maPdgcode)==fParentSelect[0][i]){
50685501 1327 isFinalOpenCharm = kTRUE;
0792aa82 1328 }
1329 }
50685501 1330 if (!isFinalOpenCharm) return;
0792aa82 1331
1332 for (Int_t i=1; i<fgkMaxIter; i++){
1333
1334 Int_t jLabel = partMother->GetMother();
1335 if (jLabel == -1){
1336 isMotherDirectCharm = kTRUE;
1337 break; // if there is no ancester
1338 }
1339 if (jLabel < 0){ // safety protection
1340 AliDebug(1, "Stack label is negative, return\n");
1341 return;
1342 }
1343
1344 // if there is an ancester
1345 AliAODMCParticle* grandMa = (AliAODMCParticle*)fMCArray->At(jLabel);
1346 Int_t grandMaPDG = grandMa->GetPdgCode();
1347
1348 for (Int_t j=0; j<fNparents; j++){
b8bec44f 1349 if (TMath::Abs(grandMaPDG)==fParentSelect[1][j]){
0792aa82 1350
1351 if (kquark == kCharm) return;
1352 // fill electron kinematics
11ff28c5 1353 fHist[iq][kElectron2nd][0][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1354 fHist[iq][kElectron2nd][0][fCentrality].fPt->Fill(mcpart->Pt());
1355 fHist[iq][kElectron2nd][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1356 fHist[iq][kElectron2nd][0][fCentrality].fEta->Fill(mcpart->Eta());
0792aa82 1357
1358 // fill mother hadron kinematics
11ff28c5 1359 fHist[iq][kDeHadron][0][fCentrality].fPdgCode->Fill(grandMaPDG);
1360 fHist[iq][kDeHadron][0][fCentrality].fPt->Fill(grandMa->Pt());
1361 fHist[iq][kDeHadron][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(grandMa));
1362 fHist[iq][kDeHadron][0][fCentrality].fEta->Fill(grandMa->Eta());
a8ef1999 1363
1364 if(eabsEta<0.9){
1365 // fill electron kinematics
11ff28c5 1366 fHist[iq][kElectron2nd][1][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1367 fHist[iq][kElectron2nd][1][fCentrality].fPt->Fill(mcpart->Pt());
1368 fHist[iq][kElectron2nd][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1369 fHist[iq][kElectron2nd][1][fCentrality].fEta->Fill(mcpart->Eta());
a8ef1999 1370
1371 // fill mother hadron kinematics
11ff28c5 1372 fHist[iq][kDeHadron][1][fCentrality].fPdgCode->Fill(grandMaPDG);
1373 fHist[iq][kDeHadron][1][fCentrality].fPt->Fill(grandMa->Pt());
1374 fHist[iq][kDeHadron][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(grandMa));
1375 fHist[iq][kDeHadron][1][fCentrality].fEta->Fill(grandMa->Eta());
a8ef1999 1376 }
1377 if(eabsY<0.5){
1378 // fill electron kinematics
11ff28c5 1379 fHist[iq][kElectron2nd][2][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1380 fHist[iq][kElectron2nd][2][fCentrality].fPt->Fill(mcpart->Pt());
1381 fHist[iq][kElectron2nd][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1382 fHist[iq][kElectron2nd][2][fCentrality].fEta->Fill(mcpart->Eta());
a8ef1999 1383
1384 // fill mother hadron kinematics
11ff28c5 1385 fHist[iq][kDeHadron][2][fCentrality].fPdgCode->Fill(grandMaPDG);
1386 fHist[iq][kDeHadron][2][fCentrality].fPt->Fill(grandMa->Pt());
1387 fHist[iq][kDeHadron][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(grandMa));
1388 fHist[iq][kDeHadron][2][fCentrality].fEta->Fill(grandMa->Eta());
a8ef1999 1389 }
0792aa82 1390
1391 return;
1392 }
1393 }
1394
1395 partMother = grandMa;
1396 } // end of iteration
1397 } // end of if
1398 if ((isMotherDirectCharm == kTRUE && kquark == kCharm) || kquark == kBeauty) {
1399 for (Int_t i=0; i<fNparents; i++){
b8bec44f 1400 if (TMath::Abs(maPdgcodeCopy)==fParentSelect[iq][i]){
0792aa82 1401
1402 // fill electron kinematics
11ff28c5 1403 fHist[iq][kElectron][0][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1404 fHist[iq][kElectron][0][fCentrality].fPt->Fill(mcpart->Pt());
1405 fHist[iq][kElectron][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1406 fHist[iq][kElectron][0][fCentrality].fEta->Fill(mcpart->Eta());
0792aa82 1407
1408 // fill mother hadron kinematics
11ff28c5 1409 fHist[iq][keHadron][0][fCentrality].fPdgCode->Fill(maPdgcodeCopy);
1410 fHist[iq][keHadron][0][fCentrality].fPt->Fill(partMotherCopy->Pt());
1411 fHist[iq][keHadron][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
1412 fHist[iq][keHadron][0][fCentrality].fEta->Fill(partMotherCopy->Eta());
a8ef1999 1413
1414 if(eabsEta<0.9){
1415 // fill electron kinematics
11ff28c5 1416 fHist[iq][kElectron][1][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1417 fHist[iq][kElectron][1][fCentrality].fPt->Fill(mcpart->Pt());
1418 fHist[iq][kElectron][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1419 fHist[iq][kElectron][1][fCentrality].fEta->Fill(mcpart->Eta());
a8ef1999 1420
1421 // fill mother hadron kinematics
11ff28c5 1422 fHist[iq][keHadron][1][fCentrality].fPdgCode->Fill(maPdgcodeCopy);
1423 fHist[iq][keHadron][1][fCentrality].fPt->Fill(partMotherCopy->Pt());
1424 fHist[iq][keHadron][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
1425 fHist[iq][keHadron][1][fCentrality].fEta->Fill(partMotherCopy->Eta());
a8ef1999 1426 }
1427 if(eabsY<0.5){
1428 // fill electron kinematics
11ff28c5 1429 fHist[iq][kElectron][2][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1430 fHist[iq][kElectron][2][fCentrality].fPt->Fill(mcpart->Pt());
1431 fHist[iq][kElectron][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1432 fHist[iq][kElectron][2][fCentrality].fEta->Fill(mcpart->Eta());
a8ef1999 1433
1434 // fill mother hadron kinematics
11ff28c5 1435 fHist[iq][keHadron][2][fCentrality].fPdgCode->Fill(maPdgcodeCopy);
1436 fHist[iq][keHadron][2][fCentrality].fPt->Fill(partMotherCopy->Pt());
1437 fHist[iq][keHadron][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
1438 fHist[iq][keHadron][2][fCentrality].fEta->Fill(partMotherCopy->Eta());
a8ef1999 1439 }
0792aa82 1440
1441 }
1442 }
1443 } // end of if
1444
1445}
259c3296 1446
1447//__________________________________________
75d81601 1448void AliHFEmcQA::IdentifyMother(Int_t motherlabel, Int_t &motherpdg, Int_t &grandmotherlabel)
259c3296 1449{
dbe3abbe 1450 // find mother pdg code and label
1451
75d81601 1452 if (motherlabel < 0) {
259c3296 1453 AliDebug(1, "Stack label is negative, return\n");
1454 return;
1455 }
faee3b18 1456 AliMCParticle *mctrack = NULL;
1457 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(motherlabel))))) return;
1458 TParticle *heavysMother = mctrack->Particle();
75d81601 1459 motherpdg = heavysMother->GetPdgCode();
1460 grandmotherlabel = heavysMother->GetFirstMother();
1461 AliDebug(1,Form("ancestor pdg code= %d\n",motherpdg));
259c3296 1462}
1463
1464//__________________________________________
5a516e0a 1465void AliHFEmcQA::HardScattering(Int_t kquark, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
259c3296 1466{
dbe3abbe 1467 // mothertype -1 means this heavy quark coming from hard vertex
1468
faee3b18 1469 AliMCParticle *mctrack1 = NULL;
1470 AliMCParticle *mctrack2 = NULL;
1471 if(!(mctrack1 = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(4))))) return;
1472 if(!(mctrack2 = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(5))))) return;
1473 TParticle *afterinitialrad1 = mctrack1->Particle();
1474 TParticle *afterinitialrad2 = mctrack2->Particle();
259c3296 1475
1476 motherlabel = -1;
1477
b8bec44f 1478 if (TMath::Abs(afterinitialrad1->GetPdgCode()) == fgkGluon && TMath::Abs(afterinitialrad2->GetPdgCode()) == fgkGluon){
259c3296 1479 AliDebug(1,"heavy from gluon gluon pair creation!\n");
1480 mothertype = -1;
1481 motherID = fgkGluon;
1482 }
b8bec44f 1483 else if (TMath::Abs(afterinitialrad1->GetPdgCode()) == kquark || TMath::Abs(afterinitialrad2->GetPdgCode()) == kquark){ // one from Q and the other from g
259c3296 1484 AliDebug(1,"heavy from flavor exitation!\n");
1485 mothertype = -1;
1486 motherID = kquark;
1487 }
b8bec44f 1488 else if (TMath::Abs(afterinitialrad1->GetPdgCode()) == TMath::Abs(afterinitialrad2->GetPdgCode())){
259c3296 1489 AliDebug(1,"heavy from q-qbar pair creation!\n");
1490 mothertype = -1;
1491 motherID = 1;
1492 }
1493 else {
1494 AliDebug(1,"something strange!\n");
1495 mothertype = -999;
1496 motherlabel = -999;
1497 motherID = -999;
1498 }
1499}
1500
1501//__________________________________________
259c3296 1502Bool_t AliHFEmcQA::IsFromInitialShower(Int_t inputmotherlabel, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
1503{
dbe3abbe 1504 // mothertype -2 means this heavy quark coming from initial state
1505
faee3b18 1506 AliMCParticle *mctrack = NULL;
259c3296 1507 if (inputmotherlabel==2 || inputmotherlabel==3){ // mother exist before initial state radiation
faee3b18 1508 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(inputmotherlabel))))) return kFALSE;
1509 TParticle *heavysMother = mctrack->Particle();
259c3296 1510 motherID = heavysMother->GetPdgCode();
1511 mothertype = -2; // there is mother before initial state radiation
1512 motherlabel = inputmotherlabel;
1513 AliDebug(1,"initial parton shower! \n");
1514
1515 return kTRUE;
1516 }
1517
1518 return kFALSE;
1519}
1520
1521//__________________________________________
259c3296 1522Bool_t AliHFEmcQA::IsFromFinalParton(Int_t inputmotherlabel, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
1523{
dbe3abbe 1524 // mothertype 2 means this heavy quark coming from final state
1525
faee3b18 1526 AliMCParticle *mctrack = NULL;
259c3296 1527 if (inputmotherlabel > 5){ // mother exist after hard scattering
faee3b18 1528 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(inputmotherlabel))))) return kFALSE;
1529 TParticle *heavysMother = mctrack->Particle();
259c3296 1530 motherID = heavysMother->GetPdgCode();
1531 mothertype = 2; //
1532 motherlabel = inputmotherlabel;
1533 AliDebug(1,Form("heavy quark from %d after hard scattering! \n",motherID));
1534
1535 return kTRUE;
1536 }
1537 return kFALSE;
1538}
1539
1540//__________________________________________
dbe3abbe 1541void AliHFEmcQA::ReportStrangeness(Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
259c3296 1542{
dbe3abbe 1543 // mark strange behavior
1544
259c3296 1545 mothertype = -888;
1546 motherlabel = -888;
1547 motherID = -888;
1548 AliDebug(1,"something strange!\n");
1549}
1550
0792aa82 1551//__________________________________________
afb48e1d 1552Int_t AliHFEmcQA::GetSource(const AliVParticle* const mcpart) const
0792aa82 1553{
9bcfd1ab 1554 // decay particle's origin
0792aa82 1555
b8bec44f 1556 //if ( TMath::Abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return -1;
0792aa82 1557
1558 Int_t origin = -1;
50685501 1559 Bool_t isFinalOpenCharm = kFALSE;
0792aa82 1560
faee3b18 1561 if(!mcpart){
1562 AliDebug(1, "Stack label is negative or no mcparticle, return\n");
1563 return -1;
1564 }
1565
0792aa82 1566 // mother
11ff28c5 1567 // Information not in the base class, cast necessary
1568 Int_t iLabel = GetMother(mcpart);
0792aa82 1569 if (iLabel<0){
1570 AliDebug(1, "Stack label is negative, return\n");
1571 return -1;
1572 }
1573
11ff28c5 1574 const AliVParticle *partMother = fMCEvent->GetTrack(iLabel);
1575 Int_t maPdgcode = partMother->PdgCode();
0792aa82 1576
1577 // if the mother is charmed hadron
b8bec44f 1578 if ( int(TMath::Abs(maPdgcode)/100.) == kCharm || int(TMath::Abs(maPdgcode)/1000.) == kCharm ) {
0792aa82 1579
1580 for (Int_t i=0; i<fNparents; i++){
b8bec44f 1581 if (TMath::Abs(maPdgcode)==fParentSelect[0][i]){
50685501 1582 isFinalOpenCharm = kTRUE;
0792aa82 1583 }
1584 }
50685501 1585 if (!isFinalOpenCharm) return -1;
0792aa82 1586
1587 // iterate until you find B hadron as a mother or become top ancester
1588 for (Int_t i=1; i<fgkMaxIter; i++){
11ff28c5 1589
1590 Int_t jLabel = GetMother(partMother);
0792aa82 1591 if (jLabel == -1){
1592 origin = kDirectCharm;
1593 return origin;
1594 }
1595 if (jLabel < 0){ // safety protection
1596 AliDebug(1, "Stack label is negative, return\n");
1597 return -1;
1598 }
1599
1600 // if there is an ancester
11ff28c5 1601 const AliVParticle* grandMa = fMCEvent->GetTrack(jLabel);
1602 Int_t grandMaPDG = grandMa->PdgCode();
0792aa82 1603
70b5ea26 1604 for (Int_t j=0; j<fNparents; j++){
b8bec44f 1605 if (TMath::Abs(grandMaPDG)==fParentSelect[1][j]){
0792aa82 1606 origin = kBeautyCharm;
1607 return origin;
1608 }
1609 }
1610
1611 partMother = grandMa;
1612 } // end of iteration
1613 } // end of if
b8bec44f 1614 else if ( int(TMath::Abs(maPdgcode)/100.) == kBeauty || int(TMath::Abs(maPdgcode)/1000.) == kBeauty ) {
0792aa82 1615 for (Int_t i=0; i<fNparents; i++){
b8bec44f 1616 if (TMath::Abs(maPdgcode)==fParentSelect[1][i]){
0792aa82 1617 origin = kDirectBeauty;
1618 return origin;
1619 }
1620 }
1621 } // end of if
b8bec44f 1622 else if ( TMath::Abs(maPdgcode) == 22 ) {
0792aa82 1623 origin = kGamma;
1624 return origin;
1625 } // end of if
b8bec44f 1626 else if ( TMath::Abs(maPdgcode) == 111 ) {
0792aa82 1627 origin = kPi0;
1628 return origin;
1629 } // end of if
1630
1631 return origin;
1632}
1633
1634//__________________________________________
afb48e1d 1635Int_t AliHFEmcQA::GetSource(const TParticle * const mcpart) const
0792aa82 1636{
9bcfd1ab 1637 // decay particle's origin
0792aa82 1638
b8bec44f 1639 //if ( TMath::Abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return -1;
0792aa82 1640
1641 Int_t origin = -1;
50685501 1642 Bool_t isFinalOpenCharm = kFALSE;
0792aa82 1643
faee3b18 1644 if(!mcpart){
1645 AliDebug(1, "no mcparticle, return\n");
1646 return -1;
1647 }
1648
0792aa82 1649 Int_t iLabel = mcpart->GetFirstMother();
1650 if (iLabel<0){
1651 AliDebug(1, "Stack label is negative, return\n");
1652 return -1;
1653 }
1654
faee3b18 1655 AliMCParticle *mctrack = NULL;
1656 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return -1;
1657 TParticle *partMother = mctrack->Particle();
0792aa82 1658 Int_t maPdgcode = partMother->GetPdgCode();
1659
1660 // if the mother is charmed hadron
b8bec44f 1661 if ( int(TMath::Abs(maPdgcode)/100.) == kCharm || int(TMath::Abs(maPdgcode)/1000.) == kCharm ) {
0792aa82 1662
1663 for (Int_t i=0; i<fNparents; i++){
b8bec44f 1664 if (TMath::Abs(maPdgcode)==fParentSelect[0][i]){
50685501 1665 isFinalOpenCharm = kTRUE;
0792aa82 1666 }
1667 }
50685501 1668 if (!isFinalOpenCharm) return -1;
0792aa82 1669
1670 // iterate until you find B hadron as a mother or become top ancester
1671 for (Int_t i=1; i<fgkMaxIter; i++){
1672
1673 Int_t jLabel = partMother->GetFirstMother();
1674 if (jLabel == -1){
1675 origin = kDirectCharm;
1676 return origin;
1677 }
1678 if (jLabel < 0){ // safety protection
1679 AliDebug(1, "Stack label is negative, return\n");
1680 return -1;
1681 }
1682
1683 // if there is an ancester
faee3b18 1684 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return -1;
1685 TParticle* grandMa = mctrack->Particle();
0792aa82 1686 Int_t grandMaPDG = grandMa->GetPdgCode();
1687
70b5ea26 1688 for (Int_t j=0; j<fNparents; j++){
b8bec44f 1689 if (TMath::Abs(grandMaPDG)==fParentSelect[1][j]){
0792aa82 1690 origin = kBeautyCharm;
1691 return origin;
1692 }
1693 }
1694
1695 partMother = grandMa;
1696 } // end of iteration
1697 } // end of if
b8bec44f 1698 else if ( int(TMath::Abs(maPdgcode)/100.) == kBeauty || int(TMath::Abs(maPdgcode)/1000.) == kBeauty ) {
0792aa82 1699 for (Int_t i=0; i<fNparents; i++){
b8bec44f 1700 if (TMath::Abs(maPdgcode)==fParentSelect[1][i]){
0792aa82 1701 origin = kDirectBeauty;
1702 return origin;
1703 }
1704 }
1705 } // end of if
b8bec44f 1706 else if ( TMath::Abs(maPdgcode) == 22 ) {
0792aa82 1707 origin = kGamma;
1708 return origin;
1709 } // end of if
b8bec44f 1710 else if ( TMath::Abs(maPdgcode) == 111 ) {
0792aa82 1711 origin = kPi0;
1712 return origin;
1713 } // end of if
1714 else origin = kElse;
1715
1716 return origin;
1717}
70da6c5a 1718
1719//__________________________________________
afb48e1d 1720Int_t AliHFEmcQA::GetElecSource(TParticle * const mcpart) const
e3fc062d 1721{
1722 // decay particle's origin
70da6c5a 1723
e3fc062d 1724 if(!mcpart){
1725 AliDebug(1, "no mcparticle, return\n");
1726 return -1;
1727 }
1728
b8bec44f 1729 if ( TMath::Abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return kMisID;
bf892a6a 1730
1731 Int_t origin = -1;
1732 Bool_t isFinalOpenCharm = kFALSE;
1733
e3fc062d 1734 Int_t iLabel = mcpart->GetFirstMother();
1735 if (iLabel<0){
1736 AliDebug(1, "Stack label is negative, return\n");
1737 return -1;
1738 }
1739
1740 AliMCParticle *mctrack = NULL;
c2690925 1741 Int_t tmpMomLabel=0;
e3fc062d 1742 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return -1;
1743 TParticle *partMother = mctrack->Particle();
c2690925 1744 TParticle *partMotherCopy = mctrack->Particle();
e3fc062d 1745 Int_t maPdgcode = partMother->GetPdgCode();
cedf0381 1746 Int_t grmaPdgcode;
1747 Int_t ggrmaPdgcode;
e3fc062d 1748
1749 // if the mother is charmed hadron
11ff28c5 1750
b8bec44f 1751 if(TMath::Abs(maPdgcode)==443){ // J/spi
11ff28c5 1752 Int_t jLabel = partMother->GetFirstMother();
cedf0381 1753 if((mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))){
1754 TParticle* grandMa = mctrack->Particle();
1755 Int_t grandMaPDG = grandMa->GetPdgCode();
b8bec44f 1756 if((int(TMath::Abs(grandMaPDG)/100.)%10) == kBeauty || (int(TMath::Abs(grandMaPDG)/1000.)%10) == kBeauty) return kB2Jpsi;
cedf0381 1757 }
1758 return kJpsi;
11ff28c5 1759 }
b8bec44f 1760 else if ( (int(TMath::Abs(maPdgcode)/100.)%10) == kCharm || (int(TMath::Abs(maPdgcode)/1000.)%10) == kCharm ) {
e3fc062d 1761
1762 for (Int_t i=0; i<fNparents; i++){
b8bec44f 1763 if (TMath::Abs(maPdgcode)==fParentSelect[0][i]){
e3fc062d 1764 isFinalOpenCharm = kTRUE;
1765 }
1766 }
1767 if (!isFinalOpenCharm) return -1;
1768
1769 // iterate until you find B hadron as a mother or become top ancester
1770 for (Int_t i=1; i<fgkMaxIter; i++){
1771
1772 Int_t jLabel = partMother->GetFirstMother();
1773 if (jLabel == -1){
11ff28c5 1774 return kDirectCharm;
e3fc062d 1775 }
1776 if (jLabel < 0){ // safety protection
1777 AliDebug(1, "Stack label is negative, return\n");
1778 return -1;
1779 }
1780
1781 // if there is an ancester
1782 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return -1;
1783 TParticle* grandMa = mctrack->Particle();
1784 Int_t grandMaPDG = grandMa->GetPdgCode();
1785
1786 for (Int_t j=0; j<fNparents; j++){
b8bec44f 1787 if (TMath::Abs(grandMaPDG)==fParentSelect[1][j]){
11ff28c5 1788 return kBeautyCharm;
e3fc062d 1789 }
1790 }
1791
1792 partMother = grandMa;
1793 } // end of iteration
1794 } // end of if
b8bec44f 1795 else if ( (int(TMath::Abs(maPdgcode)/100.)%10) == kBeauty || (int(TMath::Abs(maPdgcode)/1000.)%10) == kBeauty ) {
e3fc062d 1796 for (Int_t i=0; i<fNparents; i++){
b8bec44f 1797 if (TMath::Abs(maPdgcode)==fParentSelect[1][i]){
11ff28c5 1798 return kDirectBeauty;
e3fc062d 1799 }
1800 }
1801 } // end of if
b8bec44f 1802 else if ( TMath::Abs(maPdgcode) == 22 ) { //conversion
c2690925 1803
1804 tmpMomLabel = partMotherCopy->GetFirstMother();
1805 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tmpMomLabel))))) return -1;
1806 partMother = mctrack->Particle();
1807 maPdgcode = partMother->GetPdgCode();
cedf0381 1808
1809 // check if the ligth meson is the decay product of heavy mesons
1810 tmpMomLabel = partMother->GetFirstMother();
1811 if((mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tmpMomLabel))))) {
1812 partMother = mctrack->Particle();
1813 grmaPdgcode = partMother->GetPdgCode();
1814
b8bec44f 1815 if ( (int(TMath::Abs(grmaPdgcode)/100.)%10) == kBeauty || (int(TMath::Abs(grmaPdgcode)/1000.)%10) == kBeauty ) return kGammaB2M;
1816 if ( (int(TMath::Abs(grmaPdgcode)/100.)%10) == kCharm || (int(TMath::Abs(grmaPdgcode)/1000.)%10) == kCharm ) return kGammaD2M;
afb48e1d 1817 if ( TMath::Abs(grmaPdgcode) == 221 ) return kGammaEta2Pi0; //mj to remove secondary pi0 decay electrons from eta decays, mainly to eta signal enhance samples
1818 //if ( TMath::Abs(grmaPdgcode) == 221 ) return kElse; //mj to remove secondary pi0 decay electrons from eta decays, mainly to eta signal enhance samples
cedf0381 1819
1820 tmpMomLabel = partMother->GetFirstMother();
1821 if((mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tmpMomLabel))))) {
1822 partMother = mctrack->Particle();
1823 ggrmaPdgcode = partMother->GetPdgCode();
1824
b8bec44f 1825 if ( (int(TMath::Abs(ggrmaPdgcode)/100.)%10) == kBeauty || (int(TMath::Abs(ggrmaPdgcode)/1000.)%10) == kBeauty ) return kGammaB2M;
1826 if ( (int(TMath::Abs(ggrmaPdgcode)/100.)%10) == kCharm || (int(TMath::Abs(ggrmaPdgcode)/1000.)%10) == kCharm ) return kGammaD2M;
cedf0381 1827 }
1828 }
1829
b8bec44f 1830 if ( TMath::Abs(maPdgcode) == 111 ) {
11ff28c5 1831 return kGammaPi0;
c2690925 1832 }
b8bec44f 1833 else if ( TMath::Abs(maPdgcode) == 221 ) {
11ff28c5 1834 return kGammaEta;
c2690925 1835 }
b8bec44f 1836 else if ( TMath::Abs(maPdgcode) == 223 ) {
11ff28c5 1837 return kGammaOmega;
c2690925 1838 }
b8bec44f 1839 else if ( TMath::Abs(maPdgcode) == 333 ) {
11ff28c5 1840 return kGammaPhi;
c2690925 1841 }
b8bec44f 1842 else if ( TMath::Abs(maPdgcode) == 331 ) {
11ff28c5 1843 return kGammaEtaPrime;
c2690925 1844 }
b8bec44f 1845 else if ( TMath::Abs(maPdgcode) == 113 ) {
11ff28c5 1846 return kGammaRho0;
c2690925 1847 }
1848 else origin = kElse;
e3fc062d 1849 return origin;
c2690925 1850
11ff28c5 1851 }
cedf0381 1852 else {
1853
1854 // check if the ligth meson is the decay product of heavy mesons
1855 tmpMomLabel = partMotherCopy->GetFirstMother();
1856 if((mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tmpMomLabel))))) {
1857 partMother = mctrack->Particle();
1858 grmaPdgcode = partMother->GetPdgCode();
1859
afb48e1d 1860 if ( (int(TMath::Abs(grmaPdgcode)/100.)%10) == kBeauty || (int(TMath::Abs(grmaPdgcode)/1000.)%10) == kBeauty ) return kB2M;
1861 if ( (int(TMath::Abs(grmaPdgcode)/100.)%10) == kCharm || (int(TMath::Abs(grmaPdgcode)/1000.)%10) == kCharm ) return kD2M;
1862 if ( TMath::Abs(grmaPdgcode) == 221 ) return kEta2Pi0; //mj to remove secondary pi0 decay electrons from eta decays, mainly to eta signal enhance samples
cedf0381 1863
1864 tmpMomLabel = partMother->GetFirstMother();
1865 if((mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tmpMomLabel))))) {
1866 partMother = mctrack->Particle();
1867 ggrmaPdgcode = partMother->GetPdgCode();
1868
afb48e1d 1869 if ( (int(TMath::Abs(ggrmaPdgcode)/100.)%10) == kBeauty || (int(TMath::Abs(ggrmaPdgcode)/1000.)%10) == kBeauty ) return kB2M;
1870 if ( (int(TMath::Abs(ggrmaPdgcode)/100.)%10) == kCharm || (int(TMath::Abs(ggrmaPdgcode)/1000.)%10) == kCharm ) return kD2M;
cedf0381 1871 }
1872 }
1873
b8bec44f 1874 if ( TMath::Abs(maPdgcode) == 111 ) {
cedf0381 1875 return kPi0;
1876 }
b8bec44f 1877 else if ( TMath::Abs(maPdgcode) == 221 ) {
cedf0381 1878 return kEta;
1879 }
b8bec44f 1880 else if ( TMath::Abs(maPdgcode) == 223 ) {
cedf0381 1881 return kOmega;
1882 }
b8bec44f 1883 else if ( TMath::Abs(maPdgcode) == 333 ) {
cedf0381 1884 return kPhi;
1885 }
b8bec44f 1886 else if ( TMath::Abs(maPdgcode) == 331 ) {
cedf0381 1887 return kEtaPrime;
1888 }
b8bec44f 1889 else if ( TMath::Abs(maPdgcode) == 113 ) {
cedf0381 1890 return kRho0;
1891 }
b8bec44f 1892 else if ( TMath::Abs(maPdgcode) == 321 || TMath::Abs(maPdgcode) == 130 ) {
cedf0381 1893 return kKe3;
1894 }
1895 else{
1896 origin = kElse;
1897 }
1898
e17c1f86 1899 }
e3fc062d 1900 return origin;
70da6c5a 1901}
76d0b522 1902
1903//__________________________________________
afb48e1d 1904Int_t AliHFEmcQA::GetElecSource(const AliVParticle * const mctrack) const
76d0b522 1905{
1906 //
1907 // decay particle's origin
1908 //
1909
1910 if(!mctrack){
1911 AliDebug(1, "no mcparticle, return\n");
1912 return -1;
1913 }
1914
1915 TClass *type = mctrack->IsA();
1916
1917 if(type == AliMCParticle::Class()) {
1918 const AliMCParticle *esdmc = dynamic_cast<const AliMCParticle *>(mctrack);
1919 if(esdmc){
1920 TParticle *mcpart = esdmc->Particle();
1921 return GetElecSource(mcpart);
1922 }
1923 else return -1;
1924 }
1925 if(type == AliAODMCParticle::Class()) {
1926 const AliAODMCParticle *aodmc = dynamic_cast<const AliAODMCParticle *>(mctrack);
1927 if(aodmc){
1928 return GetElecSource(aodmc);
1929 }
1930 else return -1;
1931 }
1932 return -1;
1933}
1934//__________________________________________
afb48e1d 1935Int_t AliHFEmcQA::GetElecSource(const AliAODMCParticle * const mcpart) const
76d0b522 1936{
1937 //
1938 // Function for AliAODMCParticle
1939 //
1940
1941 if (!mcpart) return -1;
1942 if (!fMCArray) return -1;
1943
b8bec44f 1944 if ( TMath::Abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return kMisID;
76d0b522 1945
1946 Int_t origin = -1;
1947 Bool_t isFinalOpenCharm = kFALSE;
1948
1949 Int_t iLabel = mcpart->GetMother();
1950 if ((iLabel<0) || (iLabel>=fMCArray->GetEntriesFast())){
1951 AliDebug(1, "label is out of range, return\n");
1952 return -1;
1953 }
1954
1955 AliAODMCParticle *mctrack = NULL; // will change all the time
1956 Int_t tmpMomLabel=0;
1957 if(!(mctrack = dynamic_cast<AliAODMCParticle *>(fMCArray->At(TMath::Abs(iLabel))))) return -1;
1958 AliAODMCParticle *partMother = mctrack;
1959 AliAODMCParticle *partMotherCopy = mctrack;
1960 Int_t maPdgcode = mctrack->GetPdgCode();
1961 Int_t grmaPdgcode;
1962 Int_t ggrmaPdgcode;
1963
1964 // if the mother is charmed hadron
1965
b8bec44f 1966 if(TMath::Abs(maPdgcode)==443){
76d0b522 1967 //
1968 // J/spi
1969 //
1970 Int_t jLabel = partMother->GetMother();
959ea9d8 1971 if ((jLabel>=0) && (jLabel<fMCArray->GetEntriesFast())){
76d0b522 1972 if((mctrack = dynamic_cast<AliAODMCParticle *>(fMCArray->At(TMath::Abs(jLabel))))){
1973 Int_t grandMaPDG = mctrack->GetPdgCode();
b8bec44f 1974 if((int(TMath::Abs(grandMaPDG)/100.)%10) == kBeauty || (int(TMath::Abs(grandMaPDG)/1000.)%10) == kBeauty) {
76d0b522 1975 return kB2Jpsi;
1976 }
1977 }
1978 }
1979 return kJpsi;
1980 }
b8bec44f 1981 else if ( (int(TMath::Abs(maPdgcode)/100.)%10) == kCharm || (int(TMath::Abs(maPdgcode)/1000.)%10) == kCharm ) {
76d0b522 1982 //
1983 // charm
1984 //
1985 for (Int_t i=0; i<fNparents; i++){
b8bec44f 1986 if (TMath::Abs(maPdgcode)==fParentSelect[0][i]){
76d0b522 1987 isFinalOpenCharm = kTRUE;
1988 }
1989 }
1990 if (!isFinalOpenCharm) {
1991 return -1;
1992 }
1993
1994 // iterate until you find B hadron as a mother or become top ancester
1995 for (Int_t i=1; i<fgkMaxIter; i++){
1996
1997 Int_t jLabel = partMother->GetMother();
1998 if (jLabel == -1){
1999 return kDirectCharm;
2000 }
2001 if ((jLabel<0) || (jLabel>=fMCArray->GetEntriesFast())){
2002 AliDebug(1, "Stack label is negative, return\n");
2003 return -1;
2004 }
2005
2006 // if there is an ancester
2007 if(!(mctrack = dynamic_cast<AliAODMCParticle *>(fMCArray->At(TMath::Abs(jLabel))))) {
2008 return -1;
2009 }
2010 Int_t grandMaPDG = mctrack->GetPdgCode();
2011 for (Int_t j=0; j<fNparents; j++){
b8bec44f 2012 if (TMath::Abs(grandMaPDG)==fParentSelect[1][j]){
76d0b522 2013 return kBeautyCharm;
2014 }
2015 }
2016 partMother = mctrack;
2017 } // end of iteration
2018
2019 } // end of if
b8bec44f 2020 else if ( (int(TMath::Abs(maPdgcode)/100.)%10) == kBeauty || (int(TMath::Abs(maPdgcode)/1000.)%10) == kBeauty ) {
76d0b522 2021 //
2022 // beauty
2023 //
2024 for (Int_t i=0; i<fNparents; i++){
b8bec44f 2025 if (TMath::Abs(maPdgcode)==fParentSelect[1][i]){
76d0b522 2026 return kDirectBeauty;
2027 }
2028 }
2029 } // end of if
b8bec44f 2030 else if ( TMath::Abs(maPdgcode) == 22 ) {
76d0b522 2031 //
2032 //conversion
2033 //
2034 tmpMomLabel = partMotherCopy->GetMother();
2035 if((tmpMomLabel<0) || (tmpMomLabel>=fMCArray->GetEntriesFast())) {
2036 return -1;
2037 }
2038 if(!(mctrack = dynamic_cast<AliAODMCParticle *>(fMCArray->At(TMath::Abs(tmpMomLabel))))) {
2039 return -1;
2040 }
2041 partMother = mctrack;
2042 maPdgcode = partMother->GetPdgCode();
2043
2044 // check if the ligth meson is the decay product of heavy mesons
2045 tmpMomLabel = partMother->GetMother();
2046 if((tmpMomLabel>=0) && (tmpMomLabel<fMCArray->GetEntriesFast())) {
2047 if((mctrack = dynamic_cast<AliAODMCParticle *>(fMCArray->At(TMath::Abs(tmpMomLabel))))) {
2048 partMother = mctrack;
2049 grmaPdgcode = partMother->GetPdgCode();
2050
b8bec44f 2051 if ( (int(TMath::Abs(grmaPdgcode)/100.)%10) == kBeauty || (int(TMath::Abs(grmaPdgcode)/1000.)%10) == kBeauty ) {
76d0b522 2052 return kGammaB2M;
2053 }
b8bec44f 2054 if ( (int(TMath::Abs(grmaPdgcode)/100.)%10) == kCharm || (int(TMath::Abs(grmaPdgcode)/1000.)%10) == kCharm ) {
76d0b522 2055 return kGammaD2M;
2056 }
b8bec44f 2057 if ( TMath::Abs(grmaPdgcode) == 221 ) {
afb48e1d 2058 return kGammaEta2Pi0; //mj to remove secondary pi0 decay electrons from eta decays, mainly to eta signal enhance samples
2059 //return kElse; //mj to remove secondary pi0 decay electrons from eta decays, mainly to eta signal enhance samples
76d0b522 2060 }
2061
2062 tmpMomLabel = partMother->GetMother();
2063 if((tmpMomLabel>=0) && (tmpMomLabel<fMCArray->GetEntriesFast())) {
2064 if((mctrack = dynamic_cast<AliAODMCParticle *>(fMCArray->At(TMath::Abs(tmpMomLabel))))) {
2065 partMother = mctrack;
2066 ggrmaPdgcode = partMother->GetPdgCode();
2067
b8bec44f 2068 if ( (int(TMath::Abs(ggrmaPdgcode)/100.)%10) == kBeauty || (int(TMath::Abs(ggrmaPdgcode)/1000.)%10) == kBeauty ) {
76d0b522 2069 return kGammaB2M;
2070 }
b8bec44f 2071 if ( (int(TMath::Abs(ggrmaPdgcode)/100.)%10) == kCharm || (int(TMath::Abs(ggrmaPdgcode)/1000.)%10) == kCharm ) {
76d0b522 2072 return kGammaD2M;
2073 }
2074 }
2075 }
2076 }
2077 }
2078
b8bec44f 2079 if ( TMath::Abs(maPdgcode) == 111 ) {
76d0b522 2080 return kGammaPi0;
2081 }
b8bec44f 2082 else if ( TMath::Abs(maPdgcode) == 221 ) {
76d0b522 2083 return kGammaEta;
2084 }
b8bec44f 2085 else if ( TMath::Abs(maPdgcode) == 223 ) {
76d0b522 2086 return kGammaOmega;
2087 }
b8bec44f 2088 else if ( TMath::Abs(maPdgcode) == 333 ) {
76d0b522 2089 return kGammaPhi;
2090 }
b8bec44f 2091 else if ( TMath::Abs(maPdgcode) == 331 ) {
76d0b522 2092 return kGammaEtaPrime;
2093 }
b8bec44f 2094 else if ( TMath::Abs(maPdgcode) == 113 ) {
76d0b522 2095 return kGammaRho0;
2096 }
2097 else origin = kElse;
2098
2099 return origin;
2100
2101 }
2102 else {
2103 //
2104 // check if the ligth meson is the decay product of heavy mesons
2105 //
2106 tmpMomLabel = partMotherCopy->GetMother();
2107 if((tmpMomLabel>=0) && (tmpMomLabel<fMCArray->GetEntriesFast())) {
2108 if((mctrack = dynamic_cast<AliAODMCParticle *>(fMCArray->At(TMath::Abs(tmpMomLabel))))) {
2109 partMother = mctrack;
2110 grmaPdgcode = partMother->GetPdgCode();
2111
b8bec44f 2112 if ( (int(TMath::Abs(grmaPdgcode)/100.)%10) == kBeauty || (int(TMath::Abs(grmaPdgcode)/1000.)%10) == kBeauty ) {
afb48e1d 2113 return kB2M;
76d0b522 2114 }
b8bec44f 2115 if ( (int(TMath::Abs(grmaPdgcode)/100.)%10) == kCharm || (int(TMath::Abs(grmaPdgcode)/1000.)%10) == kCharm ) {
afb48e1d 2116 return kD2M;
76d0b522 2117 }
b8bec44f 2118 if ( TMath::Abs(grmaPdgcode) == 221 ) {
afb48e1d 2119 return kEta2Pi0; //mj to remove secondary pi0 decay electrons from eta decays, mainly to eta signal enhance samples
76d0b522 2120 }
2121
2122 tmpMomLabel = partMother->GetMother();
2123 if((tmpMomLabel>=0) && (tmpMomLabel<fMCArray->GetEntriesFast())) {
2124 if((mctrack = dynamic_cast<AliAODMCParticle *>(fMCArray->At(TMath::Abs(tmpMomLabel))))) {
2125 partMother = mctrack;
2126 ggrmaPdgcode = partMother->GetPdgCode();
2127
b8bec44f 2128 if ( (int(TMath::Abs(ggrmaPdgcode)/100.)%10) == kBeauty || (int(TMath::Abs(ggrmaPdgcode)/1000.)%10) == kBeauty ) {
afb48e1d 2129 return kB2M;
76d0b522 2130 }
b8bec44f 2131 if ( (int(TMath::Abs(ggrmaPdgcode)/100.)%10) == kCharm || (int(TMath::Abs(ggrmaPdgcode)/1000.)%10) == kCharm ) {
afb48e1d 2132 return kD2M;
76d0b522 2133 }
2134 }
2135 }
2136 }
2137
b8bec44f 2138 if ( TMath::Abs(maPdgcode) == 111 ) {
76d0b522 2139 return kPi0;
2140 }
b8bec44f 2141 else if ( TMath::Abs(maPdgcode) == 221 ) {
76d0b522 2142 return kEta;
2143 }
b8bec44f 2144 else if ( TMath::Abs(maPdgcode) == 223 ) {
76d0b522 2145 return kOmega;
2146 }
b8bec44f 2147 else if ( TMath::Abs(maPdgcode) == 333 ) {
76d0b522 2148 return kPhi;
2149 }
b8bec44f 2150 else if ( TMath::Abs(maPdgcode) == 331 ) {
76d0b522 2151 return kEtaPrime;
2152 }
b8bec44f 2153 else if ( TMath::Abs(maPdgcode) == 113 ) {
76d0b522 2154 return kRho0;
2155 }
b8bec44f 2156 else if ( TMath::Abs(maPdgcode) == 321 || TMath::Abs(maPdgcode) == 130 ) {
76d0b522 2157 return kKe3;
2158 }
2159 else{
2160 origin = kElse;
2161 }
2162 }
2163 }
2164 return origin;
2165}
9250ffbf 2166//__________________________________________
5a516e0a 2167Double_t AliHFEmcQA::GetWeightFactor(AliMCParticle *mctrack, Int_t iBgLevel){
9250ffbf 2168 //
8c1c76e9 2169 // Get weighting factor for the realistic background estimation, for three possible background yield levels, indicated by the argument "iLevel": the best estimate (0), the lower uncertainty level (1), and the upper uncertainty level (2)
9250ffbf 2170 //
2171 AliMCParticle *mctrackmother = NULL;
2172 Double_t weightElecBg = 0.;
2173 Double_t mesonPt = 0.;
2174 Double_t bgcategory = 0.;
2175 Int_t mArr = -1;
2176 Int_t mesonID = GetElecSource(mctrack->Particle());
2177 if(mesonID==kGammaPi0 || mesonID==kPi0) mArr=0; //pion
2178 else if(mesonID==kGammaEta || mesonID==kEta) mArr=1; //eta
2179 else if(mesonID==kGammaOmega || mesonID==kOmega) mArr=2; //omega
2180 else if(mesonID==kGammaPhi || mesonID==kPhi) mArr=3; //phi
2181 else if(mesonID==kGammaEtaPrime || mesonID==kEtaPrime) mArr=4; //etaprime
2182 else if(mesonID==kGammaRho0 || mesonID==kRho0) mArr=5; //rho
2183
11ff28c5 2184 Double_t datamc[25]={-999,-999,-999,-999,-999,-999,-999,-999,-999,-999,-999,-999,-999,-999,-999,-999, -999, -999, -999, -999, -999, -999, -999, -999, -999};
a8ef1999 2185 Double_t xr[3]={-999,-999,-999};
2186 datamc[0] = mesonID;
2187 datamc[17] = mctrack->Pt(); //electron pt
2188 datamc[18] = mctrack->Eta(); //electron eta
2189
2190 mctrack->XvYvZv(xr);
2191 datamc[9] = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
2192 datamc[10] = xr[2];
2193
2194 TParticle *mcpart = mctrack->Particle();
2195 if(mcpart){
2196 datamc[14] = mcpart->GetUniqueID();
2197 }
2198
9250ffbf 2199 if(!(mArr<0)){
2200 if(mesonID>=kGammaPi0) { // conversion electron, be careful with the enum odering
2201 Int_t glabel=TMath::Abs(mctrack->GetMother()); // gamma label
2202 if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
2203 glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's label
2204 if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
2205 mesonPt = mctrackmother->Pt(); //meson pt
2206 bgcategory = 1.;
a8ef1999 2207 datamc[1] = bgcategory;
2208 datamc[2] = mesonPt;
2209 mctrackmother->XvYvZv(xr);
2210 datamc[11] = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
2211 datamc[12] = xr[2];
2212
2213 mcpart = mctrackmother->Particle();
2214 if(mcpart){
2215 datamc[15] = mcpart->GetUniqueID();
2216 }
2217 if(glabel>fMCEvent->GetNumberOfPrimaries()) {
2218 bgcategory = 2.;
2219 datamc[1] = bgcategory;
2220 //printf("I should be gamma meson = %d mesonlabel= %d NumberOfPrimaries= %d \n",mctrackmother->PdgCode(),glabel,fMCEvent->GetNumberOfPrimaries());
2221 glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother
2222 if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
2223 datamc[3]=mctrackmother->PdgCode();
2224 datamc[4]=mctrackmother->Pt();
2225 if(TMath::Abs(mctrackmother->PdgCode())==310){
2226 bgcategory = 3.;
2227 glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother's mother
2228 if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
2229 datamc[5]=mctrackmother->PdgCode();
2230 glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother
2231 if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
2232 datamc[6]=mctrackmother->PdgCode();
2233 }
2234 }
2235 }
2236 }
2237 }
9250ffbf 2238 }
2239 }
2240 }
2241 else{ // nonHFE except for the conversion electron
2242 Int_t glabel=TMath::Abs(mctrack->GetMother());
2243 if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
2244 mesonPt=mctrackmother->Pt(); //meson pt
afb48e1d 2245 if(mesonID==kEta) bgcategory = -1.41; // to consider new branching ratio for the eta Dalitz decay
2246 else bgcategory = -1.;
a8ef1999 2247 datamc[1] = bgcategory;
2248 datamc[2] = mesonPt;
2249 mctrackmother->XvYvZv(xr);
2250 datamc[11] = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
2251 datamc[12] = xr[2];
2252
2253 mcpart = mctrackmother->Particle();
2254 if(mcpart){
2255 datamc[15] = mcpart->GetUniqueID();
2256 }
2257 if(glabel>fMCEvent->GetNumberOfPrimaries()) {
afb48e1d 2258 if(mesonID==kEta) bgcategory = -2.82; // to consider new branching ratio for the eta Dalitz decay (1.41)
2259 else bgcategory = -2.;
a8ef1999 2260 datamc[1] = bgcategory;
2261 glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother
2262 if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
2263 datamc[3]=mctrackmother->PdgCode();
2264 datamc[4]=mctrackmother->Pt();
2265 if(TMath::Abs(mctrackmother->PdgCode())==310){
2266 bgcategory = -3.;
2267 glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother's mother
2268 if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
2269 datamc[5]=mctrackmother->PdgCode();
2270 glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother
2271 if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
2272 datamc[6]=mctrackmother->PdgCode();
2273 }
2274 }
2275 }
2276 }
2277 }
9250ffbf 2278 }
2279 }
8c1c76e9 2280
a8ef1999 2281
8c1c76e9 2282 if(fIsPbPb){
11ff28c5 2283 Int_t centBin = GetWeightCentralityBin(fPerCentrality);
2284 if(centBin < 0)return 0.;
2285 weightElecBg=fElecBackgroundFactor[iBgLevel][centBin][mArr][kBgPtBins-1];
8c1c76e9 2286 for(int ii=0; ii<kBgPtBins; ii++){
2287 if((mesonPt > fBinLimit[ii]) && (mesonPt < fBinLimit[ii+1])){
11ff28c5 2288 weightElecBg = fElecBackgroundFactor[iBgLevel][centBin][mArr][ii];
8c1c76e9 2289 break;
2290 }
2291 }
9250ffbf 2292 }
8c1c76e9 2293 else{
2294 weightElecBg=fElecBackgroundFactor[iBgLevel][0][mArr][kBgPtBins-1];
9250ffbf 2295 for(int ii=0; ii<kBgPtBins; ii++){
8c1c76e9 2296 if((mesonPt > fBinLimit[ii]) && (mesonPt < fBinLimit[ii+1])){
2297 weightElecBg = fElecBackgroundFactor[iBgLevel][0][mArr][ii];
2298 break;
2299 }
9250ffbf 2300 }
8c1c76e9 2301 }
9250ffbf 2302 }
a8ef1999 2303
2304 datamc[13] = weightElecBg;
2305 datamc[16] = Double_t(fContainerStep);
2306
2307 datamc[7] = fHfeImpactR;
2308 datamc[8] = fHfeImpactnsigmaR;
2309
2310 datamc[19] = fRecPt;
2311 datamc[20] = fRecEta;
2312 datamc[21] = fRecPhi;
2313 datamc[22] = fLyrhit;
2314 datamc[23] = fLyrstat;
11ff28c5 2315 datamc[24] = fCentrality;
a8ef1999 2316
2317 if(fIsDebugStreamerON && fTreeStream){
2318 if(!iBgLevel){
2319 (*fTreeStream)<<"nonhfeQA"<<
2320 "mesonID="<<datamc[0]<<
2321 "bgcategory="<<datamc[1]<<
2322 "mesonPt="<<datamc[2]<<
2323 "mesonMomPdg="<<datamc[3]<<
2324 "mesonMomPt="<<datamc[4]<<
2325 "mesonGMomPdg="<<datamc[5]<<
2326 "mesonGGMomPdg="<<datamc[6]<<
2327 "eIPAbs="<<datamc[7]<<
2328 "eIPSig="<<datamc[8]<<
2329 "eR="<<datamc[9]<<
2330 "eZ="<<datamc[10]<<
2331 "mesonR="<<datamc[11]<<
2332 "mesonZ="<<datamc[12]<<
2333 "weightElecBg="<<datamc[13]<<
2334 "eUniqID="<<datamc[14]<<
2335 "mesonUniqID="<<datamc[15]<<
2336 "containerStep="<<datamc[16]<<
2337 "emcpt="<<datamc[17]<<
2338 "emceta="<<datamc[18]<<
2339 "erecpt="<<datamc[19]<<
2340 "ereceta="<<datamc[20]<<
2341 "erecphi="<<datamc[21]<<
2342 "itshit="<<datamc[22]<<
11ff28c5 2343 "itsstat="<<datamc[23]<<
2344 "centrality="<<datamc[24]
a8ef1999 2345 << "\n";
2346 }
2347 }
2348
2349 Double_t returnval = bgcategory*weightElecBg;
2350 if(TMath::Abs(bgcategory)>1) returnval = bgcategory/2.;
2351
2352 return returnval;
9250ffbf 2353}
2354
11ff28c5 2355//__________________________________________
5a516e0a 2356Double_t AliHFEmcQA::GetWeightFactor(const AliAODMCParticle * const mcpart, Int_t iBgLevel){
afb48e1d 2357 //
2358 // Get weighting factor for the realistic background estimation, for three possible background yield levels, indicated by the argument "iLevel": the best estimate (0), the lower uncertainty level (1), and the upper uncertainty level (2)
2359 //
2360
2361 if (!mcpart) return 0;
2362 if (!fMCArray) return 0;
2363
2364 Double_t weightElecBg = 0.;
2365 Double_t mesonPt = 0.;
2366 Double_t bgcategory = 0.;
2367 Int_t mArr = -1;
2368 Int_t mesonID = GetElecSource(mcpart);
2369 if(mesonID==kGammaPi0 || mesonID==kPi0) mArr=0; //pion
2370 else if(mesonID==kGammaEta || mesonID==kEta) mArr=1; //eta
2371 else if(mesonID==kGammaOmega || mesonID==kOmega) mArr=2; //omega
2372 else if(mesonID==kGammaPhi || mesonID==kPhi) mArr=3; //phi
2373 else if(mesonID==kGammaEtaPrime || mesonID==kEtaPrime) mArr=4; //etaprime
2374 else if(mesonID==kGammaRho0 || mesonID==kRho0) mArr=5; //rho
2375
2376 if(!(mArr<0)){
2377
2378 AliAODMCParticle *mctrackmother = NULL; // will change all the time
2379
2380 if(mesonID>=kGammaPi0) { // conversion electron, be careful with the enum odering
2381 Int_t iLabel = mcpart->GetMother(); //gamma label
2382 if(!(mctrackmother = dynamic_cast<AliAODMCParticle *>(fMCArray->At(TMath::Abs(iLabel))))) return 0;
2383 iLabel = mctrackmother->GetMother(); //gamma's mother's label
2384 if(!(mctrackmother= dynamic_cast<AliAODMCParticle *>(fMCArray->At(TMath::Abs(iLabel))))) return 0;
2385 mesonPt = mctrackmother->Pt(); //meson pt
2386 bgcategory = 1.;
2387
2388 //if(glabel>fMCEvent->GetNumberOfPrimaries()) bgcategory = 2.;
2389 }
2390 else{ // nonHFE except for the conversion electron
2391 Int_t iLabel = mcpart->GetMother(); //meson label
2392 if(!(mctrackmother = dynamic_cast<AliAODMCParticle *>(fMCArray->At(TMath::Abs(iLabel))))) return 0;
2393 mesonPt = mctrackmother->Pt(); //meson pt
2394 if(mesonID==kEta) bgcategory = -1.41; // to consider new branching ratio for the eta Dalitz decay
2395 else bgcategory = -1.;
2396
2397 //if(glabel>fMCEvent->GetNumberOfPrimaries()) bgcategory = -2.;
2398 }
2399
2400 if(fIsPbPb){
2401 Int_t centBin = GetWeightCentralityBin(fPerCentrality);
2402 if(centBin < 0)return 0.;
2403 weightElecBg=fElecBackgroundFactor[iBgLevel][centBin][mArr][kBgPtBins-1];
2404 for(int ii=0; ii<kBgPtBins; ii++){
2405 if((mesonPt > fBinLimit[ii]) && (mesonPt < fBinLimit[ii+1])){
2406 weightElecBg = fElecBackgroundFactor[iBgLevel][centBin][mArr][ii];
2407 break;
2408 }
2409 }
2410 }
2411 else{
2412 weightElecBg=fElecBackgroundFactor[iBgLevel][0][mArr][kBgPtBins-1];
2413 for(int ii=0; ii<kBgPtBins; ii++){
2414 if((mesonPt > fBinLimit[ii]) && (mesonPt < fBinLimit[ii+1])){
2415 weightElecBg = fElecBackgroundFactor[iBgLevel][0][mArr][ii];
2416 break;
2417 }
2418 }
2419 }
2420 }
2421
2422 Double_t returnval = bgcategory*weightElecBg;
2423 if(TMath::Abs(bgcategory)>1) returnval = bgcategory/2.;
2424
2425 return returnval;
2426}
2427
2428//__________________________________________
2429Int_t AliHFEmcQA::GetMother(const AliVParticle * const mcpart) const {
11ff28c5 2430 //
2431 // Wrapper to get the mother label
2432 //
2433 Int_t label = -1;
2434 TClass *type = mcpart->IsA();
2435 if(type == AliMCParticle::Class()){
2436 // ESD analysis
2437 const AliMCParticle *emcpart = static_cast<const AliMCParticle *>(mcpart);
2438 label = emcpart->GetMother();
2439 } else if(type == AliAODMCParticle::Class()){
2440 // AOD analysis
2441 const AliAODMCParticle *amcpart = static_cast<const AliAODMCParticle *>(mcpart);
2442 label = amcpart->GetMother();
2443 }
2444 return label;
2445}
2446//__________________________________________
5a516e0a 2447Int_t AliHFEmcQA::GetWeightCentralityBin(Float_t percentile) const {
11ff28c5 2448 //
2449 //translate the centrality percentile into the centrality bin of the reference weighting histograms for electron background
2450 //
2451
2452 Float_t centralityLimits[12]= {0.,5.,10., 20., 30., 40., 50., 60.,70.,80., 90., 100.};
2453 Int_t bin = -1;
2454 for(Int_t ibin = 0; ibin < 11; ibin++){
2455 if(percentile >= centralityLimits[ibin] && percentile < centralityLimits[ibin+1]){
2456 bin = ibin;
2457 break;
2458 }
2459 }
2460 return bin;
2461}
70da6c5a 2462//__________________________________________
2463void AliHFEmcQA::AliHists::FillList(TList *l) const {
2464 //
2465 // Fill Histos into a list for output
2466 //
2467 if(fPdgCode) l->Add(fPdgCode);
2468 if(fPt) l->Add(fPt);
2469 if(fY) l->Add(fY);
2470 if(fEta) l->Add(fEta);
2471}
2472
2473//__________________________________________
2474void AliHFEmcQA::AliHistsComm::FillList(TList *l) const {
2475 //
2476 // Fill Histos into a list for output
2477 //
2478 if(fNq) l->Add(fNq);
2479 if(fProcessID) l->Add(fProcessID);
2480 if(fePtRatio) l->Add(fePtRatio);
ccc37cdc 2481 if(fPtCorr) l->Add(fPtCorr);
c2690925 2482 if(fPtCorrDp) l->Add(fPtCorrDp);
2483 if(fPtCorrD0) l->Add(fPtCorrD0);
2484 if(fPtCorrDrest) l->Add(fPtCorrDrest);
cedf0381 2485
2486 if(fPtCorrDinein) l->Add(fPtCorrDinein);
2487 if(fPtCorrDineout) l->Add(fPtCorrDineout);
2488 if(fPtCorrDoutein) l->Add(fPtCorrDoutein);
2489 if(fPtCorrDouteout) l->Add(fPtCorrDouteout);
2490 if(fPtCorrDpDinein) l->Add(fPtCorrDpDinein);
2491 if(fPtCorrDpDineout) l->Add(fPtCorrDpDineout);
2492 if(fPtCorrDpDoutein) l->Add(fPtCorrDpDoutein);
2493 if(fPtCorrDpDouteout) l->Add(fPtCorrDpDouteout);
2494 if(fPtCorrD0Dinein) l->Add(fPtCorrD0Dinein);
2495 if(fPtCorrD0Dineout) l->Add(fPtCorrD0Dineout);
2496 if(fPtCorrD0Doutein) l->Add(fPtCorrD0Doutein);
2497 if(fPtCorrD0Douteout) l->Add(fPtCorrD0Douteout);
2498 if(fPtCorrDrestDinein) l->Add(fPtCorrDrestDinein);
2499 if(fPtCorrDrestDineout) l->Add(fPtCorrDrestDineout);
2500 if(fPtCorrDrestDoutein) l->Add(fPtCorrDrestDoutein);
2501 if(fPtCorrDrestDouteout) l->Add(fPtCorrDrestDouteout);
2502
2503 if(fEtaCorrD) l->Add(fEtaCorrD);
2504 if(fEtaCorrDp) l->Add(fEtaCorrDp);
2505 if(fEtaCorrD0) l->Add(fEtaCorrD0);
2506 if(fEtaCorrDrest) l->Add(fEtaCorrDrest);
2507
2508 if(fEtaCorrGD) l->Add(fEtaCorrGD);
2509 if(fEtaCorrGDp) l->Add(fEtaCorrGDp);
2510 if(fEtaCorrGD0) l->Add(fEtaCorrGD0);
2511 if(fEtaCorrGDrest) l->Add(fEtaCorrGDrest);
2512
2513 if(fEtaCorrB) l->Add(fEtaCorrB);
2514 if(fEtaCorrGB) l->Add(fEtaCorrGB);
2515 if(fPtCorrBinein) l->Add(fPtCorrBinein);
2516 if(fPtCorrBineout) l->Add(fPtCorrBineout);
2517 if(fPtCorrBoutein) l->Add(fPtCorrBoutein);
2518 if(fPtCorrBouteout) l->Add(fPtCorrBouteout);
2519
70da6c5a 2520 if(fDePtRatio) l->Add(fDePtRatio);
2521 if(feDistance) l->Add(feDistance);
2522 if(fDeDistance) l->Add(fDeDistance);
2523}
2524