reconstruction

Topics
	reconstruction data objects
	reconstruction modules

Classes
class	BeamParametersFitter
	Fitter calculating BeamParameters from CollisionBoostVector and CollisionInvariantMass. More...
class	BeamSpotAlgorithm
	Class implementing BeamSpot calibration algorithm. More...
class	BoostVectorAlgorithm
	Class implementing BoostVector calibration algorithm. More...
struct	CalibPars
	The parameters related to single calibration interval. More...
struct	CalibrationData
	Parameters and data relevant for single calibration interval. More...
class	ChebFitter
	Unbinned Maximum Likelihood fitter with a possibility to use Chebyshev interpolation. More...
class	InvariantMassAlgorithm
	Class implementing InvariantMass calibration algorithm. More...
struct	Atom
	Very small (few mins) calibration interval which cannot be further divided : Atom. More...
struct	ExpRun
	Struct containing exp number and run number. More...
struct	ExpRunEvt
	struct with expNum, runNum, evtNum More...
class	Splitter
	Class that allows to split runs into the intervals of intended properties given by the lossFunction. More...
struct	Spline
	Spline structure for zero-order & linear splines. More...
class	EventsOfDoomParameters
	DBObject containing parameters used in EventsOfDoomBuster module. More...
class	RootMergeable< T >
	Wrap a root histogram or TNtuple to make it mergeable. More...

Typedefs
typedef std::map< std::string, double >	Pars
	values of parameters in ML fit
typedef std::map< std::string, std::pair< double, double > >	Limits
	limits of parameters in ML fit

Functions
TMatrixDSym	toTMatrixDSym (Eigen::MatrixXd mIn)
	Function that converts Eigen symmetric matrix to ROOT matrix.
B2Vector3D	toB2Vector3 (Eigen::VectorXd vIn)
	Function that converts Eigen vector to ROOT vector.
int	getID (const std::vector< double > &breaks, double t)
	get id of the time point t
void	extrapolateCalibration (std::vector< CalibrationData > &calVec)
	Extrapolate calibration to intervals where it failed.
void	addShortRun (std::vector< CalibrationData > &calVec, std::pair< ExpRun, std::pair< double, double > > shortRun)
	Extrapolate calibration to the very short runs which were filtered before.
double	encodeNumber (double val, unsigned num)
	Encode integer num into double val such that val is nearly not changed (maximally by a relative shift 1e-6).
unsigned	decodeNumber (double val)
	Decode the integer number encoded in val.
template<typename Evt>
void	storePayloads (const std::vector< Evt > &evts, const std::vector< CalibrationData > &calVecConst, const std::string &objName, std::function< TObject *(Eigen::VectorXd, Eigen::MatrixXd, Eigen::MatrixXd) > getCalibObj)
	Store payloads to files.
void	storePayloadsNoIntraRun (const std::vector< CalibrationData > &calVecConst, const std::string &objName, std::function< TObject *(Eigen::VectorXd, Eigen::MatrixXd, Eigen::MatrixXd) > getCalibObj)
	Store payloads to files, where calib data have no intra-run dependence.
template<typename Evt, typename Fun>
CalibrationData	runAlgorithm (const std::vector< Evt > &evts, std::vector< std::map< ExpRun, std::pair< double, double > > > range, Fun runCalibAnalysis)
	run calibration algorithm for single calibration interval
template<typename Fun1, typename Fun2>
CalibrationAlgorithm::EResult	runCalibration (TTree *tracks, const std::string &calibName, Fun1 GetEvents, Fun2 calibAnalysis, std::function< TObject *(Eigen::VectorXd, Eigen::MatrixXd, Eigen::MatrixXd)> calibObjCreator, TString m_lossFunctionOuter, TString m_lossFunctionInner)
	Run the the calibration over the whole event sample.
std::pair< double, double >	getMinima (const std::vector< std::vector< double > > &vals, int i0, int j0)
	Get minimum inside and outside of the smaller window defined by i0, j0.
std::vector< double >	getMinimum (std::function< double(double, double)> fun, double xMin, double xMax, double yMin, double yMax)
	Get minimum of 2D function in the rectangular domain defined by xMin,xMax & yMin,yMax.
Eigen::VectorXd	getWeights (int Size)
	Get the vector of weights to calculate the integral over the Chebyshev nodes The nodes are by definition between 0 and 1, there are Size nodes To get their positions, use getNodes.
Eigen::VectorXd	getNodes (int Size)
	Get the vector of positions of the Chebyshev nodes The nodes are by definition between 0 and 1, there are Size nodes For the corresponding weights use getWeights.
Eigen::VectorXd	getPols (int Size, double x)
	Evaluate Chebyshev polynomials up to Size at point x It returns a vector of the P_i(x) for i=0..Size-1 The polynomial is defined for x between 0 and 1.
Eigen::VectorXd	getPolsSum (int Size, Eigen::VectorXd x)
	Calculate the Chebyshev polynomials of order i=0..Size-1 at points given in vector x_j and sum it over point index j It returns sum_j P_i(x_j) for i=0..Size-1 The Chebyshev polynomials are defined for x between 0 and 1.
Eigen::MatrixXd	getCoefs (int Size, bool isInverse=false)
	Transformation matrix between Cheb nodes and coefficients of the Cheb polynomials Notice, that there are two alternative ways defining polynomial interpolation:
double	evalPol (const Eigen::VectorXd &polCoef, double x)
	Evaluate Cheb.
Eigen::MatrixXd	getCoefsCheb (int oldSize)
	Transformation matrix between Cheb nodes and Cheb coefficients with better normalization of the borders.
Eigen::VectorXd	interpol (const Eigen::VectorXd &xi, double x)
	Get Interpolation vector k_i for point x from the function values at points xi (polynomial interpolation) In the second step, the function value at x can be evaluated as sum_i vals_i k_i.
double	interpol (Eigen::VectorXd xi, Eigen::VectorXd vals, double x)
	Get interpolated function value at point x when function values vals at points xi are provided.
bool	operator!= (ExpRun a, ExpRun b)
	Not equal for ExpRun.
bool	operator< (ExpRun a, ExpRun b)
	less than for ExpRun
std::map< ExpRun, std::pair< double, double > >	filter (const std::map< ExpRun, std::pair< double, double > > &runs, double cut, std::map< ExpRun, std::pair< double, double > > &runsRemoved)
	filter events to remove runs shorter than cut, it stores removed runs in runsRemoved
std::pair< int, int >	getStartEndIndexes (int nIntervals, const std::vector< int > &breaks, int indx)
	get the range of interval with nIntervals and breaks stored in a vector
std::vector< Atom >	slice (std::vector< Atom > vec, int s, int e)
	Slice the vector to contain only elements with indexes s .. e (included)
std::vector< std::map< ExpRun, std::pair< double, double > > >	breaks2intervalsSep (const std::map< ExpRun, std::pair< double, double > > &runsMap, const std::vector< Atom > &currVec, const std::vector< int > &breaks)
	Get calibration intervals according to the indexes of the breaks.
template<typename Evt>
std::map< ExpRun, std::pair< double, double > >	getRunInfo (const std::vector< Evt > &evts)
	Get the map of runs, where each run contains pair with start/end time [hours].
template<typename Evt>
ExpRunEvt	getPosition (const std::vector< Evt > &events, double tEdge)
	Get the exp-run-evt number from the event time [hours].
template<typename Evt>
std::vector< ExpRunEvt >	convertSplitPoints (const std::vector< Evt > &events, const std::vector< double > &splitPoints)
	Convert splitPoints [hours] to breakPoints in ExpRunEvt.
TString	rn ()
	Get random string.
std::vector< std::vector< double > >	merge (const std::vector< std::vector< std::vector< double > > > &toMerge)
	merge { vector<double> a, vector<double> b} into {a, b}
Eigen::VectorXd	vec2vec (const std::vector< double > &vec)
	std vector -> ROOT vector
std::vector< double >	vec2vec (const Eigen::VectorXd &v)
	ROOT vector -> std vector.
Eigen::MatrixXd	vecs2mat (const std::vector< std::vector< double > > &vecs)
	merge columns (from std::vectors) into ROOT matrix
std::vector< double >	getRangeLin (int nVals, double xMin, double xMax)
	Equidistant range between xMin and xMax for spline of the first order.
std::vector< double >	getRangeZero (int nVals, double xMin, double xMax)
	Equidistant range between xMin and xMax for spline of the zero order.
std::vector< double >	slice (const std::vector< double > &v, unsigned ind, unsigned n)
	put slice of original vector v[ind:ind+n] into new one, n is number of elements
double	eval (const std::vector< double > &spl, const std::vector< double > &vals, double x)
	Evaluate spline (zero order or first order) in point x.
VectorXd	getPolsSum (int Size, VectorXd x)
	Calculate the Chebyshev polynomials of order i=0..Size-1 at points given in vector x_j and sum it over point index j It returns sum_j P_i(x_j) for i=0..Size-1 The Chebyshev polynomials are defined for x between 0 and 1.
double	evalPol (const VectorXd &polCoef, double x)
	Evaluate Cheb.
VectorXd	interpol (const VectorXd &xi, double x)
	Get Interpolation vector k_i for point x from the function values at points xi (polynomial interpolation) In the second step, the function value at x can be evaluated as sum_i vals_i k_i.
double	interpol (VectorXd xi, VectorXd vals, double x)
	Get interpolated function value at point x when function values vals at points xi are provided.
void	plotRuns (const std::vector< std::pair< double, double > > &runs)
	plot runs on time axis
void	plotSRuns (const std::vector< std::pair< double, double > > &runs, const std::vector< int > &breaks, int offset=2)
	plot clusters or runs on time axis
void	printBySize (const std::vector< std::pair< double, double > > &runs)
	print sorted lengths of the runs
static ExpRun	getRun (const std::map< ExpRun, std::pair< double, double > > &runs, double t)
	Get exp number + run number from time.
VXDDedxTrack const *	getSVDDedxFromParticle (Particle const *particle)
	SVD dEdx value from particle.
Eigen::VectorXd	getLogFunction (Pars pars) const
	Get the -2*log(p(x)) on the Cheb nodes.
void	init (int Size, double xMin, double xMax)
	Initialize the fitter (the Chebyshev coefficients)
double	getLogLikelihoodSlow (const Pars &pars) const
	Calculate log likelihood using exact formula.
double	getLogLikelihoodFast (const Pars &pars) const
	Calculate log likelihood using approximation based on Chebyshev polynomials (typically faster)
double	operator() (const double *par) const
	Evaluate the log likelihood.
Eigen::VectorXd	getDataGrid () const
	Calculate Chebyshev coefficients for the data set.
std::pair< Pars, Eigen::MatrixXd >	fitData (const Pars &pars, Limits limits, bool UseCheb=true)
	Fit the data with specified initial values of parameters and limits on them.
double	lossFunction (const std::vector< Atom > &vec, int s, int e) const
	lossFunction of the calibration interval consisting of several "atoms" stored in vector vec The atoms included in calibration interval have indices between s and e
static std::vector< std::pair< double, double > >	splitToSmall (const std::map< ExpRun, std::pair< double, double > > &runs, double intSize=1./60)
	Split the runs into small calibration intervals (atoms) of a specified size.
double	getMinLoss (const std::vector< Atom > &vec, int e, std::vector< int > &breaks)
	Recursive function to evaluate minimal sum of the lossFuctions for the optimal clustering.
std::vector< int >	dynamicBreaks (const std::vector< Atom > &runs)
	Get optimal break points using algorithm based on dynamic programming.
static std::map< ExpRun, std::pair< double, double > >	mergeIntervals (const std::map< ExpRun, std::pair< double, double > > &I1, const std::map< ExpRun, std::pair< double, double > > &I2)
	Merge two subintervals into one subinterval.

Detailed Description

Typedef Documentation

Limits

typedef std::map<std::string, std::pair<double, double> > Limits

limits of parameters in ML fit

Definition at line 28 of file ChebFitter.h.

Pars

typedef std::map<std::string, double> Pars

values of parameters in ML fit

Definition at line 25 of file ChebFitter.h.

Function Documentation

addShortRun()

void addShortRun ( std::vector< CalibrationData > & calVec, std::pair< ExpRun, std::pair< double, double > > shortRun )

inline

Extrapolate calibration to the very short runs which were filtered before.

Definition at line 151 of file calibTools.h.

  {
    double shortStart = shortRun.second.first;
    double shortEnd   = shortRun.second.second;
 
    double distMin = 1e20;
    int iMin = -1, jMin = -1;
 
    for (unsigned i = 0; i < calVec.size(); ++i) {
      if (calVec[i].isCalibrated == false)
        continue;
      for (unsigned j = 0; j < calVec[i].subIntervals.size(); ++j) {
        for (auto I : calVec[i].subIntervals[j]) {
          double s = I.second.first;
          double e = I.second.second;
 
          double dist1 = (s - shortEnd >= 0)   ? (s - shortEnd) : 1e20;
          double dist2 = (shortStart - e >= 0) ? (shortStart - e) : 1e20;
          double dist = std::min(dist1, dist2);
 
          if (dist < distMin) {
            distMin = dist;
            iMin = i;
            jMin = j;
          }
        }
      }
    }
 
    B2ASSERT("Must be found", iMin != -1 && jMin != -1);
    calVec[iMin].subIntervals[jMin].insert(shortRun);
  }

breaks2intervalsSep()

std::vector< std::map< ExpRun, std::pair< double, double > > > breaks2intervalsSep	(	const std::map< ExpRun, std::pair< double, double > > &	runsMap,
		const std::vector< Atom > &	currVec,
		const std::vector< int > &	breaks )

Get calibration intervals according to the indexes of the breaks.

Parameters

runsMap	map defining the time ranges of the runs
currVec	vector with time intervals of the atoms (small non-divisible time intervals)
breaks	vector with integer indexes of the breaks

Returns

: a vector of the calib. intervals where each interval is a map with exp-run as a key and start- end-time as a value

Definition at line 284 of file Splitter.cc.

  {
    std::vector<std::map<ExpRun, std::pair<double, double>>> splitsNow(breaks.size() + 1);
    for (int i = 0; i < int(breaks.size()) + 1; ++i) {
      int s, e;
      std::tie(s, e) = getStartEndIndexes(currVec.size(), breaks, i);
 
      // loop over atoms in single calib interval
      for (int k = s; k <= e; ++k) {
        ExpRun r = getRun(runsMap, (currVec[k].t1 + currVec[k].t2) / 2.); //runexp of the atom
        if (splitsNow[i].count(r)) { //if already there
          splitsNow[i].at(r).first  = std::min(splitsNow[i].at(r).first, currVec[k].t1);
          splitsNow[i].at(r).second = std::max(splitsNow[i].at(r).second, currVec[k].t2);
        } else { //if new
          splitsNow[i][r].first  = currVec[k].t1;
          splitsNow[i][r].second = currVec[k].t2;
        }
      }
    }
 
    return splitsNow;
  }

convertSplitPoints()

template<typename Evt>

std::vector< ExpRunEvt > convertSplitPoints	(	const std::vector< Evt > &	events,
		const std::vector< double > &	splitPoints )

Convert splitPoints [hours] to breakPoints in ExpRunEvt.

Parameters

events	vector of events
splitPoints	the vector containing times of the edges of the calibration intervals [hours]

Returns

a vector with calibration break-points in the exp-run-evt format

Definition at line 363 of file Splitter.h.

  {
 
    std::vector<ExpRunEvt>  breakPos;
    for (auto p : splitPoints) {
      auto pos = getPosition(events, p);
      breakPos.push_back(pos);
    }
    return breakPos;
  }

decodeNumber()

unsigned decodeNumber ( double val )

inline

Decode the integer number encoded in val.

Definition at line 208 of file calibTools.h.

  {
    double factor = pow(FLT_RADIX, DBL_MANT_DIG);
    static const long long fEnc   = 4294967296; // pow(2, 32), 32 binary digits  for encoded number
 
    int e;
    double mantisa = std::frexp(val, &e);
    long long mantisaI = mantisa * factor;
 
    return (mantisaI % fEnc);
  }

dynamicBreaks()

std::vector< int > dynamicBreaks ( const std::vector< Atom > & runs )

private

Get optimal break points using algorithm based on dynamic programming.

Parameters

runs	Vector of atoms, where each atom is an intervals in time

Returns

: Optimal indexes of the break points

Definition at line 242 of file Splitter.cc.

  {
    //reset cache
    cache.resize(runs.size());
    for (auto& c : cache)
      c = std::make_pair(-1.0, std::vector<int>({}));
 
 
    std::vector<int> breaks;
    getMinLoss(runs, runs.size() - 1, breaks); //the minLoss (output) currently not used, only breaks
 
    return breaks;
  }

encodeNumber()

double encodeNumber ( double val, unsigned num )

inline

Encode integer num into double val such that val is nearly not changed (maximally by a relative shift 1e-6).

It is use to store time information to the payloads

Definition at line 186 of file calibTools.h.

  {
    double factor = pow(FLT_RADIX, DBL_MANT_DIG);
    static const long long fEnc   = 4294967296; // pow(2, 32), 32 binary digits  for encoded number
 
    int e; //exponent of the number
    double mantisa = std::frexp(val, &e);
    long long mantisaI = mantisa * factor; //mantissa as integer
 
    if (val != 0)
      mantisaI = (mantisaI / fEnc) * fEnc  + num; //adding encoded number to last digits of mantissa
    else {
      mantisaI = factor / 2 + num;
      e = -100; //if the val is zero, ensure very small number by the exponent
    }
 
    double newVal = ldexp(mantisaI / factor, e);
 
    return newVal;
  }

eval()

double eval ( const std::vector< double > & spl, const std::vector< double > & vals, double x )

inline

Evaluate spline (zero order or first order) in point x.

Definition at line 115 of file tools.h.

  {
    int order = -1;
    if (spl.size() == 0)
      order = 0;
    else if (spl.size() == vals.size() - 1)
      order = 0;
    else if (spl.size() == vals.size())
      order = 1;
    else {
      B2FATAL("Unknown order of spline");
    }
    B2ASSERT("Spline order should be zero or one", order == 0 || order == 1);
 
    if (order == 1) {
      B2ASSERT("Linear spline only meaningful for two or more nodes", spl.size() >= 2);
      B2ASSERT("As nodes as values in lin. spline", spl.size() == vals.size());
 
      if (x <= spl[0]) return vals[0];
      if (x >= spl.back()) return vals.back();
 
      // binary search for position
      int i1 = lower_bound(spl.begin(), spl.end(), x) - spl.begin() - 1;
 
      if (!(spl[i1] <= x && x <= spl[i1 + 1])) {
        B2FATAL("Wrong place founded : " <<  spl[i1] << " " << x << " " << spl[i1 + 1]);
      }
 
      // Linear interpolation between neighbouring nodes
      double v = ((spl[i1 + 1] - x) * vals[i1] + (x - spl[i1]) * vals[i1 + 1]) / (spl[i1 + 1] - spl[i1]);
      return v;
    } else if (order == 0) { //zero order polynomial
      B2ASSERT("#values vs #nodes in zero-order spline", spl.size() + 1 == vals.size());
      if (vals.size() == 1) {
        return vals[0];
      } else {
        double res = vals[0]; //by default value from lowest node
        for (unsigned i = 0; i < spl.size(); ++i) {
          if (spl[i] <= x) res = vals[i + 1];
          else break;
        }
        return res;
      }
    }
    return -99;
  }

evalPol() [1/2]

double evalPol	(	const Eigen::VectorXd &	polCoef,
		double	x )

Evaluate Cheb.

pol at point x when the coefficients of the expansion are provided

evalPol() [2/2]

double evalPol	(	const VectorXd &	polCoef,
		double	x )

Evaluate Cheb.

pol at point x when the coefficients of the expansion are provided

Definition at line 176 of file nodes.cc.

  {
    VectorXd pols = getPols(polCoef.size(), x);
 
    double s = pols.dot(polCoef);
 
    return s;
  }

extrapolateCalibration()

void extrapolateCalibration ( std::vector< CalibrationData > & calVec )

inline

Extrapolate calibration to intervals where it failed.

Definition at line 102 of file calibTools.h.

  {
    //put closest neighbor, where the statistic was low or algo failed
    for (unsigned i = 0; i < calVec.size(); ++i) {
      if (calVec[i].pars.cnt.size() != 0) continue;
      const auto& r =  calVec[i].subIntervals;
      double Start, End;
      std::tie(Start, End) = Splitter::getStartEnd(r);
 
      Eigen::Vector3d ipNow = Eigen::Vector3d::Zero();
      Eigen::MatrixXd ipeNow;
      Eigen::MatrixXd sizeMatNow;
 
      double distMin = 1e20;
      //Find the closest calibrated interval
      for (unsigned j = 0; j < calVec.size(); ++j) {
        if (calVec[j].isCalibrated == false) continue; //skip not-calibrated intervals
        const auto& rJ = calVec[j].subIntervals;
        for (unsigned jj = 0; jj < rJ.size(); ++jj) { //loop over subintervals
          const auto& rNow = rJ[jj];
          double s = rNow.begin()->second.first;
          double e = rNow.rbegin()->second.second;
 
          double dist1 = (s - End >= 0) ? (s - End) : 1e20;
          double dist2 = (Start - e >= 0) ? (Start - e) : 1e20;
          double dist = std::min(dist1, dist2);
 
          if (dist < distMin) {
            ipNow = calVec[j].pars.cnt.at(jj);
            ipeNow = calVec[j].pars.cntUnc.at(jj);
            sizeMatNow = calVec[j].pars.spreadMat;
            distMin = dist;
          }
        }
      }
 
      //Store it to vectors
      calVec[i].pars.cnt.resize(r.size());
      calVec[i].pars.cntUnc.resize(r.size());
      for (unsigned ii = 0; ii < r.size(); ++ii) {
        calVec[i].pars.cnt.at(ii) = ipNow;
        calVec[i].pars.cntUnc.at(ii) = ipeNow;
      }
      calVec[i].pars.spreadMat = sizeMatNow;
    }
 
  }

filter()

std::map< ExpRun, std::pair< double, double > > filter	(	const std::map< ExpRun, std::pair< double, double > > &	runs,
		double	cut,
		std::map< ExpRun, std::pair< double, double > > &	runsRemoved )

filter events to remove runs shorter than cut, it stores removed runs in runsRemoved

Definition at line 38 of file Splitter.cc.

  {
    std::map<ExpRun, std::pair<double, double>> runsCopy;
 
    for (auto r : runs) {
      const auto& I = r.second;
      double d =  I.second - I.first;
      if (d > cut)
        runsCopy[r.first] = r.second;
      else
        runsRemoved[r.first] = r.second;
    }
 
    return runsCopy;
  }

fitData()

std::pair< Pars, MatrixXd > fitData	(	const Pars &	pars,
		Limits	limits,
		bool	UseCheb = true )

Fit the data with specified initial values of parameters and limits on them.

Definition at line 130 of file ChebFitter.cc.

  {
    m_useCheb = UseCheb;
 
    ROOT::Math::Minimizer* minimum =
      ROOT::Math::Factory::CreateMinimizer("Minuit2", "");
 
    // set tolerance , etc...
    minimum->SetMaxFunctionCalls(10000000); // for Minuit/Minuit2
    minimum->SetMaxIterations(100000);  // for GSL
    minimum->SetTolerance(10.0);
    //minimum->SetPrecision(1e-5);
 
    //minimum->SetPrintLevel(3); //many outputs
    minimum->SetPrintLevel(0); //few outputs
    minimum->SetStrategy(2);
    minimum->SetErrorDef(1);
 
 
    // Set the free variables to be minimized !
    m_parNames.clear();
    int k = 0;
    for (auto p : pars) {
      std::string n    = p.first;
      double vCnt = p.second;
      if (limits.count(n) == 1) {
        double vMin = limits.at(n).first;
        double vMax = limits.at(n).second;
        double step = (vMax - vMin) / 100;
        minimum->SetLimitedVariable(k, n, vCnt, step, vMin, vMax);
      } else {
        double step = 1;
        minimum->SetVariable(k, n, vCnt, step);
      }
      m_parNames.push_back(n);
      ++k;
    }
 
    // create function wrapper for minimizer
    // a IMultiGenFunction type
    ROOT::Math::Functor f(*this, pars.size());
    minimum->SetFunction(f);
 
 
    // do the minimization
    minimum->Minimize();
 
 
    Pars parsF;
    for (unsigned i = 0; i < m_parNames.size(); ++i)
      parsF[m_parNames[i]] = minimum->X()[i];
 
    MatrixXd covMat(parsF.size(), parsF.size());
    for (unsigned i = 0; i < parsF.size(); ++i)
      for (unsigned j = 0; j < parsF.size(); ++j)
        covMat(i, j) = minimum->CovMatrix(i, j);
 
    // print pars
    std::stringstream log;
    log << "Minuit status : " << minimum->Status() << ", ";
    for (auto p : parsF)
      log << "\"" << p.first << "\" : " << p.second << ", ";
 
    B2INFO(log.str());
 
    delete minimum;
 
    return std::make_pair(parsF, covMat);
  }

getCoefs()

MatrixXd getCoefs	(	int	Size,
		bool	isInverse = false )

Transformation matrix between Cheb nodes and coefficients of the Cheb polynomials Notice, that there are two alternative ways defining polynomial interpolation:

• coefficients c_i in of the Cheb polynomials, i.e. f(x) = sum_i c_i P_i(x)
• Values of the f(x) in the Cheb nodes, i.e. d_j = f(x_j), where x_j are the nodes The Chebyshev polynomials are defined for x between 0 and 1

Definition at line 127 of file nodes.cc.

  {
    const int N = Size - 1;
    assert(N % 2 == 0);
 
    MatrixXd  coef(Size, Size);
 
    double mul = 1;
    double C = 1. / N;
    if (isInverse == true) {C = 1. / 2; }
 
    for (int k = 0; k <= N; ++k) {
      if (!isInverse) {
        coef(k, N) = C;
        coef(k, 0) = C * (k % 2 == 1 ? -1 : 1);
      } else {
        mul = k % 2 == 1 ? -1 : 1;
        coef(N - k, N) = C * mul;
        coef(N - k, 0) = C ;
      }
 
      for (int n = 1; n <= N - 1; ++n) {
        double el = cos(n * k * M_PI / N) * 2.*C * mul;
        if (!isInverse) coef(k, N - n) = el;
        else           coef(N - k, N - n) = el;
      }
    }
 
    return coef;
  }

getCoefsCheb()

MatrixXd getCoefsCheb

(

int

oldSize

)

Transformation matrix between Cheb nodes and Cheb coefficients with better normalization of the borders.

Definition at line 162 of file nodes.cc.

  {
    auto coef = getCoefs(oldSize);
 
    coef.row(0) *= 0.5;
    coef.row(coef.rows() - 1) *= 0.5;
 
    return coef;
  }

getDataGrid()

VectorXd getDataGrid ( ) const

private

Calculate Chebyshev coefficients for the data set.

Definition at line 109 of file ChebFitter.cc.

  {
    double a = m_nodes[0];
    double b = m_nodes[m_nodes.size() - 1];
 
 
    VectorXd polSum = VectorXd::Zero(m_nodes.size());
    for (double x : m_data) {
      double xx = (x - a) / (b - a); //normalize between 0 and 1
      polSum += getPols(m_nodes.size(), xx);
    }
 
 
    //transform to the basis of the cheb m_nodes
    VectorXd gridVals = m_coefsMat * polSum;
 
    return gridVals;
  }

getID()

int getID ( const std::vector< double > & breaks, double t )

inline

get id of the time point t

Definition at line 60 of file calibTools.h.

  {
    for (int i = 0; i < int(breaks.size()) + 1; ++i) {
      double s = (i == 0)             ? 0 : breaks[i - 1];
      double e = (i == int(breaks.size())) ? 1e20 : breaks[i];
      if (s  <= t   &&  t < e)
        return i;
    }
    return -1;
  }

getLogFunction()

VectorXd getLogFunction ( Pars pars ) const

private

Get the -2*log(p(x)) on the Cheb nodes.

Definition at line 43 of file ChebFitter.cc.

  {
    VectorXd fVals(m_nodes.size());
 
    //calc function values
    fVals = m_nodes.unaryExpr([&](double x) { return m_myFun(x, pars); });
 
    //normalize the function
    double I = fVals.dot(m_weights);
 
    fVals = -2 * log(fVals.array() / I); //normalize by integral
 
    return fVals;
 
  }

getLogLikelihoodFast()

double getLogLikelihoodFast ( const Pars & pars ) const

private

Calculate log likelihood using approximation based on Chebyshev polynomials (typically faster)

Definition at line 89 of file ChebFitter.cc.

  {
    VectorXd funVals = getLogFunction(pars);
    double LL = funVals.dot(m_dataGrid);
 
    return LL;
  }

getLogLikelihoodSlow()

double getLogLikelihoodSlow ( const Pars & pars ) const

private

Calculate log likelihood using exact formula.

Definition at line 76 of file ChebFitter.cc.

  {
 
    double L = 0;
    for (double d : m_data) {
      double v = m_myFun(d, pars);
      L += -2 * log(v);
    }
 
    return L;
  }

getMinima()

std::pair< double, double > getMinima ( const std::vector< std::vector< double > > & vals, int i0, int j0 )

inline

Get minimum inside and outside of the smaller window defined by i0, j0.

Definition at line 23 of file minimizer.h.

  {
    int N = vals.size();
    int Nzoom = (N + 1) / 2;
 
 
    double minIn  = std::numeric_limits<double>::max();
    double minOut = std::numeric_limits<double>::max();
 
    for (int i = 0; i < N; ++i)
      for (int j = 0; j < N; ++j) {
        if ((i0 <= i && i < i0 + Nzoom) &&
            (j0 <= j && j < j0 + Nzoom))
          minIn = std::min(minIn, vals[i][j]);
        else
          minOut = std::min(minOut, vals[i][j]);
      }
    return {minIn, minOut};
  }

getMinimum()

std::vector< double > getMinimum ( std::function< double(double, double)> fun, double xMin, double xMax, double yMin, double yMax )

inline

Get minimum of 2D function in the rectangular domain defined by xMin,xMax & yMin,yMax.

Definition at line 44 of file minimizer.h.

  {
    const int N = 17; //the grid has size N x N
    const int Nzoom = (N + 1) / 2; // the size of the zoomed rectangle where minimum is
 
 
    const int kMax = 35; //number of iterations
 
    std::vector<std::vector<double>> vals(N);
    for (auto& v : vals) v.resize(N);
 
    for (int k = 0; k < kMax; ++k) {
 
      // get values of the function
      for (int i = 0; i < N; ++i)
        for (int j = 0; j < N; ++j) {
          double x = xMin + i * (xMax - xMin) / (N - 1.);
          double y = yMin + j * (yMax - yMin) / (N - 1.);
          vals[i][j] = fun(x, y);
        }
 
      if (k == kMax - 1) break;
 
      double mOutMax = - std::numeric_limits<double>::max();
      int iOpt = -1, jOpt = -1;
      //find optimal rectangle
      for (int i = 0; i < N - Nzoom; ++i)
        for (int j = 0; j < N - Nzoom; ++j) {
          double mIn, mOut;
          std::tie(mIn, mOut) = getMinima(vals, i, j);
 
          if (mOut > mOutMax) {
            mOutMax = mOut;
            iOpt = i;
            jOpt = j;
          }
        }
 
      //Zoom to the optimal rectangle
      // get values of the function
 
      double xMinNow = xMin + iOpt * (xMax - xMin) / (N - 1.);
      double xMaxNow = xMin + (iOpt + Nzoom - 1) * (xMax - xMin) / (N - 1.);
 
      double yMinNow = yMin + jOpt * (yMax - yMin) / (N - 1.);
      double yMaxNow = yMin + (jOpt + Nzoom - 1) * (yMax - yMin) / (N - 1.);
 
      xMin = xMinNow;
      xMax = xMaxNow;
      yMin = yMinNow;
      yMax = yMaxNow;
 
    }
 
 
    //get the overall minimum
    double minTot = std::numeric_limits<double>::max();
    int iOpt = -1, jOpt = -1;
 
    for (int i = 0; i < N; ++i)
      for (int j = 0; j < N; ++j) {
        if (vals[i][j] < minTot) {
          minTot = vals[i][j];
          iOpt = i;
          jOpt = j;
        }
      }
 
    double xMinNow = xMin + iOpt * (xMax - xMin) / (N - 1.);
    double yMinNow = yMin + jOpt * (yMax - yMin) / (N - 1.);
 
    return {xMinNow, yMinNow};
  }

getMinLoss()

double getMinLoss ( const std::vector< Atom > & vec, int e, std::vector< int > & breaks )

private

Recursive function to evaluate minimal sum of the lossFuctions for the optimal clustering.

It return the minimum of the lossFunction and optimal break points giving such value. It acts only on atoms with indexes between 0 and e

Parameters

	vec	Vector of atoms, where each atom is an intervals in time
	e	the index of the last atom included in the optimisation problem
[out]	breaks	Output vector with indexes of optimal break points

Returns

: Minimal value of the summed loss function

Definition at line 205 of file Splitter.cc.

  {
    // If entry in cache (speed up)
    if (cache[e].first >= 0) {
      breaks = cache[e].second;
      return cache[e].first;
    }
 
 
    std::vector<int> breaksOpt;
    double minVal = 1e30;
    int iMin = -10;
    for (int i = -1; i <= e - 1; ++i) {
      auto breaksNow = breaks;
      double r1 = 0;
      if (i != -1)
        r1 = getMinLoss(vec, i, breaksNow);
      double r2 = lossFunction(vec, i + 1, e);
      double tot = r1 + r2;
 
      if (tot < minVal) { //store minimum
        minVal = tot;
        iMin = i;
        breaksOpt = breaksNow;
      }
    }
 
    if (iMin != -1)
      breaksOpt.push_back(iMin);
 
 
    breaks = breaksOpt;
    cache[e] = std::make_pair(minVal, breaks); //store solution to cache
    return minVal;
  }

getNodes()

VectorXd getNodes

(

int

Size

)

Get the vector of positions of the Chebyshev nodes The nodes are by definition between 0 and 1, there are Size nodes For the corresponding weights use getWeights.

Definition at line 65 of file nodes.cc.

  {
    assert((Size - 1) % 2 == 0);
    VectorXd xi = VectorXd::Zero(Size);
    for (int i = 0; i < Size; ++i) {
      double Cos = cos(i / (Size - 1.) * M_PI);
      xi[i] = (1 - Cos) / 2;
    }
    return xi;
  }

getPols()

VectorXd getPols	(	int	Size,
		double	x )

Evaluate Chebyshev polynomials up to Size at point x It returns a vector of the P_i(x) for i=0..Size-1 The polynomial is defined for x between 0 and 1.

Definition at line 77 of file nodes.cc.

  {
    VectorXd pol(Size);
    double C = 2 * (2 * x - 1);
 
    if (Size >= 1) pol[0] = 1;
    if (Size >= 2) pol[1] = C / 2;
 
    for (int i = 2; i < Size; ++i)
      pol[i] = C * pol[i - 1] - pol[i - 2];
    return pol;
  }

getPolsSum()

VectorXd getPolsSum	(	int	Size,
		VectorXd	x )

Calculate the Chebyshev polynomials of order i=0..Size-1 at points given in vector x_j and sum it over point index j It returns sum_j P_i(x_j) for i=0..Size-1 The Chebyshev polynomials are defined for x between 0 and 1.

Definition at line 96 of file nodes.cc.

  {
    assert(Size > 2);
 
    VectorXd polSum(Size);
 
    VectorXd pol0 = 0 * x.array() + 1;
    VectorXd pol1 = 2 * x.array() - 1;
    VectorXd C    = 2 * pol1;
 
    VectorXd pol2(x.size());
    for (int i = 2; i < Size; ++i) {
      polSum(i - 2) = pol0.sum();
 
      pol2 = C.array() * pol1.array() - pol0.array();
 
      pol0 = pol1;
      pol1 = pol2;
    }
 
    polSum(Size - 2) = pol0.sum();
    polSum(Size - 1) = pol1.sum();
 
    return polSum;
  }

getPosition()

template<typename Evt>

ExpRunEvt getPosition ( const std::vector< Evt > & events, double tEdge )

inline

Get the exp-run-evt number from the event time [hours].

Parameters

events	vector of events
tEdge	the event time of the event of interest [hours]

Returns

the position of the time point in the exp-run-evt format

Definition at line 341 of file Splitter.h.

  {
    ExpRunEvt evt(-1, -1, -1);
    double tBreak = -1e10;
    for (const auto& e : events) {
      if (e.t < tEdge) {
        if (e.t > tBreak) {
          tBreak = e.t;
          evt =  ExpRunEvt(e.exp, e.run, e.evtNo);
        }
      }
    }
    return evt;
  }

getRangeLin()

std::vector< double > getRangeLin ( int nVals, double xMin, double xMax )

inline

Equidistant range between xMin and xMax for spline of the first order.

Definition at line 83 of file tools.h.

  {
    B2ASSERT("At least one value in the spline required", nVals >= 1);
    if (nVals == 1) return {};
    std::vector<double> v(nVals);
    for (int i = 0; i < nVals; ++i)
      v[i] = xMin + i * (xMax - xMin) / (nVals - 1);
    return v;
  }

getRangeZero()

std::vector< double > getRangeZero ( int nVals, double xMin, double xMax )

inline

Equidistant range between xMin and xMax for spline of the zero order.

Definition at line 94 of file tools.h.

  {
    B2ASSERT("At least one value in the spline required", nVals >= 1);
    if (nVals == 1) return {};
    std::vector<double> v(nVals - 1);
    for (int i = 1; i < nVals; ++i)
      v[i - 1] = xMin + i * (xMax - xMin) / (nVals);
    return v;
  }

getRun()

static ExpRun getRun ( const std::map< ExpRun, std::pair< double, double > > & runs, double t )

static

Get exp number + run number from time.

Parameters

runs	map, where key contain the exp-run number and value the start- and end-time of the run
t	time of interest [hours]

Returns

: the exp-run number at the input time t

Definition at line 263 of file Splitter.cc.

  {
    ExpRun rFound(-1, -1);
    int nFound = 0;
    for (auto r : runs) { //Linear search over runs
      if (r.second.first <= t && t < r.second.second) {
        ++nFound;
        if (nFound == 2
            && rFound.run == r.first.run - 1) // to fix corner case when a time is seen in two consecutive runs due to float rounding errors
          --nFound;
        rFound = r.first;
      }
    }
 
    B2ASSERT("Exactly one interval should be found", nFound == 1);
    B2ASSERT("Assert that something was found", rFound != ExpRun(-1, -1));
    return rFound;
  }

getRunInfo()

template<typename Evt>

std::map< ExpRun, std::pair< double, double > > getRunInfo ( const std::vector< Evt > & evts )

inline

Get the map of runs, where each run contains pair with start/end time [hours].

Parameters

evts	vector of events

Returns

a map where the key is exp-run and value start/end time of the particular run [hours]

Definition at line 312 of file Splitter.h.

  {
    std::map<ExpRun, std::pair<double, double>> runsInfo;
 
    for (const auto& evt : evts) {
      int Exp = evt.exp;
      int Run = evt.run;
      double time = evt.t;
      if (runsInfo.count(ExpRun(Exp, Run))) {
        double tMin, tMax;
        std::tie(tMin, tMax) = runsInfo.at(ExpRun(Exp, Run));
        tMin = std::min(tMin, time);
        tMax = std::max(tMax, time);
        runsInfo.at(ExpRun(Exp, Run)) = {tMin, tMax};
      } else {
        runsInfo[ExpRun(Exp, Run)] = {time, time};
      }
 
    }
    return runsInfo;
  }

getStartEndIndexes()

std::pair< int, int > getStartEndIndexes	(	int	nIntervals,
		const std::vector< int > &	breaks,
		int	indx )

get the range of interval with nIntervals and breaks stored in a vector

Definition at line 83 of file Splitter.cc.

  {
    B2ASSERT("There must be at least one interval", nIntervals >= 1);
    B2ASSERT("Interval index must be positive", indx >= 0);
    B2ASSERT("Interval index must be smaller than #breaks", indx < int(breaks.size()) + 1); //There is one more interval than #breaks
    int s = (indx == 0) ? 0 : breaks[indx - 1] + 1;
    int e = (indx == int(breaks.size())) ? nIntervals - 1 : breaks[indx];
    return {s, e};
  }

getSVDDedxFromParticle()

VXDDedxTrack const * getSVDDedxFromParticle

(

Particle const *

particle

)

SVD dEdx value from particle.

Definition at line 36 of file DedxVariables.cc.

  {
    const Track* track = particle->getTrack();
    if (!track) {
      return nullptr;
    }
 
    const VXDDedxTrack* dedxTrack = track->getRelatedTo<VXDDedxTrack>();
    if (!dedxTrack) {
      return nullptr;
    }
    return dedxTrack;
  }

getWeights()

VectorXd getWeights

(

int

Size

)

Get the vector of weights to calculate the integral over the Chebyshev nodes The nodes are by definition between 0 and 1, there are Size nodes To get their positions, use getNodes.

Definition at line 29 of file nodes.cc.

  {
    const int N = Size - 1;
    assert(N % 2 == 0);
 
    std::vector<std::vector<double>> coef(Size);
    for (auto& el : coef) el.resize(Size);
 
 
    for (int k = 0; k <= N / 2; ++k) {
      coef[2 * k][N] = 1. / N;
      coef[2 * k][0] = 1. / N ;
 
      coef[2 * k][N / 2] = 2. / N * (2 * ((k + 1) % 2) - 1);
 
      for (int n = 1; n <= N / 2 - 1; ++n)
        coef[2 * k][n] = coef[2 * k][N - n] = 2. / N * cos(n * k * M_PI * 2 / N);
    }
 
    VectorXd wgt = VectorXd::Zero(Size);
 
 
    for (int i = 0; i < Size; ++i) {
      wgt[i] += coef[0][i];
      wgt[i] += coef[N][i] / (1. - N * N);
      for (int k = 1; k <= N / 2 - 1; ++k) {
        double w = 2. / (1 - 4 * k * k);
        wgt[i] += w * coef[2 * k][i];
      }
 
      wgt[i] *= 0.5; //for interval (0,1)
    }
    return wgt;
  }

init()

void init	(	int	Size,
		double	xMin,
		double	xMax )

Initialize the fitter (the Chebyshev coefficients)

Definition at line 59 of file ChebFitter.cc.

  {
    // loading the Cheb nodes
    m_nodes = (xMax - xMin) * getNodes(Size).array() + xMin;
 
    // loading the weights for integration
    m_weights = (xMax - xMin) * getWeights(Size);
 
 
    // calculate the transformation matrix from pol coefs to grid points
    m_coefsMat = getCoefsCheb(Size).transpose();
 
    m_dataGrid = getDataGrid();
 
  }

interpol() [1/3]

VectorXd interpol	(	const VectorXd &	xi,
		double	x )

Get Interpolation vector k_i for point x from the function values at points xi (polynomial interpolation) In the second step, the function value at x can be evaluated as sum_i vals_i k_i.

Definition at line 189 of file nodes.cc.

  {
    double Norm = (xi[xi.size() - 1] - xi[0]) / 2;
    VectorXd coefs(xi.size());
    for (int i = 0; i < xi.size(); ++i) {
      double num = 1, den = 1;
      for (int j = 0; j < xi.size(); ++j)
        if (j != i) {
          num *= (x     - xi(j)) / Norm;
          den *= (xi(i) - xi(j)) / Norm;
        }
      coefs(i) = num / den;
    }
    return coefs;
  }

interpol() [2/3]

double interpol	(	Eigen::VectorXd	xi,
		Eigen::VectorXd	vals,
		double	x )

Get interpolated function value at point x when function values vals at points xi are provided.

If the points xi are fixed and only vals are different between interpol calls, use interpol(xi, x) to speed up the evaluation.

interpol() [3/3]

double interpol	(	VectorXd	xi,
		VectorXd	vals,
		double	x )

Get interpolated function value at point x when function values vals at points xi are provided.

If the points xi are fixed and only vals are different between interpol calls, use interpol(xi, x) to speed up the evaluation.

Definition at line 210 of file nodes.cc.

  {
    VectorXd coefs = interpol(xi, x);
    return coefs.dot(vals);
  }

lossFunction()

double lossFunction ( const std::vector< Atom > & vec, int s, int e ) const

private

lossFunction of the calibration interval consisting of several "atoms" stored in vector vec The atoms included in calibration interval have indices between s and e

the lossFunction formula (it can be modified according to the user's taste)

Parameters

vec	Vector of atoms, where each atom is an intervals in time
s	First index of the calib. interval
e	Last index of the calib. interval

Returns

: A value of the loss function

Definition at line 135 of file Splitter.cc.

  {
 
    //raw time
    double rawTime = vec[e].t2 - vec[s].t1;
 
    //max gap
    double maxGap = 0;
    for (int i = s; i <= e - 1; ++i) {
      double d = vec[i + 1].t1 - vec[i].t2;
      maxGap = std::max(maxGap, d);
    }
 
    //net time
    double netTime = 0;
    for (int i = s; i <= e; ++i) {
      netTime += vec[i].t2 - vec[i].t1;
    }
 
    // Number of events
    double nEv = 0;
    for (int i = s; i <= e; ++i) {
      nEv += vec[i].nEv;
    }
    if (nEv == 0) nEv = 0.1;
 
    //double loss = (rawTime - tBest) * (rawTime - tBest) + gapPenalty * (maxGap * maxGap);
    //double loss = 1./nEv  + timePenalty * (rawTime * rawTime);
 
    lossFun->SetParameters(rawTime, netTime, maxGap, nEv);
    double lossNew = lossFun->Eval(0);
 
    return lossNew;
  }

merge()

std::vector< std::vector< double > > merge ( const std::vector< std::vector< std::vector< double > > > & toMerge )

inline

merge { vector<double> a, vector<double> b} into {a, b}

Definition at line 41 of file tools.h.

  {
    std::vector<std::vector<double>> allVecs;
    for (const auto& v : toMerge)
      allVecs.insert(allVecs.end(), v.begin(), v.end());
    return allVecs;
  }

mergeIntervals()

std::map< ExpRun, std::pair< double, double > > mergeIntervals ( const std::map< ExpRun, std::pair< double, double > > & I1, const std::map< ExpRun, std::pair< double, double > > & I2 )

static

Merge two subintervals into one subinterval.

Parameters

I1	First subinterval to merge
I2	Second subinterval to merge

Returns

: The resulting subinterval

Definition at line 310 of file Splitter.cc.

  {
    std::map<ExpRun, std::pair<double, double>>  I = I1;
    for (auto r : I2) {
      ExpRun run = r.first;
      if (I.count(run) == 0)
        I[run] = r.second;
      else {
        I.at(run) = std::make_pair(std::min(I1.at(run).first, I2.at(run).first),   std::max(I1.at(run).second, I2.at(run).second));
      }
    }
    return I;
  }

operator!=()

bool operator!= ( ExpRun a, ExpRun b )

inline

Not equal for ExpRun.

Definition at line 71 of file Splitter.h.

{ return (a.exp != b.exp || a.run != b.run); }

operator()()

double operator()

(

const double *

par

)

const

Evaluate the log likelihood.

Definition at line 97 of file ChebFitter.cc.

  {
    Pars pars;
    for (unsigned i = 0; i < m_parNames.size(); ++i)
      pars[m_parNames[i]] = par[i];
 
    double LL = m_useCheb ? getLogLikelihoodFast(pars) : getLogLikelihoodSlow(pars);
 
    return LL;
  }

operator<()

bool operator< ( ExpRun a, ExpRun b )

inline

less than for ExpRun

Definition at line 74 of file Splitter.h.

{return ((a.exp < b.exp) || (a.exp == b.exp && a.run < b.run));}

plotRuns()

void plotRuns

(

const std::vector< std::pair< double, double > > &

runs

)

plot runs on time axis

Definition at line 57 of file Splitter.cc.

  {
    TGraphErrors* gr = new TGraphErrors();
 
 
    for (auto r : runs) {
      double m = (r.first + r.second) / 2;
      double e = (r.second - r.first) / 2;
 
      gr->SetPoint(gr->GetN(), m, 1);
      gr->SetPointError(gr->GetN() - 1, e, 0);
    }
 
    gStyle->SetEndErrorSize(6);
 
    gr->SetLineWidth(1);
    gr->SetMarkerSize(40);
    gr->Draw("ape");
    gr->GetYaxis()->SetRangeUser(-10, 10);
    gr->GetXaxis()->SetRangeUser(0, 256);
 
    gr->GetXaxis()->SetTitle("time [hours]");
 
  }

plotSRuns()

void plotSRuns	(	const std::vector< std::pair< double, double > > &	runs,
		const std::vector< int > &	breaks,
		int	offset = 2 )

plot clusters or runs on time axis

Definition at line 94 of file Splitter.cc.

  {
    TGraphErrors* gr = new TGraphErrors();
 
    for (int i = 0; i < int(breaks.size()) + 1; ++i) {
      int s, e;
      std::tie(s, e) = getStartEndIndexes(runs.size(), breaks, i);
      double a = runs[s].first;
      double b = runs[e].second;
 
      double m = (a + b) / 2;
      double err = (b - a) / 2;
 
      gr->SetPoint(gr->GetN(), m, offset);
      gr->SetPointError(gr->GetN() - 1, err, 0);
 
    }
 
    gr->SetLineColor(kRed);
    gr->SetMarkerColor(kRed);
    gr->Draw("pe same");
  }

printBySize()

void printBySize

(

const std::vector< std::pair< double, double > > &

runs

)

print sorted lengths of the runs

Definition at line 119 of file Splitter.cc.

  {
    std::vector<double> dist;
    for (auto r : runs) {
      double d = r.second - r.first;
      dist.push_back(d);
    }
 
    sort(dist.begin(), dist.end());
 
    for (auto d : dist)
      B2INFO(d);
 
  }

rn()

TString rn ( )

inline

Get random string.

Definition at line 38 of file tools.h.

{return Form("%d", gRandom->Integer(1000000000)); }

runAlgorithm()

template<typename Evt, typename Fun>

CalibrationData runAlgorithm ( const std::vector< Evt > & evts, std::vector< std::map< ExpRun, std::pair< double, double > > > range, Fun runCalibAnalysis )

inline

run calibration algorithm for single calibration interval

Definition at line 336 of file calibTools.h.

  {
    CalibrationData calD;
    auto& r = range;
    double rStart, rEnd;
    std::tie(rStart, rEnd) = Splitter::getStartEnd(r);
    B2INFO("Start of loop startTime endTime : " << rStart << " " << rEnd);
 
    auto breaks =  Splitter::getBreaks(r);
 
    std::vector<Evt> evtsNow;
 
    std::vector<int> Counts(breaks.size() + 1, 0);
    // Select events belonging to the interval
    for (const auto& ev : evts) {
      if (rStart <= ev.t && ev.t < rEnd) {
        evtsNow.push_back(ev);
        ++Counts.at(getID(breaks, ev.t));
      }
    }
 
    B2ASSERT("Number of intervals vs number of breakPoints", r.size()  == breaks.size() + 1);
 
    //Merge smallest interval if with low stat (try it 10times)
    for (int k = 0; k < 10; ++k)  {
      int iMin = min_element(Counts.begin(), Counts.end()) - Counts.begin();
      if (Counts.size() >= 2 &&  Counts[iMin] < 50) { //merge with neighbor if possible
        auto iM = -1;
        if (iMin == 0)
          iM = iMin + 1;
        else if (iMin == int(Counts.size()) - 1)
          iM = iMin - 1;
        else {
          if (Counts[iMin + 1] < Counts[iMin - 1])
            iM = iMin + 1;
          else
            iM = iMin - 1;
        }
        B2ASSERT("Number of intervals equal to size of counters", r.size() == Counts.size());
 
        r.at(iM) = Splitter::mergeIntervals(r[iM], r[iMin]);
        r.erase(r.begin() + iMin);
        breaks =  Splitter::getBreaks(r);
        Counts[iM] += Counts[iMin];
        Counts.erase(Counts.begin() + iMin);
      }
    }
 
    B2INFO("#events "  << " : " << evtsNow.size());
    B2INFO("Breaks size " << " : " << breaks.size());
 
    calD.breakPoints = convertSplitPoints(evtsNow, breaks);
 
    calD.subIntervals = r;
 
    if (breaks.size() > 0)
      B2INFO("StartOfCalibInterval (run,evtNo,vtxIntervalsSize) " <<   calD.breakPoints.at(0).run << " " <<
             calD.breakPoints.at(0).evt << " " << calD.breakPoints.size());
 
 
    //If too few events, let have the output empty
    //Will be filled with the closest neighbor at the next stage
    if (evtsNow.size() < 50) {
      return calD;
    }
 
    // Run the calibration
    B2INFO("Start of running calibration over calibration interval");
    tie(calD.pars.cnt, calD.pars.cntUnc, calD.pars.spreadMat) = runCalibAnalysis(evtsNow, breaks);
    calD.pars.pulls.resize(calD.pars.cnt.size());
    B2INFO("End of running analysis - SpreadMatX : " << sqrt(abs(calD.pars.spreadMat(0, 0))));
    B2ASSERT("All subintervals have calibration of the mean value", calD.pars.cnt.size() == r.size());
    B2ASSERT("All subintervals have calibration of the unc. of mean", calD.pars.cntUnc.size() == r.size());
 
    calD.isCalibrated = true;
 
    return calD;
  }

runCalibration()

template<typename Fun1, typename Fun2>

CalibrationAlgorithm::EResult runCalibration	(	TTree *	tracks,
		const std::string &	calibName,
		Fun1	GetEvents,
		Fun2	calibAnalysis,
		std::function< TObject *(Eigen::VectorXd, Eigen::MatrixXd, Eigen::MatrixXd)>	calibObjCreator,
		TString	m_lossFunctionOuter,
		TString	m_lossFunctionInner )

Run the the calibration over the whole event sample.

Parameters

tracks	TTree object with mu-mu events
calibName	name of the calibration payload
GetEvents	function that transforms TTree to std::vector
calibAnalysis	function that performs the calibration on a single calibration interval
calibObjCreator	function that stores results to the payload class which inherits from TObject
m_lossFunctionOuter	Lost function for the calibration intervals of the spread parameters
m_lossFunctionInner	Lost function for the calibration subintervals (for the mean value parameters)

Returns

State of the calibration run, i.e. EResult::c_OK if everything OK

Definition at line 429 of file calibTools.h.

  {
    // Check that there are at least some data
    if (!tracks || tracks->GetEntries() < 15) {
      if (tracks)
        B2WARNING("Too few data : " << tracks->GetEntries());
      return CalibrationAlgorithm::EResult::c_NotEnoughData;
    }
    B2INFO("Number of tracks: " << tracks->GetEntries());
 
    // Tree to vector of Events
    auto evts = GetEvents(tracks);
 
    //Time range for each ExpRun
    std::map<ExpRun, std::pair<double, double>> runsInfoOrg = getRunInfo(evts);
    std::map<ExpRun, std::pair<double, double>> runsRemoved; //map with time intervals of very short runs
    auto runsInfo = filter(runsInfoOrg, 2. / 60, runsRemoved); //include only runs longer than 2mins
 
    // If nothing remains
    if (runsInfo.size() == 0) {
      B2WARNING("Too short run");
      return  CalibrationAlgorithm::EResult::c_NotEnoughData;
    }
 
    // Get intervals based on the input loss functions
    Splitter splt;
    auto splits = splt.getIntervals(runsInfo, evts, m_lossFunctionOuter, m_lossFunctionInner);
 
    //Loop over all calibration intervals
    std::vector<CalibrationData> calVec;
    for (const auto& s : splits) {
      CalibrationData calD = runAlgorithm(evts, s, calibAnalysis); // run the calibration over the interval s
      calVec.push_back(calD);
    }
 
    // extrapolate results to the low-stat intervals
    extrapolateCalibration(calVec);
 
    // Include removed short runs
    for (auto shortRun : runsRemoved) {
      addShortRun(calVec,  shortRun);
    }
 
    // Store Payloads to files
    storePayloads(evts, calVec, calibName, calibObjCreator);
 
    return CalibrationAlgorithm::EResult::c_OK;
  }

slice() [1/2]

std::vector< double > slice ( const std::vector< double > & v, unsigned ind, unsigned n )

inline

put slice of original vector v[ind:ind+n] into new one, n is number of elements

Definition at line 106 of file tools.h.

  {
    std::vector<double> vNew;
    for (unsigned i = ind; i < ind + n && i < v.size(); ++i)
      vNew.push_back(v[i]);
    return vNew;
  }

slice() [2/2]

std::vector< Atom > slice ( std::vector< Atom > vec, int s, int e )

inline

Slice the vector to contain only elements with indexes s .. e (included)

Definition at line 85 of file Splitter.h.

  {
    return std::vector<Atom>(vec.begin() + s, vec.begin() + e + 1);
  }

splitToSmall()

std::vector< std::pair< double, double > > splitToSmall ( const std::map< ExpRun, std::pair< double, double > > & runs, double intSize = 1. / 60 )

staticprivate

Split the runs into small calibration intervals (atoms) of a specified size.

By definition each of these intervals spans only over single run. These will be clustered into larger intervals in the next steps

Parameters

runs	Runs to split into the atoms
intSize	Intended size of the small intervals

Returns

: A vector with resulting time boundaries of the atoms

Definition at line 173 of file Splitter.cc.

  {
    // split into small intervals
    std::vector<std::pair<double, double>> smallRuns;
 
    for (auto r : runs) {
      const auto& I = r.second;
      if (intSize < 0) {
        smallRuns.push_back(I);
        continue;
      }
 
      double runTime = I.second - I.first;
      int nSplits = runTime / intSize; //1-m intervals
      nSplits = std::max(1, nSplits); //at least 1 interval
 
      for (int i = 0; i < nSplits; ++i) {
        double L = I.first + i * (runTime / nSplits);
        double H = I.first + (i + 1) * (runTime / nSplits);
        smallRuns.push_back({L, H});
      }
    }
    return smallRuns;
  }

storePayloads()

template<typename Evt>

void storePayloads ( const std::vector< Evt > & evts, const std::vector< CalibrationData > & calVecConst, const std::string & objName, std::function< TObject *(Eigen::VectorXd, Eigen::MatrixXd, Eigen::MatrixXd) > getCalibObj )

inline

Store payloads to files.

Definition at line 225 of file calibTools.h.

  {
    auto calVec = calVecConst;
 
    // Loop to store payloads
    ExpRun exprunLast(-1, -1); //last exprun
    EventDependency* intraRun = nullptr;
 
    // Loop over calibration intervals
    for (unsigned i = 0; i < calVec.size(); ++i) {
      const auto& r = calVec[i].subIntervals; //   splits[i];
      // Loop over calibration subintervals
      for (int k = 0; k < int(r.size()); ++k) {
 
        for (auto I : r[k]) { //interval required to be within single run
          ExpRun exprun = I.first;
 
          //Encode Start+End time in seconds of the payload
          if (calVec[i].pars.cntUnc.at(k).rows() == 3) {
            calVec[i].pars.cntUnc.at(k)(0, 1) = calVec[i].pars.cntUnc.at(k)(1, 0) = encodeNumber(calVec[i].pars.cntUnc.at(k)(0, 1),
                                                round(I.second.first  * 3600));
            calVec[i].pars.cntUnc.at(k)(0, 2) = calVec[i].pars.cntUnc.at(k)(2, 0) = encodeNumber(calVec[i].pars.cntUnc.at(k)(0, 2),
                                                round(I.second.second * 3600));
          } else {
            calVec[i].pars.cntUnc.at(k)(0, 0) = encodeNumber(calVec[i].pars.cntUnc.at(k)(0, 0), round(I.second.first  * 3600));
            calVec[i].pars.spreadMat(0, 0)    = encodeNumber(calVec[i].pars.spreadMat(0, 0),    round(I.second.second * 3600));
          }
 
          TObject* obj = getCalibObj(calVec[i].pars.cnt.at(k), calVec[i].pars.cntUnc.at(k), calVec[i].pars.spreadMat);
          if (exprun != exprunLast) { //if new run
            if (intraRun) { //if not first -> store
              auto m_iov = IntervalOfValidity(exprunLast.exp, exprunLast.run, exprunLast.exp, exprunLast.run);
              Database::Instance().storeData(objName, intraRun, m_iov);
            }
 
            intraRun = new EventDependency(obj);
          } else {
            int breakPoint;
            if (k - 1 >= 0) {
              breakPoint = calVec[i].breakPoints.at(k - 1).evt;
              B2ASSERT("Payload saving consistency", calVec[i].breakPoints.at(k - 1).run == exprun.run);
            } else {
              B2ASSERT("Payload saving consistency", i != 0);
              double rStart, rEnd;
              std::tie(rStart, rEnd) = Splitter::getStartEnd(r);
              auto pos = getPosition(evts, rStart);
              breakPoint = pos.evt;
              B2ASSERT("Payload saving consistency", pos.run == exprun.run);
            }
            intraRun->add(breakPoint, obj);
          }
          exprunLast = exprun;
        }
      } //end loop over calibration subintervals
 
    } //end loop over calibration intervals
 
    //Store the last entry
    auto m_iov = IntervalOfValidity(exprunLast.exp, exprunLast.run, exprunLast.exp, exprunLast.run);
    Database::Instance().storeData(objName, intraRun, m_iov);
  }

storePayloadsNoIntraRun()

void storePayloadsNoIntraRun ( const std::vector< CalibrationData > & calVecConst, const std::string & objName, std::function< TObject *(Eigen::VectorXd, Eigen::MatrixXd, Eigen::MatrixXd) > getCalibObj )

inline

Store payloads to files, where calib data have no intra-run dependence.

Definition at line 291 of file calibTools.h.

  {
    auto calVec = calVecConst;
 
    // Check that there is no intra-run dependence
    std::set<ExpRun> existingRuns;
    for (unsigned i = 0; i < calVec.size(); ++i) {
      const auto& r = calVec[i].subIntervals;
      // Loop over calibration subintervals
      for (int k = 0; k < int(r.size()); ++k) {
 
        for (auto I : r[k]) {
          ExpRun exprun = I.first;
          // make sure that the run isn't already in the list, to avoid duplicity
          if (existingRuns.count(exprun) != 0)
            B2FATAL("Intra-run dependence exists");
          existingRuns.insert(exprun);
        }
      }
    }
 
 
    // Loop over calibration intervals
    for (unsigned i = 0; i < calVec.size(); ++i) {
      const auto& r = calVec[i].subIntervals; //   splits[i];
      // Loop over calibration subintervals
      for (unsigned k = 0; k < r.size(); ++k) {
 
        TObject* obj = getCalibObj(calVec[i].pars.cnt.at(k), calVec[i].pars.cntUnc.at(k), calVec[i].pars.spreadMat);
 
        ExpRun start = (r[k].cbegin()->first);
        ExpRun last  = (r[k].crbegin()->first);
 
        auto iov = IntervalOfValidity(start.exp, start.run, last.exp, last.run);
        Database::Instance().storeData(objName, obj, iov);
 
 
      } //end loop over calibration subintervals
    } //end loop over calibration intervals
 
  }

toB2Vector3()

B2Vector3D toB2Vector3 ( Eigen::VectorXd vIn )

inline

Function that converts Eigen vector to ROOT vector.

Definition at line 54 of file calibTools.h.

  {
    return B2Vector3D(vIn(0), vIn(1), vIn(2));
  }

toTMatrixDSym()

TMatrixDSym toTMatrixDSym ( Eigen::MatrixXd mIn )

inline

Function that converts Eigen symmetric matrix to ROOT matrix.

Definition at line 44 of file calibTools.h.

  {
    TMatrixDSym mOut(mIn.rows());
    for (int i = 0; i < mIn.rows(); ++i)
      for (int j = 0; j < mIn.cols(); ++j)
        mOut(i, j) = (mIn(i, j) + mIn(j, i)) / 2.;
    return mOut;
  }

vec2vec() [1/2]

std::vector< double > vec2vec ( const Eigen::VectorXd & v )

inline

ROOT vector -> std vector.

Definition at line 61 of file tools.h.

  {
    std::vector<double> vNew(v.rows());
    for (int i = 0; i < v.rows(); ++i)
      vNew[i] = v(i);
    return vNew;
  }

vec2vec() [2/2]

Eigen::VectorXd vec2vec ( const std::vector< double > & vec )

inline

std vector -> ROOT vector

Definition at line 51 of file tools.h.

  {
    Eigen::VectorXd v(vec.size());
    for (unsigned i = 0; i < vec.size(); ++i) {
      v[i] = vec[i];
    }
    return v;
  }

vecs2mat()

Eigen::MatrixXd vecs2mat ( const std::vector< std::vector< double > > & vecs )

inline

merge columns (from std::vectors) into ROOT matrix

Definition at line 72 of file tools.h.

  {
    Eigen::MatrixXd m(vecs[0].size(), vecs.size());
    for (unsigned i = 0; i < vecs[0].size(); ++i)
      for (unsigned j = 0; j < vecs.size(); ++j) {
        m(i, j) = vecs[j][i];
      }
    return m;
  }