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   void SetDoTiming(Bool_t enable=true) { fDoTiming = enable; }
 139   /**
 140    * Maximum particle weight to use
 141    *
 142    * @param m
 143    */
 144   void SetMaxParticles(UShort_t m) { fMaxParticles = m; }
 145   /**
 146    * Set whether to use poisson statistics to estimate the
 147    * number of particles that has hit within a region.  If this is true,
 148    * then the average charge particle density is given by
 149    * @f[
 150    *  \lambda = -\log\left(\frac{N_e}{N_t}\right)
 151    * @f]
 152    * where $N_e$ is the number of strips within the region that has no
 153    * hits over threshold, and $N_t$ is the total number of strips in the
 154    * region/
 155    *
 156    * @param u Whether to use poisson statistics to estimate the
 157    * number of particles that has hit within a region.
 158    */
 159   void SetUsePoisson(Bool_t u) { fUsePoisson = u; }
 160   /**
 161    * In case of a displaced vertices recalculate eta and angle correction
 162    *
 163    * @param use recalculate or not
 164    *
 165    */
 166   void SetRecalculatePhi(Bool_t use) { fRecalculatePhi = use; }
 167   /**
 168    * Set whether to use the phi acceptance correction.
 169    *
 170    * How the phi acceptance is used depends on the value passed.
 171    * - 0:  No phi acceptance
 172    * - 1:  Phi acceptance correction done to estimate of particles
 173    * - 2:  Phi acceptance correction done to energy deposited
 174    *
 175    * @param u If >0, use the phi acceptance (default is false)
 176    */
 177
 178   void SetUsePhiAcceptance(UShort_t u=kPhiCorrectNch) { fUsePhiAcceptance = u; }
 179   /**
 180    * Set the luming factors used in the Poisson method
 181    *
 182    * @param eta Must be 1 or larger
 183    * @param phi Must be 1 or larger
 184    */
 185   void SetLumping(Int_t eta, Int_t phi) {
 186     fEtaLumping = (eta < 1 ? 1 : eta);
 187     fPhiLumping = (phi < 1 ? 1 : phi);
 188   }
 189   /**
 190    * Set the minimum quality of the energy loss fits
 191    *
 192    * @param cut Cut value
 193    */
 194   void SetMinQuality(UShort_t cut=10) { fMinQuality = cut; }
 195   /**
 196    * Set the maximum ratio of outlier bins to the total number of bins
 197    * with data.
 198    *
 199    * @param ratio Maximum ratio (number between 0 and 1)
 200    */
 201   void SetMaxOutliers(Double_t ratio=0.10) { fMaxOutliers = ratio; }
 202   /**
 203    * Set the maximum relative diviation between @f$N_{ch}^{Poisson}@f$
 204    * and @f$N_{ch}^{\Delta}@f$
 205    *
 206    * @param cut Relative cut (number between 0 and 1)
 207    */
 208   void SetOutlierCut(Double_t cut=0.50) { fOutlierCut = cut; }
 209   /**
 210    * Get the multiplicity cut.  If the user has set fMultCut (via
 211    * SetMultCut) then that value is used.  If not, then the lower
 212    * value of the fit range for the enery loss fits is returned.
 213    *
 214    * @param d      Detector
 215    * @param r      Ring
 216    * @param eta    Psuedo-rapidity
 217    * @param errors Factor in errors
 218    *
 219    * @return Lower cut on multiplicity
 220    */
 221   Double_t GetMultCut(UShort_t d, Char_t r, Double_t eta,
 222                       Bool_t errors=true) const;
 223   /**
 224    * Get the multiplicity cut.  If the user has set fMultCut (via
 225    * SetMultCut) then that value is used.  If not, then the lower
 226    * value of the fit range for the enery loss fits is returned.
 227    *
 228    * @param d      Detector
 229    * @param r      Ring
 230    * @param ieta   Psuedo-rapidity bin
 231    * @param errors Factor in errors
 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    * Set the minimum quality of the energy loss fits
 239    *
 240    * @return Cut value
 241    */
 242   UShort_t GetMinQuality() const { return fMinQuality; }
 243   /**
 244    * Print information
 245    *
 246    * @param option Print options
 247    *   - max  Print max weights
 248    */
 249   void Print(Option_t* option="") const;
 250   /**
 251    * Get the cuts used
 252    *
 253    * @return Reference to cuts object
 254    */
 255   AliFMDMultCuts& GetCuts() { return fCuts; }
 256   /**
 257    * Set the cuts to use
 258    *
 259    * @param c Cuts to use
 260    */
 261   void SetCuts(const AliFMDMultCuts& c) { fCuts = c; }
 262 protected:
 263   /**
 264    * Find the max weight to use for FMD<i>dr</i> in eta bin @a iEta
 265    *
 266    * @param cor   Correction
 267    * @param d     Detector
 268    * @param r     Ring
 269    * @param iEta  Eta bin
 270    *
 271    * @return The maximum weight
 272    */
 273   Int_t FindMaxWeight(const AliFMDCorrELossFit* cor,
 274                       UShort_t d, Char_t r, Int_t iEta) const;
 275
 276   /**
 277    * Find the max weight to use for FMD<i>dr</i> in eta @a eta
 278    *
 279    * @param cor   Correction
 280    * @param d     Detector
 281    * @param r     Ring
 282    * @param eta   Eta
 283    *
 284    * @return The maximum weight
 285    */
 286   Int_t FindMaxWeight(const AliFMDCorrELossFit* cor,
 287                       UShort_t d, Char_t r, Double_t iEta) const;
 288
 289   /**
 290    * Find the max weights and cache them
 291    *
 292    * @param axis Default @f$\eta@f$ axis from parent task
 293    */
 294   void CacheMaxWeights(const TAxis& axis);
 295   /**
 296    * Find the (cached) maximum weight for FMD<i>dr</i> in
 297    * @f$\eta@f$ bin @a iEta
 298    *
 299    * @param d     Detector
 300    * @param r     Ring
 301    * @param iEta  Eta bin
 302    *
 303    * @return max weight or <= 0 in case of problems
 304    */
 305   Int_t GetMaxWeight(UShort_t d, Char_t r, Int_t iEta) const;
 306   /**
 307    * Find the (cached) maximum weight for FMD<i>dr</i> iat
 308    * @f$\eta@f$
 309    *
 310    * @param d     Detector
 311    * @param r     Ring
 312    * @param eta   Eta bin
 313    *
 314    * @return max weight or <= 0 in case of problems
 315    */
 316   Int_t GetMaxWeight(UShort_t d, Char_t r, Float_t eta) const;
 317
 318   /**
 319    * Get the number of particles corresponding to the signal mult
 320    *
 321    * @param mult     Signal
 322    * @param d        Detector
 323    * @param r        Ring
 324    * @param eta      Pseudo-rapidity
 325    * @param lowFlux  Low-flux flag
 326    *
 327    * @return The number of particles
 328    */
 329   virtual Float_t NParticles(Float_t  mult,
 330                              UShort_t d,
 331                              Char_t   r,
 332                              Float_t  eta,
 333                              Bool_t   lowFlux) const;
 334   /**
 335    * Get the inverse correction factor.  This consist of
 336    *
 337    * - acceptance correction (corners of sensors)
 338    * - double hit correction (for low-flux events)
 339    * - dead strip correction
 340    *
 341    * @param d        Detector
 342    * @param r        Ring
 343    * @param t        Strip
 344    * @param eta      Pseudo-rapidity
 345    * @param lowFlux  Low-flux flag
 346    *
 347    * @return the correction factor
 348    */
 349   virtual Float_t Correction(UShort_t d, Char_t r, UShort_t t,
 350                              Float_t eta, Bool_t lowFlux) const;
 351   /**
 352    * Get the acceptance correction for strip @a t in an ring of type @a r
 353    *
 354    * @param r  Ring type ('I' or 'O')
 355    * @param t  Strip number
 356    *
 357    * @return Inverse acceptance correction
 358    */
 359   virtual Float_t AcceptanceCorrection(Char_t r, UShort_t t) const;
 360   /**
 361    * Generate the acceptance corrections
 362    *
 363    * @param r Ring to generate for
 364    *
 365    * @return Newly allocated histogram of acceptance corrections
 366    */
 367   virtual TH1D*   GenerateAcceptanceCorrection(Char_t r) const;
 368   /**
 369    * Check if, for a given region, whether this is an outlier
 370    *
 371    * The condition for an outlier event are
 372    * @f[
 373    * |N_{ch}^{Poisson} - N_{ch}^{\Delta}| / N_{ch}^{\Delta} > c
 374    * @f]
 375    *
 376    * @param eloss    @f$ N_{ch}^{\Delta}@f$ - number of charged particles
 377    * @param poisson  @f$ N_{ch}^{Poisson}@f$ - number of charged particles
 378    * @param cut      @f$ c@f$ - the cut
 379    *
 380    * @return true if the region reflects an outlier event
 381    */
 382   virtual Bool_t CheckOutlier(Double_t eloss,
 383                               Double_t poisson,
 384                               Double_t cut=0.5) const;
 385   /**
 386    * Internal data structure to keep track of the histograms
 387    */
 388   struct RingHistos : public AliForwardUtil::RingHistos
 389   {
 390     /**
 391      * Default CTOR
 392      */
 393     RingHistos();
 394     /**
 395      * Constructor
 396      *
 397      * @param d detector
 398      * @param r ring
 399      */
 400     RingHistos(UShort_t d, Char_t r);
 401     /**
 402      * Copy constructor
 403      *
 404      * @param o Object to copy from
 405      */
 406     RingHistos(const RingHistos& o);
 407     /**
 408      * Assignment operator
 409      *
 410      * @param o Object to assign from
 411      *
 412      * @return Reference to this
 413      */
 414     RingHistos& operator=(const RingHistos& o);
 415     /**
 416      * Destructor
 417      */
 418     ~RingHistos();
 419     /**
 420      * Initialize the object
 421      *
 422      * @param eAxis
 423      */
 424     void SetupForData(const TAxis& eAxis);
 425     /**
 426      * Make output
 427      *
 428      * @param dir Where to put it
 429      */
 430     void CreateOutputObjects(TList* dir);
 431     /**
 432      * Scale the histograms to the total number of events
 433      *
 434      * @param dir     Where the output is
 435      * @param nEvents Number of events
 436      */
 437     void Terminate(TList* dir, Int_t nEvents);
 438     TList*    fList;
 439     // TH2D*     fEvsN;           // Correlation of Eloss vs uncorrected Nch
 440     // TH2D*     fEvsM;           // Correlation of Eloss vs corrected Nch
 441     // TProfile* fEtaVsN;         // Average uncorrected Nch vs eta
 442     // TProfile* fEtaVsM;         // Average corrected Nch vs eta
 443     TProfile* fCorr;           // Average correction vs eta
 444     TH2D*     fDensity;        // Distribution inclusive Nch
 445     TH2D*     fELossVsPoisson; // Correlation of energy loss vs Poisson N_ch
 446     TH1D*     fDiffELossPoisson;// Relative difference to Poisson
 447     TH2D*     fELossVsPoissonOut; // Correlation of energy loss vs Poisson N_ch
 448     TH1D*     fDiffELossPoissonOut;// Relative difference to Poisson
 449     TH1D*     fOutliers;       // Fraction of outliers per event
 450     AliPoissonCalculator fPoisson; // Calculate density using Poisson method
 451     TH1D*     fELoss;          // Energy loss as seen by this
 452     TH1D*     fELossUsed;      // Energy loss in strips with signal
 453     Double_t  fMultCut;        // If set, use this
 454     TH1D*     fTotal;          // Total number of strips per eta
 455     TH1D*     fGood;           // Number of good strips per eta
 456     TH2D*     fPhiAcc;         // Phi acceptance vs IpZ
 457     TH1D*     fPhiBefore;      // Phi before re-calce
 458     TH1D*     fPhiAfter;       // Phi after re-calc
 459     ClassDef(RingHistos,10);
 460   };
 461   /**
 462    * Get the ring histogram container
 463    *
 464    * @param d Detector
 465    * @param r Ring
 466    *
 467    * @return Ring histogram container
 468    */
 469   RingHistos* GetRingHistos(UShort_t d, Char_t r) const;
 470   TList    fRingHistos;    //  List of histogram containers
 471   TH1D*    fSumOfWeights;  //  Histogram
 472   TH1D*    fWeightedSum;   //  Histogram
 473   TH1D*    fCorrections;   //  Histogram
 474   UShort_t fMaxParticles;  //  Maximum particle weight to use
 475   Bool_t   fUsePoisson;    //  If true, then use poisson statistics
 476   UShort_t fUsePhiAcceptance; // Whether to correct for corners
 477   TH1D*    fAccI;          //  Acceptance correction for inner rings
 478   TH1D*    fAccO;          //  Acceptance correction for outer rings
 479   TArrayI  fFMD1iMax;      //  Array of max weights
 480   TArrayI  fFMD2iMax;      //  Array of max weights
 481   TArrayI  fFMD2oMax;      //  Array of max weights
 482   TArrayI  fFMD3iMax;      //  Array of max weights
 483   TArrayI  fFMD3oMax;      //  Array of max weights
 484   TH2D*    fMaxWeights;    //  Histogram of max weights
 485   TH2D*    fLowCuts;       //  Histogram of low cuts
 486   Int_t    fEtaLumping;    //  How to lump eta bins for Poisson
 487   Int_t    fPhiLumping;    //  How to lump phi bins for Poisson
 488   Int_t    fDebug;         //  Debug level
 489   AliFMDMultCuts fCuts;    // Cuts
 490   Bool_t   fRecalculatePhi;  // Whether to correct for (X,Y) offset
 491   UShort_t fMinQuality;      // Least quality for fits
 492   AliForwardUtil::Histos fCache;
 493   Bool_t                 fDoTiming;
 494   TProfile*              fHTiming;
 495   Double_t               fMaxOutliers; // Maximum ratio of outlier bins
 496   Double_t               fOutlierCut;  // Maximum relative diviation
 497
 498   ClassDef(AliFMDDensityCalculator,14); // Calculate Nch density
 499 };
 500
 501 #endif
 502 // Local Variables:
 503 //   mode: C++
 504 // End:
 505