Protection to skip PYTHIA events with large jet energy compared to pTHard
[u/mrichter/AliRoot.git] / PWG4 / PartCorrBase / AliNeutralMesonSelection.cxx
CommitLineData
1c5acb87 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/* $Id: AliNeutralMesonSelection.cxx 27413 2008-07-18 13:28:12Z gconesab $ */
16
17//_________________________________________________________________________
18// Class that contains methods to select candidate pairs to neutral meson
19// 2 main selections, invariant mass around pi0 (also any other mass),
20// apperture angle to distinguish from combinatorial.
21//-- Author: Gustavo Conesa (INFN-LNF)
22
23// --- ROOT system ---
24#include <TLorentzVector.h>
25#include <TH2.h>
26#include <TList.h>
1c5acb87 27
28//---- AliRoot system ----
29#include "AliNeutralMesonSelection.h"
1c5acb87 30
31ClassImp(AliNeutralMesonSelection)
32
33
34//____________________________________________________________________________
35 AliNeutralMesonSelection::AliNeutralMesonSelection() :
36 TObject(), fM(0),
37 fInvMassMaxCut(0.), fInvMassMinCut(0.),
38 fAngleMaxParam(), fKeepNeutralMesonHistos(0),
39 fhAnglePairNoCut(0), fhAnglePairOpeningAngleCut(0),
40 fhAnglePairAllCut(0),
41 fhInvMassPairNoCut(0), fhInvMassPairOpeningAngleCut(0),
42 fhInvMassPairAllCut(0),
477d6cee 43 fHistoNEBins(0), fHistoEMax(0.), fHistoEMin(0.),
44 fHistoNPtBins(0), fHistoPtMax(0.), fHistoPtMin(0.),
45 fHistoNAngleBins(0), fHistoAngleMax(0.), fHistoAngleMin(0.),
46 fHistoNIMBins(0), fHistoIMMax(0.), fHistoIMMin(0.)
47{
1c5acb87 48 //Default Ctor
49
50 //Initialize parameters
51
52 // kGammaHadron and kGammaJet
53 fAngleMaxParam.Set(4) ;
54 fAngleMaxParam.Reset(0.);
55
56 //Initialize parameters
57 InitParameters();
58}
59
60//____________________________________________________________________________
61AliNeutralMesonSelection::AliNeutralMesonSelection(const AliNeutralMesonSelection & g) :
62 TObject(), fM(g.fM),
63 fInvMassMaxCut(g.fInvMassMaxCut), fInvMassMinCut(g.fInvMassMinCut),
64 fAngleMaxParam(g.fAngleMaxParam),
65 fKeepNeutralMesonHistos(g.fKeepNeutralMesonHistos),
66 fhAnglePairNoCut(g. fhAnglePairNoCut),
67 fhAnglePairOpeningAngleCut(g. fhAnglePairOpeningAngleCut),
68 fhAnglePairAllCut(g. fhAnglePairAllCut),
69 fhInvMassPairNoCut(g.fhInvMassPairNoCut),
70 fhInvMassPairOpeningAngleCut(g.fhInvMassPairOpeningAngleCut),
71 fhInvMassPairAllCut(g.fhInvMassPairAllCut),
72 fHistoNEBins(g.fHistoNEBins), fHistoEMax(g.fHistoEMax), fHistoEMin(g.fHistoEMin),
73 fHistoNPtBins(g.fHistoNPtBins), fHistoPtMax(g.fHistoPtMax), fHistoPtMin(g.fHistoPtMin),
74 fHistoNAngleBins(g.fHistoNAngleBins), fHistoAngleMax(g.fHistoAngleMax), fHistoAngleMin(g.fHistoAngleMin),
75 fHistoNIMBins(g.fHistoNIMBins), fHistoIMMax(g.fHistoIMMax), fHistoIMMin(g.fHistoIMMin)
76{
77 // cpy ctor
78}
79
80//_________________________________________________________________________
81AliNeutralMesonSelection & AliNeutralMesonSelection::operator = (const AliNeutralMesonSelection & source)
82{
83 // assignment operator
84
85 if(this == &source)return *this;
86 ((TObject *)this)->operator=(source);
87
88 fM = source.fM ;
89 fInvMassMaxCut = source.fInvMassMaxCut ;
90 fInvMassMinCut = source.fInvMassMinCut ;
91 fAngleMaxParam = source.fAngleMaxParam ;
92 fKeepNeutralMesonHistos = source.fKeepNeutralMesonHistos;
93
94 fhAnglePairNoCut = source. fhAnglePairNoCut ;
95 fhAnglePairOpeningAngleCut = source. fhAnglePairOpeningAngleCut ;
96 fhAnglePairAllCut = source. fhAnglePairAllCut ;
97 fhInvMassPairNoCut = source.fhInvMassPairNoCut ;
98 fhInvMassPairOpeningAngleCut = source.fhInvMassPairOpeningAngleCut ;
99 fhInvMassPairAllCut = source.fhInvMassPairAllCut ;
100
101 fHistoNEBins = source.fHistoNEBins; fHistoEMax = source.fHistoEMax; fHistoEMin = source.fHistoEMin;
102 fHistoNPtBins = source.fHistoNPtBins; fHistoPtMax = source.fHistoPtMax; fHistoPtMin = source.fHistoPtMin;
103 fHistoNAngleBins = source.fHistoNAngleBins; fHistoAngleMax = source.fHistoAngleMax; fHistoAngleMin = source.fHistoAngleMin;
104 fHistoNIMBins = source.fHistoNIMBins; fHistoIMMax = source.fHistoIMMax; fHistoIMMin = source.fHistoIMMin;
105
106 return *this;
107
108}
109
110//____________________________________________________________________________
111AliNeutralMesonSelection::~AliNeutralMesonSelection()
112{
113 //dtor
114
115 if(!fKeepNeutralMesonHistos){
116 //Histograms initialized and filled but not passed to output container
117 //delete here, I am not sure this is correct
118
119 if(fhAnglePairNoCut) delete fhAnglePairNoCut;
120 if(fhAnglePairOpeningAngleCut) delete fhAnglePairOpeningAngleCut;
121 if(fhAnglePairAllCut) delete fhAnglePairAllCut;
122 if(fhInvMassPairNoCut) delete fhInvMassPairNoCut;
123 if(fhInvMassPairOpeningAngleCut) delete fhInvMassPairOpeningAngleCut;
124 if(fhInvMassPairAllCut) delete fhInvMassPairAllCut;
125
126 }
127
128}
129//________________________________________________________________________
130TList * AliNeutralMesonSelection::GetCreateOutputObjects()
131{
477d6cee 132 // Create histograms to be saved in output file and
133 // store them in outputContainer of the analysis class that calls this class.
134
135 TList * outputContainer = new TList() ;
136 outputContainer->SetName("MesonDecayHistos") ;
137
138 fhAnglePairNoCut = new TH2F
1c5acb87 139 ("AnglePairNoCut",
140 "Angle between all #gamma pair vs E_{#pi^{0}}",fHistoNEBins,fHistoEMin,fHistoEMax,fHistoNAngleBins,fHistoAngleMin,fHistoAngleMax);
477d6cee 141 fhAnglePairNoCut->SetYTitle("Angle (rad)");
142 fhAnglePairNoCut->SetXTitle("E_{ #pi^{0}} (GeV)");
143
144 fhAnglePairOpeningAngleCut = new TH2F
1c5acb87 145 ("AnglePairOpeningAngleCut",
146 "Angle between all #gamma pair (opening angle + azimuth cut) vs E_{#pi^{0}}"
147 ,fHistoNEBins,fHistoEMin,fHistoEMax,fHistoNAngleBins,fHistoAngleMin,fHistoAngleMax);
477d6cee 148 fhAnglePairOpeningAngleCut->SetYTitle("Angle (rad)");
149 fhAnglePairOpeningAngleCut->SetXTitle("E_{ #pi^{0}} (GeV)");
150
151 fhAnglePairAllCut = new TH2F
1c5acb87 152 ("AnglePairAllCut",
153 "Angle between all #gamma pair (opening angle + inv mass cut+azimuth) vs E_{#pi^{0}}"
154 ,fHistoNEBins,fHistoEMin,fHistoEMax,fHistoNAngleBins,fHistoAngleMin,fHistoAngleMax);
477d6cee 155 fhAnglePairAllCut->SetYTitle("Angle (rad)");
156 fhAnglePairAllCut->SetXTitle("E_{ #pi^{0}} (GeV)");
157
158 //
159 fhInvMassPairNoCut = new TH2F
1c5acb87 160 ("InvMassPairNoCut","Invariant Mass of all #gamma pair vs E_{#pi^{0}}",
161 fHistoNPtBins,fHistoPtMin,fHistoPtMax,fHistoNIMBins,fHistoIMMin,fHistoIMMax);
477d6cee 162 fhInvMassPairNoCut->SetYTitle("Invariant Mass (GeV/c^{2})");
163 fhInvMassPairNoCut->SetXTitle("E_{ #pi^{0}} (GeV)");
164
165 fhInvMassPairOpeningAngleCut = new TH2F
1c5acb87 166 ("InvMassPairOpeningAngleCut",
167 "Invariant Mass of #gamma pair (angle cut) vs E_{#pi^{0}}",
168 fHistoNPtBins,fHistoPtMin,fHistoPtMax,fHistoNIMBins,fHistoIMMin,fHistoIMMax);
477d6cee 169 fhInvMassPairOpeningAngleCut->SetYTitle("Invariant Mass (GeV/c^{2})");
170 fhInvMassPairOpeningAngleCut->SetXTitle(" E_{#pi^{0}}(GeV)");
171
172 fhInvMassPairAllCut = new TH2F
1c5acb87 173 ("InvMassPairAllCut",
174 "Invariant Mass of #gamma pair (opening angle+invmass cut) vs E_{#pi^{0}}",
175 fHistoNPtBins,fHistoPtMin,fHistoPtMax,fHistoNIMBins,fHistoIMMin,fHistoIMMax);
477d6cee 176 fhInvMassPairAllCut->SetYTitle("Invariant Mass (GeV/c^{2})");
177 fhInvMassPairAllCut->SetXTitle("E_{#pi^{0}}(GeV)");
178
179 outputContainer->Add(fhAnglePairNoCut) ;
180 outputContainer->Add(fhAnglePairOpeningAngleCut) ;
181 outputContainer->Add(fhAnglePairAllCut) ;
182
183 outputContainer->Add(fhInvMassPairNoCut) ;
184 outputContainer->Add(fhInvMassPairOpeningAngleCut) ;
185 outputContainer->Add(fhInvMassPairAllCut) ;
186
187 return outputContainer;
1c5acb87 188}
189
477d6cee 190//____________________________________________________________________________
1c5acb87 191void AliNeutralMesonSelection::InitParameters()
192{
477d6cee 193
1c5acb87 194 //Initialize the parameters of the analysis.
195 fKeepNeutralMesonHistos = kFALSE ;
477d6cee 196
1c5acb87 197 fAngleMaxParam.Set(4) ;
198 fAngleMaxParam.AddAt(0.4,0);//={0.4,-0.25,0.025,-2e-4};
199 fAngleMaxParam.AddAt(-0.25,1) ;
200 fAngleMaxParam.AddAt(0.025,2) ;
201 fAngleMaxParam.AddAt(-2e-4,3) ;
202
203 fInvMassMaxCut = 0.16 ;
204 fInvMassMinCut = 0.11 ;
205
206 fM = 0.1349766;//neutralMeson mass, pi0
207
208 //Histogrammes settings
209 fHistoNEBins = 100 ;
210 fHistoEMax = 50 ;
211 fHistoEMin = 0. ;
212
213 fHistoNPtBins = 240 ;
214 fHistoPtMax = 120 ;
215 fHistoPtMin = 0. ;
216
217 fHistoNAngleBins = 200 ;
218 fHistoAngleMax = 0.2;
219 fHistoAngleMin = 0. ;
220
221 fHistoNIMBins = 300 ;
222 fHistoIMMax = 0.5 ;
223 fHistoIMMin = 0. ;
224}
225
226//__________________________________________________________________________-
227Bool_t AliNeutralMesonSelection::IsAngleInWindow(const Float_t angle,const Float_t e) const {
228 //Check if the opening angle of the candidate pairs is inside
229 //our selection windowd
230
231 Bool_t result = kFALSE;
232 Double_t max = fAngleMaxParam.At(0)*TMath::Exp(fAngleMaxParam.At(1)*e)
233 +fAngleMaxParam.At(2)+fAngleMaxParam.At(3)*e;
234 Double_t arg = (e*e-2*fM*fM)/(e*e);
235 Double_t min = 100. ;
236 if(arg>0.)
237 min = TMath::ACos(arg);
238
239 if((angle<max)&&(angle>=min))
240 result = kTRUE;
241
242 return result;
243}
244
245//____________________________________________________________________________
246Bool_t AliNeutralMesonSelection::SelectPair(TLorentzVector gammai, TLorentzVector gammaj)
247{
248
249 //Search for the neutral pion within selection cuts
250 Bool_t goodpair = kFALSE ;
251
477d6cee 252// Double_t pt = (gammai+gammaj).Pt();
1c5acb87 253 Double_t phi = (gammai+gammaj).Phi();
254 if(phi < 0)
255 phi+=TMath::TwoPi();
256 Double_t invmass = (gammai+gammaj).M();
257 Double_t angle = gammaj.Angle(gammai.Vect());
258 Double_t e = (gammai+gammaj).E();
259
260 //Fill histograms with no cuts applied.
261 fhAnglePairNoCut->Fill(e,angle);
262 fhInvMassPairNoCut->Fill(e,invmass);
263
264 //Cut on the aperture of the pair
265 if(IsAngleInWindow(angle,e)){
266 fhAnglePairOpeningAngleCut ->Fill(e,angle);
267 fhInvMassPairOpeningAngleCut->Fill(e,invmass);
477d6cee 268 //AliDebug(2,Form("Angle cut: pt %f, phi %f",pt,phi));
1c5acb87 269
270 //Cut on the invariant mass of the pair
271 if((invmass>fInvMassMinCut) && (invmass<fInvMassMaxCut)){
272 fhInvMassPairAllCut ->Fill(e,invmass);
273 fhAnglePairAllCut ->Fill(e,angle);
274 goodpair = kTRUE;
477d6cee 275 //AliDebug(2,Form("IM cut: pt %f, phi %f",pt,phi));
1c5acb87 276 }//(invmass>0.125) && (invmass<0.145)
277 }//Opening angle cut
278
279 return goodpair;
280
281}
282
283//__________________________________________________________________
284void AliNeutralMesonSelection::Print(const Option_t * opt) const
285{
286
287 //Print some relevant parameters set for the analysis
288 if(! opt)
289 return;
290
a3aebfff 291 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
1c5acb87 292
293 printf("mass : %f \n", fM );
294 printf("Invariant mass limits : %f < m < %f \n", fInvMassMinCut , fInvMassMinCut );
295 printf("Angle selection param: \n");
fbeaf916 296 printf("p0 : %f\n", fAngleMaxParam.At(0));
297 printf("p1 : %f\n", fAngleMaxParam.At(1));
298 printf("p2 : %f\n", fAngleMaxParam.At(2));
299 printf("p3 : %f\n", fAngleMaxParam.At(3));
1c5acb87 300
301 printf("Keep Neutral Meson Histos = %d\n",fKeepNeutralMesonHistos);
302
303 if(fKeepNeutralMesonHistos){
477d6cee 304 printf("Histograms: %3.1f < E < %3.1f, Nbin = %d\n", fHistoEMin, fHistoEMax, fHistoNEBins);
305 printf("Histograms: %3.1f < pT < %3.1f, Nbin = %d\n", fHistoPtMin, fHistoPtMax, fHistoNPtBins);
306 printf("Histograms: %3.1f < angle < %3.1f, Nbin = %d\n", fHistoAngleMin, fHistoAngleMax, fHistoNAngleBins);
307 printf("Histograms: %3.1f < IM < %3.1f, Nbin = %d\n", fHistoIMMin, fHistoIMMax, fHistoNIMBins);
308
1c5acb87 309 }
477d6cee 310
1c5acb87 311}