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