]> 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
Transition PWG3 --> PWGHF
[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>
35
36#include <AliLog.h>
faee3b18 37#include <AliMCEvent.h>
9bcfd1ab 38#include <AliGenEventHeader.h>
0792aa82 39#include <AliAODMCParticle.h>
c2690925 40#include <AliStack.h>
259c3296 41
42#include "AliHFEmcQA.h"
70da6c5a 43#include "AliHFEtools.h"
c2690925 44#include "AliHFEcollection.h"
259c3296 45
259c3296 46ClassImp(AliHFEmcQA)
47
48//_______________________________________________________________________________________________
8c1c76e9 49 AliHFEmcQA::AliHFEmcQA() :
50 fMCEvent(NULL)
70da6c5a 51 ,fMCHeader(NULL)
52 ,fMCArray(NULL)
53 ,fQAhistos(NULL)
c2690925 54 ,fMCQACollection(NULL)
8c1c76e9 55 ,fNparents(0)
56 ,fCentrality(0)
57 ,fIsPbPb(kFALSE)
58 ,fIsppMultiBin(kFALSE)
259c3296 59{
60 // Default constructor
e088e9f0 61 for(Int_t mom = 0; mom < 9; mom++){
62 fhD[mom] = NULL;
63 fhDLogbin[mom] = NULL;
64 }
bf892a6a 65 for(Int_t mom = 0; mom < 50; mom++){
66 fHeavyQuark[mom] = NULL;
67 }
68 for(Int_t mom = 0; mom < 2; mom++){
69 fIsHeavy[mom] = 0;
70 }
8c1c76e9 71 memset(fElecBackgroundFactor, 0, sizeof(Double_t) * kElecBgSpecies * kBgPtBins * kCentBins * kBgLevels);
9250ffbf 72 memset(fBinLimit, 0, sizeof(Double_t) * (kBgPtBins+1));
259c3296 73}
74
75//_______________________________________________________________________________________________
76AliHFEmcQA::AliHFEmcQA(const AliHFEmcQA&p):
8c1c76e9 77TObject(p)
78 ,fMCEvent(NULL)
70da6c5a 79 ,fMCHeader(NULL)
80 ,fMCArray(NULL)
81 ,fQAhistos(p.fQAhistos)
c2690925 82 ,fMCQACollection(p.fMCQACollection)
8c1c76e9 83 ,fNparents(p.fNparents)
84 ,fCentrality(0)
85 ,fIsPbPb(kFALSE)
86 ,fIsppMultiBin(kFALSE)
259c3296 87{
88 // Copy constructor
e088e9f0 89 for(Int_t mom = 0; mom < 9; mom++){
90 fhD[mom] = NULL;
91 fhDLogbin[mom] = NULL;
92 }
bf892a6a 93 for(Int_t mom = 0; mom < 50; mom++){
94 fHeavyQuark[mom] = NULL;
95 }
96 for(Int_t mom = 0; mom < 2; mom++){
97 fIsHeavy[mom] = 0;
98 }
8c1c76e9 99 memset(fElecBackgroundFactor, 0, sizeof(Double_t) * kElecBgSpecies * kBgPtBins * kCentBins * kBgLevels);
9250ffbf 100 memset(fBinLimit, 0, sizeof(Double_t) * (kBgPtBins+1));
259c3296 101}
102
103//_______________________________________________________________________________________________
104AliHFEmcQA&
105AliHFEmcQA::operator=(const AliHFEmcQA &)
106{
107 // Assignment operator
108
109 AliInfo("Not yet implemented.");
110 return *this;
111}
112
113//_______________________________________________________________________________________________
114AliHFEmcQA::~AliHFEmcQA()
115{
116 // Destructor
117
118 AliInfo("Analysis Done.");
119}
120
121//_______________________________________________________________________________________________
50685501 122void AliHFEmcQA::PostAnalyze() const
259c3296 123{
e3fc062d 124 //
125 // Post analysis
126 //
259c3296 127}
128
9250ffbf 129//_______________________________________________________________________________________________
130void AliHFEmcQA::SetBackgroundWeightFactor(Double_t *elecBackgroundFactor, Double_t *binLimit)
131{
e156c3bb 132 //
133 // copy background weighting factors into data member
134 //
8c1c76e9 135
136 memcpy(fElecBackgroundFactor,elecBackgroundFactor,sizeof(Double_t) * kElecBgSpecies * kBgPtBins * kCentBins * kBgLevels);
137 memcpy(fBinLimit,binLimit,sizeof(Double_t) * (kBgPtBins+1));
9250ffbf 138}
139
e3fc062d 140//__________________________________________
141void AliHFEmcQA::CreatDefaultHistograms(TList * const qaList)
142{
143 //
144 // make default histograms
145 //
146
147 if(!qaList) return;
148
149 fQAhistos = qaList;
150 fQAhistos->SetName("MCqa");
151
152 CreateHistograms(AliHFEmcQA::kCharm,0,"mcqa_"); // create histograms for charm
153 CreateHistograms(AliHFEmcQA::kBeauty,0,"mcqa_"); // create histograms for beauty
154 CreateHistograms(AliHFEmcQA::kOthers,0,"mcqa_"); // create histograms for beauty
155 CreateHistograms(AliHFEmcQA::kCharm,1,"mcqa_barrel_"); // create histograms for charm
156 CreateHistograms(AliHFEmcQA::kBeauty,1,"mcqa_barrel_"); // create histograms for beauty
157 CreateHistograms(AliHFEmcQA::kOthers,1,"mcqa_barrel_"); // create histograms for beauty
158 CreateHistograms(AliHFEmcQA::kCharm,2,"mcqa_unitY_"); // create histograms for charm
159 CreateHistograms(AliHFEmcQA::kBeauty,2,"mcqa_unitY_"); // create histograms for beauty
160 CreateHistograms(AliHFEmcQA::kOthers,2,"mcqa_unitY_"); // create histograms for beauty
e3fc062d 161
e97c2edf 162// check D spectra
ccc37cdc 163 TString kDspecies[9];
e97c2edf 164 kDspecies[0]="411";
165 kDspecies[1]="421";
166 kDspecies[2]="431";
167 kDspecies[3]="4122";
168 kDspecies[4]="4132";
169 kDspecies[5]="4232";
170 kDspecies[6]="4332";
ccc37cdc 171 kDspecies[7]="413";
172 kDspecies[8]="423";
173
174 const Double_t kPtbound[2] = {0.1, 20.}; //bin taken for considering inclusive e analysis binning
c2690925 175 Int_t iBin[2];
ccc37cdc 176 iBin[0] = 44; // bins in pt
c2690925 177 iBin[1] = 23; // bins in pt
ccc37cdc 178 Double_t* binEdges[1];
179 binEdges[0] = AliHFEtools::MakeLogarithmicBinning(iBin[0], kPtbound[0], kPtbound[1]);
e97c2edf 180
181 // bin size is chosen to consider ALICE D measurement
9250ffbf 182 Double_t xbins[15]={0,0.5,1.0,2.0,3.0,4.0,5.0,6.0,8.0,12,16,20,30,40,50};
e97c2edf 183 Double_t ybins[10]={-7.5,-1.0,-0.9,-0.8,-0.5,0.5,0.8,0.9,1.0,7.5};
184 TString hname;
ccc37cdc 185 for (Int_t iDmeson=0; iDmeson<9; iDmeson++){
e97c2edf 186 hname = "Dmeson"+kDspecies[iDmeson];
9250ffbf 187 fhD[iDmeson] = new TH2F(hname,hname+";p_{T} (GeV/c)",14,xbins,9,ybins);
ccc37cdc 188 hname = "DmesonLogbin"+kDspecies[iDmeson];
189 fhDLogbin[iDmeson] = new TH1F(hname,hname+";p_{T} (GeV/c)",iBin[0],binEdges[0]);
e97c2edf 190 if(fQAhistos) fQAhistos->Add(fhD[iDmeson]);
ccc37cdc 191 if(fQAhistos) fQAhistos->Add(fhDLogbin[iDmeson]);
e97c2edf 192 }
193
c2690925 194 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
195
8c1c76e9 196 Int_t kNcent;
197 if(fIsPbPb) kNcent=11;
198 else
199 {
200 if(fIsppMultiBin) kNcent=8;
201 else kNcent = 1;
202 }
203
c2690925 204 fMCQACollection = new AliHFEcollection("TaskMCQA", "MC QA histos for meason pt spectra");
8c1c76e9 205
206 for(Int_t centbin=0; centbin<kNcent; centbin++)
207 {
208 fMCQACollection->CreateTH1Farray(Form("pionspectra_centrbin%i",centbin), "pion yields: MC p_{t} ", iBin[1],kPtRange);
209 fMCQACollection->CreateTH1Farray(Form("etaspectra_centrbin%i",centbin), "eta yields: MC p_{t} ", iBin[1],kPtRange);
210 fMCQACollection->CreateTH1Farray(Form("omegaspectra_centrbin%i",centbin), "omega yields: MC p_{t} ", iBin[1],kPtRange);
211 fMCQACollection->CreateTH1Farray(Form("phispectra_centrbin%i",centbin), "phi yields: MC p_{t} ", iBin[1],kPtRange);
212 fMCQACollection->CreateTH1Farray(Form("etapspectra_centrbin%i",centbin), "etap yields: MC p_{t} ", iBin[1],kPtRange);
213 fMCQACollection->CreateTH1Farray(Form("rhospectra_centrbin%i",centbin), "rho yields: MC p_{t} ", iBin[1],kPtRange);
214
215 fMCQACollection->CreateTH1F(Form("pionspectraLog_centrbin%i",centbin), "pion yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
216 fMCQACollection->CreateTH1F(Form("etaspectraLog_centrbin%i",centbin), "eta yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
217 fMCQACollection->CreateTH1F(Form("omegaspectraLog_centrbin%i",centbin), "omega yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
218 fMCQACollection->CreateTH1F(Form("phispectraLog_centrbin%i",centbin), "phi yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
219 fMCQACollection->CreateTH1F(Form("etapspectraLog_centrbin%i",centbin), "etap yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
220 fMCQACollection->CreateTH1F(Form("rhospectraLog_centrbin%i",centbin), "rho yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
221 fMCQACollection->CreateTH1F(Form("kaonspectraLog_centrbin%i",centbin), "kaon yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
222 fMCQACollection->CreateTH1F(Form("k0LspectraLog_centrbin%i",centbin), "k0L yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
223
224 fMCQACollection->CreateTH1F(Form("piondaughters_centrbin%i",centbin), "pion yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
225 fMCQACollection->CreateTH1F(Form("etadaughters_centrbin%i",centbin), "eta yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
226 fMCQACollection->CreateTH1F(Form("omegadaughters_centrbin%i",centbin), "omega yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
227 fMCQACollection->CreateTH1F(Form("phidaughters_centrbin%i",centbin), "phi yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
228 fMCQACollection->CreateTH1F(Form("etapdaughters_centrbin%i",centbin), "etap yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
229 fMCQACollection->CreateTH1F(Form("rhodaughters_centrbin%i",centbin), "rho yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
230 }
c2690925 231
232 fQAhistos->Add(fMCQACollection->GetList());
233
e3fc062d 234}
235
259c3296 236//__________________________________________
dbe3abbe 237void AliHFEmcQA::CreateHistograms(const Int_t kquark, Int_t icut, TString hnopt)
259c3296 238{
239 // create histograms
240
faee3b18 241 if (!(kquark == kCharm || kquark == kBeauty || kquark == kOthers)) {
259c3296 242 AliDebug(1, "This task is only for heavy quark QA, return\n");
243 return;
244 }
dbe3abbe 245 Int_t iq = kquark - kCharm;
259c3296 246
247 TString kqTypeLabel[fgkqType];
dbe3abbe 248 if (kquark == kCharm){
249 kqTypeLabel[kQuark]="c";
250 kqTypeLabel[kantiQuark]="cbar";
251 kqTypeLabel[kHadron]="cHadron";
252 kqTypeLabel[keHadron]="ceHadron";
253 kqTypeLabel[kDeHadron]="nullHadron";
254 kqTypeLabel[kElectron]="ce";
255 kqTypeLabel[kElectron2nd]="nulle";
256 } else if (kquark == kBeauty){
257 kqTypeLabel[kQuark]="b";
258 kqTypeLabel[kantiQuark]="bbar";
259 kqTypeLabel[kHadron]="bHadron";
260 kqTypeLabel[keHadron]="beHadron";
261 kqTypeLabel[kDeHadron]="bDeHadron";
262 kqTypeLabel[kElectron]="be";
263 kqTypeLabel[kElectron2nd]="bce";
faee3b18 264 } else if (kquark == kOthers){
265 kqTypeLabel[kGamma-4]="gammae";
266 kqTypeLabel[kPi0-4]="pi0e";
267 kqTypeLabel[kElse-4]="elsee";
268 kqTypeLabel[kMisID-4]="miside";
259c3296 269 }
3a72645a 270
271 const Double_t kPtbound[2] = {0.1, 20.}; //bin taken for considering inclusive e analysis binning
faee3b18 272 Int_t iBin[1];
3a72645a 273 iBin[0] = 44; // bins in pt
faee3b18 274 Double_t* binEdges[1];
275 binEdges[0] = AliHFEtools::MakeLogarithmicBinning(iBin[0], kPtbound[0], kPtbound[1]);
c2690925 276
277
9250ffbf 278 Double_t xcorrbin[501];
279 for (int icorrbin = 0; icorrbin< 501; icorrbin++){
c2690925 280 xcorrbin[icorrbin]=icorrbin*0.1;
281 }
282
259c3296 283 TString hname;
faee3b18 284 if(kquark == kOthers){
285 for (Int_t iqType = 0; iqType < 4; iqType++ ){
286 hname = hnopt+"Pt_"+kqTypeLabel[iqType];
e3ae862b 287 //fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",iBin[0],binEdges[0]); //mj to compare with FONLL
288 fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",60,0.25,30.25);
faee3b18 289 hname = hnopt+"Y_"+kqTypeLabel[iqType];
290 fHist[iq][iqType][icut].fY = new TH1F(hname,hname,150,-7.5,7.5);
291 hname = hnopt+"Eta_"+kqTypeLabel[iqType];
292 fHist[iq][iqType][icut].fEta = new TH1F(hname,hname,150,-7.5,7.5);
293 // Fill List
294 if(fQAhistos) fHist[iq][iqType][icut].FillList(fQAhistos);
295 }
296 return;
297 }
259c3296 298 for (Int_t iqType = 0; iqType < fgkqType; iqType++ ){
dbe3abbe 299 if (iqType < keHadron && icut > 0) continue; // don't duplicate histogram for quark and hadron
259c3296 300 hname = hnopt+"PdgCode_"+kqTypeLabel[iqType];
dbe3abbe 301 fHist[iq][iqType][icut].fPdgCode = new TH1F(hname,hname,20001,-10000.5,10000.5);
259c3296 302 hname = hnopt+"Pt_"+kqTypeLabel[iqType];
c2690925 303 fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",iBin[0],binEdges[0]);
304 //fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",60,0.,30.); // mj to compare with FONLL
305 //fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",60,0.25,30.25); // mj to compare with FONLL
259c3296 306 hname = hnopt+"Y_"+kqTypeLabel[iqType];
dbe3abbe 307 fHist[iq][iqType][icut].fY = new TH1F(hname,hname,150,-7.5,7.5);
259c3296 308 hname = hnopt+"Eta_"+kqTypeLabel[iqType];
dbe3abbe 309 fHist[iq][iqType][icut].fEta = new TH1F(hname,hname,150,-7.5,7.5);
70da6c5a 310 // Fill List
311 if(fQAhistos) fHist[iq][iqType][icut].FillList(fQAhistos);
259c3296 312 }
313
dbe3abbe 314 if (icut == 0){
315 hname = hnopt+"Nq_"+kqTypeLabel[kQuark];
e3ae862b 316 fHistComm[iq][icut].fNq = new TH1F(hname,hname,50,-0.5,49.5);
dbe3abbe 317 hname = hnopt+"ProcessID_"+kqTypeLabel[kQuark];
318 fHistComm[iq][icut].fProcessID = new TH1F(hname,hname,21,-10.5,10.5);
319 }
320 hname = hnopt+"ePtRatio_"+kqTypeLabel[kQuark];
e088e9f0 321 fHistComm[iq][icut].fePtRatio = new TH2F(hname,hname+";p_{T} (GeV/c);momentum fraction",200,0,20,100,0,1);
ccc37cdc 322 hname = hnopt+"PtCorr_"+kqTypeLabel[kQuark];
9250ffbf 323 fHistComm[iq][icut].fPtCorr = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",500,xcorrbin,44,binEdges[0]);
c2690925 324 //fHistComm[iq][icut].fPtCorr = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",200,0,20,200,0,20);
325 hname = hnopt+"PtCorrDp_"+kqTypeLabel[kQuark];
9250ffbf 326 fHistComm[iq][icut].fPtCorrDp= new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",500,xcorrbin,44,binEdges[0]);
c2690925 327 //fHistComm[iq][icut].fPtCorrDp= new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",200,0,20,200,0,20);
328 hname = hnopt+"PtCorrD0_"+kqTypeLabel[kQuark];
9250ffbf 329 fHistComm[iq][icut].fPtCorrD0 = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",500,xcorrbin,44,binEdges[0]);
c2690925 330 //fHistComm[iq][icut].fPtCorrD0 = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",200,0,20,200,0,20);
331 hname = hnopt+"PtCorrDrest_"+kqTypeLabel[kQuark];
9250ffbf 332 fHistComm[iq][icut].fPtCorrDrest = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",500,xcorrbin,44,binEdges[0]);
c2690925 333 //fHistComm[iq][icut].fPtCorrDrest = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",200,0,20,200,0,20);
dbe3abbe 334 hname = hnopt+"DePtRatio_"+kqTypeLabel[kQuark];
e088e9f0 335 fHistComm[iq][icut].fDePtRatio = new TH2F(hname,hname+";p_{T} (GeV/c);momentum fraction",100,0,20,100,0,1);
dbe3abbe 336 hname = hnopt+"eDistance_"+kqTypeLabel[kQuark];
e088e9f0 337 fHistComm[iq][icut].feDistance= new TH2F(hname,hname+";p_{T} (GeV/c);distance (cm)",100,0,20,200,0,2);
dbe3abbe 338 hname = hnopt+"DeDistance_"+kqTypeLabel[kQuark];
e088e9f0 339 fHistComm[iq][icut].fDeDistance= new TH2F(hname,hname+";p_{T} (GeV/c);distance (cm)",100,0,20,200,0,2);
70da6c5a 340 if(fQAhistos) fHistComm[iq][icut].FillList(fQAhistos);
259c3296 341}
342
343//__________________________________________
344void AliHFEmcQA::Init()
345{
346 // called at begining every event
347
348 for (Int_t i=0; i<2; i++){
349 fIsHeavy[i] = 0;
350 }
351
352 fNparents = 7;
353
354 fParentSelect[0][0] = 411;
355 fParentSelect[0][1] = 421;
356 fParentSelect[0][2] = 431;
357 fParentSelect[0][3] = 4122;
358 fParentSelect[0][4] = 4132;
359 fParentSelect[0][5] = 4232;
360 fParentSelect[0][6] = 4332;
361
362 fParentSelect[1][0] = 511;
363 fParentSelect[1][1] = 521;
364 fParentSelect[1][2] = 531;
365 fParentSelect[1][3] = 5122;
366 fParentSelect[1][4] = 5132;
367 fParentSelect[1][5] = 5232;
368 fParentSelect[1][6] = 5332;
369
370}
371
9250ffbf 372//__________________________________________
c2690925 373void AliHFEmcQA::GetMesonKine()
374{
375 //
376 // get meson pt spectra
377 //
378
379 AliVParticle *mctrack2 = NULL;
380 AliMCParticle *mctrack0 = NULL;
381 AliVParticle *mctrackdaugt= NULL;
382 AliMCParticle *mctrackd= NULL;
383 Int_t id1=0, id2=0;
384
8c1c76e9 385
386 if(fCentrality>11) printf("warning centrality out of histogram array limits \n");
387
388
c2690925 389 for(Int_t imc = 0; imc <fMCEvent->GetNumberOfPrimaries(); imc++){
390 if(!(mctrack2 = fMCEvent->GetTrack(imc))) continue;
391 TParticle* mcpart0 = fMCEvent->Stack()->Particle(imc);
392 if(!mcpart0) continue;
393 mctrack0 = dynamic_cast<AliMCParticle *>(mctrack2);
f9f097c0 394 if(!mctrack0) continue;
8c1c76e9 395
396 if(!fIsPbPb&&!fIsppMultiBin) fCentrality=0;
397
c2690925 398 if(abs(mctrack0->PdgCode()) == 111) // pi0
399 {
400 if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
8c1c76e9 401 fMCQACollection->Fill(Form("pionspectra_centrbin%i",fCentrality),mctrack0->Pt());
402 fMCQACollection->Fill(Form("pionspectraLog_centrbin%i",fCentrality),mctrack0->Pt());
c2690925 403 }
404 id1=mctrack0->GetFirstDaughter();
405 id2=mctrack0->GetLastDaughter();
406 if(!((id2-id1)==2)) continue;
407 for(int idx=id1; idx<=id2; idx++){
408 if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
9250ffbf 409 if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
c2690925 410 if(abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
8c1c76e9 411 fMCQACollection->Fill(Form("piondaughters_centrbin%i",fCentrality),mctrackd->Pt());
c2690925 412 }
413 }
414 else if(abs(mctrack0->PdgCode()) == 221) // eta
415 {
416 if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
8c1c76e9 417 fMCQACollection->Fill(Form("etaspectra_centrbin%i",fCentrality),mctrack0->Pt());
418 fMCQACollection->Fill(Form("etaspectraLog_centrbin%i",fCentrality),mctrack0->Pt());
c2690925 419 }
420 id1=mctrack0->GetFirstDaughter();
421 id2=mctrack0->GetLastDaughter();
422 if(!((id2-id1)==2||(id2-id1)==3)) continue;
423 for(int idx=id1; idx<=id2; idx++){
424 if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
9250ffbf 425 if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
c2690925 426 if(abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
8c1c76e9 427 fMCQACollection->Fill(Form("etadaughters_centrbin%i",fCentrality),mctrackd->Pt());
c2690925 428 }
429 }
430 else if(abs(mctrack0->PdgCode()) == 223) // omega
431 {
432 if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
8c1c76e9 433 fMCQACollection->Fill(Form("omegaspectra_centrbin%i",fCentrality),mctrack0->Pt());
434 fMCQACollection->Fill(Form("omegaspectraLog_centrbin%i",fCentrality),mctrack0->Pt());
c2690925 435 }
436 id1=mctrack0->GetFirstDaughter();
437 id2=mctrack0->GetLastDaughter();
438 if(!((id2-id1)==1||(id2-id1)==2)) continue;
439 for(int idx=id1; idx<=id2; idx++){
440 if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
9250ffbf 441 if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
c2690925 442 if(abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
8c1c76e9 443 fMCQACollection->Fill(Form("omegadaughters_centrbin%i",fCentrality),mctrackd->Pt());
c2690925 444 }
445 }
446 else if(abs(mctrack0->PdgCode()) == 333) // phi
447 {
448 if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
8c1c76e9 449 fMCQACollection->Fill(Form("phispectra_centrbin%i",fCentrality),mctrack0->Pt());
450 fMCQACollection->Fill(Form("phispectraLog_centrbin%i",fCentrality),mctrack0->Pt());
c2690925 451 }
452 id1=mctrack0->GetFirstDaughter();
453 id2=mctrack0->GetLastDaughter();
454 if(!((id2-id1)==1)) continue;
455 for(int idx=id1; idx<=id2; idx++){
456 if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
9250ffbf 457 if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
c2690925 458 if(abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
8c1c76e9 459 fMCQACollection->Fill(Form("phidaughters_centrbin%i",fCentrality),mctrackd->Pt());
c2690925 460 }
461 }
462 else if(abs(mctrack0->PdgCode()) == 331) // eta prime
463 {
464 if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
8c1c76e9 465 fMCQACollection->Fill(Form("etapspectra_centrbin%i",fCentrality),mctrack0->Pt());
466 fMCQACollection->Fill(Form("etapspectraLog_centrbin%i",fCentrality),mctrack0->Pt());
c2690925 467 }
468 id1=mctrack0->GetFirstDaughter();
469 id2=mctrack0->GetLastDaughter();
470 if(!((id2-id1)==2||(id2-id1)==3)) continue;
471 for(int idx=id1; idx<=id2; idx++){
472 if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
9250ffbf 473 if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
c2690925 474 if(abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
8c1c76e9 475 fMCQACollection->Fill(Form("etapdaughters_centrbin%i",fCentrality),mctrackd->Pt());
c2690925 476 }
477 }
478 else if(abs(mctrack0->PdgCode()) == 113) // rho
479 {
480 if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
8c1c76e9 481 fMCQACollection->Fill(Form("rhospectra_centrbin%i",fCentrality),mctrack0->Pt());
482 fMCQACollection->Fill(Form("rhospectraLog_centrbin%i",fCentrality),mctrack0->Pt());
c2690925 483 }
484 id1=mctrack0->GetFirstDaughter();
485 id2=mctrack0->GetLastDaughter();
486 if(!((id2-id1)==1)) continue;
487 for(int idx=id1; idx<=id2; idx++){
488 if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
9250ffbf 489 if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
c2690925 490 if(abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
8c1c76e9 491 fMCQACollection->Fill(Form("rhodaughters_centrbin%i",fCentrality),mctrackd->Pt());
492 }
493 }
494 else if(abs(mctrack0->PdgCode()) == 321) // kaon+-
495 {
496 if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
497 fMCQACollection->Fill(Form("kaonspectraLog_centrbin%i",fCentrality),mctrack0->Pt());
498 }
499 }
500 else if(abs(mctrack0->PdgCode()) == 130) // k0L
501 {
502 if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
503 fMCQACollection->Fill(Form("k0LspectraLog_centrbin%i",fCentrality),mctrack0->Pt());
c2690925 504 }
505 }
506 }
507
508}
259c3296 509//__________________________________________
722347d8 510void AliHFEmcQA::GetQuarkKine(TParticle *part, Int_t iTrack, const Int_t kquark)
259c3296 511{
512 // get heavy quark kinematics
513
dbe3abbe 514 if (kquark != kCharm && kquark != kBeauty) {
259c3296 515 AliDebug(1, "This task is only for heavy quark QA, return\n");
516 return;
517 }
dbe3abbe 518 Int_t iq = kquark - kCharm;
259c3296 519
faee3b18 520 if (iTrack < 0 || !part) {
521 AliDebug(1, "Stack label is negative or no mcparticle, return\n");
259c3296 522 return;
523 }
524
faee3b18 525 AliMCParticle *mctrack = NULL;
259c3296 526 Int_t partPdgcode = TMath::Abs(part->GetPdgCode());
527
528 // select heavy hadron or not fragmented heavy quark
529 if ( int(partPdgcode/100.)==kquark || int(partPdgcode/1000.)==kquark || (partPdgcode==kquark && (part->GetNDaughters()==0 && iTrack>5)) ){
530
531 TParticle *partMother;
532 Int_t iLabel;
533
534 if (partPdgcode == kquark){ // in case of not fragmented heavy quark
535 partMother = part;
536 iLabel = iTrack;
537 } else{ // in case of heavy hadron, start to search for mother heavy parton
538 iLabel = part->GetFirstMother();
faee3b18 539 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return;
540 if (iLabel>-1) { partMother = mctrack->Particle(); }
259c3296 541 else {
542 AliDebug(1, "Stack label is negative, return\n");
543 return;
544 }
545 }
546
547 // heavy parton selection as a mother of heavy hadron
548 // if the heavy particle comes from string which is denoted as particle status 12|12|12...12|11,[PYTHIA p.60]
549 // in this case, the mother of heavy particle can be one of the fragmented parton of the string
550 // should I make a condition that partMother should be quark or diquark? -> not necessary
551 if ( abs(partMother->GetPdgCode()) == kquark || (partMother->GetStatusCode() == 11) ){
552 //if ( abs(partMother->GetPdgCode()) == kquark || (partMother->GetStatusCode() == 11 || partMother->GetStatusCode() == 12) ){
553
554 if ( abs(partMother->GetPdgCode()) != kquark ){
555 // search fragmented partons in the same string
dbe3abbe 556 Bool_t isSameString = kTRUE;
259c3296 557 for (Int_t i=1; i<fgkMaxIter; i++){
558 iLabel = iLabel - 1;
faee3b18 559 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return;
560 if (iLabel>-1) { partMother = mctrack->Particle(); }
259c3296 561 else {
562 AliDebug(1, "Stack label is negative, return\n");
563 return;
564 }
565 if ( abs(partMother->GetPdgCode()) == kquark ) break;
dbe3abbe 566 if ( partMother->GetStatusCode() != 12 ) isSameString = kFALSE;
567 if (!isSameString) return;
259c3296 568 }
569 }
570 AliDebug(1, "Can not find heavy parton of this heavy hadron in the string, return\n");
571 if (abs(partMother->GetPdgCode()) != kquark) return;
572
573 if (fIsHeavy[iq] >= 50) return;
574 fHeavyQuark[fIsHeavy[iq]] = partMother;
575 fIsHeavy[iq]++;
576
577 // fill kinematics for heavy parton
578 if (partMother->GetPdgCode() > 0) { // quark
dbe3abbe 579 fHist[iq][kQuark][0].fPdgCode->Fill(partMother->GetPdgCode());
580 fHist[iq][kQuark][0].fPt->Fill(partMother->Pt());
70da6c5a 581 fHist[iq][kQuark][0].fY->Fill(AliHFEtools::GetRapidity(partMother));
dbe3abbe 582 fHist[iq][kQuark][0].fEta->Fill(partMother->Eta());
259c3296 583 } else{ // antiquark
dbe3abbe 584 fHist[iq][kantiQuark][0].fPdgCode->Fill(partMother->GetPdgCode());
585 fHist[iq][kantiQuark][0].fPt->Fill(partMother->Pt());
70da6c5a 586 fHist[iq][kantiQuark][0].fY->Fill(AliHFEtools::GetRapidity(partMother));
dbe3abbe 587 fHist[iq][kantiQuark][0].fEta->Fill(partMother->Eta());
259c3296 588 }
589
590 } // end of heavy parton slection loop
591
592 } // end of heavy hadron or quark selection
593
594}
595
596//__________________________________________
597void AliHFEmcQA::EndOfEventAna(const Int_t kquark)
598{
599 // end of event analysis
600
dbe3abbe 601 if (kquark != kCharm && kquark != kBeauty) {
259c3296 602 AliDebug(1, "This task is only for heavy quark QA, return\n");
603 return;
604 }
dbe3abbe 605 Int_t iq = kquark - kCharm;
259c3296 606
607
608 // # of heavy quark per event
609 AliDebug(1,Form("Number of heavy quark in this event = %d \n",fIsHeavy[iq]));
dbe3abbe 610 fHistComm[iq][0].fNq->Fill(fIsHeavy[iq]);
259c3296 611
612 Int_t motherID[fgkMaxGener];
613 Int_t motherType[fgkMaxGener];
614 Int_t motherLabel[fgkMaxGener];
615 Int_t ancestorPdg[fgkMaxGener];
616 Int_t ancestorLabel[fgkMaxGener];
617
618 for (Int_t i = 0; i < fgkMaxGener; i++){ // initialization
619 motherID[i] = 0;
620 motherType[i] = 0;
621 motherLabel[i] = 0;
622 ancestorPdg[i] = 0;
623 ancestorLabel[i] = 0;
624 }
625
3a72645a 626
259c3296 627 // check history of found heavy quarks
628 for (Int_t i = 0; i < fIsHeavy[iq]; i++){
629
3a72645a 630 if(!fHeavyQuark[i]) return;
631
259c3296 632 ancestorLabel[0] = i;
633 ancestorPdg[0] = fHeavyQuark[i]->GetPdgCode();
634 ancestorLabel[1] = fHeavyQuark[i]->GetFirstMother();
635
636 AliDebug(1,Form("pdg code= %d\n",ancestorPdg[0]));
637 AliDebug(1,Form("ancestor label= %d\n",ancestorLabel[1]));
638
639 Int_t ig = 1;
640 while (ancestorLabel[ig] != -1){
641 // in case there is mother, get mother's pdg code and grandmother's label
642 IdentifyMother(ancestorLabel[ig], ancestorPdg[ig], ancestorLabel[ig+1]);
643 // if mother is still heavy, find again mother's ancestor
644 if (ancestorPdg[ig-1] == ancestorPdg[ig]) {
645 ig++;
646 continue; // if it is from same heavy
647 }
648 // if the heavy's mother is not heavy, check the mother's label to know if it comes from inital or final parton shower
649 if (IsFromInitialShower(ancestorLabel[ig],motherID[i],motherType[i],motherLabel[i])) break;
650 if (IsFromFinalParton(ancestorLabel[ig],motherID[i],motherType[i],motherLabel[i])) break;
651 // if it is not the above case, something is strange
dbe3abbe 652 ReportStrangeness(motherID[i],motherType[i],motherLabel[i]);
259c3296 653 break;
654 }
655 if (ancestorLabel[ig] == -1){ // from hard scattering
656 HardScattering(kquark, motherID[i],motherType[i], motherLabel[i]);
657 }
658
659 } // end of found heavy quark loop
660
661
662 // check process type
663 Int_t processID = 0;
664 for (Int_t i = 0; i < fIsHeavy[iq]; i++){
665 AliDebug(1,Form("Mother ID= %d type= %d label= %d\n",motherID[i],motherType[i],motherLabel[i]));
666 }
667
668
669 Int_t nheavypair = Int_t(fIsHeavy[iq]/2.);
670 for (Int_t ipair = 0; ipair < nheavypair; ipair++){
671
672 Int_t id1 = ipair*2;
673 Int_t id2 = ipair*2 + 1;
674
675 if (motherType[id1] == 2 && motherType[id2] == 2){
dbe3abbe 676 if (motherLabel[id1] == motherLabel[id2]) processID = kGluonSplitting; // gluon spliting
259c3296 677 else processID = -9;
678 }
679 else if (motherType[id1] == -1 && motherType[id2] == -1) {
680 if (motherLabel[id1] == -1 && motherLabel[id2] == -1) {
dbe3abbe 681 if (motherID[id1] == fgkGluon) processID = kPairCreationFromg; // gluon fusion
682 else processID = kPairCreationFromq; // q-qbar pair creation
259c3296 683 }
684 else processID = -8;
685 }
686 else if (motherType[id1] == -1 || motherType[id2] == -1) {
687 if ((motherLabel[id1] == -1 || motherLabel[id2] == -1) && (motherLabel[id1]*motherLabel[id2] == -2 || motherLabel[id1]*motherLabel[id2] == -3)) {
dbe3abbe 688 if(motherID[id1]*motherID[id2] == kquark*fgkGluon) processID = kFlavourExitation; // flavour exitation
689 else processID = kLightQuarkShower;
259c3296 690 }
691 else processID = -7;
692 }
693 else if (motherType[id1] == -2 || motherType[id2] == -2) {
dbe3abbe 694 if (motherLabel[id1] == motherLabel[id2]) processID = kInitialPartonShower; // initial parton shower
259c3296 695 else processID = -6;
696
697 }
698 else processID = -5;
699
700 if (nheavypair >1) AliDebug(1,Form("Multi pair found : process ID = %d\n",processID));
dbe3abbe 701 else fHistComm[iq][0].fProcessID->Fill(processID);
259c3296 702 AliDebug(1,Form("Process ID = %d\n",processID));
703 } // end of # heavy quark pair loop
704
705}
706
707//__________________________________________
722347d8 708void AliHFEmcQA::GetHadronKine(TParticle* mcpart, const Int_t kquark)
dbe3abbe 709{
710 // decay electron kinematics
711
712 if (kquark != kCharm && kquark != kBeauty) {
713 AliDebug(1, "This task is only for heavy quark QA, return\n");
714 return;
715 }
716 Int_t iq = kquark - kCharm;
717
faee3b18 718 if(!mcpart){
719 AliDebug(1, "no mc particle, return\n");
720 return;
721 }
dbe3abbe 722
723 Int_t iLabel = mcpart->GetFirstMother();
724 if (iLabel<0){
725 AliDebug(1, "Stack label is negative, return\n");
726 return;
727 }
728
729 TParticle *partCopy = mcpart;
730 Int_t pdgcode = mcpart->GetPdgCode();
731 Int_t pdgcodeCopy = pdgcode;
732
faee3b18 733 AliMCParticle *mctrack = NULL;
734
dbe3abbe 735 // if the mother is charmed hadron
75d81601 736 Bool_t isDirectCharm = kFALSE;
dbe3abbe 737 if ( int(abs(pdgcode)/100.) == kCharm || int(abs(pdgcode)/1000.) == kCharm ) {
738
739 // iterate until you find B hadron as a mother or become top ancester
740 for (Int_t i=1; i<fgkMaxIter; i++){
741
742 Int_t jLabel = mcpart->GetFirstMother();
743 if (jLabel == -1){
75d81601 744 isDirectCharm = kTRUE;
dbe3abbe 745 break; // if there is no ancester
746 }
747 if (jLabel < 0){ // safety protection
748 AliDebug(1, "Stack label is negative, return\n");
749 return;
750 }
751 // if there is an ancester
faee3b18 752 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return;
753 TParticle* mother = mctrack->Particle();
dbe3abbe 754 Int_t motherPDG = mother->GetPdgCode();
755
756 for (Int_t j=0; j<fNparents; j++){
757 if (abs(motherPDG)==fParentSelect[1][j]) return; // return if this hadron is originated from b
758 }
759
760 mcpart = mother;
761 } // end of iteration
762 } // end of if
75d81601 763 if((isDirectCharm == kTRUE && kquark == kCharm) || kquark == kBeauty) {
dbe3abbe 764 for (Int_t i=0; i<fNparents; i++){
765 if (abs(pdgcodeCopy)==fParentSelect[iq][i]){
766
767 // fill hadron kinematics
768 fHist[iq][kHadron][0].fPdgCode->Fill(pdgcodeCopy);
769 fHist[iq][kHadron][0].fPt->Fill(partCopy->Pt());
70da6c5a 770 fHist[iq][kHadron][0].fY->Fill(AliHFEtools::GetRapidity(partCopy));
dbe3abbe 771 fHist[iq][kHadron][0].fEta->Fill(partCopy->Eta());
772
ccc37cdc 773 if(iq==0) {
774 fhD[i]->Fill(partCopy->Pt(),AliHFEtools::GetRapidity(partCopy));
775 if(TMath::Abs(AliHFEtools::GetRapidity(partCopy))<0.5) fhDLogbin[i]->Fill(partCopy->Pt());
776 }
dbe3abbe 777 }
778 }
ccc37cdc 779 // I also want to store D* info to compare with D* measurement
780 if (abs(pdgcodeCopy)==413 && iq==0) { //D*+
781 fhD[7]->Fill(partCopy->Pt(),AliHFEtools::GetRapidity(partCopy));
782 if(TMath::Abs(AliHFEtools::GetRapidity(partCopy))<0.5) fhDLogbin[7]->Fill(partCopy->Pt());
783 }
784 if (abs(pdgcodeCopy)==423 && iq==0) { //D*0
785 fhD[8]->Fill(partCopy->Pt(),AliHFEtools::GetRapidity(partCopy));
786 if(TMath::Abs(AliHFEtools::GetRapidity(partCopy))<0.5) fhDLogbin[8]->Fill(partCopy->Pt());
787 }
dbe3abbe 788 } // end of if
789}
790
791//__________________________________________
722347d8 792void AliHFEmcQA::GetDecayedKine(TParticle* mcpart, const Int_t kquark, Int_t kdecayed, Int_t icut)
259c3296 793{
794 // decay electron kinematics
795
faee3b18 796 if (!(kquark == kCharm || kquark == kBeauty || kquark == kOthers)){
259c3296 797 AliDebug(1, "This task is only for heavy quark QA, return\n");
798 return;
799 }
dbe3abbe 800 Int_t iq = kquark - kCharm;
50685501 801 Bool_t isFinalOpenCharm = kFALSE;
259c3296 802
faee3b18 803 if(!mcpart){
804 AliDebug(1, "no mcparticle, return\n");
805 return;
259c3296 806 }
807
faee3b18 808 if(kquark==kOthers){
809 Int_t esource = -1;
810 if ( abs(mcpart->GetPdgCode()) != kdecayed ) esource = kMisID-4;
811 else esource =GetSource(mcpart)-4; // return for the cases kGamma=4, kPi0=5, kElse=6
812 if(esource==0|| esource==1 || esource==2 || esource==3){
813 fHist[iq][esource][icut].fPt->Fill(mcpart->Pt());
814 fHist[iq][esource][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart));
815 fHist[iq][esource][icut].fEta->Fill(mcpart->Eta());
816 return;
817 }
818 else {
819 AliDebug(1, "e source is out of defined ranges, return\n");
820 return;
821 }
822 }
259c3296 823
824 if ( abs(mcpart->GetPdgCode()) != kdecayed ) return;
259c3296 825
826 Int_t iLabel = mcpart->GetFirstMother();
827 if (iLabel<0){
828 AliDebug(1, "Stack label is negative, return\n");
829 return;
830 }
831
faee3b18 832 AliMCParticle *mctrack = NULL;
833 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return;
834 TParticle *partMother = mctrack->Particle();
dbe3abbe 835 TParticle *partMotherCopy = partMother;
259c3296 836 Int_t maPdgcode = partMother->GetPdgCode();
dbe3abbe 837 Int_t maPdgcodeCopy = maPdgcode;
259c3296 838
9bcfd1ab 839 // get mc primary vertex
faee3b18 840 /*
9bcfd1ab 841 TArrayF mcPrimVtx;
842 if(fMCHeader) fMCHeader->PrimaryVertex(mcPrimVtx);
843
259c3296 844 // get electron production vertex
845 TLorentzVector ePoint;
846 mcpart->ProductionVertex(ePoint);
847
9bcfd1ab 848 // calculated production vertex to primary vertex (in xy plane)
849 Float_t decayLxy = TMath::Sqrt((mcPrimVtx[0]-ePoint[0])*(mcPrimVtx[0]-ePoint[0])+(mcPrimVtx[1]-ePoint[1])*(mcPrimVtx[1]-ePoint[1]));
faee3b18 850 */
851 Float_t decayLxy = 0;
9bcfd1ab 852
259c3296 853 // if the mother is charmed hadron
dbe3abbe 854 Bool_t isMotherDirectCharm = kFALSE;
855 if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
259c3296 856
0792aa82 857 for (Int_t i=0; i<fNparents; i++){
858 if (abs(maPdgcode)==fParentSelect[0][i]){
50685501 859 isFinalOpenCharm = kTRUE;
0792aa82 860 }
861 }
50685501 862 if (!isFinalOpenCharm) return ;
0792aa82 863
259c3296 864 // iterate until you find B hadron as a mother or become top ancester
865 for (Int_t i=1; i<fgkMaxIter; i++){
866
867 Int_t jLabel = partMother->GetFirstMother();
868 if (jLabel == -1){
dbe3abbe 869 isMotherDirectCharm = kTRUE;
259c3296 870 break; // if there is no ancester
871 }
872 if (jLabel < 0){ // safety protection
873 AliDebug(1, "Stack label is negative, return\n");
874 return;
875 }
876
877 // if there is an ancester
faee3b18 878 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return;
879 TParticle* grandMa = mctrack->Particle();
259c3296 880 Int_t grandMaPDG = grandMa->GetPdgCode();
881
fc0de2a0 882 for (Int_t j=0; j<fNparents; j++){
883 if (abs(grandMaPDG)==fParentSelect[1][j]){
259c3296 884
dbe3abbe 885 if (kquark == kCharm) return;
259c3296 886 // fill electron kinematics
dbe3abbe 887 fHist[iq][kElectron2nd][icut].fPdgCode->Fill(mcpart->GetPdgCode());
888 fHist[iq][kElectron2nd][icut].fPt->Fill(mcpart->Pt());
70da6c5a 889 fHist[iq][kElectron2nd][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart));
dbe3abbe 890 fHist[iq][kElectron2nd][icut].fEta->Fill(mcpart->Eta());
259c3296 891
892 // fill mother hadron kinematics
dbe3abbe 893 fHist[iq][kDeHadron][icut].fPdgCode->Fill(grandMaPDG);
894 fHist[iq][kDeHadron][icut].fPt->Fill(grandMa->Pt());
70da6c5a 895 fHist[iq][kDeHadron][icut].fY->Fill(AliHFEtools::GetRapidity(grandMa));
dbe3abbe 896 fHist[iq][kDeHadron][icut].fEta->Fill(grandMa->Eta());
259c3296 897
898 // ratio between pT of electron and pT of mother B hadron
dbe3abbe 899 if(grandMa->Pt()) fHistComm[iq][icut].fDePtRatio->Fill(grandMa->Pt(),mcpart->Pt()/grandMa->Pt());
259c3296 900
9bcfd1ab 901 // distance between electron production point and primary vertex
902 fHistComm[iq][icut].fDeDistance->Fill(grandMa->Pt(),decayLxy);
259c3296 903 return;
904 }
dc306130 905 }
259c3296 906
907 partMother = grandMa;
908 } // end of iteration
909 } // end of if
dbe3abbe 910 if((isMotherDirectCharm == kTRUE && kquark == kCharm) || kquark == kBeauty) {
259c3296 911 for (Int_t i=0; i<fNparents; i++){
912 if (abs(maPdgcodeCopy)==fParentSelect[iq][i]){
913
e97c2edf 914//mj weighting to consider measured spectra!!!
915 Double_t mpt=partMotherCopy->Pt();
ccc37cdc 916 Double_t wfactor=2*(703.681*mpt/TMath::Power((1+TMath::Power(mpt/1.73926,2)),2.34821))/(368.608*mpt/TMath::Power((1+TMath::Power(mpt/2.74868,2)),2.34225)); // 2* considering in pythia I used particle + antiparticle differently from the measurement
917 //Double_t wfactor=(703.681*mpt/TMath::Power((1+TMath::Power(mpt/1.73926,2)),2.34821))/(368.608*mpt/TMath::Power((1+TMath::Power(mpt/2.74868,2)),2.34225));
259c3296 918 // fill electron kinematics
e97c2edf 919 if(iq==0){
920 fHist[iq][kElectron2nd][icut].fPdgCode->Fill(mcpart->GetPdgCode(),wfactor);
921 fHist[iq][kElectron2nd][icut].fPt->Fill(mcpart->Pt(),wfactor);
922 fHist[iq][kElectron2nd][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart),wfactor);
923 fHist[iq][kElectron2nd][icut].fEta->Fill(mcpart->Eta(),wfactor);
924
925 fHist[iq][kElectron][icut].fPdgCode->Fill(mcpart->GetPdgCode());
926 fHist[iq][kElectron][icut].fPt->Fill(mcpart->Pt());
927 fHist[iq][kElectron][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart));
928 fHist[iq][kElectron][icut].fEta->Fill(mcpart->Eta());
929 }
930 else{
931 fHist[iq][kElectron][icut].fPdgCode->Fill(mcpart->GetPdgCode());
932 fHist[iq][kElectron][icut].fPt->Fill(mcpart->Pt());
933 fHist[iq][kElectron][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart));
934 fHist[iq][kElectron][icut].fEta->Fill(mcpart->Eta());
935 }
936//--------------
259c3296 937 // fill mother hadron kinematics
dbe3abbe 938 fHist[iq][keHadron][icut].fPdgCode->Fill(maPdgcodeCopy);
939 fHist[iq][keHadron][icut].fPt->Fill(partMotherCopy->Pt());
70da6c5a 940 fHist[iq][keHadron][icut].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
dbe3abbe 941 fHist[iq][keHadron][icut].fEta->Fill(partMotherCopy->Eta());
259c3296 942
ccc37cdc 943 // ratio between pT of electron and pT of mother B or direct D hadron
dbe3abbe 944 if(partMotherCopy->Pt()) fHistComm[iq][icut].fePtRatio->Fill(partMotherCopy->Pt(),mcpart->Pt()/partMotherCopy->Pt());
ccc37cdc 945 fHistComm[iq][icut].fPtCorr->Fill(partMotherCopy->Pt(),mcpart->Pt());
c2690925 946 if(TMath::Abs(partMotherCopy->GetPdgCode())==411) fHistComm[iq][icut].fPtCorrDp->Fill(partMotherCopy->Pt(),mcpart->Pt());
947 else if(TMath::Abs(partMotherCopy->GetPdgCode())==421) fHistComm[iq][icut].fPtCorrD0->Fill(partMotherCopy->Pt(),mcpart->Pt());
948 else fHistComm[iq][icut].fPtCorrDrest->Fill(partMotherCopy->Pt(),mcpart->Pt());
259c3296 949
9bcfd1ab 950 // distance between electron production point and primary vertex
951 fHistComm[iq][icut].feDistance->Fill(partMotherCopy->Pt(),decayLxy);
259c3296 952 }
953 }
954 } // end of if
955}
956
0792aa82 957//____________________________________________________________________
958void AliHFEmcQA::GetDecayedKine(AliAODMCParticle *mcpart, const Int_t kquark, Int_t kdecayed, Int_t icut)
959{
960 // decay electron kinematics
961
962 if (kquark != kCharm && kquark != kBeauty) {
963 AliDebug(1, "This task is only for heavy quark QA, return\n");
964 return;
965 }
966
967 Int_t iq = kquark - kCharm;
50685501 968 Bool_t isFinalOpenCharm = kFALSE;
0792aa82 969
faee3b18 970 if(!mcpart){
971 AliDebug(1, "no mcparticle, return\n");
972 return;
973 }
974
0792aa82 975 if ( abs(mcpart->GetPdgCode()) != kdecayed ) return;
976
977 // mother
978 Int_t iLabel = mcpart->GetMother();
979 if (iLabel<0){
980 AliDebug(1, "Stack label is negative, return\n");
981 return;
982 }
983
984 AliAODMCParticle *partMother = (AliAODMCParticle*)fMCArray->At(iLabel);
985 AliAODMCParticle *partMotherCopy = partMother;
986 Int_t maPdgcode = partMother->GetPdgCode();
987 Int_t maPdgcodeCopy = maPdgcode;
988
989 Bool_t isMotherDirectCharm = kFALSE;
990 if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
991
992 for (Int_t i=0; i<fNparents; i++){
993 if (abs(maPdgcode)==fParentSelect[0][i]){
50685501 994 isFinalOpenCharm = kTRUE;
0792aa82 995 }
996 }
50685501 997 if (!isFinalOpenCharm) return;
0792aa82 998
999 for (Int_t i=1; i<fgkMaxIter; i++){
1000
1001 Int_t jLabel = partMother->GetMother();
1002 if (jLabel == -1){
1003 isMotherDirectCharm = kTRUE;
1004 break; // if there is no ancester
1005 }
1006 if (jLabel < 0){ // safety protection
1007 AliDebug(1, "Stack label is negative, return\n");
1008 return;
1009 }
1010
1011 // if there is an ancester
1012 AliAODMCParticle* grandMa = (AliAODMCParticle*)fMCArray->At(jLabel);
1013 Int_t grandMaPDG = grandMa->GetPdgCode();
1014
1015 for (Int_t j=0; j<fNparents; j++){
1016 if (abs(grandMaPDG)==fParentSelect[1][j]){
1017
1018 if (kquark == kCharm) return;
1019 // fill electron kinematics
1020 fHist[iq][kElectron2nd][icut].fPdgCode->Fill(mcpart->GetPdgCode());
1021 fHist[iq][kElectron2nd][icut].fPt->Fill(mcpart->Pt());
70da6c5a 1022 fHist[iq][kElectron2nd][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart));
0792aa82 1023 fHist[iq][kElectron2nd][icut].fEta->Fill(mcpart->Eta());
1024
1025 // fill mother hadron kinematics
1026 fHist[iq][kDeHadron][icut].fPdgCode->Fill(grandMaPDG);
1027 fHist[iq][kDeHadron][icut].fPt->Fill(grandMa->Pt());
70da6c5a 1028 fHist[iq][kDeHadron][icut].fY->Fill(AliHFEtools::GetRapidity(grandMa));
0792aa82 1029 fHist[iq][kDeHadron][icut].fEta->Fill(grandMa->Eta());
1030
1031 return;
1032 }
1033 }
1034
1035 partMother = grandMa;
1036 } // end of iteration
1037 } // end of if
1038 if ((isMotherDirectCharm == kTRUE && kquark == kCharm) || kquark == kBeauty) {
1039 for (Int_t i=0; i<fNparents; i++){
1040 if (abs(maPdgcodeCopy)==fParentSelect[iq][i]){
1041
1042 // fill electron kinematics
1043 fHist[iq][kElectron][icut].fPdgCode->Fill(mcpart->GetPdgCode());
1044 fHist[iq][kElectron][icut].fPt->Fill(mcpart->Pt());
70da6c5a 1045 fHist[iq][kElectron][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart));
0792aa82 1046 fHist[iq][kElectron][icut].fEta->Fill(mcpart->Eta());
1047
1048 // fill mother hadron kinematics
1049 fHist[iq][keHadron][icut].fPdgCode->Fill(maPdgcodeCopy);
1050 fHist[iq][keHadron][icut].fPt->Fill(partMotherCopy->Pt());
70da6c5a 1051 fHist[iq][keHadron][icut].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
0792aa82 1052 fHist[iq][keHadron][icut].fEta->Fill(partMotherCopy->Eta());
1053
1054 }
1055 }
1056 } // end of if
1057
1058}
259c3296 1059
1060//__________________________________________
75d81601 1061void AliHFEmcQA::IdentifyMother(Int_t motherlabel, Int_t &motherpdg, Int_t &grandmotherlabel)
259c3296 1062{
dbe3abbe 1063 // find mother pdg code and label
1064
75d81601 1065 if (motherlabel < 0) {
259c3296 1066 AliDebug(1, "Stack label is negative, return\n");
1067 return;
1068 }
faee3b18 1069 AliMCParticle *mctrack = NULL;
1070 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(motherlabel))))) return;
1071 TParticle *heavysMother = mctrack->Particle();
75d81601 1072 motherpdg = heavysMother->GetPdgCode();
1073 grandmotherlabel = heavysMother->GetFirstMother();
1074 AliDebug(1,Form("ancestor pdg code= %d\n",motherpdg));
259c3296 1075}
1076
1077//__________________________________________
259c3296 1078void AliHFEmcQA::HardScattering(const Int_t kquark, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
1079{
dbe3abbe 1080 // mothertype -1 means this heavy quark coming from hard vertex
1081
faee3b18 1082 AliMCParticle *mctrack1 = NULL;
1083 AliMCParticle *mctrack2 = NULL;
1084 if(!(mctrack1 = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(4))))) return;
1085 if(!(mctrack2 = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(5))))) return;
1086 TParticle *afterinitialrad1 = mctrack1->Particle();
1087 TParticle *afterinitialrad2 = mctrack2->Particle();
259c3296 1088
1089 motherlabel = -1;
1090
1091 if (abs(afterinitialrad1->GetPdgCode()) == fgkGluon && abs(afterinitialrad2->GetPdgCode()) == fgkGluon){
1092 AliDebug(1,"heavy from gluon gluon pair creation!\n");
1093 mothertype = -1;
1094 motherID = fgkGluon;
1095 }
1096 else if (abs(afterinitialrad1->GetPdgCode()) == kquark || abs(afterinitialrad2->GetPdgCode()) == kquark){ // one from Q and the other from g
1097 AliDebug(1,"heavy from flavor exitation!\n");
1098 mothertype = -1;
1099 motherID = kquark;
1100 }
1101 else if (abs(afterinitialrad1->GetPdgCode()) == abs(afterinitialrad2->GetPdgCode())){
1102 AliDebug(1,"heavy from q-qbar pair creation!\n");
1103 mothertype = -1;
1104 motherID = 1;
1105 }
1106 else {
1107 AliDebug(1,"something strange!\n");
1108 mothertype = -999;
1109 motherlabel = -999;
1110 motherID = -999;
1111 }
1112}
1113
1114//__________________________________________
259c3296 1115Bool_t AliHFEmcQA::IsFromInitialShower(Int_t inputmotherlabel, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
1116{
dbe3abbe 1117 // mothertype -2 means this heavy quark coming from initial state
1118
faee3b18 1119 AliMCParticle *mctrack = NULL;
259c3296 1120 if (inputmotherlabel==2 || inputmotherlabel==3){ // mother exist before initial state radiation
faee3b18 1121 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(inputmotherlabel))))) return kFALSE;
1122 TParticle *heavysMother = mctrack->Particle();
259c3296 1123 motherID = heavysMother->GetPdgCode();
1124 mothertype = -2; // there is mother before initial state radiation
1125 motherlabel = inputmotherlabel;
1126 AliDebug(1,"initial parton shower! \n");
1127
1128 return kTRUE;
1129 }
1130
1131 return kFALSE;
1132}
1133
1134//__________________________________________
259c3296 1135Bool_t AliHFEmcQA::IsFromFinalParton(Int_t inputmotherlabel, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
1136{
dbe3abbe 1137 // mothertype 2 means this heavy quark coming from final state
1138
faee3b18 1139 AliMCParticle *mctrack = NULL;
259c3296 1140 if (inputmotherlabel > 5){ // mother exist after hard scattering
faee3b18 1141 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(inputmotherlabel))))) return kFALSE;
1142 TParticle *heavysMother = mctrack->Particle();
259c3296 1143 motherID = heavysMother->GetPdgCode();
1144 mothertype = 2; //
1145 motherlabel = inputmotherlabel;
1146 AliDebug(1,Form("heavy quark from %d after hard scattering! \n",motherID));
1147
1148 return kTRUE;
1149 }
1150 return kFALSE;
1151}
1152
1153//__________________________________________
dbe3abbe 1154void AliHFEmcQA::ReportStrangeness(Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
259c3296 1155{
dbe3abbe 1156 // mark strange behavior
1157
259c3296 1158 mothertype = -888;
1159 motherlabel = -888;
1160 motherID = -888;
1161 AliDebug(1,"something strange!\n");
1162}
1163
0792aa82 1164//__________________________________________
c2690925 1165Int_t AliHFEmcQA::GetSource(const AliAODMCParticle * const mcpart)
0792aa82 1166{
9bcfd1ab 1167 // decay particle's origin
0792aa82 1168
9bcfd1ab 1169 //if ( abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return -1;
0792aa82 1170
1171 Int_t origin = -1;
50685501 1172 Bool_t isFinalOpenCharm = kFALSE;
0792aa82 1173
faee3b18 1174 if(!mcpart){
1175 AliDebug(1, "Stack label is negative or no mcparticle, return\n");
1176 return -1;
1177 }
1178
0792aa82 1179 // mother
1180 Int_t iLabel = mcpart->GetMother();
1181 if (iLabel<0){
1182 AliDebug(1, "Stack label is negative, return\n");
1183 return -1;
1184 }
1185
1186 AliAODMCParticle *partMother = (AliAODMCParticle*)fMCArray->At(iLabel);
1187 Int_t maPdgcode = partMother->GetPdgCode();
1188
1189 // if the mother is charmed hadron
1190 if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
1191
1192 for (Int_t i=0; i<fNparents; i++){
1193 if (abs(maPdgcode)==fParentSelect[0][i]){
50685501 1194 isFinalOpenCharm = kTRUE;
0792aa82 1195 }
1196 }
50685501 1197 if (!isFinalOpenCharm) return -1;
0792aa82 1198
1199 // iterate until you find B hadron as a mother or become top ancester
1200 for (Int_t i=1; i<fgkMaxIter; i++){
1201
1202 Int_t jLabel = partMother->GetMother();
1203 if (jLabel == -1){
1204 origin = kDirectCharm;
1205 return origin;
1206 }
1207 if (jLabel < 0){ // safety protection
1208 AliDebug(1, "Stack label is negative, return\n");
1209 return -1;
1210 }
1211
1212 // if there is an ancester
1213 AliAODMCParticle* grandMa = (AliAODMCParticle*)fMCArray->At(jLabel);
1214 Int_t grandMaPDG = grandMa->GetPdgCode();
1215
70b5ea26 1216 for (Int_t j=0; j<fNparents; j++){
1217 if (abs(grandMaPDG)==fParentSelect[1][j]){
0792aa82 1218 origin = kBeautyCharm;
1219 return origin;
1220 }
1221 }
1222
1223 partMother = grandMa;
1224 } // end of iteration
1225 } // end of if
1226 else if ( int(abs(maPdgcode)/100.) == kBeauty || int(abs(maPdgcode)/1000.) == kBeauty ) {
1227 for (Int_t i=0; i<fNparents; i++){
1228 if (abs(maPdgcode)==fParentSelect[1][i]){
1229 origin = kDirectBeauty;
1230 return origin;
1231 }
1232 }
1233 } // end of if
1234 else if ( abs(maPdgcode) == 22 ) {
1235 origin = kGamma;
1236 return origin;
1237 } // end of if
1238 else if ( abs(maPdgcode) == 111 ) {
1239 origin = kPi0;
1240 return origin;
1241 } // end of if
1242
1243 return origin;
1244}
1245
1246//__________________________________________
c2690925 1247Int_t AliHFEmcQA::GetSource(const TParticle * const mcpart)
0792aa82 1248{
9bcfd1ab 1249 // decay particle's origin
0792aa82 1250
9bcfd1ab 1251 //if ( abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return -1;
0792aa82 1252
1253 Int_t origin = -1;
50685501 1254 Bool_t isFinalOpenCharm = kFALSE;
0792aa82 1255
faee3b18 1256 if(!mcpart){
1257 AliDebug(1, "no mcparticle, return\n");
1258 return -1;
1259 }
1260
0792aa82 1261 Int_t iLabel = mcpart->GetFirstMother();
1262 if (iLabel<0){
1263 AliDebug(1, "Stack label is negative, return\n");
1264 return -1;
1265 }
1266
faee3b18 1267 AliMCParticle *mctrack = NULL;
1268 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return -1;
1269 TParticle *partMother = mctrack->Particle();
0792aa82 1270 Int_t maPdgcode = partMother->GetPdgCode();
1271
1272 // if the mother is charmed hadron
1273 if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
1274
1275 for (Int_t i=0; i<fNparents; i++){
1276 if (abs(maPdgcode)==fParentSelect[0][i]){
50685501 1277 isFinalOpenCharm = kTRUE;
0792aa82 1278 }
1279 }
50685501 1280 if (!isFinalOpenCharm) return -1;
0792aa82 1281
1282 // iterate until you find B hadron as a mother or become top ancester
1283 for (Int_t i=1; i<fgkMaxIter; i++){
1284
1285 Int_t jLabel = partMother->GetFirstMother();
1286 if (jLabel == -1){
1287 origin = kDirectCharm;
1288 return origin;
1289 }
1290 if (jLabel < 0){ // safety protection
1291 AliDebug(1, "Stack label is negative, return\n");
1292 return -1;
1293 }
1294
1295 // if there is an ancester
faee3b18 1296 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return -1;
1297 TParticle* grandMa = mctrack->Particle();
0792aa82 1298 Int_t grandMaPDG = grandMa->GetPdgCode();
1299
70b5ea26 1300 for (Int_t j=0; j<fNparents; j++){
1301 if (abs(grandMaPDG)==fParentSelect[1][j]){
0792aa82 1302 origin = kBeautyCharm;
1303 return origin;
1304 }
1305 }
1306
1307 partMother = grandMa;
1308 } // end of iteration
1309 } // end of if
1310 else if ( int(abs(maPdgcode)/100.) == kBeauty || int(abs(maPdgcode)/1000.) == kBeauty ) {
1311 for (Int_t i=0; i<fNparents; i++){
1312 if (abs(maPdgcode)==fParentSelect[1][i]){
1313 origin = kDirectBeauty;
1314 return origin;
1315 }
1316 }
1317 } // end of if
1318 else if ( abs(maPdgcode) == 22 ) {
1319 origin = kGamma;
1320 return origin;
1321 } // end of if
1322 else if ( abs(maPdgcode) == 111 ) {
1323 origin = kPi0;
1324 return origin;
1325 } // end of if
1326 else origin = kElse;
1327
1328 return origin;
1329}
70da6c5a 1330
1331//__________________________________________
e3fc062d 1332Int_t AliHFEmcQA::GetElecSource(TParticle * const mcpart)
1333{
1334 // decay particle's origin
70da6c5a 1335
e3fc062d 1336 if(!mcpart){
1337 AliDebug(1, "no mcparticle, return\n");
1338 return -1;
1339 }
1340
bf892a6a 1341 if ( abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return kMisID;
1342
1343 Int_t origin = -1;
1344 Bool_t isFinalOpenCharm = kFALSE;
1345
e3fc062d 1346 Int_t iLabel = mcpart->GetFirstMother();
1347 if (iLabel<0){
1348 AliDebug(1, "Stack label is negative, return\n");
1349 return -1;
1350 }
1351
1352 AliMCParticle *mctrack = NULL;
c2690925 1353 Int_t tmpMomLabel=0;
e3fc062d 1354 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return -1;
1355 TParticle *partMother = mctrack->Particle();
c2690925 1356 TParticle *partMotherCopy = mctrack->Particle();
e3fc062d 1357 Int_t maPdgcode = partMother->GetPdgCode();
1358
1359 // if the mother is charmed hadron
1360 if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
1361
1362 for (Int_t i=0; i<fNparents; i++){
1363 if (abs(maPdgcode)==fParentSelect[0][i]){
1364 isFinalOpenCharm = kTRUE;
1365 }
1366 }
1367 if (!isFinalOpenCharm) return -1;
1368
1369 // iterate until you find B hadron as a mother or become top ancester
1370 for (Int_t i=1; i<fgkMaxIter; i++){
1371
1372 Int_t jLabel = partMother->GetFirstMother();
1373 if (jLabel == -1){
1374 origin = kDirectCharm;
1375 return origin;
1376 }
1377 if (jLabel < 0){ // safety protection
1378 AliDebug(1, "Stack label is negative, return\n");
1379 return -1;
1380 }
1381
1382 // if there is an ancester
1383 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return -1;
1384 TParticle* grandMa = mctrack->Particle();
1385 Int_t grandMaPDG = grandMa->GetPdgCode();
1386
1387 for (Int_t j=0; j<fNparents; j++){
1388 if (abs(grandMaPDG)==fParentSelect[1][j]){
1389 origin = kBeautyCharm;
1390 return origin;
1391 }
1392 }
1393
1394 partMother = grandMa;
1395 } // end of iteration
1396 } // end of if
1397 else if ( int(abs(maPdgcode)/100.) == kBeauty || int(abs(maPdgcode)/1000.) == kBeauty ) {
1398 for (Int_t i=0; i<fNparents; i++){
1399 if (abs(maPdgcode)==fParentSelect[1][i]){
1400 origin = kDirectBeauty;
1401 return origin;
1402 }
1403 }
1404 } // end of if
c2690925 1405 else if ( abs(maPdgcode) == 22 ) { //conversion
1406
1407 tmpMomLabel = partMotherCopy->GetFirstMother();
1408 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tmpMomLabel))))) return -1;
1409 partMother = mctrack->Particle();
1410 maPdgcode = partMother->GetPdgCode();
1411 if ( abs(maPdgcode) == 111 ) {
1412 origin = kGammaPi0;
1413 return origin;
1414 }
1415 else if ( abs(maPdgcode) == 221 ) {
1416 origin = kGammaEta;
1417 return origin;
1418 }
1419 else if ( abs(maPdgcode) == 223 ) {
1420 origin = kGammaOmega;
1421 return origin;
1422 }
1423 else if ( abs(maPdgcode) == 333 ) {
1424 origin = kGammaPhi;
1425 return origin;
1426 }
1427 else if ( abs(maPdgcode) == 331 ) {
1428 origin = kGammaEtaPrime;
1429 return origin;
1430 }
1431 else if ( abs(maPdgcode) == 113 ) {
1432 origin = kGammaRho0;
1433 return origin;
1434 }
1435 else origin = kElse;
1436 //origin = kGamma; // finer category above
e3fc062d 1437 return origin;
c2690925 1438
e3fc062d 1439 } // end of if
1440 else if ( abs(maPdgcode) == 111 ) {
1441 origin = kPi0;
1442 return origin;
1443 } // end of if
1444 else if ( abs(maPdgcode) == 221 ) {
1445 origin = kEta;
1446 return origin;
1447 } // end of if
1448 else if ( abs(maPdgcode) == 223 ) {
1449 origin = kOmega;
1450 return origin;
1451 } // end of if
1452 else if ( abs(maPdgcode) == 333 ) {
1453 origin = kPhi;
1454 return origin;
1455 } // end of if
1456 else if ( abs(maPdgcode) == 331 ) {
1457 origin = kEtaPrime;
1458 return origin;
1459 } // end of if
1460 else if ( abs(maPdgcode) == 113 ) {
1461 origin = kRho0;
1462 return origin;
1463 } // end of if
1464 else origin = kElse;
1465
1466 return origin;
70da6c5a 1467}
9250ffbf 1468//__________________________________________
8c1c76e9 1469Double_t AliHFEmcQA::GetWeightFactor(AliMCParticle *mctrack, const Int_t iBgLevel){
9250ffbf 1470 //
8c1c76e9 1471 // 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 1472 //
1473 AliMCParticle *mctrackmother = NULL;
1474 Double_t weightElecBg = 0.;
1475 Double_t mesonPt = 0.;
1476 Double_t bgcategory = 0.;
1477 Int_t mArr = -1;
1478 Int_t mesonID = GetElecSource(mctrack->Particle());
1479 if(mesonID==kGammaPi0 || mesonID==kPi0) mArr=0; //pion
1480 else if(mesonID==kGammaEta || mesonID==kEta) mArr=1; //eta
1481 else if(mesonID==kGammaOmega || mesonID==kOmega) mArr=2; //omega
1482 else if(mesonID==kGammaPhi || mesonID==kPhi) mArr=3; //phi
1483 else if(mesonID==kGammaEtaPrime || mesonID==kEtaPrime) mArr=4; //etaprime
1484 else if(mesonID==kGammaRho0 || mesonID==kRho0) mArr=5; //rho
1485
1486 if(!(mArr<0)){
1487 if(mesonID>=kGammaPi0) { // conversion electron, be careful with the enum odering
1488 Int_t glabel=TMath::Abs(mctrack->GetMother()); // gamma label
1489 if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
1490 glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's label
1491 if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
1492 mesonPt = mctrackmother->Pt(); //meson pt
1493 bgcategory = 1.;
1494 }
1495 }
1496 }
1497 else{ // nonHFE except for the conversion electron
1498 Int_t glabel=TMath::Abs(mctrack->GetMother());
1499 if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
1500 mesonPt=mctrackmother->Pt(); //meson pt
1501 bgcategory = -1.;
1502 }
1503 }
8c1c76e9 1504
1505 if(fIsPbPb){
1506 if(fCentrality < 0)return 0.;
1507 weightElecBg=fElecBackgroundFactor[iBgLevel][fCentrality][mArr][kBgPtBins-1];
1508 for(int ii=0; ii<kBgPtBins; ii++){
1509 if((mesonPt > fBinLimit[ii]) && (mesonPt < fBinLimit[ii+1])){
1510 weightElecBg = fElecBackgroundFactor[iBgLevel][fCentrality][mArr][ii];
1511 break;
1512 }
1513 }
9250ffbf 1514 }
8c1c76e9 1515 else{
1516 weightElecBg=fElecBackgroundFactor[iBgLevel][0][mArr][kBgPtBins-1];
9250ffbf 1517 for(int ii=0; ii<kBgPtBins; ii++){
8c1c76e9 1518 if((mesonPt > fBinLimit[ii]) && (mesonPt < fBinLimit[ii+1])){
1519 weightElecBg = fElecBackgroundFactor[iBgLevel][0][mArr][ii];
1520 break;
1521 }
9250ffbf 1522 }
8c1c76e9 1523 }
9250ffbf 1524 }
9250ffbf 1525 return bgcategory*weightElecBg;
1526}
1527
70da6c5a 1528//__________________________________________
1529void AliHFEmcQA::AliHists::FillList(TList *l) const {
1530 //
1531 // Fill Histos into a list for output
1532 //
1533 if(fPdgCode) l->Add(fPdgCode);
1534 if(fPt) l->Add(fPt);
1535 if(fY) l->Add(fY);
1536 if(fEta) l->Add(fEta);
1537}
1538
1539//__________________________________________
1540void AliHFEmcQA::AliHistsComm::FillList(TList *l) const {
1541 //
1542 // Fill Histos into a list for output
1543 //
1544 if(fNq) l->Add(fNq);
1545 if(fProcessID) l->Add(fProcessID);
1546 if(fePtRatio) l->Add(fePtRatio);
ccc37cdc 1547 if(fPtCorr) l->Add(fPtCorr);
c2690925 1548 if(fPtCorrDp) l->Add(fPtCorrDp);
1549 if(fPtCorrD0) l->Add(fPtCorrD0);
1550 if(fPtCorrDrest) l->Add(fPtCorrDrest);
70da6c5a 1551 if(fDePtRatio) l->Add(fDePtRatio);
1552 if(feDistance) l->Add(feDistance);
1553 if(fDeDistance) l->Add(fDeDistance);
1554}
1555