New sub-structure for fitting energy loss spectra with
[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) 
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}@ is the Landau distribution, @f$ \Delta@f$ the
62    * energy loss, @f$ \xi@f the width of the Landau, and @f$
63    * \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$, @\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     /** 
125      * Constructor 
126      * 
127      * @param lowCut     Lower cut of spectrum - data below this cuts is ignored
128      * @param maxRange   Maximum range to fit to 
129      * @param minusBins  The number of bins below maximum to use 
130      */
131     ELossFitter(Double_t lowCut, Double_t maxRange, UShort_t minusBins); 
132     virtual ~ELossFitter();
133     /** 
134      * Clear internal arrays 
135      * 
136      */
137     void Clear();
138     /** 
139      * Fit a 1-particle signal to the passed energy loss distribution 
140      * 
141      * Note that this function clears the internal arrays first 
142      *
143      * @param dist    Data to fit the function to 
144      * @param sigman If larger than zero, the initial guess of the
145      *               detector induced noise. If zero or less, then this 
146      *               parameter is ignored in the fit (fixed at 0)
147      * 
148      * @return The function fitted to the data 
149      */
150     TF1* Fit1Particle(TH1* dist, Double_t sigman=-1);
151     /** 
152      * Fit a N-particle signal to the passed energy loss distribution 
153      *
154      * If there's no 1-particle fit present, it does that first 
155      *
156      * @param dist   Data to fit the function to 
157      * @param N      Number of particle signals to fit 
158      * @param sigman If larger than zero, the initial guess of the
159      *               detector induced noise. If zero or less, then this 
160      *               parameter is ignored in the fit (fixed at 0)
161      * 
162      * @return The function fitted to the data 
163      */
164     TF1* FitNParticle(TH1* dist, UShort_t n, Double_t sigman=-1);
165      
166
167     const Double_t fLowCut;     // Lower cut on data 
168     const Double_t fMaxRange;   // Maximum range to fit 
169     const UShort_t fMinusBins;  // Number of bins from maximum to fit 1st peak
170     TObjArray fFitResults;      // Array of fit results 
171     TObjArray fFunctions;       // Array of functions 
172   };
173       
174
175   //__________________________________________________________________
176   /** 
177    * Structure to hold histograms 
178    *
179    * @ingroup pwg2_forward_analysis 
180    */
181   struct Histos : public TObject
182   {     
183     /** 
184      * Constructor 
185      * 
186      * 
187      */
188     Histos() : fFMD1i(0), fFMD2i(0), fFMD2o(0), fFMD3i(0), fFMD3o(0) {}
189     /** 
190      * Copy constructor 
191      * 
192      * @param o Object to copy from 
193      */
194     Histos(const Histos& o) 
195       : TObject(o), 
196         fFMD1i(o.fFMD1i), 
197         fFMD2i(o.fFMD2i), 
198         fFMD2o(o.fFMD2o), 
199         fFMD3i(o.fFMD3i), 
200         fFMD3o(o.fFMD3o)
201     {}
202     /** 
203      * Assignement operator 
204      * 
205      * @return Reference to this 
206      */
207     Histos& operator=(const Histos&) { return *this;}
208     /** 
209      * Destructor
210      */
211     ~Histos();
212     /** 
213      * Initialize the object 
214      * 
215      * @param etaAxis Eta axis to use 
216      */
217     void Init(const TAxis& etaAxis);
218     /** 
219      * Make a histogram 
220      * 
221      * @param d        Detector
222      * @param r        Ring 
223      * @param etaAxis  Eta axis to use
224      * 
225      * @return Newly allocated histogram 
226      */
227     TH2D* Make(UShort_t d, Char_t r, const TAxis& etaAxis) const;
228     /** 
229      * Clear data 
230      * 
231      * @param option Not used 
232      */
233     void  Clear(Option_t* option="");
234     // const TH2D* Get(UShort_t d, Char_t r) const;
235     /** 
236      * Get the histogram for a particular detector,ring
237      * 
238      * @param d Detector 
239      * @param r Ring 
240      * 
241      * @return Histogram for detector,ring or nul 
242      */
243     TH2D* Get(UShort_t d, Char_t r) const;
244     TH2D* fFMD1i; // Histogram for FMD1i
245     TH2D* fFMD2i; // Histogram for FMD2i
246     TH2D* fFMD2o; // Histogram for FMD2o
247     TH2D* fFMD3i; // Histogram for FMD3i
248     TH2D* fFMD3o; // Histogram for FMD3o
249
250     ClassDef(Histos,1) 
251   };
252
253   //__________________________________________________________________
254   struct RingHistos : public TObject
255   {
256     RingHistos() : fDet(0), fRing('\0'), fName("") {}
257     RingHistos(UShort_t d, Char_t r) 
258       : fDet(d), fRing(r), fName(TString::Format("FMD%d%c", d, r)) 
259     {}
260     RingHistos(const RingHistos& o) 
261       : TObject(o), fDet(o.fDet), fRing(o.fRing), fName(o.fName)
262     {}
263     virtual ~RingHistos() {}
264     RingHistos& operator=(const RingHistos& o) 
265     {
266       TObject::operator=(o);
267       fDet  = o.fDet;
268       fRing = o.fRing;
269       fName = o.fName;
270       return *this;
271     }
272     TList* DefineOutputList(TList* d) const;
273     TList* GetOutputList(TList* d) const;
274     TH1* GetOutputHist(TList* d, const char* name) const;
275     Color_t Color() const 
276     { 
277       return ((fDet == 1 ? kRed : (fDet == 2 ? kGreen : kBlue))
278               + ((fRing == 'I' || fRing == 'i') ? 2 : -2));
279     }
280
281     UShort_t fDet; 
282     Char_t   fRing;
283     TString  fName;
284
285     ClassDef(RingHistos,1) 
286   };
287     
288 };
289
290 #endif
291 // Local Variables:
292 //  mode: C++
293 // End:
294