Fix some documentation issues
[u/mrichter/AliRoot.git] / PWGLF / FORWARD / analysis2 / AliFMDEncodedEdx.h
CommitLineData
169c3205 1#ifndef ALIFMDENCODEDEDX_H
2#define ALIFMDENCODEDEDX_H
3/**
671df6c9 4 * @file AliFMDEncodedEdx.h
169c3205 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
17class TArrayD;
18class TH1;
19class 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 */
52class AliFMDEncodedEdx
53{
54public:
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 = 60000;
350 return alirootRev >= target;
351 }
352private:
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: