AliFMDEncodedEdx.h

Go to the documentation of this file.

#ifndef ALIFMDENCODEDEDX_H
#define ALIFMDENCODEDEDX_H

#ifndef __CINT__
# include <TMath.h>
# include <TArrayD.h>
# include <TH1D.h>
# include <TH2D.h>
#else
class TArrayD;
class TH1;
class TH2;
#endif

class AliFMDEncodedEdx
{
public:
  struct Spec {
    UShort_t nBins;     // Total number of bins
    Double_t min;       // Least value 
    Double_t max;       // Largest value 
    UShort_t cutBin1;   // Last bin (plus one) of lower region
    Double_t cutVal1;   // Upper cut on lower region
    UShort_t cutBin2;   // Last bin (plus one) of middle region
    Double_t cutVal2; // Upper cut on middle region       
    Spec(UShort_t nb, Double_t l, Double_t h)
      : nBins(nb), min(l), max(h), 
  cutBin1(UShort_t(nb * 0.75 + .5)), 
  cutVal1((max-min) / 10 + min),
  cutBin2(UShort_t(nb * 0.95 + .5)),
  cutVal2((max-min) / 5  + min)
    {}
    UInt_t Encode(Double_t v) const
    {
      UInt_t   off  = 0;
      UInt_t   n    = cutBin1;
      Double_t low  = min;
      Double_t high = cutVal1;
      if (v > cutVal2) { 
  // Upper part of the plot
  off  = cutBin2;
  n    = nBins - cutBin2;
  low  = cutVal2;
  high = max;
      }
      else if (v > cutVal1) {
  // Middle part 
  off  = cutBin1;
  n    = cutBin2 - cutBin1;
  low  = cutVal1;
  high = cutVal2;
      }
      return off + UInt_t(n*(v-low)/(high-low));
    }
    Double_t Decode(UInt_t b, Double_t& w) const
    {
      Double_t off  = min;
      UInt_t   n    = cutBin1;
      Double_t high = cutVal1;
      if (b >= cutBin2) {
  // Upper part
  off  = cutVal2;
  n    = nBins - cutBin2;
  high = max;
  b    -= cutBin2;
      }
      else if (b >= cutBin1) {
  // Middle part
  off  = cutVal1;
  n    = cutBin2 - cutBin1;
  high = cutVal2;
  b    -= cutBin1;
      }
      w = (high-off)/n;
      return off + w * (b+.5);
    }
    Double_t Decode(UInt_t b) const 
    {
      Double_t w;
      return Decode(b, w);
    }
    void FillBinArray(TArrayD& a) const
    {
      a.Set(nBins+1);
      a[0] = min;
      Double_t w0 = (cutVal1 - min)     / cutBin1;
      Double_t w1 = (cutVal2 - cutVal1) / (cutBin2 - cutBin1);
      Double_t w2 = (max     - cutVal2) / (nBins   - cutBin2);
      for (UInt_t i = 1;         i <= cutBin1; i++) a[i] = a[i-1] + w0;
      for (UInt_t i = cutBin1+1; i <= cutBin2; i++) a[i] = a[i-1] + w1;
      for (UInt_t i = cutBin2+1; i <= nBins;   i++) a[i] = a[i-1] + w2;
    }
    void Print(Option_t* opt="") const 
    {
      Printf("Spec: [%8.4f,%8.4f] in %3d bins, cuts %8.4f (%3d) %8.4f (%3d)",
       min, max, nBins, cutVal1, cutBin1, cutVal2, cutBin2);
      if (opt[0] == 'T' || opt[0] == 't') {
  for (Int_t i = 0; i < nBins; i++) { 
    Double_t w = 0;
    Double_t x = Decode(i,w );
    UInt_t   j = Encode(x);
    Printf("%3d -> %8.4f (%7.4f) -> %3d", i, x, w, j);
  }
      }
    }
    static void Test() 
    {
      Spec s(125, 0, 125);
      s.Print("T");
    }
  };
  enum { 
    kNEBits = 8, 
    kNLBits = 5,
    kEMask  = 0xFF, // (1 << kNEBits) - 1
    kLMask  = 0x1F  // (1 << kNLBits) - 1
  };
  AliFMDEncodedEdx() {}
  virtual ~AliFMDEncodedEdx() {}
  static const Spec& GetdEAxis() 
  { 
    static Spec* dEAxis = 0;
    if (!dEAxis) dEAxis = new Spec((1<<kNEBits), 0 /*0.0000025*/, 11);
    return *dEAxis;
  }
  static const Spec& GetdLAxis() 
  { 
    static Spec* dLAxis = 0;
    if (!dLAxis) dLAxis = new Spec((1<<kNLBits), 0 /*0.00014*/,   0.7);
    return *dLAxis;
  }
  static UInt_t Encode(Double_t edep, Double_t length)
  {
    UInt_t uE = EncodeOne(edep,   GetdEAxis());
    UInt_t uL = EncodeOne(length, GetdLAxis());
    return (((uE & kEMask) << 0) | ((uL & kLMask) << kNEBits));
  }
  static void Decode(UInt_t bits, Double_t& edep, Double_t& length)
  {
    edep   = DecodeOne((bits >> 0)       & kEMask, GetdEAxis());
    length = DecodeOne((bits >> kNEBits) & kLMask, GetdLAxis());
  }
  static void Decode(UInt_t bits, Double_t& edep, Double_t& length, 
         Double_t& wEdep, Double_t& wLength)
  {
    edep   = DecodeOne((bits >> 0)       & kEMask, wEdep,   GetdEAxis());
    length = DecodeOne((bits >> kNEBits) & kLMask, wLength, GetdLAxis());
  }    
  static TH1* Make1D(const char* name, const char* title, UShort_t mode=true)
  {
    const Spec& a = (mode==0 || mode==2 ? GetdEAxis() : GetdLAxis());
    TArrayD     aa; a.FillBinArray(aa);

    if (mode == 2) 
      // In case we need to do dE/dx, extend the range by a factor 100
      for (Int_t i = 0; i < aa.GetSize(); i++) aa[i] *= 100;

    // Make the histogram 
    TH1* h = new TH1D(name, title, aa.GetSize()-1, aa.GetArray());
    h->SetXTitle(mode == 0 ? "d#Delta [MeV]" : 
     mode == 1 ? "dx [cm]" : 
     mode == 2 ? "d#Delta/dx [MeV/cm]" : "?");
    h->SetFillStyle(3001);
    h->SetMarkerStyle(20);
    h->Sumw2();
    
    return h;
  }
  static TH2* Make2D(const char* name, const char* title, Bool_t xedep=true)
  {
    const Spec& a1 = (xedep ? GetdEAxis() : GetdLAxis());
    const Spec& a2 = (xedep ? GetdLAxis() : GetdEAxis());
    TArrayD     aa1; a1.FillBinArray(aa1);
    TArrayD     aa2; a2.FillBinArray(aa2);
    TH2* h = new TH2D(name, title, 
          aa1.GetSize()-1, aa1.GetArray(),
          aa2.GetSize()-1, aa2.GetArray());
    h->SetXTitle(xedep ? "d#Delta [MeV]" : "dx [cm]");
    h->SetYTitle(xedep ? "dx [cm]"       : "d#Delta [MeV]");
    return h;
  }
  static Bool_t IsAvailable(UInt_t alirootRev) 
  {
    const UInt_t target = 64491;
    return alirootRev >= target;
  }
private:
  static UInt_t EncodeOne(Double_t v, const Spec& a) 
  {
    if (v < a.min) return 0;
    if (v > a.max) return a.nBins;
    return a.Encode(v);
  }
  static Double_t DecodeOne(UInt_t b, const Spec& a) 
  {
    if (b >= a.nBins) b = a.nBins-1;
    return a.Decode(b);
  }
  static Double_t DecodeOne(UInt_t b, Double_t& w, const Spec& a) 
  {
    if (b >= a.nBins) b = a.nBins-1;
    return a.Decode(b, w);
  }

  ClassDef(AliFMDEncodedEdx,1); // En-/Decode dE/dx for/from track references
};
#endif
// Local Variables:
//  mode: C++
// End: