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