move the decay bit definition from AliAnaPi0EbE to AliNeutralSelection, define there...
[u/mrichter/AliRoot.git] / PWG / CaloTrackCorrBase / 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
16 //_________________________________________________________________________
17 // Class that contains methods to select candidate pairs to neutral meson 
18 // 2 main selections, invariant mass around pi0 (also any other mass),
19 // apperture angle to distinguish from combinatorial.
20 //-- Author: Gustavo Conesa (INFN-LNF)
21
22 // --- ROOT system ---
23 #include <TLorentzVector.h>
24 #include <TH2.h>
25 #include <TList.h>
26 #include <TString.h>
27
28 //---- AliRoot system ----
29 #include "AliNeutralMesonSelection.h" 
30
31 ClassImp(AliNeutralMesonSelection)
32   
33   
34 //______________________________________________________
35   AliNeutralMesonSelection::AliNeutralMesonSelection() : 
36     TObject(),             fAsymmetryCut(1),                
37     fUseAsymmetryCut(0),   fM(0),                 
38     fInvMassMaxCut(0.),    fInvMassMinCut(0.),   
39     fLeftBandMinCut(0.),   fLeftBandMaxCut(0.),            
40     fRightBandMinCut(0.),  fRightBandMaxCut(0.),              
41     fAngleMaxParam(),      fUseAngleCut(0),       
42     fKeepNeutralMesonHistos(0),
43     fParticle(""),         fDecayBit(0),
44     // histograms
45     fhAnglePairNoCut(0),          fhAnglePairOpeningAngleCut(0),   
46     fhAnglePairAsymmetryCut(0),   fhAnglePairAllCut(0), 
47     fhInvMassPairNoCut(0),        fhInvMassPairOpeningAngleCut(0), 
48     fhInvMassPairAsymmetryCut(0), fhInvMassPairAllCut(0),
49     fhAsymmetryNoCut(0),          fhAsymmetryOpeningAngleCut(0),   
50     fhAsymmetryAllCut(0),
51     // histogram ranges and bins
52     fHistoNEBins(0),       fHistoEMax(0.),                  fHistoEMin(0.),
53     fHistoNAngleBins(0),   fHistoAngleMax(0.),              fHistoAngleMin(0.),
54     fHistoNIMBins(0),      fHistoIMMax(0.),                 fHistoIMMin(0.)
55 {
56   //Default Ctor
57   
58   //Initialize parameters
59   InitParameters();
60 }
61
62 //_________________________________________________________
63 TList *  AliNeutralMesonSelection::GetCreateOutputObjects()
64 {  
65   // Create histograms to be saved in output file and 
66   // store them in outputContainer of the analysis class that calls this class.
67   
68   TList * outputContainer = new TList() ; 
69   outputContainer->SetName("MesonDecayHistos") ; 
70   
71   if(fKeepNeutralMesonHistos){
72           
73           outputContainer->SetOwner(kFALSE);
74           
75           fhAnglePairNoCut  = new TH2F
76           ("AnglePairNoCut",
77            "Angle between all #gamma pair vs E_{#pi^{0}}",fHistoNEBins,fHistoEMin,fHistoEMax,fHistoNAngleBins,fHistoAngleMin,fHistoAngleMax); 
78           fhAnglePairNoCut->SetYTitle("Angle (rad)");
79           fhAnglePairNoCut->SetXTitle("E_{ #pi^{0}} (GeV)");
80               
81     fhAsymmetryNoCut  = new TH2F
82           ("AsymmetryNoCut","Asymmetry of all #gamma pair vs E_{#pi^{0}}",
83            fHistoNEBins,fHistoEMin,fHistoEMax,100,0,1); 
84           fhAsymmetryNoCut->SetYTitle("Asymmetry");
85           fhAsymmetryNoCut->SetXTitle("E_{ #pi^{0}} (GeV)");    
86     
87     fhInvMassPairNoCut  = new TH2F
88           ("InvMassPairNoCut","Invariant Mass of all #gamma pair vs E_{#pi^{0}}",
89            fHistoNEBins,fHistoEMin,fHistoEMax,fHistoNIMBins,fHistoIMMin,fHistoIMMax); 
90           fhInvMassPairNoCut->SetYTitle("Invariant Mass (GeV/c^{2})");
91           fhInvMassPairNoCut->SetXTitle("E_{ #pi^{0}} (GeV)");    
92     
93     outputContainer->Add(fhAnglePairNoCut) ; 
94           outputContainer->Add(fhAsymmetryNoCut) ; 
95     outputContainer->Add(fhInvMassPairNoCut) ; 
96
97     if(fUseAngleCut) {
98       fhAnglePairOpeningAngleCut  = new TH2F
99       ("AnglePairOpeningAngleCut",
100        "Angle between all #gamma pair (opening angle) vs E_{#pi^{0}}"
101        ,fHistoNEBins,fHistoEMin,fHistoEMax,fHistoNAngleBins,fHistoAngleMin,fHistoAngleMax); 
102       fhAnglePairOpeningAngleCut->SetYTitle("Angle (rad)");
103       fhAnglePairOpeningAngleCut->SetXTitle("E_{ #pi^{0}} (GeV)");
104       
105       fhAsymmetryOpeningAngleCut  = new TH2F
106       ("AsymmetryOpeningAngleCut",
107        "Asymmetry of #gamma pair (angle cut) vs E_{#pi^{0}}",
108        fHistoNEBins,fHistoEMin,fHistoEMax,100,0,1); 
109       fhAsymmetryOpeningAngleCut->SetYTitle("Asymmetry");
110       fhAsymmetryOpeningAngleCut->SetXTitle(" E_{#pi^{0}}(GeV)");   
111       
112       fhInvMassPairOpeningAngleCut  = new TH2F
113       ("InvMassPairOpeningAngleCut",
114        "Invariant Mass of #gamma pair (angle cut) vs E_{#pi^{0}}",
115        fHistoNEBins,fHistoEMin,fHistoEMax,fHistoNIMBins,fHistoIMMin,fHistoIMMax); 
116       fhInvMassPairOpeningAngleCut->SetYTitle("Invariant Mass (GeV/c^{2})");
117       fhInvMassPairOpeningAngleCut->SetXTitle(" E_{#pi^{0}}(GeV)");
118       
119       outputContainer->Add(fhAnglePairOpeningAngleCut) ;
120       outputContainer->Add(fhAsymmetryOpeningAngleCut) ;
121       outputContainer->Add(fhInvMassPairOpeningAngleCut) ;
122     }
123     
124           if(fUseAsymmetryCut) {
125       fhAnglePairAsymmetryCut  = new TH2F
126       ("AnglePairAsymmetryCut",
127        "Angle between all #gamma pair (opening angle + asymetry cut) vs E_{#pi^{0}}"
128        ,fHistoNEBins,fHistoEMin,fHistoEMax,fHistoNAngleBins,fHistoAngleMin,fHistoAngleMax); 
129       fhAnglePairAsymmetryCut->SetYTitle("Angle (rad)");
130       fhAnglePairAsymmetryCut->SetXTitle("E_{ #pi^{0}} (GeV)");
131       
132       fhInvMassPairAsymmetryCut  = new TH2F
133       ("InvMassPairAsymmetryCut",
134        "Invariant Mass of #gamma pair (opening angle + asymmetry) vs E_{#pi^{0}}",
135        fHistoNEBins,fHistoEMin,fHistoEMax,fHistoNIMBins,fHistoIMMin,fHistoIMMax); 
136       fhInvMassPairAsymmetryCut->SetYTitle("Invariant Mass (GeV/c^{2})");
137       fhInvMassPairAsymmetryCut->SetXTitle("E_{#pi^{0}}(GeV)");      
138       
139       outputContainer->Add(fhAnglePairAsymmetryCut) ;
140       outputContainer->Add(fhInvMassPairAsymmetryCut) ;
141       
142     }
143     
144           fhAnglePairAllCut  = new TH2F
145           ("AnglePairAllCut",
146            "Angle between all #gamma pair (opening angle + asymmetry + inv mass cut) vs E_{#pi^{0}}"
147            ,fHistoNEBins,fHistoEMin,fHistoEMax,fHistoNAngleBins,fHistoAngleMin,fHistoAngleMax); 
148           fhAnglePairAllCut->SetYTitle("Angle (rad)");
149           fhAnglePairAllCut->SetXTitle("E_{ #pi^{0}} (GeV)");    
150         
151           fhInvMassPairAllCut  = new TH2F
152           ("InvMassPairAllCut",
153            "Invariant Mass of #gamma pair (opening angle + asymmetry + invmass cut) vs E_{#pi^{0}}",
154            fHistoNEBins,fHistoEMin,fHistoEMax,fHistoNIMBins,fHistoIMMin,fHistoIMMax); 
155           fhInvMassPairAllCut->SetYTitle("Invariant Mass (GeV/c^{2})");
156           fhInvMassPairAllCut->SetXTitle("E_{#pi^{0}}(GeV)");
157           
158           fhAsymmetryAllCut  = new TH2F
159           ("AsymmetryAllCut",
160            "Asymmetry of #gamma pair (opening angle+invmass cut) vs E_{#pi^{0}}",
161            fHistoNEBins,fHistoEMin,fHistoEMax,100,0,1); 
162           fhAsymmetryAllCut->SetYTitle("Asymmetry");
163           fhAsymmetryAllCut->SetXTitle("E_{#pi^{0}}(GeV)");
164     
165     outputContainer->Add(fhAnglePairAllCut) ; 
166           outputContainer->Add(fhAsymmetryAllCut) ; 
167     outputContainer->Add(fhInvMassPairAllCut) ;     
168
169   }
170   
171   return outputContainer;
172
173 }
174
175 //_____________________________________________
176 void AliNeutralMesonSelection::InitParameters()
177 {
178   
179   //Initialize the parameters of the analysis.
180   fAngleMaxParam.Set(4) ;
181   fAngleMaxParam.Reset(0.);
182
183   SetParticle("Pi0");
184   
185   //Histogrammes settings
186   fHistoNEBins     = 200 ;
187   fHistoEMax       = 50  ;
188   fHistoEMin       = 0.  ;  
189
190   fHistoNAngleBins = 200 ;
191   fHistoAngleMax   = 0.5 ;
192   fHistoAngleMin   = 0.  ;
193
194 }
195
196 //______________________________________________________________________________
197 Bool_t AliNeutralMesonSelection::IsAngleInWindow(Float_t angle, Float_t e) const
198 {
199  
200   // Check if the opening angle of the candidate pairs is inside 
201   // our selection window
202   // Attention, only valid for Pi0, if needed for Eta need to revise max angle function or change parameters
203         
204   Double_t max =  fAngleMaxParam.At(0)*TMath::Exp(fAngleMaxParam.At(1)*e)
205     +fAngleMaxParam.At(2)+fAngleMaxParam.At(3)*e;
206   Double_t arg = (e*e-2*fM*fM)/(e*e);
207   Double_t min = 100. ;
208   if(arg>0.)
209     min = TMath::ACos(arg)+fShiftMinAngle[0]+fShiftMinAngle[1]*e;
210   
211   if((angle<max)&&(angle>=min)) return kTRUE  ;
212   else                          return kFALSE ;
213
214 }
215
216 //_________________________________________________________________
217 Bool_t  AliNeutralMesonSelection::SelectPair(TLorentzVector gammai, 
218                                              TLorentzVector gammaj, 
219                                              TString calo)  
220 {  
221   
222   //Search for the neutral pion within selection cuts
223   
224   //  Double_t pt  = (gammai+gammaj).Pt();
225   Double_t phi = (gammai+gammaj).Phi();
226   if(phi < 0)
227     phi+=TMath::TwoPi();
228   
229   Double_t invmass = (gammai+gammaj).M();
230   Double_t angle   = gammaj.Angle(gammai.Vect());
231   Double_t e       = (gammai+gammaj).E();
232   Double_t asy     = TMath::Abs((gammai-gammaj).E())/(gammai+gammaj).E();
233
234   //Fill histograms with no cuts applied.
235   if(fKeepNeutralMesonHistos)
236   {
237           fhAnglePairNoCut  ->Fill(e,angle);
238           fhInvMassPairNoCut->Fill(e,invmass);
239           fhAsymmetryNoCut  ->Fill(e,asy);
240   }
241   
242   //Cut on the aperture of the pair
243   if(fUseAngleCut)
244   {
245     if(IsAngleInWindow(angle,e))
246     {
247       if(fKeepNeutralMesonHistos )
248       {
249         fhAnglePairOpeningAngleCut  ->Fill(e,angle);
250         fhInvMassPairOpeningAngleCut->Fill(e,invmass);
251         fhAsymmetryOpeningAngleCut  ->Fill(e,asy);
252       }
253       //AliDebug(2,Form("Angle cut: pt %f, phi %f",pt,phi));
254     } else return kFALSE;
255   }
256   
257   // Asymmetry cut
258   if(fUseAsymmetryCut)
259   {
260     if(fAsymmetryCut > asy)
261     {
262       if(fKeepNeutralMesonHistos)
263       {
264         fhInvMassPairAsymmetryCut->Fill(e,invmass);
265         fhAnglePairAsymmetryCut  ->Fill(e,angle);
266       }
267     } else return kFALSE;
268   }
269   
270   
271   //Cut on the invariant mass of the pair
272   
273   Float_t invmassmaxcut = fInvMassMaxCut;
274   Float_t invmassRightBandMinCut = fRightBandMinCut;
275   Float_t invmassRightBandMixCut = fRightBandMaxCut;
276
277   if(calo=="EMCAL" && e > 6.)
278   { // for EMCAL, pi0s, mass depends strongly with energy for e > 6, loose max cut
279   
280     invmassmaxcut = (fInvMassMaxCutParam[0]+fInvMassMaxCut)+fInvMassMaxCutParam[1]*e+fInvMassMaxCutParam[2]*e*e;
281     invmassRightBandMinCut = (fInvMassMaxCutParam[0]+fRightBandMinCut)+fInvMassMaxCutParam[1]*e+fInvMassMaxCutParam[2]*e*e;
282     invmassRightBandMixCut = (fInvMassMaxCutParam[0]+fRightBandMaxCut)+fInvMassMaxCutParam[1]*e+fInvMassMaxCutParam[2]*e*e;
283
284     //printf("e %f, max cut %f, p00 %f,p0 %f,p1 %f,p2 %f\n",
285     //       e,invmassmaxcut,fInvMassMaxCut,fInvMassMaxCutParam[0],fInvMassMaxCutParam[1],fInvMassMaxCutParam[2]);
286   }
287   
288   // normal case, invariant mass selection around pi0/eta peak
289   if( !fParticle.Contains("SideBand") )
290   {
291     if( invmass > fInvMassMinCut && invmass < invmassmaxcut )
292     { 
293       if(fKeepNeutralMesonHistos)
294       {
295         fhInvMassPairAllCut->Fill(e,invmass);
296         fhAnglePairAllCut  ->Fill(e,angle);
297         fhAsymmetryAllCut  ->Fill(e,asy);
298       }      
299       
300       //AliDebug(2,Form("IM cut: pt %f, phi %f",pt,phi));
301       return kTRUE;
302       
303     }//(invmass>0.125) && (invmass<0.145)
304     else 
305     {
306       return kFALSE;
307     }
308   }// normal selection
309  
310   else // select a band around pi0/eta
311   {
312     if((invmass > fLeftBandMinCut  && invmass < fLeftBandMaxCut ) ||  
313        (invmass > invmassRightBandMinCut && invmass < invmassRightBandMixCut))
314     { 
315       if(fKeepNeutralMesonHistos)
316       {
317         fhInvMassPairAllCut->Fill(e,invmass);
318         fhAnglePairAllCut  ->Fill(e,angle);
319         fhAsymmetryAllCut  ->Fill(e,asy);
320       }      
321       
322       //AliDebug(2,Form("IM cut: pt %f, phi %f",pt,phi));
323       return kTRUE;
324       
325     }//(invmass>0.125) && (invmass<0.145)
326     else 
327     {
328       return kFALSE;
329     }
330   }
331   
332 }
333
334 //_______________________________________________________________
335 void  AliNeutralMesonSelection::SetParticle(TString particleName)
336 {
337   // Set some default parameters for selection of pi0 or eta
338
339   fParticle = particleName ;
340   fDecayBit = kPi0;
341
342   if(particleName.Contains("Pi0"))
343   {
344     fHistoNIMBins          = 150 ;
345     fHistoIMMax            = 0.3 ;
346     fHistoIMMin            = 0.  ;  
347     
348     fM                     = 0.135 ; // GeV
349     fInvMassMaxCut         = 0.16  ; // GeV
350     fInvMassMinCut         = 0.11  ; // GeV
351
352     fLeftBandMinCut        = 0.05  ; // GeV
353     fLeftBandMaxCut        = 0.09  ; // GeV
354     fRightBandMinCut       = 0.160 ; // GeV
355     fRightBandMaxCut       = 0.165 ; // GeV
356
357     fInvMassMaxCutParam[0] = 0.0   ;
358     fInvMassMaxCutParam[1] =-7.e-5 ;
359     fInvMassMaxCutParam[2] = 8.e-5 ;    
360          
361     fShiftMinAngle[0]      =-0.03  ;
362     fShiftMinAngle[1]      = 0.0025;
363     
364     fAngleMaxParam.AddAt( 0.8,  0) ;
365     fAngleMaxParam.AddAt(-1,    1) ;
366     fAngleMaxParam.AddAt( 0.09, 2) ; //for pi0 shift, for eta maybe 0.09 
367     fAngleMaxParam.AddAt(-2.e-3,3) ;
368     
369     fDecayBit = kPi0;
370     if(particleName.Contains("Side")) fDecayBit = kPi0Side;
371   
372   }  
373   else if(particleName.Contains("Eta"))
374   {
375     fHistoNIMBins          = 200  ; // GeV
376     fHistoIMMax            = 0.75 ; // GeV
377     fHistoIMMin            = 0.35 ; // GeV 
378   
379     fM                     = 0.547 ; // GeV
380     fInvMassMaxCut         = 0.590 ; // GeV
381     fInvMassMinCut         = 0.510 ; // GeV
382     
383     fLeftBandMinCut        = 0.450 ; // GeV
384     fLeftBandMaxCut        = 0.500 ; // GeV
385     fRightBandMinCut       = 0.600 ; // GeV
386     fRightBandMaxCut       = 0.650 ; // GeV
387     
388     fInvMassMaxCutParam[0] = 0.00 ;
389     fInvMassMaxCutParam[1] = 0.00 ;
390     fInvMassMaxCutParam[2] = 0.00 ;
391     
392     fShiftMinAngle[0]      =-0.03 ;
393     fShiftMinAngle[0]      = 0.00 ;
394     
395     fAngleMaxParam.AddAt( 0.80,  0) ; // Same as pi0
396     fAngleMaxParam.AddAt(-0.25,  1) ; // Same as pi0
397     fAngleMaxParam.AddAt( 0.12,  2) ; // Shifted with respect to pi0
398     fAngleMaxParam.AddAt(-5.e-4, 3) ; // Same as pi0
399     
400     fDecayBit = kEta;
401     if(particleName.Contains("Side")) fDecayBit = kEtaSide;
402   }
403   else 
404     printf("AliAnaNeutralMesonSelection::SetParticle(%s) *** Particle NOT defined (Pi0 or Eta), Pi0 settings by default *** \n",particleName.Data());
405   
406 }
407
408 //______________________________________________________________
409 void AliNeutralMesonSelection::Print(const Option_t * opt) const
410 {
411
412   //Print some relevant parameters set for the analysis
413   if(! opt)
414     return;
415   
416   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
417
418   printf("Particle %s, bit %d, mass : %f  \n", fParticle.Data(), fDecayBit, fM );
419   printf("Invariant mass limits : %f < m < %f \n", fInvMassMinCut , fInvMassMinCut );
420   
421   printf("Use asymmetry cut? : %d ; A < %f \n", fUseAngleCut, fAsymmetryCut );
422   
423   printf("Use angle cut? : %d  \n", fUseAngleCut );
424   if(fUseAngleCut)
425   {
426     printf("Angle selection param: \n");
427     printf("p0 :     %f\n", fAngleMaxParam.At(0));
428     printf("p1 :     %f\n", fAngleMaxParam.At(1));
429     printf("p2 :     %f\n", fAngleMaxParam.At(2));
430     printf("p3 :     %f\n", fAngleMaxParam.At(3));
431     printf("Min angle shift : p0 = %1.3f, p1 = %1.3f,\n", fShiftMinAngle[0],fShiftMinAngle[1]);
432   }
433   
434   printf("Keep Neutral Meson Histos = %d\n",fKeepNeutralMesonHistos);
435   
436   if(fKeepNeutralMesonHistos)
437   {
438     printf("Histograms: %3.1f < E  < %3.1f,  Nbin = %d\n",   fHistoEMin,     fHistoEMax,     fHistoNEBins);
439     printf("Histograms: %3.1f < angle < %3.1f, Nbin = %d\n", fHistoAngleMin, fHistoAngleMax, fHistoNAngleBins);
440     printf("Histograms: %3.1f < IM < %3.1f, Nbin = %d\n",    fHistoIMMin,    fHistoIMMax,    fHistoNIMBins);    
441   }
442   
443