AliPhysics  v5-06-11-01 (156c7f3)
 All Classes Namespaces Files Functions Variables Enumerations Enumerator Macros
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 // histograms
42 fhAnglePairNoCut(0), fhAnglePairOpeningAngleCut(0),
43 fhAnglePairAsymmetryCut(0), fhAnglePairAllCut(0),
44 fhInvMassPairNoCut(0), fhInvMassPairOpeningAngleCut(0),
45 fhInvMassPairAsymmetryCut(0), fhInvMassPairAllCut(0),
46 fhAsymmetryNoCut(0), fhAsymmetryOpeningAngleCut(0),
47 fhAsymmetryAllCut(0),
48 // histogram ranges and bins
49 fHistoNEBins(0), fHistoEMax(0.), fHistoEMin(0.),
50 fHistoNAngleBins(0), fHistoAngleMax(0.), fHistoAngleMin(0.),
51 fHistoNIMBins(0), fHistoIMMax(0.), fHistoIMMin(0.)
52 {
54 }
55 
56 //_________________________________________________________
59 //_________________________________________________________
61 {
62  TList * outputContainer = new TList() ;
63  outputContainer->SetName("MesonDecayHistos") ;
64 
66 
67  outputContainer->SetOwner(kFALSE);
68 
69  fhAnglePairNoCut = new TH2F
70  ("AnglePairNoCut",
71  "Angle between all #gamma pair vs E_{#pi^{0}}",fHistoNEBins,fHistoEMin,fHistoEMax,fHistoNAngleBins,fHistoAngleMin,fHistoAngleMax);
72  fhAnglePairNoCut->SetYTitle("Angle (rad)");
73  fhAnglePairNoCut->SetXTitle("E_{ #pi^{0}} (GeV)");
74 
75  fhAsymmetryNoCut = new TH2F
76  ("AsymmetryNoCut","Asymmetry of all #gamma pair vs E_{#pi^{0}}",
78  fhAsymmetryNoCut->SetYTitle("Asymmetry");
79  fhAsymmetryNoCut->SetXTitle("E_{ #pi^{0}} (GeV)");
80 
81  fhInvMassPairNoCut = new TH2F
82  ("InvMassPairNoCut","Invariant Mass of all #gamma pair vs E_{#pi^{0}}",
84  fhInvMassPairNoCut->SetYTitle("Invariant Mass (GeV/c^{2})");
85  fhInvMassPairNoCut->SetXTitle("E_{ #pi^{0}} (GeV)");
86 
87  outputContainer->Add(fhAnglePairNoCut) ;
88  outputContainer->Add(fhAsymmetryNoCut) ;
89  outputContainer->Add(fhInvMassPairNoCut) ;
90 
91  if(fUseAngleCut) {
93  ("AnglePairOpeningAngleCut",
94  "Angle between all #gamma pair (opening angle) vs E_{#pi^{0}}"
96  fhAnglePairOpeningAngleCut->SetYTitle("Angle (rad)");
97  fhAnglePairOpeningAngleCut->SetXTitle("E_{ #pi^{0}} (GeV)");
98 
100  ("AsymmetryOpeningAngleCut",
101  "Asymmetry of #gamma pair (angle cut) vs E_{#pi^{0}}",
103  fhAsymmetryOpeningAngleCut->SetYTitle("Asymmetry");
104  fhAsymmetryOpeningAngleCut->SetXTitle(" E_{#pi^{0}}(GeV)");
105 
107  ("InvMassPairOpeningAngleCut",
108  "Invariant Mass of #gamma pair (angle cut) vs E_{#pi^{0}}",
110  fhInvMassPairOpeningAngleCut->SetYTitle("Invariant Mass (GeV/c^{2})");
111  fhInvMassPairOpeningAngleCut->SetXTitle(" E_{#pi^{0}}(GeV)");
112 
113  outputContainer->Add(fhAnglePairOpeningAngleCut) ;
114  outputContainer->Add(fhAsymmetryOpeningAngleCut) ;
115  outputContainer->Add(fhInvMassPairOpeningAngleCut) ;
116  }
117 
118  if(fUseAsymmetryCut) {
119  fhAnglePairAsymmetryCut = new TH2F
120  ("AnglePairAsymmetryCut",
121  "Angle between all #gamma pair (opening angle + asymetry cut) vs E_{#pi^{0}}"
123  fhAnglePairAsymmetryCut->SetYTitle("Angle (rad)");
124  fhAnglePairAsymmetryCut->SetXTitle("E_{ #pi^{0}} (GeV)");
125 
126  fhInvMassPairAsymmetryCut = new TH2F
127  ("InvMassPairAsymmetryCut",
128  "Invariant Mass of #gamma pair (opening angle + asymmetry) vs E_{#pi^{0}}",
130  fhInvMassPairAsymmetryCut->SetYTitle("Invariant Mass (GeV/c^{2})");
131  fhInvMassPairAsymmetryCut->SetXTitle("E_{#pi^{0}}(GeV)");
132 
133  outputContainer->Add(fhAnglePairAsymmetryCut) ;
134  outputContainer->Add(fhInvMassPairAsymmetryCut) ;
135 
136  }
137 
138  fhAnglePairAllCut = new TH2F
139  ("AnglePairAllCut",
140  "Angle between all #gamma pair (opening angle + asymmetry + inv mass cut) vs E_{#pi^{0}}"
142  fhAnglePairAllCut->SetYTitle("Angle (rad)");
143  fhAnglePairAllCut->SetXTitle("E_{ #pi^{0}} (GeV)");
144 
145  fhInvMassPairAllCut = new TH2F
146  ("InvMassPairAllCut",
147  "Invariant Mass of #gamma pair (opening angle + asymmetry + invmass cut) vs E_{#pi^{0}}",
149  fhInvMassPairAllCut->SetYTitle("Invariant Mass (GeV/c^{2})");
150  fhInvMassPairAllCut->SetXTitle("E_{#pi^{0}}(GeV)");
151 
152  fhAsymmetryAllCut = new TH2F
153  ("AsymmetryAllCut",
154  "Asymmetry of #gamma pair (opening angle+invmass cut) vs E_{#pi^{0}}",
156  fhAsymmetryAllCut->SetYTitle("Asymmetry");
157  fhAsymmetryAllCut->SetXTitle("E_{#pi^{0}}(GeV)");
158 
159  outputContainer->Add(fhAnglePairAllCut) ;
160  outputContainer->Add(fhAsymmetryAllCut) ;
161  outputContainer->Add(fhInvMassPairAllCut) ;
162 
163  }
164 
165  return outputContainer;
166 }
167 
168 //_____________________________________________
170 //_____________________________________________
172 {
173  fAngleMaxParam.Set(4) ;
174  fAngleMaxParam.Reset(0.);
175 
176  SetParticle("Pi0");
177 
178  //Histogrammes settings
179  fHistoNEBins = 200 ;
180  fHistoEMax = 50 ;
181  fHistoEMin = 0. ;
182 
183  fHistoNAngleBins = 200 ;
184  fHistoAngleMax = 0.5 ;
185  fHistoAngleMin = 0. ;
186 }
187 
188 //______________________________________________________________________________
192 //______________________________________________________________________________
193 Bool_t AliNeutralMesonSelection::IsAngleInWindow(Float_t angle, Float_t e) const
194 {
195  Double_t max = fAngleMaxParam.At(0)*TMath::Exp(fAngleMaxParam.At(1)*e)
196  +fAngleMaxParam.At(2)+fAngleMaxParam.At(3)*e;
197  Double_t arg = (e*e-2*fM*fM)/(e*e);
198  Double_t min = 100. ;
199  if(arg>0.)
200  min = TMath::ACos(arg)+fShiftMinAngle[0]+fShiftMinAngle[1]*e;
201 
202  if((angle<max)&&(angle>=min)) return kTRUE ;
203  else return kFALSE ;
204 
205 }
206 
207 //_________________________________________________________________
213 //_________________________________________________________________
214 Bool_t AliNeutralMesonSelection::SelectPair(TLorentzVector gammai,
215  TLorentzVector gammaj,
216  Int_t calo)
217 {
218  fDecayBit = 0;
219 
220  // Double_t pt = (gammai+gammaj).Pt();
221  Double_t phi = (gammai+gammaj).Phi();
222  if(phi < 0)
223  phi+=TMath::TwoPi();
224 
225  Double_t invmass = (gammai+gammaj).M();
226  Double_t angle = gammaj.Angle(gammai.Vect());
227  Double_t e = (gammai+gammaj).E();
228  Double_t asy = TMath::Abs((gammai-gammaj).E())/(gammai+gammaj).E();
229 
230  //Fill histograms with no cuts applied.
232  {
233  fhAnglePairNoCut ->Fill(e,angle);
234  fhInvMassPairNoCut->Fill(e,invmass);
235  fhAsymmetryNoCut ->Fill(e,asy);
236  }
237 
238  //Cut on the aperture of the pair
239  if(fUseAngleCut)
240  {
241  if(IsAngleInWindow(angle,e))
242  {
244  {
245  fhAnglePairOpeningAngleCut ->Fill(e,angle);
246  fhInvMassPairOpeningAngleCut->Fill(e,invmass);
247  fhAsymmetryOpeningAngleCut ->Fill(e,asy);
248  }
249 
250  AliDebug(2,Form("Angle cut: energy %f, phi %f",e,phi));
251 
252  }
253  else return kFALSE;
254  }
255 
256  // Asymmetry cut
257  if(fUseAsymmetryCut)
258  {
259  if(fAsymmetryCut > asy)
260  {
262  {
263  fhInvMassPairAsymmetryCut->Fill(e,invmass);
264  fhAnglePairAsymmetryCut ->Fill(e,angle);
265  }
266  }
267  else
268  return kFALSE;
269  }
270 
271  //Cut on the invariant mass of the pair
272 
273  Float_t invmassmaxcut = fInvMassMaxCut;
274  Float_t invmassRightBandMinCut = fRightBandMinCut;
275  Float_t invmassRightBandMaxCut = fRightBandMaxCut;
276 
277  // kEMCAL=0, kPHOS=1
278  if(calo==0 && e > 6.)
279  {
280  // for EMCAL, pi0s, mass depends strongly with energy for e > 6, loose max cut
281 
283  invmassRightBandMinCut = (fInvMassMaxCutParam[0]+fRightBandMinCut)+fInvMassMaxCutParam[1]*e+fInvMassMaxCutParam[2]*e*e;
284  invmassRightBandMaxCut = (fInvMassMaxCutParam[0]+fRightBandMaxCut)+fInvMassMaxCutParam[1]*e+fInvMassMaxCutParam[2]*e*e;
285 
286  //printf("e %f, max cut %f, p00 %f,p0 %f,p1 %f,p2 %f\n",
287  // e,invmassmaxcut,fInvMassMaxCut,fInvMassMaxCutParam[0],fInvMassMaxCutParam[1],fInvMassMaxCutParam[2]);
288  }
289 
290  // normal case, invariant mass selection around pi0/eta peak
291  Int_t decayBit = -1;
292 
293  if( !fParticle.Contains("SideBand") )
294  {
295  if( invmass > fInvMassMinCut && invmass < invmassmaxcut )
296  {
297  if (fParticle=="Pi0") fDecayBit = kPi0;
298  else if(fParticle=="Eta") fDecayBit = kEta;
299  else AliWarning("Unexpected particle peak name!");
300 
302  {
303  fhInvMassPairAllCut->Fill(e,invmass);
304  fhAnglePairAllCut ->Fill(e,angle);
305  fhAsymmetryAllCut ->Fill(e,asy);
306  }
307 
308  AliDebug(2,Form("IM in peak cut: energy %f, im %f, bit map %d",e,invmass,fDecayBit));
309 
310  return kTRUE;
311 
312  }//(invmass>0.125) && (invmass<0.145)
313  else
314  {
315  return kFALSE;
316  }
317  }// normal selection
318 
319  else // select a band around pi0/eta
320  {
321  Bool_t ok = kFALSE;
322 
323  // Left side?
324  if( invmass > fLeftBandMinCut && invmass < fLeftBandMaxCut )
325  {
326  if (fParticle.Contains("Pi0")) fDecayBit = kPi0LeftSide;
327  else if(fParticle.Contains("Eta")) fDecayBit = kEtaLeftSide;
328  else AliWarning("Unexpected particle side_band name!");
329 
330  ok = kTRUE;
331  }
332 
333  //Right side?
334  if( invmass > invmassRightBandMinCut && invmass < invmassRightBandMaxCut )
335  {
336  if (fParticle.Contains("Pi0")) { fDecayBit = kPi0RightSide; if(ok) fDecayBit = kPi0BothSides; }
337  else if(fParticle.Contains("Eta")) { fDecayBit = kEtaRightSide; if(ok) fDecayBit = kEtaBothSides; }
338  else AliWarning("Unexpected particle side_band name!");
339 
340  ok = kTRUE;
341  }
342 
343  // Any of the sides?
344  if(ok)
345  {
347  {
348  fhInvMassPairAllCut->Fill(e,invmass);
349  fhAnglePairAllCut ->Fill(e,angle);
350  fhAsymmetryAllCut ->Fill(e,asy);
351  }
352 
353  AliDebug(2,Form("IM side band cut: energy %f, im %f, bit map %d",e,invmass,fDecayBit));
354 
355  return kTRUE;
356 
357  }//(invmass>0.125) && (invmass<0.145)
358  else
359  {
360  return kFALSE;
361  }
362  }
363 }
364 
365 //_______________________________________________________________
367 //_______________________________________________________________
368 void AliNeutralMesonSelection::SetParticle(TString particleName)
369 {
370  fParticle = particleName ;
371 
372  if(particleName.Contains("Pi0"))
373  {
374  fHistoNIMBins = 150 ;
375  fHistoIMMax = 0.3 ;
376  fHistoIMMin = 0. ;
377 
378  fM = 0.135 ; // GeV
379  fInvMassMaxCut = 0.16 ; // GeV
380  fInvMassMinCut = 0.11 ; // GeV
381 
382  fLeftBandMinCut = 0.05 ; // GeV
383  fLeftBandMaxCut = 0.09 ; // GeV
384  fRightBandMinCut = 0.160 ; // GeV
385  fRightBandMaxCut = 0.165 ; // GeV
386 
387  fInvMassMaxCutParam[0] = 0.0 ;
388  fInvMassMaxCutParam[1] =-7.e-5 ;
389  fInvMassMaxCutParam[2] = 8.e-5 ;
390 
391  fShiftMinAngle[0] =-0.03 ;
392  fShiftMinAngle[1] = 0.0025;
393 
394  fAngleMaxParam.AddAt( 0.8, 0) ;
395  fAngleMaxParam.AddAt(-1, 1) ;
396  fAngleMaxParam.AddAt( 0.09, 2) ; //for pi0 shift, for eta maybe 0.09
397  fAngleMaxParam.AddAt(-2.e-3,3) ;
398  }
399  else if(particleName.Contains("Eta"))
400  {
401  fHistoNIMBins = 200 ; // GeV
402  fHistoIMMax = 0.75 ; // GeV
403  fHistoIMMin = 0.35 ; // GeV
404 
405  fM = 0.547 ; // GeV
406  fInvMassMaxCut = 0.590 ; // GeV
407  fInvMassMinCut = 0.510 ; // GeV
408 
409  fLeftBandMinCut = 0.450 ; // GeV
410  fLeftBandMaxCut = 0.500 ; // GeV
411  fRightBandMinCut = 0.600 ; // GeV
412  fRightBandMaxCut = 0.650 ; // GeV
413 
414  fInvMassMaxCutParam[0] = 0.00 ;
415  fInvMassMaxCutParam[1] = 0.00 ;
416  fInvMassMaxCutParam[2] = 0.00 ;
417 
418  fShiftMinAngle[0] =-0.03 ;
419  fShiftMinAngle[0] = 0.00 ;
420 
421  fAngleMaxParam.AddAt( 0.80, 0) ; // Same as pi0
422  fAngleMaxParam.AddAt(-0.25, 1) ; // Same as pi0
423  fAngleMaxParam.AddAt( 0.12, 2) ; // Shifted with respect to pi0
424  fAngleMaxParam.AddAt(-5.e-4, 3) ; // Same as pi0
425  }
426  else
427  printf("AliAnaNeutralMesonSelection::SetParticle(%s) *** Particle NOT defined (Pi0 or Eta), Pi0 settings by default *** \n",particleName.Data());
428 }
429 
430 //______________________________________________________________
432 //______________________________________________________________
433 void AliNeutralMesonSelection::Print(const Option_t * opt) const
434 {
435  if(! opt)
436  return;
437 
438  printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
439 
440  printf("Particle %s, bit %d, mass : %f \n", fParticle.Data(), fDecayBit, fM );
441  printf("Invariant mass limits : %f < m < %f \n", fInvMassMinCut , fInvMassMinCut );
442 
443  printf("Use asymmetry cut? : %d ; A < %f \n", fUseAngleCut, fAsymmetryCut );
444 
445  printf("Use angle cut? : %d \n", fUseAngleCut );
446  if(fUseAngleCut)
447  {
448  printf("Angle selection param: \n");
449  printf("p0 : %f\n", fAngleMaxParam.At(0));
450  printf("p1 : %f\n", fAngleMaxParam.At(1));
451  printf("p2 : %f\n", fAngleMaxParam.At(2));
452  printf("p3 : %f\n", fAngleMaxParam.At(3));
453  printf("Min angle shift : p0 = %1.3f, p1 = %1.3f,\n", fShiftMinAngle[0],fShiftMinAngle[1]);
454  }
455 
456  printf("Keep Neutral Meson Histos = %d\n",fKeepNeutralMesonHistos);
457 
459  {
460  printf("Histograms: %3.1f < E < %3.1f, Nbin = %d\n", fHistoEMin, fHistoEMax, fHistoNEBins);
461  printf("Histograms: %3.1f < angle < %3.1f, Nbin = %d\n", fHistoAngleMin, fHistoAngleMax, fHistoNAngleBins);
462  printf("Histograms: %3.1f < IM < %3.1f, Nbin = %d\n", fHistoIMMin, fHistoIMMax, fHistoNIMBins);
463  }
464 }
Float_t fHistoAngleMin
Minimum value of angle histogram range.
Float_t fHistoAngleMax
Maximum value of angle histogram range.
ClassImp(AliAnalysisTaskTriggerRates) AliAnalysisTaskTriggerRates
Double_t fLeftBandMaxCut
Side Band selection, max left band cut.
Float_t fAsymmetryCut
Asymmetry cut.
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_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 ...
Float_t fHistoIMMax
Maximum value of Invariant Mass histogram range.
Class that contains methods to select candidate cluster pairs to neutral meson.
TH2F * fhInvMassPairAsymmetryCut
! Invariant mass of decay photons, asymmetry cut.
Bool_t fUseAsymmetryCut
Use the asymmetry cut.
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.