[u/mrichter/AliRoot.git] / FMD / AliFMDEncodedEdx.h

#ifndef ALIFMDENCODEDEDX_H
#define ALIFMDENCODEDEDX_H
/**
 * @file   AliFMDEncodedEdx.C
 * @author Christian Holm Christensen <cholm@nbi.dk>
 * @date   Thu Oct 10 11:50:59 2013
 * 
 * @brief  Class to encode the energy loss @f$d\Delta@f$ and path length
 * @f$dx@f$ of a particle into track reference bits.
 */
#ifndef __CINT__
# include <TMath.h>
# include <TArrayD.h>
# include <TH1D.h>
# include <TH2D.h>
#else
class TArrayD;
class TH1;
class TH2;
#endif

/**
 * Class to encode the energy loss @f$d\Delta@f$ and path length
 * @f$dx@f$ of a particle into track reference bits.  A total of 13
 * bits are available.  Of these 8 are used for @f$d\Delta@f$ and 5
 * for @f$dx@f$.
 *
 * The encoding is done by binning.  That is, for a given value of
 * @f$d\Delta@f$ or @f$dx@f$ we calculate a bin number and store that.
 * The reverse decoding is done by looking up the bin center of the
 * bin number stored.  Note, that the bin numbers go from 0 to 255
 * (for @f$d\Delta@f$) and 31 (for @f$dx@f$).
 *
 * The bins become progressively wider.  That is, we define 3 regions. 
 * 
 * - Lower region which spans the lower 10% of the distribution has
 *   75% of all avaible bins.
 * - Middle region which spans from the 10% to 20% of the distribtion
 *    and has 20% of the available bins.
 * - Upper region which covers the rest of the distribtion in 5% of
 *   the available bins
 *
 *  | Type          | N bins | Range | bin | value | bin | value |
 *  | ------------- | ------ | ----- | --- | ----- | --- | ----- |
 *  | @f$d\Delta@f$ |  256   | 0-11  | 192 |  1.1  | 243 |  2.2  |
 *  | @f$dx@f$      |   32   | 0-0.7 |  24 |  0.07 |  30 |  1.4  |
 *      
 * Currently there's no support of other schemas 
 * 
 * 
 */
class AliFMDEncodedEdx
{
public:
  /**
   * Specification of the bins 
   * 
   */
  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	      
    /** 
     * Constructor 
     * 
     * @param nb Total number of bins 
     * @param l  Lower value 
     * @param h  Upper value 
     */
    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)
    {}
    /** 
     * Encode a value 
     * 
     * @param v Value to encode 
     * 
     * @return Encoded (bin number) value 
     */
    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));
    }
    /** 
     * Decode a bin into a value 
     * 
     * @param b Encoded (bin number) value 
     * @param w On return, the width of the bin
     * 
     * @return Decoded value (center of bin)
     */
    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);
    }
    /** 
     * Decode a bin into a value 
     * 
     * @param b Encoded (bin number) value 
     * 
     * @return Decoded value (center of bin)
     */
    Double_t Decode(UInt_t b) const 
    {
      Double_t w;
      return Decode(b, w);
    }
    /** 
     * Fill an array with values appropriate for defining a histogram 
     * axis with the _natural_ binning of the encoding 
     * 
     * @param a On return, the modified bin-border array 
     */
    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;
    }
    /** 
     * Print information.
     * 
     * @param opt If this starts with `T' also run a test 
     */
    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);
	}
      }
    }
    /** 
     * Run a test 
     * 
     */
    static void Test() 
    {
      Spec s(125, 0, 125);
      s.Print("T");
    }
  };
  /**
   * How the 13 bits are distributed 
   */
  enum { 
    kNEBits = 8, 
    kNLBits = 5,
    kEMask  = 0xFF, // (1 << kNEBits) - 1
    kLMask  = 0x1F  // (1 << kNLBits) - 1
  };
  /** 
   * Constructor - a no-op
   */
  AliFMDEncodedEdx() {}
  /** 
   * Destructor - a no-op 
   */
  virtual ~AliFMDEncodedEdx() {}
  /** 
   * Get the @f$d\Delta@f$ bin specification. If not initialized
   * already, do so .
   * 
   * @return Constant reference to @f$d\Delta@f$ bin specification
   */
  static const Spec& GetdEAxis() 
  { 
    static Spec* dEAxis = 0;
    if (!dEAxis) dEAxis = new Spec((1<<kNEBits), 0 /*0.0000025*/, 11);
    return *dEAxis;
  }
  /** 
   * Get the @f$dx@f$ bin specification. If not initialized
   * already, do so .
   * 
   * @return Constant reference to @f$dx@f$ bin specification
   */
  static const Spec& GetdLAxis() 
  { 
    static Spec* dLAxis = 0;
    if (!dLAxis) dLAxis = new Spec((1<<kNLBits), 0 /*0.00014*/,   0.7);
    return *dLAxis;
  }
  /** 
   * Encode @f$d\Delta@f$ and @f$dx@f$ into a 13bit number.  
   * 
   * @param edep   @f$d\Delta@f$
   * @param length @f$dx@f$ 
   * 
   * @return 13-bit (lower) encoded value 
   */
  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));
  }
  /** 
   * Decode the lower 13-bit of the input into @f$d\Delta@f$ and @f$dx@f$ 
   * 
   * @param bits   Encoded 13-bit word (lower 13 bit)
   * @param edep   On return, the @f$d\Delta@f$
   * @param length On return, the @f$dx@f$ 
   */
  static void Decode(UInt_t bits, Double_t& edep, Double_t& length)
  {
    edep   = DecodeOne((bits >> 0)       & kEMask, GetdEAxis());
    length = DecodeOne((bits >> kNEBits) & kLMask, GetdLAxis());
  }
  /** 
   * Decode the lower 13-bit of the input into @f$d\Delta@f$ and @f$dx@f$ 
   * 
   * @param bits    Encoded 13-bit word (lower 13 bit)
   * @param edep    On return, the @f$d\Delta@f$
   * @param length  On return, the @f$dx@f$ 
   * @param wEdep   On return, the width of the corresponding @f$d\Delta@f$ bin
   * @param wLength On return, the width of the corresponding @f$dx@f$ bin
   */
  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());
  }    
  /** 
   * Make a 1-dimension histogram with the natural binning for the
   * encoding for either @f$d\Delta@f$ or @f$dx@f$
   * 
   * @param name   Name of produced histogram
   * @param title  Title of produced histogram 
   * @param mode   If 0, make histogram for @f$d\Delta@f$. If 1
   * for @f$dx@f$, if 2 for @f$d\Delta/dx@f$
   * 
   * @return Newly allocated histogram
   */
  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;
  }
  /** 
   * Make a 2-dimension histogram with the natural binning for the
   * encoding of @f$d\Delta@f$ versus @f$dx@f$ (or vice versa)
   * 
   * @param name   Name of produced histogram
   * @param title  Title of produced histogram 
   * @param xedep  If true, put @f$d\Delta@f$ on the X-axis, otherwise
   * for @f$dx@f$
   * 
   * @return Newly allocated histogram
   */
  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;
  }
  /** 
   * Check if the encoded @f$d\Delta@f$ and @f$dx@f$ are available in
   * the upper 13 bits of the unique ID field of track references.
   * 
   * @param alirootRev AliROOT revision of the code that _produced_
   * the track references.  Note, it cannot be the revision of AliROOT
   * running, since that can be much newer (or older) than the code
   * that made the track references.  One should get this information
   * from the object stored in the ESD tree's user objects.
   * 
   * @return true if @a alirootRev is larger than some fixed number
   * set when this class was committed to SVN.
   */
  static Bool_t IsAvailable(UInt_t alirootRev) 
  {
    const UInt_t target = 64491;
    return alirootRev >= target;
  }
private:
  /** 
   * Encode one value 
   * 
   * @param v  Value 
   * @param a  Bin specification
   * 
   * @return Encoded value 
   */
  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);
  }
  /** 
   * Decode a value 
   * 
   * @param b   Encoded value 
   * @param a   Bin specification
   * 
   * @return Decoded value 
   */  
  static Double_t DecodeOne(UInt_t b, const Spec& a) 
  {
    if (b >= a.nBins) b = a.nBins-1;
    return a.Decode(b);
  }
  /** 
   * Decode a value 
   * 
   * @param b   Encoded value 
   * @param a   Bin specification
   * @param w   On return, the bin width 
   *
   * @return Decoded value 
   */  
  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:
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: