]>
Commit | Line | Data |
---|---|---|
e6c798e6 | 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 | { | |
3a9bf5a4 | 349 | const UInt_t target = 64491; |
e6c798e6 | 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: |