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