1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 // QA class of Heavy Flavor quark and fragmeted/decayed particles
17 // -Check kinematics of Heavy Quarks/hadrons, and decayed leptons
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
27 // MinJung Kweon <minjung@physi.uni-heidelberg.de>
34 #include <TParticle.h>
35 #include "TTreeStream.h"
38 #include <AliMCEvent.h>
39 #include <AliGenEventHeader.h>
40 #include <AliAODMCParticle.h>
43 #include "AliHFEmcQA.h"
44 #include "AliHFEtools.h"
45 #include "AliHFEcollection.h"
49 //_______________________________________________________________________________________________
50 AliHFEmcQA::AliHFEmcQA() :
55 ,fMCQACollection(NULL)
60 ,fIsppMultiBin(kFALSE)
62 ,fIsDebugStreamerON(kFALSE)
69 ,fHfeImpactnsigmaR(-999)
71 ,fGetWeightHist(kFALSE)
73 // Default constructor
74 for(Int_t mom = 0; mom < 9; mom++){
77 for(Int_t mom = 0; mom < 50; mom++){
78 fHeavyQuark[mom] = NULL;
80 for(Int_t mom = 0; mom < 2; mom++){
83 memset(fElecBackgroundFactor, 0, sizeof(Double_t) * kElecBgSpecies * kBgPtBins * kCentBins * kBgLevels);
84 memset(fBinLimit, 0, sizeof(Double_t) * (kBgPtBins+1));
87 //_______________________________________________________________________________________________
88 AliHFEmcQA::AliHFEmcQA(const AliHFEmcQA&p):
93 ,fQAhistos(p.fQAhistos)
94 ,fMCQACollection(p.fMCQACollection)
95 ,fNparents(p.fNparents)
99 ,fIsppMultiBin(kFALSE)
101 ,fIsDebugStreamerON(kFALSE)
108 ,fHfeImpactnsigmaR(0)
110 ,fGetWeightHist(kFALSE)
113 for(Int_t mom = 0; mom < 9; mom++){
116 for(Int_t mom = 0; mom < 50; mom++){
117 fHeavyQuark[mom] = NULL;
119 for(Int_t mom = 0; mom < 2; mom++){
122 memset(fElecBackgroundFactor, 0, sizeof(Double_t) * kElecBgSpecies * kBgPtBins * kCentBins * kBgLevels);
123 memset(fBinLimit, 0, sizeof(Double_t) * (kBgPtBins+1));
125 //_______________________________________________________________________________________________
127 AliHFEmcQA::operator=(const AliHFEmcQA &)
129 // Assignment operator
131 AliInfo("Not yet implemented.");
135 //_______________________________________________________________________________________________
136 AliHFEmcQA::~AliHFEmcQA()
140 if(fTreeStream && fIsDebugStreamerON) delete fTreeStream;
141 AliInfo("Analysis Done.");
143 //_______________________________________________________________________________________________
144 void AliHFEmcQA::PostAnalyze() const
150 //_______________________________________________________________________________________________
151 void AliHFEmcQA::SetBackgroundWeightFactor(Double_t *elecBackgroundFactor, Double_t *binLimit)
154 // copy background weighting factors into data member
157 memcpy(fElecBackgroundFactor,elecBackgroundFactor,sizeof(Double_t) * kElecBgSpecies * kBgPtBins * kCentBins * kBgLevels);
158 memcpy(fBinLimit,binLimit,sizeof(Double_t) * (kBgPtBins+1));
160 //__________________________________________
161 void AliHFEmcQA::CreatDefaultHistograms(TList * const qaList)
164 // make default histograms
170 fQAhistos->SetName("MCqa");
172 CreateHistograms(AliHFEmcQA::kCharm); // create histograms for charm
173 CreateHistograms(AliHFEmcQA::kBeauty); // create histograms for beauty
174 CreateHistograms(AliHFEmcQA::kOthers); // create histograms for beauty
176 // prepare 2D(pt vs Y) histogram for D spectra, we consider following 9 particles
177 const Int_t nbspecies = 9;
178 TString kDspecies[nbspecies];
179 kDspecies[0]="411"; //D+
180 kDspecies[1]="421"; //D0
181 kDspecies[2]="431"; //Ds+
182 kDspecies[3]="4122"; //Lambdac+
183 kDspecies[4]="4132"; //Ksic0
184 kDspecies[5]="4232"; //Ksic+
185 kDspecies[6]="4332"; //OmegaC0
186 kDspecies[7]="413"; //D*(2010)+
187 kDspecies[8]="423"; //D*(2007)0
189 const Double_t kPtbound[2] = {0.1, 20.}; //bin taken for considering inclusive e analysis binning
191 iBin[0] = 44; // bins in pt for log binning
192 iBin[1] = 23; // bins in pt for pi0 measurement binning
193 Double_t* binEdges[1];
194 binEdges[0] = AliHFEtools::MakeLogarithmicBinning(iBin[0], kPtbound[0], kPtbound[1]);
196 // bin size is chosen to consider ALICE D measurement
197 const Int_t nptbins = 15;
198 const Int_t nybins = 9;
199 Double_t xbins[nptbins+1]={0,0.5,1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,12,16,24,32,40,50}; //pt binning for the final 7 TeV D measurement
200 Double_t ybins[nybins+1]={-7.5,-1.0,-0.9,-0.8,-0.5,0.5,0.8,0.9,1.0,7.5}; // y binning
202 for (Int_t iDmeson=0; iDmeson<nbspecies; iDmeson++){
203 hname = "Dmeson"+kDspecies[iDmeson];
204 fhD[iDmeson] = new TH2F(hname,hname+";p_{T} (GeV/c)",nptbins,xbins,nybins,ybins);
205 if(fQAhistos) fQAhistos->Add(fhD[iDmeson]);
208 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
211 if(fIsPbPb) kNcent=11;
214 if(fIsppMultiBin) kNcent=8;
218 fMCQACollection = new AliHFEcollection("TaskMCQA", "MC QA histos for meason pt spectra");
220 for(Int_t centbin=0; centbin<kNcent; centbin++)
222 fMCQACollection->CreateTH1Farray(Form("pionspectra_centrbin%i",centbin), "pion yields: MC p_{t} ", iBin[1],kPtRange);
223 fMCQACollection->CreateTH1Farray(Form("etaspectra_centrbin%i",centbin), "eta yields: MC p_{t} ", iBin[1],kPtRange);
224 fMCQACollection->CreateTH1Farray(Form("omegaspectra_centrbin%i",centbin), "omega yields: MC p_{t} ", iBin[1],kPtRange);
225 fMCQACollection->CreateTH1Farray(Form("phispectra_centrbin%i",centbin), "phi yields: MC p_{t} ", iBin[1],kPtRange);
226 fMCQACollection->CreateTH1Farray(Form("etapspectra_centrbin%i",centbin), "etap yields: MC p_{t} ", iBin[1],kPtRange);
227 fMCQACollection->CreateTH1Farray(Form("rhospectra_centrbin%i",centbin), "rho yields: MC p_{t} ", iBin[1],kPtRange);
229 fMCQACollection->CreateTH1Farray(Form("pionspectrapri_centrbin%i",centbin), "pion yields: MC p_{t} ", iBin[1],kPtRange);
230 fMCQACollection->CreateTH1Farray(Form("etaspectrapri_centrbin%i",centbin), "eta yields: MC p_{t} ", iBin[1],kPtRange);
231 fMCQACollection->CreateTH1Farray(Form("pionspectrasec_centrbin%i",centbin), "pion yields: MC p_{t} ", iBin[1],kPtRange);
232 fMCQACollection->CreateTH1Farray(Form("etaspectrasec_centrbin%i",centbin), "eta yields: MC p_{t} ", iBin[1],kPtRange);
233 fMCQACollection->CreateTH1Farray(Form("pionspectraLogpri_centrbin%i",centbin), "pion yields: MC p_{t} ", iBin[1],kPtRange);
234 fMCQACollection->CreateTH1Farray(Form("etaspectraLogpri_centrbin%i",centbin), "eta yields: MC p_{t} ", iBin[1],kPtRange);
235 fMCQACollection->CreateTH1Farray(Form("pionspectraLogsec_centrbin%i",centbin), "pion yields: MC p_{t} ", iBin[1],kPtRange);
236 fMCQACollection->CreateTH1Farray(Form("etaspectraLogsec_centrbin%i",centbin), "eta yields: MC p_{t} ", iBin[1],kPtRange);
238 fMCQACollection->CreateTH1F(Form("pionspectraLog_centrbin%i",centbin), "pion yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
239 fMCQACollection->CreateTH1F(Form("etaspectraLog_centrbin%i",centbin), "eta yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
240 fMCQACollection->CreateTH1F(Form("omegaspectraLog_centrbin%i",centbin), "omega yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
241 fMCQACollection->CreateTH1F(Form("phispectraLog_centrbin%i",centbin), "phi yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
242 fMCQACollection->CreateTH1F(Form("etapspectraLog_centrbin%i",centbin), "etap yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
243 fMCQACollection->CreateTH1F(Form("rhospectraLog_centrbin%i",centbin), "rho yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
244 fMCQACollection->CreateTH1F(Form("kaonspectraLog_centrbin%i",centbin), "kaon yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
245 fMCQACollection->CreateTH1F(Form("k0LspectraLog_centrbin%i",centbin), "k0L yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
247 fMCQACollection->CreateTH2F(Form("pionspectraLog2D_centrbin%i",centbin), "pion yields: MC p_{t} ", 32, -1.5, 30.5, iBin[0],kPtbound[0], kPtbound[1], 1);
248 fMCQACollection->CreateTH2F(Form("etaspectraLog2D_centrbin%i",centbin), "eta yields: MC p_{t} ", 32, -1.5, 30.5, iBin[0],kPtbound[0], kPtbound[1], 1);
249 fMCQACollection->CreateTH2F(Form("omegaspectraLog2D_centrbin%i",centbin), "omega yields: MC p_{t} ", 32, -1.5, 30.5, iBin[0],kPtbound[0], kPtbound[1], 1);
250 fMCQACollection->CreateTH2F(Form("phispectraLog2D_centrbin%i",centbin), "phi yields: MC p_{t} ", 32, -1.5, 30.5, iBin[0],kPtbound[0], kPtbound[1], 1);
251 fMCQACollection->CreateTH2F(Form("etapspectraLog2D_centrbin%i",centbin), "etap yields: MC p_{t} ", 32, -1.5, 30.5, iBin[0],kPtbound[0], kPtbound[1], 1);
252 fMCQACollection->CreateTH2F(Form("rhospectraLog2D_centrbin%i",centbin), "rho yields: MC p_{t} ", 32, -1.5, 30.5, iBin[0],kPtbound[0], kPtbound[1], 1);
254 fMCQACollection->CreateTH1F(Form("piondaughters_centrbin%i",centbin), "pion yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
255 fMCQACollection->CreateTH1F(Form("etadaughters_centrbin%i",centbin), "eta yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
256 fMCQACollection->CreateTH1F(Form("omegadaughters_centrbin%i",centbin), "omega yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
257 fMCQACollection->CreateTH1F(Form("phidaughters_centrbin%i",centbin), "phi yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
258 fMCQACollection->CreateTH1F(Form("etapdaughters_centrbin%i",centbin), "etap yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
259 fMCQACollection->CreateTH1F(Form("rhodaughters_centrbin%i",centbin), "rho yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
261 fMCQACollection->CreateTH1F(Form("pionspectraPrimary_centrbin%i",centbin), "pion yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
262 fMCQACollection->CreateTH1F(Form("etaspectraPrimary_centrbin%i",centbin), "eta yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
265 fQAhistos->Add(fMCQACollection->GetList());
267 if(!fTreeStream && fIsDebugStreamerON){
268 fTreeStream = new TTreeSRedirector(Form("HFEmcqadebugTree%s.root", GetName()));
273 //__________________________________________
274 void AliHFEmcQA::CreateHistograms(const Int_t kquark)
278 if (!(kquark == kCharm || kquark == kBeauty || kquark == kOthers)) {
279 AliDebug(1, "This task is only for heavy quark QA, return\n");
282 Int_t iq = kquark - kCharm;
284 TString kqTypeLabel[fgkqType];
285 if (kquark == kCharm){
286 kqTypeLabel[kQuark]="c";
287 kqTypeLabel[kantiQuark]="cbar";
288 kqTypeLabel[kHadron]="cHadron";
289 kqTypeLabel[keHadron]="ceHadron";
290 kqTypeLabel[kDeHadron]="nullHadron";
291 kqTypeLabel[kElectron]="ce";
292 kqTypeLabel[kElectron2nd]="nulle";
293 } else if (kquark == kBeauty){
294 kqTypeLabel[kQuark]="b";
295 kqTypeLabel[kantiQuark]="bbar";
296 kqTypeLabel[kHadron]="bHadron";
297 kqTypeLabel[keHadron]="beHadron";
298 kqTypeLabel[kDeHadron]="bDeHadron";
299 kqTypeLabel[kElectron]="be";
300 kqTypeLabel[kElectron2nd]="bce";
301 } else if (kquark == kOthers){
302 kqTypeLabel[kGamma-4]="gammae";
303 kqTypeLabel[kPi0-4]="pi0e";
304 kqTypeLabel[kElse-4]="elsee";
305 kqTypeLabel[kMisID-4]="miside";
308 TString kqEtaRangeLabel[fgkEtaRanges];
309 kqEtaRangeLabel[0] = "mcqa_";
310 kqEtaRangeLabel[1] = "mcqa_barrel_";
311 kqEtaRangeLabel[2] = "mcqa_unitY_";
313 const Double_t kPtbound[2] = {0.1, 20.}; //bin taken for considering inclusive e analysis binning
314 const Int_t nptbinning1 = 35;
316 iBin[0] = 44; // bins in pt
317 iBin[1] = nptbinning1; // bins in pt
318 Double_t* binEdges[1];
319 binEdges[0] = AliHFEtools::MakeLogarithmicBinning(iBin[0], kPtbound[0], kPtbound[1]);
321 // new binning for final electron analysis
322 const Double_t kPtbinning1[nptbinning1+1] = {0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4, 1.5, 1.75, 2., 2.25, 2.5, 2.75, 3., 3.5, 4., 4.5, 5., 5.5, 6., 7., 8., 10., 12., 14., 16., 18., 20.};
324 const Int_t ndptbins = 500;
325 Double_t xcorrbin[ndptbins+1];
326 for (int icorrbin = 0; icorrbin< ndptbins+1; icorrbin++){
327 xcorrbin[icorrbin]=icorrbin*0.1;
331 if(!fIsPbPb&&!fIsppMultiBin) fCentrmax=0;
332 else fCentrmax = kCentBins;
334 if(kquark == kOthers){
335 for (Int_t icut = 0; icut < fgkEtaRanges; icut++ ){
336 for (Int_t iqType = 0; iqType < 4; iqType++ ){
337 for(Int_t icentr = 0; icentr < (fCentrmax+1); icentr++)
339 hname = kqEtaRangeLabel[icut]+"Pt_"+kqTypeLabel[iqType];
340 fHist[iq][iqType][icut][icentr].fPt = new TH1F(Form("%sCentr_%i",hname.Data(),icentr),Form("%sCentr_%i;p_{T} (GeV/c)",hname.Data(),icentr),60,0.25,30.25);
341 hname = kqEtaRangeLabel[icut]+"Y_"+kqTypeLabel[iqType];
342 fHist[iq][iqType][icut][icentr].fY = new TH1F(Form("%sCentr_%i",hname.Data(),icentr),Form("%sCentr_%i;y",hname.Data(),icentr),150,-7.5,7.5);
343 hname = kqEtaRangeLabel[icut]+"Eta_"+kqTypeLabel[iqType];
344 fHist[iq][iqType][icut][icentr].fEta = new TH1F(Form("%sCentr_%i",hname.Data(),icentr),Form("%sCentr_%i;eta",hname.Data(),icentr),150,-7.5,7.5);
346 if(fQAhistos) fHist[iq][iqType][icut][icentr].FillList(fQAhistos);
352 for (Int_t icut = 0; icut < fgkEtaRanges; icut++ ){
353 for (Int_t iqType = 0; iqType < fgkqType; iqType++ ){
354 if (iqType < keHadron && icut > 0) continue; // don't duplicate histogram for quark and hadron
355 for(Int_t icentr = 0; icentr<(fCentrmax+1); icentr++)
357 hname = kqEtaRangeLabel[icut]+"PdgCode_"+kqTypeLabel[iqType];
358 fHist[iq][iqType][icut][icentr].fPdgCode = new TH1F(Form("%sCentr_%i",hname.Data(),icentr),Form("%sCentr_%i;PdgCode",hname.Data(),icentr),20001,-10000.5,10000.5);
359 hname = kqEtaRangeLabel[icut]+"Pt_"+kqTypeLabel[iqType];
360 fHist[iq][iqType][icut][icentr].fPt = new TH1F(Form("%sCentr_%i",hname.Data(),icentr),Form("%sCentr_%i;p_{T} (GeV/c)",hname.Data(),icentr),iBin[1],kPtbinning1); // new binning
361 hname = kqEtaRangeLabel[icut]+"Y_"+kqTypeLabel[iqType];
362 fHist[iq][iqType][icut][icentr].fY = new TH1F(Form("%sCentr_%i",hname.Data(),icentr),Form("%sCentr_%i;y",hname.Data(),icentr),150,-7.5,7.5);
363 hname = kqEtaRangeLabel[icut]+"Eta_"+kqTypeLabel[iqType];
364 fHist[iq][iqType][icut][icentr].fEta = new TH1F(Form("%sCentr_%i",hname.Data(),icentr),Form("%sCentr_%i;eta",hname.Data(),icentr),150,-7.5,7.5);
366 if(fQAhistos) fHist[iq][iqType][icut][icentr].FillList(fQAhistos);
371 for (Int_t icut = 0; icut < fgkEtaRanges; icut++ ){
372 hname = kqEtaRangeLabel[icut]+"PtCorr_"+kqTypeLabel[kQuark];
373 fHistComm[iq][icut].fPtCorr = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1);
374 hname = kqEtaRangeLabel[icut]+"PtCorrDp_"+kqTypeLabel[kQuark];
375 fHistComm[iq][icut].fPtCorrDp= new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1);
376 hname = kqEtaRangeLabel[icut]+"PtCorrD0_"+kqTypeLabel[kQuark];
377 fHistComm[iq][icut].fPtCorrD0 = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1);
378 hname = kqEtaRangeLabel[icut]+"PtCorrDrest_"+kqTypeLabel[kQuark];
379 fHistComm[iq][icut].fPtCorrDrest = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1);
381 hname = kqEtaRangeLabel[icut]+"ePtRatio_"+kqTypeLabel[kQuark];
382 fHistComm[iq][icut].fePtRatio = new TH2F(hname,hname+";p_{T} (GeV/c);momentum fraction",200,0,20,100,0,1);
383 hname = kqEtaRangeLabel[icut]+"DePtRatio_"+kqTypeLabel[kQuark];
384 fHistComm[iq][icut].fDePtRatio = new TH2F(hname,hname+";p_{T} (GeV/c);momentum fraction",100,0,20,100,0,1);
385 hname = kqEtaRangeLabel[icut]+"eDistance_"+kqTypeLabel[kQuark];
386 fHistComm[iq][icut].feDistance= new TH2F(hname,hname+";p_{T} (GeV/c);distance (cm)",100,0,20,200,0,2);
387 hname = kqEtaRangeLabel[icut]+"DeDistance_"+kqTypeLabel[kQuark];
388 fHistComm[iq][icut].fDeDistance= new TH2F(hname,hname+";p_{T} (GeV/c);distance (cm)",100,0,20,200,0,2);
391 hname = kqEtaRangeLabel[icut]+"PtCorrDinein_"+kqTypeLabel[kQuark];
392 fHistComm[iq][icut].fPtCorrDinein = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
393 hname = kqEtaRangeLabel[icut]+"PtCorrDineout_"+kqTypeLabel[kQuark];
394 fHistComm[iq][icut].fPtCorrDineout = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
395 hname = kqEtaRangeLabel[icut]+"PtCorrDoutein_"+kqTypeLabel[kQuark];
396 fHistComm[iq][icut].fPtCorrDoutein = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
397 hname = kqEtaRangeLabel[icut]+"PtCorrDouteout_"+kqTypeLabel[kQuark];
398 fHistComm[iq][icut].fPtCorrDouteout = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
400 hname = kqEtaRangeLabel[icut]+"PtCorrDpDinein_"+kqTypeLabel[kQuark];
401 fHistComm[iq][icut].fPtCorrDpDinein = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
402 hname = kqEtaRangeLabel[icut]+"PtCorrDpDineout_"+kqTypeLabel[kQuark];
403 fHistComm[iq][icut].fPtCorrDpDineout = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
404 hname = kqEtaRangeLabel[icut]+"PtCorrDpDoutein_"+kqTypeLabel[kQuark];
405 fHistComm[iq][icut].fPtCorrDpDoutein = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
406 hname = kqEtaRangeLabel[icut]+"PtCorrDpDouteout_"+kqTypeLabel[kQuark];
407 fHistComm[iq][icut].fPtCorrDpDouteout = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
409 hname = kqEtaRangeLabel[icut]+"PtCorrD0Dinein_"+kqTypeLabel[kQuark];
410 fHistComm[iq][icut].fPtCorrD0Dinein = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
411 hname = kqEtaRangeLabel[icut]+"PtCorrD0Dineout_"+kqTypeLabel[kQuark];
412 fHistComm[iq][icut].fPtCorrD0Dineout = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
413 hname = kqEtaRangeLabel[icut]+"PtCorrD0Doutein_"+kqTypeLabel[kQuark];
414 fHistComm[iq][icut].fPtCorrD0Doutein = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
415 hname = kqEtaRangeLabel[icut]+"PtCorrD0Douteout_"+kqTypeLabel[kQuark];
416 fHistComm[iq][icut].fPtCorrD0Douteout = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
418 hname = kqEtaRangeLabel[icut]+"PtCorrDrestDinein_"+kqTypeLabel[kQuark];
419 fHistComm[iq][icut].fPtCorrDrestDinein = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
420 hname = kqEtaRangeLabel[icut]+"PtCorrDrestDineout_"+kqTypeLabel[kQuark];
421 fHistComm[iq][icut].fPtCorrDrestDineout = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
422 hname = kqEtaRangeLabel[icut]+"PtCorrDrestDoutein_"+kqTypeLabel[kQuark];
423 fHistComm[iq][icut].fPtCorrDrestDoutein = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
424 hname = kqEtaRangeLabel[icut]+"PtCorrDrestDouteout_"+kqTypeLabel[kQuark];
425 fHistComm[iq][icut].fPtCorrDrestDouteout = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
427 hname = kqEtaRangeLabel[icut]+"fEtaCorrD_"+kqTypeLabel[kQuark];
428 fHistComm[iq][icut].fEtaCorrD = new TH2F(hname,hname+";D Y;e eta",200,-10,10,200,-10,10);
429 hname = kqEtaRangeLabel[icut]+"fEtaCorrDp_"+kqTypeLabel[kQuark];
430 fHistComm[iq][icut].fEtaCorrDp = new TH2F(hname,hname+";D Y;e eta",200,-10,10,200,-10,10);
431 hname = kqEtaRangeLabel[icut]+"fEtaCorrD0_"+kqTypeLabel[kQuark];
432 fHistComm[iq][icut].fEtaCorrD0 = new TH2F(hname,hname+";D Y;e eta",200,-10,10,200,-10,10);
433 hname = kqEtaRangeLabel[icut]+"fEtaCorrDrest_"+kqTypeLabel[kQuark];
434 fHistComm[iq][icut].fEtaCorrDrest = new TH2F(hname,hname+";D Y;e eta",200,-10,10,200,-10,10);
436 hname = kqEtaRangeLabel[icut]+"fEtaCorrGD_"+kqTypeLabel[kQuark];
437 fHistComm[iq][icut].fEtaCorrGD = new TH2F(hname,hname+";D Y;e eta",200,-10,10,200,-10,10);
438 hname = kqEtaRangeLabel[icut]+"fEtaCorrGDp_"+kqTypeLabel[kQuark];
439 fHistComm[iq][icut].fEtaCorrGDp = new TH2F(hname,hname+";D Y;e eta",200,-10,10,200,-10,10);
440 hname = kqEtaRangeLabel[icut]+"fEtaCorrGD0_"+kqTypeLabel[kQuark];
441 fHistComm[iq][icut].fEtaCorrGD0 = new TH2F(hname,hname+";D Y;e eta",200,-10,10,200,-10,10);
442 hname = kqEtaRangeLabel[icut]+"fEtaCorrGDrest_"+kqTypeLabel[kQuark];
443 fHistComm[iq][icut].fEtaCorrGDrest = new TH2F(hname,hname+";D Y;e eta",200,-10,10,200,-10,10);
445 hname = kqEtaRangeLabel[icut]+"fEtaCorrB_"+kqTypeLabel[kQuark];
446 fHistComm[iq][icut].fEtaCorrB = new TH2F(hname,hname+";B Y;e eta",200,-10,10,200,-10,10);
447 hname = kqEtaRangeLabel[icut]+"fEtaCorrGB_"+kqTypeLabel[kQuark];
448 fHistComm[iq][icut].fEtaCorrGB = new TH2F(hname,hname+";B Y;e eta",200,-10,10,200,-10,10);
450 hname = kqEtaRangeLabel[icut]+"PtCorrBinein_"+kqTypeLabel[kQuark];
451 fHistComm[iq][icut].fPtCorrBinein = new TH2F(hname,hname+";B p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
452 hname = kqEtaRangeLabel[icut]+"PtCorrBineout_"+kqTypeLabel[kQuark];
453 fHistComm[iq][icut].fPtCorrBineout = new TH2F(hname,hname+";B p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
454 hname = kqEtaRangeLabel[icut]+"PtCorrBoutein_"+kqTypeLabel[kQuark];
455 fHistComm[iq][icut].fPtCorrBoutein = new TH2F(hname,hname+";B p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
456 hname = kqEtaRangeLabel[icut]+"PtCorrBouteout_"+kqTypeLabel[kQuark];
457 fHistComm[iq][icut].fPtCorrBouteout = new TH2F(hname,hname+";B p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
459 if(fQAhistos) fHistComm[iq][icut].FillList(fQAhistos);
463 hname = kqEtaRangeLabel[0]+"Nq_"+kqTypeLabel[kQuark];
464 fHistComm[iq][0].fNq = new TH1F(hname,hname,50,-0.5,49.5);
465 hname = kqEtaRangeLabel[0]+"ProcessID_"+kqTypeLabel[kQuark];
466 fHistComm[iq][0].fProcessID = new TH1F(hname,hname,21,-10.5,10.5);
470 //__________________________________________
471 void AliHFEmcQA::Init()
473 // called at begining every event
475 for (Int_t i=0; i<2; i++){
481 fParentSelect[0][0] = 411; //D+
482 fParentSelect[0][1] = 421; //D0
483 fParentSelect[0][2] = 431; //Ds+
484 fParentSelect[0][3] = 4122; //Lambdac+
485 fParentSelect[0][4] = 4132; //Ksic0
486 fParentSelect[0][5] = 4232; //Ksic+
487 fParentSelect[0][6] = 4332; //OmegaC0
489 fParentSelect[1][0] = 511; //B0
490 fParentSelect[1][1] = 521; //B+
491 fParentSelect[1][2] = 531; //Bs0
492 fParentSelect[1][3] = 5122; //Lambdab0
493 fParentSelect[1][4] = 5132; //Ksib-
494 fParentSelect[1][5] = 5232; //Ksib0
495 fParentSelect[1][6] = 5332; //Omegab-
500 //__________________________________________
501 void AliHFEmcQA::GetMesonKine()
504 // get meson pt spectra
507 AliVParticle *mctrack2 = NULL;
508 AliMCParticle *mctrack0 = NULL;
509 AliVParticle *mctrackdaugt= NULL;
510 AliMCParticle *mctrackd= NULL;
514 if(fCentrality>=11) {
515 AliWarning(Form("Centrality out of histogram array limits: %d", fCentrality));
518 if(!fIsPbPb&&!fIsppMultiBin) fCentrality=0;
520 for(Int_t imc = 0; imc <fMCEvent->GetNumberOfPrimaries(); imc++){
521 if(!(mctrack2 = fMCEvent->GetTrack(imc))) continue;
522 TParticle* mcpart0 = fMCEvent->Stack()->Particle(imc);
523 if(!mcpart0) continue;
524 mctrack0 = dynamic_cast<AliMCParticle *>(mctrack2);
525 if(!mctrack0) continue;
527 // if(!fIsPbPb&&!fIsppMultiBin) fCentrality=0;
528 Float_t mcsource = 1;
529 if(fGetWeightHist) mcsource = GetElecSource(mctrack0, kFALSE);
531 if(TMath::Abs(mctrack0->PdgCode()) == 111) // pi0
533 if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
534 if(mcpart0->IsPrimary()){
535 fMCQACollection->Fill(Form("pionspectraPrimary_centrbin%i",fCentrality),mctrack0->Pt());
537 fMCQACollection->Fill(Form("pionspectra_centrbin%i",fCentrality),mctrack0->Pt());
538 fMCQACollection->Fill(Form("pionspectraLog_centrbin%i",fCentrality),mctrack0->Pt());
539 fMCQACollection->Fill(Form("pionspectraLog2D_centrbin%i",fCentrality),mcsource,mctrack0->Pt());
540 if(imc>fMCEvent->GetNumberOfPrimaries()) {
541 fMCQACollection->Fill(Form("pionspectrasec_centrbin%i",fCentrality),mctrack0->Pt());
542 fMCQACollection->Fill(Form("pionspectraLogsec_centrbin%i",fCentrality),mctrack0->Pt());
545 fMCQACollection->Fill(Form("pionspectrapri_centrbin%i",fCentrality),mctrack0->Pt());
546 fMCQACollection->Fill(Form("pionspectraLogpri_centrbin%i",fCentrality),mctrack0->Pt());
549 id1=mctrack0->GetFirstDaughter();
550 id2=mctrack0->GetLastDaughter();
551 if(!((id2-id1)==2)) continue;
552 for(int idx=id1; idx<=id2; idx++){
553 if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
554 if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
555 if(TMath::Abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
556 fMCQACollection->Fill(Form("piondaughters_centrbin%i",fCentrality),mctrackd->Pt());
559 else if(TMath::Abs(mctrack0->PdgCode()) == 221) // eta
561 if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
562 if(mcpart0->IsPrimary()){
563 fMCQACollection->Fill(Form("etaspectraPrimary_centrbin%i",fCentrality),mctrack0->Pt());
565 fMCQACollection->Fill(Form("etaspectra_centrbin%i",fCentrality),mctrack0->Pt());
566 fMCQACollection->Fill(Form("etaspectraLog_centrbin%i",fCentrality),mctrack0->Pt());
567 fMCQACollection->Fill(Form("etaspectraLog2D_centrbin%i",fCentrality),mcsource,mctrack0->Pt());
568 if(imc>fMCEvent->GetNumberOfPrimaries()) {
569 fMCQACollection->Fill(Form("etaspectrasec_centrbin%i",fCentrality),mctrack0->Pt());
570 fMCQACollection->Fill(Form("etaspectraLogsec_centrbin%i",fCentrality),mctrack0->Pt());
573 fMCQACollection->Fill(Form("etaspectrapri_centrbin%i",fCentrality),mctrack0->Pt());
574 fMCQACollection->Fill(Form("etaspectraLogpri_centrbin%i",fCentrality),mctrack0->Pt());
577 id1=mctrack0->GetFirstDaughter();
578 id2=mctrack0->GetLastDaughter();
579 if(!((id2-id1)==2||(id2-id1)==3)) continue;
580 for(int idx=id1; idx<=id2; idx++){
581 if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
582 if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
583 if(TMath::Abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
584 fMCQACollection->Fill(Form("etadaughters_centrbin%i",fCentrality),mctrackd->Pt());
587 else if(TMath::Abs(mctrack0->PdgCode()) == 223) // omega
589 if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
590 fMCQACollection->Fill(Form("omegaspectra_centrbin%i",fCentrality),mctrack0->Pt());
591 fMCQACollection->Fill(Form("omegaspectraLog_centrbin%i",fCentrality),mctrack0->Pt());
592 fMCQACollection->Fill(Form("omegaspectraLog2D_centrbin%i",fCentrality),mcsource,mctrack0->Pt());
594 id1=mctrack0->GetFirstDaughter();
595 id2=mctrack0->GetLastDaughter();
596 if(!((id2-id1)==1||(id2-id1)==2)) continue;
597 for(int idx=id1; idx<=id2; idx++){
598 if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
599 if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
600 if(TMath::Abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
601 fMCQACollection->Fill(Form("omegadaughters_centrbin%i",fCentrality),mctrackd->Pt());
604 else if(TMath::Abs(mctrack0->PdgCode()) == 333) // phi
606 if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
607 fMCQACollection->Fill(Form("phispectra_centrbin%i",fCentrality),mctrack0->Pt());
608 fMCQACollection->Fill(Form("phispectraLog_centrbin%i",fCentrality),mctrack0->Pt());
609 fMCQACollection->Fill(Form("phispectraLog2D_centrbin%i",fCentrality),mcsource,mctrack0->Pt());
611 id1=mctrack0->GetFirstDaughter();
612 id2=mctrack0->GetLastDaughter();
613 if(!((id2-id1)==1)) continue;
614 for(int idx=id1; idx<=id2; idx++){
615 if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
616 if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
617 if(TMath::Abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
618 fMCQACollection->Fill(Form("phidaughters_centrbin%i",fCentrality),mctrackd->Pt());
621 else if(TMath::Abs(mctrack0->PdgCode()) == 331) // eta prime
623 if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
624 fMCQACollection->Fill(Form("etapspectra_centrbin%i",fCentrality),mctrack0->Pt());
625 fMCQACollection->Fill(Form("etapspectraLog_centrbin%i",fCentrality),mctrack0->Pt());
626 fMCQACollection->Fill(Form("etapspectraLog2D_centrbin%i",fCentrality),mcsource,mctrack0->Pt());
628 id1=mctrack0->GetFirstDaughter();
629 id2=mctrack0->GetLastDaughter();
630 if(!((id2-id1)==2||(id2-id1)==3)) continue;
631 for(int idx=id1; idx<=id2; idx++){
632 if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
633 if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
634 if(TMath::Abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
635 fMCQACollection->Fill(Form("etapdaughters_centrbin%i",fCentrality),mctrackd->Pt());
638 else if(TMath::Abs(mctrack0->PdgCode()) == 113) // rho
640 if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
641 fMCQACollection->Fill(Form("rhospectra_centrbin%i",fCentrality),mctrack0->Pt());
642 fMCQACollection->Fill(Form("rhospectraLog_centrbin%i",fCentrality),mctrack0->Pt());
643 fMCQACollection->Fill(Form("rhospectraLog2D_centrbin%i",fCentrality),mcsource,mctrack0->Pt());
645 id1=mctrack0->GetFirstDaughter();
646 id2=mctrack0->GetLastDaughter();
647 if(!((id2-id1)==1)) continue;
648 for(int idx=id1; idx<=id2; idx++){
649 if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
650 if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
651 if(TMath::Abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
652 fMCQACollection->Fill(Form("rhodaughters_centrbin%i",fCentrality),mctrackd->Pt());
655 else if(TMath::Abs(mctrack0->PdgCode()) == 321) // kaon+-
657 if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
658 fMCQACollection->Fill(Form("kaonspectraLog_centrbin%i",fCentrality),mctrack0->Pt());
661 else if(TMath::Abs(mctrack0->PdgCode()) == 130) // k0L
663 if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
664 fMCQACollection->Fill(Form("k0LspectraLog_centrbin%i",fCentrality),mctrack0->Pt());
670 //__________________________________________
671 void AliHFEmcQA::GetQuarkKine(TParticle *part, Int_t iTrack, const Int_t kquark)
673 // get heavy quark kinematics
675 if (kquark != kCharm && kquark != kBeauty) {
676 AliDebug(1, "This task is only for heavy quark QA, return\n");
679 Int_t iq = kquark - kCharm;
681 if (iTrack < 0 || !part) {
682 AliDebug(1, "Stack label is negative or no mcparticle, return\n");
686 if(!fIsPbPb&&!fIsppMultiBin) fCentrality=0;
688 AliMCParticle *mctrack = NULL;
689 Int_t partPdgcode = TMath::Abs(part->GetPdgCode());
691 // select heavy hadron or not fragmented heavy quark
692 if ( int(partPdgcode/100.)==kquark || int(partPdgcode/1000.)==kquark || (partPdgcode==kquark && (part->GetNDaughters()==0 && iTrack>5)) ){
694 TParticle *partMother;
697 if (partPdgcode == kquark){ // in case of not fragmented heavy quark
700 } else{ // in case of heavy hadron, start to search for mother heavy parton
701 iLabel = part->GetFirstMother();
702 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return;
703 if (iLabel>-1) { partMother = mctrack->Particle(); }
705 AliDebug(1, "Stack label is negative, return\n");
710 // heavy parton selection as a mother of heavy hadron
711 // if the heavy particle comes from string which is denoted as particle status 12|12|12...12|11,[PYTHIA p.60]
712 // in this case, the mother of heavy particle can be one of the fragmented parton of the string
713 // should I make a condition that partMother should be quark or diquark? -> not necessary
714 if ( TMath::Abs(partMother->GetPdgCode()) == kquark || (partMother->GetStatusCode() == 11) ){
715 //if ( TMath::Abs(partMother->GetPdgCode()) == kquark || (partMother->GetStatusCode() == 11 || partMother->GetStatusCode() == 12) ){
717 if ( TMath::Abs(partMother->GetPdgCode()) != kquark ){
718 // search fragmented partons in the same string
719 Bool_t isSameString = kTRUE;
720 for (Int_t i=1; i<fgkMaxIter; i++){
722 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return;
723 if (iLabel>-1) { partMother = mctrack->Particle(); }
725 AliDebug(1, "Stack label is negative, return\n");
728 if ( TMath::Abs(partMother->GetPdgCode()) == kquark ) break;
729 if ( partMother->GetStatusCode() != 12 ) isSameString = kFALSE;
730 if (!isSameString) return;
733 AliDebug(1, "Can not find heavy parton of this heavy hadron in the string, return\n");
734 if (TMath::Abs(partMother->GetPdgCode()) != kquark) return;
736 if (fIsHeavy[iq] >= 50) return;
737 fHeavyQuark[fIsHeavy[iq]] = partMother;
740 // fill kinematics for heavy parton
741 if (partMother->GetPdgCode() > 0) { // quark
742 fHist[iq][kQuark][0][fCentrality].fPdgCode->Fill(partMother->GetPdgCode());
743 fHist[iq][kQuark][0][fCentrality].fPt->Fill(partMother->Pt());
744 fHist[iq][kQuark][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMother));
745 fHist[iq][kQuark][0][fCentrality].fEta->Fill(partMother->Eta());
747 fHist[iq][kantiQuark][0][fCentrality].fPdgCode->Fill(partMother->GetPdgCode());
748 fHist[iq][kantiQuark][0][fCentrality].fPt->Fill(partMother->Pt());
749 fHist[iq][kantiQuark][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMother));
750 fHist[iq][kantiQuark][0][fCentrality].fEta->Fill(partMother->Eta());
753 } // end of heavy parton slection loop
755 } // end of heavy hadron or quark selection
759 //__________________________________________
760 void AliHFEmcQA::EndOfEventAna(const Int_t kquark)
762 // end of event analysis
764 if (kquark != kCharm && kquark != kBeauty) {
765 AliDebug(1, "This task is only for heavy quark QA, return\n");
768 Int_t iq = kquark - kCharm;
771 // # of heavy quark per event
772 AliDebug(1,Form("Number of heavy quark in this event = %d \n",fIsHeavy[iq]));
773 fHistComm[iq][0].fNq->Fill(fIsHeavy[iq]);
775 Int_t motherID[fgkMaxGener];
776 Int_t motherType[fgkMaxGener];
777 Int_t motherLabel[fgkMaxGener];
778 Int_t ancestorPdg[fgkMaxGener];
779 Int_t ancestorLabel[fgkMaxGener];
781 for (Int_t i = 0; i < fgkMaxGener; i++){ // initialization
786 ancestorLabel[i] = 0;
790 // check history of found heavy quarks
791 for (Int_t i = 0; i < fIsHeavy[iq]; i++){
793 if(!fHeavyQuark[i]) return;
795 ancestorLabel[0] = i;
796 ancestorPdg[0] = fHeavyQuark[i]->GetPdgCode();
797 ancestorLabel[1] = fHeavyQuark[i]->GetFirstMother();
799 AliDebug(1,Form("pdg code= %d\n",ancestorPdg[0]));
800 AliDebug(1,Form("ancestor label= %d\n",ancestorLabel[1]));
803 while (ancestorLabel[ig] != -1){
804 // in case there is mother, get mother's pdg code and grandmother's label
805 IdentifyMother(ancestorLabel[ig], ancestorPdg[ig], ancestorLabel[ig+1]);
806 // if mother is still heavy, find again mother's ancestor
807 if (ancestorPdg[ig-1] == ancestorPdg[ig]) {
809 continue; // if it is from same heavy
811 // if the heavy's mother is not heavy, check the mother's label to know if it comes from inital or final parton shower
812 if (IsFromInitialShower(ancestorLabel[ig],motherID[i],motherType[i],motherLabel[i])) break;
813 if (IsFromFinalParton(ancestorLabel[ig],motherID[i],motherType[i],motherLabel[i])) break;
814 // if it is not the above case, something is strange
815 ReportStrangeness(motherID[i],motherType[i],motherLabel[i]);
818 if (ancestorLabel[ig] == -1){ // from hard scattering
819 HardScattering(kquark, motherID[i],motherType[i], motherLabel[i]);
822 } // end of found heavy quark loop
825 // check process type
827 for (Int_t i = 0; i < fIsHeavy[iq]; i++){
828 AliDebug(1,Form("Mother ID= %d type= %d label= %d\n",motherID[i],motherType[i],motherLabel[i]));
832 Int_t nheavypair = Int_t(fIsHeavy[iq]/2.);
833 for (Int_t ipair = 0; ipair < nheavypair; ipair++){
836 Int_t id2 = ipair*2 + 1;
838 if (motherType[id1] == 2 && motherType[id2] == 2){
839 if (motherLabel[id1] == motherLabel[id2]) processID = kGluonSplitting; // gluon spliting
842 else if (motherType[id1] == -1 && motherType[id2] == -1) {
843 if (motherLabel[id1] == -1 && motherLabel[id2] == -1) {
844 if (motherID[id1] == fgkGluon) processID = kPairCreationFromg; // gluon fusion
845 else processID = kPairCreationFromq; // q-qbar pair creation
849 else if (motherType[id1] == -1 || motherType[id2] == -1) {
850 if ((motherLabel[id1] == -1 || motherLabel[id2] == -1) && (motherLabel[id1]*motherLabel[id2] == -2 || motherLabel[id1]*motherLabel[id2] == -3)) {
851 if(motherID[id1]*motherID[id2] == kquark*fgkGluon) processID = kFlavourExitation; // flavour exitation
852 else processID = kLightQuarkShower;
856 else if (motherType[id1] == -2 || motherType[id2] == -2) {
857 if (motherLabel[id1] == motherLabel[id2]) processID = kInitialPartonShower; // initial parton shower
863 if (nheavypair >1) AliDebug(1,Form("Multi pair found : process ID = %d\n",processID));
864 else fHistComm[iq][0].fProcessID->Fill(processID);
865 AliDebug(1,Form("Process ID = %d\n",processID));
866 } // end of # heavy quark pair loop
870 //__________________________________________
871 void AliHFEmcQA::GetHadronKine(TParticle* mcpart, const Int_t kquark)
873 // decay electron kinematics
875 if (kquark != kCharm && kquark != kBeauty) {
876 AliDebug(1, "This task is only for heavy quark QA, return\n");
879 Int_t iq = kquark - kCharm;
882 AliDebug(1, "no mc particle, return\n");
886 Int_t iLabel = mcpart->GetFirstMother();
888 AliDebug(1, "Stack label is negative, return\n");
892 if(!fIsPbPb&&!fIsppMultiBin) fCentrality=0;
894 TParticle *partCopy = mcpart;
895 Int_t pdgcode = mcpart->GetPdgCode();
896 Int_t pdgcodeCopy = pdgcode;
898 AliMCParticle *mctrack = NULL;
900 // if the mother is charmed hadron
901 Bool_t isDirectCharm = kFALSE;
902 if ( int(TMath::Abs(pdgcode)/100.) == kCharm || int(TMath::Abs(pdgcode)/1000.) == kCharm ) {
904 // iterate until you find B hadron as a mother or become top ancester
905 for (Int_t i=1; i<fgkMaxIter; i++){
907 Int_t jLabel = mcpart->GetFirstMother();
909 isDirectCharm = kTRUE;
910 break; // if there is no ancester
912 if (jLabel < 0){ // safety protection
913 AliDebug(1, "Stack label is negative, return\n");
916 // if there is an ancester
917 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return;
918 TParticle* mother = mctrack->Particle();
919 Int_t motherPDG = mother->GetPdgCode();
921 for (Int_t j=0; j<fNparents; j++){
922 if (TMath::Abs(motherPDG)==fParentSelect[1][j]) return; // return if this hadron is originated from b
926 } // end of iteration
928 if((isDirectCharm == kTRUE && kquark == kCharm) || kquark == kBeauty) {
929 for (Int_t i=0; i<fNparents; i++){
930 if (TMath::Abs(pdgcodeCopy)==fParentSelect[iq][i]){
932 // fill hadron kinematics
933 fHist[iq][kHadron][0][fCentrality].fPdgCode->Fill(pdgcodeCopy);
934 fHist[iq][kHadron][0][fCentrality].fPt->Fill(partCopy->Pt());
935 fHist[iq][kHadron][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partCopy));
936 fHist[iq][kHadron][0][fCentrality].fEta->Fill(partCopy->Eta());
939 fhD[i]->Fill(partCopy->Pt(),AliHFEtools::GetRapidity(partCopy));
943 // I also want to store D* info to compare with D* measurement
944 if (TMath::Abs(pdgcodeCopy)==413 && iq==0) { //D*+
945 fhD[7]->Fill(partCopy->Pt(),AliHFEtools::GetRapidity(partCopy));
947 if (TMath::Abs(pdgcodeCopy)==423 && iq==0) { //D*0
948 fhD[8]->Fill(partCopy->Pt(),AliHFEtools::GetRapidity(partCopy));
953 //__________________________________________
954 void AliHFEmcQA::GetDecayedKine(TParticle* mcpart, const Int_t kquark, Int_t kdecayed)
956 // decay electron kinematics
958 if (!(kquark == kCharm || kquark == kBeauty || kquark == kOthers)){
959 AliDebug(1, "This task is only for heavy quark QA, return\n");
962 Int_t iq = kquark - kCharm;
963 Bool_t isFinalOpenCharm = kFALSE;
966 AliDebug(1, "no mcparticle, return\n");
970 if(!fIsPbPb&&!fIsppMultiBin) fCentrality=0;
972 Double_t eabsEta = TMath::Abs(mcpart->Eta());
973 Double_t eabsY = TMath::Abs(AliHFEtools::GetRapidity(mcpart));
977 if ( TMath::Abs(mcpart->GetPdgCode()) != kdecayed ) esource = kMisID-4;
978 else esource =GetSource(mcpart)-4; // return for the cases kGamma=4, kPi0=5, kElse=6
979 if(esource==0|| esource==1 || esource==2 || esource==3){
980 fHist[iq][esource][0][fCentrality].fPt->Fill(mcpart->Pt());
981 fHist[iq][esource][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
982 fHist[iq][esource][0][fCentrality].fEta->Fill(mcpart->Eta());
984 fHist[iq][esource][1][fCentrality].fPt->Fill(mcpart->Pt());
985 fHist[iq][esource][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
986 fHist[iq][esource][1][fCentrality].fEta->Fill(mcpart->Eta());
989 fHist[iq][esource][2][fCentrality].fPt->Fill(mcpart->Pt());
990 fHist[iq][esource][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
991 fHist[iq][esource][2][fCentrality].fEta->Fill(mcpart->Eta());
996 AliDebug(1, "e source is out of defined ranges, return\n");
1001 if ( TMath::Abs(mcpart->GetPdgCode()) != kdecayed ) return;
1003 Int_t iLabel = mcpart->GetFirstMother();
1005 AliDebug(1, "Stack label is negative, return\n");
1009 AliMCParticle *mctrack = NULL;
1010 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return;
1011 TParticle *partMother = mctrack->Particle();
1012 TParticle *partMotherCopy = partMother;
1013 Int_t maPdgcode = partMother->GetPdgCode();
1014 Int_t maPdgcodeCopy = maPdgcode;
1016 // get mc primary vertex
1019 if(fMCHeader) fMCHeader->PrimaryVertex(mcPrimVtx);
1021 // get electron production vertex
1022 TLorentzVector ePoint;
1023 mcpart->ProductionVertex(ePoint);
1025 // calculated production vertex to primary vertex (in xy plane)
1026 Float_t decayLxy = TMath::Sqrt((mcPrimVtx[0]-ePoint[0])*(mcPrimVtx[0]-ePoint[0])+(mcPrimVtx[1]-ePoint[1])*(mcPrimVtx[1]-ePoint[1]));
1028 Float_t decayLxy = 0;
1030 // if the mother is charmed hadron
1031 Bool_t isMotherDirectCharm = kFALSE;
1032 if ( int(TMath::Abs(maPdgcode)/100.) == kCharm || int(TMath::Abs(maPdgcode)/1000.) == kCharm ) {
1034 for (Int_t i=0; i<fNparents; i++){
1035 if (TMath::Abs(maPdgcode)==fParentSelect[0][i]){
1036 isFinalOpenCharm = kTRUE;
1039 if (!isFinalOpenCharm) return ;
1041 // iterate until you find B hadron as a mother or become top ancester
1042 for (Int_t i=1; i<fgkMaxIter; i++){
1044 Int_t jLabel = partMother->GetFirstMother();
1046 isMotherDirectCharm = kTRUE;
1047 break; // if there is no ancester
1049 if (jLabel < 0){ // safety protection
1050 AliDebug(1, "Stack label is negative, return\n");
1054 // if there is an ancester
1055 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return;
1056 TParticle* grandMa = mctrack->Particle();
1057 Int_t grandMaPDG = grandMa->GetPdgCode();
1059 for (Int_t j=0; j<fNparents; j++){
1060 if (TMath::Abs(grandMaPDG)==fParentSelect[1][j]){
1062 if (kquark == kCharm) return;
1063 // fill electron kinematics
1064 fHist[iq][kElectron2nd][0][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1065 fHist[iq][kElectron2nd][0][fCentrality].fPt->Fill(mcpart->Pt());
1066 fHist[iq][kElectron2nd][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1067 fHist[iq][kElectron2nd][0][fCentrality].fEta->Fill(mcpart->Eta());
1069 // fill mother hadron kinematics
1070 fHist[iq][kDeHadron][0][fCentrality].fPdgCode->Fill(grandMaPDG);
1071 fHist[iq][kDeHadron][0][fCentrality].fPt->Fill(grandMa->Pt());
1072 fHist[iq][kDeHadron][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(grandMa));
1073 fHist[iq][kDeHadron][0][fCentrality].fEta->Fill(grandMa->Eta());
1076 fHist[iq][kElectron2nd][1][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1077 fHist[iq][kElectron2nd][1][fCentrality].fPt->Fill(mcpart->Pt());
1078 fHist[iq][kElectron2nd][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1079 fHist[iq][kElectron2nd][1][fCentrality].fEta->Fill(mcpart->Eta());
1081 // fill mother hadron kinematics
1082 fHist[iq][kDeHadron][1][fCentrality].fPdgCode->Fill(grandMaPDG);
1083 fHist[iq][kDeHadron][1][fCentrality].fPt->Fill(grandMa->Pt());
1084 fHist[iq][kDeHadron][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(grandMa));
1085 fHist[iq][kDeHadron][1][fCentrality].fEta->Fill(grandMa->Eta());
1089 fHist[iq][kElectron2nd][2][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1090 fHist[iq][kElectron2nd][2][fCentrality].fPt->Fill(mcpart->Pt());
1091 fHist[iq][kElectron2nd][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1092 fHist[iq][kElectron2nd][2][fCentrality].fEta->Fill(mcpart->Eta());
1094 // fill mother hadron kinematics
1095 fHist[iq][kDeHadron][2][fCentrality].fPdgCode->Fill(grandMaPDG);
1096 fHist[iq][kDeHadron][2][fCentrality].fPt->Fill(grandMa->Pt());
1097 fHist[iq][kDeHadron][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(grandMa));
1098 fHist[iq][kDeHadron][2][fCentrality].fEta->Fill(grandMa->Eta());
1101 //mj: to calculate B to e eta correlation to calculate total heavy quark cross section
1102 Int_t kLabel0 = grandMa->GetFirstMother();
1103 Bool_t isGGrandmaYes = kFALSE;
1104 Double_t ggmrapidwstmp=0;
1105 if (!(kLabel0 < 0)){ // safety protection
1106 if((mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(kLabel0))))){
1107 TParticle* ggrandMatmp = mctrack->Particle();
1108 Int_t ggrandMaPDGtmp = ggrandMatmp->GetPdgCode();
1109 if ( int(TMath::Abs(ggrandMaPDGtmp)/100.) == kBeauty || int(TMath::Abs(ggrandMaPDGtmp)/1000.) == kBeauty) isGGrandmaYes = kTRUE;
1110 ggmrapidwstmp = AliHFEtools::GetRapidity(ggrandMatmp);
1114 Double_t gmrapidwstmp0 = AliHFEtools::GetRapidity(grandMa);
1115 Double_t eetawstmp0 = mcpart->Eta();
1117 Double_t gmrapidtmp0 = TMath::Abs(gmrapidwstmp0);
1118 Double_t eetatmp0 = TMath::Abs(eetawstmp0);
1120 fHistComm[iq][0].fEtaCorrB->Fill(gmrapidwstmp0,eetawstmp0);
1121 if(isGGrandmaYes) fHistComm[iq][0].fEtaCorrGB->Fill(ggmrapidwstmp,eetawstmp0);
1122 else fHistComm[iq][0].fEtaCorrGB->Fill(gmrapidwstmp0,eetawstmp0);
1124 if(gmrapidtmp0<0.5 && eetatmp0<0.5 ) fHistComm[iq][0].fPtCorrBinein->Fill(grandMa->Pt(),mcpart->Pt());
1125 else if(gmrapidtmp0<0.5 && eetatmp0>0.5 ) fHistComm[iq][0].fPtCorrBineout->Fill(grandMa->Pt(),mcpart->Pt());
1126 else if(gmrapidtmp0>0.5 && eetatmp0<0.5 ) fHistComm[iq][0].fPtCorrBoutein->Fill(grandMa->Pt(),mcpart->Pt());
1127 else if(gmrapidtmp0>0.5 && eetatmp0>0.5 ) fHistComm[iq][0].fPtCorrBouteout->Fill(grandMa->Pt(),mcpart->Pt());
1128 //======================================================================================
1130 // ratio between pT of electron and pT of mother B hadron
1132 fHistComm[iq][0].fDePtRatio->Fill(grandMa->Pt(),mcpart->Pt()/grandMa->Pt());
1134 fHistComm[iq][1].fDePtRatio->Fill(grandMa->Pt(),mcpart->Pt()/grandMa->Pt());
1137 fHistComm[iq][2].fDePtRatio->Fill(grandMa->Pt(),mcpart->Pt()/grandMa->Pt());
1141 // distance between electron production point and primary vertex
1142 fHistComm[iq][0].fDeDistance->Fill(grandMa->Pt(),decayLxy);
1144 fHistComm[iq][1].fDeDistance->Fill(grandMa->Pt(),decayLxy);
1147 fHistComm[iq][2].fDeDistance->Fill(grandMa->Pt(),decayLxy);
1153 partMother = grandMa;
1154 } // end of iteration
1156 if((isMotherDirectCharm == kTRUE && kquark == kCharm) || kquark == kBeauty) {
1157 for (Int_t i=0; i<fNparents; i++){
1158 if (TMath::Abs(maPdgcodeCopy)==fParentSelect[iq][i]){
1160 fHist[iq][kElectron][0][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1161 fHist[iq][kElectron][0][fCentrality].fPt->Fill(mcpart->Pt());
1162 fHist[iq][kElectron][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1163 fHist[iq][kElectron][0][fCentrality].fEta->Fill(mcpart->Eta());
1165 // fill mother hadron kinematics
1166 fHist[iq][keHadron][0][fCentrality].fPdgCode->Fill(maPdgcodeCopy);
1167 fHist[iq][keHadron][0][fCentrality].fPt->Fill(partMotherCopy->Pt());
1168 fHist[iq][keHadron][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
1169 fHist[iq][keHadron][0][fCentrality].fEta->Fill(partMotherCopy->Eta());
1172 fHist[iq][kElectron][1][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1173 fHist[iq][kElectron][1][fCentrality].fPt->Fill(mcpart->Pt());
1174 fHist[iq][kElectron][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1175 fHist[iq][kElectron][1][fCentrality].fEta->Fill(mcpart->Eta());
1177 // fill mother hadron kinematics
1178 fHist[iq][keHadron][1][fCentrality].fPdgCode->Fill(maPdgcodeCopy);
1179 fHist[iq][keHadron][1][fCentrality].fPt->Fill(partMotherCopy->Pt());
1180 fHist[iq][keHadron][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
1181 fHist[iq][keHadron][1][fCentrality].fEta->Fill(partMotherCopy->Eta());
1185 fHist[iq][kElectron][2][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1186 fHist[iq][kElectron][2][fCentrality].fPt->Fill(mcpart->Pt());
1187 fHist[iq][kElectron][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1188 fHist[iq][kElectron][2][fCentrality].fEta->Fill(mcpart->Eta());
1190 // fill mother hadron kinematics
1191 fHist[iq][keHadron][2][fCentrality].fPdgCode->Fill(maPdgcodeCopy);
1192 fHist[iq][keHadron][2][fCentrality].fPt->Fill(partMotherCopy->Pt());
1193 fHist[iq][keHadron][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
1194 fHist[iq][keHadron][2][fCentrality].fEta->Fill(partMotherCopy->Eta());
1197 // ratio between pT of electron and pT of mother B or direct D hadron
1198 if(partMotherCopy->Pt()) {
1199 fHistComm[iq][0].fePtRatio->Fill(partMotherCopy->Pt(),mcpart->Pt()/partMotherCopy->Pt());
1201 fHistComm[iq][1].fePtRatio->Fill(partMotherCopy->Pt(),mcpart->Pt()/partMotherCopy->Pt());
1204 fHistComm[iq][2].fePtRatio->Fill(partMotherCopy->Pt(),mcpart->Pt()/partMotherCopy->Pt());
1207 fHistComm[iq][0].fPtCorr->Fill(partMotherCopy->Pt(),mcpart->Pt());
1209 fHistComm[iq][1].fPtCorr->Fill(partMotherCopy->Pt(),mcpart->Pt());
1212 fHistComm[iq][2].fPtCorr->Fill(partMotherCopy->Pt(),mcpart->Pt());
1214 if(TMath::Abs(partMotherCopy->GetPdgCode())==411) {
1215 fHistComm[iq][0].fPtCorrDp->Fill(partMotherCopy->Pt(),mcpart->Pt());
1217 fHistComm[iq][1].fPtCorrDp->Fill(partMotherCopy->Pt(),mcpart->Pt());
1220 fHistComm[iq][2].fPtCorrDp->Fill(partMotherCopy->Pt(),mcpart->Pt());
1223 else if(TMath::Abs(partMotherCopy->GetPdgCode())==421) {
1224 fHistComm[iq][0].fPtCorrD0->Fill(partMotherCopy->Pt(),mcpart->Pt());
1226 fHistComm[iq][1].fPtCorrD0->Fill(partMotherCopy->Pt(),mcpart->Pt());
1229 fHistComm[iq][2].fPtCorrD0->Fill(partMotherCopy->Pt(),mcpart->Pt());
1233 fHistComm[iq][0].fPtCorrDrest->Fill(partMotherCopy->Pt(),mcpart->Pt());
1235 fHistComm[iq][1].fPtCorrDrest->Fill(partMotherCopy->Pt(),mcpart->Pt());
1238 fHistComm[iq][2].fPtCorrDrest->Fill(partMotherCopy->Pt(),mcpart->Pt());
1242 //mj: to calculate D to e eta correlation to calculate total heavy quark cross section
1243 Int_t kLabel = partMotherCopy->GetFirstMother();
1244 Bool_t isGrandmaYes = kFALSE;
1245 Double_t gmrapidwstmp =0;
1246 if (!(kLabel < 0)){ // safety protection
1247 if((mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(kLabel))))){
1248 TParticle* grandMatmp = mctrack->Particle();
1249 Int_t grandMaPDGtmp = grandMatmp->GetPdgCode();
1250 if ( int(TMath::Abs(grandMaPDGtmp)/100.) == kCharm || int(TMath::Abs(grandMaPDGtmp)/1000.) == kCharm ) isGrandmaYes = kTRUE;
1251 if ( int(TMath::Abs(grandMaPDGtmp)/100.) == kBeauty || int(TMath::Abs(grandMaPDGtmp)/1000.) == kBeauty) isGrandmaYes = kTRUE;
1252 gmrapidwstmp = AliHFEtools::GetRapidity(grandMatmp);
1256 Double_t mrapidwstmp = AliHFEtools::GetRapidity(partMotherCopy);
1257 Double_t eetawstmp = mcpart->Eta();
1259 Double_t mrapidtmp = TMath::Abs(mrapidwstmp);
1260 Double_t eetatmp = TMath::Abs(eetawstmp);
1262 fHistComm[iq][0].fEtaCorrD->Fill(mrapidwstmp,eetawstmp);
1263 if(isGrandmaYes) fHistComm[iq][0].fEtaCorrGD->Fill(gmrapidwstmp,eetawstmp);
1264 else fHistComm[iq][0].fEtaCorrGD->Fill(mrapidwstmp,eetawstmp);
1266 if(mrapidtmp<0.5 && eetatmp<0.5 ) fHistComm[iq][0].fPtCorrDinein->Fill(partMotherCopy->Pt(),mcpart->Pt());
1267 else if(mrapidtmp<0.5 && eetatmp>0.5 ) fHistComm[iq][0].fPtCorrDineout->Fill(partMotherCopy->Pt(),mcpart->Pt());
1268 else if(mrapidtmp>0.5 && eetatmp<0.5 ) fHistComm[iq][0].fPtCorrDoutein->Fill(partMotherCopy->Pt(),mcpart->Pt());
1269 else if(mrapidtmp>0.5 && eetatmp>0.5 ) fHistComm[iq][0].fPtCorrDouteout->Fill(partMotherCopy->Pt(),mcpart->Pt());
1270 if(TMath::Abs(partMotherCopy->GetPdgCode())==411) {
1271 fHistComm[iq][0].fEtaCorrDp->Fill(mrapidwstmp,eetawstmp);
1272 if(isGrandmaYes) fHistComm[iq][0].fEtaCorrGDp->Fill(gmrapidwstmp,eetawstmp);
1273 else fHistComm[iq][0].fEtaCorrGDp->Fill(mrapidwstmp,eetawstmp);
1274 if(mrapidtmp<0.5 && eetatmp<0.5 ) fHistComm[iq][0].fPtCorrDpDinein->Fill(partMotherCopy->Pt(),mcpart->Pt());
1275 else if(mrapidtmp<0.5 && eetatmp>0.5 ) fHistComm[iq][0].fPtCorrDpDineout->Fill(partMotherCopy->Pt(),mcpart->Pt());
1276 else if(mrapidtmp>0.5 && eetatmp<0.5 ) fHistComm[iq][0].fPtCorrDpDoutein->Fill(partMotherCopy->Pt(),mcpart->Pt());
1277 else if(mrapidtmp>0.5 && eetatmp>0.5 ) fHistComm[iq][0].fPtCorrDpDouteout->Fill(partMotherCopy->Pt(),mcpart->Pt());
1279 else if(TMath::Abs(partMotherCopy->GetPdgCode())==421) {
1280 fHistComm[iq][0].fEtaCorrD0->Fill(mrapidwstmp,eetawstmp);
1281 if(isGrandmaYes) fHistComm[iq][0].fEtaCorrGD0->Fill(gmrapidwstmp,eetawstmp);
1282 else fHistComm[iq][0].fEtaCorrGD0->Fill(mrapidwstmp,eetawstmp);
1283 if(mrapidtmp<0.5 && eetatmp<0.5 ) fHistComm[iq][0].fPtCorrD0Dinein->Fill(partMotherCopy->Pt(),mcpart->Pt());
1284 else if(mrapidtmp<0.5 && eetatmp>0.5 ) fHistComm[iq][0].fPtCorrD0Dineout->Fill(partMotherCopy->Pt(),mcpart->Pt());
1285 else if(mrapidtmp>0.5 && eetatmp<0.5 ) fHistComm[iq][0].fPtCorrD0Doutein->Fill(partMotherCopy->Pt(),mcpart->Pt());
1286 else if(mrapidtmp>0.5 && eetatmp>0.5 ) fHistComm[iq][0].fPtCorrD0Douteout->Fill(partMotherCopy->Pt(),mcpart->Pt());
1289 fHistComm[iq][0].fEtaCorrDrest->Fill(mrapidwstmp,eetawstmp);
1290 if(isGrandmaYes) fHistComm[iq][0].fEtaCorrGDrest->Fill(gmrapidwstmp,eetawstmp);
1291 else fHistComm[iq][0].fEtaCorrGDrest->Fill(mrapidwstmp,eetawstmp);
1292 if(mrapidtmp<0.5 && eetatmp<0.5 ) fHistComm[iq][0].fPtCorrDrestDinein->Fill(partMotherCopy->Pt(),mcpart->Pt());
1293 else if(mrapidtmp<0.5 && eetatmp>0.5 ) fHistComm[iq][0].fPtCorrDrestDineout->Fill(partMotherCopy->Pt(),mcpart->Pt());
1294 else if(mrapidtmp>0.5 && eetatmp<0.5 ) fHistComm[iq][0].fPtCorrDrestDoutein->Fill(partMotherCopy->Pt(),mcpart->Pt());
1295 else if(mrapidtmp>0.5 && eetatmp>0.5 ) fHistComm[iq][0].fPtCorrDrestDouteout->Fill(partMotherCopy->Pt(),mcpart->Pt());
1298 // distance between electron production point and primary vertex
1299 fHistComm[iq][0].feDistance->Fill(partMotherCopy->Pt(),decayLxy);
1301 fHistComm[iq][1].feDistance->Fill(partMotherCopy->Pt(),decayLxy);
1304 fHistComm[iq][2].feDistance->Fill(partMotherCopy->Pt(),decayLxy);
1311 //____________________________________________________________________
1312 void AliHFEmcQA::GetDecayedKine(AliAODMCParticle *mcpart, const Int_t kquark, Int_t kdecayed)
1314 // decay electron kinematics
1316 if (kquark != kCharm && kquark != kBeauty) {
1317 AliDebug(1, "This task is only for heavy quark QA, return\n");
1321 Int_t iq = kquark - kCharm;
1322 Bool_t isFinalOpenCharm = kFALSE;
1325 AliDebug(1, "no mcparticle, return\n");
1329 if ( TMath::Abs(mcpart->GetPdgCode()) != kdecayed ) return;
1331 Double_t eabsEta = TMath::Abs(mcpart->Eta());
1332 Double_t eabsY = TMath::Abs(AliHFEtools::GetRapidity(mcpart));
1335 Int_t iLabel = mcpart->GetMother();
1337 AliDebug(1, "Stack label is negative, return\n");
1341 if(!fIsPbPb&&!fIsppMultiBin) fCentrality=0;
1343 AliAODMCParticle *partMother = (AliAODMCParticle*)fMCArray->At(iLabel);
1344 AliAODMCParticle *partMotherCopy = partMother;
1345 Int_t maPdgcode = partMother->GetPdgCode();
1346 Int_t maPdgcodeCopy = maPdgcode;
1348 Bool_t isMotherDirectCharm = kFALSE;
1349 if ( int(TMath::Abs(maPdgcode)/100.) == kCharm || int(TMath::Abs(maPdgcode)/1000.) == kCharm ) {
1351 for (Int_t i=0; i<fNparents; i++){
1352 if (TMath::Abs(maPdgcode)==fParentSelect[0][i]){
1353 isFinalOpenCharm = kTRUE;
1356 if (!isFinalOpenCharm) return;
1358 for (Int_t i=1; i<fgkMaxIter; i++){
1360 Int_t jLabel = partMother->GetMother();
1362 isMotherDirectCharm = kTRUE;
1363 break; // if there is no ancester
1365 if (jLabel < 0){ // safety protection
1366 AliDebug(1, "Stack label is negative, return\n");
1370 // if there is an ancester
1371 AliAODMCParticle* grandMa = (AliAODMCParticle*)fMCArray->At(jLabel);
1372 Int_t grandMaPDG = grandMa->GetPdgCode();
1374 for (Int_t j=0; j<fNparents; j++){
1375 if (TMath::Abs(grandMaPDG)==fParentSelect[1][j]){
1377 if (kquark == kCharm) return;
1378 // fill electron kinematics
1379 fHist[iq][kElectron2nd][0][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1380 fHist[iq][kElectron2nd][0][fCentrality].fPt->Fill(mcpart->Pt());
1381 fHist[iq][kElectron2nd][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1382 fHist[iq][kElectron2nd][0][fCentrality].fEta->Fill(mcpart->Eta());
1384 // fill mother hadron kinematics
1385 fHist[iq][kDeHadron][0][fCentrality].fPdgCode->Fill(grandMaPDG);
1386 fHist[iq][kDeHadron][0][fCentrality].fPt->Fill(grandMa->Pt());
1387 fHist[iq][kDeHadron][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(grandMa));
1388 fHist[iq][kDeHadron][0][fCentrality].fEta->Fill(grandMa->Eta());
1391 // fill electron kinematics
1392 fHist[iq][kElectron2nd][1][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1393 fHist[iq][kElectron2nd][1][fCentrality].fPt->Fill(mcpart->Pt());
1394 fHist[iq][kElectron2nd][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1395 fHist[iq][kElectron2nd][1][fCentrality].fEta->Fill(mcpart->Eta());
1397 // fill mother hadron kinematics
1398 fHist[iq][kDeHadron][1][fCentrality].fPdgCode->Fill(grandMaPDG);
1399 fHist[iq][kDeHadron][1][fCentrality].fPt->Fill(grandMa->Pt());
1400 fHist[iq][kDeHadron][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(grandMa));
1401 fHist[iq][kDeHadron][1][fCentrality].fEta->Fill(grandMa->Eta());
1404 // fill electron kinematics
1405 fHist[iq][kElectron2nd][2][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1406 fHist[iq][kElectron2nd][2][fCentrality].fPt->Fill(mcpart->Pt());
1407 fHist[iq][kElectron2nd][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1408 fHist[iq][kElectron2nd][2][fCentrality].fEta->Fill(mcpart->Eta());
1410 // fill mother hadron kinematics
1411 fHist[iq][kDeHadron][2][fCentrality].fPdgCode->Fill(grandMaPDG);
1412 fHist[iq][kDeHadron][2][fCentrality].fPt->Fill(grandMa->Pt());
1413 fHist[iq][kDeHadron][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(grandMa));
1414 fHist[iq][kDeHadron][2][fCentrality].fEta->Fill(grandMa->Eta());
1421 partMother = grandMa;
1422 } // end of iteration
1424 if ((isMotherDirectCharm == kTRUE && kquark == kCharm) || kquark == kBeauty) {
1425 for (Int_t i=0; i<fNparents; i++){
1426 if (TMath::Abs(maPdgcodeCopy)==fParentSelect[iq][i]){
1428 // fill electron kinematics
1429 fHist[iq][kElectron][0][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1430 fHist[iq][kElectron][0][fCentrality].fPt->Fill(mcpart->Pt());
1431 fHist[iq][kElectron][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1432 fHist[iq][kElectron][0][fCentrality].fEta->Fill(mcpart->Eta());
1434 // fill mother hadron kinematics
1435 fHist[iq][keHadron][0][fCentrality].fPdgCode->Fill(maPdgcodeCopy);
1436 fHist[iq][keHadron][0][fCentrality].fPt->Fill(partMotherCopy->Pt());
1437 fHist[iq][keHadron][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
1438 fHist[iq][keHadron][0][fCentrality].fEta->Fill(partMotherCopy->Eta());
1441 // fill electron kinematics
1442 fHist[iq][kElectron][1][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1443 fHist[iq][kElectron][1][fCentrality].fPt->Fill(mcpart->Pt());
1444 fHist[iq][kElectron][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1445 fHist[iq][kElectron][1][fCentrality].fEta->Fill(mcpart->Eta());
1447 // fill mother hadron kinematics
1448 fHist[iq][keHadron][1][fCentrality].fPdgCode->Fill(maPdgcodeCopy);
1449 fHist[iq][keHadron][1][fCentrality].fPt->Fill(partMotherCopy->Pt());
1450 fHist[iq][keHadron][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
1451 fHist[iq][keHadron][1][fCentrality].fEta->Fill(partMotherCopy->Eta());
1454 // fill electron kinematics
1455 fHist[iq][kElectron][2][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1456 fHist[iq][kElectron][2][fCentrality].fPt->Fill(mcpart->Pt());
1457 fHist[iq][kElectron][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1458 fHist[iq][kElectron][2][fCentrality].fEta->Fill(mcpart->Eta());
1460 // fill mother hadron kinematics
1461 fHist[iq][keHadron][2][fCentrality].fPdgCode->Fill(maPdgcodeCopy);
1462 fHist[iq][keHadron][2][fCentrality].fPt->Fill(partMotherCopy->Pt());
1463 fHist[iq][keHadron][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
1464 fHist[iq][keHadron][2][fCentrality].fEta->Fill(partMotherCopy->Eta());
1473 //__________________________________________
1474 void AliHFEmcQA::IdentifyMother(Int_t motherlabel, Int_t &motherpdg, Int_t &grandmotherlabel)
1476 // find mother pdg code and label
1478 if (motherlabel < 0) {
1479 AliDebug(1, "Stack label is negative, return\n");
1482 AliMCParticle *mctrack = NULL;
1483 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(motherlabel))))) return;
1484 TParticle *heavysMother = mctrack->Particle();
1485 motherpdg = heavysMother->GetPdgCode();
1486 grandmotherlabel = heavysMother->GetFirstMother();
1487 AliDebug(1,Form("ancestor pdg code= %d\n",motherpdg));
1490 //__________________________________________
1491 void AliHFEmcQA::HardScattering(const Int_t kquark, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
1493 // mothertype -1 means this heavy quark coming from hard vertex
1495 AliMCParticle *mctrack1 = NULL;
1496 AliMCParticle *mctrack2 = NULL;
1497 if(!(mctrack1 = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(4))))) return;
1498 if(!(mctrack2 = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(5))))) return;
1499 TParticle *afterinitialrad1 = mctrack1->Particle();
1500 TParticle *afterinitialrad2 = mctrack2->Particle();
1504 if (TMath::Abs(afterinitialrad1->GetPdgCode()) == fgkGluon && TMath::Abs(afterinitialrad2->GetPdgCode()) == fgkGluon){
1505 AliDebug(1,"heavy from gluon gluon pair creation!\n");
1507 motherID = fgkGluon;
1509 else if (TMath::Abs(afterinitialrad1->GetPdgCode()) == kquark || TMath::Abs(afterinitialrad2->GetPdgCode()) == kquark){ // one from Q and the other from g
1510 AliDebug(1,"heavy from flavor exitation!\n");
1514 else if (TMath::Abs(afterinitialrad1->GetPdgCode()) == TMath::Abs(afterinitialrad2->GetPdgCode())){
1515 AliDebug(1,"heavy from q-qbar pair creation!\n");
1520 AliDebug(1,"something strange!\n");
1527 //__________________________________________
1528 Bool_t AliHFEmcQA::IsFromInitialShower(Int_t inputmotherlabel, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
1530 // mothertype -2 means this heavy quark coming from initial state
1532 AliMCParticle *mctrack = NULL;
1533 if (inputmotherlabel==2 || inputmotherlabel==3){ // mother exist before initial state radiation
1534 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(inputmotherlabel))))) return kFALSE;
1535 TParticle *heavysMother = mctrack->Particle();
1536 motherID = heavysMother->GetPdgCode();
1537 mothertype = -2; // there is mother before initial state radiation
1538 motherlabel = inputmotherlabel;
1539 AliDebug(1,"initial parton shower! \n");
1547 //__________________________________________
1548 Bool_t AliHFEmcQA::IsFromFinalParton(Int_t inputmotherlabel, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
1550 // mothertype 2 means this heavy quark coming from final state
1552 AliMCParticle *mctrack = NULL;
1553 if (inputmotherlabel > 5){ // mother exist after hard scattering
1554 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(inputmotherlabel))))) return kFALSE;
1555 TParticle *heavysMother = mctrack->Particle();
1556 motherID = heavysMother->GetPdgCode();
1558 motherlabel = inputmotherlabel;
1559 AliDebug(1,Form("heavy quark from %d after hard scattering! \n",motherID));
1566 //__________________________________________
1567 void AliHFEmcQA::ReportStrangeness(Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
1569 // mark strange behavior
1574 AliDebug(1,"something strange!\n");
1577 //__________________________________________
1578 Int_t AliHFEmcQA::GetSource(const AliVParticle* const mcpart) const
1580 // decay particle's origin
1582 //if ( TMath::Abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return -1;
1585 Bool_t isFinalOpenCharm = kFALSE;
1588 AliDebug(1, "Stack label is negative or no mcparticle, return\n");
1593 // Information not in the base class, cast necessary
1594 Int_t iLabel = GetMother(mcpart);
1596 AliDebug(1, "Stack label is negative, return\n");
1600 const AliVParticle *partMother = fMCEvent->GetTrack(iLabel);
1601 Int_t maPdgcode = partMother->PdgCode();
1603 // if the mother is charmed hadron
1604 if ( int(TMath::Abs(maPdgcode)/100.) == kCharm || int(TMath::Abs(maPdgcode)/1000.) == kCharm ) {
1606 for (Int_t i=0; i<fNparents; i++){
1607 if (TMath::Abs(maPdgcode)==fParentSelect[0][i]){
1608 isFinalOpenCharm = kTRUE;
1611 if (!isFinalOpenCharm) return -1;
1613 // iterate until you find B hadron as a mother or become top ancester
1614 for (Int_t i=1; i<fgkMaxIter; i++){
1616 Int_t jLabel = GetMother(partMother);
1618 origin = kDirectCharm;
1621 if (jLabel < 0){ // safety protection
1622 AliDebug(1, "Stack label is negative, return\n");
1626 // if there is an ancester
1627 const AliVParticle* grandMa = fMCEvent->GetTrack(jLabel);
1628 Int_t grandMaPDG = grandMa->PdgCode();
1630 for (Int_t j=0; j<fNparents; j++){
1631 if (TMath::Abs(grandMaPDG)==fParentSelect[1][j]){
1632 origin = kBeautyCharm;
1637 partMother = grandMa;
1638 } // end of iteration
1640 else if ( int(TMath::Abs(maPdgcode)/100.) == kBeauty || int(TMath::Abs(maPdgcode)/1000.) == kBeauty ) {
1641 for (Int_t i=0; i<fNparents; i++){
1642 if (TMath::Abs(maPdgcode)==fParentSelect[1][i]){
1643 origin = kDirectBeauty;
1648 else if ( TMath::Abs(maPdgcode) == 22 ) {
1652 else if ( TMath::Abs(maPdgcode) == 111 ) {
1660 //__________________________________________
1661 Int_t AliHFEmcQA::GetSource(const TParticle * const mcpart) const
1663 // decay particle's origin
1665 //if ( TMath::Abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return -1;
1668 Bool_t isFinalOpenCharm = kFALSE;
1671 AliDebug(1, "no mcparticle, return\n");
1675 Int_t iLabel = mcpart->GetFirstMother();
1677 AliDebug(1, "Stack label is negative, return\n");
1681 AliMCParticle *mctrack = NULL;
1682 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return -1;
1683 TParticle *partMother = mctrack->Particle();
1684 Int_t maPdgcode = partMother->GetPdgCode();
1686 // if the mother is charmed hadron
1687 if ( int(TMath::Abs(maPdgcode)/100.) == kCharm || int(TMath::Abs(maPdgcode)/1000.) == kCharm ) {
1689 for (Int_t i=0; i<fNparents; i++){
1690 if (TMath::Abs(maPdgcode)==fParentSelect[0][i]){
1691 isFinalOpenCharm = kTRUE;
1694 if (!isFinalOpenCharm) return -1;
1696 // iterate until you find B hadron as a mother or become top ancester
1697 for (Int_t i=1; i<fgkMaxIter; i++){
1699 Int_t jLabel = partMother->GetFirstMother();
1701 origin = kDirectCharm;
1704 if (jLabel < 0){ // safety protection
1705 AliDebug(1, "Stack label is negative, return\n");
1709 // if there is an ancester
1710 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return -1;
1711 TParticle* grandMa = mctrack->Particle();
1712 Int_t grandMaPDG = grandMa->GetPdgCode();
1714 for (Int_t j=0; j<fNparents; j++){
1715 if (TMath::Abs(grandMaPDG)==fParentSelect[1][j]){
1716 origin = kBeautyCharm;
1721 partMother = grandMa;
1722 } // end of iteration
1724 else if ( int(TMath::Abs(maPdgcode)/100.) == kBeauty || int(TMath::Abs(maPdgcode)/1000.) == kBeauty ) {
1725 for (Int_t i=0; i<fNparents; i++){
1726 if (TMath::Abs(maPdgcode)==fParentSelect[1][i]){
1727 origin = kDirectBeauty;
1732 else if ( TMath::Abs(maPdgcode) == 22 ) {
1736 else if ( TMath::Abs(maPdgcode) == 111 ) {
1740 else origin = kElse;
1745 //__________________________________________
1746 Int_t AliHFEmcQA::GetElecSource(TParticle * const mcpart, Bool_t isElec) const
1748 // decay particle's origin
1751 AliDebug(1, "no mcparticle, return\n");
1755 if(isElec) if ( TMath::Abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return kMisID;
1758 Bool_t isFinalOpenCharm = kFALSE;
1760 Int_t iLabel = mcpart->GetFirstMother();
1762 AliDebug(1, "Stack label is negative, return\n");
1766 AliMCParticle *mctrack = NULL;
1767 Int_t tmpMomLabel=0;
1768 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return -1;
1769 TParticle *partMother = mctrack->Particle();
1770 TParticle *partMotherCopy = mctrack->Particle();
1771 Int_t maPdgcode = partMother->GetPdgCode();
1772 Int_t grmaPdgcode = 0;
1775 // if the mother is charmed hadron
1777 if(TMath::Abs(maPdgcode)==443){ // J/spi
1778 Int_t jLabel = partMother->GetFirstMother();
1779 if((mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))){
1780 TParticle* grandMa = mctrack->Particle();
1781 Int_t grandMaPDG = grandMa->GetPdgCode();
1782 if((int(TMath::Abs(grandMaPDG)/100.)%10) == kBeauty || (int(TMath::Abs(grandMaPDG)/1000.)%10) == kBeauty) return kB2Jpsi;
1786 else if ( (int(TMath::Abs(maPdgcode)/100.)%10) == kCharm || (int(TMath::Abs(maPdgcode)/1000.)%10) == kCharm ) {
1788 for (Int_t i=0; i<fNparents; i++){
1789 if (TMath::Abs(maPdgcode)==fParentSelect[0][i]){
1790 isFinalOpenCharm = kTRUE;
1793 if (!isFinalOpenCharm) return -1;
1795 // iterate until you find B hadron as a mother or become top ancester
1796 for (Int_t i=1; i<fgkMaxIter; i++){
1798 Int_t jLabel = partMother->GetFirstMother();
1800 return kDirectCharm;
1802 if (jLabel < 0){ // safety protection
1803 AliDebug(1, "Stack label is negative, return\n");
1807 // if there is an ancester
1808 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return -1;
1809 TParticle* grandMa = mctrack->Particle();
1810 Int_t grandMaPDG = grandMa->GetPdgCode();
1812 for (Int_t j=0; j<fNparents; j++){
1813 if (TMath::Abs(grandMaPDG)==fParentSelect[1][j]){
1814 return kBeautyCharm;
1818 partMother = grandMa;
1819 } // end of iteration
1821 else if ( (int(TMath::Abs(maPdgcode)/100.)%10) == kBeauty || (int(TMath::Abs(maPdgcode)/1000.)%10) == kBeauty ) {
1822 for (Int_t i=0; i<fNparents; i++){
1823 if (TMath::Abs(maPdgcode)==fParentSelect[1][i]){
1824 return kDirectBeauty;
1828 else if ( TMath::Abs(maPdgcode) == 22 ) { //conversion
1830 tmpMomLabel = partMotherCopy->GetFirstMother();
1831 if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tmpMomLabel))))) return -1;
1832 partMother = mctrack->Particle();
1833 maPdgcode = partMother->GetPdgCode();
1835 // check if the ligth meson is the decay product of heavy mesons
1836 tmpMomLabel = partMother->GetFirstMother();
1837 if((mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tmpMomLabel))))) {
1838 partMother = mctrack->Particle();
1839 grmaPdgcode = partMother->GetPdgCode();
1841 if ( (int(TMath::Abs(grmaPdgcode)/100.)%10) == kBeauty || (int(TMath::Abs(grmaPdgcode)/1000.)%10) == kBeauty ) return kGammaB2M;
1842 if ( (int(TMath::Abs(grmaPdgcode)/100.)%10) == kCharm || (int(TMath::Abs(grmaPdgcode)/1000.)%10) == kCharm ) return kGammaD2M;
1844 tmpMomLabel = partMother->GetFirstMother();
1845 if((mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tmpMomLabel))))) {
1846 partMother = mctrack->Particle();
1847 ggrmaPdgcode = partMother->GetPdgCode();
1849 if ( (int(TMath::Abs(ggrmaPdgcode)/100.)%10) == kBeauty || (int(TMath::Abs(ggrmaPdgcode)/1000.)%10) == kBeauty ) return kGammaB2M;
1850 if ( (int(TMath::Abs(ggrmaPdgcode)/100.)%10) == kCharm || (int(TMath::Abs(ggrmaPdgcode)/1000.)%10) == kCharm ) return kGammaD2M;
1854 if ( TMath::Abs(maPdgcode) == 111 ) {
1855 if(grmaPdgcode == 221 || grmaPdgcode == 223 || grmaPdgcode == 333 || grmaPdgcode == 331 || grmaPdgcode == 113) return kGammaM2M;
1858 else if ( TMath::Abs(maPdgcode) == 221 ) {
1859 if(grmaPdgcode == 111 || grmaPdgcode == 223 || grmaPdgcode == 333 || grmaPdgcode == 331 || grmaPdgcode == 113) return kGammaM2M;
1862 else if ( TMath::Abs(maPdgcode) == 223 ) {
1863 if(grmaPdgcode == 111 || grmaPdgcode == 221 || grmaPdgcode == 333 || grmaPdgcode == 331 || grmaPdgcode == 113) return kGammaM2M;
1866 else if ( TMath::Abs(maPdgcode) == 333 ) {
1867 if(grmaPdgcode == 111 || grmaPdgcode == 221 || grmaPdgcode == 223 || grmaPdgcode == 331 || grmaPdgcode == 113) return kGammaM2M;
1870 else if ( TMath::Abs(maPdgcode) == 331 ) {
1871 if(grmaPdgcode == 111 || grmaPdgcode == 221 || grmaPdgcode == 223 || grmaPdgcode == 333 || grmaPdgcode == 113) return kGammaM2M;
1872 return kGammaEtaPrime;
1874 else if ( TMath::Abs(maPdgcode) == 113 ) {
1875 if(grmaPdgcode == 111 || grmaPdgcode == 221 || grmaPdgcode == 223 || grmaPdgcode == 333 || grmaPdgcode == 331) return kGammaM2M;
1878 else origin = kElse;
1884 // check if the ligth meson is the decay product of heavy mesons
1885 tmpMomLabel = partMotherCopy->GetFirstMother();
1886 if((mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tmpMomLabel))))) {
1887 partMother = mctrack->Particle();
1888 grmaPdgcode = partMother->GetPdgCode();
1890 if ( (int(TMath::Abs(grmaPdgcode)/100.)%10) == kBeauty || (int(TMath::Abs(grmaPdgcode)/1000.)%10) == kBeauty ) return kB2M;
1891 if ( (int(TMath::Abs(grmaPdgcode)/100.)%10) == kCharm || (int(TMath::Abs(grmaPdgcode)/1000.)%10) == kCharm ) return kD2M;
1893 tmpMomLabel = partMother->GetFirstMother();
1894 if((mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tmpMomLabel))))) {
1895 partMother = mctrack->Particle();
1896 ggrmaPdgcode = partMother->GetPdgCode();
1898 if ( (int(TMath::Abs(ggrmaPdgcode)/100.)%10) == kBeauty || (int(TMath::Abs(ggrmaPdgcode)/1000.)%10) == kBeauty ) return kB2M;
1899 if ( (int(TMath::Abs(ggrmaPdgcode)/100.)%10) == kCharm || (int(TMath::Abs(ggrmaPdgcode)/1000.)%10) == kCharm ) return kD2M;
1903 if ( TMath::Abs(maPdgcode) == 111 ) {
1904 if(grmaPdgcode == 221 || grmaPdgcode == 223 || grmaPdgcode == 333 || grmaPdgcode == 331 || grmaPdgcode == 113) return kM2M;
1907 else if ( TMath::Abs(maPdgcode) == 221 ) {
1908 if(grmaPdgcode == 111 || grmaPdgcode == 223 || grmaPdgcode == 333 || grmaPdgcode == 331 || grmaPdgcode == 113) return kM2M;
1911 else if ( TMath::Abs(maPdgcode) == 223 ) {
1912 if(grmaPdgcode == 111 || grmaPdgcode == 221 || grmaPdgcode == 333 || grmaPdgcode == 331 || grmaPdgcode == 113) return kM2M;
1915 else if ( TMath::Abs(maPdgcode) == 333 ) {
1916 if(grmaPdgcode == 111 || grmaPdgcode == 221 || grmaPdgcode == 223 || grmaPdgcode == 331 || grmaPdgcode == 113) return kM2M;
1919 else if ( TMath::Abs(maPdgcode) == 331 ) {
1920 if(grmaPdgcode == 111 || grmaPdgcode == 221 || grmaPdgcode == 223 || grmaPdgcode == 333 || grmaPdgcode == 113) return kM2M;
1923 else if ( TMath::Abs(maPdgcode) == 113 ) {
1924 if(grmaPdgcode == 111 || grmaPdgcode == 221 || grmaPdgcode == 223 || grmaPdgcode == 333 || grmaPdgcode == 331) return kM2M;
1927 else if ( TMath::Abs(maPdgcode) == 321 || TMath::Abs(maPdgcode) == 130 ) {
1938 //__________________________________________
1939 Int_t AliHFEmcQA::GetElecSource(const AliVParticle * const mctrack, Bool_t isElec) const
1942 // decay particle's origin
1946 AliDebug(1, "no mcparticle, return\n");
1950 TClass *type = mctrack->IsA();
1952 if(type == AliMCParticle::Class()) {
1953 const AliMCParticle *esdmc = dynamic_cast<const AliMCParticle *>(mctrack);
1955 TParticle *mcpart = esdmc->Particle();
1956 return GetElecSource(mcpart, isElec);
1960 if(type == AliAODMCParticle::Class()) {
1961 const AliAODMCParticle *aodmc = dynamic_cast<const AliAODMCParticle *>(mctrack);
1963 return GetElecSource(aodmc, isElec);
1969 //__________________________________________
1970 Int_t AliHFEmcQA::GetElecSource(const AliAODMCParticle * const mcpart, Bool_t isElec) const
1973 // Function for AliAODMCParticle
1976 if (!mcpart) return -1;
1977 if (!fMCArray) return -1;
1979 if(isElec) if ( TMath::Abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return kMisID;
1982 Bool_t isFinalOpenCharm = kFALSE;
1984 Int_t iLabel = mcpart->GetMother();
1985 if ((iLabel<0) || (iLabel>=fMCArray->GetEntriesFast())){
1986 AliDebug(1, "label is out of range, return\n");
1990 AliAODMCParticle *mctrack = NULL; // will change all the time
1991 Int_t tmpMomLabel=0;
1992 if(!(mctrack = dynamic_cast<AliAODMCParticle *>(fMCArray->At(TMath::Abs(iLabel))))) return -1;
1993 AliAODMCParticle *partMother = mctrack;
1994 AliAODMCParticle *partMotherCopy = mctrack;
1995 Int_t maPdgcode = mctrack->GetPdgCode();
1999 // if the mother is charmed hadron
2001 if(TMath::Abs(maPdgcode)==443){
2005 Int_t jLabel = partMother->GetMother();
2006 if ((jLabel>=0) && (jLabel<fMCArray->GetEntriesFast())){
2007 if((mctrack = dynamic_cast<AliAODMCParticle *>(fMCArray->At(TMath::Abs(jLabel))))){
2008 Int_t grandMaPDG = mctrack->GetPdgCode();
2009 if((int(TMath::Abs(grandMaPDG)/100.)%10) == kBeauty || (int(TMath::Abs(grandMaPDG)/1000.)%10) == kBeauty) {
2016 else if ( (int(TMath::Abs(maPdgcode)/100.)%10) == kCharm || (int(TMath::Abs(maPdgcode)/1000.)%10) == kCharm ) {
2020 for (Int_t i=0; i<fNparents; i++){
2021 if (TMath::Abs(maPdgcode)==fParentSelect[0][i]){
2022 isFinalOpenCharm = kTRUE;
2025 if (!isFinalOpenCharm) {
2029 // iterate until you find B hadron as a mother or become top ancester
2030 for (Int_t i=1; i<fgkMaxIter; i++){
2032 Int_t jLabel = partMother->GetMother();
2034 return kDirectCharm;
2036 if ((jLabel<0) || (jLabel>=fMCArray->GetEntriesFast())){
2037 AliDebug(1, "Stack label is negative, return\n");
2041 // if there is an ancester
2042 if(!(mctrack = dynamic_cast<AliAODMCParticle *>(fMCArray->At(TMath::Abs(jLabel))))) {
2045 Int_t grandMaPDG = mctrack->GetPdgCode();
2046 for (Int_t j=0; j<fNparents; j++){
2047 if (TMath::Abs(grandMaPDG)==fParentSelect[1][j]){
2048 return kBeautyCharm;
2051 partMother = mctrack;
2052 } // end of iteration
2055 else if ( (int(TMath::Abs(maPdgcode)/100.)%10) == kBeauty || (int(TMath::Abs(maPdgcode)/1000.)%10) == kBeauty ) {
2059 for (Int_t i=0; i<fNparents; i++){
2060 if (TMath::Abs(maPdgcode)==fParentSelect[1][i]){
2061 return kDirectBeauty;
2065 else if ( TMath::Abs(maPdgcode) == 22 ) {
2069 tmpMomLabel = partMotherCopy->GetMother();
2070 if((tmpMomLabel<0) || (tmpMomLabel>=fMCArray->GetEntriesFast())) {
2073 if(!(mctrack = dynamic_cast<AliAODMCParticle *>(fMCArray->At(TMath::Abs(tmpMomLabel))))) {
2076 partMother = mctrack;
2077 maPdgcode = partMother->GetPdgCode();
2079 // check if the ligth meson is the decay product of heavy mesons
2080 tmpMomLabel = partMother->GetMother();
2081 if((tmpMomLabel>=0) && (tmpMomLabel<fMCArray->GetEntriesFast())) {
2082 if((mctrack = dynamic_cast<AliAODMCParticle *>(fMCArray->At(TMath::Abs(tmpMomLabel))))) {
2083 partMother = mctrack;
2084 grmaPdgcode = partMother->GetPdgCode();
2086 if ( (int(TMath::Abs(grmaPdgcode)/100.)%10) == kBeauty || (int(TMath::Abs(grmaPdgcode)/1000.)%10) == kBeauty ) {
2089 if ( (int(TMath::Abs(grmaPdgcode)/100.)%10) == kCharm || (int(TMath::Abs(grmaPdgcode)/1000.)%10) == kCharm ) {
2093 tmpMomLabel = partMother->GetMother();
2094 if((tmpMomLabel>=0) && (tmpMomLabel<fMCArray->GetEntriesFast())) {
2095 if((mctrack = dynamic_cast<AliAODMCParticle *>(fMCArray->At(TMath::Abs(tmpMomLabel))))) {
2096 partMother = mctrack;
2097 ggrmaPdgcode = partMother->GetPdgCode();
2099 if ( (int(TMath::Abs(ggrmaPdgcode)/100.)%10) == kBeauty || (int(TMath::Abs(ggrmaPdgcode)/1000.)%10) == kBeauty ) {
2102 if ( (int(TMath::Abs(ggrmaPdgcode)/100.)%10) == kCharm || (int(TMath::Abs(ggrmaPdgcode)/1000.)%10) == kCharm ) {
2110 if ( TMath::Abs(maPdgcode) == 111 ) {
2111 if(grmaPdgcode == 221 || grmaPdgcode == 223 || grmaPdgcode == 333 || grmaPdgcode == 331 || grmaPdgcode == 113) return kGammaM2M;
2114 else if ( TMath::Abs(maPdgcode) == 221 ) {
2115 if(grmaPdgcode == 111 || grmaPdgcode == 223 || grmaPdgcode == 333 || grmaPdgcode == 331 || grmaPdgcode == 113) return kGammaM2M;
2118 else if ( TMath::Abs(maPdgcode) == 223 ) {
2119 if(grmaPdgcode == 111 || grmaPdgcode == 221 || grmaPdgcode == 333 || grmaPdgcode == 331 || grmaPdgcode == 113) return kGammaM2M;
2122 else if ( TMath::Abs(maPdgcode) == 333 ) {
2123 if(grmaPdgcode == 111 || grmaPdgcode == 221 || grmaPdgcode == 223 || grmaPdgcode == 331 || grmaPdgcode == 113) return kGammaM2M;
2126 else if ( TMath::Abs(maPdgcode) == 331 ) {
2127 if(grmaPdgcode == 111 || grmaPdgcode == 221 || grmaPdgcode == 223 || grmaPdgcode == 333 || grmaPdgcode == 113) return kGammaM2M;
2128 return kGammaEtaPrime;
2130 else if ( TMath::Abs(maPdgcode) == 113 ) {
2131 if(grmaPdgcode == 111 || grmaPdgcode == 221 || grmaPdgcode == 223 || grmaPdgcode == 333 || grmaPdgcode == 331) return kGammaM2M;
2134 else origin = kElse;
2141 // check if the ligth meson is the decay product of heavy mesons
2143 tmpMomLabel = partMotherCopy->GetMother();
2144 if((tmpMomLabel>=0) && (tmpMomLabel<fMCArray->GetEntriesFast())) {
2145 if((mctrack = dynamic_cast<AliAODMCParticle *>(fMCArray->At(TMath::Abs(tmpMomLabel))))) {
2146 partMother = mctrack;
2147 grmaPdgcode = partMother->GetPdgCode();
2149 if ( (int(TMath::Abs(grmaPdgcode)/100.)%10) == kBeauty || (int(TMath::Abs(grmaPdgcode)/1000.)%10) == kBeauty ) {
2152 if ( (int(TMath::Abs(grmaPdgcode)/100.)%10) == kCharm || (int(TMath::Abs(grmaPdgcode)/1000.)%10) == kCharm ) {
2156 tmpMomLabel = partMother->GetMother();
2157 if((tmpMomLabel>=0) && (tmpMomLabel<fMCArray->GetEntriesFast())) {
2158 if((mctrack = dynamic_cast<AliAODMCParticle *>(fMCArray->At(TMath::Abs(tmpMomLabel))))) {
2159 partMother = mctrack;
2160 ggrmaPdgcode = partMother->GetPdgCode();
2162 if ( (int(TMath::Abs(ggrmaPdgcode)/100.)%10) == kBeauty || (int(TMath::Abs(ggrmaPdgcode)/1000.)%10) == kBeauty ) {
2165 if ( (int(TMath::Abs(ggrmaPdgcode)/100.)%10) == kCharm || (int(TMath::Abs(ggrmaPdgcode)/1000.)%10) == kCharm ) {
2173 if ( TMath::Abs(maPdgcode) == 111 ) {
2174 if(grmaPdgcode == 221 || grmaPdgcode == 223 || grmaPdgcode == 333 || grmaPdgcode == 331 || grmaPdgcode == 113) return kM2M;
2177 else if ( TMath::Abs(maPdgcode) == 221 ) {
2178 if(grmaPdgcode == 111 || grmaPdgcode == 223 || grmaPdgcode == 333 || grmaPdgcode == 331 || grmaPdgcode == 113) return kM2M;
2181 else if ( TMath::Abs(maPdgcode) == 223 ) {
2182 if(grmaPdgcode == 111 || grmaPdgcode == 221 || grmaPdgcode == 333 || grmaPdgcode == 331 || grmaPdgcode == 113) return kM2M;
2185 else if ( TMath::Abs(maPdgcode) == 333 ) {
2186 if(grmaPdgcode == 111 || grmaPdgcode == 221 || grmaPdgcode == 223 || grmaPdgcode == 331 || grmaPdgcode == 113) return kM2M;
2189 else if ( TMath::Abs(maPdgcode) == 331 ) {
2190 if(grmaPdgcode == 111 || grmaPdgcode == 221 || grmaPdgcode == 223 || grmaPdgcode == 333 || grmaPdgcode == 113) return kM2M;
2193 else if ( TMath::Abs(maPdgcode) == 113 ) {
2194 if(grmaPdgcode == 111 || grmaPdgcode == 221 || grmaPdgcode == 223 || grmaPdgcode == 333 || grmaPdgcode == 331) return kM2M;
2197 else if ( TMath::Abs(maPdgcode) == 321 || TMath::Abs(maPdgcode) == 130 ) {
2207 //__________________________________________
2208 Double_t AliHFEmcQA::GetWeightFactor(AliMCParticle *mctrack, const Int_t iBgLevel){
2210 // Get weighting factor for the realistic background estimation, for three possible background yield levels, indicated by the argument "iLevel": the best estimate (0), the lower uncertainty level (1), and the upper uncertainty level (2)
2212 AliMCParticle *mctrackmother = NULL;
2213 Double_t weightElecBg = 0.;
2214 Double_t mesonPt = 0.;
2215 Double_t bgcategory = 0.;
2217 Int_t mesonID = GetElecSource(mctrack->Particle(), kTRUE);
2218 if(mesonID==kGammaPi0 || mesonID==kPi0) mArr=0; //pion
2219 else if(mesonID==kGammaEta || mesonID==kEta) mArr=1; //eta
2220 else if(mesonID==kGammaOmega || mesonID==kOmega) mArr=2; //omega
2221 else if(mesonID==kGammaPhi || mesonID==kPhi) mArr=3; //phi
2222 else if(mesonID==kGammaEtaPrime || mesonID==kEtaPrime) mArr=4; //etaprime
2223 else if(mesonID==kGammaRho0 || mesonID==kRho0) mArr=5; //rho
2225 Double_t datamc[25]={-999,-999,-999,-999,-999,-999,-999,-999,-999,-999,-999,-999,-999,-999,-999,-999, -999, -999, -999, -999, -999, -999, -999, -999, -999};
2226 Double_t xr[3]={-999,-999,-999};
2227 datamc[0] = mesonID;
2228 datamc[17] = mctrack->Pt(); //electron pt
2229 datamc[18] = mctrack->Eta(); //electron eta
2231 mctrack->XvYvZv(xr);
2232 datamc[9] = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
2235 TParticle *mcpart = mctrack->Particle();
2237 datamc[14] = mcpart->GetUniqueID();
2241 if(mesonID>=kGammaPi0) { // conversion electron, be careful with the enum odering
2242 Int_t glabel=TMath::Abs(mctrack->GetMother()); // gamma label
2243 if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
2244 glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's label
2245 if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
2246 mesonPt = mctrackmother->Pt(); //meson pt
2248 datamc[1] = bgcategory;
2249 datamc[2] = mesonPt;
2250 mctrackmother->XvYvZv(xr);
2251 datamc[11] = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
2254 mcpart = mctrackmother->Particle();
2256 datamc[15] = mcpart->GetUniqueID();
2258 if(glabel>fMCEvent->GetNumberOfPrimaries()) {
2260 datamc[1] = bgcategory;
2261 //printf("I should be gamma meson = %d mesonlabel= %d NumberOfPrimaries= %d \n",mctrackmother->PdgCode(),glabel,fMCEvent->GetNumberOfPrimaries());
2262 glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother
2263 if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
2264 datamc[3]=mctrackmother->PdgCode();
2265 datamc[4]=mctrackmother->Pt();
2266 if(TMath::Abs(mctrackmother->PdgCode())==310){
2268 glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother's mother
2269 if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
2270 datamc[5]=mctrackmother->PdgCode();
2271 glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother
2272 if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
2273 datamc[6]=mctrackmother->PdgCode();
2282 else{ // nonHFE except for the conversion electron
2283 Int_t glabel=TMath::Abs(mctrack->GetMother());
2284 if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
2285 mesonPt=mctrackmother->Pt(); //meson pt
2286 if(mesonID==kEta) bgcategory = -1.41; // to consider new branching ratio for the eta Dalitz decay
2287 else bgcategory = -1.;
2288 datamc[1] = bgcategory;
2289 datamc[2] = mesonPt;
2290 mctrackmother->XvYvZv(xr);
2291 datamc[11] = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
2294 mcpart = mctrackmother->Particle();
2296 datamc[15] = mcpart->GetUniqueID();
2298 if(glabel>fMCEvent->GetNumberOfPrimaries()) {
2299 if(mesonID==kEta) bgcategory = -2.82; // to consider new branching ratio for the eta Dalitz decay (1.41)
2300 else bgcategory = -2.;
2301 datamc[1] = bgcategory;
2302 glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother
2303 if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
2304 datamc[3]=mctrackmother->PdgCode();
2305 datamc[4]=mctrackmother->Pt();
2306 if(TMath::Abs(mctrackmother->PdgCode())==310){
2308 glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother's mother
2309 if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
2310 datamc[5]=mctrackmother->PdgCode();
2311 glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother
2312 if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
2313 datamc[6]=mctrackmother->PdgCode();
2324 Int_t centBin = GetWeightCentralityBin(fPerCentrality);
2325 if(centBin < 0)return 0.;
2326 weightElecBg=fElecBackgroundFactor[iBgLevel][centBin][mArr][kBgPtBins-1];
2327 for(int ii=0; ii<kBgPtBins; ii++){
2328 if((mesonPt > fBinLimit[ii]) && (mesonPt < fBinLimit[ii+1])){
2329 weightElecBg = fElecBackgroundFactor[iBgLevel][centBin][mArr][ii];
2335 weightElecBg=fElecBackgroundFactor[iBgLevel][0][mArr][kBgPtBins-1];
2336 for(int ii=0; ii<kBgPtBins; ii++){
2337 if((mesonPt > fBinLimit[ii]) && (mesonPt < fBinLimit[ii+1])){
2338 weightElecBg = fElecBackgroundFactor[iBgLevel][0][mArr][ii];
2345 datamc[13] = weightElecBg;
2346 datamc[16] = Double_t(fContainerStep);
2348 datamc[7] = fHfeImpactR;
2349 datamc[8] = fHfeImpactnsigmaR;
2351 datamc[19] = fRecPt;
2352 datamc[20] = fRecEta;
2353 datamc[21] = fRecPhi;
2354 datamc[22] = fLyrhit;
2355 datamc[23] = fLyrstat;
2356 datamc[24] = fCentrality;
2358 if(fIsDebugStreamerON && fTreeStream){
2360 (*fTreeStream)<<"nonhfeQA"<<
2361 "mesonID="<<datamc[0]<<
2362 "bgcategory="<<datamc[1]<<
2363 "mesonPt="<<datamc[2]<<
2364 "mesonMomPdg="<<datamc[3]<<
2365 "mesonMomPt="<<datamc[4]<<
2366 "mesonGMomPdg="<<datamc[5]<<
2367 "mesonGGMomPdg="<<datamc[6]<<
2368 "eIPAbs="<<datamc[7]<<
2369 "eIPSig="<<datamc[8]<<
2372 "mesonR="<<datamc[11]<<
2373 "mesonZ="<<datamc[12]<<
2374 "weightElecBg="<<datamc[13]<<
2375 "eUniqID="<<datamc[14]<<
2376 "mesonUniqID="<<datamc[15]<<
2377 "containerStep="<<datamc[16]<<
2378 "emcpt="<<datamc[17]<<
2379 "emceta="<<datamc[18]<<
2380 "erecpt="<<datamc[19]<<
2381 "ereceta="<<datamc[20]<<
2382 "erecphi="<<datamc[21]<<
2383 "itshit="<<datamc[22]<<
2384 "itsstat="<<datamc[23]<<
2385 "centrality="<<datamc[24]
2390 Double_t returnval = bgcategory*weightElecBg;
2395 //__________________________________________
2396 Double_t AliHFEmcQA::GetWeightFactor(const AliAODMCParticle * const mcpart, const Int_t iBgLevel){
2398 // Get weighting factor for the realistic background estimation, for three possible background yield levels, indicated by the argument "iLevel": the best estimate (0), the lower uncertainty level (1), and the upper uncertainty level (2)
2401 if (!mcpart) return 0;
2402 if (!fMCArray) return 0;
2404 Double_t weightElecBg = 0.;
2405 Double_t mesonPt = 0.;
2406 Double_t bgcategory = 0.;
2408 Int_t mesonID = GetElecSource(mcpart, kTRUE);
2409 if(mesonID==kGammaPi0 || mesonID==kPi0) mArr=0; //pion
2410 else if(mesonID==kGammaEta || mesonID==kEta) mArr=1; //eta
2411 else if(mesonID==kGammaOmega || mesonID==kOmega) mArr=2; //omega
2412 else if(mesonID==kGammaPhi || mesonID==kPhi) mArr=3; //phi
2413 else if(mesonID==kGammaEtaPrime || mesonID==kEtaPrime) mArr=4; //etaprime
2414 else if(mesonID==kGammaRho0 || mesonID==kRho0) mArr=5; //rho
2418 AliAODMCParticle *mctrackmother = NULL; // will change all the time
2420 if(mesonID>=kGammaPi0) { // conversion electron, be careful with the enum odering
2421 Int_t iLabel = mcpart->GetMother(); //gamma label
2422 if(!(mctrackmother = dynamic_cast<AliAODMCParticle *>(fMCArray->At(TMath::Abs(iLabel))))) return 0;
2423 iLabel = mctrackmother->GetMother(); //gamma's mother's label
2424 if(!(mctrackmother= dynamic_cast<AliAODMCParticle *>(fMCArray->At(TMath::Abs(iLabel))))) return 0;
2425 mesonPt = mctrackmother->Pt(); //meson pt
2428 //if(glabel>fMCEvent->GetNumberOfPrimaries()) bgcategory = 2.;
2430 else{ // nonHFE except for the conversion electron
2431 Int_t iLabel = mcpart->GetMother(); //meson label
2432 if(!(mctrackmother = dynamic_cast<AliAODMCParticle *>(fMCArray->At(TMath::Abs(iLabel))))) return 0;
2433 mesonPt = mctrackmother->Pt(); //meson pt
2434 if(mesonID==kEta) bgcategory = -1.41; // to consider new branching ratio for the eta Dalitz decay
2435 else bgcategory = -1.;
2437 //if(glabel>fMCEvent->GetNumberOfPrimaries()) bgcategory = -2.;
2441 Int_t centBin = GetWeightCentralityBin(fPerCentrality);
2442 if(centBin < 0)return 0.;
2443 weightElecBg=fElecBackgroundFactor[iBgLevel][centBin][mArr][kBgPtBins-1];
2444 for(int ii=0; ii<kBgPtBins; ii++){
2445 if((mesonPt > fBinLimit[ii]) && (mesonPt < fBinLimit[ii+1])){
2446 weightElecBg = fElecBackgroundFactor[iBgLevel][centBin][mArr][ii];
2452 weightElecBg=fElecBackgroundFactor[iBgLevel][0][mArr][kBgPtBins-1];
2453 for(int ii=0; ii<kBgPtBins; ii++){
2454 if((mesonPt > fBinLimit[ii]) && (mesonPt < fBinLimit[ii+1])){
2455 weightElecBg = fElecBackgroundFactor[iBgLevel][0][mArr][ii];
2462 Double_t returnval = bgcategory*weightElecBg;
2467 //__________________________________________
2468 Double_t AliHFEmcQA::GetWeightFactorForPrimaries(const AliAODMCParticle * const mcpart, const Int_t iBgLevel){
2470 // Get weighting factor for the realistic background estimation, for three possible background yield levels, indicated by the argument "iLevel": the best estimate (0), the lower uncertainty level (1), and the upper uncertainty level (2)
2471 // weighting will only be non zero for electrons from primary pi0 and eta mesons (via Dalitz or gamma conversions)
2474 if (!mcpart) return 0;
2475 if (!fMCArray) return 0;
2478 Double_t mesonPt = 0.;
2479 AliAODMCParticle *mctrackmother = NULL; // temp pointer
2480 Int_t mesonID = GetElecSource(mcpart, kTRUE); // get source of electron
2481 if(mesonID==kGammaPi0 || mesonID==kPi0) mArr=0; //pion
2482 else if(mesonID==kGammaEta || mesonID==kEta) mArr=1; //eta
2485 Int_t iLabel = mcpart->GetMother(); //mother label
2489 if(!(mctrackmother = dynamic_cast<AliAODMCParticle *>(fMCArray->At(TMath::Abs(iLabel))))) return 0;
2490 iLabel = mctrackmother->GetMother(); //reset label to mother of gamma
2493 if(!(mctrackmother = dynamic_cast<AliAODMCParticle *>(fMCArray->At(TMath::Abs(iLabel))))) return 0;
2494 mesonPt = mctrackmother->Pt(); //pt of (grand)mother
2495 //check mother is primary
2496 if(!(mctrackmother->IsPrimary())) return 0;
2501 Double_t weightElecBg = fElecBackgroundFactor[iBgLevel][0][mArr][kBgPtBins-1]; //set weighting for pt > max pt
2502 for(int ii=0; ii<kBgPtBins; ii++){
2503 if((mesonPt >= fBinLimit[ii]) && (mesonPt < fBinLimit[ii+1])){
2504 weightElecBg = fElecBackgroundFactor[iBgLevel][0][mArr][ii];
2509 return weightElecBg;
2512 //__________________________________________
2513 //__________________________________________
2514 Int_t AliHFEmcQA::GetMother(const AliVParticle * const mcpart) const {
2516 // Wrapper to get the mother label
2519 TClass *type = mcpart->IsA();
2520 if(type == AliMCParticle::Class()){
2522 const AliMCParticle *emcpart = static_cast<const AliMCParticle *>(mcpart);
2523 label = emcpart->GetMother();
2524 } else if(type == AliAODMCParticle::Class()){
2526 const AliAODMCParticle *amcpart = static_cast<const AliAODMCParticle *>(mcpart);
2527 label = amcpart->GetMother();
2531 //__________________________________________
2532 Int_t AliHFEmcQA::GetWeightCentralityBin(const Float_t percentile) const {
2534 //translate the centrality percentile into the centrality bin of the reference weighting histograms for electron background
2537 Float_t centralityLimits[12]= {0.,5.,10., 20., 30., 40., 50., 60.,70.,80., 90., 100.};
2539 for(Int_t ibin = 0; ibin < 11; ibin++){
2540 if(percentile >= centralityLimits[ibin] && percentile < centralityLimits[ibin+1]){
2547 //__________________________________________
2548 void AliHFEmcQA::AliHists::FillList(TList *l) const {
2550 // Fill Histos into a list for output
2552 if(fPdgCode) l->Add(fPdgCode);
2553 if(fPt) l->Add(fPt);
2555 if(fEta) l->Add(fEta);
2558 //__________________________________________
2559 void AliHFEmcQA::AliHistsComm::FillList(TList *l) const {
2561 // Fill Histos into a list for output
2563 if(fNq) l->Add(fNq);
2564 if(fProcessID) l->Add(fProcessID);
2565 if(fePtRatio) l->Add(fePtRatio);
2566 if(fPtCorr) l->Add(fPtCorr);
2567 if(fPtCorrDp) l->Add(fPtCorrDp);
2568 if(fPtCorrD0) l->Add(fPtCorrD0);
2569 if(fPtCorrDrest) l->Add(fPtCorrDrest);
2571 if(fPtCorrDinein) l->Add(fPtCorrDinein);
2572 if(fPtCorrDineout) l->Add(fPtCorrDineout);
2573 if(fPtCorrDoutein) l->Add(fPtCorrDoutein);
2574 if(fPtCorrDouteout) l->Add(fPtCorrDouteout);
2575 if(fPtCorrDpDinein) l->Add(fPtCorrDpDinein);
2576 if(fPtCorrDpDineout) l->Add(fPtCorrDpDineout);
2577 if(fPtCorrDpDoutein) l->Add(fPtCorrDpDoutein);
2578 if(fPtCorrDpDouteout) l->Add(fPtCorrDpDouteout);
2579 if(fPtCorrD0Dinein) l->Add(fPtCorrD0Dinein);
2580 if(fPtCorrD0Dineout) l->Add(fPtCorrD0Dineout);
2581 if(fPtCorrD0Doutein) l->Add(fPtCorrD0Doutein);
2582 if(fPtCorrD0Douteout) l->Add(fPtCorrD0Douteout);
2583 if(fPtCorrDrestDinein) l->Add(fPtCorrDrestDinein);
2584 if(fPtCorrDrestDineout) l->Add(fPtCorrDrestDineout);
2585 if(fPtCorrDrestDoutein) l->Add(fPtCorrDrestDoutein);
2586 if(fPtCorrDrestDouteout) l->Add(fPtCorrDrestDouteout);
2588 if(fEtaCorrD) l->Add(fEtaCorrD);
2589 if(fEtaCorrDp) l->Add(fEtaCorrDp);
2590 if(fEtaCorrD0) l->Add(fEtaCorrD0);
2591 if(fEtaCorrDrest) l->Add(fEtaCorrDrest);
2593 if(fEtaCorrGD) l->Add(fEtaCorrGD);
2594 if(fEtaCorrGDp) l->Add(fEtaCorrGDp);
2595 if(fEtaCorrGD0) l->Add(fEtaCorrGD0);
2596 if(fEtaCorrGDrest) l->Add(fEtaCorrGDrest);
2598 if(fEtaCorrB) l->Add(fEtaCorrB);
2599 if(fEtaCorrGB) l->Add(fEtaCorrGB);
2600 if(fPtCorrBinein) l->Add(fPtCorrBinein);
2601 if(fPtCorrBineout) l->Add(fPtCorrBineout);
2602 if(fPtCorrBoutein) l->Add(fPtCorrBoutein);
2603 if(fPtCorrBouteout) l->Add(fPtCorrBouteout);
2605 if(fDePtRatio) l->Add(fDePtRatio);
2606 if(feDistance) l->Add(feDistance);
2607 if(fDeDistance) l->Add(fDeDistance);