]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG4/PartCorrBase/AliMCAnalysisUtils.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
1)AliCaloPID: Posibility to recalculate PID bayesian in EMCAL
[u/mrichter/AliRoot.git] / PWG4 / PartCorrBase / AliMCAnalysisUtils.cxx
CommitLineData
7cd4e982 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/* $Id: AliMCAnalysisUtils.cxx 21839 2007-10-29 13:49:42Z gustavo $ */
16
17//_________________________________________________________________________
18// Class for analysis utils for MC data
19// stored in stack or event header.
20// Contains:
21// - method to check the origin of a given track/cluster
22// - method to obtain the generated jets
23//
24//*-- Author: Gustavo Conesa (LNF-INFN)
25//////////////////////////////////////////////////////////////////////////////
26
27
28// --- ROOT system ---
29#include <TMath.h>
30#include <TList.h>
31#include "TParticle.h"
32
33//---- ANALYSIS system ----
34#include "AliMCAnalysisUtils.h"
35#include "AliCaloTrackReader.h"
36#include "AliStack.h"
37#include "AliGenPythiaEventHeader.h"
38#include "AliAODMCParticle.h"
39
40 ClassImp(AliMCAnalysisUtils)
41
42 //________________________________________________
43 AliMCAnalysisUtils::AliMCAnalysisUtils() :
44 TObject(), fCurrentEvent(-1), fDebug(-1),
45 fJetsList(new TList), fMCGenerator("PYTHIA")
46{
47 //Ctor
48}
49
50//____________________________________________________________________________
51AliMCAnalysisUtils::AliMCAnalysisUtils(const AliMCAnalysisUtils & mcutils) :
52 TObject(mcutils), fCurrentEvent(mcutils.fCurrentEvent), fDebug(mcutils.fDebug),
53 fJetsList(mcutils.fJetsList), fMCGenerator(mcutils.fMCGenerator)
54{
55 // cpy ctor
56
57}
58
59//_________________________________________________________________________
60AliMCAnalysisUtils & AliMCAnalysisUtils::operator = (const AliMCAnalysisUtils & mcutils)
61{
62 // assignment operator
63
64 if(&mcutils == this) return *this;
65 fCurrentEvent = mcutils.fCurrentEvent ;
66 fDebug = mcutils.fDebug;
67 fJetsList = mcutils.fJetsList;
68 fMCGenerator = mcutils.fMCGenerator;
69
70 return *this;
71}
72
73//____________________________________________________________________________
74AliMCAnalysisUtils::~AliMCAnalysisUtils()
75{
76 // Remove all pointers.
77
78 if (fJetsList) {
79 fJetsList->Clear();
80 delete fJetsList ;
81 }
82}
83
84
85//_________________________________________________________________________
86Int_t AliMCAnalysisUtils::CheckOrigin(const Int_t label, AliCaloTrackReader* reader, const Int_t input = 0) {
87 //Play with the montecarlo particles if available
88 Int_t tag = 0;
89
90 //Select where the information is, ESD-galice stack or AOD mcparticles branch
91 if(reader->ReadStack()){
92 tag = CheckOriginInStack(label, reader->GetStack());
93 }
94 else if(reader->ReadAODMCParticles()){
95 tag = CheckOriginInAOD(label, reader->GetAODMCParticles(input));
96 }
97
98 return tag ;
99}
100
101//_________________________________________________________________________
102Int_t AliMCAnalysisUtils::CheckOriginInStack(const Int_t label, AliStack* stack) {
103 // Play with the MC stack if available. Tag particles depending on their origin.
104 // Do same things as in CheckOriginInAOD but different input.
105
106 if(!stack) {
107 printf("AliMCAnalysisUtils::CheckOriginInStack() - Stack is not available, check analysis settings in configuration file, STOP!!\n");
108 abort();
109 }
110 // printf("label %d, ntrack %d, nprim %d\n",label, stack->GetNtrack(), stack->GetNprimary());
111 // for(Int_t i = 0; i< stack->GetNprimary(); i++){
112 // TParticle *particle = stack->Particle(i);
113 // //particle->Print();
114 // }
115
116
117 Int_t tag = 0;
118 if(label >= 0 && label < stack->GetNtrack()){
119 //Mother
120 TParticle * mom = stack->Particle(label);
121 Int_t mPdg = TMath::Abs(mom->GetPdgCode());
122 Int_t mStatus = mom->GetStatusCode() ;
123 Int_t iParent = mom->GetFirstMother() ;
124 if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - Mother is parton %d\n",iParent);
125
126 //GrandParent
127 TParticle * parent = new TParticle ;
128 Int_t pPdg = -1;
129 Int_t pStatus =-1;
130 if(iParent > 0){
131 parent = stack->Particle(iParent);
132 pPdg = TMath::Abs(parent->GetPdgCode());
133 pStatus = parent->GetStatusCode();
134 }
135 else if(fDebug > 0 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - Parent with label %d\n",iParent);
136
137 //Check if mother is converted, if not, get the first non converted mother
138 if((mPdg == 22 || mPdg == 11) && (pPdg == 22 || pPdg == 11) && mStatus == 0){
139 SetTagBit(tag,kMCConversion);
140 //Check if the mother is photon or electron with status not stable
141 while ((pPdg == 22 || pPdg == 11) && mStatus != 1) {
142 //Mother
143 mom = stack->Particle(mom->GetFirstMother());
144 mPdg = TMath::Abs(mom->GetPdgCode());
145 mStatus = mom->GetStatusCode() ;
146 iParent = mom->GetFirstMother() ;
147 if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - Mother is parton %d\n",iParent);
148
149 //GrandParent
150 if(iParent > 0){
151 parent = stack->Particle(iParent);
152 pPdg = TMath::Abs(parent->GetPdgCode());
153 pStatus = parent->GetStatusCode();
154 }
155 }//while
156 }//mother and parent are electron or photon and have status 0
157 else if(mStatus == 0){
158 //Still a conversion but only one electron/photon generated. Just from hadrons
159 if(pPdg == 2112 || pPdg == 211 ||
160 pPdg == 321 || pPdg == 2212 || pPdg == 130 || pPdg == 13 )
161 SetTagBit(tag,kMCConversion);
162 mom = stack->Particle(mom->GetFirstMother());
163 mPdg = TMath::Abs(mom->GetPdgCode());
164 //Comment for the next lines, we do not check the parent of the hadron for the moment.
165 //iParent = mom->GetFirstMother() ;
166 //if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Mother is parton %d\n",iParent);
167
168 //GrandParent
169 //if(iParent >= 0){
170 // parent = stack->Particle(iParent);
171 // pPdg = TMath::Abs(parent->GetPdgCode());
172 //}
173 }
174 // conversion fo electrons/photons checked
175
176 //first check for typical charged particles
177 if(mPdg == 13) SetTagBit(tag,kMCMuon);
178 else if(mPdg == 211) SetTagBit(tag,kMCPion);
179 else if(mPdg == 321) SetTagBit(tag,kMCKaon);
180 else if(mPdg == 2212) SetTagBit(tag,kMCProton);
181 //check for pi0 and eta (shouldn't happen unless their decays were
182 //turned off)
183 else if(mPdg == 111) SetTagBit(tag,kMCPi0);
184 else if(mPdg == 221) SetTagBit(tag,kMCEta);
185 //Photons
186 else if(mPdg == 22){
187 SetTagBit(tag,kMCPhoton);
188 if(mStatus == 1){ //undecayed particle
189 if(fMCGenerator == "PYTHIA"){
190 if(iParent < 8 && iParent > 5) {//outgoing partons
191 if(pPdg == 22) SetTagBit(tag,kMCPrompt);
192 else SetTagBit(tag,kMCFragmentation);
193 }//Outgoing partons
194 else if(iParent <= 5) {
195 SetTagBit(tag, kMCISR); //Initial state radiation
196 }
197 else if(pStatus == 11){//Decay
198 if(pPdg == 111) SetTagBit(tag,kMCPi0Decay);
199 else if (pPdg == 221) SetTagBit(tag, kMCEtaDecay);
200 else SetTagBit(tag,kMCOtherDecay);
201 }//Decay
202 else {
203 printf("AliMCAnalysisUtils::ChecOrigingInAOD() - what is it? Mother mPdg %d, status %d \n Parent iParent %d, pPdg %d %s, status %d\n",
204 mPdg, mStatus,iParent, pPdg, parent->GetName(),pStatus);
205 SetTagBit(tag,kMCOtherDecay);//Check
206 }
207 }//PYTHIA
208
209 else if(fMCGenerator == "HERWIG"){
210 if(pStatus < 197){//Not decay
211 while(1){
212 if(parent->GetFirstMother()<=5) break;
213 iParent = parent->GetFirstMother();
214 parent=stack->Particle(iParent);
215 pStatus= parent->GetStatusCode();
216 pPdg = parent->GetPdgCode();
217 }//Look for the parton
218
219 if(iParent < 8 && iParent > 5) {
220 if(pPdg == 22) SetTagBit(tag,kMCPrompt);
221 else SetTagBit(tag,kMCFragmentation);
222 }
223 else SetTagBit(tag,kMCISR);//Initial state radiation
224 }//Not decay
225 else{//Decay
226 if(pPdg == 111) SetTagBit(tag,kMCPi0Decay);
227 else if (pPdg == 221) SetTagBit(tag,kMCEtaDecay);
228 else SetTagBit(tag,kMCOtherDecay);
229 }//Decay
230 }//HERWIG
231
232 else SetTagBit(tag,kMCUnknown);
233
234 }//Status 1 : created by event generator
235
236 else if(mStatus == 0){ // geant
237 if(pPdg == 111) SetTagBit(tag,kMCPi0Decay);
238 else if (pPdg == 221) SetTagBit(tag,kMCEtaDecay);
239 else SetTagBit(tag,kMCOtherDecay);
240 }//status 0 : geant generated
241
242 }//Mother Photon
243
244 //Electron check. Where did that electron come from?
245 else if(mPdg == 11){ //electron
246 SetTagBit(tag,kMCElectron);
247 if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInStack() - Checking ancestors of electrons");
248
249 //check first for B and C ancestry, then other possibilities.
250 //An electron from a photon parent could have other particles in
251 //its history and we would want to know that, right?
252
253 if(mStatus == 1) { //electron from event generator
254 if (pPdg == 23) SetTagBit(tag,kMCZDecay); //parent is Z-boson
255 else if (pPdg == 24) SetTagBit(tag,kMCWDecay); //parent is W-boson
256 else if (pPdg == 111) SetTagBit(tag,kMCPi0Decay); //Pi0 Dalitz decay
257 else if (pPdg == 221) SetTagBit(tag,kMCEtaDecay); //Eta Dalitz decay
258 else { //check the electron's ancestors for B/C contribution
259 Bool_t bAncestor = kFALSE;
260 Bool_t cAncestor = kFALSE;
261 TParticle * ancestors = stack->Particle(label); //start from electron
262 Int_t iAncestors = ancestors->GetFirstMother(); //get its mother index
263 if(iAncestors < 0) return tag; //no ancestor, shouldn't happen
264 ancestors = stack->Particle(iAncestors); //get electron mother
265 Int_t aPdg = TMath::Abs(ancestors->GetPdgCode()); //get electron mother pdg
266 if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOrigin: Scaning the decay chain for bottom/charm generated electron");
267 while(1){//check electron's mother's pdg and its mother(s)
268 //pdgs
269 if((499 < aPdg && aPdg < 600)||(4999 < aPdg && aPdg < 6000)) bAncestor = kTRUE;
270 if((399 < aPdg && aPdg < 500)||(3999 < aPdg && aPdg < 5000)) cAncestor = kTRUE;
271 if(bAncestor && cAncestor) break;
272 iAncestors = ancestors->GetFirstMother();
273 if(iAncestors < 0) break;
274 ancestors = stack->Particle(iAncestors);
275 aPdg = ancestors->GetPdgCode();
276 }//searching for ancestors
277 if(bAncestor && cAncestor) SetTagBit(tag,kMCEFromCFromB);//Decay chain has both B and C
278 else if(bAncestor && !cAncestor) SetTagBit(tag,kMCEFromB);//Decay chain has only B
279 else if(!bAncestor && cAncestor) SetTagBit(tag,kMCEFromC);//Decay chain has only C
280 }
281 //if it is not from any of the above, where is it from?
282 if(fDebug > 0) printf("Electron from other origin: %s (pdg = %d)\n",parent->GetName(),pPdg);
283 SetTagBit(tag,kMCOtherDecay);
284 }
285 else if (mStatus == 0) { //electron from GEANT
286
287 //Rewind ancestry and check for electron with status == 1
288 //if we find one, we'll assume that this object is from an
289 //electron but that it may have gone through some showering in
290 //material before the detector
291
292 //Not a double-counting problem because we are only accessing
293 //these histories for MC labels connected to a reco object.
294 //If you wanted to use this to sort through the kine stack
295 //directly, might it be a problem?
296 Bool_t eleFromEvGen = kFALSE;
297 Bool_t bAncestor = kFALSE;
298 Bool_t cAncestor = kFALSE;
299
300 TParticle * ancestors = stack->Particle(label); //start from electron
301 Int_t iAncestors = ancestors->GetFirstMother(); //get its mother index
302 if(iAncestors < 0) return tag; //no ancestor, shouldn't happen
303 ancestors = stack->Particle(iAncestors); //get electron mother
304 Int_t aStatus = ancestors->GetStatusCode(); //get electron mother's status
305 Int_t aPdg = TMath::Abs(ancestors->GetPdgCode()); //get electron mother pdg
306 if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOrigin: Scanning the decay chain for bottom/charm generated electron");
307 while(1){//check electron's mother's pdg and its mother(s) pdgs
308 if(aStatus == 1 && aPdg == 11) eleFromEvGen = kTRUE;
309 if(eleFromEvGen && aPdg == 23) SetTagBit(tag,kMCZDecay);
310 if(eleFromEvGen && aPdg == 24) SetTagBit(tag,kMCWDecay);
311 if(eleFromEvGen && ((499 < aPdg && aPdg < 600)||(4999 < aPdg && aPdg < 6000))) bAncestor = kTRUE;
312 if(eleFromEvGen && ((399 < aPdg && aPdg < 500)||(3999 < aPdg && aPdg < 5000))) cAncestor = kTRUE;
313 if(bAncestor && cAncestor) break;
314 iAncestors = ancestors->GetFirstMother();
315 if(iAncestors <= 5) break;
316 ancestors = stack->Particle(iAncestors);
317 aPdg = ancestors->GetPdgCode();
318 }//searching for ancestors
319 if(bAncestor && cAncestor) SetTagBit(tag,kMCEFromCFromB);//Decay chain has both B and C
320 else if(bAncestor && !cAncestor) SetTagBit(tag,kMCEFromB);//Decay chain has only B
321 else if(!bAncestor && cAncestor) SetTagBit(tag,kMCEFromC);//Decay chain has only C
322 if(pPdg == 22 || pPdg == 11 || pPdg == 2112 || pPdg == 211 ||
323 pPdg == 321 || pPdg == 2212 || pPdg == 130 || pPdg == 13)
324 SetTagBit(tag,kMCConversion);
325 else if (pPdg == 111) SetTagBit(tag,kMCPi0Decay);
326 else if (pPdg == 221) SetTagBit(tag,kMCEtaDecay);
327 else {
328 //if it is not from any of the above, where is it from?
329 if(pPdg > 10000) SetTagBit(tag,kMCUnknown);
330 if(fDebug > 0) printf("Electron from other origin: %s (pdg = %d)\n",parent->GetName(),pPdg);
331 SetTagBit(tag,kMCOtherDecay);
332 }
333 } //GEANT check
334 }//electron check
335 //Cluster was made by something else
336 else {
337 if(fDebug > 0) printf("\tSetting kMCUnknown for cluster from %s (pdg = %d, Parent pdg = %d)\n",mom->GetName(),mPdg,pPdg);
338 SetTagBit(tag,kMCUnknown);
339 }
340 }//Good label value
341 else{
342 if(label < 0 ) printf("AliMCAnalysisUtils::CheckOriginInStack() *** bad label or no stack ***: label %d \n", label);
343 if(label >= stack->GetNtrack()) printf("AliMCAnalysisUtils::CheckOriginInStack() *** large label ***: label %d, n tracks %d \n", label, stack->GetNtrack());
344 SetTagBit(tag,kMCUnknown);
345 }//Bad label
346
347 return tag;
348
349}
350
351
352//_________________________________________________________________________
353Int_t AliMCAnalysisUtils::CheckOriginInAOD(const Int_t label, TClonesArray *mcparticles) {
354 // Play with the MCParticles in AOD if available. Tag particles depending on their origin.
355 // Do same things as in CheckOriginInStack but different input.
356 if(!mcparticles) {
357 printf("AliMCAnalysisUtils::CheckOriginInAOD() - AODMCParticles is not available, check analysis settings in configuration file!!\n");
358
359 }
360
361 // printf(">>>>>>>>> Second Event, AODMCParticles, nMC %d\n", (GetAODMCParticles(1))->GetEntriesFast());
362 // for(Int_t i = 0; i< (GetAODMCParticles(1))->GetEntriesFast(); i++){
363 // AliAODMCParticle *particle = (AliAODMCParticle*) (GetAODMCParticles(1))->At(i);
364 // printf("** index %d, PDG %d, Flag %d, Mother %d, Daughter %d and %d, E %2.3f, pT %2.3f, phi %2.3f, eta %2.3f\n",
365 // i, particle->GetPdgCode(), particle->GetFlag(), particle->GetMother(), particle->GetDaughter(0),particle->GetDaughter(1),
366 // particle->E(),particle->Pt(), particle->Phi()*TMath::RadToDeg(), particle->Eta());
367 // if(particle->IsPhysicalPrimary()) printf("Is Physical Primary !\n");
368 // if(!particle->IsPrimary()) printf("Not Primary !\n");
369 // }
370
371
372 Int_t tag = 0;
373 Int_t nprimaries = mcparticles->GetEntriesFast();
374 if(label >= 0 && label < nprimaries){
375 //Mother
376 AliAODMCParticle * mom = (AliAODMCParticle *) mcparticles->At(label);
377 Int_t mPdg = TMath::Abs(mom->GetPdgCode());
378 Int_t iParent = mom->GetMother() ;
379 if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Mother is parton %d\n",iParent);
380
381 //GrandParent
382 AliAODMCParticle * parent = new AliAODMCParticle ;
383 Int_t pPdg = -1;
384 if(iParent >= 0){
385 parent = (AliAODMCParticle *) mcparticles->At(iParent);
386 pPdg = TMath::Abs(parent->GetPdgCode());
387 }
388 else if(fDebug > 0 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Parent with label %d\n",iParent);
389
390 //Check if mother is converted, if not, get the first non converted mother
391 if((mPdg == 22 || mPdg == 11) && (pPdg == 22 || pPdg == 11) && !mom->IsPrimary()){
392 SetTagBit(tag,kMCConversion);
393 //Check if the mother is photon or electron with status not stable
394 while ((pPdg == 22 || pPdg == 11) && !mom->IsPhysicalPrimary()) {
395 //Mother
396 mom = (AliAODMCParticle *) mcparticles->At(mom->GetMother());
397 mPdg = TMath::Abs(mom->GetPdgCode());
398 iParent = mom->GetMother() ;
399 if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Mother is parton %d\n",iParent);
400
401 //GrandParent
402 if(iParent >= 0){
403 parent = (AliAODMCParticle *) mcparticles->At(iParent);
404 pPdg = TMath::Abs(parent->GetPdgCode());
405 }
406 }//while
407 }//mother and parent are electron or photon and have status 0 and parent is photon or electron
408 else if(!mom->IsPrimary()){
409 //Still a conversion but only one electron/photon generated. Just from hadrons
410 if(pPdg == 2112 || pPdg == 211 ||
411 pPdg == 321 || pPdg == 2212 || pPdg == 130 || pPdg == 13 )
412 SetTagBit(tag,kMCConversion);
413 mom = (AliAODMCParticle *) mcparticles->At(mom->GetMother());
414 mPdg = TMath::Abs(mom->GetPdgCode());
415 //Comment for next lines, we do not check the parent of the hadron for the moment.
416 //iParent = mom->GetMother() ;
417 //if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Mother is parton %d\n",iParent);
418
419 //GrandParent
420 //if(iParent >= 0){
421 // parent = (AliAODMCParticle *) mcparticles->At(iParent);
422 // pPdg = TMath::Abs(parent->GetPdgCode());
423 //}
424 }
425
426 // conversion into electrons/photons checked
427
428 //first check for typical charged particles
429 if(mPdg == 13) SetTagBit(tag,kMCMuon);
430 else if(mPdg == 211) SetTagBit(tag,kMCPion);
431 else if(mPdg == 321) SetTagBit(tag,kMCKaon);
432 else if(mPdg == 2212) SetTagBit(tag,kMCProton);
433 //check for pi0 and eta (shouldn't happen unless their decays were
434 //turned off)
435 else if(mPdg == 111) SetTagBit(tag,kMCPi0);
436 else if(mPdg == 221) SetTagBit(tag,kMCEta);
437 //Photons
438 else if(mPdg == 22){
439 SetTagBit(tag,kMCPhoton);
440 if(mom->IsPhysicalPrimary()){ //undecayed particle
441 if(iParent < 8 && iParent > 5) {//outgoing partons
442 if(pPdg == 22) SetTagBit(tag,kMCPrompt);
443 else SetTagBit(tag,kMCFragmentation);
444 }//Outgoing partons
445 else if(iParent <= 5) {
446 SetTagBit(tag, kMCISR); //Initial state radiation
447 }
448 else if(parent->IsPrimary() && !parent->IsPhysicalPrimary()){//Decay
449 if(pPdg == 111) SetTagBit(tag,kMCPi0Decay);
450 else if (pPdg == 221) SetTagBit(tag, kMCEtaDecay);
451 else SetTagBit(tag,kMCOtherDecay);
452 }//Decay
453 else {
454 printf("AliMCAnalysisUtils::ChecOrigingInAOD() - what is it? Mother mPdg %d, is primary? %d, is physical %d \n Parent iParent %d, pPdg %d, is primary? %d, is physical? %d\n",
455 mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary(),iParent, pPdg,parent->IsPrimary(), parent->IsPhysicalPrimary());
456 SetTagBit(tag,kMCOtherDecay);//Check
457 }
458 }// Pythia generated
459 else if(!mom->IsPrimary()){ //Decays
460 if(pPdg == 111) SetTagBit(tag,kMCPi0Decay);
461 else if (pPdg == 221) SetTagBit(tag,kMCEtaDecay);
462 else SetTagBit(tag,kMCOtherDecay);
463 }//not primary : geant generated, decays
464 else {
465 //printf("UNKNOWN 1, mom pdg %d, primary %d, physical primary %d; parent %d, pdg %d, primary %d, physical primary %d \n",
466 //mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary(), iParent, pPdg, parent->IsPrimary(), parent->IsPhysicalPrimary());
467 SetTagBit(tag,kMCUnknown);
468 }
469 }//Mother Photon
470
471 //Electron check. Where did that electron come from?
472 else if(mPdg == 11){ //electron
473 SetTagBit(tag,kMCElectron);
474 if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Checking ancestors of electrons");
475
476 //check first for B and C ancestry, then other possibilities.
477 //An electron from a photon parent could have other particles in
478 //its history and we would want to know that, right?
479 Bool_t eleFromEvGen = kFALSE;
480 Bool_t bAncestor = kFALSE;
481 Bool_t cAncestor = kFALSE;
482
483 if(mom->IsPhysicalPrimary()) { //electron from event generator
484 if (pPdg == 23) SetTagBit(tag,kMCZDecay); //parent is Z-boson
485 else if (pPdg == 24) SetTagBit(tag,kMCWDecay); //parent is W-boson
486 else if (pPdg == 111) SetTagBit(tag,kMCPi0Decay); //Pi0 Dalitz decay
487 else if (pPdg == 221) SetTagBit(tag,kMCEtaDecay); //Eta Dalitz decay
488 else { //check the electron's ancestors for B/C contribution
489 bAncestor = kFALSE;
490 cAncestor = kFALSE;
491 AliAODMCParticle * ancestors = (AliAODMCParticle*)mcparticles->At(label); //start from electron
492 Int_t iAncestors = ancestors->GetMother(); //get its mother index
493 if(iAncestors < 0) return tag; //no ancestor, shouldn't happen
494 ancestors = (AliAODMCParticle *) mcparticles->At(iAncestors);
495 Int_t aPdg = TMath::Abs(ancestors->GetPdgCode());
496 if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOrigin: Scanning the decay chain for bottom/charm generated electron");
497 while(1){//check electron's mother's pdg and its mother(s) pdgs
498 if((499 < aPdg && aPdg < 600)||(4999 < aPdg && aPdg < 6000)) bAncestor = kTRUE;
499 if((399 < aPdg && aPdg < 500)||(3999 < aPdg && aPdg < 5000)) cAncestor = kTRUE;
500 if(bAncestor && cAncestor) break;
501 iAncestors = ancestors->GetMother();
502 if(iAncestors < 0) break;
503 ancestors = (AliAODMCParticle *) mcparticles->At(iAncestors);
504 aPdg = ancestors->GetPdgCode();
505 }//searching for ancestors
506 if(bAncestor && cAncestor) SetTagBit(tag,kMCEFromCFromB);//Decay chain has both B and C
507 else if(bAncestor && !cAncestor) SetTagBit(tag,kMCEFromB);//Decay chain has only B
508 else if(!bAncestor && cAncestor) SetTagBit(tag,kMCEFromC);//Decay chain has only C
509 }
510 //if it is not from any of the above, where is it from?
511 if(fDebug > 0) printf("Electron from other origin: %s (pdg = %d)\n",parent->GetName(),pPdg);
512 SetTagBit(tag,kMCOtherDecay);
513 }
514 else if (!mom->IsPrimary()) { //electron from GEANT
515
516 //Rewind ancestry and check for electron with status == 1
517 //if we find one, we'll assume that this object is from an
518 //electron but that it may have gone through some showering in
519 //material before the detector
520
521 //Not a double-counting problem because we are only accessing
522 //these histories for MC labels connected to a reco object.
523 //If you wanted to use this to sort through the kine stack
524 //directly, might it be a problem?
525 eleFromEvGen = kFALSE;
526 bAncestor = kFALSE;
527 cAncestor = kFALSE;
528
529 AliAODMCParticle * ancestors = (AliAODMCParticle *) mcparticles->At(label);
530 Int_t iAncestors = ancestors->GetMother();
531 if(iAncestors < 0) return tag; //no ancestor, shouldn't happen
532 ancestors = (AliAODMCParticle *) mcparticles->At(iAncestors); //get electron mother
533 Int_t aPdg = TMath::Abs(ancestors->GetPdgCode()); //get electron mother pdg
534 if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Scanning the decay chain for bottom/charm electrons");
535
536 while(1){//check electron's mother's pdg and its mother(s) pdgs
537 if(ancestors->IsPhysicalPrimary() && aPdg == 11) eleFromEvGen = kTRUE;
538 if(eleFromEvGen && aPdg == 23) SetTagBit(tag,kMCZDecay);
539 if(eleFromEvGen && aPdg == 24) SetTagBit(tag,kMCWDecay);
540 if(eleFromEvGen && ((499 < aPdg && aPdg < 600)||(4999 < aPdg && aPdg < 6000))) bAncestor = kTRUE;
541 if(eleFromEvGen && ((399 < aPdg && aPdg < 500)||(3999 < aPdg && aPdg < 5000))) cAncestor = kTRUE;
542 if(bAncestor && cAncestor) break;
543 iAncestors = ancestors->GetMother();
544 if(iAncestors <= 5) break; //don't go all the way back
545 ancestors = (AliAODMCParticle *) mcparticles->At(iAncestors);
546 aPdg = ancestors->GetPdgCode();
547 }//searching for ancestors
548 if(bAncestor && cAncestor) SetTagBit(tag,kMCEFromCFromB);//Decay chain has both B and C
549 else if(bAncestor && !cAncestor) SetTagBit(tag,kMCEFromB);//Decay chain has only B
550 else if(!bAncestor && cAncestor) SetTagBit(tag,kMCEFromC);//Decay chain has only C
551 if(pPdg == 22 || pPdg == 11 || pPdg == 2112 || pPdg == 211 ||
552 pPdg == 321 || pPdg == 2212 || pPdg == 130 || pPdg == 13)
553 SetTagBit(tag,kMCConversion);
554 else if (pPdg == 111) SetTagBit(tag,kMCPi0Decay);
555 else if (pPdg == 221) SetTagBit(tag,kMCEtaDecay);
556 else {
557 //if it is not from any of the above, where is it from?
558 if(pPdg > 10000) SetTagBit(tag,kMCUnknown);
559 if(fDebug > 0) printf("Electron from other origin: %s (pdg = %d)\n",parent->GetName(),pPdg);
560 SetTagBit(tag,kMCOtherDecay);
561 }
562 } //GEANT check
563 }//electron check
564 //cluster was made by something else
565 else {
566 if(fDebug > 0) printf("\tSetting kMCUnknown for cluster from %s (pdg = %d, Parent pdg = %d)\n",mom->GetName(),mPdg,pPdg);
567 SetTagBit(tag,kMCUnknown);
568 }
569 }//Good label value
570 else{
571 if(label < 0 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() *** bad label or no stack ***: label %d \n", label);
572 if(label >= mcparticles->GetEntriesFast()) printf("AliMCAnalysisUtils::CheckOriginInAOD() *** large label ***: label %d, n tracks %d \n", label, mcparticles->GetEntriesFast());
573 SetTagBit(tag,kMCUnknown);
574 }//Bad label
575
576 return tag;
577
578}
579
580
581
582//_________________________________________________________________________
583TList * AliMCAnalysisUtils::GetJets(AliCaloTrackReader * reader){
584 //Return list of jets (TParticles) and index of most likely parton that originated it.
585 AliStack * stack = reader->GetStack();
586 Int_t iEvent = reader->GetEventNumber();
587 AliGenEventHeader * geh = reader->GetGenEventHeader();
588 if(fCurrentEvent!=iEvent){
589 fCurrentEvent = iEvent;
590 fJetsList = new TList;
591 Int_t nTriggerJets = 0;
592 Float_t tmpjet[]={0,0,0,0};
593
594 //printf("Event %d %d\n",fCurrentEvent,iEvent);
595 //Get outgoing partons
596 if(stack->GetNtrack() < 8) return fJetsList;
597 TParticle * parton1 = stack->Particle(6);
598 TParticle * parton2 = stack->Particle(7);
599 if(fDebug > 2){
600 printf("AliMCAnalysisUtils::GetJets() - parton 6 : %s, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
601 parton1->GetName(),parton1->Pt(),parton1->Energy(),parton1->Phi()*TMath::RadToDeg(),parton1->Eta());
602 printf("AliMCAnalysisUtils::GetJets() - parton 7 : %s, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
603 parton2->GetName(),parton2->Pt(),parton2->Energy(),parton2->Phi()*TMath::RadToDeg(),parton2->Eta());
604 }
605// //Trace the jet from the mother parton
606// Float_t pt = 0;
607// Float_t pt1 = 0;
608// Float_t pt2 = 0;
609// Float_t e = 0;
610// Float_t e1 = 0;
611// Float_t e2 = 0;
612// TParticle * tmptmp = new TParticle;
613// for(Int_t i = 0; i< stack->GetNprimary(); i++){
614// tmptmp = stack->Particle(i);
615
616// if(tmptmp->GetStatusCode() == 1){
617// pt = tmptmp->Pt();
618// e = tmptmp->Energy();
619// Int_t imom = tmptmp->GetFirstMother();
620// Int_t imom1 = 0;
621// //printf("1st imom %d\n",imom);
622// while(imom > 5){
623// imom1=imom;
624// tmptmp = stack->Particle(imom);
625// imom = tmptmp->GetFirstMother();
626// //printf("imom %d \n",imom);
627// }
628// //printf("Last imom %d %d\n",imom1, imom);
629// if(imom1 == 6) {
630// pt1+=pt;
631// e1+=e;
632// }
633// else if (imom1 == 7){
634// pt2+=pt;
635// e2+=e; }
636// }// status 1
637
638// }// for
639
640// printf("JET 1, pt %2.2f, e %2.2f; JET 2, pt %2.2f, e %2.2f \n",pt1,e1,pt2,e2);
641
642 //Get the jet, different way for different generator
643 //PYTHIA
644 if(fMCGenerator == "PYTHIA"){
645 TParticle * jet = new TParticle;
646 AliGenPythiaEventHeader* pygeh= (AliGenPythiaEventHeader*) geh;
647 nTriggerJets = pygeh->NTriggerJets();
648 if(fDebug > 1)
649 printf("AliMCAnalysisUtils::GetJets() - PythiaEventHeader: Njets: %d\n",nTriggerJets);
650
651 Int_t iparton = -1;
652 for(Int_t i = 0; i< nTriggerJets; i++){
653 iparton=-1;
654 pygeh->TriggerJet(i, tmpjet);
655 jet = new TParticle(94, 21, -1, -1, -1, -1, tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3], 0,0,0,0);
656 //Assign an outgoing parton as mother
657 Float_t phidiff1 = TMath::Abs(jet->Phi()-parton1->Phi());
658 Float_t phidiff2 = TMath::Abs(jet->Phi()-parton2->Phi());
659 if(phidiff1 > phidiff2) jet->SetFirstMother(7);
660 else jet->SetFirstMother(6);
661 //jet->Print();
662 if(fDebug > 1)
663 printf("AliMCAnalysisUtils::GetJets() - PYTHIA Jet %d: mother %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
664 i, jet->GetFirstMother(),jet->Pt(),jet->Energy(),jet->Phi()*TMath::RadToDeg(),jet->Eta());
665 fJetsList->Add(jet);
666 }
667 }//Pythia triggered jets
668 //HERWIG
669 else if (fMCGenerator=="HERWIG"){
670 Int_t pdg = -1;
671 //Check parton 1
672 TParticle * tmp = parton1;
673 if(parton1->GetPdgCode()!=22){
674 while(pdg != 94){
675 if(tmp->GetFirstDaughter()==-1) return fJetsList;
676 tmp = stack->Particle(tmp->GetFirstDaughter());
677 pdg = tmp->GetPdgCode();
678 }//while
679
680 //Add found jet to list
681 TParticle *jet1 = new TParticle(*tmp);
682 jet1->SetFirstMother(6);
683 fJetsList->Add(jet1);
684 //printf("jet 1: first daughter %d, last daughter %d\n", tmp->GetFirstDaughter(), tmp->GetLastDaughter());
685 //tmp = stack->Particle(tmp->GetFirstDaughter());
686 //tmp->Print();
687 //jet1->Print();
688 if(fDebug > 1)
689 printf("AliMCAnalysisUtils::GetJets() - HERWIG Jet 1: mother %d, status %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
690 jet1->GetFirstMother(),jet1->GetStatusCode(),jet1->Pt(),jet1->Energy(),jet1->Phi()*TMath::RadToDeg(),jet1->Eta());
691 }//not photon
692
693 //Check parton 2
694 pdg = -1;
695 tmp = parton2;
696 Int_t i = -1;
697 if(parton2->GetPdgCode()!=22){
698 while(pdg != 94){
699 if(tmp->GetFirstDaughter()==-1) return fJetsList;
700 i = tmp->GetFirstDaughter();
701 tmp = stack->Particle(tmp->GetFirstDaughter());
702 pdg = tmp->GetPdgCode();
703 }//while
704 //Add found jet to list
705 TParticle *jet2 = new TParticle(*tmp);
706 jet2->SetFirstMother(7);
707 fJetsList->Add(jet2);
708 //jet2->Print();
709 if(fDebug > 1)
710 printf("AliMCAnalysisUtils::GetJets() - HERWIG Jet 2: mother %d, status %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
711 jet2->GetFirstMother(),jet2->GetStatusCode(),jet2->Pt(),jet2->Energy(),jet2->Phi()*TMath::RadToDeg(),jet2->Eta());
712 //Int_t first = tmp->GetFirstDaughter();
713 //Int_t last = tmp->GetLastDaughter();
714 //printf("jet 2: first daughter %d, last daughter %d, pdg %d\n",first, last, tmp->GetPdgCode());
715 // for(Int_t d = first ; d < last+1; d++){
716// tmp = stack->Particle(d);
717// if(i == tmp->GetFirstMother())
718// printf("Daughter n %d, Mother %d, name %s, status %d, pT %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
719// d,tmp->GetFirstMother(), tmp->GetName(), tmp->GetStatusCode(),tmp->Pt(),tmp->Energy(),tmp->Phi()*TMath::RadToDeg(),tmp->Eta());
720// }
721 //tmp->Print();
722 }//not photon
723 }//Herwig generated jets
724 }
725
726 return fJetsList;
727}
728
729
730//________________________________________________________________
731void AliMCAnalysisUtils::Print(const Option_t * opt) const
732{
733 //Print some relevant parameters set for the analysis
734
735 if(! opt)
736 return;
737
738 printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
739
740 printf("Debug level = %d\n",fDebug);
741 printf("MC Generator = %s\n",fMCGenerator.Data());
742 printf(" \n");
743
744}
745
746