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