PWG2/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 pwg2_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 pwg2_forward_algo
  53  * @ingroup pwg2_forward_aod
  54  */
  55 class AliFMDSharingFilter : public TNamed
  56 {
  57 public:
  58   enum Status {
  59     kNone             = 1,
  60     kCandidate        = 2,
  61     kMergedWithOther  = 3,
  62     kMergedInto       = 4
  63   };
  64   /**
  65    * Destructor
  66    */
  67   virtual ~AliFMDSharingFilter();
  68   /**
  69    * Default Constructor - do not use
  70    */
  71   AliFMDSharingFilter();
  72   /**
  73    * Constructor
  74    *
  75    * @param title Title of object  - not significant
  76    */
  77   AliFMDSharingFilter(const char* title);
  78   /**
  79    * Copy constructor
  80    *
  81    * @param o Object to copy from
  82    */
  83   AliFMDSharingFilter(const AliFMDSharingFilter& o);
  84   /**
  85    * Assignment operator
  86    *
  87    * @param o Object to assign from
  88    *
  89    * @return Reference to this
  90    */
  91   AliFMDSharingFilter& operator=(const AliFMDSharingFilter& o);
  92
  93   /**
  94    * Initialize
  95    *
  96    */
  97   void Init();
  98   /**
  99    * Set the debug level.  The higher the value the more output
 100    *
 101    * @param dbg Debug level
 102    */
 103   void SetDebug(Int_t dbg=1) { fDebug = dbg; }
 104
 105   /**
 106    * Enable use of angle corrected signals in the algorithm
 107    *
 108    * @param use If true, use angle corrected signals,
 109    * otherwise use de-corrected signals.  In the final output, the
 110    * signals are always angle corrected.
 111    */
 112   void SetUseAngleCorrectedSignals(Bool_t use) { fCorrectAngles = use; }
 113    /**
 114    * Enable zeroing of signals if below high cut
 115    *
 116    * @param use zero the signals if below sharing cut
 117    *
 118    */
 119   void SetZeroSharedHitsBelowThreshold(Bool_t use) { fZeroSharedHitsBelowThreshold = use; }
 120  /**
 121    * Enable a simpler merging algorithm
 122    *
 123    * @param use use the simpler algorithm
 124    *
 125    */
 126   void SetUseSimpleSharing(Bool_t use) { fUseSimpleMerging = use; }
 127   /**
 128    * In case of a simpler merging algorithm allow 3 strips to be
 129    * merged
 130    *
 131    * @param use allow three strips
 132    *
 133    */
 134   void SetAllow3Strips(Bool_t use) { fThreeStripSharing = use; }
 135   /**
 136    * Filter the input AliESDFMD object
 137    *
 138    * @param input     Input
 139    * @param lowFlux   If this is a low-flux event
 140    * @param output    Output AliESDFMD object
 141    *
 142    * @return True on success, false otherwise
 143    */
 144   Bool_t Filter(const AliESDFMD& input,
 145                 Bool_t           lowFlux,
 146                 AliESDFMD&       output);
 147   /**
 148    * Scale the histograms to the total number of events
 149    *
 150    * @param dir     Where the output is
 151    * @param nEvents Number of events
 152    */
 153   virtual void ScaleHistograms(const TList* dir, Int_t nEvents);
 154
 155   /**
 156    * Define the output histograms.  These are put in a sub list of the
 157    * passed list.   The histograms are merged before the parent task calls
 158    * AliAnalysisTaskSE::Terminate
 159    *
 160    * @param dir Directory to add to
 161    */
 162   virtual void DefineOutput(TList* dir);
 163   /**
 164    * Print information
 165    *
 166    * @param option Not used
 167    */
 168   virtual void Print(Option_t* option="") const;
 169
 170   AliFMDMultCuts& GetLCuts() { return fLCuts; }
 171   AliFMDMultCuts& GetHCuts() { return fHCuts; }
 172   const AliFMDMultCuts& GetLCuts() const { return fLCuts; }
 173   const AliFMDMultCuts& GetHCuts() const { return fHCuts; }
 174   void SetLCuts(const AliFMDMultCuts& c) { fLCuts = c; }
 175   void SetHCuts(const AliFMDMultCuts& c) { fHCuts = c; }
 176 protected:
 177   /**
 178    * Internal data structure to keep track of the histograms
 179    */
 180   struct RingHistos : public AliForwardUtil::RingHistos
 181   {
 182     /**
 183      * Default CTOR
 184      */
 185     RingHistos();
 186     /**
 187      * Constructor
 188      *
 189      * @param d detector
 190      * @param r ring
 191      */
 192     RingHistos(UShort_t d, Char_t r);
 193     /**
 194      * Copy constructor
 195      *
 196      * @param o Object to copy from
 197      */
 198     RingHistos(const RingHistos& o);
 199     /**
 200      * Assignment operator
 201      *
 202      * @param o Object to assign from
 203      *
 204      * @return Reference to this
 205      */
 206     RingHistos& operator=(const RingHistos& o);
 207     /**
 208      * Destructor
 209      */
 210     ~RingHistos();
 211     /**
 212      * Clear this object
 213      */
 214     void Clear(const Option_t* ="") { fNHits = 0; }
 215     /**
 216      * Increase number of hits
 217      *
 218      */
 219     void Incr() { fNHits++; }
 220     /**
 221      * Finish off
 222      *
 223      */
 224     void Finish();
 225     /**
 226      * Make output
 227      *
 228      * @param dir where to store
 229      */
 230     void Output(TList* dir);
 231     /**
 232      * Scale the histograms to the total number of events
 233      *
 234      * @param nEvents Number of events
 235      * @param dir     Where the output is
 236      */
 237     void ScaleHistograms(const TList* dir, Int_t nEvents);
 238     TH1D*     fBefore;       // Distribution of signals before filter
 239     TH1D*     fAfter;        // Distribution of signals after filter
 240     TH1D*     fSingle;       // Distribution of 1 signal after filter
 241     TH1D*     fDouble;       // Distribution of 2 signals after filter
 242     TH1D*     fTriple;       // Distribution of 3 signals after filter
 243     TH2D*     fSinglePerStrip;       // Distribution of 1 signal per strip
 244     TH1D*     fDistanceBefore; //Distance between signals before sharing
 245     TH1D*     fDistanceAfter; //Distance between signals after sharing
 246     TH2D*     fBeforeAfter;  // Correlation of before and after
 247     TH2D*     fNeighborsBefore; // Correlation of neighbors
 248     TH2D*     fNeighborsAfter; // Correlation of neighbors
 249     TH2D*     fSum;          // Summed signal
 250     TH1D*     fHits;         // Distribution of hit strips.
 251     Int_t     fNHits;        // Number of hit strips per event
 252     ClassDef(RingHistos,1);
 253   };
 254   /**
 255    * Get the ring histogram container
 256    *
 257    * @param d Detector
 258    * @param r Ring
 259    *
 260    * @return Ring histogram container
 261    */
 262   RingHistos* GetRingHistos(UShort_t d, Char_t r) const;
 263   /**
 264    * Get the signal in a strip
 265    *
 266    * @param fmd   ESD object
 267    * @param d     Detector
 268    * @param r     Ring
 269    * @param s     Sector
 270    * @param t     Strip
 271    *
 272    * @return The energy signal
 273    */
 274   Double_t SignalInStrip(const AliESDFMD& fmd,
 275                          UShort_t d,
 276                          Char_t   r,
 277                          UShort_t s,
 278                          UShort_t t) const;
 279   /**
 280    * The actual algorithm
 281    *
 282    * @param mult      The unfiltered signal in the strip
 283    * @param eta       Psuedo rapidity
 284    * @param prevE     Previous strip signal (or 0)
 285    * @param nextE     Next strip signal (or 0)
 286    * @param lowFlux   Whether this is a low flux event
 287    * @param d         Detector
 288    * @param r         Ring
 289    * @param s         Sector
 290    * @param t         Strip
 291    * @param usedPrev  Whether the previous strip was used in sharing or not
 292    * @param usedThis  Wether this strip was used in sharing or not.
 293    *
 294    * @return The filtered signal in the strip
 295    */
 296   Double_t MultiplicityOfStrip(Double_t mult,
 297                                Double_t eta,
 298                                Double_t prevE,
 299                                Double_t nextE,
 300                                Bool_t   lowFlux,
 301                                UShort_t d,
 302                                Char_t   r,
 303                                UShort_t s,
 304                                UShort_t t,
 305                                Bool_t&  usedPrev,
 306                                Bool_t&  usedThis) const;
 307   Double_t MultiplicityOfStrip(Double_t thisE,
 308                                Double_t prevE,
 309                                Double_t nextE,
 310                                Double_t eta,
 311                                Bool_t   lowFlux,
 312                                UShort_t d,
 313                                Char_t   r,
 314                                UShort_t s,
 315                                UShort_t t,
 316                                Status&  prevStatus,
 317                                Status&  thisStatus,
 318                                Status&  nextStatus) const;
 319   /**
 320    * Angle correct the signal
 321    *
 322    * @param mult Angle Un-corrected Signal
 323    * @param eta  Pseudo-rapidity
 324    *
 325    * @return Angle corrected signal
 326    */
 327   Double_t AngleCorrect(Double_t mult, Double_t eta) const;
 328   /**
 329    * Angle de-correct the signal
 330    *
 331    * @param mult Angle corrected Signal
 332    * @param eta  Pseudo-rapidity
 333    *
 334    * @return Angle un-corrected signal
 335    */
 336   Double_t DeAngleCorrect(Double_t mult, Double_t eta) const;
 337   /**
 338    * Get the high cut.  The high cut is defined as the
 339    * most-probably-value peak found from the energy distributions, minus
 340    * 2 times the width of the corresponding Landau.
 341    *
 342    * @param d      Detector
 343    * @param r      Ring
 344    * @param eta    Eta value
 345    * @param errors If false, do not show errors
 346    *
 347    * @return 0 or less on failure, otherwise @f$\Delta_{mp}-n\xi@f$
 348    */
 349   virtual Double_t GetHighCut(UShort_t d, Char_t r, Double_t eta,
 350                               Bool_t errors=true) const;
 351   /**
 352    * Get the low cut.  Normally, the low cut is taken to be the lower
 353    * value of the fit range used when generating the energy loss fits.
 354    * However, if fLowCut is set (using SetLowCit) to a value greater
 355    * than 0, then that value is used.
 356    *
 357    * @param d    Detector
 358    * @param r    Ring
 359    * @param eta  Eta value
 360    *
 361    * @return
 362    */
 363   virtual Double_t GetLowCut(UShort_t d, Char_t r, Double_t eta) const;
 364
 365   TList    fRingHistos;    // List of histogram containers
 366   // Double_t fLowCut;        // Low cut on sharing
 367   Bool_t   fCorrectAngles; // Whether to work on angle corrected signals
 368   TH2*     fSummed;        // Operations histogram
 369   TH2*     fHighCuts;      // High cuts used
 370   TH2*     fLowCuts;       // High cuts used
 371   AliFMDFloatMap* fOper;   // Operation done per strip
 372   Int_t    fDebug;         // Debug level
 373   Bool_t   fZeroSharedHitsBelowThreshold; //Whether to zero shared strip below cut
 374   AliFMDMultCuts fLCuts;    //Cuts object for low cuts
 375   AliFMDMultCuts fHCuts;    //Cuts object for high cuts
 376   Bool_t   fUseSimpleMerging; //enable simple sharing by HHD
 377   Bool_t   fThreeStripSharing; //In case of simple sharing allow 3 strips
 378   ClassDef(AliFMDSharingFilter,4); //
 379 };
 380
 381 #endif
 382 // Local Variables:
 383 //  mode: C++
 384 // End: