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