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