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