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