Since it contains fixes of coding rule violations, all classes are involved. Further...
[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 **************************************************************************/
15/**************************************************************************
16 * *
17 * QA class of Heavy Flavor quark and fragmeted/decayed particles *
dbe3abbe 18 * -Check kinematics of Heavy Quarks/hadrons, and decayed leptons *
19 * pT, rapidity *
20 * decay lepton kinematics w/wo acceptance *
21 * heavy hadron decay length, electron pT fraction carried from decay *
22 * -Check yield of Heavy Quarks/hadrons *
23 * Number of produced heavy quark *
24 * Number of produced hadron of given pdg code *
25 * *
259c3296 26 * *
27 * Authors: *
28 * MinJung Kweon <minjung@physi.uni-heidelberg.de> *
29 * *
30 **************************************************************************/
31
32
33#include <iostream>
34#include <TH2F.h>
35#include <TCanvas.h>
36#include <TFile.h>
37#include <TH1F.h>
38#include <TH2F.h>
39#include <TCut.h>
40#include <TRandom.h>
41#include <TDatabasePDG.h>
42#include <TROOT.h>
43#include <TParticle.h>
44
45#include <AliLog.h>
46#include <AliESDEvent.h>
47#include <AliESDtrack.h>
48#include <AliESDInputHandler.h>
49#include <AliESDVertex.h>
50#include <AliStack.h>
51
52#include "AliHFEmcQA.h"
53
54const Int_t AliHFEmcQA::fgkGluon=21;
55const Int_t AliHFEmcQA::fgkMaxGener=10;
dbe3abbe 56const Int_t AliHFEmcQA::fgkMaxIter=100;
57const Int_t AliHFEmcQA::fgkqType = 7; // number of species waiting for QA done
259c3296 58
59ClassImp(AliHFEmcQA)
60
61//_______________________________________________________________________________________________
62AliHFEmcQA::AliHFEmcQA() :
dbe3abbe 63 fStack(0x0)
259c3296 64 ,fNparents(0)
65{
66 // Default constructor
67
259c3296 68}
69
70//_______________________________________________________________________________________________
71AliHFEmcQA::AliHFEmcQA(const AliHFEmcQA&p):
72 TObject(p)
259c3296 73 ,fStack(0x0)
74 ,fNparents(p.fNparents)
75{
76 // Copy constructor
77}
78
79//_______________________________________________________________________________________________
80AliHFEmcQA&
81AliHFEmcQA::operator=(const AliHFEmcQA &)
82{
83 // Assignment operator
84
85 AliInfo("Not yet implemented.");
86 return *this;
87}
88
89//_______________________________________________________________________________________________
90AliHFEmcQA::~AliHFEmcQA()
91{
92 // Destructor
93
94 AliInfo("Analysis Done.");
95}
96
97//_______________________________________________________________________________________________
dc306130 98void AliHFEmcQA::PostAnalyze()
259c3296 99{
100}
101
102//__________________________________________
dbe3abbe 103void AliHFEmcQA::CreateHistograms(const Int_t kquark, Int_t icut, TString hnopt)
259c3296 104{
105 // create histograms
106
dbe3abbe 107 if (kquark != kCharm && kquark != kBeauty) {
259c3296 108 AliDebug(1, "This task is only for heavy quark QA, return\n");
109 return;
110 }
dbe3abbe 111 Int_t iq = kquark - kCharm;
259c3296 112
113 TString kqTypeLabel[fgkqType];
dbe3abbe 114 if (kquark == kCharm){
115 kqTypeLabel[kQuark]="c";
116 kqTypeLabel[kantiQuark]="cbar";
117 kqTypeLabel[kHadron]="cHadron";
118 kqTypeLabel[keHadron]="ceHadron";
119 kqTypeLabel[kDeHadron]="nullHadron";
120 kqTypeLabel[kElectron]="ce";
121 kqTypeLabel[kElectron2nd]="nulle";
122 } else if (kquark == kBeauty){
123 kqTypeLabel[kQuark]="b";
124 kqTypeLabel[kantiQuark]="bbar";
125 kqTypeLabel[kHadron]="bHadron";
126 kqTypeLabel[keHadron]="beHadron";
127 kqTypeLabel[kDeHadron]="bDeHadron";
128 kqTypeLabel[kElectron]="be";
129 kqTypeLabel[kElectron2nd]="bce";
259c3296 130 }
131
132
133 TString hname;
134 for (Int_t iqType = 0; iqType < fgkqType; iqType++ ){
dbe3abbe 135 if (iqType < keHadron && icut > 0) continue; // don't duplicate histogram for quark and hadron
259c3296 136 hname = hnopt+"PdgCode_"+kqTypeLabel[iqType];
dbe3abbe 137 fHist[iq][iqType][icut].fPdgCode = new TH1F(hname,hname,20001,-10000.5,10000.5);
259c3296 138 hname = hnopt+"Pt_"+kqTypeLabel[iqType];
dbe3abbe 139 fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",150,0,30);
259c3296 140 hname = hnopt+"Y_"+kqTypeLabel[iqType];
dbe3abbe 141 fHist[iq][iqType][icut].fY = new TH1F(hname,hname,150,-7.5,7.5);
259c3296 142 hname = hnopt+"Eta_"+kqTypeLabel[iqType];
dbe3abbe 143 fHist[iq][iqType][icut].fEta = new TH1F(hname,hname,150,-7.5,7.5);
259c3296 144 }
145
dbe3abbe 146 if (icut == 0){
147 hname = hnopt+"Nq_"+kqTypeLabel[kQuark];
148 fHistComm[iq][icut].fNq = new TH1F(hname,hname,10,-0.5,9.5);
149 hname = hnopt+"ProcessID_"+kqTypeLabel[kQuark];
150 fHistComm[iq][icut].fProcessID = new TH1F(hname,hname,21,-10.5,10.5);
151 }
152 hname = hnopt+"ePtRatio_"+kqTypeLabel[kQuark];
153 fHistComm[iq][icut].fePtRatio = new TH2F(hname,hname+";p_{T} (GeV/c);momentum fraction",100,0,30,100,0,1);
154 hname = hnopt+"DePtRatio_"+kqTypeLabel[kQuark];
155 fHistComm[iq][icut].fDePtRatio = new TH2F(hname,hname+";p_{T} (GeV/c);momentum fraction",100,0,30,100,0,1);
156 hname = hnopt+"eDistance_"+kqTypeLabel[kQuark];
157 fHistComm[iq][icut].feDistance= new TH2F(hname,hname+";p_{T} (GeV/c);distance (cm)",100,0,30,200,0,2);
158 hname = hnopt+"DeDistance_"+kqTypeLabel[kQuark];
159 fHistComm[iq][icut].fDeDistance= new TH2F(hname,hname+";p_{T} (GeV/c);distance (cm)",100,0,30,200,0,2);
259c3296 160
161}
162
163//__________________________________________
164void AliHFEmcQA::Init()
165{
166 // called at begining every event
167
168 for (Int_t i=0; i<2; i++){
169 fIsHeavy[i] = 0;
170 }
171
172 fNparents = 7;
173
174 fParentSelect[0][0] = 411;
175 fParentSelect[0][1] = 421;
176 fParentSelect[0][2] = 431;
177 fParentSelect[0][3] = 4122;
178 fParentSelect[0][4] = 4132;
179 fParentSelect[0][5] = 4232;
180 fParentSelect[0][6] = 4332;
181
182 fParentSelect[1][0] = 511;
183 fParentSelect[1][1] = 521;
184 fParentSelect[1][2] = 531;
185 fParentSelect[1][3] = 5122;
186 fParentSelect[1][4] = 5132;
187 fParentSelect[1][5] = 5232;
188 fParentSelect[1][6] = 5332;
189
190}
191
192//__________________________________________
193void AliHFEmcQA::GetQuarkKine(Int_t iTrack, const Int_t kquark)
194{
195 // get heavy quark kinematics
196
dbe3abbe 197 if (kquark != kCharm && kquark != kBeauty) {
259c3296 198 AliDebug(1, "This task is only for heavy quark QA, return\n");
199 return;
200 }
dbe3abbe 201 Int_t iq = kquark - kCharm;
259c3296 202
203
204 if (iTrack < 0) {
205 AliDebug(1, "Stack label is negative, return\n");
206 return;
207 }
208
209 TParticle *part = fStack->Particle(iTrack);
210 Int_t partPdgcode = TMath::Abs(part->GetPdgCode());
211
212 // select heavy hadron or not fragmented heavy quark
213 if ( int(partPdgcode/100.)==kquark || int(partPdgcode/1000.)==kquark || (partPdgcode==kquark && (part->GetNDaughters()==0 && iTrack>5)) ){
214
215 TParticle *partMother;
216 Int_t iLabel;
217
218 if (partPdgcode == kquark){ // in case of not fragmented heavy quark
219 partMother = part;
220 iLabel = iTrack;
221 } else{ // in case of heavy hadron, start to search for mother heavy parton
222 iLabel = part->GetFirstMother();
223 if (iLabel>-1) { partMother = fStack->Particle(iLabel); }
224 else {
225 AliDebug(1, "Stack label is negative, return\n");
226 return;
227 }
228 }
229
230 // heavy parton selection as a mother of heavy hadron
231 // if the heavy particle comes from string which is denoted as particle status 12|12|12...12|11,[PYTHIA p.60]
232 // in this case, the mother of heavy particle can be one of the fragmented parton of the string
233 // should I make a condition that partMother should be quark or diquark? -> not necessary
234 if ( abs(partMother->GetPdgCode()) == kquark || (partMother->GetStatusCode() == 11) ){
235 //if ( abs(partMother->GetPdgCode()) == kquark || (partMother->GetStatusCode() == 11 || partMother->GetStatusCode() == 12) ){
236
237 if ( abs(partMother->GetPdgCode()) != kquark ){
238 // search fragmented partons in the same string
dbe3abbe 239 Bool_t isSameString = kTRUE;
259c3296 240 for (Int_t i=1; i<fgkMaxIter; i++){
241 iLabel = iLabel - 1;
242 if (iLabel>-1) { partMother = fStack->Particle(iLabel); }
243 else {
244 AliDebug(1, "Stack label is negative, return\n");
245 return;
246 }
247 if ( abs(partMother->GetPdgCode()) == kquark ) break;
dbe3abbe 248 if ( partMother->GetStatusCode() != 12 ) isSameString = kFALSE;
249 if (!isSameString) return;
259c3296 250 }
251 }
252 AliDebug(1, "Can not find heavy parton of this heavy hadron in the string, return\n");
253 if (abs(partMother->GetPdgCode()) != kquark) return;
254
255 if (fIsHeavy[iq] >= 50) return;
256 fHeavyQuark[fIsHeavy[iq]] = partMother;
257 fIsHeavy[iq]++;
258
259 // fill kinematics for heavy parton
260 if (partMother->GetPdgCode() > 0) { // quark
dbe3abbe 261 fHist[iq][kQuark][0].fPdgCode->Fill(partMother->GetPdgCode());
262 fHist[iq][kQuark][0].fPt->Fill(partMother->Pt());
263 fHist[iq][kQuark][0].fY->Fill(GetRapidity(partMother));
264 fHist[iq][kQuark][0].fEta->Fill(partMother->Eta());
259c3296 265 } else{ // antiquark
dbe3abbe 266 fHist[iq][kantiQuark][0].fPdgCode->Fill(partMother->GetPdgCode());
267 fHist[iq][kantiQuark][0].fPt->Fill(partMother->Pt());
268 fHist[iq][kantiQuark][0].fY->Fill(GetRapidity(partMother));
269 fHist[iq][kantiQuark][0].fEta->Fill(partMother->Eta());
259c3296 270 }
271
272 } // end of heavy parton slection loop
273
274 } // end of heavy hadron or quark selection
275
276}
277
278//__________________________________________
279void AliHFEmcQA::EndOfEventAna(const Int_t kquark)
280{
281 // end of event analysis
282
dbe3abbe 283 if (kquark != kCharm && kquark != kBeauty) {
259c3296 284 AliDebug(1, "This task is only for heavy quark QA, return\n");
285 return;
286 }
dbe3abbe 287 Int_t iq = kquark - kCharm;
259c3296 288
289
290 // # of heavy quark per event
291 AliDebug(1,Form("Number of heavy quark in this event = %d \n",fIsHeavy[iq]));
dbe3abbe 292 fHistComm[iq][0].fNq->Fill(fIsHeavy[iq]);
259c3296 293
294 Int_t motherID[fgkMaxGener];
295 Int_t motherType[fgkMaxGener];
296 Int_t motherLabel[fgkMaxGener];
297 Int_t ancestorPdg[fgkMaxGener];
298 Int_t ancestorLabel[fgkMaxGener];
299
300 for (Int_t i = 0; i < fgkMaxGener; i++){ // initialization
301 motherID[i] = 0;
302 motherType[i] = 0;
303 motherLabel[i] = 0;
304 ancestorPdg[i] = 0;
305 ancestorLabel[i] = 0;
306 }
307
308 // check history of found heavy quarks
309 for (Int_t i = 0; i < fIsHeavy[iq]; i++){
310
311 ancestorLabel[0] = i;
312 ancestorPdg[0] = fHeavyQuark[i]->GetPdgCode();
313 ancestorLabel[1] = fHeavyQuark[i]->GetFirstMother();
314
315 AliDebug(1,Form("pdg code= %d\n",ancestorPdg[0]));
316 AliDebug(1,Form("ancestor label= %d\n",ancestorLabel[1]));
317
318 Int_t ig = 1;
319 while (ancestorLabel[ig] != -1){
320 // in case there is mother, get mother's pdg code and grandmother's label
321 IdentifyMother(ancestorLabel[ig], ancestorPdg[ig], ancestorLabel[ig+1]);
322 // if mother is still heavy, find again mother's ancestor
323 if (ancestorPdg[ig-1] == ancestorPdg[ig]) {
324 ig++;
325 continue; // if it is from same heavy
326 }
327 // if the heavy's mother is not heavy, check the mother's label to know if it comes from inital or final parton shower
328 if (IsFromInitialShower(ancestorLabel[ig],motherID[i],motherType[i],motherLabel[i])) break;
329 if (IsFromFinalParton(ancestorLabel[ig],motherID[i],motherType[i],motherLabel[i])) break;
330 // if it is not the above case, something is strange
dbe3abbe 331 ReportStrangeness(motherID[i],motherType[i],motherLabel[i]);
259c3296 332 break;
333 }
334 if (ancestorLabel[ig] == -1){ // from hard scattering
335 HardScattering(kquark, motherID[i],motherType[i], motherLabel[i]);
336 }
337
338 } // end of found heavy quark loop
339
340
341 // check process type
342 Int_t processID = 0;
343 for (Int_t i = 0; i < fIsHeavy[iq]; i++){
344 AliDebug(1,Form("Mother ID= %d type= %d label= %d\n",motherID[i],motherType[i],motherLabel[i]));
345 }
346
347
348 Int_t nheavypair = Int_t(fIsHeavy[iq]/2.);
349 for (Int_t ipair = 0; ipair < nheavypair; ipair++){
350
351 Int_t id1 = ipair*2;
352 Int_t id2 = ipair*2 + 1;
353
354 if (motherType[id1] == 2 && motherType[id2] == 2){
dbe3abbe 355 if (motherLabel[id1] == motherLabel[id2]) processID = kGluonSplitting; // gluon spliting
259c3296 356 else processID = -9;
357 }
358 else if (motherType[id1] == -1 && motherType[id2] == -1) {
359 if (motherLabel[id1] == -1 && motherLabel[id2] == -1) {
dbe3abbe 360 if (motherID[id1] == fgkGluon) processID = kPairCreationFromg; // gluon fusion
361 else processID = kPairCreationFromq; // q-qbar pair creation
259c3296 362 }
363 else processID = -8;
364 }
365 else if (motherType[id1] == -1 || motherType[id2] == -1) {
366 if ((motherLabel[id1] == -1 || motherLabel[id2] == -1) && (motherLabel[id1]*motherLabel[id2] == -2 || motherLabel[id1]*motherLabel[id2] == -3)) {
dbe3abbe 367 if(motherID[id1]*motherID[id2] == kquark*fgkGluon) processID = kFlavourExitation; // flavour exitation
368 else processID = kLightQuarkShower;
259c3296 369 }
370 else processID = -7;
371 }
372 else if (motherType[id1] == -2 || motherType[id2] == -2) {
dbe3abbe 373 if (motherLabel[id1] == motherLabel[id2]) processID = kInitialPartonShower; // initial parton shower
259c3296 374 else processID = -6;
375
376 }
377 else processID = -5;
378
379 if (nheavypair >1) AliDebug(1,Form("Multi pair found : process ID = %d\n",processID));
dbe3abbe 380 else fHistComm[iq][0].fProcessID->Fill(processID);
259c3296 381 AliDebug(1,Form("Process ID = %d\n",processID));
382 } // end of # heavy quark pair loop
383
384}
385
386//__________________________________________
dbe3abbe 387void AliHFEmcQA::GetHadronKine(Int_t iTrack, const Int_t kquark)
388{
389 // decay electron kinematics
390
391 if (kquark != kCharm && kquark != kBeauty) {
392 AliDebug(1, "This task is only for heavy quark QA, return\n");
393 return;
394 }
395 Int_t iq = kquark - kCharm;
396
397 if (iTrack < 0) {
398 AliDebug(1, "Stack label is negative, return\n");
399 return;
400 }
401
402 TParticle* mcpart = fStack->Particle(iTrack);
403
404 Int_t iLabel = mcpart->GetFirstMother();
405 if (iLabel<0){
406 AliDebug(1, "Stack label is negative, return\n");
407 return;
408 }
409
410 TParticle *partCopy = mcpart;
411 Int_t pdgcode = mcpart->GetPdgCode();
412 Int_t pdgcodeCopy = pdgcode;
413
414 // if the mother is charmed hadron
75d81601 415 Bool_t isDirectCharm = kFALSE;
dbe3abbe 416 if ( int(abs(pdgcode)/100.) == kCharm || int(abs(pdgcode)/1000.) == kCharm ) {
417
418 // iterate until you find B hadron as a mother or become top ancester
419 for (Int_t i=1; i<fgkMaxIter; i++){
420
421 Int_t jLabel = mcpart->GetFirstMother();
422 if (jLabel == -1){
75d81601 423 isDirectCharm = kTRUE;
dbe3abbe 424 break; // if there is no ancester
425 }
426 if (jLabel < 0){ // safety protection
427 AliDebug(1, "Stack label is negative, return\n");
428 return;
429 }
430 // if there is an ancester
431 TParticle* mother = fStack->Particle(jLabel);
432 Int_t motherPDG = mother->GetPdgCode();
433
434 for (Int_t j=0; j<fNparents; j++){
435 if (abs(motherPDG)==fParentSelect[1][j]) return; // return if this hadron is originated from b
436 }
437
438 mcpart = mother;
439 } // end of iteration
440 } // end of if
75d81601 441 if((isDirectCharm == kTRUE && kquark == kCharm) || kquark == kBeauty) {
dbe3abbe 442 for (Int_t i=0; i<fNparents; i++){
443 if (abs(pdgcodeCopy)==fParentSelect[iq][i]){
444
445 // fill hadron kinematics
446 fHist[iq][kHadron][0].fPdgCode->Fill(pdgcodeCopy);
447 fHist[iq][kHadron][0].fPt->Fill(partCopy->Pt());
448 fHist[iq][kHadron][0].fY->Fill(GetRapidity(partCopy));
449 fHist[iq][kHadron][0].fEta->Fill(partCopy->Eta());
450
451 }
452 }
453 } // end of if
454}
455
456//__________________________________________
457void AliHFEmcQA::GetDecayedKine(Int_t iTrack, const Int_t kquark, Int_t kdecayed, Int_t icut, Bool_t isbarrel)
259c3296 458{
459 // decay electron kinematics
460
dbe3abbe 461 if (kquark != kCharm && kquark != kBeauty) {
259c3296 462 AliDebug(1, "This task is only for heavy quark QA, return\n");
463 return;
464 }
dbe3abbe 465 Int_t iq = kquark - kCharm;
259c3296 466
467 if (iTrack < 0) {
468 AliDebug(1, "Stack label is negative, return\n");
469 return;
470 }
471
472 TParticle* mcpart = fStack->Particle(iTrack);
473
474 if ( abs(mcpart->GetPdgCode()) != kdecayed ) return;
475 if ( isbarrel && TMath::Abs(mcpart->Eta()) > 0.9 ) return;
476
477 Int_t iLabel = mcpart->GetFirstMother();
478 if (iLabel<0){
479 AliDebug(1, "Stack label is negative, return\n");
480 return;
481 }
482
483 TParticle *partMother = fStack->Particle(iLabel);
dbe3abbe 484 TParticle *partMotherCopy = partMother;
259c3296 485 Int_t maPdgcode = partMother->GetPdgCode();
dbe3abbe 486 Int_t maPdgcodeCopy = maPdgcode;
259c3296 487
488 // get electron production vertex
489 TLorentzVector ePoint;
490 mcpart->ProductionVertex(ePoint);
491
259c3296 492 // if the mother is charmed hadron
dbe3abbe 493 Bool_t isMotherDirectCharm = kFALSE;
494 if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
259c3296 495
496 // iterate until you find B hadron as a mother or become top ancester
497 for (Int_t i=1; i<fgkMaxIter; i++){
498
499 Int_t jLabel = partMother->GetFirstMother();
500 if (jLabel == -1){
dbe3abbe 501 isMotherDirectCharm = kTRUE;
259c3296 502 break; // if there is no ancester
503 }
504 if (jLabel < 0){ // safety protection
505 AliDebug(1, "Stack label is negative, return\n");
506 return;
507 }
508
509 // if there is an ancester
510 TParticle* grandMa = fStack->Particle(jLabel);
511 Int_t grandMaPDG = grandMa->GetPdgCode();
512
fc0de2a0 513 for (Int_t j=0; j<fNparents; j++){
514 if (abs(grandMaPDG)==fParentSelect[1][j]){
259c3296 515
dbe3abbe 516 if (kquark == kCharm) return;
259c3296 517 // fill electron kinematics
dbe3abbe 518 fHist[iq][kElectron2nd][icut].fPdgCode->Fill(mcpart->GetPdgCode());
519 fHist[iq][kElectron2nd][icut].fPt->Fill(mcpart->Pt());
520 fHist[iq][kElectron2nd][icut].fY->Fill(GetRapidity(mcpart));
521 fHist[iq][kElectron2nd][icut].fEta->Fill(mcpart->Eta());
259c3296 522
523 // fill mother hadron kinematics
dbe3abbe 524 fHist[iq][kDeHadron][icut].fPdgCode->Fill(grandMaPDG);
525 fHist[iq][kDeHadron][icut].fPt->Fill(grandMa->Pt());
526 fHist[iq][kDeHadron][icut].fY->Fill(GetRapidity(grandMa));
527 fHist[iq][kDeHadron][icut].fEta->Fill(grandMa->Eta());
259c3296 528
529 // ratio between pT of electron and pT of mother B hadron
dbe3abbe 530 if(grandMa->Pt()) fHistComm[iq][icut].fDePtRatio->Fill(grandMa->Pt(),mcpart->Pt()/grandMa->Pt());
259c3296 531
532 // distance between electron production point and mother hadron production point
533 TLorentzVector debPoint;
534 grandMa->ProductionVertex(debPoint);
535 TLorentzVector dedistance = ePoint - debPoint;
dbe3abbe 536 fHistComm[iq][icut].fDeDistance->Fill(grandMa->Pt(),dedistance.P());
259c3296 537 return;
538 }
dc306130 539 }
259c3296 540
541 partMother = grandMa;
542 } // end of iteration
543 } // end of if
dbe3abbe 544 if((isMotherDirectCharm == kTRUE && kquark == kCharm) || kquark == kBeauty) {
259c3296 545 for (Int_t i=0; i<fNparents; i++){
546 if (abs(maPdgcodeCopy)==fParentSelect[iq][i]){
547
548 // fill electron kinematics
dbe3abbe 549 fHist[iq][kElectron][icut].fPdgCode->Fill(mcpart->GetPdgCode());
550 fHist[iq][kElectron][icut].fPt->Fill(mcpart->Pt());
551 fHist[iq][kElectron][icut].fY->Fill(GetRapidity(mcpart));
552 fHist[iq][kElectron][icut].fEta->Fill(mcpart->Eta());
259c3296 553
554 // fill mother hadron kinematics
dbe3abbe 555 fHist[iq][keHadron][icut].fPdgCode->Fill(maPdgcodeCopy);
556 fHist[iq][keHadron][icut].fPt->Fill(partMotherCopy->Pt());
557 fHist[iq][keHadron][icut].fY->Fill(GetRapidity(partMotherCopy));
558 fHist[iq][keHadron][icut].fEta->Fill(partMotherCopy->Eta());
259c3296 559
560 // ratio between pT of electron and pT of mother B hadron
dbe3abbe 561 if(partMotherCopy->Pt()) fHistComm[iq][icut].fePtRatio->Fill(partMotherCopy->Pt(),mcpart->Pt()/partMotherCopy->Pt());
259c3296 562
563 // distance between electron production point and mother hadron production point
564 TLorentzVector ebPoint;
565 partMotherCopy->ProductionVertex(ebPoint);
566 TLorentzVector edistance = ePoint - ebPoint;
dbe3abbe 567 fHistComm[iq][icut].feDistance->Fill(partMotherCopy->Pt(),edistance.P());
259c3296 568 }
569 }
570 } // end of if
571}
572
573
574//__________________________________________
75d81601 575void AliHFEmcQA::IdentifyMother(Int_t motherlabel, Int_t &motherpdg, Int_t &grandmotherlabel)
259c3296 576{
dbe3abbe 577 // find mother pdg code and label
578
75d81601 579 if (motherlabel < 0) {
259c3296 580 AliDebug(1, "Stack label is negative, return\n");
581 return;
582 }
75d81601 583 TParticle *heavysMother = fStack->Particle(motherlabel);
584 motherpdg = heavysMother->GetPdgCode();
585 grandmotherlabel = heavysMother->GetFirstMother();
586 AliDebug(1,Form("ancestor pdg code= %d\n",motherpdg));
259c3296 587}
588
589//__________________________________________
259c3296 590void AliHFEmcQA::HardScattering(const Int_t kquark, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
591{
dbe3abbe 592 // mothertype -1 means this heavy quark coming from hard vertex
593
259c3296 594 TParticle *afterinitialrad1 = fStack->Particle(4);
595 TParticle *afterinitialrad2 = fStack->Particle(5);
596
597 motherlabel = -1;
598
599 if (abs(afterinitialrad1->GetPdgCode()) == fgkGluon && abs(afterinitialrad2->GetPdgCode()) == fgkGluon){
600 AliDebug(1,"heavy from gluon gluon pair creation!\n");
601 mothertype = -1;
602 motherID = fgkGluon;
603 }
604 else if (abs(afterinitialrad1->GetPdgCode()) == kquark || abs(afterinitialrad2->GetPdgCode()) == kquark){ // one from Q and the other from g
605 AliDebug(1,"heavy from flavor exitation!\n");
606 mothertype = -1;
607 motherID = kquark;
608 }
609 else if (abs(afterinitialrad1->GetPdgCode()) == abs(afterinitialrad2->GetPdgCode())){
610 AliDebug(1,"heavy from q-qbar pair creation!\n");
611 mothertype = -1;
612 motherID = 1;
613 }
614 else {
615 AliDebug(1,"something strange!\n");
616 mothertype = -999;
617 motherlabel = -999;
618 motherID = -999;
619 }
620}
621
622//__________________________________________
259c3296 623Bool_t AliHFEmcQA::IsFromInitialShower(Int_t inputmotherlabel, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
624{
dbe3abbe 625 // mothertype -2 means this heavy quark coming from initial state
626
259c3296 627 if (inputmotherlabel==2 || inputmotherlabel==3){ // mother exist before initial state radiation
628 TParticle *heavysMother = fStack->Particle(inputmotherlabel);
629 motherID = heavysMother->GetPdgCode();
630 mothertype = -2; // there is mother before initial state radiation
631 motherlabel = inputmotherlabel;
632 AliDebug(1,"initial parton shower! \n");
633
634 return kTRUE;
635 }
636
637 return kFALSE;
638}
639
640//__________________________________________
259c3296 641Bool_t AliHFEmcQA::IsFromFinalParton(Int_t inputmotherlabel, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
642{
dbe3abbe 643 // mothertype 2 means this heavy quark coming from final state
644
259c3296 645 if (inputmotherlabel > 5){ // mother exist after hard scattering
646 TParticle *heavysMother = fStack->Particle(inputmotherlabel);
647 motherID = heavysMother->GetPdgCode();
648 mothertype = 2; //
649 motherlabel = inputmotherlabel;
650 AliDebug(1,Form("heavy quark from %d after hard scattering! \n",motherID));
651
652 return kTRUE;
653 }
654 return kFALSE;
655}
656
657//__________________________________________
dbe3abbe 658void AliHFEmcQA::ReportStrangeness(Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
259c3296 659{
dbe3abbe 660 // mark strange behavior
661
259c3296 662 mothertype = -888;
663 motherlabel = -888;
664 motherID = -888;
665 AliDebug(1,"something strange!\n");
666}
667
668//__________________________________________
dc306130 669Float_t AliHFEmcQA::GetRapidity(TParticle *part)
259c3296 670{
dbe3abbe 671 // return rapidity
672
259c3296 673 Float_t rapidity;
674 if(part->Energy() - part->Pz() == 0 || part->Energy() + part->Pz() == 0) rapidity=-999;
675 else rapidity = 0.5*(TMath::Log((part->Energy()+part->Pz()) / (part->Energy()-part->Pz())));
676 return rapidity;
677}