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