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