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