Histo2DActions/Histo2DAction/Histo2DAction.h

/*
 * Histo2DAction.h
 *
 *  Created on: 01/set/2009
 *      Author: Nicola Mori
 */

/*! @file Histo2DAction.h The Histo2DAction class declaration file. */

#ifndef HISTO2DACTION_H_
#define HISTO2DACTION_H_

#include "../../CollectionAction/CollectionAction.h"

#include <TH2I.h>
#include <TH2F.h>
#include <TH2D.h>

/*! @brief An abstract action that fills a 2-dimensional histogram.
 *
 * This abstract class provides all the common methods needed to produce 2-dimensional histograms,
 * like reading the bins from a file, setting up the root histogram and so on. It handles histograms
 * both as a ROOT histogram and as a matrix. The OnGood() method has to be implemented in children
 * classes, with the specific physical variables to fill the histogram. Eg., a dE/dx vs. rigidity
 * histogram class derived from Histo2DAction could have this OnGood() implementation:
 *
 * Fill(event->GetTrack(0)->GetTrkTrack()->GetRigidity(), event->GetTrack(0)->GetTrkTrack()->GetDEDX());
 *
 * Fill() automatically fills the ROOT and the matrix histogram.
 * The template argument HistoType is the type of histogram (Int_t, Double_t,...). In current
 * implementation, ROOT histograms are only supported for Int_t, Float_t and Double_t. For all the
 * other types, only the the matrix histogram with text file output will be produced. If you want to
 * implement some new types of ROOT histograms, please specialize the template definitions of _CreateHisto()
 * and ~Histo2DAction() (see the specializations in the .cpp file).
 *
 * */
template<class HistoType>
class Histo2DAction: public CollectionAction {
public:

  /*! @brief Constructor.
   *
   * The histogram binning is defined by the content of the bins vector.
   *
   * @param actionName The action's name.
   * @param title The ROOT histogram title.
   * @param outFileBase The file base name for the histogram output. If "", no output file will be produced.
   * @param mode The mode of ROOT file creation (see documentation of TFile constructor
   *             in ROOT's reference guide).
   * @param outRoot If true, an output ROOT file named outFileBase + ".root" will be produced.
   * @param outText If true, an output text file named outFileBase + ".txt" will be produced. It will overwrite an
   *                eventually existing file with the same name.
   */
  Histo2DAction(const char *actionName, TString title = "", TString outFileBase = "", TString mode = "UPDATE",
      bool outRoot = false, bool outText = false);

  /*! @brief Destructor
   *
   * This has to be specialized for each kind of ROOT histogram.
   *
   * */
  ~Histo2DAction();

  /*! @brief Sets the X axis using a binning read from a a vector.
   *
   * @param label The axis' label.
   * @param bins A vector containing the histogram binning (with sign, eventually; from lowest to highest).
   */
  void SetXAxis(TString label, vector<float> &bins);
  /*! @brief Sets the X axis reading the binning from a file.
   *
   * @param label The axis' label.
   * @param binsFile The file containing the bins (with sign, eventually; from lowest to highest).
   */
  void SetXAxis(TString label, TString binsFile);
  /*! @brief Sets the X axis specifying the binning parameters.
   *
   * @param label The axis' label.
   * @param nBins The number of bins.
   * @param min The lower limit of the axis.
   * @param max The upper limit of the axis.
   * @param logBinning If true, bins will be logarithmically spaced.
   */
  void SetXAxis(TString label, unsigned int nBins, float min, float max, bool logBinning = false);

  /*! @brief Sets the Y axis using a binning read from a a vector.
   *
   * @param label The axis' label.
   * @param bins A vector containing the histogram binning (with sign, eventually; from lowest to highest).
   */
  void SetYAxis(TString label, vector<float> &bins);
  /*! @brief Sets the Y axis reading the binning from a file.
   *
   * @param label The axis' label.
   * @param binsFile The file containing the bins (with sign, eventually; from lowest to highest).
   */
  void SetYAxis(TString label, TString binsFile);

  /*! @brief Sets the Y axis specifying the binning parameters.
   *
   * @param label The axis' label.
   * @param nBins The number of bins.
   * @param min The lower limit of the axis.
   * @param max The upper limit of the axis.
   * @param logBinning If true, bins will be logarithmically spaced.
   */
  void SetYAxis(TString label, unsigned int nBins, float min, float max, bool logBinning = false);

  /*! @brief Sets the ROOT histogram's title.
   *
   * @param title The histogram title as it will appear on the histogram itself.
   */
  void SetTitle(TString &title) {
    _title = title;
  }
  /*! @brief Sets up the histogram
   *
   * This routine effectively prepares the histogram, after the desired parameters has been set by #SetXAxis() and #SetYAxis().
   *
   * @param events Pointer to PamLevel2 events (unused).
   */
  void Setup(PamLevel2 *events) {
    CollectionAction::Setup(events);
    _InitHistos();

  }

  /*! @ brief Post-analysis tasks.
   *
   * This method writes the output on a ROOT file if the proper option was set in the constructor.
   */
  void Finalize();

  /*! @brief Returns the histogram.
   *
   * Element [i][j] is the value of the i-th bin on the Y axis and on the j-th bin on the X axis. This sort
   * of index inversion is due to logically maintain the Y values on the vertical axis of the matrix, which
   * is indexed by the row number (which conventionally is the first index of a matrix element).
   *
   * @return A reference to a matrix containing the values of the bins of the histogram.
   */
  SimpleMatrix<HistoType> &GetHisto() {
    return _histo;
  }

  /*! @brief Returns a pointer to the ROOT histogram.
   *
   * @return A pointer to the root histogram
   */
  TH2 *GetRootHisto() {
    return _rootHisto;
  }

  /*! Fills the ROOT and the vector histogram.
   *
   * @param xValue The value of the X coordinate associated to the event.
   * @param yValue The value of the Y coordinate associated to the event.
   * @param weight The weight which will be applied to the event.
   */
  void Fill(double xValue, double yValue, double weight = 1.);

  /*! @brief Gets the X overflow histogram.
   *
   * The returned vector contains the histogram filled with the X overflow events, eg.,
   * those events whose X variable falls above the upper limit of the X axis. This histogram
   * is binned like the Y axis.
   *
   * @return The X overflow histogram.
   */
  vector<HistoType>& GetXOverflow() {
    return _xOverflow;
  }

  /*! @brief Gets the X underflow histogram.
   *
   * The returned vector contains the histogram filled with the X underflow events, eg.,
   * those events whose X variable falls below the lower limit of the X axis. This histogram
   * is binned like the Y axis.
   *
   * @return The X underflow histogram.
   */
  vector<HistoType>& GetXUnderflow() {
    return _xUnderflow;
  }

  /*! @brief Gets the Y overflow histogram.
   *
   * The returned vector contains the histogram filled with the Y overflow events, eg.,
   * those events whose Y variable falls above the upper limit of the Y axis. This histogram
   * is binned like the X axis.
   *
   * @return The Y overflow histogram.
   */
  vector<HistoType>& GetYOverflow() {
    return _yOverflow;
  }

  /*! @brief Gets the Y underflow histogram.
   *
   * The returned vector contains the histogram filled with the Y underflow events, eg.,
   * those events whose Y variable falls below the lower limit of the Y axis. This histogram
   * is binned like the X axis.
   *
   * @return The Y underflow histogram.
   */
  vector<HistoType>& GetYUnderflow() {
    return _yUnderflow;
  }

  /*! @brief Gets the X-Y overflow histogram.
   *
   * This single-bin histogram is built with events which fall above the lower limits of both the
   * X and Y axis.
   *
   * @return The X-Y overflow histogram.
   */
  HistoType GetXOverYOverflow() {
    return _xOverYOverflow;
  }

  /*! @brief Gets the X-Y underflow histogram.
   *
   * This single-bin histogram is built with events fall below the lower limits of both the X and Y axis.
   *
   * @return The X-Y overflow histogram.
   */
  HistoType GetXUnderYUnderflow() {
    return _xUnderYUnderflow;
  }
  /*! @brief Gets the (X over - Y under)flow histogram.
   *
   * This single-bin histogram is built with events fall above the lower limit of the Xaxis  and below
   * that of the Y axis.
   *
   * @return The (X over -Y under)flow histogram.
   */
  HistoType GetXOverYUnderflow() {
    return _xOverYUnderflow;
  }

  /*! @brief Gets the (X under - Y over)flow histogram.
   *
   * This single-bin histogram is built with events fall below the lower limit of the X axis and above
   * that of the Y axis.
   *
   * @return The (X under -Y over)flow histogram.
   */
  HistoType GetXUnderYOverflow() {
    return _xUnderYOverflow;
  }

protected:

  /*! @brief The vector containing the limits of the X bins(from lower to higher). */
  std::vector<float> _xBins;
  /*! @brief The vector containing the limits of the Y bins(from lower to higher). */
  std::vector<float> _yBins;
  /*! @brief A matrix containing the value of the histogram for each X-Y bin. */
  SimpleMatrix<HistoType> _histo;
  /*! @brief The ROOT histogram. */
  TH2 *_rootHisto;

  /*! @brief Base name of the output file. */
  TString _outFileBase;
  /*! @brief Output file open mode (UPDATE or RECREATE, see documentation of TFile). */
  TString _mode;
  /*! @brief Title for the ROOT histogram. */
  TString _title;
  /*! @brief X axis label for the ROOT histogram. */
  TString _xLabel;
  /*! @brief Y axis label for the ROOT histogram. */
  TString _yLabel;

private:

  vector<HistoType> _xUnderflow, _xOverflow, _yUnderflow, _yOverflow;
  HistoType _xUnderYUnderflow, _xOverYOverflow, _xUnderYOverflow, _xOverYUnderflow;
  bool _outRoot;
  bool _outText;
  void _CreateHisto();
  void _InitHistos();
};

template<class HistoType>
void Histo2DAction<HistoType>::_CreateHisto() {

  // No ROOT histogram for generic type; see template specializations in .cpp file.
  _rootHisto = NULL;
}

// Specializations for _CreateHistos(). See Histo2DAction.cpp
template<>
void Histo2DAction<Int_t>::_CreateHisto();

template<>
void Histo2DAction<Float_t>::_CreateHisto();

template<>
void Histo2DAction<Double_t>::_CreateHisto();

template<class HistoType>
void Histo2DAction<HistoType>::_InitHistos() {

  _CreateHisto();
  if (_xBins.size() < 2) // SetXAxis not called by the main program, or wrongly filled (only 1 bin limit)
    SetXAxis("Default X", 10, 0., 1.);
  if (_yBins.size() < 2) // SetYAxis not called by the main program, or wrongly filled (only 1 bin limit)
    SetYAxis("Default Y", 10, 0., 1.);

  if (_rootHisto) {
    Double_t *auxXArray = new Double_t[_xBins.size()];

    for (unsigned int i = 0; i < _xBins.size(); i++) {
      auxXArray[i] = _xBins[i];
    }

    Double_t *auxYArray = new Double_t[_yBins.size()];

    for (unsigned int i = 0; i < _yBins.size(); i++) {
      auxYArray[i] = _yBins[i];
    }

    _rootHisto->SetBins(_xBins.size() - 1, auxXArray, _yBins.size() - 1, auxYArray);
    _rootHisto->SetName(GetName());
    _rootHisto->SetTitle(_title);
    _rootHisto->SetXTitle(_xLabel);
    _rootHisto->SetYTitle(_yLabel);

    delete[] auxXArray;
    delete[] auxYArray;
  }

  /* The row index (first) corresponds to the position on the vertical (Y) axis. */
  _histo.Resize(_yBins.size() - 1, _xBins.size() - 1);

  _xUnderflow.resize(_yBins.size() - 1);
  _xOverflow.resize(_yBins.size() - 1);
  _yUnderflow.resize(_xBins.size() - 1);
  _yOverflow.resize(_xBins.size() - 1);

}

template<class HistoType>
Histo2DAction<HistoType>::~Histo2DAction() {

  delete _rootHisto;
  _rootHisto = NULL;
}

template<class HistoType>
Histo2DAction<HistoType>::Histo2DAction(const char *actionName, TString title, TString outFileBase, TString mode,
    bool outRoot, bool outText) :
  CollectionAction(actionName), _xBins(0), _yBins(0), _histo(0, 0), _rootHisto(NULL), _outFileBase(outFileBase), _mode(
      mode), _title(title), _xLabel(""), _yLabel(""), _xUnderflow(0), _xOverflow(0), _yUnderflow(0), _yOverflow(0),
      _xUnderYUnderflow((HistoType) 0), _xOverYOverflow((HistoType) 0), _xUnderYOverflow((HistoType) 0),
      _xOverYUnderflow((HistoType) 0), _outRoot(outRoot), _outText(outText) {

}

template<class HistoType>
void Histo2DAction<HistoType>::SetXAxis(TString label, vector<float> &bins) {

  _xBins = bins;
  _xLabel = label;

}

template<class HistoType>
void Histo2DAction<HistoType>::SetYAxis(TString label, vector<float> &bins) {

  _yBins = bins;
  _yLabel = label;

}

template<class HistoType>
void Histo2DAction<HistoType>::SetXAxis(TString label, TString binsFile) {

  // Reading the bins from file
  ifstream binListFile;
  binListFile.open(binsFile);

  TString auxString;
  _xBins.resize(0);
  while (!binListFile.eof()) {
    binListFile >> auxString;
    if (auxString != "") {
      _xBins.push_back(auxString.Atof());
    }
  }
  binListFile.close();

  _xLabel = label;

}

template<class HistoType>
void Histo2DAction<HistoType>::SetYAxis(TString label, TString binsFile) {

  // Reading the bins from file
  ifstream binListFile;
  binListFile.open(binsFile);

  TString auxString;
  _yBins.resize(0);
  while (!binListFile.eof()) {
    binListFile >> auxString;
    if (auxString != "") {
      _yBins.push_back(auxString.Atof());
    }
  }
  binListFile.close();

  _yLabel = label;

}

template<class HistoType>
void Histo2DAction<HistoType>::SetXAxis(TString label, unsigned int nBins, float min, float max, bool logBinning) {

  _xBins.resize(nBins + 1);

  if (!logBinning || (logBinning && min <= 0.)) {

#ifdef DEBUGPAMCUT
    if (logbinning && rigMin <= 0.)
    cout << "Warning: logarithmic binning was chosen for X axis but min <= 0. Using linear binning."
#endif

    float step = (max - min) / nBins;
    for (unsigned int i = 0; i < nBins + 1; i++) {
      _xBins[i] = min + i * step;
    }

  }
  else {

    double maxExp = log10(max / min);
    for (unsigned int i = 0; i < nBins + 1; i++) {
      _xBins[i] = min * pow(10., (double) i / ((double) nBins) * maxExp);
    }

  }

  _xLabel = label;
}

template<class HistoType>
void Histo2DAction<HistoType>::SetYAxis(TString label, unsigned int nBins, float min, float max, bool logBinning) {

  _yBins.resize(nBins + 1);

  if (!logBinning || (logBinning && min <= 0.)) {

#ifdef DEBUGPAMCUT
    if (logbinning && rigMin <= 0.)
    cout << "Warning: logarithmic binning was chosen for Y axis but min <= 0. Using linear binning."
#endif

    float step = (max - min) / nBins;
    for (unsigned int i = 0; i < nBins + 1; i++) {
      _yBins[i] = min + i * step;
    }

  }
  else {

    double maxExp = log10(max / min);
    for (unsigned int i = 0; i < nBins + 1; i++) {
      _yBins[i] = min * pow(10., (double) i / ((double) nBins) * maxExp);
    }

  }

  _yLabel = label;
}

template<class HistoType>
void Histo2DAction<HistoType>::Fill(double xValue, double yValue, double weight) {

  _rootHisto->Fill(xValue, yValue, weight);

  int xBin = 0, yBin = 0;

  bool UOflow = false;

  if (xValue < _xBins.front()) {
    xBin = -1;
    UOflow = true;
  }

  if (xValue >= _xBins.back()) {
    xBin = _xBins.size();
    UOflow = true;
  }

  if (yValue < _yBins.front()) {
    yBin = -1;
    UOflow = true;
  }

  if (yValue >= _yBins.back()) {
    yBin = _yBins.size();
    UOflow = true;
  }

  if (xBin == 0) {

    while (xValue >= _xBins[xBin]) {
      xBin++;
    }
    xBin--;
  }

  if (yBin == 0) {

    while (yValue >= _yBins[yBin]) {
      yBin++;
    }
    yBin--;
  }

  if (!UOflow) {
    // No under or over flow.
    _histo[yBin][xBin] += (HistoType) weight;
    return;
  }
  else {
    // An under or over flow has happened.
    if (xBin == -1) {
      if (yBin == -1) {
        _xUnderYUnderflow += (HistoType) weight;
        return;
      }
      else {
        if (yBin == (int) _yBins.size()) {
          _xUnderYOverflow += (HistoType) weight;
          return;
        }
        else {
          _xUnderflow[yBin] += (HistoType) weight;
          return;
        }
      }
    }

    if (xBin == (int) _xBins.size()) {
      if (yBin == -1) {
        _xOverYUnderflow += (HistoType) weight;
        return;
      }
      else {
        if (yBin == (int) _yBins.size()) {
          _xOverYOverflow += (HistoType) weight;
          return;
        }
        else {
          _xOverflow[yBin] += (HistoType) weight;
          return;
        }
      }
    }

    if (yBin == -1) {
      _yUnderflow[xBin] += (HistoType) weight;
      return;
    }

    if (yBin == (int) _yBins.size()) {
      _yOverflow[xBin] += (HistoType) weight;
      return;
    }

  }
}

template<class HistoType>
void Histo2DAction<HistoType>::Finalize() {

  if (_outFileBase != "") {
    // Write the ROOT file
    if (_rootHisto && _outRoot) {
      TFile outRootFile((_outFileBase + ".root"), _mode);
      outRootFile.cd();
      _rootHisto->Write();
      outRootFile.Close();
    }

    //Write the text file
    if (_outText) {
      ofstream outTextFile((_outFileBase + ".txt").Data(), ios_base::out);
      for (unsigned int i = 0; i < _histo.GetNRows(); i++) {
        for (unsigned int j = 0; j < _histo.GetNCols(); j++) {
          outTextFile << setw(7) << _histo[i][j] << "  ";
        }
        outTextFile << "\n";
      }
      outTextFile.close();
    }
  }

}
#endif /* HISTO2DACTION_H_ */
1	/*
2	* Histo2DAction.h
3	*
4	* Created on: 01/set/2009
5	* Author: Nicola Mori
6	*/
7
8	/! @file Histo2DAction.h The Histo2DAction class declaration file. /
9
10	#ifndef HISTO2DACTION_H_
11	#define HISTO2DACTION_H_
12
13	#include "../../CollectionAction/CollectionAction.h"
14
15	#include <TH2I.h>
16	#include <TH2F.h>
17	#include <TH2D.h>
18
19	/*! @brief An abstract action that fills a 2-dimensional histogram.
20	*
21	* This abstract class provides all the common methods needed to produce 2-dimensional histograms,
22	* like reading the bins from a file, setting up the root histogram and so on. It handles histograms
23	* both as a ROOT histogram and as a matrix. The OnGood() method has to be implemented in children
24	* classes, with the specific physical variables to fill the histogram. Eg., a dE/dx vs. rigidity
25	* histogram class derived from Histo2DAction could have this OnGood() implementation:
26	*
27	* Fill(event->GetTrack(0)->GetTrkTrack()->GetRigidity(), event->GetTrack(0)->GetTrkTrack()->GetDEDX());
28	*
29	* Fill() automatically fills the ROOT and the matrix histogram.
30	* The template argument HistoType is the type of histogram (Int_t, Double_t,...). In current
31	* implementation, ROOT histograms are only supported for Int_t, Float_t and Double_t. For all the
32	* other types, only the the matrix histogram with text file output will be produced. If you want to
33	* implement some new types of ROOT histograms, please specialize the template definitions of _CreateHisto()
34	* and ~Histo2DAction() (see the specializations in the .cpp file).
35	*
36	* */
37	template<class HistoType>
38	class Histo2DAction: public CollectionAction {
39	public:
40
41	/*! @brief Constructor.
42	*
43	* The histogram binning is defined by the content of the bins vector.
44	*
45	* @param actionName The action's name.
46	* @param title The ROOT histogram title.
47	* @param outFileBase The file base name for the histogram output. If "", no output file will be produced.
48	* @param mode The mode of ROOT file creation (see documentation of TFile constructor
49	* in ROOT's reference guide).
50	* @param outRoot If true, an output ROOT file named outFileBase + ".root" will be produced.
51	* @param outText If true, an output text file named outFileBase + ".txt" will be produced. It will overwrite an
52	* eventually existing file with the same name.
53	*/
54	Histo2DAction(const char *actionName, TString title = "", TString outFileBase = "", TString mode = "UPDATE",
55	bool outRoot = false, bool outText = false);
56
57	/*! @brief Destructor
58	*
59	* This has to be specialized for each kind of ROOT histogram.
60	*
61	* */
62	~Histo2DAction();
63
64	/*! @brief Sets the X axis using a binning read from a a vector.
65	*
66	* @param label The axis' label.
67	* @param bins A vector containing the histogram binning (with sign, eventually; from lowest to highest).
68	*/
69	void SetXAxis(TString label, vector<float> &bins);
70	/*! @brief Sets the X axis reading the binning from a file.
71	*
72	* @param label The axis' label.
73	* @param binsFile The file containing the bins (with sign, eventually; from lowest to highest).
74	*/
75	void SetXAxis(TString label, TString binsFile);
76	/*! @brief Sets the X axis specifying the binning parameters.
77	*
78	* @param label The axis' label.
79	* @param nBins The number of bins.
80	* @param min The lower limit of the axis.
81	* @param max The upper limit of the axis.
82	* @param logBinning If true, bins will be logarithmically spaced.
83	*/
84	void SetXAxis(TString label, unsigned int nBins, float min, float max, bool logBinning = false);
85
86	/*! @brief Sets the Y axis using a binning read from a a vector.
87	*
88	* @param label The axis' label.
89	* @param bins A vector containing the histogram binning (with sign, eventually; from lowest to highest).
90	*/
91	void SetYAxis(TString label, vector<float> &bins);
92	/*! @brief Sets the Y axis reading the binning from a file.
93	*
94	* @param label The axis' label.
95	* @param binsFile The file containing the bins (with sign, eventually; from lowest to highest).
96	*/
97	void SetYAxis(TString label, TString binsFile);
98
99	/*! @brief Sets the Y axis specifying the binning parameters.
100	*
101	* @param label The axis' label.
102	* @param nBins The number of bins.
103	* @param min The lower limit of the axis.
104	* @param max The upper limit of the axis.
105	* @param logBinning If true, bins will be logarithmically spaced.
106	*/
107	void SetYAxis(TString label, unsigned int nBins, float min, float max, bool logBinning = false);
108
109	/*! @brief Sets the ROOT histogram's title.
110	*
111	* @param title The histogram title as it will appear on the histogram itself.
112	*/
113	void SetTitle(TString &title) {
114	_title = title;
115	}
116	/*! @brief Sets up the histogram
117	*
118	* This routine effectively prepares the histogram, after the desired parameters has been set by #SetXAxis() and #SetYAxis().
119	*
120	* @param events Pointer to PamLevel2 events (unused).
121	*/
122	void Setup(PamLevel2 *events) {
123	CollectionAction::Setup(events);
124	_InitHistos();
125
126	}
127
128	/*! @ brief Post-analysis tasks.
129	*
130	* This method writes the output on a ROOT file if the proper option was set in the constructor.
131	*/
132	void Finalize();
133
134	/*! @brief Returns the histogram.
135	*
136	* Element [i][j] is the value of the i-th bin on the Y axis and on the j-th bin on the X axis. This sort
137	* of index inversion is due to logically maintain the Y values on the vertical axis of the matrix, which
138	* is indexed by the row number (which conventionally is the first index of a matrix element).
139	*
140	* @return A reference to a matrix containing the values of the bins of the histogram.
141	*/
142	SimpleMatrix<HistoType> &GetHisto() {
143	return _histo;
144	}
145
146	/*! @brief Returns a pointer to the ROOT histogram.
147	*
148	* @return A pointer to the root histogram
149	*/
150	TH2 *GetRootHisto() {
151	return _rootHisto;
152	}
153
154	/*! Fills the ROOT and the vector histogram.
155	*
156	* @param xValue The value of the X coordinate associated to the event.
157	* @param yValue The value of the Y coordinate associated to the event.
158	* @param weight The weight which will be applied to the event.
159	*/
160	void Fill(double xValue, double yValue, double weight = 1.);
161
162	/*! @brief Gets the X overflow histogram.
163	*
164	* The returned vector contains the histogram filled with the X overflow events, eg.,
165	* those events whose X variable falls above the upper limit of the X axis. This histogram
166	* is binned like the Y axis.
167	*
168	* @return The X overflow histogram.
169	*/
170	vector<HistoType>& GetXOverflow() {
171	return _xOverflow;
172	}
173
174	/*! @brief Gets the X underflow histogram.
175	*
176	* The returned vector contains the histogram filled with the X underflow events, eg.,
177	* those events whose X variable falls below the lower limit of the X axis. This histogram
178	* is binned like the Y axis.
179	*
180	* @return The X underflow histogram.
181	*/
182	vector<HistoType>& GetXUnderflow() {
183	return _xUnderflow;
184	}
185
186	/*! @brief Gets the Y overflow histogram.
187	*
188	* The returned vector contains the histogram filled with the Y overflow events, eg.,
189	* those events whose Y variable falls above the upper limit of the Y axis. This histogram
190	* is binned like the X axis.
191	*
192	* @return The Y overflow histogram.
193	*/
194	vector<HistoType>& GetYOverflow() {
195	return _yOverflow;
196	}
197
198	/*! @brief Gets the Y underflow histogram.
199	*
200	* The returned vector contains the histogram filled with the Y underflow events, eg.,
201	* those events whose Y variable falls below the lower limit of the Y axis. This histogram
202	* is binned like the X axis.
203	*
204	* @return The Y underflow histogram.
205	*/
206	vector<HistoType>& GetYUnderflow() {
207	return _yUnderflow;
208	}
209
210	/*! @brief Gets the X-Y overflow histogram.
211	*
212	* This single-bin histogram is built with events which fall above the lower limits of both the
213	* X and Y axis.
214	*
215	* @return The X-Y overflow histogram.
216	*/
217	HistoType GetXOverYOverflow() {
218	return _xOverYOverflow;
219	}
220
221	/*! @brief Gets the X-Y underflow histogram.
222	*
223	* This single-bin histogram is built with events fall below the lower limits of both the X and Y axis.
224	*
225	* @return The X-Y overflow histogram.
226	*/
227	HistoType GetXUnderYUnderflow() {
228	return _xUnderYUnderflow;
229	}
230	/*! @brief Gets the (X over - Y under)flow histogram.
231	*
232	* This single-bin histogram is built with events fall above the lower limit of the Xaxis and below
233	* that of the Y axis.
234	*
235	* @return The (X over -Y under)flow histogram.
236	*/
237	HistoType GetXOverYUnderflow() {
238	return _xOverYUnderflow;
239	}
240
241	/*! @brief Gets the (X under - Y over)flow histogram.
242	*
243	* This single-bin histogram is built with events fall below the lower limit of the X axis and above
244	* that of the Y axis.
245	*
246	* @return The (X under -Y over)flow histogram.
247	*/
248	HistoType GetXUnderYOverflow() {
249	return _xUnderYOverflow;
250	}
251
252	protected:
253
254	/! @brief The vector containing the limits of the X bins(from lower to higher). /
255	std::vector<float> _xBins;
256	/! @brief The vector containing the limits of the Y bins(from lower to higher). /
257	std::vector<float> _yBins;
258	/! @brief A matrix containing the value of the histogram for each X-Y bin. /
259	SimpleMatrix<HistoType> _histo;
260	/! @brief The ROOT histogram. /
261	TH2 *_rootHisto;
262
263	/! @brief Base name of the output file. /
264	TString _outFileBase;
265	/! @brief Output file open mode (UPDATE or RECREATE, see documentation of TFile). /
266	TString _mode;
267	/! @brief Title for the ROOT histogram. /
268	TString _title;
269	/! @brief X axis label for the ROOT histogram. /
270	TString _xLabel;
271	/! @brief Y axis label for the ROOT histogram. /
272	TString _yLabel;
273
274	private:
275
276	vector<HistoType> _xUnderflow, _xOverflow, _yUnderflow, _yOverflow;
277	HistoType _xUnderYUnderflow, _xOverYOverflow, _xUnderYOverflow, _xOverYUnderflow;
278	bool _outRoot;
279	bool _outText;
280	void _CreateHisto();
281	void _InitHistos();
282	};
283
284	template<class HistoType>
285	void Histo2DAction<HistoType>::_CreateHisto() {
286
287	// No ROOT histogram for generic type; see template specializations in .cpp file.
288	_rootHisto = NULL;
289	}
290
291	// Specializations for _CreateHistos(). See Histo2DAction.cpp
292	template<>
293	void Histo2DAction<Int_t>::_CreateHisto();
294
295	template<>
296	void Histo2DAction<Float_t>::_CreateHisto();
297
298	template<>
299	void Histo2DAction<Double_t>::_CreateHisto();
300
301	template<class HistoType>
302	void Histo2DAction<HistoType>::_InitHistos() {
303
304	_CreateHisto();
305	if (_xBins.size() < 2) // SetXAxis not called by the main program, or wrongly filled (only 1 bin limit)
306	SetXAxis("Default X", 10, 0., 1.);
307	if (_yBins.size() < 2) // SetYAxis not called by the main program, or wrongly filled (only 1 bin limit)
308	SetYAxis("Default Y", 10, 0., 1.);
309
310	if (_rootHisto) {
311	Double_t *auxXArray = new Double_t[_xBins.size()];
312
313	for (unsigned int i = 0; i < _xBins.size(); i++) {
314	auxXArray[i] = _xBins[i];
315	}
316
317	Double_t *auxYArray = new Double_t[_yBins.size()];
318
319	for (unsigned int i = 0; i < _yBins.size(); i++) {
320	auxYArray[i] = _yBins[i];
321	}
322
323	_rootHisto->SetBins(_xBins.size() - 1, auxXArray, _yBins.size() - 1, auxYArray);
324	_rootHisto->SetName(GetName());
325	_rootHisto->SetTitle(_title);
326	_rootHisto->SetXTitle(_xLabel);
327	_rootHisto->SetYTitle(_yLabel);
328
329	delete[] auxXArray;
330	delete[] auxYArray;
331	}
332
333	/* The row index (first) corresponds to the position on the vertical (Y) axis. */
334	_histo.Resize(_yBins.size() - 1, _xBins.size() - 1);
335
336	_xUnderflow.resize(_yBins.size() - 1);
337	_xOverflow.resize(_yBins.size() - 1);
338	_yUnderflow.resize(_xBins.size() - 1);
339	_yOverflow.resize(_xBins.size() - 1);
340
341	}
342
343	template<class HistoType>
344	Histo2DAction<HistoType>::~Histo2DAction() {
345
346	delete _rootHisto;
347	_rootHisto = NULL;
348	}
349
350	template<class HistoType>
351	Histo2DAction<HistoType>::Histo2DAction(const char *actionName, TString title, TString outFileBase, TString mode,
352	bool outRoot, bool outText) :
353	CollectionAction(actionName), _xBins(0), _yBins(0), _histo(0, 0), _rootHisto(NULL), _outFileBase(outFileBase), _mode(
354	mode), _title(title), _xLabel(""), _yLabel(""), _xUnderflow(0), _xOverflow(0), _yUnderflow(0), _yOverflow(0),
355	_xUnderYUnderflow((HistoType) 0), _xOverYOverflow((HistoType) 0), _xUnderYOverflow((HistoType) 0),
356	_xOverYUnderflow((HistoType) 0), _outRoot(outRoot), _outText(outText) {
357
358	}
359
360	template<class HistoType>
361	void Histo2DAction<HistoType>::SetXAxis(TString label, vector<float> &bins) {
362
363	_xBins = bins;
364	_xLabel = label;
365
366	}
367
368	template<class HistoType>
369	void Histo2DAction<HistoType>::SetYAxis(TString label, vector<float> &bins) {
370
371	_yBins = bins;
372	_yLabel = label;
373
374	}
375
376	template<class HistoType>
377	void Histo2DAction<HistoType>::SetXAxis(TString label, TString binsFile) {
378
379	// Reading the bins from file
380	ifstream binListFile;
381	binListFile.open(binsFile);
382
383	TString auxString;
384	_xBins.resize(0);
385	while (!binListFile.eof()) {
386	binListFile >> auxString;
387	if (auxString != "") {
388	_xBins.push_back(auxString.Atof());
389	}
390	}
391	binListFile.close();
392
393	_xLabel = label;
394
395	}
396
397	template<class HistoType>
398	void Histo2DAction<HistoType>::SetYAxis(TString label, TString binsFile) {
399
400	// Reading the bins from file
401	ifstream binListFile;
402	binListFile.open(binsFile);
403
404	TString auxString;
405	_yBins.resize(0);
406	while (!binListFile.eof()) {
407	binListFile >> auxString;
408	if (auxString != "") {
409	_yBins.push_back(auxString.Atof());
410	}
411	}
412	binListFile.close();
413
414	_yLabel = label;
415
416	}
417
418	template<class HistoType>
419	void Histo2DAction<HistoType>::SetXAxis(TString label, unsigned int nBins, float min, float max, bool logBinning) {
420
421	_xBins.resize(nBins + 1);
422
423	if (!logBinning \|\| (logBinning && min <= 0.)) {
424
425	#ifdef DEBUGPAMCUT
426	if (logbinning && rigMin <= 0.)
427	cout << "Warning: logarithmic binning was chosen for X axis but min <= 0. Using linear binning."
428	#endif
429
430	float step = (max - min) / nBins;
431	for (unsigned int i = 0; i < nBins + 1; i++) {
432	_xBins[i] = min + i * step;
433	}
434
435	}
436	else {
437
438	double maxExp = log10(max / min);
439	for (unsigned int i = 0; i < nBins + 1; i++) {
440	_xBins[i] = min * pow(10., (double) i / ((double) nBins) * maxExp);
441	}
442
443	}
444
445	_xLabel = label;
446	}
447
448	template<class HistoType>
449	void Histo2DAction<HistoType>::SetYAxis(TString label, unsigned int nBins, float min, float max, bool logBinning) {
450
451	_yBins.resize(nBins + 1);
452
453	if (!logBinning \|\| (logBinning && min <= 0.)) {
454
455	#ifdef DEBUGPAMCUT
456	if (logbinning && rigMin <= 0.)
457	cout << "Warning: logarithmic binning was chosen for Y axis but min <= 0. Using linear binning."
458	#endif
459
460	float step = (max - min) / nBins;
461	for (unsigned int i = 0; i < nBins + 1; i++) {
462	_yBins[i] = min + i * step;
463	}
464
465	}
466	else {
467
468	double maxExp = log10(max / min);
469	for (unsigned int i = 0; i < nBins + 1; i++) {
470	_yBins[i] = min * pow(10., (double) i / ((double) nBins) * maxExp);
471	}
472
473	}
474
475	_yLabel = label;
476	}
477
478	template<class HistoType>
479	void Histo2DAction<HistoType>::Fill(double xValue, double yValue, double weight) {
480
481	_rootHisto->Fill(xValue, yValue, weight);
482
483	int xBin = 0, yBin = 0;
484
485	bool UOflow = false;
486
487	if (xValue < _xBins.front()) {
488	xBin = -1;
489	UOflow = true;
490	}
491
492	if (xValue >= _xBins.back()) {
493	xBin = _xBins.size();
494	UOflow = true;
495	}
496
497	if (yValue < _yBins.front()) {
498	yBin = -1;
499	UOflow = true;
500	}
501
502	if (yValue >= _yBins.back()) {
503	yBin = _yBins.size();
504	UOflow = true;
505	}
506
507	if (xBin == 0) {
508
509	while (xValue >= _xBins[xBin]) {
510	xBin++;
511	}
512	xBin--;
513	}
514
515	if (yBin == 0) {
516
517	while (yValue >= _yBins[yBin]) {
518	yBin++;
519	}
520	yBin--;
521	}
522
523	if (!UOflow) {
524	// No under or over flow.
525	_histo[yBin][xBin] += (HistoType) weight;
526	return;
527	}
528	else {
529	// An under or over flow has happened.
530	if (xBin == -1) {
531	if (yBin == -1) {
532	_xUnderYUnderflow += (HistoType) weight;
533	return;
534	}
535	else {
536	if (yBin == (int) _yBins.size()) {
537	_xUnderYOverflow += (HistoType) weight;
538	return;
539	}
540	else {
541	_xUnderflow[yBin] += (HistoType) weight;
542	return;
543	}
544	}
545	}
546
547	if (xBin == (int) _xBins.size()) {
548	if (yBin == -1) {
549	_xOverYUnderflow += (HistoType) weight;
550	return;
551	}
552	else {
553	if (yBin == (int) _yBins.size()) {
554	_xOverYOverflow += (HistoType) weight;
555	return;
556	}
557	else {
558	_xOverflow[yBin] += (HistoType) weight;
559	return;
560	}
561	}
562	}
563
564	if (yBin == -1) {
565	_yUnderflow[xBin] += (HistoType) weight;
566	return;
567	}
568
569	if (yBin == (int) _yBins.size()) {
570	_yOverflow[xBin] += (HistoType) weight;
571	return;
572	}
573
574	}
575	}
576
577	template<class HistoType>
578	void Histo2DAction<HistoType>::Finalize() {
579
580	if (_outFileBase != "") {
581	// Write the ROOT file
582	if (_rootHisto && _outRoot) {
583	TFile outRootFile((_outFileBase + ".root"), _mode);
584	outRootFile.cd();
585	_rootHisto->Write();
586	outRootFile.Close();
587	}
588
589	//Write the text file
590	if (_outText) {
591	ofstream outTextFile((_outFileBase + ".txt").Data(), ios_base::out);
592	for (unsigned int i = 0; i < _histo.GetNRows(); i++) {
593	for (unsigned int j = 0; j < _histo.GetNCols(); j++) {
594	outTextFile << setw(7) << _histo[i][j] << " ";
595	}
596	outTextFile << "\n";
597	}
598	outTextFile.close();
599	}
600	}
601
602	}
603	#endif /* HISTO2DACTION_H_ */