Since it contains fixes of coding rule violations, all classes are involved. Further...
[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  *                                                                        *
17  * QA class of Heavy Flavor quark and fragmeted/decayed particles         *
18  * -Check kinematics of Heavy Quarks/hadrons, and decayed leptons         *
19  *    pT, rapidity                                                        *
20  *    decay lepton kinematics w/wo acceptance                             *
21  *    heavy hadron decay length, electron pT fraction carried from decay  *
22  * -Check yield of Heavy Quarks/hadrons                                   *
23  *    Number of produced heavy quark                                      *
24  *    Number of produced hadron of given pdg code                         *
25  *                                                                        *
26  *                                                                        *
27  * Authors:                                                               *
28  *   MinJung Kweon <minjung@physi.uni-heidelberg.de>                      *
29  *                                                                        *
30  **************************************************************************/
31
32
33 #include <iostream>
34 #include <TH2F.h>
35 #include <TCanvas.h>
36 #include <TFile.h>
37 #include <TH1F.h>
38 #include <TH2F.h>
39 #include <TCut.h>
40 #include <TRandom.h>
41 #include <TDatabasePDG.h>
42 #include <TROOT.h>
43 #include <TParticle.h>
44
45 #include <AliLog.h>
46 #include <AliESDEvent.h>
47 #include <AliESDtrack.h>
48 #include <AliESDInputHandler.h>
49 #include <AliESDVertex.h>
50 #include <AliStack.h>
51
52 #include "AliHFEmcQA.h"
53
54 const Int_t AliHFEmcQA::fgkGluon=21;
55 const Int_t AliHFEmcQA::fgkMaxGener=10;
56 const Int_t AliHFEmcQA::fgkMaxIter=100;
57 const Int_t AliHFEmcQA::fgkqType = 7;    // number of species waiting for QA done
58
59 ClassImp(AliHFEmcQA)
60
61 //_______________________________________________________________________________________________
62 AliHFEmcQA::AliHFEmcQA() : 
63         fStack(0x0) 
64         ,fNparents(0) 
65 {
66         // Default constructor
67         
68 }
69
70 //_______________________________________________________________________________________________
71 AliHFEmcQA::AliHFEmcQA(const AliHFEmcQA&p):
72         TObject(p)
73         ,fStack(0x0) 
74         ,fNparents(p.fNparents) 
75 {
76         // Copy constructor
77 }
78
79 //_______________________________________________________________________________________________
80 AliHFEmcQA&
81 AliHFEmcQA::operator=(const AliHFEmcQA &)
82 {
83   // Assignment operator
84
85   AliInfo("Not yet implemented.");
86   return *this;
87 }
88
89 //_______________________________________________________________________________________________
90 AliHFEmcQA::~AliHFEmcQA()
91 {
92         // Destructor
93
94         AliInfo("Analysis Done.");
95 }
96
97 //_______________________________________________________________________________________________
98 void AliHFEmcQA::PostAnalyze()
99 {
100 }
101
102 //__________________________________________
103 void AliHFEmcQA::CreateHistograms(const Int_t kquark, Int_t icut, TString hnopt) 
104 {
105   // create histograms
106
107   if (kquark != kCharm && kquark != kBeauty) {
108     AliDebug(1, "This task is only for heavy quark QA, return\n");
109     return; 
110   }
111   Int_t iq = kquark - kCharm; 
112
113   TString kqTypeLabel[fgkqType];
114   if (kquark == kCharm){
115     kqTypeLabel[kQuark]="c";
116     kqTypeLabel[kantiQuark]="cbar";
117     kqTypeLabel[kHadron]="cHadron";
118     kqTypeLabel[keHadron]="ceHadron";
119     kqTypeLabel[kDeHadron]="nullHadron";
120     kqTypeLabel[kElectron]="ce";
121     kqTypeLabel[kElectron2nd]="nulle";
122   } else if (kquark == kBeauty){
123     kqTypeLabel[kQuark]="b";
124     kqTypeLabel[kantiQuark]="bbar";
125     kqTypeLabel[kHadron]="bHadron";
126     kqTypeLabel[keHadron]="beHadron";
127     kqTypeLabel[kDeHadron]="bDeHadron";
128     kqTypeLabel[kElectron]="be";
129     kqTypeLabel[kElectron2nd]="bce";
130   }
131
132
133   TString hname; 
134   for (Int_t iqType = 0; iqType < fgkqType; iqType++ ){
135      if (iqType < keHadron && icut > 0) continue; // don't duplicate histogram for quark and hadron
136      hname = hnopt+"PdgCode_"+kqTypeLabel[iqType];
137      fHist[iq][iqType][icut].fPdgCode = new TH1F(hname,hname,20001,-10000.5,10000.5);
138      hname = hnopt+"Pt_"+kqTypeLabel[iqType];
139      fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",150,0,30);
140      hname = hnopt+"Y_"+kqTypeLabel[iqType];
141      fHist[iq][iqType][icut].fY = new TH1F(hname,hname,150,-7.5,7.5);
142      hname = hnopt+"Eta_"+kqTypeLabel[iqType];
143      fHist[iq][iqType][icut].fEta = new TH1F(hname,hname,150,-7.5,7.5);
144   }
145
146   if (icut == 0){ 
147     hname = hnopt+"Nq_"+kqTypeLabel[kQuark];
148     fHistComm[iq][icut].fNq = new TH1F(hname,hname,10,-0.5,9.5);
149     hname = hnopt+"ProcessID_"+kqTypeLabel[kQuark];
150     fHistComm[iq][icut].fProcessID = new TH1F(hname,hname,21,-10.5,10.5);
151   }
152   hname = hnopt+"ePtRatio_"+kqTypeLabel[kQuark];
153   fHistComm[iq][icut].fePtRatio = new TH2F(hname,hname+";p_{T} (GeV/c);momentum fraction",100,0,30,100,0,1);
154   hname = hnopt+"DePtRatio_"+kqTypeLabel[kQuark];
155   fHistComm[iq][icut].fDePtRatio = new TH2F(hname,hname+";p_{T} (GeV/c);momentum fraction",100,0,30,100,0,1);
156   hname = hnopt+"eDistance_"+kqTypeLabel[kQuark];
157   fHistComm[iq][icut].feDistance= new TH2F(hname,hname+";p_{T} (GeV/c);distance (cm)",100,0,30,200,0,2);
158   hname = hnopt+"DeDistance_"+kqTypeLabel[kQuark];
159   fHistComm[iq][icut].fDeDistance= new TH2F(hname,hname+";p_{T} (GeV/c);distance (cm)",100,0,30,200,0,2);
160
161 }
162
163 //__________________________________________
164 void AliHFEmcQA::Init()
165 {
166   // called at begining every event
167   
168   for (Int_t i=0; i<2; i++){
169      fIsHeavy[i] = 0;
170   } 
171
172   fNparents = 7;
173
174   fParentSelect[0][0] =  411;  
175   fParentSelect[0][1] =  421;
176   fParentSelect[0][2] =  431;
177   fParentSelect[0][3] = 4122;
178   fParentSelect[0][4] = 4132;
179   fParentSelect[0][5] = 4232;
180   fParentSelect[0][6] = 4332;
181
182   fParentSelect[1][0] =  511;
183   fParentSelect[1][1] =  521;
184   fParentSelect[1][2] =  531;
185   fParentSelect[1][3] = 5122;
186   fParentSelect[1][4] = 5132;
187   fParentSelect[1][5] = 5232;
188   fParentSelect[1][6] = 5332;
189
190 }
191
192 //__________________________________________
193 void AliHFEmcQA::GetQuarkKine(Int_t iTrack, const Int_t kquark) 
194 {
195   // get heavy quark kinematics
196
197     if (kquark != kCharm && kquark != kBeauty) {
198       AliDebug(1, "This task is only for heavy quark QA, return\n");
199       return; 
200     }
201     Int_t iq = kquark - kCharm; 
202
203
204     if (iTrack < 0) { 
205       AliDebug(1, "Stack label is negative, return\n");
206       return; 
207     }
208
209     TParticle *part = fStack->Particle(iTrack); 
210     Int_t partPdgcode = TMath::Abs(part->GetPdgCode());
211
212     // select heavy hadron or not fragmented heavy quark 
213     if ( int(partPdgcode/100.)==kquark || int(partPdgcode/1000.)==kquark || (partPdgcode==kquark && (part->GetNDaughters()==0 && iTrack>5)) ){ 
214
215       TParticle *partMother;
216       Int_t iLabel;
217
218       if (partPdgcode == kquark){ // in case of not fragmented heavy quark  
219         partMother = part; 
220         iLabel = iTrack;
221       } else{ // in case of heavy hadron, start to search for mother heavy parton 
222         iLabel = part->GetFirstMother(); 
223         if (iLabel>-1) { partMother = fStack->Particle(iLabel); }
224         else {
225           AliDebug(1, "Stack label is negative, return\n");
226           return; 
227         }
228       }
229
230       // heavy parton selection as a mother of heavy hadron 
231       // if the heavy particle comes from string which is denoted as particle status 12|12|12...12|11,[PYTHIA p.60]
232       // in this case, the mother of heavy particle can be one of the fragmented parton of the string
233       // should I make a condition that partMother should be quark or diquark? -> not necessary
234       if ( abs(partMother->GetPdgCode()) == kquark || (partMother->GetStatusCode() == 11) ){
235       //if ( abs(partMother->GetPdgCode()) == kquark || (partMother->GetStatusCode() == 11 || partMother->GetStatusCode() == 12) ){
236
237         if ( abs(partMother->GetPdgCode()) != kquark ){
238           // search fragmented partons in the same string
239           Bool_t isSameString = kTRUE; 
240           for (Int_t i=1; i<fgkMaxIter; i++){
241              iLabel = iLabel - 1;
242              if (iLabel>-1) { partMother = fStack->Particle(iLabel); }
243              else {
244                AliDebug(1, "Stack label is negative, return\n");
245                return; 
246              }
247              if ( abs(partMother->GetPdgCode()) == kquark ) break;
248              if ( partMother->GetStatusCode() != 12 ) isSameString = kFALSE;
249              if (!isSameString) return; 
250           }
251         }
252         AliDebug(1, "Can not find heavy parton of this heavy hadron in the string, return\n");
253         if (abs(partMother->GetPdgCode()) != kquark) return; 
254
255         if (fIsHeavy[iq] >= 50) return;  
256         fHeavyQuark[fIsHeavy[iq]] = partMother;
257         fIsHeavy[iq]++;
258
259         // fill kinematics for heavy parton
260         if (partMother->GetPdgCode() > 0) { // quark
261           fHist[iq][kQuark][0].fPdgCode->Fill(partMother->GetPdgCode());
262           fHist[iq][kQuark][0].fPt->Fill(partMother->Pt());
263           fHist[iq][kQuark][0].fY->Fill(GetRapidity(partMother));
264           fHist[iq][kQuark][0].fEta->Fill(partMother->Eta());
265         } else{ // antiquark
266           fHist[iq][kantiQuark][0].fPdgCode->Fill(partMother->GetPdgCode());
267           fHist[iq][kantiQuark][0].fPt->Fill(partMother->Pt());
268           fHist[iq][kantiQuark][0].fY->Fill(GetRapidity(partMother));
269           fHist[iq][kantiQuark][0].fEta->Fill(partMother->Eta());
270         }
271
272       } // end of heavy parton slection loop 
273
274     } // end of heavy hadron or quark selection
275
276 }
277
278 //__________________________________________
279 void AliHFEmcQA::EndOfEventAna(const Int_t kquark)
280 {
281   // end of event analysis
282
283   if (kquark != kCharm && kquark != kBeauty) {
284     AliDebug(1, "This task is only for heavy quark QA, return\n");
285     return; 
286   }
287   Int_t iq = kquark - kCharm; 
288
289
290   // # of heavy quark per event
291   AliDebug(1,Form("Number of heavy quark in this event = %d \n",fIsHeavy[iq]));
292   fHistComm[iq][0].fNq->Fill(fIsHeavy[iq]);
293
294   Int_t motherID[fgkMaxGener];
295   Int_t motherType[fgkMaxGener];
296   Int_t motherLabel[fgkMaxGener];
297   Int_t ancestorPdg[fgkMaxGener];
298   Int_t ancestorLabel[fgkMaxGener];
299
300   for (Int_t i = 0; i < fgkMaxGener; i++){ // initialization
301      motherID[i] = 0;
302      motherType[i] = 0;
303      motherLabel[i] = 0;
304      ancestorPdg[i] = 0;
305      ancestorLabel[i] = 0;
306   }
307
308   // check history of found heavy quarks
309   for (Int_t i = 0; i < fIsHeavy[iq]; i++){
310
311      ancestorLabel[0] = i;
312      ancestorPdg[0] = fHeavyQuark[i]->GetPdgCode(); 
313      ancestorLabel[1] = fHeavyQuark[i]->GetFirstMother(); 
314
315      AliDebug(1,Form("pdg code= %d\n",ancestorPdg[0]));
316      AliDebug(1,Form("ancestor label= %d\n",ancestorLabel[1]));
317
318      Int_t ig = 1;
319      while (ancestorLabel[ig] != -1){
320           // in case there is mother, get mother's pdg code and grandmother's label
321           IdentifyMother(ancestorLabel[ig], ancestorPdg[ig], ancestorLabel[ig+1]); 
322           // if mother is still heavy, find again mother's ancestor
323           if (ancestorPdg[ig-1] == ancestorPdg[ig]) {
324             ig++;
325             continue; // if it is from same heavy
326           }
327           // if the heavy's mother is not heavy, check the mother's label to know if it comes from inital or final parton shower
328           if (IsFromInitialShower(ancestorLabel[ig],motherID[i],motherType[i],motherLabel[i])) break;
329           if (IsFromFinalParton(ancestorLabel[ig],motherID[i],motherType[i],motherLabel[i])) break;
330           // if it is not the above case, something is strange
331           ReportStrangeness(motherID[i],motherType[i],motherLabel[i]);
332           break;
333      } 
334      if (ancestorLabel[ig] == -1){ // from hard scattering
335        HardScattering(kquark, motherID[i],motherType[i], motherLabel[i]);
336      }
337
338   } // end of found heavy quark loop
339
340
341   // check process type
342   Int_t processID = 0;
343   for (Int_t i = 0; i < fIsHeavy[iq]; i++){
344      AliDebug(1,Form("Mother ID= %d type= %d label= %d\n",motherID[i],motherType[i],motherLabel[i]));
345   }
346
347
348   Int_t nheavypair = Int_t(fIsHeavy[iq]/2.); 
349   for (Int_t ipair = 0; ipair < nheavypair; ipair++){
350
351      Int_t id1 = ipair*2;
352      Int_t id2 = ipair*2 + 1;
353
354      if (motherType[id1] == 2 && motherType[id2] == 2){
355        if (motherLabel[id1] == motherLabel[id2]) processID = kGluonSplitting; // gluon spliting
356        else processID = -9;
357      }
358      else if (motherType[id1] == -1 && motherType[id2] == -1) {
359        if (motherLabel[id1] == -1 && motherLabel[id2] == -1) {
360          if (motherID[id1] == fgkGluon) processID = kPairCreationFromg; // gluon fusion
361          else processID = kPairCreationFromq; // q-qbar pair creation
362        }
363        else processID = -8;
364      }
365      else if (motherType[id1] == -1 || motherType[id2] == -1) {
366        if ((motherLabel[id1] == -1 || motherLabel[id2] == -1) && (motherLabel[id1]*motherLabel[id2] == -2 || motherLabel[id1]*motherLabel[id2] == -3)) {
367          if(motherID[id1]*motherID[id2] == kquark*fgkGluon) processID = kFlavourExitation; // flavour exitation 
368          else processID = kLightQuarkShower;
369        }
370        else processID = -7;
371      }
372      else if (motherType[id1] == -2 || motherType[id2] == -2) {
373        if (motherLabel[id1] == motherLabel[id2]) processID = kInitialPartonShower; // initial parton shower
374        else processID = -6;
375        
376      }
377      else processID = -5;
378
379      if (nheavypair >1) AliDebug(1,Form("Multi pair found : process ID = %d\n",processID));
380      else fHistComm[iq][0].fProcessID->Fill(processID);
381      AliDebug(1,Form("Process ID = %d\n",processID));
382   } // end of # heavy quark pair loop
383
384 }
385
386 //__________________________________________
387 void AliHFEmcQA::GetHadronKine(Int_t iTrack, const Int_t kquark)
388 {
389     // decay electron kinematics
390
391     if (kquark != kCharm && kquark != kBeauty) {
392       AliDebug(1, "This task is only for heavy quark QA, return\n");
393       return;
394     }
395     Int_t iq = kquark - kCharm;
396
397     if (iTrack < 0) {
398       AliDebug(1, "Stack label is negative, return\n");
399       return;
400     }
401
402     TParticle* mcpart = fStack->Particle(iTrack);
403
404     Int_t iLabel = mcpart->GetFirstMother();
405     if (iLabel<0){
406       AliDebug(1, "Stack label is negative, return\n");
407       return;
408     }
409
410     TParticle *partCopy = mcpart;
411     Int_t pdgcode = mcpart->GetPdgCode();
412     Int_t pdgcodeCopy = pdgcode;
413
414     // if the mother is charmed hadron  
415     Bool_t isDirectCharm = kFALSE;
416     if ( int(abs(pdgcode)/100.) == kCharm || int(abs(pdgcode)/1000.) == kCharm ) {
417
418           // iterate until you find B hadron as a mother or become top ancester 
419           for (Int_t i=1; i<fgkMaxIter; i++){
420
421              Int_t jLabel = mcpart->GetFirstMother();
422              if (jLabel == -1){
423                isDirectCharm = kTRUE;
424                break; // if there is no ancester
425              }
426              if (jLabel < 0){ // safety protection
427                AliDebug(1, "Stack label is negative, return\n");
428                return;
429              }
430              // if there is an ancester
431              TParticle* mother = fStack->Particle(jLabel);
432              Int_t motherPDG = mother->GetPdgCode();
433     
434              for (Int_t j=0; j<fNparents; j++){
435                 if (abs(motherPDG)==fParentSelect[1][j]) return; // return if this hadron is originated from b
436              }
437
438              mcpart = mother;
439           } // end of iteration 
440     } // end of if
441     if((isDirectCharm == kTRUE && kquark == kCharm) || kquark == kBeauty) {
442          for (Int_t i=0; i<fNparents; i++){
443             if (abs(pdgcodeCopy)==fParentSelect[iq][i]){
444
445               // fill hadron kinematics
446               fHist[iq][kHadron][0].fPdgCode->Fill(pdgcodeCopy);
447               fHist[iq][kHadron][0].fPt->Fill(partCopy->Pt());
448               fHist[iq][kHadron][0].fY->Fill(GetRapidity(partCopy));
449               fHist[iq][kHadron][0].fEta->Fill(partCopy->Eta());
450
451             }
452          }
453     } // end of if
454 }
455
456 //__________________________________________
457 void AliHFEmcQA::GetDecayedKine(Int_t iTrack, const Int_t kquark, Int_t kdecayed, Int_t icut, Bool_t isbarrel) 
458 {
459     // decay electron kinematics
460     
461     if (kquark != kCharm && kquark != kBeauty) {
462       AliDebug(1, "This task is only for heavy quark QA, return\n");
463       return; 
464     }
465     Int_t iq = kquark - kCharm; 
466
467     if (iTrack < 0) { 
468       AliDebug(1, "Stack label is negative, return\n");
469       return; 
470     }
471
472     TParticle* mcpart = fStack->Particle(iTrack);
473
474     if ( abs(mcpart->GetPdgCode()) != kdecayed ) return;
475     if ( isbarrel && TMath::Abs(mcpart->Eta()) > 0.9 ) return;
476
477     Int_t iLabel = mcpart->GetFirstMother(); 
478     if (iLabel<0){
479       AliDebug(1, "Stack label is negative, return\n");
480       return; 
481     }
482
483     TParticle *partMother = fStack->Particle(iLabel);
484     TParticle *partMotherCopy = partMother;
485     Int_t maPdgcode = partMother->GetPdgCode();
486     Int_t maPdgcodeCopy = maPdgcode;
487
488     // get electron production vertex   
489     TLorentzVector ePoint;
490     mcpart->ProductionVertex(ePoint);
491
492     // if the mother is charmed hadron  
493     Bool_t isMotherDirectCharm = kFALSE;
494     if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) { 
495
496           // iterate until you find B hadron as a mother or become top ancester 
497           for (Int_t i=1; i<fgkMaxIter; i++){
498
499              Int_t jLabel = partMother->GetFirstMother(); 
500              if (jLabel == -1){
501                isMotherDirectCharm = kTRUE;
502                break; // if there is no ancester
503              }
504              if (jLabel < 0){ // safety protection
505                AliDebug(1, "Stack label is negative, return\n");
506                return; 
507              }
508
509              // if there is an ancester
510              TParticle* grandMa = fStack->Particle(jLabel);
511              Int_t grandMaPDG = grandMa->GetPdgCode();
512
513              for (Int_t j=0; j<fNparents; j++){
514                 if (abs(grandMaPDG)==fParentSelect[1][j]){
515
516                   if (kquark == kCharm) return;
517                   // fill electron kinematics
518                   fHist[iq][kElectron2nd][icut].fPdgCode->Fill(mcpart->GetPdgCode());
519                   fHist[iq][kElectron2nd][icut].fPt->Fill(mcpart->Pt());
520                   fHist[iq][kElectron2nd][icut].fY->Fill(GetRapidity(mcpart));
521                   fHist[iq][kElectron2nd][icut].fEta->Fill(mcpart->Eta());
522
523                   // fill mother hadron kinematics
524                   fHist[iq][kDeHadron][icut].fPdgCode->Fill(grandMaPDG); 
525                   fHist[iq][kDeHadron][icut].fPt->Fill(grandMa->Pt());
526                   fHist[iq][kDeHadron][icut].fY->Fill(GetRapidity(grandMa));
527                   fHist[iq][kDeHadron][icut].fEta->Fill(grandMa->Eta());
528                  
529                   // ratio between pT of electron and pT of mother B hadron 
530                   if(grandMa->Pt()) fHistComm[iq][icut].fDePtRatio->Fill(grandMa->Pt(),mcpart->Pt()/grandMa->Pt());
531
532                   // distance between electron production point and mother hadron production point
533                   TLorentzVector debPoint;
534                   grandMa->ProductionVertex(debPoint);
535                   TLorentzVector dedistance = ePoint - debPoint;
536                   fHistComm[iq][icut].fDeDistance->Fill(grandMa->Pt(),dedistance.P());
537                   return;
538                 }
539              } 
540
541              partMother = grandMa;
542           } // end of iteration 
543     } // end of if
544     if((isMotherDirectCharm == kTRUE && kquark == kCharm) || kquark == kBeauty) {
545          for (Int_t i=0; i<fNparents; i++){
546             if (abs(maPdgcodeCopy)==fParentSelect[iq][i]){
547
548               // fill electron kinematics
549               fHist[iq][kElectron][icut].fPdgCode->Fill(mcpart->GetPdgCode());
550               fHist[iq][kElectron][icut].fPt->Fill(mcpart->Pt());
551               fHist[iq][kElectron][icut].fY->Fill(GetRapidity(mcpart));
552               fHist[iq][kElectron][icut].fEta->Fill(mcpart->Eta());  
553
554               // fill mother hadron kinematics
555               fHist[iq][keHadron][icut].fPdgCode->Fill(maPdgcodeCopy); 
556               fHist[iq][keHadron][icut].fPt->Fill(partMotherCopy->Pt());
557               fHist[iq][keHadron][icut].fY->Fill(GetRapidity(partMotherCopy));
558               fHist[iq][keHadron][icut].fEta->Fill(partMotherCopy->Eta());
559
560               // ratio between pT of electron and pT of mother B hadron 
561               if(partMotherCopy->Pt()) fHistComm[iq][icut].fePtRatio->Fill(partMotherCopy->Pt(),mcpart->Pt()/partMotherCopy->Pt());
562
563               // distance between electron production point and mother hadron production point
564               TLorentzVector ebPoint;
565               partMotherCopy->ProductionVertex(ebPoint);
566               TLorentzVector edistance = ePoint - ebPoint;
567               fHistComm[iq][icut].feDistance->Fill(partMotherCopy->Pt(),edistance.P());
568             }
569          }
570     } // end of if
571 }
572
573
574 //__________________________________________
575 void AliHFEmcQA::IdentifyMother(Int_t motherlabel, Int_t &motherpdg, Int_t &grandmotherlabel)
576 {
577        // find mother pdg code and label 
578
579        if (motherlabel < 0) { 
580          AliDebug(1, "Stack label is negative, return\n");
581          return; 
582        }
583        TParticle *heavysMother = fStack->Particle(motherlabel);
584        motherpdg = heavysMother->GetPdgCode();
585        grandmotherlabel = heavysMother->GetFirstMother();
586        AliDebug(1,Form("ancestor pdg code= %d\n",motherpdg));
587 }
588
589 //__________________________________________
590 void AliHFEmcQA::HardScattering(const Int_t kquark, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
591 {
592        // mothertype -1 means this heavy quark coming from hard vertex
593
594        TParticle *afterinitialrad1  = fStack->Particle(4);
595        TParticle *afterinitialrad2  = fStack->Particle(5);
596            
597        motherlabel = -1;
598
599        if (abs(afterinitialrad1->GetPdgCode()) == fgkGluon && abs(afterinitialrad2->GetPdgCode()) == fgkGluon){
600          AliDebug(1,"heavy from gluon gluon pair creation!\n");
601          mothertype = -1;
602          motherID = fgkGluon;
603        }
604        else if (abs(afterinitialrad1->GetPdgCode()) == kquark || abs(afterinitialrad2->GetPdgCode()) == kquark){ // one from Q and the other from g
605          AliDebug(1,"heavy from flavor exitation!\n");
606          mothertype = -1;
607          motherID = kquark;
608        }
609        else if  (abs(afterinitialrad1->GetPdgCode()) == abs(afterinitialrad2->GetPdgCode())){
610          AliDebug(1,"heavy from q-qbar pair creation!\n");
611          mothertype = -1;
612          motherID = 1;
613        }
614        else {
615          AliDebug(1,"something strange!\n");
616          mothertype = -999;
617          motherlabel = -999;
618          motherID = -999;
619        }
620 }
621
622 //__________________________________________
623 Bool_t AliHFEmcQA::IsFromInitialShower(Int_t inputmotherlabel, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
624 {
625        // mothertype -2 means this heavy quark coming from initial state 
626
627        if (inputmotherlabel==2 || inputmotherlabel==3){ // mother exist before initial state radiation
628          TParticle *heavysMother = fStack->Particle(inputmotherlabel);
629          motherID = heavysMother->GetPdgCode(); 
630          mothertype = -2; // there is mother before initial state radiation
631          motherlabel = inputmotherlabel;
632          AliDebug(1,"initial parton shower! \n");
633
634          return kTRUE;
635        }
636
637        return kFALSE;
638 }
639
640 //__________________________________________
641 Bool_t AliHFEmcQA::IsFromFinalParton(Int_t inputmotherlabel, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
642 {
643        // mothertype 2 means this heavy quark coming from final state 
644
645        if (inputmotherlabel > 5){ // mother exist after hard scattering
646          TParticle *heavysMother = fStack->Particle(inputmotherlabel);
647          motherID = heavysMother->GetPdgCode(); 
648          mothertype = 2; // 
649          motherlabel = inputmotherlabel;
650          AliDebug(1,Form("heavy quark from %d after hard scattering! \n",motherID));
651
652          return kTRUE;
653        }
654        return kFALSE;
655 }
656
657 //__________________________________________
658 void AliHFEmcQA::ReportStrangeness(Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
659 {
660       // mark strange behavior  
661
662        mothertype = -888;
663        motherlabel = -888;
664        motherID = -888;
665        AliDebug(1,"something strange!\n");
666 }
667
668 //__________________________________________
669 Float_t AliHFEmcQA::GetRapidity(TParticle *part)
670 {
671       // return rapidity
672
673        Float_t rapidity;        
674        if(part->Energy() - part->Pz() == 0 || part->Energy() + part->Pz() == 0) rapidity=-999; 
675        else rapidity = 0.5*(TMath::Log((part->Energy()+part->Pz()) / (part->Energy()-part->Pz()))); 
676        return rapidity;
677 }