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