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