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