More code clean up.
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliFMDEnergyFitter.h
1 #ifndef ALIROOT_PWG2_FORWARD_ALIFMDENERGYFITTER_H
2 #define ALIROOT_PWG2_FORWARD_ALIFMDENERGYFITTER_H
3 #include <TNamed.h>
4 #include <TH1D.h>
5 #include <TAxis.h>
6 #include <TList.h>
7 #include <TObjArray.h>
8 #include <TClonesArray.h>
9 #include "AliFMDCorrELossFit.h"
10 #include "AliForwardUtil.h"
11 class AliESDFMD;
12 class TFitResult;
13 class TF1;
14 class TArrayD;
15
16 /**
17  * Class to fit the energy distribution.  
18  *
19  * @par Input: 
20  *    - AliESDFMD object  - from reconstruction
21  *
22  * @par Output: 
23  *    - Lists of histogram - one per ring.  Each list has a number of 
24  *      histograms corresponding to the number of eta bins defined.  
25  *
26  * @par Corrections used: 
27  *    - None
28  *
29  *
30  * @ingroup pwg2_forward_analysis 
31  */
32 class AliFMDEnergyFitter : public TNamed
33 {
34 public: 
35     enum { 
36       kC        = AliForwardUtil::ELossFitter::kC,
37       kDelta    = AliForwardUtil::ELossFitter::kDelta, 
38       kXi       = AliForwardUtil::ELossFitter::kXi, 
39       kSigma    = AliForwardUtil::ELossFitter::kSigma, 
40       kSigmaN   = AliForwardUtil::ELossFitter::kSigmaN,
41       kN        = AliForwardUtil::ELossFitter::kN, 
42       kA        = AliForwardUtil::ELossFitter::kA
43     };
44
45   /** 
46    * Destructor
47    */
48   virtual ~AliFMDEnergyFitter();
49   /** 
50    * Default Constructor - do not use 
51    */
52   AliFMDEnergyFitter();
53   /** 
54    * Constructor 
55    * 
56    * @param title Title of object  - not significant 
57    */
58   AliFMDEnergyFitter(const char* title);
59   /** 
60    * Copy constructor 
61    * 
62    * @param o Object to copy from 
63    */
64   AliFMDEnergyFitter(const AliFMDEnergyFitter& o);
65   /** 
66    * Assignment operator 
67    * 
68    * @param o Object to assign from 
69    * 
70    * @return Reference to this 
71    */
72   AliFMDEnergyFitter& operator=(const AliFMDEnergyFitter& o);
73
74   /** 
75    * Initialise the task
76    * 
77    * @param etaAxis The eta axis to use.  Note, that if the eta axis
78    * has already been set (using SetEtaAxis), then this parameter will be 
79    * ignored
80    */
81   void Init(const TAxis& etaAxis);
82   /** 
83    * Set the eta axis to use.  This will force the code to use this
84    * eta axis definition - irrespective of whatever axis is passed to
85    * the Init member function.  Therefore, this member function can be
86    * used to force another eta axis than one found in the correction
87    * objects. 
88    * 
89    * @param nBins  Number of bins 
90    * @param etaMin Minimum of the eta axis 
91    * @param etaMax Maximum of the eta axis 
92    */
93   void SetEtaAxis(Int_t nBins, Double_t etaMin, Double_t etaMax);
94   /** 
95    * Set the eta axis to use.  This will force the code to use this
96    * eta axis definition - irrespective of whatever axis is passed to
97    * the Init member function.  Therefore, this member function can be
98    * used to force another eta axis than one found in the correction
99    * objects. 
100    * 
101    * @param etaAxis Eta axis to use 
102    */
103   void SetEtaAxis(const TAxis& etaAxis);
104   /** 
105    * Set the low cut used for energy 
106    * 
107    * @param lowCut Low cut
108    */
109   void SetLowCut(Double_t lowCut=0.3) { fLowCut = lowCut; }
110   /** 
111    * Set the number of bins to subtract 
112    * 
113    * @param n 
114    */
115   void SetFitRangeBinWidth(UShort_t n=4) { fFitRangeBinWidth = n; }
116   /** 
117    * Whether or not to enable fitting of the final merged result.  
118    * Note, fitting takes quite a while and one should be careful not to do 
119    * this needlessly 
120    * 
121    * @param doFit Whether to do the fits or not 
122    */
123   void SetDoFits(Bool_t doFit=kTRUE) { fDoFits = doFit; }
124   /** 
125    * Set whether to make the corrections object on the output.  Note,
126    * fits should be enable for this to have any effect.
127    * 
128    * @param doMake If true (false is default), do make the corrections object. 
129    */
130   void SetDoMakeObject(Bool_t doMake=kTRUE) { fDoMakeObject = doMake; }
131   /** 
132    * Set how many particles we will try to fit at most to the data
133    * 
134    * @param n Max number of particle to try to fit 
135    */
136   void SetNParticles(UShort_t n) { fNParticles = (n < 1 ? 1 : (n > 5 ? 5 : n)); }
137   /** 
138    * Set the minimum number of entries each histogram must have 
139    * before we try to fit our response function to it
140    * 
141    * @param n Minimum number of entries
142    */
143   void SetMinEntries(UShort_t n) { fMinEntries = (n < 1 ? 1 : n); }
144   /**
145    * Set maximum energy loss to consider 
146    *
147    * @param x Maximum energy loss to consider 
148    */
149   void SetMaxE(Double_t x) { fMaxE = x; }
150   /**
151    * Set number of energy loss bins 
152    *
153    * @param x Number of energy loss bins 
154    */
155   void SetNEbins(Int_t x) { fNEbins = x; }
156   void SetMaxRelativeParameterError(Double_t e=0.2) { fMaxRelParError = e; }
157   void SetMaxChi2PerNDF(Double_t c=10) { fMaxChi2PerNDF = c; }
158   void SetMinWeight(Double_t c=1e-7) { fMinWeight = c; }
159   /**
160    * Set wheter to use increasing bin sizes 
161    *
162    * @param x Wheter to use increasing bin sizes 
163    */
164   void SetUseIncreasingBins(Bool_t x) { fUseIncreasingBins = x; }
165   /** 
166    * Fitter the input AliESDFMD object
167    * 
168    * @param input     Input 
169    * @param empty     Whether the event is 'empty'
170    * 
171    * @return True on success, false otherwise 
172    */
173   Bool_t Accumulate(const AliESDFMD& input, 
174                     Bool_t           empty);
175   /** 
176    * Scale the histograms to the total number of events 
177    * 
178    * @param dir Where the histograms are  
179    */
180   void Fit(TList* dir);
181   void MakeCorrectionsObject(TList* dir);
182   
183   /** 
184    * Define the output histograms.  These are put in a sub list of the
185    * passed list.   The histograms are merged before the parent task calls 
186    * AliAnalysisTaskSE::Terminate 
187    * 
188    * @param dir Directory to add to 
189    */
190   void DefineOutput(TList* dir);
191   /** 
192    * Set the debug level.  The higher the value the more output 
193    * 
194    * @param dbg Debug level 
195    */
196   void SetDebug(Int_t dbg=1);
197   /** 
198    * Print information
199    * 
200    * @param option Not used 
201    */
202   void Print(Option_t* option="") const;
203 protected:
204   /** 
205    * Internal data structure to keep track of the histograms
206    */
207   struct RingHistos : public AliForwardUtil::RingHistos
208   { 
209     /** 
210      * Default CTOR
211      */
212     RingHistos();
213     /** 
214      * Constructor
215      * 
216      * @param d detector
217      * @param r ring 
218      */
219     RingHistos(UShort_t d, Char_t r);
220     /** 
221      * Copy constructor 
222      * 
223      * @param o Object to copy from 
224      */
225     RingHistos(const RingHistos& o);
226     /** 
227      * Assignment operator 
228      * 
229      * @param o Object to assign from 
230      * 
231      * @return Reference to this 
232      */
233     RingHistos& operator=(const RingHistos& o);
234     /** 
235      * Destructor 
236      */
237     ~RingHistos();
238     /** 
239      * Define outputs
240      * 
241      * @param dir 
242      */
243     void Output(TList* dir);
244     /** 
245      * Initialise object 
246      * 
247      * @param eAxis      Eta axis
248      * @param maxDE      Max energy loss to consider 
249      * @param nDEbins    Number of bins 
250      * @param useIncrBin Whether to use an increasing bin size 
251      */
252     void Init(const TAxis& eAxis, 
253               Double_t     maxDE=10, 
254               Int_t        nDEbins=300, 
255               Bool_t       useIncrBin=true);
256     /** 
257      * Fill histogram 
258      * 
259      * @param empty  True if event is empty
260      * @param ieta   Eta bin
261      * @param mult   Signal 
262      */
263     void Fill(Bool_t empty, Int_t ieta, Double_t mult);
264     /** 
265      * Fit each histogram to up to @a nParticles particle responses.
266      * 
267      * @param dir         Output list 
268      * @param eta         Eta axis 
269      * @param lowCut      Lower cut 
270      * @param nParticles  Max number of convolved landaus to fit
271      * @param minEntries  Minimum number of entries 
272      * @param minusBins   Number of bins from peak to subtract to 
273      *                    get the fit range 
274      * @param relErrorCut Cut applied to relative error of parameter. 
275      *                    Note, for multi-particle weights, the cut 
276      *                    is loosend by a factor of 2 
277      * @param chi2nuCut   Cut on @f$ \chi^2/\nu@f$ - 
278      *                    the reduced @f$\chi^2@f$ 
279      */
280     TObjArray* Fit(TList* dir, 
281                    const TAxis& eta,
282                    Double_t     lowCut, 
283                    UShort_t     nParticles,
284                    UShort_t     minEntries,
285                    UShort_t     minusBins,
286                    Double_t     relErrorCut, 
287                    Double_t     chi2nuCut) const;
288     /** 
289      * Fit a signal histogram.  First, the bin @f% b_{min}@f$ with
290      * maximum bin content in the range @f$ [E_{min},\infty]@f$ is
291      * found.  Then the fit range is set to the bin range 
292      * @f$ [b_{min}-\Delta b,b_{min}+2\Delta b]@f$, and a 1 
293      * particle signal is fitted to that.  The parameters of that fit 
294      * is then used as seeds for a fit of the @f$ N@f$ particle response 
295      * to the data in the range 
296      * @f$ [b_{min}-\Delta b,N(\Delta_1+\xi_1\log(N))+2N\xi@f$
297      * 
298      * @param dist        Histogram to fit 
299      * @param lowCut      Lower cut @f$ E_{min}@f$ on signal 
300      * @param nParticles  Max number @f$ N@f$ of convolved landaus to fit
301      * @param minusBins   Number of bins @f$ \Delta b@f$ from peak to 
302      *                    subtract to get the fit range 
303      * @param relErrorCut Cut applied to relative error of parameter. 
304      *                    Note, for multi-particle weights, the cut 
305      *                    is loosend by a factor of 2 
306      * @param chi2nuCut   Cut on @f$ \chi^2/\nu@f$ - 
307      *                    the reduced @f$\chi^2@f$ 
308      * 
309      * @return The best fit function 
310      */
311     TF1* FitHist(TH1*     dist,
312                  Double_t lowCut, 
313                  UShort_t nParticles,
314                  UShort_t minusBins,
315                  Double_t relErrorCut, 
316                  Double_t chi2nuCut) const;
317     /** 
318      * Find the best fits 
319      * 
320      * @param d              Parent list
321      * @param obj            Object to add fits to
322      * @param eta            Eta axis 
323      * @param relErrorCut    Cut applied to relative error of parameter. 
324      *                       Note, for multi-particle weights, the cut 
325      *                       is loosend by a factor of 2 
326      * @param chi2nuCut      Cut on @f$ \chi^2/\nu@f$ - 
327      *                       the reduced @f$\chi^2@f$ 
328      * @param minWeightCut   Least valid @f$ a_i@f$ 
329      */
330     void FindBestFits(TList*              d, 
331                       AliFMDCorrELossFit& obj,
332                       const TAxis&        eta,     
333                       Double_t            relErrorCut, 
334                       Double_t            chi2nuCut,
335                       Double_t            minWeightCut);
336     /** 
337      * Find the best fit 
338      * 
339      * @param dist           Histogram 
340      * @param relErrorCut    Cut applied to relative error of parameter. 
341      *                       Note, for multi-particle weights, the cut 
342      *                       is loosend by a factor of 2 
343      * @param chi2nuCut      Cut on @f$ \chi^2/\nu@f$ - 
344      *                       the reduced @f$\chi^2@f$ 
345      * @param minWeightCut   Least valid @f$ a_i@f$ 
346      * 
347      * @return Best fit 
348      */
349     AliFMDCorrELossFit::ELossFit* FindBestFit(TH1* dist,
350                                               Double_t relErrorCut, 
351                                               Double_t chi2nuCut,
352                                               Double_t minWeightCut);
353     /** 
354      * Check the result of the fit. Returns true if 
355      * - @f$ \chi^2/\nu < \max{\chi^2/\nu}@f$
356      * - @f$ \Delta p_i/p_i < \delta_e@f$ for all parameters.  Note, 
357      *   for multi-particle fits, this requirement is relaxed by a 
358      *   factor of 2
359      * - @f$ a_{n} > 10^{-7}@f$ when fitting to an @f$ n@f$ 
360      *   particle response 
361      * 
362      * @param r           Result to check
363      * @param relErrorCut Cut @f$ \delta_e@f$ applied to relative error 
364      *                    of parameter.  
365      * @param chi2nuCut   Cut @f$ \max{\chi^2/\nu}@f$ 
366      * 
367      * @return true if fit is good. 
368      */
369     Bool_t CheckResult(TFitResult* r,
370                        Double_t    relErrorCut, 
371                        Double_t    chi2nuCut) const;
372     /** 
373      * Make an axis with increasing bins 
374      * 
375      * @param n    Number of bins 
376      * @param min  Minimum 
377      * @param max  Maximum
378      * 
379      * @return An axis with quadratically increasing bin size 
380      */
381     TArrayD MakeIncreasingAxis(Int_t n, Double_t min, Double_t max) const;
382     /** 
383      * Make E/E_mip histogram 
384      * 
385      * @param ieta   Eta bin
386      * @param eMin   Least signal
387      * @param eMax   Largest signal 
388      */
389     void Make(Int_t ieta, Double_t eMin, Double_t eMax, 
390               Double_t deMax=12, Int_t nDeBins=300, Bool_t incr=true);
391     /** 
392      * Make a parameter histogram
393      * 
394      * @param name   Name of histogram.
395      * @param title  Title of histogram. 
396      * @param eta    Eta axis 
397      * 
398      * @return 
399      */
400     TH1D* MakePar(const char* name, const char* title, const TAxis& eta) const;
401     /** 
402      * Make a histogram that contains the results of the fit over the full ring 
403      * 
404      * @param name  Name 
405      * @param title Title
406      * @param eta   Eta axis 
407      * @param low   Least bin
408      * @param high  Largest bin
409      * @param val   Value of parameter 
410      * @param err   Error on parameter 
411      * 
412      * @return The newly allocated histogram 
413      */
414     TH1D* MakeTotal(const char* name, 
415                     const char* title, 
416                     const TAxis& eta, 
417                     Int_t low, 
418                     Int_t high, 
419                     Double_t val, 
420                     Double_t err) const;
421     TH1D*        fEDist;        // Ring energy distribution 
422     TH1D*        fEmpty;        // Ring energy distribution for empty events
423     TList        fEtaEDists;    // Energy distributions per eta bin. 
424     TList*       fList;
425     TClonesArray fFits;
426     Int_t        fDebug;
427     ClassDef(RingHistos,1);
428   };
429   /** 
430    * Get the ring histogram container 
431    * 
432    * @param d Detector
433    * @param r Ring 
434    * 
435    * @return Ring histogram container 
436    */
437   RingHistos* GetRingHistos(UShort_t d, Char_t r) const;
438
439   TList    fRingHistos;    // List of histogram containers
440   Double_t fLowCut;        // Low cut on energy
441   UShort_t fNParticles;    // Number of landaus to try to fit 
442   UShort_t fMinEntries;    // Minimum number of entries
443   UShort_t fFitRangeBinWidth;// Number of bins to subtract from found max
444   Bool_t   fDoFits;        // Whether to actually do the fits 
445   Bool_t   fDoMakeObject;  // Whether to make corrections object
446   TAxis    fEtaAxis;       // Eta axis 
447   Double_t fMaxE;          // Maximum energy loss to consider 
448   Int_t    fNEbins;        // Number of energy loss bins 
449   Bool_t   fUseIncreasingBins; // Wheter to use increasing bin sizes 
450   Double_t fMaxRelParError;// Relative error cut
451   Double_t fMaxChi2PerNDF; // chi^2/nu cit
452   Double_t fMinWeight;     // Minimum weight value 
453   Int_t    fDebug;         // Debug level 
454   
455
456   ClassDef(AliFMDEnergyFitter,1); //
457 };
458
459 #endif
460 // Local Variables:
461 //  mode: C++ 
462 // End: