AliPhysics  vAN-20150827 (3e81cbb)
 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, TVector3 & prodVertex)
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  prodVertex.SetXYZ(daughP->Vx(),daughP->Vy(),daughP->Vz());
1484  }
1485  else
1486  {
1487  ok = kFALSE;
1488  return fDaughMom;
1489  }
1490  }
1491  else if(reader->ReadAODMCParticles())
1492  {
1493  TClonesArray* mcparticles = reader->GetAODMCParticles();
1494  if(!mcparticles)
1495  {
1496  AliWarning("AODMCParticles is not available, check analysis settings in configuration file!!");
1497 
1498  ok=kFALSE;
1499  return fDaughMom;
1500  }
1501 
1502  Int_t nprimaries = mcparticles->GetEntriesFast();
1503  if(label >= 0 && label < nprimaries)
1504  {
1505  AliAODMCParticle * momP = (AliAODMCParticle *) mcparticles->At(label);
1506  daughlabel = momP->GetDaughter(idaugh);
1507 
1508  if(daughlabel < 0 || daughlabel >= nprimaries)
1509  {
1510  ok = kFALSE;
1511  return fDaughMom;
1512  }
1513 
1514  AliAODMCParticle * daughP = (AliAODMCParticle *) mcparticles->At(daughlabel);
1515  fDaughMom.SetPxPyPzE(daughP->Px(),daughP->Py(),daughP->Pz(),daughP->E());
1516  pdg = daughP->GetPdgCode();
1517  status = daughP->GetStatus();
1518  prodVertex.SetXYZ(daughP->Xv(),daughP->Yv(),daughP->Zv());
1519  }
1520  else
1521  {
1522  ok = kFALSE;
1523  return fDaughMom;
1524  }
1525  }
1526 
1527  ok = kTRUE;
1528 
1529  return fDaughMom;
1530 }
1531 
1532 //______________________________________________________________________________________________________
1534 //______________________________________________________________________________________________________
1535 TLorentzVector AliMCAnalysisUtils::GetMother(Int_t label, const AliCaloTrackReader* reader, Bool_t & ok)
1536 {
1537  Int_t pdg = -1; Int_t status = -1; Int_t momlabel = -1;
1538 
1539  return GetMother(label,reader,pdg,status, ok,momlabel);
1540 }
1541 
1542 //_________________________________________________________________________________________
1543 // \return the kinematics of the particle that generated the signal.
1544 //_________________________________________________________________________________________
1545 TLorentzVector AliMCAnalysisUtils::GetMother(Int_t label, const AliCaloTrackReader* reader,
1546  Int_t & pdg, Int_t & status, Bool_t & ok)
1547 {
1548  Int_t momlabel = -1;
1549 
1550  return GetMother(label,reader,pdg,status, ok,momlabel);
1551 }
1552 
1553 //______________________________________________________________________________________________________
1555 //______________________________________________________________________________________________________
1556 TLorentzVector AliMCAnalysisUtils::GetMother(Int_t label, const AliCaloTrackReader* reader,
1557  Int_t & pdg, Int_t & status, Bool_t & ok, Int_t & momlabel)
1558 {
1559  fMotherMom.SetPxPyPzE(0,0,0,0);
1560 
1561  if(reader->ReadStack())
1562  {
1563  if(!reader->GetStack())
1564  {
1565  AliWarning("Stack is not available, check analysis settings in configuration file, STOP!!");
1566 
1567  ok=kFALSE;
1568  return fMotherMom;
1569  }
1570  if(label >= 0 && label < reader->GetStack()->GetNtrack())
1571  {
1572  TParticle * momP = reader->GetStack()->Particle(label);
1573  momP->Momentum(fMotherMom);
1574  pdg = momP->GetPdgCode();
1575  status = momP->GetStatusCode();
1576  momlabel = momP->GetFirstMother();
1577  }
1578  else
1579  {
1580  ok = kFALSE;
1581  return fMotherMom;
1582  }
1583  }
1584  else if(reader->ReadAODMCParticles())
1585  {
1586  TClonesArray* mcparticles = reader->GetAODMCParticles();
1587  if(!mcparticles)
1588  {
1589  AliWarning("AODMCParticles is not available, check analysis settings in configuration file!!");
1590 
1591  ok=kFALSE;
1592  return fMotherMom;
1593  }
1594 
1595  Int_t nprimaries = mcparticles->GetEntriesFast();
1596  if(label >= 0 && label < nprimaries)
1597  {
1598  AliAODMCParticle * momP = (AliAODMCParticle *) mcparticles->At(label);
1599  fMotherMom.SetPxPyPzE(momP->Px(),momP->Py(),momP->Pz(),momP->E());
1600  pdg = momP->GetPdgCode();
1601  status = momP->GetStatus();
1602  momlabel = momP->GetMother();
1603  }
1604  else
1605  {
1606  ok = kFALSE;
1607  return fMotherMom;
1608  }
1609  }
1610 
1611  ok = kTRUE;
1612 
1613  return fMotherMom;
1614 }
1615 
1616 //___________________________________________________________________________________
1618 //___________________________________________________________________________________
1619 TLorentzVector AliMCAnalysisUtils::GetMotherWithPDG(Int_t label, Int_t pdg,
1620  const AliCaloTrackReader* reader,
1621  Bool_t & ok, Int_t & momlabel)
1622 {
1623  fGMotherMom.SetPxPyPzE(0,0,0,0);
1624 
1625  if(reader->ReadStack())
1626  {
1627  if(!reader->GetStack())
1628  {
1629  AliWarning("Stack is not available, check analysis settings in configuration file!!");
1630 
1631  ok = kFALSE;
1632  return fGMotherMom;
1633  }
1634 
1635  if(label >= 0 && label < reader->GetStack()->GetNtrack())
1636  {
1637  TParticle * momP = reader->GetStack()->Particle(label);
1638 
1639  if(momP->GetPdgCode()==pdg)
1640  {
1641  AliDebug(2,"PDG of mother is already the one requested!");
1642  fGMotherMom.SetPxPyPzE(momP->Px(),momP->Py(),momP->Pz(),momP->Energy());
1643  ok=kTRUE;
1644  return fGMotherMom;
1645  }
1646 
1647  Int_t grandmomLabel = momP->GetFirstMother();
1648  Int_t grandmomPDG = -1;
1649  TParticle * grandmomP = 0x0;
1650  while (grandmomLabel >=0 )
1651  {
1652  grandmomP = reader->GetStack()->Particle(grandmomLabel);
1653  grandmomPDG = grandmomP->GetPdgCode();
1654  if(grandmomPDG==pdg)
1655  {
1656  momlabel = grandmomLabel;
1657  fGMotherMom.SetPxPyPzE(grandmomP->Px(),grandmomP->Py(),grandmomP->Pz(),grandmomP->Energy());
1658  break;
1659  }
1660 
1661  grandmomLabel = grandmomP->GetFirstMother();
1662 
1663  }
1664 
1665  if(grandmomPDG!=pdg) AliInfo(Form("Mother with PDG %d, NOT found! \n",pdg));
1666  }
1667  }
1668  else if(reader->ReadAODMCParticles())
1669  {
1670  TClonesArray* mcparticles = reader->GetAODMCParticles();
1671  if(!mcparticles)
1672  {
1673  AliWarning("AODMCParticles is not available, check analysis settings in configuration file!!");
1674 
1675  ok=kFALSE;
1676  return fGMotherMom;
1677  }
1678 
1679  Int_t nprimaries = mcparticles->GetEntriesFast();
1680  if(label >= 0 && label < nprimaries)
1681  {
1682  AliAODMCParticle * momP = (AliAODMCParticle *) mcparticles->At(label);
1683 
1684  if(momP->GetPdgCode()==pdg)
1685  {
1686  AliDebug(2,"PDG of mother is already the one requested!");
1687  fGMotherMom.SetPxPyPzE(momP->Px(),momP->Py(),momP->Pz(),momP->E());
1688  ok=kTRUE;
1689  return fGMotherMom;
1690  }
1691 
1692  Int_t grandmomLabel = momP->GetMother();
1693  Int_t grandmomPDG = -1;
1694  AliAODMCParticle * grandmomP = 0x0;
1695  while (grandmomLabel >=0 )
1696  {
1697  grandmomP = (AliAODMCParticle *) mcparticles->At(grandmomLabel);
1698  grandmomPDG = grandmomP->GetPdgCode();
1699  if(grandmomPDG==pdg)
1700  {
1701  //printf("AliMCAnalysisUtils::GetMotherWithPDG(AOD) - mother with PDG %d FOUND! \n",pdg);
1702  momlabel = grandmomLabel;
1703  fGMotherMom.SetPxPyPzE(grandmomP->Px(),grandmomP->Py(),grandmomP->Pz(),grandmomP->E());
1704  break;
1705  }
1706 
1707  grandmomLabel = grandmomP->GetMother();
1708 
1709  }
1710 
1711  if(grandmomPDG!=pdg) AliInfo(Form("Mother with PDG %d, NOT found!",pdg));
1712 
1713  }
1714  }
1715 
1716  ok = kTRUE;
1717 
1718  return fGMotherMom;
1719 }
1720 
1721 //______________________________________________________________________________________________
1723 //______________________________________________________________________________________________
1724 TLorentzVector AliMCAnalysisUtils::GetGrandMother(Int_t label, const AliCaloTrackReader* reader,
1725  Int_t & pdg, Int_t & status, Bool_t & ok,
1726  Int_t & grandMomLabel, Int_t & greatMomLabel)
1727 {
1728  fGMotherMom.SetPxPyPzE(0,0,0,0);
1729 
1730  if(reader->ReadStack())
1731  {
1732  if(!reader->GetStack())
1733  {
1734  AliWarning("Stack is not available, check analysis settings in configuration file, STOP!!");
1735 
1736  ok = kFALSE;
1737  return fGMotherMom;
1738  }
1739 
1740  if(label >= 0 && label < reader->GetStack()->GetNtrack())
1741  {
1742  TParticle * momP = reader->GetStack()->Particle(label);
1743 
1744  grandMomLabel = momP->GetFirstMother();
1745 
1746  TParticle * grandmomP = 0x0;
1747 
1748  if (grandMomLabel >=0 )
1749  {
1750  grandmomP = reader->GetStack()->Particle(grandMomLabel);
1751  pdg = grandmomP->GetPdgCode();
1752  status = grandmomP->GetStatusCode();
1753 
1754  fGMotherMom.SetPxPyPzE(grandmomP->Px(),grandmomP->Py(),grandmomP->Pz(),grandmomP->Energy());
1755  greatMomLabel = grandmomP->GetFirstMother();
1756 
1757  }
1758  }
1759  }
1760  else if(reader->ReadAODMCParticles())
1761  {
1762  TClonesArray* mcparticles = reader->GetAODMCParticles();
1763  if(!mcparticles)
1764  {
1765  AliWarning("AODMCParticles is not available, check analysis settings in configuration file!!");
1766 
1767  ok=kFALSE;
1768  return fGMotherMom;
1769  }
1770 
1771  Int_t nprimaries = mcparticles->GetEntriesFast();
1772  if(label >= 0 && label < nprimaries)
1773  {
1774  AliAODMCParticle * momP = (AliAODMCParticle *) mcparticles->At(label);
1775 
1776  grandMomLabel = momP->GetMother();
1777 
1778  AliAODMCParticle * grandmomP = 0x0;
1779 
1780  if(grandMomLabel >=0 )
1781  {
1782  grandmomP = (AliAODMCParticle *) mcparticles->At(grandMomLabel);
1783  pdg = grandmomP->GetPdgCode();
1784  status = grandmomP->GetStatus();
1785 
1786  fGMotherMom.SetPxPyPzE(grandmomP->Px(),grandmomP->Py(),grandmomP->Pz(),grandmomP->E());
1787  greatMomLabel = grandmomP->GetMother();
1788 
1789  }
1790  }
1791  }
1792 
1793  ok = kTRUE;
1794 
1795  return fGMotherMom;
1796 }
1797 
1798 //_______________________________________________________________________________________________________________
1800 //_______________________________________________________________________________________________________________
1802  Float_t & asym, Float_t & angle, Bool_t & ok)
1803 {
1804  if(reader->ReadStack())
1805  {
1806  if(!reader->GetStack())
1807  {
1808  AliWarning("Stack is not available, check analysis settings in configuration file, STOP!!");
1809 
1810  ok = kFALSE;
1811  }
1812  if(label >= 0 && label < reader->GetStack()->GetNtrack())
1813  {
1814  TParticle * momP = reader->GetStack()->Particle(label);
1815 
1816  Int_t grandmomLabel = momP->GetFirstMother();
1817  Int_t grandmomPDG = -1;
1818  TParticle * grandmomP = 0x0;
1819  while (grandmomLabel >=0 )
1820  {
1821  grandmomP = reader->GetStack()->Particle(grandmomLabel);
1822  grandmomPDG = grandmomP->GetPdgCode();
1823 
1824  if(grandmomPDG==pdg) break;
1825 
1826  grandmomLabel = grandmomP->GetFirstMother();
1827 
1828  }
1829 
1830  if(grandmomPDG==pdg && grandmomP->GetNDaughters()==2)
1831  {
1832  TParticle * d1 = reader->GetStack()->Particle(grandmomP->GetDaughter(0));
1833  TParticle * d2 = reader->GetStack()->Particle(grandmomP->GetDaughter(1));
1834  if(d1->GetPdgCode() == 22 && d1->GetPdgCode() == 22)
1835  {
1836  asym = (d1->Energy()-d2->Energy())/grandmomP->Energy();
1837  d1->Momentum(fDaughMom );
1838  d2->Momentum(fDaughMom2);
1839  angle = fDaughMom.Angle(fDaughMom2.Vect());
1840  }
1841  }
1842  else
1843  {
1844  ok=kFALSE;
1845  AliInfo(Form("Mother with PDG %d, not found!",pdg));
1846  }
1847 
1848  } // good label
1849  }
1850  else if(reader->ReadAODMCParticles())
1851  {
1852  TClonesArray* mcparticles = reader->GetAODMCParticles();
1853  if(!mcparticles)
1854  {
1855  AliWarning("AODMCParticles is not available, check analysis settings in configuration file!!");
1856 
1857  ok=kFALSE;
1858  return;
1859  }
1860 
1861  Int_t nprimaries = mcparticles->GetEntriesFast();
1862  if(label >= 0 && label < nprimaries)
1863  {
1864  AliAODMCParticle * momP = (AliAODMCParticle *) mcparticles->At(label);
1865 
1866  Int_t grandmomLabel = momP->GetMother();
1867  Int_t grandmomPDG = -1;
1868  AliAODMCParticle * grandmomP = 0x0;
1869  while (grandmomLabel >=0 )
1870  {
1871  grandmomP = (AliAODMCParticle *) mcparticles->At(grandmomLabel);
1872  grandmomPDG = grandmomP->GetPdgCode();
1873 
1874  if(grandmomPDG==pdg) break;
1875 
1876  grandmomLabel = grandmomP->GetMother();
1877 
1878  }
1879 
1880  if(grandmomPDG==pdg && grandmomP->GetNDaughters()==2)
1881  {
1882  AliAODMCParticle * d1 = (AliAODMCParticle *) mcparticles->At(grandmomP->GetDaughter(0));
1883  AliAODMCParticle * d2 = (AliAODMCParticle *) mcparticles->At(grandmomP->GetDaughter(1));
1884  if(d1->GetPdgCode() == 22 && d1->GetPdgCode() == 22)
1885  {
1886  asym = (d1->E()-d2->E())/grandmomP->E();
1887  fDaughMom .SetPxPyPzE(d1->Px(),d1->Py(),d1->Pz(),d1->E());
1888  fDaughMom2.SetPxPyPzE(d2->Px(),d2->Py(),d2->Pz(),d2->E());
1889  angle = fDaughMom.Angle(fDaughMom2.Vect());
1890  }
1891  }
1892  else
1893  {
1894  ok=kFALSE;
1895  AliInfo(Form("Mother with PDG %d, not found! \n",pdg));
1896  }
1897 
1898  } // good label
1899  }
1900 
1901  ok = kTRUE;
1902 
1903 }
1904 
1905 //_________________________________________________________________________________________________
1907 //_________________________________________________________________________________________________
1908 Int_t AliMCAnalysisUtils::GetNDaughters(Int_t label, const AliCaloTrackReader* reader, Bool_t & ok)
1909 {
1910  if(reader->ReadStack())
1911  {
1912  if(!reader->GetStack())
1913  {
1914  AliWarning("Stack is not available, check analysis settings in configuration file, STOP!!");
1915 
1916  ok=kFALSE;
1917  return -1;
1918  }
1919  if(label >= 0 && label < reader->GetStack()->GetNtrack())
1920  {
1921  TParticle * momP = reader->GetStack()->Particle(label);
1922  ok=kTRUE;
1923  return momP->GetNDaughters();
1924  }
1925  else
1926  {
1927  ok = kFALSE;
1928  return -1;
1929  }
1930  }
1931  else if(reader->ReadAODMCParticles())
1932  {
1933  TClonesArray* mcparticles = reader->GetAODMCParticles();
1934  if(!mcparticles)
1935  {
1936  AliWarning("AODMCParticles is not available, check analysis settings in configuration file!!");
1937 
1938  ok=kFALSE;
1939  return -1;
1940  }
1941 
1942  Int_t nprimaries = mcparticles->GetEntriesFast();
1943  if(label >= 0 && label < nprimaries)
1944  {
1945  AliAODMCParticle * momP = (AliAODMCParticle *) mcparticles->At(label);
1946  ok = kTRUE;
1947  return momP->GetNDaughters();
1948  }
1949  else
1950  {
1951  ok = kFALSE;
1952  return -1;
1953  }
1954  }
1955 
1956  ok = kFALSE;
1957 
1958  return -1;
1959 }
1960 
1961 //_________________________________________________________________________________
1966 //_________________________________________________________________________________
1967 Int_t AliMCAnalysisUtils::GetNOverlaps(const Int_t * label, UInt_t nlabels,
1968  Int_t mctag, Int_t mesonLabel,
1969  AliCaloTrackReader * reader, Int_t *overpdg)
1970 {
1971  Int_t ancPDG = 0, ancStatus = -1;
1972  TVector3 prodVertex;
1973  Int_t ancLabel = 0;
1974  Int_t noverlaps = 0;
1975  Bool_t ok = kFALSE;
1976 
1977  for (UInt_t ilab = 1; ilab < nlabels; ilab++ )
1978  {
1979  ancLabel = CheckCommonAncestor(label[0],label[ilab],reader,ancPDG,ancStatus,fMotherMom,prodVertex);
1980 
1981  //printf("Overlaps, i %d: Main Label %d, second label %d, ancestor: Label %d, pdg %d - tag %d \n",
1982  // ilab,label[0],label[ilab],ancLabel,ancPDG, mctag);
1983 
1984  Bool_t overlap = kFALSE;
1985 
1986  if ( ancLabel < 0 )
1987  {
1988  overlap = kTRUE;
1989  //printf("\t \t \t No Label = %d\n",ancLabel);
1990  }
1991  else if( ( ancPDG==111 || ancPDG==221 ) && ( CheckTagBit(mctag,kMCPi0) || CheckTagBit(mctag,kMCEta)) && mesonLabel != ancLabel)
1992  {
1993  //printf("\t \t meson Label %d, ancestor Label %d\n",mesonLabel,ancLabel);
1994  overlap = kTRUE;
1995  }
1996  else if( ancPDG!=22 && TMath::Abs(ancPDG)!=11 && ancPDG != 111 && ancPDG != 221 )
1997  {
1998  //printf("\t \t \t Non EM PDG = %d\n",ancPDG);
1999  overlap = kTRUE ;
2000  }
2001 
2002  if( !overlap ) continue ;
2003 
2004  // We have at least one overlap
2005 
2006  //printf("Overlap!!!!!!!!!!!!!!\n");
2007 
2008  noverlaps++;
2009 
2010  // What is the origin of the overlap?
2011  Bool_t mOK = 0, gOK = 0;
2012  Int_t mpdg = -999999, gpdg = -1;
2013  Int_t mstatus = -1, gstatus = -1;
2014  Int_t gLabel = -1, ggLabel = -1;
2015 
2016  GetMother (label[ilab],reader,mpdg,mstatus,mOK);
2017  fGMotherMom =
2018  GetGrandMother(label[ilab],reader,gpdg,gstatus,gOK, gLabel,ggLabel);
2019 
2020  //printf("\t Overlap!, mother pdg %d; grand mother pdg %d",mpdg,gpdg);
2021 
2022  if( ( mpdg == 22 || TMath::Abs(mpdg==11) ) &&
2023  ( gpdg == 22 || TMath::Abs(gpdg==11) ) &&
2024  gLabel >=0 )
2025  {
2026  Int_t labeltmp = gLabel;
2027  while( ( gpdg == 22 || TMath::Abs(gpdg==11) ) && gLabel >=0 )
2028  {
2029  mpdg=gpdg;
2030  fGMotherMom = GetGrandMother(labeltmp,reader,gpdg,gstatus,ok, gLabel,ggLabel);
2031  labeltmp=gLabel;
2032  }
2033  }
2034  overpdg[noverlaps-1] = mpdg;
2035  }
2036 
2037  return noverlaps ;
2038 }
2039 
2040 //________________________________________________________
2042 //________________________________________________________
2043 void AliMCAnalysisUtils::Print(const Option_t * opt) const
2044 {
2045  if(! opt)
2046  return;
2047 
2048  printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
2049 
2050  printf("Debug level = %d\n",fDebug);
2051  printf("MC Generator = %s\n",fMCGeneratorString.Data());
2052  printf(" \n");
2053 }
2054 
2055 //__________________________________________________
2057 //__________________________________________________
2058 void AliMCAnalysisUtils::PrintMCTag(Int_t tag) const
2059 {
2060  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",
2061  tag,
2062  CheckTagBit(tag,kMCPhoton),
2064  CheckTagBit(tag,kMCPrompt),
2066  CheckTagBit(tag,kMCISR),
2067  CheckTagBit(tag,kMCPi0Decay),
2068  CheckTagBit(tag,kMCEtaDecay),
2070  CheckTagBit(tag,kMCPi0),
2071  CheckTagBit(tag,kMCEta),
2072  CheckTagBit(tag,kMCElectron),
2073  CheckTagBit(tag,kMCMuon),
2074  CheckTagBit(tag,kMCPion),
2075  CheckTagBit(tag,kMCProton),
2077  CheckTagBit(tag,kMCKaon),
2078  CheckTagBit(tag,kMCAntiProton),
2080  CheckTagBit(tag,kMCUnknown),
2082  );
2083 }
2084 
2085 //__________________________________________________
2087 //__________________________________________________
2089 {
2090  fMCGenerator = mcgen ;
2091  if (mcgen == kPythia) fMCGeneratorString = "PYTHIA";
2092  else if(mcgen == kHerwig) fMCGeneratorString = "HERWIG";
2093  else if(mcgen == kHijing) fMCGeneratorString = "HIJING";
2094  else
2095  {
2096  fMCGeneratorString = "";
2098  }
2099 }
2100 
2101 //____________________________________________________
2103 //____________________________________________________
2105 {
2106  fMCGeneratorString = mcgen ;
2107 
2108  if (mcgen == "PYTHIA") fMCGenerator = kPythia;
2109  else if(mcgen == "HERWIG") fMCGenerator = kHerwig;
2110  else if(mcgen == "HIJING") fMCGenerator = kHijing;
2111  else
2112  {
2114  fMCGeneratorString = "" ;
2115  }
2116 }
2117 
2118 
2119 
TLorentzVector GetMotherWithPDG(Int_t label, Int_t pdg, const AliCaloTrackReader *reader, Bool_t &ok, Int_t &momLabel)
Int_t pdg
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
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 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.
TLorentzVector GetDaughter(Int_t daughter, Int_t label, const AliCaloTrackReader *reader, Int_t &pdg, Int_t &status, Bool_t &ok, Int_t &daugLabel, TVector3 &prodVertex)
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.
TString meson
Int_t CheckCommonAncestor(Int_t index1, Int_t index2, const AliCaloTrackReader *reader, Int_t &ancPDG, Int_t &ancStatus, TLorentzVector &momentum, TVector3 &prodVertex)
TList * GetJets(const AliCaloTrackReader *reader)
Bool_t CheckTagBit(Int_t tag, UInt_t test) const