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