]> 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
coverity: remove unnecesary consts in methods declaration
[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
7cd4e982 1165 Int_t iparton = -1;
1166 for(Int_t i = 0; i< nTriggerJets; i++){
f8006433 1167 iparton=-1;
1168 pygeh->TriggerJet(i, tmpjet);
1169 jet = new TParticle(94, 21, -1, -1, -1, -1, tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3], 0,0,0,0);
1170 //Assign an outgoing parton as mother
1171 Float_t phidiff1 = TMath::Abs(jet->Phi()-parton1->Phi());
1172 Float_t phidiff2 = TMath::Abs(jet->Phi()-parton2->Phi());
1173 if(phidiff1 > phidiff2) jet->SetFirstMother(7);
1174 else jet->SetFirstMother(6);
1175 //jet->Print();
1176 if(fDebug > 1)
1177 printf("AliMCAnalysisUtils::GetJets() - PYTHIA Jet %d: mother %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
1178 i, jet->GetFirstMother(),jet->Pt(),jet->Energy(),jet->Phi()*TMath::RadToDeg(),jet->Eta());
1179 fJetsList->Add(jet);
7cd4e982 1180 }
1181 }//Pythia triggered jets
1182 //HERWIG
1183 else if (fMCGenerator=="HERWIG"){
1184 Int_t pdg = -1;
1185 //Check parton 1
1186 TParticle * tmp = parton1;
1187 if(parton1->GetPdgCode()!=22){
f8006433 1188 while(pdg != 94){
1189 if(tmp->GetFirstDaughter()==-1) return fJetsList;
1190 tmp = stack->Particle(tmp->GetFirstDaughter());
1191 pdg = tmp->GetPdgCode();
1192 }//while
1193
1194 //Add found jet to list
1195 TParticle *jet1 = new TParticle(*tmp);
1196 jet1->SetFirstMother(6);
1197 fJetsList->Add(jet1);
1198 //printf("jet 1: first daughter %d, last daughter %d\n", tmp->GetFirstDaughter(), tmp->GetLastDaughter());
1199 //tmp = stack->Particle(tmp->GetFirstDaughter());
1200 //tmp->Print();
1201 //jet1->Print();
1202 if(fDebug > 1)
1203 printf("AliMCAnalysisUtils::GetJets() - HERWIG Jet 1: mother %d, status %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
1204 jet1->GetFirstMother(),jet1->GetStatusCode(),jet1->Pt(),jet1->Energy(),jet1->Phi()*TMath::RadToDeg(),jet1->Eta());
7cd4e982 1205 }//not photon
1206
1207 //Check parton 2
1208 pdg = -1;
1209 tmp = parton2;
1210 Int_t i = -1;
1211 if(parton2->GetPdgCode()!=22){
f8006433 1212 while(pdg != 94){
1213 if(tmp->GetFirstDaughter()==-1) return fJetsList;
1214 i = tmp->GetFirstDaughter();
1215 tmp = stack->Particle(tmp->GetFirstDaughter());
1216 pdg = tmp->GetPdgCode();
1217 }//while
1218 //Add found jet to list
1219 TParticle *jet2 = new TParticle(*tmp);
1220 jet2->SetFirstMother(7);
1221 fJetsList->Add(jet2);
1222 //jet2->Print();
1223 if(fDebug > 1)
1224 printf("AliMCAnalysisUtils::GetJets() - HERWIG Jet 2: mother %d, status %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
1225 jet2->GetFirstMother(),jet2->GetStatusCode(),jet2->Pt(),jet2->Energy(),jet2->Phi()*TMath::RadToDeg(),jet2->Eta());
1226 //Int_t first = tmp->GetFirstDaughter();
1227 //Int_t last = tmp->GetLastDaughter();
1228 //printf("jet 2: first daughter %d, last daughter %d, pdg %d\n",first, last, tmp->GetPdgCode());
7cd4e982 1229 // for(Int_t d = first ; d < last+1; d++){
f8006433 1230 // tmp = stack->Particle(d);
1231 // if(i == tmp->GetFirstMother())
1232 // printf("Daughter n %d, Mother %d, name %s, status %d, pT %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
1233 // d,tmp->GetFirstMother(), tmp->GetName(), tmp->GetStatusCode(),tmp->Pt(),tmp->Energy(),tmp->Phi()*TMath::RadToDeg(),tmp->Eta());
1234 // }
1235 //tmp->Print();
7cd4e982 1236 }//not photon
1237 }//Herwig generated jets
1238 }
1239
1240 return fJetsList;
1241}
1242
51a0ace5 1243
8a2dbbff 1244//__________________________________________________________________________________________________________
1245TLorentzVector AliMCAnalysisUtils::GetDaughter(Int_t idaugh, Int_t label,
64be16b4 1246 const AliCaloTrackReader* reader,
1247 Int_t & pdg, Int_t & status, Bool_t & ok, Int_t & daughlabel)
1248{
1249 //Return the kinematics of the particle that generated the signal, its pdg and its status and its label mother
1250
1251 TLorentzVector daughter(0,0,0,0);
1252
1253 if(reader->ReadStack())
1254 {
1255 if(!reader->GetStack())
1256 {
1257 if (fDebug >=0)
1258 printf("AliMCAnalysisUtils::GetDaughter() - Stack is not available, check analysis settings in configuration file, STOP!!\n");
1259
1260 ok=kFALSE;
1261 return daughter;
1262 }
1263 if(label >= 0 && label < reader->GetStack()->GetNtrack())
1264 {
1265 TParticle * momP = reader->GetStack()->Particle(label);
1266 daughlabel = momP->GetDaughter(idaugh);
1267
1268 TParticle * daughP = reader->GetStack()->Particle(daughlabel);
1269 daughP->Momentum(daughter);
1270 pdg = daughP->GetPdgCode();
1271 status = daughP->GetStatusCode();
1272 }
1273 else
1274 {
1275 ok = kFALSE;
1276 return daughter;
1277 }
1278 }
1279 else if(reader->ReadAODMCParticles())
1280 {
1281 TClonesArray* mcparticles = reader->GetAODMCParticles();
1282 if(!mcparticles)
1283 {
1284 if(fDebug >= 0)
1285 printf("AliMCAnalysisUtils::GetDaughter() - AODMCParticles is not available, check analysis settings in configuration file!!\n");
1286
1287 ok=kFALSE;
1288 return daughter;
1289 }
1290
1291 Int_t nprimaries = mcparticles->GetEntriesFast();
1292 if(label >= 0 && label < nprimaries)
1293 {
1294 AliAODMCParticle * momP = (AliAODMCParticle *) mcparticles->At(label);
1295 daughlabel = momP->GetDaughter(idaugh);
1296
1297 AliAODMCParticle * daughP = (AliAODMCParticle *) mcparticles->At(daughlabel);
1298 daughter.SetPxPyPzE(daughP->Px(),daughP->Py(),daughP->Pz(),daughP->E());
1299 pdg = daughP->GetPdgCode();
1300 status = daughP->GetStatus();
1301 }
1302 else
1303 {
1304 ok = kFALSE;
1305 return daughter;
1306 }
1307 }
1308
1309 ok = kTRUE;
1310
1311 return daughter;
1312}
1313
8a2dbbff 1314//______________________________________________________________________________________________________
1315TLorentzVector AliMCAnalysisUtils::GetMother(Int_t label, const AliCaloTrackReader* reader, Bool_t & ok)
9036b2a6 1316{
1317 //Return the kinematics of the particle that generated the signal
1318
64be16b4 1319 Int_t pdg = -1; Int_t status = -1; Int_t momlabel = -1;
1320 return GetMother(label,reader,pdg,status, ok,momlabel);
1321}
1322
8a2dbbff 1323//_________________________________________________________________________________________
1324TLorentzVector AliMCAnalysisUtils::GetMother(Int_t label, const AliCaloTrackReader* reader,
64be16b4 1325 Int_t & pdg, Int_t & status, Bool_t & ok)
1326{
1327 //Return the kinematics of the particle that generated the signal
1328
1329 Int_t momlabel = -1;
1330 return GetMother(label,reader,pdg,status, ok,momlabel);
51a0ace5 1331}
1332
8a2dbbff 1333//______________________________________________________________________________________________________
1334TLorentzVector AliMCAnalysisUtils::GetMother(Int_t label, const AliCaloTrackReader* reader,
64be16b4 1335 Int_t & pdg, Int_t & status, Bool_t & ok, Int_t & momlabel)
51a0ace5 1336{
64be16b4 1337 //Return the kinematics of the particle that generated the signal, its pdg and its status and its label mother
51a0ace5 1338
1339 TLorentzVector mom(0,0,0,0);
9036b2a6 1340
1341 if(reader->ReadStack())
1342 {
51a0ace5 1343 if(!reader->GetStack())
1344 {
9036b2a6 1345 if (fDebug >=0)
1346 printf("AliMCAnalysisUtils::GetMother() - Stack is not available, check analysis settings in configuration file, STOP!!\n");
51a0ace5 1347
1348 ok=kFALSE;
1349 return mom;
9036b2a6 1350 }
1351 if(label >= 0 && label < reader->GetStack()->GetNtrack())
1352 {
1353 TParticle * momP = reader->GetStack()->Particle(label);
1354 momP->Momentum(mom);
51a0ace5 1355 pdg = momP->GetPdgCode();
1356 status = momP->GetStatusCode();
64be16b4 1357 momlabel = momP->GetFirstMother();
51a0ace5 1358 }
1359 else
1360 {
1361 ok = kFALSE;
1362 return mom;
9036b2a6 1363 }
1364 }
1365 else if(reader->ReadAODMCParticles())
1366 {
2644ead9 1367 TClonesArray* mcparticles = reader->GetAODMCParticles();
9036b2a6 1368 if(!mcparticles)
1369 {
1370 if(fDebug >= 0)
1371 printf("AliMCAnalysisUtils::GetMother() - AODMCParticles is not available, check analysis settings in configuration file!!\n");
51a0ace5 1372
1373 ok=kFALSE;
1374 return mom;
9036b2a6 1375 }
1376
1377 Int_t nprimaries = mcparticles->GetEntriesFast();
1378 if(label >= 0 && label < nprimaries)
1379 {
1380 AliAODMCParticle * momP = (AliAODMCParticle *) mcparticles->At(label);
1381 mom.SetPxPyPzE(momP->Px(),momP->Py(),momP->Pz(),momP->E());
51a0ace5 1382 pdg = momP->GetPdgCode();
1383 status = momP->GetStatus();
64be16b4 1384 momlabel = momP->GetMother();
1385 }
51a0ace5 1386 else
1387 {
1388 ok = kFALSE;
1389 return mom;
9036b2a6 1390 }
1391 }
1392
51a0ace5 1393 ok = kTRUE;
1394
9036b2a6 1395 return mom;
1396}
1397
8a2dbbff 1398//___________________________________________________________________________________
1399TLorentzVector AliMCAnalysisUtils::GetMotherWithPDG(Int_t label, Int_t pdg,
1400 const AliCaloTrackReader* reader,
36769d30 1401 Bool_t & ok, Int_t & momlabel)
9036b2a6 1402{
1403 //Return the kinematics of the particle that generated the signal
1404
51a0ace5 1405 TLorentzVector grandmom(0,0,0,0);
9036b2a6 1406
1407
1408 if(reader->ReadStack())
1409 {
1410 if(!reader->GetStack())
1411 {
1412 if (fDebug >=0)
1413 printf("AliMCAnalysisUtils::GetMotherWithPDG() - Stack is not available, check analysis settings in configuration file, STOP!!\n");
51a0ace5 1414
1415 ok = kFALSE;
1416 return grandmom;
9036b2a6 1417 }
36769d30 1418
9036b2a6 1419 if(label >= 0 && label < reader->GetStack()->GetNtrack())
1420 {
1421 TParticle * momP = reader->GetStack()->Particle(label);
1422
1423 Int_t grandmomLabel = momP->GetFirstMother();
1424 Int_t grandmomPDG = -1;
1425 TParticle * grandmomP = 0x0;
1426 while (grandmomLabel >=0 )
1427 {
1428 grandmomP = reader->GetStack()->Particle(grandmomLabel);
1429 grandmomPDG = grandmomP->GetPdgCode();
1430 if(grandmomPDG==pdg)
1431 {
36769d30 1432 momlabel = grandmomLabel;
9036b2a6 1433 grandmom.SetPxPyPzE(grandmomP->Px(),grandmomP->Py(),grandmomP->Pz(),grandmomP->Energy());
1434 break;
1435 }
1436
1437 grandmomLabel = grandmomP->GetFirstMother();
1438
1439 }
1440
1441 if(grandmomPDG!=pdg) printf("AliMCAnalysisUtils::GetMotherWithPDG(ESD) - mother with PDG %d, not found! \n",pdg);
1442 }
1443 }
1444 else if(reader->ReadAODMCParticles())
1445 {
2644ead9 1446 TClonesArray* mcparticles = reader->GetAODMCParticles();
9036b2a6 1447 if(!mcparticles)
1448 {
1449 if(fDebug >= 0)
1450 printf("AliMCAnalysisUtils::GetMotherWithPDG() - AODMCParticles is not available, check analysis settings in configuration file!!\n");
51a0ace5 1451
1452 ok=kFALSE;
1453 return grandmom;
9036b2a6 1454 }
1455
1456 Int_t nprimaries = mcparticles->GetEntriesFast();
1457 if(label >= 0 && label < nprimaries)
1458 {
1459 AliAODMCParticle * momP = (AliAODMCParticle *) mcparticles->At(label);
1460
1461 Int_t grandmomLabel = momP->GetMother();
1462 Int_t grandmomPDG = -1;
1463 AliAODMCParticle * grandmomP = 0x0;
1464 while (grandmomLabel >=0 )
1465 {
1466 grandmomP = (AliAODMCParticle *) mcparticles->At(grandmomLabel);
1467 grandmomPDG = grandmomP->GetPdgCode();
1468 if(grandmomPDG==pdg)
1469 {
1470 //printf("AliMCAnalysisUtils::GetMotherWithPDG(AOD) - mother with PDG %d FOUND! \n",pdg);
36769d30 1471 momlabel = grandmomLabel;
9036b2a6 1472 grandmom.SetPxPyPzE(grandmomP->Px(),grandmomP->Py(),grandmomP->Pz(),grandmomP->E());
1473 break;
1474 }
1475
1476 grandmomLabel = grandmomP->GetMother();
1477
1478 }
1479
1480 if(grandmomPDG!=pdg) printf("AliMCAnalysisUtils::GetMotherWithPDG(AOD) - mother with PDG %d, NOT found! \n",pdg);
1481
1482 }
1483 }
1484
51a0ace5 1485 ok = kTRUE;
9036b2a6 1486
1487 return grandmom;
1488}
7cd4e982 1489
8a2dbbff 1490//______________________________________________________________________________________________
1491TLorentzVector AliMCAnalysisUtils::GetGrandMother(Int_t label, const AliCaloTrackReader* reader,
914b9fe7 1492 Int_t & pdg, Int_t & status, Bool_t & ok,
1493 Int_t & grandMomLabel, Int_t & greatMomLabel)
1494{
1495 //Return the kinematics of the particle that generated the signal
1496
1497 TLorentzVector grandmom(0,0,0,0);
1498
1499 if(reader->ReadStack())
1500 {
1501 if(!reader->GetStack())
1502 {
1503 if (fDebug >=0)
64be16b4 1504 printf("AliMCAnalysisUtils::GetGrandMother() - Stack is not available, check analysis settings in configuration file, STOP!!\n");
914b9fe7 1505
1506 ok = kFALSE;
1507 return grandmom;
1508 }
1509
1510 if(label >= 0 && label < reader->GetStack()->GetNtrack())
1511 {
1512 TParticle * momP = reader->GetStack()->Particle(label);
1513
1514 grandMomLabel = momP->GetFirstMother();
1515
1516 TParticle * grandmomP = 0x0;
1517
1518 if (grandMomLabel >=0 )
1519 {
1520 grandmomP = reader->GetStack()->Particle(grandMomLabel);
1521 pdg = grandmomP->GetPdgCode();
1522 status = grandmomP->GetStatusCode();
1523
1524 grandmom.SetPxPyPzE(grandmomP->Px(),grandmomP->Py(),grandmomP->Pz(),grandmomP->Energy());
1525 greatMomLabel = grandmomP->GetFirstMother();
1526
1527 }
1528 }
1529 }
1530 else if(reader->ReadAODMCParticles())
1531 {
1532 TClonesArray* mcparticles = reader->GetAODMCParticles();
1533 if(!mcparticles)
1534 {
1535 if(fDebug >= 0)
64be16b4 1536 printf("AliMCAnalysisUtils::GetGrandMother() - AODMCParticles is not available, check analysis settings in configuration file!!\n");
914b9fe7 1537
1538 ok=kFALSE;
1539 return grandmom;
1540 }
1541
1542 Int_t nprimaries = mcparticles->GetEntriesFast();
1543 if(label >= 0 && label < nprimaries)
1544 {
1545 AliAODMCParticle * momP = (AliAODMCParticle *) mcparticles->At(label);
1546
1547 grandMomLabel = momP->GetMother();
1548
1549 AliAODMCParticle * grandmomP = 0x0;
1550
1551 if(grandMomLabel >=0 )
1552 {
1553 grandmomP = (AliAODMCParticle *) mcparticles->At(grandMomLabel);
1554 pdg = grandmomP->GetPdgCode();
1555 status = grandmomP->GetStatus();
1556
1557 grandmom.SetPxPyPzE(grandmomP->Px(),grandmomP->Py(),grandmomP->Pz(),grandmomP->E());
1558 greatMomLabel = grandmomP->GetMother();
1559
1560 }
1561 }
1562 }
1563
1564 ok = kTRUE;
1565
1566 return grandmom;
1567}
1568
8a2dbbff 1569//_______________________________________________________________________________________________________________
1570void AliMCAnalysisUtils::GetMCDecayAsymmetryAngleForPDG(Int_t label, Int_t pdg, const AliCaloTrackReader* reader,
4efb6957 1571 Float_t & asym, Float_t & angle, Bool_t & ok)
8e81c2cf 1572{
1573 //In case of an eta or pi0 decay into 2 photons, get the asymmetry in the energy of the photons
1574
8e81c2cf 1575
1576 if(reader->ReadStack())
1577 {
1578 if(!reader->GetStack())
1579 {
1580 if (fDebug >=0)
64be16b4 1581 printf("AliMCAnalysisUtils::GetMCDecayAsymmetryForPDG() - Stack is not available, check analysis settings in configuration file, STOP!!\n");
8e81c2cf 1582
1583 ok = kFALSE;
8e81c2cf 1584 }
1585 if(label >= 0 && label < reader->GetStack()->GetNtrack())
1586 {
1587 TParticle * momP = reader->GetStack()->Particle(label);
1588
1589 Int_t grandmomLabel = momP->GetFirstMother();
1590 Int_t grandmomPDG = -1;
1591 TParticle * grandmomP = 0x0;
1592 while (grandmomLabel >=0 )
1593 {
1594 grandmomP = reader->GetStack()->Particle(grandmomLabel);
1595 grandmomPDG = grandmomP->GetPdgCode();
1596
1597 if(grandmomPDG==pdg) break;
1598
1599 grandmomLabel = grandmomP->GetFirstMother();
1600
1601 }
1602
1603 if(grandmomPDG==pdg && grandmomP->GetNDaughters()==2)
1604 {
1605 TParticle * d1 = reader->GetStack()->Particle(grandmomP->GetDaughter(0));
1606 TParticle * d2 = reader->GetStack()->Particle(grandmomP->GetDaughter(1));
1607 if(d1->GetPdgCode() == 22 && d1->GetPdgCode() == 22)
1608 {
1609 asym = (d1->Energy()-d2->Energy())/grandmomP->Energy();
4efb6957 1610 TLorentzVector lvd1, lvd2;
1611 d1->Momentum(lvd1);
1612 d2->Momentum(lvd2);
1613 angle = lvd1.Angle(lvd2.Vect());
8e81c2cf 1614 }
1615 }
1616 else
1617 {
1618 ok=kFALSE;
64be16b4 1619 printf("AliMCAnalysisUtils::GetMCDecayAsymmetryForPDG(ESD) - mother with PDG %d, not found! \n",pdg);
8e81c2cf 1620 }
1621
1622 } // good label
1623 }
1624 else if(reader->ReadAODMCParticles())
1625 {
2644ead9 1626 TClonesArray* mcparticles = reader->GetAODMCParticles();
8e81c2cf 1627 if(!mcparticles)
1628 {
1629 if(fDebug >= 0)
64be16b4 1630 printf("AliMCAnalysisUtils::GetMCDecayAsymmetryForPDG() - AODMCParticles is not available, check analysis settings in configuration file!!\n");
8e81c2cf 1631
1632 ok=kFALSE;
c31af35c 1633 return;
8e81c2cf 1634 }
1635
1636 Int_t nprimaries = mcparticles->GetEntriesFast();
1637 if(label >= 0 && label < nprimaries)
1638 {
1639 AliAODMCParticle * momP = (AliAODMCParticle *) mcparticles->At(label);
1640
1641 Int_t grandmomLabel = momP->GetMother();
1642 Int_t grandmomPDG = -1;
1643 AliAODMCParticle * grandmomP = 0x0;
1644 while (grandmomLabel >=0 )
1645 {
1646 grandmomP = (AliAODMCParticle *) mcparticles->At(grandmomLabel);
1647 grandmomPDG = grandmomP->GetPdgCode();
1648
1649 if(grandmomPDG==pdg) break;
1650
1651 grandmomLabel = grandmomP->GetMother();
1652
1653 }
1654
1655 if(grandmomPDG==pdg && grandmomP->GetNDaughters()==2)
1656 {
1657 AliAODMCParticle * d1 = (AliAODMCParticle *) mcparticles->At(grandmomP->GetDaughter(0));
1658 AliAODMCParticle * d2 = (AliAODMCParticle *) mcparticles->At(grandmomP->GetDaughter(1));
1659 if(d1->GetPdgCode() == 22 && d1->GetPdgCode() == 22)
1660 {
1661 asym = (d1->E()-d2->E())/grandmomP->E();
4efb6957 1662 TLorentzVector lvd1, lvd2;
1663 lvd1.SetPxPyPzE(d1->Px(),d1->Py(),d1->Pz(),d1->E());
1664 lvd2.SetPxPyPzE(d2->Px(),d2->Py(),d2->Pz(),d2->E());
1665 angle = lvd1.Angle(lvd2.Vect());
8e81c2cf 1666 }
1667 }
1668 else
1669 {
1670 ok=kFALSE;
64be16b4 1671 printf("AliMCAnalysisUtils::GetMCDecayAsymmetryForPDG(AOD) - mother with PDG %d, not found! \n",pdg);
8e81c2cf 1672 }
1673
1674 } // good label
1675 }
1676
1677 ok = kTRUE;
1678
8e81c2cf 1679}
1680
1681
8a2dbbff 1682//_________________________________________________________________________________________________
1683Int_t AliMCAnalysisUtils::GetNDaughters(Int_t label, const AliCaloTrackReader* reader, Bool_t & ok)
64be16b4 1684{
1685 // Return the the number of daughters of a given MC particle
1686
1687
1688 if(reader->ReadStack())
1689 {
1690 if(!reader->GetStack())
1691 {
1692 if (fDebug >=0)
1693 printf("AliMCAnalysisUtils::GetNDaughters() - Stack is not available, check analysis settings in configuration file, STOP!!\n");
1694
1695 ok=kFALSE;
1696 return -1;
1697 }
1698 if(label >= 0 && label < reader->GetStack()->GetNtrack())
1699 {
1700 TParticle * momP = reader->GetStack()->Particle(label);
1701 ok=kTRUE;
1702 return momP->GetNDaughters();
1703 }
1704 else
1705 {
1706 ok = kFALSE;
1707 return -1;
1708 }
1709 }
1710 else if(reader->ReadAODMCParticles())
1711 {
1712 TClonesArray* mcparticles = reader->GetAODMCParticles();
1713 if(!mcparticles)
1714 {
1715 if(fDebug >= 0)
1716 printf("AliMCAnalysisUtils::GetNDaughters() - AODMCParticles is not available, check analysis settings in configuration file!!\n");
1717
1718 ok=kFALSE;
1719 return -1;
1720 }
1721
1722 Int_t nprimaries = mcparticles->GetEntriesFast();
1723 if(label >= 0 && label < nprimaries)
1724 {
1725 AliAODMCParticle * momP = (AliAODMCParticle *) mcparticles->At(label);
1726 ok = kTRUE;
1727 return momP->GetNDaughters();
1728 }
1729 else
1730 {
1731 ok = kFALSE;
1732 return -1;
1733 }
1734 }
1735
1736 ok = kFALSE;
1737
1738 return -1;
1739}
1740
8a2dbbff 1741//_________________________________________________________________________________
1742Int_t AliMCAnalysisUtils::GetNOverlaps(const Int_t * label, UInt_t nlabels,
1743 Int_t mctag, Int_t mesonLabel,
4914e781 1744 AliCaloTrackReader * reader, Int_t *overpdg)
1745{
1746 // Compare the primary depositing more energy with the rest,
1747 // if no photon/electron (conversion) or neutral meson as comon ancestor, consider it as other particle contributing
1748 // Give as input the meson label in case it was a pi0 or eta merged cluster
1749 // Init overpdg with nlabels
1750
1751 Int_t ancPDG = 0, ancStatus = -1;
1752 TLorentzVector momentum; TVector3 prodVertex;
1753 Int_t ancLabel = 0;
1754 Int_t noverlaps = 0;
1755 Bool_t ok = kFALSE;
1756
1757 for (UInt_t ilab = 1; ilab < nlabels; ilab++ )
1758 {
1759 ancLabel = CheckCommonAncestor(label[0],label[ilab],reader,ancPDG,ancStatus,momentum,prodVertex);
1760
1761 //printf("Overlaps, i %d: Main Label %d, second label %d, ancestor: Label %d, pdg %d - tag %d \n",
1762 // ilab,label[0],label[ilab],ancLabel,ancPDG, mctag);
1763
1764 Bool_t overlap = kFALSE;
1765
1766 if ( ancLabel < 0 )
1767 {
1768 overlap = kTRUE;
1769 //printf("\t \t \t No Label = %d\n",ancLabel);
1770 }
1771 else if( ( ancPDG==111 || ancPDG==221 ) && ( CheckTagBit(mctag,kMCPi0) || CheckTagBit(mctag,kMCEta)) && mesonLabel != ancLabel)
1772 {
1773 //printf("\t \t meson Label %d, ancestor Label %d\n",mesonLabel,ancLabel);
1774 overlap = kTRUE;
1775 }
1776 else if( ancPDG!=22 && TMath::Abs(ancPDG)!=11 && ancPDG != 111 && ancPDG != 221 )
1777 {
1778 //printf("\t \t \t Non EM PDG = %d\n",ancPDG);
1779 overlap = kTRUE ;
1780 }
1781
1782 if( !overlap ) continue ;
1783
1784 // We have at least one overlap
1785
1786 //printf("Overlap!!!!!!!!!!!!!!\n");
1787
1788 noverlaps++;
1789
1790 // What is the origin of the overlap?
1791 Bool_t mOK = 0, gOK = 0;
1792 Int_t mpdg = -999999, gpdg = -1;
1793 Int_t mstatus = -1, gstatus = -1;
1794 Int_t gLabel = -1, ggLabel = -1;
1795 TLorentzVector mother = GetMother (label[ilab],reader,mpdg,mstatus,mOK);
1796 TLorentzVector grandmother = GetGrandMother(label[ilab],reader,gpdg,gstatus,gOK, gLabel,ggLabel);
1797
1798 //printf("\t Overlap!, mother pdg %d; grand mother pdg %d",mpdg,gpdg);
1799
1800 if( ( mpdg == 22 || TMath::Abs(mpdg==11) ) &&
1801 ( gpdg == 22 || TMath::Abs(gpdg==11) ) &&
1802 gLabel >=0 )
1803 {
1804 Int_t labeltmp = gLabel;
1805 while( ( gpdg == 22 || TMath::Abs(gpdg==11) ) && gLabel >=0 )
1806 {
1807 mpdg=gpdg;
1808 grandmother = GetGrandMother(labeltmp,reader,gpdg,gstatus,ok, gLabel,ggLabel);
1809 labeltmp=gLabel;
1810 }
1811 }
1812 overpdg[noverlaps-1] = mpdg;
1813 }
1814
1815 return noverlaps ;
1816
1817}
1818
c5693f62 1819//________________________________________________________
7cd4e982 1820void AliMCAnalysisUtils::Print(const Option_t * opt) const
1821{
1822 //Print some relevant parameters set for the analysis
f8006433 1823
1824 if(! opt)
1825 return;
1826
1827 printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
1828
1829 printf("Debug level = %d\n",fDebug);
1830 printf("MC Generator = %s\n",fMCGenerator.Data());
1831 printf(" \n");
1832
7cd4e982 1833}
1834
8a2dbbff 1835//__________________________________________________
1836void AliMCAnalysisUtils::PrintMCTag(Int_t tag) const
bbb79e66 1837{
1838 // print the assigned origins to this particle
1839
bfc290d0 1840 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 1841 tag,
1842 CheckTagBit(tag,kMCPhoton),
bfc290d0 1843 CheckTagBit(tag,kMCConversion),
bbb79e66 1844 CheckTagBit(tag,kMCPrompt),
1845 CheckTagBit(tag,kMCFragmentation),
1846 CheckTagBit(tag,kMCISR),
1847 CheckTagBit(tag,kMCPi0Decay),
1848 CheckTagBit(tag,kMCEtaDecay),
1849 CheckTagBit(tag,kMCOtherDecay),
bbb79e66 1850 CheckTagBit(tag,kMCPi0),
bfc290d0 1851 CheckTagBit(tag,kMCEta),
bbb79e66 1852 CheckTagBit(tag,kMCElectron),
bbb79e66 1853 CheckTagBit(tag,kMCMuon),
1854 CheckTagBit(tag,kMCPion),
1855 CheckTagBit(tag,kMCProton),
1856 CheckTagBit(tag,kMCAntiNeutron),
1857 CheckTagBit(tag,kMCKaon),
1858 CheckTagBit(tag,kMCAntiProton),
1859 CheckTagBit(tag,kMCAntiNeutron),
bfc290d0 1860 CheckTagBit(tag,kMCOther),
1861 CheckTagBit(tag,kMCUnknown),
bbb79e66 1862 CheckTagBit(tag,kMCBadLabel)
1863 );
bbb79e66 1864}
1865
7cd4e982 1866