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      * @return List of fits
 333      */
 334     TObjArray* Fit(TList* dir,
 335                    const TAxis& eta,
 336                    Double_t     lowCut,
 337                    UShort_t     nParticles,
 338                    UShort_t     minEntries,
 339                    UShort_t     minusBins,
 340                    Double_t     relErrorCut,
 341                    Double_t     chi2nuCut) const;
 342     /**
 343      * Fit a signal histogram.  First, the bin @f$ b_{min}@f$ with
 344      * maximum bin content in the range @f$ [E_{min},\infty]@f$ is
 345      * found.  Then the fit range is set to the bin range
 346      * @f$ [b_{min}-\Delta b,b_{min}+2\Delta b]@f$, and a 1
 347      * particle signal is fitted to that.  The parameters of that fit
 348      * is then used as seeds for a fit of the @f$ N@f$ particle response
 349      * to the data in the range
 350      * @f$ [b_{min}-\Delta b,N(\Delta_1+\xi_1\log(N))+2N\xi@f$
 351      *
 352      * @param dist        Histogram to fit
 353      * @param lowCut      Lower cut @f$ E_{min}@f$ on signal
 354      * @param nParticles  Max number @f$ N@f$ of convolved landaus to fit
 355      * @param minusBins   Number of bins @f$ \Delta b@f$ from peak to
 356      *                    subtract to get the fit range
 357      * @param relErrorCut Cut applied to relative error of parameter.
 358      *                    Note, for multi-particle weights, the cut
 359      *                    is loosend by a factor of 2
 360      * @param chi2nuCut   Cut on @f$ \chi^2/\nu@f$ -
 361      *                    the reduced @f$\chi^2@f$
 362      *
 363      * @return The best fit function
 364      */
 365     TF1* FitHist(TH1*     dist,
 366                  Double_t lowCut,
 367                  UShort_t nParticles,
 368                  UShort_t minusBins,
 369                  Double_t relErrorCut,
 370                  Double_t chi2nuCut) const;
 371     /**
 372      * Find the best fits
 373      *
 374      * @param d              Parent list
 375      * @param obj            Object to add fits to
 376      * @param eta            Eta axis
 377      * @param relErrorCut    Cut applied to relative error of parameter.
 378      *                       Note, for multi-particle weights, the cut
 379      *                       is loosend by a factor of 2
 380      * @param chi2nuCut      Cut on @f$ \chi^2/\nu@f$ -
 381      *                       the reduced @f$\chi^2@f$
 382      * @param minWeightCut   Least valid @f$ a_i@f$
 383      */
 384     void FindBestFits(const TList*        d,
 385                       AliFMDCorrELossFit& obj,
 386                       const TAxis&        eta,
 387                       Double_t            relErrorCut,
 388                       Double_t            chi2nuCut,
 389                       Double_t            minWeightCut);
 390     /**
 391      * Find the best fit
 392      *
 393      * @param dist           Histogram
 394      * @param relErrorCut    Cut applied to relative error of parameter.
 395      *                       Note, for multi-particle weights, the cut
 396      *                       is loosend by a factor of 2
 397      * @param chi2nuCut      Cut on @f$ \chi^2/\nu@f$ -
 398      *                       the reduced @f$\chi^2@f$
 399      * @param minWeightCut   Least valid @f$ a_i@f$
 400      *
 401      * @return Best fit
 402      */
 403     AliFMDCorrELossFit::ELossFit* FindBestFit(const TH1* dist,
 404                                               Double_t relErrorCut,
 405                                               Double_t chi2nuCut,
 406                                               Double_t minWeightCut);
 407     /**
 408      * Check the result of the fit. Returns true if
 409      * - @f$ \chi^2/\nu < \max{\chi^2/\nu}@f$
 410      * - @f$ \Delta p_i/p_i < \delta_e@f$ for all parameters.  Note,
 411      *   for multi-particle fits, this requirement is relaxed by a
 412      *   factor of 2
 413      * - @f$ a_{n} > 10^{-7}@f$ when fitting to an @f$ n@f$
 414      *   particle response
 415      *
 416      * @param r           Result to check
 417      * @param relErrorCut Cut @f$ \delta_e@f$ applied to relative error
 418      *                    of parameter.
 419      * @param chi2nuCut   Cut @f$ \max{\chi^2/\nu}@f$
 420      *
 421      * @return true if fit is good.
 422      */
 423     Bool_t CheckResult(TFitResult* r,
 424                        Double_t    relErrorCut,
 425                        Double_t    chi2nuCut) const;
 426     /**
 427      * Make an axis with increasing bins
 428      *
 429      * @param n    Number of bins
 430      * @param min  Minimum
 431      * @param max  Maximum
 432      *
 433      * @return An axis with quadratically increasing bin size
 434      */
 435     TArrayD MakeIncreasingAxis(Int_t n, Double_t min, Double_t max) const;
 436     /**
 437      * Make E/E_mip histogram
 438      *
 439      * @param ieta    Eta bin
 440      * @param eMin    Least signal
 441      * @param eMax    Largest signal
 442      * @param deMax   Maximum energy loss
 443      * @param nDeBins Number energy loss bins
 444      * @param incr    Whether to make bins of increasing size
 445      */
 446     void Make(Int_t ieta, Double_t eMin, Double_t eMax,
 447               Double_t deMax=12, Int_t nDeBins=300, Bool_t incr=true);
 448     /**
 449      * Make a parameter histogram
 450      *
 451      * @param name   Name of histogram.
 452      * @param title  Title of histogram.
 453      * @param eta    Eta axis
 454      *
 455      * @return
 456      */
 457     TH1D* MakePar(const char* name, const char* title, const TAxis& eta) const;
 458     /**
 459      * Make a histogram that contains the results of the fit over the full ring
 460      *
 461      * @param name  Name
 462      * @param title Title
 463      * @param eta   Eta axis
 464      * @param low   Least bin
 465      * @param high  Largest bin
 466      * @param val   Value of parameter
 467      * @param err   Error on parameter
 468      *
 469      * @return The newly allocated histogram
 470      */
 471     TH1D* MakeTotal(const char* name,
 472                     const char* title,
 473                     const TAxis& eta,
 474                     Int_t low,
 475                     Int_t high,
 476                     Double_t val,
 477                     Double_t err) const;
 478     TH1D*        fEDist;        // Ring energy distribution
 479     TH1D*        fEmpty;        // Ring energy distribution for empty events
 480     TList        fEtaEDists;    // Energy distributions per eta bin.
 481     TList*       fList;
 482     TClonesArray fFits;
 483     Int_t        fDebug;
 484     ClassDef(RingHistos,1);
 485   };
 486   /**
 487    * Get the ring histogram container
 488    *
 489    * @param d Detector
 490    * @param r Ring
 491    *
 492    * @return Ring histogram container
 493    */
 494   RingHistos* GetRingHistos(UShort_t d, Char_t r) const;
 495
 496   TList    fRingHistos;    // List of histogram containers
 497   Double_t fLowCut;        // Low cut on energy
 498   UShort_t fNParticles;    // Number of landaus to try to fit
 499   UShort_t fMinEntries;    // Minimum number of entries
 500   UShort_t fFitRangeBinWidth;// Number of bins to subtract from found max
 501   Bool_t   fDoFits;        // Whether to actually do the fits
 502   Bool_t   fDoMakeObject;  // Whether to make corrections object
 503   TAxis    fEtaAxis;       // Eta axis
 504   TAxis    fCentralityAxis;// Centrality axis
 505   Double_t fMaxE;          // Maximum energy loss to consider
 506   Int_t    fNEbins;        // Number of energy loss bins
 507   Bool_t   fUseIncreasingBins; // Wheter to use increasing bin sizes
 508   Double_t fMaxRelParError;// Relative error cut
 509   Double_t fMaxChi2PerNDF; // chi^2/nu cit
 510   Double_t fMinWeight;     // Minimum weight value
 511   Int_t    fDebug;         // Debug level
 512
 513
 514   ClassDef(AliFMDEnergyFitter,2); //
 515 };
 516
 517 #endif
 518 // Local Variables:
 519 //  mode: C++
 520 // End: