JRECONSTRUCTION::JMuonPrefit Struct Reference

Wrapper class to make pre-fit of muon trajectory. More...

#include <JMuonPrefit.hh>

Inheritance diagram for JRECONSTRUCTION::JMuonPrefit:

Classes
struct	cmz
	Auxiliary data structure for sorting of hits. More...
struct	input_type
	Input data type. More...

Public Types
typedef JEstimator< JLine1Z >	JEstimator_t
typedef JTRIGGER::JHitR1	hit_type
typedef std::vector< hit_type >	buffer_type

Public Member Functions
	JMuonPrefit (const JMuonPrefitParameters_t &parameters, const int debug=0)
	Constructor.
	JMuonPrefit (const JMuonPrefitParameters_t &parameters, const JOmega3D &omega, const int debug=0)
	Constructor.
input_type	getInput (const JModuleRouter &router, const JDAQEvent &event, const coverage_type &coverage) const
	Get input data.
JEvt	operator() (const input_type &input)
	Fit function.
void	reset ()
	Reset fit parameters.
bool	equals (const JMuonPrefitParameters_t &parameters) const
	Equality.
	ClassDef (JMuonPrefitParameters_t, 1)

Public Attributes
JOmega3D	omega
int	debug
int	factoryLimit
	factory limit for combinatorics
int	NMaxHits
	maximal number of hits
double	sigma_ns
	time resolution [ns]
double	gridAngle_deg
	grid angle for directions [deg]
int	useL0
	option for L0 hit use
int	numberOfOutliers
	maximum number of outliers
size_t	numberOfPrefits
	number of prefits
size_t	numberOfPostfits
	number of postfits
double	TMaxLocal_ns
	time window for local coincidences [ns]
double	ctMin
	minimal cosine space angle between PMT axes
double	roadWidth_m
	road width [m]
double	Qwatershed
	minimal quality step

Static Public Attributes
static const struct JRECONSTRUCTION::JMuonPrefit::cmz	cmz
static constexpr double	STANDARD_DEVIATIONS = 3.0
	number of standard deviations for outlier removal [unit]

Private Attributes
buffer_type	data
JMatrixNZ	V
JVectorNZ	Y

Detailed Description

Wrapper class to make pre-fit of muon trajectory.

The JMuonPrefit fit is used to generate start values for subsequent fits (usually JMuonSimplex and JMuonGandalf).
To this end, a scan of directions is made and the time and transverse positions of the track are fitted for each direction (JFIT::JEstimator<JLine1Z>).
The directions are spaced by the parameters JMuonPrefitParameters_t::gridAngle_deg.
This angle corresponds to the envisaged angular accuracy of the result.
The probability that one of the results is less than this angle away from the correct value, multiple start values should be considered (JMuonPrefitParameters_t::numberOfPrefits).
Note that the CPU time scales with the inverse of the square of this angle.
The chi-squared is based on the time residuals.

Definition at line 79 of file JMuonPrefit.hh.

Member Typedef Documentation

JEstimator_t

JEstimator<JLine1Z> JRECONSTRUCTION::JMuonPrefit::JEstimator_t

Definition at line 83 of file JMuonPrefit.hh.

hit_type

JTRIGGER::JHitR1 JRECONSTRUCTION::JMuonPrefit::hit_type

Definition at line 84 of file JMuonPrefit.hh.

buffer_type

std::vector<hit_type> JRECONSTRUCTION::JMuonPrefit::buffer_type

Definition at line 85 of file JMuonPrefit.hh.

Constructor & Destructor Documentation

JMuonPrefit() [1/2]

JRECONSTRUCTION::JMuonPrefit::JMuonPrefit ( const JMuonPrefitParameters_t & parameters, const int debug = 0 )

inline

Constructor.

Parameters

parameters	parameters
debug	debug

Definition at line 126 of file JMuonPrefit.hh.

                                                          :
      JMuonPrefitParameters_t(parameters),
      omega (parameters.gridAngle_deg * JMATH::PI/180.0),
      debug (debug)
    {}

JMuonPrefit() [2/2]

JRECONSTRUCTION::JMuonPrefit::JMuonPrefit ( const JMuonPrefitParameters_t & parameters, const JOmega3D & omega, const int debug = 0 )

inline

Constructor.

Parameters

parameters	parameters
omega	directions
debug	debug

Definition at line 141 of file JMuonPrefit.hh.

                                                          :
      JMuonPrefitParameters_t(parameters),
      omega (omega),
      debug (debug)
    {}

Member Function Documentation

getInput()

input_type JRECONSTRUCTION::JMuonPrefit::getInput ( const JModuleRouter & router, const JDAQEvent & event, const coverage_type & coverage ) const

inline

Get input data.

Parameters

router	module router
event	event
coverage	coverage

Returns

input data

Definition at line 158 of file JMuonPrefit.hh.

    {
      using namespace std;
      using namespace JTRIGGER;
      using namespace KM3NETDAQ;
 
      const JBuildL0<hit_type> buildL0;
      const JBuildL2<hit_type> buildL2(JL2Parameters(2, TMaxLocal_ns, ctMin));
 
      input_type   input(event.getDAQEventHeader(), coverage);
    
      buffer_type& dataL0 = input.dataL0;
      buffer_type& dataL1 = input.dataL1;
 
      buildL2(event, router, !useL0, back_inserter(dataL1));
 
      // 3D cluster of unique optical modules
 
      const JMatch3B<hit_type> match3B(roadWidth_m, TMaxLocal_ns);
 
      sort(dataL1.begin(), dataL1.end(), hit_type::compare);
 
      buffer_type::iterator __end = dataL1.end();
 
      __end = unique(dataL1.begin(), __end, equal_to<JDAQModuleIdentifier>());
    
      __end = clusterizeWeight(dataL1.begin(), __end, match3B);
 
      dataL1.erase(__end, dataL1.end());
 
 
      if (useL0) {
 
        buildL0(event, router, true, back_inserter(dataL0));
 
        __end = dataL0.end();
 
        for (buffer_type::iterator i = dataL0.begin(); i != __end; ) {
 
          if (match3B.count(*i, dataL1.begin(), dataL1.end()) != 0)
            ++i;
          else
            swap(*i, *--__end);
        }
 
        dataL0.erase(__end, dataL0.end());
      }
 
      return input;
    }

operator()()

JEvt JRECONSTRUCTION::JMuonPrefit::operator() ( const input_type & input )

inline

Fit function.

Parameters

input

input data

Returns

fit results

Definition at line 219 of file JMuonPrefit.hh.

    {
      using namespace std;
      using namespace JPP;
 
      const double HIT_OFF  =  1.0e3 * sigma_ns * sigma_ns;  // [ns^2]
 
      const JMatch1D<hit_type> match1D(roadWidth_m, TMaxLocal_ns);
 
      const buffer_type& dataL0 = input.dataL0;
      const buffer_type& dataL1 = input.dataL1;
 
      data.reserve(dataL0.size() + 
                   dataL1.size());
 
      JEvent event(JMUONPREFIT);
 
      JEvt out;
 
      for (JOmega3D_t::const_iterator dir = omega.begin(); dir != omega.end(); ++dir) {
 
        const JRotation3D R(*dir);
 
 
        buffer_type::iterator __end = copy(dataL1.begin(), dataL1.end(), data.begin());
 
        for (buffer_type::iterator i = data.begin(); i != __end; ++i) {
          i->rotate(R);
        }
 
 
        // reduce data
 
        if (distance(data.begin(), __end) > NMaxHits) {
 
          buffer_type::iterator __p = data.begin();
 
          advance(__p, NMaxHits);
 
          partial_sort(data.begin(), __p, __end, cmz);
 
          __end = __p;
        }
 
 
        // 1D cluster
 
        __end = clusterizeWeight(data.begin(), __end, match1D);
 
        if (useL0 && 
            distance(data.begin(), __end) < useL0   &&
            distance(data.begin(), __end) < NMaxHits) {
 
          buffer_type::iterator p = __end;                                     // begin L0 data
          buffer_type::iterator q = copy(dataL0.begin(), dataL0.end(), p);     // end   L0 data
 
          for (buffer_type::iterator i = p; i != q; ++i) {
 
            if (find_if(data.begin(), __end, make_predicate(&hit_type::getModuleID, i->getModuleID())) == __end) {
 
              i->rotate(R);
 
              if (match1D.count(*i, data.begin(), __end) != 0) {
                *p = *i;
                ++p;
              }
            }
          }
 
          __end = clusterize(__end, p, match1D);
        }
 
          
        if (distance(data.begin(), __end) <= NUMBER_OF_PARAMETERS) {
          continue;
        }
 
 
        // 1D fit
 
        JLine1Z  tz;
        double   chi2 = numeric_limits<double>::max();
        int      NDF  = distance(data.begin(), __end) - NUMBER_OF_PARAMETERS;
        int      N    = getCount(data.begin(), __end);
 
        try {
 
          sort(data.begin(), __end, hit_type::compare);
 
          (*this)(data.begin(), __end);
 
          V.set(*this, data.begin(), __end, gridAngle_deg, sigma_ns);
          Y.set(*this, data.begin(), __end);
 
          V.invert();
 
          tz   = *this;
          chi2 = getChi2(Y, V);
        }
        catch(const exception& error) {
          continue;
        }
 
 
        // outlier removal
 
        if (distance(data.begin(), __end) <= factoryLimit) {
 
          int number_of_outliers = numberOfOutliers;
 
          if (number_of_outliers > NDF - 1) {
            number_of_outliers = NDF - 1;
          }
 
          double ymin = chi2 - STANDARD_DEVIATIONS * STANDARD_DEVIATIONS;
            
          for (int n = 1; n <= number_of_outliers; ++n, ymin -= STANDARD_DEVIATIONS * STANDARD_DEVIATIONS) {
              
            sort(data.begin(), __end, hit_type::compare);
              
            buffer_type::iterator __end1 = prev(__end, n);
 
            do {
 
              try {
                
                (*this)(data.begin(), __end1);
                
                V.set(*this, data.begin(), __end1, gridAngle_deg, sigma_ns);
                Y.set(*this, data.begin(), __end1);
                
                V.invert();
                
                double y = getChi2(Y, V);
                
                if (y <= -(STANDARD_DEVIATIONS * STANDARD_DEVIATIONS)) {
 
                  WARNING(endl << "JMuonPrefit: invalid chi2 " << y << endl);
 
                } else {
 
                  if (y < 0.0) {
                    y = 0.0;
                  }
 
                  if (y < ymin) {
                    ymin = y;
                    tz   = *this;
                    chi2 = ymin;
                    NDF  = distance(data.begin(), __end1) - NUMBER_OF_PARAMETERS;
                    N    = getCount(data.begin(), __end1);
                  }
                }
              }
              catch(const exception& error) {}
 
            } while (next_permutation(data.begin(), __end1, __end, hit_type::compare));
          }
 
        } else {
 
          const int number_of_outliers = NDF - 1;
 
          for (int n = 1; n <= number_of_outliers; ++n) {
              
            double ymax =  0.0;
            int    k    = -1;
                
            for (size_t i = 0; i != Y.size(); ++i) {
                  
              double y = getChi2(Y, V, i);
                  
              if (y > ymax) {
                ymax = y;
                k    = i;
              }
            }
                
            if (ymax < STANDARD_DEVIATIONS * STANDARD_DEVIATIONS) {
              break;
            }
 
            try {
 
              V.update(k, +HIT_OFF);
 
              this->update(data.begin(), __end, V);
 
              Y.set(*this, data.begin(), __end);
 
              tz   = *this;
              chi2 = getChi2(Y, V);
              NDF -= 1;
              N   -= getCount(data[k]);
            }
            catch(const exception& error) {
              break;
            }
          }
        }
 
        tz.rotate_back(R);
 
        out.push_back(getFit(event(), tz, *dir, getQuality(chi2, N, NDF), NDF));
 
        // set additional values
 
        out.rbegin()->setW(JPP_COVERAGE_ORIENTATION, input.coverage.orientation);
        out.rbegin()->setW(JPP_COVERAGE_POSITION,    input.coverage.position);
      }
 
 
      if (numberOfPrefits > 0) {
 
        JEvt::iterator __end = out.end();
 
        if (Qwatershed > 0) {                         // sort distinct maxima
 
          __end = out.begin();                        // output
 
          for (JEvt::iterator p1 = out.begin(); p1 != out.end() && (size_t) distance(out.begin(), __end) < numberOfPrefits; ) {
 
            JEvt::iterator p2 = p1;                   // p2 -> best-of-the-rest
 
            for (JEvt::iterator i = p1; i != out.end(); ++i) {
              if (qualitySorter(*i, *p2)) {
                p2 = i;
              }
            }
 
            swap(*p2, *p1);                           // put best-of-the-rest at head
 
            p2 = p1++;                                // p2 = head; p1 = next to head
 
            const JPointing pointing(*p2);            // reference direction
 
            sort(p1, out.end(), pointing);            // sort according space angle with respect to reference direction
 
            for (double Q = p2->getQ();
                 p1 != out.end() && (p1->getQ() >= p2->getQ() - Qwatershed ||   // bottom
                                     p1->getQ() <= Q);                          // descend
                 ++p1) {
              Q = p1->getQ();
            }
 
            swap(*(__end++), *p2);                    // output head
          }
 
        } else {
 
          advance(__end = out.begin(), min(numberOfPrefits, out.size()));
 
          partial_sort(out.begin(), __end, out.end(), qualitySorter);
        }
 
        if (numberOfPostfits > 0) {
 
          __end = gridify(__end, out.end(), numberOfPostfits);
        }
 
        out.erase(__end, out.end());
 
      } else {
 
        sort(out.begin(), out.end(), qualitySorter);
      }
 
      return out;
    }

reset()

void JRECONSTRUCTION::JMuonPrefitParameters_t::reset ( )

inlineinherited

Reset fit parameters.

Definition at line 39 of file JMuonPrefitParameters_t.hh.

    {
      factoryLimit     = 8;
      NMaxHits         = 50;
      sigma_ns         = 5;
      gridAngle_deg    = 1;
      useL0            = 0;
      numberOfOutliers = 3;
      numberOfPrefits  = 12;
      numberOfPostfits = 0;
      TMaxLocal_ns     = 18.0;
      ctMin            = 0;
      roadWidth_m      = 200;
      Qwatershed       = 0.0;
    }

equals()

bool JRECONSTRUCTION::JMuonPrefitParameters_t::equals ( const JMuonPrefitParameters_t & parameters ) const

inlineinherited

Equality.

Parameters

parameters

fit parameters

Returns

true if equals; else false

Definition at line 61 of file JMuonPrefitParameters_t.hh.

    {
      return (this->factoryLimit     == parameters.factoryLimit     &&
              this->NMaxHits         == parameters.NMaxHits         &&
              this->sigma_ns         == parameters.sigma_ns         &&
              this->gridAngle_deg    == parameters.gridAngle_deg    &&
              this->useL0            == parameters.useL0            &&
              this->numberOfOutliers == parameters.numberOfOutliers &&         
              this->numberOfPrefits  == parameters.numberOfPrefits  &&
              this->numberOfPostfits == parameters.numberOfPostfits &&
              this->TMaxLocal_ns     == parameters.TMaxLocal_ns     &&
              this->roadWidth_m      == parameters.roadWidth_m      &&
              this->Qwatershed       == parameters.Qwatershed);
    }

ClassDef()

JRECONSTRUCTION::JMuonPrefitParameters_t::ClassDef ( JMuonPrefitParameters_t , 1  )

inherited

Member Data Documentation

cmz

const struct JRECONSTRUCTION::JMuonPrefit::cmz JRECONSTRUCTION::JMuonPrefit::cmz

static

omega

JOmega3D JRECONSTRUCTION::JMuonPrefit::omega

Definition at line 512 of file JMuonPrefit.hh.

debug

int JRECONSTRUCTION::JMuonPrefit::debug

Definition at line 513 of file JMuonPrefit.hh.

STANDARD_DEVIATIONS

double JRECONSTRUCTION::JMuonPrefit::STANDARD_DEVIATIONS = 3.0

staticconstexpr

number of standard deviations for outlier removal [unit]

Definition at line 515 of file JMuonPrefit.hh.

data

buffer_type JRECONSTRUCTION::JMuonPrefit::data

private

Definition at line 518 of file JMuonPrefit.hh.

V

JMatrixNZ JRECONSTRUCTION::JMuonPrefit::V

private

Definition at line 519 of file JMuonPrefit.hh.

Y

JVectorNZ JRECONSTRUCTION::JMuonPrefit::Y

private

Definition at line 520 of file JMuonPrefit.hh.

factoryLimit

int JRECONSTRUCTION::JMuonPrefitParameters_t::factoryLimit

inherited

factory limit for combinatorics

Definition at line 78 of file JMuonPrefitParameters_t.hh.

NMaxHits

int JRECONSTRUCTION::JMuonPrefitParameters_t::NMaxHits

inherited

maximal number of hits

Definition at line 79 of file JMuonPrefitParameters_t.hh.

sigma_ns

double JRECONSTRUCTION::JMuonPrefitParameters_t::sigma_ns

inherited

time resolution [ns]

Definition at line 80 of file JMuonPrefitParameters_t.hh.

gridAngle_deg

double JRECONSTRUCTION::JMuonPrefitParameters_t::gridAngle_deg

inherited

grid angle for directions [deg]

Definition at line 81 of file JMuonPrefitParameters_t.hh.

useL0

int JRECONSTRUCTION::JMuonPrefitParameters_t::useL0

inherited

option for L0 hit use

Definition at line 82 of file JMuonPrefitParameters_t.hh.

numberOfOutliers

int JRECONSTRUCTION::JMuonPrefitParameters_t::numberOfOutliers

inherited

maximum number of outliers

Definition at line 83 of file JMuonPrefitParameters_t.hh.

numberOfPrefits

size_t JRECONSTRUCTION::JMuonPrefitParameters_t::numberOfPrefits

inherited

number of prefits

Definition at line 84 of file JMuonPrefitParameters_t.hh.

numberOfPostfits

size_t JRECONSTRUCTION::JMuonPrefitParameters_t::numberOfPostfits

inherited

number of postfits

Definition at line 85 of file JMuonPrefitParameters_t.hh.

TMaxLocal_ns

double JRECONSTRUCTION::JMuonPrefitParameters_t::TMaxLocal_ns

inherited

time window for local coincidences [ns]

Definition at line 86 of file JMuonPrefitParameters_t.hh.

ctMin

double JRECONSTRUCTION::JMuonPrefitParameters_t::ctMin

inherited

minimal cosine space angle between PMT axes

Definition at line 87 of file JMuonPrefitParameters_t.hh.

roadWidth_m

double JRECONSTRUCTION::JMuonPrefitParameters_t::roadWidth_m

inherited

road width [m]

Definition at line 88 of file JMuonPrefitParameters_t.hh.

Qwatershed

double JRECONSTRUCTION::JMuonPrefitParameters_t::Qwatershed

inherited

minimal quality step

Definition at line 89 of file JMuonPrefitParameters_t.hh.

---

The documentation for this struct was generated from the following file:

• JMuonPrefit.hh