]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG/CaloTrackCorrBase/AliMCAnalysisUtils.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
avoid generator name string comparisons, change to int; avoid TLorentzVectors declara...
[u/mrichter/AliRoot.git] / PWG / CaloTrackCorrBase / 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),
0d360687 49fMCGenerator(kPythia),
50fMCGeneratorString("PYTHIA"),
51fDaughMom(), fDaughMom2(),
52fMotherMom(), fGMotherMom()
7cd4e982 53{
54 //Ctor
55}
7cd4e982 56
c5693f62 57//_______________________________________
7cd4e982 58AliMCAnalysisUtils::~AliMCAnalysisUtils()
59{
60 // Remove all pointers.
61
0d360687 62 if (fJetsList)
63 {
7cd4e982 64 fJetsList->Clear();
65 delete fJetsList ;
66 }
67}
68
c5693f62 69//_____________________________________________________________________________________________
8a2dbbff 70Int_t AliMCAnalysisUtils::CheckCommonAncestor(Int_t index1, Int_t index2,
c5693f62 71 const AliCaloTrackReader* reader,
72 Int_t & ancPDG, Int_t & ancStatus,
73 TLorentzVector & momentum, TVector3 & prodVertex)
840124fe 74{
0d360687 75 // Check the first common ancestor of 2 clusters, given the most likely labels of the primaries generating such clusters.
76
48a0baf4 77 Int_t label1[100];
78 Int_t label2[100];
79 label1[0]= index1;
80 label2[0]= index2;
81 Int_t counter1 = 0;
82 Int_t counter2 = 0;
840124fe 83
2644ead9 84 if(label1[0]==label2[0])
85 {
48a0baf4 86 //printf("AliMCAnalysisUtils::CheckCommonAncestor() - Already the same label: %d\n",label1[0]);
87 counter1=1;
88 counter2=1;
89 }
2644ead9 90 else
91 {
92 if(reader->ReadAODMCParticles())
93 {
94 TClonesArray * mcparticles = reader->GetAODMCParticles();
48a0baf4 95
96 Int_t label=label1[0];
2644ead9 97 if(label < 0) return -1;
98
36769d30 99 while(label > -1 && counter1 < 99)
100 {
48a0baf4 101 counter1++;
102 AliAODMCParticle * mom = (AliAODMCParticle *) mcparticles->At(label);
36769d30 103 if(mom)
104 {
840124fe 105 label = mom->GetMother() ;
106 label1[counter1]=label;
48a0baf4 107 }
108 //printf("\t counter %d, label %d\n", counter1,label);
109 }
36769d30 110
48a0baf4 111 //printf("Org label2=%d,\n",label2[0]);
112 label=label2[0];
2644ead9 113 if(label < 0) return -1;
36769d30 114
115 while(label > -1 && counter2 < 99)
116 {
48a0baf4 117 counter2++;
118 AliAODMCParticle * mom = (AliAODMCParticle *) mcparticles->At(label);
36769d30 119 if(mom)
120 {
48a0baf4 121 label = mom->GetMother() ;
122 label2[counter2]=label;
123 }
124 //printf("\t counter %d, label %d\n", counter2,label);
125 }
126 }//AOD MC
36769d30 127 else
128 { //Kine stack from ESDs
48a0baf4 129 AliStack * stack = reader->GetStack();
36769d30 130
48a0baf4 131 Int_t label=label1[0];
36769d30 132 while(label > -1 && counter1 < 99)
133 {
48a0baf4 134 counter1++;
135 TParticle * mom = stack->Particle(label);
36769d30 136 if(mom)
137 {
48a0baf4 138 label = mom->GetFirstMother() ;
139 label1[counter1]=label;
140 }
141 //printf("\t counter %d, label %d\n", counter1,label);
142 }
36769d30 143
48a0baf4 144 //printf("Org label2=%d,\n",label2[0]);
36769d30 145
48a0baf4 146 label=label2[0];
36769d30 147 while(label > -1 && counter2 < 99)
148 {
48a0baf4 149 counter2++;
150 TParticle * mom = stack->Particle(label);
36769d30 151 if(mom)
152 {
48a0baf4 153 label = mom->GetFirstMother() ;
154 label2[counter2]=label;
155 }
156 //printf("\t counter %d, label %d\n", counter2,label);
157 }
158 }// Kine stack from ESDs
159 }//First labels not the same
160
36769d30 161 if((counter1==99 || counter2==99) && fDebug >=0)
162 printf("AliMCAnalysisUtils::CheckCommonAncestor() - Genealogy too large c1: %d, c2= %d\n", counter1, counter2);
48a0baf4 163 //printf("CheckAncestor:\n");
36769d30 164
48a0baf4 165 Int_t commonparents = 0;
166 Int_t ancLabel = -1;
167 //printf("counters %d %d \n",counter1, counter2);
36769d30 168 for (Int_t c1 = 0; c1 < counter1; c1++)
169 {
170 for (Int_t c2 = 0; c2 < counter2; c2++)
171 {
172 if(label1[c1]==label2[c2] && label1[c1]>-1)
173 {
48a0baf4 174 ancLabel = label1[c1];
175 commonparents++;
36769d30 176
177 if(reader->ReadAODMCParticles())
178 {
2644ead9 179 AliAODMCParticle * mom = (AliAODMCParticle *) reader->GetAODMCParticles()->At(label1[c1]);
36769d30 180
181 if (mom)
182 {
48a0baf4 183 ancPDG = mom->GetPdgCode();
184 ancStatus = mom->GetStatus();
952615e5 185 momentum.SetPxPyPzE(mom->Px(),mom->Py(),mom->Pz(),mom->E());
c4a7d28a 186 prodVertex.SetXYZ(mom->Xv(),mom->Yv(),mom->Zv());
48a0baf4 187 }
188 }
36769d30 189 else
190 {
48a0baf4 191 TParticle * mom = (reader->GetStack())->Particle(label1[c1]);
36769d30 192
193 if (mom)
194 {
48a0baf4 195 ancPDG = mom->GetPdgCode();
196 ancStatus = mom->GetStatusCode();
952615e5 197 mom->Momentum(momentum);
c4a7d28a 198 prodVertex.SetXYZ(mom->Vx(),mom->Vy(),mom->Vz());
48a0baf4 199 }
200 }
36769d30 201
48a0baf4 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
36769d30 209 if(ancLabel < 0)
210 {
211 ancPDG = -10000;
212 ancStatus = -10000;
213 momentum.SetXYZT(0,0,0,0);
214 prodVertex.SetXYZ(-10,-10,-10);
215 }
216
48a0baf4 217 return ancLabel;
218}
219
3e57baab 220//________________________________________________________________________________________
8a2dbbff 221Int_t AliMCAnalysisUtils::CheckOrigin(const Int_t * label, Int_t nlabels,
3e57baab 222 const AliCaloTrackReader* reader, TString calorimeter)
840124fe 223{
0d360687 224 // Play with the montecarlo particles if available.
225
3c75ddf7 226 Int_t tag = 0;
227
840124fe 228 if(nlabels<=0) {
229 printf("AliMCAnalysisUtils::CheckOrigin(nlabel<=0) - No MC labels available, please check!!!\n");
230 return kMCBadLabel;
231 }
232
3e57baab 233 TObjArray* arrayCluster = 0;
234 if ( calorimeter == "EMCAL" ) arrayCluster = reader->GetEMCALClusters();
235 else if ( calorimeter == "PHOS" ) arrayCluster= reader->GetPHOSClusters();
236
3c75ddf7 237 //Select where the information is, ESD-galice stack or AOD mcparticles branch
238 if(reader->ReadStack()){
3e57baab 239 tag = CheckOriginInStack(label, nlabels, reader->GetStack(), arrayCluster);
3c75ddf7 240 }
241 else if(reader->ReadAODMCParticles()){
3e57baab 242 tag = CheckOriginInAOD(label, nlabels, reader->GetAODMCParticles(),arrayCluster);
3c75ddf7 243 }
244
245 return tag ;
840124fe 246}
247
3e57baab 248//____________________________________________________________________________________________________
249Int_t AliMCAnalysisUtils::CheckOrigin(Int_t label, const AliCaloTrackReader* reader, TString calorimeter)
840124fe 250{
0d360687 251 // Play with the montecarlo particles if available.
252
3c75ddf7 253 Int_t tag = 0;
b5ef1905 254
255 if(label<0) {
898c9d44 256 printf("AliMCAnalysisUtils::CheckOrigin(label<0) - No MC labels available, please check!!!\n");
257 return kMCBadLabel;
b5ef1905 258 }
259
3e57baab 260 TObjArray* arrayCluster = 0;
261 if ( calorimeter == "EMCAL" ) arrayCluster = reader->GetEMCALClusters();
262 else if ( calorimeter == "PHOS" ) arrayCluster= reader->GetPHOSClusters();
263
3c75ddf7 264 Int_t labels[]={label};
265
266 //Select where the information is, ESD-galice stack or AOD mcparticles branch
267 if(reader->ReadStack()){
3e57baab 268 tag = CheckOriginInStack(labels, 1,reader->GetStack(),arrayCluster);
3c75ddf7 269 }
270 else if(reader->ReadAODMCParticles()){
3e57baab 271 tag = CheckOriginInAOD(labels, 1,reader->GetAODMCParticles(),arrayCluster);
3c75ddf7 272 }
273
274 return tag ;
7cd4e982 275}
276
3e57baab 277//__________________________________________________________________________________________
278Int_t AliMCAnalysisUtils::CheckOriginInStack(const Int_t *labels, Int_t nlabels,
279 AliStack* stack, const TObjArray* arrayCluster)
840124fe 280{
7cd4e982 281 // Play with the MC stack if available. Tag particles depending on their origin.
282 // Do same things as in CheckOriginInAOD but different input.
556e55b0 283
284 //generally speaking, label is the MC label of a reconstructed
285 //entity (track, cluster, etc) for which we want to know something
286 //about its heritage, but one can also use it directly with stack
287 //particles not connected to reconstructed entities
f8006433 288
36769d30 289 if(!stack)
290 {
1cd71065 291 if (fDebug >=0)
f8006433 292 printf("AliMCAnalysisUtils::CheckOriginInStack() - Stack is not available, check analysis settings in configuration file, STOP!!\n");
293 return -1;
7cd4e982 294 }
f8006433 295
7cd4e982 296 Int_t tag = 0;
902aa95c 297 Int_t label=labels[0];//Most significant particle contributing to the cluster
3c75ddf7 298
36769d30 299 if(label >= 0 && label < stack->GetNtrack())
300 {
556e55b0 301 //MC particle of interest is the "mom" of the entity
7cd4e982 302 TParticle * mom = stack->Particle(label);
eb3e1b50 303 Int_t iMom = label;
304 Int_t mPdgSign = mom->GetPdgCode();
305 Int_t mPdg = TMath::Abs(mPdgSign);
306 Int_t mStatus = mom->GetStatusCode() ;
307 Int_t iParent = mom->GetFirstMother() ;
36769d30 308
0d360687 309 if(fDebug > 0 && label < 8 && fMCGenerator != kBoxLike) printf("AliMCAnalysisUtils::CheckOriginInStack() - Mother is parton %d\n",iParent);
7cd4e982 310
556e55b0 311 //GrandParent of the entity
d7c10d78 312 TParticle * parent = NULL;
7cd4e982 313 Int_t pPdg = -1;
314 Int_t pStatus =-1;
36769d30 315 if(iParent >= 0)
316 {
7cd4e982 317 parent = stack->Particle(iParent);
36769d30 318
319 if(parent)
320 {
95fb4ab7 321 pPdg = TMath::Abs(parent->GetPdgCode());
322 pStatus = parent->GetStatusCode();
323 }
7cd4e982 324 }
325 else if(fDebug > 0 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - Parent with label %d\n",iParent);
8dacfd76 326
36769d30 327 if(fDebug > 2 )
328 {
3c75ddf7 329 printf("AliMCAnalysisUtils::CheckOriginInStack() - Cluster most contributing mother and its parent: \n");
330 printf("\t Mother label %d, pdg %d, status %d\n",iMom, mPdg, mStatus);
331 printf("\t Parent label %d, pdg %d, status %d\n",iParent, pPdg, pStatus);
f8006433 332 }
902aa95c 333
556e55b0 334 //Check if "mother" of entity is converted, if not, get the first non converted mother
36769d30 335 if((mPdg == 22 || mPdg == 11) && (pPdg == 22 || pPdg == 11) && mStatus == 0)
336 {
8dacfd76 337 SetTagBit(tag,kMCConversion);
36769d30 338
8dacfd76 339 //Check if the mother is photon or electron with status not stable
36769d30 340 while ((pPdg == 22 || pPdg == 11) && mStatus != 1)
341 {
f8006433 342 //Mother
eb3e1b50 343 iMom = mom->GetFirstMother();
344 mom = stack->Particle(iMom);
345 mPdgSign = mom->GetPdgCode();
346 mPdg = TMath::Abs(mPdgSign);
347 mStatus = mom->GetStatusCode() ;
348 iParent = mom->GetFirstMother() ;
36769d30 349
f8006433 350 if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - Mother is parton %d\n",iParent);
351
352 //GrandParent
36769d30 353 if(iParent >= 0)
354 {
f8006433 355 parent = stack->Particle(iParent);
36769d30 356 if(parent)
357 {
95fb4ab7 358 pPdg = TMath::Abs(parent->GetPdgCode());
359 pStatus = parent->GetStatusCode();
360 }
f8006433 361 }
0e453907 362 else {// in case of gun/box simulations
363 pPdg = 0;
364 pStatus = 0;
365 break;
366 }
36769d30 367 }//while
368
369 if(fDebug > 2 )
370 {
f8006433 371 printf("AliMCAnalysisUtils::CheckOriginInStack() - Converted photon/electron: \n");
372 printf("\t Mother label %d, pdg %d, status %d\n",iMom, mPdg, mStatus);
373 printf("\t Parent label %d, pdg %d, status %d\n",iParent, pPdg, pStatus);
374 }
375
7cd4e982 376 }//mother and parent are electron or photon and have status 0
36769d30 377 else if((mPdg == 22 || mPdg == 11) && mStatus == 0)
378 {
902aa95c 379 //Still a conversion but only one electron/photon generated. Just from hadrons but not decays.
380 if(pPdg == 2112 || pPdg == 211 || pPdg == 321 ||
36769d30 381 pPdg == 2212 || pPdg == 130 || pPdg == 13 )
382 {
f8006433 383 SetTagBit(tag,kMCConversion);
eb3e1b50 384 iMom = mom->GetFirstMother();
385 mom = stack->Particle(iMom);
386 mPdgSign = mom->GetPdgCode();
387 mPdg = TMath::Abs(mPdgSign);
f8006433 388
36769d30 389 if(fDebug > 2 )
390 {
f8006433 391 printf("AliMCAnalysisUtils::CheckOriginInStack() - Converted hadron: \n");
392 printf("\t Mother label %d, pdg %d, status %d\n",iMom, mPdg, mStatus);
393 }
394 }//hadron converted
395
8dacfd76 396 //Comment for the next lines, we do not check the parent of the hadron for the moment.
397 //iParent = mom->GetFirstMother() ;
398 //if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - Mother is parton %d\n",iParent);
399
400 //GrandParent
401 //if(iParent >= 0){
402 // parent = stack->Particle(iParent);
403 // pPdg = TMath::Abs(parent->GetPdgCode());
404 //}
405 }
902aa95c 406 // conversion into electrons/photons checked
8dacfd76 407
7cd4e982 408 //first check for typical charged particles
eb3e1b50 409 if (mPdg == 13) SetTagBit(tag,kMCMuon);
410 else if(mPdg == 211) SetTagBit(tag,kMCPion);
411 else if(mPdg == 321) SetTagBit(tag,kMCKaon);
412 else if(mPdgSign == 2212) SetTagBit(tag,kMCProton);
413 else if(mPdgSign == -2212) SetTagBit(tag,kMCAntiProton);
414 else if(mPdgSign == 2112) SetTagBit(tag,kMCNeutron);
415 else if(mPdgSign == -2112) SetTagBit(tag,kMCAntiNeutron);
840124fe 416
8dacfd76 417 //check for pi0 and eta (shouldn't happen unless their decays were turned off)
36769d30 418 else if(mPdg == 111)
419 {
4d52b003 420
f8006433 421 SetTagBit(tag,kMCPi0Decay);
36769d30 422
4d52b003 423 if(fDebug > 2 )
424 printf("AliMCAnalysisUtils::CheckOriginInStack() - First mother is directly pi0, not decayed by generator \n");
36769d30 425
f8006433 426 CheckOverlapped2GammaDecay(labels,nlabels, iMom, stack, tag); //set to kMCPi0 if 2 gammas in same cluster
3e57baab 427
f8006433 428 }
36769d30 429 else if(mPdg == 221)
430 {
f8006433 431 SetTagBit(tag,kMCEtaDecay);
36769d30 432
4d52b003 433 if(fDebug > 2 )
434 printf("AliMCAnalysisUtils::CheckOriginInStack() - First mother is directly eta, not decayed by generator \n");
36769d30 435
f8006433 436 CheckOverlapped2GammaDecay(labels,nlabels, iMom, stack, tag); //set to kMCEta if 2 gammas in same cluster
437 }
8dacfd76 438 //Photons
36769d30 439 else if(mPdg == 22)
440 {
7cd4e982 441 SetTagBit(tag,kMCPhoton);
4d52b003 442
443 if(pPdg == 111)
444 {
445 SetTagBit(tag,kMCPi0Decay);
446
447 if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - PYTHIA pi0 decay photon, parent pi0 with status 11 \n");
448
449 CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag); //set to kMCPi0 if 2 gammas in same cluster
3e57baab 450 // In case it did not merge, check if the decay companion is lost
451 if(!CheckTagBit(tag, kMCPi0) && !CheckTagBit(tag,kMCDecayPairInCalo))
452 CheckLostDecayPair(arrayCluster,iMom, iParent, stack, tag);
453
454 //printf("Bit set is Merged %d, Pair in calo %d, Lost %d\n",CheckTagBit(tag, kMCPi0),CheckTagBit(tag,kMCDecayPairInCalo),CheckTagBit(tag,kMCDecayPairLost));
4d52b003 455 }
456 else if (pPdg == 221)
457 {
458 SetTagBit(tag, kMCEtaDecay);
459
460 if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - PYTHIA eta decay photon, parent pi0 with status 11 \n");
461
462 CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag);//set to kMCEta if 2 gammas in same cluster
3e57baab 463 // In case it did not merge, check if the decay companion is lost
464 if(!CheckTagBit(tag, kMCEta) && !CheckTagBit(tag,kMCDecayPairInCalo))
465 CheckLostDecayPair(arrayCluster,iMom, iParent, stack, tag);
4d52b003 466 }
467 else if(mStatus == 1)
36769d30 468 { //undecayed particle
0d360687 469 if(fMCGenerator == kPythia)
36769d30 470 {
471 if(iParent < 8 && iParent > 5)
472 {//outgoing partons
f8006433 473 if(pPdg == 22) SetTagBit(tag,kMCPrompt);
36769d30 474 else SetTagBit(tag,kMCFragmentation);
f8006433 475 }//Outgoing partons
36769d30 476 else if(iParent <= 5)
477 {
f8006433 478 SetTagBit(tag, kMCISR); //Initial state radiation
479 }
4d52b003 480 else SetTagBit(tag,kMCUnknown);
481 }//PYTHIA
482
0d360687 483 else if(fMCGenerator == kHerwig)
36769d30 484 {
485 if(pStatus < 197)
486 {//Not decay
487 while(1)
488 {
489 if(parent)
490 {
95fb4ab7 491 if(parent->GetFirstMother()<=5) break;
492 iParent = parent->GetFirstMother();
493 parent=stack->Particle(iParent);
494 pStatus= parent->GetStatusCode();
495 pPdg = TMath::Abs(parent->GetPdgCode());
496 } else break;
f8006433 497 }//Look for the parton
498
36769d30 499 if(iParent < 8 && iParent > 5)
500 {
f8006433 501 if(pPdg == 22) SetTagBit(tag,kMCPrompt);
36769d30 502 else SetTagBit(tag,kMCFragmentation);
f8006433 503 }
504 else SetTagBit(tag,kMCISR);//Initial state radiation
505 }//Not decay
4d52b003 506 else SetTagBit(tag,kMCUnknown);
f8006433 507 }//HERWIG
4d52b003 508 }
509 else SetTagBit(tag,kMCOtherDecay);
510
7cd4e982 511 }//Mother Photon
8dacfd76 512
7cd4e982 513 //Electron check. Where did that electron come from?
514 else if(mPdg == 11){ //electron
36769d30 515 if(pPdg == 11 && parent)
516 {
e495fc0d 517 Int_t iGrandma = parent->GetFirstMother();
36769d30 518 if(iGrandma >= 0)
519 {
e495fc0d 520 TParticle* gma = (TParticle*)stack->Particle(iGrandma); //get mother
521 Int_t gPdg = TMath::Abs(gma->GetPdgCode());
f8006433 522
36769d30 523 if (gPdg == 23) { SetTagBit(tag,kMCZDecay); } //parent is Z-boson
e495fc0d 524 else if (gPdg == 24) { SetTagBit(tag,kMCWDecay); } //parent is W-boson
525 }
526 }
36769d30 527
528 SetTagBit(tag,kMCElectron);
529
902aa95c 530 if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInStack() - Checking ancestors of electrons\n");
36769d30 531
3e57baab 532 if (pPdg == 111) { SetTagBit(tag,kMCPi0Decay); SetTagBit(tag, kMCDecayDalitz) ; } //Pi0 Dalitz decay
533 else if (pPdg == 221) { SetTagBit(tag,kMCEtaDecay); SetTagBit(tag, kMCDecayDalitz) ; } //Eta Dalitz decay
d409d6b9 534 else if((499 < pPdg && pPdg < 600)||(4999 < pPdg && pPdg < 6000)) { SetTagBit(tag,kMCEFromB); } //b-->e decay
b3fdfed5 535 else if((399 < pPdg && pPdg < 500)||(3999 < pPdg && pPdg < 5000)) { //check charm decay
36769d30 536 if(parent)
537 {
b3fdfed5 538 Int_t iGrandma = parent->GetFirstMother();
36769d30 539 if(iGrandma >= 0)
540 {
b3fdfed5 541 TParticle* gma = (TParticle*)stack->Particle(iGrandma); //get mother of charm
542 Int_t gPdg = TMath::Abs(gma->GetPdgCode());
543 if((499 < gPdg && gPdg < 600)||(4999 < gPdg && gPdg < 6000)) SetTagBit(tag,kMCEFromCFromB); //b-->c-->e
544 else SetTagBit(tag,kMCEFromC); //c-->e
36769d30 545 }
546 else SetTagBit(tag,kMCEFromC); //c-->e
b3fdfed5 547 }//parent
36769d30 548 }
549 else
550 {
f8006433 551 //if it is not from any of the above, where is it from?
552 if(pPdg > 10000) SetTagBit(tag,kMCUnknown);
36769d30 553
f8006433 554 else SetTagBit(tag,kMCOtherDecay);
36769d30 555
95fb4ab7 556 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 557 }
7cd4e982 558 }//electron check
559 //Cluster was made by something else
36769d30 560 else
561 {
562 if(fDebug > 0)
563 printf("AliMCAnalysisUtils::CheckOriginInStack() - \tSetting kMCUnknown for cluster from %s (pdg = %d, Parent pdg = %d)\n",
564 mom->GetName(),mPdg,pPdg);
565
7cd4e982 566 SetTagBit(tag,kMCUnknown);
567 }
568 }//Good label value
36769d30 569 else
570 {// Bad label
902aa95c 571 if(label < 0 && (fDebug >= 0))
f8006433 572 printf("AliMCAnalysisUtils::CheckOriginInStack() *** bad label or no stack ***: label %d \n", label);
36769d30 573
574 if(label >= stack->GetNtrack() && (fDebug >= 0))
f8006433 575 printf("AliMCAnalysisUtils::CheckOriginInStack() *** large label ***: label %d, n tracks %d \n", label, stack->GetNtrack());
36769d30 576
7cd4e982 577 SetTagBit(tag,kMCUnknown);
578 }//Bad label
579
580 return tag;
581
582}
583
584
3e57baab 585//________________________________________________________________________________________________________
8a2dbbff 586Int_t AliMCAnalysisUtils::CheckOriginInAOD(const Int_t *labels, Int_t nlabels,
3e57baab 587 const TClonesArray *mcparticles, const TObjArray* arrayCluster)
840124fe 588{
7cd4e982 589 // Play with the MCParticles in AOD if available. Tag particles depending on their origin.
590 // Do same things as in CheckOriginInStack but different input.
0d360687 591
36769d30 592 if(!mcparticles)
593 {
1cd71065 594 if(fDebug >= 0)
f8006433 595 printf("AliMCAnalysisUtils::CheckOriginInAOD() - AODMCParticles is not available, check analysis settings in configuration file!!\n");
596 return -1;
7cd4e982 597 }
7cd4e982 598
599 Int_t tag = 0;
902aa95c 600 Int_t label=labels[0];//Most significant particle contributing to the cluster
f8006433 601
7cd4e982 602 Int_t nprimaries = mcparticles->GetEntriesFast();
36769d30 603 if(label >= 0 && label < nprimaries)
604 {
7cd4e982 605 //Mother
606 AliAODMCParticle * mom = (AliAODMCParticle *) mcparticles->At(label);
eb3e1b50 607 Int_t iMom = label;
608 Int_t mPdgSign = mom->GetPdgCode();
609 Int_t mPdg = TMath::Abs(mPdgSign);
610 Int_t iParent = mom->GetMother() ;
36769d30 611
0d360687 612 if(fDebug > 0 && label < 8 && fMCGenerator != kBoxLike) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Mother is parton %d\n",iParent);
7cd4e982 613
614 //GrandParent
95fb4ab7 615 AliAODMCParticle * parent = NULL ;
7cd4e982 616 Int_t pPdg = -1;
36769d30 617 if(iParent >= 0)
618 {
7cd4e982 619 parent = (AliAODMCParticle *) mcparticles->At(iParent);
620 pPdg = TMath::Abs(parent->GetPdgCode());
621 }
622 else if(fDebug > 0 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Parent with label %d\n",iParent);
f8006433 623
36769d30 624 if(fDebug > 2 )
625 {
3c75ddf7 626 printf("AliMCAnalysisUtils::CheckOriginInAOD() - Cluster most contributing mother and its parent: \n");
36769d30 627
628 printf("\t Mother label %d, pdg %d, Primary? %d, Physical Primary? %d\n",
629 iMom, mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary());
630
3c75ddf7 631 if(parent)
36769d30 632 printf("\t Parent label %d, pdg %d, Primary? %d, Physical Primary? %d\n",
633 iParent, pPdg, parent->IsPrimary(), parent->IsPhysicalPrimary());
f8006433 634 }
902aa95c 635
8dacfd76 636 //Check if mother is converted, if not, get the first non converted mother
36769d30 637 if((mPdg == 22 || mPdg == 11) && (pPdg == 22 || pPdg == 11) && !mom->IsPrimary())
638 {
8dacfd76 639 SetTagBit(tag,kMCConversion);
36769d30 640
8dacfd76 641 //Check if the mother is photon or electron with status not stable
36769d30 642 while ((pPdg == 22 || pPdg == 11) && !mom->IsPhysicalPrimary())
643 {
f8006433 644 //Mother
eb3e1b50 645 iMom = mom->GetMother();
646 mom = (AliAODMCParticle *) mcparticles->At(iMom);
647 mPdgSign = mom->GetPdgCode();
648 mPdg = TMath::Abs(mPdgSign);
649 iParent = mom->GetMother() ;
f8006433 650 if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Mother is parton %d\n",iParent);
651
652 //GrandParent
36769d30 653 if(iParent >= 0 && parent)
654 {
f8006433 655 parent = (AliAODMCParticle *) mcparticles->At(iParent);
656 pPdg = TMath::Abs(parent->GetPdgCode());
657 }
658 // printf("\t While Mother label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iMom, mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary());
659 // printf("\t While Parent label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iParent, pPdg, parent->IsPrimary(), parent->IsPhysicalPrimary());
660
902aa95c 661 }//while
f8006433 662
36769d30 663 if(fDebug > 2 )
664 {
f8006433 665 printf("AliMCAnalysisUtils::CheckOriginInAOD() - Converted photon/electron : \n");
36769d30 666
667 printf("\t Mother label %d, pdg %d, Primary? %d, Physical Primary? %d\n",
668 iMom, mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary());
669
95fb4ab7 670 if(parent)
36769d30 671 printf("\t Parent label %d, pdg %d, Primary? %d, Physical Primary? %d\n",
672 iParent, pPdg, parent->IsPrimary(), parent->IsPhysicalPrimary());
f8006433 673 }
674
8dacfd76 675 }//mother and parent are electron or photon and have status 0 and parent is photon or electron
36769d30 676 else if((mPdg == 22 || mPdg == 11) && !mom->IsPrimary())
677 {
8dacfd76 678 //Still a conversion but only one electron/photon generated. Just from hadrons
902aa95c 679 if(pPdg == 2112 || pPdg == 211 || pPdg == 321 ||
36769d30 680 pPdg == 2212 || pPdg == 130 || pPdg == 13 )
681 {
f8006433 682 SetTagBit(tag,kMCConversion);
eb3e1b50 683 iMom = mom->GetMother();
684 mom = (AliAODMCParticle *) mcparticles->At(iMom);
685 mPdgSign = mom->GetPdgCode();
686 mPdg = TMath::Abs(mPdgSign);
f8006433 687
36769d30 688 if(fDebug > 2 )
689 {
f8006433 690 printf("AliMCAnalysisUtils::CheckOriginInAOD() - Converted hadron : \n");
691 printf("\t Mother label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iMom, mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary());
692 }
693 }//hadron converted
694
8dacfd76 695 //Comment for next lines, we do not check the parent of the hadron for the moment.
696 //iParent = mom->GetMother() ;
697 //if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Mother is parton %d\n",iParent);
698
699 //GrandParent
700 //if(iParent >= 0){
701 // parent = (AliAODMCParticle *) mcparticles->At(iParent);
702 // pPdg = TMath::Abs(parent->GetPdgCode());
703 //}
704 }
705
3c75ddf7 706 //printf("Final mother mPDG %d\n",mPdg);
c5693f62 707
8dacfd76 708 // conversion into electrons/photons checked
709
710 //first check for typical charged particles
eb3e1b50 711 if (mPdg == 13) SetTagBit(tag,kMCMuon);
712 else if(mPdg == 211) SetTagBit(tag,kMCPion);
713 else if(mPdg == 321) SetTagBit(tag,kMCKaon);
714 else if(mPdgSign == 2212) SetTagBit(tag,kMCProton);
715 else if(mPdgSign == 2112) SetTagBit(tag,kMCNeutron);
716 else if(mPdgSign == -2212) SetTagBit(tag,kMCAntiProton);
717 else if(mPdgSign == -2112) SetTagBit(tag,kMCAntiNeutron);
718
8dacfd76 719 //check for pi0 and eta (shouldn't happen unless their decays were turned off)
36769d30 720 else if(mPdg == 111)
721 {
f8006433 722 SetTagBit(tag,kMCPi0Decay);
36769d30 723
f8006433 724 if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - First mother is directly pi0, not decayed by generator \n");
36769d30 725
f8006433 726 CheckOverlapped2GammaDecay(labels,nlabels, iMom, mcparticles, tag); //set to kMCPi0 if 2 gammas in same cluster
727 }
36769d30 728 else if(mPdg == 221)
729 {
f8006433 730 SetTagBit(tag,kMCEtaDecay);
36769d30 731
f8006433 732 if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - First mother is directly eta, not decayed by generator \n");
36769d30 733
f8006433 734 CheckOverlapped2GammaDecay(labels,nlabels, iMom, mcparticles, tag); //set to kMCEta if 2 gammas in same cluster
735 }
8dacfd76 736 //Photons
36769d30 737 else if(mPdg == 22)
738 {
7cd4e982 739 SetTagBit(tag,kMCPhoton);
4d52b003 740
741 if(pPdg == 111)
742 {
743 SetTagBit(tag,kMCPi0Decay);
744
745 if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Generator pi0 decay photon \n");
746
747 CheckOverlapped2GammaDecay(labels,nlabels, iParent, mcparticles, tag); //set to kMCPi0 if 2 gammas in same cluster
3e57baab 748 // In case it did not merge, check if the decay companion is lost
749 if(!CheckTagBit(tag, kMCPi0) && !CheckTagBit(tag,kMCDecayPairInCalo) && !CheckTagBit(tag,kMCDecayPairLost))
750 CheckLostDecayPair(arrayCluster,iMom, iParent, mcparticles, tag);
751
752 //printf("Bit set is Merged %d, Pair in calo %d, Lost %d\n",CheckTagBit(tag, kMCPi0),CheckTagBit(tag,kMCDecayPairInCalo),CheckTagBit(tag,kMCDecayPairLost));
4d52b003 753 }
754 else if (pPdg == 221)
755 {
756 SetTagBit(tag, kMCEtaDecay);
757
758 if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Generator eta decay photon \n");
759
760 CheckOverlapped2GammaDecay(labels,nlabels, iParent, mcparticles, tag); //set to kMCEta if 2 gammas in same cluster
3e57baab 761 // In case it did not merge, check if the decay companion is lost
762 if(!CheckTagBit(tag, kMCEta) && !CheckTagBit(tag,kMCDecayPairInCalo))
763 CheckLostDecayPair(arrayCluster,iMom, iParent, mcparticles, tag);
4d52b003 764 }
0d360687 765 else if( mom->IsPhysicalPrimary() && ( fMCGenerator == kPythia || fMCGenerator == kHerwig ) ) //undecayed particle
3c75ddf7 766 {
36769d30 767 if(iParent < 8 && iParent > 5 )
768 {//outgoing partons
3c75ddf7 769 if(pPdg == 22) SetTagBit(tag,kMCPrompt);
770 else SetTagBit(tag,kMCFragmentation);
771 }//Outgoing partons
0d360687 772 else if( iParent <= 5 && ( fMCGenerator == kPythia || fMCGenerator == kHerwig ) )
36769d30 773 {
3c75ddf7 774 SetTagBit(tag, kMCISR); //Initial state radiation
775 }
4d52b003 776 else SetTagBit(tag,kMCUnknown);
840124fe 777 }//Physical primary
4d52b003 778 else SetTagBit(tag,kMCOtherDecay);
779
7cd4e982 780 }//Mother Photon
781
782 //Electron check. Where did that electron come from?
36769d30 783 else if(mPdg == 11)
784 { //electron
785 if(pPdg == 11 && parent)
786 {
e495fc0d 787 Int_t iGrandma = parent->GetMother();
36769d30 788 if(iGrandma >= 0)
789 {
e495fc0d 790 AliAODMCParticle* gma = (AliAODMCParticle*)mcparticles->At(iGrandma);
791 Int_t gPdg = TMath::Abs(gma->GetPdgCode());
f8006433 792
36769d30 793 if (gPdg == 23) { SetTagBit(tag,kMCZDecay); } //parent is Z-boson
e495fc0d 794 else if (gPdg == 24) { SetTagBit(tag,kMCWDecay); } //parent is W-boson
795 }
796 }
36769d30 797
798 SetTagBit(tag,kMCElectron);
799
3e57baab 800 if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Checking ancestors of electrons");
801
802 if (pPdg == 111) { SetTagBit(tag,kMCPi0Decay); SetTagBit(tag,kMCDecayDalitz);} //Pi0 Dalitz decay
803 else if (pPdg == 221) { SetTagBit(tag,kMCEtaDecay); SetTagBit(tag,kMCDecayDalitz);} //Eta Dalitz decay
d409d6b9 804 else if((499 < pPdg && pPdg < 600)||(4999 < pPdg && pPdg < 6000)) { SetTagBit(tag,kMCEFromB);} //b-hadron decay
36769d30 805 else if((399 < pPdg && pPdg < 500)||(3999 < pPdg && pPdg < 5000))
806 { //c-hadron decay check
807 if(parent)
808 {
b3fdfed5 809 Int_t iGrandma = parent->GetMother();
36769d30 810 if(iGrandma >= 0)
811 {
b3fdfed5 812 AliAODMCParticle* gma = (AliAODMCParticle*)mcparticles->At(iGrandma); //charm's mother
813 Int_t gPdg = TMath::Abs(gma->GetPdgCode());
814 if((499 < gPdg && gPdg < 600)||(4999 < gPdg && gPdg < 6000)) SetTagBit(tag,kMCEFromCFromB); //b-->c-->e decay
815 else SetTagBit(tag,kMCEFromC); //c-hadron decay
36769d30 816 }
817 else SetTagBit(tag,kMCEFromC); //c-hadron decay
b3fdfed5 818 }//parent
36769d30 819 } else
820 { //prompt or other decay
f8006433 821 TParticlePDG* foo = TDatabasePDG::Instance()->GetParticle(pPdg);
822 TParticlePDG* foo1 = TDatabasePDG::Instance()->GetParticle(mPdg);
36769d30 823
824 if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Electron from other origin: %s (pPdg = %d) %s (mPdg = %d)\n",
825 foo->GetName(), pPdg,foo1->GetName(),mPdg);
826
f8006433 827 if(pPdg > 10000) SetTagBit(tag,kMCUnknown);
36769d30 828 else SetTagBit(tag,kMCOtherDecay);
d409d6b9 829 }
8dacfd76 830 }//electron check
7cd4e982 831 //cluster was made by something else
36769d30 832 else
833 {
834 if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInAOD() - \tSetting kMCUnknown for cluster with pdg = %d, Parent pdg = %d\n",
835 mPdg,pPdg);
7cd4e982 836 SetTagBit(tag,kMCUnknown);
837 }
838 }//Good label value
36769d30 839 else
840 {//Bad label
902aa95c 841 if(label < 0 && (fDebug >= 0) )
f8006433 842 printf("AliMCAnalysisUtils::CheckOriginInAOD() *** bad label or no mcparticles ***: label %d \n", label);
36769d30 843
844 if(label >= mcparticles->GetEntriesFast() && (fDebug >= 0) )
f8006433 845 printf("AliMCAnalysisUtils::CheckOriginInAOD() *** large label ***: label %d, n tracks %d \n", label, mcparticles->GetEntriesFast());
36769d30 846
7cd4e982 847 SetTagBit(tag,kMCUnknown);
f8006433 848
7cd4e982 849 }//Bad label
850
851 return tag;
852
853}
854
8a2dbbff 855//_________________________________________________________________________________________
856void AliMCAnalysisUtils::CheckOverlapped2GammaDecay(const Int_t *labels, Int_t nlabels,
857 Int_t mesonIndex, AliStack *stack,
c5693f62 858 Int_t &tag)
840124fe 859{
0d360687 860 // Check if cluster is formed from the contribution of 2 decay photons from pi0 or eta. Input in stack.
3c75ddf7 861
862 if(labels[0] < 0 || labels[0] > stack->GetNtrack() || nlabels <= 1) {
36769d30 863 if(fDebug > 2) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) - Exit : label[0] %d, n primaries %d, nlabels %d \n",
f8006433 864 labels[0],stack->GetNtrack(), nlabels);
3c75ddf7 865 return;
866 }
867
868 TParticle * meson = stack->Particle(mesonIndex);
869 Int_t mesonPdg = meson->GetPdgCode();
36769d30 870 if(mesonPdg!=111 && mesonPdg!=221)
871 {
872 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) - Wrong pi0/eta PDG : %d \n",mesonPdg);
3c75ddf7 873 return;
874 }
875
36769d30 876 if(fDebug > 2) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) - %s, label %d\n",meson->GetName(), mesonIndex);
3c75ddf7 877
878 //Check if meson decayed into 2 daughters or if both were kept.
36769d30 879 if(meson->GetNDaughters() != 2)
880 {
3c75ddf7 881 if(fDebug > 2)
36769d30 882 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) - Not overalapped. Number of daughters is %d, not 2 \n",meson->GetNDaughters());
3c75ddf7 883 return;
884 }
885
886 //Get the daughters
887 Int_t iPhoton0 = meson->GetDaughter(0);
888 Int_t iPhoton1 = meson->GetDaughter(1);
889 TParticle *photon0 = stack->Particle(iPhoton0);
890 TParticle *photon1 = stack->Particle(iPhoton1);
891
892 //Check if both daughters are photons
36769d30 893 if(photon0->GetPdgCode() != 22 || photon1->GetPdgCode()!=22)
894 {
3c75ddf7 895 if(fDebug > 2)
36769d30 896 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) - Not overalapped. PDG: daughter 1 = %d, of daughter 2 = %d \n",photon0->GetPdgCode(),photon1->GetPdgCode());
3c75ddf7 897 return;
898 }
899
900 if(fDebug > 2)
36769d30 901 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) - Daughter labels : photon0 = %d, photon1 = %d \n",iPhoton0,iPhoton1);
3c75ddf7 902
903 //Check if both photons contribute to the cluster
904 Bool_t okPhoton0 = kFALSE;
905 Bool_t okPhoton1 = kFALSE;
906
36769d30 907 if(fDebug > 3) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) - Labels loop:\n");
908
909 Bool_t conversion = kFALSE;
3c75ddf7 910
36769d30 911 for(Int_t i = 0; i < nlabels; i++)
912 {
3c75ddf7 913 if(fDebug > 3) printf("\t at begin:label %d/%d: %d, ok? photon1 %d, photon2 %d\n", i+1, nlabels, labels[i], okPhoton0, okPhoton1);
914
915 //If we already found both, break the loop
916 if(okPhoton0 && okPhoton1) break;
917
36769d30 918 Int_t index = labels[i];
919 if (iPhoton0 == index)
920 {
3c75ddf7 921 okPhoton0 = kTRUE;
922 continue;
923 }
36769d30 924 else if (iPhoton1 == index)
925 {
3c75ddf7 926 okPhoton1 = kTRUE;
927 continue;
928 }
929
930 //Trace back the mother in case it was a conversion
931
36769d30 932 if(index >= stack->GetNtrack())
933 {
934 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) Particle index %d larger than size of list %d !!\n",
935 index,stack->GetNtrack());
3c75ddf7 936 continue;
937 }
938
939 TParticle * daught = stack->Particle(index);
940 Int_t tmpindex = daught->GetFirstMother();
941 if(fDebug > 3) printf("\t Conversion? : mother %d\n",tmpindex);
36769d30 942 while(tmpindex>=0)
943 {
3c75ddf7 944 //MC particle of interest is the mother
945 if(fDebug > 3) printf("\t \t parent index %d\n",tmpindex);
946 daught = stack->Particle(tmpindex);
36769d30 947 if (iPhoton0 == tmpindex)
948 {
949 conversion = kTRUE;
950 okPhoton0 = kTRUE;
c5693f62 951 break;
3c75ddf7 952 }
36769d30 953 else if (iPhoton1 == tmpindex)
954 {
955 conversion = kTRUE;
956 okPhoton1 = kTRUE;
c5693f62 957 break;
3c75ddf7 958 }
36769d30 959
3c75ddf7 960 tmpindex = daught->GetFirstMother();
36769d30 961
3c75ddf7 962 }//While to check if pi0/eta daughter was one of these contributors to the cluster
963
964 if(i == 0 && (!okPhoton0 && !okPhoton1) && fDebug>=0)
36769d30 965 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) - Something happens, first label should be from a photon decay!\n");
3c75ddf7 966
967 }//loop on list of labels
968
969 //If both photons contribute tag as the corresponding meson.
36769d30 970 if(okPhoton0 && okPhoton1)
971 {
3c75ddf7 972 if(fDebug > 2)
36769d30 973 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) - %s OVERLAPPED DECAY \n", meson->GetName());
974
975 if(!CheckTagBit(tag,kMCConversion) && conversion) SetTagBit(tag,kMCConversion) ;
3c75ddf7 976
977 if(mesonPdg == 111) SetTagBit(tag,kMCPi0);
36769d30 978 else SetTagBit(tag,kMCEta);
3c75ddf7 979 }
980
902aa95c 981}
982
8a2dbbff 983//________________________________________________________________________________________________________
984void AliMCAnalysisUtils::CheckOverlapped2GammaDecay(const Int_t *labels, Int_t nlabels, Int_t mesonIndex,
985 const TClonesArray *mcparticles, Int_t & tag )
840124fe 986{
0d360687 987 // Check if cluster is formed from the contribution of 2 decay photons from pi0 or eta. Input in AODMCParticles.
f8006433 988
2644ead9 989 if(labels[0] < 0 || labels[0] > mcparticles->GetEntriesFast() || nlabels <= 1)
990 {
3c75ddf7 991 if(fDebug > 2)
992 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Exit : label[0] %d, n primaries %d, nlabels %d \n",
993 labels[0],mcparticles->GetEntriesFast(), nlabels);
994 return;
995 }
c5693f62 996
f8006433 997 AliAODMCParticle * meson = (AliAODMCParticle *) mcparticles->At(mesonIndex);
3c75ddf7 998 Int_t mesonPdg = meson->GetPdgCode();
2644ead9 999 if(mesonPdg != 111 && mesonPdg != 221)
1000 {
3c75ddf7 1001 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Wrong pi0/eta PDG : %d \n",mesonPdg);
1002 return;
1003 }
1004
1005 if(fDebug > 2) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - pdg %d, label %d, ndaughters %d\n", mesonPdg, mesonIndex, meson->GetNDaughters());
1006
1007
1008 //Get the daughters
2644ead9 1009 if(meson->GetNDaughters() != 2)
1010 {
3c75ddf7 1011 if(fDebug > 2)
36769d30 1012 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) - Not overalapped. Number of daughters is %d, not 2 \n",meson->GetNDaughters());
3c75ddf7 1013 return;
1014 }
2644ead9 1015
3c75ddf7 1016 Int_t iPhoton0 = meson->GetDaughter(0);
1017 Int_t iPhoton1 = meson->GetDaughter(1);
1018 //if((iPhoton0 == -1) || (iPhoton1 == -1)){
1019 // if(fDebug > 2)
1020 // printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Exit : Not overlapped. At least a daughter do not exists : d1 %d, d2 %d \n", iPhoton0, iPhoton1);
1021 // return;
1022 //}
1023 AliAODMCParticle *photon0 = (AliAODMCParticle *) mcparticles->At(iPhoton0);
1024 AliAODMCParticle *photon1 = (AliAODMCParticle *) mcparticles->At(iPhoton1);
2644ead9 1025
3c75ddf7 1026 //Check if both daughters are photons
2644ead9 1027 if(photon0->GetPdgCode() != 22 && photon1->GetPdgCode()!=22)
1028 {
3c75ddf7 1029 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Not overlapped. PDG: daughter 1 = %d, of daughter 2 = %d \n",photon0->GetPdgCode(),photon1->GetPdgCode());
1030 return;
1031 }
1032
1033 if(fDebug > 2)
1034 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Daughter labels : photon0 = %d, photon1 = %d \n",iPhoton0,iPhoton1);
1035
1036 //Check if both photons contribute to the cluster
1037 Bool_t okPhoton0 = kFALSE;
1038 Bool_t okPhoton1 = kFALSE;
1039
1040 if(fDebug > 3)
1041 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Labels loop:\n");
1042
36769d30 1043 Bool_t conversion = kFALSE;
1044
2644ead9 1045 for(Int_t i = 0; i < nlabels; i++)
1046 {
1047 if(fDebug > 3)
3c75ddf7 1048 printf("\t label %d/%d: %d, ok? %d, %d\n", i, nlabels, labels[i], okPhoton0, okPhoton1);
f8006433 1049
2644ead9 1050 if(labels[i]<0) continue;
1051
3c75ddf7 1052 //If we already found both, break the loop
1053 if(okPhoton0 && okPhoton1) break;
f8006433 1054
3c75ddf7 1055 Int_t index = labels[i];
36769d30 1056 if (iPhoton0 == index)
1057 {
3c75ddf7 1058 okPhoton0 = kTRUE;
1059 continue;
1060 }
36769d30 1061 else if (iPhoton1 == index)
1062 {
3c75ddf7 1063 okPhoton1 = kTRUE;
1064 continue;
1065 }
1066
1067 //Trace back the mother in case it was a conversion
1068
36769d30 1069 if(index >= mcparticles->GetEntriesFast())
1070 {
3c75ddf7 1071 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) Particle index %d larger than size of list %d !!\n",index,mcparticles->GetEntriesFast());
1072 continue;
1073 }
1074
1075 AliAODMCParticle * daught = (AliAODMCParticle*) mcparticles->At(index);
1076 Int_t tmpindex = daught->GetMother();
1077 if(fDebug > 3)
1078 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Conversion? : mother %d\n",tmpindex);
1079
1080 while(tmpindex>=0){
1081
1082 //MC particle of interest is the mother
1083 if(fDebug > 3)
1084 printf("\t parent index %d\n",tmpindex);
1085 daught = (AliAODMCParticle*) mcparticles->At(tmpindex);
1086 //printf("tmpindex %d\n",tmpindex);
36769d30 1087 if (iPhoton0 == tmpindex)
1088 {
1089 conversion = kTRUE;
1090 okPhoton0 = kTRUE;
3c75ddf7 1091 break;
1092 }
36769d30 1093 else if (iPhoton1 == tmpindex)
1094 {
1095 conversion = kTRUE;
1096 okPhoton1 = kTRUE;
3c75ddf7 1097 break;
36769d30 1098 }
1099
3c75ddf7 1100 tmpindex = daught->GetMother();
36769d30 1101
3c75ddf7 1102 }//While to check if pi0/eta daughter was one of these contributors to the cluster
1103
2644ead9 1104 if(i == 0 && (!okPhoton0 && !okPhoton1) && fDebug>=0 )
3c75ddf7 1105 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Something happens, first label should be from a photon decay!\n");
1106
1107 }//loop on list of labels
1108
1109 //If both photons contribute tag as the corresponding meson.
36769d30 1110 if(okPhoton0 && okPhoton1)
1111 {
3c75ddf7 1112 if(fDebug > 2)
36769d30 1113 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - %s OVERLAPPED DECAY \n",
1114 (TDatabasePDG::Instance()->GetParticle(mesonPdg))->GetName());
1115
1116 if(!CheckTagBit(tag,kMCConversion) && conversion)
1117 {
1118 if(fDebug > 2)
1119 printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) Second decay photon produced a conversion \n");
1120 SetTagBit(tag,kMCConversion) ;
1121 }
3c75ddf7 1122
1123 if(mesonPdg == 111) SetTagBit(tag,kMCPi0);
1124 else SetTagBit(tag,kMCEta);
1125 }
1126
902aa95c 1127}
7cd4e982 1128
3e57baab 1129//______________________________________________________________________________________________________
1130void AliMCAnalysisUtils::CheckLostDecayPair(const TObjArray* arrayCluster, Int_t iMom, Int_t iParent,
1131 AliStack * stack, Int_t & tag)
1132{
0d360687 1133 // Check on ESDs if the current decay photon has the second photon companion lost.
3e57baab 1134
1135 if(!arrayCluster || iMom < 0 || iParent < 0|| !stack) return;
1136
1137 TParticle * parent= stack->Particle(iParent);
1138
1139 if(parent->GetNDaughters()!=2)
1140 {
1141 SetTagBit(tag, kMCDecayPairLost);
1142 return ;
1143 }
1144
1145 Int_t pairLabel = -1;
1146 if ( iMom != parent->GetDaughter(0) ) pairLabel = parent->GetDaughter(0);
1147 else if( iMom != parent->GetDaughter(1) ) pairLabel = parent->GetDaughter(1);
1148
1149 if(pairLabel<0)
1150 {
1151 SetTagBit(tag, kMCDecayPairLost);
1152 return ;
1153 }
1154
1155 for(Int_t iclus = 0; iclus < arrayCluster->GetEntriesFast(); iclus++)
1156 {
1157 AliVCluster * cluster = (AliVCluster*) arrayCluster->At(iclus);
1158 for(UInt_t ilab = 0; ilab< cluster->GetNLabels(); ilab++)
1159 {
1160 Int_t label = cluster->GetLabels()[ilab];
1161 if(label==pairLabel)
1162 {
1163 SetTagBit(tag, kMCDecayPairInCalo);
1164 return ;
1165 }
1166 else if(label== iParent || label== iMom)
1167 {
1168 continue;
1169 }
1170 else // check the ancestry
1171 {
1172 TParticle * mother = stack->Particle(label);
1173 Int_t momPDG = TMath::Abs(mother->GetPdgCode());
1174 if(momPDG!=11 && momPDG!=22) continue;
1175
1176 //Check if "mother" of entity is converted, if not, get the first non converted mother
1177 Int_t iParentClus = mother->GetFirstMother();
1178 if(iParentClus < 0) continue;
1179
1180 TParticle * parentClus = stack->Particle(iParentClus);
1181 if(!parentClus) continue;
1182
1183 Int_t parentClusPDG = TMath::Abs(parentClus->GetPdgCode());
1184 Int_t parentClusStatus = parentClus->GetStatusCode();
1185
1186 if( parentClusPDG != 22 && parentClusPDG != 11 && parentClusStatus != 0) continue;
1187
1188 //printf("Conversion\n");
1189
1190 //Check if the mother is photon or electron with status not stable
1191 while ((parentClusPDG == 22 || parentClusPDG == 11) && parentClusStatus != 1)
1192 {
1193 //New Mother
1194 label = iParentClus;
1195 momPDG = parentClusPDG;
1196
1197 iParentClus = parentClus->GetFirstMother();
1198 if(iParentClus < 0) break;
1199
1200 parentClus = stack->Particle(iParentClus);
1201 if(!parentClus) break;
1202
1203 parentClusPDG = TMath::Abs(parentClus->GetPdgCode());
1204 parentClusStatus = parentClus->GetStatusCode() ;
1205 }//while
1206
1207 if((momPDG == 22 || parentClusPDG ==22) && (label==pairLabel || iParentClus == pairLabel))
1208 {
1209 SetTagBit(tag, kMCDecayPairInCalo);
1210 //printf("Conversion is paired\n");
1211 return ;
1212 }
1213 else continue;
1214
1215 }
1216 }
1217 } // cluster loop
1218
1219 SetTagBit(tag, kMCDecayPairLost);
1220
1221}
1222
1223//______________________________________________________________________________________________________
1224void AliMCAnalysisUtils::CheckLostDecayPair(const TObjArray * arrayCluster,Int_t iMom, Int_t iParent,
1225 const TClonesArray* mcparticles, Int_t & tag)
1226{
0d360687 1227 // Check on AODs if the current decay photon has the second photon companion lost.
3e57baab 1228
1229 if(!arrayCluster || iMom < 0 || iParent < 0|| !mcparticles) return;
1230
1231 AliAODMCParticle * parent = (AliAODMCParticle*) mcparticles->At(iParent);
1232
1233 //printf("*** Check label %d with parent %d\n",iMom, iParent);
1234
1235 if(parent->GetNDaughters()!=2)
1236 {
1237 SetTagBit(tag, kMCDecayPairLost);
1238 //printf("\t ndaugh = %d\n",parent->GetNDaughters());
1239 return ;
1240 }
1241
1242 Int_t pairLabel = -1;
1243 if ( iMom != parent->GetDaughter(0) ) pairLabel = parent->GetDaughter(0);
1244 else if( iMom != parent->GetDaughter(1) ) pairLabel = parent->GetDaughter(1);
1245
1246 if(pairLabel<0)
1247 {
1248 //printf("\t pair Label not found = %d\n",pairLabel);
1249 SetTagBit(tag, kMCDecayPairLost);
1250 return ;
1251 }
1252
1253 //printf("\t *** find pair %d\n",pairLabel);
1254
1255 for(Int_t iclus = 0; iclus < arrayCluster->GetEntriesFast(); iclus++)
1256 {
1257 AliVCluster * cluster = (AliVCluster*) arrayCluster->At(iclus);
1258 //printf("\t \t ** Cluster %d, nlabels %d\n",iclus,cluster->GetNLabels());
1259 for(UInt_t ilab = 0; ilab< cluster->GetNLabels(); ilab++)
1260 {
1261 Int_t label = cluster->GetLabels()[ilab];
1262
1263 //printf("\t \t label %d\n",label);
1264
1265 if(label==pairLabel)
1266 {
1267 //printf("\t \t Pair found\n");
1268 SetTagBit(tag, kMCDecayPairInCalo);
1269 return ;
1270 }
1271 else if(label== iParent || label== iMom)
1272 {
1273 //printf("\t \t skip\n");
1274 continue;
1275 }
1276 else // check the ancestry
1277 {
1278 AliAODMCParticle * mother = (AliAODMCParticle*) mcparticles->At(label);
1279 Int_t momPDG = TMath::Abs(mother->GetPdgCode());
1280 if(momPDG!=11 && momPDG!=22)
1281 {
1282 //printf("\t \t skip, pdg %d\n",momPDG);
1283 continue;
1284 }
1285
1286 //Check if "mother" of entity is converted, if not, get the first non converted mother
1287 Int_t iParentClus = mother->GetMother();
1288 if(iParentClus < 0) continue;
1289
1290 AliAODMCParticle * parentClus = (AliAODMCParticle*) mcparticles->At(iParentClus);
1291 if(!parentClus) continue;
1292
1293 Int_t parentClusPDG = TMath::Abs(parentClus->GetPdgCode());
1294 Int_t parentClusStatus = parentClus->GetStatus();
1295
1296 if( parentClusPDG != 22 && parentClusPDG != 11 && parentClusStatus != 0)
1297 {
1298 //printf("\t \t skip, not a conversion, parent: pdg %d, status %d\n",parentClusPDG,parentClusStatus);
1299 continue;
1300 }
1301
1302 //printf("\t \t Conversion\n");
1303
1304 //Check if the mother is photon or electron with status not stable
1305 while ((parentClusPDG == 22 || parentClusPDG == 11) && parentClusStatus != 1)
1306 {
1307 //New Mother
1308 label = iParentClus;
1309 momPDG = parentClusPDG;
1310
1311 iParentClus = parentClus->GetMother();
1312 if(iParentClus < 0) break;
1313
1314 parentClus = (AliAODMCParticle*) mcparticles->At(iParentClus);
1315 if(!parentClus) break;
1316
1317 parentClusPDG = TMath::Abs(parentClus->GetPdgCode());
1318 parentClusStatus = parentClus->GetStatus() ;
1319 }//while
1320
1321 if((momPDG == 22 || parentClusPDG ==22) && (label==pairLabel || iParentClus == pairLabel))
1322 {
1323 SetTagBit(tag, kMCDecayPairInCalo);
1324 //printf("\t \t Conversion is paired: mom %d, parent %d\n",label,iParentClus);
1325 return ;
1326 }
1327 else
1328 {
1329 //printf("\t \t Skip, finally label %d, pdg %d, parent label %d, pdg %d, status %d\n",label,momPDG,iParentClus,parentClusPDG,parentClusStatus);
1330 continue;
1331 }
1332
1333 }
1334 }
1335 } // cluster loop
1336
1337
1338 SetTagBit(tag, kMCDecayPairLost);
1339
1340
1341}
1342
8a2dbbff 1343//_____________________________________________________________________
c5693f62 1344TList * AliMCAnalysisUtils::GetJets(const AliCaloTrackReader * reader)
840124fe 1345{
0d360687 1346 // Return list of jets (TParticles) and index of most likely parton that originated it.
1347
7cd4e982 1348 AliStack * stack = reader->GetStack();
1349 Int_t iEvent = reader->GetEventNumber();
1350 AliGenEventHeader * geh = reader->GetGenEventHeader();
1351 if(fCurrentEvent!=iEvent){
1352 fCurrentEvent = iEvent;
1353 fJetsList = new TList;
1354 Int_t nTriggerJets = 0;
1355 Float_t tmpjet[]={0,0,0,0};
1356
1357 //printf("Event %d %d\n",fCurrentEvent,iEvent);
1358 //Get outgoing partons
1359 if(stack->GetNtrack() < 8) return fJetsList;
1360 TParticle * parton1 = stack->Particle(6);
1361 TParticle * parton2 = stack->Particle(7);
1362 if(fDebug > 2){
1363 printf("AliMCAnalysisUtils::GetJets() - parton 6 : %s, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
f8006433 1364 parton1->GetName(),parton1->Pt(),parton1->Energy(),parton1->Phi()*TMath::RadToDeg(),parton1->Eta());
7cd4e982 1365 printf("AliMCAnalysisUtils::GetJets() - parton 7 : %s, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
f8006433 1366 parton2->GetName(),parton2->Pt(),parton2->Energy(),parton2->Phi()*TMath::RadToDeg(),parton2->Eta());
7cd4e982 1367 }
f8006433 1368 // //Trace the jet from the mother parton
1369 // Float_t pt = 0;
1370 // Float_t pt1 = 0;
1371 // Float_t pt2 = 0;
1372 // Float_t e = 0;
1373 // Float_t e1 = 0;
1374 // Float_t e2 = 0;
1375 // TParticle * tmptmp = new TParticle;
1376 // for(Int_t i = 0; i< stack->GetNprimary(); i++){
1377 // tmptmp = stack->Particle(i);
7cd4e982 1378
f8006433 1379 // if(tmptmp->GetStatusCode() == 1){
1380 // pt = tmptmp->Pt();
1381 // e = tmptmp->Energy();
1382 // Int_t imom = tmptmp->GetFirstMother();
1383 // Int_t imom1 = 0;
1384 // //printf("1st imom %d\n",imom);
1385 // while(imom > 5){
1386 // imom1=imom;
1387 // tmptmp = stack->Particle(imom);
1388 // imom = tmptmp->GetFirstMother();
1389 // //printf("imom %d \n",imom);
1390 // }
1391 // //printf("Last imom %d %d\n",imom1, imom);
1392 // if(imom1 == 6) {
1393 // pt1+=pt;
1394 // e1+=e;
1395 // }
1396 // else if (imom1 == 7){
1397 // pt2+=pt;
1398 // e2+=e; }
1399 // }// status 1
1400
1401 // }// for
7cd4e982 1402
f8006433 1403 // printf("JET 1, pt %2.2f, e %2.2f; JET 2, pt %2.2f, e %2.2f \n",pt1,e1,pt2,e2);
7cd4e982 1404
1405 //Get the jet, different way for different generator
1406 //PYTHIA
0d360687 1407 if(fMCGenerator == kPythia)
4d52b003 1408 {
f8006433 1409 TParticle * jet = 0x0;
7cd4e982 1410 AliGenPythiaEventHeader* pygeh= (AliGenPythiaEventHeader*) geh;
1411 nTriggerJets = pygeh->NTriggerJets();
1412 if(fDebug > 1)
f8006433 1413 printf("AliMCAnalysisUtils::GetJets() - PythiaEventHeader: Njets: %d\n",nTriggerJets);
1414
265cab75 1415 for(Int_t i = 0; i< nTriggerJets; i++)
1416 {
f8006433 1417 pygeh->TriggerJet(i, tmpjet);
1418 jet = new TParticle(94, 21, -1, -1, -1, -1, tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3], 0,0,0,0);
1419 //Assign an outgoing parton as mother
1420 Float_t phidiff1 = TMath::Abs(jet->Phi()-parton1->Phi());
1421 Float_t phidiff2 = TMath::Abs(jet->Phi()-parton2->Phi());
1422 if(phidiff1 > phidiff2) jet->SetFirstMother(7);
1423 else jet->SetFirstMother(6);
1424 //jet->Print();
1425 if(fDebug > 1)
1426 printf("AliMCAnalysisUtils::GetJets() - PYTHIA Jet %d: mother %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
1427 i, jet->GetFirstMother(),jet->Pt(),jet->Energy(),jet->Phi()*TMath::RadToDeg(),jet->Eta());
1428 fJetsList->Add(jet);
7cd4e982 1429 }
1430 }//Pythia triggered jets
1431 //HERWIG
0d360687 1432 else if (fMCGenerator == kHerwig)
265cab75 1433 {
7cd4e982 1434 Int_t pdg = -1;
1435 //Check parton 1
1436 TParticle * tmp = parton1;
265cab75 1437 if(parton1->GetPdgCode()!=22)
1438 {
f8006433 1439 while(pdg != 94){
1440 if(tmp->GetFirstDaughter()==-1) return fJetsList;
1441 tmp = stack->Particle(tmp->GetFirstDaughter());
1442 pdg = tmp->GetPdgCode();
1443 }//while
1444
1445 //Add found jet to list
1446 TParticle *jet1 = new TParticle(*tmp);
1447 jet1->SetFirstMother(6);
1448 fJetsList->Add(jet1);
1449 //printf("jet 1: first daughter %d, last daughter %d\n", tmp->GetFirstDaughter(), tmp->GetLastDaughter());
1450 //tmp = stack->Particle(tmp->GetFirstDaughter());
1451 //tmp->Print();
1452 //jet1->Print();
1453 if(fDebug > 1)
1454 printf("AliMCAnalysisUtils::GetJets() - HERWIG Jet 1: mother %d, status %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
1455 jet1->GetFirstMother(),jet1->GetStatusCode(),jet1->Pt(),jet1->Energy(),jet1->Phi()*TMath::RadToDeg(),jet1->Eta());
7cd4e982 1456 }//not photon
1457
1458 //Check parton 2
1459 pdg = -1;
1460 tmp = parton2;
265cab75 1461 if(parton2->GetPdgCode()!=22)
1462 {
1463 while(pdg != 94)
1464 {
f8006433 1465 if(tmp->GetFirstDaughter()==-1) return fJetsList;
f8006433 1466 tmp = stack->Particle(tmp->GetFirstDaughter());
1467 pdg = tmp->GetPdgCode();
1468 }//while
265cab75 1469
f8006433 1470 //Add found jet to list
1471 TParticle *jet2 = new TParticle(*tmp);
1472 jet2->SetFirstMother(7);
1473 fJetsList->Add(jet2);
1474 //jet2->Print();
1475 if(fDebug > 1)
1476 printf("AliMCAnalysisUtils::GetJets() - HERWIG Jet 2: mother %d, status %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
1477 jet2->GetFirstMother(),jet2->GetStatusCode(),jet2->Pt(),jet2->Energy(),jet2->Phi()*TMath::RadToDeg(),jet2->Eta());
1478 //Int_t first = tmp->GetFirstDaughter();
1479 //Int_t last = tmp->GetLastDaughter();
1480 //printf("jet 2: first daughter %d, last daughter %d, pdg %d\n",first, last, tmp->GetPdgCode());
7cd4e982 1481 // for(Int_t d = first ; d < last+1; d++){
f8006433 1482 // tmp = stack->Particle(d);
1483 // if(i == tmp->GetFirstMother())
1484 // printf("Daughter n %d, Mother %d, name %s, status %d, pT %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
1485 // d,tmp->GetFirstMother(), tmp->GetName(), tmp->GetStatusCode(),tmp->Pt(),tmp->Energy(),tmp->Phi()*TMath::RadToDeg(),tmp->Eta());
1486 // }
1487 //tmp->Print();
7cd4e982 1488 }//not photon
1489 }//Herwig generated jets
1490 }
1491
1492 return fJetsList;
1493}
1494
51a0ace5 1495
8a2dbbff 1496//__________________________________________________________________________________________________________
1497TLorentzVector AliMCAnalysisUtils::GetDaughter(Int_t idaugh, Int_t label,
64be16b4 1498 const AliCaloTrackReader* reader,
1499 Int_t & pdg, Int_t & status, Bool_t & ok, Int_t & daughlabel)
1500{
0d360687 1501 // Return the kinematics of the particle that generated the signal, its pdg and its status and its label mother.
1502 fDaughMom.SetPxPyPzE(0,0,0,0);
64be16b4 1503
1504 if(reader->ReadStack())
1505 {
1506 if(!reader->GetStack())
1507 {
1508 if (fDebug >=0)
1509 printf("AliMCAnalysisUtils::GetDaughter() - Stack is not available, check analysis settings in configuration file, STOP!!\n");
1510
1511 ok=kFALSE;
0d360687 1512 return fDaughMom;
64be16b4 1513 }
4d6b53b6 1514
1515 Int_t nprimaries = reader->GetStack()->GetNtrack();
1516 if(label >= 0 && label < nprimaries)
64be16b4 1517 {
1518 TParticle * momP = reader->GetStack()->Particle(label);
1519 daughlabel = momP->GetDaughter(idaugh);
1520
4d6b53b6 1521 if(daughlabel < 0 || daughlabel >= nprimaries)
1522 {
1523 ok = kFALSE;
0d360687 1524 return fDaughMom;
4d6b53b6 1525 }
1526
64be16b4 1527 TParticle * daughP = reader->GetStack()->Particle(daughlabel);
0d360687 1528 daughP->Momentum(fDaughMom);
64be16b4 1529 pdg = daughP->GetPdgCode();
1530 status = daughP->GetStatusCode();
1531 }
1532 else
1533 {
1534 ok = kFALSE;
0d360687 1535 return fDaughMom;
64be16b4 1536 }
1537 }
1538 else if(reader->ReadAODMCParticles())
1539 {
1540 TClonesArray* mcparticles = reader->GetAODMCParticles();
1541 if(!mcparticles)
1542 {
1543 if(fDebug >= 0)
1544 printf("AliMCAnalysisUtils::GetDaughter() - AODMCParticles is not available, check analysis settings in configuration file!!\n");
1545
1546 ok=kFALSE;
0d360687 1547 return fDaughMom;
64be16b4 1548 }
1549
1550 Int_t nprimaries = mcparticles->GetEntriesFast();
1551 if(label >= 0 && label < nprimaries)
1552 {
1553 AliAODMCParticle * momP = (AliAODMCParticle *) mcparticles->At(label);
1554 daughlabel = momP->GetDaughter(idaugh);
1555
4d6b53b6 1556 if(daughlabel < 0 || daughlabel >= nprimaries)
1557 {
1558 ok = kFALSE;
0d360687 1559 return fDaughMom;
4d6b53b6 1560 }
1561
64be16b4 1562 AliAODMCParticle * daughP = (AliAODMCParticle *) mcparticles->At(daughlabel);
0d360687 1563 fDaughMom.SetPxPyPzE(daughP->Px(),daughP->Py(),daughP->Pz(),daughP->E());
64be16b4 1564 pdg = daughP->GetPdgCode();
1565 status = daughP->GetStatus();
1566 }
1567 else
1568 {
1569 ok = kFALSE;
0d360687 1570 return fDaughMom;
64be16b4 1571 }
1572 }
1573
1574 ok = kTRUE;
1575
0d360687 1576 return fDaughMom;
64be16b4 1577}
1578
8a2dbbff 1579//______________________________________________________________________________________________________
1580TLorentzVector AliMCAnalysisUtils::GetMother(Int_t label, const AliCaloTrackReader* reader, Bool_t & ok)
9036b2a6 1581{
0d360687 1582 // Return the kinematics of the particle that generated the signal.
9036b2a6 1583
64be16b4 1584 Int_t pdg = -1; Int_t status = -1; Int_t momlabel = -1;
1585 return GetMother(label,reader,pdg,status, ok,momlabel);
1586}
1587
8a2dbbff 1588//_________________________________________________________________________________________
1589TLorentzVector AliMCAnalysisUtils::GetMother(Int_t label, const AliCaloTrackReader* reader,
64be16b4 1590 Int_t & pdg, Int_t & status, Bool_t & ok)
1591{
0d360687 1592 // Return the kinematics of the particle that generated the signal.
64be16b4 1593
1594 Int_t momlabel = -1;
1595 return GetMother(label,reader,pdg,status, ok,momlabel);
51a0ace5 1596}
1597
8a2dbbff 1598//______________________________________________________________________________________________________
1599TLorentzVector AliMCAnalysisUtils::GetMother(Int_t label, const AliCaloTrackReader* reader,
64be16b4 1600 Int_t & pdg, Int_t & status, Bool_t & ok, Int_t & momlabel)
51a0ace5 1601{
0d360687 1602 // Return the kinematics of the particle that generated the signal, its pdg and its status and its label mother.
51a0ace5 1603
0d360687 1604 fMotherMom.SetPxPyPzE(0,0,0,0);
9036b2a6 1605
1606 if(reader->ReadStack())
1607 {
51a0ace5 1608 if(!reader->GetStack())
1609 {
9036b2a6 1610 if (fDebug >=0)
1611 printf("AliMCAnalysisUtils::GetMother() - Stack is not available, check analysis settings in configuration file, STOP!!\n");
51a0ace5 1612
1613 ok=kFALSE;
0d360687 1614 return fMotherMom;
9036b2a6 1615 }
1616 if(label >= 0 && label < reader->GetStack()->GetNtrack())
1617 {
1618 TParticle * momP = reader->GetStack()->Particle(label);
0d360687 1619 momP->Momentum(fMotherMom);
1620 pdg = momP->GetPdgCode();
1621 status = momP->GetStatusCode();
64be16b4 1622 momlabel = momP->GetFirstMother();
51a0ace5 1623 }
1624 else
1625 {
1626 ok = kFALSE;
0d360687 1627 return fMotherMom;
9036b2a6 1628 }
1629 }
1630 else if(reader->ReadAODMCParticles())
1631 {
2644ead9 1632 TClonesArray* mcparticles = reader->GetAODMCParticles();
9036b2a6 1633 if(!mcparticles)
1634 {
1635 if(fDebug >= 0)
1636 printf("AliMCAnalysisUtils::GetMother() - AODMCParticles is not available, check analysis settings in configuration file!!\n");
51a0ace5 1637
1638 ok=kFALSE;
0d360687 1639 return fMotherMom;
9036b2a6 1640 }
1641
1642 Int_t nprimaries = mcparticles->GetEntriesFast();
1643 if(label >= 0 && label < nprimaries)
1644 {
1645 AliAODMCParticle * momP = (AliAODMCParticle *) mcparticles->At(label);
0d360687 1646 fMotherMom.SetPxPyPzE(momP->Px(),momP->Py(),momP->Pz(),momP->E());
1647 pdg = momP->GetPdgCode();
1648 status = momP->GetStatus();
64be16b4 1649 momlabel = momP->GetMother();
1650 }
51a0ace5 1651 else
1652 {
1653 ok = kFALSE;
0d360687 1654 return fMotherMom;
9036b2a6 1655 }
1656 }
1657
51a0ace5 1658 ok = kTRUE;
1659
0d360687 1660 return fMotherMom;
9036b2a6 1661}
1662
8a2dbbff 1663//___________________________________________________________________________________
1664TLorentzVector AliMCAnalysisUtils::GetMotherWithPDG(Int_t label, Int_t pdg,
1665 const AliCaloTrackReader* reader,
36769d30 1666 Bool_t & ok, Int_t & momlabel)
9036b2a6 1667{
0d360687 1668 // Return the kinematics of the particle that generated the signal.
9036b2a6 1669
0d360687 1670 fGMotherMom.SetPxPyPzE(0,0,0,0);
9036b2a6 1671
1672 if(reader->ReadStack())
1673 {
1674 if(!reader->GetStack())
1675 {
1676 if (fDebug >=0)
1677 printf("AliMCAnalysisUtils::GetMotherWithPDG() - Stack is not available, check analysis settings in configuration file, STOP!!\n");
51a0ace5 1678
1679 ok = kFALSE;
0d360687 1680 return fGMotherMom;
9036b2a6 1681 }
36769d30 1682
9036b2a6 1683 if(label >= 0 && label < reader->GetStack()->GetNtrack())
1684 {
1685 TParticle * momP = reader->GetStack()->Particle(label);
1686
1687 Int_t grandmomLabel = momP->GetFirstMother();
1688 Int_t grandmomPDG = -1;
1689 TParticle * grandmomP = 0x0;
1690 while (grandmomLabel >=0 )
1691 {
1692 grandmomP = reader->GetStack()->Particle(grandmomLabel);
1693 grandmomPDG = grandmomP->GetPdgCode();
1694 if(grandmomPDG==pdg)
1695 {
36769d30 1696 momlabel = grandmomLabel;
0d360687 1697 fGMotherMom.SetPxPyPzE(grandmomP->Px(),grandmomP->Py(),grandmomP->Pz(),grandmomP->Energy());
9036b2a6 1698 break;
1699 }
1700
1701 grandmomLabel = grandmomP->GetFirstMother();
1702
1703 }
1704
1705 if(grandmomPDG!=pdg) printf("AliMCAnalysisUtils::GetMotherWithPDG(ESD) - mother with PDG %d, not found! \n",pdg);
1706 }
1707 }
1708 else if(reader->ReadAODMCParticles())
1709 {
2644ead9 1710 TClonesArray* mcparticles = reader->GetAODMCParticles();
9036b2a6 1711 if(!mcparticles)
1712 {
1713 if(fDebug >= 0)
1714 printf("AliMCAnalysisUtils::GetMotherWithPDG() - AODMCParticles is not available, check analysis settings in configuration file!!\n");
51a0ace5 1715
1716 ok=kFALSE;
0d360687 1717 return fGMotherMom;
9036b2a6 1718 }
1719
1720 Int_t nprimaries = mcparticles->GetEntriesFast();
1721 if(label >= 0 && label < nprimaries)
1722 {
1723 AliAODMCParticle * momP = (AliAODMCParticle *) mcparticles->At(label);
1724
1725 Int_t grandmomLabel = momP->GetMother();
1726 Int_t grandmomPDG = -1;
1727 AliAODMCParticle * grandmomP = 0x0;
1728 while (grandmomLabel >=0 )
1729 {
1730 grandmomP = (AliAODMCParticle *) mcparticles->At(grandmomLabel);
1731 grandmomPDG = grandmomP->GetPdgCode();
1732 if(grandmomPDG==pdg)
1733 {
1734 //printf("AliMCAnalysisUtils::GetMotherWithPDG(AOD) - mother with PDG %d FOUND! \n",pdg);
36769d30 1735 momlabel = grandmomLabel;
0d360687 1736 fGMotherMom.SetPxPyPzE(grandmomP->Px(),grandmomP->Py(),grandmomP->Pz(),grandmomP->E());
9036b2a6 1737 break;
1738 }
1739
1740 grandmomLabel = grandmomP->GetMother();
1741
1742 }
1743
1744 if(grandmomPDG!=pdg) printf("AliMCAnalysisUtils::GetMotherWithPDG(AOD) - mother with PDG %d, NOT found! \n",pdg);
1745
1746 }
1747 }
1748
51a0ace5 1749 ok = kTRUE;
9036b2a6 1750
0d360687 1751 return fGMotherMom;
9036b2a6 1752}
7cd4e982 1753
8a2dbbff 1754//______________________________________________________________________________________________
1755TLorentzVector AliMCAnalysisUtils::GetGrandMother(Int_t label, const AliCaloTrackReader* reader,
914b9fe7 1756 Int_t & pdg, Int_t & status, Bool_t & ok,
1757 Int_t & grandMomLabel, Int_t & greatMomLabel)
1758{
0d360687 1759 // Return the kinematics of the particle that generated the signal.
914b9fe7 1760
0d360687 1761 fGMotherMom.SetPxPyPzE(0,0,0,0);
914b9fe7 1762
1763 if(reader->ReadStack())
1764 {
1765 if(!reader->GetStack())
1766 {
1767 if (fDebug >=0)
64be16b4 1768 printf("AliMCAnalysisUtils::GetGrandMother() - Stack is not available, check analysis settings in configuration file, STOP!!\n");
914b9fe7 1769
1770 ok = kFALSE;
0d360687 1771 return fGMotherMom;
914b9fe7 1772 }
1773
1774 if(label >= 0 && label < reader->GetStack()->GetNtrack())
1775 {
1776 TParticle * momP = reader->GetStack()->Particle(label);
1777
1778 grandMomLabel = momP->GetFirstMother();
1779
1780 TParticle * grandmomP = 0x0;
1781
1782 if (grandMomLabel >=0 )
1783 {
1784 grandmomP = reader->GetStack()->Particle(grandMomLabel);
1785 pdg = grandmomP->GetPdgCode();
1786 status = grandmomP->GetStatusCode();
1787
0d360687 1788 fGMotherMom.SetPxPyPzE(grandmomP->Px(),grandmomP->Py(),grandmomP->Pz(),grandmomP->Energy());
914b9fe7 1789 greatMomLabel = grandmomP->GetFirstMother();
1790
1791 }
1792 }
1793 }
1794 else if(reader->ReadAODMCParticles())
1795 {
1796 TClonesArray* mcparticles = reader->GetAODMCParticles();
1797 if(!mcparticles)
1798 {
1799 if(fDebug >= 0)
64be16b4 1800 printf("AliMCAnalysisUtils::GetGrandMother() - AODMCParticles is not available, check analysis settings in configuration file!!\n");
914b9fe7 1801
1802 ok=kFALSE;
0d360687 1803 return fGMotherMom;
914b9fe7 1804 }
1805
1806 Int_t nprimaries = mcparticles->GetEntriesFast();
1807 if(label >= 0 && label < nprimaries)
1808 {
1809 AliAODMCParticle * momP = (AliAODMCParticle *) mcparticles->At(label);
1810
1811 grandMomLabel = momP->GetMother();
1812
1813 AliAODMCParticle * grandmomP = 0x0;
1814
1815 if(grandMomLabel >=0 )
1816 {
1817 grandmomP = (AliAODMCParticle *) mcparticles->At(grandMomLabel);
1818 pdg = grandmomP->GetPdgCode();
1819 status = grandmomP->GetStatus();
1820
0d360687 1821 fGMotherMom.SetPxPyPzE(grandmomP->Px(),grandmomP->Py(),grandmomP->Pz(),grandmomP->E());
914b9fe7 1822 greatMomLabel = grandmomP->GetMother();
1823
1824 }
1825 }
1826 }
1827
1828 ok = kTRUE;
1829
0d360687 1830 return fGMotherMom;
914b9fe7 1831}
1832
8a2dbbff 1833//_______________________________________________________________________________________________________________
1834void AliMCAnalysisUtils::GetMCDecayAsymmetryAngleForPDG(Int_t label, Int_t pdg, const AliCaloTrackReader* reader,
4efb6957 1835 Float_t & asym, Float_t & angle, Bool_t & ok)
8e81c2cf 1836{
0d360687 1837 // In case of an eta or pi0 decay into 2 photons, get the asymmetry in the energy of the photons.
8e81c2cf 1838
1839 if(reader->ReadStack())
1840 {
1841 if(!reader->GetStack())
1842 {
1843 if (fDebug >=0)
64be16b4 1844 printf("AliMCAnalysisUtils::GetMCDecayAsymmetryForPDG() - Stack is not available, check analysis settings in configuration file, STOP!!\n");
8e81c2cf 1845
1846 ok = kFALSE;
8e81c2cf 1847 }
1848 if(label >= 0 && label < reader->GetStack()->GetNtrack())
1849 {
1850 TParticle * momP = reader->GetStack()->Particle(label);
1851
1852 Int_t grandmomLabel = momP->GetFirstMother();
1853 Int_t grandmomPDG = -1;
1854 TParticle * grandmomP = 0x0;
1855 while (grandmomLabel >=0 )
1856 {
1857 grandmomP = reader->GetStack()->Particle(grandmomLabel);
1858 grandmomPDG = grandmomP->GetPdgCode();
1859
1860 if(grandmomPDG==pdg) break;
1861
1862 grandmomLabel = grandmomP->GetFirstMother();
1863
1864 }
1865
1866 if(grandmomPDG==pdg && grandmomP->GetNDaughters()==2)
1867 {
1868 TParticle * d1 = reader->GetStack()->Particle(grandmomP->GetDaughter(0));
1869 TParticle * d2 = reader->GetStack()->Particle(grandmomP->GetDaughter(1));
1870 if(d1->GetPdgCode() == 22 && d1->GetPdgCode() == 22)
1871 {
1872 asym = (d1->Energy()-d2->Energy())/grandmomP->Energy();
0d360687 1873 d1->Momentum(fDaughMom );
1874 d2->Momentum(fDaughMom2);
1875 angle = fDaughMom.Angle(fDaughMom2.Vect());
8e81c2cf 1876 }
1877 }
1878 else
1879 {
1880 ok=kFALSE;
64be16b4 1881 printf("AliMCAnalysisUtils::GetMCDecayAsymmetryForPDG(ESD) - mother with PDG %d, not found! \n",pdg);
8e81c2cf 1882 }
1883
1884 } // good label
1885 }
1886 else if(reader->ReadAODMCParticles())
1887 {
2644ead9 1888 TClonesArray* mcparticles = reader->GetAODMCParticles();
8e81c2cf 1889 if(!mcparticles)
1890 {
1891 if(fDebug >= 0)
64be16b4 1892 printf("AliMCAnalysisUtils::GetMCDecayAsymmetryForPDG() - AODMCParticles is not available, check analysis settings in configuration file!!\n");
8e81c2cf 1893
1894 ok=kFALSE;
c31af35c 1895 return;
8e81c2cf 1896 }
1897
1898 Int_t nprimaries = mcparticles->GetEntriesFast();
1899 if(label >= 0 && label < nprimaries)
1900 {
1901 AliAODMCParticle * momP = (AliAODMCParticle *) mcparticles->At(label);
1902
1903 Int_t grandmomLabel = momP->GetMother();
1904 Int_t grandmomPDG = -1;
1905 AliAODMCParticle * grandmomP = 0x0;
1906 while (grandmomLabel >=0 )
1907 {
1908 grandmomP = (AliAODMCParticle *) mcparticles->At(grandmomLabel);
1909 grandmomPDG = grandmomP->GetPdgCode();
1910
1911 if(grandmomPDG==pdg) break;
1912
1913 grandmomLabel = grandmomP->GetMother();
1914
1915 }
1916
1917 if(grandmomPDG==pdg && grandmomP->GetNDaughters()==2)
1918 {
1919 AliAODMCParticle * d1 = (AliAODMCParticle *) mcparticles->At(grandmomP->GetDaughter(0));
1920 AliAODMCParticle * d2 = (AliAODMCParticle *) mcparticles->At(grandmomP->GetDaughter(1));
1921 if(d1->GetPdgCode() == 22 && d1->GetPdgCode() == 22)
1922 {
1923 asym = (d1->E()-d2->E())/grandmomP->E();
0d360687 1924 fDaughMom .SetPxPyPzE(d1->Px(),d1->Py(),d1->Pz(),d1->E());
1925 fDaughMom2.SetPxPyPzE(d2->Px(),d2->Py(),d2->Pz(),d2->E());
1926 angle = fDaughMom.Angle(fDaughMom2.Vect());
8e81c2cf 1927 }
1928 }
1929 else
1930 {
1931 ok=kFALSE;
64be16b4 1932 printf("AliMCAnalysisUtils::GetMCDecayAsymmetryForPDG(AOD) - mother with PDG %d, not found! \n",pdg);
8e81c2cf 1933 }
1934
1935 } // good label
1936 }
1937
1938 ok = kTRUE;
1939
8e81c2cf 1940}
1941
8a2dbbff 1942//_________________________________________________________________________________________________
1943Int_t AliMCAnalysisUtils::GetNDaughters(Int_t label, const AliCaloTrackReader* reader, Bool_t & ok)
64be16b4 1944{
0d360687 1945 // Return the the number of daughters of a given MC particle.
64be16b4 1946
1947 if(reader->ReadStack())
1948 {
1949 if(!reader->GetStack())
1950 {
1951 if (fDebug >=0)
1952 printf("AliMCAnalysisUtils::GetNDaughters() - Stack is not available, check analysis settings in configuration file, STOP!!\n");
1953
1954 ok=kFALSE;
1955 return -1;
1956 }
1957 if(label >= 0 && label < reader->GetStack()->GetNtrack())
1958 {
1959 TParticle * momP = reader->GetStack()->Particle(label);
1960 ok=kTRUE;
1961 return momP->GetNDaughters();
1962 }
1963 else
1964 {
1965 ok = kFALSE;
1966 return -1;
1967 }
1968 }
1969 else if(reader->ReadAODMCParticles())
1970 {
1971 TClonesArray* mcparticles = reader->GetAODMCParticles();
1972 if(!mcparticles)
1973 {
1974 if(fDebug >= 0)
1975 printf("AliMCAnalysisUtils::GetNDaughters() - AODMCParticles is not available, check analysis settings in configuration file!!\n");
1976
1977 ok=kFALSE;
1978 return -1;
1979 }
1980
1981 Int_t nprimaries = mcparticles->GetEntriesFast();
1982 if(label >= 0 && label < nprimaries)
1983 {
1984 AliAODMCParticle * momP = (AliAODMCParticle *) mcparticles->At(label);
1985 ok = kTRUE;
1986 return momP->GetNDaughters();
1987 }
1988 else
1989 {
1990 ok = kFALSE;
1991 return -1;
1992 }
1993 }
1994
1995 ok = kFALSE;
1996
1997 return -1;
1998}
1999
8a2dbbff 2000//_________________________________________________________________________________
2001Int_t AliMCAnalysisUtils::GetNOverlaps(const Int_t * label, UInt_t nlabels,
2002 Int_t mctag, Int_t mesonLabel,
4914e781 2003 AliCaloTrackReader * reader, Int_t *overpdg)
2004{
2005 // Compare the primary depositing more energy with the rest,
0d360687 2006 // if no photon/electron (conversion) or neutral meson as comon ancestor, consider it as other particle contributing.
2007 // Give as input the meson label in case it was a pi0 or eta merged cluster.
2008 // Init overpdg with nlabels.
4914e781 2009
2010 Int_t ancPDG = 0, ancStatus = -1;
0d360687 2011 TVector3 prodVertex;
4914e781 2012 Int_t ancLabel = 0;
2013 Int_t noverlaps = 0;
2014 Bool_t ok = kFALSE;
2015
2016 for (UInt_t ilab = 1; ilab < nlabels; ilab++ )
2017 {
0d360687 2018 ancLabel = CheckCommonAncestor(label[0],label[ilab],reader,ancPDG,ancStatus,fMotherMom,prodVertex);
4914e781 2019
2020 //printf("Overlaps, i %d: Main Label %d, second label %d, ancestor: Label %d, pdg %d - tag %d \n",
2021 // ilab,label[0],label[ilab],ancLabel,ancPDG, mctag);
2022
2023 Bool_t overlap = kFALSE;
2024
2025 if ( ancLabel < 0 )
2026 {
2027 overlap = kTRUE;
2028 //printf("\t \t \t No Label = %d\n",ancLabel);
2029 }
2030 else if( ( ancPDG==111 || ancPDG==221 ) && ( CheckTagBit(mctag,kMCPi0) || CheckTagBit(mctag,kMCEta)) && mesonLabel != ancLabel)
2031 {
2032 //printf("\t \t meson Label %d, ancestor Label %d\n",mesonLabel,ancLabel);
2033 overlap = kTRUE;
2034 }
2035 else if( ancPDG!=22 && TMath::Abs(ancPDG)!=11 && ancPDG != 111 && ancPDG != 221 )
2036 {
2037 //printf("\t \t \t Non EM PDG = %d\n",ancPDG);
2038 overlap = kTRUE ;
2039 }
2040
2041 if( !overlap ) continue ;
2042
2043 // We have at least one overlap
2044
2045 //printf("Overlap!!!!!!!!!!!!!!\n");
2046
2047 noverlaps++;
2048
2049 // What is the origin of the overlap?
2050 Bool_t mOK = 0, gOK = 0;
2051 Int_t mpdg = -999999, gpdg = -1;
2052 Int_t mstatus = -1, gstatus = -1;
2053 Int_t gLabel = -1, ggLabel = -1;
0d360687 2054
2055 GetMother (label[ilab],reader,mpdg,mstatus,mOK);
2056 fGMotherMom =
2057 GetGrandMother(label[ilab],reader,gpdg,gstatus,gOK, gLabel,ggLabel);
4914e781 2058
2059 //printf("\t Overlap!, mother pdg %d; grand mother pdg %d",mpdg,gpdg);
2060
2061 if( ( mpdg == 22 || TMath::Abs(mpdg==11) ) &&
2062 ( gpdg == 22 || TMath::Abs(gpdg==11) ) &&
2063 gLabel >=0 )
2064 {
2065 Int_t labeltmp = gLabel;
2066 while( ( gpdg == 22 || TMath::Abs(gpdg==11) ) && gLabel >=0 )
2067 {
2068 mpdg=gpdg;
0d360687 2069 fGMotherMom = GetGrandMother(labeltmp,reader,gpdg,gstatus,ok, gLabel,ggLabel);
4914e781 2070 labeltmp=gLabel;
2071 }
2072 }
2073 overpdg[noverlaps-1] = mpdg;
2074 }
2075
2076 return noverlaps ;
2077
2078}
2079
c5693f62 2080//________________________________________________________
7cd4e982 2081void AliMCAnalysisUtils::Print(const Option_t * opt) const
2082{
0d360687 2083 // Print some relevant parameters set for the analysis.
f8006433 2084
2085 if(! opt)
2086 return;
2087
2088 printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
2089
2090 printf("Debug level = %d\n",fDebug);
0d360687 2091 printf("MC Generator = %s\n",fMCGeneratorString.Data());
f8006433 2092 printf(" \n");
2093
7cd4e982 2094}
2095
8a2dbbff 2096//__________________________________________________
2097void AliMCAnalysisUtils::PrintMCTag(Int_t tag) const
bbb79e66 2098{
0d360687 2099 // Print the assigned origins to this particle.
bbb79e66 2100
4250012f 2101 printf("AliMCAnalysisUtils::PrintMCTag() - tag %d \n photon %d, conv %d, prompt %d, frag %d, isr %d, \n pi0 decay %d, eta decay %d, other decay %d pi0 %d, eta %d \n electron %d, muon %d,pion %d, proton %d, neutron %d, \n kaon %d, a-proton %d, a-neutron %d, unk %d, bad %d\n",
bbb79e66 2102 tag,
2103 CheckTagBit(tag,kMCPhoton),
bfc290d0 2104 CheckTagBit(tag,kMCConversion),
bbb79e66 2105 CheckTagBit(tag,kMCPrompt),
2106 CheckTagBit(tag,kMCFragmentation),
2107 CheckTagBit(tag,kMCISR),
2108 CheckTagBit(tag,kMCPi0Decay),
2109 CheckTagBit(tag,kMCEtaDecay),
2110 CheckTagBit(tag,kMCOtherDecay),
bbb79e66 2111 CheckTagBit(tag,kMCPi0),
bfc290d0 2112 CheckTagBit(tag,kMCEta),
bbb79e66 2113 CheckTagBit(tag,kMCElectron),
bbb79e66 2114 CheckTagBit(tag,kMCMuon),
2115 CheckTagBit(tag,kMCPion),
2116 CheckTagBit(tag,kMCProton),
2117 CheckTagBit(tag,kMCAntiNeutron),
2118 CheckTagBit(tag,kMCKaon),
2119 CheckTagBit(tag,kMCAntiProton),
2120 CheckTagBit(tag,kMCAntiNeutron),
bfc290d0 2121 CheckTagBit(tag,kMCUnknown),
bbb79e66 2122 CheckTagBit(tag,kMCBadLabel)
2123 );
bbb79e66 2124}
2125
0d360687 2126//__________________________________________________
2127void AliMCAnalysisUtils::SetMCGenerator(Int_t mcgen)
2128{
2129 // Set the generator type.
2130
2131 fMCGenerator = mcgen ;
2132 if (mcgen == kPythia) fMCGeneratorString = "PYTHIA";
2133 else if(mcgen == kHerwig) fMCGeneratorString = "HERWIG";
2134 else if(mcgen == kHijing) fMCGeneratorString = "HIJING";
2135 else
2136 {
2137 fMCGeneratorString = "";
2138 fMCGenerator = kBoxLike ;
2139 }
2140
2141}
2142
2143//__________________________________________________
2144void AliMCAnalysisUtils::SetMCGenerator(TString mcgen)
2145{
2146 // Set the generator type.
2147
2148 fMCGeneratorString = mcgen ;
2149
2150 if (mcgen == "PYTHIA") fMCGenerator = kPythia;
2151 else if(mcgen == "HERWIG") fMCGenerator = kHerwig;
2152 else if(mcgen == "HIJING") fMCGenerator = kHijing;
2153 else
2154 {
2155 fMCGenerator = kBoxLike;
2156 fMCGeneratorString = "" ;
2157 }
2158}
2159
2160
7cd4e982 2161