]> 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
adds possiblilty to use filelists as imput to run on PBS farm
[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>
259c3296 37#include <AliStack.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() :
70da6c5a 53 fStack(NULL)
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)
70da6c5a 66 ,fStack(NULL)
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{
96}
97
98//__________________________________________
dbe3abbe 99void AliHFEmcQA::CreateHistograms(const Int_t kquark, Int_t icut, TString hnopt)
259c3296 100{
101 // create histograms
102
dbe3abbe 103 if (kquark != kCharm && kquark != kBeauty) {
259c3296 104 AliDebug(1, "This task is only for heavy quark QA, return\n");
105 return;
106 }
dbe3abbe 107 Int_t iq = kquark - kCharm;
259c3296 108
109 TString kqTypeLabel[fgkqType];
dbe3abbe 110 if (kquark == kCharm){
111 kqTypeLabel[kQuark]="c";
112 kqTypeLabel[kantiQuark]="cbar";
113 kqTypeLabel[kHadron]="cHadron";
114 kqTypeLabel[keHadron]="ceHadron";
115 kqTypeLabel[kDeHadron]="nullHadron";
116 kqTypeLabel[kElectron]="ce";
117 kqTypeLabel[kElectron2nd]="nulle";
118 } else if (kquark == kBeauty){
119 kqTypeLabel[kQuark]="b";
120 kqTypeLabel[kantiQuark]="bbar";
121 kqTypeLabel[kHadron]="bHadron";
122 kqTypeLabel[keHadron]="beHadron";
123 kqTypeLabel[kDeHadron]="bDeHadron";
124 kqTypeLabel[kElectron]="be";
125 kqTypeLabel[kElectron2nd]="bce";
259c3296 126 }
127
128
129 TString hname;
130 for (Int_t iqType = 0; iqType < fgkqType; iqType++ ){
dbe3abbe 131 if (iqType < keHadron && icut > 0) continue; // don't duplicate histogram for quark and hadron
259c3296 132 hname = hnopt+"PdgCode_"+kqTypeLabel[iqType];
dbe3abbe 133 fHist[iq][iqType][icut].fPdgCode = new TH1F(hname,hname,20001,-10000.5,10000.5);
259c3296 134 hname = hnopt+"Pt_"+kqTypeLabel[iqType];
78ea5ef4 135 fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",250,0,50);
259c3296 136 hname = hnopt+"Y_"+kqTypeLabel[iqType];
dbe3abbe 137 fHist[iq][iqType][icut].fY = new TH1F(hname,hname,150,-7.5,7.5);
259c3296 138 hname = hnopt+"Eta_"+kqTypeLabel[iqType];
dbe3abbe 139 fHist[iq][iqType][icut].fEta = new TH1F(hname,hname,150,-7.5,7.5);
70da6c5a 140 // Fill List
141 if(fQAhistos) fHist[iq][iqType][icut].FillList(fQAhistos);
259c3296 142 }
143
dbe3abbe 144 if (icut == 0){
145 hname = hnopt+"Nq_"+kqTypeLabel[kQuark];
146 fHistComm[iq][icut].fNq = new TH1F(hname,hname,10,-0.5,9.5);
147 hname = hnopt+"ProcessID_"+kqTypeLabel[kQuark];
148 fHistComm[iq][icut].fProcessID = new TH1F(hname,hname,21,-10.5,10.5);
149 }
150 hname = hnopt+"ePtRatio_"+kqTypeLabel[kQuark];
151 fHistComm[iq][icut].fePtRatio = new TH2F(hname,hname+";p_{T} (GeV/c);momentum fraction",100,0,30,100,0,1);
152 hname = hnopt+"DePtRatio_"+kqTypeLabel[kQuark];
153 fHistComm[iq][icut].fDePtRatio = new TH2F(hname,hname+";p_{T} (GeV/c);momentum fraction",100,0,30,100,0,1);
154 hname = hnopt+"eDistance_"+kqTypeLabel[kQuark];
155 fHistComm[iq][icut].feDistance= new TH2F(hname,hname+";p_{T} (GeV/c);distance (cm)",100,0,30,200,0,2);
156 hname = hnopt+"DeDistance_"+kqTypeLabel[kQuark];
157 fHistComm[iq][icut].fDeDistance= new TH2F(hname,hname+";p_{T} (GeV/c);distance (cm)",100,0,30,200,0,2);
70da6c5a 158 if(fQAhistos) fHistComm[iq][icut].FillList(fQAhistos);
259c3296 159}
160
161//__________________________________________
162void AliHFEmcQA::Init()
163{
164 // called at begining every event
165
166 for (Int_t i=0; i<2; i++){
167 fIsHeavy[i] = 0;
168 }
169
170 fNparents = 7;
171
172 fParentSelect[0][0] = 411;
173 fParentSelect[0][1] = 421;
174 fParentSelect[0][2] = 431;
175 fParentSelect[0][3] = 4122;
176 fParentSelect[0][4] = 4132;
177 fParentSelect[0][5] = 4232;
178 fParentSelect[0][6] = 4332;
179
180 fParentSelect[1][0] = 511;
181 fParentSelect[1][1] = 521;
182 fParentSelect[1][2] = 531;
183 fParentSelect[1][3] = 5122;
184 fParentSelect[1][4] = 5132;
185 fParentSelect[1][5] = 5232;
186 fParentSelect[1][6] = 5332;
187
188}
189
190//__________________________________________
722347d8 191void AliHFEmcQA::GetQuarkKine(TParticle *part, Int_t iTrack, const Int_t kquark)
259c3296 192{
193 // get heavy quark kinematics
194
dbe3abbe 195 if (kquark != kCharm && kquark != kBeauty) {
259c3296 196 AliDebug(1, "This task is only for heavy quark QA, return\n");
197 return;
198 }
dbe3abbe 199 Int_t iq = kquark - kCharm;
259c3296 200
259c3296 201 if (iTrack < 0) {
202 AliDebug(1, "Stack label is negative, return\n");
203 return;
204 }
205
722347d8 206 //TParticle *part = fStack->Particle(iTrack);
259c3296 207 Int_t partPdgcode = TMath::Abs(part->GetPdgCode());
208
209 // select heavy hadron or not fragmented heavy quark
210 if ( int(partPdgcode/100.)==kquark || int(partPdgcode/1000.)==kquark || (partPdgcode==kquark && (part->GetNDaughters()==0 && iTrack>5)) ){
211
212 TParticle *partMother;
213 Int_t iLabel;
214
215 if (partPdgcode == kquark){ // in case of not fragmented heavy quark
216 partMother = part;
217 iLabel = iTrack;
218 } else{ // in case of heavy hadron, start to search for mother heavy parton
219 iLabel = part->GetFirstMother();
220 if (iLabel>-1) { partMother = fStack->Particle(iLabel); }
221 else {
222 AliDebug(1, "Stack label is negative, return\n");
223 return;
224 }
225 }
226
227 // heavy parton selection as a mother of heavy hadron
228 // if the heavy particle comes from string which is denoted as particle status 12|12|12...12|11,[PYTHIA p.60]
229 // in this case, the mother of heavy particle can be one of the fragmented parton of the string
230 // should I make a condition that partMother should be quark or diquark? -> not necessary
231 if ( abs(partMother->GetPdgCode()) == kquark || (partMother->GetStatusCode() == 11) ){
232 //if ( abs(partMother->GetPdgCode()) == kquark || (partMother->GetStatusCode() == 11 || partMother->GetStatusCode() == 12) ){
233
234 if ( abs(partMother->GetPdgCode()) != kquark ){
235 // search fragmented partons in the same string
dbe3abbe 236 Bool_t isSameString = kTRUE;
259c3296 237 for (Int_t i=1; i<fgkMaxIter; i++){
238 iLabel = iLabel - 1;
239 if (iLabel>-1) { partMother = fStack->Particle(iLabel); }
240 else {
241 AliDebug(1, "Stack label is negative, return\n");
242 return;
243 }
244 if ( abs(partMother->GetPdgCode()) == kquark ) break;
dbe3abbe 245 if ( partMother->GetStatusCode() != 12 ) isSameString = kFALSE;
246 if (!isSameString) return;
259c3296 247 }
248 }
249 AliDebug(1, "Can not find heavy parton of this heavy hadron in the string, return\n");
250 if (abs(partMother->GetPdgCode()) != kquark) return;
251
252 if (fIsHeavy[iq] >= 50) return;
253 fHeavyQuark[fIsHeavy[iq]] = partMother;
254 fIsHeavy[iq]++;
255
256 // fill kinematics for heavy parton
257 if (partMother->GetPdgCode() > 0) { // quark
dbe3abbe 258 fHist[iq][kQuark][0].fPdgCode->Fill(partMother->GetPdgCode());
259 fHist[iq][kQuark][0].fPt->Fill(partMother->Pt());
70da6c5a 260 fHist[iq][kQuark][0].fY->Fill(AliHFEtools::GetRapidity(partMother));
dbe3abbe 261 fHist[iq][kQuark][0].fEta->Fill(partMother->Eta());
259c3296 262 } else{ // antiquark
dbe3abbe 263 fHist[iq][kantiQuark][0].fPdgCode->Fill(partMother->GetPdgCode());
264 fHist[iq][kantiQuark][0].fPt->Fill(partMother->Pt());
70da6c5a 265 fHist[iq][kantiQuark][0].fY->Fill(AliHFEtools::GetRapidity(partMother));
dbe3abbe 266 fHist[iq][kantiQuark][0].fEta->Fill(partMother->Eta());
259c3296 267 }
268
269 } // end of heavy parton slection loop
270
271 } // end of heavy hadron or quark selection
272
273}
274
275//__________________________________________
276void AliHFEmcQA::EndOfEventAna(const Int_t kquark)
277{
278 // end of event analysis
279
dbe3abbe 280 if (kquark != kCharm && kquark != kBeauty) {
259c3296 281 AliDebug(1, "This task is only for heavy quark QA, return\n");
282 return;
283 }
dbe3abbe 284 Int_t iq = kquark - kCharm;
259c3296 285
286
287 // # of heavy quark per event
288 AliDebug(1,Form("Number of heavy quark in this event = %d \n",fIsHeavy[iq]));
dbe3abbe 289 fHistComm[iq][0].fNq->Fill(fIsHeavy[iq]);
259c3296 290
291 Int_t motherID[fgkMaxGener];
292 Int_t motherType[fgkMaxGener];
293 Int_t motherLabel[fgkMaxGener];
294 Int_t ancestorPdg[fgkMaxGener];
295 Int_t ancestorLabel[fgkMaxGener];
296
297 for (Int_t i = 0; i < fgkMaxGener; i++){ // initialization
298 motherID[i] = 0;
299 motherType[i] = 0;
300 motherLabel[i] = 0;
301 ancestorPdg[i] = 0;
302 ancestorLabel[i] = 0;
303 }
304
305 // check history of found heavy quarks
306 for (Int_t i = 0; i < fIsHeavy[iq]; i++){
307
308 ancestorLabel[0] = i;
309 ancestorPdg[0] = fHeavyQuark[i]->GetPdgCode();
310 ancestorLabel[1] = fHeavyQuark[i]->GetFirstMother();
311
312 AliDebug(1,Form("pdg code= %d\n",ancestorPdg[0]));
313 AliDebug(1,Form("ancestor label= %d\n",ancestorLabel[1]));
314
315 Int_t ig = 1;
316 while (ancestorLabel[ig] != -1){
317 // in case there is mother, get mother's pdg code and grandmother's label
318 IdentifyMother(ancestorLabel[ig], ancestorPdg[ig], ancestorLabel[ig+1]);
319 // if mother is still heavy, find again mother's ancestor
320 if (ancestorPdg[ig-1] == ancestorPdg[ig]) {
321 ig++;
322 continue; // if it is from same heavy
323 }
324 // if the heavy's mother is not heavy, check the mother's label to know if it comes from inital or final parton shower
325 if (IsFromInitialShower(ancestorLabel[ig],motherID[i],motherType[i],motherLabel[i])) break;
326 if (IsFromFinalParton(ancestorLabel[ig],motherID[i],motherType[i],motherLabel[i])) break;
327 // if it is not the above case, something is strange
dbe3abbe 328 ReportStrangeness(motherID[i],motherType[i],motherLabel[i]);
259c3296 329 break;
330 }
331 if (ancestorLabel[ig] == -1){ // from hard scattering
332 HardScattering(kquark, motherID[i],motherType[i], motherLabel[i]);
333 }
334
335 } // end of found heavy quark loop
336
337
338 // check process type
339 Int_t processID = 0;
340 for (Int_t i = 0; i < fIsHeavy[iq]; i++){
341 AliDebug(1,Form("Mother ID= %d type= %d label= %d\n",motherID[i],motherType[i],motherLabel[i]));
342 }
343
344
345 Int_t nheavypair = Int_t(fIsHeavy[iq]/2.);
346 for (Int_t ipair = 0; ipair < nheavypair; ipair++){
347
348 Int_t id1 = ipair*2;
349 Int_t id2 = ipair*2 + 1;
350
351 if (motherType[id1] == 2 && motherType[id2] == 2){
dbe3abbe 352 if (motherLabel[id1] == motherLabel[id2]) processID = kGluonSplitting; // gluon spliting
259c3296 353 else processID = -9;
354 }
355 else if (motherType[id1] == -1 && motherType[id2] == -1) {
356 if (motherLabel[id1] == -1 && motherLabel[id2] == -1) {
dbe3abbe 357 if (motherID[id1] == fgkGluon) processID = kPairCreationFromg; // gluon fusion
358 else processID = kPairCreationFromq; // q-qbar pair creation
259c3296 359 }
360 else processID = -8;
361 }
362 else if (motherType[id1] == -1 || motherType[id2] == -1) {
363 if ((motherLabel[id1] == -1 || motherLabel[id2] == -1) && (motherLabel[id1]*motherLabel[id2] == -2 || motherLabel[id1]*motherLabel[id2] == -3)) {
dbe3abbe 364 if(motherID[id1]*motherID[id2] == kquark*fgkGluon) processID = kFlavourExitation; // flavour exitation
365 else processID = kLightQuarkShower;
259c3296 366 }
367 else processID = -7;
368 }
369 else if (motherType[id1] == -2 || motherType[id2] == -2) {
dbe3abbe 370 if (motherLabel[id1] == motherLabel[id2]) processID = kInitialPartonShower; // initial parton shower
259c3296 371 else processID = -6;
372
373 }
374 else processID = -5;
375
376 if (nheavypair >1) AliDebug(1,Form("Multi pair found : process ID = %d\n",processID));
dbe3abbe 377 else fHistComm[iq][0].fProcessID->Fill(processID);
259c3296 378 AliDebug(1,Form("Process ID = %d\n",processID));
379 } // end of # heavy quark pair loop
380
381}
382
383//__________________________________________
722347d8 384void AliHFEmcQA::GetHadronKine(TParticle* mcpart, const Int_t kquark)
dbe3abbe 385{
386 // decay electron kinematics
387
388 if (kquark != kCharm && kquark != kBeauty) {
389 AliDebug(1, "This task is only for heavy quark QA, return\n");
390 return;
391 }
392 Int_t iq = kquark - kCharm;
393
722347d8 394 //TParticle* mcpart = fStack->Particle(iTrack);
dbe3abbe 395
396 Int_t iLabel = mcpart->GetFirstMother();
397 if (iLabel<0){
398 AliDebug(1, "Stack label is negative, return\n");
399 return;
400 }
401
402 TParticle *partCopy = mcpart;
403 Int_t pdgcode = mcpart->GetPdgCode();
404 Int_t pdgcodeCopy = pdgcode;
405
406 // if the mother is charmed hadron
75d81601 407 Bool_t isDirectCharm = kFALSE;
dbe3abbe 408 if ( int(abs(pdgcode)/100.) == kCharm || int(abs(pdgcode)/1000.) == kCharm ) {
409
410 // iterate until you find B hadron as a mother or become top ancester
411 for (Int_t i=1; i<fgkMaxIter; i++){
412
413 Int_t jLabel = mcpart->GetFirstMother();
414 if (jLabel == -1){
75d81601 415 isDirectCharm = kTRUE;
dbe3abbe 416 break; // if there is no ancester
417 }
418 if (jLabel < 0){ // safety protection
419 AliDebug(1, "Stack label is negative, return\n");
420 return;
421 }
422 // if there is an ancester
423 TParticle* mother = fStack->Particle(jLabel);
424 Int_t motherPDG = mother->GetPdgCode();
425
426 for (Int_t j=0; j<fNparents; j++){
427 if (abs(motherPDG)==fParentSelect[1][j]) return; // return if this hadron is originated from b
428 }
429
430 mcpart = mother;
431 } // end of iteration
432 } // end of if
75d81601 433 if((isDirectCharm == kTRUE && kquark == kCharm) || kquark == kBeauty) {
dbe3abbe 434 for (Int_t i=0; i<fNparents; i++){
435 if (abs(pdgcodeCopy)==fParentSelect[iq][i]){
436
437 // fill hadron kinematics
438 fHist[iq][kHadron][0].fPdgCode->Fill(pdgcodeCopy);
439 fHist[iq][kHadron][0].fPt->Fill(partCopy->Pt());
70da6c5a 440 fHist[iq][kHadron][0].fY->Fill(AliHFEtools::GetRapidity(partCopy));
dbe3abbe 441 fHist[iq][kHadron][0].fEta->Fill(partCopy->Eta());
442
443 }
444 }
445 } // end of if
446}
447
448//__________________________________________
722347d8 449void AliHFEmcQA::GetDecayedKine(TParticle* mcpart, const Int_t kquark, Int_t kdecayed, Int_t icut)
259c3296 450{
451 // decay electron kinematics
452
dbe3abbe 453 if (kquark != kCharm && kquark != kBeauty) {
259c3296 454 AliDebug(1, "This task is only for heavy quark QA, return\n");
455 return;
456 }
dbe3abbe 457 Int_t iq = kquark - kCharm;
50685501 458 Bool_t isFinalOpenCharm = kFALSE;
259c3296 459
722347d8 460/*
259c3296 461 if (iTrack < 0) {
462 AliDebug(1, "Stack label is negative, return\n");
463 return;
464 }
722347d8 465 */
259c3296 466
722347d8 467 //TParticle* mcpart = fStack->Particle(iTrack);
259c3296 468
469 if ( abs(mcpart->GetPdgCode()) != kdecayed ) return;
259c3296 470
471 Int_t iLabel = mcpart->GetFirstMother();
472 if (iLabel<0){
473 AliDebug(1, "Stack label is negative, return\n");
474 return;
475 }
476
477 TParticle *partMother = fStack->Particle(iLabel);
dbe3abbe 478 TParticle *partMotherCopy = partMother;
259c3296 479 Int_t maPdgcode = partMother->GetPdgCode();
dbe3abbe 480 Int_t maPdgcodeCopy = maPdgcode;
259c3296 481
9bcfd1ab 482 // get mc primary vertex
483 TArrayF mcPrimVtx;
484 if(fMCHeader) fMCHeader->PrimaryVertex(mcPrimVtx);
485
259c3296 486 // get electron production vertex
487 TLorentzVector ePoint;
488 mcpart->ProductionVertex(ePoint);
489
9bcfd1ab 490 // calculated production vertex to primary vertex (in xy plane)
491 Float_t decayLxy = TMath::Sqrt((mcPrimVtx[0]-ePoint[0])*(mcPrimVtx[0]-ePoint[0])+(mcPrimVtx[1]-ePoint[1])*(mcPrimVtx[1]-ePoint[1]));
492
259c3296 493 // if the mother is charmed hadron
dbe3abbe 494 Bool_t isMotherDirectCharm = kFALSE;
495 if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
259c3296 496
0792aa82 497 for (Int_t i=0; i<fNparents; i++){
498 if (abs(maPdgcode)==fParentSelect[0][i]){
50685501 499 isFinalOpenCharm = kTRUE;
0792aa82 500 }
501 }
50685501 502 if (!isFinalOpenCharm) return ;
0792aa82 503
259c3296 504 // iterate until you find B hadron as a mother or become top ancester
505 for (Int_t i=1; i<fgkMaxIter; i++){
506
507 Int_t jLabel = partMother->GetFirstMother();
508 if (jLabel == -1){
dbe3abbe 509 isMotherDirectCharm = kTRUE;
259c3296 510 break; // if there is no ancester
511 }
512 if (jLabel < 0){ // safety protection
513 AliDebug(1, "Stack label is negative, return\n");
514 return;
515 }
516
517 // if there is an ancester
518 TParticle* grandMa = fStack->Particle(jLabel);
519 Int_t grandMaPDG = grandMa->GetPdgCode();
520
fc0de2a0 521 for (Int_t j=0; j<fNparents; j++){
522 if (abs(grandMaPDG)==fParentSelect[1][j]){
259c3296 523
dbe3abbe 524 if (kquark == kCharm) return;
259c3296 525 // fill electron kinematics
dbe3abbe 526 fHist[iq][kElectron2nd][icut].fPdgCode->Fill(mcpart->GetPdgCode());
527 fHist[iq][kElectron2nd][icut].fPt->Fill(mcpart->Pt());
70da6c5a 528 fHist[iq][kElectron2nd][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart));
dbe3abbe 529 fHist[iq][kElectron2nd][icut].fEta->Fill(mcpart->Eta());
259c3296 530
531 // fill mother hadron kinematics
dbe3abbe 532 fHist[iq][kDeHadron][icut].fPdgCode->Fill(grandMaPDG);
533 fHist[iq][kDeHadron][icut].fPt->Fill(grandMa->Pt());
70da6c5a 534 fHist[iq][kDeHadron][icut].fY->Fill(AliHFEtools::GetRapidity(grandMa));
dbe3abbe 535 fHist[iq][kDeHadron][icut].fEta->Fill(grandMa->Eta());
259c3296 536
537 // ratio between pT of electron and pT of mother B hadron
dbe3abbe 538 if(grandMa->Pt()) fHistComm[iq][icut].fDePtRatio->Fill(grandMa->Pt(),mcpart->Pt()/grandMa->Pt());
259c3296 539
9bcfd1ab 540 // distance between electron production point and primary vertex
541 fHistComm[iq][icut].fDeDistance->Fill(grandMa->Pt(),decayLxy);
259c3296 542 return;
543 }
dc306130 544 }
259c3296 545
546 partMother = grandMa;
547 } // end of iteration
548 } // end of if
dbe3abbe 549 if((isMotherDirectCharm == kTRUE && kquark == kCharm) || kquark == kBeauty) {
259c3296 550 for (Int_t i=0; i<fNparents; i++){
551 if (abs(maPdgcodeCopy)==fParentSelect[iq][i]){
552
553 // fill electron kinematics
dbe3abbe 554 fHist[iq][kElectron][icut].fPdgCode->Fill(mcpart->GetPdgCode());
555 fHist[iq][kElectron][icut].fPt->Fill(mcpart->Pt());
70da6c5a 556 fHist[iq][kElectron][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart));
dbe3abbe 557 fHist[iq][kElectron][icut].fEta->Fill(mcpart->Eta());
259c3296 558
559 // fill mother hadron kinematics
dbe3abbe 560 fHist[iq][keHadron][icut].fPdgCode->Fill(maPdgcodeCopy);
561 fHist[iq][keHadron][icut].fPt->Fill(partMotherCopy->Pt());
70da6c5a 562 fHist[iq][keHadron][icut].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
dbe3abbe 563 fHist[iq][keHadron][icut].fEta->Fill(partMotherCopy->Eta());
259c3296 564
565 // ratio between pT of electron and pT of mother B hadron
dbe3abbe 566 if(partMotherCopy->Pt()) fHistComm[iq][icut].fePtRatio->Fill(partMotherCopy->Pt(),mcpart->Pt()/partMotherCopy->Pt());
259c3296 567
9bcfd1ab 568 // distance between electron production point and primary vertex
569 fHistComm[iq][icut].feDistance->Fill(partMotherCopy->Pt(),decayLxy);
259c3296 570 }
571 }
572 } // end of if
573}
574
0792aa82 575//____________________________________________________________________
576void AliHFEmcQA::GetDecayedKine(AliAODMCParticle *mcpart, const Int_t kquark, Int_t kdecayed, Int_t icut)
577{
578 // decay electron kinematics
579
580 if (kquark != kCharm && kquark != kBeauty) {
581 AliDebug(1, "This task is only for heavy quark QA, return\n");
582 return;
583 }
584
585 Int_t iq = kquark - kCharm;
50685501 586 Bool_t isFinalOpenCharm = kFALSE;
0792aa82 587
588 if ( abs(mcpart->GetPdgCode()) != kdecayed ) return;
589
590 // mother
591 Int_t iLabel = mcpart->GetMother();
592 if (iLabel<0){
593 AliDebug(1, "Stack label is negative, return\n");
594 return;
595 }
596
597 AliAODMCParticle *partMother = (AliAODMCParticle*)fMCArray->At(iLabel);
598 AliAODMCParticle *partMotherCopy = partMother;
599 Int_t maPdgcode = partMother->GetPdgCode();
600 Int_t maPdgcodeCopy = maPdgcode;
601
602 Bool_t isMotherDirectCharm = kFALSE;
603 if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
604
605 for (Int_t i=0; i<fNparents; i++){
606 if (abs(maPdgcode)==fParentSelect[0][i]){
50685501 607 isFinalOpenCharm = kTRUE;
0792aa82 608 }
609 }
50685501 610 if (!isFinalOpenCharm) return;
0792aa82 611
612 for (Int_t i=1; i<fgkMaxIter; i++){
613
614 Int_t jLabel = partMother->GetMother();
615 if (jLabel == -1){
616 isMotherDirectCharm = kTRUE;
617 break; // if there is no ancester
618 }
619 if (jLabel < 0){ // safety protection
620 AliDebug(1, "Stack label is negative, return\n");
621 return;
622 }
623
624 // if there is an ancester
625 AliAODMCParticle* grandMa = (AliAODMCParticle*)fMCArray->At(jLabel);
626 Int_t grandMaPDG = grandMa->GetPdgCode();
627
628 for (Int_t j=0; j<fNparents; j++){
629 if (abs(grandMaPDG)==fParentSelect[1][j]){
630
631 if (kquark == kCharm) return;
632 // fill electron kinematics
633 fHist[iq][kElectron2nd][icut].fPdgCode->Fill(mcpart->GetPdgCode());
634 fHist[iq][kElectron2nd][icut].fPt->Fill(mcpart->Pt());
70da6c5a 635 fHist[iq][kElectron2nd][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart));
0792aa82 636 fHist[iq][kElectron2nd][icut].fEta->Fill(mcpart->Eta());
637
638 // fill mother hadron kinematics
639 fHist[iq][kDeHadron][icut].fPdgCode->Fill(grandMaPDG);
640 fHist[iq][kDeHadron][icut].fPt->Fill(grandMa->Pt());
70da6c5a 641 fHist[iq][kDeHadron][icut].fY->Fill(AliHFEtools::GetRapidity(grandMa));
0792aa82 642 fHist[iq][kDeHadron][icut].fEta->Fill(grandMa->Eta());
643
644 return;
645 }
646 }
647
648 partMother = grandMa;
649 } // end of iteration
650 } // end of if
651 if ((isMotherDirectCharm == kTRUE && kquark == kCharm) || kquark == kBeauty) {
652 for (Int_t i=0; i<fNparents; i++){
653 if (abs(maPdgcodeCopy)==fParentSelect[iq][i]){
654
655 // fill electron kinematics
656 fHist[iq][kElectron][icut].fPdgCode->Fill(mcpart->GetPdgCode());
657 fHist[iq][kElectron][icut].fPt->Fill(mcpart->Pt());
70da6c5a 658 fHist[iq][kElectron][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart));
0792aa82 659 fHist[iq][kElectron][icut].fEta->Fill(mcpart->Eta());
660
661 // fill mother hadron kinematics
662 fHist[iq][keHadron][icut].fPdgCode->Fill(maPdgcodeCopy);
663 fHist[iq][keHadron][icut].fPt->Fill(partMotherCopy->Pt());
70da6c5a 664 fHist[iq][keHadron][icut].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
0792aa82 665 fHist[iq][keHadron][icut].fEta->Fill(partMotherCopy->Eta());
666
667 }
668 }
669 } // end of if
670
671}
259c3296 672
673//__________________________________________
75d81601 674void AliHFEmcQA::IdentifyMother(Int_t motherlabel, Int_t &motherpdg, Int_t &grandmotherlabel)
259c3296 675{
dbe3abbe 676 // find mother pdg code and label
677
75d81601 678 if (motherlabel < 0) {
259c3296 679 AliDebug(1, "Stack label is negative, return\n");
680 return;
681 }
75d81601 682 TParticle *heavysMother = fStack->Particle(motherlabel);
683 motherpdg = heavysMother->GetPdgCode();
684 grandmotherlabel = heavysMother->GetFirstMother();
685 AliDebug(1,Form("ancestor pdg code= %d\n",motherpdg));
259c3296 686}
687
688//__________________________________________
259c3296 689void AliHFEmcQA::HardScattering(const Int_t kquark, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
690{
dbe3abbe 691 // mothertype -1 means this heavy quark coming from hard vertex
692
259c3296 693 TParticle *afterinitialrad1 = fStack->Particle(4);
694 TParticle *afterinitialrad2 = fStack->Particle(5);
695
696 motherlabel = -1;
697
698 if (abs(afterinitialrad1->GetPdgCode()) == fgkGluon && abs(afterinitialrad2->GetPdgCode()) == fgkGluon){
699 AliDebug(1,"heavy from gluon gluon pair creation!\n");
700 mothertype = -1;
701 motherID = fgkGluon;
702 }
703 else if (abs(afterinitialrad1->GetPdgCode()) == kquark || abs(afterinitialrad2->GetPdgCode()) == kquark){ // one from Q and the other from g
704 AliDebug(1,"heavy from flavor exitation!\n");
705 mothertype = -1;
706 motherID = kquark;
707 }
708 else if (abs(afterinitialrad1->GetPdgCode()) == abs(afterinitialrad2->GetPdgCode())){
709 AliDebug(1,"heavy from q-qbar pair creation!\n");
710 mothertype = -1;
711 motherID = 1;
712 }
713 else {
714 AliDebug(1,"something strange!\n");
715 mothertype = -999;
716 motherlabel = -999;
717 motherID = -999;
718 }
719}
720
721//__________________________________________
259c3296 722Bool_t AliHFEmcQA::IsFromInitialShower(Int_t inputmotherlabel, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
723{
dbe3abbe 724 // mothertype -2 means this heavy quark coming from initial state
725
259c3296 726 if (inputmotherlabel==2 || inputmotherlabel==3){ // mother exist before initial state radiation
727 TParticle *heavysMother = fStack->Particle(inputmotherlabel);
728 motherID = heavysMother->GetPdgCode();
729 mothertype = -2; // there is mother before initial state radiation
730 motherlabel = inputmotherlabel;
731 AliDebug(1,"initial parton shower! \n");
732
733 return kTRUE;
734 }
735
736 return kFALSE;
737}
738
739//__________________________________________
259c3296 740Bool_t AliHFEmcQA::IsFromFinalParton(Int_t inputmotherlabel, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
741{
dbe3abbe 742 // mothertype 2 means this heavy quark coming from final state
743
259c3296 744 if (inputmotherlabel > 5){ // mother exist after hard scattering
745 TParticle *heavysMother = fStack->Particle(inputmotherlabel);
746 motherID = heavysMother->GetPdgCode();
747 mothertype = 2; //
748 motherlabel = inputmotherlabel;
749 AliDebug(1,Form("heavy quark from %d after hard scattering! \n",motherID));
750
751 return kTRUE;
752 }
753 return kFALSE;
754}
755
756//__________________________________________
dbe3abbe 757void AliHFEmcQA::ReportStrangeness(Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
259c3296 758{
dbe3abbe 759 // mark strange behavior
760
259c3296 761 mothertype = -888;
762 motherlabel = -888;
763 motherID = -888;
764 AliDebug(1,"something strange!\n");
765}
766
0792aa82 767//__________________________________________
9bcfd1ab 768Int_t AliHFEmcQA::GetSource(AliAODMCParticle * const mcpart)
0792aa82 769{
9bcfd1ab 770 // decay particle's origin
0792aa82 771
9bcfd1ab 772 //if ( abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return -1;
0792aa82 773
774 Int_t origin = -1;
50685501 775 Bool_t isFinalOpenCharm = kFALSE;
0792aa82 776
777 // mother
778 Int_t iLabel = mcpart->GetMother();
779 if (iLabel<0){
780 AliDebug(1, "Stack label is negative, return\n");
781 return -1;
782 }
783
784 AliAODMCParticle *partMother = (AliAODMCParticle*)fMCArray->At(iLabel);
785 Int_t maPdgcode = partMother->GetPdgCode();
786
787 // if the mother is charmed hadron
788 if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
789
790 for (Int_t i=0; i<fNparents; i++){
791 if (abs(maPdgcode)==fParentSelect[0][i]){
50685501 792 isFinalOpenCharm = kTRUE;
0792aa82 793 }
794 }
50685501 795 if (!isFinalOpenCharm) return -1;
0792aa82 796
797 // iterate until you find B hadron as a mother or become top ancester
798 for (Int_t i=1; i<fgkMaxIter; i++){
799
800 Int_t jLabel = partMother->GetMother();
801 if (jLabel == -1){
802 origin = kDirectCharm;
803 return origin;
804 }
805 if (jLabel < 0){ // safety protection
806 AliDebug(1, "Stack label is negative, return\n");
807 return -1;
808 }
809
810 // if there is an ancester
811 AliAODMCParticle* grandMa = (AliAODMCParticle*)fMCArray->At(jLabel);
812 Int_t grandMaPDG = grandMa->GetPdgCode();
813
70b5ea26 814 for (Int_t j=0; j<fNparents; j++){
815 if (abs(grandMaPDG)==fParentSelect[1][j]){
0792aa82 816 origin = kBeautyCharm;
817 return origin;
818 }
819 }
820
821 partMother = grandMa;
822 } // end of iteration
823 } // end of if
824 else if ( int(abs(maPdgcode)/100.) == kBeauty || int(abs(maPdgcode)/1000.) == kBeauty ) {
825 for (Int_t i=0; i<fNparents; i++){
826 if (abs(maPdgcode)==fParentSelect[1][i]){
827 origin = kDirectBeauty;
828 return origin;
829 }
830 }
831 } // end of if
832 else if ( abs(maPdgcode) == 22 ) {
833 origin = kGamma;
834 return origin;
835 } // end of if
836 else if ( abs(maPdgcode) == 111 ) {
837 origin = kPi0;
838 return origin;
839 } // end of if
840
841 return origin;
842}
843
844//__________________________________________
9bcfd1ab 845Int_t AliHFEmcQA::GetSource(TParticle * const mcpart)
0792aa82 846{
9bcfd1ab 847 // decay particle's origin
0792aa82 848
9bcfd1ab 849 //if ( abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return -1;
0792aa82 850
851 Int_t origin = -1;
50685501 852 Bool_t isFinalOpenCharm = kFALSE;
0792aa82 853
854 Int_t iLabel = mcpart->GetFirstMother();
855 if (iLabel<0){
856 AliDebug(1, "Stack label is negative, return\n");
857 return -1;
858 }
859
860 TParticle *partMother = fStack->Particle(iLabel);
861 Int_t maPdgcode = partMother->GetPdgCode();
862
863 // if the mother is charmed hadron
864 if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
865
866 for (Int_t i=0; i<fNparents; i++){
867 if (abs(maPdgcode)==fParentSelect[0][i]){
50685501 868 isFinalOpenCharm = kTRUE;
0792aa82 869 }
870 }
50685501 871 if (!isFinalOpenCharm) return -1;
0792aa82 872
873 // iterate until you find B hadron as a mother or become top ancester
874 for (Int_t i=1; i<fgkMaxIter; i++){
875
876 Int_t jLabel = partMother->GetFirstMother();
877 if (jLabel == -1){
878 origin = kDirectCharm;
879 return origin;
880 }
881 if (jLabel < 0){ // safety protection
882 AliDebug(1, "Stack label is negative, return\n");
883 return -1;
884 }
885
886 // if there is an ancester
887 TParticle* grandMa = fStack->Particle(jLabel);
888 Int_t grandMaPDG = grandMa->GetPdgCode();
889
70b5ea26 890 for (Int_t j=0; j<fNparents; j++){
891 if (abs(grandMaPDG)==fParentSelect[1][j]){
0792aa82 892 origin = kBeautyCharm;
893 return origin;
894 }
895 }
896
897 partMother = grandMa;
898 } // end of iteration
899 } // end of if
900 else if ( int(abs(maPdgcode)/100.) == kBeauty || int(abs(maPdgcode)/1000.) == kBeauty ) {
901 for (Int_t i=0; i<fNparents; i++){
902 if (abs(maPdgcode)==fParentSelect[1][i]){
903 origin = kDirectBeauty;
904 return origin;
905 }
906 }
907 } // end of if
908 else if ( abs(maPdgcode) == 22 ) {
909 origin = kGamma;
910 return origin;
911 } // end of if
912 else if ( abs(maPdgcode) == 111 ) {
913 origin = kPi0;
914 return origin;
915 } // end of if
916 else origin = kElse;
917
918 return origin;
919}
70da6c5a 920
921//__________________________________________
922void AliHFEmcQA::MakeHistograms(){
923 //
924 // Create the QA histograms
925 //
926 fQAhistos = new TList;
927 fQAhistos->SetName("MCqa");
928
929 CreateHistograms(AliHFEmcQA::kCharm,0,"mcqa_"); // create histograms for charm
930 CreateHistograms(AliHFEmcQA::kBeauty,0,"mcqa_"); // create histograms for beauty
931 CreateHistograms(AliHFEmcQA::kCharm,1,"mcqa_barrel_"); // create histograms for charm
932 CreateHistograms(AliHFEmcQA::kBeauty,1,"mcqa_barrel_"); // create histograms for beauty
933 CreateHistograms(AliHFEmcQA::kCharm,2,"mcqa_unitY_"); // create histograms for charm
934 CreateHistograms(AliHFEmcQA::kBeauty,2,"mcqa_unitY_"); // create histograms for beauty
935 CreateHistograms(AliHFEmcQA::kCharm,3,"mcqa_reccut_"); // create histograms for charm
936 CreateHistograms(AliHFEmcQA::kBeauty,3,"mcqa_reccut_"); // create histograms for beauty
937 CreateHistograms(AliHFEmcQA::kCharm,4,"mcqa_recpidcut_"); // create histograms for charm
938 CreateHistograms(AliHFEmcQA::kBeauty,4,"mcqa_recpidcut_"); // create histograms for beauty
939}
940
941//__________________________________________
942void AliHFEmcQA::AliHists::FillList(TList *l) const {
943 //
944 // Fill Histos into a list for output
945 //
946 if(fPdgCode) l->Add(fPdgCode);
947 if(fPt) l->Add(fPt);
948 if(fY) l->Add(fY);
949 if(fEta) l->Add(fEta);
950}
951
952//__________________________________________
953void AliHFEmcQA::AliHistsComm::FillList(TList *l) const {
954 //
955 // Fill Histos into a list for output
956 //
957 if(fNq) l->Add(fNq);
958 if(fProcessID) l->Add(fProcessID);
959 if(fePtRatio) l->Add(fePtRatio);
960 if(fDePtRatio) l->Add(fDePtRatio);
961 if(feDistance) l->Add(feDistance);
962 if(fDeDistance) l->Add(fDeDistance);
963}
964