Line data Source code
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 4 : : nBins(nb), min(l), max(h),
76 2 : cutBin1(UShort_t(nb * 0.75 + .5)),
77 2 : cutVal1((max-min) / 10 + min),
78 2 : cutBin2(UShort_t(nb * 0.95 + .5)),
79 2 : cutVal2((max-min) / 5 + min)
80 6 : {}
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 440 : UInt_t n = cutBin1;
92 220 : Double_t low = min;
93 220 : Double_t high = cutVal1;
94 220 : if (v > cutVal2) {
95 : // Upper part of the plot
96 0 : off = cutBin2;
97 0 : n = nBins - cutBin2;
98 : low = cutVal2;
99 0 : high = max;
100 0 : }
101 220 : else if (v > cutVal1) {
102 : // Middle part
103 : off = cutBin1;
104 0 : n = cutBin2 - cutBin1;
105 : low = cutVal1;
106 : high = cutVal2;
107 0 : }
108 220 : 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 222 : if (!dEAxis) dEAxis = new Spec((1<<kNEBits), 0 /*0.0000025*/, 11);
224 110 : return *dEAxis;
225 0 : }
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 222 : if (!dLAxis) dLAxis = new Spec((1<<kNLBits), 0 /*0.00014*/, 0.7);
236 110 : return *dLAxis;
237 0 : }
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 220 : UInt_t uE = EncodeOne(edep, GetdEAxis());
249 110 : UInt_t uL = EncodeOne(length, GetdLAxis());
250 110 : 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 440 : if (v < a.min) return 0;
364 220 : if (v > a.max) return a.nBins;
365 220 : return a.Encode(v);
366 220 : }
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:
|
