PWGLF/FORWARD/analysis2/AliFMDSharingFilter.h

   1 //
   2 // Class to do the sharing correction.  That is, a filter that merges
   3 // adjacent strip signals presumably originating from a single particle
   4 // that impinges on the detector in such a way that it deposite energy
   5 // into two or more strips.
   6 //
   7 #ifndef ALIFMDSHARINGFILTER_H
   8 #define ALIFMDSHARINGFILTER_H
   9 /**
  10  * @file   AliFMDSharingFilter.h
  11  * @author Christian Holm Christensen <cholm@dalsgaard.hehi.nbi.dk>
  12  * @date   Wed Mar 23 14:03:57 2011
  13  *
  14  * @brief
  15  *
  16  *
  17  * @ingroup pwglf_forward_aod
  18  */
  19 #include <TNamed.h>
  20 #include <TH2.h>
  21 #include <TList.h>
  22 #include "AliForwardUtil.h"
  23 #include "AliFMDMultCuts.h"
  24 #include <TBits.h>
  25 class AliESDFMD;
  26 class TAxis;
  27 class TList;
  28 class TH2;
  29 class AliFMDFloatMap;
  30
  31 /**
  32  * Class to do the sharing correction.  That is, a filter that merges
  33  * adjacent strip signals presumably originating from a single particle
  34  * that impinges on the detector in such a way that it deposite energy
  35  * into two or more strips.
  36  *
  37  * @image html alice-int-2012-040-share_fraction.png "Energy loss sharing"
  38  *
  39  * @par Input:
  40  *    - AliESDFMD object  - from reconstruction
  41  *
  42  * @par Output:
  43  *    - AliESDFMD object  - copy of input, but with signals merged
  44  *
  45  * @par Corrections used:
  46  *    - AliFMDCorrELossFit
  47  *
  48  * @par Histograms:
  49  *    - For each ring (FMD1i, FMD2i, FMD2o, FMD3i, FMD3o) the distribution of
  50  *      signals before and after the filter.
  51  *    - For each ring (see above), an array of distributions of number of
  52  *      hit strips for each vertex bin (if enabled - see SetupForData method)
  53  *
  54  *
  55  * @ingroup pwglf_forward_algo
  56  * @ingroup pwglf_forward_aod
  57  */
  58 class AliFMDSharingFilter : public TNamed
  59 {
  60 public:
  61   /**
  62    * Status of a strip
  63    * @deprecated Not used
  64    */
  65   enum Status {
  66     /** Nothing yet */
  67     kNone             = 1,
  68     /** Candidate for merging */
  69     kCandidate        = 2,
  70     /** This was merged into other strip */
  71     kMergedWithOther  = 3,
  72     /** Other strips was merged into this */
  73     kMergedInto       = 4
  74   };
  75   /**
  76    * Destructor
  77    */
  78   virtual ~AliFMDSharingFilter();
  79   /**
  80    * Default Constructor - do not use
  81    */
  82   AliFMDSharingFilter();
  83   /**
  84    * Constructor
  85    *
  86    * @param title Title of object  - not significant
  87    */
  88   AliFMDSharingFilter(const char* title);
  89   /**
  90    * Copy constructor
  91    *
  92    * @param o Object to copy from
  93    */
  94   AliFMDSharingFilter(const AliFMDSharingFilter& o);
  95   /**
  96    * Assignment operator
  97    *
  98    * @param o Object to assign from
  99    *
 100    * @return Reference to this
 101    */
 102   AliFMDSharingFilter& operator=(const AliFMDSharingFilter& o);
 103
 104   /**
 105    * Initialize
 106    *
 107    * @param axis Default eta axis from parent task
 108    */
 109   void SetupForData(const TAxis& axis);
 110   /**
 111    * Set the debug level.  The higher the value the more output
 112    *
 113    * @param dbg Debug level
 114    */
 115   virtual void SetDebug(Int_t dbg=1) { fDebug = dbg; }
 116
 117   /**
 118    * Enable use of angle corrected signals in the algorithm
 119    *
 120    * @param use If true, use angle corrected signals,
 121    * otherwise use de-corrected signals.  In the final output, the
 122    * signals are always angle corrected.
 123    */
 124   void SetUseAngleCorrectedSignals(Bool_t use) { fCorrectAngles = use; }
 125    /**
 126    * Enable zeroing of signals if below high cut
 127    *
 128    * @param use zero the signals if below sharing cut
 129    *
 130    */
 131   void SetZeroSharedHitsBelowThreshold(Bool_t use) { fZeroSharedHitsBelowThreshold = use; }
 132  /**
 133    * Enable a simpler merging algorithm
 134    *
 135    * @param use use the simpler algorithm
 136    *
 137    */
 138   void SetUseSimpleSharing(Bool_t use) { fUseSimpleMerging = use; }
 139   /**
 140    * In case of a simpler merging algorithm allow 3 strips to be
 141    * merged
 142    *
 143    * @param use allow three strips
 144    *
 145    */
 146   void SetAllow3Strips(Bool_t use) { fThreeStripSharing = use; }
 147
 148   /**
 149    * In case of a displaced vertices recalculate eta and angle correction
 150    *
 151    * @param use recalculate or not
 152    *
 153    */
 154   void SetRecalculateEta(Bool_t use) { fRecalculateEta = use; }
 155   /**
 156    * Set whether to consider invalid multiplicities as null (or empty)
 157    * signal.
 158    *
 159    * @param flag If true, count invalids as empty
 160    */
 161   void SetInvalidIsEmpty(Bool_t flag) { fInvalidIsEmpty = flag; }
 162
 163   /**
 164    * Filter the input AliESDFMD object
 165    *
 166    * @param input     Input
 167    * @param lowFlux   If this is a low-flux event
 168    * @param output    Output AliESDFMD object
 169    * @param zvtx      Vertex position
 170    *
 171    * @return True on success, false otherwise
 172    */
 173   Bool_t Filter(const AliESDFMD& input,
 174                 Bool_t           lowFlux,
 175                 AliESDFMD&       output,
 176                 Double_t         zvtx);
 177   /**
 178    * Scale the histograms to the total number of events
 179    *
 180    * @param dir     Where the output is
 181    * @param output  Output list
 182    * @param nEvents Number of events
 183    */
 184   virtual void Terminate(const TList* dir, TList* output, Int_t nEvents);
 185
 186   /**
 187    * Define the output histograms.  These are put in a sub list of the
 188    * passed list.   The histograms are merged before the parent task calls
 189    * AliAnalysisTaskSE::Terminate
 190    *
 191    * @param dir Directory to add to
 192    */
 193   virtual void CreateOutputObjects(TList* dir);
 194   /**
 195    * Print information
 196    *
 197    * @param option Not used
 198    */
 199   virtual void Print(Option_t* option="") const;
 200
 201   /**
 202    * Get the low cuts
 203    *
 204    * @return Reference to low cuts
 205    */
 206   AliFMDMultCuts& GetLCuts() { return fLCuts; }
 207   /**
 208    * Get the high cuts
 209    *
 210    * @return Reference to high cuts
 211    */
 212   AliFMDMultCuts& GetHCuts() { return fHCuts; }
 213   /**
 214    * Get the low cuts
 215    *
 216    * @return Reference to low cuts
 217    */
 218   const AliFMDMultCuts& GetLCuts() const { return fLCuts; }
 219   /**
 220    * Get the high cuts
 221    *
 222    * @return Reference to high cuts
 223    */
 224   const AliFMDMultCuts& GetHCuts() const { return fHCuts; }
 225   /**
 226    * Set the low cuts
 227    *
 228    * @param c Cuts object
 229    */
 230   void SetLCuts(const AliFMDMultCuts& c) { fLCuts = c; }
 231   /**
 232    * Set the high cuts
 233    *
 234    * @param c Cuts object
 235    */
 236   void SetHCuts(const AliFMDMultCuts& c) { fHCuts = c; }
 237   /**
 238    * Add a dead strip
 239    *
 240    * @param d  Detector
 241    * @param r  Ring
 242    * @param s  Sector
 243    * @param t  Strip
 244    */
 245   void AddDead(UShort_t d, Char_t r, UShort_t s, UShort_t t);
 246   /**
 247    * Add a dead region in a detector ring
 248    *
 249    * @param d   Detector
 250    * @param r   Ring
 251    * @param s1  First sector (inclusive)
 252    * @param s2  Last sector (inclusive)
 253    * @param t1  First strip (inclusive)
 254    * @param t2  Last strip (inclusive)
 255    */
 256   void AddDeadRegion(UShort_t d, Char_t r, UShort_t s1, UShort_t s2,
 257                      UShort_t t1, UShort_t t2);
 258   /**
 259    * Add dead strips from a script.  The script is supposed to accept
 260    * a pointer to this object (AliFMDSharingFilter) and then call
 261    * AddDead or AddDeadRegion as needed.
 262    *
 263    * @code
 264    * void deadstrips(AliFMDSharingFilter* filter)
 265    * {
 266    *   filter->AddDead(...);
 267    *   // ... and so on
 268    * }
 269    * @endcode
 270    *
 271    * @param script The script to read dead strips from.
 272    */
 273   void AddDead(const Char_t* script);
 274 protected:
 275   /**
 276    * Internal data structure to keep track of the histograms
 277    */
 278   struct RingHistos : public AliForwardUtil::RingHistos
 279   {
 280     /**
 281      * Default CTOR
 282      */
 283     RingHistos();
 284     /**
 285      * Constructor
 286      *
 287      * @param d detector
 288      * @param r ring
 289      */
 290     RingHistos(UShort_t d, Char_t r);
 291     /**
 292      * Copy constructor
 293      *
 294      * @param o Object to copy from
 295      */
 296     RingHistos(const RingHistos& o);
 297     /**
 298      * Assignment operator
 299      *
 300      * @param o Object to assign from
 301      *
 302      * @return Reference to this
 303      */
 304     RingHistos& operator=(const RingHistos& o);
 305     /**
 306      * Destructor
 307      */
 308     ~RingHistos();
 309     /**
 310      * Clear this object
 311      */
 312     // void Clear(const Option_t* ="") { fNHits = 0; }
 313     /**
 314      * Increase number of hits
 315      *
 316      */
 317     // void Incr() { fNHits++; }
 318     /**
 319      * Finish off
 320      *
 321      */
 322     // void Finish();
 323     /**
 324      * Make output
 325      *
 326      * @param dir where to store
 327      */
 328     void CreateOutputObjects(TList* dir);
 329     /**
 330      * Scale the histograms to the total number of events
 331      *
 332      * @param nEvents Number of events
 333      * @param dir     Where the output is
 334      */
 335     void Terminate(const TList* dir, Int_t nEvents);
 336     TH1D*     fBefore;       // Distribution of signals before filter
 337     TH1D*     fAfter;        // Distribution of signals after filter
 338     TH1D*     fSingle;       // Distribution of 1 signal after filter
 339     TH1D*     fDouble;       // Distribution of 2 signals after filter
 340     TH1D*     fTriple;       // Distribution of 3 signals after filter
 341     TH2D*     fSinglePerStrip;       // Distribution of 1 signal per strip
 342     // TH1D*     fDistanceBefore; //Distance between signals before sharing
 343     // TH1D*     fDistanceAfter; //Distance between signals after sharing
 344     TH2D*     fBeforeAfter;  // Correlation of before and after
 345     TH2D*     fNeighborsBefore; // Correlation of neighbors
 346     TH2D*     fNeighborsAfter; // Correlation of neighbors
 347     TH2D*     fSum;          // Summed signal
 348     // TH1D*     fHits;         // Distribution of hit strips.
 349     // Int_t     fNHits;        // Number of hit strips per event
 350     ClassDef(RingHistos,3);
 351   };
 352   /**
 353    * Get the ring histogram container
 354    *
 355    * @param d Detector
 356    * @param r Ring
 357    *
 358    * @return Ring histogram container
 359    */
 360   RingHistos* GetRingHistos(UShort_t d, Char_t r) const;
 361   /**
 362    * Get the signal in a strip
 363    *
 364    * @param fmd   ESD object
 365    * @param d     Detector
 366    * @param r     Ring
 367    * @param s     Sector
 368    * @param t     Strip
 369    *
 370    * @return The energy signal
 371    */
 372   Double_t SignalInStrip(const AliESDFMD& fmd,
 373                          UShort_t d,
 374                          Char_t   r,
 375                          UShort_t s,
 376                          UShort_t t) const;
 377   /**
 378    * Angle correct the signal
 379    *
 380    * @param mult Angle Un-corrected Signal
 381    * @param eta  Pseudo-rapidity
 382    *
 383    * @return Angle corrected signal
 384    */
 385   Double_t AngleCorrect(Double_t mult, Double_t eta) const;
 386   /**
 387    * Angle de-correct the signal
 388    *
 389    * @param mult Angle corrected Signal
 390    * @param eta  Pseudo-rapidity
 391    *
 392    * @return Angle un-corrected signal
 393    */
 394   Double_t DeAngleCorrect(Double_t mult, Double_t eta) const;
 395   /**
 396    * Get the high cut.  The high cut is defined as the
 397    * most-probably-value peak found from the energy distributions, minus
 398    * 2 times the width of the corresponding Landau.
 399    *
 400    * @param d      Detector
 401    * @param r      Ring
 402    * @param eta    Eta value
 403    * @param errors If false, do not show errors
 404    *
 405    * @return 0 or less on failure, otherwise @f$\Delta_{mp}-n\xi@f$
 406    */
 407   virtual Double_t GetHighCut(UShort_t d, Char_t r, Double_t eta,
 408                               Bool_t errors=true) const;
 409   /**
 410    * Get the low cut.  Normally, the low cut is taken to be the lower
 411    * value of the fit range used when generating the energy loss fits.
 412    * However, if fLowCut is set (using SetLowCit) to a value greater
 413    * than 0, then that value is used.
 414    *
 415    * @param d    Detector
 416    * @param r    Ring
 417    * @param eta  Eta value
 418    *
 419    * @return
 420    */
 421   virtual Double_t GetLowCut(UShort_t d, Char_t r, Double_t eta) const;
 422
 423   virtual Bool_t IsDead(UShort_t d, Char_t r, UShort_t s, UShort_t t) const;
 424   TList    fRingHistos;    // List of histogram containers
 425   // Double_t fLowCut;        // Low cut on sharing
 426   Bool_t   fCorrectAngles; // Whether to work on angle corrected signals
 427   // TH2*     fSummed;        // Operations histogram
 428   TH2*     fHighCuts;      // High cuts used
 429   TH2*     fLowCuts;       // High cuts used
 430   // AliFMDFloatMap* fOper;   // Operation done per strip
 431   Int_t    fDebug;         // Debug level
 432   Bool_t   fZeroSharedHitsBelowThreshold; //Whether to zero shared strip below cut
 433   AliFMDMultCuts fLCuts;    //Cuts object for low cuts
 434   AliFMDMultCuts fHCuts;    //Cuts object for high cuts
 435   Bool_t   fUseSimpleMerging; //enable simple sharing by HHD
 436   Bool_t   fThreeStripSharing; //In case of simple sharing allow 3 strips
 437   Bool_t   fRecalculateEta; //Whether to recalculate eta and angle correction (disp vtx)
 438   // TArrayI  fExtraDead;      // List of extra dead channels
 439   TBits    fXtraDead;
 440   Bool_t   fInvalidIsEmpty;  // Consider kInvalidMult as zero
 441
 442   ClassDef(AliFMDSharingFilter,8); //
 443 };
 444
 445 #endif
 446 // Local Variables:
 447 //  mode: C++
 448 // End: