Major update of the HFE package (comments inside the code
[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
41 #include "AliHFEmcQA.h"
42 #include "AliHFEtools.h"
43
44 const Int_t AliHFEmcQA::fgkGluon=21;
45 const Int_t AliHFEmcQA::fgkMaxGener=10;
46 const Int_t AliHFEmcQA::fgkMaxIter=100;
47 const Int_t AliHFEmcQA::fgkqType = 7;    // number of species waiting for QA done
48
49 ClassImp(AliHFEmcQA)
50
51 //_______________________________________________________________________________________________
52 AliHFEmcQA::AliHFEmcQA() : 
53         fMCEvent(NULL) 
54         ,fMCHeader(NULL)
55         ,fMCArray(NULL)
56         ,fQAhistos(NULL)
57         ,fNparents(0) 
58 {
59         // Default constructor
60         
61 }
62
63 //_______________________________________________________________________________________________
64 AliHFEmcQA::AliHFEmcQA(const AliHFEmcQA&p):
65         TObject(p)
66         ,fMCEvent(NULL) 
67         ,fMCHeader(NULL)
68         ,fMCArray(NULL)
69         ,fQAhistos(p.fQAhistos)
70         ,fNparents(p.fNparents) 
71 {
72         // Copy constructor
73 }
74
75 //_______________________________________________________________________________________________
76 AliHFEmcQA&
77 AliHFEmcQA::operator=(const AliHFEmcQA &)
78 {
79   // Assignment operator
80
81   AliInfo("Not yet implemented.");
82   return *this;
83 }
84
85 //_______________________________________________________________________________________________
86 AliHFEmcQA::~AliHFEmcQA()
87 {
88         // Destructor
89
90         AliInfo("Analysis Done.");
91 }
92
93 //_______________________________________________________________________________________________
94 void AliHFEmcQA::PostAnalyze() const
95 {
96         //
97         // Post analysis
98         //
99 }
100
101 //__________________________________________
102 void AliHFEmcQA::CreatDefaultHistograms(TList * const qaList)
103 {      
104   //
105   // make default histograms
106   //
107   
108   if(!qaList) return;
109
110   fQAhistos = qaList;
111   fQAhistos->SetName("MCqa");
112
113   CreateHistograms(AliHFEmcQA::kCharm,0,"mcqa_");               // create histograms for charm
114   CreateHistograms(AliHFEmcQA::kBeauty,0,"mcqa_");              // create histograms for beauty
115   CreateHistograms(AliHFEmcQA::kOthers,0,"mcqa_");              // create histograms for beauty
116   CreateHistograms(AliHFEmcQA::kCharm,1,"mcqa_barrel_");        // create histograms for charm 
117   CreateHistograms(AliHFEmcQA::kBeauty,1,"mcqa_barrel_");       // create histograms for beauty
118   CreateHistograms(AliHFEmcQA::kOthers,1,"mcqa_barrel_");       // create histograms for beauty
119   CreateHistograms(AliHFEmcQA::kCharm,2,"mcqa_unitY_");         // create histograms for charm 
120   CreateHistograms(AliHFEmcQA::kBeauty,2,"mcqa_unitY_");        // create histograms for beauty
121   CreateHistograms(AliHFEmcQA::kOthers,2,"mcqa_unitY_");        // create histograms for beauty
122   CreateHistograms(AliHFEmcQA::kCharm,3,"mcqa_reccut_");        // create histograms for charm 
123   CreateHistograms(AliHFEmcQA::kBeauty,3,"mcqa_reccut_");       // create histograms for beauty
124   CreateHistograms(AliHFEmcQA::kOthers,3,"mcqa_reccut_");       // create histograms for beauty
125   CreateHistograms(AliHFEmcQA::kCharm,4,"mcqa_recpidcut_");     // create histograms for charm 
126   CreateHistograms(AliHFEmcQA::kBeauty,4,"mcqa_recpidcut_");    // create histograms for beauty
127   CreateHistograms(AliHFEmcQA::kOthers,4,"mcqa_recpidcut_");    // create histograms for beauty
128   CreateHistograms(AliHFEmcQA::kCharm,5,"mcqa_secvtxcut_");     // create histograms for charm 
129   CreateHistograms(AliHFEmcQA::kBeauty,5,"mcqa_secvtxcut_");    // create histograms for beauty
130   CreateHistograms(AliHFEmcQA::kOthers,5,"mcqa_secvtxcut_");    // create histograms for beauty
131  
132 }
133   
134 //__________________________________________
135 void AliHFEmcQA::CreateHistograms(const Int_t kquark, Int_t icut, TString hnopt) 
136 {
137   // create histograms
138
139   if (!(kquark == kCharm || kquark == kBeauty || kquark == kOthers)) {
140     AliDebug(1, "This task is only for heavy quark QA, return\n");
141     return; 
142   }
143   Int_t iq = kquark - kCharm; 
144
145   TString kqTypeLabel[fgkqType];
146   if (kquark == kCharm){
147     kqTypeLabel[kQuark]="c";
148     kqTypeLabel[kantiQuark]="cbar";
149     kqTypeLabel[kHadron]="cHadron";
150     kqTypeLabel[keHadron]="ceHadron";
151     kqTypeLabel[kDeHadron]="nullHadron";
152     kqTypeLabel[kElectron]="ce";
153     kqTypeLabel[kElectron2nd]="nulle";
154   } else if (kquark == kBeauty){
155     kqTypeLabel[kQuark]="b";
156     kqTypeLabel[kantiQuark]="bbar";
157     kqTypeLabel[kHadron]="bHadron";
158     kqTypeLabel[keHadron]="beHadron";
159     kqTypeLabel[kDeHadron]="bDeHadron";
160     kqTypeLabel[kElectron]="be";
161     kqTypeLabel[kElectron2nd]="bce";
162   } else if (kquark == kOthers){
163     kqTypeLabel[kGamma-4]="gammae";
164     kqTypeLabel[kPi0-4]="pi0e";
165     kqTypeLabel[kElse-4]="elsee";
166     kqTypeLabel[kMisID-4]="miside";
167   }
168
169   const Double_t kPtbound[2] = {0.1, 20.}; //bin taken for considering inclusive e analysis binning
170   Int_t iBin[1];
171   iBin[0] = 44; // bins in pt
172   Double_t* binEdges[1];
173   binEdges[0] =  AliHFEtools::MakeLogarithmicBinning(iBin[0], kPtbound[0], kPtbound[1]);
174   
175   
176   TString hname; 
177   if(kquark == kOthers){
178     for (Int_t iqType = 0; iqType < 4; iqType++ ){
179        hname = hnopt+"Pt_"+kqTypeLabel[iqType];
180        fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",iBin[0],binEdges[0]);
181        //fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",500,0,50);
182        hname = hnopt+"Y_"+kqTypeLabel[iqType];
183        fHist[iq][iqType][icut].fY = new TH1F(hname,hname,150,-7.5,7.5);
184        hname = hnopt+"Eta_"+kqTypeLabel[iqType];
185        fHist[iq][iqType][icut].fEta = new TH1F(hname,hname,150,-7.5,7.5);
186        // Fill List
187        if(fQAhistos) fHist[iq][iqType][icut].FillList(fQAhistos);
188     }
189     return;
190   }
191   for (Int_t iqType = 0; iqType < fgkqType; iqType++ ){
192      if (iqType < keHadron && icut > 0) continue; // don't duplicate histogram for quark and hadron
193      hname = hnopt+"PdgCode_"+kqTypeLabel[iqType];
194      fHist[iq][iqType][icut].fPdgCode = new TH1F(hname,hname,20001,-10000.5,10000.5);
195      hname = hnopt+"Pt_"+kqTypeLabel[iqType];
196      fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",iBin[0],binEdges[0]);
197      //fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",500,0,50);
198      hname = hnopt+"Y_"+kqTypeLabel[iqType];
199      fHist[iq][iqType][icut].fY = new TH1F(hname,hname,150,-7.5,7.5);
200      hname = hnopt+"Eta_"+kqTypeLabel[iqType];
201      fHist[iq][iqType][icut].fEta = new TH1F(hname,hname,150,-7.5,7.5);
202      // Fill List
203      if(fQAhistos) fHist[iq][iqType][icut].FillList(fQAhistos);
204   }
205
206   if (icut == 0){ 
207     hname = hnopt+"Nq_"+kqTypeLabel[kQuark];
208     fHistComm[iq][icut].fNq = new TH1F(hname,hname,10,-0.5,9.5);
209     hname = hnopt+"ProcessID_"+kqTypeLabel[kQuark];
210     fHistComm[iq][icut].fProcessID = new TH1F(hname,hname,21,-10.5,10.5);
211   }
212   hname = hnopt+"ePtRatio_"+kqTypeLabel[kQuark];
213   fHistComm[iq][icut].fePtRatio = new TH2F(hname,hname+";p_{T} (GeV/c);momentum fraction",100,0,30,100,0,1);
214   hname = hnopt+"DePtRatio_"+kqTypeLabel[kQuark];
215   fHistComm[iq][icut].fDePtRatio = new TH2F(hname,hname+";p_{T} (GeV/c);momentum fraction",100,0,30,100,0,1);
216   hname = hnopt+"eDistance_"+kqTypeLabel[kQuark];
217   fHistComm[iq][icut].feDistance= new TH2F(hname,hname+";p_{T} (GeV/c);distance (cm)",100,0,30,200,0,2);
218   hname = hnopt+"DeDistance_"+kqTypeLabel[kQuark];
219   fHistComm[iq][icut].fDeDistance= new TH2F(hname,hname+";p_{T} (GeV/c);distance (cm)",100,0,30,200,0,2);
220   if(fQAhistos) fHistComm[iq][icut].FillList(fQAhistos);
221 }
222
223 //__________________________________________
224 void AliHFEmcQA::Init()
225 {
226   // called at begining every event
227   
228   for (Int_t i=0; i<2; i++){
229      fIsHeavy[i] = 0;
230   } 
231
232   fNparents = 7;
233
234   fParentSelect[0][0] =  411;  
235   fParentSelect[0][1] =  421;
236   fParentSelect[0][2] =  431;
237   fParentSelect[0][3] = 4122;
238   fParentSelect[0][4] = 4132;
239   fParentSelect[0][5] = 4232;
240   fParentSelect[0][6] = 4332;
241
242   fParentSelect[1][0] =  511;
243   fParentSelect[1][1] =  521;
244   fParentSelect[1][2] =  531;
245   fParentSelect[1][3] = 5122;
246   fParentSelect[1][4] = 5132;
247   fParentSelect[1][5] = 5232;
248   fParentSelect[1][6] = 5332;
249
250 }
251
252 //__________________________________________
253 void AliHFEmcQA::GetQuarkKine(TParticle *part, Int_t iTrack, const Int_t kquark) 
254 {
255   // get heavy quark kinematics
256
257     if (kquark != kCharm && kquark != kBeauty) {
258       AliDebug(1, "This task is only for heavy quark QA, return\n");
259       return; 
260     }
261     Int_t iq = kquark - kCharm; 
262
263     if (iTrack < 0 || !part) { 
264       AliDebug(1, "Stack label is negative or no mcparticle, return\n");
265       return; 
266     }
267
268     AliMCParticle *mctrack = NULL;
269     Int_t partPdgcode = TMath::Abs(part->GetPdgCode());
270
271     // select heavy hadron or not fragmented heavy quark 
272     if ( int(partPdgcode/100.)==kquark || int(partPdgcode/1000.)==kquark || (partPdgcode==kquark && (part->GetNDaughters()==0 && iTrack>5)) ){ 
273
274       TParticle *partMother;
275       Int_t iLabel;
276
277       if (partPdgcode == kquark){ // in case of not fragmented heavy quark  
278         partMother = part; 
279         iLabel = iTrack;
280       } else{ // in case of heavy hadron, start to search for mother heavy parton 
281         iLabel = part->GetFirstMother(); 
282         if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return; 
283         if (iLabel>-1) { partMother = mctrack->Particle(); }
284         else {
285           AliDebug(1, "Stack label is negative, return\n");
286           return; 
287         }
288       }
289
290       // heavy parton selection as a mother of heavy hadron 
291       // if the heavy particle comes from string which is denoted as particle status 12|12|12...12|11,[PYTHIA p.60]
292       // in this case, the mother of heavy particle can be one of the fragmented parton of the string
293       // should I make a condition that partMother should be quark or diquark? -> not necessary
294       if ( abs(partMother->GetPdgCode()) == kquark || (partMother->GetStatusCode() == 11) ){
295       //if ( abs(partMother->GetPdgCode()) == kquark || (partMother->GetStatusCode() == 11 || partMother->GetStatusCode() == 12) ){
296
297         if ( abs(partMother->GetPdgCode()) != kquark ){
298           // search fragmented partons in the same string
299           Bool_t isSameString = kTRUE; 
300           for (Int_t i=1; i<fgkMaxIter; i++){
301              iLabel = iLabel - 1;
302              if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return; 
303              if (iLabel>-1) { partMother = mctrack->Particle(); }
304              else {
305                AliDebug(1, "Stack label is negative, return\n");
306                return; 
307              }
308              if ( abs(partMother->GetPdgCode()) == kquark ) break;
309              if ( partMother->GetStatusCode() != 12 ) isSameString = kFALSE;
310              if (!isSameString) return; 
311           }
312         }
313         AliDebug(1, "Can not find heavy parton of this heavy hadron in the string, return\n");
314         if (abs(partMother->GetPdgCode()) != kquark) return; 
315
316         if (fIsHeavy[iq] >= 50) return;  
317         fHeavyQuark[fIsHeavy[iq]] = partMother;
318         fIsHeavy[iq]++;
319
320         // fill kinematics for heavy parton
321         if (partMother->GetPdgCode() > 0) { // quark
322           fHist[iq][kQuark][0].fPdgCode->Fill(partMother->GetPdgCode());
323           fHist[iq][kQuark][0].fPt->Fill(partMother->Pt());
324           fHist[iq][kQuark][0].fY->Fill(AliHFEtools::GetRapidity(partMother));
325           fHist[iq][kQuark][0].fEta->Fill(partMother->Eta());
326         } else{ // antiquark
327           fHist[iq][kantiQuark][0].fPdgCode->Fill(partMother->GetPdgCode());
328           fHist[iq][kantiQuark][0].fPt->Fill(partMother->Pt());
329           fHist[iq][kantiQuark][0].fY->Fill(AliHFEtools::GetRapidity(partMother));
330           fHist[iq][kantiQuark][0].fEta->Fill(partMother->Eta());
331         }
332
333       } // end of heavy parton slection loop 
334
335     } // end of heavy hadron or quark selection
336
337 }
338
339 //__________________________________________
340 void AliHFEmcQA::EndOfEventAna(const Int_t kquark)
341 {
342   // end of event analysis
343
344   if (kquark != kCharm && kquark != kBeauty) {
345     AliDebug(1, "This task is only for heavy quark QA, return\n");
346     return; 
347   }
348   Int_t iq = kquark - kCharm; 
349
350
351   // # of heavy quark per event
352   AliDebug(1,Form("Number of heavy quark in this event = %d \n",fIsHeavy[iq]));
353   fHistComm[iq][0].fNq->Fill(fIsHeavy[iq]);
354
355   Int_t motherID[fgkMaxGener];
356   Int_t motherType[fgkMaxGener];
357   Int_t motherLabel[fgkMaxGener];
358   Int_t ancestorPdg[fgkMaxGener];
359   Int_t ancestorLabel[fgkMaxGener];
360
361   for (Int_t i = 0; i < fgkMaxGener; i++){ // initialization
362      motherID[i] = 0;
363      motherType[i] = 0;
364      motherLabel[i] = 0;
365      ancestorPdg[i] = 0;
366      ancestorLabel[i] = 0;
367   }
368
369
370   // check history of found heavy quarks
371   for (Int_t i = 0; i < fIsHeavy[iq]; i++){
372
373      if(!fHeavyQuark[i]) return;
374
375      ancestorLabel[0] = i;
376      ancestorPdg[0] = fHeavyQuark[i]->GetPdgCode(); 
377      ancestorLabel[1] = fHeavyQuark[i]->GetFirstMother(); 
378
379      AliDebug(1,Form("pdg code= %d\n",ancestorPdg[0]));
380      AliDebug(1,Form("ancestor label= %d\n",ancestorLabel[1]));
381
382      Int_t ig = 1;
383      while (ancestorLabel[ig] != -1){
384           // in case there is mother, get mother's pdg code and grandmother's label
385           IdentifyMother(ancestorLabel[ig], ancestorPdg[ig], ancestorLabel[ig+1]); 
386           // if mother is still heavy, find again mother's ancestor
387           if (ancestorPdg[ig-1] == ancestorPdg[ig]) {
388             ig++;
389             continue; // if it is from same heavy
390           }
391           // if the heavy's mother is not heavy, check the mother's label to know if it comes from inital or final parton shower
392           if (IsFromInitialShower(ancestorLabel[ig],motherID[i],motherType[i],motherLabel[i])) break;
393           if (IsFromFinalParton(ancestorLabel[ig],motherID[i],motherType[i],motherLabel[i])) break;
394           // if it is not the above case, something is strange
395           ReportStrangeness(motherID[i],motherType[i],motherLabel[i]);
396           break;
397      } 
398      if (ancestorLabel[ig] == -1){ // from hard scattering
399        HardScattering(kquark, motherID[i],motherType[i], motherLabel[i]);
400      }
401
402   } // end of found heavy quark loop
403
404
405   // check process type
406   Int_t processID = 0;
407   for (Int_t i = 0; i < fIsHeavy[iq]; i++){
408      AliDebug(1,Form("Mother ID= %d type= %d label= %d\n",motherID[i],motherType[i],motherLabel[i]));
409   }
410
411
412   Int_t nheavypair = Int_t(fIsHeavy[iq]/2.); 
413   for (Int_t ipair = 0; ipair < nheavypair; ipair++){
414
415      Int_t id1 = ipair*2;
416      Int_t id2 = ipair*2 + 1;
417
418      if (motherType[id1] == 2 && motherType[id2] == 2){
419        if (motherLabel[id1] == motherLabel[id2]) processID = kGluonSplitting; // gluon spliting
420        else processID = -9;
421      }
422      else if (motherType[id1] == -1 && motherType[id2] == -1) {
423        if (motherLabel[id1] == -1 && motherLabel[id2] == -1) {
424          if (motherID[id1] == fgkGluon) processID = kPairCreationFromg; // gluon fusion
425          else processID = kPairCreationFromq; // q-qbar pair creation
426        }
427        else processID = -8;
428      }
429      else if (motherType[id1] == -1 || motherType[id2] == -1) {
430        if ((motherLabel[id1] == -1 || motherLabel[id2] == -1) && (motherLabel[id1]*motherLabel[id2] == -2 || motherLabel[id1]*motherLabel[id2] == -3)) {
431          if(motherID[id1]*motherID[id2] == kquark*fgkGluon) processID = kFlavourExitation; // flavour exitation 
432          else processID = kLightQuarkShower;
433        }
434        else processID = -7;
435      }
436      else if (motherType[id1] == -2 || motherType[id2] == -2) {
437        if (motherLabel[id1] == motherLabel[id2]) processID = kInitialPartonShower; // initial parton shower
438        else processID = -6;
439        
440      }
441      else processID = -5;
442
443      if (nheavypair >1) AliDebug(1,Form("Multi pair found : process ID = %d\n",processID));
444      else fHistComm[iq][0].fProcessID->Fill(processID);
445      AliDebug(1,Form("Process ID = %d\n",processID));
446   } // end of # heavy quark pair loop
447
448 }
449
450 //__________________________________________
451 void AliHFEmcQA::GetHadronKine(TParticle* mcpart, const Int_t kquark)
452 {
453     // decay electron kinematics
454
455     if (kquark != kCharm && kquark != kBeauty) {
456       AliDebug(1, "This task is only for heavy quark QA, return\n");
457       return;
458     }
459     Int_t iq = kquark - kCharm;
460
461     if(!mcpart){
462       AliDebug(1, "no mc particle, return\n");
463       return;
464     }
465
466     Int_t iLabel = mcpart->GetFirstMother();
467     if (iLabel<0){
468       AliDebug(1, "Stack label is negative, return\n");
469       return;
470     }
471
472     TParticle *partCopy = mcpart;
473     Int_t pdgcode = mcpart->GetPdgCode();
474     Int_t pdgcodeCopy = pdgcode;
475
476     AliMCParticle *mctrack = NULL;
477
478     // if the mother is charmed hadron  
479     Bool_t isDirectCharm = kFALSE;
480     if ( int(abs(pdgcode)/100.) == kCharm || int(abs(pdgcode)/1000.) == kCharm ) {
481
482           // iterate until you find B hadron as a mother or become top ancester 
483           for (Int_t i=1; i<fgkMaxIter; i++){
484
485              Int_t jLabel = mcpart->GetFirstMother();
486              if (jLabel == -1){
487                isDirectCharm = kTRUE;
488                break; // if there is no ancester
489              }
490              if (jLabel < 0){ // safety protection
491                AliDebug(1, "Stack label is negative, return\n");
492                return;
493              }
494              // if there is an ancester
495              if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return; 
496              TParticle* mother = mctrack->Particle();
497              Int_t motherPDG = mother->GetPdgCode();
498     
499              for (Int_t j=0; j<fNparents; j++){
500                 if (abs(motherPDG)==fParentSelect[1][j]) return; // return if this hadron is originated from b
501              }
502
503              mcpart = mother;
504           } // end of iteration 
505     } // end of if
506     if((isDirectCharm == kTRUE && kquark == kCharm) || kquark == kBeauty) {
507          for (Int_t i=0; i<fNparents; i++){
508             if (abs(pdgcodeCopy)==fParentSelect[iq][i]){
509
510               // fill hadron kinematics
511               fHist[iq][kHadron][0].fPdgCode->Fill(pdgcodeCopy);
512               fHist[iq][kHadron][0].fPt->Fill(partCopy->Pt());
513               fHist[iq][kHadron][0].fY->Fill(AliHFEtools::GetRapidity(partCopy));
514               fHist[iq][kHadron][0].fEta->Fill(partCopy->Eta());
515
516             }
517          }
518     } // end of if
519 }
520
521 //__________________________________________
522 void AliHFEmcQA::GetDecayedKine(TParticle* mcpart, const Int_t kquark, Int_t kdecayed, Int_t icut) 
523 {
524     // decay electron kinematics
525     
526     if (!(kquark == kCharm || kquark == kBeauty || kquark == kOthers)){
527       AliDebug(1, "This task is only for heavy quark QA, return\n");
528       return; 
529     }
530     Int_t iq = kquark - kCharm; 
531     Bool_t isFinalOpenCharm = kFALSE;
532
533     if(!mcpart){
534       AliDebug(1, "no mcparticle, return\n");
535       return;
536     }
537
538     if(kquark==kOthers){
539       Int_t esource = -1;
540       if ( abs(mcpart->GetPdgCode()) != kdecayed ) esource = kMisID-4;
541       else esource =GetSource(mcpart)-4; // return for the cases kGamma=4, kPi0=5, kElse=6
542       if(esource==0|| esource==1 || esource==2 || esource==3){
543         fHist[iq][esource][icut].fPt->Fill(mcpart->Pt());
544         fHist[iq][esource][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart));
545         fHist[iq][esource][icut].fEta->Fill(mcpart->Eta());
546         return; 
547       }
548       else {
549         AliDebug(1, "e source is out of defined ranges, return\n");
550         return;
551       }
552     }
553
554     if ( abs(mcpart->GetPdgCode()) != kdecayed ) return;
555
556     Int_t iLabel = mcpart->GetFirstMother(); 
557     if (iLabel<0){
558       AliDebug(1, "Stack label is negative, return\n");
559       return; 
560     }
561
562     AliMCParticle *mctrack = NULL;
563     if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return; 
564     TParticle *partMother = mctrack->Particle();
565     TParticle *partMotherCopy = partMother;
566     Int_t maPdgcode = partMother->GetPdgCode();
567     Int_t maPdgcodeCopy = maPdgcode;
568
569     // get mc primary vertex
570     /*
571     TArrayF mcPrimVtx;
572     if(fMCHeader) fMCHeader->PrimaryVertex(mcPrimVtx);
573
574     // get electron production vertex   
575     TLorentzVector ePoint;
576     mcpart->ProductionVertex(ePoint);
577
578     // calculated production vertex to primary vertex (in xy plane)
579     Float_t decayLxy = TMath::Sqrt((mcPrimVtx[0]-ePoint[0])*(mcPrimVtx[0]-ePoint[0])+(mcPrimVtx[1]-ePoint[1])*(mcPrimVtx[1]-ePoint[1]));
580     */ 
581     Float_t decayLxy = 0;
582
583     // if the mother is charmed hadron  
584     Bool_t isMotherDirectCharm = kFALSE;
585     if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) { 
586
587          for (Int_t i=0; i<fNparents; i++){
588             if (abs(maPdgcode)==fParentSelect[0][i]){
589               isFinalOpenCharm = kTRUE;
590             } 
591          }  
592          if (!isFinalOpenCharm) return ;
593
594           // iterate until you find B hadron as a mother or become top ancester 
595           for (Int_t i=1; i<fgkMaxIter; i++){
596
597              Int_t jLabel = partMother->GetFirstMother(); 
598              if (jLabel == -1){
599                isMotherDirectCharm = kTRUE;
600                break; // if there is no ancester
601              }
602              if (jLabel < 0){ // safety protection
603                AliDebug(1, "Stack label is negative, return\n");
604                return; 
605              }
606
607              // if there is an ancester
608              if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return; 
609              TParticle* grandMa = mctrack->Particle();
610              Int_t grandMaPDG = grandMa->GetPdgCode();
611
612              for (Int_t j=0; j<fNparents; j++){
613                 if (abs(grandMaPDG)==fParentSelect[1][j]){
614
615                   if (kquark == kCharm) return;
616                   // fill electron kinematics
617                   fHist[iq][kElectron2nd][icut].fPdgCode->Fill(mcpart->GetPdgCode());
618                   fHist[iq][kElectron2nd][icut].fPt->Fill(mcpart->Pt());
619                   fHist[iq][kElectron2nd][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart));
620                   fHist[iq][kElectron2nd][icut].fEta->Fill(mcpart->Eta());
621
622                   // fill mother hadron kinematics
623                   fHist[iq][kDeHadron][icut].fPdgCode->Fill(grandMaPDG); 
624                   fHist[iq][kDeHadron][icut].fPt->Fill(grandMa->Pt());
625                   fHist[iq][kDeHadron][icut].fY->Fill(AliHFEtools::GetRapidity(grandMa));
626                   fHist[iq][kDeHadron][icut].fEta->Fill(grandMa->Eta());
627                  
628                   // ratio between pT of electron and pT of mother B hadron 
629                   if(grandMa->Pt()) fHistComm[iq][icut].fDePtRatio->Fill(grandMa->Pt(),mcpart->Pt()/grandMa->Pt());
630
631                   // distance between electron production point and primary vertex
632                   fHistComm[iq][icut].fDeDistance->Fill(grandMa->Pt(),decayLxy);
633                   return;
634                 }
635              } 
636
637              partMother = grandMa;
638           } // end of iteration 
639     } // end of if
640     if((isMotherDirectCharm == kTRUE && kquark == kCharm) || kquark == kBeauty) {
641          for (Int_t i=0; i<fNparents; i++){
642             if (abs(maPdgcodeCopy)==fParentSelect[iq][i]){
643
644               // fill electron kinematics
645               fHist[iq][kElectron][icut].fPdgCode->Fill(mcpart->GetPdgCode());
646               fHist[iq][kElectron][icut].fPt->Fill(mcpart->Pt());
647               fHist[iq][kElectron][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart));
648               fHist[iq][kElectron][icut].fEta->Fill(mcpart->Eta());  
649
650               // fill mother hadron kinematics
651               fHist[iq][keHadron][icut].fPdgCode->Fill(maPdgcodeCopy); 
652               fHist[iq][keHadron][icut].fPt->Fill(partMotherCopy->Pt());
653               fHist[iq][keHadron][icut].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
654               fHist[iq][keHadron][icut].fEta->Fill(partMotherCopy->Eta());
655
656               // ratio between pT of electron and pT of mother B hadron 
657               if(partMotherCopy->Pt()) fHistComm[iq][icut].fePtRatio->Fill(partMotherCopy->Pt(),mcpart->Pt()/partMotherCopy->Pt());
658
659               // distance between electron production point and primary vertex
660               fHistComm[iq][icut].feDistance->Fill(partMotherCopy->Pt(),decayLxy);
661             }
662          }
663     } // end of if
664 }
665
666 //____________________________________________________________________
667 void  AliHFEmcQA::GetDecayedKine(AliAODMCParticle *mcpart, const Int_t kquark, Int_t kdecayed, Int_t icut)
668 {
669   // decay electron kinematics
670
671   if (kquark != kCharm && kquark != kBeauty) {
672     AliDebug(1, "This task is only for heavy quark QA, return\n");
673     return;
674   }
675
676   Int_t iq = kquark - kCharm;
677   Bool_t isFinalOpenCharm = kFALSE;
678
679   if(!mcpart){
680     AliDebug(1, "no mcparticle, return\n");
681     return;
682   }
683
684   if ( abs(mcpart->GetPdgCode()) != kdecayed ) return;
685
686   // mother
687   Int_t iLabel = mcpart->GetMother();
688   if (iLabel<0){
689     AliDebug(1, "Stack label is negative, return\n");
690     return;
691   }
692
693   AliAODMCParticle *partMother = (AliAODMCParticle*)fMCArray->At(iLabel);
694   AliAODMCParticle *partMotherCopy = partMother;
695   Int_t maPdgcode = partMother->GetPdgCode();
696   Int_t maPdgcodeCopy = maPdgcode;
697
698   Bool_t isMotherDirectCharm = kFALSE;
699   if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
700
701     for (Int_t i=0; i<fNparents; i++){
702        if (abs(maPdgcode)==fParentSelect[0][i]){
703          isFinalOpenCharm = kTRUE;
704        }
705     } 
706     if (!isFinalOpenCharm) return;
707
708     for (Int_t i=1; i<fgkMaxIter; i++){
709
710        Int_t jLabel = partMother->GetMother();
711        if (jLabel == -1){
712          isMotherDirectCharm = kTRUE;
713          break; // if there is no ancester
714        }
715        if (jLabel < 0){ // safety protection
716          AliDebug(1, "Stack label is negative, return\n");
717          return;
718        }
719
720        // if there is an ancester
721        AliAODMCParticle* grandMa = (AliAODMCParticle*)fMCArray->At(jLabel);
722        Int_t grandMaPDG = grandMa->GetPdgCode();
723
724        for (Int_t j=0; j<fNparents; j++){
725           if (abs(grandMaPDG)==fParentSelect[1][j]){
726
727             if (kquark == kCharm) return;
728             // fill electron kinematics
729             fHist[iq][kElectron2nd][icut].fPdgCode->Fill(mcpart->GetPdgCode());
730             fHist[iq][kElectron2nd][icut].fPt->Fill(mcpart->Pt());
731             fHist[iq][kElectron2nd][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart));
732             fHist[iq][kElectron2nd][icut].fEta->Fill(mcpart->Eta());
733
734             // fill mother hadron kinematics
735             fHist[iq][kDeHadron][icut].fPdgCode->Fill(grandMaPDG);
736             fHist[iq][kDeHadron][icut].fPt->Fill(grandMa->Pt());
737             fHist[iq][kDeHadron][icut].fY->Fill(AliHFEtools::GetRapidity(grandMa));
738             fHist[iq][kDeHadron][icut].fEta->Fill(grandMa->Eta());
739
740             return;
741           }
742        }
743
744        partMother = grandMa;
745     } // end of iteration 
746   } // end of if
747   if ((isMotherDirectCharm == kTRUE && kquark == kCharm) || kquark == kBeauty) {
748     for (Int_t i=0; i<fNparents; i++){
749        if (abs(maPdgcodeCopy)==fParentSelect[iq][i]){
750
751          // fill electron kinematics
752          fHist[iq][kElectron][icut].fPdgCode->Fill(mcpart->GetPdgCode());
753          fHist[iq][kElectron][icut].fPt->Fill(mcpart->Pt());
754          fHist[iq][kElectron][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart));
755          fHist[iq][kElectron][icut].fEta->Fill(mcpart->Eta());
756
757          // fill mother hadron kinematics
758          fHist[iq][keHadron][icut].fPdgCode->Fill(maPdgcodeCopy);
759          fHist[iq][keHadron][icut].fPt->Fill(partMotherCopy->Pt());
760          fHist[iq][keHadron][icut].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
761          fHist[iq][keHadron][icut].fEta->Fill(partMotherCopy->Eta());
762
763        }
764     }
765   } // end of if
766
767 }
768
769 //__________________________________________
770 void AliHFEmcQA::IdentifyMother(Int_t motherlabel, Int_t &motherpdg, Int_t &grandmotherlabel)
771 {
772        // find mother pdg code and label 
773
774        if (motherlabel < 0) { 
775          AliDebug(1, "Stack label is negative, return\n");
776          return; 
777        }
778        AliMCParticle *mctrack = NULL;
779        if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(motherlabel))))) return; 
780        TParticle *heavysMother = mctrack->Particle();
781        motherpdg = heavysMother->GetPdgCode();
782        grandmotherlabel = heavysMother->GetFirstMother();
783        AliDebug(1,Form("ancestor pdg code= %d\n",motherpdg));
784 }
785
786 //__________________________________________
787 void AliHFEmcQA::HardScattering(const Int_t kquark, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
788 {
789        // mothertype -1 means this heavy quark coming from hard vertex
790
791        AliMCParticle *mctrack1 = NULL;
792        AliMCParticle *mctrack2 = NULL;
793        if(!(mctrack1 = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(4))))) return; 
794        if(!(mctrack2 = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(5))))) return; 
795        TParticle *afterinitialrad1  = mctrack1->Particle();
796        TParticle *afterinitialrad2  = mctrack2->Particle();
797            
798        motherlabel = -1;
799
800        if (abs(afterinitialrad1->GetPdgCode()) == fgkGluon && abs(afterinitialrad2->GetPdgCode()) == fgkGluon){
801          AliDebug(1,"heavy from gluon gluon pair creation!\n");
802          mothertype = -1;
803          motherID = fgkGluon;
804        }
805        else if (abs(afterinitialrad1->GetPdgCode()) == kquark || abs(afterinitialrad2->GetPdgCode()) == kquark){ // one from Q and the other from g
806          AliDebug(1,"heavy from flavor exitation!\n");
807          mothertype = -1;
808          motherID = kquark;
809        }
810        else if  (abs(afterinitialrad1->GetPdgCode()) == abs(afterinitialrad2->GetPdgCode())){
811          AliDebug(1,"heavy from q-qbar pair creation!\n");
812          mothertype = -1;
813          motherID = 1;
814        }
815        else {
816          AliDebug(1,"something strange!\n");
817          mothertype = -999;
818          motherlabel = -999;
819          motherID = -999;
820        }
821 }
822
823 //__________________________________________
824 Bool_t AliHFEmcQA::IsFromInitialShower(Int_t inputmotherlabel, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
825 {
826        // mothertype -2 means this heavy quark coming from initial state 
827
828        AliMCParticle *mctrack = NULL;
829        if (inputmotherlabel==2 || inputmotherlabel==3){ // mother exist before initial state radiation
830          if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(inputmotherlabel))))) return kFALSE; 
831          TParticle *heavysMother = mctrack->Particle();
832          motherID = heavysMother->GetPdgCode(); 
833          mothertype = -2; // there is mother before initial state radiation
834          motherlabel = inputmotherlabel;
835          AliDebug(1,"initial parton shower! \n");
836
837          return kTRUE;
838        }
839
840        return kFALSE;
841 }
842
843 //__________________________________________
844 Bool_t AliHFEmcQA::IsFromFinalParton(Int_t inputmotherlabel, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
845 {
846        // mothertype 2 means this heavy quark coming from final state 
847
848        AliMCParticle *mctrack = NULL;
849        if (inputmotherlabel > 5){ // mother exist after hard scattering
850          if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(inputmotherlabel))))) return kFALSE; 
851          TParticle *heavysMother = mctrack->Particle();
852          motherID = heavysMother->GetPdgCode(); 
853          mothertype = 2; // 
854          motherlabel = inputmotherlabel;
855          AliDebug(1,Form("heavy quark from %d after hard scattering! \n",motherID));
856
857          return kTRUE;
858        }
859        return kFALSE;
860 }
861
862 //__________________________________________
863 void AliHFEmcQA::ReportStrangeness(Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
864 {
865       // mark strange behavior  
866
867        mothertype = -888;
868        motherlabel = -888;
869        motherID = -888;
870        AliDebug(1,"something strange!\n");
871 }
872
873 //__________________________________________
874 Int_t AliHFEmcQA::GetSource(AliAODMCParticle * const mcpart)
875 {        
876   // decay particle's origin 
877
878   //if ( abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return -1;
879        
880   Int_t origin = -1;
881   Bool_t isFinalOpenCharm = kFALSE;
882
883   if(!mcpart){
884     AliDebug(1, "Stack label is negative or no mcparticle, return\n");
885     return -1;
886   }
887
888   // mother
889   Int_t iLabel = mcpart->GetMother();
890   if (iLabel<0){
891     AliDebug(1, "Stack label is negative, return\n");
892     return -1;
893   } 
894        
895   AliAODMCParticle *partMother = (AliAODMCParticle*)fMCArray->At(iLabel);
896   Int_t maPdgcode = partMother->GetPdgCode();
897   
898   // if the mother is charmed hadron  
899   if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
900     
901     for (Int_t i=0; i<fNparents; i++){
902        if (abs(maPdgcode)==fParentSelect[0][i]){
903          isFinalOpenCharm = kTRUE;
904        }
905     }
906     if (!isFinalOpenCharm) return -1;
907
908     // iterate until you find B hadron as a mother or become top ancester 
909     for (Int_t i=1; i<fgkMaxIter; i++){
910
911        Int_t jLabel = partMother->GetMother();
912        if (jLabel == -1){
913          origin = kDirectCharm;
914          return origin;
915        }
916        if (jLabel < 0){ // safety protection
917          AliDebug(1, "Stack label is negative, return\n");
918          return -1;
919        }
920
921        // if there is an ancester
922        AliAODMCParticle* grandMa = (AliAODMCParticle*)fMCArray->At(jLabel);
923        Int_t grandMaPDG = grandMa->GetPdgCode();
924
925        for (Int_t j=0; j<fNparents; j++){
926           if (abs(grandMaPDG)==fParentSelect[1][j]){
927             origin = kBeautyCharm;
928             return origin;
929           }
930        }
931
932        partMother = grandMa;
933     } // end of iteration 
934   } // end of if
935   else if ( int(abs(maPdgcode)/100.) == kBeauty || int(abs(maPdgcode)/1000.) == kBeauty ) {
936     for (Int_t i=0; i<fNparents; i++){
937        if (abs(maPdgcode)==fParentSelect[1][i]){
938          origin = kDirectBeauty;
939          return origin;
940        }
941     }
942   } // end of if
943   else if ( abs(maPdgcode) == 22 ) {
944     origin = kGamma;
945     return origin;
946   } // end of if
947   else if ( abs(maPdgcode) == 111 ) {
948     origin = kPi0;
949     return origin;
950   } // end of if
951
952   return origin;
953 }
954
955 //__________________________________________
956 Int_t AliHFEmcQA::GetSource(TParticle * const mcpart)
957 {
958   // decay particle's origin 
959
960   //if ( abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return -1;
961
962   Int_t origin = -1;
963   Bool_t isFinalOpenCharm = kFALSE;
964
965   if(!mcpart){
966     AliDebug(1, "no mcparticle, return\n");
967     return -1;
968   }
969
970   Int_t iLabel = mcpart->GetFirstMother();
971   if (iLabel<0){
972     AliDebug(1, "Stack label is negative, return\n");
973     return -1;
974   }
975
976   AliMCParticle *mctrack = NULL;
977   if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return -1; 
978   TParticle *partMother = mctrack->Particle();
979   Int_t maPdgcode = partMother->GetPdgCode();
980
981    // if the mother is charmed hadron  
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 -1;
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           origin = kDirectCharm;
997           return origin;
998         }
999         if (jLabel < 0){ // safety protection
1000           AliDebug(1, "Stack label is negative, return\n");
1001           return -1;
1002         }
1003
1004         // if there is an ancester
1005         if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return -1; 
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              origin = kBeautyCharm;
1012              return origin;
1013            }
1014         }
1015
1016         partMother = grandMa;
1017      } // end of iteration 
1018    } // end of if
1019    else if ( int(abs(maPdgcode)/100.) == kBeauty || int(abs(maPdgcode)/1000.) == kBeauty ) {
1020      for (Int_t i=0; i<fNparents; i++){
1021         if (abs(maPdgcode)==fParentSelect[1][i]){
1022           origin = kDirectBeauty;
1023           return origin;
1024         }
1025      }
1026    } // end of if
1027    else if ( abs(maPdgcode) == 22 ) {
1028      origin = kGamma;
1029      return origin;
1030    } // end of if
1031    else if ( abs(maPdgcode) == 111 ) {
1032      origin = kPi0;
1033      return origin;
1034    } // end of if
1035    else origin = kElse;
1036
1037    return origin;
1038 }
1039
1040 //__________________________________________
1041 Int_t AliHFEmcQA::GetElecSource(TParticle * const mcpart)
1042 {
1043   // decay particle's origin 
1044
1045   if ( abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return kMisID;
1046
1047   Int_t origin = -1;
1048   Bool_t isFinalOpenCharm = kFALSE;
1049
1050   if(!mcpart){
1051     AliDebug(1, "no mcparticle, return\n");
1052     return -1;
1053   }
1054
1055   Int_t iLabel = mcpart->GetFirstMother();
1056   if (iLabel<0){
1057     AliDebug(1, "Stack label is negative, return\n");
1058     return -1;
1059   }
1060
1061   AliMCParticle *mctrack = NULL;
1062   if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return -1; 
1063   TParticle *partMother = mctrack->Particle();
1064   Int_t maPdgcode = partMother->GetPdgCode();
1065
1066    // if the mother is charmed hadron  
1067    if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
1068
1069      for (Int_t i=0; i<fNparents; i++){
1070         if (abs(maPdgcode)==fParentSelect[0][i]){
1071           isFinalOpenCharm = kTRUE;
1072         }
1073      }
1074      if (!isFinalOpenCharm) return -1;
1075
1076      // iterate until you find B hadron as a mother or become top ancester 
1077      for (Int_t i=1; i<fgkMaxIter; i++){
1078
1079         Int_t jLabel = partMother->GetFirstMother();
1080         if (jLabel == -1){
1081           origin = kDirectCharm;
1082           return origin;
1083         }
1084         if (jLabel < 0){ // safety protection
1085           AliDebug(1, "Stack label is negative, return\n");
1086           return -1;
1087         }
1088
1089         // if there is an ancester
1090         if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return -1; 
1091         TParticle* grandMa = mctrack->Particle();
1092         Int_t grandMaPDG = grandMa->GetPdgCode();
1093
1094         for (Int_t j=0; j<fNparents; j++){
1095            if (abs(grandMaPDG)==fParentSelect[1][j]){
1096              origin = kBeautyCharm;
1097              return origin;
1098            }
1099         }
1100
1101         partMother = grandMa;
1102      } // end of iteration 
1103    } // end of if
1104    else if ( int(abs(maPdgcode)/100.) == kBeauty || int(abs(maPdgcode)/1000.) == kBeauty ) {
1105      for (Int_t i=0; i<fNparents; i++){
1106         if (abs(maPdgcode)==fParentSelect[1][i]){
1107           origin = kDirectBeauty;
1108           return origin;
1109         }
1110      }
1111    } // end of if
1112    else if ( abs(maPdgcode) == 22 ) {
1113      origin = kGamma;
1114      return origin;
1115    } // end of if
1116    else if ( abs(maPdgcode) == 111 ) {
1117      origin = kPi0;
1118      return origin;
1119    } // end of if
1120    else if ( abs(maPdgcode) == 221 ) {
1121      origin = kEta;
1122      return origin;
1123    } // end of if
1124    else if ( abs(maPdgcode) == 223 ) {
1125      origin = kOmega;
1126      return origin;
1127    } // end of if
1128    else if ( abs(maPdgcode) == 333 ) {
1129      origin = kPhi;
1130      return origin;
1131    } // end of if
1132    else if ( abs(maPdgcode) == 331 ) {
1133      origin = kEtaPrime;
1134      return origin;
1135    } // end of if
1136    else if ( abs(maPdgcode) == 113 ) {
1137      origin = kRho0;
1138      return origin;
1139    } // end of if
1140    else origin = kElse;
1141
1142    return origin;
1143 }
1144
1145 //__________________________________________
1146 void AliHFEmcQA::AliHists::FillList(TList *l) const {
1147   //
1148   // Fill Histos into a list for output
1149   //
1150   if(fPdgCode) l->Add(fPdgCode);
1151   if(fPt) l->Add(fPt);
1152   if(fY) l->Add(fY);
1153   if(fEta) l->Add(fEta);
1154 }
1155
1156 //__________________________________________
1157 void AliHFEmcQA::AliHistsComm::FillList(TList *l) const { 
1158   //
1159   // Fill Histos into a list for output
1160   //
1161   if(fNq) l->Add(fNq);
1162   if(fProcessID) l->Add(fProcessID);
1163   if(fePtRatio) l->Add(fePtRatio);
1164   if(fDePtRatio) l->Add(fDePtRatio);
1165   if(feDistance) l->Add(feDistance);
1166   if(fDeDistance) l->Add(fDeDistance);
1167 }
1168