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