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