PWG2/FORWARD/analysis2/AliFMDEnergyFitter.h

   1 #ifndef ALIFMDENERGYFITTER_H
   2 #define 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_algo
  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   /**
 157    * Set the maximum relative error
 158    *
 159    * @param e Maximum relative error
 160    */
 161   void SetMaxRelativeParameterError(Double_t e=0.2) { fMaxRelParError = e; }
 162   /**
 163    * Set the maximum @f$ \chi^2/\nu@f$
 164    *
 165    * @param c Maximum @f$ \chi^2/\nu@f$
 166    */
 167   void SetMaxChi2PerNDF(Double_t c=10) { fMaxChi2PerNDF = c; }
 168   /**
 169    * Set the least weight
 170    *
 171    * @param c Least weight
 172    */
 173   void SetMinWeight(Double_t c=1e-7) { fMinWeight = c; }
 174   /**
 175    * Set wheter to use increasing bin sizes
 176    *
 177    * @param x Wheter to use increasing bin sizes
 178    */
 179   void SetUseIncreasingBins(Bool_t x) { fUseIncreasingBins = x; }
 180   /**
 181    * Fitter the input AliESDFMD object
 182    *
 183    * @param input     Input
 184    * @param empty     Whether the event is 'empty'
 185    *
 186    * @return True on success, false otherwise
 187    */
 188   Bool_t Accumulate(const AliESDFMD& input,
 189                     Bool_t           empty);
 190   /**
 191    * Scale the histograms to the total number of events
 192    *
 193    * @param dir Where the histograms are
 194    */
 195   void Fit(const TList* dir);
 196   /**
 197    * Generate the corrections object
 198    *
 199    * @param dir List to analyse
 200    */
 201   void MakeCorrectionsObject(TList* dir);
 202
 203   /**
 204    * Define the output histograms.  These are put in a sub list of the
 205    * passed list.   The histograms are merged before the parent task calls
 206    * AliAnalysisTaskSE::Terminate
 207    *
 208    * @param dir Directory to add to
 209    */
 210   void DefineOutput(TList* dir);
 211   /**
 212    * Set the debug level.  The higher the value the more output
 213    *
 214    * @param dbg Debug level
 215    */
 216   void SetDebug(Int_t dbg=1);
 217   /**
 218    * Print information
 219    *
 220    * @param option Not used
 221    */
 222   void Print(Option_t* option="") const;
 223 protected:
 224   /**
 225    * Internal data structure to keep track of the histograms
 226    */
 227   struct RingHistos : public AliForwardUtil::RingHistos
 228   {
 229     /**
 230      * Default CTOR
 231      */
 232     RingHistos();
 233     /**
 234      * Constructor
 235      *
 236      * @param d detector
 237      * @param r ring
 238      */
 239     RingHistos(UShort_t d, Char_t r);
 240     /**
 241      * Copy constructor
 242      *
 243      * @param o Object to copy from
 244      */
 245     RingHistos(const RingHistos& o);
 246     /**
 247      * Assignment operator
 248      *
 249      * @param o Object to assign from
 250      *
 251      * @return Reference to this
 252      */
 253     RingHistos& operator=(const RingHistos& o);
 254     /**
 255      * Destructor
 256      */
 257     ~RingHistos();
 258     /**
 259      * Define outputs
 260      *
 261      * @param dir
 262      */
 263     void Output(TList* dir);
 264     /**
 265      * Initialise object
 266      *
 267      * @param eAxis      Eta axis
 268      * @param maxDE      Max energy loss to consider
 269      * @param nDEbins    Number of bins
 270      * @param useIncrBin Whether to use an increasing bin size
 271      */
 272     void Init(const TAxis& eAxis,
 273               Double_t     maxDE=10,
 274               Int_t        nDEbins=300,
 275               Bool_t       useIncrBin=true);
 276     /**
 277      * Fill histogram
 278      *
 279      * @param empty  True if event is empty
 280      * @param ieta   Eta bin
 281      * @param mult   Signal
 282      */
 283     void Fill(Bool_t empty, Int_t ieta, Double_t mult);
 284     /**
 285      * Fit each histogram to up to @a nParticles particle responses.
 286      *
 287      * @param dir         Output list
 288      * @param eta         Eta axis
 289      * @param lowCut      Lower cut
 290      * @param nParticles  Max number of convolved landaus to fit
 291      * @param minEntries  Minimum number of entries
 292      * @param minusBins   Number of bins from peak to subtract to
 293      *                    get the fit range
 294      * @param relErrorCut Cut applied to relative error of parameter.
 295      *                    Note, for multi-particle weights, the cut
 296      *                    is loosend by a factor of 2
 297      * @param chi2nuCut   Cut on @f$ \chi^2/\nu@f$ -
 298      *                    the reduced @f$\chi^2@f$
 299      */
 300     TObjArray* Fit(TList* dir,
 301                    const TAxis& eta,
 302                    Double_t     lowCut,
 303                    UShort_t     nParticles,
 304                    UShort_t     minEntries,
 305                    UShort_t     minusBins,
 306                    Double_t     relErrorCut,
 307                    Double_t     chi2nuCut) const;
 308     /**
 309      * Fit a signal histogram.  First, the bin @f$ b_{min}@f$ with
 310      * maximum bin content in the range @f$ [E_{min},\infty]@f$ is
 311      * found.  Then the fit range is set to the bin range
 312      * @f$ [b_{min}-\Delta b,b_{min}+2\Delta b]@f$, and a 1
 313      * particle signal is fitted to that.  The parameters of that fit
 314      * is then used as seeds for a fit of the @f$ N@f$ particle response
 315      * to the data in the range
 316      * @f$ [b_{min}-\Delta b,N(\Delta_1+\xi_1\log(N))+2N\xi@f$
 317      *
 318      * @param dist        Histogram to fit
 319      * @param lowCut      Lower cut @f$ E_{min}@f$ on signal
 320      * @param nParticles  Max number @f$ N@f$ of convolved landaus to fit
 321      * @param minusBins   Number of bins @f$ \Delta b@f$ from peak to
 322      *                    subtract to get the fit range
 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      *
 329      * @return The best fit function
 330      */
 331     TF1* FitHist(TH1*     dist,
 332                  Double_t lowCut,
 333                  UShort_t nParticles,
 334                  UShort_t minusBins,
 335                  Double_t relErrorCut,
 336                  Double_t chi2nuCut) const;
 337     /**
 338      * Find the best fits
 339      *
 340      * @param d              Parent list
 341      * @param obj            Object to add fits to
 342      * @param eta            Eta axis
 343      * @param relErrorCut    Cut applied to relative error of parameter.
 344      *                       Note, for multi-particle weights, the cut
 345      *                       is loosend by a factor of 2
 346      * @param chi2nuCut      Cut on @f$ \chi^2/\nu@f$ -
 347      *                       the reduced @f$\chi^2@f$
 348      * @param minWeightCut   Least valid @f$ a_i@f$
 349      */
 350     void FindBestFits(TList*              d,
 351                       AliFMDCorrELossFit& obj,
 352                       const TAxis&        eta,
 353                       Double_t            relErrorCut,
 354                       Double_t            chi2nuCut,
 355                       Double_t            minWeightCut);
 356     /**
 357      * Find the best fit
 358      *
 359      * @param dist           Histogram
 360      * @param relErrorCut    Cut applied to relative error of parameter.
 361      *                       Note, for multi-particle weights, the cut
 362      *                       is loosend by a factor of 2
 363      * @param chi2nuCut      Cut on @f$ \chi^2/\nu@f$ -
 364      *                       the reduced @f$\chi^2@f$
 365      * @param minWeightCut   Least valid @f$ a_i@f$
 366      *
 367      * @return Best fit
 368      */
 369     AliFMDCorrELossFit::ELossFit* FindBestFit(TH1* dist,
 370                                               Double_t relErrorCut,
 371                                               Double_t chi2nuCut,
 372                                               Double_t minWeightCut);
 373     /**
 374      * Check the result of the fit. Returns true if
 375      * - @f$ \chi^2/\nu < \max{\chi^2/\nu}@f$
 376      * - @f$ \Delta p_i/p_i < \delta_e@f$ for all parameters.  Note,
 377      *   for multi-particle fits, this requirement is relaxed by a
 378      *   factor of 2
 379      * - @f$ a_{n} > 10^{-7}@f$ when fitting to an @f$ n@f$
 380      *   particle response
 381      *
 382      * @param r           Result to check
 383      * @param relErrorCut Cut @f$ \delta_e@f$ applied to relative error
 384      *                    of parameter.
 385      * @param chi2nuCut   Cut @f$ \max{\chi^2/\nu}@f$
 386      *
 387      * @return true if fit is good.
 388      */
 389     Bool_t CheckResult(TFitResult* r,
 390                        Double_t    relErrorCut,
 391                        Double_t    chi2nuCut) const;
 392     /**
 393      * Make an axis with increasing bins
 394      *
 395      * @param n    Number of bins
 396      * @param min  Minimum
 397      * @param max  Maximum
 398      *
 399      * @return An axis with quadratically increasing bin size
 400      */
 401     TArrayD MakeIncreasingAxis(Int_t n, Double_t min, Double_t max) const;
 402     /**
 403      * Make E/E_mip histogram
 404      *
 405      * @param ieta    Eta bin
 406      * @param eMin    Least signal
 407      * @param eMax    Largest signal
 408      * @param deMax   Maximum energy loss
 409      * @param nDeBins Number energy loss bins
 410      * @param incr    Whether to make bins of increasing size
 411      */
 412     void Make(Int_t ieta, Double_t eMin, Double_t eMax,
 413               Double_t deMax=12, Int_t nDeBins=300, Bool_t incr=true);
 414     /**
 415      * Make a parameter histogram
 416      *
 417      * @param name   Name of histogram.
 418      * @param title  Title of histogram.
 419      * @param eta    Eta axis
 420      *
 421      * @return
 422      */
 423     TH1D* MakePar(const char* name, const char* title, const TAxis& eta) const;
 424     /**
 425      * Make a histogram that contains the results of the fit over the full ring
 426      *
 427      * @param name  Name
 428      * @param title Title
 429      * @param eta   Eta axis
 430      * @param low   Least bin
 431      * @param high  Largest bin
 432      * @param val   Value of parameter
 433      * @param err   Error on parameter
 434      *
 435      * @return The newly allocated histogram
 436      */
 437     TH1D* MakeTotal(const char* name,
 438                     const char* title,
 439                     const TAxis& eta,
 440                     Int_t low,
 441                     Int_t high,
 442                     Double_t val,
 443                     Double_t err) const;
 444     TH1D*        fEDist;        // Ring energy distribution
 445     TH1D*        fEmpty;        // Ring energy distribution for empty events
 446     TList        fEtaEDists;    // Energy distributions per eta bin.
 447     TList*       fList;
 448     TClonesArray fFits;
 449     Int_t        fDebug;
 450     ClassDef(RingHistos,1);
 451   };
 452   /**
 453    * Get the ring histogram container
 454    *
 455    * @param d Detector
 456    * @param r Ring
 457    *
 458    * @return Ring histogram container
 459    */
 460   RingHistos* GetRingHistos(UShort_t d, Char_t r) const;
 461
 462   TList    fRingHistos;    // List of histogram containers
 463   Double_t fLowCut;        // Low cut on energy
 464   UShort_t fNParticles;    // Number of landaus to try to fit
 465   UShort_t fMinEntries;    // Minimum number of entries
 466   UShort_t fFitRangeBinWidth;// Number of bins to subtract from found max
 467   Bool_t   fDoFits;        // Whether to actually do the fits
 468   Bool_t   fDoMakeObject;  // Whether to make corrections object
 469   TAxis    fEtaAxis;       // Eta axis
 470   Double_t fMaxE;          // Maximum energy loss to consider
 471   Int_t    fNEbins;        // Number of energy loss bins
 472   Bool_t   fUseIncreasingBins; // Wheter to use increasing bin sizes
 473   Double_t fMaxRelParError;// Relative error cut
 474   Double_t fMaxChi2PerNDF; // chi^2/nu cit
 475   Double_t fMinWeight;     // Minimum weight value
 476   Int_t    fDebug;         // Debug level
 477
 478
 479   ClassDef(AliFMDEnergyFitter,1); //
 480 };
 481
 482 #endif
 483 // Local Variables:
 484 //  mode: C++
 485 // End: