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