PWGLF/FORWARD/analysis2/AliFMDDensityCalculator.h

   1 // This class calculates the inclusive charged particle density
   2 // in each for the 5 FMD rings.
   3 //
   4 #ifndef ALIFMDDENSITYCALCULATOR_H
   5 #define ALIFMDDENSITYCALCULATOR_H
   6 /**
   7  * @file   AliFMDDensityCalculator.h
   8  * @author Christian Holm Christensen <cholm@dalsgaard.hehi.nbi.dk>
   9  * @date   Wed Mar 23 14:02:09 2011
  10  *
  11  * @brief
  12  *
  13  *
  14  * @ingroup pwglf_forward_aod
  15  */
  16 #include <TNamed.h>
  17 #include <TList.h>
  18 #include <TArrayI.h>
  19 #include "AliForwardUtil.h"
  20 #include "AliFMDMultCuts.h"
  21 #include "AliPoissonCalculator.h"
  22 class AliESDFMD;
  23 class TH2D;
  24 class TH1D;
  25 class TProfile;
  26 class AliFMDCorrELossFit;
  27
  28 /**
  29  * This class calculates the inclusive charged particle density
  30  * in each for the 5 FMD rings.
  31  *
  32  * @par Input:
  33  *   - AliESDFMD object possibly corrected for sharing
  34  *
  35  * @par Output:
  36  *   - 5 RingHistos objects - each with a number of vertex dependent
  37  *     2D histograms of the inclusive charge particle density
  38  *
  39  * @par Corrections used:
  40  *   - AliFMDAnaCalibEnergyDistribution
  41  *   - AliFMDDoubleHitCorrection
  42  *   - AliFMDDeadCorrection
  43  *
  44  * @ingroup pwglf_forward_algo
  45  * @ingroup pwglf_forward_aod
  46  */
  47 class AliFMDDensityCalculator : public TNamed
  48 {
  49 public:
  50   /**
  51    * How to correct for the missing phi coverage at the corners of the
  52    * sensors
  53    *
  54    */
  55   enum {
  56     /** No correction */
  57     kPhiNoCorrect,
  58     /** Correct the calculated number charged particles */
  59     kPhiCorrectNch,
  60     /** Correct the energy loss */
  61     kPhiCorrectELoss
  62   };
  63   /**
  64    * Constructor
  65    */
  66   AliFMDDensityCalculator();
  67   /**
  68    * Constructor
  69    *
  70    * @param name Name of object
  71    */
  72   AliFMDDensityCalculator(const char* name);
  73   /**
  74    * Copy constructor
  75    *
  76    * @param o Object to copy from
  77    */
  78   AliFMDDensityCalculator(const AliFMDDensityCalculator& o);
  79   /**
  80    * Destructor
  81    */
  82   virtual ~AliFMDDensityCalculator();
  83   /**
  84    * Assignement operator
  85    *
  86    * @param o Object to assign from
  87    *
  88    * @return Reference to this object
  89    */
  90   AliFMDDensityCalculator& operator=(const AliFMDDensityCalculator& o);
  91   /**
  92    * Initialize this sub-algorithm
  93    *
  94    * @param etaAxis Not used
  95    */
  96   virtual void Init(const TAxis& etaAxis);
  97   /**
  98    * Do the calculations
  99    *
 100    * @param fmd      AliESDFMD object (possibly) corrected for sharing
 101    * @param hists    Histogram cache
 102    * @param vtxBin   Vertex bin
 103    * @param lowFlux  Low flux flag.
 104    *
 105    * @return true on successs
 106    */
 107   virtual Bool_t Calculate(const AliESDFMD& fmd,
 108                            AliForwardUtil::Histos& hists,
 109                            UShort_t vtxBin, Bool_t lowFlux, Double_t cent=-1, Double_t vz=0);
 110   /**
 111    * Scale the histograms to the total number of events
 112    *
 113    * @param dir     where to put the output
 114    * @param nEvents Number of events
 115    */
 116   virtual void ScaleHistograms(const TList* dir, Int_t nEvents);
 117   /**
 118    * Output diagnostic histograms to directory
 119    *
 120    * @param dir List to write in
 121    */
 122   virtual void DefineOutput(TList* dir);
 123   /**
 124    * Set the debug level.  The higher the value the more output
 125    *
 126    * @param dbg Debug level
 127    */
 128   void SetDebug(Int_t dbg=1) { fDebug = dbg; }
 129   /**
 130    * Maximum particle weight to use
 131    *
 132    * @param m
 133    */
 134   void SetMaxParticles(UShort_t m) { fMaxParticles = m; }
 135   /**
 136    * Set whether to use poisson statistics to estimate the
 137    * number of particles that has hit within a region.  If this is true,
 138    * then the average charge particle density is given by
 139    * @f[
 140    *  \lambda = -\log\left(\frac{N_e}{N_t}\right)
 141    * @f]
 142    * where $N_e$ is the number of strips within the region that has no
 143    * hits over threshold, and $N_t$ is the total number of strips in the
 144    * region/
 145    *
 146    * @param u Whether to use poisson statistics to estimate the
 147    * number of particles that has hit within a region.
 148    */
 149   void SetUsePoisson(Bool_t u) { fUsePoisson = u; }
 150   /**
 151    * In case of a displaced vertices recalculate eta and angle correction
 152    *
 153    * @param use recalculate or not
 154    *
 155    */
 156   void SetRecalculateEta(Bool_t use) { fRecalculateEta = use; }
 157   /**
 158    * Set whether to use the phi acceptance correction.
 159    *
 160    * How the phi acceptance is used depends on the value passed.
 161    * - 0:  No phi acceptance
 162    * - 1:  Phi acceptance correction done to estimate of particles
 163    * - 2:  Phi acceptance correction done to energy deposited
 164    *
 165    * @param u If >0, use the phi acceptance (default is false)
 166    */
 167
 168   void SetUsePhiAcceptance(UShort_t u=kPhiCorrectNch) { fUsePhiAcceptance = u; }
 169   /**
 170    * Set the luming factors used in the Poisson method
 171    *
 172    * @param eta Must be 1 or larger
 173    * @param phi Must be 1 or larger
 174    */
 175   void SetLumping(Int_t eta, Int_t phi) {
 176     fEtaLumping = (eta < 1 ? 1 : eta);
 177     fPhiLumping = (phi < 1 ? 1 : phi);
 178   }
 179   /**
 180    * Get the multiplicity cut.  If the user has set fMultCut (via
 181    * SetMultCut) then that value is used.  If not, then the lower
 182    * value of the fit range for the enery loss fits is returned.
 183    *
 184    * @return Lower cut on multiplicity
 185    */
 186   Double_t GetMultCut(UShort_t d, Char_t r, Double_t eta,
 187                       Bool_t errors=true) const;
 188   /**
 189    * Get the multiplicity cut.  If the user has set fMultCut (via
 190    * SetMultCut) then that value is used.  If not, then the lower
 191    * value of the fit range for the enery loss fits is returned.
 192    *
 193    * @return Lower cut on multiplicity
 194    */
 195   Double_t GetMultCut(UShort_t d, Char_t r, Int_t ieta,
 196                       Bool_t errors=true) const;
 197   /**
 198    * Print information
 199    *
 200    * @param option Print options
 201    *   - max  Print max weights
 202    */
 203   void Print(Option_t* option="") const;
 204   /**
 205    * Get the cuts used
 206    *
 207    * @return Reference to cuts object
 208    */
 209   AliFMDMultCuts& GetCuts() { return fCuts; }
 210   /**
 211    * Set the cuts to use
 212    *
 213    * @param c Cuts to use
 214    */
 215   void SetCuts(const AliFMDMultCuts& c) { fCuts = c; }
 216 protected:
 217   /**
 218    * Find the max weight to use for FMD<i>dr</i> in eta bin @a iEta
 219    *
 220    * @param cor   Correction
 221    * @param d     Detector
 222    * @param r     Ring
 223    * @param iEta  Eta bin
 224    */
 225   Int_t FindMaxWeight(const AliFMDCorrELossFit* cor,
 226                       UShort_t d, Char_t r, Int_t iEta) const;
 227
 228   /**
 229    * Find the max weights and cache them
 230    *
 231    */
 232   void CacheMaxWeights();
 233   /**
 234    * Find the (cached) maximum weight for FMD<i>dr</i> in
 235    * @f$\eta@f$ bin @a iEta
 236    *
 237    * @param d     Detector
 238    * @param r     Ring
 239    * @param iEta  Eta bin
 240    *
 241    * @return max weight or <= 0 in case of problems
 242    */
 243   Int_t GetMaxWeight(UShort_t d, Char_t r, Int_t iEta) const;
 244   /**
 245    * Find the (cached) maximum weight for FMD<i>dr</i> iat
 246    * @f$\eta@f$
 247    *
 248    * @param d     Detector
 249    * @param r     Ring
 250    * @param eta   Eta bin
 251    *
 252    * @return max weight or <= 0 in case of problems
 253    */
 254   Int_t GetMaxWeight(UShort_t d, Char_t r, Float_t eta) const;
 255
 256   /**
 257    * Get the number of particles corresponding to the signal mult
 258    *
 259    * @param mult     Signal
 260    * @param d        Detector
 261    * @param r        Ring
 262    * @param s        Sector
 263    * @param t        Strip (not used)
 264    * @param v        Vertex bin
 265    * @param eta      Pseudo-rapidity
 266    * @param lowFlux  Low-flux flag
 267    *
 268    * @return The number of particles
 269    */
 270   virtual Float_t NParticles(Float_t mult,
 271                              UShort_t d, Char_t r, UShort_t s, UShort_t t,
 272                              UShort_t v, Float_t eta, Bool_t lowFlux) const;
 273   /**
 274    * Get the inverse correction factor.  This consist of
 275    *
 276    * - acceptance correction (corners of sensors)
 277    * - double hit correction (for low-flux events)
 278    * - dead strip correction
 279    *
 280    * @param d        Detector
 281    * @param r        Ring
 282    * @param s        Sector
 283    * @param t        Strip (not used)
 284    * @param v        Vertex bin
 285    * @param eta      Pseudo-rapidity
 286    * @param lowFlux  Low-flux flag
 287    *
 288    * @return
 289    */
 290   virtual Float_t Correction(UShort_t d, Char_t r, UShort_t s, UShort_t t,
 291                              UShort_t v, Float_t eta, Bool_t lowFlux) const;
 292   /**
 293    * Get the acceptance correction for strip @a t in an ring of type @a r
 294    *
 295    * @param r  Ring type ('I' or 'O')
 296    * @param t  Strip number
 297    *
 298    * @return Inverse acceptance correction
 299    */
 300   virtual Float_t AcceptanceCorrection(Char_t r, UShort_t t) const;
 301   /**
 302    * Generate the acceptance corrections
 303    *
 304    * @param r Ring to generate for
 305    *
 306    * @return Newly allocated histogram of acceptance corrections
 307    */
 308   virtual TH1D*   GenerateAcceptanceCorrection(Char_t r) const;
 309   /**
 310    * Internal data structure to keep track of the histograms
 311    */
 312   struct RingHistos : public AliForwardUtil::RingHistos
 313   {
 314     /**
 315      * Default CTOR
 316      */
 317     RingHistos();
 318     /**
 319      * Constructor
 320      *
 321      * @param d detector
 322      * @param r ring
 323      */
 324     RingHistos(UShort_t d, Char_t r);
 325     /**
 326      * Copy constructor
 327      *
 328      * @param o Object to copy from
 329      */
 330     RingHistos(const RingHistos& o);
 331     /**
 332      * Assignment operator
 333      *
 334      * @param o Object to assign from
 335      *
 336      * @return Reference to this
 337      */
 338     RingHistos& operator=(const RingHistos& o);
 339     /**
 340      * Destructor
 341      */
 342     ~RingHistos();
 343     /**
 344      * Initialize the object
 345      *
 346      * @param eAxis
 347      */
 348     void Init(const TAxis& eAxis);
 349     /**
 350      * Make output
 351      *
 352      * @param dir Where to put it
 353      */
 354     void Output(TList* dir);
 355     /**
 356      * Scale the histograms to the total number of events
 357      *
 358      * @param dir     Where the output is
 359      * @param nEvents Number of events
 360      */
 361     void ScaleHistograms(TList* dir, Int_t nEvents);
 362     TH2D*     fEvsN;           // Correlation of Eloss vs uncorrected Nch
 363     TH2D*     fEvsM;           // Correlation of Eloss vs corrected Nch
 364     TProfile* fEtaVsN;         // Average uncorrected Nch vs eta
 365     TProfile* fEtaVsM;         // Average corrected Nch vs eta
 366     TProfile* fCorr;           // Average correction vs eta
 367     TH2D*     fDensity;        // Distribution inclusive Nch
 368     TH2D*     fELossVsPoisson; // Correlation of energy loss vs Poisson N_ch
 369     AliPoissonCalculator fPoisson; // Calculate density using Poisson method
 370     TH1D*     fELoss;          // Energy loss as seen by this
 371     TH1D*     fELossUsed;      // Energy loss in strips with signal
 372     Double_t  fMultCut;        // If set, use this
 373
 374     ClassDef(RingHistos,6);
 375   };
 376   /**
 377    * Get the ring histogram container
 378    *
 379    * @param d Detector
 380    * @param r Ring
 381    *
 382    * @return Ring histogram container
 383    */
 384   RingHistos* GetRingHistos(UShort_t d, Char_t r) const;
 385   TList    fRingHistos;    //  List of histogram containers
 386   TH1D*    fSumOfWeights;  //  Histogram
 387   TH1D*    fWeightedSum;   //  Histogram
 388   TH1D*    fCorrections;   //  Histogram
 389   UShort_t fMaxParticles;  //  Maximum particle weight to use
 390   Bool_t   fUsePoisson;    //  If true, then use poisson statistics
 391   UShort_t fUsePhiAcceptance; // Whether to correct for corners
 392   TH1D*    fAccI;          //  Acceptance correction for inner rings
 393   TH1D*    fAccO;          //  Acceptance correction for outer rings
 394   TArrayI  fFMD1iMax;      //  Array of max weights
 395   TArrayI  fFMD2iMax;      //  Array of max weights
 396   TArrayI  fFMD2oMax;      //  Array of max weights
 397   TArrayI  fFMD3iMax;      //  Array of max weights
 398   TArrayI  fFMD3oMax;      //  Array of max weights
 399   TH2D*    fMaxWeights;    //  Histogram of max weights
 400   TH2D*    fLowCuts;       //  Histogram of low cuts
 401   Int_t    fEtaLumping;    //  How to lump eta bins for Poisson
 402   Int_t    fPhiLumping;    //  How to lump phi bins for Poisson
 403   Int_t    fDebug;         //  Debug level
 404   AliFMDMultCuts fCuts;    // Cuts
 405   Bool_t fRecalculateEta;  //  //Whether to recalculate eta and angle correction (disp vtx)
 406
 407   ClassDef(AliFMDDensityCalculator,7); // Calculate Nch density
 408 };
 409
 410 #endif
 411 // Local Variables:
 412 //   mode: C++
 413 // End:
 414