]> 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
corrected (again) the setting of WDecay parent for electrons
[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"
8dacfd76 32#include "TDatabasePDG.h"
7cd4e982 33
34//---- ANALYSIS system ----
35#include "AliMCAnalysisUtils.h"
36#include "AliCaloTrackReader.h"
37#include "AliStack.h"
38#include "AliGenPythiaEventHeader.h"
39#include "AliAODMCParticle.h"
40
41 ClassImp(AliMCAnalysisUtils)
42
43 //________________________________________________
44 AliMCAnalysisUtils::AliMCAnalysisUtils() :
45 TObject(), fCurrentEvent(-1), fDebug(-1),
46 fJetsList(new TList), fMCGenerator("PYTHIA")
47{
48 //Ctor
49}
50
51//____________________________________________________________________________
52AliMCAnalysisUtils::AliMCAnalysisUtils(const AliMCAnalysisUtils & mcutils) :
53 TObject(mcutils), fCurrentEvent(mcutils.fCurrentEvent), fDebug(mcutils.fDebug),
54 fJetsList(mcutils.fJetsList), fMCGenerator(mcutils.fMCGenerator)
55{
56 // cpy ctor
57
58}
59
60//_________________________________________________________________________
61AliMCAnalysisUtils & AliMCAnalysisUtils::operator = (const AliMCAnalysisUtils & mcutils)
62{
63 // assignment operator
64
65 if(&mcutils == this) return *this;
66 fCurrentEvent = mcutils.fCurrentEvent ;
67 fDebug = mcutils.fDebug;
68 fJetsList = mcutils.fJetsList;
69 fMCGenerator = mcutils.fMCGenerator;
70
71 return *this;
72}
73
74//____________________________________________________________________________
75AliMCAnalysisUtils::~AliMCAnalysisUtils()
76{
77 // Remove all pointers.
78
79 if (fJetsList) {
80 fJetsList->Clear();
81 delete fJetsList ;
82 }
83}
84
85
86//_________________________________________________________________________
87Int_t AliMCAnalysisUtils::CheckOrigin(const Int_t label, AliCaloTrackReader* reader, const Int_t input = 0) {
88 //Play with the montecarlo particles if available
89 Int_t tag = 0;
90
91 //Select where the information is, ESD-galice stack or AOD mcparticles branch
92 if(reader->ReadStack()){
93 tag = CheckOriginInStack(label, reader->GetStack());
94 }
95 else if(reader->ReadAODMCParticles()){
96 tag = CheckOriginInAOD(label, reader->GetAODMCParticles(input));
97 }
98
99 return tag ;
100}
101
102//_________________________________________________________________________
103Int_t AliMCAnalysisUtils::CheckOriginInStack(const Int_t label, AliStack* stack) {
104 // Play with the MC stack if available. Tag particles depending on their origin.
105 // Do same things as in CheckOriginInAOD but different input.
556e55b0 106
107 //generally speaking, label is the MC label of a reconstructed
108 //entity (track, cluster, etc) for which we want to know something
109 //about its heritage, but one can also use it directly with stack
110 //particles not connected to reconstructed entities
111
7cd4e982 112 if(!stack) {
113 printf("AliMCAnalysisUtils::CheckOriginInStack() - Stack is not available, check analysis settings in configuration file, STOP!!\n");
114 abort();
115 }
7cd4e982 116
117 Int_t tag = 0;
556e55b0 118 if(label >= 0 && label < stack->GetNtrack()){
119 //MC particle of interest is the "mom" of the entity
7cd4e982 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
556e55b0 126 //GrandParent of the entity
7cd4e982 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);
8dacfd76 136
556e55b0 137 //Check if "mother" of entity is converted, if not, get the first non converted mother
7cd4e982 138 if((mPdg == 22 || mPdg == 11) && (pPdg == 22 || pPdg == 11) && mStatus == 0){
8dacfd76 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
7cd4e982 156 }//mother and parent are electron or photon and have status 0
8dacfd76 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::CheckOriginInStack() - 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
7cd4e982 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);
8dacfd76 181 //check for pi0 and eta (shouldn't happen unless their decays were turned off)
7cd4e982 182 else if(mPdg == 111) SetTagBit(tag,kMCPi0);
183 else if(mPdg == 221) SetTagBit(tag,kMCEta);
8dacfd76 184 //Photons
7cd4e982 185 else if(mPdg == 22){
186 SetTagBit(tag,kMCPhoton);
187 if(mStatus == 1){ //undecayed particle
8dacfd76 188 if(fMCGenerator == "PYTHIA"){
189 if(iParent < 8 && iParent > 5) {//outgoing partons
190 if(pPdg == 22) SetTagBit(tag,kMCPrompt);
191 else SetTagBit(tag,kMCFragmentation);
192 }//Outgoing partons
193 else if(iParent <= 5) {
194 SetTagBit(tag, kMCISR); //Initial state radiation
195 }
196 else if(pStatus == 11){//Decay
197 if(pPdg == 111) SetTagBit(tag,kMCPi0Decay);
198 else if (pPdg == 221) SetTagBit(tag, kMCEtaDecay);
199 else SetTagBit(tag,kMCOtherDecay);
200 }//Decay
201 else {
202 printf("AliMCAnalysisUtils::ChecOrigingInAOD() - what is it? Mother mPdg %d, status %d \n Parent iParent %d, pPdg %d %s, status %d\n",
203 mPdg, mStatus,iParent, pPdg, parent->GetName(),pStatus);
204 SetTagBit(tag,kMCOtherDecay);//Check
205 }
206 }//PYTHIA
207
208 else if(fMCGenerator == "HERWIG"){
209 if(pStatus < 197){//Not decay
210 while(1){
211 if(parent->GetFirstMother()<=5) break;
212 iParent = parent->GetFirstMother();
213 parent=stack->Particle(iParent);
214 pStatus= parent->GetStatusCode();
215 pPdg = parent->GetPdgCode();
216 }//Look for the parton
7cd4e982 217
8dacfd76 218 if(iParent < 8 && iParent > 5) {
219 if(pPdg == 22) SetTagBit(tag,kMCPrompt);
220 else SetTagBit(tag,kMCFragmentation);
221 }
222 else SetTagBit(tag,kMCISR);//Initial state radiation
223 }//Not decay
224 else{//Decay
225 if(pPdg == 111) SetTagBit(tag,kMCPi0Decay);
226 else if (pPdg == 221) SetTagBit(tag,kMCEtaDecay);
227 else SetTagBit(tag,kMCOtherDecay);
228 }//Decay
229 }//HERWIG
230
231 else SetTagBit(tag,kMCUnknown);
232
7cd4e982 233 }//Status 1 : created by event generator
8dacfd76 234
7cd4e982 235 else if(mStatus == 0){ // geant
8dacfd76 236 if(pPdg == 111) SetTagBit(tag,kMCPi0Decay);
237 else if (pPdg == 221) SetTagBit(tag,kMCEtaDecay);
238 else SetTagBit(tag,kMCOtherDecay);
7cd4e982 239 }//status 0 : geant generated
8dacfd76 240
7cd4e982 241 }//Mother Photon
8dacfd76 242
7cd4e982 243 //Electron check. Where did that electron come from?
244 else if(mPdg == 11){ //electron
e495fc0d 245 if(pPdg == 11){
246 Int_t iGrandma = parent->GetFirstMother();
247 if(iGrandma >= 0) {
248 TParticle* gma = (TParticle*)stack->Particle(iGrandma); //get mother
249 Int_t gPdg = TMath::Abs(gma->GetPdgCode());
250
251 if (gPdg == 23) { SetTagBit(tag,kMCZDecay); } //parent is Z-boson
252 else if (gPdg == 24) { SetTagBit(tag,kMCWDecay); } //parent is W-boson
253 }
254 }
7cd4e982 255 SetTagBit(tag,kMCElectron);
256 if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInStack() - Checking ancestors of electrons");
e495fc0d 257 if (pPdg == 111) { SetTagBit(tag,kMCPi0Decay); } //Pi0 Dalitz decay
d409d6b9 258 else if (pPdg == 221) { SetTagBit(tag,kMCEtaDecay); } //Eta Dalitz decay
259 else if((499 < pPdg && pPdg < 600)||(4999 < pPdg && pPdg < 6000)) { SetTagBit(tag,kMCEFromB); } //b-->e decay
260 else if((399 < pPdg && pPdg < 500)||(3999 < pPdg && pPdg < 5000)) { //check charm decay
261 Int_t iGrandma = parent->GetFirstMother();
262 if(iGrandma >= 0) {
263 TParticle* gma = (TParticle*)stack->Particle(iGrandma); //get mother of charm
264 Int_t gPdg = TMath::Abs(gma->GetPdgCode());
265 if((499 < gPdg && gPdg < 600)||(4999 < gPdg && gPdg < 6000)) SetTagBit(tag,kMCEFromCFromB); //b-->c-->e
266 else SetTagBit(tag,kMCEFromC); //c-->e
267 } else SetTagBit(tag,kMCEFromC); //c-->e
268 } else {
8dacfd76 269 //if it is not from any of the above, where is it from?
d409d6b9 270 if(pPdg > 10000) SetTagBit(tag,kMCUnknown);
271 else SetTagBit(tag,kMCOtherDecay);
272 if(fDebug > 0) printf("STACK Status %d Electron from other origin: %s (pPdg = %d) %s (mpdg = %d)\n",mStatus,parent->GetName(),pPdg,mom->GetName(),mPdg);
273 }
7cd4e982 274 }//electron check
275 //Cluster was made by something else
276 else {
8dacfd76 277 if(fDebug > 0) printf("\tSetting kMCUnknown for cluster from %s (pdg = %d, Parent pdg = %d)\n",mom->GetName(),mPdg,pPdg);
7cd4e982 278 SetTagBit(tag,kMCUnknown);
279 }
280 }//Good label value
281 else{
282 if(label < 0 ) printf("AliMCAnalysisUtils::CheckOriginInStack() *** bad label or no stack ***: label %d \n", label);
283 if(label >= stack->GetNtrack()) printf("AliMCAnalysisUtils::CheckOriginInStack() *** large label ***: label %d, n tracks %d \n", label, stack->GetNtrack());
284 SetTagBit(tag,kMCUnknown);
285 }//Bad label
286
287 return tag;
288
289}
290
291
292//_________________________________________________________________________
293Int_t AliMCAnalysisUtils::CheckOriginInAOD(const Int_t label, TClonesArray *mcparticles) {
294 // Play with the MCParticles in AOD if available. Tag particles depending on their origin.
295 // Do same things as in CheckOriginInStack but different input.
296 if(!mcparticles) {
297 printf("AliMCAnalysisUtils::CheckOriginInAOD() - AODMCParticles is not available, check analysis settings in configuration file!!\n");
298
299 }
7cd4e982 300
301 Int_t tag = 0;
302 Int_t nprimaries = mcparticles->GetEntriesFast();
303 if(label >= 0 && label < nprimaries){
304 //Mother
305 AliAODMCParticle * mom = (AliAODMCParticle *) mcparticles->At(label);
306 Int_t mPdg = TMath::Abs(mom->GetPdgCode());
307 Int_t iParent = mom->GetMother() ;
308 if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Mother is parton %d\n",iParent);
309
310 //GrandParent
311 AliAODMCParticle * parent = new AliAODMCParticle ;
312 Int_t pPdg = -1;
313 if(iParent >= 0){
314 parent = (AliAODMCParticle *) mcparticles->At(iParent);
315 pPdg = TMath::Abs(parent->GetPdgCode());
316 }
317 else if(fDebug > 0 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Parent with label %d\n",iParent);
318
8dacfd76 319 //Check if mother is converted, if not, get the first non converted mother
320 if((mPdg == 22 || mPdg == 11) && (pPdg == 22 || pPdg == 11) && !mom->IsPrimary()){
321 SetTagBit(tag,kMCConversion);
322 //Check if the mother is photon or electron with status not stable
323 while ((pPdg == 22 || pPdg == 11) && !mom->IsPhysicalPrimary()) {
324 //Mother
325 mom = (AliAODMCParticle *) mcparticles->At(mom->GetMother());
326 mPdg = TMath::Abs(mom->GetPdgCode());
327 iParent = mom->GetMother() ;
328 if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Mother is parton %d\n",iParent);
329
330 //GrandParent
331 if(iParent >= 0){
332 parent = (AliAODMCParticle *) mcparticles->At(iParent);
333 pPdg = TMath::Abs(parent->GetPdgCode());
334 }
335 }//while
336 }//mother and parent are electron or photon and have status 0 and parent is photon or electron
337 else if(!mom->IsPrimary()){
338 //Still a conversion but only one electron/photon generated. Just from hadrons
339 if(pPdg == 2112 || pPdg == 211 ||
340 pPdg == 321 || pPdg == 2212 || pPdg == 130 || pPdg == 13 )
341 SetTagBit(tag,kMCConversion);
342 mom = (AliAODMCParticle *) mcparticles->At(mom->GetMother());
343 mPdg = TMath::Abs(mom->GetPdgCode());
344 //Comment for next lines, we do not check the parent of the hadron for the moment.
345 //iParent = mom->GetMother() ;
346 //if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Mother is parton %d\n",iParent);
347
348 //GrandParent
349 //if(iParent >= 0){
350 // parent = (AliAODMCParticle *) mcparticles->At(iParent);
351 // pPdg = TMath::Abs(parent->GetPdgCode());
352 //}
353 }
354
355 // conversion into electrons/photons checked
356
357 //first check for typical charged particles
7cd4e982 358 if(mPdg == 13) SetTagBit(tag,kMCMuon);
359 else if(mPdg == 211) SetTagBit(tag,kMCPion);
360 else if(mPdg == 321) SetTagBit(tag,kMCKaon);
361 else if(mPdg == 2212) SetTagBit(tag,kMCProton);
8dacfd76 362 //check for pi0 and eta (shouldn't happen unless their decays were turned off)
7cd4e982 363 else if(mPdg == 111) SetTagBit(tag,kMCPi0);
364 else if(mPdg == 221) SetTagBit(tag,kMCEta);
8dacfd76 365 //Photons
7cd4e982 366 else if(mPdg == 22){
367 SetTagBit(tag,kMCPhoton);
368 if(mom->IsPhysicalPrimary()){ //undecayed particle
8dacfd76 369 if(iParent < 8 && iParent > 5) {//outgoing partons
370 if(pPdg == 22) SetTagBit(tag,kMCPrompt);
371 else SetTagBit(tag,kMCFragmentation);
372 }//Outgoing partons
373 else if(iParent <= 5) {
374 SetTagBit(tag, kMCISR); //Initial state radiation
375 }
376 else if(parent->IsPrimary() && !parent->IsPhysicalPrimary()){//Decay
377 if(pPdg == 111) SetTagBit(tag,kMCPi0Decay);
378 else if (pPdg == 221) SetTagBit(tag, kMCEtaDecay);
379 else SetTagBit(tag,kMCOtherDecay);
380 }//Decay
381 else {
382 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",
383 mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary(),iParent, pPdg,parent->IsPrimary(), parent->IsPhysicalPrimary());
384 SetTagBit(tag,kMCOtherDecay);//Check
385 }
386 }// Pythia generated
387 else if(!mom->IsPrimary()){ //Decays
388 if(pPdg == 111) SetTagBit(tag,kMCPi0Decay);
389 else if (pPdg == 221) SetTagBit(tag,kMCEtaDecay);
390 else SetTagBit(tag,kMCOtherDecay);
7cd4e982 391 }//not primary : geant generated, decays
392 else {
8dacfd76 393 //printf("UNKNOWN 1, mom pdg %d, primary %d, physical primary %d; parent %d, pdg %d, primary %d, physical primary %d \n",
394 //mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary(), iParent, pPdg, parent->IsPrimary(), parent->IsPhysicalPrimary());
395 SetTagBit(tag,kMCUnknown);
7cd4e982 396 }
397 }//Mother Photon
398
399 //Electron check. Where did that electron come from?
400 else if(mPdg == 11){ //electron
e495fc0d 401 if(pPdg == 11){
402 Int_t iGrandma = parent->GetMother();
403 if(iGrandma >= 0) {
404 AliAODMCParticle* gma = (AliAODMCParticle*)mcparticles->At(iGrandma);
405 Int_t gPdg = TMath::Abs(gma->GetPdgCode());
406
407 if (gPdg == 23) { SetTagBit(tag,kMCZDecay); } //parent is Z-boson
408 else if (gPdg == 24) { SetTagBit(tag,kMCWDecay); } //parent is W-boson
409 }
410 }
7cd4e982 411 SetTagBit(tag,kMCElectron);
412 if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Checking ancestors of electrons");
e495fc0d 413 if (pPdg == 111) { SetTagBit(tag,kMCPi0Decay); } //Pi0 Dalitz decay
d409d6b9 414 else if (pPdg == 221) { SetTagBit(tag,kMCEtaDecay); } //Eta Dalitz decay
415 else if((499 < pPdg && pPdg < 600)||(4999 < pPdg && pPdg < 6000)) { SetTagBit(tag,kMCEFromB);} //b-hadron decay
416 else if((399 < pPdg && pPdg < 500)||(3999 < pPdg && pPdg < 5000)) { //c-hadron decay check
417 Int_t iGrandma = parent->GetMother();
418 if(iGrandma >= 0) {
419 AliAODMCParticle* gma = (AliAODMCParticle*)mcparticles->At(iGrandma); //charm's mother
420 Int_t gPdg = TMath::Abs(gma->GetPdgCode());
421 if((499 < gPdg && gPdg < 600)||(4999 < gPdg && gPdg < 6000)) SetTagBit(tag,kMCEFromCFromB); //b-->c-->e decay
422 else SetTagBit(tag,kMCEFromC); //c-hadron decay
423 } else SetTagBit(tag,kMCEFromC); //c-hadron decay
424 } else { //prompt or other decay
425 TParticlePDG* foo = TDatabasePDG::Instance()->GetParticle(pPdg);
426 TParticlePDG* foo1 = TDatabasePDG::Instance()->GetParticle(mPdg);
427 if(fDebug > 0) printf("AOD Electron from other origin: %s (pPdg = %d) %s (mPdg = %d)\n",foo->GetName(),pPdg,foo1->GetName(),mPdg);
428 if(pPdg > 10000) SetTagBit(tag,kMCUnknown);
429 else SetTagBit(tag,kMCOtherDecay);
430 }
8dacfd76 431 }//electron check
7cd4e982 432 //cluster was made by something else
433 else {
8dacfd76 434 if(fDebug > 0) printf("\tSetting kMCUnknown for cluster from %s (pdg = %d, Parent pdg = %d)\n",mom->GetName(),mPdg,pPdg);
7cd4e982 435 SetTagBit(tag,kMCUnknown);
436 }
437 }//Good label value
438 else{
439 if(label < 0 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() *** bad label or no stack ***: label %d \n", label);
440 if(label >= mcparticles->GetEntriesFast()) printf("AliMCAnalysisUtils::CheckOriginInAOD() *** large label ***: label %d, n tracks %d \n", label, mcparticles->GetEntriesFast());
441 SetTagBit(tag,kMCUnknown);
442 }//Bad label
443
444 return tag;
445
446}
447
448
449
450//_________________________________________________________________________
8dacfd76 451TList * AliMCAnalysisUtils::GetJets(AliCaloTrackReader * const reader){
7cd4e982 452 //Return list of jets (TParticles) and index of most likely parton that originated it.
453 AliStack * stack = reader->GetStack();
454 Int_t iEvent = reader->GetEventNumber();
455 AliGenEventHeader * geh = reader->GetGenEventHeader();
456 if(fCurrentEvent!=iEvent){
457 fCurrentEvent = iEvent;
458 fJetsList = new TList;
459 Int_t nTriggerJets = 0;
460 Float_t tmpjet[]={0,0,0,0};
461
462 //printf("Event %d %d\n",fCurrentEvent,iEvent);
463 //Get outgoing partons
464 if(stack->GetNtrack() < 8) return fJetsList;
465 TParticle * parton1 = stack->Particle(6);
466 TParticle * parton2 = stack->Particle(7);
467 if(fDebug > 2){
468 printf("AliMCAnalysisUtils::GetJets() - parton 6 : %s, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
469 parton1->GetName(),parton1->Pt(),parton1->Energy(),parton1->Phi()*TMath::RadToDeg(),parton1->Eta());
470 printf("AliMCAnalysisUtils::GetJets() - parton 7 : %s, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
471 parton2->GetName(),parton2->Pt(),parton2->Energy(),parton2->Phi()*TMath::RadToDeg(),parton2->Eta());
472 }
473// //Trace the jet from the mother parton
474// Float_t pt = 0;
475// Float_t pt1 = 0;
476// Float_t pt2 = 0;
477// Float_t e = 0;
478// Float_t e1 = 0;
479// Float_t e2 = 0;
480// TParticle * tmptmp = new TParticle;
481// for(Int_t i = 0; i< stack->GetNprimary(); i++){
482// tmptmp = stack->Particle(i);
483
484// if(tmptmp->GetStatusCode() == 1){
485// pt = tmptmp->Pt();
486// e = tmptmp->Energy();
487// Int_t imom = tmptmp->GetFirstMother();
488// Int_t imom1 = 0;
489// //printf("1st imom %d\n",imom);
490// while(imom > 5){
491// imom1=imom;
492// tmptmp = stack->Particle(imom);
493// imom = tmptmp->GetFirstMother();
494// //printf("imom %d \n",imom);
495// }
496// //printf("Last imom %d %d\n",imom1, imom);
497// if(imom1 == 6) {
498// pt1+=pt;
499// e1+=e;
500// }
501// else if (imom1 == 7){
502// pt2+=pt;
503// e2+=e; }
504// }// status 1
505
506// }// for
507
508// printf("JET 1, pt %2.2f, e %2.2f; JET 2, pt %2.2f, e %2.2f \n",pt1,e1,pt2,e2);
509
510 //Get the jet, different way for different generator
511 //PYTHIA
512 if(fMCGenerator == "PYTHIA"){
513 TParticle * jet = new TParticle;
514 AliGenPythiaEventHeader* pygeh= (AliGenPythiaEventHeader*) geh;
515 nTriggerJets = pygeh->NTriggerJets();
516 if(fDebug > 1)
517 printf("AliMCAnalysisUtils::GetJets() - PythiaEventHeader: Njets: %d\n",nTriggerJets);
518
519 Int_t iparton = -1;
520 for(Int_t i = 0; i< nTriggerJets; i++){
521 iparton=-1;
522 pygeh->TriggerJet(i, tmpjet);
523 jet = new TParticle(94, 21, -1, -1, -1, -1, tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3], 0,0,0,0);
524 //Assign an outgoing parton as mother
525 Float_t phidiff1 = TMath::Abs(jet->Phi()-parton1->Phi());
526 Float_t phidiff2 = TMath::Abs(jet->Phi()-parton2->Phi());
527 if(phidiff1 > phidiff2) jet->SetFirstMother(7);
528 else jet->SetFirstMother(6);
529 //jet->Print();
530 if(fDebug > 1)
531 printf("AliMCAnalysisUtils::GetJets() - PYTHIA Jet %d: mother %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
532 i, jet->GetFirstMother(),jet->Pt(),jet->Energy(),jet->Phi()*TMath::RadToDeg(),jet->Eta());
533 fJetsList->Add(jet);
534 }
535 }//Pythia triggered jets
536 //HERWIG
537 else if (fMCGenerator=="HERWIG"){
538 Int_t pdg = -1;
539 //Check parton 1
540 TParticle * tmp = parton1;
541 if(parton1->GetPdgCode()!=22){
542 while(pdg != 94){
543 if(tmp->GetFirstDaughter()==-1) return fJetsList;
544 tmp = stack->Particle(tmp->GetFirstDaughter());
545 pdg = tmp->GetPdgCode();
546 }//while
547
548 //Add found jet to list
549 TParticle *jet1 = new TParticle(*tmp);
550 jet1->SetFirstMother(6);
551 fJetsList->Add(jet1);
552 //printf("jet 1: first daughter %d, last daughter %d\n", tmp->GetFirstDaughter(), tmp->GetLastDaughter());
553 //tmp = stack->Particle(tmp->GetFirstDaughter());
554 //tmp->Print();
555 //jet1->Print();
556 if(fDebug > 1)
557 printf("AliMCAnalysisUtils::GetJets() - HERWIG Jet 1: mother %d, status %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
558 jet1->GetFirstMother(),jet1->GetStatusCode(),jet1->Pt(),jet1->Energy(),jet1->Phi()*TMath::RadToDeg(),jet1->Eta());
559 }//not photon
560
561 //Check parton 2
562 pdg = -1;
563 tmp = parton2;
564 Int_t i = -1;
565 if(parton2->GetPdgCode()!=22){
566 while(pdg != 94){
567 if(tmp->GetFirstDaughter()==-1) return fJetsList;
568 i = tmp->GetFirstDaughter();
569 tmp = stack->Particle(tmp->GetFirstDaughter());
570 pdg = tmp->GetPdgCode();
571 }//while
572 //Add found jet to list
573 TParticle *jet2 = new TParticle(*tmp);
574 jet2->SetFirstMother(7);
575 fJetsList->Add(jet2);
576 //jet2->Print();
577 if(fDebug > 1)
578 printf("AliMCAnalysisUtils::GetJets() - HERWIG Jet 2: mother %d, status %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
579 jet2->GetFirstMother(),jet2->GetStatusCode(),jet2->Pt(),jet2->Energy(),jet2->Phi()*TMath::RadToDeg(),jet2->Eta());
580 //Int_t first = tmp->GetFirstDaughter();
581 //Int_t last = tmp->GetLastDaughter();
582 //printf("jet 2: first daughter %d, last daughter %d, pdg %d\n",first, last, tmp->GetPdgCode());
583 // for(Int_t d = first ; d < last+1; d++){
584// tmp = stack->Particle(d);
585// if(i == tmp->GetFirstMother())
586// printf("Daughter n %d, Mother %d, name %s, status %d, pT %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
587// d,tmp->GetFirstMother(), tmp->GetName(), tmp->GetStatusCode(),tmp->Pt(),tmp->Energy(),tmp->Phi()*TMath::RadToDeg(),tmp->Eta());
588// }
589 //tmp->Print();
590 }//not photon
591 }//Herwig generated jets
592 }
593
594 return fJetsList;
595}
596
597
598//________________________________________________________________
599void AliMCAnalysisUtils::Print(const Option_t * opt) const
600{
601 //Print some relevant parameters set for the analysis
602
603 if(! opt)
604 return;
605
606 printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
607
608 printf("Debug level = %d\n",fDebug);
609 printf("MC Generator = %s\n",fMCGenerator.Data());
610 printf(" \n");
611
612}
613
614