-
Pere Mato Vila authoredPere Mato Vila authored
TrackMathCore.h 7.78 KiB
#ifndef INCLUDE_TRACKMATHCORE_H
#define INCLUDE_TRACKMATHCORE_H
// dummy track class for testing I/o of matric
//#ifdef USE_REFLEX
//#define G__DICTIONARY
//#endif
#include "Math/Point3D.h"
#include "Math/Vector4D.h"
#include "Math/SMatrix.h"
#include "Rtypes.h" // for Double32_t
#include "TError.h"
#include <vector>
#include <string>
#include <iostream>
#include <cassert>
typedef ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > Vector4D;
typedef ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<Double32_t> > Vector4D32;
typedef ROOT::Math::DisplacementVector3D<ROOT::Math::Cartesian3D<double> > Vector3D;
typedef ROOT::Math::PositionVector3D<ROOT::Math::Cartesian3D<double> > Point3D;
typedef ROOT::Math::PositionVector3D<ROOT::Math::Cartesian3D<Double32_t> > Point3D32;
typedef ROOT::Math::SMatrix<double,4,4,ROOT::Math::MatRepStd<double,4,4> > Matrix4D;
typedef ROOT::Math::SMatrix<Double32_t,4,4,ROOT::Math::MatRepStd<Double32_t,4,4> > Matrix4D32;
typedef ROOT::Math::SMatrix<double,6,6,ROOT::Math::MatRepSym<double,6> > SymMatrix6D;
typedef ROOT::Math::SMatrix<Double32_t,6,6,ROOT::Math::MatRepSym<Double32_t,6> > SymMatrix6D32;;
// track class containing a vector and a point
class TrackD {
public:
typedef Matrix4D::const_iterator const_iterator;
TrackD() {}
TrackD(double * begin, double * end) {
double * itr = begin;
fPos.SetCoordinates(itr, itr+3); itr +=3;
fVec.SetCoordinates(itr,itr+4); itr+=4;
R__ASSERT(itr == end);
}
enum { kSize = 3 + 4 };
static std::string Type() { return "TrackD"; }
static bool IsD32() { return false; }
TrackD & operator += (const TrackD & t) {
fPos += Vector3D(t.fPos);
fVec += t.fVec;
return *this;
}
double Sum() const {
double s = 0;
double d[4];
fPos.GetCoordinates(d,d+3); for (int i=0; i<3;++i)
s+= d[i];
fVec.GetCoordinates(d,d+4);
for (int i=0; i<4;++i)
s+= d[i];
return s;
}
void Print() const {
std::cout << "Point " << fPos << std::endl;
std::cout << "Vec " << fVec << std::endl;
}
private:
Point3D fPos;
Vector4D fVec;
};
// track class based on of Double32
class TrackD32 {
public:
typedef Matrix4D::const_iterator const_iterator;
TrackD32() {}
enum { kSize = 3 + 4 };
static std::string Type() { return "TrackD32"; }
static bool IsD32() { return true; }
TrackD32(double * begin, double * end) {
double * itr = begin;
fPos.SetCoordinates(itr, itr+3); itr +=3;
fVec.SetCoordinates(itr,itr+4); itr+=4;
R__ASSERT(itr == end);
}
TrackD32 & operator += (const TrackD32 & t) {
fPos += Vector3D(t.fPos);
fVec += t.fVec;
return *this;
}
double Sum() const {
double s = 0;
double d[4];
fPos.GetCoordinates(d,d+3);
for (int i=0; i<3;++i)
s+= d[i];
fVec.GetCoordinates(d,d+4);
for (int i=0; i<4;++i)
s+= d[i];
return s;
}
void Print() const {
std::cout << "Point " << fPos << std::endl;
std::cout << "Vec " << fVec << std::endl;
}
private:
Point3D32 fPos;
Vector4D32 fVec;
};
// track class (containing a vector and a point) and matrices
class TrackErrD {
public:
typedef Matrix4D::const_iterator const_iterator;
TrackErrD() {}
TrackErrD(double * begin, double * end) {
double * itr = begin;
fPos.SetCoordinates(itr, itr+3); itr +=3;
fVec.SetCoordinates(itr,itr+4); itr+=4;
fMat = Matrix4D(itr,itr+16); itr += 16;
fSymMat = SymMatrix6D(itr,itr+21);
R__ASSERT(itr+21 == end);
}
enum { kSize =
3 + 4 + Matrix4D32::kSize + SymMatrix6D32::rep_type::kSize
};
static std::string Type() { return "TrackErrD"; }
static bool IsD32() { return false; }
TrackErrD & operator += (const TrackErrD & t) {
fPos += Vector3D(t.fPos);
fVec += t.fVec;
fMat += t.fMat;
fSymMat += t.fSymMat;
return *this;
}
double Sum() const {
double s = 0;
double d[4];
fPos.GetCoordinates(d,d+3);
for (int i=0; i<3;++i)
s+= d[i];
fVec.GetCoordinates(d,d+4);
for (int i=0; i<4;++i)
s+= d[i];
for (const_iterator itr = fMat.begin(); itr != fMat.end(); ++itr)
s+= *itr;
for (const_iterator itr = fSymMat.begin(); itr != fSymMat.end(); ++itr)
s+= *itr;
return s;
}
void Print() const {
std::cout << "Point " << fPos << std::endl;
std::cout << "Vec " << fVec << std::endl;
std::cout << "Mat " << fMat << std::endl;
std::cout << "SymMat " << fSymMat << std::endl;
}
private:
Point3D fPos;
Vector4D fVec; Matrix4D fMat;
SymMatrix6D fSymMat;
};
// track class based on of Double32
class TrackErrD32 {
public:
typedef Matrix4D::const_iterator const_iterator;
TrackErrD32() {}
enum { kSize =
3 + 4 + Matrix4D32::kSize + SymMatrix6D32::rep_type::kSize
};
static std::string Type() { return "TrackErrD32"; }
static bool IsD32() { return true; }
TrackErrD32(double * begin, double * end) {
double * itr = begin;
fPos.SetCoordinates(itr, itr+3); itr +=3;
fVec.SetCoordinates(itr,itr+4); itr+=4;
fMat = Matrix4D32(itr,itr+16); itr += 16;
fSymMat = SymMatrix6D32(itr,itr+21);
R__ASSERT(itr+21 == end);
}
TrackErrD32 & operator += (const TrackErrD32 & t) {
fPos += Vector3D(t.fPos);
fVec += t.fVec;
fMat += t.fMat;
fSymMat += t.fSymMat;
return *this;
}
double Sum() const {
double s = 0;
double d[4];
fPos.GetCoordinates(d,d+3);
for (int i=0; i<3;++i)
s+= d[i];
fVec.GetCoordinates(d,d+4);
for (int i=0; i<4;++i)
s+= d[i];
for (const_iterator itr = fMat.begin(); itr != fMat.end(); ++itr)
s+= *itr;
for (const_iterator itr = fSymMat.begin(); itr != fSymMat.end(); ++itr)
s+= *itr;
return s;
}
void Print() const {
std::cout << "Point " << fPos << std::endl;
std::cout << "Vec " << fVec << std::endl;
std::cout << "Mat " << fMat << std::endl;
std::cout << "SymMat " << fSymMat << std::endl;
}
private:
Point3D32 fPos;
Vector4D32 fVec; Matrix4D32 fMat;
SymMatrix6D32 fSymMat;
};
// class containning a vector of tracks
template<class T>
class VecTrack {
public:
typedef typename std::vector<T>::const_iterator It;
VecTrack() {
// create klen empty trackD
fTrks.reserve(kLen);
for (int i = 0; i < kLen; ++i) {
fTrks.push_back(T() );
}
}
VecTrack(double * ibegin, double * iend)
{
fTrks.reserve(kLen);
double * itr = ibegin;
for (int i = 0; i < kLen; ++i) {
fTrks.push_back(T(itr, itr + T::kSize) );
itr += T::kSize;
}
R__ASSERT( itr == iend);
}
enum { kLen = 3, kSize = kLen*T::kSize };
static std::string Type() { return "VecTrack<"+T::Type()+">"; }
static bool IsD32() { return false; }
VecTrack & operator += (const VecTrack & v) {
for (unsigned int i = 0; i < fTrks.size() ; ++i)
fTrks[i] += v.fTrks[i];
return *this;
}
It begin() const { return fTrks.begin(); }
It end() const { return fTrks.end(); }
double Sum() const {
double s = 0;
for (unsigned int i = 0; i < fTrks.size() ; ++i)
s += fTrks[i].Sum();
return s;
}
void Print() const {
for (unsigned int i = 0; i < fTrks.size() ; ++i) {
std::cout << "\n======> Track ========<" << i << std::endl;
fTrks[i].Print();
}
}
private:
std::vector<T> fTrks;
};
// for instantiating the template VecTrackD class for reflex
struct Dummy {
VecTrack<TrackD> v1;
VecTrack<TrackErrD> v2;
};
#endif