making the container name parameterized
[u/mrichter/AliRoot.git] / PWG / CaloTrackCorrBase / 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 **************************************************************************/
1c5acb87 15
16//_________________________________________________________________________
17// Class that contains methods to select candidate pairs to neutral meson
18// 2 main selections, invariant mass around pi0 (also any other mass),
19// apperture angle to distinguish from combinatorial.
20//-- Author: Gustavo Conesa (INFN-LNF)
21
22// --- ROOT system ---
23#include <TLorentzVector.h>
24#include <TH2.h>
25#include <TList.h>
3f150b4b 26#include <TString.h>
27
1c5acb87 28//---- AliRoot system ----
29#include "AliNeutralMesonSelection.h"
1c5acb87 30
31ClassImp(AliNeutralMesonSelection)
32
33
b960c7eb 34//______________________________________________________
1c5acb87 35 AliNeutralMesonSelection::AliNeutralMesonSelection() :
3f150b4b 36 TObject(), fAsymmetryCut(1),
37 fUseAsymmetryCut(0), fM(0),
38 fInvMassMaxCut(0.), fInvMassMinCut(0.),
39 fLeftBandMinCut(0.), fLeftBandMaxCut(0.),
40 fRightBandMinCut(0.), fRightBandMaxCut(0.),
41 fAngleMaxParam(), fUseAngleCut(0),
42 fKeepNeutralMesonHistos(0), fParticle(""),
43 // histograms
44 fhAnglePairNoCut(0), fhAnglePairOpeningAngleCut(0),
45 fhAnglePairAsymmetryCut(0), fhAnglePairAllCut(0),
46 fhInvMassPairNoCut(0), fhInvMassPairOpeningAngleCut(0),
47 fhInvMassPairAsymmetryCut(0), fhInvMassPairAllCut(0),
48 fhAsymmetryNoCut(0), fhAsymmetryOpeningAngleCut(0),
49 fhAsymmetryAllCut(0),
50 // histogram ranges and bins
9542ac75 51 fHistoNEBins(0), fHistoEMax(0.), fHistoEMin(0.),
9542ac75 52 fHistoNAngleBins(0), fHistoAngleMax(0.), fHistoAngleMin(0.),
53 fHistoNIMBins(0), fHistoIMMax(0.), fHistoIMMin(0.)
477d6cee 54{
1c5acb87 55 //Default Ctor
56
57 //Initialize parameters
1c5acb87 58 InitParameters();
59}
1c5acb87 60
b960c7eb 61//_________________________________________________________
1c5acb87 62TList * AliNeutralMesonSelection::GetCreateOutputObjects()
63{
477d6cee 64 // Create histograms to be saved in output file and
65 // store them in outputContainer of the analysis class that calls this class.
66
67 TList * outputContainer = new TList() ;
68 outputContainer->SetName("MesonDecayHistos") ;
9542ac75 69
fad7739b 70 if(fKeepNeutralMesonHistos){
71
72 outputContainer->SetOwner(kFALSE);
73
74 fhAnglePairNoCut = new TH2F
75 ("AnglePairNoCut",
76 "Angle between all #gamma pair vs E_{#pi^{0}}",fHistoNEBins,fHistoEMin,fHistoEMax,fHistoNAngleBins,fHistoAngleMin,fHistoAngleMax);
77 fhAnglePairNoCut->SetYTitle("Angle (rad)");
78 fhAnglePairNoCut->SetXTitle("E_{ #pi^{0}} (GeV)");
9542ac75 79
80 fhAsymmetryNoCut = new TH2F
81 ("AsymmetryNoCut","Asymmetry of all #gamma pair vs E_{#pi^{0}}",
3bfcb597 82 fHistoNEBins,fHistoEMin,fHistoEMax,100,0,1);
9542ac75 83 fhAsymmetryNoCut->SetYTitle("Asymmetry");
84 fhAsymmetryNoCut->SetXTitle("E_{ #pi^{0}} (GeV)");
85
86 fhInvMassPairNoCut = new TH2F
87 ("InvMassPairNoCut","Invariant Mass of all #gamma pair vs E_{#pi^{0}}",
3bfcb597 88 fHistoNEBins,fHistoEMin,fHistoEMax,fHistoNIMBins,fHistoIMMin,fHistoIMMax);
9542ac75 89 fhInvMassPairNoCut->SetYTitle("Invariant Mass (GeV/c^{2})");
90 fhInvMassPairNoCut->SetXTitle("E_{ #pi^{0}} (GeV)");
91
92 outputContainer->Add(fhAnglePairNoCut) ;
93 outputContainer->Add(fhAsymmetryNoCut) ;
94 outputContainer->Add(fhInvMassPairNoCut) ;
95
96 if(fUseAngleCut) {
97 fhAnglePairOpeningAngleCut = new TH2F
98 ("AnglePairOpeningAngleCut",
99 "Angle between all #gamma pair (opening angle) vs E_{#pi^{0}}"
100 ,fHistoNEBins,fHistoEMin,fHistoEMax,fHistoNAngleBins,fHistoAngleMin,fHistoAngleMax);
101 fhAnglePairOpeningAngleCut->SetYTitle("Angle (rad)");
102 fhAnglePairOpeningAngleCut->SetXTitle("E_{ #pi^{0}} (GeV)");
103
104 fhAsymmetryOpeningAngleCut = new TH2F
105 ("AsymmetryOpeningAngleCut",
106 "Asymmetry of #gamma pair (angle cut) vs E_{#pi^{0}}",
3bfcb597 107 fHistoNEBins,fHistoEMin,fHistoEMax,100,0,1);
9542ac75 108 fhAsymmetryOpeningAngleCut->SetYTitle("Asymmetry");
109 fhAsymmetryOpeningAngleCut->SetXTitle(" E_{#pi^{0}}(GeV)");
110
111 fhInvMassPairOpeningAngleCut = new TH2F
112 ("InvMassPairOpeningAngleCut",
113 "Invariant Mass of #gamma pair (angle cut) vs E_{#pi^{0}}",
3bfcb597 114 fHistoNEBins,fHistoEMin,fHistoEMax,fHistoNIMBins,fHistoIMMin,fHistoIMMax);
9542ac75 115 fhInvMassPairOpeningAngleCut->SetYTitle("Invariant Mass (GeV/c^{2})");
116 fhInvMassPairOpeningAngleCut->SetXTitle(" E_{#pi^{0}}(GeV)");
117
118 outputContainer->Add(fhAnglePairOpeningAngleCut) ;
119 outputContainer->Add(fhAsymmetryOpeningAngleCut) ;
120 outputContainer->Add(fhInvMassPairOpeningAngleCut) ;
121 }
122
123 if(fUseAsymmetryCut) {
124 fhAnglePairAsymmetryCut = new TH2F
125 ("AnglePairAsymmetryCut",
126 "Angle between all #gamma pair (opening angle + asymetry cut) vs E_{#pi^{0}}"
127 ,fHistoNEBins,fHistoEMin,fHistoEMax,fHistoNAngleBins,fHistoAngleMin,fHistoAngleMax);
128 fhAnglePairAsymmetryCut->SetYTitle("Angle (rad)");
129 fhAnglePairAsymmetryCut->SetXTitle("E_{ #pi^{0}} (GeV)");
130
131 fhInvMassPairAsymmetryCut = new TH2F
132 ("InvMassPairAsymmetryCut",
133 "Invariant Mass of #gamma pair (opening angle + asymmetry) vs E_{#pi^{0}}",
3bfcb597 134 fHistoNEBins,fHistoEMin,fHistoEMax,fHistoNIMBins,fHistoIMMin,fHistoIMMax);
9542ac75 135 fhInvMassPairAsymmetryCut->SetYTitle("Invariant Mass (GeV/c^{2})");
136 fhInvMassPairAsymmetryCut->SetXTitle("E_{#pi^{0}}(GeV)");
137
138 outputContainer->Add(fhAnglePairAsymmetryCut) ;
139 outputContainer->Add(fhInvMassPairAsymmetryCut) ;
140
141 }
142
fad7739b 143 fhAnglePairAllCut = new TH2F
144 ("AnglePairAllCut",
9542ac75 145 "Angle between all #gamma pair (opening angle + asymmetry + inv mass cut) vs E_{#pi^{0}}"
fad7739b 146 ,fHistoNEBins,fHistoEMin,fHistoEMax,fHistoNAngleBins,fHistoAngleMin,fHistoAngleMax);
147 fhAnglePairAllCut->SetYTitle("Angle (rad)");
148 fhAnglePairAllCut->SetXTitle("E_{ #pi^{0}} (GeV)");
9542ac75 149
fad7739b 150 fhInvMassPairAllCut = new TH2F
151 ("InvMassPairAllCut",
9542ac75 152 "Invariant Mass of #gamma pair (opening angle + asymmetry + invmass cut) vs E_{#pi^{0}}",
3bfcb597 153 fHistoNEBins,fHistoEMin,fHistoEMax,fHistoNIMBins,fHistoIMMin,fHistoIMMax);
fad7739b 154 fhInvMassPairAllCut->SetYTitle("Invariant Mass (GeV/c^{2})");
155 fhInvMassPairAllCut->SetXTitle("E_{#pi^{0}}(GeV)");
156
9542ac75 157 fhAsymmetryAllCut = new TH2F
158 ("AsymmetryAllCut",
159 "Asymmetry of #gamma pair (opening angle+invmass cut) vs E_{#pi^{0}}",
3bfcb597 160 fHistoNEBins,fHistoEMin,fHistoEMax,100,0,1);
9542ac75 161 fhAsymmetryAllCut->SetYTitle("Asymmetry");
162 fhAsymmetryAllCut->SetXTitle("E_{#pi^{0}}(GeV)");
163
164 outputContainer->Add(fhAnglePairAllCut) ;
165 outputContainer->Add(fhAsymmetryAllCut) ;
166 outputContainer->Add(fhInvMassPairAllCut) ;
167
fad7739b 168 }
9542ac75 169
477d6cee 170 return outputContainer;
9542ac75 171
1c5acb87 172}
173
b960c7eb 174//_____________________________________________
1c5acb87 175void AliNeutralMesonSelection::InitParameters()
176{
477d6cee 177
1c5acb87 178 //Initialize the parameters of the analysis.
1c5acb87 179 fAngleMaxParam.Set(4) ;
9542ac75 180 fAngleMaxParam.Reset(0.);
1c5acb87 181
3bfcb597 182 SetParticle("Pi0");
1c5acb87 183
3bfcb597 184 //Histogrammes settings
185 fHistoNEBins = 200 ;
186 fHistoEMax = 50 ;
a5fb4114 187 fHistoEMin = 0. ;
1c5acb87 188
3bfcb597 189 fHistoNAngleBins = 200 ;
190 fHistoAngleMax = 0.5 ;
1c5acb87 191 fHistoAngleMin = 0. ;
192
1c5acb87 193}
194
b960c7eb 195//_____________________________________________________________________
196Bool_t AliNeutralMesonSelection::IsAngleInWindow(const Float_t angle,
197 const Float_t e) const
a5fb4114 198{
199
200 // Check if the opening angle of the candidate pairs is inside
201 // our selection window
3bfcb597 202 // Attention, only valid for Pi0, if needed for Eta need to revise max angle function or change parameters
9542ac75 203
1c5acb87 204 Double_t max = fAngleMaxParam.At(0)*TMath::Exp(fAngleMaxParam.At(1)*e)
205 +fAngleMaxParam.At(2)+fAngleMaxParam.At(3)*e;
206 Double_t arg = (e*e-2*fM*fM)/(e*e);
207 Double_t min = 100. ;
208 if(arg>0.)
7287751b 209 min = TMath::ACos(arg)+fShiftMinAngle[0]+fShiftMinAngle[1]*e;
3bfcb597 210
9542ac75 211 if((angle<max)&&(angle>=min)) return kTRUE ;
212 else return kFALSE ;
213
1c5acb87 214}
215
b960c7eb 216//_________________________________________________________________
217Bool_t AliNeutralMesonSelection::SelectPair(TLorentzVector gammai,
218 TLorentzVector gammaj,
219 TString calo)
1c5acb87 220{
221
222 //Search for the neutral pion within selection cuts
1c5acb87 223
9542ac75 224 // Double_t pt = (gammai+gammaj).Pt();
1c5acb87 225 Double_t phi = (gammai+gammaj).Phi();
226 if(phi < 0)
227 phi+=TMath::TwoPi();
b960c7eb 228
1c5acb87 229 Double_t invmass = (gammai+gammaj).M();
230 Double_t angle = gammaj.Angle(gammai.Vect());
231 Double_t e = (gammai+gammaj).E();
9542ac75 232 Double_t asy = TMath::Abs((gammai-gammaj).E())/(gammai+gammaj).E();
233
1c5acb87 234 //Fill histograms with no cuts applied.
b960c7eb 235 if(fKeepNeutralMesonHistos)
236 {
fad7739b 237 fhAnglePairNoCut ->Fill(e,angle);
238 fhInvMassPairNoCut->Fill(e,invmass);
9542ac75 239 fhAsymmetryNoCut ->Fill(e,asy);
fad7739b 240 }
9542ac75 241
1c5acb87 242 //Cut on the aperture of the pair
b960c7eb 243 if(fUseAngleCut)
244 {
245 if(IsAngleInWindow(angle,e))
246 {
247 if(fKeepNeutralMesonHistos )
248 {
9542ac75 249 fhAnglePairOpeningAngleCut ->Fill(e,angle);
250 fhInvMassPairOpeningAngleCut->Fill(e,invmass);
251 fhAsymmetryOpeningAngleCut ->Fill(e,asy);
252 }
253 //AliDebug(2,Form("Angle cut: pt %f, phi %f",pt,phi));
254 } else return kFALSE;
255 }
256
257 // Asymmetry cut
b960c7eb 258 if(fUseAsymmetryCut)
259 {
260 if(fAsymmetryCut > asy)
261 {
262 if(fKeepNeutralMesonHistos)
263 {
9542ac75 264 fhInvMassPairAsymmetryCut->Fill(e,invmass);
265 fhAnglePairAsymmetryCut ->Fill(e,angle);
266 }
267 } else return kFALSE;
268 }
1c5acb87 269
9542ac75 270
271 //Cut on the invariant mass of the pair
3bfcb597 272
273 Float_t invmassmaxcut = fInvMassMaxCut;
80515912 274 Float_t invmassRightBandMinCut = fRightBandMinCut;
275 Float_t invmassRightBandMixCut = fRightBandMaxCut;
276
b960c7eb 277 if(calo=="EMCAL" && e > 6.)
278 { // for EMCAL, pi0s, mass depends strongly with energy for e > 6, loose max cut
3bfcb597 279
280 invmassmaxcut = (fInvMassMaxCutParam[0]+fInvMassMaxCut)+fInvMassMaxCutParam[1]*e+fInvMassMaxCutParam[2]*e*e;
80515912 281 invmassRightBandMinCut = (fInvMassMaxCutParam[0]+fRightBandMinCut)+fInvMassMaxCutParam[1]*e+fInvMassMaxCutParam[2]*e*e;
282 invmassRightBandMixCut = (fInvMassMaxCutParam[0]+fRightBandMaxCut)+fInvMassMaxCutParam[1]*e+fInvMassMaxCutParam[2]*e*e;
283
3bfcb597 284 //printf("e %f, max cut %f, p00 %f,p0 %f,p1 %f,p2 %f\n",
285 // e,invmassmaxcut,fInvMassMaxCut,fInvMassMaxCutParam[0],fInvMassMaxCutParam[1],fInvMassMaxCutParam[2]);
286 }
287
3f150b4b 288 // normal case, invariant mass selection around pi0/eta peak
289 if( !fParticle.Contains("SideBand") )
290 {
291 if( invmass > fInvMassMinCut && invmass < invmassmaxcut )
292 {
293 if(fKeepNeutralMesonHistos)
294 {
295 fhInvMassPairAllCut->Fill(e,invmass);
296 fhAnglePairAllCut ->Fill(e,angle);
297 fhAsymmetryAllCut ->Fill(e,asy);
298 }
299
300 //AliDebug(2,Form("IM cut: pt %f, phi %f",pt,phi));
301 return kTRUE;
302
303 }//(invmass>0.125) && (invmass<0.145)
304 else
b960c7eb 305 {
3f150b4b 306 return kFALSE;
307 }
308 }// normal selection
309
310 else // select a band around pi0/eta
311 {
312 if((invmass > fLeftBandMinCut && invmass < fLeftBandMaxCut ) ||
80515912 313 (invmass > invmassRightBandMinCut && invmass < invmassRightBandMixCut))
3f150b4b 314 {
315 if(fKeepNeutralMesonHistos)
316 {
317 fhInvMassPairAllCut->Fill(e,invmass);
318 fhAnglePairAllCut ->Fill(e,angle);
319 fhAsymmetryAllCut ->Fill(e,asy);
320 }
321
322 //AliDebug(2,Form("IM cut: pt %f, phi %f",pt,phi));
323 return kTRUE;
324
325 }//(invmass>0.125) && (invmass<0.145)
326 else
327 {
328 return kFALSE;
329 }
330 }
1c5acb87 331
332}
333
b960c7eb 334//_______________________________________________________________
335void AliNeutralMesonSelection::SetParticle(TString particleName)
336{
3bfcb597 337 // Set some default parameters for selection of pi0 or eta
3f150b4b 338
339 fParticle = particleName ;
3bfcb597 340
3f150b4b 341 if(particleName.Contains("Pi0"))
b960c7eb 342 {
343 fHistoNIMBins = 150 ;
344 fHistoIMMax = 0.3 ;
345 fHistoIMMin = 0. ;
3bfcb597 346
b960c7eb 347 fM = 0.135 ; // GeV
348 fInvMassMaxCut = 0.16 ; // GeV
349 fInvMassMinCut = 0.11 ; // GeV
3bfcb597 350
3f150b4b 351 fLeftBandMinCut = 0.05 ; // GeV
352 fLeftBandMaxCut = 0.09 ; // GeV
353 fRightBandMinCut = 0.160 ; // GeV
354 fRightBandMaxCut = 0.165 ; // GeV
355
3bfcb597 356 fInvMassMaxCutParam[0] = 0.0 ;
357 fInvMassMaxCutParam[1] =-7.e-5 ;
358 fInvMassMaxCutParam[2] = 8.e-5 ;
359
b960c7eb 360 fShiftMinAngle[0] =-0.03 ;
361 fShiftMinAngle[1] = 0.0025;
7287751b 362
363 fAngleMaxParam.AddAt( 0.8, 0) ;
364 fAngleMaxParam.AddAt(-1, 1) ;
365 fAngleMaxParam.AddAt( 0.09, 2) ; //for pi0 shift, for eta maybe 0.09
366 fAngleMaxParam.AddAt(-2.e-3,3) ;
3bfcb597 367
3f150b4b 368 }
369 else if(particleName.Contains("Eta"))
b960c7eb 370 {
371 fHistoNIMBins = 200 ; // GeV
372 fHistoIMMax = 0.75 ; // GeV
373 fHistoIMMin = 0.35 ; // GeV
374
375 fM = 0.547 ; // GeV
376 fInvMassMaxCut = 0.590 ; // GeV
377 fInvMassMinCut = 0.510 ; // GeV
3bfcb597 378
3f150b4b 379 fLeftBandMinCut = 0.450 ; // GeV
380 fLeftBandMaxCut = 0.500 ; // GeV
381 fRightBandMinCut = 0.600 ; // GeV
382 fRightBandMaxCut = 0.650 ; // GeV
7287751b 383
3f150b4b 384 fInvMassMaxCutParam[0] = 0.00 ;
385 fInvMassMaxCutParam[1] = 0.00 ;
386 fInvMassMaxCutParam[2] = 0.00 ;
387
388 fShiftMinAngle[0] =-0.03 ;
389 fShiftMinAngle[0] = 0.00 ;
3bfcb597 390
7287751b 391 fAngleMaxParam.AddAt( 0.80, 0) ; // Same as pi0
392 fAngleMaxParam.AddAt(-0.25, 1) ; // Same as pi0
393 fAngleMaxParam.AddAt( 0.12, 2) ; // Shifted with respect to pi0
394 fAngleMaxParam.AddAt(-5.e-4, 3) ; // Same as pi0
3bfcb597 395
396 }
397 else
398 printf("AliAnaNeutralMesonSelection::SetParticle(%s) *** Particle NOT defined (Pi0 or Eta), Pi0 settings by default *** \n",particleName.Data());
399
400
401}
402
b960c7eb 403//______________________________________________________________
1c5acb87 404void AliNeutralMesonSelection::Print(const Option_t * opt) const
405{
406
407 //Print some relevant parameters set for the analysis
408 if(! opt)
409 return;
410
a3aebfff 411 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
1c5acb87 412
3f150b4b 413 printf("Particle %s, mass : %f \n", fParticle.Data(), fM );
1c5acb87 414 printf("Invariant mass limits : %f < m < %f \n", fInvMassMinCut , fInvMassMinCut );
9542ac75 415
416 printf("Use asymmetry cut? : %d ; A < %f \n", fUseAngleCut, fAsymmetryCut );
417
418 printf("Use angle cut? : %d \n", fUseAngleCut );
3f150b4b 419 if(fUseAngleCut)
420 {
9542ac75 421 printf("Angle selection param: \n");
422 printf("p0 : %f\n", fAngleMaxParam.At(0));
423 printf("p1 : %f\n", fAngleMaxParam.At(1));
424 printf("p2 : %f\n", fAngleMaxParam.At(2));
425 printf("p3 : %f\n", fAngleMaxParam.At(3));
7287751b 426 printf("Min angle shift : p0 = %1.3f, p1 = %1.3f,\n", fShiftMinAngle[0],fShiftMinAngle[1]);
9542ac75 427 }
428
1c5acb87 429 printf("Keep Neutral Meson Histos = %d\n",fKeepNeutralMesonHistos);
430
3f150b4b 431 if(fKeepNeutralMesonHistos)
432 {
a5fb4114 433 printf("Histograms: %3.1f < E < %3.1f, Nbin = %d\n", fHistoEMin, fHistoEMax, fHistoNEBins);
477d6cee 434 printf("Histograms: %3.1f < angle < %3.1f, Nbin = %d\n", fHistoAngleMin, fHistoAngleMax, fHistoNAngleBins);
3bfcb597 435 printf("Histograms: %3.1f < IM < %3.1f, Nbin = %d\n", fHistoIMMin, fHistoIMMax, fHistoNIMBins);
1c5acb87 436 }
477d6cee 437
1c5acb87 438}