JRECONSTRUCTION::JMuonEnergy Class Reference

Auxiliary class to to determine muon energy. More...

#include <JMuonEnergy.hh>

Inheritance diagram for JRECONSTRUCTION::JMuonEnergy:

Classes
struct	input_type
	Input data type. More...
struct	JResult
	Auxiliary class for energy estimation. More...

Public Types
typedef JRegressor< JEnergy >	JRegressor_t
typedef JModuleL0	module_type
typedef std::vector< module_type >	detector_type

Public Member Functions
	JMuonEnergy (const JMuonEnergyParameters_t &parameters, const storage_type &storage, const JPMTParametersMap &pmtParameters, const JEnergyCorrection &correct, const int debug=0)
	Constructor.
input_type	getInput (const JModuleRouter &router, const JSummaryRouter &summary, const JDAQEvent &event, const JEvt &in, const coverage_type &coverage) const
	Get input data.
JEvt	operator() (const input_type &input)
	Fit function.
void	reset ()
	Reset fit parameters.
bool	equals (const JMuonEnergyParameters_t &parameters) const
	Equality.
	ClassDef (JMuonEnergyParameters_t, 2)

Public Attributes
const JPMTParametersMap &	pmtParameters
JEnergyCorrection	correct
double	roadWidth_m
	road width [m]
double	R_Hz
	default rate [Hz]
size_t	numberOfPrefits
	number of prefits
double	EMin_log
	minimal energy [log10(GeV)]
double	EMax_log
	maximal energy [log10(GeV)]
double	TMin_ns
	minimal time w.r.t. Cherenkov hypothesis [ns]
double	TMax_ns
	maximal time w.r.t. Cherenkov hypothesis [ns]
double	ZMin_m
	minimal z-position [m]
double	resolution
	energy resolution [log10(GeV)]
int	mestimator
	M-estimator (see JFIT::JMEstimator_t)

Detailed Description

Auxiliary class to to determine muon energy.

Definition at line 69 of file JMuonEnergy.hh.

Member Typedef Documentation

JRegressor_t

JRegressor<JEnergy> JRECONSTRUCTION::JMuonEnergy::JRegressor_t

Definition at line 74 of file JMuonEnergy.hh.

module_type

JModuleL0 JRECONSTRUCTION::JMuonEnergy::module_type

Definition at line 75 of file JMuonEnergy.hh.

detector_type

std::vector<module_type> JRECONSTRUCTION::JMuonEnergy::detector_type

Definition at line 76 of file JMuonEnergy.hh.

Constructor & Destructor Documentation

JMuonEnergy()

JRECONSTRUCTION::JMuonEnergy::JMuonEnergy ( const JMuonEnergyParameters_t & parameters, const storage_type & storage, const JPMTParametersMap & pmtParameters, const JEnergyCorrection & correct, const int debug = 0 )

inline

Constructor.

Parameters

parameters	parameters
storage	storage
pmtParameters	PMT parameters
correct	energy correction
debug	debug

Definition at line 123 of file JMuonEnergy.hh.

                                                         :
      JMuonEnergyParameters_t(parameters),
      JRegressor_t (storage),
      pmtParameters(pmtParameters),
      correct      (correct)
    {
      using namespace JPP;
      
      if (this->getRmax() < roadWidth_m) {
        roadWidth_m = this->getRmax();
      }
 
      JRegressor_t::debug  = debug;
 
      this->estimator.reset(getMEstimator(mestimator));
    }

Member Function Documentation

getInput()

input_type JRECONSTRUCTION::JMuonEnergy::getInput ( const JModuleRouter & router, const JSummaryRouter & summary, const JDAQEvent & event, const JEvt & in, const coverage_type & coverage ) const

inline

Get input data.

Parameters

router	module router
summary	summary data
event	event
in	start values
coverage	coverage

Returns

input data

Definition at line 155 of file JMuonEnergy.hh.

    {
      using namespace std;
      using namespace JTRIGGER;
 
      input_type   input(event.getDAQEventHeader(), in, coverage);
 
      const JBuildL0 <JHitR0>   buildL0;
      map<int, vector<JHitR0> > data;
 
      const JDAQTimeslice timeslice(event, true);
 
      JSuperFrame2D<JHit> buffer;
 
      for (JDAQTimeslice::const_iterator i = timeslice.begin(); i != timeslice.end(); ++i) {
 
        if (router.hasModule(i->getModuleID())) {
 
          buffer(*i, router.getModule(i->getModuleID()));
 
          buildL0(buffer, back_inserter(data[i->getModuleID()]));
        }
      }
 
      for (const auto& module : router.getReference()) {
        if (!module.empty()) {
          input.data.push_back(module_type(module, summary.getSummaryFrame(module.getID(), R_Hz), data[module.getID()]));
        }
      }
 
      return input;
    }

operator()()

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

inline

Fit function.

Parameters

input

input data

Returns

fit results

Definition at line 231 of file JMuonEnergy.hh.

    {
      using namespace std;
      using namespace JPP;
 
      JEvent event(JMUONENERGY);
 
      JEvt out;
 
      // select start values
 
      JEvt in  = input.in;
 
      in.select(numberOfPrefits, qualitySorter);
 
      if (!in.empty()) {
        in.select(JHistory::is_event(in.begin()->getHistory()));
      }
 
      for (JEvt::const_iterator track = in.begin(); track != in.end(); ++track) {
 
        JFit fit(*track);
 
        // move track
 
        if (fit.getW(JSTART_LENGTH_METRES, 0.0) > 0.0) {
 
          fit.move(fit.getW(JSTART_ZMIN_M), getSpeedOfLight());
 
          fit.setW(JSTART_ZMIN_M, 0.0);
          fit.setW(JSTART_ZMAX_M, fit.getW(JSTART_LENGTH_METRES));
        }
 
        const JRotation3D  R (getDirection(fit));
        const JLine1Z      tz(getPosition (fit).rotate(R), fit.getT());
 
        double zmin = numeric_limits<double>::lowest();
 
        if (fit.getW(JSTART_LENGTH_METRES, 0.0) > 0.0) {
          zmin = this->ZMin_m;
        }
 
        vector<JNPEHit> data;
 
        for (const auto& module : input.data) {
 
          JPosition3D pos(module->getPosition());
 
          pos.transform(R, tz.getPosition());
 
          if (pos.getX() <= roadWidth_m) {
 
            const double z1 = pos.getZ()  -  pos.getX() / getTanThetaC();
            const double t1 = tz .getT()  + (pos.getZ() + pos.getX() * getKappaC()) * getInverseSpeedOfLight();
 
            if (z1 >= zmin) {
 
              for (size_t i = 0; i != module->size(); ++i) {
 
                if (module.getStatus(i)) {
 
                  const struct {
 
                    bool operator()(const JHitR0& hit) const
                    {
                      return (hit.getPMT() == pmt && T_ns(hit.getT()));
                    }
                    
                    const JTimeRange T_ns;
                    const size_t     pmt;
 
                  } match = { T_ns + t1, i };
 
                  const JPMTIdentifier id(module->getID(), i);
 
                  const double ps = getSurvivalProbability(pmtParameters.getPMTParameters(id));
 
                  JPMT pmt = module->getPMT(i);
 
                  pmt.transform(R, tz.getPosition());
 
                  const JNPEHit hit(this->getNPE(pmt, module.frame.getRate(i)), count_if(module.begin(), module.end(), match), ps);
 
                  DEBUG("hit: " << setw(8) << module->getID() << '.' << FILL(2,'0') << i << ' ' 
                        << FIXED(7,3) << hypot(pmt.getX(), pmt.getY())                   << ' '
                        << FIXED(7,3) << module.frame.getRate(i) * 1.0e-3                << ' '
                        << SCIENTIFIC(9,3) << hit.getY0() << ' '
                        << SCIENTIFIC(9,3) << hit.getY1() << ' '
                        << SCIENTIFIC(9,3) << hit.getYA() << ' '
                        << SCIENTIFIC(9,3) << hit.getYB() << ' '
                        << setw(2)         << hit.getN()  << endl);
 
                  data.push_back(hit);
                }
              }
            }
          }
        }
 
        const int NDF = distance(data.begin(), data.end()) - 1;
 
        if (NDF >= 0) {
          
          // 5-point search between given limits
 
          const int      N  =  5;
          
          JResult result[N];
          
          for (int i = 0; i != N; ++i) {
            result[i].x = EMin_log + i * (EMax_log - EMin_log) / (N-1);
          }
          
          map<double, double> buffer;
          
          do {
            
            int j = 0;
            
            for (int i = 0; i != N; ++i) {
              
              if (!result[i]) {
                
                const JEnergy x    = result[i].x;
                const double  chi2 = (*this)(x, data.begin(), data.end());
                
                result[i].chi2 = chi2;
                buffer[chi2]   = x.getE();
              }
              
              if (result[i].chi2 < result[j].chi2) {
                j = i;
              }
            }  
            
            
            for (int i = 0; i != N; ++i) {
              DEBUG(' ' << FIXED(5,2) << result[i].x << ' ' << FIXED(9,3) << result[i].chi2);
            }
            DEBUG(endl);
            
            // squeeze range
            
            switch (j) {
              
            case 0:
              result[4] = result[1];
              result[2] = result[0];
              result[0] = JResult(2*result[2].x - result[4].x);
              break;
              
            case 1:
              result[4] = result[2];
              result[2] = result[1];
              break;
              
            case 2:
              result[0] = result[1];
              result[4] = result[3];
              break;
              
            case 3:
              result[0] = result[2];
              result[2] = result[3];
              break;
              
            case 4:
              result[0] = result[3];
              result[2] = result[4];
              result[4] = JResult(2*result[2].x - result[0].x);
              break;
            }
            
            result[1] = JResult(0.5 * (result[0].x + result[2].x));
            result[3] = JResult(0.5 * (result[2].x + result[4].x));
            
          } while (result[4].x - result[0].x > resolution);
          
          
          if (result[1].chi2 != result[3].chi2) {
            
            result[2].x    += 0.25 * (result[3].x - result[1].x) * (result[1].chi2 - result[3].chi2) / (result[1].chi2 + result[3].chi2 - 2*result[2].chi2);
            result[2].chi2  = (*this)(result[2].x, data.begin(), data.end());
            
          }
          
          const double chi2 = result[2].chi2;
          const double E    = result[2].x.getE();
          
          // calculate additional variables
 
          double Emin = numeric_limits<double>::max();
          double Emax = numeric_limits<double>::lowest();
          
          for (map<double, double>::const_iterator i = buffer.begin(); i != buffer.end() && i->first <= chi2 + 1.0; ++i) {
            if (i->second < Emin) { Emin = i->second; }
            if (i->second > Emax) { Emax = i->second; }
          }
          
          const double mu_range   = gWater(E);           // range of a muon with the reconstructed energy
          
          double noise_likelihood = 0.0;                 // log-likelihood of every hit being from K40   
          int    number_of_hits   = 0;                   // number of hits selected for JEnergy
          
          for (vector<JNPEHit>::const_iterator i = data.begin(); i != data.end(); ++i) {
            noise_likelihood += log10(getP(i->getY0(), i->getN()));       // probability of getting the observed multiplicity with K40 
            number_of_hits   += i->getN();                                // hit multiplicity
          }
 
          fit.push_back(event());
 
          // set corrected energy
 
          fit.setE(pow(10.0, correct(log10(E))));
 
          out.push_back(fit);
 
          // set additional values
 
          out.rbegin()->setW(JENERGY_ENERGY,            E);
          out.rbegin()->setW(JENERGY_CHI2,              chi2);
          out.rbegin()->setW(JENERGY_MUON_RANGE_METRES, mu_range);
          out.rbegin()->setW(JENERGY_NOISE_LIKELIHOOD,  noise_likelihood);
          out.rbegin()->setW(JENERGY_NDF,               NDF);
          out.rbegin()->setW(JENERGY_NUMBER_OF_HITS,    number_of_hits);
          out.rbegin()->setW(JENERGY_MINIMAL_ENERGY,    Emin);
          out.rbegin()->setW(JENERGY_MAXIMAL_ENERGY,    Emax);
          out.rbegin()->setW(JPP_COVERAGE_ORIENTATION,  input.coverage.orientation);
          out.rbegin()->setW(JPP_COVERAGE_POSITION,     input.coverage.position);
        }
      }
 
      // apply default sorter
 
      sort(out.begin(), out.end(), qualitySorter);
 
      copy(input.in.begin(), input.in.end(), back_inserter(out));
 
      return out;
    }

reset()

void JRECONSTRUCTION::JMuonEnergyParameters_t::reset ( )

inlineinherited

Reset fit parameters.

Definition at line 44 of file JMuonEnergyParameters_t.hh.

    {
      roadWidth_m     = std::numeric_limits<double>::max();
      R_Hz            = 6.0e3;
      numberOfPrefits = 0;
      EMin_log        = 1.0;
      EMax_log        = 8.0;
      TMin_ns         =  -50.0;     
      TMax_ns         = +450.0;     
      ZMin_m          = std::numeric_limits<double>::lowest();
      resolution      = 0.01;
      mestimator      = JFIT::EM_NORMAL;
    }

equals()

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

inlineinherited

Equality.

Parameters

parameters

fit parameters

Returns

true if equals; else false

Definition at line 64 of file JMuonEnergyParameters_t.hh.

    {
      return (this->roadWidth_m      == parameters.roadWidth_m     &&
              this->R_Hz             == parameters.R_Hz            &&
              this->numberOfPrefits  == parameters.numberOfPrefits &&
              this->EMin_log         == parameters.EMin_log        &&
              this->EMax_log         == parameters.EMax_log        &&
              this->TMin_ns          == parameters.TMin_ns         &&
              this->TMax_ns          == parameters.TMax_ns         &&
              this->ZMin_m           == parameters.ZMin_m          &&
              this->resolution       == parameters.resolution      &&
              this->mestimator       == parameters.mestimator);
    }

ClassDef()

JRECONSTRUCTION::JMuonEnergyParameters_t::ClassDef ( JMuonEnergyParameters_t , 2  )

inherited

Member Data Documentation

pmtParameters

const JPMTParametersMap& JRECONSTRUCTION::JMuonEnergy::pmtParameters

Definition at line 472 of file JMuonEnergy.hh.

correct

JEnergyCorrection JRECONSTRUCTION::JMuonEnergy::correct

Definition at line 474 of file JMuonEnergy.hh.

roadWidth_m

double JRECONSTRUCTION::JMuonEnergyParameters_t::roadWidth_m

inherited

road width [m]

Definition at line 80 of file JMuonEnergyParameters_t.hh.

R_Hz

double JRECONSTRUCTION::JMuonEnergyParameters_t::R_Hz

inherited

default rate [Hz]

Definition at line 81 of file JMuonEnergyParameters_t.hh.

numberOfPrefits

size_t JRECONSTRUCTION::JMuonEnergyParameters_t::numberOfPrefits

inherited

number of prefits

Definition at line 82 of file JMuonEnergyParameters_t.hh.

EMin_log

double JRECONSTRUCTION::JMuonEnergyParameters_t::EMin_log

inherited

minimal energy [log10(GeV)]

Definition at line 83 of file JMuonEnergyParameters_t.hh.

EMax_log

double JRECONSTRUCTION::JMuonEnergyParameters_t::EMax_log

inherited

maximal energy [log10(GeV)]

Definition at line 84 of file JMuonEnergyParameters_t.hh.

TMin_ns

double JRECONSTRUCTION::JMuonEnergyParameters_t::TMin_ns

inherited

minimal time w.r.t. Cherenkov hypothesis [ns]

Definition at line 85 of file JMuonEnergyParameters_t.hh.

TMax_ns

double JRECONSTRUCTION::JMuonEnergyParameters_t::TMax_ns

inherited

maximal time w.r.t. Cherenkov hypothesis [ns]

Definition at line 86 of file JMuonEnergyParameters_t.hh.

ZMin_m

double JRECONSTRUCTION::JMuonEnergyParameters_t::ZMin_m

inherited

minimal z-position [m]

Definition at line 87 of file JMuonEnergyParameters_t.hh.

resolution

double JRECONSTRUCTION::JMuonEnergyParameters_t::resolution

inherited

energy resolution [log10(GeV)]

Definition at line 88 of file JMuonEnergyParameters_t.hh.

mestimator

int JRECONSTRUCTION::JMuonEnergyParameters_t::mestimator

inherited

M-estimator (see JFIT::JMEstimator_t)

Definition at line 89 of file JMuonEnergyParameters_t.hh.

---

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

• JMuonEnergy.hh