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