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      * Just sum the signals
  83      */
  84     kSum
  85   };
  86   /**
  87    * How to obtain the fiducial cuts
  88    */
  89   enum FiducialMethod {
  90     /**
  91      * Select bins by fixed cut.  Bins with a secondary correction
  92      * less than the cut is considered as non-valid
  93      */
  94     kByCut,
  95     /**
  96      * A bin is considered non-valid, if it is less then twice as
  97      * large as it's neighbors (in eta)
  98      */
  99     kDistance
 100   };
 101   /**
 102    * FMD ring bits for skipping
 103    */
 104    enum FMDRingBits {
 105      kFMD1I=0x11,
 106      kFMD1 =kFMD1I,
 107      kFMD2I=0x21,
 108      kFMD2O=0x22,
 109      kFMD2 =kFMD2I|kFMD2O,
 110      kFMD3I=0x31,
 111      kFMD3O=0x32,
 112      kFMD3 =kFMD2I|kFMD2O
 113   };
 114   /**
 115    * Constructor
 116    */
 117   AliFMDHistCollector();
 118   /**
 119    * Constructor
 120    *
 121    * @param title Name of object
 122    */
 123   AliFMDHistCollector(const char* title);
 124   /**
 125    * Copy constructor
 126    *
 127    * @param o Object to copy from
 128    */
 129   AliFMDHistCollector(const AliFMDHistCollector& o);
 130
 131   /**
 132    * Destructor
 133    */
 134   virtual ~AliFMDHistCollector();
 135   /**
 136    * Assignement operator
 137    *
 138    * @param o Object to assign from
 139    *
 140    * @return Reference to this object
 141    */
 142   AliFMDHistCollector& operator=(const AliFMDHistCollector&);
 143   /**
 144    * Intialise
 145    *
 146    * @param vtxAxis  @f$ v_z@f$ axis
 147    * @param etaAxis  @f$ \eta@f$ axis
 148    */
 149   virtual void Init(const TAxis& vtxAxis,
 150                     const TAxis& etaAxis);
 151   /**
 152    * Do the calculations
 153    *
 154    * @param hists    Cache of histograms
 155    * @param sums     Cache to sum ring histograms in
 156    * @param vtxBin   Vertex bin (1 based)
 157    * @param out      Output histogram
 158    *
 159    * @return true on successs
 160    */
 161   virtual Bool_t Collect(const AliForwardUtil::Histos& hists,
 162                          AliForwardUtil::Histos&       sums,
 163                          UShort_t                      vtxBin,
 164                          TH2D&                         out);
 165   /**
 166    * Output diagnostic histograms to directory
 167    *
 168    * @param dir List to write in
 169    */
 170   virtual void DefineOutput(TList* dir);
 171   /**
 172    * Set the merge method
 173    *
 174    * @param m Method
 175    */
 176   void SetMergeMethod(MergeMethod m) { fMergeMethod = m; }
 177   /**
 178    * Set the method for finding the fidicual area of the secondary maps
 179    *
 180    * @param m Method
 181    */
 182   void SetFiducialMethod(FiducialMethod m) { fFiducialMethod = m; }
 183   /**
 184    * Set the number of extra bins (beyond the secondary map border)
 185    * to cut away.
 186    *
 187    * @param n Number of bins
 188    */
 189   void SetNCutBins(UInt_t n=2) { fNCutBins = n; }
 190   /**
 191    * Set the correction cut, that is, when bins in the secondary
 192    * correction maps have a value lower than this cut, they are
 193    * considered uncertain and not used
 194    *
 195    * @param cut Cut-off
 196    */
 197   void SetCorrectionCut(Float_t cut=0.5) { fCorrectionCut = cut; }
 198   /**
 199    * Set FMD rings to skip. Argument should be
 200    * kFirstRingToSkip|kSecondRingToSkip...
 201    *
 202    * @param mask bit pattern
 203    */
 204   void SetFMDRingsToSkip(UShort_t mask) { fSkipFMDRings = mask; }
 205  /**
 206    * Set whether to make bg maps or not
 207    *
 208    * @param use make them
 209    */
 210   void SetMakeBGHitMaps(Bool_t use) { fBgAndHitMaps = use; }
 211
 212   /**
 213    * Set the debug level. The higher the value the more output
 214    *
 215    * @param dbg Debug level
 216    */
 217
 218   void SetDebug(Int_t dbg=1) { fDebug = dbg; }
 219   /**
 220    * Print information
 221    *
 222    * @param option Not used
 223    */
 224   void Print(Option_t* option="") const;
 225 protected:
 226   /**
 227    * Get the first and last eta bin to use for a given ring and vertex
 228    *
 229    * @param d        Detector
 230    * @param r        Ring
 231    * @param vtxBin   Vertex bin (1 based)
 232    * @param first    On return, the first eta bin to use
 233    * @param last     On return, the last eta bin to use
 234    */
 235   virtual void GetFirstAndLast(UShort_t d, Char_t r, UShort_t vtxBin,
 236                                Int_t& first, Int_t& last) const;
 237   /**
 238    * Get the first and last eta bin to use for a given ring and vertex
 239    *
 240    * @param idx      Ring index as given by GetIdx
 241    * @param vtxBin   Vertex bin (1 based)
 242    * @param first    On return, the first eta bin to use
 243    * @param last     On return, the last eta bin to use
 244    */
 245   virtual void GetFirstAndLast(Int_t idx, UShort_t vtxBin,
 246                                Int_t& first, Int_t& last) const;
 247   /**
 248    * Get the first eta bin to use for a given ring and vertex
 249    *
 250    * @param d Detector
 251    * @param r Ring
 252    * @param v vertex bin (1 based)
 253    *
 254    * @return First eta bin to use, or -1 in case of problems
 255    */
 256   Int_t GetFirst(UShort_t d, Char_t r, UShort_t v) const;
 257   /**
 258    * Get the first eta bin to use for a given ring and vertex
 259    *
 260    * @param idx Ring index as given by GetIdx
 261    * @param v vertex bin (1 based)
 262    *
 263    * @return First eta bin to use, or -1 in case of problems
 264    */
 265   Int_t GetFirst(Int_t idx, UShort_t v) const;
 266   /**
 267    * Get the last eta bin to use for a given ring and vertex
 268    *
 269    * @param d Detector
 270    * @param r Ring
 271    * @param v vertex bin (1 based)
 272    *
 273    * @return Last eta bin to use, or -1 in case of problems
 274    */
 275   Int_t GetLast(UShort_t d, Char_t r, UShort_t v) const;
 276   /**
 277    * Get the last eta bin to use for a given ring and vertex
 278    *
 279    * @param idx Ring index as given by GetIdx
 280    * @param v vertex bin (1 based)
 281    *
 282    * @return Last eta bin to use, or -1 in case of problems
 283    */
 284   Int_t GetLast(Int_t idx, UShort_t v) const;
 285   /**
 286    * Get the detector and ring from the ring index
 287    *
 288    * @param idx Ring index
 289    * @param d   On return, the detector or 0 in case of errors
 290    * @param r   On return, the ring id or '0' in case of errors
 291    */
 292   void GetDetRing(Int_t idx, UShort_t& d, Char_t& r) const;
 293   /**
 294    * Get the ring index from detector number and ring identifier
 295    *
 296    * @param d Detector
 297    * @param r Ring identifier
 298    *
 299    * @return ring index or -1 in case of problems
 300    */
 301   Int_t GetIdx(UShort_t d, Char_t r) const;
 302   /**
 303    * Get the possibly overlapping histogram of eta bin @a e in
 304    * detector and ring
 305    *
 306    * @param d Detector
 307    * @param r Ring
 308    * @param e Eta bin
 309    * @param v Vertex bin (1 based)
 310    *
 311    * @return Overlapping histogram index or -1
 312    */
 313   Int_t GetOverlap(UShort_t d, Char_t r, Int_t e, UShort_t v) const;
 314   /**
 315    * Get the possibly overlapping histogram of eta bin @a e in
 316    * detector and ring
 317    *
 318    * @param i Ring index
 319    * @param e Eta bin
 320    * @param v Vertex bin (1 based)
 321    *
 322    * @return Overlapping histogram index or -1
 323    */
 324   Int_t GetOverlap(Int_t i, Int_t e, UShort_t v) const;
 325   /**
 326    * Check if there's an overlapping histogram with this eta bin of
 327    * the detector and ring
 328    *
 329    * @param d Detector
 330    * @param r Ring
 331    * @param e eta bin
 332    * @param v Vertex bin (1 based)
 333    *
 334    * @return True if there's an overlapping histogram
 335    */
 336   Bool_t HasOverlap(UShort_t d, Char_t r, Int_t e, UShort_t v) const;
 337   /**
 338    * Check if there's an overlapping histogram with this eta bin of
 339    * ring
 340    *
 341    * @param i Ring index
 342    * @param e eta bin
 343    * @param v Vertex bin
 344    *
 345    * @return True if there's an overlapping histogram
 346    */
 347   Bool_t HasOverlap(Int_t i, Int_t e, UShort_t v) const;
 348   /**
 349    * Check if we should include the bin in the data range
 350    *
 351    * @param bg Secondary map histogram
 352    * @param ie Eta bin
 353    * @param ip Phi bin
 354    *
 355    * @return True if to be used
 356    */
 357   Bool_t CheckCorrection(const TH2D* bg, Int_t ie, Int_t ip) const;
 358   /**
 359    * Merge bins accoring to set method
 360    *
 361    * @param c   Current content
 362    * @param e   Current error
 363    * @param oc  Old content
 364    * @param oe  Old error
 365    * @param rc  On return, the new content
 366    * @param re  On return, tne new error
 367    */
 368   void MergeBins(Double_t c,   Double_t e,
 369                  Double_t oc,  Double_t oe,
 370                  Double_t& rc, Double_t& re) const;
 371
 372
 373   Int_t       fNCutBins;        // Number of additional bins to cut away
 374   Float_t     fCorrectionCut;   // Cut-off on secondary corrections
 375   TArrayI     fFirstBins;       // Array of first eta bins
 376   TArrayI     fLastBins;        // Array of last eta bins
 377   Int_t       fDebug;           // Debug level
 378   TList*      fList;            // Output list
 379   TH2D*       fSumRings;        // Sum per ring (on y-axis)
 380   TH2D*       fCoverage;        // Sum per ring (on y-axis)
 381   MergeMethod fMergeMethod;     // Merge methiod for overlapping bins
 382   FiducialMethod fFiducialMethod; // Fidicual method
 383   UShort_t    fSkipFMDRings;    // FMD rings to ignore
 384   Bool_t      fBgAndHitMaps;    // Make hit/bg maps or not
 385
 386   ClassDef(AliFMDHistCollector,3); // Calculate Nch density
 387 };
 388
 389 //____________________________________________________________________
 390 inline void
 391 AliFMDHistCollector::GetFirstAndLast(UShort_t d, Char_t r, UShort_t vtxbin,
 392                                      Int_t& first, Int_t& last) const
 393 {
 394   GetFirstAndLast(GetIdx(d,r), vtxbin, first, last);
 395 }
 396 //____________________________________________________________________
 397 inline Int_t
 398 AliFMDHistCollector::GetFirst(UShort_t d, Char_t r, UShort_t v) const
 399 {
 400   return GetFirst(GetIdx(d,r), v);
 401 }
 402 //____________________________________________________________________
 403 inline Int_t
 404 AliFMDHistCollector::GetLast(UShort_t d, Char_t r, UShort_t v) const
 405 {
 406   return GetLast(GetIdx(d, r), v);
 407 }
 408 //____________________________________________________________________
 409 inline Bool_t
 410 AliFMDHistCollector::HasOverlap(UShort_t d, Char_t r, Int_t e, UShort_t v) const
 411 {
 412   return GetOverlap(d,r,e,v) >= 0;
 413 }
 414 //____________________________________________________________________
 415 inline Bool_t
 416 AliFMDHistCollector::HasOverlap(Int_t i, Int_t e, UShort_t v) const
 417 {
 418   return GetOverlap(i,e,v) >= 0;
 419 }
 420
 421 #endif
 422 // Local Variables:
 423 //   mode: C++
 424 // End:
 425