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