AliPhysics  2797316 (2797316)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
1134  if ( !mother )
1135  {
1136  AliInfo(Form("MC Mother not available for label %d",label));
1137  continue;
1138  }
1139 
1140  Int_t momPDG = TMath::Abs(mother->GetPdgCode());
1141  if ( momPDG!=11 && momPDG!=22 ) continue;
1142 
1143  // Check if "mother" of entity is converted, if not, get the first non converted mother
1144  Int_t iParentClus = mother->GetFirstMother();
1145  if ( iParentClus < 0 ) continue;
1146 
1147  TParticle * parentClus = stack->Particle(iParentClus);
1148  if ( !parentClus ) continue;
1149 
1150  Int_t parentClusPDG = TMath::Abs(parentClus->GetPdgCode());
1151  Int_t parentClusStatus = parentClus->GetStatusCode();
1152 
1153  if( parentClusPDG != 22 && parentClusPDG != 11 && parentClusStatus != 0 ) continue;
1154 
1155  //printf("Conversion\n");
1156 
1157  // Check if the mother is photon or electron with status not stable
1158  while ( (parentClusPDG == 22 || parentClusPDG == 11) && parentClusStatus != 1 )
1159  {
1160  //New Mother
1161  label = iParentClus;
1162  momPDG = parentClusPDG;
1163 
1164  iParentClus = parentClus->GetFirstMother();
1165  if(iParentClus < 0) break;
1166 
1167  parentClus = stack->Particle(iParentClus);
1168  if(!parentClus) break;
1169 
1170  parentClusPDG = TMath::Abs(parentClus->GetPdgCode());
1171  parentClusStatus = parentClus->GetStatusCode() ;
1172  }//while
1173 
1174  if ( (momPDG == 22 || parentClusPDG ==22) && (label==pairLabel || iParentClus == pairLabel) )
1175  {
1177  //printf("Conversion is paired\n");
1178  return ;
1179  }
1180  else continue;
1181 
1182  }
1183  }
1184  } // cluster loop
1185 
1187 }
1188 
1189 //________________________________________________________________________________________________________
1191 //________________________________________________________________________________________________________
1192 void AliMCAnalysisUtils::CheckLostDecayPair(const TObjArray * arrayCluster,Int_t iMom, Int_t iParent,
1193  const TClonesArray* mcparticles, Int_t & tag)
1194 {
1195  if(!arrayCluster || iMom < 0 || iParent < 0|| !mcparticles) return;
1196 
1197  AliAODMCParticle * parent = (AliAODMCParticle*) mcparticles->At(iParent);
1198 
1199  //printf("*** Check label %d with parent %d\n",iMom, iParent);
1200 
1201  if(parent->GetNDaughters()!=2)
1202  {
1204  //printf("\t ndaugh = %d\n",parent->GetNDaughters());
1205  return ;
1206  }
1207 
1208  Int_t pairLabel = -1;
1209  if ( iMom != parent->GetDaughter(0) ) pairLabel = parent->GetDaughter(0);
1210  else if( iMom != parent->GetDaughter(1) ) pairLabel = parent->GetDaughter(1);
1211 
1212  if(pairLabel<0)
1213  {
1214  //printf("\t pair Label not found = %d\n",pairLabel);
1216  return ;
1217  }
1218 
1219  //printf("\t *** find pair %d\n",pairLabel);
1220 
1221  for(Int_t iclus = 0; iclus < arrayCluster->GetEntriesFast(); iclus++)
1222  {
1223  AliVCluster * cluster = (AliVCluster*) arrayCluster->At(iclus);
1224  //printf("\t \t ** Cluster %d, nlabels %d\n",iclus,cluster->GetNLabels());
1225  for(UInt_t ilab = 0; ilab< cluster->GetNLabels(); ilab++)
1226  {
1227  Int_t label = cluster->GetLabels()[ilab];
1228 
1229  //printf("\t \t label %d\n",label);
1230 
1231  if ( label==pairLabel )
1232  {
1233  //printf("\t \t Pair found\n");
1235  return ;
1236  }
1237  else if ( label== iParent || label== iMom )
1238  {
1239  //printf("\t \t skip\n");
1240  continue;
1241  }
1242  else // check the ancestry
1243  {
1244  AliAODMCParticle * mother = (AliAODMCParticle*) mcparticles->At(label);
1245 
1246  if ( !mother )
1247  {
1248  AliInfo(Form("MC Mother not available for label %d",label));
1249  continue;
1250  }
1251 
1252  Int_t momPDG = TMath::Abs(mother->GetPdgCode());
1253  if ( momPDG!=11 && momPDG!=22 ) continue;
1254 
1255  // Check if "mother" of entity is converted, if not, get the first non converted mother
1256  Int_t iParentClus = mother->GetMother();
1257  if(iParentClus < 0) continue;
1258 
1259  AliAODMCParticle * parentClus = (AliAODMCParticle*) mcparticles->At(iParentClus);
1260  if(!parentClus) continue;
1261 
1262  Int_t parentClusPDG = TMath::Abs(parentClus->GetPdgCode());
1263  Int_t parentClusStatus = parentClus->GetStatus();
1264 
1265  if ( parentClusPDG != 22 && parentClusPDG != 11 && parentClusStatus != 0 )
1266  {
1267  //printf("\t \t skip, not a conversion, parent: pdg %d, status %d\n",parentClusPDG,parentClusStatus);
1268  continue;
1269  }
1270 
1271  //printf("\t \t Conversion\n");
1272 
1273  // Check if the mother is photon or electron with status not stable
1274  while ((parentClusPDG == 22 || parentClusPDG == 11) && parentClusStatus != 1)
1275  {
1276  //New Mother
1277  label = iParentClus;
1278  momPDG = parentClusPDG;
1279 
1280  iParentClus = parentClus->GetMother();
1281  if ( iParentClus < 0 ) break;
1282 
1283  parentClus = (AliAODMCParticle*) mcparticles->At(iParentClus);
1284  if ( !parentClus ) break;
1285 
1286  parentClusPDG = TMath::Abs(parentClus->GetPdgCode());
1287  parentClusStatus = parentClus->GetStatus() ;
1288  }//while
1289 
1290  if ( (momPDG == 22 || parentClusPDG ==22) && (label==pairLabel || iParentClus == pairLabel) )
1291  {
1293  //printf("\t \t Conversion is paired: mom %d, parent %d\n",label,iParentClus);
1294  return ;
1295  }
1296  else
1297  {
1298  //printf("\t \t Skip, finally label %d, pdg %d, parent label %d, pdg %d, status %d\n",label,momPDG,iParentClus,parentClusPDG,parentClusStatus);
1299  continue;
1300  }
1301 
1302  }
1303  }
1304  } // cluster loop
1305 
1307 }
1308 
1309 //_____________________________________________________________________
1311 //_____________________________________________________________________
1313 {
1314  AliStack * stack = reader->GetStack();
1315  Int_t iEvent = reader->GetEventNumber();
1316  AliGenEventHeader * geh = reader->GetGenEventHeader();
1317  if(fCurrentEvent!=iEvent){
1318  fCurrentEvent = iEvent;
1319  fJetsList = new TList;
1320  Int_t nTriggerJets = 0;
1321  Float_t tmpjet[]={0,0,0,0};
1322 
1323  //printf("Event %d %d\n",fCurrentEvent,iEvent);
1324  //Get outgoing partons
1325  if(stack->GetNtrack() < 8) return fJetsList;
1326  TParticle * parton1 = stack->Particle(6);
1327  TParticle * parton2 = stack->Particle(7);
1328 
1329  AliDebug(2,Form("Parton 6 : %s, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f",
1330  parton1->GetName(),parton1->Pt(),parton1->Energy(),parton1->Phi()*TMath::RadToDeg(),parton1->Eta()));
1331  AliDebug(2,Form("Parton 7 : %s, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f",
1332  parton2->GetName(),parton2->Pt(),parton2->Energy(),parton2->Phi()*TMath::RadToDeg(),parton2->Eta()));
1333 
1334  // //Trace the jet from the mother parton
1335  // Float_t pt = 0;
1336  // Float_t pt1 = 0;
1337  // Float_t pt2 = 0;
1338  // Float_t e = 0;
1339  // Float_t e1 = 0;
1340  // Float_t e2 = 0;
1341  // TParticle * tmptmp = new TParticle;
1342  // for(Int_t i = 0; i< stack->GetNprimary(); i++){
1343  // tmptmp = stack->Particle(i);
1344 
1345  // if(tmptmp->GetStatusCode() == 1){
1346  // pt = tmptmp->Pt();
1347  // e = tmptmp->Energy();
1348  // Int_t imom = tmptmp->GetFirstMother();
1349  // Int_t imom1 = 0;
1350  // //printf("1st imom %d\n",imom);
1351  // while(imom > 5){
1352  // imom1=imom;
1353  // tmptmp = stack->Particle(imom);
1354  // imom = tmptmp->GetFirstMother();
1355  // //printf("imom %d \n",imom);
1356  // }
1357  // //printf("Last imom %d %d\n",imom1, imom);
1358  // if(imom1 == 6) {
1359  // pt1+=pt;
1360  // e1+=e;
1361  // }
1362  // else if (imom1 == 7){
1363  // pt2+=pt;
1364  // e2+=e; }
1365  // }// status 1
1366 
1367  // }// for
1368 
1369  // printf("JET 1, pt %2.2f, e %2.2f; JET 2, pt %2.2f, e %2.2f \n",pt1,e1,pt2,e2);
1370 
1371  //Get the jet, different way for different generator
1372  //PYTHIA
1373  if(fMCGenerator == kPythia)
1374  {
1375  TParticle * jet = 0x0;
1376  AliGenPythiaEventHeader* pygeh= (AliGenPythiaEventHeader*) geh;
1377  nTriggerJets = pygeh->NTriggerJets();
1378  AliDebug(2,Form("PythiaEventHeader: Njets: %d",nTriggerJets));
1379 
1380  for(Int_t i = 0; i< nTriggerJets; i++)
1381  {
1382  pygeh->TriggerJet(i, tmpjet);
1383  jet = new TParticle(94, 21, -1, -1, -1, -1, tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3], 0,0,0,0);
1384  //Assign an outgoing parton as mother
1385  Float_t phidiff1 = TMath::Abs(jet->Phi()-parton1->Phi());
1386  Float_t phidiff2 = TMath::Abs(jet->Phi()-parton2->Phi());
1387  if(phidiff1 > phidiff2) jet->SetFirstMother(7);
1388  else jet->SetFirstMother(6);
1389  //jet->Print();
1390  AliDebug(1,Form("PYTHIA Jet %d: mother %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f",
1391  i, jet->GetFirstMother(),jet->Pt(),jet->Energy(),jet->Phi()*TMath::RadToDeg(),jet->Eta()));
1392  fJetsList->Add(jet);
1393  }
1394  }//Pythia triggered jets
1395  //HERWIG
1396  else if (fMCGenerator == kHerwig)
1397  {
1398  Int_t pdg = -1;
1399  //Check parton 1
1400  TParticle * tmp = parton1;
1401  if(parton1->GetPdgCode()!=22)
1402  {
1403  while(pdg != 94){
1404  if(tmp->GetFirstDaughter()==-1) return fJetsList;
1405  tmp = stack->Particle(tmp->GetFirstDaughter());
1406  pdg = tmp->GetPdgCode();
1407  }//while
1408 
1409  //Add found jet to list
1410  TParticle *jet1 = new TParticle(*tmp);
1411  jet1->SetFirstMother(6);
1412  fJetsList->Add(jet1);
1413  //printf("jet 1: first daughter %d, last daughter %d\n", tmp->GetFirstDaughter(), tmp->GetLastDaughter());
1414  //tmp = stack->Particle(tmp->GetFirstDaughter());
1415  //tmp->Print();
1416  //jet1->Print();
1417  AliDebug(1,Form("HERWIG Jet 1: mother %d, status %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f",
1418  jet1->GetFirstMother(),jet1->GetStatusCode(),jet1->Pt(),jet1->Energy(),jet1->Phi()*TMath::RadToDeg(),jet1->Eta()));
1419  }//not photon
1420 
1421  //Check parton 2
1422  pdg = -1;
1423  tmp = parton2;
1424  if(parton2->GetPdgCode()!=22)
1425  {
1426  while(pdg != 94)
1427  {
1428  if(tmp->GetFirstDaughter()==-1) return fJetsList;
1429  tmp = stack->Particle(tmp->GetFirstDaughter());
1430  pdg = tmp->GetPdgCode();
1431  }//while
1432 
1433  //Add found jet to list
1434  TParticle *jet2 = new TParticle(*tmp);
1435  jet2->SetFirstMother(7);
1436  fJetsList->Add(jet2);
1437  //jet2->Print();
1438  AliDebug(2,Form("HERWIG Jet 2: mother %d, status %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f",
1439  jet2->GetFirstMother(),jet2->GetStatusCode(),jet2->Pt(),jet2->Energy(),jet2->Phi()*TMath::RadToDeg(),jet2->Eta()));
1440  //Int_t first = tmp->GetFirstDaughter();
1441  //Int_t last = tmp->GetLastDaughter();
1442  //printf("jet 2: first daughter %d, last daughter %d, pdg %d\n",first, last, tmp->GetPdgCode());
1443  // for(Int_t d = first ; d < last+1; d++){
1444  // tmp = stack->Particle(d);
1445  // if(i == tmp->GetFirstMother())
1446  // printf("Daughter n %d, Mother %d, name %s, status %d, pT %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
1447  // d,tmp->GetFirstMother(), tmp->GetName(), tmp->GetStatusCode(),tmp->Pt(),tmp->Energy(),tmp->Phi()*TMath::RadToDeg(),tmp->Eta());
1448  // }
1449  //tmp->Print();
1450  }//not photon
1451  }//Herwig generated jets
1452  }
1453 
1454  return fJetsList;
1455 }
1456 
1457 //______________________________________________________________________________
1459 //______________________________________________________________________________
1460 TLorentzVector AliMCAnalysisUtils::GetDaughter(Int_t idaugh, Int_t label,
1461  const AliCaloTrackReader* reader,
1462  Int_t & pdg, Int_t & status,
1463  Bool_t & ok, Int_t & daughlabel, TVector3 & prodVertex)
1464 {
1465  fDaughMom.SetPxPyPzE(0,0,0,0);
1466 
1467  if(reader->ReadStack())
1468  {
1469  if(!reader->GetStack())
1470  {
1471  AliWarning("Stack is not available, check analysis settings in configuration file!!");
1472 
1473  ok=kFALSE;
1474  return fDaughMom;
1475  }
1476 
1477  Int_t nprimaries = reader->GetStack()->GetNtrack();
1478  if(label >= 0 && label < nprimaries)
1479  {
1480  TParticle * momP = reader->GetStack()->Particle(label);
1481  daughlabel = momP->GetDaughter(idaugh);
1482 
1483  if(daughlabel < 0 || daughlabel >= nprimaries)
1484  {
1485  ok = kFALSE;
1486  return fDaughMom;
1487  }
1488 
1489  TParticle * daughP = reader->GetStack()->Particle(daughlabel);
1490  daughP->Momentum(fDaughMom);
1491  pdg = daughP->GetPdgCode();
1492  status = daughP->GetStatusCode();
1493  prodVertex.SetXYZ(daughP->Vx(),daughP->Vy(),daughP->Vz());
1494  }
1495  else
1496  {
1497  ok = kFALSE;
1498  return fDaughMom;
1499  }
1500  }
1501  else if(reader->ReadAODMCParticles())
1502  {
1503  TClonesArray* mcparticles = reader->GetAODMCParticles();
1504  if(!mcparticles)
1505  {
1506  AliWarning("AODMCParticles is not available, check analysis settings in configuration file!!");
1507 
1508  ok=kFALSE;
1509  return fDaughMom;
1510  }
1511 
1512  Int_t nprimaries = mcparticles->GetEntriesFast();
1513  if(label >= 0 && label < nprimaries)
1514  {
1515  AliAODMCParticle * momP = (AliAODMCParticle *) mcparticles->At(label);
1516  daughlabel = momP->GetDaughter(idaugh);
1517 
1518  if(daughlabel < 0 || daughlabel >= nprimaries)
1519  {
1520  ok = kFALSE;
1521  return fDaughMom;
1522  }
1523 
1524  AliAODMCParticle * daughP = (AliAODMCParticle *) mcparticles->At(daughlabel);
1525  fDaughMom.SetPxPyPzE(daughP->Px(),daughP->Py(),daughP->Pz(),daughP->E());
1526  pdg = daughP->GetPdgCode();
1527  status = daughP->GetStatus();
1528  prodVertex.SetXYZ(daughP->Xv(),daughP->Yv(),daughP->Zv());
1529  }
1530  else
1531  {
1532  ok = kFALSE;
1533  return fDaughMom;
1534  }
1535  }
1536 
1537  ok = kTRUE;
1538 
1539  return fDaughMom;
1540 }
1541 
1542 //______________________________________________________________________________________________________
1544 //______________________________________________________________________________________________________
1545 TLorentzVector AliMCAnalysisUtils::GetMother(Int_t label, const AliCaloTrackReader* reader, Bool_t & ok)
1546 {
1547  Int_t pdg = -1; Int_t status = -1; Int_t momlabel = -1;
1548 
1549  return GetMother(label,reader,pdg,status, ok,momlabel);
1550 }
1551 
1552 //_________________________________________________________________________________________
1553 // \return the kinematics of the particle that generated the signal.
1554 //_________________________________________________________________________________________
1555 TLorentzVector AliMCAnalysisUtils::GetMother(Int_t label, const AliCaloTrackReader* reader,
1556  Int_t & pdg, Int_t & status, Bool_t & ok)
1557 {
1558  Int_t momlabel = -1;
1559 
1560  return GetMother(label,reader,pdg,status, ok,momlabel);
1561 }
1562 
1563 //______________________________________________________________________________________________________
1565 //______________________________________________________________________________________________________
1566 TLorentzVector AliMCAnalysisUtils::GetMother(Int_t label, const AliCaloTrackReader* reader,
1567  Int_t & pdg, Int_t & status, Bool_t & ok, Int_t & momlabel)
1568 {
1569  fMotherMom.SetPxPyPzE(0,0,0,0);
1570 
1571  if(reader->ReadStack())
1572  {
1573  if(!reader->GetStack())
1574  {
1575  AliWarning("Stack is not available, check analysis settings in configuration file, STOP!!");
1576 
1577  ok=kFALSE;
1578  return fMotherMom;
1579  }
1580  if(label >= 0 && label < reader->GetStack()->GetNtrack())
1581  {
1582  TParticle * momP = reader->GetStack()->Particle(label);
1583  momP->Momentum(fMotherMom);
1584  pdg = momP->GetPdgCode();
1585  status = momP->GetStatusCode();
1586  momlabel = momP->GetFirstMother();
1587  }
1588  else
1589  {
1590  ok = kFALSE;
1591  return fMotherMom;
1592  }
1593  }
1594  else if(reader->ReadAODMCParticles())
1595  {
1596  TClonesArray* mcparticles = reader->GetAODMCParticles();
1597  if(!mcparticles)
1598  {
1599  AliWarning("AODMCParticles is not available, check analysis settings in configuration file!!");
1600 
1601  ok=kFALSE;
1602  return fMotherMom;
1603  }
1604 
1605  Int_t nprimaries = mcparticles->GetEntriesFast();
1606  if(label >= 0 && label < nprimaries)
1607  {
1608  AliAODMCParticle * momP = (AliAODMCParticle *) mcparticles->At(label);
1609  fMotherMom.SetPxPyPzE(momP->Px(),momP->Py(),momP->Pz(),momP->E());
1610  pdg = momP->GetPdgCode();
1611  status = momP->GetStatus();
1612  momlabel = momP->GetMother();
1613  }
1614  else
1615  {
1616  ok = kFALSE;
1617  return fMotherMom;
1618  }
1619  }
1620 
1621  ok = kTRUE;
1622 
1623  return fMotherMom;
1624 }
1625 
1626 //___________________________________________________________________________________
1628 //___________________________________________________________________________________
1629 TLorentzVector AliMCAnalysisUtils::GetMotherWithPDG(Int_t label, Int_t pdg,
1630  const AliCaloTrackReader* reader,
1631  Bool_t & ok, Int_t & momlabel)
1632 {
1633  fGMotherMom.SetPxPyPzE(0,0,0,0);
1634 
1635  if(reader->ReadStack())
1636  {
1637  if(!reader->GetStack())
1638  {
1639  AliWarning("Stack is not available, check analysis settings in configuration file!!");
1640 
1641  ok = kFALSE;
1642  return fGMotherMom;
1643  }
1644 
1645  if(label >= 0 && label < reader->GetStack()->GetNtrack())
1646  {
1647  TParticle * momP = reader->GetStack()->Particle(label);
1648 
1649  if(momP->GetPdgCode()==pdg)
1650  {
1651  AliDebug(2,"PDG of mother is already the one requested!");
1652  fGMotherMom.SetPxPyPzE(momP->Px(),momP->Py(),momP->Pz(),momP->Energy());
1653  ok=kTRUE;
1654  return fGMotherMom;
1655  }
1656 
1657  Int_t grandmomLabel = momP->GetFirstMother();
1658  Int_t grandmomPDG = -1;
1659  TParticle * grandmomP = 0x0;
1660  while (grandmomLabel >=0 )
1661  {
1662  grandmomP = reader->GetStack()->Particle(grandmomLabel);
1663  grandmomPDG = grandmomP->GetPdgCode();
1664  if(grandmomPDG==pdg)
1665  {
1666  momlabel = grandmomLabel;
1667  fGMotherMom.SetPxPyPzE(grandmomP->Px(),grandmomP->Py(),grandmomP->Pz(),grandmomP->Energy());
1668  break;
1669  }
1670 
1671  grandmomLabel = grandmomP->GetFirstMother();
1672 
1673  }
1674 
1675  if(grandmomPDG!=pdg) AliInfo(Form("Mother with PDG %d, NOT found! \n",pdg));
1676  }
1677  }
1678  else if(reader->ReadAODMCParticles())
1679  {
1680  TClonesArray* mcparticles = reader->GetAODMCParticles();
1681  if(!mcparticles)
1682  {
1683  AliWarning("AODMCParticles is not available, check analysis settings in configuration file!!");
1684 
1685  ok=kFALSE;
1686  return fGMotherMom;
1687  }
1688 
1689  Int_t nprimaries = mcparticles->GetEntriesFast();
1690  if(label >= 0 && label < nprimaries)
1691  {
1692  AliAODMCParticle * momP = (AliAODMCParticle *) mcparticles->At(label);
1693 
1694  if(momP->GetPdgCode()==pdg)
1695  {
1696  AliDebug(2,"PDG of mother is already the one requested!");
1697  fGMotherMom.SetPxPyPzE(momP->Px(),momP->Py(),momP->Pz(),momP->E());
1698  ok=kTRUE;
1699  return fGMotherMom;
1700  }
1701 
1702  Int_t grandmomLabel = momP->GetMother();
1703  Int_t grandmomPDG = -1;
1704  AliAODMCParticle * grandmomP = 0x0;
1705  while (grandmomLabel >=0 )
1706  {
1707  grandmomP = (AliAODMCParticle *) mcparticles->At(grandmomLabel);
1708  grandmomPDG = grandmomP->GetPdgCode();
1709  if(grandmomPDG==pdg)
1710  {
1711  //printf("AliMCAnalysisUtils::GetMotherWithPDG(AOD) - mother with PDG %d FOUND! \n",pdg);
1712  momlabel = grandmomLabel;
1713  fGMotherMom.SetPxPyPzE(grandmomP->Px(),grandmomP->Py(),grandmomP->Pz(),grandmomP->E());
1714  break;
1715  }
1716 
1717  grandmomLabel = grandmomP->GetMother();
1718 
1719  }
1720 
1721  if(grandmomPDG!=pdg) AliInfo(Form("Mother with PDG %d, NOT found!",pdg));
1722 
1723  }
1724  }
1725 
1726  ok = kTRUE;
1727 
1728  return fGMotherMom;
1729 }
1730 
1731 //______________________________________________________________________________________________
1733 //______________________________________________________________________________________________
1734 TLorentzVector AliMCAnalysisUtils::GetGrandMother(Int_t label, const AliCaloTrackReader* reader,
1735  Int_t & pdg, Int_t & status, Bool_t & ok,
1736  Int_t & grandMomLabel, Int_t & greatMomLabel)
1737 {
1738  fGMotherMom.SetPxPyPzE(0,0,0,0);
1739 
1740  if(reader->ReadStack())
1741  {
1742  if(!reader->GetStack())
1743  {
1744  AliWarning("Stack is not available, check analysis settings in configuration file, STOP!!");
1745 
1746  ok = kFALSE;
1747  return fGMotherMom;
1748  }
1749 
1750  if(label >= 0 && label < reader->GetStack()->GetNtrack())
1751  {
1752  TParticle * momP = reader->GetStack()->Particle(label);
1753 
1754  grandMomLabel = momP->GetFirstMother();
1755 
1756  TParticle * grandmomP = 0x0;
1757 
1758  if (grandMomLabel >=0 )
1759  {
1760  grandmomP = reader->GetStack()->Particle(grandMomLabel);
1761  pdg = grandmomP->GetPdgCode();
1762  status = grandmomP->GetStatusCode();
1763 
1764  fGMotherMom.SetPxPyPzE(grandmomP->Px(),grandmomP->Py(),grandmomP->Pz(),grandmomP->Energy());
1765  greatMomLabel = grandmomP->GetFirstMother();
1766 
1767  }
1768  }
1769  }
1770  else if(reader->ReadAODMCParticles())
1771  {
1772  TClonesArray* mcparticles = reader->GetAODMCParticles();
1773  if(!mcparticles)
1774  {
1775  AliWarning("AODMCParticles is not available, check analysis settings in configuration file!!");
1776 
1777  ok=kFALSE;
1778  return fGMotherMom;
1779  }
1780 
1781  Int_t nprimaries = mcparticles->GetEntriesFast();
1782  if(label >= 0 && label < nprimaries)
1783  {
1784  AliAODMCParticle * momP = (AliAODMCParticle *) mcparticles->At(label);
1785 
1786  grandMomLabel = momP->GetMother();
1787 
1788  AliAODMCParticle * grandmomP = 0x0;
1789 
1790  if(grandMomLabel >=0 )
1791  {
1792  grandmomP = (AliAODMCParticle *) mcparticles->At(grandMomLabel);
1793  pdg = grandmomP->GetPdgCode();
1794  status = grandmomP->GetStatus();
1795 
1796  fGMotherMom.SetPxPyPzE(grandmomP->Px(),grandmomP->Py(),grandmomP->Pz(),grandmomP->E());
1797  greatMomLabel = grandmomP->GetMother();
1798 
1799  }
1800  }
1801  }
1802 
1803  ok = kTRUE;
1804 
1805  return fGMotherMom;
1806 }
1807 
1808 //_______________________________________________________________________________________________________________
1810 //_______________________________________________________________________________________________________________
1812  Float_t & asym, Float_t & angle, Bool_t & ok)
1813 {
1814  if(reader->ReadStack())
1815  {
1816  if(!reader->GetStack())
1817  {
1818  AliWarning("Stack is not available, check analysis settings in configuration file, STOP!!");
1819 
1820  ok = kFALSE;
1821  }
1822  if(label >= 0 && label < reader->GetStack()->GetNtrack())
1823  {
1824  TParticle * momP = reader->GetStack()->Particle(label);
1825 
1826  Int_t grandmomLabel = momP->GetFirstMother();
1827  Int_t grandmomPDG = -1;
1828  TParticle * grandmomP = 0x0;
1829  while (grandmomLabel >=0 )
1830  {
1831  grandmomP = reader->GetStack()->Particle(grandmomLabel);
1832  grandmomPDG = grandmomP->GetPdgCode();
1833 
1834  if(grandmomPDG==pdg) break;
1835 
1836  grandmomLabel = grandmomP->GetFirstMother();
1837 
1838  }
1839 
1840  if(grandmomPDG==pdg && grandmomP->GetNDaughters()==2)
1841  {
1842  TParticle * d1 = reader->GetStack()->Particle(grandmomP->GetDaughter(0));
1843  TParticle * d2 = reader->GetStack()->Particle(grandmomP->GetDaughter(1));
1844  if(d1->GetPdgCode() == 22 && d1->GetPdgCode() == 22)
1845  {
1846  asym = (d1->Energy()-d2->Energy())/grandmomP->Energy();
1847  d1->Momentum(fDaughMom );
1848  d2->Momentum(fDaughMom2);
1849  angle = fDaughMom.Angle(fDaughMom2.Vect());
1850  }
1851  }
1852  else
1853  {
1854  ok=kFALSE;
1855  AliInfo(Form("Mother with PDG %d, not found!",pdg));
1856  }
1857 
1858  } // good label
1859  }
1860  else if(reader->ReadAODMCParticles())
1861  {
1862  TClonesArray* mcparticles = reader->GetAODMCParticles();
1863  if(!mcparticles)
1864  {
1865  AliWarning("AODMCParticles is not available, check analysis settings in configuration file!!");
1866 
1867  ok=kFALSE;
1868  return;
1869  }
1870 
1871  Int_t nprimaries = mcparticles->GetEntriesFast();
1872  if(label >= 0 && label < nprimaries)
1873  {
1874  AliAODMCParticle * momP = (AliAODMCParticle *) mcparticles->At(label);
1875 
1876  Int_t grandmomLabel = momP->GetMother();
1877  Int_t grandmomPDG = -1;
1878  AliAODMCParticle * grandmomP = 0x0;
1879  while (grandmomLabel >=0 )
1880  {
1881  grandmomP = (AliAODMCParticle *) mcparticles->At(grandmomLabel);
1882  grandmomPDG = grandmomP->GetPdgCode();
1883 
1884  if(grandmomPDG==pdg) break;
1885 
1886  grandmomLabel = grandmomP->GetMother();
1887 
1888  }
1889 
1890  if(grandmomPDG==pdg && grandmomP->GetNDaughters()==2)
1891  {
1892  AliAODMCParticle * d1 = (AliAODMCParticle *) mcparticles->At(grandmomP->GetDaughter(0));
1893  AliAODMCParticle * d2 = (AliAODMCParticle *) mcparticles->At(grandmomP->GetDaughter(1));
1894  if(d1->GetPdgCode() == 22 && d1->GetPdgCode() == 22)
1895  {
1896  asym = (d1->E()-d2->E())/grandmomP->E();
1897  fDaughMom .SetPxPyPzE(d1->Px(),d1->Py(),d1->Pz(),d1->E());
1898  fDaughMom2.SetPxPyPzE(d2->Px(),d2->Py(),d2->Pz(),d2->E());
1899  angle = fDaughMom.Angle(fDaughMom2.Vect());
1900  }
1901  }
1902  else
1903  {
1904  ok=kFALSE;
1905  AliInfo(Form("Mother with PDG %d, not found! \n",pdg));
1906  }
1907 
1908  } // good label
1909  }
1910 
1911  ok = kTRUE;
1912 
1913 }
1914 
1915 //_________________________________________________________________________________________________
1917 //_________________________________________________________________________________________________
1918 Int_t AliMCAnalysisUtils::GetNDaughters(Int_t label, const AliCaloTrackReader* reader, Bool_t & ok)
1919 {
1920  if(reader->ReadStack())
1921  {
1922  if(!reader->GetStack())
1923  {
1924  AliWarning("Stack is not available, check analysis settings in configuration file, STOP!!");
1925 
1926  ok=kFALSE;
1927  return -1;
1928  }
1929  if(label >= 0 && label < reader->GetStack()->GetNtrack())
1930  {
1931  TParticle * momP = reader->GetStack()->Particle(label);
1932  ok=kTRUE;
1933  return momP->GetNDaughters();
1934  }
1935  else
1936  {
1937  ok = kFALSE;
1938  return -1;
1939  }
1940  }
1941  else if(reader->ReadAODMCParticles())
1942  {
1943  TClonesArray* mcparticles = reader->GetAODMCParticles();
1944  if(!mcparticles)
1945  {
1946  AliWarning("AODMCParticles is not available, check analysis settings in configuration file!!");
1947 
1948  ok=kFALSE;
1949  return -1;
1950  }
1951 
1952  Int_t nprimaries = mcparticles->GetEntriesFast();
1953  if(label >= 0 && label < nprimaries)
1954  {
1955  AliAODMCParticle * momP = (AliAODMCParticle *) mcparticles->At(label);
1956  ok = kTRUE;
1957  return momP->GetNDaughters();
1958  }
1959  else
1960  {
1961  ok = kFALSE;
1962  return -1;
1963  }
1964  }
1965 
1966  ok = kFALSE;
1967 
1968  return -1;
1969 }
1970 
1971 //_________________________________________________________________________________
1976 //_________________________________________________________________________________
1977 Int_t AliMCAnalysisUtils::GetNOverlaps(const Int_t * label, UInt_t nlabels,
1978  Int_t mctag, Int_t mesonLabel,
1979  AliCaloTrackReader * reader, Int_t *overpdg)
1980 {
1981  Int_t ancPDG = 0, ancStatus = -1;
1982  TVector3 prodVertex;
1983  Int_t ancLabel = 0;
1984  Int_t noverlaps = 0;
1985  Bool_t ok = kFALSE;
1986 
1987  for (UInt_t ilab = 1; ilab < nlabels; ilab++ )
1988  {
1989  ancLabel = CheckCommonAncestor(label[0],label[ilab],reader,ancPDG,ancStatus,fMotherMom,prodVertex);
1990 
1991  //printf("Overlaps, i %d: Main Label %d, second label %d, ancestor: Label %d, pdg %d - tag %d \n",
1992  // ilab,label[0],label[ilab],ancLabel,ancPDG, mctag);
1993 
1994  Bool_t overlap = kFALSE;
1995 
1996  if ( ancLabel < 0 )
1997  {
1998  overlap = kTRUE;
1999  //printf("\t \t \t No Label = %d\n",ancLabel);
2000  }
2001  else if( ( ancPDG==111 || ancPDG==221 ) && ( CheckTagBit(mctag,kMCPi0) || CheckTagBit(mctag,kMCEta)) && mesonLabel != ancLabel)
2002  {
2003  //printf("\t \t meson Label %d, ancestor Label %d\n",mesonLabel,ancLabel);
2004  overlap = kTRUE;
2005  }
2006  else if( ancPDG!=22 && TMath::Abs(ancPDG)!=11 && ancPDG != 111 && ancPDG != 221 )
2007  {
2008  //printf("\t \t \t Non EM PDG = %d\n",ancPDG);
2009  overlap = kTRUE ;
2010  }
2011 
2012  if( !overlap ) continue ;
2013 
2014  // We have at least one overlap
2015 
2016  //printf("Overlap!!!!!!!!!!!!!!\n");
2017 
2018  noverlaps++;
2019 
2020  // What is the origin of the overlap?
2021  Bool_t mOK = 0, gOK = 0;
2022  Int_t mpdg = -999999, gpdg = -1;
2023  Int_t mstatus = -1, gstatus = -1;
2024  Int_t gLabel = -1, ggLabel = -1;
2025 
2026  GetMother (label[ilab],reader,mpdg,mstatus,mOK);
2027  fGMotherMom =
2028  GetGrandMother(label[ilab],reader,gpdg,gstatus,gOK, gLabel,ggLabel);
2029 
2030  //printf("\t Overlap!, mother pdg %d; grand mother pdg %d",mpdg,gpdg);
2031 
2032  if( ( mpdg == 22 || TMath::Abs(mpdg==11) ) &&
2033  ( gpdg == 22 || TMath::Abs(gpdg==11) ) &&
2034  gLabel >=0 )
2035  {
2036  Int_t labeltmp = gLabel;
2037  while( ( gpdg == 22 || TMath::Abs(gpdg==11) ) && gLabel >=0 )
2038  {
2039  mpdg=gpdg;
2040  fGMotherMom = GetGrandMother(labeltmp,reader,gpdg,gstatus,ok, gLabel,ggLabel);
2041  labeltmp=gLabel;
2042  }
2043  }
2044  overpdg[noverlaps-1] = mpdg;
2045  }
2046 
2047  return noverlaps ;
2048 }
2049 
2050 //________________________________________________________
2052 //________________________________________________________
2053 void AliMCAnalysisUtils::Print(const Option_t * opt) const
2054 {
2055  if(! opt)
2056  return;
2057 
2058  printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
2059 
2060  printf("Debug level = %d\n",fDebug);
2061  printf("MC Generator = %s\n",fMCGeneratorString.Data());
2062  printf(" \n");
2063 }
2064 
2065 //__________________________________________________
2067 //__________________________________________________
2068 void AliMCAnalysisUtils::PrintMCTag(Int_t tag) const
2069 {
2070  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",
2071  tag,
2072  CheckTagBit(tag,kMCPhoton),
2074  CheckTagBit(tag,kMCPrompt),
2076  CheckTagBit(tag,kMCISR),
2077  CheckTagBit(tag,kMCPi0Decay),
2078  CheckTagBit(tag,kMCEtaDecay),
2080  CheckTagBit(tag,kMCPi0),
2081  CheckTagBit(tag,kMCEta),
2082  CheckTagBit(tag,kMCElectron),
2083  CheckTagBit(tag,kMCMuon),
2084  CheckTagBit(tag,kMCPion),
2085  CheckTagBit(tag,kMCProton),
2087  CheckTagBit(tag,kMCKaon),
2088  CheckTagBit(tag,kMCAntiProton),
2090  CheckTagBit(tag,kMCUnknown),
2092  );
2093 }
2094 
2095 //__________________________________________________
2097 //__________________________________________________
2099 {
2100  fMCGenerator = mcgen ;
2101  if (mcgen == kPythia) fMCGeneratorString = "PYTHIA";
2102  else if(mcgen == kHerwig) fMCGeneratorString = "HERWIG";
2103  else if(mcgen == kHijing) fMCGeneratorString = "HIJING";
2104  else
2105  {
2106  fMCGeneratorString = "";
2108  }
2109 }
2110 
2111 //____________________________________________________
2113 //____________________________________________________
2115 {
2116  fMCGeneratorString = mcgen ;
2117 
2118  if (mcgen == "PYTHIA") fMCGenerator = kPythia;
2119  else if(mcgen == "HERWIG") fMCGenerator = kHerwig;
2120  else if(mcgen == "HIJING") fMCGenerator = kHijing;
2121  else
2122  {
2124  fMCGeneratorString = "" ;
2125  }
2126 }
2127 
2128 
2129 
TLorentzVector GetMotherWithPDG(Int_t label, Int_t pdg, const AliCaloTrackReader *reader, Bool_t &ok, Int_t &momLabel)
Int_t pdg
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)
ClassImp(AliAnalysisTaskCRC) AliAnalysisTaskCRC
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