PWGLF/FORWARD/analysis2/AliFMDEnergyFitter.h

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