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