PWGLF/FORWARD/analysis2/AliFMDHistCollector.h

   1 //
   2 // This class collects the event histograms into single histograms,
   3 // one for each ring in each vertex bin.
   4 //
   5 #ifndef ALIFMDHISTCOLLECTOR_H
   6 #define ALIFMDHISTCOLLECTOR_H
   7 /**
   8  * @file   AliFMDHistCollector.h
   9  * @author Christian Holm Christensen <cholm@dalsgaard.hehi.nbi.dk>
  10  * @date   Wed Mar 23 14:03:01 2011
  11  *
  12  * @brief
  13  *
  14  *
  15  * @ingroup pwglf_forward_aod
  16  */
  17 #include <TNamed.h>
  18 #include <TList.h>
  19 #include <TArrayI.h>
  20 #include "AliForwardUtil.h"
  21 class AliESDFMD;
  22 class TH2D;
  23
  24 /**
  25  * This class collects the event histograms into single histograms,
  26  * one for each ring in each vertex bin.
  27  *
  28  * @par Input:
  29  *   - AliESDFMD object possibly corrected for sharing
  30  *
  31  * @par Output:
  32  *   - 5 RingHistos objects - each with a number of vertex dependent
  33  *     2D histograms of the inclusive charge particle density
  34  *
  35  * @par HistCollector used:
  36  *   - AliFMDCorrSecondaryMap
  37  *
  38  * @ingroup pwglf_forward_algo
  39  * @ingroup pwglf_forward_aod
  40  */
  41 class AliFMDHistCollector : public TNamed
  42 {
  43 public:
  44   /**
  45    * Methods to use when merging overlapping bins @f$b_1@f$, @f$b_2@f$
  46    * with content @f$c_1@f$, @f$c_2@f$, and errors @f$e_1@f$,
  47    * @f$e_2@f$ into bin @f$ b@f$ with content @f$c@f$ and error @f$e@f$
  48    */
  49   enum MergeMethod {
  50     /**
  51      * @f[
  52      *   c = \frac{1}{2}(c_1+c_2)
  53      * @f]
  54      * @f[
  55      *   e = \sqrt{e_1^2+e_2^2}
  56      * @f]
  57      */
  58     kStraightMean,
  59     /**
  60      * As above, exept zero's are ignored
  61      */
  62     kStraightMeanNoZero,
  63     /**
  64      * @f[
  65      *   c = \frac{\frac{c_1}{e_1^2}+\frac{c_2}{e_2^2}}{
  66      *             \frac{1}{e_1^2}+\frac{1}{e_2^2}}
  67      * @f]
  68      * @f[
  69      *   e = \sqrt{\frac{1}{\frac{1}{e_1^2}+\frac{1}{e_2^2}}}
  70      * @f]
  71      */
  72     kWeightedMean,
  73     /**
  74      * @f[
  75      *     c = \left\{\begin{array}{cl}
  76      *          c_1 & \mbox{if $e_1 < e_2$} \\
  77      *          c_2 & \mbox{otherwise}\end{array}\right.
  78      * @f]
  79      */
  80     kLeastError
  81   };
  82   /**
  83    * How to obtain the fiducial cuts
  84    */
  85   enum FiducialMethod {
  86     /**
  87      * Select bins by fixed cut.  Bins with a secondary correction
  88      * less than the cut is considered as non-valid
  89      */
  90     kByCut,
  91     /**
  92      * A bin is considered non-valid, if it is less then twice as
  93      * large as it's neighbors (in eta)
  94      */
  95     kDistance
  96   };
  97   /**
  98    * Constructor
  99    */
 100   AliFMDHistCollector();
 101   /**
 102    * Constructor
 103    *
 104    * @param title Name of object
 105    */
 106   AliFMDHistCollector(const char* title);
 107   /**
 108    * Copy constructor
 109    *
 110    * @param o Object to copy from
 111    */
 112   AliFMDHistCollector(const AliFMDHistCollector& o);
 113
 114   /**
 115    * Destructor
 116    */
 117   virtual ~AliFMDHistCollector();
 118   /**
 119    * Assignement operator
 120    *
 121    * @param o Object to assign from
 122    *
 123    * @return Reference to this object
 124    */
 125   AliFMDHistCollector& operator=(const AliFMDHistCollector&);
 126   /**
 127    * Intialise
 128    *
 129    * @param vtxAxis  @f$ v_z@f$ axis
 130    * @param etaAxis  @f$ \eta@f$ axis
 131    */
 132   virtual void Init(const TAxis& vtxAxis,
 133                     const TAxis& etaAxis);
 134   /**
 135    * Do the calculations
 136    *
 137    * @param hists    Cache of histograms
 138    * @param sums     Cache to sum ring histograms in
 139    * @param vtxBin   Vertex bin (1 based)
 140    * @param out      Output histogram
 141    *
 142    * @return true on successs
 143    */
 144   virtual Bool_t Collect(const AliForwardUtil::Histos& hists,
 145                          AliForwardUtil::Histos&       sums,
 146                          UShort_t                      vtxBin,
 147                          TH2D&                         out);
 148   /**
 149    * Output diagnostic histograms to directory
 150    *
 151    * @param dir List to write in
 152    */
 153   virtual void DefineOutput(TList* dir);
 154   /**
 155    * Set the merge method
 156    *
 157    * @param m Method
 158    */
 159   void SetMergeMethod(MergeMethod m) { fMergeMethod = m; }
 160   /**
 161    * Set the method for finding the fidicual area of the secondary maps
 162    *
 163    * @param m Method
 164    */
 165   void SetFiducialMethod(FiducialMethod m) { fFiducialMethod = m; }
 166   /**
 167    * Set the number of extra bins (beyond the secondary map border)
 168    * to cut away.
 169    *
 170    * @param n Number of bins
 171    */
 172   void SetNCutBins(UInt_t n=2) { fNCutBins = n; }
 173   /**
 174    * Set the correction cut, that is, when bins in the secondary
 175    * correction maps have a value lower than this cut, they are
 176    * considered uncertain and not used
 177    *
 178    * @param cut Cut-off
 179    */
 180   void SetCorrectionCut(Float_t cut=0.5) { fCorrectionCut = cut; }
 181   /**
 182    * Set the debug level.  The higher the value the more output
 183    *
 184    * @param dbg Debug level
 185    */
 186   void SetDebug(Int_t dbg=1) { fDebug = dbg; }
 187   /**
 188    * Print information
 189    *
 190    * @param option Not used
 191    */
 192   void Print(Option_t* option="") const;
 193 protected:
 194   /**
 195    * Get the first and last eta bin to use for a given ring and vertex
 196    *
 197    * @param d        Detector
 198    * @param r        Ring
 199    * @param vtxBin   Vertex bin (1 based)
 200    * @param first    On return, the first eta bin to use
 201    * @param last     On return, the last eta bin to use
 202    */
 203   virtual void GetFirstAndLast(UShort_t d, Char_t r, UShort_t vtxBin,
 204                                Int_t& first, Int_t& last) const;
 205   /**
 206    * Get the first and last eta bin to use for a given ring and vertex
 207    *
 208    * @param idx      Ring index as given by GetIdx
 209    * @param vtxBin   Vertex bin (1 based)
 210    * @param first    On return, the first eta bin to use
 211    * @param last     On return, the last eta bin to use
 212    */
 213   virtual void GetFirstAndLast(Int_t idx, UShort_t vtxBin,
 214                                Int_t& first, Int_t& last) const;
 215   /**
 216    * Get the first eta bin to use for a given ring and vertex
 217    *
 218    * @param d Detector
 219    * @param r Ring
 220    * @param v vertex bin (1 based)
 221    *
 222    * @return First eta bin to use, or -1 in case of problems
 223    */
 224   Int_t GetFirst(UShort_t d, Char_t r, UShort_t v) const;
 225   /**
 226    * Get the first eta bin to use for a given ring and vertex
 227    *
 228    * @param idx Ring index as given by GetIdx
 229    * @param v vertex bin (1 based)
 230    *
 231    * @return First eta bin to use, or -1 in case of problems
 232    */
 233   Int_t GetFirst(Int_t idx, UShort_t v) const;
 234   /**
 235    * Get the last eta bin to use for a given ring and vertex
 236    *
 237    * @param d Detector
 238    * @param r Ring
 239    * @param v vertex bin (1 based)
 240    *
 241    * @return Last eta bin to use, or -1 in case of problems
 242    */
 243   Int_t GetLast(UShort_t d, Char_t r, UShort_t v) const;
 244   /**
 245    * Get the last eta bin to use for a given ring and vertex
 246    *
 247    * @param idx Ring index as given by GetIdx
 248    * @param v vertex bin (1 based)
 249    *
 250    * @return Last eta bin to use, or -1 in case of problems
 251    */
 252   Int_t GetLast(Int_t idx, UShort_t v) const;
 253   /**
 254    * Get the detector and ring from the ring index
 255    *
 256    * @param idx Ring index
 257    * @param d   On return, the detector or 0 in case of errors
 258    * @param r   On return, the ring id or '0' in case of errors
 259    */
 260   void GetDetRing(Int_t idx, UShort_t& d, Char_t& r) const;
 261   /**
 262    * Get the ring index from detector number and ring identifier
 263    *
 264    * @param d Detector
 265    * @param r Ring identifier
 266    *
 267    * @return ring index or -1 in case of problems
 268    */
 269   Int_t GetIdx(UShort_t d, Char_t r) const;
 270   /**
 271    * Get the possibly overlapping histogram of eta bin @a e in
 272    * detector and ring
 273    *
 274    * @param d Detector
 275    * @param r Ring
 276    * @param e Eta bin
 277    * @param v Vertex bin (1 based)
 278    *
 279    * @return Overlapping histogram index or -1
 280    */
 281   Int_t GetOverlap(UShort_t d, Char_t r, Int_t e, UShort_t v) const;
 282   /**
 283    * Get the possibly overlapping histogram of eta bin @a e in
 284    * detector and ring
 285    *
 286    * @param i Ring index
 287    * @param e Eta bin
 288    * @param v Vertex bin (1 based)
 289    *
 290    * @return Overlapping histogram index or -1
 291    */
 292   Int_t GetOverlap(Int_t i, Int_t e, UShort_t v) const;
 293   /**
 294    * Check if there's an overlapping histogram with this eta bin of
 295    * the detector and ring
 296    *
 297    * @param d Detector
 298    * @param r Ring
 299    * @param e eta bin
 300    * @param v Vertex bin (1 based)
 301    *
 302    * @return True if there's an overlapping histogram
 303    */
 304   Bool_t HasOverlap(UShort_t d, Char_t r, Int_t e, UShort_t v) const;
 305   /**
 306    * Check if there's an overlapping histogram with this eta bin of
 307    * ring
 308    *
 309    * @param i Ring index
 310    * @param e eta bin
 311    * @param v Vertex bin
 312    *
 313    * @return True if there's an overlapping histogram
 314    */
 315   Bool_t HasOverlap(Int_t i, Int_t e, UShort_t v) const;
 316   /**
 317    * Check if we should include the bin in the data range
 318    *
 319    * @param bg Secondary map histogram
 320    * @param ie Eta bin
 321    * @param ip Phi bin
 322    *
 323    * @return True if to be used
 324    */
 325   Bool_t CheckCorrection(const TH2D* bg, Int_t ie, Int_t ip) const;
 326   /**
 327    * Merge bins accoring to set method
 328    *
 329    * @param c   Current content
 330    * @param e   Current error
 331    * @param oc  Old content
 332    * @param oe  Old error
 333    * @param rc  On return, the new content
 334    * @param re  On return, tne new error
 335    */
 336   void MergeBins(Double_t c,   Double_t e,
 337                  Double_t oc,  Double_t oe,
 338                  Double_t& rc, Double_t& re) const;
 339
 340
 341   Int_t       fNCutBins;        // Number of additional bins to cut away
 342   Float_t     fCorrectionCut;   // Cut-off on secondary corrections
 343   TArrayI     fFirstBins;       // Array of first eta bins
 344   TArrayI     fLastBins;        // Array of last eta bins
 345   Int_t       fDebug;           // Debug level
 346   TList*      fList;            // Output list
 347   TH2D*       fSumRings;        // Sum per ring (on y-axis)
 348   TH2D*       fCoverage;        // Sum per ring (on y-axis)
 349   MergeMethod fMergeMethod;     // Merge methiod for overlapping bins
 350   FiducialMethod fFiducialMethod; // Fidicual method
 351
 352   ClassDef(AliFMDHistCollector,1); // Calculate Nch density
 353 };
 354
 355 //____________________________________________________________________
 356 inline void
 357 AliFMDHistCollector::GetFirstAndLast(UShort_t d, Char_t r, UShort_t vtxbin,
 358                                      Int_t& first, Int_t& last) const
 359 {
 360   GetFirstAndLast(GetIdx(d,r), vtxbin, first, last);
 361 }
 362 //____________________________________________________________________
 363 inline Int_t
 364 AliFMDHistCollector::GetFirst(UShort_t d, Char_t r, UShort_t v) const
 365 {
 366   return GetFirst(GetIdx(d,r), v);
 367 }
 368 //____________________________________________________________________
 369 inline Int_t
 370 AliFMDHistCollector::GetLast(UShort_t d, Char_t r, UShort_t v) const
 371 {
 372   return GetLast(GetIdx(d, r), v);
 373 }
 374 //____________________________________________________________________
 375 inline Bool_t
 376 AliFMDHistCollector::HasOverlap(UShort_t d, Char_t r, Int_t e, UShort_t v) const
 377 {
 378   return GetOverlap(d,r,e,v) >= 0;
 379 }
 380 //____________________________________________________________________
 381 inline Bool_t
 382 AliFMDHistCollector::HasOverlap(Int_t i, Int_t e, UShort_t v) const
 383 {
 384   return GetOverlap(i,e,v) >= 0;
 385 }
 386
 387 #endif
 388 // Local Variables:
 389 //   mode: C++
 390 // End:
 391