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