f0cc48a9b20e1eef6a8e46ba6ce060c609faa549
[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 // There is a 3rd cut based on the gamma correlation on phi or pt.
37 //-- Author: Gustavo Conesa (INFN-LNF)
38
39 // --- ROOT system ---
40 #include <TParticle.h>
41 #include <TLorentzVector.h>
42 #include <TH2.h>
43 #include <TList.h>
44 #include <TArrayD.h>
45  
46 //---- AliRoot system ----
47 #include "AliNeutralMesonSelection.h" 
48 #include "AliLog.h"
49
50 ClassImp(AliNeutralMesonSelection)
51   
52   
53 //____________________________________________________________________________
54   AliNeutralMesonSelection::AliNeutralMesonSelection() : 
55     TObject(), fSelect(0), fM(0),
56     fInvMassMaxCut(0.), fInvMassMinCut(0.),
57     fAngleMaxParam(),  fMinPt(0),
58     fDeltaPhiMaxCut(0.), fDeltaPhiMinCut(0.),
59     fRatioMaxCut(0), fRatioMinCut(0),  fKeepNeutralMesonHistos(0),
60     fhAnglePairNoCut(0),  fhAnglePairCorrelationCut(0),
61     fhAnglePairOpeningAngleCut(0), fhAnglePairAllCut(0), 
62     fhInvMassPairNoCut(0),   fhInvMassPairCorrelationCut(0), 
63     fhInvMassPairOpeningAngleCut(0), fhInvMassPairAllCut(0) 
64 {
65   //Default Ctor
66   
67   //Initialize parameters
68   
69   // kGammaHadron and kGammaJet 
70   fAngleMaxParam.Set(4) ;
71   fAngleMaxParam.Reset(0.);
72   
73   //Initialize parameters
74   InitParameters();
75 }
76
77 //____________________________________________________________________________
78 AliNeutralMesonSelection::AliNeutralMesonSelection(const AliNeutralMesonSelection & g) :   
79   TObject(),  
80   fSelect(g.fSelect), fM(g.fM),
81   fInvMassMaxCut(g.fInvMassMaxCut), fInvMassMinCut(g.fInvMassMinCut),
82   fAngleMaxParam(g.fAngleMaxParam), fMinPt(g.fMinPt),
83   fDeltaPhiMaxCut(g.fDeltaPhiMaxCut), fDeltaPhiMinCut(g.fDeltaPhiMinCut), 
84   fRatioMaxCut(g.fRatioMaxCut), fRatioMinCut(g.fRatioMinCut), 
85   fKeepNeutralMesonHistos(g.fKeepNeutralMesonHistos),
86   fhAnglePairNoCut(g. fhAnglePairNoCut), 
87   fhAnglePairCorrelationCut(g. fhAnglePairCorrelationCut), 
88   fhAnglePairOpeningAngleCut(g. fhAnglePairOpeningAngleCut), 
89   fhAnglePairAllCut(g. fhAnglePairAllCut), 
90   fhInvMassPairNoCut(g.fhInvMassPairNoCut),  
91   fhInvMassPairCorrelationCut(g.fhInvMassPairCorrelationCut), 
92   fhInvMassPairOpeningAngleCut(g.fhInvMassPairOpeningAngleCut), 
93   fhInvMassPairAllCut(g.fhInvMassPairAllCut)
94 {
95   // cpy ctor
96 }
97
98 //_________________________________________________________________________
99 AliNeutralMesonSelection & AliNeutralMesonSelection::operator = (const AliNeutralMesonSelection & source)
100 {
101   // assignment operator
102   
103   if(this == &source)return *this;
104   ((TObject *)this)->operator=(source);
105
106   fSelect = source.fSelect ;
107   fM = source.fM ;
108   fInvMassMaxCut = source.fInvMassMaxCut ; 
109   fInvMassMinCut = source.fInvMassMinCut ;
110   fAngleMaxParam = source.fAngleMaxParam ;
111   fMinPt = source.fMinPt ;
112   fDeltaPhiMaxCut = source.fDeltaPhiMaxCut ; 
113   fDeltaPhiMinCut = source.fDeltaPhiMinCut ; 
114   fRatioMaxCut = source.fRatioMaxCut ; 
115   fRatioMinCut = source.fRatioMinCut ;
116   fKeepNeutralMesonHistos = source.fKeepNeutralMesonHistos;
117  
118   fhAnglePairNoCut = source. fhAnglePairNoCut ; 
119   fhAnglePairCorrelationCut = source. fhAnglePairCorrelationCut ; 
120   fhAnglePairOpeningAngleCut = source. fhAnglePairOpeningAngleCut ; 
121   fhAnglePairAllCut = source. fhAnglePairAllCut ; 
122   fhInvMassPairNoCut = source.fhInvMassPairNoCut ; 
123   fhInvMassPairCorrelationCut = source.fhInvMassPairCorrelationCut ; 
124   fhInvMassPairOpeningAngleCut = source.fhInvMassPairOpeningAngleCut ; 
125   fhInvMassPairAllCut = source.fhInvMassPairAllCut ; 
126   
127   return *this;
128   
129 }
130
131 //____________________________________________________________________________
132 AliNeutralMesonSelection::~AliNeutralMesonSelection() 
133 {
134  // Remove all pointers except analysis output pointers.
135
136 }
137
138
139
140 //________________________________________________________________________
141 TList *  AliNeutralMesonSelection::GetCreateOutputObjects()
142 {  
143
144   // Create histograms to be saved in output file and 
145   // store them in outputContainer
146   TList * outputContainer = new TList() ; 
147   outputContainer->SetName("MesonDecayHistos") ; 
148   
149   fhAnglePairNoCut  = new TH2F
150     ("AnglePairNoCut",
151      "Angle between all #gamma pair vs E_{#pi^{0}}",200,0,50,200,0,0.2); 
152   fhAnglePairNoCut->SetYTitle("Angle (rad)");
153   fhAnglePairNoCut->SetXTitle("E_{ #pi^{0}} (GeV)");
154   
155   fhAnglePairOpeningAngleCut  = new TH2F
156     ("AnglePairOpeningAngleCut",
157      "Angle between all #gamma pair (opening angle + azimuth cut) vs E_{#pi^{0}}"
158      ,200,0,50,200,0,0.2); 
159   fhAnglePairOpeningAngleCut->SetYTitle("Angle (rad)");
160   fhAnglePairOpeningAngleCut->SetXTitle("E_{ #pi^{0}} (GeV)");
161   
162   fhAnglePairAllCut  = new TH2F
163     ("AnglePairAllCut",
164      "Angle between all #gamma pair (opening angle + inv mass cut+azimuth) vs E_{#pi^{0}}"
165      ,200,0,50,200,0,0.2); 
166   fhAnglePairAllCut->SetYTitle("Angle (rad)");
167   fhAnglePairAllCut->SetXTitle("E_{ #pi^{0}} (GeV)");    
168   
169   //
170   fhInvMassPairNoCut  = new TH2F
171     ("InvMassPairNoCut","Invariant Mass of all #gamma pair vs E_{#pi^{0}}",
172      120,0,120,360,0,0.5); 
173   fhInvMassPairNoCut->SetYTitle("Invariant Mass (GeV/c^{2})");
174   fhInvMassPairNoCut->SetXTitle("E_{ #pi^{0}} (GeV)");
175   
176   fhInvMassPairOpeningAngleCut  = new TH2F
177     ("InvMassPairOpeningAngleCut",
178      "Invariant Mass of #gamma pair (angle cut) vs E_{#pi^{0}}",
179      120,0,120,360,0,0.5); 
180   fhInvMassPairOpeningAngleCut->SetYTitle("Invariant Mass (GeV/c^{2})");
181   fhInvMassPairOpeningAngleCut->SetXTitle(" E_{#pi^{0}}(GeV)");
182   
183   fhInvMassPairAllCut  = new TH2F
184     ("InvMassPairAllCut",
185      "Invariant Mass of #gamma pair (opening angle+invmass cut) vs E_{#pi^{0}}",
186      120,0,120,360,0,0.5); 
187   fhInvMassPairAllCut->SetYTitle("Invariant Mass (GeV/c^{2})");
188   fhInvMassPairAllCut->SetXTitle("E_{#pi^{0}}(GeV)");
189
190   fhAnglePairCorrelationCut  = new TH2F
191     ("AnglePairCorrelationCut",
192      "Angle between correlated #gamma pair vs E_{#pi^{0}}",200,0,50,200,0,0.2); 
193   fhAnglePairCorrelationCut->SetYTitle("Angle (rad)");
194   fhAnglePairCorrelationCut->SetXTitle("E_{ #pi^{0}} (GeV)");
195
196   fhInvMassPairCorrelationCut  = new TH2F
197     ("InvMassPairCorrelationCut","Invariant Mass of correlated #gamma pair vs E_{#pi^{0}}",
198      120,0,120,360,0,0.5); 
199   fhInvMassPairCorrelationCut->SetYTitle("Invariant Mass (GeV/c^{2})");
200   fhInvMassPairCorrelationCut->SetXTitle("E_{ #pi^{0}} (GeV)");
201   
202   outputContainer->Add(fhAnglePairNoCut) ; 
203   outputContainer->Add(fhAnglePairOpeningAngleCut) ;
204   outputContainer->Add(fhAnglePairAllCut) ; 
205   
206   outputContainer->Add(fhInvMassPairNoCut) ; 
207   outputContainer->Add(fhInvMassPairOpeningAngleCut) ; 
208   outputContainer->Add(fhInvMassPairAllCut) ; 
209
210   outputContainer->Add(fhAnglePairCorrelationCut) ; 
211   outputContainer->Add(fhInvMassPairCorrelationCut) ; 
212   
213   return outputContainer;
214 }
215
216  //____________________________________________________________________________
217 void AliNeutralMesonSelection::InitParameters()
218 {
219  
220   //Initialize the parameters of the analysis.
221   fKeepNeutralMesonHistos = kTRUE ;
222
223   //-------------kHadron, kJetLeadCone-----------------
224   fAngleMaxParam.Set(4) ;
225   fAngleMaxParam.AddAt(0.4,0);//={0.4,-0.25,0.025,-2e-4};
226   fAngleMaxParam.AddAt(-0.25,1) ;
227   fAngleMaxParam.AddAt(0.025,2) ;
228   fAngleMaxParam.AddAt(-2e-4,3) ;
229
230   fInvMassMaxCut  = 0.16 ;
231   fInvMassMinCut  = 0.11 ;
232
233   fM = 0.1349766;//neutralMeson mass
234
235   fMinPt = 0.   ;
236   fDeltaPhiMaxCut      = 4.5;
237   fDeltaPhiMinCut      = 1.5 ;
238   fRatioMaxCut    = 1.0 ;
239   fRatioMinCut    = 0.1 ; 
240 }
241
242 //__________________________________________________________________________-
243 Bool_t AliNeutralMesonSelection::IsAngleInWindow(const Float_t angle,const Float_t e) const {
244   //Check if the opening angle of the candidate pairs is inside 
245   //our selection windowd
246
247   Bool_t result = kFALSE;
248   Double_t max =  fAngleMaxParam.At(0)*TMath::Exp(fAngleMaxParam.At(1)*e)
249     +fAngleMaxParam.At(2)+fAngleMaxParam.At(3)*e;
250   Double_t arg = (e*e-2*fM*fM)/(e*e);
251   Double_t min = 100. ;
252   if(arg>0.)
253     min = TMath::ACos(arg);
254
255   if((angle<max)&&(angle>=min))
256     result = kTRUE;
257  
258   return result;
259 }
260
261 //____________________________________________________________________________
262 Bool_t  AliNeutralMesonSelection::CutPtPhi(Double_t ptg, Double_t phig, Double_t pt, Double_t phi)  const
263
264   //Select pair if delta
265   Bool_t cut = kFALSE ;
266  
267   if(fSelect == kNoSelectPhiPt) cut = kTRUE ;
268   else if((phig-phi) > fDeltaPhiMinCut && ((phig-phi) < fDeltaPhiMaxCut)){
269     //Cut on pt
270     if((fSelect == kSelectPhiPtRatio && ptg > 0. && pt/ptg  > fRatioMinCut &&  pt/ptg  < fRatioMaxCut) ||
271         (fSelect == kSelectPhiMinPt && pt > fMinPt)  )  cut = kTRUE ;
272   }
273   else cut = kFALSE ;
274   
275   return cut ;
276   
277 }
278
279 //____________________________________________________________________________
280 Bool_t  AliNeutralMesonSelection::SelectPair(TParticle * pGamma, TLorentzVector gammai, TLorentzVector gammaj)  
281 {  
282   
283   //Search for the neutral pion within selection cuts
284   
285   Double_t ptg = pGamma->Pt();
286   Double_t phig = pGamma->Phi() ;
287   Bool_t goodpair = kFALSE ;
288   
289   Double_t pt  = (gammai+gammaj).Pt();
290   Double_t phi = (gammai+gammaj).Phi();
291   if(phi < 0)
292     phi+=TMath::TwoPi();
293   Double_t invmass = (gammai+gammaj).M();
294   Double_t angle   = gammaj.Angle(gammai.Vect());
295   Double_t e       = (gammai+gammaj).E();
296   
297   //Fill histograms with no cuts applied.
298   fhAnglePairNoCut->Fill(e,angle);
299   fhInvMassPairNoCut->Fill(e,invmass);
300   
301   //Cut on phig-phi meson
302   if(CutPtPhi(ptg, phig, pt, phi)){
303     
304     fhAnglePairCorrelationCut     ->Fill(e,angle);
305     fhInvMassPairCorrelationCut->Fill(e,invmass);
306     
307     //Cut on the aperture of the pair
308     if(IsAngleInWindow(angle,e)){
309       fhAnglePairOpeningAngleCut     ->Fill(e,angle);
310       fhInvMassPairOpeningAngleCut->Fill(e,invmass);
311       AliDebug(2,Form("Angle cut: pt %f, phi %f",pt,phi));
312       
313       //Cut on the invariant mass of the pair
314       if((invmass>fInvMassMinCut) && (invmass<fInvMassMaxCut)){ 
315         fhInvMassPairAllCut  ->Fill(e,invmass);
316         fhAnglePairAllCut       ->Fill(e,angle);
317         goodpair = kTRUE;
318         AliDebug(2,Form("IM cut: pt %f, phi %f",pt,phi));
319       }//(invmass>0.125) && (invmass<0.145)
320     }//Opening angle cut
321   } // cut on pt and phi
322   
323   
324   return goodpair; 
325   
326 }
327
328 //__________________________________________________________________
329 void AliNeutralMesonSelection::Print(const Option_t * opt) const
330 {
331
332   //Print some relevant parameters set for the analysis
333   if(! opt)
334     return;
335   
336   Info("Print", "%s %s", GetName(), GetTitle() ) ;
337
338   printf("mass : %f  \n", fM );
339   printf("Invariant mass limits : %f < m < %f \n", fInvMassMinCut , fInvMassMinCut );
340   printf("Angle selection param: \n");
341   printf("p0 :     %f", fAngleMaxParam.At(0));
342   printf("p1 :     %f", fAngleMaxParam.At(1));
343   printf("p2 :     %f", fAngleMaxParam.At(2));
344   printf("p3 :     %f", fAngleMaxParam.At(3));
345
346   printf("pT meson       >    %f\n", fMinPt) ; 
347   printf("Phi gamma-meson      <     %f\n", fDeltaPhiMaxCut) ; 
348   printf("Phi gamma-meson      >     %f\n", fDeltaPhiMinCut) ;
349   printf("pT meson / pT Gamma             <     %f\n", fRatioMaxCut) ; 
350   printf("pT meson / pT Gamma             >     %f\n", fRatioMinCut) ;
351   printf("Keep Neutral Meson Histos = %d\n",fKeepNeutralMesonHistos);
352
353