AliPhysics  master (68a99cd)
AliAODConversionPhoton.cxx
Go to the documentation of this file.
1 #include "AliAODMCParticle.h"
2 #include "AliAODMCHeader.h"
5 
6 using namespace std;
7 
9 
13  fLeadingNeutralPionIndex(-1),
14  fLeadingNeutralPionDaughterIndex(-1),
15  fCaloClusterRef(-1),
16  fDCArPrimVtx(0),
17  fDCAzPrimVtx(0),
18  fInvMassPair(0),
19  fNCaloPhotonMCLabels(0),
20  fNCaloPhotonMotherMCLabels(0),
21  fNNeutralPionLabels(0),
22  fCaloPhotonMCFlags(0),
23  fPairedId(-1),
24  fCaloPhoton(kFALSE),
25  fUseForMesonPair(kTRUE)
26 {
27  // initialize calo photon MC labels
28  for (Int_t i =0; i<50; i++){
29  fCaloPhotonMCLabels[i]=-1;
30  }
31  for (Int_t i =0; i<20; i++){
32  fCaloPhotonMotherMCLabels[i]=-1;
33  fNeutralPionLabels[i]=-1;
34  fNeutralPionEnergyFraction[i]=0;
35  }
36  //Standard constructor
37 }
38 
40 AliAODConversionParticle(kfphoton),
42  fLeadingNeutralPionIndex(-1),
43  fLeadingNeutralPionDaughterIndex(-1),
44  fCaloClusterRef(-1),
45  fDCArPrimVtx(0),
46  fDCAzPrimVtx(0),
47  fInvMassPair(0),
48  fNCaloPhotonMCLabels(0),
49  fNCaloPhotonMotherMCLabels(0),
50  fNNeutralPionLabels(0),
51  fCaloPhotonMCFlags(0),
52  fPairedId(-1),
53  fCaloPhoton(kFALSE),
54  fUseForMesonPair(kTRUE)
55 {
56  //Constructor from kfphoton
57  // puts the mass to zero and store dilepton mass
58  SetMass(kfphoton->M());
59 
60  //SetE(P());
61  // initialize calo photon MC labels
62  for (Int_t i =0; i<50; i++){
63  fCaloPhotonMCLabels[i]=-1;
64  }
65  for (Int_t i =0; i<20; i++){
67  fNeutralPionLabels[i]=-1;
69  }
70 
71 }
72 
78  fCaloClusterRef(-1),
79  fDCArPrimVtx(0),
80  fDCAzPrimVtx(0),
81  fInvMassPair(0),
86  fPairedId(-1),
87  fCaloPhoton(kFALSE),
88  fUseForMesonPair(kTRUE)
89 {
90  //Constructor from TLorentzVector
91 
92  // initialize calo photon MC labels
93  for (Int_t i =0; i<50; i++){
94  fCaloPhotonMCLabels[i]=-1;
95  }
96  for (Int_t i =0; i<20; i++){
98  fNeutralPionLabels[i]=-1;
100  }
101 }
102 
103 
104 
106 AliAODConversionParticle(original),
107 AliConversionPhotonBase(original),
111  fDCArPrimVtx(original.fDCArPrimVtx),
112  fDCAzPrimVtx(original.fDCAzPrimVtx),
113  fInvMassPair(original.fInvMassPair),
118  fPairedId(original.fPairedId),
119  fCaloPhoton(original.fCaloPhoton),
121 {
122  //Copy constructor
123 
124  // initialize calo photon MC labels
125  for (Int_t i =0; i<50; i++){
127  }
128  for (Int_t i =0; i<20; i++){
130  fNeutralPionLabels[i]=original.fNeutralPionLabels[i];
132  }
133 }
134 
136 {
137  // empty standard destructor
138 }
139 
141 {
142  // assignment operator
143  return *this;
144 }
145 
148 
149  Double_t primCo[3] = {primVertex->GetX(),primVertex->GetY(),primVertex->GetZ()};
150 
151  Double_t absoluteP = TMath::Sqrt(TMath::Power(GetPx(),2) + TMath::Power(GetPy(),2) + TMath::Power(GetPz(),2));
152  Double_t p[3] = {GetPx()/absoluteP,GetPy()/absoluteP,GetPz()/absoluteP};
153  Double_t CP[3];
154 
155  CP[0] = fConversionPoint[0] - primCo[0];
156  CP[1] = fConversionPoint[1] - primCo[1];
157  CP[2] = fConversionPoint[2] - primCo[2];
158 
159  Double_t Lambda = - (CP[0]*p[0]+CP[1]*p[1]+CP[2]*p[2])/(p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);
160 
161  Double_t S[3];
162  S[0] = fConversionPoint[0] + p[0]*Lambda;
163  S[1] = fConversionPoint[1] + p[1]*Lambda;
164  S[2] = fConversionPoint[2] + p[2]*Lambda;
165 
166  fDCArPrimVtx = TMath::Sqrt( TMath::Power(primCo[0]-S[0],2) + TMath::Power(primCo[1]-S[1],2));
167  fDCAzPrimVtx = primCo[2]-S[2];
168 
169  return;
170 }
171 
172 
173 void AliAODConversionPhoton::SetCaloPhotonMCFlags(AliMCEvent *mcEvent, Bool_t enableSort, Bool_t mergedAnalysis, AliVCluster* cluster){
174 
175 
176  TParticle* PhotonDummyMerged;
177  TParticle* PhotonDummyMergedMother;
178  TParticle* PhotonDummyMergedGrandMother;
179  TParticle* PhotonDummyMergedGGMother;
180  Int_t photonDummyMergedPDG= -1;
181  Int_t photonDummyMergedMotherPDG= -1;
182  Int_t photonDummyMergedGrandMotherPDG= -1;
183  Int_t photonDummyMergedGGMotherPDG= -1;
184  Bool_t foundNeutralPion = kFALSE;
185  Int_t neutralPionLabel = -1;
186  if(mergedAnalysis){
187  for (Int_t j = 0; j< fNCaloPhotonMCLabels; j++){
188  neutralPionLabel = -1;
189  foundNeutralPion = kFALSE;
190  PhotonDummyMerged = mcEvent->Particle(GetCaloPhotonMCLabel(j)); // main particle
191  photonDummyMergedPDG = PhotonDummyMerged->GetPdgCode();
192  if(TMath::Abs(photonDummyMergedPDG)==111){
193  foundNeutralPion = kTRUE;
194  neutralPionLabel = GetCaloPhotonMCLabel(j);
195  }
196  if(PhotonDummyMerged->GetMother(0)>-1 && !foundNeutralPion){
197  PhotonDummyMergedMother = mcEvent->Particle(PhotonDummyMerged->GetMother(0)); // mother
198  photonDummyMergedMotherPDG = PhotonDummyMergedMother->GetPdgCode();
199  if(TMath::Abs(photonDummyMergedMotherPDG)==111){
200  foundNeutralPion = kTRUE;
201  neutralPionLabel = PhotonDummyMerged->GetMother(0);
202  }
203  if(PhotonDummyMergedMother->GetMother(0)>-1 && !foundNeutralPion){
204  PhotonDummyMergedGrandMother = mcEvent->Particle(PhotonDummyMergedMother->GetMother(0)); // grandmother
205  photonDummyMergedGrandMotherPDG = PhotonDummyMergedGrandMother->GetPdgCode();
206  if(TMath::Abs(photonDummyMergedGrandMotherPDG)==111){
207  foundNeutralPion = kTRUE;
208  neutralPionLabel = PhotonDummyMergedMother->GetMother(0);
209  }
210  if(PhotonDummyMergedGrandMother->GetMother(0)>-1 && !foundNeutralPion){
211  PhotonDummyMergedGGMother = mcEvent->Particle(PhotonDummyMergedGrandMother->GetMother(0)); // grand-grandmother
212  photonDummyMergedGGMotherPDG = PhotonDummyMergedGGMother->GetPdgCode();
213  if(TMath::Abs(photonDummyMergedGGMotherPDG)==111){
214  foundNeutralPion = kTRUE;
215  neutralPionLabel = PhotonDummyMergedGrandMother->GetMother(0);
216  }
217  }
218  }
219  }
220  // check if current found pi0 was already found before in the cluster
221  // if it was found before, add the cluster energy fraction of the current label to this pi0
222  if(foundNeutralPion){
223  Bool_t newNeutralPion = true;
224  for(Int_t k=0; k<fNNeutralPionLabels; k++){
225  if (fNeutralPionLabels[k] == neutralPionLabel){ //check if mother is already contained in fNeutralPionLabels
226  newNeutralPion = false;
227  fNeutralPionEnergyFraction[k] += cluster->GetClusterMCEdepFraction(j);
228  }
229  }
230  // if the pi0 is new, then increase pi0 count by one and set E-fraction and label to array
231  if (newNeutralPion){
232  fNeutralPionLabels[fNNeutralPionLabels] = neutralPionLabel;
233  fNeutralPionEnergyFraction[fNNeutralPionLabels] = cluster->GetClusterMCEdepFraction(j);
234  fNNeutralPionLabels++; //only if particle label is not yet contained in array, count up fNeutralPionLabels
235  }
236  }
237  }
238 
239  // search for the largest contributing pi0 in the cluster
240  Double_t largestEfraction = 0;
241  for(Int_t k=0; k<fNNeutralPionLabels; k++){
242  if(fNeutralPionEnergyFraction[k]>largestEfraction){
243  largestEfraction = fNeutralPionEnergyFraction[k];
245  }
246  }
247  // if the cluster contains a pi0 contribution, find the MC label/index of the largest contributing daughter
248  // for this, go upwards three mothers from each MC label to find leading daughter index
250  for(Int_t l=0; l<fNCaloPhotonMCLabels; l++){
252  continue;
253  PhotonDummyMerged = mcEvent->Particle(GetCaloPhotonMCLabel(l)); // main particle
254  if(PhotonDummyMerged->GetMother(0)>-1){
255  if(PhotonDummyMerged->GetMother(0)==fNeutralPionLabels[fLeadingNeutralPionIndex]){
257  }else{
258  PhotonDummyMergedMother = mcEvent->Particle(PhotonDummyMerged->GetMother(0));
259  if(PhotonDummyMergedMother->GetMother(0)>-1){
260  if(PhotonDummyMergedMother->GetMother(0)==fNeutralPionLabels[fLeadingNeutralPionIndex]){
262  }else{
263  PhotonDummyMergedGrandMother = mcEvent->Particle(PhotonDummyMergedMother->GetMother(0));
264  if(PhotonDummyMergedGrandMother->GetMother(0)>-1){
265  if(PhotonDummyMergedGrandMother->GetMother(0)==fNeutralPionLabels[fLeadingNeutralPionIndex]){
267  }else{
268  PhotonDummyMergedGGMother = mcEvent->Particle(PhotonDummyMergedGrandMother->GetMother(0));
269  if(PhotonDummyMergedGGMother->GetMother(0)>-1){
270  if(PhotonDummyMergedGGMother->GetMother(0)==fNeutralPionLabels[fLeadingNeutralPionIndex]){
272  }
273  }
274  }
275  }
276  }
277  }
278  }
279  }
280  }
281  }
282  }
283 
284  Bool_t isPhoton = kFALSE; // largest contribution to cluster is photon
285  Bool_t isElectron = kFALSE; // largest contribution to cluster is electron
286  Bool_t isConversion = kFALSE; // largest contribution to cluster is converted electron
287  Bool_t isConversionFullyContained = kFALSE; // largest contribution to cluster is converted electron, second electron has been found in same cluster
288  Bool_t isMerged = kFALSE; // cluster contains more than one particle from the same decay
289  Bool_t isMergedPartConv = kFALSE; // cluster contains more than one particle from the same decay and at least one of the particles came from a conversion
290  Bool_t isDalitz = kFALSE; // this cluster was created by a particle stemming from a dalitz decay
291  Bool_t isDalitzMerged = kFALSE; // this cluster was created by a particle stemming from a dalitz decay and more than one particle of the dalitz decay is contained in the cluster
292  Bool_t isPhotonWithElecMother = kFALSE; // this cluster is from a photon with an electron as mother
293  Bool_t isShower = kFALSE; // this cluster contains as a largest contribution a particle from a shower or radiative process
294  Bool_t isSubLeadingEM = kFALSE; // cluster contains at least one electron or photon from a pi0, eta or eta_prime in subleading contribution
295  Bool_t isElectronFromFragPhoton = kFALSE; // largest contribution to cluster is from converted electron, but photon stems from fragmentation photon ( q -> q gamma)
296 
297  TParticle* Photon = 0x0;
298  TParticle* PhotonMother = 0x0;
299  TParticle* PhotonGrandMother = 0x0;
300  if (fNCaloPhotonMCLabels==0) return;
301 
302  if (enableSort){
303  // sort the array according to the energy of contributing particles
304  if (fNCaloPhotonMCLabels>1){
305  // cout << "start sorting" << endl;
306  Int_t* sortIdx = new Int_t[fNCaloPhotonMCLabels];
307  Double_t* energyPerPart = new Double_t[fNCaloPhotonMCLabels];
308  Long_t* orginalContrib = new Long_t[fNCaloPhotonMCLabels];
309  for (Int_t i = 0; i < fNCaloPhotonMCLabels; i++){
310  orginalContrib[i] = fCaloPhotonMCLabels[i];
311  if (fCaloPhotonMCLabels[i]> -1){
312  TParticle* dummy = mcEvent->Particle(fCaloPhotonMCLabels[i]);
313  if (dummy){
314  energyPerPart[i] = dummy->Energy();
315  // suppress energy of hadrons !!! DIRTY hack !!!
316  if (!(TMath::Abs(dummy->GetPdgCode())== 11 || TMath::Abs(dummy->GetPdgCode())== 22)){
317  // cout << "suppressed hadron energy for:" << dummy->GetPdgCode() << endl;
318  energyPerPart[i]= 0.25; // put energy to mip
319  }
320  }
321  } else {
322  energyPerPart[i] = 0;
323  }
324  }
325 
326  TMath::Sort(fNCaloPhotonMCLabels,energyPerPart,sortIdx);
327  for(Int_t index = 0; index < fNCaloPhotonMCLabels; index++) {
328  fCaloPhotonMCLabels[index] = orginalContrib [sortIdx[index]] ;
329 
330  }
331  delete [] sortIdx;
332  delete [] energyPerPart;
333  delete [] orginalContrib;
334  }
335  }
336 
337  if(mergedAnalysis && fNNeutralPionLabels>0 && fLeadingNeutralPionDaughterIndex!=0 && fNeutralPionEnergyFraction[fLeadingNeutralPionIndex]>cluster->GetClusterMCEdepFraction(0)){
338  // check if leading pi0 comes not from label 0 in cluster
339  // for this do:
340  // -> check if neutral pions were found in cluster
341  // -> if the leading daughter index is not 0
342  // -> the leading neutral pion has a larger cluster energy fraction than the cluster label 0
343  // load particle corresponding to largest daughter of leading pi0
344  Photon = mcEvent->Particle(GetCaloPhotonMCLabel(fLeadingNeutralPionDaughterIndex));
345  } else {
346  // load particle corresponding to MC label 0 in cluster
347  if(GetCaloPhotonMCLabel(0)>-1)
348  Photon = mcEvent->Particle(GetCaloPhotonMCLabel(0));
349  }
350 
351 
352  if(Photon == NULL){
353  return;
354  }
355 
356  Int_t particleMotherLabel = Photon->GetMother(0);
357  Int_t particleGrandMotherLabel = -1;
358  Int_t particleGrandMotherX2Label = -1;
359  Int_t particleGrandMotherX3Label = -1;
360  Int_t particleGrandMotherX4Label = -1;
361  Int_t particleGrandMotherX5Label = -1;
362  Int_t particleGrandMotherX6Label = -1;
363  Int_t particleGrandMotherX7Label = -1;
364 
365  Int_t particleMotherPDG = -1;
366  Int_t particleGrandMotherPDG = -1;
367  Int_t particleGrandMotherX2PDG = -1;
368  Int_t particleGrandMotherX3PDG = -1;
369  Int_t particleGrandMotherX4PDG = -1;
370  Int_t particleGrandMotherX5PDG = -1;
371  Int_t particleGrandMotherX6PDG = -1;
372 
373  Int_t particleMotherNDaugthers = 0;
374  Int_t particleGrandMotherNDaugthers = 0;
375 
376  if (particleMotherLabel > -1){
377  PhotonMother = mcEvent->Particle(particleMotherLabel);
378  particleMotherNDaugthers = PhotonMother->GetNDaughters();
379  particleGrandMotherLabel = PhotonMother->GetMother(0);
380  particleMotherPDG = PhotonMother->GetPdgCode();
381  if (particleGrandMotherLabel > -1){
382  PhotonGrandMother = mcEvent->Particle(particleGrandMotherLabel);
383  particleGrandMotherPDG = PhotonGrandMother->GetPdgCode();
384  particleGrandMotherNDaugthers = PhotonGrandMother->GetNDaughters();
385  particleGrandMotherX2Label = PhotonGrandMother->GetMother(0);
386  TParticle* dummyGMM = 0x0;
387  if (particleGrandMotherX2Label > -1){
388  dummyGMM = mcEvent->Particle(particleGrandMotherX2Label);
389  particleGrandMotherX2PDG = dummyGMM->GetPdgCode();
390  particleGrandMotherX3Label = dummyGMM->GetMother(0);
391  if (particleGrandMotherX3Label > -1){
392  dummyGMM = mcEvent->Particle(particleGrandMotherX3Label);
393  particleGrandMotherX3PDG = dummyGMM->GetPdgCode();
394  particleGrandMotherX4Label = dummyGMM->GetMother(0);
395  if (particleGrandMotherX4Label > -1){
396  dummyGMM = mcEvent->Particle(particleGrandMotherX4Label);
397  particleGrandMotherX4PDG = dummyGMM->GetPdgCode();
398  particleGrandMotherX5Label = dummyGMM->GetMother(0);
399  if (particleGrandMotherX5Label > -1){
400  dummyGMM = mcEvent->Particle(particleGrandMotherX5Label);
401  particleGrandMotherX5PDG = dummyGMM->GetPdgCode();
402  particleGrandMotherX6Label = dummyGMM->GetMother(0);
403  if (particleGrandMotherX6Label > -1){
404  dummyGMM = mcEvent->Particle(particleGrandMotherX6Label);
405  particleGrandMotherX6PDG = dummyGMM->GetPdgCode();
406  particleGrandMotherX7Label = dummyGMM->GetMother(0);
407  }
408  }
409  }
410  }
411  }
412  }
413  }
414 
415  //determine mother/grandmother of leading particle and if it is pion/eta/eta_prime: fill array fCaloPhotonMotherMCLabels at position 0
416  if (particleMotherLabel > -1){
417  if( TMath::Abs(particleMotherPDG) == 111 || TMath::Abs(particleMotherPDG) == 221 || TMath::Abs(particleMotherPDG) == 331 ){
418  fCaloPhotonMotherMCLabels[0] = particleMotherLabel;
420  } else if (particleGrandMotherLabel > -1 && TMath::Abs(particleMotherPDG) == 22 &&
421  (TMath::Abs(particleGrandMotherPDG) == 111 || TMath::Abs(particleGrandMotherPDG) == 221 || TMath::Abs(particleGrandMotherPDG) == 331) ){
422  fCaloPhotonMotherMCLabels[0] = particleGrandMotherLabel;
424  } else {
425  fCaloPhotonMotherMCLabels[0] = particleMotherLabel;
427  }
428  }
429 
430  // Check whether the first contribution was photon
431  if(TMath::Abs(Photon->GetPdgCode()) == 22){
432  isPhoton = kTRUE;
433  // did it decay via the dalitz channel
434  if (particleMotherLabel > -1 && particleMotherNDaugthers == 3)
435  isDalitz = kTRUE;
436  // Test whether particle stems from a shower or radiation
437  if (TMath::Abs(particleMotherPDG) == 11){ // check whether photon stems from electron
438  isPhotonWithElecMother = kTRUE;
439  if (particleGrandMotherLabel > -1){ // test whether first particle has a grandmother
440  if (TMath::Abs(particleGrandMotherPDG) == 22 )
441  isShower = kTRUE; // check whether grandmother is a photon (meaning this is most likely a shower)
442  }
443  }
444  }
445 
446  // Check whether the first contribution was electron
447  if( TMath::Abs(Photon->GetPdgCode()) == 11 ){
448  isElectron = kTRUE;
449  if (particleMotherLabel > -1) {
450  // was it a conversion
451  if (TMath::Abs(particleMotherPDG) == 22)
452  isConversion = kTRUE;
453  // did it decay via the dalitz channel
454  if (particleGrandMotherLabel > -1 && particleGrandMotherNDaugthers == 3 )
455  isDalitz = kTRUE;
456  }
457  if (particleGrandMotherLabel > -1){ // check whether electron has a grandmother
458  if (TMath::Abs(particleGrandMotherPDG) == 11 || TMath::Abs(particleGrandMotherPDG) == 22){ // test whether electron has photon or electron as grandmother (meaning will most likely be a shower)
459  isShower = kTRUE;
460  }
461  }
462  // consider the rare case, where the conversion electron stems from photon, which stems from electron, which stems from photon, ...
463  if (isConversion && TMath::Abs(particleGrandMotherPDG) == 11 && particleGrandMotherX2PDG == 22){
464  SetCaloPhotonMCLabel(0,particleGrandMotherLabel);
465  fCaloPhotonMotherMCLabels[0] = particleGrandMotherX3Label;
466  if (TMath::Abs(particleGrandMotherX3PDG) == 11 && particleGrandMotherX4PDG == 22 ){
467  SetCaloPhotonMCLabel(0,particleGrandMotherX3Label);
468  fCaloPhotonMotherMCLabels[0] = particleGrandMotherX5Label;
469  if (TMath::Abs(particleGrandMotherX5PDG) == 11 && particleGrandMotherX6PDG == 22 ){
470  SetCaloPhotonMCLabel(0,particleGrandMotherX5Label);
471  fCaloPhotonMotherMCLabels[0] = particleGrandMotherX7Label;
472  }
473  }
474  }
475 
476  // consider the case, where a photon stems from a photon which stems from a photon etc...which is common for frag. photons
477  if (particleMotherLabel > -1){
478  TParticle *dummyMother = mcEvent->Particle(particleMotherLabel);
479  Bool_t originReached = kFALSE;
480  if (dummyMother){
481  while (dummyMother->GetPdgCode() == 22 && !originReached){ // follow conversion photon's history, as long as the mother is a photon
482  if (dummyMother->GetMother(0) > -1){
483  dummyMother = mcEvent->Particle(dummyMother->GetMother(0));
484  if ((TMath::Abs(dummyMother->GetPdgCode()) == 11) || (TMath::Abs(dummyMother->GetPdgCode()) == 22)){ // in case of additional conversion skip to photon's grandma, which should be a photon
485  if (dummyMother->GetMother(0) > -1){ // also mother of photon could be a photon (can happen for fragmentation)
486  dummyMother = mcEvent->Particle(dummyMother->GetMother(0));
487  } else {
488  originReached = kTRUE;
489  }
490  } else {
491  originReached = kTRUE;
492  }
493  isElectronFromFragPhoton = (TMath::Abs(dummyMother->GetPdgCode()) < 6);// photon stems from quark = fragmentation photon
494  } else {
495  originReached = kTRUE;
496  }
497  }
498  }
499  }
500  }
501 
502  Bool_t enablePrintOuts = kFALSE;
503  // if (fNCaloPhotonMCLabels>2)
504  // enablePrintOuts = kTRUE;
505 
506  // check whether there were other contributions to the cluster
507  if (fNCaloPhotonMCLabels>1){
508  TParticle* dummyPart =NULL;
509  TParticle* dummyPartMother =NULL;
510  TParticle* dummyPartGrandMother =NULL;
511 
512  for (Int_t i = 1; i< fNCaloPhotonMCLabels; i++){
513  if (i > 49) continue; // abort if more than 50 entries to the cluster have been checked (more are not stored in these objects)
514  if (enablePrintOuts) cout << "checking particle: " << i << endl;
515  if (GetCaloPhotonMCLabel(i) < 0) continue;
516  dummyPart = mcEvent->Particle(GetCaloPhotonMCLabel(i));
517  Int_t dummyPartMotherLabel = dummyPart->GetMother(0);
518  Int_t dummyPartGrandMotherLabel = -1;
519  Int_t dummyPartMotherPDG = -1;
520  Int_t dummyPartGrandMotherPDG = -1;
521  // check whether this particle has a mother & obtain the pdg code
522  if (dummyPartMotherLabel > -1){
523  dummyPartMother = mcEvent->Particle(dummyPartMotherLabel);
524  dummyPartGrandMotherLabel = dummyPartMother->GetMother(0);
525  dummyPartMotherPDG = dummyPartMother->GetPdgCode();
526  // check whether this particle has a grandmother & obtain its pdg code
527  if (dummyPartGrandMotherLabel > -1){
528  dummyPartGrandMother = mcEvent->Particle(dummyPartGrandMotherLabel);
529  dummyPartGrandMotherPDG = dummyPartGrandMother->GetPdgCode();
530  }
531  }
532  // largest contribution was from photon and is not from shower or electron mother
533  if (isPhoton && (!isShower || !isPhotonWithElecMother )){
534  if (enablePrintOuts) cout << "lead gamma" << endl;
535  if (dummyPartMotherLabel > -1 && particleMotherLabel > -1){ // test whether first particle has a mother
536  if (dummyPartMotherLabel == particleMotherLabel)
537  isMerged = kTRUE; // test whether current and first particle have the same mother => i.e. other gamma from decay or dalitz electron
538  if (dummyPartGrandMotherLabel > -1){ // test whether first particle has a grandmother
539  // check whether particle is an electron from a conversion of a photon from the original mother
540  if (TMath::Abs(dummyPart->GetPdgCode()) == 11 && dummyPartGrandMotherLabel == particleMotherLabel )
541  isMergedPartConv = kTRUE;
542  // check whether particle is an electron from a dalitz decay from the original mother
543  if (TMath::Abs(dummyPart->GetPdgCode()) == 11 && dummyPartMotherLabel == particleMotherLabel )
544  isDalitzMerged = kTRUE;
545  }
546  }
547  }
548 
549  // largest contribution was from electron & not a from a shower
550  if (isElectron && !isShower){
551  if (enablePrintOuts) cout << "lead electron" << endl;
552  if (dummyPartMotherLabel > -1 && particleMotherLabel > -1){ // test whether first particle has a mother
553  if ( TMath::Abs(dummyPart->GetPdgCode()) == 11 && isConversion && dummyPartMotherLabel == particleMotherLabel ) {
554  isConversionFullyContained = kTRUE; // test whether conversion is fully contained in cluster
555  if (enablePrintOuts) cout << "found full conversion" << endl;
556  } else if (TMath::Abs(dummyPart->GetPdgCode()) == 11) { // current particle is an electron
557  if (enablePrintOuts) cout << "current is electron" << endl;
558  // check whether current particle is a shower
559  Bool_t isShowerCurrent = kFALSE;
560  if (TMath::Abs(dummyPartGrandMotherPDG) == 11)
561  isShowerCurrent = kTRUE;
562  if (TMath::Abs(dummyPartGrandMotherPDG) == 22)
563  isShowerCurrent = kTRUE;
564 
565  // check whether current particle is a conversion
566  Bool_t isConvCurrent = kFALSE;
567  if (TMath::Abs(dummyPartMotherPDG) == 22 && !isShowerCurrent)
568  isConvCurrent = kTRUE;
569 
570  if (enablePrintOuts) cout << "conv current: " << isConvCurrent << "\t shower current: " << isShowerCurrent << endl;
571  // check whether orginal electron and this electron stem from the same particle and electron originated in dalitz decay
572  if( dummyPartMotherLabel == particleMotherLabel && TMath::Abs(particleMotherPDG) != 22 && !isShowerCurrent ) {
573  isDalitzMerged = kTRUE;
574  isMerged = kTRUE;
575  // both particles are from conversions
576  } else if (isConversion && isConvCurrent){
577  if (particleGrandMotherLabel > -1 && dummyPartGrandMotherLabel > -1){ // test whether first particle has a grandmother and current particle has one as well
578  // check whether orginal electron and this electron stem from the same particle and electron stems from conversion
579  if( dummyPartGrandMotherLabel == particleGrandMotherLabel && !isShowerCurrent )
580  isMergedPartConv = kTRUE;
581  }
582  // check whether current electron is from conversion & orginal e are from same mother
583  } else if (dummyPartGrandMotherLabel > -1 && isConvCurrent){
584  if (dummyPartGrandMotherLabel == particleMotherLabel && (TMath::Abs(particleMotherPDG) == 111 || TMath::Abs(particleMotherPDG) == 221 || TMath::Abs(particleMotherPDG) == 331) ){
585  isDalitzMerged = kTRUE;
586  isMerged = kTRUE;
587  }
588  // check whether current electron & orginal e is from conversion are from same mother
589  } else if (isConversion && particleGrandMotherLabel > -1){
590  if (dummyPartMotherLabel == particleGrandMotherLabel && (TMath::Abs(dummyPartMotherPDG) == 111 || TMath::Abs(dummyPartMotherPDG) == 221 || TMath::Abs(dummyPartMotherPDG) == 331) ){
591  isDalitzMerged = kTRUE;
592  isMerged = kTRUE;
593  }
594  }
595  } else if (TMath::Abs(dummyPart->GetPdgCode()) == 22){ // test whether this particle is a photon
596  // check whether current particle is a shower
597  Bool_t isShowerCurrent = kFALSE;
598  if (TMath::Abs(dummyPartMotherPDG) == 11)
599  isShowerCurrent = kTRUE;
600  if (TMath::Abs(dummyPartMotherPDG) == 22)
601  isShowerCurrent = kTRUE;
602 
603  // check whether orginal electron and this photon stem from the same particle and electron originated in dalitz
604  if( dummyPartMotherLabel == particleMotherLabel && !isShowerCurrent ) {
605  isDalitzMerged = kTRUE;
606  isMerged = kTRUE;
607  } else if (particleGrandMotherLabel > -1){ // test whether first particle has a grandmother
608  // check whether orginal electron and this photon stem from the same particle and electron stems from conversion
609  if( dummyPartMotherLabel == particleGrandMotherLabel && !isShowerCurrent )
610  isMergedPartConv = kTRUE;
611  }
612  }
613  }
614  }
615 
616 
617  if (dummyPartMotherLabel > -1){ // test whether particle has a mother
618  if (TMath::Abs(dummyPart->GetPdgCode()) == 22){ // test whether particle is a photon
619  //check if photon directly comes from a pion/eta/eta_prime decay
620  Bool_t helpN = true;
621  for(Int_t j=0; j<fNCaloPhotonMotherMCLabels; j++){
622  if (fCaloPhotonMotherMCLabels[j] == dummyPartMotherLabel){ //check if mother is already contained in fCaloPhotonMotherMCLabels
623  helpN = false;
624  }
625  }
626  if (helpN){
628  fNCaloPhotonMotherMCLabels++; //only if particle label is not yet contained in array, count up fNCaloPhotonMotherMCLabels
629  }
630  if ((TMath::Abs(dummyPartMotherPDG) == 111 || TMath::Abs(dummyPartMotherPDG) == 221 || TMath::Abs(dummyPartMotherPDG) == 331) && !isPhoton && !isElectron)
631  isSubLeadingEM = kTRUE;
632  } else if (TMath::Abs(dummyPart->GetPdgCode()) == 11){ //test whether particle is an electron
633  //check if electron comes from a pion decay
634  if ( TMath::Abs(dummyPartMotherPDG) != 22 ){
635  Bool_t helpN = true;
636  for(Int_t j=0; j<fNCaloPhotonMotherMCLabels; j++){
637  if (fCaloPhotonMotherMCLabels[j] == dummyPartMotherLabel){ //check if mother is already contained in fCaloPhotonMotherMCLabels
638  helpN = false;
639  }
640  }
641  if (helpN){
643  fNCaloPhotonMotherMCLabels++; //only if particle label is not yet contained in array, count up fNCaloPhotonMotherMCLabels
644  }
645  if ((TMath::Abs(dummyPartMotherPDG) == 111 || TMath::Abs(dummyPartMotherPDG) == 221 || TMath::Abs(dummyPartMotherPDG) == 331) && !isPhoton && !isElectron)
646  isSubLeadingEM = kTRUE;
647  } else if (dummyPartGrandMotherLabel > -1){ //if it is not a dalitz decay, test whether particle has a grandmother
648  //check if it is a conversion electron that has pion/eta/eta_prime as grandmother
649  if ( (TMath::Abs(dummyPartGrandMotherPDG) == 111 || TMath::Abs(dummyPartGrandMotherPDG) == 221 || TMath::Abs(dummyPartGrandMotherPDG) == 331)){
650 
651  Bool_t helpN = true;
652  for(Int_t j=0; j<fNCaloPhotonMotherMCLabels; j++){
653  if (fCaloPhotonMotherMCLabels[j] == dummyPartGrandMotherLabel){ //check if grandmother is already contained in fCaloPhotonMotherMCLabels
654  helpN = false;
655  }
656  }
657  if (helpN){
658  fCaloPhotonMotherMCLabels[fNCaloPhotonMotherMCLabels] = dummyPartGrandMotherLabel;
659  fNCaloPhotonMotherMCLabels++; //only if particle label is not yet contained in array, count up fNCaloPhotonMotherMCLabels
660  }
661  if (!isPhoton && !isElectron)
662  isSubLeadingEM = kTRUE;
663  }
664  }
665  }
666  }
667  }
668  }
669  fCaloPhotonMCFlags = isPhoton *1 + isElectron *2 + isConversion*4+ isConversionFullyContained *8 + isMerged *16 + isMergedPartConv*32 + isDalitz *64 + isDalitzMerged *128 + isPhotonWithElecMother *256 + isShower * 512 + isSubLeadingEM * 1024 + isElectronFromFragPhoton * 2048;
670 }
671 
672 void AliAODConversionPhoton::SetCaloPhotonMCFlagsAOD(AliVEvent* event, Bool_t enableSort, Bool_t mergedAnalysis, AliVCluster* cluster){
673 
674  TClonesArray *AODMCTrackArray = dynamic_cast<TClonesArray*>(event->FindListObject(AliAODMCParticle::StdBranchName()));
675  if (!AODMCTrackArray) return;
676 
677  AliAODMCParticle* PhotonDummyMerged;
678  AliAODMCParticle* PhotonDummyMergedMother;
679  AliAODMCParticle* PhotonDummyMergedGrandMother;
680  AliAODMCParticle* PhotonDummyMergedGGMother;
681  Int_t photonDummyMergedPDG= -1;
682  Int_t photonDummyMergedMotherPDG= -1;
683  Int_t photonDummyMergedGrandMotherPDG= -1;
684  Int_t photonDummyMergedGGMotherPDG= -1;
685  Bool_t foundNeutralPion = kFALSE;
686  Int_t neutralPionLabel = -1;
687  if(mergedAnalysis){
688  // this part searches the MC labels of the cluster/photon candidate for all contributions from pi0s and save all pi0s that contributed
689  for (Int_t j = 0; j< fNCaloPhotonMCLabels; j++){
690  neutralPionLabel = -1;
691  foundNeutralPion = kFALSE;
692  PhotonDummyMerged = (AliAODMCParticle*) AODMCTrackArray->At(GetCaloPhotonMCLabel(j)); // main particle
693  photonDummyMergedPDG = PhotonDummyMerged->GetPdgCode();
694  if(TMath::Abs(photonDummyMergedPDG)==111){
695  foundNeutralPion = kTRUE;
696  neutralPionLabel = GetCaloPhotonMCLabel(j);
697  }
698  if(PhotonDummyMerged->GetMother()>-1 && !foundNeutralPion){
699  PhotonDummyMergedMother = (AliAODMCParticle*) AODMCTrackArray->At(PhotonDummyMerged->GetMother()); // mother
700  photonDummyMergedMotherPDG = PhotonDummyMergedMother->GetPdgCode();
701  if(TMath::Abs(photonDummyMergedMotherPDG)==111){
702  foundNeutralPion = kTRUE;
703  neutralPionLabel = PhotonDummyMerged->GetMother();
704  }
705  if(PhotonDummyMergedMother->GetMother()>-1 && !foundNeutralPion){
706  PhotonDummyMergedGrandMother = (AliAODMCParticle*) AODMCTrackArray->At(PhotonDummyMergedMother->GetMother()); // grandmother
707  photonDummyMergedGrandMotherPDG = PhotonDummyMergedGrandMother->GetPdgCode();
708  if(TMath::Abs(photonDummyMergedGrandMotherPDG)==111){
709  foundNeutralPion = kTRUE;
710  neutralPionLabel = PhotonDummyMergedMother->GetMother();
711  }
712  if(PhotonDummyMergedGrandMother->GetMother()>-1 && !foundNeutralPion){
713  PhotonDummyMergedGGMother = (AliAODMCParticle*) AODMCTrackArray->At(PhotonDummyMergedGrandMother->GetMother()); // grand-grandmother
714  photonDummyMergedGGMotherPDG = PhotonDummyMergedGGMother->GetPdgCode();
715  if(TMath::Abs(photonDummyMergedGGMotherPDG)==111){
716  foundNeutralPion = kTRUE;
717  neutralPionLabel = PhotonDummyMergedGrandMother->GetMother();
718  }
719  }
720  }
721  }
722  // check if current found pi0 was already found before in the cluster
723  // if it was found before, add the cluster energy fraction of the current label to this pi0
724  if(foundNeutralPion){
725  Bool_t newNeutralPion = true;
726  for(Int_t k=0; k<fNNeutralPionLabels; k++){
727  if (fNeutralPionLabels[k] == neutralPionLabel){ //check if mother is already contained in fNeutralPionLabels
728  newNeutralPion = false;
729  fNeutralPionEnergyFraction[k] += cluster->GetClusterMCEdepFraction(j);
730  }
731  }
732  // if the pi0 is new, then increase pi0 count by one and set E-fraction and label to array
733  if (newNeutralPion){
734  fNeutralPionLabels[fNNeutralPionLabels] = neutralPionLabel;
735  fNeutralPionEnergyFraction[fNNeutralPionLabels] = cluster->GetClusterMCEdepFraction(j);
736  fNNeutralPionLabels++; //only if particle label is not yet contained in array, count up fNeutralPionLabels
737  }
738  }
739  }
740 
741  // search for the largest contributing pi0 in the cluster
742  Double_t largestEfraction = 0;
743  for(Int_t k=0; k<fNNeutralPionLabels; k++){
744  if(fNeutralPionEnergyFraction[k]>largestEfraction){
745  largestEfraction = fNeutralPionEnergyFraction[k];
747  }
748  }
749  // if the cluster contains a pi0 contribution, find the MC label/index of the largest contributing daughter
750  // for this, go upwards three mothers from each MC label to find leading daughter index
752  for(Int_t l=0; l<fNCaloPhotonMCLabels; l++){
754  continue;
755  PhotonDummyMerged = (AliAODMCParticle*) AODMCTrackArray->At(GetCaloPhotonMCLabel(l)); // main particle
756  if(PhotonDummyMerged->GetMother()>-1){
757  if(PhotonDummyMerged->GetMother()==fNeutralPionLabels[fLeadingNeutralPionIndex]){
759  }else{
760  PhotonDummyMergedMother = (AliAODMCParticle*) AODMCTrackArray->At(PhotonDummyMerged->GetMother());
761  if(PhotonDummyMergedMother->GetMother()>-1){
762  if(PhotonDummyMergedMother->GetMother()==fNeutralPionLabels[fLeadingNeutralPionIndex]){
764  }else{
765  PhotonDummyMergedGrandMother = (AliAODMCParticle*) AODMCTrackArray->At(PhotonDummyMergedMother->GetMother());
766  if(PhotonDummyMergedGrandMother->GetMother()>-1){
767  if(PhotonDummyMergedGrandMother->GetMother()==fNeutralPionLabels[fLeadingNeutralPionIndex]){
769  }else{
770  PhotonDummyMergedGGMother = (AliAODMCParticle*) AODMCTrackArray->At(PhotonDummyMergedGrandMother->GetMother());
771  if(PhotonDummyMergedGGMother->GetMother()>-1){
772  if(PhotonDummyMergedGGMother->GetMother()==fNeutralPionLabels[fLeadingNeutralPionIndex]){
774  }
775  }
776  }
777  }
778  }
779  }
780  }
781  }
782  }
783  }
784  }
785 
786  Bool_t isPhoton = kFALSE; // largest contribution to cluster is photon
787  Bool_t isElectron = kFALSE; // largest contribution to cluster is electron
788  Bool_t isConversion = kFALSE; // largest contribution to cluster is converted electron
789  Bool_t isConversionFullyContained = kFALSE; // largest contribution to cluster is converted electron, second electron has been found in same cluster
790  Bool_t isMerged = kFALSE; // largest contribution to cluster is photon, second photon or electron from dalitz has been found in same cluster
791  Bool_t isMergedPartConv = kFALSE; // cluster contains more than one particle from the same decay and at least one of the particles came from a conversion
792  Bool_t isDalitz = kFALSE; // this cluster was created by a particle stemming from a dality decay
793  Bool_t isDalitzMerged = kFALSE; // this cluster was created by a particle stemming from a dality decay and more than one particle of the dalitz decay is contained in the cluster
794  Bool_t isPhotonWithElecMother = kFALSE; // this cluster is from a photon with an electron as mother
795  Bool_t isShower = kFALSE; // this cluster contains as a largest contribution a particle from a shower or radiative process
796  Bool_t isSubLeadingEM = kFALSE; // cluster contains at least one electron or photon from a pi0 or eta in subleading contribution
797  Bool_t isElectronFromFragPhoton = kFALSE; // largest contribution to cluster is from converted electron, but photon stems from fragmentation photon ( q -> q gamma)
798 
799  AliAODMCParticle* Photon = NULL;
800  AliAODMCParticle* PhotonMother = NULL;
801  AliAODMCParticle* PhotonGrandMother = NULL;
802 
803  if (fNCaloPhotonMCLabels==0) return;
804 
805  if (enableSort){
806  // sort the array according to the energy of contributing particles
807  if (fNCaloPhotonMCLabels>1){
808  // cout << "start sorting" << endl;
809  Int_t* sortIdx = new Int_t[fNCaloPhotonMCLabels];
810  Double_t* energyPerPart = new Double_t[fNCaloPhotonMCLabels];
811  Long_t* orginalContrib = new Long_t[fNCaloPhotonMCLabels];
812  for (Int_t i = 0; i < fNCaloPhotonMCLabels; i++){
813  orginalContrib[i] = fCaloPhotonMCLabels[i];
814  if (fCaloPhotonMCLabels[i]> -1){
815  AliAODMCParticle* dummy = (AliAODMCParticle*) AODMCTrackArray->At(GetCaloPhotonMCLabel(i));
816  energyPerPart[i] = dummy->E();
817  // suppress energy of hadrons !!! DIRTY hack !!!
818  if (!(TMath::Abs(dummy->GetPdgCode())== 11 || TMath::Abs(dummy->GetPdgCode())== 22)){
819  // cout << "suppressed hadron energy" << endl;
820  energyPerPart[i]= 0.25; // put energy to mip
821  }
822  } else {
823  energyPerPart[i] = 0;
824  }
825  }
826 
827  TMath::Sort(fNCaloPhotonMCLabels,energyPerPart,sortIdx);
828  for(Int_t index = 0; index < fNCaloPhotonMCLabels; index++) {
829  fCaloPhotonMCLabels[index] = orginalContrib [sortIdx[index]] ;
830 
831  }
832  delete [] sortIdx;
833  delete [] energyPerPart;
834  delete [] orginalContrib;
835  }
836  }
837 
838  if(mergedAnalysis && fNNeutralPionLabels>0 && fLeadingNeutralPionDaughterIndex!=0 && fNeutralPionEnergyFraction[fLeadingNeutralPionIndex]>cluster->GetClusterMCEdepFraction(0)){
839  // check if leading pi0 comes not from label 0 in cluster
840  // for this do:
841  // -> check if neutral pions were found in cluster
842  // -> if the leading daughter index is not 0
843  // -> the leading neutral pion has a larger cluster energy fraction than the cluster label 0
844  // load particle corresponding to largest daughter of leading pi0
845  Photon = (AliAODMCParticle*) AODMCTrackArray->At(GetCaloPhotonMCLabel(fLeadingNeutralPionDaughterIndex));
846  } else {
847  // load particle corresponding to MC label 0 in cluster
848  if(GetCaloPhotonMCLabel(0)>-1){
849  Photon = (AliAODMCParticle*) AODMCTrackArray->At(GetCaloPhotonMCLabel(0));
850  }
851  }
852 
853 
854 
855  if(Photon == NULL){
856  return;
857  }
858 
859  Int_t particleMotherLabel = Photon->GetMother();
860  Int_t particleGrandMotherLabel = -1;
861  Int_t particleMotherPDG = -1;
862  Int_t particleGrandMotherPDG = -1;
863  Int_t particleMotherNDaugthers = 0;
864  Int_t particleGrandMotherNDaugthers = 0;
865  Int_t particleGrandMotherX2Label = -1;
866  Int_t particleGrandMotherX3Label = -1;
867  Int_t particleGrandMotherX4Label = -1;
868  Int_t particleGrandMotherX5Label = -1;
869  Int_t particleGrandMotherX6Label = -1;
870  Int_t particleGrandMotherX7Label = -1;
871  Int_t particleGrandMotherX2PDG = -1;
872  Int_t particleGrandMotherX3PDG = -1;
873  Int_t particleGrandMotherX4PDG = -1;
874  Int_t particleGrandMotherX5PDG = -1;
875  Int_t particleGrandMotherX6PDG = -1;
876 
877  if (particleMotherLabel > -1){
878  PhotonMother = (AliAODMCParticle*) AODMCTrackArray->At(Photon->GetMother());
879  particleMotherNDaugthers = PhotonMother->GetNDaughters();
880  particleGrandMotherLabel = PhotonMother->GetMother();
881  particleMotherPDG = PhotonMother->GetPdgCode();
882  if (particleGrandMotherLabel > -1){
883  PhotonGrandMother = (AliAODMCParticle*) AODMCTrackArray->At(PhotonMother->GetMother());
884  particleGrandMotherPDG = PhotonGrandMother->GetPdgCode();
885  particleGrandMotherNDaugthers = PhotonGrandMother->GetNDaughters();
886  AliAODMCParticle* dummyGMM = NULL;
887 
888  particleGrandMotherX2Label = PhotonGrandMother->GetMother();
889  if (particleGrandMotherX2Label > -1){
890  dummyGMM = (AliAODMCParticle*) AODMCTrackArray->At(particleGrandMotherX2Label);
891  particleGrandMotherX2PDG = dummyGMM->GetPdgCode();
892  particleGrandMotherX3Label = dummyGMM->GetMother();
893  if (particleGrandMotherX3Label > -1){
894  dummyGMM = (AliAODMCParticle*) AODMCTrackArray->At(particleGrandMotherX3Label);
895  particleGrandMotherX3PDG = dummyGMM->GetPdgCode();
896  particleGrandMotherX4Label = dummyGMM->GetMother();
897  if (particleGrandMotherX4Label > -1){
898  dummyGMM = (AliAODMCParticle*) AODMCTrackArray->At(particleGrandMotherX4Label);
899  particleGrandMotherX4PDG = dummyGMM->GetPdgCode();
900  particleGrandMotherX5Label = dummyGMM->GetMother();
901  if (particleGrandMotherX5Label > -1){
902  dummyGMM = (AliAODMCParticle*) AODMCTrackArray->At(particleGrandMotherX5Label);
903  particleGrandMotherX5PDG = dummyGMM->GetPdgCode();
904  particleGrandMotherX6Label = dummyGMM->GetMother();
905  if (particleGrandMotherX6Label > -1){
906  dummyGMM = (AliAODMCParticle*) AODMCTrackArray->At(particleGrandMotherX6Label);
907  particleGrandMotherX6PDG = dummyGMM->GetPdgCode();
908  particleGrandMotherX7Label = dummyGMM->GetMother();
909  }
910  }
911  }
912  }
913  }
914  }
915  }
916 
917  //determine mother/grandmother of leading particle and if it is pion/eta/eta_prime: fill array fCaloPhotonMotherMCLabels at position 0
918  if (particleMotherLabel > -1){
919  if( TMath::Abs(particleMotherPDG) == 111 || TMath::Abs(particleMotherPDG) == 221 || TMath::Abs(particleMotherPDG) == 331 ){
920  fCaloPhotonMotherMCLabels[0] = particleMotherLabel;
922  } else if (TMath::Abs(particleMotherPDG) == 22 && particleGrandMotherLabel > -1){
923  if ( TMath::Abs(particleMotherPDG) == 22 && (TMath::Abs(particleGrandMotherPDG) == 111 || TMath::Abs(particleGrandMotherPDG) == 221 || TMath::Abs(particleGrandMotherPDG) == 331) ){
924  fCaloPhotonMotherMCLabels[0] = particleGrandMotherLabel;
926  }
927  } else {
928  fCaloPhotonMotherMCLabels[0] = particleMotherLabel;
930  }
931  }
932 
933  // Check whether the first contribution was photon
934  if(TMath::Abs(Photon->GetPdgCode()) == 22){
935  isPhoton = kTRUE;
936  // did it decay via the dalitz channel
937  if (particleMotherLabel > -1 && particleMotherNDaugthers == 3)
938  isDalitz = kTRUE;
939  // Test whether particle stems from a shower or radiation
940  if (TMath::Abs(particleMotherPDG) == 11){ // check whether photon stems from electron
941  isPhotonWithElecMother = kTRUE;
942  if (particleGrandMotherLabel > -1){ // test whether first particle has a grandmother
943  if (TMath::Abs(particleGrandMotherPDG) == 22 )
944  isShower = kTRUE; // check whether grandmother is a photon (meaning this is most likely a shower)
945  }
946  }
947  }
948  // Check whether the first contribution was electron
949  if(TMath::Abs(Photon->GetPdgCode()) == 11 ){
950  isElectron = kTRUE;
951  if (particleMotherLabel > -1) {
952  // was it a conversion
953  if (TMath::Abs(particleMotherPDG) == 22)
954  isConversion = kTRUE;
955  // did it decay via the dalitz channel
956  if (particleGrandMotherLabel > -1 && particleGrandMotherNDaugthers == 3 )
957  isDalitz = kTRUE;
958  }
959  if (particleGrandMotherLabel > -1){ // check whether electron has a grandmother
960  if (TMath::Abs(particleGrandMotherPDG) == 11 || TMath::Abs(particleGrandMotherPDG) == 22){ // test whether electron has photon or electron as grandmother (meaning will most likely be a shower)
961  isShower = kTRUE;
962  }
963  }
964  // consider the rare case, where the conversion electron stems from photon, which stems from electron, which stems from photon, ...
965  if (isConversion && TMath::Abs(particleGrandMotherPDG) == 11 && particleGrandMotherX2PDG == 22){
966  SetCaloPhotonMCLabel(0,particleGrandMotherLabel);
967  fCaloPhotonMotherMCLabels[0] = particleGrandMotherX3Label;
968  if (TMath::Abs(particleGrandMotherX3PDG) == 11 && particleGrandMotherX4PDG == 22 ){
969  SetCaloPhotonMCLabel(0,particleGrandMotherX3Label);
970  fCaloPhotonMotherMCLabels[0] = particleGrandMotherX5Label;
971  if (TMath::Abs(particleGrandMotherX5PDG) == 11 && particleGrandMotherX6PDG == 22 ){
972  SetCaloPhotonMCLabel(0,particleGrandMotherX5Label);
973  fCaloPhotonMotherMCLabels[0] = particleGrandMotherX7Label;
974  }
975  }
976  }
977 
978  // consider the case, where a photon stems from a photon which stems from a photon etc...which is common for frag. photons
979  if (particleMotherLabel > -1){
980  AliAODMCParticle* dummyMother = (AliAODMCParticle*) AODMCTrackArray->At(Photon->GetMother());
981  if (dummyMother){
982  Bool_t originReached = kFALSE;
983  while (dummyMother->GetPdgCode() == 22 && !originReached){ // follow conversion photon's history, as long as the mother is a photon
984  if (dummyMother->GetMother() > -1){
985  dummyMother = (AliAODMCParticle*) AODMCTrackArray->At(dummyMother->GetMother());
986  if ((TMath::Abs(dummyMother->GetPdgCode()) == 11) || (TMath::Abs(dummyMother->GetPdgCode()) == 22)){ // in case of additional conversion skip to photon's grandma, which should be a photon
987  if (dummyMother->GetMother() > -1){
988  dummyMother = (AliAODMCParticle*) AODMCTrackArray->At(dummyMother->GetMother());
989  } else {
990  originReached = kTRUE;
991  }
992  } else {
993  originReached = kTRUE;
994  }
995  isElectronFromFragPhoton = (TMath::Abs(dummyMother->GetPdgCode()) < 6);// photon stems from quark = fragmentation photon
996  } else {
997  originReached = kTRUE;
998  }
999  }
1000  }
1001  }
1002  }
1003 
1004 
1005 
1006  // check whether there were other contributions to the cluster
1007  if (fNCaloPhotonMCLabels>1){
1008  AliAODMCParticle* dummyPart = NULL;
1009  AliAODMCParticle* dummyPartMother = NULL;
1010  AliAODMCParticle* dummyPartGrandMother = NULL;
1011  for (Int_t i = 1; i< fNCaloPhotonMCLabels; i++){
1012  if (i > 49) continue; // abort if more than 50 entries to the cluster have been checked (more are not stored in these objects)
1013  dummyPart = (AliAODMCParticle*) AODMCTrackArray->At(GetCaloPhotonMCLabel(i));
1014  Int_t dummyPartMotherLabel = dummyPart->GetMother();
1015  Int_t dummyPartGrandMotherLabel = -1;
1016  Int_t dummyPartMotherPDG = -1;
1017  Int_t dummyPartGrandMotherPDG = -1;
1018 
1019  // check whether this particle has a mother & obtain the pdg code
1020  if (dummyPartMotherLabel > -1){
1021  dummyPartMother = (AliAODMCParticle*) AODMCTrackArray->At(dummyPart->GetMother());
1022  dummyPartGrandMotherLabel = dummyPartMother->GetMother();
1023  dummyPartMotherPDG = dummyPartMother->GetPdgCode();
1024  // check whether this particle has a grandmother & obtain its pdg code
1025  if (dummyPartGrandMotherLabel > -1){
1026  dummyPartGrandMother = (AliAODMCParticle*) AODMCTrackArray->At(dummyPartMother->GetMother());
1027  dummyPartGrandMotherPDG = dummyPartGrandMother->GetPdgCode();
1028  }
1029  }
1030 
1031 
1032  // largest contribution was from photon and is not from shower or electron mother
1033  if (isPhoton && (!isShower || !isPhotonWithElecMother )){
1034  if (dummyPartMotherLabel > -1 && particleMotherLabel > -1){ // test whether first particle has a mother
1035  if (dummyPartMotherLabel == particleMotherLabel)
1036  isMerged = kTRUE; // test whether current and first particle have the same mother => i.e. other gamma from decay or dalitz electron
1037  if (dummyPartGrandMotherLabel > -1){ // test whether first particle has a grandmother
1038  // check whether particle is an electron from a conversion of a photon from the original mother
1039  if (TMath::Abs(dummyPart->GetPdgCode()) == 11 && dummyPartGrandMotherLabel == particleMotherLabel )
1040  isMergedPartConv = kTRUE;
1041  // check whether particle is an electron from a dalitz decay from the original mother
1042  if (TMath::Abs(dummyPart->GetPdgCode()) == 11 && dummyPartMotherLabel == particleMotherLabel )
1043  isDalitzMerged = kTRUE;
1044  }
1045  }
1046  }
1047 
1048  // largest contribution was from electron & not a from a shower
1049  if (isElectron && !isShower){
1050  if (dummyPartMotherLabel > -1 && particleMotherLabel > -1){ // test whether first particle has a mother
1051  if ( TMath::Abs(dummyPart->GetPdgCode()) == 11 && isConversion && dummyPartMotherLabel == particleMotherLabel ) {
1052  isConversionFullyContained = kTRUE; // test whether conversion is fully contained in cluster
1053  } else if (TMath::Abs(dummyPart->GetPdgCode()) == 11) { // current particle is an electron
1054 
1055  // check whether current particle is a shower
1056  Bool_t isShowerCurrent = kFALSE;
1057  if (TMath::Abs(dummyPartGrandMotherPDG) == 11)
1058  isShowerCurrent = kTRUE;
1059  if (TMath::Abs(dummyPartGrandMotherPDG) == 22)
1060  isShowerCurrent = kTRUE;
1061 
1062  // check whether current particle is a conversion
1063  Bool_t isConvCurrent = kFALSE;
1064  if (TMath::Abs(dummyPartMotherPDG) == 22 && !isShowerCurrent)
1065  isConvCurrent = kTRUE;
1066 
1067  // check whether orginal electron and this electron stem from the same particle and electron originated in dalitz decay
1068  if( dummyPartMotherLabel == particleMotherLabel && TMath::Abs(particleMotherPDG) != 22 && !isShowerCurrent ) {
1069  isDalitzMerged = kTRUE;
1070  isMerged = kTRUE;
1071  // both particles are from conversions
1072  } else if (isConversion && isConvCurrent){
1073  if (particleGrandMotherLabel > -1 && dummyPartGrandMotherLabel > -1){ // test whether first particle has a grandmother and current particle has one as well
1074  // check whether orginal electron and this electron stem from the same particle and electron stems from conversion
1075  if( dummyPartGrandMotherLabel == particleGrandMotherLabel && !isShowerCurrent )
1076  isMergedPartConv = kTRUE;
1077  }
1078  // check whether current electron is from conversion & orginal e are from same mother
1079  } else if (dummyPartGrandMotherLabel > -1 && isConvCurrent){
1080  if (dummyPartGrandMotherLabel == particleMotherLabel && (TMath::Abs(particleMotherPDG) == 111 || TMath::Abs(particleMotherPDG) == 221 || TMath::Abs(particleMotherPDG) == 331) ){
1081  isDalitzMerged = kTRUE;
1082  isMerged = kTRUE;
1083  }
1084  // check whether current electron & orginal e is from conversion are from same mother
1085  } else if (isConversion && particleGrandMotherLabel > -1){
1086  if (dummyPartMotherLabel == particleGrandMotherLabel && (TMath::Abs(dummyPartMotherPDG) == 111 || TMath::Abs(dummyPartMotherPDG) == 221 || TMath::Abs(dummyPartMotherPDG) == 331) ){
1087  isDalitzMerged = kTRUE;
1088  isMerged = kTRUE;
1089  }
1090  }
1091  } else if (TMath::Abs(dummyPart->GetPdgCode()) == 22){ // test whether this particle is a photon
1092  // check whether current particle is a shower
1093  Bool_t isShowerCurrent = kFALSE;
1094  if (TMath::Abs(dummyPartMotherPDG) == 11)
1095  isShowerCurrent = kTRUE;
1096  if (TMath::Abs(dummyPartMotherPDG) == 22)
1097  isShowerCurrent = kTRUE;
1098 
1099  // check whether orginal electron and this photon stem from the same particle and electron originated in dalitz
1100  if( dummyPartMotherLabel == particleMotherLabel && !isShowerCurrent ) {
1101  isDalitzMerged = kTRUE;
1102  isMerged = kTRUE;
1103  } else if (particleGrandMotherLabel > -1){ // test whether first particle has a grandmother
1104  // check whether orginal electron and this photon stem from the same particle and electron stems from conversion
1105  if( dummyPartMotherLabel == particleGrandMotherLabel && !isShowerCurrent )
1106  isMergedPartConv = kTRUE;
1107  }
1108  }
1109  }
1110  }
1111 
1112  if (dummyPartMotherLabel > -1){ // test whether particle has a mother
1113  if (TMath::Abs(dummyPart->GetPdgCode()) == 22){ // test whether particle is a photon
1114  //check if photon directly comes from a pion/eta/eta_prime decay
1115  Bool_t helpN = true;
1116  for(Int_t j=0; j<fNCaloPhotonMotherMCLabels; j++){
1117  if (fCaloPhotonMotherMCLabels[j] == dummyPartMotherLabel){ //check if mother is already contained in fCaloPhotonMotherMCLabels
1118  helpN = false;
1119  }
1120  }
1121  if (helpN){
1122  fCaloPhotonMotherMCLabels[fNCaloPhotonMotherMCLabels] = dummyPartMotherLabel;
1123  fNCaloPhotonMotherMCLabels++; //only if particle label is not yet contained in array, count up fNCaloPhotonMotherMCLabels
1124  }
1125  if ((TMath::Abs(dummyPartMotherPDG) == 111 || TMath::Abs(dummyPartMotherPDG) == 221 || TMath::Abs(dummyPartMotherPDG) == 331) && !isPhoton && !isElectron)
1126  isSubLeadingEM = kTRUE;
1127  } else if (TMath::Abs(dummyPart->GetPdgCode()) == 11){ //test whether particle is an electron
1128  //check if electron comes from a pion decay
1129  if ( TMath::Abs(dummyPartMotherPDG) != 22 ){
1130  Bool_t helpN = true;
1131  for(Int_t j=0; j<fNCaloPhotonMotherMCLabels; j++){
1132  if (fCaloPhotonMotherMCLabels[j] == dummyPartMotherLabel){ //check if mother is already contained in fCaloPhotonMotherMCLabels
1133  helpN = false;
1134  }
1135  }
1136  if (helpN){
1137  fCaloPhotonMotherMCLabels[fNCaloPhotonMotherMCLabels] = dummyPartMotherLabel;
1138  fNCaloPhotonMotherMCLabels++; //only if particle label is not yet contained in array, count up fNCaloPhotonMotherMCLabels
1139  }
1140  if ((TMath::Abs(dummyPartMotherPDG) == 111 || TMath::Abs(dummyPartMotherPDG) == 221 || TMath::Abs(dummyPartMotherPDG) == 331) && !isPhoton && !isElectron)
1141  isSubLeadingEM = kTRUE;
1142  } else if (dummyPartGrandMotherLabel > -1){ //if it is not a dalitz decay, test whether particle has a grandmother
1143  //check if it is a conversion electron that has pion/eta/eta_prime as grandmother
1144  if ( (TMath::Abs(dummyPartGrandMotherPDG) == 111 || TMath::Abs(dummyPartGrandMotherPDG) == 221 || TMath::Abs(dummyPartGrandMotherPDG) == 331)){
1145 
1146  Bool_t helpN = true;
1147  for(Int_t j=0; j<fNCaloPhotonMotherMCLabels; j++){
1148  if (fCaloPhotonMotherMCLabels[j] == dummyPartGrandMotherLabel){ //check if grandmother is already contained in fCaloPhotonMotherMCLabels
1149  helpN = false;
1150  }
1151  }
1152  if (helpN){
1153  fCaloPhotonMotherMCLabels[fNCaloPhotonMotherMCLabels] = dummyPartGrandMotherLabel;
1154  fNCaloPhotonMotherMCLabels++; //only if particle label is not yet contained in array, count up fNCaloPhotonMotherMCLabels
1155  }
1156  if (!isPhoton && !isElectron)
1157  isSubLeadingEM = kTRUE;
1158  }
1159  }
1160  }
1161  }
1162  }
1163  }
1164  fCaloPhotonMCFlags = isPhoton *1 + isElectron *2 + isConversion*4+ isConversionFullyContained *8 + isMerged *16 + isMergedPartConv*32 + isDalitz *64 + isDalitzMerged *128 + isPhotonWithElecMother *256 + isShower * 512 + isSubLeadingEM * 1024 + isElectronFromFragPhoton * 2048;;
1165 }
1166 
1167 //_____________________________________________________________________________
1168 // prints information of MC particles contributing to cluster
1169 //_____________________________________________________________________________
1171  cout << endl << endl << "particles contributing: " << endl;
1172  for (Int_t i =0 ; i < fNCaloPhotonMCLabels; i++ ){
1173  TParticle *dummy = NULL;
1174  if (fCaloPhotonMCLabels[i]>0){
1175  dummy = (TParticle*)mcEvent->Particle(fCaloPhotonMCLabels[i]);
1176  cout << i << "\t"<< fCaloPhotonMCLabels[i] << "\t pdg code: " <<dummy->GetPdgCode() << "\t prod radius: "<< dummy->R() << "\t energy: " << dummy->Energy() ;
1177  if (dummy->GetMother(0) > -1){
1178  TParticle* dummyMother = (TParticle*)mcEvent->Particle(dummy->GetMother(0));
1179  cout << "\t mother part: " << dummy->GetMother(0) << "\t mother pdg code: " << dummyMother->GetPdgCode() << "\t energy: " << dummyMother->Energy() << endl;
1180  if (dummyMother->GetMother(0) > -1){
1181  TParticle* dummyGrandMother = (TParticle*)mcEvent->Particle(dummyMother->GetMother(0));
1182  cout << "\t grandmother part: " << dummyMother->GetMother(0) << "\t grandmother pdg code: " << dummyGrandMother->GetPdgCode() << "\t energy: " << dummyGrandMother->Energy() << endl;
1183  } else {
1184  cout << endl;
1185  }
1186  } else {
1187  cout << endl;
1188  }
1189  }
1190 // if (dummy) delete dummy;
1191  }
1192  cout << "mothers contributing: " << endl;
1193  for (Int_t i =0 ; i < fNCaloPhotonMotherMCLabels; i++ ){
1194  TParticle *dummy = NULL;
1195  if (fCaloPhotonMotherMCLabels[i]>0){
1196  dummy = (TParticle*)mcEvent->Particle(fCaloPhotonMotherMCLabels[i]);
1197  cout << i << "\t"<< fCaloPhotonMotherMCLabels[i] << "\t pdg code: " <<dummy->GetPdgCode() << "\t prod radius: "<< dummy->R() << "\t energy: " << dummy->Energy() << endl;
1198 
1199 
1200  }
1201 // if (dummy) delete dummy;
1202  }
1203 }
1204 
1205 //_____________________________________________________________________________
1206 // prints information of cluster as it is flagged
1207 //_____________________________________________________________________________
1209  cout << fCaloPhotonMCFlags << "\t photon: " << IsLargestComponentPhoton() << "\t electron: " << IsLargestComponentElectron() << "\t conversion: " << IsConversion() << "\t conversion full: "
1210  << IsConversionFullyContained() << "\t merged: " << IsMerged() << "\t merged p.conv: " << IsMergedPartConv() << "\t Dalitz: " << IsDalitz() << "\t Dalitz merged: " << IsDalitzMerged()
1211  << "\t rad: " << IsPhotonWithElecMother() << "\t shower: " << IsShower() << "\t no EM lead: " << IsEMNonLeading() << "\t EM sub lead: " << IsSubLeadingEM() << endl;
1212 }
virtual Double_t GetPy() const
double Double_t
Definition: External.C:58
virtual Double_t GetPz() const
void SetCaloPhotonMCLabel(Int_t i, Int_t labelCaloPhoton)
AliAODConversionPhoton & operator=(const AliAODConversionPhoton &g)
void PrintCaloMCLabelsAndInfo(AliMCEvent *mcEvent)
void SetCaloPhotonMCFlagsAOD(AliVEvent *event, Bool_t enableSort, Bool_t mergedAnalysis=kFALSE, AliVCluster *cluster=NULL)
int Int_t
Definition: External.C:63
virtual Double_t GetPx() const
bool Bool_t
Definition: External.C:53
void CalculateDistanceOfClossetApproachToPrimVtx(const AliVVertex *primVertex)
void SetCaloPhotonMCFlags(AliMCEvent *mcEvent, Bool_t enableSort, Bool_t mergedAnalysis=kFALSE, AliVCluster *cluster=NULL)