]> 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
Fixing bug caused by UUID not being generated properly.
[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),
d151d2ee 54 fJetsList(new TList), fMCGenerator(mcutils.fMCGenerator)
7cd4e982 55{
56 // cpy ctor
57
58}
59
60//_________________________________________________________________________
d151d2ee 61//AliMCAnalysisUtils & 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//}
7cd4e982 73
74//____________________________________________________________________________
75AliMCAnalysisUtils::~AliMCAnalysisUtils()
76{
77 // Remove all pointers.
78
79 if (fJetsList) {
80 fJetsList->Clear();
81 delete fJetsList ;
82 }
83}
84
85
902aa95c 86//_________________________________________________________________________
87Int_t AliMCAnalysisUtils::CheckOrigin(const Int_t * label, const Int_t nlabels, 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, nlabels, reader->GetStack());
94 }
95 else if(reader->ReadAODMCParticles()){
96 tag = CheckOriginInAOD(label, nlabels, reader->GetAODMCParticles(input));
97 }
98
99 return tag ;
100}
101
7cd4e982 102//_________________________________________________________________________
103Int_t AliMCAnalysisUtils::CheckOrigin(const Int_t label, AliCaloTrackReader* reader, const Int_t input = 0) {
104 //Play with the montecarlo particles if available
105 Int_t tag = 0;
902aa95c 106 Int_t labels[]={label};
7cd4e982 107
108 //Select where the information is, ESD-galice stack or AOD mcparticles branch
109 if(reader->ReadStack()){
902aa95c 110 tag = CheckOriginInStack(labels, 1,reader->GetStack());
7cd4e982 111 }
112 else if(reader->ReadAODMCParticles()){
902aa95c 113 tag = CheckOriginInAOD(labels, 1,reader->GetAODMCParticles(input));
7cd4e982 114 }
115
116 return tag ;
117}
118
119//_________________________________________________________________________
902aa95c 120Int_t AliMCAnalysisUtils::CheckOriginInStack(const Int_t *labels, const Int_t nlabels, AliStack* stack) {
7cd4e982 121 // Play with the MC stack if available. Tag particles depending on their origin.
122 // Do same things as in CheckOriginInAOD but different input.
556e55b0 123
124 //generally speaking, label is the MC label of a reconstructed
125 //entity (track, cluster, etc) for which we want to know something
126 //about its heritage, but one can also use it directly with stack
127 //particles not connected to reconstructed entities
128
7cd4e982 129 if(!stack) {
1cd71065 130 if (fDebug >=0)
131 printf("AliMCAnalysisUtils::CheckOriginInStack() - Stack is not available, check analysis settings in configuration file, STOP!!\n");
132 return -1;
7cd4e982 133 }
7cd4e982 134
135 Int_t tag = 0;
902aa95c 136 Int_t label=labels[0];//Most significant particle contributing to the cluster
137
556e55b0 138 if(label >= 0 && label < stack->GetNtrack()){
139 //MC particle of interest is the "mom" of the entity
7cd4e982 140 TParticle * mom = stack->Particle(label);
902aa95c 141 Int_t iMom = label;
7cd4e982 142 Int_t mPdg = TMath::Abs(mom->GetPdgCode());
143 Int_t mStatus = mom->GetStatusCode() ;
144 Int_t iParent = mom->GetFirstMother() ;
145 if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - Mother is parton %d\n",iParent);
146
556e55b0 147 //GrandParent of the entity
7cd4e982 148 TParticle * parent = new TParticle ;
149 Int_t pPdg = -1;
150 Int_t pStatus =-1;
151 if(iParent > 0){
152 parent = stack->Particle(iParent);
153 pPdg = TMath::Abs(parent->GetPdgCode());
154 pStatus = parent->GetStatusCode();
155 }
156 else if(fDebug > 0 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - Parent with label %d\n",iParent);
8dacfd76 157
902aa95c 158 if(fDebug > 2 ) {
159 printf("AliMCAnalysisUtils::CheckOriginInStack() - Cluster most contributing mother and its parent: \n");
160 printf("\t Mother label %d, pdg %d, status %d\n",iMom, mPdg, mStatus);
161 printf("\t Parent label %d, pdg %d, status %d\n",iParent, pPdg, pStatus);
162 }
163
556e55b0 164 //Check if "mother" of entity is converted, if not, get the first non converted mother
7cd4e982 165 if((mPdg == 22 || mPdg == 11) && (pPdg == 22 || pPdg == 11) && mStatus == 0){
8dacfd76 166 SetTagBit(tag,kMCConversion);
167 //Check if the mother is photon or electron with status not stable
168 while ((pPdg == 22 || pPdg == 11) && mStatus != 1) {
169 //Mother
902aa95c 170 iMom = mom->GetFirstMother();
171 mom = stack->Particle(iMom);
8dacfd76 172 mPdg = TMath::Abs(mom->GetPdgCode());
173 mStatus = mom->GetStatusCode() ;
174 iParent = mom->GetFirstMother() ;
175 if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - Mother is parton %d\n",iParent);
176
177 //GrandParent
902aa95c 178 if(iParent >= 0){
8dacfd76 179 parent = stack->Particle(iParent);
180 pPdg = TMath::Abs(parent->GetPdgCode());
181 pStatus = parent->GetStatusCode();
182 }
183 }//while
902aa95c 184 if(fDebug > 2 ) {
185 printf("AliMCAnalysisUtils::CheckOriginInStack() - Converted photon/electron: \n");
186 printf("\t Mother label %d, pdg %d, status %d\n",iMom, mPdg, mStatus);
187 printf("\t Parent label %d, pdg %d, status %d\n",iParent, pPdg, pStatus);
188 }
189
7cd4e982 190 }//mother and parent are electron or photon and have status 0
902aa95c 191 else if((mPdg == 22 || mPdg == 11) && mStatus == 0){
192 //Still a conversion but only one electron/photon generated. Just from hadrons but not decays.
193 if(pPdg == 2112 || pPdg == 211 || pPdg == 321 ||
194 pPdg == 2212 || pPdg == 130 || pPdg == 13 ) {
195 SetTagBit(tag,kMCConversion);
196 iMom = mom->GetFirstMother();
197 mom = stack->Particle(iMom);
198 mPdg = TMath::Abs(mom->GetPdgCode());
199
200 if(fDebug > 2 ) {
201 printf("AliMCAnalysisUtils::CheckOriginInStack() - Converted hadron: \n");
202 printf("\t Mother label %d, pdg %d, status %d\n",iMom, mPdg, mStatus);
203 }
204 }//hadron converted
205
8dacfd76 206 //Comment for the next lines, we do not check the parent of the hadron for the moment.
207 //iParent = mom->GetFirstMother() ;
208 //if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - Mother is parton %d\n",iParent);
209
210 //GrandParent
211 //if(iParent >= 0){
212 // parent = stack->Particle(iParent);
213 // pPdg = TMath::Abs(parent->GetPdgCode());
214 //}
215 }
902aa95c 216 // conversion into electrons/photons checked
8dacfd76 217
7cd4e982 218 //first check for typical charged particles
219 if(mPdg == 13) SetTagBit(tag,kMCMuon);
220 else if(mPdg == 211) SetTagBit(tag,kMCPion);
221 else if(mPdg == 321) SetTagBit(tag,kMCKaon);
222 else if(mPdg == 2212) SetTagBit(tag,kMCProton);
8dacfd76 223 //check for pi0 and eta (shouldn't happen unless their decays were turned off)
902aa95c 224 else if(mPdg == 111) {
225 SetTagBit(tag,kMCPi0Decay);
226 if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - First mother is directly pi0, not decayed by generator \n");
227 CheckOverlapped2GammaDecay(labels,nlabels, iMom, stack, tag); //set to kMCPi0 if 2 gammas in same cluster
228 }
229 else if(mPdg == 221) {
230 SetTagBit(tag,kMCEtaDecay);
231 if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - First mother is directly eta, not decayed by generator \n");
232 CheckOverlapped2GammaDecay(labels,nlabels, iMom, stack, tag); //set to kMCEta if 2 gammas in same cluster
233 }
8dacfd76 234 //Photons
7cd4e982 235 else if(mPdg == 22){
236 SetTagBit(tag,kMCPhoton);
237 if(mStatus == 1){ //undecayed particle
8dacfd76 238 if(fMCGenerator == "PYTHIA"){
239 if(iParent < 8 && iParent > 5) {//outgoing partons
240 if(pPdg == 22) SetTagBit(tag,kMCPrompt);
241 else SetTagBit(tag,kMCFragmentation);
242 }//Outgoing partons
243 else if(iParent <= 5) {
244 SetTagBit(tag, kMCISR); //Initial state radiation
245 }
246 else if(pStatus == 11){//Decay
902aa95c 247 if(pPdg == 111) {
248 SetTagBit(tag,kMCPi0Decay);
249 if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - PYTHIA pi0 decay photon, parent pi0 with status 11 \n");
250 CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag); //set to kMCPi0 if 2 gammas in same cluster
251 }
252 else if (pPdg == 221) {
253 SetTagBit(tag, kMCEtaDecay);
254 if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - PYTHIA eta decay photon, parent pi0 with status 11 \n");
255 CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag);//set to kMCEta if 2 gammas in same cluster
256 }
257 else SetTagBit(tag,kMCOtherDecay);
8dacfd76 258 }//Decay
259 else {
74f449c1 260 if(fDebug > 1 ) printf("AliMCAnalysisUtils::CheckOrigingInStack() - what is it in PYTHIA? Wrong generator setting? Mother mPdg %d, status %d \n Parent iParent %d, pPdg %d %s, status %d\n",
8dacfd76 261 mPdg, mStatus,iParent, pPdg, parent->GetName(),pStatus);
74f449c1 262 if(pPdg == 111) {
263 SetTagBit(tag,kMCPi0Decay);
264 if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - PYTHIA pi0 decay photon, parent pi0 with status 11 \n");
265 CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag); //set to kMCPi0 if 2 gammas in same cluster
266 }
267 else if (pPdg == 221) {
268 SetTagBit(tag, kMCEtaDecay);
269 if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - PYTHIA eta decay photon, parent pi0 with status 11 \n");
270 CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag);//set to kMCEta if 2 gammas in same cluster
271 }
272 else SetTagBit(tag,kMCOtherDecay);
8dacfd76 273 }
274 }//PYTHIA
275
276 else if(fMCGenerator == "HERWIG"){
277 if(pStatus < 197){//Not decay
278 while(1){
279 if(parent->GetFirstMother()<=5) break;
280 iParent = parent->GetFirstMother();
281 parent=stack->Particle(iParent);
282 pStatus= parent->GetStatusCode();
902aa95c 283 pPdg = TMath::Abs(parent->GetPdgCode());
8dacfd76 284 }//Look for the parton
7cd4e982 285
8dacfd76 286 if(iParent < 8 && iParent > 5) {
287 if(pPdg == 22) SetTagBit(tag,kMCPrompt);
288 else SetTagBit(tag,kMCFragmentation);
289 }
290 else SetTagBit(tag,kMCISR);//Initial state radiation
291 }//Not decay
292 else{//Decay
902aa95c 293 if(pPdg == 111) {
294 SetTagBit(tag,kMCPi0Decay);
295 if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - HERWIG pi0 decay photon \n");
296 CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag); //set to kMCPi0 if 2 gammas in same cluster
297 }
298 else if (pPdg == 221) {
299 SetTagBit(tag,kMCEtaDecay);
300 if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - HERWIG eta decay photon \n");
301 CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag); //set to kMCEta if 2 gammas in same cluster
302 }
8dacfd76 303 else SetTagBit(tag,kMCOtherDecay);
304 }//Decay
305 }//HERWIG
306
307 else SetTagBit(tag,kMCUnknown);
308
7cd4e982 309 }//Status 1 : created by event generator
8dacfd76 310
7cd4e982 311 else if(mStatus == 0){ // geant
902aa95c 312 if(pPdg == 111) {
313 SetTagBit(tag,kMCPi0Decay);
314 if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - Transport MC pi0 decay photon \n");
315 CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag); //set to kMCPi0 if 2 gammas in same cluster
316 }
317 else if (pPdg == 221) {
318 SetTagBit(tag,kMCEtaDecay);
319 if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - Transport MC eta decay photon \n");
320 CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag); //set to kMCEta if 2 gammas in same cluster
321 }
8dacfd76 322 else SetTagBit(tag,kMCOtherDecay);
7cd4e982 323 }//status 0 : geant generated
8dacfd76 324
7cd4e982 325 }//Mother Photon
8dacfd76 326
7cd4e982 327 //Electron check. Where did that electron come from?
328 else if(mPdg == 11){ //electron
e495fc0d 329 if(pPdg == 11){
330 Int_t iGrandma = parent->GetFirstMother();
331 if(iGrandma >= 0) {
332 TParticle* gma = (TParticle*)stack->Particle(iGrandma); //get mother
333 Int_t gPdg = TMath::Abs(gma->GetPdgCode());
334
335 if (gPdg == 23) { SetTagBit(tag,kMCZDecay); } //parent is Z-boson
336 else if (gPdg == 24) { SetTagBit(tag,kMCWDecay); } //parent is W-boson
337 }
338 }
7cd4e982 339 SetTagBit(tag,kMCElectron);
902aa95c 340 if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInStack() - Checking ancestors of electrons\n");
e495fc0d 341 if (pPdg == 111) { SetTagBit(tag,kMCPi0Decay); } //Pi0 Dalitz decay
d409d6b9 342 else if (pPdg == 221) { SetTagBit(tag,kMCEtaDecay); } //Eta Dalitz decay
343 else if((499 < pPdg && pPdg < 600)||(4999 < pPdg && pPdg < 6000)) { SetTagBit(tag,kMCEFromB); } //b-->e decay
344 else if((399 < pPdg && pPdg < 500)||(3999 < pPdg && pPdg < 5000)) { //check charm decay
345 Int_t iGrandma = parent->GetFirstMother();
346 if(iGrandma >= 0) {
347 TParticle* gma = (TParticle*)stack->Particle(iGrandma); //get mother of charm
348 Int_t gPdg = TMath::Abs(gma->GetPdgCode());
349 if((499 < gPdg && gPdg < 600)||(4999 < gPdg && gPdg < 6000)) SetTagBit(tag,kMCEFromCFromB); //b-->c-->e
350 else SetTagBit(tag,kMCEFromC); //c-->e
351 } else SetTagBit(tag,kMCEFromC); //c-->e
352 } else {
8dacfd76 353 //if it is not from any of the above, where is it from?
d409d6b9 354 if(pPdg > 10000) SetTagBit(tag,kMCUnknown);
355 else SetTagBit(tag,kMCOtherDecay);
902aa95c 356 if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInStack() - Status %d Electron from other origin: %s (pPdg = %d) %s (mpdg = %d)\n",mStatus,parent->GetName(),pPdg,mom->GetName(),mPdg);
d409d6b9 357 }
7cd4e982 358 }//electron check
359 //Cluster was made by something else
360 else {
902aa95c 361 if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInStack() - \tSetting kMCUnknown for cluster from %s (pdg = %d, Parent pdg = %d)\n",mom->GetName(),mPdg,pPdg);
7cd4e982 362 SetTagBit(tag,kMCUnknown);
363 }
364 }//Good label value
902aa95c 365 else{// Bad label
366
367 if(label < 0 && (fDebug >= 0))
368 printf("AliMCAnalysisUtils::CheckOriginInStack() *** bad label or no stack ***: label %d \n", label);
369 if(label >= stack->GetNtrack() && (fDebug >= 0))
370 printf("AliMCAnalysisUtils::CheckOriginInStack() *** large label ***: label %d, n tracks %d \n", label, stack->GetNtrack());
7cd4e982 371 SetTagBit(tag,kMCUnknown);
372 }//Bad label
373
374 return tag;
375
376}
377
378
379//_________________________________________________________________________
902aa95c 380Int_t AliMCAnalysisUtils::CheckOriginInAOD(const Int_t *labels, const Int_t nlabels, TClonesArray *mcparticles) {
7cd4e982 381 // Play with the MCParticles in AOD if available. Tag particles depending on their origin.
382 // Do same things as in CheckOriginInStack but different input.
383 if(!mcparticles) {
1cd71065 384 if(fDebug >= 0)
385 printf("AliMCAnalysisUtils::CheckOriginInAOD() - AODMCParticles is not available, check analysis settings in configuration file!!\n");
386 return -1;
7cd4e982 387 }
7cd4e982 388
389 Int_t tag = 0;
902aa95c 390 Int_t label=labels[0];//Most significant particle contributing to the cluster
391
7cd4e982 392 Int_t nprimaries = mcparticles->GetEntriesFast();
393 if(label >= 0 && label < nprimaries){
394 //Mother
395 AliAODMCParticle * mom = (AliAODMCParticle *) mcparticles->At(label);
902aa95c 396 Int_t iMom = label;
7cd4e982 397 Int_t mPdg = TMath::Abs(mom->GetPdgCode());
398 Int_t iParent = mom->GetMother() ;
399 if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Mother is parton %d\n",iParent);
400
401 //GrandParent
402 AliAODMCParticle * parent = new AliAODMCParticle ;
403 Int_t pPdg = -1;
404 if(iParent >= 0){
405 parent = (AliAODMCParticle *) mcparticles->At(iParent);
406 pPdg = TMath::Abs(parent->GetPdgCode());
407 }
408 else if(fDebug > 0 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Parent with label %d\n",iParent);
409
902aa95c 410 if(fDebug > 2 ) {
411 printf("AliMCAnalysisUtils::CheckOriginInAOD() - Cluster most contributing mother and its parent: \n");
412 printf("\t Mother label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iMom, mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary());
413 printf("\t Parent label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iParent, pPdg, parent->IsPrimary(), parent->IsPhysicalPrimary());
414 }
415
8dacfd76 416 //Check if mother is converted, if not, get the first non converted mother
417 if((mPdg == 22 || mPdg == 11) && (pPdg == 22 || pPdg == 11) && !mom->IsPrimary()){
418 SetTagBit(tag,kMCConversion);
419 //Check if the mother is photon or electron with status not stable
420 while ((pPdg == 22 || pPdg == 11) && !mom->IsPhysicalPrimary()) {
421 //Mother
902aa95c 422 iMom = mom->GetMother();
423 mom = (AliAODMCParticle *) mcparticles->At(iMom);
8dacfd76 424 mPdg = TMath::Abs(mom->GetPdgCode());
425 iParent = mom->GetMother() ;
426 if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Mother is parton %d\n",iParent);
427
428 //GrandParent
429 if(iParent >= 0){
430 parent = (AliAODMCParticle *) mcparticles->At(iParent);
431 pPdg = TMath::Abs(parent->GetPdgCode());
432 }
902aa95c 433 // printf("\t While Mother label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iMom, mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary());
434 // printf("\t While Parent label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iParent, pPdg, parent->IsPrimary(), parent->IsPhysicalPrimary());
435
436 }//while
437
438 if(fDebug > 2 ) {
439 printf("AliMCAnalysisUtils::CheckOriginInAOD() - Converted photon/electron : \n");
440 printf("\t Mother label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iMom, mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary());
441 printf("\t Parent label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iParent, pPdg, parent->IsPrimary(), parent->IsPhysicalPrimary());
442 }
443
8dacfd76 444 }//mother and parent are electron or photon and have status 0 and parent is photon or electron
902aa95c 445 else if((mPdg == 22 || mPdg == 11) && !mom->IsPrimary()){
8dacfd76 446 //Still a conversion but only one electron/photon generated. Just from hadrons
902aa95c 447 if(pPdg == 2112 || pPdg == 211 || pPdg == 321 ||
448 pPdg == 2212 || pPdg == 130 || pPdg == 13 ) {
449 SetTagBit(tag,kMCConversion);
450 iMom = mom->GetMother();
451 mom = (AliAODMCParticle *) mcparticles->At(iMom);
452 mPdg = TMath::Abs(mom->GetPdgCode());
453
454 if(fDebug > 2 ) {
455 printf("AliMCAnalysisUtils::CheckOriginInAOD() - Converted hadron : \n");
456 printf("\t Mother label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iMom, mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary());
457 }
458 }//hadron converted
459
8dacfd76 460 //Comment for next lines, we do not check the parent of the hadron for the moment.
461 //iParent = mom->GetMother() ;
462 //if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Mother is parton %d\n",iParent);
463
464 //GrandParent
465 //if(iParent >= 0){
466 // parent = (AliAODMCParticle *) mcparticles->At(iParent);
467 // pPdg = TMath::Abs(parent->GetPdgCode());
468 //}
469 }
470
471 // conversion into electrons/photons checked
472
473 //first check for typical charged particles
7cd4e982 474 if(mPdg == 13) SetTagBit(tag,kMCMuon);
475 else if(mPdg == 211) SetTagBit(tag,kMCPion);
476 else if(mPdg == 321) SetTagBit(tag,kMCKaon);
477 else if(mPdg == 2212) SetTagBit(tag,kMCProton);
8dacfd76 478 //check for pi0 and eta (shouldn't happen unless their decays were turned off)
902aa95c 479 else if(mPdg == 111) {
480 SetTagBit(tag,kMCPi0Decay);
481 if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - First mother is directly pi0, not decayed by generator \n");
482 CheckOverlapped2GammaDecay(labels,nlabels, iMom, mcparticles, tag); //set to kMCPi0 if 2 gammas in same cluster
483 }
484 else if(mPdg == 221) {
485 SetTagBit(tag,kMCEtaDecay);
486 if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - First mother is directly eta, not decayed by generator \n");
487 CheckOverlapped2GammaDecay(labels,nlabels, iMom, mcparticles, tag); //set to kMCEta if 2 gammas in same cluster
488 }
8dacfd76 489 //Photons
7cd4e982 490 else if(mPdg == 22){
491 SetTagBit(tag,kMCPhoton);
492 if(mom->IsPhysicalPrimary()){ //undecayed particle
8dacfd76 493 if(iParent < 8 && iParent > 5) {//outgoing partons
494 if(pPdg == 22) SetTagBit(tag,kMCPrompt);
495 else SetTagBit(tag,kMCFragmentation);
496 }//Outgoing partons
497 else if(iParent <= 5) {
498 SetTagBit(tag, kMCISR); //Initial state radiation
499 }
500 else if(parent->IsPrimary() && !parent->IsPhysicalPrimary()){//Decay
902aa95c 501 if(pPdg == 111){
502 SetTagBit(tag,kMCPi0Decay);
503 if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Generator pi0 decay photon \n");
504 CheckOverlapped2GammaDecay(labels,nlabels, iParent, mcparticles, tag); //set to kMCPi0 if 2 gammas in same cluster
505 }
506 else if (pPdg == 221) {
507 SetTagBit(tag, kMCEtaDecay);
508 if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Generator eta decay photon \n");
509 CheckOverlapped2GammaDecay(labels,nlabels, iParent, mcparticles, tag); //set to kMCEta if 2 gammas in same cluster
510 }
8dacfd76 511 else SetTagBit(tag,kMCOtherDecay);
512 }//Decay
513 else {
514 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",
515 mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary(),iParent, pPdg,parent->IsPrimary(), parent->IsPhysicalPrimary());
516 SetTagBit(tag,kMCOtherDecay);//Check
517 }
518 }// Pythia generated
519 else if(!mom->IsPrimary()){ //Decays
902aa95c 520 if(pPdg == 111){
521 SetTagBit(tag,kMCPi0Decay);
522 if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Transport MC pi0 decay photon \n");
523 CheckOverlapped2GammaDecay(labels,nlabels, iParent, mcparticles, tag); //set to kMCPi0 if 2 gammas in same cluster
524 }
525 else if (pPdg == 221) {
526 SetTagBit(tag,kMCEtaDecay);
527 if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Transport MC eta decay photon \n");
528 CheckOverlapped2GammaDecay(labels,nlabels, iParent, mcparticles, tag); //set to kMCEta if 2 gammas in same cluster
529 }
8dacfd76 530 else SetTagBit(tag,kMCOtherDecay);
7cd4e982 531 }//not primary : geant generated, decays
532 else {
8dacfd76 533 //printf("UNKNOWN 1, mom pdg %d, primary %d, physical primary %d; parent %d, pdg %d, primary %d, physical primary %d \n",
534 //mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary(), iParent, pPdg, parent->IsPrimary(), parent->IsPhysicalPrimary());
535 SetTagBit(tag,kMCUnknown);
7cd4e982 536 }
537 }//Mother Photon
538
539 //Electron check. Where did that electron come from?
540 else if(mPdg == 11){ //electron
e495fc0d 541 if(pPdg == 11){
542 Int_t iGrandma = parent->GetMother();
543 if(iGrandma >= 0) {
544 AliAODMCParticle* gma = (AliAODMCParticle*)mcparticles->At(iGrandma);
545 Int_t gPdg = TMath::Abs(gma->GetPdgCode());
546
547 if (gPdg == 23) { SetTagBit(tag,kMCZDecay); } //parent is Z-boson
548 else if (gPdg == 24) { SetTagBit(tag,kMCWDecay); } //parent is W-boson
549 }
550 }
7cd4e982 551 SetTagBit(tag,kMCElectron);
552 if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Checking ancestors of electrons");
e495fc0d 553 if (pPdg == 111) { SetTagBit(tag,kMCPi0Decay); } //Pi0 Dalitz decay
d409d6b9 554 else if (pPdg == 221) { SetTagBit(tag,kMCEtaDecay); } //Eta Dalitz decay
555 else if((499 < pPdg && pPdg < 600)||(4999 < pPdg && pPdg < 6000)) { SetTagBit(tag,kMCEFromB);} //b-hadron decay
556 else if((399 < pPdg && pPdg < 500)||(3999 < pPdg && pPdg < 5000)) { //c-hadron decay check
557 Int_t iGrandma = parent->GetMother();
558 if(iGrandma >= 0) {
559 AliAODMCParticle* gma = (AliAODMCParticle*)mcparticles->At(iGrandma); //charm's mother
560 Int_t gPdg = TMath::Abs(gma->GetPdgCode());
561 if((499 < gPdg && gPdg < 600)||(4999 < gPdg && gPdg < 6000)) SetTagBit(tag,kMCEFromCFromB); //b-->c-->e decay
562 else SetTagBit(tag,kMCEFromC); //c-hadron decay
563 } else SetTagBit(tag,kMCEFromC); //c-hadron decay
564 } else { //prompt or other decay
565 TParticlePDG* foo = TDatabasePDG::Instance()->GetParticle(pPdg);
566 TParticlePDG* foo1 = TDatabasePDG::Instance()->GetParticle(mPdg);
902aa95c 567 if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Electron from other origin: %s (pPdg = %d) %s (mPdg = %d)\n",foo->GetName(), pPdg,foo1->GetName(),mPdg);
d409d6b9 568 if(pPdg > 10000) SetTagBit(tag,kMCUnknown);
569 else SetTagBit(tag,kMCOtherDecay);
570 }
8dacfd76 571 }//electron check
7cd4e982 572 //cluster was made by something else
573 else {
902aa95c 574 if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInAOD() - \tSetting kMCUnknown for cluster with pdg = %d, Parent pdg = %d\n",mPdg,pPdg);
7cd4e982 575 SetTagBit(tag,kMCUnknown);
576 }
577 }//Good label value
902aa95c 578 else{//Bad label
579
580 if(label < 0 && (fDebug >= 0) )
581 printf("AliMCAnalysisUtils::CheckOriginInAOD() *** bad label or no mcparticles ***: label %d \n", label);
582 if(label >= mcparticles->GetEntriesFast() && (fDebug >= 0) )
583 printf("AliMCAnalysisUtils::CheckOriginInAOD() *** large label ***: label %d, n tracks %d \n", label, mcparticles->GetEntriesFast());
7cd4e982 584 SetTagBit(tag,kMCUnknown);
902aa95c 585
7cd4e982 586 }//Bad label
587
588 return tag;
589
590}
591
902aa95c 592//_________________________________________________________________________
593void AliMCAnalysisUtils::CheckOverlapped2GammaDecay(const Int_t *labels, const Int_t nlabels, const Int_t mesonIndex,
594 AliStack *stack, Int_t &tag) {
595 //Check if cluster is formed from the contribution of 2 decay photons from pi0 or eta. Input in stack
596
597 if(labels[0] < 0 || labels[0] > stack->GetNtrack() || nlabels <= 1) {
598 if(fDebug > 2) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Exit : label[0] %d, n primaries %d, nlabels %d \n",
599 labels[0],stack->GetNtrack(), nlabels);
600 return;
601 }
602
603 TParticle * meson = stack->Particle(mesonIndex);
604 Int_t mesonPdg = meson->GetPdgCode();
605 if(mesonPdg!=111 && mesonPdg!=221){
606 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Wrong pi0/eta PDG : %d \n",mesonPdg);
607 return;
608 }
609
610 if(fDebug > 2) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - %s, label %d\n",meson->GetName(), mesonIndex);
611
612 //Check if meson decayed into 2 daughters or if both were kept.
613 if(meson->GetNDaughters() != 2){
614 if(fDebug > 2)
615 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Not overalapped. Number of daughters is %d, not 2 \n",meson->GetNDaughters());
616 return;
617 }
618
619 //Get the daughters
620 Int_t iPhoton0 = meson->GetDaughter(0);
621 Int_t iPhoton1 = meson->GetDaughter(1);
622 TParticle *photon0 = stack->Particle(iPhoton0);
623 TParticle *photon1 = stack->Particle(iPhoton1);
624
625 //Check if both daughters are photons
626 if(photon0->GetPdgCode() != 22 || photon1->GetPdgCode()!=22){
627 if(fDebug > 2)
628 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Not overalapped. PDG: daughter 1 = %d, of daughter 2 = %d \n",photon0->GetPdgCode(),photon1->GetPdgCode());
629 return;
630 }
631
632 if(fDebug > 2)
633 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Daughter labels : photon0 = %d, photon1 = %d \n",iPhoton0,iPhoton1);
634
635 //Check if both photons contribute to the cluster
636 Bool_t okPhoton0 = kFALSE;
637 Bool_t okPhoton1 = kFALSE;
638
639 if(fDebug > 3) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Labels loop:\n");
640
641 for(Int_t i = 0; i < nlabels; i++){
642 if(fDebug > 3) printf("\t at begin:label %d/%d: %d, ok? photon1 %d, photon2 %d\n", i+1, nlabels, labels[i], okPhoton0, okPhoton1);
643
644 //If we already found both, break the loop
645 if(okPhoton0 && okPhoton1) break;
646
647 Int_t index = labels[i];
648 if (iPhoton0 == index) {
649 okPhoton0 = kTRUE;
650 continue;
651 }
652 else if (iPhoton1 == index) {
653 okPhoton1 = kTRUE;
654 continue;
655 }
656
657 //Trace back the mother in case it was a conversion
658 TParticle * daught = stack->Particle(index);
659 Int_t tmpindex = daught->GetFirstMother();
660 if(fDebug > 3) printf("\t Conversion? : mother %d\n",tmpindex);
661 while(tmpindex>=0){
662 //MC particle of interest is the mother
663 if(fDebug > 3) printf("\t \t parent index %d\n",tmpindex);
664 daught = stack->Particle(tmpindex);
665 if (iPhoton0 == tmpindex) {
666 okPhoton0 = kTRUE;
667 break;
668 }
669 else if (iPhoton1 == tmpindex) {
670 okPhoton1 = kTRUE;
671 break;
672 }
673 tmpindex = daught->GetFirstMother();
674 }//While to check if pi0/eta daughter was one of these contributors to the cluster
675
676 if(i == 0 && (!okPhoton0 && !okPhoton1) && fDebug>=0)
677 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Something happens, first label should be from a photon decay!\n");
678
679 }//loop on list of labels
680
681 //If both photons contribute tag as the corresponding meson.
682 if(okPhoton0 && okPhoton1){
683 if(fDebug > 2)
684 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - %s OVERLAPPED DECAY \n", meson->GetName());
685
686 if(mesonPdg == 111) SetTagBit(tag,kMCPi0);
687 else SetTagBit(tag,kMCEta);
688 }
689
690}
691
692//_________________________________________________________________________
693void AliMCAnalysisUtils::CheckOverlapped2GammaDecay(const Int_t *labels, const Int_t nlabels, const Int_t mesonIndex,
694 TClonesArray *mcparticles, Int_t &tag) {
695 //Check if cluster is formed from the contribution of 2 decay photons from pi0 or eta. Input in AODMCParticles
696
697 if(labels[0] < 0 || labels[0] > mcparticles->GetEntriesFast() || nlabels <= 1) {
698 if(fDebug > 2) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Exit : label[0] %d, n primaries %d, nlabels %d \n",
699 labels[0],mcparticles->GetEntriesFast(), nlabels);
700 return;
701 }
702
703
704 AliAODMCParticle * meson = (AliAODMCParticle *) mcparticles->At(mesonIndex);
705 Int_t mesonPdg = meson->GetPdgCode();
706 if(mesonPdg != 111 && mesonPdg != 221) {
707 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Wrong pi0/eta PDG : %d \n",mesonPdg);
708 return;
709 }
710
711 if(fDebug > 2) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - pdg %d, label %d, ndaughters %d\n", mesonPdg, mesonIndex, meson->GetNDaughters());
712
713
714 //Get the daughters
715 if(meson->GetNDaughters() != 2){
716 if(fDebug > 2)
717 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Not overalapped. Number of daughters is %d, not 2 \n",meson->GetNDaughters());
718 return;
719 }
720 Int_t iPhoton0 = meson->GetDaughter(0);
721 Int_t iPhoton1 = meson->GetDaughter(1);
722 //if((iPhoton0 == -1) || (iPhoton1 == -1)){
723 // if(fDebug > 2)
724 // printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Exit : Not overalapped. At least a daughter do not exists : d1 %d, d2 %d \n", iPhoton0, iPhoton1);
725 // return;
726 //}
727 AliAODMCParticle *photon0 = (AliAODMCParticle *) mcparticles->At(iPhoton0);
728 AliAODMCParticle *photon1 = (AliAODMCParticle *) mcparticles->At(iPhoton1);
729
730 //Check if both daughters are photons
731 if(photon0->GetPdgCode() != 22 && photon1->GetPdgCode()!=22){
732 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Not overalapped. PDG: daughter 1 = %d, of daughter 2 = %d \n",photon0->GetPdgCode(),photon1->GetPdgCode());
733 return;
734 }
735
736 if(fDebug > 2)
737 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Daughter labels : photon0 = %d, photon1 = %d \n",iPhoton0,iPhoton1);
738
739 //Check if both photons contribute to the cluster
740 Bool_t okPhoton0 = kFALSE;
741 Bool_t okPhoton1 = kFALSE;
742
743 if(fDebug > 3) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Labels loop:\n");
744
745 for(Int_t i = 0; i < nlabels; i++){
746 if(fDebug > 3) printf("\t label %d/%d: %d, ok? %d, %d\n", i, nlabels, labels[i], okPhoton0, okPhoton1);
747
748 //If we already found both, break the loop
749 if(okPhoton0 && okPhoton1) break;
750
751 Int_t index = labels[i];
752 if (iPhoton0 == index) {
753 okPhoton0 = kTRUE;
754 continue;
755 }
756 else if (iPhoton1 == index) {
757 okPhoton1 = kTRUE;
758 continue;
759 }
760
761 //Trace back the mother in case it was a conversion
762 AliAODMCParticle * daught = (AliAODMCParticle*) mcparticles->At(index);
763 Int_t tmpindex = daught->GetMother();
764 if(fDebug > 3) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Conversion? : mother %d\n",tmpindex);
7cd4e982 765
902aa95c 766 while(tmpindex>=0){
767
768 //MC particle of interest is the mother
769 if(fDebug > 3) printf("\t parent index %d\n",tmpindex);
770 daught = (AliAODMCParticle*) mcparticles->At(tmpindex);
771 //printf("tmpindex %d\n",tmpindex);
772 if (iPhoton0 == tmpindex) {
773 okPhoton0 = kTRUE;
774 break;
775 }
776 else if (iPhoton1 == tmpindex) {
777 okPhoton1 = kTRUE;
778 break;
779 }
780 tmpindex = daught->GetMother();
781 }//While to check if pi0/eta daughter was one of these contributors to the cluster
782
783 if(i == 0 && (!okPhoton0 && !okPhoton1) && fDebug>=0 )
784 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Something happens, first label should be from a photon decay!\n");
785
786 }//loop on list of labels
787
788 //If both photons contribute tag as the corresponding meson.
789 if(okPhoton0 && okPhoton1){
790 if(fDebug > 2)
791 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - %s OVERLAPPED DECAY \n",(TDatabasePDG::Instance()->GetParticle(mesonPdg))->GetName());
792
793 if(mesonPdg == 111) SetTagBit(tag,kMCPi0);
794 else SetTagBit(tag,kMCEta);
795 }
796
797}
7cd4e982 798
799//_________________________________________________________________________
8dacfd76 800TList * AliMCAnalysisUtils::GetJets(AliCaloTrackReader * const reader){
7cd4e982 801 //Return list of jets (TParticles) and index of most likely parton that originated it.
802 AliStack * stack = reader->GetStack();
803 Int_t iEvent = reader->GetEventNumber();
804 AliGenEventHeader * geh = reader->GetGenEventHeader();
805 if(fCurrentEvent!=iEvent){
806 fCurrentEvent = iEvent;
807 fJetsList = new TList;
808 Int_t nTriggerJets = 0;
809 Float_t tmpjet[]={0,0,0,0};
810
811 //printf("Event %d %d\n",fCurrentEvent,iEvent);
812 //Get outgoing partons
813 if(stack->GetNtrack() < 8) return fJetsList;
814 TParticle * parton1 = stack->Particle(6);
815 TParticle * parton2 = stack->Particle(7);
816 if(fDebug > 2){
817 printf("AliMCAnalysisUtils::GetJets() - parton 6 : %s, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
818 parton1->GetName(),parton1->Pt(),parton1->Energy(),parton1->Phi()*TMath::RadToDeg(),parton1->Eta());
819 printf("AliMCAnalysisUtils::GetJets() - parton 7 : %s, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
820 parton2->GetName(),parton2->Pt(),parton2->Energy(),parton2->Phi()*TMath::RadToDeg(),parton2->Eta());
821 }
822// //Trace the jet from the mother parton
823// Float_t pt = 0;
824// Float_t pt1 = 0;
825// Float_t pt2 = 0;
826// Float_t e = 0;
827// Float_t e1 = 0;
828// Float_t e2 = 0;
829// TParticle * tmptmp = new TParticle;
830// for(Int_t i = 0; i< stack->GetNprimary(); i++){
831// tmptmp = stack->Particle(i);
832
833// if(tmptmp->GetStatusCode() == 1){
834// pt = tmptmp->Pt();
835// e = tmptmp->Energy();
836// Int_t imom = tmptmp->GetFirstMother();
837// Int_t imom1 = 0;
838// //printf("1st imom %d\n",imom);
839// while(imom > 5){
840// imom1=imom;
841// tmptmp = stack->Particle(imom);
842// imom = tmptmp->GetFirstMother();
843// //printf("imom %d \n",imom);
844// }
845// //printf("Last imom %d %d\n",imom1, imom);
846// if(imom1 == 6) {
847// pt1+=pt;
848// e1+=e;
849// }
850// else if (imom1 == 7){
851// pt2+=pt;
852// e2+=e; }
853// }// status 1
854
855// }// for
856
857// printf("JET 1, pt %2.2f, e %2.2f; JET 2, pt %2.2f, e %2.2f \n",pt1,e1,pt2,e2);
858
859 //Get the jet, different way for different generator
860 //PYTHIA
861 if(fMCGenerator == "PYTHIA"){
862 TParticle * jet = new TParticle;
863 AliGenPythiaEventHeader* pygeh= (AliGenPythiaEventHeader*) geh;
864 nTriggerJets = pygeh->NTriggerJets();
865 if(fDebug > 1)
866 printf("AliMCAnalysisUtils::GetJets() - PythiaEventHeader: Njets: %d\n",nTriggerJets);
867
868 Int_t iparton = -1;
869 for(Int_t i = 0; i< nTriggerJets; i++){
870 iparton=-1;
871 pygeh->TriggerJet(i, tmpjet);
872 jet = new TParticle(94, 21, -1, -1, -1, -1, tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3], 0,0,0,0);
873 //Assign an outgoing parton as mother
874 Float_t phidiff1 = TMath::Abs(jet->Phi()-parton1->Phi());
875 Float_t phidiff2 = TMath::Abs(jet->Phi()-parton2->Phi());
876 if(phidiff1 > phidiff2) jet->SetFirstMother(7);
877 else jet->SetFirstMother(6);
878 //jet->Print();
879 if(fDebug > 1)
880 printf("AliMCAnalysisUtils::GetJets() - PYTHIA Jet %d: mother %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
881 i, jet->GetFirstMother(),jet->Pt(),jet->Energy(),jet->Phi()*TMath::RadToDeg(),jet->Eta());
882 fJetsList->Add(jet);
883 }
884 }//Pythia triggered jets
885 //HERWIG
886 else if (fMCGenerator=="HERWIG"){
887 Int_t pdg = -1;
888 //Check parton 1
889 TParticle * tmp = parton1;
890 if(parton1->GetPdgCode()!=22){
891 while(pdg != 94){
892 if(tmp->GetFirstDaughter()==-1) return fJetsList;
893 tmp = stack->Particle(tmp->GetFirstDaughter());
894 pdg = tmp->GetPdgCode();
895 }//while
896
897 //Add found jet to list
898 TParticle *jet1 = new TParticle(*tmp);
899 jet1->SetFirstMother(6);
900 fJetsList->Add(jet1);
901 //printf("jet 1: first daughter %d, last daughter %d\n", tmp->GetFirstDaughter(), tmp->GetLastDaughter());
902 //tmp = stack->Particle(tmp->GetFirstDaughter());
903 //tmp->Print();
904 //jet1->Print();
905 if(fDebug > 1)
906 printf("AliMCAnalysisUtils::GetJets() - HERWIG Jet 1: mother %d, status %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
907 jet1->GetFirstMother(),jet1->GetStatusCode(),jet1->Pt(),jet1->Energy(),jet1->Phi()*TMath::RadToDeg(),jet1->Eta());
908 }//not photon
909
910 //Check parton 2
911 pdg = -1;
912 tmp = parton2;
913 Int_t i = -1;
914 if(parton2->GetPdgCode()!=22){
915 while(pdg != 94){
916 if(tmp->GetFirstDaughter()==-1) return fJetsList;
917 i = tmp->GetFirstDaughter();
918 tmp = stack->Particle(tmp->GetFirstDaughter());
919 pdg = tmp->GetPdgCode();
920 }//while
921 //Add found jet to list
922 TParticle *jet2 = new TParticle(*tmp);
923 jet2->SetFirstMother(7);
924 fJetsList->Add(jet2);
925 //jet2->Print();
926 if(fDebug > 1)
927 printf("AliMCAnalysisUtils::GetJets() - HERWIG Jet 2: mother %d, status %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
928 jet2->GetFirstMother(),jet2->GetStatusCode(),jet2->Pt(),jet2->Energy(),jet2->Phi()*TMath::RadToDeg(),jet2->Eta());
929 //Int_t first = tmp->GetFirstDaughter();
930 //Int_t last = tmp->GetLastDaughter();
931 //printf("jet 2: first daughter %d, last daughter %d, pdg %d\n",first, last, tmp->GetPdgCode());
932 // for(Int_t d = first ; d < last+1; d++){
933// tmp = stack->Particle(d);
934// if(i == tmp->GetFirstMother())
935// printf("Daughter n %d, Mother %d, name %s, status %d, pT %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
936// d,tmp->GetFirstMother(), tmp->GetName(), tmp->GetStatusCode(),tmp->Pt(),tmp->Energy(),tmp->Phi()*TMath::RadToDeg(),tmp->Eta());
937// }
938 //tmp->Print();
939 }//not photon
940 }//Herwig generated jets
941 }
942
943 return fJetsList;
944}
945
946
947//________________________________________________________________
948void AliMCAnalysisUtils::Print(const Option_t * opt) const
949{
950 //Print some relevant parameters set for the analysis
951
952 if(! opt)
953 return;
954
955 printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
956
957 printf("Debug level = %d\n",fDebug);
958 printf("MC Generator = %s\n",fMCGenerator.Data());
959 printf(" \n");
960
961}
962
963