PWG2/FORWARD/analysis2/AliFMDEnergyFitter.h

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