fbea85a91f323d2efb740e69bbcd63954f618e78
[u/mrichter/AliRoot.git] / PWGHF / hfe / AliHFEmcQA.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15 //
16 // QA class of Heavy Flavor quark and fragmeted/decayed particles
17 // -Check kinematics of Heavy Quarks/hadrons, and decayed leptons
18 //    pT, rapidity
19 //    decay lepton kinematics w/wo acceptance
20 //    heavy hadron decay length, electron pT fraction carried from decay
21 // -Check yield of Heavy Quarks/hadrons
22 //    Number of produced heavy quark
23 //    Number of produced hadron of given pdg code
24 //
25 //
26 // Authors:
27 //   MinJung Kweon <minjung@physi.uni-heidelberg.de>
28 //
29
30 #include <TH2F.h>
31 #include <TH1F.h>
32 #include <TH2F.h>
33 #include <TList.h>
34 #include <TParticle.h>
35 #include "TTreeStream.h"
36
37 #include <AliLog.h>
38 #include <AliMCEvent.h>
39 #include <AliGenEventHeader.h>
40 #include <AliAODMCParticle.h>
41 #include <AliStack.h>
42
43 #include "AliHFEmcQA.h"
44 #include "AliHFEtools.h"
45 #include "AliHFEcollection.h"
46
47 ClassImp(AliHFEmcQA)
48
49 //_______________________________________________________________________________________________
50     AliHFEmcQA::AliHFEmcQA() :
51          fMCEvent(NULL)
52         ,fMCHeader(NULL)
53         ,fMCArray(NULL)
54         ,fQAhistos(NULL)
55         ,fMCQACollection(NULL)
56         ,fNparents(0)
57         ,fCentrality(0)
58         ,fPerCentrality(-1)
59         ,fIsPbPb(kFALSE)
60         ,fIsppMultiBin(kFALSE)
61         ,fContainerStep(0)
62         ,fIsDebugStreamerON(kFALSE)
63         ,fRecPt(-999)
64         ,fRecEta(-999)
65         ,fRecPhi(-999)
66         ,fLyrhit(0)
67         ,fLyrstat(0)
68         ,fHfeImpactR(-999)
69         ,fHfeImpactnsigmaR(-999)
70         ,fTreeStream(NULL)
71 {
72         // Default constructor
73   for(Int_t mom = 0; mom < 9; mom++){
74     fhD[mom] = NULL;
75   }
76   for(Int_t mom = 0; mom < 50; mom++){
77     fHeavyQuark[mom] = NULL;
78   }
79   for(Int_t mom = 0; mom < 2; mom++){
80     fIsHeavy[mom] = 0;
81   }
82   memset(fElecBackgroundFactor, 0, sizeof(Double_t) * kElecBgSpecies * kBgPtBins * kCentBins * kBgLevels);
83   memset(fBinLimit, 0, sizeof(Double_t) * (kBgPtBins+1));
84 }
85
86 //_______________________________________________________________________________________________
87 AliHFEmcQA::AliHFEmcQA(const AliHFEmcQA&p):
88 TObject(p)
89         ,fMCEvent(NULL)
90         ,fMCHeader(NULL)
91         ,fMCArray(NULL)
92         ,fQAhistos(p.fQAhistos)
93         ,fMCQACollection(p.fMCQACollection)
94         ,fNparents(p.fNparents)
95         ,fCentrality(0)
96         ,fPerCentrality(-1)
97         ,fIsPbPb(kFALSE)
98         ,fIsppMultiBin(kFALSE)
99         ,fContainerStep(0)
100         ,fIsDebugStreamerON(kFALSE)
101         ,fRecPt(-999)
102         ,fRecEta(-999)
103         ,fRecPhi(-999)
104         ,fLyrhit(0)
105         ,fLyrstat(0)
106         ,fHfeImpactR(0)
107         ,fHfeImpactnsigmaR(0)
108         ,fTreeStream(NULL)
109 {
110         // Copy constructor
111   for(Int_t mom = 0; mom < 9; mom++){
112     fhD[mom] = NULL;
113   }
114   for(Int_t mom = 0; mom < 50; mom++){
115     fHeavyQuark[mom] = NULL;
116   }
117   for(Int_t mom = 0; mom < 2; mom++){
118     fIsHeavy[mom] = 0;
119   }
120   memset(fElecBackgroundFactor, 0, sizeof(Double_t) * kElecBgSpecies * kBgPtBins * kCentBins * kBgLevels);
121   memset(fBinLimit, 0, sizeof(Double_t) * (kBgPtBins+1));
122 }
123
124 //_______________________________________________________________________________________________
125 AliHFEmcQA&
126 AliHFEmcQA::operator=(const AliHFEmcQA &)
127 {
128   // Assignment operator
129
130   AliInfo("Not yet implemented.");
131   return *this;
132 }
133
134 //_______________________________________________________________________________________________
135 AliHFEmcQA::~AliHFEmcQA()
136 {
137         // Destructor
138
139         if(fTreeStream && fIsDebugStreamerON) delete fTreeStream;
140         AliInfo("Analysis Done.");
141 }
142
143 //_______________________________________________________________________________________________
144 void AliHFEmcQA::PostAnalyze() const
145 {
146         //
147         // Post analysis
148         //
149 }
150
151 //_______________________________________________________________________________________________
152 void AliHFEmcQA::SetBackgroundWeightFactor(Double_t *elecBackgroundFactor, Double_t *binLimit)
153 {
154    //
155    // copy background weighting factors into data member
156    //
157   
158   memcpy(fElecBackgroundFactor,elecBackgroundFactor,sizeof(Double_t) * kElecBgSpecies * kBgPtBins * kCentBins * kBgLevels);
159   memcpy(fBinLimit,binLimit,sizeof(Double_t) * (kBgPtBins+1));
160 }
161
162 //__________________________________________
163 void AliHFEmcQA::CreatDefaultHistograms(TList * const qaList)
164 {      
165   //
166   // make default histograms
167   //
168   
169   if(!qaList) return;
170
171   fQAhistos = qaList;
172   fQAhistos->SetName("MCqa");
173
174   CreateHistograms(AliHFEmcQA::kCharm);               // create histograms for charm
175   CreateHistograms(AliHFEmcQA::kBeauty);              // create histograms for beauty
176   CreateHistograms(AliHFEmcQA::kOthers);              // create histograms for beauty
177  
178 // prepare 2D(pt vs Y) histogram for D spectra, we consider following 9 particles
179   const Int_t nbspecies = 9;
180   TString kDspecies[nbspecies];
181   kDspecies[0]="411";   //D+
182   kDspecies[1]="421";   //D0
183   kDspecies[2]="431";   //Ds+
184   kDspecies[3]="4122";  //Lambdac+
185   kDspecies[4]="4132";  //Ksic0
186   kDspecies[5]="4232";  //Ksic+
187   kDspecies[6]="4332";  //OmegaC0
188   kDspecies[7]="413";   //D*(2010)+
189   kDspecies[8]="423";   //D*(2007)0
190
191   const Double_t kPtbound[2] = {0.1, 20.}; //bin taken for considering inclusive e analysis binning
192   Int_t iBin[2];
193   iBin[0] = 44; // bins in pt for log binning
194   iBin[1] = 23; // bins in pt for pi0 measurement binning
195   Double_t* binEdges[1];
196   binEdges[0] =  AliHFEtools::MakeLogarithmicBinning(iBin[0], kPtbound[0], kPtbound[1]);
197
198   // bin size is chosen to consider ALICE D measurement
199   const Int_t nptbins = 15;
200   const Int_t nybins = 9;
201   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 
202   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
203   TString hname;
204   for (Int_t iDmeson=0; iDmeson<nbspecies; iDmeson++){
205      hname = "Dmeson"+kDspecies[iDmeson];
206      fhD[iDmeson] = new TH2F(hname,hname+";p_{T} (GeV/c)",nptbins,xbins,nybins,ybins);
207      if(fQAhistos) fQAhistos->Add(fhD[iDmeson]);
208   }
209
210   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
212   Int_t kNcent;
213   if(fIsPbPb) kNcent=11;
214   else
215   {
216       if(fIsppMultiBin) kNcent=8;
217       else kNcent = 1;
218   }
219
220   fMCQACollection = new AliHFEcollection("TaskMCQA", "MC QA histos for meason pt spectra");
221
222   for(Int_t centbin=0; centbin<kNcent; centbin++)
223   {
224       fMCQACollection->CreateTH1Farray(Form("pionspectra_centrbin%i",centbin), "pion yields: MC p_{t} ", iBin[1],kPtRange);
225       fMCQACollection->CreateTH1Farray(Form("etaspectra_centrbin%i",centbin), "eta yields: MC p_{t} ", iBin[1],kPtRange);
226       fMCQACollection->CreateTH1Farray(Form("omegaspectra_centrbin%i",centbin), "omega yields: MC p_{t} ", iBin[1],kPtRange);
227       fMCQACollection->CreateTH1Farray(Form("phispectra_centrbin%i",centbin), "phi yields: MC p_{t} ", iBin[1],kPtRange);
228       fMCQACollection->CreateTH1Farray(Form("etapspectra_centrbin%i",centbin), "etap yields: MC p_{t} ", iBin[1],kPtRange);
229       fMCQACollection->CreateTH1Farray(Form("rhospectra_centrbin%i",centbin), "rho yields: MC p_{t} ", iBin[1],kPtRange);
230
231       fMCQACollection->CreateTH1F(Form("pionspectraLog_centrbin%i",centbin), "pion yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
232       fMCQACollection->CreateTH1F(Form("etaspectraLog_centrbin%i",centbin), "eta yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
233       fMCQACollection->CreateTH1F(Form("omegaspectraLog_centrbin%i",centbin), "omega yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
234       fMCQACollection->CreateTH1F(Form("phispectraLog_centrbin%i",centbin), "phi yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
235       fMCQACollection->CreateTH1F(Form("etapspectraLog_centrbin%i",centbin), "etap yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
236       fMCQACollection->CreateTH1F(Form("rhospectraLog_centrbin%i",centbin), "rho yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
237       fMCQACollection->CreateTH1F(Form("kaonspectraLog_centrbin%i",centbin), "kaon yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
238       fMCQACollection->CreateTH1F(Form("k0LspectraLog_centrbin%i",centbin), "k0L yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
239
240       fMCQACollection->CreateTH1F(Form("piondaughters_centrbin%i",centbin), "pion yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
241       fMCQACollection->CreateTH1F(Form("etadaughters_centrbin%i",centbin), "eta yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
242       fMCQACollection->CreateTH1F(Form("omegadaughters_centrbin%i",centbin), "omega yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
243       fMCQACollection->CreateTH1F(Form("phidaughters_centrbin%i",centbin), "phi yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
244       fMCQACollection->CreateTH1F(Form("etapdaughters_centrbin%i",centbin), "etap yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
245       fMCQACollection->CreateTH1F(Form("rhodaughters_centrbin%i",centbin), "rho yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
246   }
247
248   fQAhistos->Add(fMCQACollection->GetList());
249
250   if(!fTreeStream && fIsDebugStreamerON){
251    fTreeStream = new TTreeSRedirector(Form("HFEmcqadebugTree%s.root", GetName()));
252   }
253
254 }
255   
256 //__________________________________________
257 void AliHFEmcQA::CreateHistograms(const Int_t kquark) 
258 {
259   // create histograms
260
261   if (!(kquark == kCharm || kquark == kBeauty || kquark == kOthers)) {
262     AliDebug(1, "This task is only for heavy quark QA, return\n");
263     return; 
264   }
265   Int_t iq = kquark - kCharm; 
266
267   TString kqTypeLabel[fgkqType];
268   if (kquark == kCharm){
269     kqTypeLabel[kQuark]="c";
270     kqTypeLabel[kantiQuark]="cbar";
271     kqTypeLabel[kHadron]="cHadron";
272     kqTypeLabel[keHadron]="ceHadron";
273     kqTypeLabel[kDeHadron]="nullHadron";
274     kqTypeLabel[kElectron]="ce";
275     kqTypeLabel[kElectron2nd]="nulle";
276   } else if (kquark == kBeauty){
277     kqTypeLabel[kQuark]="b";
278     kqTypeLabel[kantiQuark]="bbar";
279     kqTypeLabel[kHadron]="bHadron";
280     kqTypeLabel[keHadron]="beHadron";
281     kqTypeLabel[kDeHadron]="bDeHadron";
282     kqTypeLabel[kElectron]="be";
283     kqTypeLabel[kElectron2nd]="bce";
284   } else if (kquark == kOthers){
285     kqTypeLabel[kGamma-4]="gammae";
286     kqTypeLabel[kPi0-4]="pi0e";
287     kqTypeLabel[kElse-4]="elsee";
288     kqTypeLabel[kMisID-4]="miside";
289   }
290
291   TString kqEtaRangeLabel[fgkEtaRanges];
292   kqEtaRangeLabel[0] = "mcqa_";
293   kqEtaRangeLabel[1] = "mcqa_barrel_";
294   kqEtaRangeLabel[2] = "mcqa_unitY_";
295
296   const Double_t kPtbound[2] = {0.1, 20.}; //bin taken for considering inclusive e analysis binning
297   const Int_t nptbinning1 = 35;
298   Int_t iBin[2];
299   iBin[0] = 44; // bins in pt
300   iBin[1] = nptbinning1; // bins in pt
301   Double_t* binEdges[1];
302   binEdges[0] =  AliHFEtools::MakeLogarithmicBinning(iBin[0], kPtbound[0], kPtbound[1]);
303
304   // new binning for final electron analysis
305   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.};
306
307   const Int_t ndptbins = 500;
308   Double_t xcorrbin[ndptbins+1];
309   for (int icorrbin = 0; icorrbin< ndptbins+1; icorrbin++){
310     xcorrbin[icorrbin]=icorrbin*0.1;
311   }
312
313   Int_t fCentrmax = 0;
314   if(!fIsPbPb&&!fIsppMultiBin) fCentrmax=0;
315   else fCentrmax = kCentBins;
316   TString hname; 
317   if(kquark == kOthers){
318    for (Int_t icut = 0; icut < fgkEtaRanges; icut++ ){
319        for (Int_t iqType = 0; iqType < 4; iqType++ ){
320            for(Int_t icentr = 0; icentr < (fCentrmax+1); icentr++)
321            {
322                hname = kqEtaRangeLabel[icut]+"Pt_"+kqTypeLabel[iqType];
323                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);
324                hname = kqEtaRangeLabel[icut]+"Y_"+kqTypeLabel[iqType];
325                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);
326                hname = kqEtaRangeLabel[icut]+"Eta_"+kqTypeLabel[iqType];
327                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);
328                // Fill List
329                if(fQAhistos) fHist[iq][iqType][icut][icentr].FillList(fQAhistos);
330            }
331        }
332    }
333    return;
334   }
335   for (Int_t icut = 0; icut < fgkEtaRanges; icut++ ){
336    for (Int_t iqType = 0; iqType < fgkqType; iqType++ ){
337        if (iqType < keHadron && icut > 0) continue; // don't duplicate histogram for quark and hadron
338         for(Int_t icentr = 0; icentr<(fCentrmax+1); icentr++)
339            {
340                hname = kqEtaRangeLabel[icut]+"PdgCode_"+kqTypeLabel[iqType];
341                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);
342                hname = kqEtaRangeLabel[icut]+"Pt_"+kqTypeLabel[iqType];
343                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
344                hname = kqEtaRangeLabel[icut]+"Y_"+kqTypeLabel[iqType];
345                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);
346                hname = kqEtaRangeLabel[icut]+"Eta_"+kqTypeLabel[iqType];
347                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);
348                // Fill List
349                if(fQAhistos) fHist[iq][iqType][icut][icentr].FillList(fQAhistos);
350            }
351    }
352   }
353
354   for (Int_t icut = 0; icut < fgkEtaRanges; icut++ ){
355     hname = kqEtaRangeLabel[icut]+"PtCorr_"+kqTypeLabel[kQuark];
356     fHistComm[iq][icut].fPtCorr = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); 
357     hname = kqEtaRangeLabel[icut]+"PtCorrDp_"+kqTypeLabel[kQuark];
358     fHistComm[iq][icut].fPtCorrDp= new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1);
359     hname = kqEtaRangeLabel[icut]+"PtCorrD0_"+kqTypeLabel[kQuark];
360     fHistComm[iq][icut].fPtCorrD0 = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1);
361     hname = kqEtaRangeLabel[icut]+"PtCorrDrest_"+kqTypeLabel[kQuark];
362     fHistComm[iq][icut].fPtCorrDrest = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1);
363   
364     hname = kqEtaRangeLabel[icut]+"ePtRatio_"+kqTypeLabel[kQuark];
365     fHistComm[iq][icut].fePtRatio = new TH2F(hname,hname+";p_{T} (GeV/c);momentum fraction",200,0,20,100,0,1);
366     hname = kqEtaRangeLabel[icut]+"DePtRatio_"+kqTypeLabel[kQuark];
367     fHistComm[iq][icut].fDePtRatio = new TH2F(hname,hname+";p_{T} (GeV/c);momentum fraction",100,0,20,100,0,1);
368     hname = kqEtaRangeLabel[icut]+"eDistance_"+kqTypeLabel[kQuark];
369     fHistComm[iq][icut].feDistance= new TH2F(hname,hname+";p_{T} (GeV/c);distance (cm)",100,0,20,200,0,2);
370     hname = kqEtaRangeLabel[icut]+"DeDistance_"+kqTypeLabel[kQuark];
371     fHistComm[iq][icut].fDeDistance= new TH2F(hname,hname+";p_{T} (GeV/c);distance (cm)",100,0,20,200,0,2);
372     if(fQAhistos) fHistComm[iq][icut].FillList(fQAhistos);
373   }
374
375   hname = kqEtaRangeLabel[0]+"Nq_"+kqTypeLabel[kQuark];
376   fHistComm[iq][0].fNq = new TH1F(hname,hname,50,-0.5,49.5);
377   hname = kqEtaRangeLabel[0]+"ProcessID_"+kqTypeLabel[kQuark];
378   fHistComm[iq][0].fProcessID = new TH1F(hname,hname,21,-10.5,10.5);
379   
380 }
381
382 //__________________________________________
383 void AliHFEmcQA::Init()
384 {
385   // called at begining every event
386   
387   for (Int_t i=0; i<2; i++){
388      fIsHeavy[i] = 0;
389   } 
390
391   fNparents = 7;
392
393   fParentSelect[0][0] =  411; //D+  
394   fParentSelect[0][1] =  421; //D0
395   fParentSelect[0][2] =  431; //Ds+
396   fParentSelect[0][3] = 4122; //Lambdac+
397   fParentSelect[0][4] = 4132; //Ksic0
398   fParentSelect[0][5] = 4232; //Ksic+
399   fParentSelect[0][6] = 4332; //OmegaC0
400
401   fParentSelect[1][0] =  511; //B0
402   fParentSelect[1][1] =  521; //B+
403   fParentSelect[1][2] =  531; //Bs0
404   fParentSelect[1][3] = 5122; //Lambdab0
405   fParentSelect[1][4] = 5132; //Ksib-
406   fParentSelect[1][5] = 5232; //Ksib0
407   fParentSelect[1][6] = 5332; //Omegab-
408
409
410 }
411
412 //__________________________________________
413 void AliHFEmcQA::GetMesonKine() 
414 {
415   //
416   // get meson pt spectra
417   //
418
419   AliVParticle *mctrack2 = NULL;
420   AliMCParticle *mctrack0 = NULL;
421   AliVParticle *mctrackdaugt= NULL;
422   AliMCParticle *mctrackd= NULL;
423   Int_t id1=0, id2=0;
424
425
426   if(fCentrality>11) printf("warning centrality out of histogram array limits \n");
427   if(!fIsPbPb&&!fIsppMultiBin) fCentrality=0;
428
429   for(Int_t imc = 0; imc <fMCEvent->GetNumberOfPrimaries(); imc++){
430      if(!(mctrack2 = fMCEvent->GetTrack(imc))) continue;
431      TParticle* mcpart0 = fMCEvent->Stack()->Particle(imc);
432      if(!mcpart0) continue;
433      mctrack0 = dynamic_cast<AliMCParticle *>(mctrack2);
434      if(!mctrack0) continue;
435
436 //     if(!fIsPbPb&&!fIsppMultiBin) fCentrality=0;
437
438      if(abs(mctrack0->PdgCode()) == 111) // pi0 
439        {
440           if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
441             fMCQACollection->Fill(Form("pionspectra_centrbin%i",fCentrality),mctrack0->Pt());
442             fMCQACollection->Fill(Form("pionspectraLog_centrbin%i",fCentrality),mctrack0->Pt());
443           }
444           id1=mctrack0->GetFirstDaughter();
445           id2=mctrack0->GetLastDaughter();
446           if(!((id2-id1)==2)) continue;
447           for(int idx=id1; idx<=id2; idx++){
448             if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
449             if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
450             if(abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
451              fMCQACollection->Fill(Form("piondaughters_centrbin%i",fCentrality),mctrackd->Pt());
452           }
453        }
454      else if(abs(mctrack0->PdgCode()) == 221) // eta 
455        {
456           if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
457             fMCQACollection->Fill(Form("etaspectra_centrbin%i",fCentrality),mctrack0->Pt());
458             fMCQACollection->Fill(Form("etaspectraLog_centrbin%i",fCentrality),mctrack0->Pt());
459           } 
460           id1=mctrack0->GetFirstDaughter();
461           id2=mctrack0->GetLastDaughter();
462           if(!((id2-id1)==2||(id2-id1)==3)) continue;
463           for(int idx=id1; idx<=id2; idx++){
464             if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
465             if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
466             if(abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
467              fMCQACollection->Fill(Form("etadaughters_centrbin%i",fCentrality),mctrackd->Pt());
468           }
469        }
470      else if(abs(mctrack0->PdgCode()) == 223) // omega
471        {
472           if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
473             fMCQACollection->Fill(Form("omegaspectra_centrbin%i",fCentrality),mctrack0->Pt());
474             fMCQACollection->Fill(Form("omegaspectraLog_centrbin%i",fCentrality),mctrack0->Pt());
475           }
476           id1=mctrack0->GetFirstDaughter();
477           id2=mctrack0->GetLastDaughter();
478           if(!((id2-id1)==1||(id2-id1)==2)) continue;
479           for(int idx=id1; idx<=id2; idx++){
480             if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
481             if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
482             if(abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
483              fMCQACollection->Fill(Form("omegadaughters_centrbin%i",fCentrality),mctrackd->Pt());
484           }
485        }
486      else if(abs(mctrack0->PdgCode()) == 333) // phi 
487        {
488           if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
489             fMCQACollection->Fill(Form("phispectra_centrbin%i",fCentrality),mctrack0->Pt());
490             fMCQACollection->Fill(Form("phispectraLog_centrbin%i",fCentrality),mctrack0->Pt());
491           } 
492           id1=mctrack0->GetFirstDaughter();
493           id2=mctrack0->GetLastDaughter();
494           if(!((id2-id1)==1)) continue;
495           for(int idx=id1; idx<=id2; idx++){
496             if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
497             if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
498             if(abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
499              fMCQACollection->Fill(Form("phidaughters_centrbin%i",fCentrality),mctrackd->Pt());
500           }
501        }
502      else if(abs(mctrack0->PdgCode()) == 331) // eta prime
503        {
504           if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
505             fMCQACollection->Fill(Form("etapspectra_centrbin%i",fCentrality),mctrack0->Pt());
506             fMCQACollection->Fill(Form("etapspectraLog_centrbin%i",fCentrality),mctrack0->Pt());
507           }
508           id1=mctrack0->GetFirstDaughter();
509           id2=mctrack0->GetLastDaughter();
510           if(!((id2-id1)==2||(id2-id1)==3)) continue;
511           for(int idx=id1; idx<=id2; idx++){
512             if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
513             if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
514             if(abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
515              fMCQACollection->Fill(Form("etapdaughters_centrbin%i",fCentrality),mctrackd->Pt());
516           }
517        }
518      else if(abs(mctrack0->PdgCode()) == 113) // rho
519        {
520           if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
521             fMCQACollection->Fill(Form("rhospectra_centrbin%i",fCentrality),mctrack0->Pt());
522             fMCQACollection->Fill(Form("rhospectraLog_centrbin%i",fCentrality),mctrack0->Pt());
523           }
524           id1=mctrack0->GetFirstDaughter();
525           id2=mctrack0->GetLastDaughter();
526           if(!((id2-id1)==1)) continue;
527           for(int idx=id1; idx<=id2; idx++){
528             if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
529             if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
530             if(abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
531              fMCQACollection->Fill(Form("rhodaughters_centrbin%i",fCentrality),mctrackd->Pt());
532           }
533        }
534      else if(abs(mctrack0->PdgCode()) == 321) // kaon+-
535        {
536           if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
537             fMCQACollection->Fill(Form("kaonspectraLog_centrbin%i",fCentrality),mctrack0->Pt());
538           }
539        }
540      else if(abs(mctrack0->PdgCode()) == 130) // k0L
541        {
542           if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
543             fMCQACollection->Fill(Form("k0LspectraLog_centrbin%i",fCentrality),mctrack0->Pt());
544           }
545        }
546      }
547
548 }
549 //__________________________________________
550 void AliHFEmcQA::GetQuarkKine(TParticle *part, Int_t iTrack, const Int_t kquark) 
551 {
552   // get heavy quark kinematics
553
554     if (kquark != kCharm && kquark != kBeauty) {
555       AliDebug(1, "This task is only for heavy quark QA, return\n");
556       return; 
557     }
558     Int_t iq = kquark - kCharm; 
559
560     if (iTrack < 0 || !part) { 
561       AliDebug(1, "Stack label is negative or no mcparticle, return\n");
562       return; 
563     }
564
565     if(!fIsPbPb&&!fIsppMultiBin) fCentrality=0;
566
567     AliMCParticle *mctrack = NULL;
568     Int_t partPdgcode = TMath::Abs(part->GetPdgCode());
569
570     // select heavy hadron or not fragmented heavy quark 
571     if ( int(partPdgcode/100.)==kquark || int(partPdgcode/1000.)==kquark || (partPdgcode==kquark && (part->GetNDaughters()==0 && iTrack>5)) ){ 
572
573       TParticle *partMother;
574       Int_t iLabel;
575
576       if (partPdgcode == kquark){ // in case of not fragmented heavy quark  
577         partMother = part; 
578         iLabel = iTrack;
579       } else{ // in case of heavy hadron, start to search for mother heavy parton 
580         iLabel = part->GetFirstMother(); 
581         if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return; 
582         if (iLabel>-1) { partMother = mctrack->Particle(); }
583         else {
584           AliDebug(1, "Stack label is negative, return\n");
585           return; 
586         }
587       }
588
589       // heavy parton selection as a mother of heavy hadron 
590       // if the heavy particle comes from string which is denoted as particle status 12|12|12...12|11,[PYTHIA p.60]
591       // in this case, the mother of heavy particle can be one of the fragmented parton of the string
592       // should I make a condition that partMother should be quark or diquark? -> not necessary
593       if ( abs(partMother->GetPdgCode()) == kquark || (partMother->GetStatusCode() == 11) ){
594       //if ( abs(partMother->GetPdgCode()) == kquark || (partMother->GetStatusCode() == 11 || partMother->GetStatusCode() == 12) ){
595
596         if ( abs(partMother->GetPdgCode()) != kquark ){
597           // search fragmented partons in the same string
598           Bool_t isSameString = kTRUE; 
599           for (Int_t i=1; i<fgkMaxIter; i++){
600              iLabel = iLabel - 1;
601              if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return; 
602              if (iLabel>-1) { partMother = mctrack->Particle(); }
603              else {
604                AliDebug(1, "Stack label is negative, return\n");
605                return; 
606              }
607              if ( abs(partMother->GetPdgCode()) == kquark ) break;
608              if ( partMother->GetStatusCode() != 12 ) isSameString = kFALSE;
609              if (!isSameString) return; 
610           }
611         }
612         AliDebug(1, "Can not find heavy parton of this heavy hadron in the string, return\n");
613         if (abs(partMother->GetPdgCode()) != kquark) return; 
614
615         if (fIsHeavy[iq] >= 50) return;  
616         fHeavyQuark[fIsHeavy[iq]] = partMother;
617         fIsHeavy[iq]++;
618
619         // fill kinematics for heavy parton
620         if (partMother->GetPdgCode() > 0) { // quark
621           fHist[iq][kQuark][0][fCentrality].fPdgCode->Fill(partMother->GetPdgCode());
622           fHist[iq][kQuark][0][fCentrality].fPt->Fill(partMother->Pt());
623           fHist[iq][kQuark][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMother));
624           fHist[iq][kQuark][0][fCentrality].fEta->Fill(partMother->Eta());
625         } else{ // antiquark
626           fHist[iq][kantiQuark][0][fCentrality].fPdgCode->Fill(partMother->GetPdgCode());
627           fHist[iq][kantiQuark][0][fCentrality].fPt->Fill(partMother->Pt());
628           fHist[iq][kantiQuark][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMother));
629           fHist[iq][kantiQuark][0][fCentrality].fEta->Fill(partMother->Eta());
630         }
631
632       } // end of heavy parton slection loop 
633
634     } // end of heavy hadron or quark selection
635
636 }
637
638 //__________________________________________
639 void AliHFEmcQA::EndOfEventAna(const Int_t kquark)
640 {
641   // end of event analysis
642
643   if (kquark != kCharm && kquark != kBeauty) {
644     AliDebug(1, "This task is only for heavy quark QA, return\n");
645     return; 
646   }
647   Int_t iq = kquark - kCharm; 
648
649
650   // # of heavy quark per event
651   AliDebug(1,Form("Number of heavy quark in this event = %d \n",fIsHeavy[iq]));
652   fHistComm[iq][0].fNq->Fill(fIsHeavy[iq]);
653
654   Int_t motherID[fgkMaxGener];
655   Int_t motherType[fgkMaxGener];
656   Int_t motherLabel[fgkMaxGener];
657   Int_t ancestorPdg[fgkMaxGener];
658   Int_t ancestorLabel[fgkMaxGener];
659
660   for (Int_t i = 0; i < fgkMaxGener; i++){ // initialization
661      motherID[i] = 0;
662      motherType[i] = 0;
663      motherLabel[i] = 0;
664      ancestorPdg[i] = 0;
665      ancestorLabel[i] = 0;
666   }
667
668
669   // check history of found heavy quarks
670   for (Int_t i = 0; i < fIsHeavy[iq]; i++){
671
672      if(!fHeavyQuark[i]) return;
673
674      ancestorLabel[0] = i;
675      ancestorPdg[0] = fHeavyQuark[i]->GetPdgCode(); 
676      ancestorLabel[1] = fHeavyQuark[i]->GetFirstMother(); 
677
678      AliDebug(1,Form("pdg code= %d\n",ancestorPdg[0]));
679      AliDebug(1,Form("ancestor label= %d\n",ancestorLabel[1]));
680
681      Int_t ig = 1;
682      while (ancestorLabel[ig] != -1){
683           // in case there is mother, get mother's pdg code and grandmother's label
684           IdentifyMother(ancestorLabel[ig], ancestorPdg[ig], ancestorLabel[ig+1]); 
685           // if mother is still heavy, find again mother's ancestor
686           if (ancestorPdg[ig-1] == ancestorPdg[ig]) {
687             ig++;
688             continue; // if it is from same heavy
689           }
690           // if the heavy's mother is not heavy, check the mother's label to know if it comes from inital or final parton shower
691           if (IsFromInitialShower(ancestorLabel[ig],motherID[i],motherType[i],motherLabel[i])) break;
692           if (IsFromFinalParton(ancestorLabel[ig],motherID[i],motherType[i],motherLabel[i])) break;
693           // if it is not the above case, something is strange
694           ReportStrangeness(motherID[i],motherType[i],motherLabel[i]);
695           break;
696      } 
697      if (ancestorLabel[ig] == -1){ // from hard scattering
698        HardScattering(kquark, motherID[i],motherType[i], motherLabel[i]);
699      }
700
701   } // end of found heavy quark loop
702
703
704   // check process type
705   Int_t processID = 0;
706   for (Int_t i = 0; i < fIsHeavy[iq]; i++){
707      AliDebug(1,Form("Mother ID= %d type= %d label= %d\n",motherID[i],motherType[i],motherLabel[i]));
708   }
709
710
711   Int_t nheavypair = Int_t(fIsHeavy[iq]/2.); 
712   for (Int_t ipair = 0; ipair < nheavypair; ipair++){
713
714      Int_t id1 = ipair*2;
715      Int_t id2 = ipair*2 + 1;
716
717      if (motherType[id1] == 2 && motherType[id2] == 2){
718        if (motherLabel[id1] == motherLabel[id2]) processID = kGluonSplitting; // gluon spliting
719        else processID = -9;
720      }
721      else if (motherType[id1] == -1 && motherType[id2] == -1) {
722        if (motherLabel[id1] == -1 && motherLabel[id2] == -1) {
723          if (motherID[id1] == fgkGluon) processID = kPairCreationFromg; // gluon fusion
724          else processID = kPairCreationFromq; // q-qbar pair creation
725        }
726        else processID = -8;
727      }
728      else if (motherType[id1] == -1 || motherType[id2] == -1) {
729        if ((motherLabel[id1] == -1 || motherLabel[id2] == -1) && (motherLabel[id1]*motherLabel[id2] == -2 || motherLabel[id1]*motherLabel[id2] == -3)) {
730          if(motherID[id1]*motherID[id2] == kquark*fgkGluon) processID = kFlavourExitation; // flavour exitation 
731          else processID = kLightQuarkShower;
732        }
733        else processID = -7;
734      }
735      else if (motherType[id1] == -2 || motherType[id2] == -2) {
736        if (motherLabel[id1] == motherLabel[id2]) processID = kInitialPartonShower; // initial parton shower
737        else processID = -6;
738        
739      }
740      else processID = -5;
741
742      if (nheavypair >1) AliDebug(1,Form("Multi pair found : process ID = %d\n",processID));
743      else fHistComm[iq][0].fProcessID->Fill(processID);
744      AliDebug(1,Form("Process ID = %d\n",processID));
745   } // end of # heavy quark pair loop
746
747 }
748
749 //__________________________________________
750 void AliHFEmcQA::GetHadronKine(TParticle* mcpart, const Int_t kquark)
751 {
752     // decay electron kinematics
753
754     if (kquark != kCharm && kquark != kBeauty) {
755       AliDebug(1, "This task is only for heavy quark QA, return\n");
756       return;
757     }
758     Int_t iq = kquark - kCharm;
759
760     if(!mcpart){
761       AliDebug(1, "no mc particle, return\n");
762       return;
763     }
764
765     Int_t iLabel = mcpart->GetFirstMother();
766     if (iLabel<0){
767       AliDebug(1, "Stack label is negative, return\n");
768       return;
769     }
770
771     if(!fIsPbPb&&!fIsppMultiBin) fCentrality=0;
772
773     TParticle *partCopy = mcpart;
774     Int_t pdgcode = mcpart->GetPdgCode();
775     Int_t pdgcodeCopy = pdgcode;
776
777     AliMCParticle *mctrack = NULL;
778
779     // if the mother is charmed hadron  
780     Bool_t isDirectCharm = kFALSE;
781     if ( int(abs(pdgcode)/100.) == kCharm || int(abs(pdgcode)/1000.) == kCharm ) {
782
783           // iterate until you find B hadron as a mother or become top ancester 
784           for (Int_t i=1; i<fgkMaxIter; i++){
785
786              Int_t jLabel = mcpart->GetFirstMother();
787              if (jLabel == -1){
788                isDirectCharm = kTRUE;
789                break; // if there is no ancester
790              }
791              if (jLabel < 0){ // safety protection
792                AliDebug(1, "Stack label is negative, return\n");
793                return;
794              }
795              // if there is an ancester
796              if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return; 
797              TParticle* mother = mctrack->Particle();
798              Int_t motherPDG = mother->GetPdgCode();
799     
800              for (Int_t j=0; j<fNparents; j++){
801                 if (abs(motherPDG)==fParentSelect[1][j]) return; // return if this hadron is originated from b
802              }
803
804              mcpart = mother;
805           } // end of iteration 
806     } // end of if
807     if((isDirectCharm == kTRUE && kquark == kCharm) || kquark == kBeauty) {
808          for (Int_t i=0; i<fNparents; i++){
809             if (abs(pdgcodeCopy)==fParentSelect[iq][i]){
810
811               // fill hadron kinematics
812               fHist[iq][kHadron][0][fCentrality].fPdgCode->Fill(pdgcodeCopy);
813               fHist[iq][kHadron][0][fCentrality].fPt->Fill(partCopy->Pt());
814               fHist[iq][kHadron][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partCopy));
815               fHist[iq][kHadron][0][fCentrality].fEta->Fill(partCopy->Eta());
816
817               if(iq==0) {
818                fhD[i]->Fill(partCopy->Pt(),AliHFEtools::GetRapidity(partCopy));
819               }
820             }
821          }
822          // I also want to store D* info to compare with D* measurement 
823          if (abs(pdgcodeCopy)==413 && iq==0) { //D*+
824                fhD[7]->Fill(partCopy->Pt(),AliHFEtools::GetRapidity(partCopy));
825          }
826          if (abs(pdgcodeCopy)==423 && iq==0) { //D*0
827                fhD[8]->Fill(partCopy->Pt(),AliHFEtools::GetRapidity(partCopy));
828          }
829     } // end of if
830 }
831
832 //__________________________________________
833 void AliHFEmcQA::GetDecayedKine(TParticle* mcpart, const Int_t kquark, Int_t kdecayed) 
834 {
835     // decay electron kinematics
836     
837     if (!(kquark == kCharm || kquark == kBeauty || kquark == kOthers)){
838       AliDebug(1, "This task is only for heavy quark QA, return\n");
839       return; 
840     }
841     Int_t iq = kquark - kCharm; 
842     Bool_t isFinalOpenCharm = kFALSE;
843
844     if(!mcpart){
845       AliDebug(1, "no mcparticle, return\n");
846       return;
847     }
848
849     if(!fIsPbPb&&!fIsppMultiBin) fCentrality=0;
850
851     Double_t eabsEta = TMath::Abs(mcpart->Eta());
852     Double_t eabsY = TMath::Abs(AliHFEtools::GetRapidity(mcpart));
853
854     if(kquark==kOthers){
855       Int_t esource = -1;
856       if ( abs(mcpart->GetPdgCode()) != kdecayed ) esource = kMisID-4;
857       else esource =GetSource(mcpart)-4; // return for the cases kGamma=4, kPi0=5, kElse=6
858       if(esource==0|| esource==1 || esource==2 || esource==3){
859         fHist[iq][esource][0][fCentrality].fPt->Fill(mcpart->Pt());
860         fHist[iq][esource][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
861         fHist[iq][esource][0][fCentrality].fEta->Fill(mcpart->Eta());
862         if(eabsEta<0.9){
863           fHist[iq][esource][1][fCentrality].fPt->Fill(mcpart->Pt());
864           fHist[iq][esource][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
865           fHist[iq][esource][1][fCentrality].fEta->Fill(mcpart->Eta());
866         }
867         if(eabsY<0.5){
868           fHist[iq][esource][2][fCentrality].fPt->Fill(mcpart->Pt());
869           fHist[iq][esource][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
870           fHist[iq][esource][2][fCentrality].fEta->Fill(mcpart->Eta());
871         }
872         return; 
873       }
874       else {
875         AliDebug(1, "e source is out of defined ranges, return\n");
876         return;
877       }
878     }
879
880     if ( abs(mcpart->GetPdgCode()) != kdecayed ) return;
881
882     Int_t iLabel = mcpart->GetFirstMother(); 
883     if (iLabel<0){
884       AliDebug(1, "Stack label is negative, return\n");
885       return; 
886     }
887
888     AliMCParticle *mctrack = NULL;
889     if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return; 
890     TParticle *partMother = mctrack->Particle();
891     TParticle *partMotherCopy = partMother;
892     Int_t maPdgcode = partMother->GetPdgCode();
893     Int_t maPdgcodeCopy = maPdgcode;
894
895     // get mc primary vertex
896     /*
897     TArrayF mcPrimVtx;
898     if(fMCHeader) fMCHeader->PrimaryVertex(mcPrimVtx);
899
900     // get electron production vertex   
901     TLorentzVector ePoint;
902     mcpart->ProductionVertex(ePoint);
903
904     // calculated production vertex to primary vertex (in xy plane)
905     Float_t decayLxy = TMath::Sqrt((mcPrimVtx[0]-ePoint[0])*(mcPrimVtx[0]-ePoint[0])+(mcPrimVtx[1]-ePoint[1])*(mcPrimVtx[1]-ePoint[1]));
906     */ 
907     Float_t decayLxy = 0;
908
909     // if the mother is charmed hadron  
910     Bool_t isMotherDirectCharm = kFALSE;
911     if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) { 
912
913          for (Int_t i=0; i<fNparents; i++){
914             if (abs(maPdgcode)==fParentSelect[0][i]){
915               isFinalOpenCharm = kTRUE;
916             } 
917          }  
918          if (!isFinalOpenCharm) return ;
919
920           // iterate until you find B hadron as a mother or become top ancester 
921           for (Int_t i=1; i<fgkMaxIter; i++){
922
923              Int_t jLabel = partMother->GetFirstMother(); 
924              if (jLabel == -1){
925                isMotherDirectCharm = kTRUE;
926                break; // if there is no ancester
927              }
928              if (jLabel < 0){ // safety protection
929                AliDebug(1, "Stack label is negative, return\n");
930                return; 
931              }
932
933              // if there is an ancester
934              if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return; 
935              TParticle* grandMa = mctrack->Particle();
936              Int_t grandMaPDG = grandMa->GetPdgCode();
937
938              for (Int_t j=0; j<fNparents; j++){
939                 if (abs(grandMaPDG)==fParentSelect[1][j]){
940
941                   if (kquark == kCharm) return;
942                   // fill electron kinematics
943                   fHist[iq][kElectron2nd][0][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
944                   fHist[iq][kElectron2nd][0][fCentrality].fPt->Fill(mcpart->Pt());
945                   fHist[iq][kElectron2nd][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
946                   fHist[iq][kElectron2nd][0][fCentrality].fEta->Fill(mcpart->Eta());
947
948                   // fill mother hadron kinematics
949                   fHist[iq][kDeHadron][0][fCentrality].fPdgCode->Fill(grandMaPDG); 
950                   fHist[iq][kDeHadron][0][fCentrality].fPt->Fill(grandMa->Pt());
951                   fHist[iq][kDeHadron][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(grandMa));
952                   fHist[iq][kDeHadron][0][fCentrality].fEta->Fill(grandMa->Eta());
953
954                   if(eabsEta<0.9){
955                     fHist[iq][kElectron2nd][1][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
956                     fHist[iq][kElectron2nd][1][fCentrality].fPt->Fill(mcpart->Pt());
957                     fHist[iq][kElectron2nd][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
958                     fHist[iq][kElectron2nd][1][fCentrality].fEta->Fill(mcpart->Eta());
959
960                     // fill mother hadron kinematics
961                     fHist[iq][kDeHadron][1][fCentrality].fPdgCode->Fill(grandMaPDG); 
962                     fHist[iq][kDeHadron][1][fCentrality].fPt->Fill(grandMa->Pt());
963                     fHist[iq][kDeHadron][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(grandMa));
964                     fHist[iq][kDeHadron][1][fCentrality].fEta->Fill(grandMa->Eta());
965                   }
966
967                   if(eabsY<0.5){
968                     fHist[iq][kElectron2nd][2][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
969                     fHist[iq][kElectron2nd][2][fCentrality].fPt->Fill(mcpart->Pt());
970                     fHist[iq][kElectron2nd][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
971                     fHist[iq][kElectron2nd][2][fCentrality].fEta->Fill(mcpart->Eta());
972
973                     // fill mother hadron kinematics
974                     fHist[iq][kDeHadron][2][fCentrality].fPdgCode->Fill(grandMaPDG); 
975                     fHist[iq][kDeHadron][2][fCentrality].fPt->Fill(grandMa->Pt());
976                     fHist[iq][kDeHadron][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(grandMa));
977                     fHist[iq][kDeHadron][2][fCentrality].fEta->Fill(grandMa->Eta());
978                   }
979
980                   // ratio between pT of electron and pT of mother B hadron 
981                   if(grandMa->Pt()) {
982                     fHistComm[iq][0].fDePtRatio->Fill(grandMa->Pt(),mcpart->Pt()/grandMa->Pt());
983                     if(eabsEta<0.9){
984                       fHistComm[iq][1].fDePtRatio->Fill(grandMa->Pt(),mcpart->Pt()/grandMa->Pt());
985                     }
986                     if(eabsY<0.5){
987                       fHistComm[iq][2].fDePtRatio->Fill(grandMa->Pt(),mcpart->Pt()/grandMa->Pt());
988                     }
989                   }
990
991                   // distance between electron production point and primary vertex
992                   fHistComm[iq][0].fDeDistance->Fill(grandMa->Pt(),decayLxy);
993                   if(eabsEta<0.9){
994                     fHistComm[iq][1].fDeDistance->Fill(grandMa->Pt(),decayLxy);
995                   }
996                   if(eabsY<0.5){
997                     fHistComm[iq][2].fDeDistance->Fill(grandMa->Pt(),decayLxy);
998                   }
999                   return;
1000                 }
1001              } 
1002
1003              partMother = grandMa;
1004           } // end of iteration 
1005     } // end of if
1006     if((isMotherDirectCharm == kTRUE && kquark == kCharm) || kquark == kBeauty) {
1007          for (Int_t i=0; i<fNparents; i++){
1008             if (abs(maPdgcodeCopy)==fParentSelect[iq][i]){
1009
1010               fHist[iq][kElectron][0][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1011               fHist[iq][kElectron][0][fCentrality].fPt->Fill(mcpart->Pt());
1012               fHist[iq][kElectron][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1013               fHist[iq][kElectron][0][fCentrality].fEta->Fill(mcpart->Eta());  
1014
1015               // fill mother hadron kinematics
1016               fHist[iq][keHadron][0][fCentrality].fPdgCode->Fill(maPdgcodeCopy); 
1017               fHist[iq][keHadron][0][fCentrality].fPt->Fill(partMotherCopy->Pt());
1018               fHist[iq][keHadron][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
1019               fHist[iq][keHadron][0][fCentrality].fEta->Fill(partMotherCopy->Eta());
1020
1021               if(eabsEta<0.9){
1022                 fHist[iq][kElectron][1][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1023                 fHist[iq][kElectron][1][fCentrality].fPt->Fill(mcpart->Pt());
1024                 fHist[iq][kElectron][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1025                 fHist[iq][kElectron][1][fCentrality].fEta->Fill(mcpart->Eta());  
1026
1027                 // fill mother hadron kinematics
1028                 fHist[iq][keHadron][1][fCentrality].fPdgCode->Fill(maPdgcodeCopy); 
1029                 fHist[iq][keHadron][1][fCentrality].fPt->Fill(partMotherCopy->Pt());
1030                 fHist[iq][keHadron][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
1031                 fHist[iq][keHadron][1][fCentrality].fEta->Fill(partMotherCopy->Eta());
1032               }
1033
1034               if(eabsY<0.5){
1035                 fHist[iq][kElectron][2][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1036                 fHist[iq][kElectron][2][fCentrality].fPt->Fill(mcpart->Pt());
1037                 fHist[iq][kElectron][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1038                 fHist[iq][kElectron][2][fCentrality].fEta->Fill(mcpart->Eta());  
1039
1040                 // fill mother hadron kinematics
1041                 fHist[iq][keHadron][2][fCentrality].fPdgCode->Fill(maPdgcodeCopy); 
1042                 fHist[iq][keHadron][2][fCentrality].fPt->Fill(partMotherCopy->Pt());
1043                 fHist[iq][keHadron][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
1044                 fHist[iq][keHadron][2][fCentrality].fEta->Fill(partMotherCopy->Eta());
1045               }
1046
1047               // ratio between pT of electron and pT of mother B or direct D hadron 
1048               if(partMotherCopy->Pt()) {
1049                  fHistComm[iq][0].fePtRatio->Fill(partMotherCopy->Pt(),mcpart->Pt()/partMotherCopy->Pt());
1050                  if(eabsEta<0.9){
1051                    fHistComm[iq][1].fePtRatio->Fill(partMotherCopy->Pt(),mcpart->Pt()/partMotherCopy->Pt());
1052                  }
1053                  if(eabsY<0.5){
1054                    fHistComm[iq][2].fePtRatio->Fill(partMotherCopy->Pt(),mcpart->Pt()/partMotherCopy->Pt());
1055                  }
1056               }
1057               fHistComm[iq][0].fPtCorr->Fill(partMotherCopy->Pt(),mcpart->Pt());
1058               if(eabsEta<0.9){
1059                 fHistComm[iq][1].fPtCorr->Fill(partMotherCopy->Pt(),mcpart->Pt());
1060               }
1061               if(eabsY<0.5){
1062                 fHistComm[iq][2].fPtCorr->Fill(partMotherCopy->Pt(),mcpart->Pt());
1063               }
1064               if(TMath::Abs(partMotherCopy->GetPdgCode())==411) {
1065                 fHistComm[iq][0].fPtCorrDp->Fill(partMotherCopy->Pt(),mcpart->Pt());
1066                 if(eabsEta<0.9){
1067                   fHistComm[iq][1].fPtCorrDp->Fill(partMotherCopy->Pt(),mcpart->Pt());
1068                 }
1069                 if(eabsY<0.5){
1070                   fHistComm[iq][2].fPtCorrDp->Fill(partMotherCopy->Pt(),mcpart->Pt());
1071                 }
1072               }
1073               else if(TMath::Abs(partMotherCopy->GetPdgCode())==421) {
1074                 fHistComm[iq][0].fPtCorrD0->Fill(partMotherCopy->Pt(),mcpart->Pt());
1075                 if(eabsEta<0.9){
1076                   fHistComm[iq][1].fPtCorrD0->Fill(partMotherCopy->Pt(),mcpart->Pt());
1077                 }
1078                 if(eabsY<0.5){
1079                   fHistComm[iq][2].fPtCorrD0->Fill(partMotherCopy->Pt(),mcpart->Pt());
1080                 }
1081               }
1082               else {
1083                 fHistComm[iq][0].fPtCorrDrest->Fill(partMotherCopy->Pt(),mcpart->Pt());
1084                 if(eabsEta<0.9){
1085                   fHistComm[iq][1].fPtCorrDrest->Fill(partMotherCopy->Pt(),mcpart->Pt());
1086                 }
1087                 if(eabsY<0.5){
1088                   fHistComm[iq][2].fPtCorrDrest->Fill(partMotherCopy->Pt(),mcpart->Pt());
1089                 }
1090               }
1091
1092               // distance between electron production point and primary vertex
1093               fHistComm[iq][0].feDistance->Fill(partMotherCopy->Pt(),decayLxy);
1094               if(eabsEta<0.9){
1095                 fHistComm[iq][1].feDistance->Fill(partMotherCopy->Pt(),decayLxy);
1096               }
1097               if(eabsY<0.5){
1098                 fHistComm[iq][2].feDistance->Fill(partMotherCopy->Pt(),decayLxy);
1099               }
1100             }
1101          }
1102     } // end of if
1103 }
1104
1105 //____________________________________________________________________
1106 void  AliHFEmcQA::GetDecayedKine(AliAODMCParticle *mcpart, const Int_t kquark, Int_t kdecayed)
1107 {
1108   // decay electron kinematics
1109
1110   if (kquark != kCharm && kquark != kBeauty) {
1111     AliDebug(1, "This task is only for heavy quark QA, return\n");
1112     return;
1113   }
1114
1115   Int_t iq = kquark - kCharm;
1116   Bool_t isFinalOpenCharm = kFALSE;
1117
1118   if(!mcpart){
1119     AliDebug(1, "no mcparticle, return\n");
1120     return;
1121   }
1122
1123   if ( abs(mcpart->GetPdgCode()) != kdecayed ) return;
1124
1125   Double_t eabsEta = TMath::Abs(mcpart->Eta());
1126   Double_t eabsY = TMath::Abs(AliHFEtools::GetRapidity(mcpart));
1127
1128   // mother
1129   Int_t iLabel = mcpart->GetMother();
1130   if (iLabel<0){
1131     AliDebug(1, "Stack label is negative, return\n");
1132     return;
1133   }
1134
1135   if(!fIsPbPb&&!fIsppMultiBin) fCentrality=0;
1136
1137   AliAODMCParticle *partMother = (AliAODMCParticle*)fMCArray->At(iLabel);
1138   AliAODMCParticle *partMotherCopy = partMother;
1139   Int_t maPdgcode = partMother->GetPdgCode();
1140   Int_t maPdgcodeCopy = maPdgcode;
1141
1142   Bool_t isMotherDirectCharm = kFALSE;
1143   if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
1144
1145     for (Int_t i=0; i<fNparents; i++){
1146        if (abs(maPdgcode)==fParentSelect[0][i]){
1147          isFinalOpenCharm = kTRUE;
1148        }
1149     } 
1150     if (!isFinalOpenCharm) return;
1151
1152     for (Int_t i=1; i<fgkMaxIter; i++){
1153
1154        Int_t jLabel = partMother->GetMother();
1155        if (jLabel == -1){
1156          isMotherDirectCharm = kTRUE;
1157          break; // if there is no ancester
1158        }
1159        if (jLabel < 0){ // safety protection
1160          AliDebug(1, "Stack label is negative, return\n");
1161          return;
1162        }
1163
1164        // if there is an ancester
1165        AliAODMCParticle* grandMa = (AliAODMCParticle*)fMCArray->At(jLabel);
1166        Int_t grandMaPDG = grandMa->GetPdgCode();
1167
1168        for (Int_t j=0; j<fNparents; j++){
1169           if (abs(grandMaPDG)==fParentSelect[1][j]){
1170
1171             if (kquark == kCharm) return;
1172             // fill electron kinematics
1173             fHist[iq][kElectron2nd][0][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1174             fHist[iq][kElectron2nd][0][fCentrality].fPt->Fill(mcpart->Pt());
1175             fHist[iq][kElectron2nd][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1176             fHist[iq][kElectron2nd][0][fCentrality].fEta->Fill(mcpart->Eta());
1177
1178             // fill mother hadron kinematics
1179             fHist[iq][kDeHadron][0][fCentrality].fPdgCode->Fill(grandMaPDG);
1180             fHist[iq][kDeHadron][0][fCentrality].fPt->Fill(grandMa->Pt());
1181             fHist[iq][kDeHadron][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(grandMa));
1182             fHist[iq][kDeHadron][0][fCentrality].fEta->Fill(grandMa->Eta());
1183
1184             if(eabsEta<0.9){
1185               // fill electron kinematics
1186               fHist[iq][kElectron2nd][1][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1187               fHist[iq][kElectron2nd][1][fCentrality].fPt->Fill(mcpart->Pt());
1188               fHist[iq][kElectron2nd][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1189               fHist[iq][kElectron2nd][1][fCentrality].fEta->Fill(mcpart->Eta());
1190
1191               // fill mother hadron kinematics
1192               fHist[iq][kDeHadron][1][fCentrality].fPdgCode->Fill(grandMaPDG);
1193               fHist[iq][kDeHadron][1][fCentrality].fPt->Fill(grandMa->Pt());
1194               fHist[iq][kDeHadron][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(grandMa));
1195               fHist[iq][kDeHadron][1][fCentrality].fEta->Fill(grandMa->Eta());
1196             }
1197             if(eabsY<0.5){
1198               // fill electron kinematics
1199               fHist[iq][kElectron2nd][2][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1200               fHist[iq][kElectron2nd][2][fCentrality].fPt->Fill(mcpart->Pt());
1201               fHist[iq][kElectron2nd][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1202               fHist[iq][kElectron2nd][2][fCentrality].fEta->Fill(mcpart->Eta());
1203
1204               // fill mother hadron kinematics
1205               fHist[iq][kDeHadron][2][fCentrality].fPdgCode->Fill(grandMaPDG);
1206               fHist[iq][kDeHadron][2][fCentrality].fPt->Fill(grandMa->Pt());
1207               fHist[iq][kDeHadron][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(grandMa));
1208               fHist[iq][kDeHadron][2][fCentrality].fEta->Fill(grandMa->Eta());
1209             }
1210
1211             return;
1212           }
1213        }
1214
1215        partMother = grandMa;
1216     } // end of iteration 
1217   } // end of if
1218   if ((isMotherDirectCharm == kTRUE && kquark == kCharm) || kquark == kBeauty) {
1219     for (Int_t i=0; i<fNparents; i++){
1220        if (abs(maPdgcodeCopy)==fParentSelect[iq][i]){
1221
1222          // fill electron kinematics
1223          fHist[iq][kElectron][0][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1224          fHist[iq][kElectron][0][fCentrality].fPt->Fill(mcpart->Pt());
1225          fHist[iq][kElectron][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1226          fHist[iq][kElectron][0][fCentrality].fEta->Fill(mcpart->Eta());
1227
1228          // fill mother hadron kinematics
1229          fHist[iq][keHadron][0][fCentrality].fPdgCode->Fill(maPdgcodeCopy);
1230          fHist[iq][keHadron][0][fCentrality].fPt->Fill(partMotherCopy->Pt());
1231          fHist[iq][keHadron][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
1232          fHist[iq][keHadron][0][fCentrality].fEta->Fill(partMotherCopy->Eta());
1233
1234          if(eabsEta<0.9){
1235            // fill electron kinematics
1236            fHist[iq][kElectron][1][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1237            fHist[iq][kElectron][1][fCentrality].fPt->Fill(mcpart->Pt());
1238            fHist[iq][kElectron][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1239            fHist[iq][kElectron][1][fCentrality].fEta->Fill(mcpart->Eta());
1240
1241            // fill mother hadron kinematics
1242            fHist[iq][keHadron][1][fCentrality].fPdgCode->Fill(maPdgcodeCopy);
1243            fHist[iq][keHadron][1][fCentrality].fPt->Fill(partMotherCopy->Pt());
1244            fHist[iq][keHadron][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
1245            fHist[iq][keHadron][1][fCentrality].fEta->Fill(partMotherCopy->Eta());
1246          }
1247          if(eabsY<0.5){
1248            // fill electron kinematics
1249            fHist[iq][kElectron][2][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1250            fHist[iq][kElectron][2][fCentrality].fPt->Fill(mcpart->Pt());
1251            fHist[iq][kElectron][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1252            fHist[iq][kElectron][2][fCentrality].fEta->Fill(mcpart->Eta());
1253
1254            // fill mother hadron kinematics
1255            fHist[iq][keHadron][2][fCentrality].fPdgCode->Fill(maPdgcodeCopy);
1256            fHist[iq][keHadron][2][fCentrality].fPt->Fill(partMotherCopy->Pt());
1257            fHist[iq][keHadron][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
1258            fHist[iq][keHadron][2][fCentrality].fEta->Fill(partMotherCopy->Eta());
1259          }
1260
1261        }
1262     }
1263   } // end of if
1264
1265 }
1266
1267 //__________________________________________
1268 void AliHFEmcQA::IdentifyMother(Int_t motherlabel, Int_t &motherpdg, Int_t &grandmotherlabel)
1269 {
1270        // find mother pdg code and label 
1271
1272        if (motherlabel < 0) { 
1273          AliDebug(1, "Stack label is negative, return\n");
1274          return; 
1275        }
1276        AliMCParticle *mctrack = NULL;
1277        if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(motherlabel))))) return; 
1278        TParticle *heavysMother = mctrack->Particle();
1279        motherpdg = heavysMother->GetPdgCode();
1280        grandmotherlabel = heavysMother->GetFirstMother();
1281        AliDebug(1,Form("ancestor pdg code= %d\n",motherpdg));
1282 }
1283
1284 //__________________________________________
1285 void AliHFEmcQA::HardScattering(const Int_t kquark, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
1286 {
1287        // mothertype -1 means this heavy quark coming from hard vertex
1288
1289        AliMCParticle *mctrack1 = NULL;
1290        AliMCParticle *mctrack2 = NULL;
1291        if(!(mctrack1 = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(4))))) return; 
1292        if(!(mctrack2 = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(5))))) return; 
1293        TParticle *afterinitialrad1  = mctrack1->Particle();
1294        TParticle *afterinitialrad2  = mctrack2->Particle();
1295            
1296        motherlabel = -1;
1297
1298        if (abs(afterinitialrad1->GetPdgCode()) == fgkGluon && abs(afterinitialrad2->GetPdgCode()) == fgkGluon){
1299          AliDebug(1,"heavy from gluon gluon pair creation!\n");
1300          mothertype = -1;
1301          motherID = fgkGluon;
1302        }
1303        else if (abs(afterinitialrad1->GetPdgCode()) == kquark || abs(afterinitialrad2->GetPdgCode()) == kquark){ // one from Q and the other from g
1304          AliDebug(1,"heavy from flavor exitation!\n");
1305          mothertype = -1;
1306          motherID = kquark;
1307        }
1308        else if  (abs(afterinitialrad1->GetPdgCode()) == abs(afterinitialrad2->GetPdgCode())){
1309          AliDebug(1,"heavy from q-qbar pair creation!\n");
1310          mothertype = -1;
1311          motherID = 1;
1312        }
1313        else {
1314          AliDebug(1,"something strange!\n");
1315          mothertype = -999;
1316          motherlabel = -999;
1317          motherID = -999;
1318        }
1319 }
1320
1321 //__________________________________________
1322 Bool_t AliHFEmcQA::IsFromInitialShower(Int_t inputmotherlabel, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
1323 {
1324        // mothertype -2 means this heavy quark coming from initial state 
1325
1326        AliMCParticle *mctrack = NULL;
1327        if (inputmotherlabel==2 || inputmotherlabel==3){ // mother exist before initial state radiation
1328          if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(inputmotherlabel))))) return kFALSE; 
1329          TParticle *heavysMother = mctrack->Particle();
1330          motherID = heavysMother->GetPdgCode(); 
1331          mothertype = -2; // there is mother before initial state radiation
1332          motherlabel = inputmotherlabel;
1333          AliDebug(1,"initial parton shower! \n");
1334
1335          return kTRUE;
1336        }
1337
1338        return kFALSE;
1339 }
1340
1341 //__________________________________________
1342 Bool_t AliHFEmcQA::IsFromFinalParton(Int_t inputmotherlabel, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
1343 {
1344        // mothertype 2 means this heavy quark coming from final state 
1345
1346        AliMCParticle *mctrack = NULL;
1347        if (inputmotherlabel > 5){ // mother exist after hard scattering
1348          if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(inputmotherlabel))))) return kFALSE; 
1349          TParticle *heavysMother = mctrack->Particle();
1350          motherID = heavysMother->GetPdgCode(); 
1351          mothertype = 2; // 
1352          motherlabel = inputmotherlabel;
1353          AliDebug(1,Form("heavy quark from %d after hard scattering! \n",motherID));
1354
1355          return kTRUE;
1356        }
1357        return kFALSE;
1358 }
1359
1360 //__________________________________________
1361 void AliHFEmcQA::ReportStrangeness(Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
1362 {
1363       // mark strange behavior  
1364
1365        mothertype = -888;
1366        motherlabel = -888;
1367        motherID = -888;
1368        AliDebug(1,"something strange!\n");
1369 }
1370
1371 //__________________________________________
1372 Int_t AliHFEmcQA::GetSource(const AliVParticle* const mcpart)
1373 {        
1374   // decay particle's origin 
1375
1376   //if ( abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return -1;
1377        
1378   Int_t origin = -1;
1379   Bool_t isFinalOpenCharm = kFALSE;
1380
1381   if(!mcpart){
1382     AliDebug(1, "Stack label is negative or no mcparticle, return\n");
1383     return -1;
1384   }
1385
1386   // mother
1387   // Information not in the base class, cast necessary
1388   Int_t iLabel = GetMother(mcpart);
1389   if (iLabel<0){
1390     AliDebug(1, "Stack label is negative, return\n");
1391     return -1;
1392   } 
1393        
1394   const AliVParticle *partMother = fMCEvent->GetTrack(iLabel);
1395   Int_t maPdgcode = partMother->PdgCode();
1396   
1397   // if the mother is charmed hadron  
1398   if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
1399     
1400     for (Int_t i=0; i<fNparents; i++){
1401        if (abs(maPdgcode)==fParentSelect[0][i]){
1402          isFinalOpenCharm = kTRUE;
1403        }
1404     }
1405     if (!isFinalOpenCharm) return -1;
1406
1407     // iterate until you find B hadron as a mother or become top ancester 
1408     for (Int_t i=1; i<fgkMaxIter; i++){
1409       
1410        Int_t jLabel = GetMother(partMother);
1411        if (jLabel == -1){
1412          origin = kDirectCharm;
1413          return origin;
1414        }
1415        if (jLabel < 0){ // safety protection
1416          AliDebug(1, "Stack label is negative, return\n");
1417          return -1;
1418        }
1419
1420        // if there is an ancester
1421        const AliVParticle* grandMa = fMCEvent->GetTrack(jLabel);
1422        Int_t grandMaPDG = grandMa->PdgCode();
1423
1424        for (Int_t j=0; j<fNparents; j++){
1425           if (abs(grandMaPDG)==fParentSelect[1][j]){
1426             origin = kBeautyCharm;
1427             return origin;
1428           }
1429        }
1430
1431        partMother = grandMa;
1432     } // end of iteration 
1433   } // end of if
1434   else if ( int(abs(maPdgcode)/100.) == kBeauty || int(abs(maPdgcode)/1000.) == kBeauty ) {
1435     for (Int_t i=0; i<fNparents; i++){
1436        if (abs(maPdgcode)==fParentSelect[1][i]){
1437          origin = kDirectBeauty;
1438          return origin;
1439        }
1440     }
1441   } // end of if
1442   else if ( abs(maPdgcode) == 22 ) {
1443     origin = kGamma;
1444     return origin;
1445   } // end of if
1446   else if ( abs(maPdgcode) == 111 ) {
1447     origin = kPi0;
1448     return origin;
1449   } // end of if
1450
1451   return origin;
1452 }
1453
1454 //__________________________________________
1455 Int_t AliHFEmcQA::GetSource(const TParticle * const mcpart)
1456 {
1457   // decay particle's origin 
1458
1459   //if ( abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return -1;
1460
1461   Int_t origin = -1;
1462   Bool_t isFinalOpenCharm = kFALSE;
1463
1464   if(!mcpart){
1465     AliDebug(1, "no mcparticle, return\n");
1466     return -1;
1467   }
1468
1469   Int_t iLabel = mcpart->GetFirstMother();
1470   if (iLabel<0){
1471     AliDebug(1, "Stack label is negative, return\n");
1472     return -1;
1473   }
1474
1475   AliMCParticle *mctrack = NULL;
1476   if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return -1; 
1477   TParticle *partMother = mctrack->Particle();
1478   Int_t maPdgcode = partMother->GetPdgCode();
1479
1480    // if the mother is charmed hadron  
1481    if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
1482
1483      for (Int_t i=0; i<fNparents; i++){
1484         if (abs(maPdgcode)==fParentSelect[0][i]){
1485           isFinalOpenCharm = kTRUE;
1486         }
1487      }
1488      if (!isFinalOpenCharm) return -1;
1489
1490      // iterate until you find B hadron as a mother or become top ancester 
1491      for (Int_t i=1; i<fgkMaxIter; i++){
1492
1493         Int_t jLabel = partMother->GetFirstMother();
1494         if (jLabel == -1){
1495           origin = kDirectCharm;
1496           return origin;
1497         }
1498         if (jLabel < 0){ // safety protection
1499           AliDebug(1, "Stack label is negative, return\n");
1500           return -1;
1501         }
1502
1503         // if there is an ancester
1504         if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return -1; 
1505         TParticle* grandMa = mctrack->Particle();
1506         Int_t grandMaPDG = grandMa->GetPdgCode();
1507
1508         for (Int_t j=0; j<fNparents; j++){
1509            if (abs(grandMaPDG)==fParentSelect[1][j]){
1510              origin = kBeautyCharm;
1511              return origin;
1512            }
1513         }
1514
1515         partMother = grandMa;
1516      } // end of iteration 
1517    } // end of if
1518    else if ( int(abs(maPdgcode)/100.) == kBeauty || int(abs(maPdgcode)/1000.) == kBeauty ) {
1519      for (Int_t i=0; i<fNparents; i++){
1520         if (abs(maPdgcode)==fParentSelect[1][i]){
1521           origin = kDirectBeauty;
1522           return origin;
1523         }
1524      }
1525    } // end of if
1526    else if ( abs(maPdgcode) == 22 ) {
1527      origin = kGamma;
1528      return origin;
1529    } // end of if
1530    else if ( abs(maPdgcode) == 111 ) {
1531      origin = kPi0;
1532      return origin;
1533    } // end of if
1534    else origin = kElse;
1535
1536    return origin;
1537 }
1538
1539 //__________________________________________
1540 Int_t AliHFEmcQA::GetElecSource(TParticle * const mcpart)
1541 {
1542   // decay particle's origin 
1543
1544   if(!mcpart){
1545     AliDebug(1, "no mcparticle, return\n");
1546     return -1;
1547   }
1548
1549   if ( abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return kMisID;
1550
1551   Int_t origin = -1;
1552   Bool_t isFinalOpenCharm = kFALSE;
1553
1554   Int_t iLabel = mcpart->GetFirstMother();
1555   if (iLabel<0){
1556     AliDebug(1, "Stack label is negative, return\n");
1557     return -1;
1558   }
1559
1560   AliMCParticle *mctrack = NULL;
1561   Int_t tmpMomLabel=0;
1562   if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return -1; 
1563   TParticle *partMother = mctrack->Particle();
1564   TParticle *partMotherCopy = mctrack->Particle();
1565   Int_t maPdgcode = partMother->GetPdgCode();
1566
1567    // if the mother is charmed hadron  
1568
1569    if(abs(maPdgcode)==443){ // J/spi
1570       Int_t jLabel = partMother->GetFirstMother();
1571       if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return kJpsi; 
1572       TParticle* grandMa = mctrack->Particle();
1573       Int_t grandMaPDG = grandMa->GetPdgCode();
1574       if((int(abs(grandMaPDG)/100.)%10) == kBeauty || (int(abs(grandMaPDG)/1000.)%10) == kBeauty) return kB2Jpsi;
1575       else return -1;   
1576    } 
1577    else if ( (int(abs(maPdgcode)/100.)%10) == kCharm || (int(abs(maPdgcode)/1000.)%10) == kCharm ) {
1578
1579      for (Int_t i=0; i<fNparents; i++){
1580         if (abs(maPdgcode)==fParentSelect[0][i]){
1581           isFinalOpenCharm = kTRUE;
1582         }
1583      }
1584      if (!isFinalOpenCharm) return -1;
1585
1586      // iterate until you find B hadron as a mother or become top ancester 
1587      for (Int_t i=1; i<fgkMaxIter; i++){
1588
1589         Int_t jLabel = partMother->GetFirstMother();
1590         if (jLabel == -1){
1591           return kDirectCharm;
1592         }
1593         if (jLabel < 0){ // safety protection
1594           AliDebug(1, "Stack label is negative, return\n");
1595           return -1;
1596         }
1597
1598         // if there is an ancester
1599         if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return -1; 
1600         TParticle* grandMa = mctrack->Particle();
1601         Int_t grandMaPDG = grandMa->GetPdgCode();
1602
1603         for (Int_t j=0; j<fNparents; j++){
1604            if (abs(grandMaPDG)==fParentSelect[1][j]){
1605              return kBeautyCharm;
1606            }
1607         }
1608
1609         partMother = grandMa;
1610      } // end of iteration 
1611    } // end of if
1612    else if ( (int(abs(maPdgcode)/100.)%10) == kBeauty || (int(abs(maPdgcode)/1000.)%10) == kBeauty ) {
1613      for (Int_t i=0; i<fNparents; i++){
1614         if (abs(maPdgcode)==fParentSelect[1][i]){
1615           return kDirectBeauty;
1616         }
1617      }
1618    } // end of if
1619    else if ( abs(maPdgcode) == 22 ) { //conversion
1620
1621      tmpMomLabel = partMotherCopy->GetFirstMother();
1622      if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tmpMomLabel))))) return -1;
1623      partMother = mctrack->Particle();
1624      maPdgcode = partMother->GetPdgCode();
1625      if ( abs(maPdgcode) == 111 ) {
1626        return kGammaPi0;
1627      } 
1628      else if ( abs(maPdgcode) == 221 ) {
1629        return kGammaEta;
1630      } 
1631      else if ( abs(maPdgcode) == 223 ) {
1632        return kGammaOmega;
1633      } 
1634      else if ( abs(maPdgcode) == 333 ) {
1635        return kGammaPhi;
1636      }
1637      else if ( abs(maPdgcode) == 331 ) {
1638        return kGammaEtaPrime; 
1639      }
1640      else if ( abs(maPdgcode) == 113 ) {
1641        return kGammaRho0;
1642      }
1643      else origin = kElse;
1644      return origin;
1645
1646    } 
1647    else if ( abs(maPdgcode) == 111 ) {
1648      return kPi0;
1649    } 
1650    else if ( abs(maPdgcode) == 221 ) {
1651      return kEta;
1652    } 
1653    else if ( abs(maPdgcode) == 223 ) {
1654      return kOmega;
1655    } 
1656    else if ( abs(maPdgcode) == 333 ) {
1657      return kPhi;
1658    } 
1659    else if ( abs(maPdgcode) == 331 ) {
1660      return kEtaPrime;
1661    } 
1662    else if ( abs(maPdgcode) == 113 ) {
1663      return kRho0;
1664    } 
1665    else if ( abs(maPdgcode) == 321 || abs(maPdgcode) == 130 ) {
1666      return kKe3;
1667    }
1668    else{ 
1669     origin = kElse;
1670    }
1671    return origin;
1672 }
1673 //__________________________________________
1674 Double_t AliHFEmcQA::GetWeightFactor(AliMCParticle *mctrack, const Int_t iBgLevel){
1675   //
1676   // 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)
1677   //
1678   AliMCParticle *mctrackmother = NULL;  
1679   Double_t weightElecBg = 0.;
1680   Double_t mesonPt = 0.;
1681   Double_t bgcategory = 0.;
1682   Int_t mArr = -1;  
1683   Int_t mesonID = GetElecSource(mctrack->Particle());
1684   if(mesonID==kGammaPi0 || mesonID==kPi0) mArr=0;                //pion
1685   else if(mesonID==kGammaEta || mesonID==kEta) mArr=1;           //eta
1686   else if(mesonID==kGammaOmega || mesonID==kOmega) mArr=2;       //omega
1687   else if(mesonID==kGammaPhi || mesonID==kPhi) mArr=3;           //phi
1688   else if(mesonID==kGammaEtaPrime || mesonID==kEtaPrime) mArr=4; //etaprime
1689   else if(mesonID==kGammaRho0 || mesonID==kRho0) mArr=5;         //rho
1690
1691   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};
1692   Double_t xr[3]={-999,-999,-999};
1693   datamc[0] = mesonID;
1694   datamc[17] = mctrack->Pt(); //electron pt
1695   datamc[18] = mctrack->Eta(); //electron eta
1696
1697   mctrack->XvYvZv(xr);
1698   datamc[9] = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
1699   datamc[10] = xr[2];
1700
1701   TParticle *mcpart = mctrack->Particle();
1702   if(mcpart){
1703     datamc[14] = mcpart->GetUniqueID();
1704   }
1705
1706   if(!(mArr<0)){
1707      if(mesonID>=kGammaPi0) {  // conversion electron, be careful with the enum odering 
1708         Int_t glabel=TMath::Abs(mctrack->GetMother()); // gamma label
1709         if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
1710           glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's label
1711           if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
1712             mesonPt = mctrackmother->Pt(); //meson pt
1713             bgcategory = 1.;
1714             datamc[1] = bgcategory;
1715             datamc[2] = mesonPt;
1716             mctrackmother->XvYvZv(xr);
1717             datamc[11] = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
1718             datamc[12] = xr[2];
1719
1720             mcpart = mctrackmother->Particle();
1721             if(mcpart){
1722               datamc[15] = mcpart->GetUniqueID();
1723             }
1724             if(glabel>fMCEvent->GetNumberOfPrimaries()) {
1725               bgcategory = 2.;
1726               datamc[1] = bgcategory;
1727               //printf("I should be gamma meson = %d  mesonlabel= %d  NumberOfPrimaries= %d \n",mctrackmother->PdgCode(),glabel,fMCEvent->GetNumberOfPrimaries()); 
1728               glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother
1729               if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
1730                 datamc[3]=mctrackmother->PdgCode();
1731                 datamc[4]=mctrackmother->Pt();
1732                 if(TMath::Abs(mctrackmother->PdgCode())==310){
1733                   bgcategory = 3.;
1734                   glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother's mother
1735                   if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
1736                     datamc[5]=mctrackmother->PdgCode();
1737                     glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother
1738                     if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
1739                        datamc[6]=mctrackmother->PdgCode();
1740                     }
1741                   }
1742                 }
1743               }
1744             } 
1745           } 
1746         }
1747      }
1748      else{ // nonHFE except for the conversion electron
1749         Int_t glabel=TMath::Abs(mctrack->GetMother()); 
1750         if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
1751           mesonPt=mctrackmother->Pt(); //meson pt
1752           bgcategory = -1.;
1753           datamc[1] = bgcategory;
1754           datamc[2] = mesonPt;
1755           mctrackmother->XvYvZv(xr);
1756           datamc[11] = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
1757           datamc[12] = xr[2];
1758
1759           mcpart = mctrackmother->Particle();
1760           if(mcpart){
1761             datamc[15] = mcpart->GetUniqueID();
1762           }
1763           if(glabel>fMCEvent->GetNumberOfPrimaries()) {
1764             bgcategory = -2.;
1765             datamc[1] = bgcategory;
1766             glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother
1767             if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
1768               datamc[3]=mctrackmother->PdgCode();
1769               datamc[4]=mctrackmother->Pt();
1770               if(TMath::Abs(mctrackmother->PdgCode())==310){
1771                bgcategory = -3.;
1772                glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother's mother
1773                if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
1774                  datamc[5]=mctrackmother->PdgCode();
1775                  glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother
1776                  if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
1777                    datamc[6]=mctrackmother->PdgCode();
1778                  }
1779                }
1780               }  
1781             }
1782           }
1783         }
1784      }
1785
1786
1787      if(fIsPbPb){
1788        Int_t centBin = GetWeightCentralityBin(fPerCentrality);
1789        if(centBin < 0)return 0.;
1790        weightElecBg=fElecBackgroundFactor[iBgLevel][centBin][mArr][kBgPtBins-1];                        
1791        for(int ii=0; ii<kBgPtBins; ii++){              
1792          if((mesonPt > fBinLimit[ii]) && (mesonPt < fBinLimit[ii+1])){
1793            weightElecBg = fElecBackgroundFactor[iBgLevel][centBin][mArr][ii];
1794            break;
1795          }
1796        }
1797      }
1798      else{
1799        weightElecBg=fElecBackgroundFactor[iBgLevel][0][mArr][kBgPtBins-1];                         
1800        for(int ii=0; ii<kBgPtBins; ii++){              
1801          if((mesonPt > fBinLimit[ii]) && (mesonPt < fBinLimit[ii+1])){
1802            weightElecBg = fElecBackgroundFactor[iBgLevel][0][mArr][ii];
1803            break;
1804          }
1805        }
1806      }    
1807   }
1808
1809   datamc[13] = weightElecBg;
1810   datamc[16] = Double_t(fContainerStep);
1811
1812   datamc[7] = fHfeImpactR;
1813   datamc[8] = fHfeImpactnsigmaR;
1814
1815   datamc[19] = fRecPt;
1816   datamc[20] = fRecEta;
1817   datamc[21] = fRecPhi;
1818   datamc[22] = fLyrhit;
1819   datamc[23] = fLyrstat;
1820   datamc[24] = fCentrality;
1821
1822   if(fIsDebugStreamerON && fTreeStream){
1823    if(!iBgLevel){
1824     (*fTreeStream)<<"nonhfeQA"<<
1825         "mesonID="<<datamc[0]<<
1826         "bgcategory="<<datamc[1]<<
1827         "mesonPt="<<datamc[2]<<
1828         "mesonMomPdg="<<datamc[3]<<
1829         "mesonMomPt="<<datamc[4]<<
1830         "mesonGMomPdg="<<datamc[5]<<
1831         "mesonGGMomPdg="<<datamc[6]<<
1832         "eIPAbs="<<datamc[7]<<
1833         "eIPSig="<<datamc[8]<<
1834         "eR="<<datamc[9]<<
1835         "eZ="<<datamc[10]<<
1836         "mesonR="<<datamc[11]<<
1837         "mesonZ="<<datamc[12]<<
1838         "weightElecBg="<<datamc[13]<< 
1839         "eUniqID="<<datamc[14]<<
1840         "mesonUniqID="<<datamc[15]<<
1841         "containerStep="<<datamc[16]<<
1842         "emcpt="<<datamc[17]<<
1843         "emceta="<<datamc[18]<<
1844         "erecpt="<<datamc[19]<<
1845         "ereceta="<<datamc[20]<<
1846         "erecphi="<<datamc[21]<< 
1847         "itshit="<<datamc[22]<<
1848         "itsstat="<<datamc[23]<<
1849         "centrality="<<datamc[24]
1850         << "\n";
1851    }
1852   }
1853
1854   Double_t returnval = bgcategory*weightElecBg;
1855   if(TMath::Abs(bgcategory)>1) returnval = bgcategory/2.;
1856
1857   return returnval;
1858 }
1859
1860 //__________________________________________
1861 Int_t AliHFEmcQA::GetMother(const AliVParticle * const mcpart){
1862   //
1863   // Wrapper to get the mother label
1864   //
1865   Int_t label = -1; 
1866   TClass *type = mcpart->IsA();
1867   if(type == AliMCParticle::Class()){
1868     // ESD analysis
1869     const AliMCParticle *emcpart = static_cast<const AliMCParticle *>(mcpart);
1870     label = emcpart->GetMother();
1871   } else if(type == AliAODMCParticle::Class()){
1872     // AOD analysis
1873     const AliAODMCParticle *amcpart = static_cast<const AliAODMCParticle *>(mcpart);
1874     label = amcpart->GetMother();
1875   }
1876   return label;
1877 }
1878 //__________________________________________
1879 Int_t AliHFEmcQA::GetWeightCentralityBin(const Float_t percentile) const {
1880   //
1881   //translate the centrality percentile into the centrality bin of the reference weighting histograms for electron background
1882   //
1883
1884   Float_t centralityLimits[12]= {0.,5.,10., 20., 30., 40., 50., 60.,70.,80., 90., 100.};
1885   Int_t bin = -1;
1886   for(Int_t ibin = 0; ibin < 11; ibin++){
1887     if(percentile >= centralityLimits[ibin] && percentile < centralityLimits[ibin+1]){
1888       bin = ibin;
1889       break;
1890     }
1891   } 
1892   return bin;
1893 }
1894 //__________________________________________
1895 void AliHFEmcQA::AliHists::FillList(TList *l) const {
1896   //
1897   // Fill Histos into a list for output
1898   //
1899   if(fPdgCode) l->Add(fPdgCode);
1900   if(fPt) l->Add(fPt);
1901   if(fY) l->Add(fY);
1902   if(fEta) l->Add(fEta);
1903 }
1904
1905 //__________________________________________
1906 void AliHFEmcQA::AliHistsComm::FillList(TList *l) const { 
1907   //
1908   // Fill Histos into a list for output
1909   //
1910   if(fNq) l->Add(fNq);
1911   if(fProcessID) l->Add(fProcessID);
1912   if(fePtRatio) l->Add(fePtRatio);
1913   if(fPtCorr) l->Add(fPtCorr);
1914   if(fPtCorrDp) l->Add(fPtCorrDp);
1915   if(fPtCorrD0) l->Add(fPtCorrD0);
1916   if(fPtCorrDrest) l->Add(fPtCorrDrest);
1917   if(fDePtRatio) l->Add(fDePtRatio);
1918   if(feDistance) l->Add(feDistance);
1919   if(fDeDistance) l->Add(fDeDistance);
1920 }
1921