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