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