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 #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 pwg2_forward_algo
  45  * @ingroup pwg2_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);
 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    * Set whether to use the phi acceptance correction.
 152    *
 153    * How the phi acceptance is used depends on the value passed.
 154    * - 0:  No phi acceptance
 155    * - 1:  Phi acceptance correction done to estimate of particles
 156    * - 2:  Phi acceptance correction done to energy deposited
 157    *
 158    * @param u If >0, use the phi acceptance (default is false)
 159    */
 160   void SetUsePhiAcceptance(UShort_t u=kPhiCorrectNch) { fUsePhiAcceptance = u; }
 161   /**
 162    * Set the luming factors used in the Poisson method
 163    *
 164    * @param eta Must be 1 or larger
 165    * @param phi Must be 1 or larger
 166    */
 167   void SetLumping(Int_t eta, Int_t phi) {
 168     fEtaLumping = (eta < 1 ? 1 : eta);
 169     fPhiLumping = (phi < 1 ? 1 : phi);
 170   }
 171   /**
 172    * Get the multiplicity cut.  If the user has set fMultCut (via
 173    * SetMultCut) then that value is used.  If not, then the lower
 174    * value of the fit range for the enery loss fits is returned.
 175    *
 176    * @return Lower cut on multiplicity
 177    */
 178   Double_t GetMultCut(UShort_t d, Char_t r, Double_t eta,
 179                       Bool_t errors=true) const;
 180   /**
 181    * Get the multiplicity cut.  If the user has set fMultCut (via
 182    * SetMultCut) then that value is used.  If not, then the lower
 183    * value of the fit range for the enery loss fits is returned.
 184    *
 185    * @return Lower cut on multiplicity
 186    */
 187   Double_t GetMultCut(UShort_t d, Char_t r, Int_t ieta,
 188                       Bool_t errors=true) const;
 189   /**
 190    * Print information
 191    *
 192    * @param option Print options
 193    *   - max  Print max weights
 194    */
 195   void Print(Option_t* option="") const;
 196   /**
 197    * Get the cuts used
 198    *
 199    * @return Reference to cuts object
 200    */
 201   AliFMDMultCuts& GetCuts() { return fCuts; }
 202   /**
 203    * Set the cuts to use
 204    *
 205    * @param c Cuts to use
 206    */
 207   void SetCuts(const AliFMDMultCuts& c) { fCuts = c; }
 208 protected:
 209   /**
 210    * Find the max weight to use for FMD<i>dr</i> in eta bin @a iEta
 211    *
 212    * @param cor   Correction
 213    * @param d     Detector
 214    * @param r     Ring
 215    * @param iEta  Eta bin
 216    */
 217   Int_t FindMaxWeight(const AliFMDCorrELossFit* cor,
 218                       UShort_t d, Char_t r, Int_t iEta) const;
 219
 220   /**
 221    * Find the max weights and cache them
 222    *
 223    */
 224   void CacheMaxWeights();
 225   /**
 226    * Find the (cached) maximum weight for FMD<i>dr</i> in
 227    * @f$\eta@f$ bin @a iEta
 228    *
 229    * @param d     Detector
 230    * @param r     Ring
 231    * @param iEta  Eta bin
 232    *
 233    * @return max weight or <= 0 in case of problems
 234    */
 235   Int_t GetMaxWeight(UShort_t d, Char_t r, Int_t iEta) const;
 236   /**
 237    * Find the (cached) maximum weight for FMD<i>dr</i> iat
 238    * @f$\eta@f$
 239    *
 240    * @param d     Detector
 241    * @param r     Ring
 242    * @param eta   Eta bin
 243    *
 244    * @return max weight or <= 0 in case of problems
 245    */
 246   Int_t GetMaxWeight(UShort_t d, Char_t r, Float_t eta) const;
 247
 248   /**
 249    * Get the number of particles corresponding to the signal mult
 250    *
 251    * @param mult     Signal
 252    * @param d        Detector
 253    * @param r        Ring
 254    * @param s        Sector
 255    * @param t        Strip (not used)
 256    * @param v        Vertex bin
 257    * @param eta      Pseudo-rapidity
 258    * @param lowFlux  Low-flux flag
 259    *
 260    * @return The number of particles
 261    */
 262   virtual Float_t NParticles(Float_t mult,
 263                              UShort_t d, Char_t r, UShort_t s, UShort_t t,
 264                              UShort_t v, Float_t eta, Bool_t lowFlux) const;
 265   /**
 266    * Get the inverse correction factor.  This consist of
 267    *
 268    * - acceptance correction (corners of sensors)
 269    * - double hit correction (for low-flux events)
 270    * - dead strip correction
 271    *
 272    * @param d        Detector
 273    * @param r        Ring
 274    * @param s        Sector
 275    * @param t        Strip (not used)
 276    * @param v        Vertex bin
 277    * @param eta      Pseudo-rapidity
 278    * @param lowFlux  Low-flux flag
 279    *
 280    * @return
 281    */
 282   virtual Float_t Correction(UShort_t d, Char_t r, UShort_t s, UShort_t t,
 283                              UShort_t v, Float_t eta, Bool_t lowFlux) const;
 284   /**
 285    * Get the acceptance correction for strip @a t in an ring of type @a r
 286    *
 287    * @param r  Ring type ('I' or 'O')
 288    * @param t  Strip number
 289    *
 290    * @return Inverse acceptance correction
 291    */
 292   virtual Float_t AcceptanceCorrection(Char_t r, UShort_t t) const;
 293   /**
 294    * Generate the acceptance corrections
 295    *
 296    * @param r Ring to generate for
 297    *
 298    * @return Newly allocated histogram of acceptance corrections
 299    */
 300   virtual TH1D*   GenerateAcceptanceCorrection(Char_t r) const;
 301   /**
 302    * Internal data structure to keep track of the histograms
 303    */
 304   struct RingHistos : public AliForwardUtil::RingHistos
 305   {
 306     /**
 307      * Default CTOR
 308      */
 309     RingHistos();
 310     /**
 311      * Constructor
 312      *
 313      * @param d detector
 314      * @param r ring
 315      */
 316     RingHistos(UShort_t d, Char_t r);
 317     /**
 318      * Copy constructor
 319      *
 320      * @param o Object to copy from
 321      */
 322     RingHistos(const RingHistos& o);
 323     /**
 324      * Assignment operator
 325      *
 326      * @param o Object to assign from
 327      *
 328      * @return Reference to this
 329      */
 330     RingHistos& operator=(const RingHistos& o);
 331     /**
 332      * Destructor
 333      */
 334     ~RingHistos();
 335     /**
 336      * Initialize the object
 337      *
 338      * @param eAxis
 339      */
 340     void Init(const TAxis& eAxis);
 341     /**
 342      * Make output
 343      *
 344      * @param dir Where to put it
 345      */
 346     void Output(TList* dir);
 347     /**
 348      * Scale the histograms to the total number of events
 349      *
 350      * @param dir     Where the output is
 351      * @param nEvents Number of events
 352      */
 353     void ScaleHistograms(TList* dir, Int_t nEvents);
 354     TH2D*     fEvsN;           // Correlation of Eloss vs uncorrected Nch
 355     TH2D*     fEvsM;           // Correlation of Eloss vs corrected Nch
 356     TProfile* fEtaVsN;         // Average uncorrected Nch vs eta
 357     TProfile* fEtaVsM;         // Average corrected Nch vs eta
 358     TProfile* fCorr;           // Average correction vs eta
 359     TH2D*     fDensity;        // Distribution inclusive Nch
 360     TH2D*     fELossVsPoisson; // Correlation of energy loss vs Poisson N_ch
 361     AliPoissonCalculator fPoisson; // Calculate density using Poisson method
 362     TH1D*     fELoss;          // Energy loss as seen by this
 363     TH1D*     fELossUsed;      // Energy loss in strips with signal
 364     Double_t  fMultCut;        // If set, use this
 365
 366     ClassDef(RingHistos,6);
 367   };
 368   /**
 369    * Get the ring histogram container
 370    *
 371    * @param d Detector
 372    * @param r Ring
 373    *
 374    * @return Ring histogram container
 375    */
 376   RingHistos* GetRingHistos(UShort_t d, Char_t r) const;
 377   TList    fRingHistos;    //  List of histogram containers
 378   TH1D*    fSumOfWeights;  //  Histogram
 379   TH1D*    fWeightedSum;   //  Histogram
 380   TH1D*    fCorrections;   //  Histogram
 381   UShort_t fMaxParticles;  //  Maximum particle weight to use
 382   Bool_t   fUsePoisson;    //  If true, then use poisson statistics
 383   UShort_t fUsePhiAcceptance; // Whether to correct for corners
 384   TH1D*    fAccI;          //  Acceptance correction for inner rings
 385   TH1D*    fAccO;          //  Acceptance correction for outer rings
 386   TArrayI  fFMD1iMax;      //  Array of max weights
 387   TArrayI  fFMD2iMax;      //  Array of max weights
 388   TArrayI  fFMD2oMax;      //  Array of max weights
 389   TArrayI  fFMD3iMax;      //  Array of max weights
 390   TArrayI  fFMD3oMax;      //  Array of max weights
 391   TH2D*    fMaxWeights;    //  Histogram of max weights
 392   TH2D*    fLowCuts;       //  Histogram of low cuts
 393   Int_t    fEtaLumping;    //  How to lump eta bins for Poisson
 394   Int_t    fPhiLumping;    //  How to lump phi bins for Poisson
 395   Int_t    fDebug;         //  Debug level
 396   AliFMDMultCuts fCuts;    // Cuts
 397
 398   ClassDef(AliFMDDensityCalculator,7); // Calculate Nch density
 399 };
 400
 401 #endif
 402 // Local Variables:
 403 //   mode: C++
 404 // End:
 405