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 TH2;
  23 class TH2D;
  24 class TH1D;
  25 class TObjArray;
  26
  27 /**
  28  * This class collects the event histograms into single histograms,
  29  * one for each ring in each vertex bin.
  30  *
  31  * @par Input:
  32  *   - AliESDFMD object possibly corrected for sharing
  33  *
  34  * @par Output:
  35  *   - 5 RingHistos objects - each with a number of vertex dependent
  36  *     2D histograms of the inclusive charge particle density
  37  *
  38  * @par HistCollector used:
  39  *   - AliFMDCorrSecondaryMap
  40  *
  41  * @ingroup pwglf_forward_algo
  42  * @ingroup pwglf_forward_aod
  43  */
  44 class AliFMDHistCollector : public TNamed
  45 {
  46 public:
  47   /**
  48    * Methods to use when merging overlapping bins @f$b_1@f$, @f$b_2@f$
  49    * with content @f$c_1@f$, @f$c_2@f$, and errors @f$e_1@f$,
  50    * @f$e_2@f$ into bin @f$ b@f$ with content @f$c@f$ and error @f$e@f$
  51    */
  52   enum MergeMethod {
  53     /**
  54      * @f[
  55      *   c = \frac{1}{2}(c_1+c_2)
  56      * @f]
  57      * @f[
  58      *   e = \sqrt{e_1^2+e_2^2}
  59      * @f]
  60      */
  61     kStraightMean,
  62     /**
  63      * As above, exept zero's are ignored
  64      */
  65     kStraightMeanNoZero,
  66     /**
  67      * @f[
  68      *   c = \frac{\frac{c_1}{e_1^2}+\frac{c_2}{e_2^2}}{
  69      *             \frac{1}{e_1^2}+\frac{1}{e_2^2}}
  70      * @f]
  71      * @f[
  72      *   e = \sqrt{\frac{1}{\frac{1}{e_1^2}+\frac{1}{e_2^2}}}
  73      * @f]
  74      */
  75     kWeightedMean,
  76     /**
  77      * @f[
  78      *     c = \left\{\begin{array}{cl}
  79      *          c_1 & \mbox{if $e_1 < e_2$} \\
  80      *          c_2 & \mbox{otherwise}\end{array}\right.
  81      * @f]
  82      */
  83     kLeastError,
  84     /**
  85      * Just sum the signals
  86      */
  87     kSum,
  88     /**
  89      * In overlaps, prefer inners, or if both are inners,
  90      * do the straight mean
  91      */
  92     kPreferInner,
  93     /**
  94      * In overlaps, prefer outers, or if both are outers (doesn't happen),
  95      * do the straight mean
  96      */
  97     kPreferOuter
  98
  99   };
 100   /**
 101    * How to obtain the fiducial cuts
 102    */
 103   enum FiducialMethod {
 104     /**
 105      * Select bins by fixed cut.  Bins with a secondary correction
 106      * less than the cut is considered as non-valid
 107      */
 108     kByCut,
 109     /**
 110      * A bin is considered non-valid, if it is less then twice as
 111      * large as it's neighbors (in eta)
 112      */
 113     kDistance
 114   };
 115   /**
 116    * FMD ring bits for skipping
 117    */
 118    enum FMDRingBits {
 119      kFMD1I=0x01,
 120      kFMD1 =kFMD1I,
 121      kFMD2I=0x02,
 122      kFMD2O=0x04,
 123      kFMD2 =kFMD2I|kFMD2O,
 124      kFMD3I=0x08,
 125      kFMD3O=0x10,
 126      kFMD3 =kFMD3I|kFMD3O
 127   };
 128   /**
 129    * Constructor
 130    */
 131   AliFMDHistCollector();
 132   /**
 133    * Constructor
 134    *
 135    * @param title Name of object
 136    */
 137   AliFMDHistCollector(const char* title);
 138   /**
 139    * Copy constructor
 140    *
 141    * @param o Object to copy from
 142    */
 143   AliFMDHistCollector(const AliFMDHistCollector& o);
 144
 145   /**
 146    * Destructor
 147    */
 148   virtual ~AliFMDHistCollector();
 149   /**
 150    * Assignement operator
 151    *
 152    * @param o Object to assign from
 153    *
 154    * @return Reference to this object
 155    */
 156   AliFMDHistCollector& operator=(const AliFMDHistCollector&);
 157   /**
 158    * Intialise
 159    *
 160    * @param vtxAxis  @f$ v_z@f$ axis
 161    * @param etaAxis  @f$ \eta@f$ axis
 162    */
 163   virtual void SetupForData(const TAxis& vtxAxis,
 164                     const TAxis& etaAxis);
 165   /**
 166    * Do the calculations
 167    *
 168    * @param hists    Cache of histograms
 169    * @param sums     Cache to sum ring histograms in
 170    * @param vtxBin   Vertex bin (1 based)
 171    * @param out      Output histogram
 172    * @param cent     Centrality
 173    * @param eta2phi    Copy eta coverage to phi acceptance
 174    *
 175    * @return true on successs
 176    */
 177   virtual Bool_t Collect(const AliForwardUtil::Histos& hists,
 178                          AliForwardUtil::Histos&       sums,
 179                          UShort_t                      vtxBin,
 180                          TH2D&                         out,
 181                          Double_t                      cent=-1.0,
 182                          Bool_t                        eta2phi=false);
 183   /**
 184    * Output diagnostic histograms to directory
 185    *
 186    * @param dir List to write in
 187    */
 188   virtual void CreateOutputObjects(TList* dir);
 189   /**
 190    * Set the merge method
 191    *
 192    * @param m Method
 193    */
 194   void SetMergeMethod(MergeMethod m) { fMergeMethod = m; }
 195   MergeMethod GetMergeMethod() const { return fMergeMethod; }
 196   /**
 197    * Set the method for finding the fidicual area of the secondary maps
 198    *
 199    * @param m Method
 200    */
 201   void SetFiducialMethod(FiducialMethod m) { fFiducialMethod = m; }
 202   /**
 203    * Set the number of extra bins (beyond the secondary map border)
 204    * to cut away.
 205    *
 206    * @param n Number of bins
 207    */
 208   void SetNCutBins(UInt_t n=2) { fNCutBins = n; }
 209   /**
 210    * Set the correction cut, that is, when bins in the secondary
 211    * correction maps have a value lower than this cut, they are
 212    * considered uncertain and not used
 213    *
 214    * @param cut Cut-off
 215    */
 216   void SetCorrectionCut(Float_t cut=0.5) { fCorrectionCut = cut; }
 217   /**
 218    * Set FMD rings to skip. Argument should be
 219    * kFirstRingToSkip|kSecondRingToSkip...
 220    *
 221    * @param mask bit pattern
 222    */
 223   void SetFMDRingsToSkip(UShort_t mask) { fSkipFMDRings = mask; }
 224  /**
 225    * Set whether to make bg maps or not
 226    *
 227    * @param use make them
 228    */
 229   void SetMakeBGHitMaps(Bool_t use) { fBgAndHitMaps = use; }
 230   /**
 231    * Set whether to make by-centrality sums for each ring
 232    *
 233    * @param use If true, make by-centrality sums
 234    */
 235   void SetMakeCentralitySums(Bool_t use) { fDoByCent = use; }
 236   /**
 237    * Set the debug level. The higher the value the more output
 238    *
 239    * @param dbg Debug level
 240    */
 241   void SetDebug(Int_t dbg=1) { fDebug = dbg; }
 242   /**
 243    * Print information
 244    *
 245    * @param option Not used
 246    */
 247   void Print(Option_t* option="") const;
 248 protected:
 249   /**
 250    * Get the detector and ring from the ring index
 251    *
 252    * @param idx Ring index
 253    * @param d   On return, the detector or 0 in case of errors
 254    * @param r   On return, the ring id or '0' in case of errors
 255    */
 256   static void GetDetRing(Int_t idx, UShort_t& d, Char_t& r);
 257   /**
 258    * Get the ring index from detector number and ring identifier
 259    *
 260    * @param d Detector
 261    * @param r Ring identifier
 262    *
 263    * @return ring index or -1 in case of problems
 264    */
 265   static Int_t GetIdx(UShort_t d, Char_t r);
 266   /**
 267    * Check if the detector @a d, ring @a r is listed <i>in</i> the @a
 268    * skips bit mask.  If the detector/ring is in the mask, return true.
 269    *
 270    * That is, use case is
 271    * @code
 272    *  for (UShort_t d=1. d<=3, d++) {
 273    *    UShort_t nr = (d == 1 ? 1 : 2);
 274    *    for (UShort_t q = 0; q < nr; q++) {
 275    *      Char_t r = (q == 0 ? 'I' : 'O');
 276    *      if (CheckSkips(d, r, skips)) continue;
 277    *      // Process detector/ring
 278    *    }
 279    *  }
 280    * @endcode
 281    *
 282    * @param d      Detector
 283    * @param r      Ring
 284    * @param skips  Mask of detector/rings to skip
 285    *
 286    * @return True if detector @a d, ring @a r is in the mask @a skips
 287    */
 288   static Bool_t CheckSkip(UShort_t d, Char_t r, UShort_t skips);
 289   /**
 290    * Check the correction
 291    *
 292    * @param m   Fiducial method used
 293    * @param cut Cut value
 294    * @param bg  Secondary map
 295    * @param ie  @f$\eta@f$ bin
 296    * @param ip  @f$\varphi@f$ bin
 297    *
 298    * @return true if OK.
 299    */
 300   static Bool_t CheckCorrection(FiducialMethod m, Double_t cut,
 301                                 const TH2D* bg, Int_t ie, Int_t ip);
 302
 303   /**
 304    * Merge bins accoring to set method
 305    *
 306    * @param m   Merging method
 307    * @param c   Current content
 308    * @param e   Current error
 309    * @param oc  Old content
 310    * @param oe  Old error
 311    * @param rc  On return, the new content
 312    * @param re  On return, tne new error
 313    */
 314   static void MergeBins(MergeMethod   m,
 315                         Double_t c,   Double_t e,
 316                         Double_t oc,  Double_t oe,
 317                         Double_t& rc, Double_t& re);
 318
 319   //==================================================================
 320   /**
 321    * Structure to hold per-vertex bin cache of per-ring histograms
 322    */
 323   struct VtxBin : public TObject
 324   {
 325     /**
 326      * Constructor
 327      *
 328      * @param index   Index number
 329      * @param minIpZ  Least @f$IP_{z}@f$
 330      * @param maxIpZ  Largest @f$IP_{z}@f$
 331      * @param nCut    Cut on n
 332      */
 333     VtxBin(Int_t index=0, Double_t minIpZ=999, Double_t maxIpZ=-999,
 334            Int_t nCut=0);
 335     /**
 336      * Copy constructor
 337      *
 338      * @param o Object to copy from
 339      */
 340     VtxBin(const VtxBin& o);
 341     /**
 342      * Assignment operator
 343      *
 344      * @param o Object to assign from
 345      *
 346      * @return Reference to this object
 347      */
 348     VtxBin& operator=(const VtxBin& o);
 349     /**
 350      * Override to give name based on cuts
 351      *
 352      * @return Name
 353      */
 354     const Char_t* GetName() const;
 355     /**
 356      * Set up for data
 357      *
 358      * @param coverage    Diagnostics histogram to be filled
 359      * @param skip        Skip flags
 360      * @param fiducial    Fiducial cut method
 361      * @param cut         Fiducial cut
 362      * @param l           Parent output list
 363      * @param etaAxis     @f$\eta@f$ axis used
 364      * @param doHitMap    If true, also do a per-ring sum
 365      * @param storeSecMap If true, store used secondary map
 366      */
 367     void SetupForData(TH2*           coverage,
 368                       UShort_t       skip,
 369                       FiducialMethod fiducial,
 370                       Double_t       cut,
 371                       TList*         l,
 372                       const TAxis&   etaAxis,
 373                       Bool_t         doHitMap,
 374                       Bool_t         storeSecMap);
 375     /**
 376      * Process one event in this vertex bin
 377      *
 378      * @param hists      Histograms
 379      * @param sums       Sum histograms
 380      * @param out        Per-event output histogram
 381      * @param sumRings   Sum per ring
 382      * @param skipped    Histogram of skipped rings
 383      * @param cent       Event centrality
 384      * @param m          Merging method
 385      * @param skips      Which rings to skip
 386      * @param byCent     List (or null) of per centrality sums
 387      * @param eta2phi    Copy eta coverage to phi acceptance
 388      *
 389      * @return true on success
 390      */
 391     Bool_t Collect(const AliForwardUtil::Histos& hists,
 392                    AliForwardUtil::Histos&       sums,
 393                    TH2D&                         out,
 394                    TH2D*                         sumRings,
 395                    TH1D*                         skipped,
 396                    Double_t                      cent,
 397                    MergeMethod                   m,
 398                    UShort_t                      skips,
 399                    TList*                        byCent,
 400                    Bool_t                        eta2phi);
 401     /**
 402      * Check if there's an overlap between detector @a d, ring @a r
 403      * and some other ring for the given @f$\eta@f$ @a bin.  If so,
 404      * return the ring index.  If not, return -1.
 405      *
 406      * @param d    Current detector
 407      * @param r    Current ring
 408      * @param bin  Current @f$\eta@f$ bin
 409      *
 410      * @return Index of overlapping ring, or -1
 411      */
 412     Int_t GetOverlap(UShort_t d, Char_t r, Int_t bin) const;
 413     /**
 414      * Get the first and last @f$\eta@f$ bin for a detector
 415      *
 416      * @param d      Current detector
 417      * @param r      Current ring
 418      * @param first  On return, the first @f$\eta@f$ bin
 419      * @param last   On return, the last @f$\eta@f$ bin
 420      */
 421     void GetFirstAndLast(UShort_t d, UShort_t r,
 422                          Int_t& first, Int_t& last) const {
 423       GetFirstAndLast(GetIdx(d,r), first, last);
 424     }
 425     /**
 426      * Get the first and last @f$\eta@f$ bin for a detector
 427      *
 428      * @param idx    Current ring index
 429      * @param first  On return, the first @f$\eta@f$ bin
 430      * @param last   On return, the last @f$\eta@f$ bin
 431      */
 432     void GetFirstAndLast(Int_t idx,Int_t& first, Int_t& last) const;
 433     /**
 434      * Get the first @f$\eta@f$ bin
 435      *
 436      * @param idx Ring index (0-4)
 437      *
 438      * @return bin number
 439      */
 440     Int_t GetFirst(Int_t idx) const;
 441     /**
 442      * Get the last @f$\eta@f$ bin
 443      *
 444      * @param idx Ring index (0-4)
 445      *
 446      * @return bin number
 447      */
 448     Int_t GetLast(Int_t idx) const;
 449     /**
 450      * Get the first @f$\eta@f$ bin
 451      *
 452      * @param d  Detector
 453      * @param r  Ring
 454      *
 455      * @return bin number
 456      */
 457     Int_t GetFirst(UShort_t d, Char_t r) const { return GetFirst(GetIdx(d,r));}
 458     /**
 459      * Get the last @f$\eta@f$ bin
 460      *
 461      * @param d  Detector
 462      * @param r  Ring
 463      *
 464      * @return bin number
 465      */
 466     Int_t GetLast(UShort_t d, Char_t r) const { return GetLast(GetIdx(d,r));}
 467
 468     Int_t                   fIndex;     // Vertex bin index
 469     Double_t                fLow;       // Low @f$ ip_z @f$
 470     Double_t                fHigh;      // High @f$ ip_z @f$
 471     AliForwardUtil::Histos* fHitMap;    // Hit map (optional)
 472     TArrayI                 fFirstBin;  // Per-ring first bin
 473     TArrayI                 fLastBin;   // Per-ring last bin
 474     Int_t                   fNCutBins;  // Number of bins to cut
 475
 476     // ClassDef(VtxBin,1); // Vertex bin in histogram collector
 477   };
 478   /**
 479    * Get a vertex bin
 480    *
 481    * @param ivtx Bin number (1-nVz)
 482    *
 483    * @return Bin or null
 484    */
 485   VtxBin* GetVtxBin(Int_t ivtx);
 486   /**
 487    * Get a vertex bin
 488    *
 489    * @param ivtx Bin number (1-nVz)
 490    *
 491    * @return Bin or null
 492    */
 493   const VtxBin* GetVtxBin(Int_t ivtx) const;
 494
 495   Int_t       fNCutBins;        // Number of additional bins to cut away
 496   Float_t     fCorrectionCut;   // Cut-off on secondary corrections
 497   Int_t       fDebug;           // Debug level
 498   TList*      fList;            // Output list
 499   TH2D*       fSumRings;        // Sum per ring (on y-axis)
 500   TH2D*       fCoverage;        // Sum per ring (on y-axis)
 501   TH1D*       fSkipped;         // Skipped rings
 502   MergeMethod fMergeMethod;     // Merge methiod for overlapping bins
 503   FiducialMethod fFiducialMethod; // Fidicual method
 504   UShort_t    fSkipFMDRings;    // FMD rings to ignore
 505   Bool_t      fBgAndHitMaps;    // Make hit/bg maps or not
 506   TObjArray*  fVtxList;         //! Per-vertex list
 507   TList*      fByCent;          // By centrality sums
 508   Bool_t      fDoByCent;        // Whether to do by centrality sum
 509   ClassDef(AliFMDHistCollector,6); // Calculate Nch density
 510 };
 511
 512
 513 #endif
 514 // Local Variables:
 515 //   mode: C++
 516 // End:
 517