Update hfe package
[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
373     if(icut <1){
374       hname = kqEtaRangeLabel[icut]+"PtCorrDinein_"+kqTypeLabel[kQuark];
375       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
376       hname = kqEtaRangeLabel[icut]+"PtCorrDineout_"+kqTypeLabel[kQuark];
377       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
378       hname = kqEtaRangeLabel[icut]+"PtCorrDoutein_"+kqTypeLabel[kQuark];
379       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
380       hname = kqEtaRangeLabel[icut]+"PtCorrDouteout_"+kqTypeLabel[kQuark];
381       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
382
383       hname = kqEtaRangeLabel[icut]+"PtCorrDpDinein_"+kqTypeLabel[kQuark];
384       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
385       hname = kqEtaRangeLabel[icut]+"PtCorrDpDineout_"+kqTypeLabel[kQuark];
386       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
387       hname = kqEtaRangeLabel[icut]+"PtCorrDpDoutein_"+kqTypeLabel[kQuark];
388       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
389       hname = kqEtaRangeLabel[icut]+"PtCorrDpDouteout_"+kqTypeLabel[kQuark];
390       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
391
392       hname = kqEtaRangeLabel[icut]+"PtCorrD0Dinein_"+kqTypeLabel[kQuark];
393       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
394       hname = kqEtaRangeLabel[icut]+"PtCorrD0Dineout_"+kqTypeLabel[kQuark];
395       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
396       hname = kqEtaRangeLabel[icut]+"PtCorrD0Doutein_"+kqTypeLabel[kQuark];
397       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
398       hname = kqEtaRangeLabel[icut]+"PtCorrD0Douteout_"+kqTypeLabel[kQuark];
399       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
400
401       hname = kqEtaRangeLabel[icut]+"PtCorrDrestDinein_"+kqTypeLabel[kQuark];
402       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
403       hname = kqEtaRangeLabel[icut]+"PtCorrDrestDineout_"+kqTypeLabel[kQuark];
404       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
405       hname = kqEtaRangeLabel[icut]+"PtCorrDrestDoutein_"+kqTypeLabel[kQuark];
406       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
407       hname = kqEtaRangeLabel[icut]+"PtCorrDrestDouteout_"+kqTypeLabel[kQuark];
408       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
409
410       hname = kqEtaRangeLabel[icut]+"fEtaCorrD_"+kqTypeLabel[kQuark];
411       fHistComm[iq][icut].fEtaCorrD = new TH2F(hname,hname+";D Y;e eta",200,-10,10,200,-10,10); 
412       hname = kqEtaRangeLabel[icut]+"fEtaCorrDp_"+kqTypeLabel[kQuark];
413       fHistComm[iq][icut].fEtaCorrDp = new TH2F(hname,hname+";D Y;e eta",200,-10,10,200,-10,10); 
414       hname = kqEtaRangeLabel[icut]+"fEtaCorrD0_"+kqTypeLabel[kQuark];
415       fHistComm[iq][icut].fEtaCorrD0 = new TH2F(hname,hname+";D Y;e eta",200,-10,10,200,-10,10); 
416       hname = kqEtaRangeLabel[icut]+"fEtaCorrDrest_"+kqTypeLabel[kQuark];
417       fHistComm[iq][icut].fEtaCorrDrest = new TH2F(hname,hname+";D Y;e eta",200,-10,10,200,-10,10); 
418
419       hname = kqEtaRangeLabel[icut]+"fEtaCorrGD_"+kqTypeLabel[kQuark];
420       fHistComm[iq][icut].fEtaCorrGD = new TH2F(hname,hname+";D Y;e eta",200,-10,10,200,-10,10); 
421       hname = kqEtaRangeLabel[icut]+"fEtaCorrGDp_"+kqTypeLabel[kQuark];
422       fHistComm[iq][icut].fEtaCorrGDp = new TH2F(hname,hname+";D Y;e eta",200,-10,10,200,-10,10); 
423       hname = kqEtaRangeLabel[icut]+"fEtaCorrGD0_"+kqTypeLabel[kQuark];
424       fHistComm[iq][icut].fEtaCorrGD0 = new TH2F(hname,hname+";D Y;e eta",200,-10,10,200,-10,10); 
425       hname = kqEtaRangeLabel[icut]+"fEtaCorrGDrest_"+kqTypeLabel[kQuark];
426       fHistComm[iq][icut].fEtaCorrGDrest = new TH2F(hname,hname+";D Y;e eta",200,-10,10,200,-10,10); 
427
428       hname = kqEtaRangeLabel[icut]+"fEtaCorrB_"+kqTypeLabel[kQuark];
429       fHistComm[iq][icut].fEtaCorrB = new TH2F(hname,hname+";B Y;e eta",200,-10,10,200,-10,10); 
430       hname = kqEtaRangeLabel[icut]+"fEtaCorrGB_"+kqTypeLabel[kQuark];
431       fHistComm[iq][icut].fEtaCorrGB = new TH2F(hname,hname+";B Y;e eta",200,-10,10,200,-10,10); 
432
433       hname = kqEtaRangeLabel[icut]+"PtCorrBinein_"+kqTypeLabel[kQuark];
434       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
435       hname = kqEtaRangeLabel[icut]+"PtCorrBineout_"+kqTypeLabel[kQuark];
436       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
437       hname = kqEtaRangeLabel[icut]+"PtCorrBoutein_"+kqTypeLabel[kQuark];
438       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
439       hname = kqEtaRangeLabel[icut]+"PtCorrBouteout_"+kqTypeLabel[kQuark];
440       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
441     }
442     if(fQAhistos) fHistComm[iq][icut].FillList(fQAhistos);
443   }
444
445
446   hname = kqEtaRangeLabel[0]+"Nq_"+kqTypeLabel[kQuark];
447   fHistComm[iq][0].fNq = new TH1F(hname,hname,50,-0.5,49.5);
448   hname = kqEtaRangeLabel[0]+"ProcessID_"+kqTypeLabel[kQuark];
449   fHistComm[iq][0].fProcessID = new TH1F(hname,hname,21,-10.5,10.5);
450   
451 }
452
453 //__________________________________________
454 void AliHFEmcQA::Init()
455 {
456   // called at begining every event
457   
458   for (Int_t i=0; i<2; i++){
459      fIsHeavy[i] = 0;
460   } 
461
462   fNparents = 7;
463
464   fParentSelect[0][0] =  411; //D+  
465   fParentSelect[0][1] =  421; //D0
466   fParentSelect[0][2] =  431; //Ds+
467   fParentSelect[0][3] = 4122; //Lambdac+
468   fParentSelect[0][4] = 4132; //Ksic0
469   fParentSelect[0][5] = 4232; //Ksic+
470   fParentSelect[0][6] = 4332; //OmegaC0
471
472   fParentSelect[1][0] =  511; //B0
473   fParentSelect[1][1] =  521; //B+
474   fParentSelect[1][2] =  531; //Bs0
475   fParentSelect[1][3] = 5122; //Lambdab0
476   fParentSelect[1][4] = 5132; //Ksib-
477   fParentSelect[1][5] = 5232; //Ksib0
478   fParentSelect[1][6] = 5332; //Omegab-
479
480
481 }
482
483 //__________________________________________
484 void AliHFEmcQA::GetMesonKine() 
485 {
486   //
487   // get meson pt spectra
488   //
489
490   AliVParticle *mctrack2 = NULL;
491   AliMCParticle *mctrack0 = NULL;
492   AliVParticle *mctrackdaugt= NULL;
493   AliMCParticle *mctrackd= NULL;
494   Int_t id1=0, id2=0;
495
496
497   if(fCentrality>11) printf("warning centrality out of histogram array limits \n");
498   if(!fIsPbPb&&!fIsppMultiBin) fCentrality=0;
499
500   for(Int_t imc = 0; imc <fMCEvent->GetNumberOfPrimaries(); imc++){
501      if(!(mctrack2 = fMCEvent->GetTrack(imc))) continue;
502      TParticle* mcpart0 = fMCEvent->Stack()->Particle(imc);
503      if(!mcpart0) continue;
504      mctrack0 = dynamic_cast<AliMCParticle *>(mctrack2);
505      if(!mctrack0) continue;
506
507 //     if(!fIsPbPb&&!fIsppMultiBin) fCentrality=0;
508
509      if(abs(mctrack0->PdgCode()) == 111) // pi0 
510        {
511           if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
512             fMCQACollection->Fill(Form("pionspectra_centrbin%i",fCentrality),mctrack0->Pt());
513             fMCQACollection->Fill(Form("pionspectraLog_centrbin%i",fCentrality),mctrack0->Pt());
514           }
515           id1=mctrack0->GetFirstDaughter();
516           id2=mctrack0->GetLastDaughter();
517           if(!((id2-id1)==2)) continue;
518           for(int idx=id1; idx<=id2; idx++){
519             if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
520             if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
521             if(abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
522              fMCQACollection->Fill(Form("piondaughters_centrbin%i",fCentrality),mctrackd->Pt());
523           }
524        }
525      else if(abs(mctrack0->PdgCode()) == 221) // eta 
526        {
527           if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
528             fMCQACollection->Fill(Form("etaspectra_centrbin%i",fCentrality),mctrack0->Pt());
529             fMCQACollection->Fill(Form("etaspectraLog_centrbin%i",fCentrality),mctrack0->Pt());
530           } 
531           id1=mctrack0->GetFirstDaughter();
532           id2=mctrack0->GetLastDaughter();
533           if(!((id2-id1)==2||(id2-id1)==3)) continue;
534           for(int idx=id1; idx<=id2; idx++){
535             if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
536             if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
537             if(abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
538              fMCQACollection->Fill(Form("etadaughters_centrbin%i",fCentrality),mctrackd->Pt());
539           }
540        }
541      else if(abs(mctrack0->PdgCode()) == 223) // omega
542        {
543           if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
544             fMCQACollection->Fill(Form("omegaspectra_centrbin%i",fCentrality),mctrack0->Pt());
545             fMCQACollection->Fill(Form("omegaspectraLog_centrbin%i",fCentrality),mctrack0->Pt());
546           }
547           id1=mctrack0->GetFirstDaughter();
548           id2=mctrack0->GetLastDaughter();
549           if(!((id2-id1)==1||(id2-id1)==2)) continue;
550           for(int idx=id1; idx<=id2; idx++){
551             if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
552             if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
553             if(abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
554              fMCQACollection->Fill(Form("omegadaughters_centrbin%i",fCentrality),mctrackd->Pt());
555           }
556        }
557      else if(abs(mctrack0->PdgCode()) == 333) // phi 
558        {
559           if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
560             fMCQACollection->Fill(Form("phispectra_centrbin%i",fCentrality),mctrack0->Pt());
561             fMCQACollection->Fill(Form("phispectraLog_centrbin%i",fCentrality),mctrack0->Pt());
562           } 
563           id1=mctrack0->GetFirstDaughter();
564           id2=mctrack0->GetLastDaughter();
565           if(!((id2-id1)==1)) continue;
566           for(int idx=id1; idx<=id2; idx++){
567             if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
568             if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
569             if(abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
570              fMCQACollection->Fill(Form("phidaughters_centrbin%i",fCentrality),mctrackd->Pt());
571           }
572        }
573      else if(abs(mctrack0->PdgCode()) == 331) // eta prime
574        {
575           if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
576             fMCQACollection->Fill(Form("etapspectra_centrbin%i",fCentrality),mctrack0->Pt());
577             fMCQACollection->Fill(Form("etapspectraLog_centrbin%i",fCentrality),mctrack0->Pt());
578           }
579           id1=mctrack0->GetFirstDaughter();
580           id2=mctrack0->GetLastDaughter();
581           if(!((id2-id1)==2||(id2-id1)==3)) continue;
582           for(int idx=id1; idx<=id2; idx++){
583             if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
584             if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
585             if(abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
586              fMCQACollection->Fill(Form("etapdaughters_centrbin%i",fCentrality),mctrackd->Pt());
587           }
588        }
589      else if(abs(mctrack0->PdgCode()) == 113) // rho
590        {
591           if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
592             fMCQACollection->Fill(Form("rhospectra_centrbin%i",fCentrality),mctrack0->Pt());
593             fMCQACollection->Fill(Form("rhospectraLog_centrbin%i",fCentrality),mctrack0->Pt());
594           }
595           id1=mctrack0->GetFirstDaughter();
596           id2=mctrack0->GetLastDaughter();
597           if(!((id2-id1)==1)) continue;
598           for(int idx=id1; idx<=id2; idx++){
599             if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
600             if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
601             if(abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
602              fMCQACollection->Fill(Form("rhodaughters_centrbin%i",fCentrality),mctrackd->Pt());
603           }
604        }
605      else if(abs(mctrack0->PdgCode()) == 321) // kaon+-
606        {
607           if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
608             fMCQACollection->Fill(Form("kaonspectraLog_centrbin%i",fCentrality),mctrack0->Pt());
609           }
610        }
611      else if(abs(mctrack0->PdgCode()) == 130) // k0L
612        {
613           if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
614             fMCQACollection->Fill(Form("k0LspectraLog_centrbin%i",fCentrality),mctrack0->Pt());
615           }
616        }
617      }
618
619 }
620 //__________________________________________
621 void AliHFEmcQA::GetQuarkKine(TParticle *part, Int_t iTrack, const Int_t kquark) 
622 {
623   // get heavy quark kinematics
624
625     if (kquark != kCharm && kquark != kBeauty) {
626       AliDebug(1, "This task is only for heavy quark QA, return\n");
627       return; 
628     }
629     Int_t iq = kquark - kCharm; 
630
631     if (iTrack < 0 || !part) { 
632       AliDebug(1, "Stack label is negative or no mcparticle, return\n");
633       return; 
634     }
635
636     if(!fIsPbPb&&!fIsppMultiBin) fCentrality=0;
637
638     AliMCParticle *mctrack = NULL;
639     Int_t partPdgcode = TMath::Abs(part->GetPdgCode());
640
641     // select heavy hadron or not fragmented heavy quark 
642     if ( int(partPdgcode/100.)==kquark || int(partPdgcode/1000.)==kquark || (partPdgcode==kquark && (part->GetNDaughters()==0 && iTrack>5)) ){ 
643
644       TParticle *partMother;
645       Int_t iLabel;
646
647       if (partPdgcode == kquark){ // in case of not fragmented heavy quark  
648         partMother = part; 
649         iLabel = iTrack;
650       } else{ // in case of heavy hadron, start to search for mother heavy parton 
651         iLabel = part->GetFirstMother(); 
652         if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return; 
653         if (iLabel>-1) { partMother = mctrack->Particle(); }
654         else {
655           AliDebug(1, "Stack label is negative, return\n");
656           return; 
657         }
658       }
659
660       // heavy parton selection as a mother of heavy hadron 
661       // if the heavy particle comes from string which is denoted as particle status 12|12|12...12|11,[PYTHIA p.60]
662       // in this case, the mother of heavy particle can be one of the fragmented parton of the string
663       // should I make a condition that partMother should be quark or diquark? -> not necessary
664       if ( abs(partMother->GetPdgCode()) == kquark || (partMother->GetStatusCode() == 11) ){
665       //if ( abs(partMother->GetPdgCode()) == kquark || (partMother->GetStatusCode() == 11 || partMother->GetStatusCode() == 12) ){
666
667         if ( abs(partMother->GetPdgCode()) != kquark ){
668           // search fragmented partons in the same string
669           Bool_t isSameString = kTRUE; 
670           for (Int_t i=1; i<fgkMaxIter; i++){
671              iLabel = iLabel - 1;
672              if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return; 
673              if (iLabel>-1) { partMother = mctrack->Particle(); }
674              else {
675                AliDebug(1, "Stack label is negative, return\n");
676                return; 
677              }
678              if ( abs(partMother->GetPdgCode()) == kquark ) break;
679              if ( partMother->GetStatusCode() != 12 ) isSameString = kFALSE;
680              if (!isSameString) return; 
681           }
682         }
683         AliDebug(1, "Can not find heavy parton of this heavy hadron in the string, return\n");
684         if (abs(partMother->GetPdgCode()) != kquark) return; 
685
686         if (fIsHeavy[iq] >= 50) return;  
687         fHeavyQuark[fIsHeavy[iq]] = partMother;
688         fIsHeavy[iq]++;
689
690         // fill kinematics for heavy parton
691         if (partMother->GetPdgCode() > 0) { // quark
692           fHist[iq][kQuark][0][fCentrality].fPdgCode->Fill(partMother->GetPdgCode());
693           fHist[iq][kQuark][0][fCentrality].fPt->Fill(partMother->Pt());
694           fHist[iq][kQuark][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMother));
695           fHist[iq][kQuark][0][fCentrality].fEta->Fill(partMother->Eta());
696         } else{ // antiquark
697           fHist[iq][kantiQuark][0][fCentrality].fPdgCode->Fill(partMother->GetPdgCode());
698           fHist[iq][kantiQuark][0][fCentrality].fPt->Fill(partMother->Pt());
699           fHist[iq][kantiQuark][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMother));
700           fHist[iq][kantiQuark][0][fCentrality].fEta->Fill(partMother->Eta());
701         }
702
703       } // end of heavy parton slection loop 
704
705     } // end of heavy hadron or quark selection
706
707 }
708
709 //__________________________________________
710 void AliHFEmcQA::EndOfEventAna(const Int_t kquark)
711 {
712   // end of event analysis
713
714   if (kquark != kCharm && kquark != kBeauty) {
715     AliDebug(1, "This task is only for heavy quark QA, return\n");
716     return; 
717   }
718   Int_t iq = kquark - kCharm; 
719
720
721   // # of heavy quark per event
722   AliDebug(1,Form("Number of heavy quark in this event = %d \n",fIsHeavy[iq]));
723   fHistComm[iq][0].fNq->Fill(fIsHeavy[iq]);
724
725   Int_t motherID[fgkMaxGener];
726   Int_t motherType[fgkMaxGener];
727   Int_t motherLabel[fgkMaxGener];
728   Int_t ancestorPdg[fgkMaxGener];
729   Int_t ancestorLabel[fgkMaxGener];
730
731   for (Int_t i = 0; i < fgkMaxGener; i++){ // initialization
732      motherID[i] = 0;
733      motherType[i] = 0;
734      motherLabel[i] = 0;
735      ancestorPdg[i] = 0;
736      ancestorLabel[i] = 0;
737   }
738
739
740   // check history of found heavy quarks
741   for (Int_t i = 0; i < fIsHeavy[iq]; i++){
742
743      if(!fHeavyQuark[i]) return;
744
745      ancestorLabel[0] = i;
746      ancestorPdg[0] = fHeavyQuark[i]->GetPdgCode(); 
747      ancestorLabel[1] = fHeavyQuark[i]->GetFirstMother(); 
748
749      AliDebug(1,Form("pdg code= %d\n",ancestorPdg[0]));
750      AliDebug(1,Form("ancestor label= %d\n",ancestorLabel[1]));
751
752      Int_t ig = 1;
753      while (ancestorLabel[ig] != -1){
754           // in case there is mother, get mother's pdg code and grandmother's label
755           IdentifyMother(ancestorLabel[ig], ancestorPdg[ig], ancestorLabel[ig+1]); 
756           // if mother is still heavy, find again mother's ancestor
757           if (ancestorPdg[ig-1] == ancestorPdg[ig]) {
758             ig++;
759             continue; // if it is from same heavy
760           }
761           // if the heavy's mother is not heavy, check the mother's label to know if it comes from inital or final parton shower
762           if (IsFromInitialShower(ancestorLabel[ig],motherID[i],motherType[i],motherLabel[i])) break;
763           if (IsFromFinalParton(ancestorLabel[ig],motherID[i],motherType[i],motherLabel[i])) break;
764           // if it is not the above case, something is strange
765           ReportStrangeness(motherID[i],motherType[i],motherLabel[i]);
766           break;
767      } 
768      if (ancestorLabel[ig] == -1){ // from hard scattering
769        HardScattering(kquark, motherID[i],motherType[i], motherLabel[i]);
770      }
771
772   } // end of found heavy quark loop
773
774
775   // check process type
776   Int_t processID = 0;
777   for (Int_t i = 0; i < fIsHeavy[iq]; i++){
778      AliDebug(1,Form("Mother ID= %d type= %d label= %d\n",motherID[i],motherType[i],motherLabel[i]));
779   }
780
781
782   Int_t nheavypair = Int_t(fIsHeavy[iq]/2.); 
783   for (Int_t ipair = 0; ipair < nheavypair; ipair++){
784
785      Int_t id1 = ipair*2;
786      Int_t id2 = ipair*2 + 1;
787
788      if (motherType[id1] == 2 && motherType[id2] == 2){
789        if (motherLabel[id1] == motherLabel[id2]) processID = kGluonSplitting; // gluon spliting
790        else processID = -9;
791      }
792      else if (motherType[id1] == -1 && motherType[id2] == -1) {
793        if (motherLabel[id1] == -1 && motherLabel[id2] == -1) {
794          if (motherID[id1] == fgkGluon) processID = kPairCreationFromg; // gluon fusion
795          else processID = kPairCreationFromq; // q-qbar pair creation
796        }
797        else processID = -8;
798      }
799      else if (motherType[id1] == -1 || motherType[id2] == -1) {
800        if ((motherLabel[id1] == -1 || motherLabel[id2] == -1) && (motherLabel[id1]*motherLabel[id2] == -2 || motherLabel[id1]*motherLabel[id2] == -3)) {
801          if(motherID[id1]*motherID[id2] == kquark*fgkGluon) processID = kFlavourExitation; // flavour exitation 
802          else processID = kLightQuarkShower;
803        }
804        else processID = -7;
805      }
806      else if (motherType[id1] == -2 || motherType[id2] == -2) {
807        if (motherLabel[id1] == motherLabel[id2]) processID = kInitialPartonShower; // initial parton shower
808        else processID = -6;
809        
810      }
811      else processID = -5;
812
813      if (nheavypair >1) AliDebug(1,Form("Multi pair found : process ID = %d\n",processID));
814      else fHistComm[iq][0].fProcessID->Fill(processID);
815      AliDebug(1,Form("Process ID = %d\n",processID));
816   } // end of # heavy quark pair loop
817
818 }
819
820 //__________________________________________
821 void AliHFEmcQA::GetHadronKine(TParticle* mcpart, const Int_t kquark)
822 {
823     // decay electron kinematics
824
825     if (kquark != kCharm && kquark != kBeauty) {
826       AliDebug(1, "This task is only for heavy quark QA, return\n");
827       return;
828     }
829     Int_t iq = kquark - kCharm;
830
831     if(!mcpart){
832       AliDebug(1, "no mc particle, return\n");
833       return;
834     }
835
836     Int_t iLabel = mcpart->GetFirstMother();
837     if (iLabel<0){
838       AliDebug(1, "Stack label is negative, return\n");
839       return;
840     }
841
842     if(!fIsPbPb&&!fIsppMultiBin) fCentrality=0;
843
844     TParticle *partCopy = mcpart;
845     Int_t pdgcode = mcpart->GetPdgCode();
846     Int_t pdgcodeCopy = pdgcode;
847
848     AliMCParticle *mctrack = NULL;
849
850     // if the mother is charmed hadron  
851     Bool_t isDirectCharm = kFALSE;
852     if ( int(abs(pdgcode)/100.) == kCharm || int(abs(pdgcode)/1000.) == kCharm ) {
853
854           // iterate until you find B hadron as a mother or become top ancester 
855           for (Int_t i=1; i<fgkMaxIter; i++){
856
857              Int_t jLabel = mcpart->GetFirstMother();
858              if (jLabel == -1){
859                isDirectCharm = kTRUE;
860                break; // if there is no ancester
861              }
862              if (jLabel < 0){ // safety protection
863                AliDebug(1, "Stack label is negative, return\n");
864                return;
865              }
866              // if there is an ancester
867              if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return; 
868              TParticle* mother = mctrack->Particle();
869              Int_t motherPDG = mother->GetPdgCode();
870     
871              for (Int_t j=0; j<fNparents; j++){
872                 if (abs(motherPDG)==fParentSelect[1][j]) return; // return if this hadron is originated from b
873              }
874
875              mcpart = mother;
876           } // end of iteration 
877     } // end of if
878     if((isDirectCharm == kTRUE && kquark == kCharm) || kquark == kBeauty) {
879          for (Int_t i=0; i<fNparents; i++){
880             if (abs(pdgcodeCopy)==fParentSelect[iq][i]){
881
882               // fill hadron kinematics
883               fHist[iq][kHadron][0][fCentrality].fPdgCode->Fill(pdgcodeCopy);
884               fHist[iq][kHadron][0][fCentrality].fPt->Fill(partCopy->Pt());
885               fHist[iq][kHadron][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partCopy));
886               fHist[iq][kHadron][0][fCentrality].fEta->Fill(partCopy->Eta());
887
888               if(iq==0) {
889                fhD[i]->Fill(partCopy->Pt(),AliHFEtools::GetRapidity(partCopy));
890               }
891             }
892          }
893          // I also want to store D* info to compare with D* measurement 
894          if (abs(pdgcodeCopy)==413 && iq==0) { //D*+
895                fhD[7]->Fill(partCopy->Pt(),AliHFEtools::GetRapidity(partCopy));
896          }
897          if (abs(pdgcodeCopy)==423 && iq==0) { //D*0
898                fhD[8]->Fill(partCopy->Pt(),AliHFEtools::GetRapidity(partCopy));
899          }
900     } // end of if
901 }
902
903 //__________________________________________
904 void AliHFEmcQA::GetDecayedKine(TParticle* mcpart, const Int_t kquark, Int_t kdecayed) 
905 {
906     // decay electron kinematics
907     
908     if (!(kquark == kCharm || kquark == kBeauty || kquark == kOthers)){
909       AliDebug(1, "This task is only for heavy quark QA, return\n");
910       return; 
911     }
912     Int_t iq = kquark - kCharm; 
913     Bool_t isFinalOpenCharm = kFALSE;
914
915     if(!mcpart){
916       AliDebug(1, "no mcparticle, return\n");
917       return;
918     }
919
920     if(!fIsPbPb&&!fIsppMultiBin) fCentrality=0;
921
922     Double_t eabsEta = TMath::Abs(mcpart->Eta());
923     Double_t eabsY = TMath::Abs(AliHFEtools::GetRapidity(mcpart));
924
925     if(kquark==kOthers){
926       Int_t esource = -1;
927       if ( abs(mcpart->GetPdgCode()) != kdecayed ) esource = kMisID-4;
928       else esource =GetSource(mcpart)-4; // return for the cases kGamma=4, kPi0=5, kElse=6
929       if(esource==0|| esource==1 || esource==2 || esource==3){
930         fHist[iq][esource][0][fCentrality].fPt->Fill(mcpart->Pt());
931         fHist[iq][esource][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
932         fHist[iq][esource][0][fCentrality].fEta->Fill(mcpart->Eta());
933         if(eabsEta<0.9){
934           fHist[iq][esource][1][fCentrality].fPt->Fill(mcpart->Pt());
935           fHist[iq][esource][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
936           fHist[iq][esource][1][fCentrality].fEta->Fill(mcpart->Eta());
937         }
938         if(eabsY<0.5){
939           fHist[iq][esource][2][fCentrality].fPt->Fill(mcpart->Pt());
940           fHist[iq][esource][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
941           fHist[iq][esource][2][fCentrality].fEta->Fill(mcpart->Eta());
942         }
943         return; 
944       }
945       else {
946         AliDebug(1, "e source is out of defined ranges, return\n");
947         return;
948       }
949     }
950
951     if ( abs(mcpart->GetPdgCode()) != kdecayed ) return;
952
953     Int_t iLabel = mcpart->GetFirstMother(); 
954     if (iLabel<0){
955       AliDebug(1, "Stack label is negative, return\n");
956       return; 
957     }
958
959     AliMCParticle *mctrack = NULL;
960     if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return; 
961     TParticle *partMother = mctrack->Particle();
962     TParticle *partMotherCopy = partMother;
963     Int_t maPdgcode = partMother->GetPdgCode();
964     Int_t maPdgcodeCopy = maPdgcode;
965
966     // get mc primary vertex
967     /*
968     TArrayF mcPrimVtx;
969     if(fMCHeader) fMCHeader->PrimaryVertex(mcPrimVtx);
970
971     // get electron production vertex   
972     TLorentzVector ePoint;
973     mcpart->ProductionVertex(ePoint);
974
975     // calculated production vertex to primary vertex (in xy plane)
976     Float_t decayLxy = TMath::Sqrt((mcPrimVtx[0]-ePoint[0])*(mcPrimVtx[0]-ePoint[0])+(mcPrimVtx[1]-ePoint[1])*(mcPrimVtx[1]-ePoint[1]));
977     */ 
978     Float_t decayLxy = 0;
979
980     // if the mother is charmed hadron  
981     Bool_t isMotherDirectCharm = kFALSE;
982     if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) { 
983
984          for (Int_t i=0; i<fNparents; i++){
985             if (abs(maPdgcode)==fParentSelect[0][i]){
986               isFinalOpenCharm = kTRUE;
987             } 
988          }  
989          if (!isFinalOpenCharm) return ;
990
991           // iterate until you find B hadron as a mother or become top ancester 
992           for (Int_t i=1; i<fgkMaxIter; i++){
993
994              Int_t jLabel = partMother->GetFirstMother(); 
995              if (jLabel == -1){
996                isMotherDirectCharm = kTRUE;
997                break; // if there is no ancester
998              }
999              if (jLabel < 0){ // safety protection
1000                AliDebug(1, "Stack label is negative, return\n");
1001                return; 
1002              }
1003
1004              // if there is an ancester
1005              if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return; 
1006              TParticle* grandMa = mctrack->Particle();
1007              Int_t grandMaPDG = grandMa->GetPdgCode();
1008
1009              for (Int_t j=0; j<fNparents; j++){
1010                 if (abs(grandMaPDG)==fParentSelect[1][j]){
1011
1012                   if (kquark == kCharm) return;
1013                   // fill electron kinematics
1014                   fHist[iq][kElectron2nd][0][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1015                   fHist[iq][kElectron2nd][0][fCentrality].fPt->Fill(mcpart->Pt());
1016                   fHist[iq][kElectron2nd][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1017                   fHist[iq][kElectron2nd][0][fCentrality].fEta->Fill(mcpart->Eta());
1018
1019                   // fill mother hadron kinematics
1020                   fHist[iq][kDeHadron][0][fCentrality].fPdgCode->Fill(grandMaPDG); 
1021                   fHist[iq][kDeHadron][0][fCentrality].fPt->Fill(grandMa->Pt());
1022                   fHist[iq][kDeHadron][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(grandMa));
1023                   fHist[iq][kDeHadron][0][fCentrality].fEta->Fill(grandMa->Eta());
1024
1025                   if(eabsEta<0.9){
1026                     fHist[iq][kElectron2nd][1][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1027                     fHist[iq][kElectron2nd][1][fCentrality].fPt->Fill(mcpart->Pt());
1028                     fHist[iq][kElectron2nd][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1029                     fHist[iq][kElectron2nd][1][fCentrality].fEta->Fill(mcpart->Eta());
1030
1031                     // fill mother hadron kinematics
1032                     fHist[iq][kDeHadron][1][fCentrality].fPdgCode->Fill(grandMaPDG); 
1033                     fHist[iq][kDeHadron][1][fCentrality].fPt->Fill(grandMa->Pt());
1034                     fHist[iq][kDeHadron][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(grandMa));
1035                     fHist[iq][kDeHadron][1][fCentrality].fEta->Fill(grandMa->Eta());
1036                   }
1037
1038                   if(eabsY<0.5){
1039                     fHist[iq][kElectron2nd][2][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1040                     fHist[iq][kElectron2nd][2][fCentrality].fPt->Fill(mcpart->Pt());
1041                     fHist[iq][kElectron2nd][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1042                     fHist[iq][kElectron2nd][2][fCentrality].fEta->Fill(mcpart->Eta());
1043
1044                     // fill mother hadron kinematics
1045                     fHist[iq][kDeHadron][2][fCentrality].fPdgCode->Fill(grandMaPDG); 
1046                     fHist[iq][kDeHadron][2][fCentrality].fPt->Fill(grandMa->Pt());
1047                     fHist[iq][kDeHadron][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(grandMa));
1048                     fHist[iq][kDeHadron][2][fCentrality].fEta->Fill(grandMa->Eta());
1049                   }
1050
1051                   //mj: to calculate B to e eta correlation to calculate total heavy quark cross section
1052                   Int_t kLabel0 = grandMa->GetFirstMother();
1053                   Bool_t isGGrandmaYes = kFALSE;
1054                   Double_t ggmrapidwstmp=0;
1055                   if (!(kLabel0 < 0)){ // safety protection
1056                     if((mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(kLabel0))))){
1057                       TParticle* ggrandMatmp = mctrack->Particle();
1058                       Int_t ggrandMaPDGtmp = ggrandMatmp->GetPdgCode();
1059                       if ( int(abs(ggrandMaPDGtmp)/100.) == kBeauty || int(abs(ggrandMaPDGtmp)/1000.) == kBeauty) isGGrandmaYes = kTRUE;
1060                       ggmrapidwstmp = AliHFEtools::GetRapidity(ggrandMatmp);
1061                     }
1062                   }
1063
1064                   Double_t gmrapidwstmp0 = AliHFEtools::GetRapidity(grandMa);
1065                   Double_t eetawstmp0 = mcpart->Eta();
1066   
1067                   Double_t gmrapidtmp0 = TMath::Abs(gmrapidwstmp0);
1068                   Double_t eetatmp0 = TMath::Abs(eetawstmp0);
1069
1070                   fHistComm[iq][0].fEtaCorrB->Fill(gmrapidwstmp0,eetawstmp0);
1071                   if(isGGrandmaYes) fHistComm[iq][0].fEtaCorrGB->Fill(ggmrapidwstmp,eetawstmp0);
1072                   else fHistComm[iq][0].fEtaCorrGB->Fill(gmrapidwstmp0,eetawstmp0);
1073
1074                   if(gmrapidtmp0<0.5 && eetatmp0<0.5 ) fHistComm[iq][0].fPtCorrBinein->Fill(grandMa->Pt(),mcpart->Pt());
1075                   else if(gmrapidtmp0<0.5 && eetatmp0>0.5 ) fHistComm[iq][0].fPtCorrBineout->Fill(grandMa->Pt(),mcpart->Pt());
1076                   else if(gmrapidtmp0>0.5 && eetatmp0<0.5 ) fHistComm[iq][0].fPtCorrBoutein->Fill(grandMa->Pt(),mcpart->Pt());
1077                   else if(gmrapidtmp0>0.5 && eetatmp0>0.5 ) fHistComm[iq][0].fPtCorrBouteout->Fill(grandMa->Pt(),mcpart->Pt());
1078                   //======================================================================================
1079
1080                   // ratio between pT of electron and pT of mother B hadron 
1081                   if(grandMa->Pt()) {
1082                     fHistComm[iq][0].fDePtRatio->Fill(grandMa->Pt(),mcpart->Pt()/grandMa->Pt());
1083                     if(eabsEta<0.9){
1084                       fHistComm[iq][1].fDePtRatio->Fill(grandMa->Pt(),mcpart->Pt()/grandMa->Pt());
1085                     }
1086                     if(eabsY<0.5){
1087                       fHistComm[iq][2].fDePtRatio->Fill(grandMa->Pt(),mcpart->Pt()/grandMa->Pt());
1088                     }
1089                   }
1090
1091                   // distance between electron production point and primary vertex
1092                   fHistComm[iq][0].fDeDistance->Fill(grandMa->Pt(),decayLxy);
1093                   if(eabsEta<0.9){
1094                     fHistComm[iq][1].fDeDistance->Fill(grandMa->Pt(),decayLxy);
1095                   }
1096                   if(eabsY<0.5){
1097                     fHistComm[iq][2].fDeDistance->Fill(grandMa->Pt(),decayLxy);
1098                   }
1099                   return;
1100                 }
1101              } 
1102
1103              partMother = grandMa;
1104           } // end of iteration 
1105     } // end of if
1106     if((isMotherDirectCharm == kTRUE && kquark == kCharm) || kquark == kBeauty) {
1107          for (Int_t i=0; i<fNparents; i++){
1108             if (abs(maPdgcodeCopy)==fParentSelect[iq][i]){
1109
1110               fHist[iq][kElectron][0][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1111               fHist[iq][kElectron][0][fCentrality].fPt->Fill(mcpart->Pt());
1112               fHist[iq][kElectron][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1113               fHist[iq][kElectron][0][fCentrality].fEta->Fill(mcpart->Eta());  
1114
1115               // fill mother hadron kinematics
1116               fHist[iq][keHadron][0][fCentrality].fPdgCode->Fill(maPdgcodeCopy); 
1117               fHist[iq][keHadron][0][fCentrality].fPt->Fill(partMotherCopy->Pt());
1118               fHist[iq][keHadron][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
1119               fHist[iq][keHadron][0][fCentrality].fEta->Fill(partMotherCopy->Eta());
1120
1121               if(eabsEta<0.9){
1122                 fHist[iq][kElectron][1][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1123                 fHist[iq][kElectron][1][fCentrality].fPt->Fill(mcpart->Pt());
1124                 fHist[iq][kElectron][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1125                 fHist[iq][kElectron][1][fCentrality].fEta->Fill(mcpart->Eta());  
1126
1127                 // fill mother hadron kinematics
1128                 fHist[iq][keHadron][1][fCentrality].fPdgCode->Fill(maPdgcodeCopy); 
1129                 fHist[iq][keHadron][1][fCentrality].fPt->Fill(partMotherCopy->Pt());
1130                 fHist[iq][keHadron][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
1131                 fHist[iq][keHadron][1][fCentrality].fEta->Fill(partMotherCopy->Eta());
1132               }
1133
1134               if(eabsY<0.5){
1135                 fHist[iq][kElectron][2][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1136                 fHist[iq][kElectron][2][fCentrality].fPt->Fill(mcpart->Pt());
1137                 fHist[iq][kElectron][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1138                 fHist[iq][kElectron][2][fCentrality].fEta->Fill(mcpart->Eta());  
1139
1140                 // fill mother hadron kinematics
1141                 fHist[iq][keHadron][2][fCentrality].fPdgCode->Fill(maPdgcodeCopy); 
1142                 fHist[iq][keHadron][2][fCentrality].fPt->Fill(partMotherCopy->Pt());
1143                 fHist[iq][keHadron][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
1144                 fHist[iq][keHadron][2][fCentrality].fEta->Fill(partMotherCopy->Eta());
1145               }
1146
1147               // ratio between pT of electron and pT of mother B or direct D hadron 
1148               if(partMotherCopy->Pt()) {
1149                  fHistComm[iq][0].fePtRatio->Fill(partMotherCopy->Pt(),mcpart->Pt()/partMotherCopy->Pt());
1150                  if(eabsEta<0.9){
1151                    fHistComm[iq][1].fePtRatio->Fill(partMotherCopy->Pt(),mcpart->Pt()/partMotherCopy->Pt());
1152                  }
1153                  if(eabsY<0.5){
1154                    fHistComm[iq][2].fePtRatio->Fill(partMotherCopy->Pt(),mcpart->Pt()/partMotherCopy->Pt());
1155                  }
1156               }
1157               fHistComm[iq][0].fPtCorr->Fill(partMotherCopy->Pt(),mcpart->Pt());
1158               if(eabsEta<0.9){
1159                 fHistComm[iq][1].fPtCorr->Fill(partMotherCopy->Pt(),mcpart->Pt());
1160               }
1161               if(eabsY<0.5){
1162                 fHistComm[iq][2].fPtCorr->Fill(partMotherCopy->Pt(),mcpart->Pt());
1163               }
1164               if(TMath::Abs(partMotherCopy->GetPdgCode())==411) {
1165                 fHistComm[iq][0].fPtCorrDp->Fill(partMotherCopy->Pt(),mcpart->Pt());
1166                 if(eabsEta<0.9){
1167                   fHistComm[iq][1].fPtCorrDp->Fill(partMotherCopy->Pt(),mcpart->Pt());
1168                 }
1169                 if(eabsY<0.5){
1170                   fHistComm[iq][2].fPtCorrDp->Fill(partMotherCopy->Pt(),mcpart->Pt());
1171                 }
1172               }
1173               else if(TMath::Abs(partMotherCopy->GetPdgCode())==421) {
1174                 fHistComm[iq][0].fPtCorrD0->Fill(partMotherCopy->Pt(),mcpart->Pt());
1175                 if(eabsEta<0.9){
1176                   fHistComm[iq][1].fPtCorrD0->Fill(partMotherCopy->Pt(),mcpart->Pt());
1177                 }
1178                 if(eabsY<0.5){
1179                   fHistComm[iq][2].fPtCorrD0->Fill(partMotherCopy->Pt(),mcpart->Pt());
1180                 }
1181               }
1182               else {
1183                 fHistComm[iq][0].fPtCorrDrest->Fill(partMotherCopy->Pt(),mcpart->Pt());
1184                 if(eabsEta<0.9){
1185                   fHistComm[iq][1].fPtCorrDrest->Fill(partMotherCopy->Pt(),mcpart->Pt());
1186                 }
1187                 if(eabsY<0.5){
1188                   fHistComm[iq][2].fPtCorrDrest->Fill(partMotherCopy->Pt(),mcpart->Pt());
1189                 }
1190               }
1191
1192               //mj: to calculate D to e eta correlation to calculate total heavy quark cross section
1193               Int_t kLabel = partMotherCopy->GetFirstMother();
1194               Bool_t isGrandmaYes = kFALSE;
1195               Double_t gmrapidwstmp =0;
1196               if (!(kLabel < 0)){ // safety protection
1197                 if((mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(kLabel))))){ 
1198                   TParticle* grandMatmp = mctrack->Particle();
1199                   Int_t grandMaPDGtmp = grandMatmp->GetPdgCode();
1200                   if ( int(abs(grandMaPDGtmp)/100.) == kCharm || int(abs(grandMaPDGtmp)/1000.) == kCharm ) isGrandmaYes = kTRUE;
1201                   if ( int(abs(grandMaPDGtmp)/100.) == kBeauty || int(abs(grandMaPDGtmp)/1000.) == kBeauty) isGrandmaYes = kTRUE;
1202                   gmrapidwstmp = AliHFEtools::GetRapidity(grandMatmp);
1203                 }
1204               }
1205
1206               Double_t mrapidwstmp = AliHFEtools::GetRapidity(partMotherCopy);
1207               Double_t eetawstmp = mcpart->Eta();
1208
1209               Double_t mrapidtmp = TMath::Abs(mrapidwstmp);
1210               Double_t eetatmp = TMath::Abs(eetawstmp);
1211
1212               fHistComm[iq][0].fEtaCorrD->Fill(mrapidwstmp,eetawstmp);
1213               if(isGrandmaYes) fHistComm[iq][0].fEtaCorrGD->Fill(gmrapidwstmp,eetawstmp);
1214               else fHistComm[iq][0].fEtaCorrGD->Fill(mrapidwstmp,eetawstmp);
1215
1216               if(mrapidtmp<0.5 && eetatmp<0.5 ) fHistComm[iq][0].fPtCorrDinein->Fill(partMotherCopy->Pt(),mcpart->Pt());
1217               else if(mrapidtmp<0.5 && eetatmp>0.5 ) fHistComm[iq][0].fPtCorrDineout->Fill(partMotherCopy->Pt(),mcpart->Pt());
1218               else if(mrapidtmp>0.5 && eetatmp<0.5 ) fHistComm[iq][0].fPtCorrDoutein->Fill(partMotherCopy->Pt(),mcpart->Pt());
1219               else if(mrapidtmp>0.5 && eetatmp>0.5 ) fHistComm[iq][0].fPtCorrDouteout->Fill(partMotherCopy->Pt(),mcpart->Pt());
1220               if(TMath::Abs(partMotherCopy->GetPdgCode())==411) {
1221                 fHistComm[iq][0].fEtaCorrDp->Fill(mrapidwstmp,eetawstmp);
1222                 if(isGrandmaYes) fHistComm[iq][0].fEtaCorrGDp->Fill(gmrapidwstmp,eetawstmp);
1223                 else fHistComm[iq][0].fEtaCorrGDp->Fill(mrapidwstmp,eetawstmp);
1224                 if(mrapidtmp<0.5 && eetatmp<0.5 ) fHistComm[iq][0].fPtCorrDpDinein->Fill(partMotherCopy->Pt(),mcpart->Pt());
1225                 else if(mrapidtmp<0.5 && eetatmp>0.5 ) fHistComm[iq][0].fPtCorrDpDineout->Fill(partMotherCopy->Pt(),mcpart->Pt());
1226                 else if(mrapidtmp>0.5 && eetatmp<0.5 ) fHistComm[iq][0].fPtCorrDpDoutein->Fill(partMotherCopy->Pt(),mcpart->Pt());
1227                 else if(mrapidtmp>0.5 && eetatmp>0.5 ) fHistComm[iq][0].fPtCorrDpDouteout->Fill(partMotherCopy->Pt(),mcpart->Pt());
1228               }
1229               else if(TMath::Abs(partMotherCopy->GetPdgCode())==421) {
1230                 fHistComm[iq][0].fEtaCorrD0->Fill(mrapidwstmp,eetawstmp);
1231                 if(isGrandmaYes) fHistComm[iq][0].fEtaCorrGD0->Fill(gmrapidwstmp,eetawstmp);
1232                 else fHistComm[iq][0].fEtaCorrGD0->Fill(mrapidwstmp,eetawstmp);
1233                 if(mrapidtmp<0.5 && eetatmp<0.5 ) fHistComm[iq][0].fPtCorrD0Dinein->Fill(partMotherCopy->Pt(),mcpart->Pt());
1234                 else if(mrapidtmp<0.5 && eetatmp>0.5 ) fHistComm[iq][0].fPtCorrD0Dineout->Fill(partMotherCopy->Pt(),mcpart->Pt());
1235                 else if(mrapidtmp>0.5 && eetatmp<0.5 ) fHistComm[iq][0].fPtCorrD0Doutein->Fill(partMotherCopy->Pt(),mcpart->Pt());
1236                 else if(mrapidtmp>0.5 && eetatmp>0.5 ) fHistComm[iq][0].fPtCorrD0Douteout->Fill(partMotherCopy->Pt(),mcpart->Pt());
1237               }
1238               else {
1239                 fHistComm[iq][0].fEtaCorrDrest->Fill(mrapidwstmp,eetawstmp);
1240                 if(isGrandmaYes) fHistComm[iq][0].fEtaCorrGDrest->Fill(gmrapidwstmp,eetawstmp);
1241                 else fHistComm[iq][0].fEtaCorrGDrest->Fill(mrapidwstmp,eetawstmp);
1242                 if(mrapidtmp<0.5 && eetatmp<0.5 ) fHistComm[iq][0].fPtCorrDrestDinein->Fill(partMotherCopy->Pt(),mcpart->Pt());
1243                 else if(mrapidtmp<0.5 && eetatmp>0.5 ) fHistComm[iq][0].fPtCorrDrestDineout->Fill(partMotherCopy->Pt(),mcpart->Pt());
1244                 else if(mrapidtmp>0.5 && eetatmp<0.5 ) fHistComm[iq][0].fPtCorrDrestDoutein->Fill(partMotherCopy->Pt(),mcpart->Pt());
1245                 else if(mrapidtmp>0.5 && eetatmp>0.5 ) fHistComm[iq][0].fPtCorrDrestDouteout->Fill(partMotherCopy->Pt(),mcpart->Pt());
1246               }
1247
1248               // distance between electron production point and primary vertex
1249               fHistComm[iq][0].feDistance->Fill(partMotherCopy->Pt(),decayLxy);
1250               if(eabsEta<0.9){
1251                 fHistComm[iq][1].feDistance->Fill(partMotherCopy->Pt(),decayLxy);
1252               }
1253               if(eabsY<0.5){
1254                 fHistComm[iq][2].feDistance->Fill(partMotherCopy->Pt(),decayLxy);
1255               }
1256             }
1257          }
1258     } // end of if
1259 }
1260
1261 //____________________________________________________________________
1262 void  AliHFEmcQA::GetDecayedKine(AliAODMCParticle *mcpart, const Int_t kquark, Int_t kdecayed)
1263 {
1264   // decay electron kinematics
1265
1266   if (kquark != kCharm && kquark != kBeauty) {
1267     AliDebug(1, "This task is only for heavy quark QA, return\n");
1268     return;
1269   }
1270
1271   Int_t iq = kquark - kCharm;
1272   Bool_t isFinalOpenCharm = kFALSE;
1273
1274   if(!mcpart){
1275     AliDebug(1, "no mcparticle, return\n");
1276     return;
1277   }
1278
1279   if ( abs(mcpart->GetPdgCode()) != kdecayed ) return;
1280
1281   Double_t eabsEta = TMath::Abs(mcpart->Eta());
1282   Double_t eabsY = TMath::Abs(AliHFEtools::GetRapidity(mcpart));
1283
1284   // mother
1285   Int_t iLabel = mcpart->GetMother();
1286   if (iLabel<0){
1287     AliDebug(1, "Stack label is negative, return\n");
1288     return;
1289   }
1290
1291   if(!fIsPbPb&&!fIsppMultiBin) fCentrality=0;
1292
1293   AliAODMCParticle *partMother = (AliAODMCParticle*)fMCArray->At(iLabel);
1294   AliAODMCParticle *partMotherCopy = partMother;
1295   Int_t maPdgcode = partMother->GetPdgCode();
1296   Int_t maPdgcodeCopy = maPdgcode;
1297
1298   Bool_t isMotherDirectCharm = kFALSE;
1299   if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
1300
1301     for (Int_t i=0; i<fNparents; i++){
1302        if (abs(maPdgcode)==fParentSelect[0][i]){
1303          isFinalOpenCharm = kTRUE;
1304        }
1305     } 
1306     if (!isFinalOpenCharm) return;
1307
1308     for (Int_t i=1; i<fgkMaxIter; i++){
1309
1310        Int_t jLabel = partMother->GetMother();
1311        if (jLabel == -1){
1312          isMotherDirectCharm = kTRUE;
1313          break; // if there is no ancester
1314        }
1315        if (jLabel < 0){ // safety protection
1316          AliDebug(1, "Stack label is negative, return\n");
1317          return;
1318        }
1319
1320        // if there is an ancester
1321        AliAODMCParticle* grandMa = (AliAODMCParticle*)fMCArray->At(jLabel);
1322        Int_t grandMaPDG = grandMa->GetPdgCode();
1323
1324        for (Int_t j=0; j<fNparents; j++){
1325           if (abs(grandMaPDG)==fParentSelect[1][j]){
1326
1327             if (kquark == kCharm) return;
1328             // fill electron kinematics
1329             fHist[iq][kElectron2nd][0][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1330             fHist[iq][kElectron2nd][0][fCentrality].fPt->Fill(mcpart->Pt());
1331             fHist[iq][kElectron2nd][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1332             fHist[iq][kElectron2nd][0][fCentrality].fEta->Fill(mcpart->Eta());
1333
1334             // fill mother hadron kinematics
1335             fHist[iq][kDeHadron][0][fCentrality].fPdgCode->Fill(grandMaPDG);
1336             fHist[iq][kDeHadron][0][fCentrality].fPt->Fill(grandMa->Pt());
1337             fHist[iq][kDeHadron][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(grandMa));
1338             fHist[iq][kDeHadron][0][fCentrality].fEta->Fill(grandMa->Eta());
1339
1340             if(eabsEta<0.9){
1341               // fill electron kinematics
1342               fHist[iq][kElectron2nd][1][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1343               fHist[iq][kElectron2nd][1][fCentrality].fPt->Fill(mcpart->Pt());
1344               fHist[iq][kElectron2nd][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1345               fHist[iq][kElectron2nd][1][fCentrality].fEta->Fill(mcpart->Eta());
1346
1347               // fill mother hadron kinematics
1348               fHist[iq][kDeHadron][1][fCentrality].fPdgCode->Fill(grandMaPDG);
1349               fHist[iq][kDeHadron][1][fCentrality].fPt->Fill(grandMa->Pt());
1350               fHist[iq][kDeHadron][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(grandMa));
1351               fHist[iq][kDeHadron][1][fCentrality].fEta->Fill(grandMa->Eta());
1352             }
1353             if(eabsY<0.5){
1354               // fill electron kinematics
1355               fHist[iq][kElectron2nd][2][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1356               fHist[iq][kElectron2nd][2][fCentrality].fPt->Fill(mcpart->Pt());
1357               fHist[iq][kElectron2nd][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1358               fHist[iq][kElectron2nd][2][fCentrality].fEta->Fill(mcpart->Eta());
1359
1360               // fill mother hadron kinematics
1361               fHist[iq][kDeHadron][2][fCentrality].fPdgCode->Fill(grandMaPDG);
1362               fHist[iq][kDeHadron][2][fCentrality].fPt->Fill(grandMa->Pt());
1363               fHist[iq][kDeHadron][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(grandMa));
1364               fHist[iq][kDeHadron][2][fCentrality].fEta->Fill(grandMa->Eta());
1365             }
1366
1367             return;
1368           }
1369        }
1370
1371        partMother = grandMa;
1372     } // end of iteration 
1373   } // end of if
1374   if ((isMotherDirectCharm == kTRUE && kquark == kCharm) || kquark == kBeauty) {
1375     for (Int_t i=0; i<fNparents; i++){
1376        if (abs(maPdgcodeCopy)==fParentSelect[iq][i]){
1377
1378          // fill electron kinematics
1379          fHist[iq][kElectron][0][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1380          fHist[iq][kElectron][0][fCentrality].fPt->Fill(mcpart->Pt());
1381          fHist[iq][kElectron][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1382          fHist[iq][kElectron][0][fCentrality].fEta->Fill(mcpart->Eta());
1383
1384          // fill mother hadron kinematics
1385          fHist[iq][keHadron][0][fCentrality].fPdgCode->Fill(maPdgcodeCopy);
1386          fHist[iq][keHadron][0][fCentrality].fPt->Fill(partMotherCopy->Pt());
1387          fHist[iq][keHadron][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
1388          fHist[iq][keHadron][0][fCentrality].fEta->Fill(partMotherCopy->Eta());
1389
1390          if(eabsEta<0.9){
1391            // fill electron kinematics
1392            fHist[iq][kElectron][1][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1393            fHist[iq][kElectron][1][fCentrality].fPt->Fill(mcpart->Pt());
1394            fHist[iq][kElectron][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1395            fHist[iq][kElectron][1][fCentrality].fEta->Fill(mcpart->Eta());
1396
1397            // fill mother hadron kinematics
1398            fHist[iq][keHadron][1][fCentrality].fPdgCode->Fill(maPdgcodeCopy);
1399            fHist[iq][keHadron][1][fCentrality].fPt->Fill(partMotherCopy->Pt());
1400            fHist[iq][keHadron][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
1401            fHist[iq][keHadron][1][fCentrality].fEta->Fill(partMotherCopy->Eta());
1402          }
1403          if(eabsY<0.5){
1404            // fill electron kinematics
1405            fHist[iq][kElectron][2][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1406            fHist[iq][kElectron][2][fCentrality].fPt->Fill(mcpart->Pt());
1407            fHist[iq][kElectron][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1408            fHist[iq][kElectron][2][fCentrality].fEta->Fill(mcpart->Eta());
1409
1410            // fill mother hadron kinematics
1411            fHist[iq][keHadron][2][fCentrality].fPdgCode->Fill(maPdgcodeCopy);
1412            fHist[iq][keHadron][2][fCentrality].fPt->Fill(partMotherCopy->Pt());
1413            fHist[iq][keHadron][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
1414            fHist[iq][keHadron][2][fCentrality].fEta->Fill(partMotherCopy->Eta());
1415          }
1416
1417        }
1418     }
1419   } // end of if
1420
1421 }
1422
1423 //__________________________________________
1424 void AliHFEmcQA::IdentifyMother(Int_t motherlabel, Int_t &motherpdg, Int_t &grandmotherlabel)
1425 {
1426        // find mother pdg code and label 
1427
1428        if (motherlabel < 0) { 
1429          AliDebug(1, "Stack label is negative, return\n");
1430          return; 
1431        }
1432        AliMCParticle *mctrack = NULL;
1433        if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(motherlabel))))) return; 
1434        TParticle *heavysMother = mctrack->Particle();
1435        motherpdg = heavysMother->GetPdgCode();
1436        grandmotherlabel = heavysMother->GetFirstMother();
1437        AliDebug(1,Form("ancestor pdg code= %d\n",motherpdg));
1438 }
1439
1440 //__________________________________________
1441 void AliHFEmcQA::HardScattering(const Int_t kquark, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
1442 {
1443        // mothertype -1 means this heavy quark coming from hard vertex
1444
1445        AliMCParticle *mctrack1 = NULL;
1446        AliMCParticle *mctrack2 = NULL;
1447        if(!(mctrack1 = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(4))))) return; 
1448        if(!(mctrack2 = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(5))))) return; 
1449        TParticle *afterinitialrad1  = mctrack1->Particle();
1450        TParticle *afterinitialrad2  = mctrack2->Particle();
1451            
1452        motherlabel = -1;
1453
1454        if (abs(afterinitialrad1->GetPdgCode()) == fgkGluon && abs(afterinitialrad2->GetPdgCode()) == fgkGluon){
1455          AliDebug(1,"heavy from gluon gluon pair creation!\n");
1456          mothertype = -1;
1457          motherID = fgkGluon;
1458        }
1459        else if (abs(afterinitialrad1->GetPdgCode()) == kquark || abs(afterinitialrad2->GetPdgCode()) == kquark){ // one from Q and the other from g
1460          AliDebug(1,"heavy from flavor exitation!\n");
1461          mothertype = -1;
1462          motherID = kquark;
1463        }
1464        else if  (abs(afterinitialrad1->GetPdgCode()) == abs(afterinitialrad2->GetPdgCode())){
1465          AliDebug(1,"heavy from q-qbar pair creation!\n");
1466          mothertype = -1;
1467          motherID = 1;
1468        }
1469        else {
1470          AliDebug(1,"something strange!\n");
1471          mothertype = -999;
1472          motherlabel = -999;
1473          motherID = -999;
1474        }
1475 }
1476
1477 //__________________________________________
1478 Bool_t AliHFEmcQA::IsFromInitialShower(Int_t inputmotherlabel, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
1479 {
1480        // mothertype -2 means this heavy quark coming from initial state 
1481
1482        AliMCParticle *mctrack = NULL;
1483        if (inputmotherlabel==2 || inputmotherlabel==3){ // mother exist before initial state radiation
1484          if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(inputmotherlabel))))) return kFALSE; 
1485          TParticle *heavysMother = mctrack->Particle();
1486          motherID = heavysMother->GetPdgCode(); 
1487          mothertype = -2; // there is mother before initial state radiation
1488          motherlabel = inputmotherlabel;
1489          AliDebug(1,"initial parton shower! \n");
1490
1491          return kTRUE;
1492        }
1493
1494        return kFALSE;
1495 }
1496
1497 //__________________________________________
1498 Bool_t AliHFEmcQA::IsFromFinalParton(Int_t inputmotherlabel, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
1499 {
1500        // mothertype 2 means this heavy quark coming from final state 
1501
1502        AliMCParticle *mctrack = NULL;
1503        if (inputmotherlabel > 5){ // mother exist after hard scattering
1504          if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(inputmotherlabel))))) return kFALSE; 
1505          TParticle *heavysMother = mctrack->Particle();
1506          motherID = heavysMother->GetPdgCode(); 
1507          mothertype = 2; // 
1508          motherlabel = inputmotherlabel;
1509          AliDebug(1,Form("heavy quark from %d after hard scattering! \n",motherID));
1510
1511          return kTRUE;
1512        }
1513        return kFALSE;
1514 }
1515
1516 //__________________________________________
1517 void AliHFEmcQA::ReportStrangeness(Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
1518 {
1519       // mark strange behavior  
1520
1521        mothertype = -888;
1522        motherlabel = -888;
1523        motherID = -888;
1524        AliDebug(1,"something strange!\n");
1525 }
1526
1527 //__________________________________________
1528 Int_t AliHFEmcQA::GetSource(const AliVParticle* const mcpart)
1529 {        
1530   // decay particle's origin 
1531
1532   //if ( abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return -1;
1533        
1534   Int_t origin = -1;
1535   Bool_t isFinalOpenCharm = kFALSE;
1536
1537   if(!mcpart){
1538     AliDebug(1, "Stack label is negative or no mcparticle, return\n");
1539     return -1;
1540   }
1541
1542   // mother
1543   // Information not in the base class, cast necessary
1544   Int_t iLabel = GetMother(mcpart);
1545   if (iLabel<0){
1546     AliDebug(1, "Stack label is negative, return\n");
1547     return -1;
1548   } 
1549        
1550   const AliVParticle *partMother = fMCEvent->GetTrack(iLabel);
1551   Int_t maPdgcode = partMother->PdgCode();
1552   
1553   // if the mother is charmed hadron  
1554   if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
1555     
1556     for (Int_t i=0; i<fNparents; i++){
1557        if (abs(maPdgcode)==fParentSelect[0][i]){
1558          isFinalOpenCharm = kTRUE;
1559        }
1560     }
1561     if (!isFinalOpenCharm) return -1;
1562
1563     // iterate until you find B hadron as a mother or become top ancester 
1564     for (Int_t i=1; i<fgkMaxIter; i++){
1565       
1566        Int_t jLabel = GetMother(partMother);
1567        if (jLabel == -1){
1568          origin = kDirectCharm;
1569          return origin;
1570        }
1571        if (jLabel < 0){ // safety protection
1572          AliDebug(1, "Stack label is negative, return\n");
1573          return -1;
1574        }
1575
1576        // if there is an ancester
1577        const AliVParticle* grandMa = fMCEvent->GetTrack(jLabel);
1578        Int_t grandMaPDG = grandMa->PdgCode();
1579
1580        for (Int_t j=0; j<fNparents; j++){
1581           if (abs(grandMaPDG)==fParentSelect[1][j]){
1582             origin = kBeautyCharm;
1583             return origin;
1584           }
1585        }
1586
1587        partMother = grandMa;
1588     } // end of iteration 
1589   } // end of if
1590   else if ( int(abs(maPdgcode)/100.) == kBeauty || int(abs(maPdgcode)/1000.) == kBeauty ) {
1591     for (Int_t i=0; i<fNparents; i++){
1592        if (abs(maPdgcode)==fParentSelect[1][i]){
1593          origin = kDirectBeauty;
1594          return origin;
1595        }
1596     }
1597   } // end of if
1598   else if ( abs(maPdgcode) == 22 ) {
1599     origin = kGamma;
1600     return origin;
1601   } // end of if
1602   else if ( abs(maPdgcode) == 111 ) {
1603     origin = kPi0;
1604     return origin;
1605   } // end of if
1606
1607   return origin;
1608 }
1609
1610 //__________________________________________
1611 Int_t AliHFEmcQA::GetSource(const TParticle * const mcpart)
1612 {
1613   // decay particle's origin 
1614
1615   //if ( abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return -1;
1616
1617   Int_t origin = -1;
1618   Bool_t isFinalOpenCharm = kFALSE;
1619
1620   if(!mcpart){
1621     AliDebug(1, "no mcparticle, return\n");
1622     return -1;
1623   }
1624
1625   Int_t iLabel = mcpart->GetFirstMother();
1626   if (iLabel<0){
1627     AliDebug(1, "Stack label is negative, return\n");
1628     return -1;
1629   }
1630
1631   AliMCParticle *mctrack = NULL;
1632   if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return -1; 
1633   TParticle *partMother = mctrack->Particle();
1634   Int_t maPdgcode = partMother->GetPdgCode();
1635
1636    // if the mother is charmed hadron  
1637    if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
1638
1639      for (Int_t i=0; i<fNparents; i++){
1640         if (abs(maPdgcode)==fParentSelect[0][i]){
1641           isFinalOpenCharm = kTRUE;
1642         }
1643      }
1644      if (!isFinalOpenCharm) return -1;
1645
1646      // iterate until you find B hadron as a mother or become top ancester 
1647      for (Int_t i=1; i<fgkMaxIter; i++){
1648
1649         Int_t jLabel = partMother->GetFirstMother();
1650         if (jLabel == -1){
1651           origin = kDirectCharm;
1652           return origin;
1653         }
1654         if (jLabel < 0){ // safety protection
1655           AliDebug(1, "Stack label is negative, return\n");
1656           return -1;
1657         }
1658
1659         // if there is an ancester
1660         if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return -1; 
1661         TParticle* grandMa = mctrack->Particle();
1662         Int_t grandMaPDG = grandMa->GetPdgCode();
1663
1664         for (Int_t j=0; j<fNparents; j++){
1665            if (abs(grandMaPDG)==fParentSelect[1][j]){
1666              origin = kBeautyCharm;
1667              return origin;
1668            }
1669         }
1670
1671         partMother = grandMa;
1672      } // end of iteration 
1673    } // end of if
1674    else if ( int(abs(maPdgcode)/100.) == kBeauty || int(abs(maPdgcode)/1000.) == kBeauty ) {
1675      for (Int_t i=0; i<fNparents; i++){
1676         if (abs(maPdgcode)==fParentSelect[1][i]){
1677           origin = kDirectBeauty;
1678           return origin;
1679         }
1680      }
1681    } // end of if
1682    else if ( abs(maPdgcode) == 22 ) {
1683      origin = kGamma;
1684      return origin;
1685    } // end of if
1686    else if ( abs(maPdgcode) == 111 ) {
1687      origin = kPi0;
1688      return origin;
1689    } // end of if
1690    else origin = kElse;
1691
1692    return origin;
1693 }
1694
1695 //__________________________________________
1696 Int_t AliHFEmcQA::GetElecSource(TParticle * const mcpart)
1697 {
1698   // decay particle's origin 
1699
1700   if(!mcpart){
1701     AliDebug(1, "no mcparticle, return\n");
1702     return -1;
1703   }
1704
1705   if ( abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return kMisID;
1706
1707   Int_t origin = -1;
1708   Bool_t isFinalOpenCharm = kFALSE;
1709
1710   Int_t iLabel = mcpart->GetFirstMother();
1711   if (iLabel<0){
1712     AliDebug(1, "Stack label is negative, return\n");
1713     return -1;
1714   }
1715
1716   AliMCParticle *mctrack = NULL;
1717   Int_t tmpMomLabel=0;
1718   if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return -1; 
1719   TParticle *partMother = mctrack->Particle();
1720   TParticle *partMotherCopy = mctrack->Particle();
1721   Int_t maPdgcode = partMother->GetPdgCode();
1722   Int_t grmaPdgcode;
1723   Int_t ggrmaPdgcode;
1724
1725    // if the mother is charmed hadron  
1726
1727    if(abs(maPdgcode)==443){ // J/spi
1728       Int_t jLabel = partMother->GetFirstMother();
1729       if((mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))){
1730         TParticle* grandMa = mctrack->Particle();
1731         Int_t grandMaPDG = grandMa->GetPdgCode();
1732         if((int(abs(grandMaPDG)/100.)%10) == kBeauty || (int(abs(grandMaPDG)/1000.)%10) == kBeauty) return kB2Jpsi;
1733       }
1734       return kJpsi;   
1735    } 
1736    else if ( (int(abs(maPdgcode)/100.)%10) == kCharm || (int(abs(maPdgcode)/1000.)%10) == kCharm ) {
1737
1738      for (Int_t i=0; i<fNparents; i++){
1739         if (abs(maPdgcode)==fParentSelect[0][i]){
1740           isFinalOpenCharm = kTRUE;
1741         }
1742      }
1743      if (!isFinalOpenCharm) return -1;
1744
1745      // iterate until you find B hadron as a mother or become top ancester 
1746      for (Int_t i=1; i<fgkMaxIter; i++){
1747
1748         Int_t jLabel = partMother->GetFirstMother();
1749         if (jLabel == -1){
1750           return kDirectCharm;
1751         }
1752         if (jLabel < 0){ // safety protection
1753           AliDebug(1, "Stack label is negative, return\n");
1754           return -1;
1755         }
1756
1757         // if there is an ancester
1758         if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return -1; 
1759         TParticle* grandMa = mctrack->Particle();
1760         Int_t grandMaPDG = grandMa->GetPdgCode();
1761
1762         for (Int_t j=0; j<fNparents; j++){
1763            if (abs(grandMaPDG)==fParentSelect[1][j]){
1764              return kBeautyCharm;
1765            }
1766         }
1767
1768         partMother = grandMa;
1769      } // end of iteration 
1770    } // end of if
1771    else if ( (int(abs(maPdgcode)/100.)%10) == kBeauty || (int(abs(maPdgcode)/1000.)%10) == kBeauty ) {
1772      for (Int_t i=0; i<fNparents; i++){
1773         if (abs(maPdgcode)==fParentSelect[1][i]){
1774           return kDirectBeauty;
1775         }
1776      }
1777    } // end of if
1778    else if ( abs(maPdgcode) == 22 ) { //conversion
1779
1780      tmpMomLabel = partMotherCopy->GetFirstMother();
1781      if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tmpMomLabel))))) return -1;
1782      partMother = mctrack->Particle();
1783      maPdgcode = partMother->GetPdgCode();
1784
1785      // check if the ligth meson is the decay product of heavy mesons
1786      tmpMomLabel = partMother->GetFirstMother();
1787      if((mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tmpMomLabel))))) {
1788       partMother = mctrack->Particle();
1789       grmaPdgcode = partMother->GetPdgCode();
1790
1791       if ( (int(abs(grmaPdgcode)/100.)%10) == kBeauty || (int(abs(grmaPdgcode)/1000.)%10) == kBeauty ) return kGammaB2M;
1792       if ( (int(abs(grmaPdgcode)/100.)%10) == kCharm || (int(abs(grmaPdgcode)/1000.)%10) == kCharm ) return kGammaD2M;
1793       if ( abs(grmaPdgcode) == 221 ) return kElse; //mj to remove secondary pi0 decay electrons from eta decays, mainly to eta signal enhance samples
1794
1795       tmpMomLabel = partMother->GetFirstMother();
1796       if((mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tmpMomLabel))))) {
1797        partMother = mctrack->Particle();
1798        ggrmaPdgcode = partMother->GetPdgCode();
1799
1800        if ( (int(abs(ggrmaPdgcode)/100.)%10) == kBeauty || (int(abs(ggrmaPdgcode)/1000.)%10) == kBeauty ) return kGammaB2M;
1801        if ( (int(abs(ggrmaPdgcode)/100.)%10) == kCharm || (int(abs(ggrmaPdgcode)/1000.)%10) == kCharm ) return kGammaD2M;
1802       }
1803      }
1804
1805      if ( abs(maPdgcode) == 111 ) {
1806        return kGammaPi0;
1807      } 
1808      else if ( abs(maPdgcode) == 221 ) {
1809        return kGammaEta;
1810      } 
1811      else if ( abs(maPdgcode) == 223 ) {
1812        return kGammaOmega;
1813      } 
1814      else if ( abs(maPdgcode) == 333 ) {
1815        return kGammaPhi;
1816      }
1817      else if ( abs(maPdgcode) == 331 ) {
1818        return kGammaEtaPrime; 
1819      }
1820      else if ( abs(maPdgcode) == 113 ) {
1821        return kGammaRho0;
1822      }
1823      else origin = kElse;
1824      return origin;
1825
1826    } 
1827    else {
1828
1829      // check if the ligth meson is the decay product of heavy mesons
1830      tmpMomLabel = partMotherCopy->GetFirstMother();
1831      if((mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tmpMomLabel))))) {
1832       partMother = mctrack->Particle();
1833       grmaPdgcode = partMother->GetPdgCode();
1834
1835       if ( (int(abs(grmaPdgcode)/100.)%10) == kBeauty || (int(abs(grmaPdgcode)/1000.)%10) == kBeauty ) return kGammaB2M;
1836       if ( (int(abs(grmaPdgcode)/100.)%10) == kCharm || (int(abs(grmaPdgcode)/1000.)%10) == kCharm ) return kGammaD2M;
1837       if ( abs(grmaPdgcode) == 221 ) return kElse; //mj to remove secondary pi0 decay electrons from eta decays, mainly to eta signal enhance samples
1838
1839       tmpMomLabel = partMother->GetFirstMother();
1840       if((mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tmpMomLabel))))) {
1841        partMother = mctrack->Particle();
1842        ggrmaPdgcode = partMother->GetPdgCode();
1843
1844        if ( (int(abs(ggrmaPdgcode)/100.)%10) == kBeauty || (int(abs(ggrmaPdgcode)/1000.)%10) == kBeauty ) return kGammaB2M;
1845        if ( (int(abs(ggrmaPdgcode)/100.)%10) == kCharm || (int(abs(ggrmaPdgcode)/1000.)%10) == kCharm ) return kGammaD2M;
1846       }
1847      }
1848
1849      if ( abs(maPdgcode) == 111 ) {
1850        return kPi0;
1851      } 
1852      else if ( abs(maPdgcode) == 221 ) {
1853        return kEta;
1854      } 
1855      else if ( abs(maPdgcode) == 223 ) {
1856        return kOmega;
1857      } 
1858      else if ( abs(maPdgcode) == 333 ) {
1859        return kPhi;
1860      } 
1861      else if ( abs(maPdgcode) == 331 ) {
1862        return kEtaPrime;
1863      } 
1864      else if ( abs(maPdgcode) == 113 ) {
1865        return kRho0;
1866      } 
1867      else if ( abs(maPdgcode) == 321 || abs(maPdgcode) == 130 ) {
1868        return kKe3;
1869      }
1870      else{ 
1871       origin = kElse;
1872      }
1873
1874    }
1875    return origin;
1876 }
1877 //__________________________________________
1878 Double_t AliHFEmcQA::GetWeightFactor(AliMCParticle *mctrack, const Int_t iBgLevel){
1879   //
1880   // 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)
1881   //
1882   AliMCParticle *mctrackmother = NULL;  
1883   Double_t weightElecBg = 0.;
1884   Double_t mesonPt = 0.;
1885   Double_t bgcategory = 0.;
1886   Int_t mArr = -1;  
1887   Int_t mesonID = GetElecSource(mctrack->Particle());
1888   if(mesonID==kGammaPi0 || mesonID==kPi0) mArr=0;                //pion
1889   else if(mesonID==kGammaEta || mesonID==kEta) mArr=1;           //eta
1890   else if(mesonID==kGammaOmega || mesonID==kOmega) mArr=2;       //omega
1891   else if(mesonID==kGammaPhi || mesonID==kPhi) mArr=3;           //phi
1892   else if(mesonID==kGammaEtaPrime || mesonID==kEtaPrime) mArr=4; //etaprime
1893   else if(mesonID==kGammaRho0 || mesonID==kRho0) mArr=5;         //rho
1894
1895   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};
1896   Double_t xr[3]={-999,-999,-999};
1897   datamc[0] = mesonID;
1898   datamc[17] = mctrack->Pt(); //electron pt
1899   datamc[18] = mctrack->Eta(); //electron eta
1900
1901   mctrack->XvYvZv(xr);
1902   datamc[9] = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
1903   datamc[10] = xr[2];
1904
1905   TParticle *mcpart = mctrack->Particle();
1906   if(mcpart){
1907     datamc[14] = mcpart->GetUniqueID();
1908   }
1909
1910   if(!(mArr<0)){
1911      if(mesonID>=kGammaPi0) {  // conversion electron, be careful with the enum odering 
1912         Int_t glabel=TMath::Abs(mctrack->GetMother()); // gamma label
1913         if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
1914           glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's label
1915           if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
1916             mesonPt = mctrackmother->Pt(); //meson pt
1917             bgcategory = 1.;
1918             datamc[1] = bgcategory;
1919             datamc[2] = mesonPt;
1920             mctrackmother->XvYvZv(xr);
1921             datamc[11] = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
1922             datamc[12] = xr[2];
1923
1924             mcpart = mctrackmother->Particle();
1925             if(mcpart){
1926               datamc[15] = mcpart->GetUniqueID();
1927             }
1928             if(glabel>fMCEvent->GetNumberOfPrimaries()) {
1929               bgcategory = 2.;
1930               datamc[1] = bgcategory;
1931               //printf("I should be gamma meson = %d  mesonlabel= %d  NumberOfPrimaries= %d \n",mctrackmother->PdgCode(),glabel,fMCEvent->GetNumberOfPrimaries()); 
1932               glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother
1933               if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
1934                 datamc[3]=mctrackmother->PdgCode();
1935                 datamc[4]=mctrackmother->Pt();
1936                 if(TMath::Abs(mctrackmother->PdgCode())==310){
1937                   bgcategory = 3.;
1938                   glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother's mother
1939                   if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
1940                     datamc[5]=mctrackmother->PdgCode();
1941                     glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother
1942                     if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
1943                        datamc[6]=mctrackmother->PdgCode();
1944                     }
1945                   }
1946                 }
1947               }
1948             } 
1949           } 
1950         }
1951      }
1952      else{ // nonHFE except for the conversion electron
1953         Int_t glabel=TMath::Abs(mctrack->GetMother()); 
1954         if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
1955           mesonPt=mctrackmother->Pt(); //meson pt
1956           bgcategory = -1.;
1957           datamc[1] = bgcategory;
1958           datamc[2] = mesonPt;
1959           mctrackmother->XvYvZv(xr);
1960           datamc[11] = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
1961           datamc[12] = xr[2];
1962
1963           mcpart = mctrackmother->Particle();
1964           if(mcpart){
1965             datamc[15] = mcpart->GetUniqueID();
1966           }
1967           if(glabel>fMCEvent->GetNumberOfPrimaries()) {
1968             bgcategory = -2.;
1969             datamc[1] = bgcategory;
1970             glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother
1971             if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
1972               datamc[3]=mctrackmother->PdgCode();
1973               datamc[4]=mctrackmother->Pt();
1974               if(TMath::Abs(mctrackmother->PdgCode())==310){
1975                bgcategory = -3.;
1976                glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother's mother
1977                if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
1978                  datamc[5]=mctrackmother->PdgCode();
1979                  glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother
1980                  if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
1981                    datamc[6]=mctrackmother->PdgCode();
1982                  }
1983                }
1984               }  
1985             }
1986           }
1987         }
1988      }
1989
1990
1991      if(fIsPbPb){
1992        Int_t centBin = GetWeightCentralityBin(fPerCentrality);
1993        if(centBin < 0)return 0.;
1994        weightElecBg=fElecBackgroundFactor[iBgLevel][centBin][mArr][kBgPtBins-1];                        
1995        for(int ii=0; ii<kBgPtBins; ii++){              
1996          if((mesonPt > fBinLimit[ii]) && (mesonPt < fBinLimit[ii+1])){
1997            weightElecBg = fElecBackgroundFactor[iBgLevel][centBin][mArr][ii];
1998            break;
1999          }
2000        }
2001      }
2002      else{
2003        weightElecBg=fElecBackgroundFactor[iBgLevel][0][mArr][kBgPtBins-1];                         
2004        for(int ii=0; ii<kBgPtBins; ii++){              
2005          if((mesonPt > fBinLimit[ii]) && (mesonPt < fBinLimit[ii+1])){
2006            weightElecBg = fElecBackgroundFactor[iBgLevel][0][mArr][ii];
2007            break;
2008          }
2009        }
2010      }    
2011   }
2012
2013   datamc[13] = weightElecBg;
2014   datamc[16] = Double_t(fContainerStep);
2015
2016   datamc[7] = fHfeImpactR;
2017   datamc[8] = fHfeImpactnsigmaR;
2018
2019   datamc[19] = fRecPt;
2020   datamc[20] = fRecEta;
2021   datamc[21] = fRecPhi;
2022   datamc[22] = fLyrhit;
2023   datamc[23] = fLyrstat;
2024   datamc[24] = fCentrality;
2025
2026   if(fIsDebugStreamerON && fTreeStream){
2027    if(!iBgLevel){
2028     (*fTreeStream)<<"nonhfeQA"<<
2029         "mesonID="<<datamc[0]<<
2030         "bgcategory="<<datamc[1]<<
2031         "mesonPt="<<datamc[2]<<
2032         "mesonMomPdg="<<datamc[3]<<
2033         "mesonMomPt="<<datamc[4]<<
2034         "mesonGMomPdg="<<datamc[5]<<
2035         "mesonGGMomPdg="<<datamc[6]<<
2036         "eIPAbs="<<datamc[7]<<
2037         "eIPSig="<<datamc[8]<<
2038         "eR="<<datamc[9]<<
2039         "eZ="<<datamc[10]<<
2040         "mesonR="<<datamc[11]<<
2041         "mesonZ="<<datamc[12]<<
2042         "weightElecBg="<<datamc[13]<< 
2043         "eUniqID="<<datamc[14]<<
2044         "mesonUniqID="<<datamc[15]<<
2045         "containerStep="<<datamc[16]<<
2046         "emcpt="<<datamc[17]<<
2047         "emceta="<<datamc[18]<<
2048         "erecpt="<<datamc[19]<<
2049         "ereceta="<<datamc[20]<<
2050         "erecphi="<<datamc[21]<< 
2051         "itshit="<<datamc[22]<<
2052         "itsstat="<<datamc[23]<<
2053         "centrality="<<datamc[24]
2054         << "\n";
2055    }
2056   }
2057
2058   Double_t returnval = bgcategory*weightElecBg;
2059   if(TMath::Abs(bgcategory)>1) returnval = bgcategory/2.;
2060
2061   return returnval;
2062 }
2063
2064 //__________________________________________
2065 Int_t AliHFEmcQA::GetMother(const AliVParticle * const mcpart){
2066   //
2067   // Wrapper to get the mother label
2068   //
2069   Int_t label = -1; 
2070   TClass *type = mcpart->IsA();
2071   if(type == AliMCParticle::Class()){
2072     // ESD analysis
2073     const AliMCParticle *emcpart = static_cast<const AliMCParticle *>(mcpart);
2074     label = emcpart->GetMother();
2075   } else if(type == AliAODMCParticle::Class()){
2076     // AOD analysis
2077     const AliAODMCParticle *amcpart = static_cast<const AliAODMCParticle *>(mcpart);
2078     label = amcpart->GetMother();
2079   }
2080   return label;
2081 }
2082 //__________________________________________
2083 Int_t AliHFEmcQA::GetWeightCentralityBin(const Float_t percentile) const {
2084   //
2085   //translate the centrality percentile into the centrality bin of the reference weighting histograms for electron background
2086   //
2087
2088   Float_t centralityLimits[12]= {0.,5.,10., 20., 30., 40., 50., 60.,70.,80., 90., 100.};
2089   Int_t bin = -1;
2090   for(Int_t ibin = 0; ibin < 11; ibin++){
2091     if(percentile >= centralityLimits[ibin] && percentile < centralityLimits[ibin+1]){
2092       bin = ibin;
2093       break;
2094     }
2095   } 
2096   return bin;
2097 }
2098 //__________________________________________
2099 void AliHFEmcQA::AliHists::FillList(TList *l) const {
2100   //
2101   // Fill Histos into a list for output
2102   //
2103   if(fPdgCode) l->Add(fPdgCode);
2104   if(fPt) l->Add(fPt);
2105   if(fY) l->Add(fY);
2106   if(fEta) l->Add(fEta);
2107 }
2108
2109 //__________________________________________
2110 void AliHFEmcQA::AliHistsComm::FillList(TList *l) const { 
2111   //
2112   // Fill Histos into a list for output
2113   //
2114   if(fNq) l->Add(fNq);
2115   if(fProcessID) l->Add(fProcessID);
2116   if(fePtRatio) l->Add(fePtRatio);
2117   if(fPtCorr) l->Add(fPtCorr);
2118   if(fPtCorrDp) l->Add(fPtCorrDp);
2119   if(fPtCorrD0) l->Add(fPtCorrD0);
2120   if(fPtCorrDrest) l->Add(fPtCorrDrest);
2121
2122   if(fPtCorrDinein) l->Add(fPtCorrDinein);
2123   if(fPtCorrDineout) l->Add(fPtCorrDineout);
2124   if(fPtCorrDoutein) l->Add(fPtCorrDoutein);
2125   if(fPtCorrDouteout) l->Add(fPtCorrDouteout);
2126   if(fPtCorrDpDinein) l->Add(fPtCorrDpDinein);
2127   if(fPtCorrDpDineout) l->Add(fPtCorrDpDineout);
2128   if(fPtCorrDpDoutein) l->Add(fPtCorrDpDoutein);
2129   if(fPtCorrDpDouteout) l->Add(fPtCorrDpDouteout);
2130   if(fPtCorrD0Dinein) l->Add(fPtCorrD0Dinein);
2131   if(fPtCorrD0Dineout) l->Add(fPtCorrD0Dineout);
2132   if(fPtCorrD0Doutein) l->Add(fPtCorrD0Doutein);
2133   if(fPtCorrD0Douteout) l->Add(fPtCorrD0Douteout);
2134   if(fPtCorrDrestDinein) l->Add(fPtCorrDrestDinein);
2135   if(fPtCorrDrestDineout) l->Add(fPtCorrDrestDineout);
2136   if(fPtCorrDrestDoutein) l->Add(fPtCorrDrestDoutein);
2137   if(fPtCorrDrestDouteout) l->Add(fPtCorrDrestDouteout);
2138
2139   if(fEtaCorrD) l->Add(fEtaCorrD);
2140   if(fEtaCorrDp) l->Add(fEtaCorrDp);
2141   if(fEtaCorrD0) l->Add(fEtaCorrD0);
2142   if(fEtaCorrDrest) l->Add(fEtaCorrDrest);
2143
2144   if(fEtaCorrGD) l->Add(fEtaCorrGD);
2145   if(fEtaCorrGDp) l->Add(fEtaCorrGDp);
2146   if(fEtaCorrGD0) l->Add(fEtaCorrGD0);
2147   if(fEtaCorrGDrest) l->Add(fEtaCorrGDrest);
2148
2149   if(fEtaCorrB) l->Add(fEtaCorrB);
2150   if(fEtaCorrGB) l->Add(fEtaCorrGB);
2151   if(fPtCorrBinein) l->Add(fPtCorrBinein);
2152   if(fPtCorrBineout) l->Add(fPtCorrBineout);
2153   if(fPtCorrBoutein) l->Add(fPtCorrBoutein);
2154   if(fPtCorrBouteout) l->Add(fPtCorrBouteout);
2155
2156   if(fDePtRatio) l->Add(fDePtRatio);
2157   if(feDistance) l->Add(feDistance);
2158   if(fDeDistance) l->Add(fDeDistance);
2159 }
2160