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