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