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