1 #ifndef ALIFMDENCODEDEDX_H
2 #define ALIFMDENCODEDEDX_H
4 * @file AliFMDEncodedEdx.C
5 * @author Christian Holm Christensen <cholm@nbi.dk>
6 * @date Thu Oct 10 11:50:59 2013
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.
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
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$).
34 * The bins become progressively wider. That is, we define 3 regions.
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
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 |
48 * Currently there's no support of other schemas
52 class AliFMDEncodedEdx
56 * Specification of the bins
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
70 * @param nb Total number of bins
71 * @param l Lower value
72 * @param h Upper value
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)
84 * @param v Value to encode
86 * @return Encoded (bin number) value
88 UInt_t Encode(Double_t v) const
93 Double_t high = cutVal1;
95 // Upper part of the plot
101 else if (v > cutVal1) {
104 n = cutBin2 - cutBin1;
108 return off + UInt_t(n*(v-low)/(high-low));
111 * Decode a bin into a value
113 * @param b Encoded (bin number) value
114 * @param w On return, the width of the bin
116 * @return Decoded value (center of bin)
118 Double_t Decode(UInt_t b, Double_t& w) const
122 Double_t high = cutVal1;
130 else if (b >= cutBin1) {
133 n = cutBin2 - cutBin1;
138 return off + w * (b+.5);
141 * Decode a bin into a value
143 * @param b Encoded (bin number) value
145 * @return Decoded value (center of bin)
147 Double_t Decode(UInt_t b) const
153 * Fill an array with values appropriate for defining a histogram
154 * axis with the _natural_ binning of the encoding
156 * @param a On return, the modified bin-border array
158 void FillBinArray(TArrayD& a) const
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;
172 * @param opt If this starts with `T' also run a test
174 void Print(Option_t* opt="") const
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++) {
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);
198 * How the 13 bits are distributed
203 kEMask = 0xFF, // (1 << kNEBits) - 1
204 kLMask = 0x1F // (1 << kNLBits) - 1
207 * Constructor - a no-op
209 AliFMDEncodedEdx() {}
211 * Destructor - a no-op
213 virtual ~AliFMDEncodedEdx() {}
215 * Get the @f$d\Delta@f$ bin specification. If not initialized
218 * @return Constant reference to @f$d\Delta@f$ bin specification
220 static const Spec& GetdEAxis()
222 static Spec* dEAxis = 0;
223 if (!dEAxis) dEAxis = new Spec((1<<kNEBits), 0 /*0.0000025*/, 11);
227 * Get the @f$dx@f$ bin specification. If not initialized
230 * @return Constant reference to @f$dx@f$ bin specification
232 static const Spec& GetdLAxis()
234 static Spec* dLAxis = 0;
235 if (!dLAxis) dLAxis = new Spec((1<<kNLBits), 0 /*0.00014*/, 0.7);
239 * Encode @f$d\Delta@f$ and @f$dx@f$ into a 13bit number.
241 * @param edep @f$d\Delta@f$
242 * @param length @f$dx@f$
244 * @return 13-bit (lower) encoded value
246 static UInt_t Encode(Double_t edep, Double_t length)
248 UInt_t uE = EncodeOne(edep, GetdEAxis());
249 UInt_t uL = EncodeOne(length, GetdLAxis());
250 return (((uE & kEMask) << 0) | ((uL & kLMask) << kNEBits));
253 * Decode the lower 13-bit of the input into @f$d\Delta@f$ and @f$dx@f$
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$
259 static void Decode(UInt_t bits, Double_t& edep, Double_t& length)
261 edep = DecodeOne((bits >> 0) & kEMask, GetdEAxis());
262 length = DecodeOne((bits >> kNEBits) & kLMask, GetdLAxis());
265 * Decode the lower 13-bit of the input into @f$d\Delta@f$ and @f$dx@f$
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
273 static void Decode(UInt_t bits, Double_t& edep, Double_t& length,
274 Double_t& wEdep, Double_t& wLength)
276 edep = DecodeOne((bits >> 0) & kEMask, wEdep, GetdEAxis());
277 length = DecodeOne((bits >> kNEBits) & kLMask, wLength, GetdLAxis());
280 * Make a 1-dimension histogram with the natural binning for the
281 * encoding for either @f$d\Delta@f$ or @f$dx@f$
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$
288 * @return Newly allocated histogram
290 static TH1* Make1D(const char* name, const char* title, UShort_t mode=true)
292 const Spec& a = (mode==0 || mode==2 ? GetdEAxis() : GetdLAxis());
293 TArrayD aa; a.FillBinArray(aa);
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;
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);
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)
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
319 * @return Newly allocated histogram
321 static TH2* Make2D(const char* name, const char* title, Bool_t xedep=true)
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]");
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.
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.
344 * @return true if @a alirootRev is larger than some fixed number
345 * set when this class was committed to SVN.
347 static Bool_t IsAvailable(UInt_t alirootRev)
349 const UInt_t target = 64491;
350 return alirootRev >= target;
357 * @param a Bin specification
359 * @return Encoded value
361 static UInt_t EncodeOne(Double_t v, const Spec& a)
363 if (v < a.min) return 0;
364 if (v > a.max) return a.nBins;
370 * @param b Encoded value
371 * @param a Bin specification
373 * @return Decoded value
375 static Double_t DecodeOne(UInt_t b, const Spec& a)
377 if (b >= a.nBins) b = a.nBins-1;
383 * @param b Encoded value
384 * @param a Bin specification
385 * @param w On return, the bin width
387 * @return Decoded value
389 static Double_t DecodeOne(UInt_t b, Double_t& w, const Spec& a)
391 if (b >= a.nBins) b = a.nBins-1;
392 return a.Decode(b, w);
395 ClassDef(AliFMDEncodedEdx,1); // En-/Decode dE/dx for/from track references