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