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