[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliFMDEnergyFitter.h

#ifndef ALIROOT_PWG2_FORWARD_ALIFMDENERGYFITTER_H
#define ALIROOT_PWG2_FORWARD_ALIFMDENERGYFITTER_H
#include <TNamed.h>
#include <TH1D.h>
#include <TAxis.h>
#include <TList.h>
#include <TObjArray.h>
#include "AliForwardUtil.h"
class AliESDFMD;
class TFitResult;
class TF1;
class TArrayD;

/**
 * Class to fit the energy distribution.  
 *
 * @par Input: 
 *    - AliESDFMD object  - from reconstruction
 *
 * @par Output: 
 *    - Lists of histogram - one per ring.  Each list has a number of 
 *      histograms corresponding to the number of eta bins defined.  
 *
 * @par Corrections used: 
 *    - None
 *
 *
 * @ingroup pwg2_forward_analysis 
 */
class AliFMDEnergyFitter : public TNamed
{
public: 
  /** 
   * Destructor
   */
  virtual ~AliFMDEnergyFitter();
  /** 
   * Default Constructor - do not use 
   */
  AliFMDEnergyFitter();
  /** 
   * Constructor 
   * 
   * @param title Title of object  - not significant 
   */
  AliFMDEnergyFitter(const char* title);
  /** 
   * Copy constructor 
   * 
   * @param o Object to copy from 
   */
  AliFMDEnergyFitter(const AliFMDEnergyFitter& o);
  /** 
   * Assignment operator 
   * 
   * @param o Object to assign from 
   * 
   * @return Reference to this 
   */
  AliFMDEnergyFitter& operator=(const AliFMDEnergyFitter& o);

  /** 
   * Initialise the task
   * 
   * @param etaAxis The eta axis to use.  Note, that if the eta axis
   * has already been set (using SetEtaAxis), then this parameter will be 
   * ignored
   */
  void Init(const TAxis& etaAxis);
  /** 
   * Set the eta axis to use.  This will force the code to use this
   * eta axis definition - irrespective of whatever axis is passed to
   * the Init member function.  Therefore, this member function can be
   * used to force another eta axis than one found in the correction
   * objects. 
   * 
   * @param nBins  Number of bins 
   * @param etaMin Minimum of the eta axis 
   * @param etaMax Maximum of the eta axis 
   */
  void SetEtaAxis(Int_t nBins, Double_t etaMin, Double_t etaMax);
  /** 
   * Set the eta axis to use.  This will force the code to use this
   * eta axis definition - irrespective of whatever axis is passed to
   * the Init member function.  Therefore, this member function can be
   * used to force another eta axis than one found in the correction
   * objects. 
   * 
   * @param etaAxis Eta axis to use 
   */
  void SetEtaAxis(const TAxis& etaAxis);
  /** 
   * Set the low cut used for energy 
   * 
   * @param lowCut Low cut
   */
  void SetLowCut(Double_t lowCut=0.3) { fLowCut = lowCut; }
  /** 
   * Set the number of bins to subtract 
   * 
   * @param n 
   */
  void SetBinsToSubtract(UShort_t n=4) { fBinsToSubtract = n; }
  /** 
   * Whether or not to enable fitting of the final merged result.  
   * Note, fitting takes quite a while and one should be careful not to do 
   * this needlessly 
   * 
   * @param doFit Whether to do the fits or not 
   */
  void DoFits(Bool_t doFit=kTRUE) { fDoFits = doFit; }
  /** 
   * Whether or not to enable fitting of the final merged result.  
   * Note, fitting takes quite a while and one should be careful not to do 
   * this needlessly 
   * 
   * @param doFit Whether to do the fits or not 
   */
  void SetNLandau(UShort_t n) { fNLandau = (n < 1 ? 1 : (n > 5 ? 5 : n)); }
  /** 
   * Whether or not to enable fitting of the final merged result.  
   * Note, fitting takes quite a while and one should be careful not to do 
   * this needlessly 
   * 
   * @param doFit Whether to do the fits or not 
   */
  void SetMinEntries(UShort_t n) { fMinEntries = (n < 1 ? 1 : n); }
  /**
   * Set maximum energy loss to consider 
   *
   * @param x Maximum energy loss to consider 
   */
  void SetMaxE(Double_t x) { fMaxE = x; }
  /**
   * Set number of energy loss bins 
   *
   * @param x Number of energy loss bins 
   */
  void SetNEbins(Int_t x) { fNEbins = x; }
  /**
   * Set wheter to use increasing bin sizes 
   *
   * @param x Wheter to use increasing bin sizes 
   */
  void SetUseIncreasingBins(Bool_t x) { fUseIncreasingBins = x; }
  /** 
   * Fitter the input AliESDFMD object
   * 
   * @param input     Input 
   * @param empty     Whether the event is 'empty'
   * 
   * @return True on success, false otherwise 
   */
  Bool_t Accumulate(const AliESDFMD& input, 
		    Bool_t           empty);
  /** 
   * Scale the histograms to the total number of events 
   * 
   * @param nEvents Number of events 
   */
  void Fit(TList* dir);
  
  /** 
   * Define the output histograms.  These are put in a sub list of the
   * passed list.   The histograms are merged before the parent task calls 
   * AliAnalysisTaskSE::Terminate 
   * 
   * @param dir Directory to add to 
   */
  void DefineOutput(TList* dir);
  /** 
   * Set the debug level.  The higher the value the more output 
   * 
   * @param dbg Debug level 
   */
  void SetDebug(Int_t dbg=1);
protected:
  /** 
   * Internal data structure to keep track of the histograms
   */
  struct RingHistos : public AliForwardUtil::RingHistos
  { 
    /** 
     * Default CTOR
     */
    RingHistos();
    /** 
     * Constructor
     * 
     * @param d detector
     * @param r ring 
     */
    RingHistos(UShort_t d, Char_t r);
    /** 
     * Copy constructor 
     * 
     * @param o Object to copy from 
     */
    RingHistos(const RingHistos& o);
    /** 
     * Assignment operator 
     * 
     * @param o Object to assign from 
     * 
     * @return Reference to this 
     */
    RingHistos& operator=(const RingHistos& o);
    /** 
     * Destructor 
     */
    ~RingHistos();
    /** 
     * Define outputs
     * 
     * @param dir 
     */
    void Output(TList* dir);
    /** 
     * Initialise object 
     * 
     * @param eAxis 
     */
    void Init(const TAxis& eAxis, 
	      Double_t     maxDE=10, 
	      Int_t        nDEbins=300, 
	      Bool_t       useIncrBin=true);
    /** 
     * Fill histogram 
     * 
     * @param empty  True if event is empty
     * @param ieta   Eta bin
     * @param mult   Signal 
     */
    void Fill(Bool_t empty, Int_t ieta, Double_t mult);
    /** 
     * Scale the histograms to the total number of events 
     * 
     * @param dir     Output list 
     * @param eta     Eta axis 
     * @param lowCut  Lower cut 
     * @param nLandau Max number of convolved landaus to fit
     */
    TObjArray* Fit(TList* dir, 
		   const TAxis& eta,
		   Double_t     lowCut, 
		   UShort_t     nLandau,
		   UShort_t     minEntries,
		   UShort_t     minusBins) const;
    /** 
     * Fit a signal histogram 
     * 
     * @param dist     Historgam to fit 
     * @param lowCut   Lower cut on signal 
     * @param nLandau  Max number of convolved landaus to fit
     * 
     * @return The best fit function 
     */
    TF1* FitHist(TH1*     dist,
		 Double_t lowCut, 
		 UShort_t nLandau,
		 UShort_t minusBins) const;
    /** 
     * Check the result of the fit. Returns true if 
     * - the reduced  @f$ \chi^2/\nu@f$ is less than 5
     * - and that the relative error @f$ \Delta p_i/p_i@f$ on each
     *   parameter is less than 20 percent.
     * - If this is a fit to N particles if the N particle weight is 
     *   larger than @$f 10^{-7}@f$ 
     * 
     * @param r Result to check
     * 
     * @return true if fit is good. 
     */
    Bool_t CheckResult(TFitResult* r) const;
    /** 
     * Make an axis with increasing bins 
     * 
     * @param n    Number of bins 
     * @param min  Minimum 
     * @param max  Maximum
     * 
     * @return An axis with quadratically increasing bin size 
     */
    TArrayD MakeIncreasingAxis(Int_t n, Double_t min, Double_t max) const;
    /** 
     * Make E/E_mip histogram 
     * 
     * @param ieta   Eta bin
     * @param eMin   Least signal
     * @param eMax   Largest signal 
     */
    void Make(Int_t ieta, Double_t eMin, Double_t eMax, 
	      Double_t deMax=12, Int_t nDeBins=300, Bool_t incr=true);
    /** 
     * Make a parameter histogram
     * 
     * @param name   Name of histogram.
     * @param title  Title of histogram. 
     * @param eta    Eta axis 
     * 
     * @return 
     */
    TH1D* MakePar(const char* name, const char* title, const TAxis& eta) const;
    /** 
     * Make a histogram that contains the results of the fit over the full ring 
     * 
     * @param name  Name 
     * @param title Title
     * @param eta   Eta axis 
     * @param low   Least bin
     * @param high  Largest bin
     * @param val   Value of parameter 
     * @param err   Error on parameter 
     * 
     * @return The newly allocated histogram 
     */
    TH1D* MakeTotal(const char* name, 
		    const char* title, 
		    const TAxis& eta, 
		    Int_t low, 
		    Int_t high, 
		    Double_t val, 
		    Double_t err) const;
    TH1D*     fEDist;        // Ring energy distribution 
    TH1D*     fEmpty;        // Ring energy distribution for empty events
    TList     fEtaEDists;    // Energy distributions per eta bin. 
    TList*    fList;
    Int_t     fDebug;
    ClassDef(RingHistos,1);
  };
  /** 
   * Get the ring histogram container 
   * 
   * @param d Detector
   * @param r Ring 
   * 
   * @return Ring histogram container 
   */
  RingHistos* GetRingHistos(UShort_t d, Char_t r) const;

  TList    fRingHistos;    // List of histogram containers
  Double_t fLowCut;        // Low cut on energy
  UShort_t fNLandau;       // Number of landaus to try to fit 
  UShort_t fMinEntries;    // Minimum number of entries
  UShort_t fBinsToSubtract;// Number of bins to subtract from found max
  Bool_t   fDoFits;        // Wheter to actually do the fits 
  TAxis    fEtaAxis;       // Eta axis 
  Double_t fMaxE;          // Maximum energy loss to consider 
  Int_t    fNEbins;        // Number of energy loss bins 
  Bool_t   fUseIncreasingBins; // Wheter to use increasing bin sizes 
  Int_t    fDebug;         // Debug level 

  ClassDef(AliFMDEnergyFitter,1); //
};

#endif
// Local Variables:
//  mode: C++ 
// End:
Commit	Line	Data
f8715167	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 "AliForwardUtil.h"
	9	class AliESDFMD;
	10	class TFitResult;
	11	class TF1;
	12	class TArrayD;
	13
	14	/**
	15	* Class to fit the energy distribution.
	16	*
	17	* @par Input:
	18	* - AliESDFMD object - from reconstruction
	19	*
	20	* @par Output:
	21	* - Lists of histogram - one per ring. Each list has a number of
	22	* histograms corresponding to the number of eta bins defined.
	23	*
	24	* @par Corrections used:
	25	* - None
	26	*
	27	*
	28	* @ingroup pwg2_forward_analysis
	29	*/
	30	class AliFMDEnergyFitter : public TNamed
	31	{
	32	public:
	33	/**
	34	* Destructor
	35	*/
	36	virtual ~AliFMDEnergyFitter();
	37	/**
	38	* Default Constructor - do not use
	39	*/
	40	AliFMDEnergyFitter();
	41	/**
	42	* Constructor
	43	*
	44	* @param title Title of object - not significant
	45	*/
	46	AliFMDEnergyFitter(const char* title);
	47	/**
	48	* Copy constructor
	49	*
	50	* @param o Object to copy from
	51	*/
	52	AliFMDEnergyFitter(const AliFMDEnergyFitter& o);
	53	/**
	54	* Assignment operator
	55	*
	56	* @param o Object to assign from
	57	*
	58	* @return Reference to this
	59	*/
	60	AliFMDEnergyFitter& operator=(const AliFMDEnergyFitter& o);
	61
4b9857f3	62	/**
	63	* Initialise the task
	64	*
	65	* @param etaAxis The eta axis to use. Note, that if the eta axis
	66	* has already been set (using SetEtaAxis), then this parameter will be
	67	* ignored
	68	*/
f8715167	69	void Init(const TAxis& etaAxis);
4b9857f3	70	/**
	71	* Set the eta axis to use. This will force the code to use this
	72	* eta axis definition - irrespective of whatever axis is passed to
	73	* the Init member function. Therefore, this member function can be
	74	* used to force another eta axis than one found in the correction
	75	* objects.
	76	*
	77	* @param nBins Number of bins
	78	* @param etaMin Minimum of the eta axis
	79	* @param etaMax Maximum of the eta axis
	80	*/
	81	void SetEtaAxis(Int_t nBins, Double_t etaMin, Double_t etaMax);
	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 etaAxis Eta axis to use
	90	*/
	91	void SetEtaAxis(const TAxis& etaAxis);
f8715167	92	/**
	93	* Set the low cut used for energy
	94	*
	95	* @param lowCut Low cut
	96	*/
	97	void SetLowCut(Double_t lowCut=0.3) { fLowCut = lowCut; }
4b9857f3	98	/**
	99	* Set the number of bins to subtract
	100	*
	101	* @param n
	102	*/
f8715167	103	void SetBinsToSubtract(UShort_t n=4) { fBinsToSubtract = n; }
	104	/**
	105	* Whether or not to enable fitting of the final merged result.
	106	* Note, fitting takes quite a while and one should be careful not to do
	107	* this needlessly
	108	*
	109	* @param doFit Whether to do the fits or not
	110	*/
	111	void DoFits(Bool_t doFit=kTRUE) { fDoFits = doFit; }
	112	/**
	113	* Whether or not to enable fitting of the final merged result.
	114	* Note, fitting takes quite a while and one should be careful not to do
	115	* this needlessly
	116	*
	117	* @param doFit Whether to do the fits or not
	118	*/
	119	void SetNLandau(UShort_t n) { fNLandau = (n < 1 ? 1 : (n > 5 ? 5 : n)); }
	120	/**
	121	* Whether or not to enable fitting of the final merged result.
	122	* Note, fitting takes quite a while and one should be careful not to do
	123	* this needlessly
	124	*
	125	* @param doFit Whether to do the fits or not
	126	*/
	127	void SetMinEntries(UShort_t n) { fMinEntries = (n < 1 ? 1 : n); }
4b9857f3	128	/**
	129	* Set maximum energy loss to consider
	130	*
	131	* @param x Maximum energy loss to consider
	132	*/
	133	void SetMaxE(Double_t x) { fMaxE = x; }
	134	/**
	135	* Set number of energy loss bins
	136	*
	137	* @param x Number of energy loss bins
	138	*/
	139	void SetNEbins(Int_t x) { fNEbins = x; }
	140	/**
	141	* Set wheter to use increasing bin sizes
	142	*
	143	* @param x Wheter to use increasing bin sizes
	144	*/
	145	void SetUseIncreasingBins(Bool_t x) { fUseIncreasingBins = x; }
f8715167	146	/**
	147	* Fitter the input AliESDFMD object
	148	*
	149	* @param input Input
	150	* @param empty Whether the event is 'empty'
	151	*
	152	* @return True on success, false otherwise
	153	*/
	154	Bool_t Accumulate(const AliESDFMD& input,
	155	Bool_t empty);
	156	/**
	157	* Scale the histograms to the total number of events
	158	*
	159	* @param nEvents Number of events
	160	*/
	161	void Fit(TList* dir);
	162
	163	/**
	164	* Define the output histograms. These are put in a sub list of the
	165	* passed list. The histograms are merged before the parent task calls
	166	* AliAnalysisTaskSE::Terminate
	167	*
	168	* @param dir Directory to add to
	169	*/
	170	void DefineOutput(TList* dir);
	171	/**
	172	* Set the debug level. The higher the value the more output
	173	*
	174	* @param dbg Debug level
	175	*/
	176	void SetDebug(Int_t dbg=1);
	177	protected:
	178	/**
	179	* Internal data structure to keep track of the histograms
	180	*/
	181	struct RingHistos : public AliForwardUtil::RingHistos
	182	{
	183	/**
	184	* Default CTOR
	185	*/
	186	RingHistos();
	187	/**
	188	* Constructor
	189	*
	190	* @param d detector
	191	* @param r ring
	192	*/
	193	RingHistos(UShort_t d, Char_t r);
	194	/**
	195	* Copy constructor
	196	*
	197	* @param o Object to copy from
	198	*/
	199	RingHistos(const RingHistos& o);
	200	/**
	201	* Assignment operator
	202	*
	203	* @param o Object to assign from
	204	*
	205	* @return Reference to this
	206	*/
	207	RingHistos& operator=(const RingHistos& o);
	208	/**
	209	* Destructor
210	*/
211	~RingHistos();
212	/**
213	* Define outputs
214	*
215	* @param dir
216	*/
217	void Output(TList* dir);
218	/**
219	* Initialise object
220	*
221	* @param eAxis
222	*/
66b34429	223	void Init(const TAxis& eAxis,
4b9857f3	224	Double_t maxDE=10,
	225	Int_t nDEbins=300,
	226	Bool_t useIncrBin=true);
f8715167	227	/**
	228	* Fill histogram
	229	*
	230	* @param empty True if event is empty
	231	* @param ieta Eta bin
	232	* @param mult Signal
	233	*/
	234	void Fill(Bool_t empty, Int_t ieta, Double_t mult);
	235	/**
	236	* Scale the histograms to the total number of events
	237	*
	238	* @param dir Output list
	239	* @param eta Eta axis
	240	* @param lowCut Lower cut
	241	* @param nLandau Max number of convolved landaus to fit
	242	*/
4b9857f3	243	TObjArray* Fit(TList* dir,
	244	const TAxis& eta,
	245	Double_t lowCut,
	246	UShort_t nLandau,
	247	UShort_t minEntries,
	248	UShort_t minusBins) const;
f8715167	249	/**
	250	* Fit a signal histogram
	251	*
	252	* @param dist Historgam to fit
	253	* @param lowCut Lower cut on signal
	254	* @param nLandau Max number of convolved landaus to fit
	255	*
	256	* @return The best fit function
	257	*/
4b9857f3	258	TF1* FitHist(TH1* dist,
	259	Double_t lowCut,
	260	UShort_t nLandau,
	261	UShort_t minusBins) const;
f8715167	262	/**
4b9857f3	263	* Check the result of the fit. Returns true if
	264	* - the reduced @f$ \chi^2/\nu@f$ is less than 5
	265	* - and that the relative error @f$ \Delta p_i/p_i@f$ on each
	266	* parameter is less than 20 percent.
	267	* - If this is a fit to N particles if the N particle weight is
	268	* larger than @$f 10^{-7}@f$
f8715167	269	*
4b9857f3	270	* @param r Result to check
f8715167	271	*
4b9857f3	272	* @return true if fit is good.
f8715167	273	*/
4b9857f3	274	Bool_t CheckResult(TFitResult* r) const;
f8715167	275	/**
4b9857f3	276	* Make an axis with increasing bins
f8715167	277	*
4b9857f3	278	* @param n Number of bins
	279	* @param min Minimum
	280	* @param max Maximum
f8715167	281	*
4b9857f3	282	* @return An axis with quadratically increasing bin size
f8715167	283	*/
66b34429	284	TArrayD MakeIncreasingAxis(Int_t n, Double_t min, Double_t max) const;
f8715167	285	/**
	286	* Make E/E_mip histogram
	287	*
	288	* @param ieta Eta bin
	289	* @param eMin Least signal
	290	* @param eMax Largest signal
	291	*/
66b34429	292	void Make(Int_t ieta, Double_t eMin, Double_t eMax,
66b34429	293	Double_t deMax=12, Int_t nDeBins=300, Bool_t incr=true);
f8715167	294	/**
	295	* Make a parameter histogram
	296	*
	297	* @param name Name of histogram.
	298	* @param title Title of histogram.
	299	* @param eta Eta axis
	300	*
	301	* @return
	302	*/
	303	TH1D* MakePar(const char* name, const char* title, const TAxis& eta) const;
4b9857f3	304	/**
	305	* Make a histogram that contains the results of the fit over the full ring
	306	*
	307	* @param name Name
	308	* @param title Title
	309	* @param eta Eta axis
	310	* @param low Least bin
	311	* @param high Largest bin
	312	* @param val Value of parameter
	313	* @param err Error on parameter
	314	*
	315	* @return The newly allocated histogram
	316	*/
f8715167	317	TH1D* MakeTotal(const char* name,
	318	const char* title,
	319	const TAxis& eta,
	320	Int_t low,
	321	Int_t high,
	322	Double_t val,
	323	Double_t err) const;
	324	TH1D* fEDist; // Ring energy distribution
	325	TH1D* fEmpty; // Ring energy distribution for empty events
	326	TList fEtaEDists; // Energy distributions per eta bin.
	327	TList* fList;
	328	Int_t fDebug;
	329	ClassDef(RingHistos,1);
	330	};
	331	/**
	332	* Get the ring histogram container
	333	*
	334	* @param d Detector
	335	* @param r Ring
	336	*
	337	* @return Ring histogram container
	338	*/
	339	RingHistos* GetRingHistos(UShort_t d, Char_t r) const;
	340
	341	TList fRingHistos; // List of histogram containers
	342	Double_t fLowCut; // Low cut on energy
	343	UShort_t fNLandau; // Number of landaus to try to fit
	344	UShort_t fMinEntries; // Minimum number of entries
	345	UShort_t fBinsToSubtract;// Number of bins to subtract from found max
	346	Bool_t fDoFits; // Wheter to actually do the fits
	347	TAxis fEtaAxis; // Eta axis
4b9857f3	348	Double_t fMaxE; // Maximum energy loss to consider
	349	Int_t fNEbins; // Number of energy loss bins
	350	Bool_t fUseIncreasingBins; // Wheter to use increasing bin sizes
f8715167	351	Int_t fDebug; // Debug level
	352
	353	ClassDef(AliFMDEnergyFitter,1); //
	354	};
	355
	356	#endif
	357	// Local Variables:
	358	// mode: C++
	359	// End: