Handle deletion of AliFMDAltroMapping object. If we make the object
[u/mrichter/AliRoot.git] / FMD / AliFMDEncodedEdx.h
1 #ifndef ALIFMDENCODEDEDX_H
2 #define ALIFMDENCODEDEDX_H
3 /**
4  * @file   AliFMDEncodedEdx.C
5  * @author Christian Holm Christensen <cholm@nbi.dk>
6  * @date   Thu Oct 10 11:50:59 2013
7  * 
8  * @brief  Class to encode the energy loss @f$d\Delta@f$ and path length
9  * @f$dx@f$ of a particle into track reference bits.
10  */
11 #ifndef __CINT__
12 # include <TMath.h>
13 # include <TArrayD.h>
14 # include <TH1D.h>
15 # include <TH2D.h>
16 #else
17 class TArrayD;
18 class TH1;
19 class TH2;
20 #endif
21
22 /**
23  * Class to encode the energy loss @f$d\Delta@f$ and path length
24  * @f$dx@f$ of a particle into track reference bits.  A total of 13
25  * bits are available.  Of these 8 are used for @f$d\Delta@f$ and 5
26  * for @f$dx@f$.
27  *
28  * The encoding is done by binning.  That is, for a given value of
29  * @f$d\Delta@f$ or @f$dx@f$ we calculate a bin number and store that.
30  * The reverse decoding is done by looking up the bin center of the
31  * bin number stored.  Note, that the bin numbers go from 0 to 255
32  * (for @f$d\Delta@f$) and 31 (for @f$dx@f$).
33  *
34  * The bins become progressively wider.  That is, we define 3 regions. 
35  * 
36  * - Lower region which spans the lower 10% of the distribution has
37  *   75% of all avaible bins.
38  * - Middle region which spans from the 10% to 20% of the distribtion
39  *    and has 20% of the available bins.
40  * - Upper region which covers the rest of the distribtion in 5% of
41  *   the available bins
42  *
43  *  | Type          | N bins | Range | bin | value | bin | value |
44  *  | ------------- | ------ | ----- | --- | ----- | --- | ----- |
45  *  | @f$d\Delta@f$ |  256   | 0-11  | 192 |  1.1  | 243 |  2.2  |
46  *  | @f$dx@f$      |   32   | 0-0.7 |  24 |  0.07 |  30 |  1.4  |
47  *      
48  * Currently there's no support of other schemas 
49  * 
50  * 
51  */
52 class AliFMDEncodedEdx
53 {
54 public:
55   /**
56    * Specification of the bins 
57    * 
58    */
59   struct Spec {
60     UShort_t nBins;     // Total number of bins
61     Double_t min;       // Least value 
62     Double_t max;       // Largest value 
63     UShort_t cutBin1;   // Last bin (plus one) of lower region
64     Double_t cutVal1;   // Upper cut on lower region
65     UShort_t cutBin2;   // Last bin (plus one) of middle region
66     Double_t cutVal2;   // Upper cut on middle region         
67     /** 
68      * Constructor 
69      * 
70      * @param nb Total number of bins 
71      * @param l  Lower value 
72      * @param h  Upper value 
73      */
74     Spec(UShort_t nb, Double_t l, Double_t h)
75       : nBins(nb), min(l), max(h), 
76         cutBin1(UShort_t(nb * 0.75 + .5)), 
77         cutVal1((max-min) / 10 + min),
78         cutBin2(UShort_t(nb * 0.95 + .5)),
79         cutVal2((max-min) / 5  + min)
80     {}
81     /** 
82      * Encode a value 
83      * 
84      * @param v Value to encode 
85      * 
86      * @return Encoded (bin number) value 
87      */
88     UInt_t Encode(Double_t v) const
89     {
90       UInt_t   off  = 0;
91       UInt_t   n    = cutBin1;
92       Double_t low  = min;
93       Double_t high = cutVal1;
94       if (v > cutVal2) { 
95         // Upper part of the plot
96         off  = cutBin2;
97         n    = nBins - cutBin2;
98         low  = cutVal2;
99         high = max;
100       }
101       else if (v > cutVal1) {
102         // Middle part 
103         off  = cutBin1;
104         n    = cutBin2 - cutBin1;
105         low  = cutVal1;
106         high = cutVal2;
107       }
108       return off + UInt_t(n*(v-low)/(high-low));
109     }
110     /** 
111      * Decode a bin into a value 
112      * 
113      * @param b Encoded (bin number) value 
114      * @param w On return, the width of the bin
115      * 
116      * @return Decoded value (center of bin)
117      */
118     Double_t Decode(UInt_t b, Double_t& w) const
119     {
120       Double_t off  = min;
121       UInt_t   n    = cutBin1;
122       Double_t high = cutVal1;
123       if (b >= cutBin2) {
124         // Upper part
125         off  = cutVal2;
126         n    = nBins - cutBin2;
127         high = max;
128         b    -= cutBin2;
129       }
130       else if (b >= cutBin1) {
131         // Middle part
132         off  = cutVal1;
133         n    = cutBin2 - cutBin1;
134         high = cutVal2;
135         b    -= cutBin1;
136       }
137       w = (high-off)/n;
138       return off + w * (b+.5);
139     }
140     /** 
141      * Decode a bin into a value 
142      * 
143      * @param b Encoded (bin number) value 
144      * 
145      * @return Decoded value (center of bin)
146      */
147     Double_t Decode(UInt_t b) const 
148     {
149       Double_t w;
150       return Decode(b, w);
151     }
152     /** 
153      * Fill an array with values appropriate for defining a histogram 
154      * axis with the _natural_ binning of the encoding 
155      * 
156      * @param a On return, the modified bin-border array 
157      */
158     void FillBinArray(TArrayD& a) const
159     {
160       a.Set(nBins+1);
161       a[0] = min;
162       Double_t w0 = (cutVal1 - min)     / cutBin1;
163       Double_t w1 = (cutVal2 - cutVal1) / (cutBin2 - cutBin1);
164       Double_t w2 = (max     - cutVal2) / (nBins   - cutBin2);
165       for (UInt_t i = 1;         i <= cutBin1; i++) a[i] = a[i-1] + w0;
166       for (UInt_t i = cutBin1+1; i <= cutBin2; i++) a[i] = a[i-1] + w1;
167       for (UInt_t i = cutBin2+1; i <= nBins;   i++) a[i] = a[i-1] + w2;
168     }
169     /** 
170      * Print information.
171      * 
172      * @param opt If this starts with `T' also run a test 
173      */
174     void Print(Option_t* opt="") const 
175     {
176       Printf("Spec: [%8.4f,%8.4f] in %3d bins, cuts %8.4f (%3d) %8.4f (%3d)",
177              min, max, nBins, cutVal1, cutBin1, cutVal2, cutBin2);
178       if (opt[0] == 'T' || opt[0] == 't') {
179         for (Int_t i = 0; i < nBins; i++) { 
180           Double_t w = 0;
181           Double_t x = Decode(i,w );
182           UInt_t   j = Encode(x);
183           Printf("%3d -> %8.4f (%7.4f) -> %3d", i, x, w, j);
184         }
185       }
186     }
187     /** 
188      * Run a test 
189      * 
190      */
191     static void Test() 
192     {
193       Spec s(125, 0, 125);
194       s.Print("T");
195     }
196   };
197   /**
198    * How the 13 bits are distributed 
199    */
200   enum { 
201     kNEBits = 8, 
202     kNLBits = 5,
203     kEMask  = 0xFF, // (1 << kNEBits) - 1
204     kLMask  = 0x1F  // (1 << kNLBits) - 1
205   };
206   /** 
207    * Constructor - a no-op
208    */
209   AliFMDEncodedEdx() {}
210   /** 
211    * Destructor - a no-op 
212    */
213   virtual ~AliFMDEncodedEdx() {}
214   /** 
215    * Get the @f$d\Delta@f$ bin specification. If not initialized
216    * already, do so .
217    * 
218    * @return Constant reference to @f$d\Delta@f$ bin specification
219    */
220   static const Spec& GetdEAxis() 
221   { 
222     static Spec* dEAxis = 0;
223     if (!dEAxis) dEAxis = new Spec((1<<kNEBits), 0 /*0.0000025*/, 11);
224     return *dEAxis;
225   }
226   /** 
227    * Get the @f$dx@f$ bin specification. If not initialized
228    * already, do so .
229    * 
230    * @return Constant reference to @f$dx@f$ bin specification
231    */
232   static const Spec& GetdLAxis() 
233   { 
234     static Spec* dLAxis = 0;
235     if (!dLAxis) dLAxis = new Spec((1<<kNLBits), 0 /*0.00014*/,   0.7);
236     return *dLAxis;
237   }
238   /** 
239    * Encode @f$d\Delta@f$ and @f$dx@f$ into a 13bit number.  
240    * 
241    * @param edep   @f$d\Delta@f$
242    * @param length @f$dx@f$ 
243    * 
244    * @return 13-bit (lower) encoded value 
245    */
246   static UInt_t Encode(Double_t edep, Double_t length)
247   {
248     UInt_t uE = EncodeOne(edep,   GetdEAxis());
249     UInt_t uL = EncodeOne(length, GetdLAxis());
250     return (((uE & kEMask) << 0) | ((uL & kLMask) << kNEBits));
251   }
252   /** 
253    * Decode the lower 13-bit of the input into @f$d\Delta@f$ and @f$dx@f$ 
254    * 
255    * @param bits   Encoded 13-bit word (lower 13 bit)
256    * @param edep   On return, the @f$d\Delta@f$
257    * @param length On return, the @f$dx@f$ 
258    */
259   static void Decode(UInt_t bits, Double_t& edep, Double_t& length)
260   {
261     edep   = DecodeOne((bits >> 0)       & kEMask, GetdEAxis());
262     length = DecodeOne((bits >> kNEBits) & kLMask, GetdLAxis());
263   }
264   /** 
265    * Decode the lower 13-bit of the input into @f$d\Delta@f$ and @f$dx@f$ 
266    * 
267    * @param bits    Encoded 13-bit word (lower 13 bit)
268    * @param edep    On return, the @f$d\Delta@f$
269    * @param length  On return, the @f$dx@f$ 
270    * @param wEdep   On return, the width of the corresponding @f$d\Delta@f$ bin
271    * @param wLength On return, the width of the corresponding @f$dx@f$ bin
272    */
273   static void Decode(UInt_t bits, Double_t& edep, Double_t& length, 
274                      Double_t& wEdep, Double_t& wLength)
275   {
276     edep   = DecodeOne((bits >> 0)       & kEMask, wEdep,   GetdEAxis());
277     length = DecodeOne((bits >> kNEBits) & kLMask, wLength, GetdLAxis());
278   }    
279   /** 
280    * Make a 1-dimension histogram with the natural binning for the
281    * encoding for either @f$d\Delta@f$ or @f$dx@f$
282    * 
283    * @param name   Name of produced histogram
284    * @param title  Title of produced histogram 
285    * @param mode   If 0, make histogram for @f$d\Delta@f$. If 1
286    * for @f$dx@f$, if 2 for @f$d\Delta/dx@f$
287    * 
288    * @return Newly allocated histogram
289    */
290   static TH1* Make1D(const char* name, const char* title, UShort_t mode=true)
291   {
292     const Spec& a = (mode==0 || mode==2 ? GetdEAxis() : GetdLAxis());
293     TArrayD     aa; a.FillBinArray(aa);
294
295     if (mode == 2) 
296       // In case we need to do dE/dx, extend the range by a factor 100
297       for (Int_t i = 0; i < aa.GetSize(); i++) aa[i] *= 100;
298
299     // Make the histogram 
300     TH1* h = new TH1D(name, title, aa.GetSize()-1, aa.GetArray());
301     h->SetXTitle(mode == 0 ? "d#Delta [MeV]" : 
302                  mode == 1 ? "dx [cm]" : 
303                  mode == 2 ? "d#Delta/dx [MeV/cm]" : "?");
304     h->SetFillStyle(3001);
305     h->SetMarkerStyle(20);
306     h->Sumw2();
307     
308     return h;
309   }
310   /** 
311    * Make a 2-dimension histogram with the natural binning for the
312    * encoding of @f$d\Delta@f$ versus @f$dx@f$ (or vice versa)
313    * 
314    * @param name   Name of produced histogram
315    * @param title  Title of produced histogram 
316    * @param xedep  If true, put @f$d\Delta@f$ on the X-axis, otherwise
317    * for @f$dx@f$
318    * 
319    * @return Newly allocated histogram
320    */
321   static TH2* Make2D(const char* name, const char* title, Bool_t xedep=true)
322   {
323     const Spec& a1 = (xedep ? GetdEAxis() : GetdLAxis());
324     const Spec& a2 = (xedep ? GetdLAxis() : GetdEAxis());
325     TArrayD     aa1; a1.FillBinArray(aa1);
326     TArrayD     aa2; a2.FillBinArray(aa2);
327     TH2* h = new TH2D(name, title, 
328                       aa1.GetSize()-1, aa1.GetArray(),
329                       aa2.GetSize()-1, aa2.GetArray());
330     h->SetXTitle(xedep ? "d#Delta [MeV]" : "dx [cm]");
331     h->SetYTitle(xedep ? "dx [cm]"       : "d#Delta [MeV]");
332     return h;
333   }
334   /** 
335    * Check if the encoded @f$d\Delta@f$ and @f$dx@f$ are available in
336    * the upper 13 bits of the unique ID field of track references.
337    * 
338    * @param alirootRev AliROOT revision of the code that _produced_
339    * the track references.  Note, it cannot be the revision of AliROOT
340    * running, since that can be much newer (or older) than the code
341    * that made the track references.  One should get this information
342    * from the object stored in the ESD tree's user objects.
343    * 
344    * @return true if @a alirootRev is larger than some fixed number
345    * set when this class was committed to SVN.
346    */
347   static Bool_t IsAvailable(UInt_t alirootRev) 
348   {
349     const UInt_t target = 64491;
350     return alirootRev >= target;
351   }
352 private:
353   /** 
354    * Encode one value 
355    * 
356    * @param v  Value 
357    * @param a  Bin specification
358    * 
359    * @return Encoded value 
360    */
361   static UInt_t EncodeOne(Double_t v, const Spec& a) 
362   {
363     if (v < a.min) return 0;
364     if (v > a.max) return a.nBins;
365     return a.Encode(v);
366   }
367   /** 
368    * Decode a value 
369    * 
370    * @param b   Encoded value 
371    * @param a   Bin specification
372    * 
373    * @return Decoded value 
374    */  
375   static Double_t DecodeOne(UInt_t b, const Spec& a) 
376   {
377     if (b >= a.nBins) b = a.nBins-1;
378     return a.Decode(b);
379   }
380   /** 
381    * Decode a value 
382    * 
383    * @param b   Encoded value 
384    * @param a   Bin specification
385    * @param w   On return, the bin width 
386    *
387    * @return Decoded value 
388    */  
389   static Double_t DecodeOne(UInt_t b, Double_t& w, const Spec& a) 
390   {
391     if (b >= a.nBins) b = a.nBins-1;
392     return a.Decode(b, w);
393   }
394
395   ClassDef(AliFMDEncodedEdx,1); // En-/Decode dE/dx for/from track references
396 };
397 #endif
398 // Local Variables:
399 //  mode: C++
400 // End: