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