Change binning of histograms, correct problem with delete of PHOS geometry pointer...
[u/mrichter/AliRoot.git] / PWG4 / PartCorrBase / AliNeutralMesonSelection.cxx
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>
27  
28 //---- AliRoot system ----
29 #include "AliNeutralMesonSelection.h" 
30
31 ClassImp(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),
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 {
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 //____________________________________________________________________________
61 AliNeutralMesonSelection::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 //_________________________________________________________________________
81 AliNeutralMesonSelection & 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 //____________________________________________________________________________
111 AliNeutralMesonSelection::~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 //________________________________________________________________________
130 TList *  AliNeutralMesonSelection::GetCreateOutputObjects()
131 {  
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         if(fKeepNeutralMesonHistos) outputContainer->SetOwner(kFALSE);
138         
139   fhAnglePairNoCut  = new TH2F
140     ("AnglePairNoCut",
141      "Angle between all #gamma pair vs E_{#pi^{0}}",fHistoNEBins,fHistoEMin,fHistoEMax,fHistoNAngleBins,fHistoAngleMin,fHistoAngleMax); 
142   fhAnglePairNoCut->SetYTitle("Angle (rad)");
143   fhAnglePairNoCut->SetXTitle("E_{ #pi^{0}} (GeV)");
144   
145   fhAnglePairOpeningAngleCut  = new TH2F
146     ("AnglePairOpeningAngleCut",
147      "Angle between all #gamma pair (opening angle + azimuth cut) vs E_{#pi^{0}}"
148      ,fHistoNEBins,fHistoEMin,fHistoEMax,fHistoNAngleBins,fHistoAngleMin,fHistoAngleMax); 
149   fhAnglePairOpeningAngleCut->SetYTitle("Angle (rad)");
150   fhAnglePairOpeningAngleCut->SetXTitle("E_{ #pi^{0}} (GeV)");
151   
152   fhAnglePairAllCut  = new TH2F
153     ("AnglePairAllCut",
154      "Angle between all #gamma pair (opening angle + inv mass cut+azimuth) vs E_{#pi^{0}}"
155      ,fHistoNEBins,fHistoEMin,fHistoEMax,fHistoNAngleBins,fHistoAngleMin,fHistoAngleMax); 
156   fhAnglePairAllCut->SetYTitle("Angle (rad)");
157   fhAnglePairAllCut->SetXTitle("E_{ #pi^{0}} (GeV)");    
158   
159   //
160   fhInvMassPairNoCut  = new TH2F
161     ("InvMassPairNoCut","Invariant Mass of all #gamma pair vs E_{#pi^{0}}",
162      fHistoNPtBins,fHistoPtMin,fHistoPtMax,fHistoNIMBins,fHistoIMMin,fHistoIMMax); 
163   fhInvMassPairNoCut->SetYTitle("Invariant Mass (GeV/c^{2})");
164   fhInvMassPairNoCut->SetXTitle("E_{ #pi^{0}} (GeV)");
165   
166   fhInvMassPairOpeningAngleCut  = new TH2F
167     ("InvMassPairOpeningAngleCut",
168      "Invariant Mass of #gamma pair (angle cut) vs E_{#pi^{0}}",
169      fHistoNPtBins,fHistoPtMin,fHistoPtMax,fHistoNIMBins,fHistoIMMin,fHistoIMMax); 
170   fhInvMassPairOpeningAngleCut->SetYTitle("Invariant Mass (GeV/c^{2})");
171   fhInvMassPairOpeningAngleCut->SetXTitle(" E_{#pi^{0}}(GeV)");
172   
173   fhInvMassPairAllCut  = new TH2F
174     ("InvMassPairAllCut",
175      "Invariant Mass of #gamma pair (opening angle+invmass cut) vs E_{#pi^{0}}",
176      fHistoNPtBins,fHistoPtMin,fHistoPtMax,fHistoNIMBins,fHistoIMMin,fHistoIMMax); 
177   fhInvMassPairAllCut->SetYTitle("Invariant Mass (GeV/c^{2})");
178   fhInvMassPairAllCut->SetXTitle("E_{#pi^{0}}(GeV)");
179   
180   outputContainer->Add(fhAnglePairNoCut) ; 
181   outputContainer->Add(fhAnglePairOpeningAngleCut) ;
182   outputContainer->Add(fhAnglePairAllCut) ; 
183   
184   outputContainer->Add(fhInvMassPairNoCut) ; 
185   outputContainer->Add(fhInvMassPairOpeningAngleCut) ; 
186   outputContainer->Add(fhInvMassPairAllCut) ; 
187   
188   return outputContainer;
189 }
190
191 //____________________________________________________________________________
192 void AliNeutralMesonSelection::InitParameters()
193 {
194   
195   //Initialize the parameters of the analysis.
196   fKeepNeutralMesonHistos = kFALSE ;
197   
198   fAngleMaxParam.Set(4) ;
199   fAngleMaxParam.AddAt(0.4,0);//={0.4,-0.25,0.025,-2e-4};
200   fAngleMaxParam.AddAt(-0.25,1) ;
201   fAngleMaxParam.AddAt(0.025,2) ;
202   fAngleMaxParam.AddAt(-2e-4,3) ;
203
204   fInvMassMaxCut  = 0.16 ;
205   fInvMassMinCut  = 0.11 ;
206
207   fM = 0.1349766;//neutralMeson mass, pi0
208   
209  //Histogrammes settings
210   fHistoNEBins = 100 ;
211   fHistoEMax   = 20 ;
212   fHistoEMin   = 0.  ;  
213   
214   fHistoNPtBins = 100 ;
215   fHistoPtMax   = 20 ;
216   fHistoPtMin   = 0.  ;
217
218   fHistoNAngleBins = 100 ;
219   fHistoAngleMax   = 0.2;
220   fHistoAngleMin   = 0.  ;
221
222   fHistoNIMBins = 100 ;
223   fHistoIMMax   = 0.3   ;
224   fHistoIMMin   = 0.  ;  
225 }
226
227 //__________________________________________________________________________-
228 Bool_t AliNeutralMesonSelection::IsAngleInWindow(const Float_t angle,const Float_t e) const {
229   //Check if the opening angle of the candidate pairs is inside 
230   //our selection windowd
231
232   Bool_t result = kFALSE;
233   Double_t max =  fAngleMaxParam.At(0)*TMath::Exp(fAngleMaxParam.At(1)*e)
234     +fAngleMaxParam.At(2)+fAngleMaxParam.At(3)*e;
235   Double_t arg = (e*e-2*fM*fM)/(e*e);
236   Double_t min = 100. ;
237   if(arg>0.)
238     min = TMath::ACos(arg);
239
240   if((angle<max)&&(angle>=min))
241     result = kTRUE;
242  
243   return result;
244 }
245
246 //____________________________________________________________________________
247 Bool_t  AliNeutralMesonSelection::SelectPair(TLorentzVector gammai, TLorentzVector gammaj)  
248 {  
249   
250   //Search for the neutral pion within selection cuts
251   Bool_t goodpair = kFALSE ;
252   
253 //  Double_t pt  = (gammai+gammaj).Pt();
254   Double_t phi = (gammai+gammaj).Phi();
255   if(phi < 0)
256     phi+=TMath::TwoPi();
257   Double_t invmass = (gammai+gammaj).M();
258   Double_t angle   = gammaj.Angle(gammai.Vect());
259   Double_t e       = (gammai+gammaj).E();
260   
261   //Fill histograms with no cuts applied.
262   fhAnglePairNoCut->Fill(e,angle);
263   fhInvMassPairNoCut->Fill(e,invmass);
264     
265   //Cut on the aperture of the pair
266   if(IsAngleInWindow(angle,e)){
267     fhAnglePairOpeningAngleCut     ->Fill(e,angle);
268     fhInvMassPairOpeningAngleCut->Fill(e,invmass);
269     //AliDebug(2,Form("Angle cut: pt %f, phi %f",pt,phi));
270     
271     //Cut on the invariant mass of the pair
272     if((invmass>fInvMassMinCut) && (invmass<fInvMassMaxCut)){ 
273       fhInvMassPairAllCut  ->Fill(e,invmass);
274       fhAnglePairAllCut       ->Fill(e,angle);
275       goodpair = kTRUE;
276       //AliDebug(2,Form("IM cut: pt %f, phi %f",pt,phi));
277     }//(invmass>0.125) && (invmass<0.145)
278   }//Opening angle cut
279   
280   return goodpair; 
281   
282 }
283
284 //__________________________________________________________________
285 void AliNeutralMesonSelection::Print(const Option_t * opt) const
286 {
287
288   //Print some relevant parameters set for the analysis
289   if(! opt)
290     return;
291   
292   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
293
294   printf("mass : %f  \n", fM );
295   printf("Invariant mass limits : %f < m < %f \n", fInvMassMinCut , fInvMassMinCut );
296   printf("Angle selection param: \n");
297   printf("p0 :     %f\n", fAngleMaxParam.At(0));
298   printf("p1 :     %f\n", fAngleMaxParam.At(1));
299   printf("p2 :     %f\n", fAngleMaxParam.At(2));
300   printf("p3 :     %f\n", fAngleMaxParam.At(3));
301
302   printf("Keep Neutral Meson Histos = %d\n",fKeepNeutralMesonHistos);
303   
304   if(fKeepNeutralMesonHistos){
305     printf("Histograms: %3.1f < E  < %3.1f,  Nbin = %d\n", fHistoEMin,  fHistoEMax,  fHistoNEBins);
306     printf("Histograms: %3.1f < pT < %3.1f,  Nbin = %d\n", fHistoPtMin,  fHistoPtMax,  fHistoNPtBins);
307     printf("Histograms: %3.1f < angle < %3.1f, Nbin = %d\n", fHistoAngleMin, fHistoAngleMax, fHistoNAngleBins);
308     printf("Histograms: %3.1f < IM < %3.1f, Nbin = %d\n", fHistoIMMin, fHistoIMMax, fHistoNIMBins);
309     
310   }
311   
312