]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGHF/hfe/AliHFEmcQA.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
updated
[u/mrichter/AliRoot.git] / PWGHF / hfe / AliHFEmcQA.cxx
CommitLineData
259c3296 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
50685501 15//
16// QA class of Heavy Flavor quark and fragmeted/decayed particles
17// -Check kinematics of Heavy Quarks/hadrons, and decayed leptons
18// pT, rapidity
19// decay lepton kinematics w/wo acceptance
20// heavy hadron decay length, electron pT fraction carried from decay
21// -Check yield of Heavy Quarks/hadrons
22// Number of produced heavy quark
23// Number of produced hadron of given pdg code
24//
25//
26// Authors:
27// MinJung Kweon <minjung@physi.uni-heidelberg.de>
28//
259c3296 29
259c3296 30#include <TH2F.h>
259c3296 31#include <TH1F.h>
32#include <TH2F.h>
70da6c5a 33#include <TList.h>
259c3296 34#include <TParticle.h>
a8ef1999 35#include "TTreeStream.h"
259c3296 36
37#include <AliLog.h>
faee3b18 38#include <AliMCEvent.h>
9bcfd1ab 39#include <AliGenEventHeader.h>
0792aa82 40#include <AliAODMCParticle.h>
c2690925 41#include <AliStack.h>
259c3296 42
43#include "AliHFEmcQA.h"
70da6c5a 44#include "AliHFEtools.h"
c2690925 45#include "AliHFEcollection.h"
259c3296 46
259c3296 47ClassImp(AliHFEmcQA)
48
49//_______________________________________________________________________________________________
8c1c76e9 50 AliHFEmcQA::AliHFEmcQA() :
51 fMCEvent(NULL)
70da6c5a 52 ,fMCHeader(NULL)
53 ,fMCArray(NULL)
54 ,fQAhistos(NULL)
c2690925 55 ,fMCQACollection(NULL)
8c1c76e9 56 ,fNparents(0)
57 ,fCentrality(0)
11ff28c5 58 ,fPerCentrality(-1)
8c1c76e9 59 ,fIsPbPb(kFALSE)
60 ,fIsppMultiBin(kFALSE)
a8ef1999 61 ,fContainerStep(0)
62 ,fIsDebugStreamerON(kFALSE)
63 ,fRecPt(-999)
64 ,fRecEta(-999)
65 ,fRecPhi(-999)
66 ,fLyrhit(0)
67 ,fLyrstat(0)
68 ,fHfeImpactR(-999)
69 ,fHfeImpactnsigmaR(-999)
70 ,fTreeStream(NULL)
259c3296 71{
72 // Default constructor
e088e9f0 73 for(Int_t mom = 0; mom < 9; mom++){
74 fhD[mom] = NULL;
e088e9f0 75 }
bf892a6a 76 for(Int_t mom = 0; mom < 50; mom++){
77 fHeavyQuark[mom] = NULL;
78 }
79 for(Int_t mom = 0; mom < 2; mom++){
80 fIsHeavy[mom] = 0;
81 }
8c1c76e9 82 memset(fElecBackgroundFactor, 0, sizeof(Double_t) * kElecBgSpecies * kBgPtBins * kCentBins * kBgLevels);
9250ffbf 83 memset(fBinLimit, 0, sizeof(Double_t) * (kBgPtBins+1));
259c3296 84}
85
86//_______________________________________________________________________________________________
87AliHFEmcQA::AliHFEmcQA(const AliHFEmcQA&p):
8c1c76e9 88TObject(p)
89 ,fMCEvent(NULL)
70da6c5a 90 ,fMCHeader(NULL)
91 ,fMCArray(NULL)
92 ,fQAhistos(p.fQAhistos)
c2690925 93 ,fMCQACollection(p.fMCQACollection)
8c1c76e9 94 ,fNparents(p.fNparents)
95 ,fCentrality(0)
11ff28c5 96 ,fPerCentrality(-1)
8c1c76e9 97 ,fIsPbPb(kFALSE)
98 ,fIsppMultiBin(kFALSE)
a8ef1999 99 ,fContainerStep(0)
100 ,fIsDebugStreamerON(kFALSE)
101 ,fRecPt(-999)
102 ,fRecEta(-999)
103 ,fRecPhi(-999)
104 ,fLyrhit(0)
105 ,fLyrstat(0)
106 ,fHfeImpactR(0)
107 ,fHfeImpactnsigmaR(0)
108 ,fTreeStream(NULL)
259c3296 109{
110 // Copy constructor
e088e9f0 111 for(Int_t mom = 0; mom < 9; mom++){
112 fhD[mom] = NULL;
e088e9f0 113 }
bf892a6a 114 for(Int_t mom = 0; mom < 50; mom++){
115 fHeavyQuark[mom] = NULL;
116 }
117 for(Int_t mom = 0; mom < 2; mom++){
118 fIsHeavy[mom] = 0;
119 }
8c1c76e9 120 memset(fElecBackgroundFactor, 0, sizeof(Double_t) * kElecBgSpecies * kBgPtBins * kCentBins * kBgLevels);
9250ffbf 121 memset(fBinLimit, 0, sizeof(Double_t) * (kBgPtBins+1));
259c3296 122}
123
124//_______________________________________________________________________________________________
125AliHFEmcQA&
126AliHFEmcQA::operator=(const AliHFEmcQA &)
127{
128 // Assignment operator
129
130 AliInfo("Not yet implemented.");
131 return *this;
132}
133
134//_______________________________________________________________________________________________
135AliHFEmcQA::~AliHFEmcQA()
136{
137 // Destructor
138
a8ef1999 139 if(fTreeStream && fIsDebugStreamerON) delete fTreeStream;
259c3296 140 AliInfo("Analysis Done.");
141}
142
143//_______________________________________________________________________________________________
50685501 144void AliHFEmcQA::PostAnalyze() const
259c3296 145{
e3fc062d 146 //
147 // Post analysis
148 //
259c3296 149}
150
9250ffbf 151//_______________________________________________________________________________________________
152void AliHFEmcQA::SetBackgroundWeightFactor(Double_t *elecBackgroundFactor, Double_t *binLimit)
153{
e156c3bb 154 //
155 // copy background weighting factors into data member
156 //
8c1c76e9 157
158 memcpy(fElecBackgroundFactor,elecBackgroundFactor,sizeof(Double_t) * kElecBgSpecies * kBgPtBins * kCentBins * kBgLevels);
159 memcpy(fBinLimit,binLimit,sizeof(Double_t) * (kBgPtBins+1));
9250ffbf 160}
161
e3fc062d 162//__________________________________________
163void AliHFEmcQA::CreatDefaultHistograms(TList * const qaList)
164{
165 //
166 // make default histograms
167 //
168
169 if(!qaList) return;
170
171 fQAhistos = qaList;
172 fQAhistos->SetName("MCqa");
173
a8ef1999 174 CreateHistograms(AliHFEmcQA::kCharm); // create histograms for charm
175 CreateHistograms(AliHFEmcQA::kBeauty); // create histograms for beauty
176 CreateHistograms(AliHFEmcQA::kOthers); // create histograms for beauty
e3fc062d 177
e17c1f86 178// prepare 2D(pt vs Y) histogram for D spectra, we consider following 9 particles
179 const Int_t nbspecies = 9;
180 TString kDspecies[nbspecies];
181 kDspecies[0]="411"; //D+
182 kDspecies[1]="421"; //D0
183 kDspecies[2]="431"; //Ds+
184 kDspecies[3]="4122"; //Lambdac+
185 kDspecies[4]="4132"; //Ksic0
186 kDspecies[5]="4232"; //Ksic+
187 kDspecies[6]="4332"; //OmegaC0
188 kDspecies[7]="413"; //D*(2010)+
189 kDspecies[8]="423"; //D*(2007)0
ccc37cdc 190
191 const Double_t kPtbound[2] = {0.1, 20.}; //bin taken for considering inclusive e analysis binning
c2690925 192 Int_t iBin[2];
e17c1f86 193 iBin[0] = 44; // bins in pt for log binning
194 iBin[1] = 23; // bins in pt for pi0 measurement binning
ccc37cdc 195 Double_t* binEdges[1];
196 binEdges[0] = AliHFEtools::MakeLogarithmicBinning(iBin[0], kPtbound[0], kPtbound[1]);
e97c2edf 197
198 // bin size is chosen to consider ALICE D measurement
e17c1f86 199 const Int_t nptbins = 15;
200 const Int_t nybins = 9;
201 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
202 Double_t ybins[nybins+1]={-7.5,-1.0,-0.9,-0.8,-0.5,0.5,0.8,0.9,1.0,7.5}; // y binning
e97c2edf 203 TString hname;
e17c1f86 204 for (Int_t iDmeson=0; iDmeson<nbspecies; iDmeson++){
e97c2edf 205 hname = "Dmeson"+kDspecies[iDmeson];
e17c1f86 206 fhD[iDmeson] = new TH2F(hname,hname+";p_{T} (GeV/c)",nptbins,xbins,nybins,ybins);
e97c2edf 207 if(fQAhistos) fQAhistos->Add(fhD[iDmeson]);
208 }
209
c2690925 210 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
211
8c1c76e9 212 Int_t kNcent;
213 if(fIsPbPb) kNcent=11;
214 else
215 {
216 if(fIsppMultiBin) kNcent=8;
217 else kNcent = 1;
218 }
219
c2690925 220 fMCQACollection = new AliHFEcollection("TaskMCQA", "MC QA histos for meason pt spectra");
8c1c76e9 221
222 for(Int_t centbin=0; centbin<kNcent; centbin++)
223 {
224 fMCQACollection->CreateTH1Farray(Form("pionspectra_centrbin%i",centbin), "pion yields: MC p_{t} ", iBin[1],kPtRange);
225 fMCQACollection->CreateTH1Farray(Form("etaspectra_centrbin%i",centbin), "eta yields: MC p_{t} ", iBin[1],kPtRange);
226 fMCQACollection->CreateTH1Farray(Form("omegaspectra_centrbin%i",centbin), "omega yields: MC p_{t} ", iBin[1],kPtRange);
227 fMCQACollection->CreateTH1Farray(Form("phispectra_centrbin%i",centbin), "phi yields: MC p_{t} ", iBin[1],kPtRange);
228 fMCQACollection->CreateTH1Farray(Form("etapspectra_centrbin%i",centbin), "etap yields: MC p_{t} ", iBin[1],kPtRange);
229 fMCQACollection->CreateTH1Farray(Form("rhospectra_centrbin%i",centbin), "rho yields: MC p_{t} ", iBin[1],kPtRange);
230
231 fMCQACollection->CreateTH1F(Form("pionspectraLog_centrbin%i",centbin), "pion yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
232 fMCQACollection->CreateTH1F(Form("etaspectraLog_centrbin%i",centbin), "eta yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
233 fMCQACollection->CreateTH1F(Form("omegaspectraLog_centrbin%i",centbin), "omega yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
234 fMCQACollection->CreateTH1F(Form("phispectraLog_centrbin%i",centbin), "phi yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
235 fMCQACollection->CreateTH1F(Form("etapspectraLog_centrbin%i",centbin), "etap yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
236 fMCQACollection->CreateTH1F(Form("rhospectraLog_centrbin%i",centbin), "rho yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
237 fMCQACollection->CreateTH1F(Form("kaonspectraLog_centrbin%i",centbin), "kaon yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
238 fMCQACollection->CreateTH1F(Form("k0LspectraLog_centrbin%i",centbin), "k0L yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
239
240 fMCQACollection->CreateTH1F(Form("piondaughters_centrbin%i",centbin), "pion yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
241 fMCQACollection->CreateTH1F(Form("etadaughters_centrbin%i",centbin), "eta yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
242 fMCQACollection->CreateTH1F(Form("omegadaughters_centrbin%i",centbin), "omega yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
243 fMCQACollection->CreateTH1F(Form("phidaughters_centrbin%i",centbin), "phi yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
244 fMCQACollection->CreateTH1F(Form("etapdaughters_centrbin%i",centbin), "etap yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
245 fMCQACollection->CreateTH1F(Form("rhodaughters_centrbin%i",centbin), "rho yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
246 }
c2690925 247
248 fQAhistos->Add(fMCQACollection->GetList());
249
a8ef1999 250 if(!fTreeStream && fIsDebugStreamerON){
251 fTreeStream = new TTreeSRedirector(Form("HFEmcqadebugTree%s.root", GetName()));
252 }
253
e3fc062d 254}
255
259c3296 256//__________________________________________
a8ef1999 257void AliHFEmcQA::CreateHistograms(const Int_t kquark)
259c3296 258{
259 // create histograms
260
faee3b18 261 if (!(kquark == kCharm || kquark == kBeauty || kquark == kOthers)) {
259c3296 262 AliDebug(1, "This task is only for heavy quark QA, return\n");
263 return;
264 }
dbe3abbe 265 Int_t iq = kquark - kCharm;
259c3296 266
267 TString kqTypeLabel[fgkqType];
dbe3abbe 268 if (kquark == kCharm){
269 kqTypeLabel[kQuark]="c";
270 kqTypeLabel[kantiQuark]="cbar";
271 kqTypeLabel[kHadron]="cHadron";
272 kqTypeLabel[keHadron]="ceHadron";
273 kqTypeLabel[kDeHadron]="nullHadron";
274 kqTypeLabel[kElectron]="ce";
275 kqTypeLabel[kElectron2nd]="nulle";
276 } else if (kquark == kBeauty){
277 kqTypeLabel[kQuark]="b";
278 kqTypeLabel[kantiQuark]="bbar";
279 kqTypeLabel[kHadron]="bHadron";
280 kqTypeLabel[keHadron]="beHadron";
281 kqTypeLabel[kDeHadron]="bDeHadron";
282 kqTypeLabel[kElectron]="be";
283 kqTypeLabel[kElectron2nd]="bce";
faee3b18 284 } else if (kquark == kOthers){
285 kqTypeLabel[kGamma-4]="gammae";
286 kqTypeLabel[kPi0-4]="pi0e";
287 kqTypeLabel[kElse-4]="elsee";
288 kqTypeLabel[kMisID-4]="miside";
259c3296 289 }
3a72645a 290
a8ef1999 291 TString kqEtaRangeLabel[fgkEtaRanges];
292 kqEtaRangeLabel[0] = "mcqa_";
293 kqEtaRangeLabel[1] = "mcqa_barrel_";
294 kqEtaRangeLabel[2] = "mcqa_unitY_";
295
3a72645a 296 const Double_t kPtbound[2] = {0.1, 20.}; //bin taken for considering inclusive e analysis binning
e17c1f86 297 const Int_t nptbinning1 = 35;
298 Int_t iBin[2];
3a72645a 299 iBin[0] = 44; // bins in pt
e17c1f86 300 iBin[1] = nptbinning1; // bins in pt
faee3b18 301 Double_t* binEdges[1];
302 binEdges[0] = AliHFEtools::MakeLogarithmicBinning(iBin[0], kPtbound[0], kPtbound[1]);
c2690925 303
e17c1f86 304 // new binning for final electron analysis
305 const Double_t kPtbinning1[nptbinning1+1] = {0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4, 1.5, 1.75, 2., 2.25, 2.5, 2.75, 3., 3.5, 4., 4.5, 5., 5.5, 6., 7., 8., 10., 12., 14., 16., 18., 20.};
c2690925 306
e17c1f86 307 const Int_t ndptbins = 500;
308 Double_t xcorrbin[ndptbins+1];
309 for (int icorrbin = 0; icorrbin< ndptbins+1; icorrbin++){
c2690925 310 xcorrbin[icorrbin]=icorrbin*0.1;
311 }
312
11ff28c5 313 Int_t fCentrmax = 0;
314 if(!fIsPbPb&&!fIsppMultiBin) fCentrmax=0;
315 else fCentrmax = kCentBins;
259c3296 316 TString hname;
faee3b18 317 if(kquark == kOthers){
a8ef1999 318 for (Int_t icut = 0; icut < fgkEtaRanges; icut++ ){
11ff28c5 319 for (Int_t iqType = 0; iqType < 4; iqType++ ){
320 for(Int_t icentr = 0; icentr < (fCentrmax+1); icentr++)
321 {
322 hname = kqEtaRangeLabel[icut]+"Pt_"+kqTypeLabel[iqType];
323 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);
324 hname = kqEtaRangeLabel[icut]+"Y_"+kqTypeLabel[iqType];
325 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);
326 hname = kqEtaRangeLabel[icut]+"Eta_"+kqTypeLabel[iqType];
327 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);
328 // Fill List
329 if(fQAhistos) fHist[iq][iqType][icut][icentr].FillList(fQAhistos);
330 }
331 }
a8ef1999 332 }
333 return;
faee3b18 334 }
a8ef1999 335 for (Int_t icut = 0; icut < fgkEtaRanges; icut++ ){
336 for (Int_t iqType = 0; iqType < fgkqType; iqType++ ){
11ff28c5 337 if (iqType < keHadron && icut > 0) continue; // don't duplicate histogram for quark and hadron
338 for(Int_t icentr = 0; icentr<(fCentrmax+1); icentr++)
339 {
340 hname = kqEtaRangeLabel[icut]+"PdgCode_"+kqTypeLabel[iqType];
341 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);
342 hname = kqEtaRangeLabel[icut]+"Pt_"+kqTypeLabel[iqType];
343 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
344 hname = kqEtaRangeLabel[icut]+"Y_"+kqTypeLabel[iqType];
345 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);
346 hname = kqEtaRangeLabel[icut]+"Eta_"+kqTypeLabel[iqType];
347 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);
348 // Fill List
349 if(fQAhistos) fHist[iq][iqType][icut][icentr].FillList(fQAhistos);
350 }
a8ef1999 351 }
259c3296 352 }
353
a8ef1999 354 for (Int_t icut = 0; icut < fgkEtaRanges; icut++ ){
355 hname = kqEtaRangeLabel[icut]+"PtCorr_"+kqTypeLabel[kQuark];
356 fHistComm[iq][icut].fPtCorr = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1);
357 hname = kqEtaRangeLabel[icut]+"PtCorrDp_"+kqTypeLabel[kQuark];
358 fHistComm[iq][icut].fPtCorrDp= new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1);
359 hname = kqEtaRangeLabel[icut]+"PtCorrD0_"+kqTypeLabel[kQuark];
360 fHistComm[iq][icut].fPtCorrD0 = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1);
361 hname = kqEtaRangeLabel[icut]+"PtCorrDrest_"+kqTypeLabel[kQuark];
362 fHistComm[iq][icut].fPtCorrDrest = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1);
363
364 hname = kqEtaRangeLabel[icut]+"ePtRatio_"+kqTypeLabel[kQuark];
365 fHistComm[iq][icut].fePtRatio = new TH2F(hname,hname+";p_{T} (GeV/c);momentum fraction",200,0,20,100,0,1);
366 hname = kqEtaRangeLabel[icut]+"DePtRatio_"+kqTypeLabel[kQuark];
367 fHistComm[iq][icut].fDePtRatio = new TH2F(hname,hname+";p_{T} (GeV/c);momentum fraction",100,0,20,100,0,1);
368 hname = kqEtaRangeLabel[icut]+"eDistance_"+kqTypeLabel[kQuark];
369 fHistComm[iq][icut].feDistance= new TH2F(hname,hname+";p_{T} (GeV/c);distance (cm)",100,0,20,200,0,2);
370 hname = kqEtaRangeLabel[icut]+"DeDistance_"+kqTypeLabel[kQuark];
371 fHistComm[iq][icut].fDeDistance= new TH2F(hname,hname+";p_{T} (GeV/c);distance (cm)",100,0,20,200,0,2);
372 if(fQAhistos) fHistComm[iq][icut].FillList(fQAhistos);
dbe3abbe 373 }
a8ef1999 374
375 hname = kqEtaRangeLabel[0]+"Nq_"+kqTypeLabel[kQuark];
376 fHistComm[iq][0].fNq = new TH1F(hname,hname,50,-0.5,49.5);
377 hname = kqEtaRangeLabel[0]+"ProcessID_"+kqTypeLabel[kQuark];
378 fHistComm[iq][0].fProcessID = new TH1F(hname,hname,21,-10.5,10.5);
e17c1f86 379
259c3296 380}
381
382//__________________________________________
383void AliHFEmcQA::Init()
384{
385 // called at begining every event
386
387 for (Int_t i=0; i<2; i++){
388 fIsHeavy[i] = 0;
389 }
390
391 fNparents = 7;
392
e17c1f86 393 fParentSelect[0][0] = 411; //D+
394 fParentSelect[0][1] = 421; //D0
395 fParentSelect[0][2] = 431; //Ds+
396 fParentSelect[0][3] = 4122; //Lambdac+
397 fParentSelect[0][4] = 4132; //Ksic0
398 fParentSelect[0][5] = 4232; //Ksic+
399 fParentSelect[0][6] = 4332; //OmegaC0
400
401 fParentSelect[1][0] = 511; //B0
402 fParentSelect[1][1] = 521; //B+
403 fParentSelect[1][2] = 531; //Bs0
404 fParentSelect[1][3] = 5122; //Lambdab0
405 fParentSelect[1][4] = 5132; //Ksib-
406 fParentSelect[1][5] = 5232; //Ksib0
407 fParentSelect[1][6] = 5332; //Omegab-
259c3296 408
11ff28c5 409
259c3296 410}
411
9250ffbf 412//__________________________________________
c2690925 413void AliHFEmcQA::GetMesonKine()
414{
415 //
416 // get meson pt spectra
417 //
418
419 AliVParticle *mctrack2 = NULL;
420 AliMCParticle *mctrack0 = NULL;
421 AliVParticle *mctrackdaugt= NULL;
422 AliMCParticle *mctrackd= NULL;
423 Int_t id1=0, id2=0;
424
8c1c76e9 425
426 if(fCentrality>11) printf("warning centrality out of histogram array limits \n");
11ff28c5 427 if(!fIsPbPb&&!fIsppMultiBin) fCentrality=0;
8c1c76e9 428
c2690925 429 for(Int_t imc = 0; imc <fMCEvent->GetNumberOfPrimaries(); imc++){
430 if(!(mctrack2 = fMCEvent->GetTrack(imc))) continue;
431 TParticle* mcpart0 = fMCEvent->Stack()->Particle(imc);
432 if(!mcpart0) continue;
433 mctrack0 = dynamic_cast<AliMCParticle *>(mctrack2);
f9f097c0 434 if(!mctrack0) continue;
8c1c76e9 435
11ff28c5 436// if(!fIsPbPb&&!fIsppMultiBin) fCentrality=0;
8c1c76e9 437
c2690925 438 if(abs(mctrack0->PdgCode()) == 111) // pi0
439 {
440 if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
8c1c76e9 441 fMCQACollection->Fill(Form("pionspectra_centrbin%i",fCentrality),mctrack0->Pt());
442 fMCQACollection->Fill(Form("pionspectraLog_centrbin%i",fCentrality),mctrack0->Pt());
c2690925 443 }
444 id1=mctrack0->GetFirstDaughter();
445 id2=mctrack0->GetLastDaughter();
446 if(!((id2-id1)==2)) continue;
447 for(int idx=id1; idx<=id2; idx++){
448 if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
9250ffbf 449 if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
c2690925 450 if(abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
8c1c76e9 451 fMCQACollection->Fill(Form("piondaughters_centrbin%i",fCentrality),mctrackd->Pt());
c2690925 452 }
453 }
454 else if(abs(mctrack0->PdgCode()) == 221) // eta
455 {
456 if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
8c1c76e9 457 fMCQACollection->Fill(Form("etaspectra_centrbin%i",fCentrality),mctrack0->Pt());
458 fMCQACollection->Fill(Form("etaspectraLog_centrbin%i",fCentrality),mctrack0->Pt());
c2690925 459 }
460 id1=mctrack0->GetFirstDaughter();
461 id2=mctrack0->GetLastDaughter();
462 if(!((id2-id1)==2||(id2-id1)==3)) continue;
463 for(int idx=id1; idx<=id2; idx++){
464 if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
9250ffbf 465 if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
c2690925 466 if(abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
8c1c76e9 467 fMCQACollection->Fill(Form("etadaughters_centrbin%i",fCentrality),mctrackd->Pt());
c2690925 468 }
469 }
470 else if(abs(mctrack0->PdgCode()) == 223) // omega
471 {
472 if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
8c1c76e9 473 fMCQACollection->Fill(Form("omegaspectra_centrbin%i",fCentrality),mctrack0->Pt());
474 fMCQACollection->Fill(Form("omegaspectraLog_centrbin%i",fCentrality),mctrack0->Pt());
c2690925 475 }
476 id1=mctrack0->GetFirstDaughter();
477 id2=mctrack0->GetLastDaughter();
478 if(!((id2-id1)==1||(id2-id1)==2)) continue;
479 for(int idx=id1; idx<=id2; idx++){
480 if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
9250ffbf 481 if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
c2690925 482 if(abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
8c1c76e9 483 fMCQACollection->Fill(Form("omegadaughters_centrbin%i",fCentrality),mctrackd->Pt());
c2690925 484 }
485 }
486 else if(abs(mctrack0->PdgCode()) == 333) // phi
487 {
488 if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
8c1c76e9 489 fMCQACollection->Fill(Form("phispectra_centrbin%i",fCentrality),mctrack0->Pt());
490 fMCQACollection->Fill(Form("phispectraLog_centrbin%i",fCentrality),mctrack0->Pt());
c2690925 491 }
492 id1=mctrack0->GetFirstDaughter();
493 id2=mctrack0->GetLastDaughter();
494 if(!((id2-id1)==1)) continue;
495 for(int idx=id1; idx<=id2; idx++){
496 if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
9250ffbf 497 if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
c2690925 498 if(abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
8c1c76e9 499 fMCQACollection->Fill(Form("phidaughters_centrbin%i",fCentrality),mctrackd->Pt());
c2690925 500 }
501 }
502 else if(abs(mctrack0->PdgCode()) == 331) // eta prime
503 {
504 if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
8c1c76e9 505 fMCQACollection->Fill(Form("etapspectra_centrbin%i",fCentrality),mctrack0->Pt());
506 fMCQACollection->Fill(Form("etapspectraLog_centrbin%i",fCentrality),mctrack0->Pt());
c2690925 507 }
508 id1=mctrack0->GetFirstDaughter();
509 id2=mctrack0->GetLastDaughter();
510 if(!((id2-id1)==2||(id2-id1)==3)) continue;
511 for(int idx=id1; idx<=id2; idx++){
512 if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
9250ffbf 513 if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
c2690925 514 if(abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
8c1c76e9 515 fMCQACollection->Fill(Form("etapdaughters_centrbin%i",fCentrality),mctrackd->Pt());
c2690925 516 }
517 }
518 else if(abs(mctrack0->PdgCode()) == 113) // rho
519 {
520 if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
8c1c76e9 521 fMCQACollection->Fill(Form("rhospectra_centrbin%i",fCentrality),mctrack0->Pt());
522 fMCQACollection->Fill(Form("rhospectraLog_centrbin%i",fCentrality),mctrack0->Pt());
c2690925 523 }
524 id1=mctrack0->GetFirstDaughter();
525 id2=mctrack0->GetLastDaughter();
526 if(!((id2-id1)==1)) continue;
527 for(int idx=id1; idx<=id2; idx++){
528 if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
9250ffbf 529 if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
c2690925 530 if(abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
8c1c76e9 531 fMCQACollection->Fill(Form("rhodaughters_centrbin%i",fCentrality),mctrackd->Pt());
532 }
533 }
534 else if(abs(mctrack0->PdgCode()) == 321) // kaon+-
535 {
536 if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
537 fMCQACollection->Fill(Form("kaonspectraLog_centrbin%i",fCentrality),mctrack0->Pt());
538 }
539 }
540 else if(abs(mctrack0->PdgCode()) == 130) // k0L
541 {
542 if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
543 fMCQACollection->Fill(Form("k0LspectraLog_centrbin%i",fCentrality),mctrack0->Pt());
c2690925 544 }
545 }
546 }
547
548}
259c3296 549//__________________________________________
722347d8 550void AliHFEmcQA::GetQuarkKine(TParticle *part, Int_t iTrack, const Int_t kquark)
259c3296 551{
552 // get heavy quark kinematics
553
dbe3abbe 554 if (kquark != kCharm && kquark != kBeauty) {
259c3296 555 AliDebug(1, "This task is only for heavy quark QA, return\n");
556 return;
557 }
dbe3abbe 558 Int_t iq = kquark - kCharm;
259c3296 559
faee3b18 560 if (iTrack < 0 || !part) {
561 AliDebug(1, "Stack label is negative or no mcparticle, return\n");
259c3296 562 return;
563 }
564
11ff28c5 565 if(!fIsPbPb&&!fIsppMultiBin) fCentrality=0;
566
faee3b18 567 AliMCParticle *mctrack = NULL;
259c3296 568 Int_t partPdgcode = TMath::Abs(part->GetPdgCode());
569
570 // select heavy hadron or not fragmented heavy quark
571 if ( int(partPdgcode/100.)==kquark || int(partPdgcode/1000.)==kquark || (partPdgcode==kquark && (part->GetNDaughters()==0 && iTrack>5)) ){
572
573 TParticle *partMother;
574 Int_t iLabel;
575
576 if (partPdgcode == kquark){ // in case of not fragmented heavy quark
577 partMother = part;
578 iLabel = iTrack;
579 } else{ // in case of heavy hadron, start to search for mother heavy parton
580 iLabel = part->GetFirstMother();
faee3b18 581 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return;
582 if (iLabel>-1) { partMother = mctrack->Particle(); }
259c3296 583 else {
584 AliDebug(1, "Stack label is negative, return\n");
585 return;
586 }
587 }
588
589 // heavy parton selection as a mother of heavy hadron
590 // if the heavy particle comes from string which is denoted as particle status 12|12|12...12|11,[PYTHIA p.60]
591 // in this case, the mother of heavy particle can be one of the fragmented parton of the string
592 // should I make a condition that partMother should be quark or diquark? -> not necessary
593 if ( abs(partMother->GetPdgCode()) == kquark || (partMother->GetStatusCode() == 11) ){
594 //if ( abs(partMother->GetPdgCode()) == kquark || (partMother->GetStatusCode() == 11 || partMother->GetStatusCode() == 12) ){
595
596 if ( abs(partMother->GetPdgCode()) != kquark ){
597 // search fragmented partons in the same string
dbe3abbe 598 Bool_t isSameString = kTRUE;
259c3296 599 for (Int_t i=1; i<fgkMaxIter; i++){
600 iLabel = iLabel - 1;
faee3b18 601 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return;
602 if (iLabel>-1) { partMother = mctrack->Particle(); }
259c3296 603 else {
604 AliDebug(1, "Stack label is negative, return\n");
605 return;
606 }
607 if ( abs(partMother->GetPdgCode()) == kquark ) break;
dbe3abbe 608 if ( partMother->GetStatusCode() != 12 ) isSameString = kFALSE;
609 if (!isSameString) return;
259c3296 610 }
611 }
612 AliDebug(1, "Can not find heavy parton of this heavy hadron in the string, return\n");
613 if (abs(partMother->GetPdgCode()) != kquark) return;
614
615 if (fIsHeavy[iq] >= 50) return;
616 fHeavyQuark[fIsHeavy[iq]] = partMother;
617 fIsHeavy[iq]++;
618
619 // fill kinematics for heavy parton
620 if (partMother->GetPdgCode() > 0) { // quark
11ff28c5 621 fHist[iq][kQuark][0][fCentrality].fPdgCode->Fill(partMother->GetPdgCode());
622 fHist[iq][kQuark][0][fCentrality].fPt->Fill(partMother->Pt());
623 fHist[iq][kQuark][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMother));
624 fHist[iq][kQuark][0][fCentrality].fEta->Fill(partMother->Eta());
259c3296 625 } else{ // antiquark
11ff28c5 626 fHist[iq][kantiQuark][0][fCentrality].fPdgCode->Fill(partMother->GetPdgCode());
627 fHist[iq][kantiQuark][0][fCentrality].fPt->Fill(partMother->Pt());
628 fHist[iq][kantiQuark][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMother));
629 fHist[iq][kantiQuark][0][fCentrality].fEta->Fill(partMother->Eta());
259c3296 630 }
631
632 } // end of heavy parton slection loop
633
634 } // end of heavy hadron or quark selection
635
636}
637
638//__________________________________________
639void AliHFEmcQA::EndOfEventAna(const Int_t kquark)
640{
641 // end of event analysis
642
dbe3abbe 643 if (kquark != kCharm && kquark != kBeauty) {
259c3296 644 AliDebug(1, "This task is only for heavy quark QA, return\n");
645 return;
646 }
dbe3abbe 647 Int_t iq = kquark - kCharm;
259c3296 648
649
650 // # of heavy quark per event
651 AliDebug(1,Form("Number of heavy quark in this event = %d \n",fIsHeavy[iq]));
dbe3abbe 652 fHistComm[iq][0].fNq->Fill(fIsHeavy[iq]);
259c3296 653
654 Int_t motherID[fgkMaxGener];
655 Int_t motherType[fgkMaxGener];
656 Int_t motherLabel[fgkMaxGener];
657 Int_t ancestorPdg[fgkMaxGener];
658 Int_t ancestorLabel[fgkMaxGener];
659
660 for (Int_t i = 0; i < fgkMaxGener; i++){ // initialization
661 motherID[i] = 0;
662 motherType[i] = 0;
663 motherLabel[i] = 0;
664 ancestorPdg[i] = 0;
665 ancestorLabel[i] = 0;
666 }
667
3a72645a 668
259c3296 669 // check history of found heavy quarks
670 for (Int_t i = 0; i < fIsHeavy[iq]; i++){
671
3a72645a 672 if(!fHeavyQuark[i]) return;
673
259c3296 674 ancestorLabel[0] = i;
675 ancestorPdg[0] = fHeavyQuark[i]->GetPdgCode();
676 ancestorLabel[1] = fHeavyQuark[i]->GetFirstMother();
677
678 AliDebug(1,Form("pdg code= %d\n",ancestorPdg[0]));
679 AliDebug(1,Form("ancestor label= %d\n",ancestorLabel[1]));
680
681 Int_t ig = 1;
682 while (ancestorLabel[ig] != -1){
683 // in case there is mother, get mother's pdg code and grandmother's label
684 IdentifyMother(ancestorLabel[ig], ancestorPdg[ig], ancestorLabel[ig+1]);
685 // if mother is still heavy, find again mother's ancestor
686 if (ancestorPdg[ig-1] == ancestorPdg[ig]) {
687 ig++;
688 continue; // if it is from same heavy
689 }
690 // if the heavy's mother is not heavy, check the mother's label to know if it comes from inital or final parton shower
691 if (IsFromInitialShower(ancestorLabel[ig],motherID[i],motherType[i],motherLabel[i])) break;
692 if (IsFromFinalParton(ancestorLabel[ig],motherID[i],motherType[i],motherLabel[i])) break;
693 // if it is not the above case, something is strange
dbe3abbe 694 ReportStrangeness(motherID[i],motherType[i],motherLabel[i]);
259c3296 695 break;
696 }
697 if (ancestorLabel[ig] == -1){ // from hard scattering
698 HardScattering(kquark, motherID[i],motherType[i], motherLabel[i]);
699 }
700
701 } // end of found heavy quark loop
702
703
704 // check process type
705 Int_t processID = 0;
706 for (Int_t i = 0; i < fIsHeavy[iq]; i++){
707 AliDebug(1,Form("Mother ID= %d type= %d label= %d\n",motherID[i],motherType[i],motherLabel[i]));
708 }
709
710
711 Int_t nheavypair = Int_t(fIsHeavy[iq]/2.);
712 for (Int_t ipair = 0; ipair < nheavypair; ipair++){
713
714 Int_t id1 = ipair*2;
715 Int_t id2 = ipair*2 + 1;
716
717 if (motherType[id1] == 2 && motherType[id2] == 2){
dbe3abbe 718 if (motherLabel[id1] == motherLabel[id2]) processID = kGluonSplitting; // gluon spliting
259c3296 719 else processID = -9;
720 }
721 else if (motherType[id1] == -1 && motherType[id2] == -1) {
722 if (motherLabel[id1] == -1 && motherLabel[id2] == -1) {
dbe3abbe 723 if (motherID[id1] == fgkGluon) processID = kPairCreationFromg; // gluon fusion
724 else processID = kPairCreationFromq; // q-qbar pair creation
259c3296 725 }
726 else processID = -8;
727 }
728 else if (motherType[id1] == -1 || motherType[id2] == -1) {
729 if ((motherLabel[id1] == -1 || motherLabel[id2] == -1) && (motherLabel[id1]*motherLabel[id2] == -2 || motherLabel[id1]*motherLabel[id2] == -3)) {
dbe3abbe 730 if(motherID[id1]*motherID[id2] == kquark*fgkGluon) processID = kFlavourExitation; // flavour exitation
731 else processID = kLightQuarkShower;
259c3296 732 }
733 else processID = -7;
734 }
735 else if (motherType[id1] == -2 || motherType[id2] == -2) {
dbe3abbe 736 if (motherLabel[id1] == motherLabel[id2]) processID = kInitialPartonShower; // initial parton shower
259c3296 737 else processID = -6;
738
739 }
740 else processID = -5;
741
742 if (nheavypair >1) AliDebug(1,Form("Multi pair found : process ID = %d\n",processID));
dbe3abbe 743 else fHistComm[iq][0].fProcessID->Fill(processID);
259c3296 744 AliDebug(1,Form("Process ID = %d\n",processID));
745 } // end of # heavy quark pair loop
746
747}
748
749//__________________________________________
722347d8 750void AliHFEmcQA::GetHadronKine(TParticle* mcpart, const Int_t kquark)
dbe3abbe 751{
752 // decay electron kinematics
753
754 if (kquark != kCharm && kquark != kBeauty) {
755 AliDebug(1, "This task is only for heavy quark QA, return\n");
756 return;
757 }
758 Int_t iq = kquark - kCharm;
759
faee3b18 760 if(!mcpart){
761 AliDebug(1, "no mc particle, return\n");
762 return;
763 }
dbe3abbe 764
765 Int_t iLabel = mcpart->GetFirstMother();
766 if (iLabel<0){
767 AliDebug(1, "Stack label is negative, return\n");
768 return;
769 }
770
11ff28c5 771 if(!fIsPbPb&&!fIsppMultiBin) fCentrality=0;
772
dbe3abbe 773 TParticle *partCopy = mcpart;
774 Int_t pdgcode = mcpart->GetPdgCode();
775 Int_t pdgcodeCopy = pdgcode;
776
faee3b18 777 AliMCParticle *mctrack = NULL;
778
dbe3abbe 779 // if the mother is charmed hadron
75d81601 780 Bool_t isDirectCharm = kFALSE;
dbe3abbe 781 if ( int(abs(pdgcode)/100.) == kCharm || int(abs(pdgcode)/1000.) == kCharm ) {
782
783 // iterate until you find B hadron as a mother or become top ancester
784 for (Int_t i=1; i<fgkMaxIter; i++){
785
786 Int_t jLabel = mcpart->GetFirstMother();
787 if (jLabel == -1){
75d81601 788 isDirectCharm = kTRUE;
dbe3abbe 789 break; // if there is no ancester
790 }
791 if (jLabel < 0){ // safety protection
792 AliDebug(1, "Stack label is negative, return\n");
793 return;
794 }
795 // if there is an ancester
faee3b18 796 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return;
797 TParticle* mother = mctrack->Particle();
dbe3abbe 798 Int_t motherPDG = mother->GetPdgCode();
799
800 for (Int_t j=0; j<fNparents; j++){
801 if (abs(motherPDG)==fParentSelect[1][j]) return; // return if this hadron is originated from b
802 }
803
804 mcpart = mother;
805 } // end of iteration
806 } // end of if
75d81601 807 if((isDirectCharm == kTRUE && kquark == kCharm) || kquark == kBeauty) {
dbe3abbe 808 for (Int_t i=0; i<fNparents; i++){
809 if (abs(pdgcodeCopy)==fParentSelect[iq][i]){
810
811 // fill hadron kinematics
11ff28c5 812 fHist[iq][kHadron][0][fCentrality].fPdgCode->Fill(pdgcodeCopy);
813 fHist[iq][kHadron][0][fCentrality].fPt->Fill(partCopy->Pt());
814 fHist[iq][kHadron][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partCopy));
815 fHist[iq][kHadron][0][fCentrality].fEta->Fill(partCopy->Eta());
dbe3abbe 816
ccc37cdc 817 if(iq==0) {
818 fhD[i]->Fill(partCopy->Pt(),AliHFEtools::GetRapidity(partCopy));
ccc37cdc 819 }
dbe3abbe 820 }
821 }
ccc37cdc 822 // I also want to store D* info to compare with D* measurement
823 if (abs(pdgcodeCopy)==413 && iq==0) { //D*+
824 fhD[7]->Fill(partCopy->Pt(),AliHFEtools::GetRapidity(partCopy));
ccc37cdc 825 }
826 if (abs(pdgcodeCopy)==423 && iq==0) { //D*0
827 fhD[8]->Fill(partCopy->Pt(),AliHFEtools::GetRapidity(partCopy));
ccc37cdc 828 }
dbe3abbe 829 } // end of if
830}
831
832//__________________________________________
a8ef1999 833void AliHFEmcQA::GetDecayedKine(TParticle* mcpart, const Int_t kquark, Int_t kdecayed)
259c3296 834{
835 // decay electron kinematics
836
faee3b18 837 if (!(kquark == kCharm || kquark == kBeauty || kquark == kOthers)){
259c3296 838 AliDebug(1, "This task is only for heavy quark QA, return\n");
839 return;
840 }
dbe3abbe 841 Int_t iq = kquark - kCharm;
50685501 842 Bool_t isFinalOpenCharm = kFALSE;
259c3296 843
faee3b18 844 if(!mcpart){
845 AliDebug(1, "no mcparticle, return\n");
846 return;
259c3296 847 }
848
11ff28c5 849 if(!fIsPbPb&&!fIsppMultiBin) fCentrality=0;
850
a8ef1999 851 Double_t eabsEta = TMath::Abs(mcpart->Eta());
852 Double_t eabsY = TMath::Abs(AliHFEtools::GetRapidity(mcpart));
853
faee3b18 854 if(kquark==kOthers){
855 Int_t esource = -1;
856 if ( abs(mcpart->GetPdgCode()) != kdecayed ) esource = kMisID-4;
857 else esource =GetSource(mcpart)-4; // return for the cases kGamma=4, kPi0=5, kElse=6
858 if(esource==0|| esource==1 || esource==2 || esource==3){
11ff28c5 859 fHist[iq][esource][0][fCentrality].fPt->Fill(mcpart->Pt());
860 fHist[iq][esource][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
861 fHist[iq][esource][0][fCentrality].fEta->Fill(mcpart->Eta());
a8ef1999 862 if(eabsEta<0.9){
11ff28c5 863 fHist[iq][esource][1][fCentrality].fPt->Fill(mcpart->Pt());
864 fHist[iq][esource][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
865 fHist[iq][esource][1][fCentrality].fEta->Fill(mcpart->Eta());
a8ef1999 866 }
867 if(eabsY<0.5){
11ff28c5 868 fHist[iq][esource][2][fCentrality].fPt->Fill(mcpart->Pt());
869 fHist[iq][esource][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
870 fHist[iq][esource][2][fCentrality].fEta->Fill(mcpart->Eta());
a8ef1999 871 }
faee3b18 872 return;
873 }
874 else {
875 AliDebug(1, "e source is out of defined ranges, return\n");
876 return;
877 }
878 }
259c3296 879
880 if ( abs(mcpart->GetPdgCode()) != kdecayed ) return;
259c3296 881
882 Int_t iLabel = mcpart->GetFirstMother();
883 if (iLabel<0){
884 AliDebug(1, "Stack label is negative, return\n");
885 return;
886 }
887
faee3b18 888 AliMCParticle *mctrack = NULL;
889 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return;
890 TParticle *partMother = mctrack->Particle();
dbe3abbe 891 TParticle *partMotherCopy = partMother;
259c3296 892 Int_t maPdgcode = partMother->GetPdgCode();
dbe3abbe 893 Int_t maPdgcodeCopy = maPdgcode;
259c3296 894
9bcfd1ab 895 // get mc primary vertex
faee3b18 896 /*
9bcfd1ab 897 TArrayF mcPrimVtx;
898 if(fMCHeader) fMCHeader->PrimaryVertex(mcPrimVtx);
899
259c3296 900 // get electron production vertex
901 TLorentzVector ePoint;
902 mcpart->ProductionVertex(ePoint);
903
9bcfd1ab 904 // calculated production vertex to primary vertex (in xy plane)
905 Float_t decayLxy = TMath::Sqrt((mcPrimVtx[0]-ePoint[0])*(mcPrimVtx[0]-ePoint[0])+(mcPrimVtx[1]-ePoint[1])*(mcPrimVtx[1]-ePoint[1]));
faee3b18 906 */
907 Float_t decayLxy = 0;
9bcfd1ab 908
259c3296 909 // if the mother is charmed hadron
dbe3abbe 910 Bool_t isMotherDirectCharm = kFALSE;
911 if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
259c3296 912
0792aa82 913 for (Int_t i=0; i<fNparents; i++){
914 if (abs(maPdgcode)==fParentSelect[0][i]){
50685501 915 isFinalOpenCharm = kTRUE;
0792aa82 916 }
917 }
50685501 918 if (!isFinalOpenCharm) return ;
0792aa82 919
259c3296 920 // iterate until you find B hadron as a mother or become top ancester
921 for (Int_t i=1; i<fgkMaxIter; i++){
922
923 Int_t jLabel = partMother->GetFirstMother();
924 if (jLabel == -1){
dbe3abbe 925 isMotherDirectCharm = kTRUE;
259c3296 926 break; // if there is no ancester
927 }
928 if (jLabel < 0){ // safety protection
929 AliDebug(1, "Stack label is negative, return\n");
930 return;
931 }
932
933 // if there is an ancester
faee3b18 934 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return;
935 TParticle* grandMa = mctrack->Particle();
259c3296 936 Int_t grandMaPDG = grandMa->GetPdgCode();
937
fc0de2a0 938 for (Int_t j=0; j<fNparents; j++){
939 if (abs(grandMaPDG)==fParentSelect[1][j]){
259c3296 940
dbe3abbe 941 if (kquark == kCharm) return;
259c3296 942 // fill electron kinematics
11ff28c5 943 fHist[iq][kElectron2nd][0][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
944 fHist[iq][kElectron2nd][0][fCentrality].fPt->Fill(mcpart->Pt());
945 fHist[iq][kElectron2nd][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
946 fHist[iq][kElectron2nd][0][fCentrality].fEta->Fill(mcpart->Eta());
259c3296 947
948 // fill mother hadron kinematics
11ff28c5 949 fHist[iq][kDeHadron][0][fCentrality].fPdgCode->Fill(grandMaPDG);
950 fHist[iq][kDeHadron][0][fCentrality].fPt->Fill(grandMa->Pt());
951 fHist[iq][kDeHadron][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(grandMa));
952 fHist[iq][kDeHadron][0][fCentrality].fEta->Fill(grandMa->Eta());
a8ef1999 953
954 if(eabsEta<0.9){
11ff28c5 955 fHist[iq][kElectron2nd][1][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
956 fHist[iq][kElectron2nd][1][fCentrality].fPt->Fill(mcpart->Pt());
957 fHist[iq][kElectron2nd][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
958 fHist[iq][kElectron2nd][1][fCentrality].fEta->Fill(mcpart->Eta());
a8ef1999 959
960 // fill mother hadron kinematics
11ff28c5 961 fHist[iq][kDeHadron][1][fCentrality].fPdgCode->Fill(grandMaPDG);
962 fHist[iq][kDeHadron][1][fCentrality].fPt->Fill(grandMa->Pt());
963 fHist[iq][kDeHadron][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(grandMa));
964 fHist[iq][kDeHadron][1][fCentrality].fEta->Fill(grandMa->Eta());
a8ef1999 965 }
966
967 if(eabsY<0.5){
11ff28c5 968 fHist[iq][kElectron2nd][2][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
969 fHist[iq][kElectron2nd][2][fCentrality].fPt->Fill(mcpart->Pt());
970 fHist[iq][kElectron2nd][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
971 fHist[iq][kElectron2nd][2][fCentrality].fEta->Fill(mcpart->Eta());
a8ef1999 972
973 // fill mother hadron kinematics
11ff28c5 974 fHist[iq][kDeHadron][2][fCentrality].fPdgCode->Fill(grandMaPDG);
975 fHist[iq][kDeHadron][2][fCentrality].fPt->Fill(grandMa->Pt());
976 fHist[iq][kDeHadron][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(grandMa));
977 fHist[iq][kDeHadron][2][fCentrality].fEta->Fill(grandMa->Eta());
a8ef1999 978 }
979
259c3296 980 // ratio between pT of electron and pT of mother B hadron
a8ef1999 981 if(grandMa->Pt()) {
982 fHistComm[iq][0].fDePtRatio->Fill(grandMa->Pt(),mcpart->Pt()/grandMa->Pt());
983 if(eabsEta<0.9){
984 fHistComm[iq][1].fDePtRatio->Fill(grandMa->Pt(),mcpart->Pt()/grandMa->Pt());
985 }
986 if(eabsY<0.5){
987 fHistComm[iq][2].fDePtRatio->Fill(grandMa->Pt(),mcpart->Pt()/grandMa->Pt());
988 }
989 }
259c3296 990
9bcfd1ab 991 // distance between electron production point and primary vertex
a8ef1999 992 fHistComm[iq][0].fDeDistance->Fill(grandMa->Pt(),decayLxy);
993 if(eabsEta<0.9){
994 fHistComm[iq][1].fDeDistance->Fill(grandMa->Pt(),decayLxy);
995 }
996 if(eabsY<0.5){
997 fHistComm[iq][2].fDeDistance->Fill(grandMa->Pt(),decayLxy);
998 }
259c3296 999 return;
1000 }
dc306130 1001 }
259c3296 1002
1003 partMother = grandMa;
1004 } // end of iteration
1005 } // end of if
dbe3abbe 1006 if((isMotherDirectCharm == kTRUE && kquark == kCharm) || kquark == kBeauty) {
259c3296 1007 for (Int_t i=0; i<fNparents; i++){
1008 if (abs(maPdgcodeCopy)==fParentSelect[iq][i]){
1009
11ff28c5 1010 fHist[iq][kElectron][0][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1011 fHist[iq][kElectron][0][fCentrality].fPt->Fill(mcpart->Pt());
1012 fHist[iq][kElectron][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1013 fHist[iq][kElectron][0][fCentrality].fEta->Fill(mcpart->Eta());
a8ef1999 1014
259c3296 1015 // fill mother hadron kinematics
11ff28c5 1016 fHist[iq][keHadron][0][fCentrality].fPdgCode->Fill(maPdgcodeCopy);
1017 fHist[iq][keHadron][0][fCentrality].fPt->Fill(partMotherCopy->Pt());
1018 fHist[iq][keHadron][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
1019 fHist[iq][keHadron][0][fCentrality].fEta->Fill(partMotherCopy->Eta());
a8ef1999 1020
1021 if(eabsEta<0.9){
11ff28c5 1022 fHist[iq][kElectron][1][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1023 fHist[iq][kElectron][1][fCentrality].fPt->Fill(mcpart->Pt());
1024 fHist[iq][kElectron][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1025 fHist[iq][kElectron][1][fCentrality].fEta->Fill(mcpart->Eta());
a8ef1999 1026
1027 // fill mother hadron kinematics
11ff28c5 1028 fHist[iq][keHadron][1][fCentrality].fPdgCode->Fill(maPdgcodeCopy);
1029 fHist[iq][keHadron][1][fCentrality].fPt->Fill(partMotherCopy->Pt());
1030 fHist[iq][keHadron][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
1031 fHist[iq][keHadron][1][fCentrality].fEta->Fill(partMotherCopy->Eta());
a8ef1999 1032 }
1033
1034 if(eabsY<0.5){
11ff28c5 1035 fHist[iq][kElectron][2][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1036 fHist[iq][kElectron][2][fCentrality].fPt->Fill(mcpart->Pt());
1037 fHist[iq][kElectron][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1038 fHist[iq][kElectron][2][fCentrality].fEta->Fill(mcpart->Eta());
a8ef1999 1039
1040 // fill mother hadron kinematics
11ff28c5 1041 fHist[iq][keHadron][2][fCentrality].fPdgCode->Fill(maPdgcodeCopy);
1042 fHist[iq][keHadron][2][fCentrality].fPt->Fill(partMotherCopy->Pt());
1043 fHist[iq][keHadron][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
1044 fHist[iq][keHadron][2][fCentrality].fEta->Fill(partMotherCopy->Eta());
a8ef1999 1045 }
259c3296 1046
ccc37cdc 1047 // ratio between pT of electron and pT of mother B or direct D hadron
a8ef1999 1048 if(partMotherCopy->Pt()) {
1049 fHistComm[iq][0].fePtRatio->Fill(partMotherCopy->Pt(),mcpart->Pt()/partMotherCopy->Pt());
1050 if(eabsEta<0.9){
1051 fHistComm[iq][1].fePtRatio->Fill(partMotherCopy->Pt(),mcpart->Pt()/partMotherCopy->Pt());
1052 }
1053 if(eabsY<0.5){
1054 fHistComm[iq][2].fePtRatio->Fill(partMotherCopy->Pt(),mcpart->Pt()/partMotherCopy->Pt());
1055 }
1056 }
1057 fHistComm[iq][0].fPtCorr->Fill(partMotherCopy->Pt(),mcpart->Pt());
1058 if(eabsEta<0.9){
1059 fHistComm[iq][1].fPtCorr->Fill(partMotherCopy->Pt(),mcpart->Pt());
1060 }
1061 if(eabsY<0.5){
1062 fHistComm[iq][2].fPtCorr->Fill(partMotherCopy->Pt(),mcpart->Pt());
1063 }
1064 if(TMath::Abs(partMotherCopy->GetPdgCode())==411) {
1065 fHistComm[iq][0].fPtCorrDp->Fill(partMotherCopy->Pt(),mcpart->Pt());
1066 if(eabsEta<0.9){
1067 fHistComm[iq][1].fPtCorrDp->Fill(partMotherCopy->Pt(),mcpart->Pt());
1068 }
1069 if(eabsY<0.5){
1070 fHistComm[iq][2].fPtCorrDp->Fill(partMotherCopy->Pt(),mcpart->Pt());
1071 }
1072 }
1073 else if(TMath::Abs(partMotherCopy->GetPdgCode())==421) {
1074 fHistComm[iq][0].fPtCorrD0->Fill(partMotherCopy->Pt(),mcpart->Pt());
1075 if(eabsEta<0.9){
1076 fHistComm[iq][1].fPtCorrD0->Fill(partMotherCopy->Pt(),mcpart->Pt());
1077 }
1078 if(eabsY<0.5){
1079 fHistComm[iq][2].fPtCorrD0->Fill(partMotherCopy->Pt(),mcpart->Pt());
1080 }
1081 }
1082 else {
1083 fHistComm[iq][0].fPtCorrDrest->Fill(partMotherCopy->Pt(),mcpart->Pt());
1084 if(eabsEta<0.9){
1085 fHistComm[iq][1].fPtCorrDrest->Fill(partMotherCopy->Pt(),mcpart->Pt());
1086 }
1087 if(eabsY<0.5){
1088 fHistComm[iq][2].fPtCorrDrest->Fill(partMotherCopy->Pt(),mcpart->Pt());
1089 }
1090 }
259c3296 1091
9bcfd1ab 1092 // distance between electron production point and primary vertex
a8ef1999 1093 fHistComm[iq][0].feDistance->Fill(partMotherCopy->Pt(),decayLxy);
1094 if(eabsEta<0.9){
1095 fHistComm[iq][1].feDistance->Fill(partMotherCopy->Pt(),decayLxy);
1096 }
1097 if(eabsY<0.5){
1098 fHistComm[iq][2].feDistance->Fill(partMotherCopy->Pt(),decayLxy);
1099 }
259c3296 1100 }
1101 }
1102 } // end of if
1103}
1104
0792aa82 1105//____________________________________________________________________
a8ef1999 1106void AliHFEmcQA::GetDecayedKine(AliAODMCParticle *mcpart, const Int_t kquark, Int_t kdecayed)
0792aa82 1107{
1108 // decay electron kinematics
1109
1110 if (kquark != kCharm && kquark != kBeauty) {
1111 AliDebug(1, "This task is only for heavy quark QA, return\n");
1112 return;
1113 }
1114
1115 Int_t iq = kquark - kCharm;
50685501 1116 Bool_t isFinalOpenCharm = kFALSE;
0792aa82 1117
faee3b18 1118 if(!mcpart){
1119 AliDebug(1, "no mcparticle, return\n");
1120 return;
1121 }
1122
0792aa82 1123 if ( abs(mcpart->GetPdgCode()) != kdecayed ) return;
1124
a8ef1999 1125 Double_t eabsEta = TMath::Abs(mcpart->Eta());
1126 Double_t eabsY = TMath::Abs(AliHFEtools::GetRapidity(mcpart));
1127
0792aa82 1128 // mother
1129 Int_t iLabel = mcpart->GetMother();
1130 if (iLabel<0){
1131 AliDebug(1, "Stack label is negative, return\n");
1132 return;
1133 }
1134
11ff28c5 1135 if(!fIsPbPb&&!fIsppMultiBin) fCentrality=0;
1136
0792aa82 1137 AliAODMCParticle *partMother = (AliAODMCParticle*)fMCArray->At(iLabel);
1138 AliAODMCParticle *partMotherCopy = partMother;
1139 Int_t maPdgcode = partMother->GetPdgCode();
1140 Int_t maPdgcodeCopy = maPdgcode;
1141
1142 Bool_t isMotherDirectCharm = kFALSE;
1143 if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
1144
1145 for (Int_t i=0; i<fNparents; i++){
1146 if (abs(maPdgcode)==fParentSelect[0][i]){
50685501 1147 isFinalOpenCharm = kTRUE;
0792aa82 1148 }
1149 }
50685501 1150 if (!isFinalOpenCharm) return;
0792aa82 1151
1152 for (Int_t i=1; i<fgkMaxIter; i++){
1153
1154 Int_t jLabel = partMother->GetMother();
1155 if (jLabel == -1){
1156 isMotherDirectCharm = kTRUE;
1157 break; // if there is no ancester
1158 }
1159 if (jLabel < 0){ // safety protection
1160 AliDebug(1, "Stack label is negative, return\n");
1161 return;
1162 }
1163
1164 // if there is an ancester
1165 AliAODMCParticle* grandMa = (AliAODMCParticle*)fMCArray->At(jLabel);
1166 Int_t grandMaPDG = grandMa->GetPdgCode();
1167
1168 for (Int_t j=0; j<fNparents; j++){
1169 if (abs(grandMaPDG)==fParentSelect[1][j]){
1170
1171 if (kquark == kCharm) return;
1172 // fill electron kinematics
11ff28c5 1173 fHist[iq][kElectron2nd][0][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1174 fHist[iq][kElectron2nd][0][fCentrality].fPt->Fill(mcpart->Pt());
1175 fHist[iq][kElectron2nd][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1176 fHist[iq][kElectron2nd][0][fCentrality].fEta->Fill(mcpart->Eta());
0792aa82 1177
1178 // fill mother hadron kinematics
11ff28c5 1179 fHist[iq][kDeHadron][0][fCentrality].fPdgCode->Fill(grandMaPDG);
1180 fHist[iq][kDeHadron][0][fCentrality].fPt->Fill(grandMa->Pt());
1181 fHist[iq][kDeHadron][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(grandMa));
1182 fHist[iq][kDeHadron][0][fCentrality].fEta->Fill(grandMa->Eta());
a8ef1999 1183
1184 if(eabsEta<0.9){
1185 // fill electron kinematics
11ff28c5 1186 fHist[iq][kElectron2nd][1][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1187 fHist[iq][kElectron2nd][1][fCentrality].fPt->Fill(mcpart->Pt());
1188 fHist[iq][kElectron2nd][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1189 fHist[iq][kElectron2nd][1][fCentrality].fEta->Fill(mcpart->Eta());
a8ef1999 1190
1191 // fill mother hadron kinematics
11ff28c5 1192 fHist[iq][kDeHadron][1][fCentrality].fPdgCode->Fill(grandMaPDG);
1193 fHist[iq][kDeHadron][1][fCentrality].fPt->Fill(grandMa->Pt());
1194 fHist[iq][kDeHadron][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(grandMa));
1195 fHist[iq][kDeHadron][1][fCentrality].fEta->Fill(grandMa->Eta());
a8ef1999 1196 }
1197 if(eabsY<0.5){
1198 // fill electron kinematics
11ff28c5 1199 fHist[iq][kElectron2nd][2][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1200 fHist[iq][kElectron2nd][2][fCentrality].fPt->Fill(mcpart->Pt());
1201 fHist[iq][kElectron2nd][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1202 fHist[iq][kElectron2nd][2][fCentrality].fEta->Fill(mcpart->Eta());
a8ef1999 1203
1204 // fill mother hadron kinematics
11ff28c5 1205 fHist[iq][kDeHadron][2][fCentrality].fPdgCode->Fill(grandMaPDG);
1206 fHist[iq][kDeHadron][2][fCentrality].fPt->Fill(grandMa->Pt());
1207 fHist[iq][kDeHadron][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(grandMa));
1208 fHist[iq][kDeHadron][2][fCentrality].fEta->Fill(grandMa->Eta());
a8ef1999 1209 }
0792aa82 1210
1211 return;
1212 }
1213 }
1214
1215 partMother = grandMa;
1216 } // end of iteration
1217 } // end of if
1218 if ((isMotherDirectCharm == kTRUE && kquark == kCharm) || kquark == kBeauty) {
1219 for (Int_t i=0; i<fNparents; i++){
1220 if (abs(maPdgcodeCopy)==fParentSelect[iq][i]){
1221
1222 // fill electron kinematics
11ff28c5 1223 fHist[iq][kElectron][0][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1224 fHist[iq][kElectron][0][fCentrality].fPt->Fill(mcpart->Pt());
1225 fHist[iq][kElectron][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1226 fHist[iq][kElectron][0][fCentrality].fEta->Fill(mcpart->Eta());
0792aa82 1227
1228 // fill mother hadron kinematics
11ff28c5 1229 fHist[iq][keHadron][0][fCentrality].fPdgCode->Fill(maPdgcodeCopy);
1230 fHist[iq][keHadron][0][fCentrality].fPt->Fill(partMotherCopy->Pt());
1231 fHist[iq][keHadron][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
1232 fHist[iq][keHadron][0][fCentrality].fEta->Fill(partMotherCopy->Eta());
a8ef1999 1233
1234 if(eabsEta<0.9){
1235 // fill electron kinematics
11ff28c5 1236 fHist[iq][kElectron][1][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1237 fHist[iq][kElectron][1][fCentrality].fPt->Fill(mcpart->Pt());
1238 fHist[iq][kElectron][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1239 fHist[iq][kElectron][1][fCentrality].fEta->Fill(mcpart->Eta());
a8ef1999 1240
1241 // fill mother hadron kinematics
11ff28c5 1242 fHist[iq][keHadron][1][fCentrality].fPdgCode->Fill(maPdgcodeCopy);
1243 fHist[iq][keHadron][1][fCentrality].fPt->Fill(partMotherCopy->Pt());
1244 fHist[iq][keHadron][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
1245 fHist[iq][keHadron][1][fCentrality].fEta->Fill(partMotherCopy->Eta());
a8ef1999 1246 }
1247 if(eabsY<0.5){
1248 // fill electron kinematics
11ff28c5 1249 fHist[iq][kElectron][2][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1250 fHist[iq][kElectron][2][fCentrality].fPt->Fill(mcpart->Pt());
1251 fHist[iq][kElectron][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1252 fHist[iq][kElectron][2][fCentrality].fEta->Fill(mcpart->Eta());
a8ef1999 1253
1254 // fill mother hadron kinematics
11ff28c5 1255 fHist[iq][keHadron][2][fCentrality].fPdgCode->Fill(maPdgcodeCopy);
1256 fHist[iq][keHadron][2][fCentrality].fPt->Fill(partMotherCopy->Pt());
1257 fHist[iq][keHadron][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
1258 fHist[iq][keHadron][2][fCentrality].fEta->Fill(partMotherCopy->Eta());
a8ef1999 1259 }
0792aa82 1260
1261 }
1262 }
1263 } // end of if
1264
1265}
259c3296 1266
1267//__________________________________________
75d81601 1268void AliHFEmcQA::IdentifyMother(Int_t motherlabel, Int_t &motherpdg, Int_t &grandmotherlabel)
259c3296 1269{
dbe3abbe 1270 // find mother pdg code and label
1271
75d81601 1272 if (motherlabel < 0) {
259c3296 1273 AliDebug(1, "Stack label is negative, return\n");
1274 return;
1275 }
faee3b18 1276 AliMCParticle *mctrack = NULL;
1277 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(motherlabel))))) return;
1278 TParticle *heavysMother = mctrack->Particle();
75d81601 1279 motherpdg = heavysMother->GetPdgCode();
1280 grandmotherlabel = heavysMother->GetFirstMother();
1281 AliDebug(1,Form("ancestor pdg code= %d\n",motherpdg));
259c3296 1282}
1283
1284//__________________________________________
259c3296 1285void AliHFEmcQA::HardScattering(const Int_t kquark, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
1286{
dbe3abbe 1287 // mothertype -1 means this heavy quark coming from hard vertex
1288
faee3b18 1289 AliMCParticle *mctrack1 = NULL;
1290 AliMCParticle *mctrack2 = NULL;
1291 if(!(mctrack1 = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(4))))) return;
1292 if(!(mctrack2 = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(5))))) return;
1293 TParticle *afterinitialrad1 = mctrack1->Particle();
1294 TParticle *afterinitialrad2 = mctrack2->Particle();
259c3296 1295
1296 motherlabel = -1;
1297
1298 if (abs(afterinitialrad1->GetPdgCode()) == fgkGluon && abs(afterinitialrad2->GetPdgCode()) == fgkGluon){
1299 AliDebug(1,"heavy from gluon gluon pair creation!\n");
1300 mothertype = -1;
1301 motherID = fgkGluon;
1302 }
1303 else if (abs(afterinitialrad1->GetPdgCode()) == kquark || abs(afterinitialrad2->GetPdgCode()) == kquark){ // one from Q and the other from g
1304 AliDebug(1,"heavy from flavor exitation!\n");
1305 mothertype = -1;
1306 motherID = kquark;
1307 }
1308 else if (abs(afterinitialrad1->GetPdgCode()) == abs(afterinitialrad2->GetPdgCode())){
1309 AliDebug(1,"heavy from q-qbar pair creation!\n");
1310 mothertype = -1;
1311 motherID = 1;
1312 }
1313 else {
1314 AliDebug(1,"something strange!\n");
1315 mothertype = -999;
1316 motherlabel = -999;
1317 motherID = -999;
1318 }
1319}
1320
1321//__________________________________________
259c3296 1322Bool_t AliHFEmcQA::IsFromInitialShower(Int_t inputmotherlabel, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
1323{
dbe3abbe 1324 // mothertype -2 means this heavy quark coming from initial state
1325
faee3b18 1326 AliMCParticle *mctrack = NULL;
259c3296 1327 if (inputmotherlabel==2 || inputmotherlabel==3){ // mother exist before initial state radiation
faee3b18 1328 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(inputmotherlabel))))) return kFALSE;
1329 TParticle *heavysMother = mctrack->Particle();
259c3296 1330 motherID = heavysMother->GetPdgCode();
1331 mothertype = -2; // there is mother before initial state radiation
1332 motherlabel = inputmotherlabel;
1333 AliDebug(1,"initial parton shower! \n");
1334
1335 return kTRUE;
1336 }
1337
1338 return kFALSE;
1339}
1340
1341//__________________________________________
259c3296 1342Bool_t AliHFEmcQA::IsFromFinalParton(Int_t inputmotherlabel, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
1343{
dbe3abbe 1344 // mothertype 2 means this heavy quark coming from final state
1345
faee3b18 1346 AliMCParticle *mctrack = NULL;
259c3296 1347 if (inputmotherlabel > 5){ // mother exist after hard scattering
faee3b18 1348 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(inputmotherlabel))))) return kFALSE;
1349 TParticle *heavysMother = mctrack->Particle();
259c3296 1350 motherID = heavysMother->GetPdgCode();
1351 mothertype = 2; //
1352 motherlabel = inputmotherlabel;
1353 AliDebug(1,Form("heavy quark from %d after hard scattering! \n",motherID));
1354
1355 return kTRUE;
1356 }
1357 return kFALSE;
1358}
1359
1360//__________________________________________
dbe3abbe 1361void AliHFEmcQA::ReportStrangeness(Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
259c3296 1362{
dbe3abbe 1363 // mark strange behavior
1364
259c3296 1365 mothertype = -888;
1366 motherlabel = -888;
1367 motherID = -888;
1368 AliDebug(1,"something strange!\n");
1369}
1370
0792aa82 1371//__________________________________________
11ff28c5 1372Int_t AliHFEmcQA::GetSource(const AliVParticle* const mcpart)
0792aa82 1373{
9bcfd1ab 1374 // decay particle's origin
0792aa82 1375
9bcfd1ab 1376 //if ( abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return -1;
0792aa82 1377
1378 Int_t origin = -1;
50685501 1379 Bool_t isFinalOpenCharm = kFALSE;
0792aa82 1380
faee3b18 1381 if(!mcpart){
1382 AliDebug(1, "Stack label is negative or no mcparticle, return\n");
1383 return -1;
1384 }
1385
0792aa82 1386 // mother
11ff28c5 1387 // Information not in the base class, cast necessary
1388 Int_t iLabel = GetMother(mcpart);
0792aa82 1389 if (iLabel<0){
1390 AliDebug(1, "Stack label is negative, return\n");
1391 return -1;
1392 }
1393
11ff28c5 1394 const AliVParticle *partMother = fMCEvent->GetTrack(iLabel);
1395 Int_t maPdgcode = partMother->PdgCode();
0792aa82 1396
1397 // if the mother is charmed hadron
1398 if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
1399
1400 for (Int_t i=0; i<fNparents; i++){
1401 if (abs(maPdgcode)==fParentSelect[0][i]){
50685501 1402 isFinalOpenCharm = kTRUE;
0792aa82 1403 }
1404 }
50685501 1405 if (!isFinalOpenCharm) return -1;
0792aa82 1406
1407 // iterate until you find B hadron as a mother or become top ancester
1408 for (Int_t i=1; i<fgkMaxIter; i++){
11ff28c5 1409
1410 Int_t jLabel = GetMother(partMother);
0792aa82 1411 if (jLabel == -1){
1412 origin = kDirectCharm;
1413 return origin;
1414 }
1415 if (jLabel < 0){ // safety protection
1416 AliDebug(1, "Stack label is negative, return\n");
1417 return -1;
1418 }
1419
1420 // if there is an ancester
11ff28c5 1421 const AliVParticle* grandMa = fMCEvent->GetTrack(jLabel);
1422 Int_t grandMaPDG = grandMa->PdgCode();
0792aa82 1423
70b5ea26 1424 for (Int_t j=0; j<fNparents; j++){
1425 if (abs(grandMaPDG)==fParentSelect[1][j]){
0792aa82 1426 origin = kBeautyCharm;
1427 return origin;
1428 }
1429 }
1430
1431 partMother = grandMa;
1432 } // end of iteration
1433 } // end of if
1434 else if ( int(abs(maPdgcode)/100.) == kBeauty || int(abs(maPdgcode)/1000.) == kBeauty ) {
1435 for (Int_t i=0; i<fNparents; i++){
1436 if (abs(maPdgcode)==fParentSelect[1][i]){
1437 origin = kDirectBeauty;
1438 return origin;
1439 }
1440 }
1441 } // end of if
1442 else if ( abs(maPdgcode) == 22 ) {
1443 origin = kGamma;
1444 return origin;
1445 } // end of if
1446 else if ( abs(maPdgcode) == 111 ) {
1447 origin = kPi0;
1448 return origin;
1449 } // end of if
1450
1451 return origin;
1452}
1453
1454//__________________________________________
c2690925 1455Int_t AliHFEmcQA::GetSource(const TParticle * const mcpart)
0792aa82 1456{
9bcfd1ab 1457 // decay particle's origin
0792aa82 1458
9bcfd1ab 1459 //if ( abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return -1;
0792aa82 1460
1461 Int_t origin = -1;
50685501 1462 Bool_t isFinalOpenCharm = kFALSE;
0792aa82 1463
faee3b18 1464 if(!mcpart){
1465 AliDebug(1, "no mcparticle, return\n");
1466 return -1;
1467 }
1468
0792aa82 1469 Int_t iLabel = mcpart->GetFirstMother();
1470 if (iLabel<0){
1471 AliDebug(1, "Stack label is negative, return\n");
1472 return -1;
1473 }
1474
faee3b18 1475 AliMCParticle *mctrack = NULL;
1476 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return -1;
1477 TParticle *partMother = mctrack->Particle();
0792aa82 1478 Int_t maPdgcode = partMother->GetPdgCode();
1479
1480 // if the mother is charmed hadron
1481 if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
1482
1483 for (Int_t i=0; i<fNparents; i++){
1484 if (abs(maPdgcode)==fParentSelect[0][i]){
50685501 1485 isFinalOpenCharm = kTRUE;
0792aa82 1486 }
1487 }
50685501 1488 if (!isFinalOpenCharm) return -1;
0792aa82 1489
1490 // iterate until you find B hadron as a mother or become top ancester
1491 for (Int_t i=1; i<fgkMaxIter; i++){
1492
1493 Int_t jLabel = partMother->GetFirstMother();
1494 if (jLabel == -1){
1495 origin = kDirectCharm;
1496 return origin;
1497 }
1498 if (jLabel < 0){ // safety protection
1499 AliDebug(1, "Stack label is negative, return\n");
1500 return -1;
1501 }
1502
1503 // if there is an ancester
faee3b18 1504 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return -1;
1505 TParticle* grandMa = mctrack->Particle();
0792aa82 1506 Int_t grandMaPDG = grandMa->GetPdgCode();
1507
70b5ea26 1508 for (Int_t j=0; j<fNparents; j++){
1509 if (abs(grandMaPDG)==fParentSelect[1][j]){
0792aa82 1510 origin = kBeautyCharm;
1511 return origin;
1512 }
1513 }
1514
1515 partMother = grandMa;
1516 } // end of iteration
1517 } // end of if
1518 else if ( int(abs(maPdgcode)/100.) == kBeauty || int(abs(maPdgcode)/1000.) == kBeauty ) {
1519 for (Int_t i=0; i<fNparents; i++){
1520 if (abs(maPdgcode)==fParentSelect[1][i]){
1521 origin = kDirectBeauty;
1522 return origin;
1523 }
1524 }
1525 } // end of if
1526 else if ( abs(maPdgcode) == 22 ) {
1527 origin = kGamma;
1528 return origin;
1529 } // end of if
1530 else if ( abs(maPdgcode) == 111 ) {
1531 origin = kPi0;
1532 return origin;
1533 } // end of if
1534 else origin = kElse;
1535
1536 return origin;
1537}
70da6c5a 1538
1539//__________________________________________
e3fc062d 1540Int_t AliHFEmcQA::GetElecSource(TParticle * const mcpart)
1541{
1542 // decay particle's origin
70da6c5a 1543
e3fc062d 1544 if(!mcpart){
1545 AliDebug(1, "no mcparticle, return\n");
1546 return -1;
1547 }
1548
bf892a6a 1549 if ( abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return kMisID;
1550
1551 Int_t origin = -1;
1552 Bool_t isFinalOpenCharm = kFALSE;
1553
e3fc062d 1554 Int_t iLabel = mcpart->GetFirstMother();
1555 if (iLabel<0){
1556 AliDebug(1, "Stack label is negative, return\n");
1557 return -1;
1558 }
1559
1560 AliMCParticle *mctrack = NULL;
c2690925 1561 Int_t tmpMomLabel=0;
e3fc062d 1562 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return -1;
1563 TParticle *partMother = mctrack->Particle();
c2690925 1564 TParticle *partMotherCopy = mctrack->Particle();
e3fc062d 1565 Int_t maPdgcode = partMother->GetPdgCode();
1566
1567 // if the mother is charmed hadron
11ff28c5 1568
1569 if(abs(maPdgcode)==443){ // J/spi
1570 Int_t jLabel = partMother->GetFirstMother();
1571 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return kJpsi;
1572 TParticle* grandMa = mctrack->Particle();
1573 Int_t grandMaPDG = grandMa->GetPdgCode();
1574 if((int(abs(grandMaPDG)/100.)%10) == kBeauty || (int(abs(grandMaPDG)/1000.)%10) == kBeauty) return kB2Jpsi;
1575 else return -1;
1576 }
1577 else if ( (int(abs(maPdgcode)/100.)%10) == kCharm || (int(abs(maPdgcode)/1000.)%10) == kCharm ) {
e3fc062d 1578
1579 for (Int_t i=0; i<fNparents; i++){
1580 if (abs(maPdgcode)==fParentSelect[0][i]){
1581 isFinalOpenCharm = kTRUE;
1582 }
1583 }
1584 if (!isFinalOpenCharm) return -1;
1585
1586 // iterate until you find B hadron as a mother or become top ancester
1587 for (Int_t i=1; i<fgkMaxIter; i++){
1588
1589 Int_t jLabel = partMother->GetFirstMother();
1590 if (jLabel == -1){
11ff28c5 1591 return kDirectCharm;
e3fc062d 1592 }
1593 if (jLabel < 0){ // safety protection
1594 AliDebug(1, "Stack label is negative, return\n");
1595 return -1;
1596 }
1597
1598 // if there is an ancester
1599 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return -1;
1600 TParticle* grandMa = mctrack->Particle();
1601 Int_t grandMaPDG = grandMa->GetPdgCode();
1602
1603 for (Int_t j=0; j<fNparents; j++){
1604 if (abs(grandMaPDG)==fParentSelect[1][j]){
11ff28c5 1605 return kBeautyCharm;
e3fc062d 1606 }
1607 }
1608
1609 partMother = grandMa;
1610 } // end of iteration
1611 } // end of if
e17c1f86 1612 else if ( (int(abs(maPdgcode)/100.)%10) == kBeauty || (int(abs(maPdgcode)/1000.)%10) == kBeauty ) {
e3fc062d 1613 for (Int_t i=0; i<fNparents; i++){
1614 if (abs(maPdgcode)==fParentSelect[1][i]){
11ff28c5 1615 return kDirectBeauty;
e3fc062d 1616 }
1617 }
1618 } // end of if
c2690925 1619 else if ( abs(maPdgcode) == 22 ) { //conversion
1620
1621 tmpMomLabel = partMotherCopy->GetFirstMother();
1622 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tmpMomLabel))))) return -1;
1623 partMother = mctrack->Particle();
1624 maPdgcode = partMother->GetPdgCode();
1625 if ( abs(maPdgcode) == 111 ) {
11ff28c5 1626 return kGammaPi0;
c2690925 1627 }
1628 else if ( abs(maPdgcode) == 221 ) {
11ff28c5 1629 return kGammaEta;
c2690925 1630 }
1631 else if ( abs(maPdgcode) == 223 ) {
11ff28c5 1632 return kGammaOmega;
c2690925 1633 }
1634 else if ( abs(maPdgcode) == 333 ) {
11ff28c5 1635 return kGammaPhi;
c2690925 1636 }
1637 else if ( abs(maPdgcode) == 331 ) {
11ff28c5 1638 return kGammaEtaPrime;
c2690925 1639 }
1640 else if ( abs(maPdgcode) == 113 ) {
11ff28c5 1641 return kGammaRho0;
c2690925 1642 }
1643 else origin = kElse;
e3fc062d 1644 return origin;
c2690925 1645
11ff28c5 1646 }
e3fc062d 1647 else if ( abs(maPdgcode) == 111 ) {
11ff28c5 1648 return kPi0;
1649 }
e3fc062d 1650 else if ( abs(maPdgcode) == 221 ) {
11ff28c5 1651 return kEta;
1652 }
e3fc062d 1653 else if ( abs(maPdgcode) == 223 ) {
11ff28c5 1654 return kOmega;
1655 }
e3fc062d 1656 else if ( abs(maPdgcode) == 333 ) {
11ff28c5 1657 return kPhi;
1658 }
e3fc062d 1659 else if ( abs(maPdgcode) == 331 ) {
11ff28c5 1660 return kEtaPrime;
1661 }
e3fc062d 1662 else if ( abs(maPdgcode) == 113 ) {
11ff28c5 1663 return kRho0;
1664 }
1665 else if ( abs(maPdgcode) == 321 || abs(maPdgcode) == 130 ) {
1666 return kKe3;
1667 }
e17c1f86 1668 else{
1669 origin = kElse;
1670 }
e3fc062d 1671 return origin;
70da6c5a 1672}
9250ffbf 1673//__________________________________________
8c1c76e9 1674Double_t AliHFEmcQA::GetWeightFactor(AliMCParticle *mctrack, const Int_t iBgLevel){
9250ffbf 1675 //
8c1c76e9 1676 // Get weighting factor for the realistic background estimation, for three possible background yield levels, indicated by the argument "iLevel": the best estimate (0), the lower uncertainty level (1), and the upper uncertainty level (2)
9250ffbf 1677 //
1678 AliMCParticle *mctrackmother = NULL;
1679 Double_t weightElecBg = 0.;
1680 Double_t mesonPt = 0.;
1681 Double_t bgcategory = 0.;
1682 Int_t mArr = -1;
1683 Int_t mesonID = GetElecSource(mctrack->Particle());
1684 if(mesonID==kGammaPi0 || mesonID==kPi0) mArr=0; //pion
1685 else if(mesonID==kGammaEta || mesonID==kEta) mArr=1; //eta
1686 else if(mesonID==kGammaOmega || mesonID==kOmega) mArr=2; //omega
1687 else if(mesonID==kGammaPhi || mesonID==kPhi) mArr=3; //phi
1688 else if(mesonID==kGammaEtaPrime || mesonID==kEtaPrime) mArr=4; //etaprime
1689 else if(mesonID==kGammaRho0 || mesonID==kRho0) mArr=5; //rho
1690
11ff28c5 1691 Double_t datamc[25]={-999,-999,-999,-999,-999,-999,-999,-999,-999,-999,-999,-999,-999,-999,-999,-999, -999, -999, -999, -999, -999, -999, -999, -999, -999};
a8ef1999 1692 Double_t xr[3]={-999,-999,-999};
1693 datamc[0] = mesonID;
1694 datamc[17] = mctrack->Pt(); //electron pt
1695 datamc[18] = mctrack->Eta(); //electron eta
1696
1697 mctrack->XvYvZv(xr);
1698 datamc[9] = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
1699 datamc[10] = xr[2];
1700
1701 TParticle *mcpart = mctrack->Particle();
1702 if(mcpart){
1703 datamc[14] = mcpart->GetUniqueID();
1704 }
1705
9250ffbf 1706 if(!(mArr<0)){
1707 if(mesonID>=kGammaPi0) { // conversion electron, be careful with the enum odering
1708 Int_t glabel=TMath::Abs(mctrack->GetMother()); // gamma label
1709 if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
1710 glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's label
1711 if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
1712 mesonPt = mctrackmother->Pt(); //meson pt
1713 bgcategory = 1.;
a8ef1999 1714 datamc[1] = bgcategory;
1715 datamc[2] = mesonPt;
1716 mctrackmother->XvYvZv(xr);
1717 datamc[11] = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
1718 datamc[12] = xr[2];
1719
1720 mcpart = mctrackmother->Particle();
1721 if(mcpart){
1722 datamc[15] = mcpart->GetUniqueID();
1723 }
1724 if(glabel>fMCEvent->GetNumberOfPrimaries()) {
1725 bgcategory = 2.;
1726 datamc[1] = bgcategory;
1727 //printf("I should be gamma meson = %d mesonlabel= %d NumberOfPrimaries= %d \n",mctrackmother->PdgCode(),glabel,fMCEvent->GetNumberOfPrimaries());
1728 glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother
1729 if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
1730 datamc[3]=mctrackmother->PdgCode();
1731 datamc[4]=mctrackmother->Pt();
1732 if(TMath::Abs(mctrackmother->PdgCode())==310){
1733 bgcategory = 3.;
1734 glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother's mother
1735 if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
1736 datamc[5]=mctrackmother->PdgCode();
1737 glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother
1738 if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
1739 datamc[6]=mctrackmother->PdgCode();
1740 }
1741 }
1742 }
1743 }
1744 }
9250ffbf 1745 }
1746 }
1747 }
1748 else{ // nonHFE except for the conversion electron
1749 Int_t glabel=TMath::Abs(mctrack->GetMother());
1750 if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
1751 mesonPt=mctrackmother->Pt(); //meson pt
1752 bgcategory = -1.;
a8ef1999 1753 datamc[1] = bgcategory;
1754 datamc[2] = mesonPt;
1755 mctrackmother->XvYvZv(xr);
1756 datamc[11] = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
1757 datamc[12] = xr[2];
1758
1759 mcpart = mctrackmother->Particle();
1760 if(mcpart){
1761 datamc[15] = mcpart->GetUniqueID();
1762 }
1763 if(glabel>fMCEvent->GetNumberOfPrimaries()) {
1764 bgcategory = -2.;
1765 datamc[1] = bgcategory;
1766 glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother
1767 if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
1768 datamc[3]=mctrackmother->PdgCode();
1769 datamc[4]=mctrackmother->Pt();
1770 if(TMath::Abs(mctrackmother->PdgCode())==310){
1771 bgcategory = -3.;
1772 glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother's mother
1773 if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
1774 datamc[5]=mctrackmother->PdgCode();
1775 glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother
1776 if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
1777 datamc[6]=mctrackmother->PdgCode();
1778 }
1779 }
1780 }
1781 }
1782 }
9250ffbf 1783 }
1784 }
8c1c76e9 1785
a8ef1999 1786
8c1c76e9 1787 if(fIsPbPb){
11ff28c5 1788 Int_t centBin = GetWeightCentralityBin(fPerCentrality);
1789 if(centBin < 0)return 0.;
1790 weightElecBg=fElecBackgroundFactor[iBgLevel][centBin][mArr][kBgPtBins-1];
8c1c76e9 1791 for(int ii=0; ii<kBgPtBins; ii++){
1792 if((mesonPt > fBinLimit[ii]) && (mesonPt < fBinLimit[ii+1])){
11ff28c5 1793 weightElecBg = fElecBackgroundFactor[iBgLevel][centBin][mArr][ii];
8c1c76e9 1794 break;
1795 }
1796 }
9250ffbf 1797 }
8c1c76e9 1798 else{
1799 weightElecBg=fElecBackgroundFactor[iBgLevel][0][mArr][kBgPtBins-1];
9250ffbf 1800 for(int ii=0; ii<kBgPtBins; ii++){
8c1c76e9 1801 if((mesonPt > fBinLimit[ii]) && (mesonPt < fBinLimit[ii+1])){
1802 weightElecBg = fElecBackgroundFactor[iBgLevel][0][mArr][ii];
1803 break;
1804 }
9250ffbf 1805 }
8c1c76e9 1806 }
9250ffbf 1807 }
a8ef1999 1808
1809 datamc[13] = weightElecBg;
1810 datamc[16] = Double_t(fContainerStep);
1811
1812 datamc[7] = fHfeImpactR;
1813 datamc[8] = fHfeImpactnsigmaR;
1814
1815 datamc[19] = fRecPt;
1816 datamc[20] = fRecEta;
1817 datamc[21] = fRecPhi;
1818 datamc[22] = fLyrhit;
1819 datamc[23] = fLyrstat;
11ff28c5 1820 datamc[24] = fCentrality;
a8ef1999 1821
1822 if(fIsDebugStreamerON && fTreeStream){
1823 if(!iBgLevel){
1824 (*fTreeStream)<<"nonhfeQA"<<
1825 "mesonID="<<datamc[0]<<
1826 "bgcategory="<<datamc[1]<<
1827 "mesonPt="<<datamc[2]<<
1828 "mesonMomPdg="<<datamc[3]<<
1829 "mesonMomPt="<<datamc[4]<<
1830 "mesonGMomPdg="<<datamc[5]<<
1831 "mesonGGMomPdg="<<datamc[6]<<
1832 "eIPAbs="<<datamc[7]<<
1833 "eIPSig="<<datamc[8]<<
1834 "eR="<<datamc[9]<<
1835 "eZ="<<datamc[10]<<
1836 "mesonR="<<datamc[11]<<
1837 "mesonZ="<<datamc[12]<<
1838 "weightElecBg="<<datamc[13]<<
1839 "eUniqID="<<datamc[14]<<
1840 "mesonUniqID="<<datamc[15]<<
1841 "containerStep="<<datamc[16]<<
1842 "emcpt="<<datamc[17]<<
1843 "emceta="<<datamc[18]<<
1844 "erecpt="<<datamc[19]<<
1845 "ereceta="<<datamc[20]<<
1846 "erecphi="<<datamc[21]<<
1847 "itshit="<<datamc[22]<<
11ff28c5 1848 "itsstat="<<datamc[23]<<
1849 "centrality="<<datamc[24]
a8ef1999 1850 << "\n";
1851 }
1852 }
1853
1854 Double_t returnval = bgcategory*weightElecBg;
1855 if(TMath::Abs(bgcategory)>1) returnval = bgcategory/2.;
1856
1857 return returnval;
9250ffbf 1858}
1859
11ff28c5 1860//__________________________________________
1861Int_t AliHFEmcQA::GetMother(const AliVParticle * const mcpart){
1862 //
1863 // Wrapper to get the mother label
1864 //
1865 Int_t label = -1;
1866 TClass *type = mcpart->IsA();
1867 if(type == AliMCParticle::Class()){
1868 // ESD analysis
1869 const AliMCParticle *emcpart = static_cast<const AliMCParticle *>(mcpart);
1870 label = emcpart->GetMother();
1871 } else if(type == AliAODMCParticle::Class()){
1872 // AOD analysis
1873 const AliAODMCParticle *amcpart = static_cast<const AliAODMCParticle *>(mcpart);
1874 label = amcpart->GetMother();
1875 }
1876 return label;
1877}
1878//__________________________________________
1879Int_t AliHFEmcQA::GetWeightCentralityBin(const Float_t percentile) const {
1880 //
1881 //translate the centrality percentile into the centrality bin of the reference weighting histograms for electron background
1882 //
1883
1884 Float_t centralityLimits[12]= {0.,5.,10., 20., 30., 40., 50., 60.,70.,80., 90., 100.};
1885 Int_t bin = -1;
1886 for(Int_t ibin = 0; ibin < 11; ibin++){
1887 if(percentile >= centralityLimits[ibin] && percentile < centralityLimits[ibin+1]){
1888 bin = ibin;
1889 break;
1890 }
1891 }
1892 return bin;
1893}
70da6c5a 1894//__________________________________________
1895void AliHFEmcQA::AliHists::FillList(TList *l) const {
1896 //
1897 // Fill Histos into a list for output
1898 //
1899 if(fPdgCode) l->Add(fPdgCode);
1900 if(fPt) l->Add(fPt);
1901 if(fY) l->Add(fY);
1902 if(fEta) l->Add(fEta);
1903}
1904
1905//__________________________________________
1906void AliHFEmcQA::AliHistsComm::FillList(TList *l) const {
1907 //
1908 // Fill Histos into a list for output
1909 //
1910 if(fNq) l->Add(fNq);
1911 if(fProcessID) l->Add(fProcessID);
1912 if(fePtRatio) l->Add(fePtRatio);
ccc37cdc 1913 if(fPtCorr) l->Add(fPtCorr);
c2690925 1914 if(fPtCorrDp) l->Add(fPtCorrDp);
1915 if(fPtCorrD0) l->Add(fPtCorrD0);
1916 if(fPtCorrDrest) l->Add(fPtCorrDrest);
70da6c5a 1917 if(fDePtRatio) l->Add(fDePtRatio);
1918 if(feDistance) l->Add(feDistance);
1919 if(fDeDistance) l->Add(fDeDistance);
1920}
1921