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: