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;

  TString _outFileBase;
  TString _mode;
  TString _title, _xLabel, _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;
          cout << "_xUnderYOverflow: " << GetXUnderYOverflow() << ", weight: " << weight << ", (HistoType) weight: "
              << (HistoType) weight << endl;
          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	pam-fi	1.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	pam-fi	1.2	/*! @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	pam-fi	1.1	/*! @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	pam-fi	1.5	/*! @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	pam-fi	1.2	/*! Fills the ROOT and the vector histogram.
155			*
156	pam-fi	1.3	* @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	pam-fi	1.2	* @param weight The weight which will be applied to the event.
159			*/
160	pam-fi	1.1	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	pam-fi	1.2	/! @brief The vector containing the limits of the X bins(from lower to higher). /
255	pam-fi	1.1	std::vector<float> _xBins;
256	pam-fi	1.2	/! @brief The vector containing the limits of the Y bins(from lower to higher). /
257	pam-fi	1.1	std::vector<float> _yBins;
258	pam-fi	1.2	/! @brief A matrix containing the value of the histogram for each X-Y bin. /
259	pam-fi	1.1	SimpleMatrix<HistoType> _histo;
260	pam-fi	1.2	/! @brief The ROOT histogram. /
261	pam-fi	1.1	TH2 *_rootHisto;
262
263	pam-fi	1.6	TString _outFileBase;
264			TString _mode;
265			TString _title, _xLabel, _yLabel;
266
267	pam-fi	1.1	private:
268
269			vector<HistoType> _xUnderflow, _xOverflow, _yUnderflow, _yOverflow;
270			HistoType _xUnderYUnderflow, _xOverYOverflow, _xUnderYOverflow, _xOverYUnderflow;
271			bool _outRoot;
272			bool _outText;
273			void _CreateHisto();
274			void _InitHistos();
275			};
276
277			template<class HistoType>
278			void Histo2DAction<HistoType>::_CreateHisto() {
279
280			// No ROOT histogram for generic type; see template specializations in .cpp file.
281			_rootHisto = NULL;
282			}
283
284	pam-fi	1.4	// Specializations for _CreateHistos(). See Histo2DAction.cpp
285			template<>
286			void Histo2DAction<Int_t>::_CreateHisto();
287
288			template<>
289			void Histo2DAction<Float_t>::_CreateHisto();
290
291			template<>
292			void Histo2DAction<Double_t>::_CreateHisto();
293
294	pam-fi	1.1	template<class HistoType>
295			void Histo2DAction<HistoType>::_InitHistos() {
296
297			_CreateHisto();
298	pam-fi	1.5	if (_xBins.size() < 2) // SetXAxis not called by the main program, or wrongly filled (only 1 bin limit)
299			SetXAxis("Default X", 10, 0., 1.);
300			if (_yBins.size() < 2) // SetYAxis not called by the main program, or wrongly filled (only 1 bin limit)
301			SetYAxis("Default Y", 10, 0., 1.);
302	pam-fi	1.1
303			if (_rootHisto) {
304			Double_t *auxXArray = new Double_t[_xBins.size()];
305
306			for (unsigned int i = 0; i < _xBins.size(); i++) {
307			auxXArray[i] = _xBins[i];
308			}
309
310			Double_t *auxYArray = new Double_t[_yBins.size()];
311
312			for (unsigned int i = 0; i < _yBins.size(); i++) {
313			auxYArray[i] = _yBins[i];
314			}
315
316			_rootHisto->SetBins(_xBins.size() - 1, auxXArray, _yBins.size() - 1, auxYArray);
317			_rootHisto->SetName(GetName());
318			_rootHisto->SetTitle(_title);
319			_rootHisto->SetXTitle(_xLabel);
320			_rootHisto->SetYTitle(_yLabel);
321
322			delete[] auxXArray;
323			delete[] auxYArray;
324			}
325
326			/* The row index (first) corresponds to the position on the vertical (Y) axis. */
327			_histo.Resize(_yBins.size() - 1, _xBins.size() - 1);
328
329	pam-fi	1.6	_xUnderflow.resize(_yBins.size() - 1);
330			_xOverflow.resize(_yBins.size() - 1);
331			_yUnderflow.resize(_xBins.size() - 1);
332			_yOverflow.resize(_xBins.size() - 1);
333	pam-fi	1.1
334			}
335
336			template<class HistoType>
337			Histo2DAction<HistoType>::~Histo2DAction() {
338	pam-fi	1.3
339			delete _rootHisto;
340			_rootHisto = NULL;
341	pam-fi	1.1	}
342
343			template<class HistoType>
344			Histo2DAction<HistoType>::Histo2DAction(const char *actionName, TString title, TString outFileBase, TString mode,
345			bool outRoot, bool outText) :
346	pam-fi	1.6	CollectionAction(actionName), _xBins(0), _yBins(0), _histo(0, 0), _rootHisto(NULL), _outFileBase(outFileBase), _mode(
347			mode), _title(title), _xLabel(""), _yLabel(""), _xUnderflow(0), _xOverflow(0), _yUnderflow(0), _yOverflow(0),
348			_xUnderYUnderflow((HistoType) 0), _xOverYOverflow((HistoType) 0), _xUnderYOverflow((HistoType) 0),
349			_xOverYUnderflow((HistoType) 0), _outRoot(outRoot), _outText(outText) {
350	pam-fi	1.1
351			}
352
353			template<class HistoType>
354			void Histo2DAction<HistoType>::SetXAxis(TString label, vector<float> &bins) {
355
356			_xBins = bins;
357			_xLabel = label;
358
359			}
360
361			template<class HistoType>
362			void Histo2DAction<HistoType>::SetYAxis(TString label, vector<float> &bins) {
363
364			_yBins = bins;
365			_yLabel = label;
366
367			}
368
369			template<class HistoType>
370			void Histo2DAction<HistoType>::SetXAxis(TString label, TString binsFile) {
371
372			// Reading the bins from file
373			ifstream binListFile;
374			binListFile.open(binsFile);
375
376			TString auxString;
377			_xBins.resize(0);
378			while (!binListFile.eof()) {
379			binListFile >> auxString;
380			if (auxString != "") {
381			_xBins.push_back(auxString.Atof());
382			}
383			}
384			binListFile.close();
385
386			_xLabel = label;
387
388			}
389
390			template<class HistoType>
391			void Histo2DAction<HistoType>::SetYAxis(TString label, TString binsFile) {
392
393			// Reading the bins from file
394			ifstream binListFile;
395			binListFile.open(binsFile);
396
397			TString auxString;
398			_yBins.resize(0);
399			while (!binListFile.eof()) {
400			binListFile >> auxString;
401			if (auxString != "") {
402			_yBins.push_back(auxString.Atof());
403			}
404			}
405			binListFile.close();
406
407			_yLabel = label;
408
409			}
410
411			template<class HistoType>
412			void Histo2DAction<HistoType>::SetXAxis(TString label, unsigned int nBins, float min, float max, bool logBinning) {
413
414			_xBins.resize(nBins + 1);
415
416			if (!logBinning \|\| (logBinning && min <= 0.)) {
417
418			#ifdef DEBUGPAMCUT
419			if (logbinning && rigMin <= 0.)
420			cout << "Warning: logarithmic binning was chosen for X axis but min <= 0. Using linear binning."
421			#endif
422
423			float step = (max - min) / nBins;
424			for (unsigned int i = 0; i < nBins + 1; i++) {
425			_xBins[i] = min + i * step;
426			}
427
428			}
429			else {
430
431			double maxExp = log10(max / min);
432			for (unsigned int i = 0; i < nBins + 1; i++) {
433			_xBins[i] = min * pow(10., (double) i / ((double) nBins) * maxExp);
434			}
435
436			}
437
438			_xLabel = label;
439			}
440
441			template<class HistoType>
442			void Histo2DAction<HistoType>::SetYAxis(TString label, unsigned int nBins, float min, float max, bool logBinning) {
443
444			_yBins.resize(nBins + 1);
445
446			if (!logBinning \|\| (logBinning && min <= 0.)) {
447
448			#ifdef DEBUGPAMCUT
449			if (logbinning && rigMin <= 0.)
450			cout << "Warning: logarithmic binning was chosen for Y axis but min <= 0. Using linear binning."
451			#endif
452
453			float step = (max - min) / nBins;
454			for (unsigned int i = 0; i < nBins + 1; i++) {
455			_yBins[i] = min + i * step;
456			}
457
458			}
459			else {
460
461			double maxExp = log10(max / min);
462			for (unsigned int i = 0; i < nBins + 1; i++) {
463			_yBins[i] = min * pow(10., (double) i / ((double) nBins) * maxExp);
464			}
465
466			}
467
468			_yLabel = label;
469			}
470
471			template<class HistoType>
472	pam-fi	1.5	void Histo2DAction<HistoType>::Fill(double xValue, double yValue, double weight) {
473	pam-fi	1.1
474			_rootHisto->Fill(xValue, yValue, weight);
475
476			int xBin = 0, yBin = 0;
477
478			bool UOflow = false;
479
480			if (xValue < _xBins.front()) {
481			xBin = -1;
482			UOflow = true;
483			}
484
485	pam-fi	1.5	if (xValue >= _xBins.back()) {
486	pam-fi	1.1	xBin = _xBins.size();
487			UOflow = true;
488			}
489
490			if (yValue < _yBins.front()) {
491			yBin = -1;
492			UOflow = true;
493			}
494
495	pam-fi	1.5	if (yValue >= _yBins.back()) {
496	pam-fi	1.1	yBin = _yBins.size();
497			UOflow = true;
498			}
499
500			if (xBin == 0) {
501
502			while (xValue >= _xBins[xBin]) {
503			xBin++;
504			}
505			xBin--;
506			}
507
508			if (yBin == 0) {
509
510			while (yValue >= _yBins[yBin]) {
511			yBin++;
512			}
513			yBin--;
514			}
515
516			if (!UOflow) {
517			// No under or over flow.
518			_histo[yBin][xBin] += (HistoType) weight;
519			return;
520			}
521			else {
522			// An under or over flow has happened.
523			if (xBin == -1) {
524			if (yBin == -1) {
525			_xUnderYUnderflow += (HistoType) weight;
526			return;
527			}
528			else {
529	pam-fi	1.2	if (yBin == (int) _yBins.size()) {
530	pam-fi	1.1	_xUnderYOverflow += (HistoType) weight;
531	pam-fi	1.6	cout << "_xUnderYOverflow: " << GetXUnderYOverflow() << ", weight: " << weight << ", (HistoType) weight: "
532			<< (HistoType) weight << endl;
533	pam-fi	1.1	return;
534			}
535			else {
536			_xUnderflow[yBin] += (HistoType) weight;
537			return;
538			}
539			}
540			}
541
542	pam-fi	1.2	if (xBin == (int) _xBins.size()) {
543	pam-fi	1.1	if (yBin == -1) {
544			_xOverYUnderflow += (HistoType) weight;
545			return;
546			}
547			else {
548	pam-fi	1.2	if (yBin == (int) _yBins.size()) {
549	pam-fi	1.1	_xOverYOverflow += (HistoType) weight;
550			return;
551			}
552			else {
553			_xOverflow[yBin] += (HistoType) weight;
554			return;
555			}
556			}
557			}
558
559			if (yBin == -1) {
560			_yUnderflow[xBin] += (HistoType) weight;
561			return;
562			}
563
564	pam-fi	1.2	if (yBin == (int) _yBins.size()) {
565	pam-fi	1.1	_yOverflow[xBin] += (HistoType) weight;
566			return;
567			}
568
569			}
570			}
571
572			template<class HistoType>
573			void Histo2DAction<HistoType>::Finalize() {
574
575			if (_outFileBase != "") {
576			// Write the ROOT file
577			if (_rootHisto && _outRoot) {
578			TFile outRootFile((_outFileBase + ".root"), _mode);
579			outRootFile.cd();
580			_rootHisto->Write();
581			outRootFile.Close();
582			}
583
584			//Write the text file
585			if (_outText) {
586			ofstream outTextFile((_outFileBase + ".txt").Data(), ios_base::out);
587			for (unsigned int i = 0; i < _histo.GetNRows(); i++) {
588			for (unsigned int j = 0; j < _histo.GetNCols(); j++) {
589	pam-fi	1.6	outTextFile << setw(7) << _histo[i][j] << " ";
590	pam-fi	1.1	}
591			outTextFile << "\n";
592			}
593			outTextFile.close();
594			}
595			}
596
597			}
598			#endif /* HISTO2DACTION_H_ */