#ifndef __JRECONSTRUCTION__JEVT__
#define __JRECONSTRUCTION__JEVT__
#include <istream>
#include <ostream>
#include <iomanip>
#include <limits>
#include <vector>
#include <algorithm>
#include <cmath>
#include <TROOT.h>
#include <TObject.h>
#include "JLang/JType.hh"
#include "JROOT/JRoot.hh"
#include "JROOT/JTreeParameters.hh"
#include "JReconstruction/JFitStatus.hh"
#include "JReconstruction/JHistory.hh"
/**
* \author mdejong
*/
namespace JRECONSTRUCTION {}
namespace JPP { using namespace JRECONSTRUCTION; }
namespace JFIT {
/**
* Data structure for track fit results with history and optional associated values.
*/
class JFit :
public TObject,
public JHistory
{
public:
/**
* Default constructor.
*
* Parameters are initialized with non physical values.
*/
JFit():
__x(std::numeric_limits<double>::lowest()),
__y(std::numeric_limits<double>::lowest()),
__z(std::numeric_limits<double>::lowest()),
__dx(0.0),
__dy(0.0),
__dz(0.0),
__t(std::numeric_limits<double>::lowest()),
__Q(std::numeric_limits<double>::lowest()),
__NDF(-1),
__E(0.0),
__status(JRECONSTRUCTION::UNDEFINED)
{}
/**
* Constructor.
*
* \param history history
* \param x X-position
* \param y Y-position
* \param z Z-position
* \param dx X-slope
* \param dy Y-slope
* \param dz Z-slope
* \param t time
* \param Q quality
* \param NDF number of degrees of freedom
* \param E energy
* \param status status
*/
JFit(const JHistory& history,
const double x,
const double y,
const double z,
const double dx,
const double dy,
const double dz,
const double t,
const double Q,
const int NDF,
const double E = 0.0,
const int status = JRECONSTRUCTION::SINGLE_STAGE) :
JHistory(history)
{
__x = x;
__y = y;
__z = z;
__dx = dx;
__dy = dy;
__dz = dz;
__t = t;
__Q = Q;
__NDF = NDF;
__E = E;
__status = status;
}
/**
* Constructor for storing position only.
*
* Note that the type of fit can be obtained via the history information.
*
* \param history history
* \param x X-position
* \param y Y-position
* \param z Z-position
* \param status status
*/
JFit(const JHistory& history,
const double x,
const double y,
const double z,
const int status = JRECONSTRUCTION::SINGLE_STAGE):
JHistory(history)
{
__x = x;
__y = y;
__z = z;
__dx = 0.0;
__dy = 0.0;
__dz = 0.0;
__t = 0.0;
__Q = 0.0;
__NDF = -1;
__E = 0.0;
__status = status;
}
/**
* Add event to history.
*
* \param type application type
* \return this fit
*/
JFit& add(const int type)
{
getHistory().add(type);
return *this;
}
double getX() const { return __x; } //!< Get X-position.
double getY() const { return __y; } //!< Get Y-position.
double getZ() const { return __z; } //!< Get Z-position.
double getDX() const { return __dx; } //!< Get X-slope.
double getDY() const { return __dy; } //!< Get Y-slope.
double getDZ() const { return __dz; } //!< Get Z-slope.
double getT() const { return __t; } //!< Get time.
double getQ() const { return __Q; } //!< Get quality.
int getNDF() const { return __NDF; } //!< Get number of degrees of freedom.
double getE() const { return __E; } //!< Get energy.
int getStatus() const { return __status; } //!< Get status of the fit; negative values should refer to a bad fit.
/**
* Set status of the fit.
*
* \param status status
*/
void setStatus(const int status)
{
this->__status = status;
}
/**
* Move vertex along this track with given velocity.
*
* \param step step
* \param velocity velocity
*/
void move(const double step, const double velocity)
{
__x += step * __dx;
__y += step * __dy;
__z += step * __dz;
__t += step / velocity;
}
/**
* Set energy.
*
* \param E energy
*/
void setE(const double E)
{
__E = E;
}
/**
* Get dimension of error matrix.
*
* \return dimension
*/
size_t getDimensionOfErrorMatrix() const
{
return ((size_t) sqrt(1.0 + 8.0*V.size()) - 1) / 2;
}
/**
* Get error matrix.
*
* Note that only the lower-half of the matrix is returned.
*
* \return error matrix
*/
const std::vector<double>& getV() const
{
return V;
}
/**
* Get element of error matrix.
*
* \param row row (<tt>0 <= row < dimension</tt>)
* \param col col (<tt>0 <= col < dimension</tt>)
* \return value
*/
double getV(const size_t row, const size_t col) const
{
using namespace std;
const size_t __row = max(row, col) + 1;
const size_t __col = min(row, col);
return V.at((__row * (__row + 1)) / 2 - __row + __col);
}
/**
* Set error matrix.
*
* The given size corresponds to the dimension of a 2D-array of which
* the elements should be accessible via the usual array operators.\n
* Note that only the lower-half of the given matrix is stored.
*
* \param size size
* \param data matrix
*/
template<class T>
void setV(const size_t size, const T& data)
{
V.clear();
for (size_t row = 0; row != size; ++row) {
for (size_t col = 0; col <= row; ++col) {
V.push_back(data[row][col]);
}
}
}
/**
* Get associated values.
*
* \return values
*/
const std::vector<double>& getW() const
{
return this->W;
}
/**
* Set associated values.
*
* \param W values
*/
void setW(const std::vector<double>& W)
{
this->W = W;
}
/**
* Get number of associated values.
*
* \return number of values
*/
int getN() const
{
return W.size();
}
/**
* Check availability of value.
*
* \param i index
* \return true if available; else false
*/
bool hasW(const int i) const
{
return (i >= 0 && i < (int) W.size());
}
/**
* Get associated value.
*
* \param i index
* \return value
*/
double getW(const int i) const
{
return W.at(i);
}
/**
* Get associated value.
*
* \param i index
* \param value default value
* \return value
*/
double getW(const int i, const double value) const
{
if (hasW(i))
return W.at(i);
else
return value;
}
/**
* Set associated value.
*
* \param i index
* \param value value
*/
void setW(const int i, const double value)
{
if (i >= (int) W.size()) {
W.resize(i + 1, 0.0);
}
W[i] = value;
}
ClassDef(JFit, 7);
protected:
double __x;
double __y;
double __z;
double __dx;
double __dy;
double __dz;
double __t;
double __Q;
int __NDF;
std::vector<double> V;
std::vector<double> W;
double __E;
int __status;
};
/**
* Data structure for set of track fit results.
*/
class JEvt :
public TObject,
public std::vector<JFit>
{
public:
/**
* Default constructor.
*/
JEvt()
{}
/**
* Select fits.
*
* \param selector fit selector
*/
template<class JSelector_t>
void select(const JSelector_t& selector)
{
using namespace std;
if (!empty()) {
iterator __end = partition(this->begin(), this->end(), selector);
this->erase(__end, this->end());
}
}
/**
* Select fits.
*
* \param number_of_fits maximal number of best fits to select
* \param comparator quality comparator
*/
template<class JComparator_t>
void select(const size_t number_of_fits,
const JComparator_t& comparator)
{
using namespace std;
if (!empty()) {
iterator __end = this->end();
if (number_of_fits > 0 && number_of_fits < this->size()) {
advance(__end = this->begin(), number_of_fits);
partial_sort(this->begin(), __end, this->end(), comparator);
} else {
sort(this->begin(), __end, comparator);
}
this->erase(__end, this->end());
}
}
/**
* Write event to output.
*
* \param out output stream
* \param event event
* \return output stream
*/
friend inline std::ostream& operator<<(std::ostream& out, const JEvt& event)
{
using namespace std;
out << "Event: " << endl;
for (JEvt::const_iterator fit = event.begin(); fit != event.end(); ++fit) {
out << *fit;
}
return out;
}
ClassDef(JEvt, 5);
};
}
namespace JRECONSTRUCTION {
typedef JFIT::JFit JFit;
typedef JFIT::JEvt JEvt;
}
/**
* Write fit results to output.
*
* \param out output stream
* \param fit fit results
* \return output stream
*/
std::ostream& operator<<(std::ostream& out, const JRECONSTRUCTION::JFit& fit);
/**
* Get TTree parameters for given data type.
*
* \return TTree parameters
*/
inline JROOT::JTreeParameters getTreeParameters(JLANG::JType<JRECONSTRUCTION::JEvt>)
{
return JROOT::JTreeParameters("EVT", "evt", "", 0);
}
#endif