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