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