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