Update of the HFE package
[u/mrichter/AliRoot.git] / PWG3 / hfe / AliHFEmcQA.cxx
CommitLineData
259c3296 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
50685501 15//
16// QA class of Heavy Flavor quark and fragmeted/decayed particles
17// -Check kinematics of Heavy Quarks/hadrons, and decayed leptons
18// pT, rapidity
19// decay lepton kinematics w/wo acceptance
20// heavy hadron decay length, electron pT fraction carried from decay
21// -Check yield of Heavy Quarks/hadrons
22// Number of produced heavy quark
23// Number of produced hadron of given pdg code
24//
25//
26// Authors:
27// MinJung Kweon <minjung@physi.uni-heidelberg.de>
28//
259c3296 29
259c3296 30#include <TH2F.h>
259c3296 31#include <TH1F.h>
32#include <TH2F.h>
70da6c5a 33#include <TList.h>
259c3296 34#include <TParticle.h>
35
36#include <AliLog.h>
faee3b18 37#include <AliMCEvent.h>
9bcfd1ab 38#include <AliGenEventHeader.h>
0792aa82 39#include <AliAODMCParticle.h>
c2690925 40#include <AliStack.h>
259c3296 41
42#include "AliHFEmcQA.h"
70da6c5a 43#include "AliHFEtools.h"
c2690925 44#include "AliHFEcollection.h"
259c3296 45
259c3296 46
47ClassImp(AliHFEmcQA)
48
49//_______________________________________________________________________________________________
50AliHFEmcQA::AliHFEmcQA() :
bf892a6a 51 fMCEvent(NULL)
70da6c5a 52 ,fMCHeader(NULL)
53 ,fMCArray(NULL)
54 ,fQAhistos(NULL)
c2690925 55 ,fMCQACollection(NULL)
259c3296 56 ,fNparents(0)
57{
58 // Default constructor
e088e9f0 59 for(Int_t mom = 0; mom < 9; mom++){
60 fhD[mom] = NULL;
61 fhDLogbin[mom] = NULL;
62 }
bf892a6a 63 for(Int_t mom = 0; mom < 50; mom++){
64 fHeavyQuark[mom] = NULL;
65 }
66 for(Int_t mom = 0; mom < 2; mom++){
67 fIsHeavy[mom] = 0;
68 }
9250ffbf 69 memset(fElecBackgroundFactor, 0, sizeof(Double_t) * kElecBgSpecies * kBgPtBins);
70 memset(fBinLimit, 0, sizeof(Double_t) * (kBgPtBins+1));
259c3296 71}
72
73//_______________________________________________________________________________________________
74AliHFEmcQA::AliHFEmcQA(const AliHFEmcQA&p):
75 TObject(p)
faee3b18 76 ,fMCEvent(NULL)
70da6c5a 77 ,fMCHeader(NULL)
78 ,fMCArray(NULL)
79 ,fQAhistos(p.fQAhistos)
c2690925 80 ,fMCQACollection(p.fMCQACollection)
259c3296 81 ,fNparents(p.fNparents)
82{
83 // Copy constructor
e088e9f0 84 for(Int_t mom = 0; mom < 9; mom++){
85 fhD[mom] = NULL;
86 fhDLogbin[mom] = NULL;
87 }
bf892a6a 88 for(Int_t mom = 0; mom < 50; mom++){
89 fHeavyQuark[mom] = NULL;
90 }
91 for(Int_t mom = 0; mom < 2; mom++){
92 fIsHeavy[mom] = 0;
93 }
9250ffbf 94 memset(fElecBackgroundFactor, 0, sizeof(Double_t) * kElecBgSpecies * kBgPtBins);
95 memset(fBinLimit, 0, sizeof(Double_t) * (kBgPtBins+1));
259c3296 96}
97
98//_______________________________________________________________________________________________
99AliHFEmcQA&
100AliHFEmcQA::operator=(const AliHFEmcQA &)
101{
102 // Assignment operator
103
104 AliInfo("Not yet implemented.");
105 return *this;
106}
107
108//_______________________________________________________________________________________________
109AliHFEmcQA::~AliHFEmcQA()
110{
111 // Destructor
112
113 AliInfo("Analysis Done.");
114}
115
116//_______________________________________________________________________________________________
50685501 117void AliHFEmcQA::PostAnalyze() const
259c3296 118{
e3fc062d 119 //
120 // Post analysis
121 //
259c3296 122}
123
9250ffbf 124//_______________________________________________________________________________________________
125void AliHFEmcQA::SetBackgroundWeightFactor(Double_t *elecBackgroundFactor, Double_t *binLimit)
126{
e156c3bb 127 //
128 // copy background weighting factors into data member
129 //
9250ffbf 130 memcpy(fElecBackgroundFactor,elecBackgroundFactor,sizeof(Double_t) * kElecBgSpecies * kBgPtBins);
131 memcpy(fBinLimit,binLimit,sizeof(Double_t) * (kBgPtBins+1));
132/*
133 for(Int_t j=0;j < 6;j++){
134 for(Int_t i=0;i < 44;i++){
135 fElecBackgroundFactor[j][i] = elecBackgroundFactor[j][i];
136 }
137 }
138 for(Int_t i=0;i < 45;i++){
139 fBinLimit[i]=binLimit[i];
140 }
141 */
142}
143
259c3296 144//__________________________________________
e3fc062d 145void AliHFEmcQA::CreatDefaultHistograms(TList * const qaList)
146{
147 //
148 // make default histograms
149 //
150
151 if(!qaList) return;
152
153 fQAhistos = qaList;
154 fQAhistos->SetName("MCqa");
155
156 CreateHistograms(AliHFEmcQA::kCharm,0,"mcqa_"); // create histograms for charm
157 CreateHistograms(AliHFEmcQA::kBeauty,0,"mcqa_"); // create histograms for beauty
158 CreateHistograms(AliHFEmcQA::kOthers,0,"mcqa_"); // create histograms for beauty
159 CreateHistograms(AliHFEmcQA::kCharm,1,"mcqa_barrel_"); // create histograms for charm
160 CreateHistograms(AliHFEmcQA::kBeauty,1,"mcqa_barrel_"); // create histograms for beauty
161 CreateHistograms(AliHFEmcQA::kOthers,1,"mcqa_barrel_"); // create histograms for beauty
162 CreateHistograms(AliHFEmcQA::kCharm,2,"mcqa_unitY_"); // create histograms for charm
163 CreateHistograms(AliHFEmcQA::kBeauty,2,"mcqa_unitY_"); // create histograms for beauty
164 CreateHistograms(AliHFEmcQA::kOthers,2,"mcqa_unitY_"); // create histograms for beauty
e3fc062d 165
e97c2edf 166// check D spectra
ccc37cdc 167 TString kDspecies[9];
e97c2edf 168 kDspecies[0]="411";
169 kDspecies[1]="421";
170 kDspecies[2]="431";
171 kDspecies[3]="4122";
172 kDspecies[4]="4132";
173 kDspecies[5]="4232";
174 kDspecies[6]="4332";
ccc37cdc 175 kDspecies[7]="413";
176 kDspecies[8]="423";
177
178 const Double_t kPtbound[2] = {0.1, 20.}; //bin taken for considering inclusive e analysis binning
c2690925 179 Int_t iBin[2];
ccc37cdc 180 iBin[0] = 44; // bins in pt
c2690925 181 iBin[1] = 23; // bins in pt
ccc37cdc 182 Double_t* binEdges[1];
183 binEdges[0] = AliHFEtools::MakeLogarithmicBinning(iBin[0], kPtbound[0], kPtbound[1]);
e97c2edf 184
185 // bin size is chosen to consider ALICE D measurement
9250ffbf 186 Double_t xbins[15]={0,0.5,1.0,2.0,3.0,4.0,5.0,6.0,8.0,12,16,20,30,40,50};
e97c2edf 187 Double_t ybins[10]={-7.5,-1.0,-0.9,-0.8,-0.5,0.5,0.8,0.9,1.0,7.5};
188 TString hname;
ccc37cdc 189 for (Int_t iDmeson=0; iDmeson<9; iDmeson++){
e97c2edf 190 hname = "Dmeson"+kDspecies[iDmeson];
9250ffbf 191 fhD[iDmeson] = new TH2F(hname,hname+";p_{T} (GeV/c)",14,xbins,9,ybins);
ccc37cdc 192 hname = "DmesonLogbin"+kDspecies[iDmeson];
193 fhDLogbin[iDmeson] = new TH1F(hname,hname+";p_{T} (GeV/c)",iBin[0],binEdges[0]);
e97c2edf 194 if(fQAhistos) fQAhistos->Add(fhD[iDmeson]);
ccc37cdc 195 if(fQAhistos) fQAhistos->Add(fhDLogbin[iDmeson]);
e97c2edf 196 }
197
c2690925 198 const Double_t kPtRange[24] = {0.,0.3,0.4,0.5,0.6,0.8,1.,1.2,1.4,1.6,1.8,2.,2.2,2.4,2.6,2.8,3.,3.5,4.,5.,6.,7.,20.,30.}; // to cope with Ana's bin
199
200 fMCQACollection = new AliHFEcollection("TaskMCQA", "MC QA histos for meason pt spectra");
201 fMCQACollection->CreateTH1Farray("pionspectra", "pion yields: MC p_{t} ", iBin[1],kPtRange);
202 fMCQACollection->CreateTH1Farray("etaspectra", "eta yields: MC p_{t} ", iBin[1],kPtRange);
203 fMCQACollection->CreateTH1Farray("omegaspectra", "omega yields: MC p_{t} ", iBin[1],kPtRange);
204 fMCQACollection->CreateTH1Farray("phispectra", "phi yields: MC p_{t} ", iBin[1],kPtRange);
205 fMCQACollection->CreateTH1Farray("etapspectra", "etap yields: MC p_{t} ", iBin[1],kPtRange);
206 fMCQACollection->CreateTH1Farray("rhospectra", "rho yields: MC p_{t} ", iBin[1],kPtRange);
207
208 fMCQACollection->CreateTH1F("pionspectraLog", "pion yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
209 fMCQACollection->CreateTH1F("etaspectraLog", "eta yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
210 fMCQACollection->CreateTH1F("omegaspectraLog", "omega yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
211 fMCQACollection->CreateTH1F("phispectraLog", "phi yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
212 fMCQACollection->CreateTH1F("etapspectraLog", "etap yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
213 fMCQACollection->CreateTH1F("rhospectraLog", "rho yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
214
215 fMCQACollection->CreateTH1F("piondaughters", "pion yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
216 fMCQACollection->CreateTH1F("etadaughters", "eta yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
217 fMCQACollection->CreateTH1F("omegadaughters", "omega yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
218 fMCQACollection->CreateTH1F("phidaughters", "phi yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
219 fMCQACollection->CreateTH1F("etapdaughters", "etap yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
220 fMCQACollection->CreateTH1F("rhodaughters", "rho yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
221
222 fQAhistos->Add(fMCQACollection->GetList());
223
e3fc062d 224}
225
226//__________________________________________
dbe3abbe 227void AliHFEmcQA::CreateHistograms(const Int_t kquark, Int_t icut, TString hnopt)
259c3296 228{
229 // create histograms
230
faee3b18 231 if (!(kquark == kCharm || kquark == kBeauty || kquark == kOthers)) {
259c3296 232 AliDebug(1, "This task is only for heavy quark QA, return\n");
233 return;
234 }
dbe3abbe 235 Int_t iq = kquark - kCharm;
259c3296 236
237 TString kqTypeLabel[fgkqType];
dbe3abbe 238 if (kquark == kCharm){
239 kqTypeLabel[kQuark]="c";
240 kqTypeLabel[kantiQuark]="cbar";
241 kqTypeLabel[kHadron]="cHadron";
242 kqTypeLabel[keHadron]="ceHadron";
243 kqTypeLabel[kDeHadron]="nullHadron";
244 kqTypeLabel[kElectron]="ce";
245 kqTypeLabel[kElectron2nd]="nulle";
246 } else if (kquark == kBeauty){
247 kqTypeLabel[kQuark]="b";
248 kqTypeLabel[kantiQuark]="bbar";
249 kqTypeLabel[kHadron]="bHadron";
250 kqTypeLabel[keHadron]="beHadron";
251 kqTypeLabel[kDeHadron]="bDeHadron";
252 kqTypeLabel[kElectron]="be";
253 kqTypeLabel[kElectron2nd]="bce";
faee3b18 254 } else if (kquark == kOthers){
255 kqTypeLabel[kGamma-4]="gammae";
256 kqTypeLabel[kPi0-4]="pi0e";
257 kqTypeLabel[kElse-4]="elsee";
258 kqTypeLabel[kMisID-4]="miside";
259c3296 259 }
3a72645a 260
261 const Double_t kPtbound[2] = {0.1, 20.}; //bin taken for considering inclusive e analysis binning
faee3b18 262 Int_t iBin[1];
3a72645a 263 iBin[0] = 44; // bins in pt
faee3b18 264 Double_t* binEdges[1];
265 binEdges[0] = AliHFEtools::MakeLogarithmicBinning(iBin[0], kPtbound[0], kPtbound[1]);
c2690925 266
267
9250ffbf 268 Double_t xcorrbin[501];
269 for (int icorrbin = 0; icorrbin< 501; icorrbin++){
c2690925 270 xcorrbin[icorrbin]=icorrbin*0.1;
271 }
272
259c3296 273 TString hname;
faee3b18 274 if(kquark == kOthers){
275 for (Int_t iqType = 0; iqType < 4; iqType++ ){
276 hname = hnopt+"Pt_"+kqTypeLabel[iqType];
e3ae862b 277 //fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",iBin[0],binEdges[0]); //mj to compare with FONLL
278 fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",60,0.25,30.25);
faee3b18 279 hname = hnopt+"Y_"+kqTypeLabel[iqType];
280 fHist[iq][iqType][icut].fY = new TH1F(hname,hname,150,-7.5,7.5);
281 hname = hnopt+"Eta_"+kqTypeLabel[iqType];
282 fHist[iq][iqType][icut].fEta = new TH1F(hname,hname,150,-7.5,7.5);
283 // Fill List
284 if(fQAhistos) fHist[iq][iqType][icut].FillList(fQAhistos);
285 }
286 return;
287 }
259c3296 288 for (Int_t iqType = 0; iqType < fgkqType; iqType++ ){
dbe3abbe 289 if (iqType < keHadron && icut > 0) continue; // don't duplicate histogram for quark and hadron
259c3296 290 hname = hnopt+"PdgCode_"+kqTypeLabel[iqType];
dbe3abbe 291 fHist[iq][iqType][icut].fPdgCode = new TH1F(hname,hname,20001,-10000.5,10000.5);
259c3296 292 hname = hnopt+"Pt_"+kqTypeLabel[iqType];
c2690925 293 fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",iBin[0],binEdges[0]);
294 //fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",60,0.,30.); // mj to compare with FONLL
295 //fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",60,0.25,30.25); // mj to compare with FONLL
259c3296 296 hname = hnopt+"Y_"+kqTypeLabel[iqType];
dbe3abbe 297 fHist[iq][iqType][icut].fY = new TH1F(hname,hname,150,-7.5,7.5);
259c3296 298 hname = hnopt+"Eta_"+kqTypeLabel[iqType];
dbe3abbe 299 fHist[iq][iqType][icut].fEta = new TH1F(hname,hname,150,-7.5,7.5);
70da6c5a 300 // Fill List
301 if(fQAhistos) fHist[iq][iqType][icut].FillList(fQAhistos);
259c3296 302 }
303
dbe3abbe 304 if (icut == 0){
305 hname = hnopt+"Nq_"+kqTypeLabel[kQuark];
e3ae862b 306 fHistComm[iq][icut].fNq = new TH1F(hname,hname,50,-0.5,49.5);
dbe3abbe 307 hname = hnopt+"ProcessID_"+kqTypeLabel[kQuark];
308 fHistComm[iq][icut].fProcessID = new TH1F(hname,hname,21,-10.5,10.5);
309 }
310 hname = hnopt+"ePtRatio_"+kqTypeLabel[kQuark];
e088e9f0 311 fHistComm[iq][icut].fePtRatio = new TH2F(hname,hname+";p_{T} (GeV/c);momentum fraction",200,0,20,100,0,1);
ccc37cdc 312 hname = hnopt+"PtCorr_"+kqTypeLabel[kQuark];
9250ffbf 313 fHistComm[iq][icut].fPtCorr = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",500,xcorrbin,44,binEdges[0]);
c2690925 314 //fHistComm[iq][icut].fPtCorr = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",200,0,20,200,0,20);
315 hname = hnopt+"PtCorrDp_"+kqTypeLabel[kQuark];
9250ffbf 316 fHistComm[iq][icut].fPtCorrDp= new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",500,xcorrbin,44,binEdges[0]);
c2690925 317 //fHistComm[iq][icut].fPtCorrDp= new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",200,0,20,200,0,20);
318 hname = hnopt+"PtCorrD0_"+kqTypeLabel[kQuark];
9250ffbf 319 fHistComm[iq][icut].fPtCorrD0 = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",500,xcorrbin,44,binEdges[0]);
c2690925 320 //fHistComm[iq][icut].fPtCorrD0 = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",200,0,20,200,0,20);
321 hname = hnopt+"PtCorrDrest_"+kqTypeLabel[kQuark];
9250ffbf 322 fHistComm[iq][icut].fPtCorrDrest = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",500,xcorrbin,44,binEdges[0]);
c2690925 323 //fHistComm[iq][icut].fPtCorrDrest = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",200,0,20,200,0,20);
dbe3abbe 324 hname = hnopt+"DePtRatio_"+kqTypeLabel[kQuark];
e088e9f0 325 fHistComm[iq][icut].fDePtRatio = new TH2F(hname,hname+";p_{T} (GeV/c);momentum fraction",100,0,20,100,0,1);
dbe3abbe 326 hname = hnopt+"eDistance_"+kqTypeLabel[kQuark];
e088e9f0 327 fHistComm[iq][icut].feDistance= new TH2F(hname,hname+";p_{T} (GeV/c);distance (cm)",100,0,20,200,0,2);
dbe3abbe 328 hname = hnopt+"DeDistance_"+kqTypeLabel[kQuark];
e088e9f0 329 fHistComm[iq][icut].fDeDistance= new TH2F(hname,hname+";p_{T} (GeV/c);distance (cm)",100,0,20,200,0,2);
70da6c5a 330 if(fQAhistos) fHistComm[iq][icut].FillList(fQAhistos);
259c3296 331}
332
333//__________________________________________
334void AliHFEmcQA::Init()
335{
336 // called at begining every event
337
338 for (Int_t i=0; i<2; i++){
339 fIsHeavy[i] = 0;
340 }
341
342 fNparents = 7;
343
344 fParentSelect[0][0] = 411;
345 fParentSelect[0][1] = 421;
346 fParentSelect[0][2] = 431;
347 fParentSelect[0][3] = 4122;
348 fParentSelect[0][4] = 4132;
349 fParentSelect[0][5] = 4232;
350 fParentSelect[0][6] = 4332;
351
352 fParentSelect[1][0] = 511;
353 fParentSelect[1][1] = 521;
354 fParentSelect[1][2] = 531;
355 fParentSelect[1][3] = 5122;
356 fParentSelect[1][4] = 5132;
357 fParentSelect[1][5] = 5232;
358 fParentSelect[1][6] = 5332;
359
360}
361
9250ffbf 362//__________________________________________
c2690925 363void AliHFEmcQA::GetMesonKine()
364{
365 //
366 // get meson pt spectra
367 //
368
369 AliVParticle *mctrack2 = NULL;
370 AliMCParticle *mctrack0 = NULL;
371 AliVParticle *mctrackdaugt= NULL;
372 AliMCParticle *mctrackd= NULL;
373 Int_t id1=0, id2=0;
374
375 for(Int_t imc = 0; imc <fMCEvent->GetNumberOfPrimaries(); imc++){
376 if(!(mctrack2 = fMCEvent->GetTrack(imc))) continue;
377 TParticle* mcpart0 = fMCEvent->Stack()->Particle(imc);
378 if(!mcpart0) continue;
379 mctrack0 = dynamic_cast<AliMCParticle *>(mctrack2);
f9f097c0 380 if(!mctrack0) continue;
c2690925 381 if(abs(mctrack0->PdgCode()) == 111) // pi0
382 {
383 if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
384 fMCQACollection->Fill("pionspectra",mctrack0->Pt());
385 fMCQACollection->Fill("pionspectraLog",mctrack0->Pt());
386 }
387 id1=mctrack0->GetFirstDaughter();
388 id2=mctrack0->GetLastDaughter();
389 if(!((id2-id1)==2)) continue;
390 for(int idx=id1; idx<=id2; idx++){
391 if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
9250ffbf 392 if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
c2690925 393 if(abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
394 fMCQACollection->Fill("piondaughters",mctrackd->Pt());
395 }
396 }
397 else if(abs(mctrack0->PdgCode()) == 221) // eta
398 {
399 if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
400 fMCQACollection->Fill("etaspectra",mctrack0->Pt());
401 fMCQACollection->Fill("etaspectraLog",mctrack0->Pt());
402 }
403 id1=mctrack0->GetFirstDaughter();
404 id2=mctrack0->GetLastDaughter();
405 if(!((id2-id1)==2||(id2-id1)==3)) continue;
406 for(int idx=id1; idx<=id2; idx++){
407 if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
9250ffbf 408 if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
c2690925 409 if(abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
410 fMCQACollection->Fill("etadaughters",mctrackd->Pt());
411 }
412 }
413 else if(abs(mctrack0->PdgCode()) == 223) // omega
414 {
415 if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
416 fMCQACollection->Fill("omegaspectra",mctrack0->Pt());
417 fMCQACollection->Fill("omegaspectraLog",mctrack0->Pt());
418 }
419 id1=mctrack0->GetFirstDaughter();
420 id2=mctrack0->GetLastDaughter();
421 if(!((id2-id1)==1||(id2-id1)==2)) continue;
422 for(int idx=id1; idx<=id2; idx++){
423 if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
9250ffbf 424 if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
c2690925 425 if(abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
426 fMCQACollection->Fill("omegadaughters",mctrackd->Pt());
427 }
428 }
429 else if(abs(mctrack0->PdgCode()) == 333) // phi
430 {
431 if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
432 fMCQACollection->Fill("phispectra",mctrack0->Pt());
433 fMCQACollection->Fill("phispectraLog",mctrack0->Pt());
434 }
435 id1=mctrack0->GetFirstDaughter();
436 id2=mctrack0->GetLastDaughter();
437 if(!((id2-id1)==1)) continue;
438 for(int idx=id1; idx<=id2; idx++){
439 if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
9250ffbf 440 if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
c2690925 441 if(abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
442 fMCQACollection->Fill("phidaughters",mctrackd->Pt());
443 }
444 }
445 else if(abs(mctrack0->PdgCode()) == 331) // eta prime
446 {
447 if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
448 fMCQACollection->Fill("etapspectra",mctrack0->Pt());
449 fMCQACollection->Fill("etapspectraLog",mctrack0->Pt());
450 }
451 id1=mctrack0->GetFirstDaughter();
452 id2=mctrack0->GetLastDaughter();
453 if(!((id2-id1)==2||(id2-id1)==3)) continue;
454 for(int idx=id1; idx<=id2; idx++){
455 if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
9250ffbf 456 if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
c2690925 457 if(abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
458 fMCQACollection->Fill("etapdaughters",mctrackd->Pt());
459 }
460 }
461 else if(abs(mctrack0->PdgCode()) == 113) // rho
462 {
463 if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
464 fMCQACollection->Fill("rhospectra",mctrack0->Pt());
465 fMCQACollection->Fill("rhospectraLog",mctrack0->Pt());
466 }
467 id1=mctrack0->GetFirstDaughter();
468 id2=mctrack0->GetLastDaughter();
469 if(!((id2-id1)==1)) continue;
470 for(int idx=id1; idx<=id2; idx++){
471 if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
9250ffbf 472 if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
c2690925 473 if(abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
474 fMCQACollection->Fill("rhodaughters",mctrackd->Pt());
475 }
476 }
477 }
478
479}
259c3296 480//__________________________________________
722347d8 481void AliHFEmcQA::GetQuarkKine(TParticle *part, Int_t iTrack, const Int_t kquark)
259c3296 482{
483 // get heavy quark kinematics
484
dbe3abbe 485 if (kquark != kCharm && kquark != kBeauty) {
259c3296 486 AliDebug(1, "This task is only for heavy quark QA, return\n");
487 return;
488 }
dbe3abbe 489 Int_t iq = kquark - kCharm;
259c3296 490
faee3b18 491 if (iTrack < 0 || !part) {
492 AliDebug(1, "Stack label is negative or no mcparticle, return\n");
259c3296 493 return;
494 }
495
faee3b18 496 AliMCParticle *mctrack = NULL;
259c3296 497 Int_t partPdgcode = TMath::Abs(part->GetPdgCode());
498
499 // select heavy hadron or not fragmented heavy quark
500 if ( int(partPdgcode/100.)==kquark || int(partPdgcode/1000.)==kquark || (partPdgcode==kquark && (part->GetNDaughters()==0 && iTrack>5)) ){
501
502 TParticle *partMother;
503 Int_t iLabel;
504
505 if (partPdgcode == kquark){ // in case of not fragmented heavy quark
506 partMother = part;
507 iLabel = iTrack;
508 } else{ // in case of heavy hadron, start to search for mother heavy parton
509 iLabel = part->GetFirstMother();
faee3b18 510 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return;
511 if (iLabel>-1) { partMother = mctrack->Particle(); }
259c3296 512 else {
513 AliDebug(1, "Stack label is negative, return\n");
514 return;
515 }
516 }
517
518 // heavy parton selection as a mother of heavy hadron
519 // if the heavy particle comes from string which is denoted as particle status 12|12|12...12|11,[PYTHIA p.60]
520 // in this case, the mother of heavy particle can be one of the fragmented parton of the string
521 // should I make a condition that partMother should be quark or diquark? -> not necessary
522 if ( abs(partMother->GetPdgCode()) == kquark || (partMother->GetStatusCode() == 11) ){
523 //if ( abs(partMother->GetPdgCode()) == kquark || (partMother->GetStatusCode() == 11 || partMother->GetStatusCode() == 12) ){
524
525 if ( abs(partMother->GetPdgCode()) != kquark ){
526 // search fragmented partons in the same string
dbe3abbe 527 Bool_t isSameString = kTRUE;
259c3296 528 for (Int_t i=1; i<fgkMaxIter; i++){
529 iLabel = iLabel - 1;
faee3b18 530 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return;
531 if (iLabel>-1) { partMother = mctrack->Particle(); }
259c3296 532 else {
533 AliDebug(1, "Stack label is negative, return\n");
534 return;
535 }
536 if ( abs(partMother->GetPdgCode()) == kquark ) break;
dbe3abbe 537 if ( partMother->GetStatusCode() != 12 ) isSameString = kFALSE;
538 if (!isSameString) return;
259c3296 539 }
540 }
541 AliDebug(1, "Can not find heavy parton of this heavy hadron in the string, return\n");
542 if (abs(partMother->GetPdgCode()) != kquark) return;
543
544 if (fIsHeavy[iq] >= 50) return;
545 fHeavyQuark[fIsHeavy[iq]] = partMother;
546 fIsHeavy[iq]++;
547
548 // fill kinematics for heavy parton
549 if (partMother->GetPdgCode() > 0) { // quark
dbe3abbe 550 fHist[iq][kQuark][0].fPdgCode->Fill(partMother->GetPdgCode());
551 fHist[iq][kQuark][0].fPt->Fill(partMother->Pt());
70da6c5a 552 fHist[iq][kQuark][0].fY->Fill(AliHFEtools::GetRapidity(partMother));
dbe3abbe 553 fHist[iq][kQuark][0].fEta->Fill(partMother->Eta());
259c3296 554 } else{ // antiquark
dbe3abbe 555 fHist[iq][kantiQuark][0].fPdgCode->Fill(partMother->GetPdgCode());
556 fHist[iq][kantiQuark][0].fPt->Fill(partMother->Pt());
70da6c5a 557 fHist[iq][kantiQuark][0].fY->Fill(AliHFEtools::GetRapidity(partMother));
dbe3abbe 558 fHist[iq][kantiQuark][0].fEta->Fill(partMother->Eta());
259c3296 559 }
560
561 } // end of heavy parton slection loop
562
563 } // end of heavy hadron or quark selection
564
565}
566
567//__________________________________________
568void AliHFEmcQA::EndOfEventAna(const Int_t kquark)
569{
570 // end of event analysis
571
dbe3abbe 572 if (kquark != kCharm && kquark != kBeauty) {
259c3296 573 AliDebug(1, "This task is only for heavy quark QA, return\n");
574 return;
575 }
dbe3abbe 576 Int_t iq = kquark - kCharm;
259c3296 577
578
579 // # of heavy quark per event
580 AliDebug(1,Form("Number of heavy quark in this event = %d \n",fIsHeavy[iq]));
dbe3abbe 581 fHistComm[iq][0].fNq->Fill(fIsHeavy[iq]);
259c3296 582
583 Int_t motherID[fgkMaxGener];
584 Int_t motherType[fgkMaxGener];
585 Int_t motherLabel[fgkMaxGener];
586 Int_t ancestorPdg[fgkMaxGener];
587 Int_t ancestorLabel[fgkMaxGener];
588
589 for (Int_t i = 0; i < fgkMaxGener; i++){ // initialization
590 motherID[i] = 0;
591 motherType[i] = 0;
592 motherLabel[i] = 0;
593 ancestorPdg[i] = 0;
594 ancestorLabel[i] = 0;
595 }
596
3a72645a 597
259c3296 598 // check history of found heavy quarks
599 for (Int_t i = 0; i < fIsHeavy[iq]; i++){
600
3a72645a 601 if(!fHeavyQuark[i]) return;
602
259c3296 603 ancestorLabel[0] = i;
604 ancestorPdg[0] = fHeavyQuark[i]->GetPdgCode();
605 ancestorLabel[1] = fHeavyQuark[i]->GetFirstMother();
606
607 AliDebug(1,Form("pdg code= %d\n",ancestorPdg[0]));
608 AliDebug(1,Form("ancestor label= %d\n",ancestorLabel[1]));
609
610 Int_t ig = 1;
611 while (ancestorLabel[ig] != -1){
612 // in case there is mother, get mother's pdg code and grandmother's label
613 IdentifyMother(ancestorLabel[ig], ancestorPdg[ig], ancestorLabel[ig+1]);
614 // if mother is still heavy, find again mother's ancestor
615 if (ancestorPdg[ig-1] == ancestorPdg[ig]) {
616 ig++;
617 continue; // if it is from same heavy
618 }
619 // if the heavy's mother is not heavy, check the mother's label to know if it comes from inital or final parton shower
620 if (IsFromInitialShower(ancestorLabel[ig],motherID[i],motherType[i],motherLabel[i])) break;
621 if (IsFromFinalParton(ancestorLabel[ig],motherID[i],motherType[i],motherLabel[i])) break;
622 // if it is not the above case, something is strange
dbe3abbe 623 ReportStrangeness(motherID[i],motherType[i],motherLabel[i]);
259c3296 624 break;
625 }
626 if (ancestorLabel[ig] == -1){ // from hard scattering
627 HardScattering(kquark, motherID[i],motherType[i], motherLabel[i]);
628 }
629
630 } // end of found heavy quark loop
631
632
633 // check process type
634 Int_t processID = 0;
635 for (Int_t i = 0; i < fIsHeavy[iq]; i++){
636 AliDebug(1,Form("Mother ID= %d type= %d label= %d\n",motherID[i],motherType[i],motherLabel[i]));
637 }
638
639
640 Int_t nheavypair = Int_t(fIsHeavy[iq]/2.);
641 for (Int_t ipair = 0; ipair < nheavypair; ipair++){
642
643 Int_t id1 = ipair*2;
644 Int_t id2 = ipair*2 + 1;
645
646 if (motherType[id1] == 2 && motherType[id2] == 2){
dbe3abbe 647 if (motherLabel[id1] == motherLabel[id2]) processID = kGluonSplitting; // gluon spliting
259c3296 648 else processID = -9;
649 }
650 else if (motherType[id1] == -1 && motherType[id2] == -1) {
651 if (motherLabel[id1] == -1 && motherLabel[id2] == -1) {
dbe3abbe 652 if (motherID[id1] == fgkGluon) processID = kPairCreationFromg; // gluon fusion
653 else processID = kPairCreationFromq; // q-qbar pair creation
259c3296 654 }
655 else processID = -8;
656 }
657 else if (motherType[id1] == -1 || motherType[id2] == -1) {
658 if ((motherLabel[id1] == -1 || motherLabel[id2] == -1) && (motherLabel[id1]*motherLabel[id2] == -2 || motherLabel[id1]*motherLabel[id2] == -3)) {
dbe3abbe 659 if(motherID[id1]*motherID[id2] == kquark*fgkGluon) processID = kFlavourExitation; // flavour exitation
660 else processID = kLightQuarkShower;
259c3296 661 }
662 else processID = -7;
663 }
664 else if (motherType[id1] == -2 || motherType[id2] == -2) {
dbe3abbe 665 if (motherLabel[id1] == motherLabel[id2]) processID = kInitialPartonShower; // initial parton shower
259c3296 666 else processID = -6;
667
668 }
669 else processID = -5;
670
671 if (nheavypair >1) AliDebug(1,Form("Multi pair found : process ID = %d\n",processID));
dbe3abbe 672 else fHistComm[iq][0].fProcessID->Fill(processID);
259c3296 673 AliDebug(1,Form("Process ID = %d\n",processID));
674 } // end of # heavy quark pair loop
675
676}
677
678//__________________________________________
722347d8 679void AliHFEmcQA::GetHadronKine(TParticle* mcpart, const Int_t kquark)
dbe3abbe 680{
681 // decay electron kinematics
682
683 if (kquark != kCharm && kquark != kBeauty) {
684 AliDebug(1, "This task is only for heavy quark QA, return\n");
685 return;
686 }
687 Int_t iq = kquark - kCharm;
688
faee3b18 689 if(!mcpart){
690 AliDebug(1, "no mc particle, return\n");
691 return;
692 }
dbe3abbe 693
694 Int_t iLabel = mcpart->GetFirstMother();
695 if (iLabel<0){
696 AliDebug(1, "Stack label is negative, return\n");
697 return;
698 }
699
700 TParticle *partCopy = mcpart;
701 Int_t pdgcode = mcpart->GetPdgCode();
702 Int_t pdgcodeCopy = pdgcode;
703
faee3b18 704 AliMCParticle *mctrack = NULL;
705
dbe3abbe 706 // if the mother is charmed hadron
75d81601 707 Bool_t isDirectCharm = kFALSE;
dbe3abbe 708 if ( int(abs(pdgcode)/100.) == kCharm || int(abs(pdgcode)/1000.) == kCharm ) {
709
710 // iterate until you find B hadron as a mother or become top ancester
711 for (Int_t i=1; i<fgkMaxIter; i++){
712
713 Int_t jLabel = mcpart->GetFirstMother();
714 if (jLabel == -1){
75d81601 715 isDirectCharm = kTRUE;
dbe3abbe 716 break; // if there is no ancester
717 }
718 if (jLabel < 0){ // safety protection
719 AliDebug(1, "Stack label is negative, return\n");
720 return;
721 }
722 // if there is an ancester
faee3b18 723 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return;
724 TParticle* mother = mctrack->Particle();
dbe3abbe 725 Int_t motherPDG = mother->GetPdgCode();
726
727 for (Int_t j=0; j<fNparents; j++){
728 if (abs(motherPDG)==fParentSelect[1][j]) return; // return if this hadron is originated from b
729 }
730
731 mcpart = mother;
732 } // end of iteration
733 } // end of if
75d81601 734 if((isDirectCharm == kTRUE && kquark == kCharm) || kquark == kBeauty) {
dbe3abbe 735 for (Int_t i=0; i<fNparents; i++){
736 if (abs(pdgcodeCopy)==fParentSelect[iq][i]){
737
738 // fill hadron kinematics
739 fHist[iq][kHadron][0].fPdgCode->Fill(pdgcodeCopy);
740 fHist[iq][kHadron][0].fPt->Fill(partCopy->Pt());
70da6c5a 741 fHist[iq][kHadron][0].fY->Fill(AliHFEtools::GetRapidity(partCopy));
dbe3abbe 742 fHist[iq][kHadron][0].fEta->Fill(partCopy->Eta());
743
ccc37cdc 744 if(iq==0) {
745 fhD[i]->Fill(partCopy->Pt(),AliHFEtools::GetRapidity(partCopy));
746 if(TMath::Abs(AliHFEtools::GetRapidity(partCopy))<0.5) fhDLogbin[i]->Fill(partCopy->Pt());
747 }
dbe3abbe 748 }
749 }
ccc37cdc 750 // I also want to store D* info to compare with D* measurement
751 if (abs(pdgcodeCopy)==413 && iq==0) { //D*+
752 fhD[7]->Fill(partCopy->Pt(),AliHFEtools::GetRapidity(partCopy));
753 if(TMath::Abs(AliHFEtools::GetRapidity(partCopy))<0.5) fhDLogbin[7]->Fill(partCopy->Pt());
754 }
755 if (abs(pdgcodeCopy)==423 && iq==0) { //D*0
756 fhD[8]->Fill(partCopy->Pt(),AliHFEtools::GetRapidity(partCopy));
757 if(TMath::Abs(AliHFEtools::GetRapidity(partCopy))<0.5) fhDLogbin[8]->Fill(partCopy->Pt());
758 }
dbe3abbe 759 } // end of if
760}
761
762//__________________________________________
722347d8 763void AliHFEmcQA::GetDecayedKine(TParticle* mcpart, const Int_t kquark, Int_t kdecayed, Int_t icut)
259c3296 764{
765 // decay electron kinematics
766
faee3b18 767 if (!(kquark == kCharm || kquark == kBeauty || kquark == kOthers)){
259c3296 768 AliDebug(1, "This task is only for heavy quark QA, return\n");
769 return;
770 }
dbe3abbe 771 Int_t iq = kquark - kCharm;
50685501 772 Bool_t isFinalOpenCharm = kFALSE;
259c3296 773
faee3b18 774 if(!mcpart){
775 AliDebug(1, "no mcparticle, return\n");
776 return;
259c3296 777 }
778
faee3b18 779 if(kquark==kOthers){
780 Int_t esource = -1;
781 if ( abs(mcpart->GetPdgCode()) != kdecayed ) esource = kMisID-4;
782 else esource =GetSource(mcpart)-4; // return for the cases kGamma=4, kPi0=5, kElse=6
783 if(esource==0|| esource==1 || esource==2 || esource==3){
784 fHist[iq][esource][icut].fPt->Fill(mcpart->Pt());
785 fHist[iq][esource][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart));
786 fHist[iq][esource][icut].fEta->Fill(mcpart->Eta());
787 return;
788 }
789 else {
790 AliDebug(1, "e source is out of defined ranges, return\n");
791 return;
792 }
793 }
259c3296 794
795 if ( abs(mcpart->GetPdgCode()) != kdecayed ) return;
259c3296 796
797 Int_t iLabel = mcpart->GetFirstMother();
798 if (iLabel<0){
799 AliDebug(1, "Stack label is negative, return\n");
800 return;
801 }
802
faee3b18 803 AliMCParticle *mctrack = NULL;
804 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return;
805 TParticle *partMother = mctrack->Particle();
dbe3abbe 806 TParticle *partMotherCopy = partMother;
259c3296 807 Int_t maPdgcode = partMother->GetPdgCode();
dbe3abbe 808 Int_t maPdgcodeCopy = maPdgcode;
259c3296 809
9bcfd1ab 810 // get mc primary vertex
faee3b18 811 /*
9bcfd1ab 812 TArrayF mcPrimVtx;
813 if(fMCHeader) fMCHeader->PrimaryVertex(mcPrimVtx);
814
259c3296 815 // get electron production vertex
816 TLorentzVector ePoint;
817 mcpart->ProductionVertex(ePoint);
818
9bcfd1ab 819 // calculated production vertex to primary vertex (in xy plane)
820 Float_t decayLxy = TMath::Sqrt((mcPrimVtx[0]-ePoint[0])*(mcPrimVtx[0]-ePoint[0])+(mcPrimVtx[1]-ePoint[1])*(mcPrimVtx[1]-ePoint[1]));
faee3b18 821 */
822 Float_t decayLxy = 0;
9bcfd1ab 823
259c3296 824 // if the mother is charmed hadron
dbe3abbe 825 Bool_t isMotherDirectCharm = kFALSE;
826 if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
259c3296 827
0792aa82 828 for (Int_t i=0; i<fNparents; i++){
829 if (abs(maPdgcode)==fParentSelect[0][i]){
50685501 830 isFinalOpenCharm = kTRUE;
0792aa82 831 }
832 }
50685501 833 if (!isFinalOpenCharm) return ;
0792aa82 834
259c3296 835 // iterate until you find B hadron as a mother or become top ancester
836 for (Int_t i=1; i<fgkMaxIter; i++){
837
838 Int_t jLabel = partMother->GetFirstMother();
839 if (jLabel == -1){
dbe3abbe 840 isMotherDirectCharm = kTRUE;
259c3296 841 break; // if there is no ancester
842 }
843 if (jLabel < 0){ // safety protection
844 AliDebug(1, "Stack label is negative, return\n");
845 return;
846 }
847
848 // if there is an ancester
faee3b18 849 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return;
850 TParticle* grandMa = mctrack->Particle();
259c3296 851 Int_t grandMaPDG = grandMa->GetPdgCode();
852
fc0de2a0 853 for (Int_t j=0; j<fNparents; j++){
854 if (abs(grandMaPDG)==fParentSelect[1][j]){
259c3296 855
dbe3abbe 856 if (kquark == kCharm) return;
259c3296 857 // fill electron kinematics
dbe3abbe 858 fHist[iq][kElectron2nd][icut].fPdgCode->Fill(mcpart->GetPdgCode());
859 fHist[iq][kElectron2nd][icut].fPt->Fill(mcpart->Pt());
70da6c5a 860 fHist[iq][kElectron2nd][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart));
dbe3abbe 861 fHist[iq][kElectron2nd][icut].fEta->Fill(mcpart->Eta());
259c3296 862
863 // fill mother hadron kinematics
dbe3abbe 864 fHist[iq][kDeHadron][icut].fPdgCode->Fill(grandMaPDG);
865 fHist[iq][kDeHadron][icut].fPt->Fill(grandMa->Pt());
70da6c5a 866 fHist[iq][kDeHadron][icut].fY->Fill(AliHFEtools::GetRapidity(grandMa));
dbe3abbe 867 fHist[iq][kDeHadron][icut].fEta->Fill(grandMa->Eta());
259c3296 868
869 // ratio between pT of electron and pT of mother B hadron
dbe3abbe 870 if(grandMa->Pt()) fHistComm[iq][icut].fDePtRatio->Fill(grandMa->Pt(),mcpart->Pt()/grandMa->Pt());
259c3296 871
9bcfd1ab 872 // distance between electron production point and primary vertex
873 fHistComm[iq][icut].fDeDistance->Fill(grandMa->Pt(),decayLxy);
259c3296 874 return;
875 }
dc306130 876 }
259c3296 877
878 partMother = grandMa;
879 } // end of iteration
880 } // end of if
dbe3abbe 881 if((isMotherDirectCharm == kTRUE && kquark == kCharm) || kquark == kBeauty) {
259c3296 882 for (Int_t i=0; i<fNparents; i++){
883 if (abs(maPdgcodeCopy)==fParentSelect[iq][i]){
884
e97c2edf 885//mj weighting to consider measured spectra!!!
886 Double_t mpt=partMotherCopy->Pt();
ccc37cdc 887 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
888 //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 889 // fill electron kinematics
e97c2edf 890 if(iq==0){
891 fHist[iq][kElectron2nd][icut].fPdgCode->Fill(mcpart->GetPdgCode(),wfactor);
892 fHist[iq][kElectron2nd][icut].fPt->Fill(mcpart->Pt(),wfactor);
893 fHist[iq][kElectron2nd][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart),wfactor);
894 fHist[iq][kElectron2nd][icut].fEta->Fill(mcpart->Eta(),wfactor);
895
896 fHist[iq][kElectron][icut].fPdgCode->Fill(mcpart->GetPdgCode());
897 fHist[iq][kElectron][icut].fPt->Fill(mcpart->Pt());
898 fHist[iq][kElectron][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart));
899 fHist[iq][kElectron][icut].fEta->Fill(mcpart->Eta());
900 }
901 else{
902 fHist[iq][kElectron][icut].fPdgCode->Fill(mcpart->GetPdgCode());
903 fHist[iq][kElectron][icut].fPt->Fill(mcpart->Pt());
904 fHist[iq][kElectron][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart));
905 fHist[iq][kElectron][icut].fEta->Fill(mcpart->Eta());
906 }
907//--------------
259c3296 908 // fill mother hadron kinematics
dbe3abbe 909 fHist[iq][keHadron][icut].fPdgCode->Fill(maPdgcodeCopy);
910 fHist[iq][keHadron][icut].fPt->Fill(partMotherCopy->Pt());
70da6c5a 911 fHist[iq][keHadron][icut].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
dbe3abbe 912 fHist[iq][keHadron][icut].fEta->Fill(partMotherCopy->Eta());
259c3296 913
ccc37cdc 914 // ratio between pT of electron and pT of mother B or direct D hadron
dbe3abbe 915 if(partMotherCopy->Pt()) fHistComm[iq][icut].fePtRatio->Fill(partMotherCopy->Pt(),mcpart->Pt()/partMotherCopy->Pt());
ccc37cdc 916 fHistComm[iq][icut].fPtCorr->Fill(partMotherCopy->Pt(),mcpart->Pt());
c2690925 917 if(TMath::Abs(partMotherCopy->GetPdgCode())==411) fHistComm[iq][icut].fPtCorrDp->Fill(partMotherCopy->Pt(),mcpart->Pt());
918 else if(TMath::Abs(partMotherCopy->GetPdgCode())==421) fHistComm[iq][icut].fPtCorrD0->Fill(partMotherCopy->Pt(),mcpart->Pt());
919 else fHistComm[iq][icut].fPtCorrDrest->Fill(partMotherCopy->Pt(),mcpart->Pt());
259c3296 920
9bcfd1ab 921 // distance between electron production point and primary vertex
922 fHistComm[iq][icut].feDistance->Fill(partMotherCopy->Pt(),decayLxy);
259c3296 923 }
924 }
925 } // end of if
926}
927
0792aa82 928//____________________________________________________________________
929void AliHFEmcQA::GetDecayedKine(AliAODMCParticle *mcpart, const Int_t kquark, Int_t kdecayed, Int_t icut)
930{
931 // decay electron kinematics
932
933 if (kquark != kCharm && kquark != kBeauty) {
934 AliDebug(1, "This task is only for heavy quark QA, return\n");
935 return;
936 }
937
938 Int_t iq = kquark - kCharm;
50685501 939 Bool_t isFinalOpenCharm = kFALSE;
0792aa82 940
faee3b18 941 if(!mcpart){
942 AliDebug(1, "no mcparticle, return\n");
943 return;
944 }
945
0792aa82 946 if ( abs(mcpart->GetPdgCode()) != kdecayed ) return;
947
948 // mother
949 Int_t iLabel = mcpart->GetMother();
950 if (iLabel<0){
951 AliDebug(1, "Stack label is negative, return\n");
952 return;
953 }
954
955 AliAODMCParticle *partMother = (AliAODMCParticle*)fMCArray->At(iLabel);
956 AliAODMCParticle *partMotherCopy = partMother;
957 Int_t maPdgcode = partMother->GetPdgCode();
958 Int_t maPdgcodeCopy = maPdgcode;
959
960 Bool_t isMotherDirectCharm = kFALSE;
961 if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
962
963 for (Int_t i=0; i<fNparents; i++){
964 if (abs(maPdgcode)==fParentSelect[0][i]){
50685501 965 isFinalOpenCharm = kTRUE;
0792aa82 966 }
967 }
50685501 968 if (!isFinalOpenCharm) return;
0792aa82 969
970 for (Int_t i=1; i<fgkMaxIter; i++){
971
972 Int_t jLabel = partMother->GetMother();
973 if (jLabel == -1){
974 isMotherDirectCharm = kTRUE;
975 break; // if there is no ancester
976 }
977 if (jLabel < 0){ // safety protection
978 AliDebug(1, "Stack label is negative, return\n");
979 return;
980 }
981
982 // if there is an ancester
983 AliAODMCParticle* grandMa = (AliAODMCParticle*)fMCArray->At(jLabel);
984 Int_t grandMaPDG = grandMa->GetPdgCode();
985
986 for (Int_t j=0; j<fNparents; j++){
987 if (abs(grandMaPDG)==fParentSelect[1][j]){
988
989 if (kquark == kCharm) return;
990 // fill electron kinematics
991 fHist[iq][kElectron2nd][icut].fPdgCode->Fill(mcpart->GetPdgCode());
992 fHist[iq][kElectron2nd][icut].fPt->Fill(mcpart->Pt());
70da6c5a 993 fHist[iq][kElectron2nd][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart));
0792aa82 994 fHist[iq][kElectron2nd][icut].fEta->Fill(mcpart->Eta());
995
996 // fill mother hadron kinematics
997 fHist[iq][kDeHadron][icut].fPdgCode->Fill(grandMaPDG);
998 fHist[iq][kDeHadron][icut].fPt->Fill(grandMa->Pt());
70da6c5a 999 fHist[iq][kDeHadron][icut].fY->Fill(AliHFEtools::GetRapidity(grandMa));
0792aa82 1000 fHist[iq][kDeHadron][icut].fEta->Fill(grandMa->Eta());
1001
1002 return;
1003 }
1004 }
1005
1006 partMother = grandMa;
1007 } // end of iteration
1008 } // end of if
1009 if ((isMotherDirectCharm == kTRUE && kquark == kCharm) || kquark == kBeauty) {
1010 for (Int_t i=0; i<fNparents; i++){
1011 if (abs(maPdgcodeCopy)==fParentSelect[iq][i]){
1012
1013 // fill electron kinematics
1014 fHist[iq][kElectron][icut].fPdgCode->Fill(mcpart->GetPdgCode());
1015 fHist[iq][kElectron][icut].fPt->Fill(mcpart->Pt());
70da6c5a 1016 fHist[iq][kElectron][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart));
0792aa82 1017 fHist[iq][kElectron][icut].fEta->Fill(mcpart->Eta());
1018
1019 // fill mother hadron kinematics
1020 fHist[iq][keHadron][icut].fPdgCode->Fill(maPdgcodeCopy);
1021 fHist[iq][keHadron][icut].fPt->Fill(partMotherCopy->Pt());
70da6c5a 1022 fHist[iq][keHadron][icut].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
0792aa82 1023 fHist[iq][keHadron][icut].fEta->Fill(partMotherCopy->Eta());
1024
1025 }
1026 }
1027 } // end of if
1028
1029}
259c3296 1030
1031//__________________________________________
75d81601 1032void AliHFEmcQA::IdentifyMother(Int_t motherlabel, Int_t &motherpdg, Int_t &grandmotherlabel)
259c3296 1033{
dbe3abbe 1034 // find mother pdg code and label
1035
75d81601 1036 if (motherlabel < 0) {
259c3296 1037 AliDebug(1, "Stack label is negative, return\n");
1038 return;
1039 }
faee3b18 1040 AliMCParticle *mctrack = NULL;
1041 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(motherlabel))))) return;
1042 TParticle *heavysMother = mctrack->Particle();
75d81601 1043 motherpdg = heavysMother->GetPdgCode();
1044 grandmotherlabel = heavysMother->GetFirstMother();
1045 AliDebug(1,Form("ancestor pdg code= %d\n",motherpdg));
259c3296 1046}
1047
1048//__________________________________________
259c3296 1049void AliHFEmcQA::HardScattering(const Int_t kquark, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
1050{
dbe3abbe 1051 // mothertype -1 means this heavy quark coming from hard vertex
1052
faee3b18 1053 AliMCParticle *mctrack1 = NULL;
1054 AliMCParticle *mctrack2 = NULL;
1055 if(!(mctrack1 = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(4))))) return;
1056 if(!(mctrack2 = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(5))))) return;
1057 TParticle *afterinitialrad1 = mctrack1->Particle();
1058 TParticle *afterinitialrad2 = mctrack2->Particle();
259c3296 1059
1060 motherlabel = -1;
1061
1062 if (abs(afterinitialrad1->GetPdgCode()) == fgkGluon && abs(afterinitialrad2->GetPdgCode()) == fgkGluon){
1063 AliDebug(1,"heavy from gluon gluon pair creation!\n");
1064 mothertype = -1;
1065 motherID = fgkGluon;
1066 }
1067 else if (abs(afterinitialrad1->GetPdgCode()) == kquark || abs(afterinitialrad2->GetPdgCode()) == kquark){ // one from Q and the other from g
1068 AliDebug(1,"heavy from flavor exitation!\n");
1069 mothertype = -1;
1070 motherID = kquark;
1071 }
1072 else if (abs(afterinitialrad1->GetPdgCode()) == abs(afterinitialrad2->GetPdgCode())){
1073 AliDebug(1,"heavy from q-qbar pair creation!\n");
1074 mothertype = -1;
1075 motherID = 1;
1076 }
1077 else {
1078 AliDebug(1,"something strange!\n");
1079 mothertype = -999;
1080 motherlabel = -999;
1081 motherID = -999;
1082 }
1083}
1084
1085//__________________________________________
259c3296 1086Bool_t AliHFEmcQA::IsFromInitialShower(Int_t inputmotherlabel, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
1087{
dbe3abbe 1088 // mothertype -2 means this heavy quark coming from initial state
1089
faee3b18 1090 AliMCParticle *mctrack = NULL;
259c3296 1091 if (inputmotherlabel==2 || inputmotherlabel==3){ // mother exist before initial state radiation
faee3b18 1092 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(inputmotherlabel))))) return kFALSE;
1093 TParticle *heavysMother = mctrack->Particle();
259c3296 1094 motherID = heavysMother->GetPdgCode();
1095 mothertype = -2; // there is mother before initial state radiation
1096 motherlabel = inputmotherlabel;
1097 AliDebug(1,"initial parton shower! \n");
1098
1099 return kTRUE;
1100 }
1101
1102 return kFALSE;
1103}
1104
1105//__________________________________________
259c3296 1106Bool_t AliHFEmcQA::IsFromFinalParton(Int_t inputmotherlabel, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
1107{
dbe3abbe 1108 // mothertype 2 means this heavy quark coming from final state
1109
faee3b18 1110 AliMCParticle *mctrack = NULL;
259c3296 1111 if (inputmotherlabel > 5){ // mother exist after hard scattering
faee3b18 1112 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(inputmotherlabel))))) return kFALSE;
1113 TParticle *heavysMother = mctrack->Particle();
259c3296 1114 motherID = heavysMother->GetPdgCode();
1115 mothertype = 2; //
1116 motherlabel = inputmotherlabel;
1117 AliDebug(1,Form("heavy quark from %d after hard scattering! \n",motherID));
1118
1119 return kTRUE;
1120 }
1121 return kFALSE;
1122}
1123
1124//__________________________________________
dbe3abbe 1125void AliHFEmcQA::ReportStrangeness(Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
259c3296 1126{
dbe3abbe 1127 // mark strange behavior
1128
259c3296 1129 mothertype = -888;
1130 motherlabel = -888;
1131 motherID = -888;
1132 AliDebug(1,"something strange!\n");
1133}
1134
1135//__________________________________________
c2690925 1136Int_t AliHFEmcQA::GetSource(const AliAODMCParticle * const mcpart)
0792aa82 1137{
9bcfd1ab 1138 // decay particle's origin
0792aa82 1139
9bcfd1ab 1140 //if ( abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return -1;
0792aa82 1141
1142 Int_t origin = -1;
50685501 1143 Bool_t isFinalOpenCharm = kFALSE;
0792aa82 1144
faee3b18 1145 if(!mcpart){
1146 AliDebug(1, "Stack label is negative or no mcparticle, return\n");
1147 return -1;
1148 }
1149
0792aa82 1150 // mother
1151 Int_t iLabel = mcpart->GetMother();
1152 if (iLabel<0){
1153 AliDebug(1, "Stack label is negative, return\n");
1154 return -1;
1155 }
1156
1157 AliAODMCParticle *partMother = (AliAODMCParticle*)fMCArray->At(iLabel);
1158 Int_t maPdgcode = partMother->GetPdgCode();
1159
1160 // if the mother is charmed hadron
1161 if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
1162
1163 for (Int_t i=0; i<fNparents; i++){
1164 if (abs(maPdgcode)==fParentSelect[0][i]){
50685501 1165 isFinalOpenCharm = kTRUE;
0792aa82 1166 }
1167 }
50685501 1168 if (!isFinalOpenCharm) return -1;
0792aa82 1169
1170 // iterate until you find B hadron as a mother or become top ancester
1171 for (Int_t i=1; i<fgkMaxIter; i++){
1172
1173 Int_t jLabel = partMother->GetMother();
1174 if (jLabel == -1){
1175 origin = kDirectCharm;
1176 return origin;
1177 }
1178 if (jLabel < 0){ // safety protection
1179 AliDebug(1, "Stack label is negative, return\n");
1180 return -1;
1181 }
1182
1183 // if there is an ancester
1184 AliAODMCParticle* grandMa = (AliAODMCParticle*)fMCArray->At(jLabel);
1185 Int_t grandMaPDG = grandMa->GetPdgCode();
1186
70b5ea26 1187 for (Int_t j=0; j<fNparents; j++){
1188 if (abs(grandMaPDG)==fParentSelect[1][j]){
0792aa82 1189 origin = kBeautyCharm;
1190 return origin;
1191 }
1192 }
1193
1194 partMother = grandMa;
1195 } // end of iteration
1196 } // end of if
1197 else if ( int(abs(maPdgcode)/100.) == kBeauty || int(abs(maPdgcode)/1000.) == kBeauty ) {
1198 for (Int_t i=0; i<fNparents; i++){
1199 if (abs(maPdgcode)==fParentSelect[1][i]){
1200 origin = kDirectBeauty;
1201 return origin;
1202 }
1203 }
1204 } // end of if
1205 else if ( abs(maPdgcode) == 22 ) {
1206 origin = kGamma;
1207 return origin;
1208 } // end of if
1209 else if ( abs(maPdgcode) == 111 ) {
1210 origin = kPi0;
1211 return origin;
1212 } // end of if
1213
1214 return origin;
1215}
1216
1217//__________________________________________
c2690925 1218Int_t AliHFEmcQA::GetSource(const TParticle * const mcpart)
0792aa82 1219{
9bcfd1ab 1220 // decay particle's origin
0792aa82 1221
9bcfd1ab 1222 //if ( abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return -1;
0792aa82 1223
1224 Int_t origin = -1;
50685501 1225 Bool_t isFinalOpenCharm = kFALSE;
0792aa82 1226
faee3b18 1227 if(!mcpart){
1228 AliDebug(1, "no mcparticle, return\n");
1229 return -1;
1230 }
1231
0792aa82 1232 Int_t iLabel = mcpart->GetFirstMother();
1233 if (iLabel<0){
1234 AliDebug(1, "Stack label is negative, return\n");
1235 return -1;
1236 }
1237
faee3b18 1238 AliMCParticle *mctrack = NULL;
1239 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return -1;
1240 TParticle *partMother = mctrack->Particle();
0792aa82 1241 Int_t maPdgcode = partMother->GetPdgCode();
1242
1243 // if the mother is charmed hadron
1244 if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
1245
1246 for (Int_t i=0; i<fNparents; i++){
1247 if (abs(maPdgcode)==fParentSelect[0][i]){
50685501 1248 isFinalOpenCharm = kTRUE;
0792aa82 1249 }
1250 }
50685501 1251 if (!isFinalOpenCharm) return -1;
0792aa82 1252
1253 // iterate until you find B hadron as a mother or become top ancester
1254 for (Int_t i=1; i<fgkMaxIter; i++){
1255
1256 Int_t jLabel = partMother->GetFirstMother();
1257 if (jLabel == -1){
1258 origin = kDirectCharm;
1259 return origin;
1260 }
1261 if (jLabel < 0){ // safety protection
1262 AliDebug(1, "Stack label is negative, return\n");
1263 return -1;
1264 }
1265
1266 // if there is an ancester
faee3b18 1267 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return -1;
1268 TParticle* grandMa = mctrack->Particle();
0792aa82 1269 Int_t grandMaPDG = grandMa->GetPdgCode();
1270
70b5ea26 1271 for (Int_t j=0; j<fNparents; j++){
1272 if (abs(grandMaPDG)==fParentSelect[1][j]){
0792aa82 1273 origin = kBeautyCharm;
1274 return origin;
1275 }
1276 }
1277
1278 partMother = grandMa;
1279 } // end of iteration
1280 } // end of if
1281 else if ( int(abs(maPdgcode)/100.) == kBeauty || int(abs(maPdgcode)/1000.) == kBeauty ) {
1282 for (Int_t i=0; i<fNparents; i++){
1283 if (abs(maPdgcode)==fParentSelect[1][i]){
1284 origin = kDirectBeauty;
1285 return origin;
1286 }
1287 }
1288 } // end of if
1289 else if ( abs(maPdgcode) == 22 ) {
1290 origin = kGamma;
1291 return origin;
1292 } // end of if
1293 else if ( abs(maPdgcode) == 111 ) {
1294 origin = kPi0;
1295 return origin;
1296 } // end of if
1297 else origin = kElse;
1298
1299 return origin;
1300}
70da6c5a 1301
1302//__________________________________________
e3fc062d 1303Int_t AliHFEmcQA::GetElecSource(TParticle * const mcpart)
1304{
1305 // decay particle's origin
70da6c5a 1306
e3fc062d 1307 if(!mcpart){
1308 AliDebug(1, "no mcparticle, return\n");
1309 return -1;
1310 }
1311
bf892a6a 1312 if ( abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return kMisID;
1313
1314 Int_t origin = -1;
1315 Bool_t isFinalOpenCharm = kFALSE;
1316
e3fc062d 1317 Int_t iLabel = mcpart->GetFirstMother();
1318 if (iLabel<0){
1319 AliDebug(1, "Stack label is negative, return\n");
1320 return -1;
1321 }
1322
1323 AliMCParticle *mctrack = NULL;
c2690925 1324 Int_t tmpMomLabel=0;
e3fc062d 1325 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return -1;
1326 TParticle *partMother = mctrack->Particle();
c2690925 1327 TParticle *partMotherCopy = mctrack->Particle();
e3fc062d 1328 Int_t maPdgcode = partMother->GetPdgCode();
1329
1330 // if the mother is charmed hadron
1331 if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
1332
1333 for (Int_t i=0; i<fNparents; i++){
1334 if (abs(maPdgcode)==fParentSelect[0][i]){
1335 isFinalOpenCharm = kTRUE;
1336 }
1337 }
1338 if (!isFinalOpenCharm) return -1;
1339
1340 // iterate until you find B hadron as a mother or become top ancester
1341 for (Int_t i=1; i<fgkMaxIter; i++){
1342
1343 Int_t jLabel = partMother->GetFirstMother();
1344 if (jLabel == -1){
1345 origin = kDirectCharm;
1346 return origin;
1347 }
1348 if (jLabel < 0){ // safety protection
1349 AliDebug(1, "Stack label is negative, return\n");
1350 return -1;
1351 }
1352
1353 // if there is an ancester
1354 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return -1;
1355 TParticle* grandMa = mctrack->Particle();
1356 Int_t grandMaPDG = grandMa->GetPdgCode();
1357
1358 for (Int_t j=0; j<fNparents; j++){
1359 if (abs(grandMaPDG)==fParentSelect[1][j]){
1360 origin = kBeautyCharm;
1361 return origin;
1362 }
1363 }
1364
1365 partMother = grandMa;
1366 } // end of iteration
1367 } // end of if
1368 else if ( int(abs(maPdgcode)/100.) == kBeauty || int(abs(maPdgcode)/1000.) == kBeauty ) {
1369 for (Int_t i=0; i<fNparents; i++){
1370 if (abs(maPdgcode)==fParentSelect[1][i]){
1371 origin = kDirectBeauty;
1372 return origin;
1373 }
1374 }
1375 } // end of if
c2690925 1376 else if ( abs(maPdgcode) == 22 ) { //conversion
1377
1378 tmpMomLabel = partMotherCopy->GetFirstMother();
1379 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tmpMomLabel))))) return -1;
1380 partMother = mctrack->Particle();
1381 maPdgcode = partMother->GetPdgCode();
1382 if ( abs(maPdgcode) == 111 ) {
1383 origin = kGammaPi0;
1384 return origin;
1385 }
1386 else if ( abs(maPdgcode) == 221 ) {
1387 origin = kGammaEta;
1388 return origin;
1389 }
1390 else if ( abs(maPdgcode) == 223 ) {
1391 origin = kGammaOmega;
1392 return origin;
1393 }
1394 else if ( abs(maPdgcode) == 333 ) {
1395 origin = kGammaPhi;
1396 return origin;
1397 }
1398 else if ( abs(maPdgcode) == 331 ) {
1399 origin = kGammaEtaPrime;
1400 return origin;
1401 }
1402 else if ( abs(maPdgcode) == 113 ) {
1403 origin = kGammaRho0;
1404 return origin;
1405 }
1406 else origin = kElse;
1407 //origin = kGamma; // finer category above
e3fc062d 1408 return origin;
c2690925 1409
e3fc062d 1410 } // end of if
1411 else if ( abs(maPdgcode) == 111 ) {
1412 origin = kPi0;
1413 return origin;
1414 } // end of if
1415 else if ( abs(maPdgcode) == 221 ) {
1416 origin = kEta;
1417 return origin;
1418 } // end of if
1419 else if ( abs(maPdgcode) == 223 ) {
1420 origin = kOmega;
1421 return origin;
1422 } // end of if
1423 else if ( abs(maPdgcode) == 333 ) {
1424 origin = kPhi;
1425 return origin;
1426 } // end of if
1427 else if ( abs(maPdgcode) == 331 ) {
1428 origin = kEtaPrime;
1429 return origin;
1430 } // end of if
1431 else if ( abs(maPdgcode) == 113 ) {
1432 origin = kRho0;
1433 return origin;
1434 } // end of if
1435 else origin = kElse;
1436
1437 return origin;
70da6c5a 1438}
1439
1440//__________________________________________
9250ffbf 1441Double_t AliHFEmcQA::GetWeightFactor(AliMCParticle *mctrack){
1442 //
1443 // Get weighting factor for the realistic background estimation
1444 //
1445 AliMCParticle *mctrackmother = NULL;
1446 Double_t weightElecBg = 0.;
1447 Double_t mesonPt = 0.;
1448 Double_t bgcategory = 0.;
1449 Int_t mArr = -1;
1450 Int_t mesonID = GetElecSource(mctrack->Particle());
1451 if(mesonID==kGammaPi0 || mesonID==kPi0) mArr=0; //pion
1452 else if(mesonID==kGammaEta || mesonID==kEta) mArr=1; //eta
1453 else if(mesonID==kGammaOmega || mesonID==kOmega) mArr=2; //omega
1454 else if(mesonID==kGammaPhi || mesonID==kPhi) mArr=3; //phi
1455 else if(mesonID==kGammaEtaPrime || mesonID==kEtaPrime) mArr=4; //etaprime
1456 else if(mesonID==kGammaRho0 || mesonID==kRho0) mArr=5; //rho
1457
1458 if(!(mArr<0)){
1459 if(mesonID>=kGammaPi0) { // conversion electron, be careful with the enum odering
1460 Int_t glabel=TMath::Abs(mctrack->GetMother()); // gamma label
1461 if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
1462 glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's label
1463 if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
1464 mesonPt = mctrackmother->Pt(); //meson pt
1465 bgcategory = 1.;
1466 }
1467 }
1468 }
1469 else{ // nonHFE except for the conversion electron
1470 Int_t glabel=TMath::Abs(mctrack->GetMother());
1471 if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
1472 mesonPt=mctrackmother->Pt(); //meson pt
1473 bgcategory = -1.;
1474 }
1475 }
1476 weightElecBg=fElecBackgroundFactor[mArr][kBgPtBins-1];
1477 for(int ii=0; ii<kBgPtBins; ii++){
1478 if( (mesonPt > fBinLimit[ii]) && (mesonPt < fBinLimit[ii+1])){
1479 weightElecBg = fElecBackgroundFactor[mArr][ii];
1480 }
1481 }
1482 /*for(int jj=0; jj<6; jj++){
1483 for(int ii=0; ii<kBgPtBins; ii++){
1484 printf("species= %d ptbin= %d wfactor= %lf\n",jj,ii,fElecBackgroundFactor[jj][ii]);
1485 }
1486 }*/
1487 }
1488
1489 return bgcategory*weightElecBg;
1490}
1491
1492//__________________________________________
70da6c5a 1493void AliHFEmcQA::AliHists::FillList(TList *l) const {
1494 //
1495 // Fill Histos into a list for output
1496 //
1497 if(fPdgCode) l->Add(fPdgCode);
1498 if(fPt) l->Add(fPt);
1499 if(fY) l->Add(fY);
1500 if(fEta) l->Add(fEta);
1501}
1502
1503//__________________________________________
1504void AliHFEmcQA::AliHistsComm::FillList(TList *l) const {
1505 //
1506 // Fill Histos into a list for output
1507 //
1508 if(fNq) l->Add(fNq);
1509 if(fProcessID) l->Add(fProcessID);
1510 if(fePtRatio) l->Add(fePtRatio);
ccc37cdc 1511 if(fPtCorr) l->Add(fPtCorr);
c2690925 1512 if(fPtCorrDp) l->Add(fPtCorrDp);
1513 if(fPtCorrD0) l->Add(fPtCorrD0);
1514 if(fPtCorrDrest) l->Add(fPtCorrDrest);
70da6c5a 1515 if(fDePtRatio) l->Add(fDePtRatio);
1516 if(feDistance) l->Add(feDistance);
1517 if(fDeDistance) l->Add(fDeDistance);
1518}
1519