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