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