Improved eloss fitting - see NIM B1, 16
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliForwardUtil.h
1 #ifndef ALIROOT_PWG2_FORWARD_ALIFORWARDUTIL_H
2 #define ALIROOT_PWG2_FORWARD_ALIFORWARDUTIL_H
3 #include <TObject.h>
4 #include <TString.h>
5 #include <TObjArray.h>
6 class TH2D;
7 class TH1I;
8 class TH1;
9 class TF1;
10 class TAxis;
11 class AliESDEvent;
12
13 /** 
14  * Utilities used in the forward multiplcity analysis 
15  * 
16  * @ingroup pwg2_forward_analysis 
17  */
18 class AliForwardUtil : public TObject
19 {
20 public:
21   //__________________________________________________________________
22   /**
23    * Number of steps to do in the Landau, Gaussiam convolution 
24    */
25   static Int_t fgConvolutionSteps;
26   //------------------------------------------------------------------
27   /** 
28    * How many sigma's of the Gaussian in the Landau, Gaussian
29    * convolution to integrate over
30    */
31   static Double_t fgConvolutionNSigma;
32   //------------------------------------------------------------------
33   /** 
34    * Calculate the shifted Landau
35    * @f[
36    *    f'_{L}(x;\Delta,\xi) = f_L(x;\Delta+0.22278298\xi)
37    * @f]
38    *
39    * where @f$ f_{L}@f$ is the ROOT implementation of the Landau
40    * distribution (known to have @f$ \Delta_{p}=-0.22278298@f$ for
41    * @f$\Delta=0,\xi=1@f$. 
42    *
43    * @param x      Where to evaluate @f$ f'_{L}@f$ 
44    * @param delta  Most probable value 
45    * @param xi     The 'width' of the distribution 
46    *
47    * @return @f$ f'_{L}(x;\Delta,\xi) @f$
48    */
49   static Double_t Landau(Double_t x, Double_t delta, Double_t xi);
50   
51   //------------------------------------------------------------------
52   /** 
53    * Calculate the value of a Landau convolved with a Gaussian 
54    * 
55    * @f[ 
56    * f(x;\Delta,\xi,\sigma') = \frac{1}{\sigma' \sqrt{2 \pi}}
57    *    \int_{-\infty}^{+\infty} d\Delta' f'_{L}(x;\Delta',\xi)
58    *    \exp{-\frac{(\Delta-\Delta')^2}{2\sigma'^2}}
59    * @f]
60    * 
61    * where @f$ f'_{L}@f$ is the Landau distribution, @f$ \Delta@f$ the
62    * energy loss, @f$ \xi@f$ the width of the Landau, and 
63    * @f$ \sigma'^2=\sigma^2-\sigma_n^2 @f$.  Here, @f$\sigma@f$ is the
64    * variance of the Gaussian, and @f$\sigma_n@f$ is a parameter modelling 
65    * noise in the detector.  
66    *
67    * Note that this function uses the constants fgConvolutionSteps and
68    * fgConvolutionNSigma
69    * 
70    * References: 
71    *  - <a href="http://dx.doi.org/10.1016/0168-583X(84)90472-5">Nucl.Instrum.Meth.B1:16</a>
72    *  - <a href="http://dx.doi.org/10.1103/PhysRevA.28.615">Phys.Rev.A28:615</a>
73    *  - <a href="http://root.cern.ch/root/htmldoc/tutorials/fit/langaus.C.html">ROOT implementation</a>
74    * 
75    * @param x         where to evaluate @f$ f@f$
76    * @param delta     @f$ \Delta@f$ of @f$ f(x;\Delta,\xi,\sigma')@f$
77    * @param xi        @f$ \xi@f$ of @f$ f(x;\Delta,\xi,\sigma')@f$
78    * @param sigma     @f$ \sigma@f$ of @f$\sigma'^2=\sigma^2-\sigma_n^2 @f$
79    * @param sigma_n   @f$ \sigma_n@f$ of @f$\sigma'^2=\sigma^2-\sigma_n^2 @f$
80    * 
81    * @return @f$ f@f$ evaluated at @f$ x@f$.  
82    */
83   static Double_t LandauGaus(Double_t x, Double_t delta, Double_t xi, 
84                              Double_t sigma, Double_t sigma_n);
85   
86   //------------------------------------------------------------------
87   /** 
88    * Evaluate 
89    * @f[ 
90       f_N(x;\Delta,\xi,\sigma') = \sum_{i=1}^N a_i f(x;\Delta_i,\xi_i,\sigma'_i)
91      @f] 
92    * 
93    * where @f$ f(x;\Delta,\xi,\sigma')@f$ is the convolution of a
94    * Landau with a Gaussian (see LandauGaus).  Note that 
95    * @f$ a_1 = 1@f$, @f$\Delta_i = i(\Delta_1 + \xi\log(i))@f$, 
96    * @f$\xi_i=i\xi_1@f$, and @f$\sigma_i'^2 = \sigma_n^2 + i\sigma_1^2@f$. 
97    *  
98    * References: 
99    *  - <a href="http://dx.doi.org/10.1016/0168-583X(84)90472-5">Nucl.Instrum.Meth.B1:16</a>
100    *  - <a href="http://dx.doi.org/10.1103/PhysRevA.28.615">Phys.Rev.A28:615</a>
101    *  - <a href="http://root.cern.ch/root/htmldoc/tutorials/fit/langaus.C.html">ROOT implementation</a>
102    * 
103    * @param x        Where to evaluate @f$ f_N@f$
104    * @param delta    @f$ \Delta_1@f$ 
105    * @param xi       @f$ \xi_1@f$
106    * @param sigma    @f$ \sigma_1@f$ 
107    * @param sigma_n  @f$ \sigma_n@f$ 
108    * @param n        @f$ N@f$ in the sum above.
109    * @param a        Array of size @f$ N-1@f$ of the weights @f$ a_i@f$ for 
110    *                 @f$ i > 1@f$ 
111    * 
112    * @return @f$ f_N(x;\Delta,\xi,\sigma')@f$ 
113    */
114   static Double_t NLandauGaus(Double_t x, Double_t delta, Double_t xi, 
115                               Double_t sigma, Double_t sigma_n, Int_t n, 
116                               Double_t* a);
117   //__________________________________________________________________
118   /** 
119    * Structure to do fits to the energy loss spectrum 
120    * 
121    */
122   struct ELossFitter 
123   {
124     enum { 
125       kC     = 0,
126       kDelta, 
127       kXi, 
128       kSigma, 
129       kSigmaN, 
130       kN, 
131       kA
132     };
133     /** 
134      * Constructor 
135      * 
136      * @param lowCut     Lower cut of spectrum - data below this cuts is ignored
137      * @param maxRange   Maximum range to fit to 
138      * @param minusBins  The number of bins below maximum to use 
139      */
140     ELossFitter(Double_t lowCut, Double_t maxRange, UShort_t minusBins); 
141     virtual ~ELossFitter();
142     /** 
143      * Clear internal arrays 
144      * 
145      */
146     void Clear();
147     /** 
148      * Fit a 1-particle signal to the passed energy loss distribution 
149      * 
150      * Note that this function clears the internal arrays first 
151      *
152      * @param dist    Data to fit the function to 
153      * @param sigman If larger than zero, the initial guess of the
154      *               detector induced noise. If zero or less, then this 
155      *               parameter is ignored in the fit (fixed at 0)
156      * 
157      * @return The function fitted to the data 
158      */
159     TF1* Fit1Particle(TH1* dist, Double_t sigman=-1);
160     /** 
161      * Fit a N-particle signal to the passed energy loss distribution 
162      *
163      * If there's no 1-particle fit present, it does that first 
164      *
165      * @param dist   Data to fit the function to 
166      * @param n      Number of particle signals to fit 
167      * @param sigman If larger than zero, the initial guess of the
168      *               detector induced noise. If zero or less, then this 
169      *               parameter is ignored in the fit (fixed at 0)
170      * 
171      * @return The function fitted to the data 
172      */
173     TF1* FitNParticle(TH1* dist, UShort_t n, Double_t sigman=-1);
174      
175
176     const Double_t fLowCut;     // Lower cut on data 
177     const Double_t fMaxRange;   // Maximum range to fit 
178     const UShort_t fMinusBins;  // Number of bins from maximum to fit 1st peak
179     TObjArray fFitResults;      // Array of fit results 
180     TObjArray fFunctions;       // Array of functions 
181   };
182       
183
184   //__________________________________________________________________
185   /** 
186    * Structure to hold histograms 
187    *
188    * @ingroup pwg2_forward_analysis 
189    */
190   struct Histos : public TObject
191   {     
192     /** 
193      * Constructor 
194      * 
195      * 
196      */
197     Histos() : fFMD1i(0), fFMD2i(0), fFMD2o(0), fFMD3i(0), fFMD3o(0) {}
198     /** 
199      * Copy constructor 
200      * 
201      * @param o Object to copy from 
202      */
203     Histos(const Histos& o) 
204       : TObject(o), 
205         fFMD1i(o.fFMD1i), 
206         fFMD2i(o.fFMD2i), 
207         fFMD2o(o.fFMD2o), 
208         fFMD3i(o.fFMD3i), 
209         fFMD3o(o.fFMD3o)
210     {}
211     /** 
212      * Assignement operator 
213      * 
214      * @return Reference to this 
215      */
216     Histos& operator=(const Histos&) { return *this;}
217     /** 
218      * Destructor
219      */
220     ~Histos();
221     /** 
222      * Initialize the object 
223      * 
224      * @param etaAxis Eta axis to use 
225      */
226     void Init(const TAxis& etaAxis);
227     /** 
228      * Make a histogram 
229      * 
230      * @param d        Detector
231      * @param r        Ring 
232      * @param etaAxis  Eta axis to use
233      * 
234      * @return Newly allocated histogram 
235      */
236     TH2D* Make(UShort_t d, Char_t r, const TAxis& etaAxis) const;
237     /** 
238      * Clear data 
239      * 
240      * @param option Not used 
241      */
242     void  Clear(Option_t* option="");
243     // const TH2D* Get(UShort_t d, Char_t r) const;
244     /** 
245      * Get the histogram for a particular detector,ring
246      * 
247      * @param d Detector 
248      * @param r Ring 
249      * 
250      * @return Histogram for detector,ring or nul 
251      */
252     TH2D* Get(UShort_t d, Char_t r) const;
253     TH2D* fFMD1i; // Histogram for FMD1i
254     TH2D* fFMD2i; // Histogram for FMD2i
255     TH2D* fFMD2o; // Histogram for FMD2o
256     TH2D* fFMD3i; // Histogram for FMD3i
257     TH2D* fFMD3o; // Histogram for FMD3o
258
259     ClassDef(Histos,1) 
260   };
261
262   //__________________________________________________________________
263   struct RingHistos : public TObject
264   {
265     RingHistos() : fDet(0), fRing('\0'), fName("") {}
266     RingHistos(UShort_t d, Char_t r) 
267       : fDet(d), fRing(r), fName(TString::Format("FMD%d%c", d, r)) 
268     {}
269     RingHistos(const RingHistos& o) 
270       : TObject(o), fDet(o.fDet), fRing(o.fRing), fName(o.fName)
271     {}
272     virtual ~RingHistos() {}
273     RingHistos& operator=(const RingHistos& o) 
274     {
275       TObject::operator=(o);
276       fDet  = o.fDet;
277       fRing = o.fRing;
278       fName = o.fName;
279       return *this;
280     }
281     TList* DefineOutputList(TList* d) const;
282     TList* GetOutputList(TList* d) const;
283     TH1* GetOutputHist(TList* d, const char* name) const;
284     Color_t Color() const 
285     { 
286       return ((fDet == 1 ? kRed : (fDet == 2 ? kGreen : kBlue))
287               + ((fRing == 'I' || fRing == 'i') ? 2 : -2));
288     }
289
290     UShort_t fDet; 
291     Char_t   fRing;
292     TString  fName;
293
294     ClassDef(RingHistos,1) 
295   };
296     
297 };
298
299 #endif
300 // Local Variables:
301 //  mode: C++
302 // End:
303