AliPhysics  5364b50 (5364b50)
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 fDCArPrimVtx(0),
14 fDCAzPrimVtx(0),
15 fInvMassPair(0),
16 fCaloPhoton(0),
17 fCaloClusterRef(-1),
18 fNCaloPhotonMCLabels(0),
19 fNCaloPhotonMotherMCLabels(0),
20 fCaloPhotonMCFlags(0)
21 {
22  // initialize calo photon MC labels
23  for (Int_t i =0; i<50; i++){
24  fCaloPhotonMCLabels[i]=-1;
25  }
26  for (Int_t i =0; i<20; i++){
27  fCaloPhotonMotherMCLabels[i]=-1;
28  }
29  //Standard constructor
30 }
31 
33 AliAODConversionParticle(kfphoton),
35 fDCArPrimVtx(0),
36 fDCAzPrimVtx(0),
37 fInvMassPair(0),
38 fCaloPhoton(0),
39 fCaloClusterRef(-1),
40 fNCaloPhotonMCLabels(0),
41 fNCaloPhotonMotherMCLabels(0),
42 fCaloPhotonMCFlags(0)
43 {
44  //Constructor from kfphoton
45 
46  // puts the mass to zero and store dilepton mass
47  SetMass(kfphoton->M());
48 
49  //SetE(P());
50 
51  // initialize calo photon MC labels
52  for (Int_t i =0; i<50; i++){
53  fCaloPhotonMCLabels[i]=-1;
54  }
55  for (Int_t i =0; i<20; i++){
57  }
58 
59 }
60 
64 fDCArPrimVtx(0),
65 fDCAzPrimVtx(0),
66 fInvMassPair(0),
67 fCaloPhoton(0),
68 fCaloClusterRef(-1),
72 {
73  //Constructor from TLorentzVector
74 
75  // initialize calo photon MC labels
76  for (Int_t i =0; i<50; i++){
77  fCaloPhotonMCLabels[i]=-1;
78  }
79  for (Int_t i =0; i<20; i++){
81  }
82 }
83 
84 
85 
87 AliAODConversionParticle(original),
88 AliConversionPhotonBase(original),
89 fDCArPrimVtx(original.fDCArPrimVtx),
90 fDCAzPrimVtx(original.fDCAzPrimVtx),
91 fInvMassPair(original.fInvMassPair),
92 fCaloPhoton(original.fCaloPhoton),
97 {
98  //Copy constructor
99 
100  // initialize calo photon MC labels
101  for (Int_t i =0; i<50; i++){
103  }
104  for (Int_t i =0; i<20; i++){
106  }
107 }
108 
110 {
111  // empty standard destructor
112 }
113 
115 {
116  // assignment operator
117  return *this;
118 }
119 
122 
123  Double_t primCo[3] = {primVertex->GetX(),primVertex->GetY(),primVertex->GetZ()};
124 
125  Double_t absoluteP = TMath::Sqrt(TMath::Power(GetPx(),2) + TMath::Power(GetPy(),2) + TMath::Power(GetPz(),2));
126  Double_t p[3] = {GetPx()/absoluteP,GetPy()/absoluteP,GetPz()/absoluteP};
127  Double_t CP[3];
128 
129  CP[0] = fConversionPoint[0] - primCo[0];
130  CP[1] = fConversionPoint[1] - primCo[1];
131  CP[2] = fConversionPoint[2] - primCo[2];
132 
133  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]);
134 
135  Double_t S[3];
136  S[0] = fConversionPoint[0] + p[0]*Lambda;
137  S[1] = fConversionPoint[1] + p[1]*Lambda;
138  S[2] = fConversionPoint[2] + p[2]*Lambda;
139 
140  fDCArPrimVtx = TMath::Sqrt( TMath::Power(primCo[0]-S[0],2) + TMath::Power(primCo[1]-S[1],2));
141  fDCAzPrimVtx = primCo[2]-S[2];
142 
143  return;
144 }
145 
146 
147 void AliAODConversionPhoton::SetCaloPhotonMCFlags(AliMCEvent *mcEvent, Bool_t enableSort){
148 
149  Bool_t isPhoton = kFALSE; // largest contribution to cluster is photon
150  Bool_t isElectron = kFALSE; // largest contribution to cluster is electron
151  Bool_t isConversion = kFALSE; // largest contribution to cluster is converted electron
152  Bool_t isConversionFullyContained = kFALSE; // largest contribution to cluster is converted electron, second electron has been found in same cluster
153  Bool_t isMerged = kFALSE; // cluster contains more than one particle from the same decay
154  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
155  Bool_t isDalitz = kFALSE; // this cluster was created by a particle stemming from a dalitz decay
156  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
157  Bool_t isPhotonWithElecMother = kFALSE; // this cluster is from a photon with an electron as mother
158  Bool_t isShower = kFALSE; // this cluster contains as a largest contribution a particle from a shower or radiative process
159  Bool_t isSubLeadingEM = kFALSE; // cluster contains at least one electron or photon from a pi0, eta or eta_prime in subleading contribution
160  Bool_t isElectronFromFragPhoton = kFALSE; // largest contribution to cluster is from converted electron, but photon stems from fragmentation photon ( q -> q gamma)
161 
162 
163 
164  TParticle* Photon = 0x0;
165  if (fNCaloPhotonMCLabels==0) return;
166 
167  if (enableSort){
168  // sort the array according to the energy of contributing particles
169  if (fNCaloPhotonMCLabels>1){
170  // cout << "start sorting" << endl;
171  Int_t* sortIdx = new Int_t[fNCaloPhotonMCLabels];
172  Double_t* energyPerPart = new Double_t[fNCaloPhotonMCLabels];
173  Long_t* orginalContrib = new Long_t[fNCaloPhotonMCLabels];
174  for (Int_t i = 0; i < fNCaloPhotonMCLabels; i++){
175  orginalContrib[i] = fCaloPhotonMCLabels[i];
176  if (fCaloPhotonMCLabels[i]> -1){
177  TParticle* dummy = mcEvent->Particle(fCaloPhotonMCLabels[i]);
178  energyPerPart[i] = dummy->Energy();
179  // suppress energy of hadrons !!! DIRTY hack !!!
180  if (!(TMath::Abs(dummy->GetPdgCode())== 11 || TMath::Abs(dummy->GetPdgCode())== 22)){
181  // cout << "suppressed hadron energy for:" << dummy->GetPdgCode() << endl;
182  energyPerPart[i]= 0.25; // put energy to mip
183  }
184  } else {
185  energyPerPart[i] = 0;
186  }
187  }
188 
189  TMath::Sort(fNCaloPhotonMCLabels,energyPerPart,sortIdx);
190  for(Int_t index = 0; index < fNCaloPhotonMCLabels; index++) {
191  fCaloPhotonMCLabels[index] = orginalContrib [sortIdx[index]] ;
192 
193  }
194  delete [] sortIdx;
195  delete [] energyPerPart;
196  delete [] orginalContrib;
197  }
198  }
199 
200  if(GetCaloPhotonMCLabel(0)>-1) Photon = mcEvent->Particle(GetCaloPhotonMCLabel(0));
201 
202  if(Photon == NULL){
203  return;
204  }
205 
206  Int_t particleMotherLabel = Photon->GetMother(0);
207  Int_t particleGrandMotherLabel = -1;
208  Int_t particleGrandMotherX2Label = -1;
209  Int_t particleGrandMotherX3Label = -1;
210  Int_t particleGrandMotherX4Label = -1;
211  Int_t particleGrandMotherX5Label = -1;
212  Int_t particleGrandMotherX6Label = -1;
213  Int_t particleGrandMotherX7Label = -1;
214 
215  Int_t particleMotherPDG = -1;
216  Int_t particleGrandMotherPDG = -1;
217  Int_t particleGrandMotherX2PDG = -1;
218  Int_t particleGrandMotherX3PDG = -1;
219  Int_t particleGrandMotherX4PDG = -1;
220  Int_t particleGrandMotherX5PDG = -1;
221  Int_t particleGrandMotherX6PDG = -1;
222 
223  Int_t particleMotherNDaugthers = 0;
224  Int_t particleGrandMotherNDaugthers = 0;
225 
226  if (particleMotherLabel > -1){
227  particleMotherNDaugthers = mcEvent->Particle(Photon->GetMother(0))->GetNDaughters();
228  particleGrandMotherLabel = mcEvent->Particle(Photon->GetMother(0))->GetMother(0);
229  particleMotherPDG = mcEvent->Particle(Photon->GetMother(0))->GetPdgCode();
230  if (particleGrandMotherLabel > -1){
231  particleGrandMotherPDG = mcEvent->Particle(mcEvent->Particle(Photon->GetMother(0))->GetMother(0))->GetPdgCode();
232  particleGrandMotherNDaugthers = mcEvent->Particle(mcEvent->Particle(Photon->GetMother(0))->GetMother(0))->GetNDaughters();
233  particleGrandMotherX2Label = mcEvent->Particle(mcEvent->Particle(Photon->GetMother(0))->GetMother(0))->GetMother(0);
234  if (particleGrandMotherX2Label > -1){
235  particleGrandMotherX2PDG = mcEvent->Particle(mcEvent->Particle(mcEvent->Particle(Photon->GetMother(0))->GetMother(0))->GetMother(0))->GetPdgCode();
236  particleGrandMotherX3Label = mcEvent->Particle(mcEvent->Particle(mcEvent->Particle(Photon->GetMother(0))->GetMother(0))->GetMother(0))->GetMother(0);
237  if (particleGrandMotherX3Label > -1){
238  particleGrandMotherX3PDG = mcEvent->Particle(mcEvent->Particle(mcEvent->Particle(mcEvent->Particle(Photon->GetMother(0))->GetMother(0))->GetMother(0))->GetMother(0))->GetPdgCode();
239  particleGrandMotherX4Label = mcEvent->Particle(mcEvent->Particle(mcEvent->Particle(mcEvent->Particle(Photon->GetMother(0))->GetMother(0))->GetMother(0))->GetMother(0))->GetMother(0);
240  if (particleGrandMotherX4Label > -1){
241  particleGrandMotherX4PDG = mcEvent->Particle(mcEvent->Particle(mcEvent->Particle(mcEvent->Particle(mcEvent->Particle(Photon->GetMother(0))->GetMother(0))->GetMother(0))->GetMother(0))->GetMother(0))->GetPdgCode();
242  particleGrandMotherX5Label = mcEvent->Particle(mcEvent->Particle(mcEvent->Particle(mcEvent->Particle(mcEvent->Particle(Photon->GetMother(0))->GetMother(0))->GetMother(0))->GetMother(0))->GetMother(0))->GetMother(0);
243  if (particleGrandMotherX5Label > -1){
244  particleGrandMotherX5PDG = mcEvent->Particle(mcEvent->Particle(mcEvent->Particle(mcEvent->Particle(mcEvent->Particle(mcEvent->Particle(Photon->GetMother(0))->GetMother(0))->GetMother(0))->GetMother(0))->GetMother(0))->GetMother(0))->GetPdgCode();
245  particleGrandMotherX6Label = mcEvent->Particle(mcEvent->Particle(mcEvent->Particle(mcEvent->Particle(mcEvent->Particle(mcEvent->Particle(Photon->GetMother(0))->GetMother(0))->GetMother(0))->GetMother(0))->GetMother(0))->GetMother(0))->GetMother(0);
246  if (particleGrandMotherX6Label > -1){
247  particleGrandMotherX6PDG = mcEvent->Particle(mcEvent->Particle(mcEvent->Particle(mcEvent->Particle(mcEvent->Particle(mcEvent->Particle(mcEvent->Particle(Photon->GetMother(0))->GetMother(0))->GetMother(0))->GetMother(0))->GetMother(0))->GetMother(0))->GetMother(0))->GetPdgCode();
248  particleGrandMotherX7Label = mcEvent->Particle(mcEvent->Particle(mcEvent->Particle(mcEvent->Particle(mcEvent->Particle(mcEvent->Particle(mcEvent->Particle(Photon->GetMother(0))->GetMother(0))->GetMother(0))->GetMother(0))->GetMother(0))->GetMother(0))->GetMother(0))->GetMother(0);
249  }
250  }
251  }
252  }
253  }
254  }
255  }
256 
257  //determine mother/grandmother of leading particle and if it is pion/eta/eta_prime: fill array fCaloPhotonMotherMCLabels at position 0
258  if (particleMotherLabel > -1){
259  if( TMath::Abs(particleMotherPDG) == 111 || TMath::Abs(particleMotherPDG) == 221 || TMath::Abs(particleMotherPDG) == 331 ){
260  fCaloPhotonMotherMCLabels[0] = particleMotherLabel;
262  } else if (particleGrandMotherLabel > -1 && TMath::Abs(particleMotherPDG) == 22 && (TMath::Abs(particleGrandMotherPDG) == 111 || TMath::Abs(particleGrandMotherPDG) == 221 || TMath::Abs(particleGrandMotherPDG) == 331) ){
263  fCaloPhotonMotherMCLabels[0] = particleGrandMotherLabel;
265  } else {
266  fCaloPhotonMotherMCLabels[0] = particleMotherLabel;
268  }
269  }
270 
271  // Check whether the first contribution was photon
272  if(TMath::Abs(mcEvent->Particle(GetCaloPhotonMCLabel(0))->GetPdgCode()) == 22){
273  isPhoton = kTRUE;
274  // did it decay via the dalitz channel
275  if (particleMotherLabel > -1 && particleMotherNDaugthers == 3)
276  isDalitz = kTRUE;
277  // Test whether particle stems from a shower or radiation
278  if (TMath::Abs(particleMotherPDG) == 11){ // check whether photon stems from electron
279  isPhotonWithElecMother = kTRUE;
280  if (particleGrandMotherLabel > -1){ // test whether first particle has a grandmother
281  if (TMath::Abs(particleGrandMotherPDG) == 22 )
282  isShower = kTRUE; // check whether grandmother is a photon (meaning this is most likely a shower)
283  }
284  }
285  }
286 
287  // Check whether the first contribution was electron
288  if( TMath::Abs(mcEvent->Particle(GetCaloPhotonMCLabel(0))->GetPdgCode()) == 11 ){
289  isElectron = kTRUE;
290  if (particleMotherLabel > -1) {
291  // was it a conversion
292  if (TMath::Abs(particleMotherPDG) == 22)
293  isConversion = kTRUE;
294  // did it decay via the dalitz channel
295  if (particleGrandMotherLabel > -1 && particleGrandMotherNDaugthers == 3 )
296  isDalitz = kTRUE;
297  }
298  if (particleGrandMotherLabel > -1){ // check whether electron has a grandmother
299  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)
300  isShower = kTRUE;
301  }
302  }
303  // consider the rare case, where the conversion electron stems from photon, which stems from electron, which stems from photon, ...
304  if (isConversion && TMath::Abs(particleGrandMotherPDG) == 11 && particleGrandMotherX2PDG == 22){
305 // printf("gamma -> electron -> gamma -> electron!\n");
306  SetCaloPhotonMCLabel(0,particleGrandMotherLabel);
307  fCaloPhotonMotherMCLabels[0] = particleGrandMotherX3Label;
308  if (TMath::Abs(particleGrandMotherX3PDG) == 11 && particleGrandMotherX4PDG == 22 ){
309 // printf("gamma -> electron -> gamma -> electron -> gamma -> electron!\n");
310  SetCaloPhotonMCLabel(0,particleGrandMotherX3Label);
311  fCaloPhotonMotherMCLabels[0] = particleGrandMotherX5Label;
312  if (TMath::Abs(particleGrandMotherX5PDG) == 11 && particleGrandMotherX6PDG == 22 ){
313 // printf("gamma -> electron -> gamma -> electron -> gamma -> electron -> gamma -> electron!\n");
314  SetCaloPhotonMCLabel(0,particleGrandMotherX5Label);
315  fCaloPhotonMotherMCLabels[0] = particleGrandMotherX7Label;
316  }
317  }
318  }
319 
320  // consider the case, where a photon stems from a photon which stems from a photon etc...which is common for frag. photons
321  TParticle *dummyMother = mcEvent->Particle(particleMotherLabel);
322  while (dummyMother->GetPdgCode() == 22){ // follow conversion photon's history, as long as the mother is a photon
323  dummyMother = mcEvent->Particle(dummyMother->GetMother(0));
324  if (TMath::Abs(dummyMother->GetPdgCode()) == 11) // in case of additional conversion skip to photon's grandma, which should be a photon
325  dummyMother = mcEvent->Particle(dummyMother->GetMother(0));
326  isElectronFromFragPhoton = (TMath::Abs(dummyMother->GetPdgCode()) < 6);// photon stems from quark = fragmentation photon
327  }
328 
329  }
330 
331  Bool_t enablePrintOuts = kFALSE;
332 // if (fNCaloPhotonMCLabels>2)
333 // enablePrintOuts = kTRUE;
334 
335  // check whether there were other contributions to the cluster
336  if (fNCaloPhotonMCLabels>1){
337  TParticle* dummyPart =NULL;
338 
339  for (Int_t i = 1; i< fNCaloPhotonMCLabels; i++){
340  if (i > 49) continue; // abort if more than 50 entries to the cluster have been checked (more are not stored in these objects)
341  if (enablePrintOuts) cout << "checking particle: " << i << endl;
342  if (GetCaloPhotonMCLabel(i) < 0) continue;
343  dummyPart = mcEvent->Particle(GetCaloPhotonMCLabel(i));
344  Int_t dummyPartMotherLabel = dummyPart->GetMother(0);
345  Int_t dummyPartGrandMotherLabel = -1;
346  Int_t dummyPartMotherPDG = -1;
347  Int_t dummyPartGrandMotherPDG = -1;
348  // check whether this particle has a mother & obtain the pdg code
349  if (dummyPartMotherLabel > -1){
350  dummyPartGrandMotherLabel = mcEvent->Particle(dummyPart->GetMother(0))->GetMother(0);
351  dummyPartMotherPDG = mcEvent->Particle(dummyPart->GetMother(0))->GetPdgCode();
352  // check whether this particle has a grandmother & obtain its pdg code
353  if (dummyPartGrandMotherLabel > -1){
354  dummyPartGrandMotherPDG = mcEvent->Particle(mcEvent->Particle(dummyPart->GetMother(0))->GetMother(0))->GetPdgCode();
355  }
356  }
357  // largest contribution was from photon and is not from shower or electron mother
358  if (isPhoton && (!isShower || !isPhotonWithElecMother )){
359  if (enablePrintOuts) cout << "lead gamma" << endl;
360  if (dummyPartMotherLabel > -1 && particleMotherLabel > -1){ // test whether first particle has a mother
361  if (dummyPartMotherLabel == particleMotherLabel)
362  isMerged = kTRUE; // test whether current and first particle have the same mother => i.e. other gamma from decay or dalitz electron
363  if (dummyPartGrandMotherLabel > -1){ // test whether first particle has a grandmother
364  // check whether particle is an electron from a conversion of a photon from the original mother
365  if (TMath::Abs(dummyPart->GetPdgCode()) == 11 && dummyPartGrandMotherLabel == particleMotherLabel )
366  isMergedPartConv = kTRUE;
367  // check whether particle is an electron from a dalitz decay from the original mother
368  if (TMath::Abs(dummyPart->GetPdgCode()) == 11 && dummyPartMotherLabel == particleMotherLabel )
369  isDalitzMerged = kTRUE;
370  }
371  }
372  }
373 
374  // largest contribution was from electron & not a from a shower
375  if (isElectron && !isShower){
376  if (enablePrintOuts) cout << "lead electron" << endl;
377  if (dummyPartMotherLabel > -1 && particleMotherLabel > -1){ // test whether first particle has a mother
378  if ( TMath::Abs(dummyPart->GetPdgCode()) == 11 && isConversion && dummyPartMotherLabel == particleMotherLabel ) {
379  isConversionFullyContained = kTRUE; // test whether conversion is fully contained in cluster
380  if (enablePrintOuts) cout << "found full conversion" << endl;
381  } else if (TMath::Abs(dummyPart->GetPdgCode()) == 11) { // current particle is an electron
382  if (enablePrintOuts) cout << "current is electron" << endl;
383  // check whether current particle is a shower
384  Bool_t isShowerCurrent = kFALSE;
385  if (TMath::Abs(dummyPartGrandMotherPDG) == 11)
386  isShowerCurrent = kTRUE;
387  if (TMath::Abs(dummyPartGrandMotherPDG) == 22)
388  isShowerCurrent = kTRUE;
389 
390  // check whether current particle is a conversion
391  Bool_t isConvCurrent = kFALSE;
392  if (TMath::Abs(dummyPartMotherPDG) == 22 && !isShowerCurrent)
393  isConvCurrent = kTRUE;
394 
395  if (enablePrintOuts) cout << "conv current: " << isConvCurrent << "\t shower current: " << isShowerCurrent << endl;
396  // check whether orginal electron and this electron stem from the same particle and electron originated in dalitz decay
397  if( dummyPartMotherLabel == particleMotherLabel && TMath::Abs(particleMotherPDG) != 22 && !isShowerCurrent ) {
398  isDalitzMerged = kTRUE;
399  isMerged = kTRUE;
400  // both particles are from conversions
401  } else if (isConversion && isConvCurrent){
402  if (particleGrandMotherLabel > -1 && dummyPartGrandMotherLabel > -1){ // test whether first particle has a grandmother and current particle has one as well
403  // check whether orginal electron and this electron stem from the same particle and electron stems from conversion
404  if( dummyPartGrandMotherLabel == particleGrandMotherLabel && !isShowerCurrent )
405  isMergedPartConv = kTRUE;
406  }
407  // check whether current electron is from conversion & orginal e are from same mother
408  } else if (dummyPartGrandMotherLabel > -1 && isConvCurrent){
409  if (dummyPartGrandMotherLabel == particleMotherLabel && (TMath::Abs(particleMotherPDG) == 111 || TMath::Abs(particleMotherPDG) == 221 || TMath::Abs(particleMotherPDG) == 331) ){
410  isDalitzMerged = kTRUE;
411  isMerged = kTRUE;
412  }
413  // check whether current electron & orginal e is from conversion are from same mother
414  } else if (isConversion && particleGrandMotherLabel > -1){
415  if (dummyPartMotherLabel == particleGrandMotherLabel && (TMath::Abs(dummyPartMotherPDG) == 111 || TMath::Abs(dummyPartMotherPDG) == 221 || TMath::Abs(dummyPartMotherPDG) == 331) ){
416  isDalitzMerged = kTRUE;
417  isMerged = kTRUE;
418  }
419  }
420  } else if (TMath::Abs(dummyPart->GetPdgCode()) == 22){ // test whether this particle is a photon
421  // check whether current particle is a shower
422  Bool_t isShowerCurrent = kFALSE;
423  if (TMath::Abs(dummyPartMotherPDG) == 11)
424  isShowerCurrent = kTRUE;
425  if (TMath::Abs(dummyPartMotherPDG) == 22)
426  isShowerCurrent = kTRUE;
427 
428  // check whether orginal electron and this photon stem from the same particle and electron originated in dalitz
429  if( dummyPartMotherLabel == particleMotherLabel && !isShowerCurrent ) {
430  isDalitzMerged = kTRUE;
431  isMerged = kTRUE;
432  } else if (particleGrandMotherLabel > -1){ // test whether first particle has a grandmother
433  // check whether orginal electron and this photon stem from the same particle and electron stems from conversion
434  if( dummyPartMotherLabel == particleGrandMotherLabel && !isShowerCurrent )
435  isMergedPartConv = kTRUE;
436  }
437  }
438  }
439  }
440 
441 
442  if (dummyPartMotherLabel > -1){ // test whether particle has a mother
443  if (TMath::Abs(dummyPart->GetPdgCode()) == 22){ // test whether particle is a photon
444  //check if photon directly comes from a pion/eta/eta_prime decay
445  Bool_t helpN = true;
446  for(Int_t j=0; j<fNCaloPhotonMotherMCLabels; j++){
447  if (fCaloPhotonMotherMCLabels[j] == dummyPartMotherLabel){ //check if mother is already contained in fCaloPhotonMotherMCLabels
448  helpN = false;
449  }
450  }
451  if (helpN){
453  fNCaloPhotonMotherMCLabels++; //only if particle label is not yet contained in array, count up fNCaloPhotonMotherMCLabels
454  }
455  if ((TMath::Abs(dummyPartMotherPDG) == 111 || TMath::Abs(dummyPartMotherPDG) == 221 || TMath::Abs(dummyPartMotherPDG) == 331) && !isPhoton && !isElectron)
456  isSubLeadingEM = kTRUE;
457  } else if (TMath::Abs(dummyPart->GetPdgCode()) == 11){ //test whether particle is an electron
458  //check if electron comes from a pion decay
459  if ( TMath::Abs(dummyPartMotherPDG) != 22 ){
460  Bool_t helpN = true;
461  for(Int_t j=0; j<fNCaloPhotonMotherMCLabels; j++){
462  if (fCaloPhotonMotherMCLabels[j] == dummyPartMotherLabel){ //check if mother is already contained in fCaloPhotonMotherMCLabels
463  helpN = false;
464  }
465  }
466  if (helpN){
468  fNCaloPhotonMotherMCLabels++; //only if particle label is not yet contained in array, count up fNCaloPhotonMotherMCLabels
469  }
470  if ((TMath::Abs(dummyPartMotherPDG) == 111 || TMath::Abs(dummyPartMotherPDG) == 221 || TMath::Abs(dummyPartMotherPDG) == 331) && !isPhoton && !isElectron)
471  isSubLeadingEM = kTRUE;
472  } else if (dummyPartGrandMotherLabel > -1){ //if it is not a dalitz decay, test whether particle has a grandmother
473  //check if it is a conversion electron that has pion/eta/eta_prime as grandmother
474  if ( (TMath::Abs(dummyPartGrandMotherPDG) == 111 || TMath::Abs(dummyPartGrandMotherPDG) == 221 || TMath::Abs(dummyPartGrandMotherPDG) == 331)){
475 
476  Bool_t helpN = true;
477  for(Int_t j=0; j<fNCaloPhotonMotherMCLabels; j++){
478  if (fCaloPhotonMotherMCLabels[j] == dummyPartGrandMotherLabel){ //check if grandmother is already contained in fCaloPhotonMotherMCLabels
479  helpN = false;
480  }
481  }
482  if (helpN){
483  fCaloPhotonMotherMCLabels[fNCaloPhotonMotherMCLabels] = dummyPartGrandMotherLabel;
484  fNCaloPhotonMotherMCLabels++; //only if particle label is not yet contained in array, count up fNCaloPhotonMotherMCLabels
485  }
486  if (!isPhoton && !isElectron)
487  isSubLeadingEM = kTRUE;
488  }
489  }
490  }
491  }
492  }
493  }
494  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;
495 }
496 
497 void AliAODConversionPhoton::SetCaloPhotonMCFlagsAOD(AliVEvent* event, Bool_t enableSort){
498 
499  TClonesArray *AODMCTrackArray = dynamic_cast<TClonesArray*>(event->FindListObject(AliAODMCParticle::StdBranchName()));
500  if (!AODMCTrackArray) return;
501 
502  Bool_t isPhoton = kFALSE; // largest contribution to cluster is photon
503  Bool_t isElectron = kFALSE; // largest contribution to cluster is electron
504  Bool_t isConversion = kFALSE; // largest contribution to cluster is converted electron
505  Bool_t isConversionFullyContained = kFALSE; // largest contribution to cluster is converted electron, second electron has been found in same cluster
506  Bool_t isMerged = kFALSE; // largest contribution to cluster is photon, second photon or electron from dalitz has been found in same cluster
507  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
508  Bool_t isDalitz = kFALSE; // this cluster was created by a particle stemming from a dality decay
509  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
510  Bool_t isPhotonWithElecMother = kFALSE; // this cluster is from a photon with an electron as mother
511  Bool_t isShower = kFALSE; // this cluster contains as a largest contribution a particle from a shower or radiative process
512  Bool_t isSubLeadingEM = kFALSE; // cluster contains at least one electron or photon from a pi0 or eta in subleading contribution
513 
514  AliAODMCParticle* Photon;
515  AliAODMCParticle* PhotonMother;
516  AliAODMCParticle* PhotonGrandMother;
517 
518  if (fNCaloPhotonMCLabels==0) return;
519  Photon = (AliAODMCParticle*) AODMCTrackArray->At(GetCaloPhotonMCLabel(0));
520 
521  if (enableSort){
522  // sort the array according to the energy of contributing particles
523  if (fNCaloPhotonMCLabels>1){
524  // cout << "start sorting" << endl;
525  Int_t* sortIdx = new Int_t[fNCaloPhotonMCLabels];
526  Double_t* energyPerPart = new Double_t[fNCaloPhotonMCLabels];
527  Long_t* orginalContrib = new Long_t[fNCaloPhotonMCLabels];
528  for (Int_t i = 0; i < fNCaloPhotonMCLabels; i++){
529  orginalContrib[i] = fCaloPhotonMCLabels[i];
530  if (fCaloPhotonMCLabels[i]> -1){
531  AliAODMCParticle* dummy = (AliAODMCParticle*) AODMCTrackArray->At(GetCaloPhotonMCLabel(i));
532  energyPerPart[i] = dummy->E();
533  // suppress energy of hadrons !!! DIRTY hack !!!
534  if (!(TMath::Abs(dummy->GetPdgCode())== 11 || TMath::Abs(dummy->GetPdgCode())== 22)){
535  // cout << "suppressed hadron energy" << endl;
536  energyPerPart[i]= 0.25; // put energy to mip
537  }
538  } else {
539  energyPerPart[i] = 0;
540  }
541  }
542 
543  TMath::Sort(fNCaloPhotonMCLabels,energyPerPart,sortIdx);
544  for(Int_t index = 0; index < fNCaloPhotonMCLabels; index++) {
545  fCaloPhotonMCLabels[index] = orginalContrib [sortIdx[index]] ;
546 
547  }
548  delete [] sortIdx;
549  delete [] energyPerPart;
550  delete [] orginalContrib;
551  }
552  }
553 
554  if(Photon == NULL){
555  return;
556  }
557 
558  Int_t particleMotherLabel = Photon->GetMother();
559  Int_t particleGrandMotherLabel = -1;
560  Int_t particleMotherPDG = -1;
561  Int_t particleGrandMotherPDG = -1;
562  Int_t particleMotherNDaugthers = 0;
563  Int_t particleGrandMotherNDaugthers = 0;
564  if (particleMotherLabel > -1){
565  PhotonMother = (AliAODMCParticle*) AODMCTrackArray->At(Photon->GetMother());
566  particleMotherNDaugthers = PhotonMother->GetNDaughters();
567  particleGrandMotherLabel = PhotonMother->GetMother();
568  particleMotherPDG = PhotonMother->GetPdgCode();
569  if (particleGrandMotherLabel > -1){
570  PhotonGrandMother = (AliAODMCParticle*) AODMCTrackArray->At(PhotonMother->GetMother());
571  particleGrandMotherPDG = PhotonGrandMother->GetPdgCode();
572  particleGrandMotherNDaugthers = PhotonGrandMother->GetNDaughters();
573  }
574  }
575 
576  //determine mother/grandmother of leading particle and if it is pion/eta/eta_prime: fill array fCaloPhotonMotherMCLabels at position 0
577  if (particleMotherLabel > -1){
578  if( TMath::Abs(particleMotherPDG) == 111 || TMath::Abs(particleMotherPDG) == 221 || TMath::Abs(particleMotherPDG) == 331 ){
579  fCaloPhotonMotherMCLabels[0] = particleMotherLabel;
581  } else if (TMath::Abs(particleMotherPDG) == 22 && particleGrandMotherLabel > -1){
582  if ( TMath::Abs(particleMotherPDG) == 22 && (TMath::Abs(particleGrandMotherPDG) == 111 || TMath::Abs(particleGrandMotherPDG) == 221 || TMath::Abs(particleGrandMotherPDG) == 331) ){
583  fCaloPhotonMotherMCLabels[0] = particleGrandMotherLabel;
585  }
586  } else {
587  fCaloPhotonMotherMCLabels[0] = particleMotherLabel;
589  }
590  }
591 
592  // Check whether the first contribution was photon
593  if(TMath::Abs(Photon->GetPdgCode()) == 22){
594  isPhoton = kTRUE;
595  // did it decay via the dalitz channel
596  if (particleMotherLabel > -1 && particleMotherNDaugthers == 3)
597  isDalitz = kTRUE;
598  // Test whether particle stems from a shower or radiation
599  if (TMath::Abs(particleMotherPDG) == 11){ // check whether photon stems from electron
600  isPhotonWithElecMother = kTRUE;
601  if (particleGrandMotherLabel > -1){ // test whether first particle has a grandmother
602  if (TMath::Abs(particleGrandMotherPDG) == 22 )
603  isShower = kTRUE; // check whether grandmother is a photon (meaning this is most likely a shower)
604  }
605  }
606  }
607  // Check whether the first contribution was electron
608  if(TMath::Abs(Photon->GetPdgCode()) == 11 ){
609  isElectron = kTRUE;
610  if (particleMotherLabel > -1) {
611  // was it a conversion
612  if (TMath::Abs(particleMotherPDG) == 22)
613  isConversion = kTRUE;
614  // did it decay via the dalitz channel
615  if (particleGrandMotherLabel > -1 && particleGrandMotherNDaugthers == 3 )
616  isDalitz = kTRUE;
617  }
618  if (particleGrandMotherLabel > -1){ // check whether electron has a grandmother
619  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)
620  isShower = kTRUE;
621  }
622  }
623  }
624 
625 
626  // check whether there were other contributions to the cluster
627  if (fNCaloPhotonMCLabels>1){
628  AliAODMCParticle* dummyPart = NULL;
629  AliAODMCParticle* dummyPartMother = NULL;
630  AliAODMCParticle* dummyPartGrandMother = NULL;
631  for (Int_t i = 1; i< fNCaloPhotonMCLabels; i++){
632  if (i > 49) continue; // abort if more than 50 entries to the cluster have been checked (more are not stored in these objects)
633  dummyPart = (AliAODMCParticle*) AODMCTrackArray->At(GetCaloPhotonMCLabel(i));
634  Int_t dummyPartMotherLabel = dummyPart->GetMother();
635  Int_t dummyPartGrandMotherLabel = -1;
636  Int_t dummyPartMotherPDG = -1;
637  Int_t dummyPartGrandMotherPDG = -1;
638 
639  // check whether this particle has a mother & obtain the pdg code
640  if (dummyPartMotherLabel > -1){
641  dummyPartMother = (AliAODMCParticle*) AODMCTrackArray->At(dummyPart->GetMother());
642  dummyPartGrandMotherLabel = dummyPartMother->GetMother();
643  dummyPartMotherPDG = dummyPartMother->GetPdgCode();
644  // check whether this particle has a grandmother & obtain its pdg code
645  if (dummyPartGrandMotherLabel > -1){
646  dummyPartGrandMother = (AliAODMCParticle*) AODMCTrackArray->At(dummyPartMother->GetMother());
647  dummyPartGrandMotherPDG = dummyPartGrandMother->GetPdgCode();
648  }
649  }
650 
651 
652  // largest contribution was from photon and is not from shower or electron mother
653  if (isPhoton && (!isShower || !isPhotonWithElecMother )){
654  if (dummyPartMotherLabel > -1 && particleMotherLabel > -1){ // test whether first particle has a mother
655  if (dummyPartMotherLabel == particleMotherLabel)
656  isMerged = kTRUE; // test whether current and first particle have the same mother => i.e. other gamma from decay or dalitz electron
657  if (dummyPartGrandMotherLabel > -1){ // test whether first particle has a grandmother
658  // check whether particle is an electron from a conversion of a photon from the original mother
659  if (TMath::Abs(dummyPart->GetPdgCode()) == 11 && dummyPartGrandMotherLabel == particleMotherLabel )
660  isMergedPartConv = kTRUE;
661  // check whether particle is an electron from a dalitz decay from the original mother
662  if (TMath::Abs(dummyPart->GetPdgCode()) == 11 && dummyPartMotherLabel == particleMotherLabel )
663  isDalitzMerged = kTRUE;
664  }
665  }
666  }
667 
668  // largest contribution was from electron & not a from a shower
669  if (isElectron && !isShower){
670  if (dummyPartMotherLabel > -1 && particleMotherLabel > -1){ // test whether first particle has a mother
671  if ( TMath::Abs(dummyPart->GetPdgCode()) == 11 && isConversion && dummyPartMotherLabel == particleMotherLabel ) {
672  isConversionFullyContained = kTRUE; // test whether conversion is fully contained in cluster
673  } else if (TMath::Abs(dummyPart->GetPdgCode()) == 11) { // current particle is an electron
674 
675  // check whether current particle is a shower
676  Bool_t isShowerCurrent = kFALSE;
677  if (TMath::Abs(dummyPartGrandMotherPDG) == 11)
678  isShowerCurrent = kTRUE;
679  if (TMath::Abs(dummyPartGrandMotherPDG) == 22)
680  isShowerCurrent = kTRUE;
681 
682  // check whether current particle is a conversion
683  Bool_t isConvCurrent = kFALSE;
684  if (TMath::Abs(dummyPartMotherPDG) == 22 && !isShowerCurrent)
685  isConvCurrent = kTRUE;
686 
687  // check whether orginal electron and this electron stem from the same particle and electron originated in dalitz decay
688  if( dummyPartMotherLabel == particleMotherLabel && TMath::Abs(particleMotherPDG) != 22 && !isShowerCurrent ) {
689  isDalitzMerged = kTRUE;
690  isMerged = kTRUE;
691  // both particles are from conversions
692  } else if (isConversion && isConvCurrent){
693  if (particleGrandMotherLabel > -1 && dummyPartGrandMotherLabel > -1){ // test whether first particle has a grandmother and current particle has one as well
694  // check whether orginal electron and this electron stem from the same particle and electron stems from conversion
695  if( dummyPartGrandMotherLabel == particleGrandMotherLabel && !isShowerCurrent )
696  isMergedPartConv = kTRUE;
697  }
698  // check whether current electron is from conversion & orginal e are from same mother
699  } else if (dummyPartGrandMotherLabel > -1 && isConvCurrent){
700  if (dummyPartGrandMotherLabel == particleMotherLabel && (TMath::Abs(particleMotherPDG) == 111 || TMath::Abs(particleMotherPDG) == 221 || TMath::Abs(particleMotherPDG) == 331) ){
701  isDalitzMerged = kTRUE;
702  isMerged = kTRUE;
703  }
704  // check whether current electron & orginal e is from conversion are from same mother
705  } else if (isConversion && particleGrandMotherLabel > -1){
706  if (dummyPartMotherLabel == particleGrandMotherLabel && (TMath::Abs(dummyPartMotherPDG) == 111 || TMath::Abs(dummyPartMotherPDG) == 221 || TMath::Abs(dummyPartMotherPDG) == 331) ){
707  isDalitzMerged = kTRUE;
708  isMerged = kTRUE;
709  }
710  }
711  } else if (TMath::Abs(dummyPart->GetPdgCode()) == 22){ // test whether this particle is a photon
712  // check whether current particle is a shower
713  Bool_t isShowerCurrent = kFALSE;
714  if (TMath::Abs(dummyPartMotherPDG) == 11)
715  isShowerCurrent = kTRUE;
716  if (TMath::Abs(dummyPartMotherPDG) == 22)
717  isShowerCurrent = kTRUE;
718 
719  // check whether orginal electron and this photon stem from the same particle and electron originated in dalitz
720  if( dummyPartMotherLabel == particleMotherLabel && !isShowerCurrent ) {
721  isDalitzMerged = kTRUE;
722  isMerged = kTRUE;
723  } else if (particleGrandMotherLabel > -1){ // test whether first particle has a grandmother
724  // check whether orginal electron and this photon stem from the same particle and electron stems from conversion
725  if( dummyPartMotherLabel == particleGrandMotherLabel && !isShowerCurrent )
726  isMergedPartConv = kTRUE;
727  }
728  }
729  }
730  }
731 
732  if (dummyPartMotherLabel > -1){ // test whether particle has a mother
733  if (TMath::Abs(dummyPart->GetPdgCode()) == 22){ // test whether particle is a photon
734  //check if photon directly comes from a pion/eta/eta_prime decay
735  Bool_t helpN = true;
736  for(Int_t j=0; j<fNCaloPhotonMotherMCLabels; j++){
737  if (fCaloPhotonMotherMCLabels[j] == dummyPartMotherLabel){ //check if mother is already contained in fCaloPhotonMotherMCLabels
738  helpN = false;
739  }
740  }
741  if (helpN){
743  fNCaloPhotonMotherMCLabels++; //only if particle label is not yet contained in array, count up fNCaloPhotonMotherMCLabels
744  }
745  if ((TMath::Abs(dummyPartMotherPDG) == 111 || TMath::Abs(dummyPartMotherPDG) == 221 || TMath::Abs(dummyPartMotherPDG) == 331) && !isPhoton && !isElectron)
746  isSubLeadingEM = kTRUE;
747  } else if (TMath::Abs(dummyPart->GetPdgCode()) == 11){ //test whether particle is an electron
748  //check if electron comes from a pion decay
749  if ( TMath::Abs(dummyPartMotherPDG) != 22 ){
750  Bool_t helpN = true;
751  for(Int_t j=0; j<fNCaloPhotonMotherMCLabels; j++){
752  if (fCaloPhotonMotherMCLabels[j] == dummyPartMotherLabel){ //check if mother is already contained in fCaloPhotonMotherMCLabels
753  helpN = false;
754  }
755  }
756  if (helpN){
758  fNCaloPhotonMotherMCLabels++; //only if particle label is not yet contained in array, count up fNCaloPhotonMotherMCLabels
759  }
760  if ((TMath::Abs(dummyPartMotherPDG) == 111 || TMath::Abs(dummyPartMotherPDG) == 221 || TMath::Abs(dummyPartMotherPDG) == 331) && !isPhoton && !isElectron)
761  isSubLeadingEM = kTRUE;
762  } else if (dummyPartGrandMotherLabel > -1){ //if it is not a dalitz decay, test whether particle has a grandmother
763  //check if it is a conversion electron that has pion/eta/eta_prime as grandmother
764  if ( (TMath::Abs(dummyPartGrandMotherPDG) == 111 || TMath::Abs(dummyPartGrandMotherPDG) == 221 || TMath::Abs(dummyPartGrandMotherPDG) == 331)){
765 
766  Bool_t helpN = true;
767  for(Int_t j=0; j<fNCaloPhotonMotherMCLabels; j++){
768  if (fCaloPhotonMotherMCLabels[j] == dummyPartGrandMotherLabel){ //check if grandmother is already contained in fCaloPhotonMotherMCLabels
769  helpN = false;
770  }
771  }
772  if (helpN){
773  fCaloPhotonMotherMCLabels[fNCaloPhotonMotherMCLabels] = dummyPartGrandMotherLabel;
774  fNCaloPhotonMotherMCLabels++; //only if particle label is not yet contained in array, count up fNCaloPhotonMotherMCLabels
775  }
776  if (!isPhoton && !isElectron)
777  isSubLeadingEM = kTRUE;
778  }
779  }
780  }
781  }
782  }
783  }
784  fCaloPhotonMCFlags = isPhoton *1 + isElectron *2 + isConversion*4+ isConversionFullyContained *8 + isMerged *16 + isMergedPartConv*32 + isDalitz *64 + isDalitzMerged *128 + isPhotonWithElecMother *256 + isShower * 512 + isSubLeadingEM * 1024;
785 }
786 
787 //_____________________________________________________________________________
788 // prints information of MC particles contributing to cluster
789 //_____________________________________________________________________________
791  cout << endl << endl << "particles contributing: " << endl;
792  for (Int_t i =0 ; i < fNCaloPhotonMCLabels; i++ ){
793  TParticle *dummy = NULL;
794  if (fCaloPhotonMCLabels[i]>0){
795  dummy = (TParticle*)mcEvent->Particle(fCaloPhotonMCLabels[i]);
796  cout << i << "\t"<< fCaloPhotonMCLabels[i] << "\t pdg code: " <<dummy->GetPdgCode() << "\t prod radius: "<< dummy->R() << "\t energy: " << dummy->Energy() ;
797  if (dummy->GetMother(0) > -1){
798  TParticle* dummyMother = (TParticle*)mcEvent->Particle(dummy->GetMother(0));
799  cout << "\t mother part: " << dummy->GetMother(0) << "\t mother pdg code: " << dummyMother->GetPdgCode() << "\t energy: " << dummyMother->Energy() << endl;
800  if (dummyMother->GetMother(0) > -1){
801  TParticle* dummyGrandMother = (TParticle*)mcEvent->Particle(dummyMother->GetMother(0));
802  cout << "\t grandmother part: " << dummyMother->GetMother(0) << "\t grandmother pdg code: " << dummyGrandMother->GetPdgCode() << "\t energy: " << dummyGrandMother->Energy() << endl;
803  } else {
804  cout << endl;
805  }
806  } else {
807  cout << endl;
808  }
809  }
810 // if (dummy) delete dummy;
811  }
812  cout << "mothers contributing: " << endl;
813  for (Int_t i =0 ; i < fNCaloPhotonMotherMCLabels; i++ ){
814  TParticle *dummy = NULL;
815  if (fCaloPhotonMotherMCLabels[i]>0){
816  dummy = (TParticle*)mcEvent->Particle(fCaloPhotonMotherMCLabels[i]);
817  cout << i << "\t"<< fCaloPhotonMotherMCLabels[i] << "\t pdg code: " <<dummy->GetPdgCode() << "\t prod radius: "<< dummy->R() << "\t energy: " << dummy->Energy() << endl;
818 
819 
820  }
821 // if (dummy) delete dummy;
822  }
823 }
824 
825 //_____________________________________________________________________________
826 // prints information of cluster as it is flagged
827 //_____________________________________________________________________________
829  cout << fCaloPhotonMCFlags << "\t photon: " << IsLargestComponentPhoton() << "\t electron: " << IsLargestComponentElectron() << "\t conversion: " << IsConversion() << "\t conversion full: "
830  << IsConversionFullyContained() << "\t merged: " << IsMerged() << "\t merged p.conv: " << IsMergedPartConv() << "\t Dalitz: " << IsDalitz() << "\t Dalitz merged: " << IsDalitzMerged()
831  << "\t rad: " << IsPhotonWithElecMother() << "\t shower: " << IsShower() << "\t no EM lead: " << IsEMNonLeading() << "\t EM sub lead: " << IsSubLeadingEM() << endl;
832 }
virtual Double_t GetPy() const
double Double_t
Definition: External.C:58
void SetCaloPhotonMCFlags(AliMCEvent *mcEvent, Bool_t enableSort)
virtual Double_t GetPz() const
void SetCaloPhotonMCLabel(Int_t i, Int_t labelCaloPhoton)
AliAODConversionPhoton & operator=(const AliAODConversionPhoton &g)
void PrintCaloMCLabelsAndInfo(AliMCEvent *mcEvent)
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 SetCaloPhotonMCFlagsAOD(AliVEvent *event, Bool_t enableSort)