AliPhysics  6cf2591 (6cf2591)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliNeutralMesonSelection.cxx
Go to the documentation of this file.
1  /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  * *
4  * Author: The ALICE Off-line Project. *
5  * Contributors are mentioned in the code where appropriate. *
6  * *
7  * Permission to use, copy, modify and distribute this software and its *
8  * documentation strictly for non-commercial purposes is hereby granted *
9  * without fee, provided that the above copyright notice appears in all *
10  * copies and that both the copyright notice and this permission notice *
11  * appear in the supporting documentation. The authors make no claims *
12  * about the suitability of this software for any purpose. It is *
13  * provided "as is" without express or implied warranty. *
14  **************************************************************************/
15 
16 // --- ROOT system ---
17 #include <TLorentzVector.h>
18 #include <TH2.h>
19 #include <TList.h>
20 
21 //---- AliRoot system ----
23 #include "AliLog.h"
24 
28 
29 //____________________________________________________
31 //____________________________________________________
33 TObject(), fAsymmetryCut(1),
34 fUseAsymmetryCut(0), fM(0),
35 fInvMassMaxCut(0.), fInvMassMinCut(0.),
36 fLeftBandMinCut(0.), fLeftBandMaxCut(0.),
37 fRightBandMinCut(0.), fRightBandMaxCut(0.),
38 fAngleMaxParam(), fUseAngleCut(0),
39 fKeepNeutralMesonHistos(0),
40 fParticle(""), fDecayBit(0),
41 fDebug(0),
42 // histograms
43 fhAnglePairNoCut(0), fhAnglePairOpeningAngleCut(0),
44 fhAnglePairAsymmetryCut(0), fhAnglePairAllCut(0),
45 fhInvMassPairNoCut(0), fhInvMassPairOpeningAngleCut(0),
46 fhInvMassPairAsymmetryCut(0), fhInvMassPairAllCut(0),
47 fhAsymmetryNoCut(0), fhAsymmetryOpeningAngleCut(0),
48 fhAsymmetryAllCut(0),
49 // histogram ranges and bins
50 fHistoNEBins(0), fHistoEMax(0.), fHistoEMin(0.),
51 fHistoNAngleBins(0), fHistoAngleMax(0.), fHistoAngleMin(0.),
52 fHistoNIMBins(0), fHistoIMMax(0.), fHistoIMMin(0.)
53 {
55 }
56 
57 //_________________________________________________________
60 //_________________________________________________________
62 {
63  TList * outputContainer = new TList() ;
64  outputContainer->SetName("MesonDecayHistos") ;
65 
67 
68  outputContainer->SetOwner(kFALSE);
69 
71  ("AnglePairNoCut",
72  "Angle between all #gamma pair vs E_{#pi^{0}}",fHistoNEBins,fHistoEMin,fHistoEMax,fHistoNAngleBins,fHistoAngleMin,fHistoAngleMax);
73  fhAnglePairNoCut->SetYTitle("Angle (rad)");
74  fhAnglePairNoCut->SetXTitle("E_{ #pi^{0}} (GeV)");
75 
77  ("AsymmetryNoCut","Asymmetry of all #gamma pair vs E_{#pi^{0}}",
79  fhAsymmetryNoCut->SetYTitle("Asymmetry");
80  fhAsymmetryNoCut->SetXTitle("E_{ #pi^{0}} (GeV)");
81 
83  ("InvMassPairNoCut","Invariant Mass of all #gamma pair vs E_{#pi^{0}}",
85  fhInvMassPairNoCut->SetYTitle("Invariant Mass (GeV/c^{2})");
86  fhInvMassPairNoCut->SetXTitle("E_{ #pi^{0}} (GeV)");
87 
88  outputContainer->Add(fhAnglePairNoCut) ;
89  outputContainer->Add(fhAsymmetryNoCut) ;
90  outputContainer->Add(fhInvMassPairNoCut) ;
91 
92  if(fUseAngleCut) {
94  ("AnglePairOpeningAngleCut",
95  "Angle between all #gamma pair (opening angle) vs E_{#pi^{0}}"
97  fhAnglePairOpeningAngleCut->SetYTitle("Angle (rad)");
98  fhAnglePairOpeningAngleCut->SetXTitle("E_{ #pi^{0}} (GeV)");
99 
101  ("AsymmetryOpeningAngleCut",
102  "Asymmetry of #gamma pair (angle cut) vs E_{#pi^{0}}",
104  fhAsymmetryOpeningAngleCut->SetYTitle("Asymmetry");
105  fhAsymmetryOpeningAngleCut->SetXTitle(" E_{#pi^{0}}(GeV)");
106 
108  ("InvMassPairOpeningAngleCut",
109  "Invariant Mass of #gamma pair (angle cut) vs E_{#pi^{0}}",
111  fhInvMassPairOpeningAngleCut->SetYTitle("Invariant Mass (GeV/c^{2})");
112  fhInvMassPairOpeningAngleCut->SetXTitle(" E_{#pi^{0}}(GeV)");
113 
114  outputContainer->Add(fhAnglePairOpeningAngleCut) ;
115  outputContainer->Add(fhAsymmetryOpeningAngleCut) ;
116  outputContainer->Add(fhInvMassPairOpeningAngleCut) ;
117  }
118 
119  if(fUseAsymmetryCut) {
121  ("AnglePairAsymmetryCut",
122  "Angle between all #gamma pair (opening angle + asymetry cut) vs E_{#pi^{0}}"
124  fhAnglePairAsymmetryCut->SetYTitle("Angle (rad)");
125  fhAnglePairAsymmetryCut->SetXTitle("E_{ #pi^{0}} (GeV)");
126 
128  ("InvMassPairAsymmetryCut",
129  "Invariant Mass of #gamma pair (opening angle + asymmetry) vs E_{#pi^{0}}",
131  fhInvMassPairAsymmetryCut->SetYTitle("Invariant Mass (GeV/c^{2})");
132  fhInvMassPairAsymmetryCut->SetXTitle("E_{#pi^{0}}(GeV)");
133 
134  outputContainer->Add(fhAnglePairAsymmetryCut) ;
135  outputContainer->Add(fhInvMassPairAsymmetryCut) ;
136 
137  }
138 
139  fhAnglePairAllCut = new TH2F
140  ("AnglePairAllCut",
141  "Angle between all #gamma pair (opening angle + asymmetry + inv mass cut) vs E_{#pi^{0}}"
143  fhAnglePairAllCut->SetYTitle("Angle (rad)");
144  fhAnglePairAllCut->SetXTitle("E_{ #pi^{0}} (GeV)");
145 
147  ("InvMassPairAllCut",
148  "Invariant Mass of #gamma pair (opening angle + asymmetry + invmass cut) vs E_{#pi^{0}}",
150  fhInvMassPairAllCut->SetYTitle("Invariant Mass (GeV/c^{2})");
151  fhInvMassPairAllCut->SetXTitle("E_{#pi^{0}}(GeV)");
152 
153  fhAsymmetryAllCut = new TH2F
154  ("AsymmetryAllCut",
155  "Asymmetry of #gamma pair (opening angle+invmass cut) vs E_{#pi^{0}}",
157  fhAsymmetryAllCut->SetYTitle("Asymmetry");
158  fhAsymmetryAllCut->SetXTitle("E_{#pi^{0}}(GeV)");
159 
160  outputContainer->Add(fhAnglePairAllCut) ;
161  outputContainer->Add(fhAsymmetryAllCut) ;
162  outputContainer->Add(fhInvMassPairAllCut) ;
163 
164  }
165 
166  return outputContainer;
167 }
168 
169 //_____________________________________________
171 //_____________________________________________
173 {
174  fAngleMaxParam.Set(4) ;
175  fAngleMaxParam.Reset(0.);
176 
177  SetParticle("Pi0");
178 
179  //Histogrammes settings
180  fHistoNEBins = 200 ;
181  fHistoEMax = 50 ;
182  fHistoEMin = 0. ;
183 
184  fHistoNAngleBins = 200 ;
185  fHistoAngleMax = 0.5 ;
186  fHistoAngleMin = 0. ;
187 }
188 
189 //______________________________________________________________________________
193 //______________________________________________________________________________
195 {
196  Double_t max = fAngleMaxParam.At(0)*TMath::Exp(fAngleMaxParam.At(1)*e)
197  +fAngleMaxParam.At(2)+fAngleMaxParam.At(3)*e;
198  Double_t arg = (e*e-2*fM*fM)/(e*e);
199  Double_t min = 100. ;
200  if(arg>0.)
201  min = TMath::ACos(arg)+fShiftMinAngle[0]+fShiftMinAngle[1]*e;
202 
203  if((angle<max)&&(angle>=min)) return kTRUE ;
204  else return kFALSE ;
205 
206 }
207 
208 //_________________________________________________________________
214 //_________________________________________________________________
216  TLorentzVector gammaj,
217  Int_t calo)
218 {
219  fDecayBit = 0;
220 
221  // Double_t pt = (gammai+gammaj).Pt();
222  Double_t phi = (gammai+gammaj).Phi();
223  if(phi < 0)
224  phi+=TMath::TwoPi();
225 
226  Double_t invmass = (gammai+gammaj).M();
227  Double_t angle = gammaj.Angle(gammai.Vect());
228  Double_t e = (gammai+gammaj).E();
229  Double_t asy = TMath::Abs((gammai-gammaj).E())/(gammai+gammaj).E();
230 
231  //Fill histograms with no cuts applied.
233  {
234  fhAnglePairNoCut ->Fill(e,angle);
235  fhInvMassPairNoCut->Fill(e,invmass);
236  fhAsymmetryNoCut ->Fill(e,asy);
237  }
238 
239  //Cut on the aperture of the pair
240  if(fUseAngleCut)
241  {
242  if(IsAngleInWindow(angle,e))
243  {
245  {
246  fhAnglePairOpeningAngleCut ->Fill(e,angle);
247  fhInvMassPairOpeningAngleCut->Fill(e,invmass);
248  fhAsymmetryOpeningAngleCut ->Fill(e,asy);
249  }
250 
251  AliDebug(2,Form("Angle cut: energy %f, phi %f",e,phi));
252 
253  }
254  else return kFALSE;
255  }
256 
257  // Asymmetry cut
258  if(fUseAsymmetryCut)
259  {
260  if(fAsymmetryCut > asy)
261  {
263  {
264  fhInvMassPairAsymmetryCut->Fill(e,invmass);
265  fhAnglePairAsymmetryCut ->Fill(e,angle);
266  }
267  }
268  else
269  return kFALSE;
270  }
271 
272  //Cut on the invariant mass of the pair
273 
274  Float_t invmassmaxcut = fInvMassMaxCut;
275  Float_t invmassRightBandMinCut = fRightBandMinCut;
276  Float_t invmassRightBandMaxCut = fRightBandMaxCut;
277 
278  // kEMCAL=0, kPHOS=1
279  if(calo==0 && e > 6.)
280  {
281  // for EMCAL, pi0s, mass depends strongly with energy for e > 6, loose max cut
282 
284  invmassRightBandMinCut = (fInvMassMaxCutParam[0]+fRightBandMinCut)+fInvMassMaxCutParam[1]*e+fInvMassMaxCutParam[2]*e*e;
285  invmassRightBandMaxCut = (fInvMassMaxCutParam[0]+fRightBandMaxCut)+fInvMassMaxCutParam[1]*e+fInvMassMaxCutParam[2]*e*e;
286 
287  //printf("e %f, max cut %f, p00 %f,p0 %f,p1 %f,p2 %f\n",
288  // e,invmassmaxcut,fInvMassMaxCut,fInvMassMaxCutParam[0],fInvMassMaxCutParam[1],fInvMassMaxCutParam[2]);
289  }
290 
291  // normal case, invariant mass selection around pi0/eta peak
292  if( !fParticle.Contains("SideBand") )
293  {
294  if( invmass > fInvMassMinCut && invmass < invmassmaxcut )
295  {
296  if (fParticle=="Pi0") fDecayBit = kPi0;
297  else if(fParticle=="Eta") fDecayBit = kEta;
298  else AliWarning("Unexpected particle peak name!");
299 
301  {
302  fhInvMassPairAllCut->Fill(e,invmass);
303  fhAnglePairAllCut ->Fill(e,angle);
304  fhAsymmetryAllCut ->Fill(e,asy);
305  }
306 
307  AliDebug(2,Form("IM in peak cut: energy %f, im %f, bit map %d",e,invmass,fDecayBit));
308 
309  return kTRUE;
310 
311  }//(invmass>0.125) && (invmass<0.145)
312  else
313  {
314  return kFALSE;
315  }
316  }// normal selection
317 
318  else // select a band around pi0/eta
319  {
320  Bool_t ok = kFALSE;
321 
322  // Left side?
323  if( invmass > fLeftBandMinCut && invmass < fLeftBandMaxCut )
324  {
325  if (fParticle.Contains("Pi0")) fDecayBit = kPi0LeftSide;
326  else if(fParticle.Contains("Eta")) fDecayBit = kEtaLeftSide;
327  else AliWarning("Unexpected particle side_band name!");
328 
329  ok = kTRUE;
330  }
331 
332  //Right side?
333  if( invmass > invmassRightBandMinCut && invmass < invmassRightBandMaxCut )
334  {
335  if (fParticle.Contains("Pi0")) { fDecayBit = kPi0RightSide; if(ok) fDecayBit = kPi0BothSides; }
336  else if(fParticle.Contains("Eta")) { fDecayBit = kEtaRightSide; if(ok) fDecayBit = kEtaBothSides; }
337  else AliWarning("Unexpected particle side_band name!");
338 
339  ok = kTRUE;
340  }
341 
342  // Any of the sides?
343  if(ok)
344  {
346  {
347  fhInvMassPairAllCut->Fill(e,invmass);
348  fhAnglePairAllCut ->Fill(e,angle);
349  fhAsymmetryAllCut ->Fill(e,asy);
350  }
351 
352  AliDebug(2,Form("IM side band cut: energy %f, im %f, bit map %d",e,invmass,fDecayBit));
353 
354  return kTRUE;
355 
356  }//(invmass>0.125) && (invmass<0.145)
357  else
358  {
359  return kFALSE;
360  }
361  }
362 }
363 
364 //_______________________________________________________________
366 //_______________________________________________________________
368 {
369  fParticle = particleName ;
370 
371  if(particleName.Contains("Pi0"))
372  {
373  fHistoNIMBins = 150 ;
374  fHistoIMMax = 0.3 ;
375  fHistoIMMin = 0. ;
376 
377  fM = 0.135 ; // GeV
378  fInvMassMaxCut = 0.16 ; // GeV
379  fInvMassMinCut = 0.11 ; // GeV
380 
381  fLeftBandMinCut = 0.05 ; // GeV
382  fLeftBandMaxCut = 0.09 ; // GeV
383  fRightBandMinCut = 0.160 ; // GeV
384  fRightBandMaxCut = 0.165 ; // GeV
385 
386  fInvMassMaxCutParam[0] = 0.0 ;
387  fInvMassMaxCutParam[1] =-7.e-5 ;
388  fInvMassMaxCutParam[2] = 8.e-5 ;
389 
390  fShiftMinAngle[0] =-0.03 ;
391  fShiftMinAngle[1] = 0.0025;
392 
393  fAngleMaxParam.AddAt( 0.8, 0) ;
394  fAngleMaxParam.AddAt(-1, 1) ;
395  fAngleMaxParam.AddAt( 0.09, 2) ; //for pi0 shift, for eta maybe 0.09
396  fAngleMaxParam.AddAt(-2.e-3,3) ;
397  }
398  else if(particleName.Contains("Eta"))
399  {
400  fHistoNIMBins = 200 ; // GeV
401  fHistoIMMax = 0.75 ; // GeV
402  fHistoIMMin = 0.35 ; // GeV
403 
404  fM = 0.547 ; // GeV
405  fInvMassMaxCut = 0.590 ; // GeV
406  fInvMassMinCut = 0.510 ; // GeV
407 
408  fLeftBandMinCut = 0.450 ; // GeV
409  fLeftBandMaxCut = 0.500 ; // GeV
410  fRightBandMinCut = 0.600 ; // GeV
411  fRightBandMaxCut = 0.650 ; // GeV
412 
413  fInvMassMaxCutParam[0] = 0.00 ;
414  fInvMassMaxCutParam[1] = 0.00 ;
415  fInvMassMaxCutParam[2] = 0.00 ;
416 
417  fShiftMinAngle[0] =-0.03 ;
418  fShiftMinAngle[0] = 0.00 ;
419 
420  fAngleMaxParam.AddAt( 0.80, 0) ; // Same as pi0
421  fAngleMaxParam.AddAt(-0.25, 1) ; // Same as pi0
422  fAngleMaxParam.AddAt( 0.12, 2) ; // Shifted with respect to pi0
423  fAngleMaxParam.AddAt(-5.e-4, 3) ; // Same as pi0
424  }
425  else
426  printf("AliAnaNeutralMesonSelection::SetParticle(%s) *** Particle NOT defined (Pi0 or Eta), Pi0 settings by default *** \n",particleName.Data());
427 }
428 
429 //______________________________________________________________
431 //______________________________________________________________
433 {
434  if(! opt)
435  return;
436 
437  printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
438 
439  printf("Particle %s, bit %d, mass : %f \n", fParticle.Data(), fDecayBit, fM );
440  printf("Invariant mass limits : %f < m < %f \n", fInvMassMinCut , fInvMassMinCut );
441 
442  printf("Use asymmetry cut? : %d ; A < %f \n", fUseAngleCut, fAsymmetryCut );
443 
444  printf("Use angle cut? : %d \n", fUseAngleCut );
445  if(fUseAngleCut)
446  {
447  printf("Angle selection param: \n");
448  printf("p0 : %f\n", fAngleMaxParam.At(0));
449  printf("p1 : %f\n", fAngleMaxParam.At(1));
450  printf("p2 : %f\n", fAngleMaxParam.At(2));
451  printf("p3 : %f\n", fAngleMaxParam.At(3));
452  printf("Min angle shift : p0 = %1.3f, p1 = %1.3f,\n", fShiftMinAngle[0],fShiftMinAngle[1]);
453  }
454 
455  printf("Keep Neutral Meson Histos = %d\n",fKeepNeutralMesonHistos);
456 
458  {
459  printf("Histograms: %3.1f < E < %3.1f, Nbin = %d\n", fHistoEMin, fHistoEMax, fHistoNEBins);
460  printf("Histograms: %3.1f < angle < %3.1f, Nbin = %d\n", fHistoAngleMin, fHistoAngleMax, fHistoNAngleBins);
461  printf("Histograms: %3.1f < IM < %3.1f, Nbin = %d\n", fHistoIMMin, fHistoIMMax, fHistoNIMBins);
462  }
463 }
Float_t fHistoAngleMin
Minimum value of angle histogram range.
Float_t fHistoAngleMax
Maximum value of angle histogram range.
Double_t fLeftBandMaxCut
Side Band selection, max left band cut.
double Double_t
Definition: External.C:58
Float_t fAsymmetryCut
Asymmetry cut.
Definition: External.C:236
void InitParameters()
Initialize the parameters of the analysis.
TH2F * fhInvMassPairAllCut
! Invariant mass of decay photons, all cuts.
AliNeutralMesonSelection()
Default constructor. Initialize parameters.
Double_t fRightBandMinCut
Side Band selection, min right band cut.
void Print(const Option_t *opt) const
Print some relevant parameters set for the analysis.
Double_t fLeftBandMinCut
Side Band selection, min left band cut.
TH2F * fhAsymmetryOpeningAngleCut
! Asymmetry of decay photons, cut on opening angle.
TH2F * fhInvMassPairNoCut
! Invariant mass of decay photons, no cuts.
Bool_t fKeepNeutralMesonHistos
Keep neutral meson selection histograms.
Bool_t IsAngleInWindow(Float_t angle, Float_t e) const
Bool_t fUseAngleCut
Select pairs depending on their opening angle.
Double_t fRightBandMaxCut
Side Band selection, max right band cut.
int Int_t
Definition: External.C:63
float Float_t
Definition: External.C:68
Int_t fHistoNIMBins
Number of bins in Invariant Mass axis.
Float_t fHistoEMax
Maximum value of pi0 E histogram range.
Double_t fInvMassMaxCutParam[3]
Variable invariant mass max cut, for pi0 in EMCAL.
Double_t fInvMassMaxCut
Invariant Mass cut maximum.
TArrayD fAngleMaxParam
Maximum opening angle selection parameters.
void SetParticle(TString particleName)
Set some default parameters for selection of pi0 or eta.
Float_t fShiftMinAngle[2]
Correction shift for min angle from true kinematic limit, resolution effects.
Bool_t SelectPair(TLorentzVector particlei, TLorentzVector particlej, Int_t calo)
TH2F * fhAnglePairAllCut
! Aperture angle of decay photons, all cuts.
Double_t fM
Mass of the neutral meson.
Int_t fHistoNEBins
Number of bins in pi0 E axis.
TH2F * fhAnglePairNoCut
! Aperture angle of decay photons, no cuts.
TH2F * fhAsymmetryNoCut
! Asymmetry of decay photons, no cuts.
Float_t fHistoEMin
Minimum value of pi0 E histogram range.
Float_t fHistoIMMin
Minimum value of Invariant Mass histogram range.
TH2F * fhAnglePairOpeningAngleCut
! Aperture angle of decay photons, cut on opening angle.
UInt_t fDecayBit
Decay type flag, set while selecting, depending on fParticle and side range. See enum decayTypes for ...
ClassImp(AliAnalysisTaskCRC) AliAnalysisTaskCRC
Float_t fHistoIMMax
Maximum value of Invariant Mass histogram range.
Class that contains methods to select candidate cluster pairs to neutral meson.
const char Option_t
Definition: External.C:48
TH2F * fhInvMassPairAsymmetryCut
! Invariant mass of decay photons, asymmetry cut.
Bool_t fUseAsymmetryCut
Use the asymmetry cut.
bool Bool_t
Definition: External.C:53
Double_t fInvMassMinCut
Invariant Masscut minimun.
Int_t fHistoNAngleBins
Number of bins in angle axis.
TH2F * fhAsymmetryAllCut
! Asymmetry of decay photons, all cuts.
TH2F * fhAnglePairAsymmetryCut
! Aperture angle of decay photons, asymmetry cut.
TString fParticle
Meutral meson name (Pi0, Eta, +SideBand).
TH2F * fhInvMassPairOpeningAngleCut
! Invariant mass of decay photons, cut on opening angle.