]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG3/hfe/AliHFEmcQA.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Major update of the HFE package (comments inside the code
[u/mrichter/AliRoot.git] / PWG3 / 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>
259c3296 40
41#include "AliHFEmcQA.h"
70da6c5a 42#include "AliHFEtools.h"
259c3296 43
44const Int_t AliHFEmcQA::fgkGluon=21;
45const Int_t AliHFEmcQA::fgkMaxGener=10;
dbe3abbe 46const Int_t AliHFEmcQA::fgkMaxIter=100;
47const Int_t AliHFEmcQA::fgkqType = 7; // number of species waiting for QA done
259c3296 48
49ClassImp(AliHFEmcQA)
50
51//_______________________________________________________________________________________________
52AliHFEmcQA::AliHFEmcQA() :
faee3b18 53 fMCEvent(NULL)
70da6c5a 54 ,fMCHeader(NULL)
55 ,fMCArray(NULL)
56 ,fQAhistos(NULL)
259c3296 57 ,fNparents(0)
58{
59 // Default constructor
60
259c3296 61}
62
63//_______________________________________________________________________________________________
64AliHFEmcQA::AliHFEmcQA(const AliHFEmcQA&p):
65 TObject(p)
faee3b18 66 ,fMCEvent(NULL)
70da6c5a 67 ,fMCHeader(NULL)
68 ,fMCArray(NULL)
69 ,fQAhistos(p.fQAhistos)
259c3296 70 ,fNparents(p.fNparents)
71{
72 // Copy constructor
73}
74
75//_______________________________________________________________________________________________
76AliHFEmcQA&
77AliHFEmcQA::operator=(const AliHFEmcQA &)
78{
79 // Assignment operator
80
81 AliInfo("Not yet implemented.");
82 return *this;
83}
84
85//_______________________________________________________________________________________________
86AliHFEmcQA::~AliHFEmcQA()
87{
88 // Destructor
89
90 AliInfo("Analysis Done.");
91}
92
93//_______________________________________________________________________________________________
50685501 94void AliHFEmcQA::PostAnalyze() const
259c3296 95{
e3fc062d 96 //
97 // Post analysis
98 //
259c3296 99}
100
e3fc062d 101//__________________________________________
102void AliHFEmcQA::CreatDefaultHistograms(TList * const qaList)
103{
104 //
105 // make default histograms
106 //
107
108 if(!qaList) return;
109
110 fQAhistos = qaList;
111 fQAhistos->SetName("MCqa");
112
113 CreateHistograms(AliHFEmcQA::kCharm,0,"mcqa_"); // create histograms for charm
114 CreateHistograms(AliHFEmcQA::kBeauty,0,"mcqa_"); // create histograms for beauty
115 CreateHistograms(AliHFEmcQA::kOthers,0,"mcqa_"); // create histograms for beauty
116 CreateHistograms(AliHFEmcQA::kCharm,1,"mcqa_barrel_"); // create histograms for charm
117 CreateHistograms(AliHFEmcQA::kBeauty,1,"mcqa_barrel_"); // create histograms for beauty
118 CreateHistograms(AliHFEmcQA::kOthers,1,"mcqa_barrel_"); // create histograms for beauty
119 CreateHistograms(AliHFEmcQA::kCharm,2,"mcqa_unitY_"); // create histograms for charm
120 CreateHistograms(AliHFEmcQA::kBeauty,2,"mcqa_unitY_"); // create histograms for beauty
121 CreateHistograms(AliHFEmcQA::kOthers,2,"mcqa_unitY_"); // create histograms for beauty
122 CreateHistograms(AliHFEmcQA::kCharm,3,"mcqa_reccut_"); // create histograms for charm
123 CreateHistograms(AliHFEmcQA::kBeauty,3,"mcqa_reccut_"); // create histograms for beauty
124 CreateHistograms(AliHFEmcQA::kOthers,3,"mcqa_reccut_"); // create histograms for beauty
125 CreateHistograms(AliHFEmcQA::kCharm,4,"mcqa_recpidcut_"); // create histograms for charm
126 CreateHistograms(AliHFEmcQA::kBeauty,4,"mcqa_recpidcut_"); // create histograms for beauty
127 CreateHistograms(AliHFEmcQA::kOthers,4,"mcqa_recpidcut_"); // create histograms for beauty
128 CreateHistograms(AliHFEmcQA::kCharm,5,"mcqa_secvtxcut_"); // create histograms for charm
129 CreateHistograms(AliHFEmcQA::kBeauty,5,"mcqa_secvtxcut_"); // create histograms for beauty
130 CreateHistograms(AliHFEmcQA::kOthers,5,"mcqa_secvtxcut_"); // create histograms for beauty
131
132}
133
259c3296 134//__________________________________________
dbe3abbe 135void AliHFEmcQA::CreateHistograms(const Int_t kquark, Int_t icut, TString hnopt)
259c3296 136{
137 // create histograms
138
faee3b18 139 if (!(kquark == kCharm || kquark == kBeauty || kquark == kOthers)) {
259c3296 140 AliDebug(1, "This task is only for heavy quark QA, return\n");
141 return;
142 }
dbe3abbe 143 Int_t iq = kquark - kCharm;
259c3296 144
145 TString kqTypeLabel[fgkqType];
dbe3abbe 146 if (kquark == kCharm){
147 kqTypeLabel[kQuark]="c";
148 kqTypeLabel[kantiQuark]="cbar";
149 kqTypeLabel[kHadron]="cHadron";
150 kqTypeLabel[keHadron]="ceHadron";
151 kqTypeLabel[kDeHadron]="nullHadron";
152 kqTypeLabel[kElectron]="ce";
153 kqTypeLabel[kElectron2nd]="nulle";
154 } else if (kquark == kBeauty){
155 kqTypeLabel[kQuark]="b";
156 kqTypeLabel[kantiQuark]="bbar";
157 kqTypeLabel[kHadron]="bHadron";
158 kqTypeLabel[keHadron]="beHadron";
159 kqTypeLabel[kDeHadron]="bDeHadron";
160 kqTypeLabel[kElectron]="be";
161 kqTypeLabel[kElectron2nd]="bce";
faee3b18 162 } else if (kquark == kOthers){
163 kqTypeLabel[kGamma-4]="gammae";
164 kqTypeLabel[kPi0-4]="pi0e";
165 kqTypeLabel[kElse-4]="elsee";
166 kqTypeLabel[kMisID-4]="miside";
259c3296 167 }
3a72645a 168
169 const Double_t kPtbound[2] = {0.1, 20.}; //bin taken for considering inclusive e analysis binning
faee3b18 170 Int_t iBin[1];
3a72645a 171 iBin[0] = 44; // bins in pt
faee3b18 172 Double_t* binEdges[1];
173 binEdges[0] = AliHFEtools::MakeLogarithmicBinning(iBin[0], kPtbound[0], kPtbound[1]);
3a72645a 174
faee3b18 175
259c3296 176 TString hname;
faee3b18 177 if(kquark == kOthers){
178 for (Int_t iqType = 0; iqType < 4; iqType++ ){
179 hname = hnopt+"Pt_"+kqTypeLabel[iqType];
3a72645a 180 fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",iBin[0],binEdges[0]);
181 //fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",500,0,50);
faee3b18 182 hname = hnopt+"Y_"+kqTypeLabel[iqType];
183 fHist[iq][iqType][icut].fY = new TH1F(hname,hname,150,-7.5,7.5);
184 hname = hnopt+"Eta_"+kqTypeLabel[iqType];
185 fHist[iq][iqType][icut].fEta = new TH1F(hname,hname,150,-7.5,7.5);
186 // Fill List
187 if(fQAhistos) fHist[iq][iqType][icut].FillList(fQAhistos);
188 }
189 return;
190 }
259c3296 191 for (Int_t iqType = 0; iqType < fgkqType; iqType++ ){
dbe3abbe 192 if (iqType < keHadron && icut > 0) continue; // don't duplicate histogram for quark and hadron
259c3296 193 hname = hnopt+"PdgCode_"+kqTypeLabel[iqType];
dbe3abbe 194 fHist[iq][iqType][icut].fPdgCode = new TH1F(hname,hname,20001,-10000.5,10000.5);
259c3296 195 hname = hnopt+"Pt_"+kqTypeLabel[iqType];
3a72645a 196 fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",iBin[0],binEdges[0]);
197 //fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",500,0,50);
259c3296 198 hname = hnopt+"Y_"+kqTypeLabel[iqType];
dbe3abbe 199 fHist[iq][iqType][icut].fY = new TH1F(hname,hname,150,-7.5,7.5);
259c3296 200 hname = hnopt+"Eta_"+kqTypeLabel[iqType];
dbe3abbe 201 fHist[iq][iqType][icut].fEta = new TH1F(hname,hname,150,-7.5,7.5);
70da6c5a 202 // Fill List
203 if(fQAhistos) fHist[iq][iqType][icut].FillList(fQAhistos);
259c3296 204 }
205
dbe3abbe 206 if (icut == 0){
207 hname = hnopt+"Nq_"+kqTypeLabel[kQuark];
208 fHistComm[iq][icut].fNq = new TH1F(hname,hname,10,-0.5,9.5);
209 hname = hnopt+"ProcessID_"+kqTypeLabel[kQuark];
210 fHistComm[iq][icut].fProcessID = new TH1F(hname,hname,21,-10.5,10.5);
211 }
212 hname = hnopt+"ePtRatio_"+kqTypeLabel[kQuark];
213 fHistComm[iq][icut].fePtRatio = new TH2F(hname,hname+";p_{T} (GeV/c);momentum fraction",100,0,30,100,0,1);
214 hname = hnopt+"DePtRatio_"+kqTypeLabel[kQuark];
215 fHistComm[iq][icut].fDePtRatio = new TH2F(hname,hname+";p_{T} (GeV/c);momentum fraction",100,0,30,100,0,1);
216 hname = hnopt+"eDistance_"+kqTypeLabel[kQuark];
217 fHistComm[iq][icut].feDistance= new TH2F(hname,hname+";p_{T} (GeV/c);distance (cm)",100,0,30,200,0,2);
218 hname = hnopt+"DeDistance_"+kqTypeLabel[kQuark];
219 fHistComm[iq][icut].fDeDistance= new TH2F(hname,hname+";p_{T} (GeV/c);distance (cm)",100,0,30,200,0,2);
70da6c5a 220 if(fQAhistos) fHistComm[iq][icut].FillList(fQAhistos);
259c3296 221}
222
223//__________________________________________
224void AliHFEmcQA::Init()
225{
226 // called at begining every event
227
228 for (Int_t i=0; i<2; i++){
229 fIsHeavy[i] = 0;
230 }
231
232 fNparents = 7;
233
234 fParentSelect[0][0] = 411;
235 fParentSelect[0][1] = 421;
236 fParentSelect[0][2] = 431;
237 fParentSelect[0][3] = 4122;
238 fParentSelect[0][4] = 4132;
239 fParentSelect[0][5] = 4232;
240 fParentSelect[0][6] = 4332;
241
242 fParentSelect[1][0] = 511;
243 fParentSelect[1][1] = 521;
244 fParentSelect[1][2] = 531;
245 fParentSelect[1][3] = 5122;
246 fParentSelect[1][4] = 5132;
247 fParentSelect[1][5] = 5232;
248 fParentSelect[1][6] = 5332;
249
250}
251
252//__________________________________________
722347d8 253void AliHFEmcQA::GetQuarkKine(TParticle *part, Int_t iTrack, const Int_t kquark)
259c3296 254{
255 // get heavy quark kinematics
256
dbe3abbe 257 if (kquark != kCharm && kquark != kBeauty) {
259c3296 258 AliDebug(1, "This task is only for heavy quark QA, return\n");
259 return;
260 }
dbe3abbe 261 Int_t iq = kquark - kCharm;
259c3296 262
faee3b18 263 if (iTrack < 0 || !part) {
264 AliDebug(1, "Stack label is negative or no mcparticle, return\n");
259c3296 265 return;
266 }
267
faee3b18 268 AliMCParticle *mctrack = NULL;
259c3296 269 Int_t partPdgcode = TMath::Abs(part->GetPdgCode());
270
271 // select heavy hadron or not fragmented heavy quark
272 if ( int(partPdgcode/100.)==kquark || int(partPdgcode/1000.)==kquark || (partPdgcode==kquark && (part->GetNDaughters()==0 && iTrack>5)) ){
273
274 TParticle *partMother;
275 Int_t iLabel;
276
277 if (partPdgcode == kquark){ // in case of not fragmented heavy quark
278 partMother = part;
279 iLabel = iTrack;
280 } else{ // in case of heavy hadron, start to search for mother heavy parton
281 iLabel = part->GetFirstMother();
faee3b18 282 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return;
283 if (iLabel>-1) { partMother = mctrack->Particle(); }
259c3296 284 else {
285 AliDebug(1, "Stack label is negative, return\n");
286 return;
287 }
288 }
289
290 // heavy parton selection as a mother of heavy hadron
291 // if the heavy particle comes from string which is denoted as particle status 12|12|12...12|11,[PYTHIA p.60]
292 // in this case, the mother of heavy particle can be one of the fragmented parton of the string
293 // should I make a condition that partMother should be quark or diquark? -> not necessary
294 if ( abs(partMother->GetPdgCode()) == kquark || (partMother->GetStatusCode() == 11) ){
295 //if ( abs(partMother->GetPdgCode()) == kquark || (partMother->GetStatusCode() == 11 || partMother->GetStatusCode() == 12) ){
296
297 if ( abs(partMother->GetPdgCode()) != kquark ){
298 // search fragmented partons in the same string
dbe3abbe 299 Bool_t isSameString = kTRUE;
259c3296 300 for (Int_t i=1; i<fgkMaxIter; i++){
301 iLabel = iLabel - 1;
faee3b18 302 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return;
303 if (iLabel>-1) { partMother = mctrack->Particle(); }
259c3296 304 else {
305 AliDebug(1, "Stack label is negative, return\n");
306 return;
307 }
308 if ( abs(partMother->GetPdgCode()) == kquark ) break;
dbe3abbe 309 if ( partMother->GetStatusCode() != 12 ) isSameString = kFALSE;
310 if (!isSameString) return;
259c3296 311 }
312 }
313 AliDebug(1, "Can not find heavy parton of this heavy hadron in the string, return\n");
314 if (abs(partMother->GetPdgCode()) != kquark) return;
315
316 if (fIsHeavy[iq] >= 50) return;
317 fHeavyQuark[fIsHeavy[iq]] = partMother;
318 fIsHeavy[iq]++;
319
320 // fill kinematics for heavy parton
321 if (partMother->GetPdgCode() > 0) { // quark
dbe3abbe 322 fHist[iq][kQuark][0].fPdgCode->Fill(partMother->GetPdgCode());
323 fHist[iq][kQuark][0].fPt->Fill(partMother->Pt());
70da6c5a 324 fHist[iq][kQuark][0].fY->Fill(AliHFEtools::GetRapidity(partMother));
dbe3abbe 325 fHist[iq][kQuark][0].fEta->Fill(partMother->Eta());
259c3296 326 } else{ // antiquark
dbe3abbe 327 fHist[iq][kantiQuark][0].fPdgCode->Fill(partMother->GetPdgCode());
328 fHist[iq][kantiQuark][0].fPt->Fill(partMother->Pt());
70da6c5a 329 fHist[iq][kantiQuark][0].fY->Fill(AliHFEtools::GetRapidity(partMother));
dbe3abbe 330 fHist[iq][kantiQuark][0].fEta->Fill(partMother->Eta());
259c3296 331 }
332
333 } // end of heavy parton slection loop
334
335 } // end of heavy hadron or quark selection
336
337}
338
339//__________________________________________
340void AliHFEmcQA::EndOfEventAna(const Int_t kquark)
341{
342 // end of event analysis
343
dbe3abbe 344 if (kquark != kCharm && kquark != kBeauty) {
259c3296 345 AliDebug(1, "This task is only for heavy quark QA, return\n");
346 return;
347 }
dbe3abbe 348 Int_t iq = kquark - kCharm;
259c3296 349
350
351 // # of heavy quark per event
352 AliDebug(1,Form("Number of heavy quark in this event = %d \n",fIsHeavy[iq]));
dbe3abbe 353 fHistComm[iq][0].fNq->Fill(fIsHeavy[iq]);
259c3296 354
355 Int_t motherID[fgkMaxGener];
356 Int_t motherType[fgkMaxGener];
357 Int_t motherLabel[fgkMaxGener];
358 Int_t ancestorPdg[fgkMaxGener];
359 Int_t ancestorLabel[fgkMaxGener];
360
361 for (Int_t i = 0; i < fgkMaxGener; i++){ // initialization
362 motherID[i] = 0;
363 motherType[i] = 0;
364 motherLabel[i] = 0;
365 ancestorPdg[i] = 0;
366 ancestorLabel[i] = 0;
367 }
368
3a72645a 369
259c3296 370 // check history of found heavy quarks
371 for (Int_t i = 0; i < fIsHeavy[iq]; i++){
372
3a72645a 373 if(!fHeavyQuark[i]) return;
374
259c3296 375 ancestorLabel[0] = i;
376 ancestorPdg[0] = fHeavyQuark[i]->GetPdgCode();
377 ancestorLabel[1] = fHeavyQuark[i]->GetFirstMother();
378
379 AliDebug(1,Form("pdg code= %d\n",ancestorPdg[0]));
380 AliDebug(1,Form("ancestor label= %d\n",ancestorLabel[1]));
381
382 Int_t ig = 1;
383 while (ancestorLabel[ig] != -1){
384 // in case there is mother, get mother's pdg code and grandmother's label
385 IdentifyMother(ancestorLabel[ig], ancestorPdg[ig], ancestorLabel[ig+1]);
386 // if mother is still heavy, find again mother's ancestor
387 if (ancestorPdg[ig-1] == ancestorPdg[ig]) {
388 ig++;
389 continue; // if it is from same heavy
390 }
391 // if the heavy's mother is not heavy, check the mother's label to know if it comes from inital or final parton shower
392 if (IsFromInitialShower(ancestorLabel[ig],motherID[i],motherType[i],motherLabel[i])) break;
393 if (IsFromFinalParton(ancestorLabel[ig],motherID[i],motherType[i],motherLabel[i])) break;
394 // if it is not the above case, something is strange
dbe3abbe 395 ReportStrangeness(motherID[i],motherType[i],motherLabel[i]);
259c3296 396 break;
397 }
398 if (ancestorLabel[ig] == -1){ // from hard scattering
399 HardScattering(kquark, motherID[i],motherType[i], motherLabel[i]);
400 }
401
402 } // end of found heavy quark loop
403
404
405 // check process type
406 Int_t processID = 0;
407 for (Int_t i = 0; i < fIsHeavy[iq]; i++){
408 AliDebug(1,Form("Mother ID= %d type= %d label= %d\n",motherID[i],motherType[i],motherLabel[i]));
409 }
410
411
412 Int_t nheavypair = Int_t(fIsHeavy[iq]/2.);
413 for (Int_t ipair = 0; ipair < nheavypair; ipair++){
414
415 Int_t id1 = ipair*2;
416 Int_t id2 = ipair*2 + 1;
417
418 if (motherType[id1] == 2 && motherType[id2] == 2){
dbe3abbe 419 if (motherLabel[id1] == motherLabel[id2]) processID = kGluonSplitting; // gluon spliting
259c3296 420 else processID = -9;
421 }
422 else if (motherType[id1] == -1 && motherType[id2] == -1) {
423 if (motherLabel[id1] == -1 && motherLabel[id2] == -1) {
dbe3abbe 424 if (motherID[id1] == fgkGluon) processID = kPairCreationFromg; // gluon fusion
425 else processID = kPairCreationFromq; // q-qbar pair creation
259c3296 426 }
427 else processID = -8;
428 }
429 else if (motherType[id1] == -1 || motherType[id2] == -1) {
430 if ((motherLabel[id1] == -1 || motherLabel[id2] == -1) && (motherLabel[id1]*motherLabel[id2] == -2 || motherLabel[id1]*motherLabel[id2] == -3)) {
dbe3abbe 431 if(motherID[id1]*motherID[id2] == kquark*fgkGluon) processID = kFlavourExitation; // flavour exitation
432 else processID = kLightQuarkShower;
259c3296 433 }
434 else processID = -7;
435 }
436 else if (motherType[id1] == -2 || motherType[id2] == -2) {
dbe3abbe 437 if (motherLabel[id1] == motherLabel[id2]) processID = kInitialPartonShower; // initial parton shower
259c3296 438 else processID = -6;
439
440 }
441 else processID = -5;
442
443 if (nheavypair >1) AliDebug(1,Form("Multi pair found : process ID = %d\n",processID));
dbe3abbe 444 else fHistComm[iq][0].fProcessID->Fill(processID);
259c3296 445 AliDebug(1,Form("Process ID = %d\n",processID));
446 } // end of # heavy quark pair loop
447
448}
449
450//__________________________________________
722347d8 451void AliHFEmcQA::GetHadronKine(TParticle* mcpart, const Int_t kquark)
dbe3abbe 452{
453 // decay electron kinematics
454
455 if (kquark != kCharm && kquark != kBeauty) {
456 AliDebug(1, "This task is only for heavy quark QA, return\n");
457 return;
458 }
459 Int_t iq = kquark - kCharm;
460
faee3b18 461 if(!mcpart){
462 AliDebug(1, "no mc particle, return\n");
463 return;
464 }
dbe3abbe 465
466 Int_t iLabel = mcpart->GetFirstMother();
467 if (iLabel<0){
468 AliDebug(1, "Stack label is negative, return\n");
469 return;
470 }
471
472 TParticle *partCopy = mcpart;
473 Int_t pdgcode = mcpart->GetPdgCode();
474 Int_t pdgcodeCopy = pdgcode;
475
faee3b18 476 AliMCParticle *mctrack = NULL;
477
dbe3abbe 478 // if the mother is charmed hadron
75d81601 479 Bool_t isDirectCharm = kFALSE;
dbe3abbe 480 if ( int(abs(pdgcode)/100.) == kCharm || int(abs(pdgcode)/1000.) == kCharm ) {
481
482 // iterate until you find B hadron as a mother or become top ancester
483 for (Int_t i=1; i<fgkMaxIter; i++){
484
485 Int_t jLabel = mcpart->GetFirstMother();
486 if (jLabel == -1){
75d81601 487 isDirectCharm = kTRUE;
dbe3abbe 488 break; // if there is no ancester
489 }
490 if (jLabel < 0){ // safety protection
491 AliDebug(1, "Stack label is negative, return\n");
492 return;
493 }
494 // if there is an ancester
faee3b18 495 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return;
496 TParticle* mother = mctrack->Particle();
dbe3abbe 497 Int_t motherPDG = mother->GetPdgCode();
498
499 for (Int_t j=0; j<fNparents; j++){
500 if (abs(motherPDG)==fParentSelect[1][j]) return; // return if this hadron is originated from b
501 }
502
503 mcpart = mother;
504 } // end of iteration
505 } // end of if
75d81601 506 if((isDirectCharm == kTRUE && kquark == kCharm) || kquark == kBeauty) {
dbe3abbe 507 for (Int_t i=0; i<fNparents; i++){
508 if (abs(pdgcodeCopy)==fParentSelect[iq][i]){
509
510 // fill hadron kinematics
511 fHist[iq][kHadron][0].fPdgCode->Fill(pdgcodeCopy);
512 fHist[iq][kHadron][0].fPt->Fill(partCopy->Pt());
70da6c5a 513 fHist[iq][kHadron][0].fY->Fill(AliHFEtools::GetRapidity(partCopy));
dbe3abbe 514 fHist[iq][kHadron][0].fEta->Fill(partCopy->Eta());
515
516 }
517 }
518 } // end of if
519}
520
521//__________________________________________
722347d8 522void AliHFEmcQA::GetDecayedKine(TParticle* mcpart, const Int_t kquark, Int_t kdecayed, Int_t icut)
259c3296 523{
524 // decay electron kinematics
525
faee3b18 526 if (!(kquark == kCharm || kquark == kBeauty || kquark == kOthers)){
259c3296 527 AliDebug(1, "This task is only for heavy quark QA, return\n");
528 return;
529 }
dbe3abbe 530 Int_t iq = kquark - kCharm;
50685501 531 Bool_t isFinalOpenCharm = kFALSE;
259c3296 532
faee3b18 533 if(!mcpart){
534 AliDebug(1, "no mcparticle, return\n");
535 return;
259c3296 536 }
537
faee3b18 538 if(kquark==kOthers){
539 Int_t esource = -1;
540 if ( abs(mcpart->GetPdgCode()) != kdecayed ) esource = kMisID-4;
541 else esource =GetSource(mcpart)-4; // return for the cases kGamma=4, kPi0=5, kElse=6
542 if(esource==0|| esource==1 || esource==2 || esource==3){
543 fHist[iq][esource][icut].fPt->Fill(mcpart->Pt());
544 fHist[iq][esource][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart));
545 fHist[iq][esource][icut].fEta->Fill(mcpart->Eta());
546 return;
547 }
548 else {
549 AliDebug(1, "e source is out of defined ranges, return\n");
550 return;
551 }
552 }
259c3296 553
554 if ( abs(mcpart->GetPdgCode()) != kdecayed ) return;
259c3296 555
556 Int_t iLabel = mcpart->GetFirstMother();
557 if (iLabel<0){
558 AliDebug(1, "Stack label is negative, return\n");
559 return;
560 }
561
faee3b18 562 AliMCParticle *mctrack = NULL;
563 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return;
564 TParticle *partMother = mctrack->Particle();
dbe3abbe 565 TParticle *partMotherCopy = partMother;
259c3296 566 Int_t maPdgcode = partMother->GetPdgCode();
dbe3abbe 567 Int_t maPdgcodeCopy = maPdgcode;
259c3296 568
9bcfd1ab 569 // get mc primary vertex
faee3b18 570 /*
9bcfd1ab 571 TArrayF mcPrimVtx;
572 if(fMCHeader) fMCHeader->PrimaryVertex(mcPrimVtx);
573
259c3296 574 // get electron production vertex
575 TLorentzVector ePoint;
576 mcpart->ProductionVertex(ePoint);
577
9bcfd1ab 578 // calculated production vertex to primary vertex (in xy plane)
579 Float_t decayLxy = TMath::Sqrt((mcPrimVtx[0]-ePoint[0])*(mcPrimVtx[0]-ePoint[0])+(mcPrimVtx[1]-ePoint[1])*(mcPrimVtx[1]-ePoint[1]));
faee3b18 580 */
581 Float_t decayLxy = 0;
9bcfd1ab 582
259c3296 583 // if the mother is charmed hadron
dbe3abbe 584 Bool_t isMotherDirectCharm = kFALSE;
585 if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
259c3296 586
0792aa82 587 for (Int_t i=0; i<fNparents; i++){
588 if (abs(maPdgcode)==fParentSelect[0][i]){
50685501 589 isFinalOpenCharm = kTRUE;
0792aa82 590 }
591 }
50685501 592 if (!isFinalOpenCharm) return ;
0792aa82 593
259c3296 594 // iterate until you find B hadron as a mother or become top ancester
595 for (Int_t i=1; i<fgkMaxIter; i++){
596
597 Int_t jLabel = partMother->GetFirstMother();
598 if (jLabel == -1){
dbe3abbe 599 isMotherDirectCharm = kTRUE;
259c3296 600 break; // if there is no ancester
601 }
602 if (jLabel < 0){ // safety protection
603 AliDebug(1, "Stack label is negative, return\n");
604 return;
605 }
606
607 // if there is an ancester
faee3b18 608 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return;
609 TParticle* grandMa = mctrack->Particle();
259c3296 610 Int_t grandMaPDG = grandMa->GetPdgCode();
611
fc0de2a0 612 for (Int_t j=0; j<fNparents; j++){
613 if (abs(grandMaPDG)==fParentSelect[1][j]){
259c3296 614
dbe3abbe 615 if (kquark == kCharm) return;
259c3296 616 // fill electron kinematics
dbe3abbe 617 fHist[iq][kElectron2nd][icut].fPdgCode->Fill(mcpart->GetPdgCode());
618 fHist[iq][kElectron2nd][icut].fPt->Fill(mcpart->Pt());
70da6c5a 619 fHist[iq][kElectron2nd][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart));
dbe3abbe 620 fHist[iq][kElectron2nd][icut].fEta->Fill(mcpart->Eta());
259c3296 621
622 // fill mother hadron kinematics
dbe3abbe 623 fHist[iq][kDeHadron][icut].fPdgCode->Fill(grandMaPDG);
624 fHist[iq][kDeHadron][icut].fPt->Fill(grandMa->Pt());
70da6c5a 625 fHist[iq][kDeHadron][icut].fY->Fill(AliHFEtools::GetRapidity(grandMa));
dbe3abbe 626 fHist[iq][kDeHadron][icut].fEta->Fill(grandMa->Eta());
259c3296 627
628 // ratio between pT of electron and pT of mother B hadron
dbe3abbe 629 if(grandMa->Pt()) fHistComm[iq][icut].fDePtRatio->Fill(grandMa->Pt(),mcpart->Pt()/grandMa->Pt());
259c3296 630
9bcfd1ab 631 // distance between electron production point and primary vertex
632 fHistComm[iq][icut].fDeDistance->Fill(grandMa->Pt(),decayLxy);
259c3296 633 return;
634 }
dc306130 635 }
259c3296 636
637 partMother = grandMa;
638 } // end of iteration
639 } // end of if
dbe3abbe 640 if((isMotherDirectCharm == kTRUE && kquark == kCharm) || kquark == kBeauty) {
259c3296 641 for (Int_t i=0; i<fNparents; i++){
642 if (abs(maPdgcodeCopy)==fParentSelect[iq][i]){
643
644 // fill electron kinematics
dbe3abbe 645 fHist[iq][kElectron][icut].fPdgCode->Fill(mcpart->GetPdgCode());
646 fHist[iq][kElectron][icut].fPt->Fill(mcpart->Pt());
70da6c5a 647 fHist[iq][kElectron][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart));
dbe3abbe 648 fHist[iq][kElectron][icut].fEta->Fill(mcpart->Eta());
259c3296 649
650 // fill mother hadron kinematics
dbe3abbe 651 fHist[iq][keHadron][icut].fPdgCode->Fill(maPdgcodeCopy);
652 fHist[iq][keHadron][icut].fPt->Fill(partMotherCopy->Pt());
70da6c5a 653 fHist[iq][keHadron][icut].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
dbe3abbe 654 fHist[iq][keHadron][icut].fEta->Fill(partMotherCopy->Eta());
259c3296 655
656 // ratio between pT of electron and pT of mother B hadron
dbe3abbe 657 if(partMotherCopy->Pt()) fHistComm[iq][icut].fePtRatio->Fill(partMotherCopy->Pt(),mcpart->Pt()/partMotherCopy->Pt());
259c3296 658
9bcfd1ab 659 // distance between electron production point and primary vertex
660 fHistComm[iq][icut].feDistance->Fill(partMotherCopy->Pt(),decayLxy);
259c3296 661 }
662 }
663 } // end of if
664}
665
0792aa82 666//____________________________________________________________________
667void AliHFEmcQA::GetDecayedKine(AliAODMCParticle *mcpart, const Int_t kquark, Int_t kdecayed, Int_t icut)
668{
669 // decay electron kinematics
670
671 if (kquark != kCharm && kquark != kBeauty) {
672 AliDebug(1, "This task is only for heavy quark QA, return\n");
673 return;
674 }
675
676 Int_t iq = kquark - kCharm;
50685501 677 Bool_t isFinalOpenCharm = kFALSE;
0792aa82 678
faee3b18 679 if(!mcpart){
680 AliDebug(1, "no mcparticle, return\n");
681 return;
682 }
683
0792aa82 684 if ( abs(mcpart->GetPdgCode()) != kdecayed ) return;
685
686 // mother
687 Int_t iLabel = mcpart->GetMother();
688 if (iLabel<0){
689 AliDebug(1, "Stack label is negative, return\n");
690 return;
691 }
692
693 AliAODMCParticle *partMother = (AliAODMCParticle*)fMCArray->At(iLabel);
694 AliAODMCParticle *partMotherCopy = partMother;
695 Int_t maPdgcode = partMother->GetPdgCode();
696 Int_t maPdgcodeCopy = maPdgcode;
697
698 Bool_t isMotherDirectCharm = kFALSE;
699 if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
700
701 for (Int_t i=0; i<fNparents; i++){
702 if (abs(maPdgcode)==fParentSelect[0][i]){
50685501 703 isFinalOpenCharm = kTRUE;
0792aa82 704 }
705 }
50685501 706 if (!isFinalOpenCharm) return;
0792aa82 707
708 for (Int_t i=1; i<fgkMaxIter; i++){
709
710 Int_t jLabel = partMother->GetMother();
711 if (jLabel == -1){
712 isMotherDirectCharm = kTRUE;
713 break; // if there is no ancester
714 }
715 if (jLabel < 0){ // safety protection
716 AliDebug(1, "Stack label is negative, return\n");
717 return;
718 }
719
720 // if there is an ancester
721 AliAODMCParticle* grandMa = (AliAODMCParticle*)fMCArray->At(jLabel);
722 Int_t grandMaPDG = grandMa->GetPdgCode();
723
724 for (Int_t j=0; j<fNparents; j++){
725 if (abs(grandMaPDG)==fParentSelect[1][j]){
726
727 if (kquark == kCharm) return;
728 // fill electron kinematics
729 fHist[iq][kElectron2nd][icut].fPdgCode->Fill(mcpart->GetPdgCode());
730 fHist[iq][kElectron2nd][icut].fPt->Fill(mcpart->Pt());
70da6c5a 731 fHist[iq][kElectron2nd][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart));
0792aa82 732 fHist[iq][kElectron2nd][icut].fEta->Fill(mcpart->Eta());
733
734 // fill mother hadron kinematics
735 fHist[iq][kDeHadron][icut].fPdgCode->Fill(grandMaPDG);
736 fHist[iq][kDeHadron][icut].fPt->Fill(grandMa->Pt());
70da6c5a 737 fHist[iq][kDeHadron][icut].fY->Fill(AliHFEtools::GetRapidity(grandMa));
0792aa82 738 fHist[iq][kDeHadron][icut].fEta->Fill(grandMa->Eta());
739
740 return;
741 }
742 }
743
744 partMother = grandMa;
745 } // end of iteration
746 } // end of if
747 if ((isMotherDirectCharm == kTRUE && kquark == kCharm) || kquark == kBeauty) {
748 for (Int_t i=0; i<fNparents; i++){
749 if (abs(maPdgcodeCopy)==fParentSelect[iq][i]){
750
751 // fill electron kinematics
752 fHist[iq][kElectron][icut].fPdgCode->Fill(mcpart->GetPdgCode());
753 fHist[iq][kElectron][icut].fPt->Fill(mcpart->Pt());
70da6c5a 754 fHist[iq][kElectron][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart));
0792aa82 755 fHist[iq][kElectron][icut].fEta->Fill(mcpart->Eta());
756
757 // fill mother hadron kinematics
758 fHist[iq][keHadron][icut].fPdgCode->Fill(maPdgcodeCopy);
759 fHist[iq][keHadron][icut].fPt->Fill(partMotherCopy->Pt());
70da6c5a 760 fHist[iq][keHadron][icut].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
0792aa82 761 fHist[iq][keHadron][icut].fEta->Fill(partMotherCopy->Eta());
762
763 }
764 }
765 } // end of if
766
767}
259c3296 768
769//__________________________________________
75d81601 770void AliHFEmcQA::IdentifyMother(Int_t motherlabel, Int_t &motherpdg, Int_t &grandmotherlabel)
259c3296 771{
dbe3abbe 772 // find mother pdg code and label
773
75d81601 774 if (motherlabel < 0) {
259c3296 775 AliDebug(1, "Stack label is negative, return\n");
776 return;
777 }
faee3b18 778 AliMCParticle *mctrack = NULL;
779 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(motherlabel))))) return;
780 TParticle *heavysMother = mctrack->Particle();
75d81601 781 motherpdg = heavysMother->GetPdgCode();
782 grandmotherlabel = heavysMother->GetFirstMother();
783 AliDebug(1,Form("ancestor pdg code= %d\n",motherpdg));
259c3296 784}
785
786//__________________________________________
259c3296 787void AliHFEmcQA::HardScattering(const Int_t kquark, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
788{
dbe3abbe 789 // mothertype -1 means this heavy quark coming from hard vertex
790
faee3b18 791 AliMCParticle *mctrack1 = NULL;
792 AliMCParticle *mctrack2 = NULL;
793 if(!(mctrack1 = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(4))))) return;
794 if(!(mctrack2 = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(5))))) return;
795 TParticle *afterinitialrad1 = mctrack1->Particle();
796 TParticle *afterinitialrad2 = mctrack2->Particle();
259c3296 797
798 motherlabel = -1;
799
800 if (abs(afterinitialrad1->GetPdgCode()) == fgkGluon && abs(afterinitialrad2->GetPdgCode()) == fgkGluon){
801 AliDebug(1,"heavy from gluon gluon pair creation!\n");
802 mothertype = -1;
803 motherID = fgkGluon;
804 }
805 else if (abs(afterinitialrad1->GetPdgCode()) == kquark || abs(afterinitialrad2->GetPdgCode()) == kquark){ // one from Q and the other from g
806 AliDebug(1,"heavy from flavor exitation!\n");
807 mothertype = -1;
808 motherID = kquark;
809 }
810 else if (abs(afterinitialrad1->GetPdgCode()) == abs(afterinitialrad2->GetPdgCode())){
811 AliDebug(1,"heavy from q-qbar pair creation!\n");
812 mothertype = -1;
813 motherID = 1;
814 }
815 else {
816 AliDebug(1,"something strange!\n");
817 mothertype = -999;
818 motherlabel = -999;
819 motherID = -999;
820 }
821}
822
823//__________________________________________
259c3296 824Bool_t AliHFEmcQA::IsFromInitialShower(Int_t inputmotherlabel, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
825{
dbe3abbe 826 // mothertype -2 means this heavy quark coming from initial state
827
faee3b18 828 AliMCParticle *mctrack = NULL;
259c3296 829 if (inputmotherlabel==2 || inputmotherlabel==3){ // mother exist before initial state radiation
faee3b18 830 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(inputmotherlabel))))) return kFALSE;
831 TParticle *heavysMother = mctrack->Particle();
259c3296 832 motherID = heavysMother->GetPdgCode();
833 mothertype = -2; // there is mother before initial state radiation
834 motherlabel = inputmotherlabel;
835 AliDebug(1,"initial parton shower! \n");
836
837 return kTRUE;
838 }
839
840 return kFALSE;
841}
842
843//__________________________________________
259c3296 844Bool_t AliHFEmcQA::IsFromFinalParton(Int_t inputmotherlabel, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
845{
dbe3abbe 846 // mothertype 2 means this heavy quark coming from final state
847
faee3b18 848 AliMCParticle *mctrack = NULL;
259c3296 849 if (inputmotherlabel > 5){ // mother exist after hard scattering
faee3b18 850 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(inputmotherlabel))))) return kFALSE;
851 TParticle *heavysMother = mctrack->Particle();
259c3296 852 motherID = heavysMother->GetPdgCode();
853 mothertype = 2; //
854 motherlabel = inputmotherlabel;
855 AliDebug(1,Form("heavy quark from %d after hard scattering! \n",motherID));
856
857 return kTRUE;
858 }
859 return kFALSE;
860}
861
862//__________________________________________
dbe3abbe 863void AliHFEmcQA::ReportStrangeness(Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
259c3296 864{
dbe3abbe 865 // mark strange behavior
866
259c3296 867 mothertype = -888;
868 motherlabel = -888;
869 motherID = -888;
870 AliDebug(1,"something strange!\n");
871}
872
0792aa82 873//__________________________________________
9bcfd1ab 874Int_t AliHFEmcQA::GetSource(AliAODMCParticle * const mcpart)
0792aa82 875{
9bcfd1ab 876 // decay particle's origin
0792aa82 877
9bcfd1ab 878 //if ( abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return -1;
0792aa82 879
880 Int_t origin = -1;
50685501 881 Bool_t isFinalOpenCharm = kFALSE;
0792aa82 882
faee3b18 883 if(!mcpart){
884 AliDebug(1, "Stack label is negative or no mcparticle, return\n");
885 return -1;
886 }
887
0792aa82 888 // mother
889 Int_t iLabel = mcpart->GetMother();
890 if (iLabel<0){
891 AliDebug(1, "Stack label is negative, return\n");
892 return -1;
893 }
894
895 AliAODMCParticle *partMother = (AliAODMCParticle*)fMCArray->At(iLabel);
896 Int_t maPdgcode = partMother->GetPdgCode();
897
898 // if the mother is charmed hadron
899 if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
900
901 for (Int_t i=0; i<fNparents; i++){
902 if (abs(maPdgcode)==fParentSelect[0][i]){
50685501 903 isFinalOpenCharm = kTRUE;
0792aa82 904 }
905 }
50685501 906 if (!isFinalOpenCharm) return -1;
0792aa82 907
908 // iterate until you find B hadron as a mother or become top ancester
909 for (Int_t i=1; i<fgkMaxIter; i++){
910
911 Int_t jLabel = partMother->GetMother();
912 if (jLabel == -1){
913 origin = kDirectCharm;
914 return origin;
915 }
916 if (jLabel < 0){ // safety protection
917 AliDebug(1, "Stack label is negative, return\n");
918 return -1;
919 }
920
921 // if there is an ancester
922 AliAODMCParticle* grandMa = (AliAODMCParticle*)fMCArray->At(jLabel);
923 Int_t grandMaPDG = grandMa->GetPdgCode();
924
70b5ea26 925 for (Int_t j=0; j<fNparents; j++){
926 if (abs(grandMaPDG)==fParentSelect[1][j]){
0792aa82 927 origin = kBeautyCharm;
928 return origin;
929 }
930 }
931
932 partMother = grandMa;
933 } // end of iteration
934 } // end of if
935 else if ( int(abs(maPdgcode)/100.) == kBeauty || int(abs(maPdgcode)/1000.) == kBeauty ) {
936 for (Int_t i=0; i<fNparents; i++){
937 if (abs(maPdgcode)==fParentSelect[1][i]){
938 origin = kDirectBeauty;
939 return origin;
940 }
941 }
942 } // end of if
943 else if ( abs(maPdgcode) == 22 ) {
944 origin = kGamma;
945 return origin;
946 } // end of if
947 else if ( abs(maPdgcode) == 111 ) {
948 origin = kPi0;
949 return origin;
950 } // end of if
951
952 return origin;
953}
954
955//__________________________________________
9bcfd1ab 956Int_t AliHFEmcQA::GetSource(TParticle * const mcpart)
0792aa82 957{
9bcfd1ab 958 // decay particle's origin
0792aa82 959
9bcfd1ab 960 //if ( abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return -1;
0792aa82 961
962 Int_t origin = -1;
50685501 963 Bool_t isFinalOpenCharm = kFALSE;
0792aa82 964
faee3b18 965 if(!mcpart){
966 AliDebug(1, "no mcparticle, return\n");
967 return -1;
968 }
969
0792aa82 970 Int_t iLabel = mcpart->GetFirstMother();
971 if (iLabel<0){
972 AliDebug(1, "Stack label is negative, return\n");
973 return -1;
974 }
975
faee3b18 976 AliMCParticle *mctrack = NULL;
977 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return -1;
978 TParticle *partMother = mctrack->Particle();
0792aa82 979 Int_t maPdgcode = partMother->GetPdgCode();
980
981 // if the mother is charmed hadron
982 if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
983
984 for (Int_t i=0; i<fNparents; i++){
985 if (abs(maPdgcode)==fParentSelect[0][i]){
50685501 986 isFinalOpenCharm = kTRUE;
0792aa82 987 }
988 }
50685501 989 if (!isFinalOpenCharm) return -1;
0792aa82 990
991 // iterate until you find B hadron as a mother or become top ancester
992 for (Int_t i=1; i<fgkMaxIter; i++){
993
994 Int_t jLabel = partMother->GetFirstMother();
995 if (jLabel == -1){
996 origin = kDirectCharm;
997 return origin;
998 }
999 if (jLabel < 0){ // safety protection
1000 AliDebug(1, "Stack label is negative, return\n");
1001 return -1;
1002 }
1003
1004 // if there is an ancester
faee3b18 1005 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return -1;
1006 TParticle* grandMa = mctrack->Particle();
0792aa82 1007 Int_t grandMaPDG = grandMa->GetPdgCode();
1008
70b5ea26 1009 for (Int_t j=0; j<fNparents; j++){
1010 if (abs(grandMaPDG)==fParentSelect[1][j]){
0792aa82 1011 origin = kBeautyCharm;
1012 return origin;
1013 }
1014 }
1015
1016 partMother = grandMa;
1017 } // end of iteration
1018 } // end of if
1019 else if ( int(abs(maPdgcode)/100.) == kBeauty || int(abs(maPdgcode)/1000.) == kBeauty ) {
1020 for (Int_t i=0; i<fNparents; i++){
1021 if (abs(maPdgcode)==fParentSelect[1][i]){
1022 origin = kDirectBeauty;
1023 return origin;
1024 }
1025 }
1026 } // end of if
1027 else if ( abs(maPdgcode) == 22 ) {
1028 origin = kGamma;
1029 return origin;
1030 } // end of if
1031 else if ( abs(maPdgcode) == 111 ) {
1032 origin = kPi0;
1033 return origin;
1034 } // end of if
1035 else origin = kElse;
1036
1037 return origin;
1038}
70da6c5a 1039
1040//__________________________________________
e3fc062d 1041Int_t AliHFEmcQA::GetElecSource(TParticle * const mcpart)
1042{
1043 // decay particle's origin
70da6c5a 1044
3a72645a 1045 if ( abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return kMisID;
e3fc062d 1046
1047 Int_t origin = -1;
1048 Bool_t isFinalOpenCharm = kFALSE;
1049
1050 if(!mcpart){
1051 AliDebug(1, "no mcparticle, return\n");
1052 return -1;
1053 }
1054
1055 Int_t iLabel = mcpart->GetFirstMother();
1056 if (iLabel<0){
1057 AliDebug(1, "Stack label is negative, return\n");
1058 return -1;
1059 }
1060
1061 AliMCParticle *mctrack = NULL;
1062 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return -1;
1063 TParticle *partMother = mctrack->Particle();
1064 Int_t maPdgcode = partMother->GetPdgCode();
1065
1066 // if the mother is charmed hadron
1067 if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
1068
1069 for (Int_t i=0; i<fNparents; i++){
1070 if (abs(maPdgcode)==fParentSelect[0][i]){
1071 isFinalOpenCharm = kTRUE;
1072 }
1073 }
1074 if (!isFinalOpenCharm) return -1;
1075
1076 // iterate until you find B hadron as a mother or become top ancester
1077 for (Int_t i=1; i<fgkMaxIter; i++){
1078
1079 Int_t jLabel = partMother->GetFirstMother();
1080 if (jLabel == -1){
1081 origin = kDirectCharm;
1082 return origin;
1083 }
1084 if (jLabel < 0){ // safety protection
1085 AliDebug(1, "Stack label is negative, return\n");
1086 return -1;
1087 }
1088
1089 // if there is an ancester
1090 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return -1;
1091 TParticle* grandMa = mctrack->Particle();
1092 Int_t grandMaPDG = grandMa->GetPdgCode();
1093
1094 for (Int_t j=0; j<fNparents; j++){
1095 if (abs(grandMaPDG)==fParentSelect[1][j]){
1096 origin = kBeautyCharm;
1097 return origin;
1098 }
1099 }
1100
1101 partMother = grandMa;
1102 } // end of iteration
1103 } // end of if
1104 else if ( int(abs(maPdgcode)/100.) == kBeauty || int(abs(maPdgcode)/1000.) == kBeauty ) {
1105 for (Int_t i=0; i<fNparents; i++){
1106 if (abs(maPdgcode)==fParentSelect[1][i]){
1107 origin = kDirectBeauty;
1108 return origin;
1109 }
1110 }
1111 } // end of if
1112 else if ( abs(maPdgcode) == 22 ) {
1113 origin = kGamma;
1114 return origin;
1115 } // end of if
1116 else if ( abs(maPdgcode) == 111 ) {
1117 origin = kPi0;
1118 return origin;
1119 } // end of if
1120 else if ( abs(maPdgcode) == 221 ) {
1121 origin = kEta;
1122 return origin;
1123 } // end of if
1124 else if ( abs(maPdgcode) == 223 ) {
1125 origin = kOmega;
1126 return origin;
1127 } // end of if
1128 else if ( abs(maPdgcode) == 333 ) {
1129 origin = kPhi;
1130 return origin;
1131 } // end of if
1132 else if ( abs(maPdgcode) == 331 ) {
1133 origin = kEtaPrime;
1134 return origin;
1135 } // end of if
1136 else if ( abs(maPdgcode) == 113 ) {
1137 origin = kRho0;
1138 return origin;
1139 } // end of if
1140 else origin = kElse;
1141
1142 return origin;
70da6c5a 1143}
1144
1145//__________________________________________
1146void AliHFEmcQA::AliHists::FillList(TList *l) const {
1147 //
1148 // Fill Histos into a list for output
1149 //
1150 if(fPdgCode) l->Add(fPdgCode);
1151 if(fPt) l->Add(fPt);
1152 if(fY) l->Add(fY);
1153 if(fEta) l->Add(fEta);
1154}
1155
1156//__________________________________________
1157void AliHFEmcQA::AliHistsComm::FillList(TList *l) const {
1158 //
1159 // Fill Histos into a list for output
1160 //
1161 if(fNq) l->Add(fNq);
1162 if(fProcessID) l->Add(fProcessID);
1163 if(fePtRatio) l->Add(fePtRatio);
1164 if(fDePtRatio) l->Add(fDePtRatio);
1165 if(feDistance) l->Add(feDistance);
1166 if(fDeDistance) l->Add(fDeDistance);
1167}
1168