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