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