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