Jpp 19.3.0-rc.3
the software that should make you happy
Loading...
Searching...
No Matches
JRECONSTRUCTION::JMuonEnergy Class Reference

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

#include <JMuonEnergy.hh>

Inheritance diagram for JRECONSTRUCTION::JMuonEnergy:
JRECONSTRUCTION::JMuonEnergyParameters_t JFIT::JRegressor< JModel_t, JMinimiser_t > TObject

Classes

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

Public Types

typedef JRegressor< JEnergyJRegressor_t
 
typedef JModuleL0 module_type
 
typedef std::vector< module_typedetector_type
 

Public Member Functions

 JMuonEnergy (const JMuonEnergyParameters_t &parameters, const storage_type &storage, 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

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.
 

Detailed Description

Auxiliary class to to determine muon energy.

Definition at line 66 of file JMuonEnergy.hh.

Member Typedef Documentation

◆ JRegressor_t

◆ module_type

◆ detector_type

Constructor & Destructor Documentation

◆ JMuonEnergy()

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

Constructor.

Parameters
parametersparameters
storagestorage
correctenergy correction
debugdebug

Definition at line 118 of file JMuonEnergy.hh.

121 :
122 JMuonEnergyParameters_t(parameters),
123 JRegressor_t(storage),
125 {
126 using namespace JPP;
127
128 if (this->getRmax() < roadWidth_m) {
129 roadWidth_m = this->getRmax();
130 }
131
132 JRegressor_t::debug = debug;
133
134 this->estimator.reset(getMEstimator(mestimator));
135 }
int debug
debug level
Definition JSirene.cc:72
JRegressor< JEnergy > JRegressor_t
JMEstimator * getMEstimator(const int type)
Get M-Estimator.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).

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
routermodule router
summarysummary data
eventevent
instart values
coveragecoverage
Returns
input data

Definition at line 148 of file JMuonEnergy.hh.

153 {
154 using namespace std;
155 using namespace JTRIGGER;
156
157 input_type input(event.getDAQEventHeader(), in, coverage);
158
159 const JBuildL0 <JHitR0> buildL0;
161
162 const JDAQTimeslice timeslice(event, true);
163
164 JSuperFrame2D<JHit> buffer;
165
166 for (JDAQTimeslice::const_iterator i = timeslice.begin(); i != timeslice.end(); ++i) {
167
168 if (router.hasModule(i->getModuleID())) {
169
170 buffer(*i, router.getModule(i->getModuleID()));
171
172 buildL0(buffer, back_inserter(data[i->getModuleID()]));
173 }
174 }
175
176 for (const auto& module : router.getReference()) {
177 if (!module.empty()) {
178 input.data.push_back(module_type(module, summary.getSummaryFrame(module.getID(), R_Hz), data[module.getID()]));
179 }
180 }
181
182 return input;
183 }
bool hasModule(const JObjectID &id) const
Has module.
const JModule & getModule(const JObjectID &id) const
Get module parameters.
const JClass_t & getReference() const
Get reference to object.
Definition JReference.hh:38
const JDAQSummaryFrame & getSummaryFrame(const JDAQModuleIdentifier &module) const
Get summary frame.
2-dimensional frame with time calibrated data from one optical module.
const JDAQEventHeader & getDAQEventHeader() const
Get DAQ event header.
Auxiliary classes and methods for triggering.

◆ operator()()

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

Fit function.

Parameters
inputinput data
Returns
fit results

Definition at line 224 of file JMuonEnergy.hh.

225 {
226 using namespace std;
227 using namespace JPP;
228
229 JEvent event(JMUONENERGY);
230
231 JEvt out;
232
233 // select start values
234
235 JEvt in = input.in;
236
237 in.select(numberOfPrefits, qualitySorter);
238
239 if (!in.empty()) {
240 in.select(JHistory::is_event(in.begin()->getHistory()));
241 }
242
243 for (JEvt::const_iterator track = in.begin(); track != in.end(); ++track) {
244
245 const JRotation3D R (getDirection(*track));
246 const JLine1Z tz(getPosition (*track).rotate(R), track->getT());
247
248 double zmin = numeric_limits<double>::lowest();
249
250 if (track->hasW(JSTART_LENGTH_METRES) && track->getW(JSTART_LENGTH_METRES) > 0.0) {
251 zmin = this->ZMin_m;
252 }
253
255
256 for (const auto& module : input.data) {
257
258 JPosition3D pos(module->getPosition());
259
260 pos.transform(R, tz.getPosition());
261
262 if (pos.getX() <= roadWidth_m) {
263
264 const double z1 = pos.getZ() - pos.getX() / getTanThetaC();
265 const double t1 = tz .getT() + (pos.getZ() + pos.getX() * getKappaC()) * getInverseSpeedOfLight();
266
267 if (z1 >= zmin) {
268
269 for (size_t i = 0; i != module->size(); ++i) {
270
271 if (module.getStatus(i)) {
272
273 const struct {
274
275 bool operator()(const JHitR0& hit) const
276 {
277 return (hit.getPMT() == pmt && T_ns(hit.getT()));
278 }
279
280 const JTimeRange T_ns;
281 const size_t pmt;
282
283 } match = { T_ns + t1, i };
284
285 JPMT pmt = module->getPMT(i);
286
287 pmt.transform(R, tz.getPosition());
288
289 const JNPEHit hit(this->getNPE(pmt, module.frame.getRate(i)), count_if(module.begin(), module.end(), match));
290
291 DEBUG("hit: " << setw(8) << module->getID() << '.' << FILL(2,'0') << i << ' '
292 << FIXED(7,3) << hypot(pmt.getX(), pmt.getY()) << ' '
293 << FIXED(7,3) << module.frame.getRate(i) * 1.0e-3 << ' '
294 << SCIENTIFIC(9,3) << hit.getY0() << ' '
295 << SCIENTIFIC(9,3) << hit.getY1() << ' '
296 << SCIENTIFIC(9,3) << hit.getYA() << ' '
297 << SCIENTIFIC(9,3) << hit.getYB() << ' '
298 << setw(2) << hit.getN() << endl);
299
300 data.push_back(hit);
301 }
302 }
303 }
304 }
305 }
306
307 const int NDF = distance(data.begin(), data.end()) - 1;
308
309 if (NDF >= 0) {
310
311 // 5-point search between given limits
312
313 const int N = 5;
314
315 JResult result[N];
316
317 for (int i = 0; i != N; ++i) {
318 result[i].x = EMin_log + i * (EMax_log - EMin_log) / (N-1);
319 }
320
321 map<double, double> buffer;
322
323 do {
324
325 int j = 0;
326
327 for (int i = 0; i != N; ++i) {
328
329 if (!result[i]) {
330
331 const JEnergy x = result[i].x;
332 const double chi2 = (*this)(x, data.begin(), data.end());
333
334 result[i].chi2 = chi2;
335 buffer[chi2] = x.getE();
336 }
337
338 if (result[i].chi2 < result[j].chi2) {
339 j = i;
340 }
341 }
342
343
344 for (int i = 0; i != N; ++i) {
345 DEBUG(' ' << FIXED(5,2) << result[i].x << ' ' << FIXED(9,3) << result[i].chi2);
346 }
347 DEBUG(endl);
348
349 // squeeze range
350
351 switch (j) {
352
353 case 0:
354 result[4] = result[1];
355 result[2] = result[0];
356 result[0] = JResult(2*result[2].x - result[4].x);
357 break;
358
359 case 1:
360 result[4] = result[2];
361 result[2] = result[1];
362 break;
363
364 case 2:
365 result[0] = result[1];
366 result[4] = result[3];
367 break;
368
369 case 3:
370 result[0] = result[2];
371 result[2] = result[3];
372 break;
373
374 case 4:
375 result[0] = result[3];
376 result[2] = result[4];
377 result[4] = JResult(2*result[2].x - result[0].x);
378 break;
379 }
380
381 result[1] = JResult(0.5 * (result[0].x + result[2].x));
382 result[3] = JResult(0.5 * (result[2].x + result[4].x));
383
384 } while (result[4].x - result[0].x > resolution);
385
386
387 if (result[1].chi2 != result[3].chi2) {
388
389 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);
390 result[2].chi2 = (*this)(result[2].x, data.begin(), data.end());
391
392 }
393
394 const double chi2 = result[2].chi2;
395 const double E = result[2].x.getE();
396
397 // calculate additional variables
398
399 double Emin = numeric_limits<double>::max();
400 double Emax = numeric_limits<double>::lowest();
401
402 for (map<double, double>::const_iterator i = buffer.begin(); i != buffer.end() && i->first <= chi2 + 1.0; ++i) {
403 if (i->second < Emin) { Emin = i->second; }
404 if (i->second > Emax) { Emax = i->second; }
405 }
406
407 const double mu_range = gWater(E); // range of a muon with the reconstructed energy
408
409 double noise_likelihood = 0.0; // log-likelihood of every hit being from K40
410 int number_of_hits = 0; // number of hits selected for JEnergy
411
412 for (vector<JNPEHit>::const_iterator i = data.begin(); i != data.end(); ++i) {
413 noise_likelihood += log10(getP(i->getY0(), i->getN())); // probability of getting the observed multiplicity with K40
414 number_of_hits += i->getN(); // hit multiplicity
415 }
416
417 JFit fit = *track;
418
419 fit.push_back(event());
420
421 // set corrected energy
422
423 fit.setE(correct(E));
424
425 out.push_back(fit);
426
427 // set additional values
428
429 out.rbegin()->setW(track->getW());
430 out.rbegin()->setW(JENERGY_ENERGY, E);
431 out.rbegin()->setW(JENERGY_CHI2, chi2);
432 out.rbegin()->setW(JENERGY_MUON_RANGE_METRES, mu_range);
433 out.rbegin()->setW(JENERGY_NOISE_LIKELIHOOD, noise_likelihood);
434 out.rbegin()->setW(JENERGY_NDF, NDF);
435 out.rbegin()->setW(JENERGY_NUMBER_OF_HITS, number_of_hits);
436 out.rbegin()->setW(JENERGY_MINIMAL_ENERGY, Emin);
437 out.rbegin()->setW(JENERGY_MAXIMAL_ENERGY, Emax);
438 out.rbegin()->setW(JPP_COVERAGE_ORIENTATION, input.coverage.orientation);
439 out.rbegin()->setW(JPP_COVERAGE_POSITION, input.coverage.position);
440 }
441 }
442
443 // apply default sorter
444
445 sort(out.begin(), out.end(), qualitySorter);
446
447 copy(input.in.begin(), input.in.end(), back_inserter(out));
448
449 return out;
450 }
#define DEBUG(A)
Message macros.
Definition JMessage.hh:62
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
Data structure for PMT geometry, calibration and status.
Definition JPMT.hh:49
Data structure for fit of energy.
Definition JEnergy.hh:31
Data structure for fit of straight line paralel to z-axis.
Definition JLine1Z.hh:29
void transform(const JAxis3D &axis)
Transform axis to reference frame of given axis.
Definition JAxis3D.hh:359
Data structure for position in three dimensions.
JPosition3D & rotate(const JRotation3D &R)
Rotate.
double getY() const
Get y position.
Definition JVector3D.hh:104
double getX() const
Get x position.
Definition JVector3D.hh:94
JEvt operator()(const input_type &input)
Fit function.
Reduced data structure for L0 hit.
Definition JHitR0.hh:27
JPMT_t getPMT() const
Get PMT.
Definition JHitR0.hh:60
int getN() const
Get count.
Definition JHitR0.hh:71
double getT() const
Get calibrated time of hit.
static const int JMUONENERGY
static const int JENERGY_NDF
number of degrees of freedom from JEnergy.cc
static const int JENERGY_ENERGY
uncorrected energy [GeV] from JEnergy.cc
static const int JENERGY_NOISE_LIKELIHOOD
log likelihood of every hit being K40 from JEnergy.cc
static const int JENERGY_CHI2
chi2 from JEnergy.cc
static const int JENERGY_MUON_RANGE_METRES
range of a muon with the reconstructed energy [m] from JEnergy.cc
static const int JPP_COVERAGE_POSITION
coverage of dynamic position calibration from any Jpp application
static const int JENERGY_NUMBER_OF_HITS
number of hits from JEnergy.cc
static const int JSTART_LENGTH_METRES
distance between projected positions on the track of optical modules for which the response does not ...
static const int JPP_COVERAGE_ORIENTATION
coverage of dynamic orientation calibration from any Jpp application
static const int JENERGY_MINIMAL_ENERGY
minimal energy [GeV] from JEnergy.cc
static const int JENERGY_MAXIMAL_ENERGY
maximal energy [GeV] from JEnergy.cc
void copy(const Head &from, JHead &to)
Copy header from from to to.
Definition JHead.cc:163
double getNPE(const Hit &hit)
Get true charge of hit.
double getP(const double expval, bool hit)
Get Poisson probability to observe a hit or not for given expectation value for the number of hits.
double getKappaC()
Get average R-dependence of arrival time of Cherenkov light (a.k.a.
const double getInverseSpeedOfLight()
Get inverse speed of light.
static const JGeaneWater gWater
Function object for energy loss of muon in sea water.
Definition JGeane.hh:381
double getTanThetaC()
Get average tangent of Cherenkov angle of water corresponding to group velocity.
JPosition3D getPosition(const JFit &fit)
Get position.
bool qualitySorter(const JFit &first, const JFit &second)
Comparison of fit results.
JDirection3D getDirection(const JFit &fit)
Get direction.
return result
Definition JPolint.hh:862
int j
Definition JPolint.hh:801
Auxiliary data structure for sequence of same character.
Definition JManip.hh:330
Auxiliary data structure for floating point format specification.
Definition JManip.hh:448
Acoustic event fit.
Acoustic single fit.
Auxiliary class to test history.
Definition JHistory.hh:136
Auxiliary class for simultaneously handling light yields and response of PMT.
Definition JNPEHit.hh:21
double resolution
energy resolution [log10(GeV)]
double EMin_log
minimal energy [log10(GeV)]
double EMax_log
maximal energy [log10(GeV)]
Auxiliary data structure for floating point format specification.
Definition JManip.hh:488

◆ reset()

void JRECONSTRUCTION::JMuonEnergyParameters_t::reset ( )
inlineinherited

Reset fit parameters.

Definition at line 44 of file JMuonEnergyParameters_t.hh.

45 {
46 roadWidth_m = std::numeric_limits<double>::max();
47 R_Hz = 6.0e3;
49 EMin_log = 1.0;
50 EMax_log = 8.0;
51 TMin_ns = -50.0;
52 TMax_ns = +450.0;
53 ZMin_m = std::numeric_limits<double>::lowest();
54 resolution = 0.01;
56 }
@ EM_NORMAL
double TMin_ns
minimal time w.r.t. Cherenkov hypothesis [ns]
double TMax_ns
maximal time w.r.t. Cherenkov hypothesis [ns]

◆ equals()

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

Equality.

Parameters
parametersfit parameters
Returns
true if equals; else false

Definition at line 64 of file JMuonEnergyParameters_t.hh.

65 {
66 return (this->roadWidth_m == parameters.roadWidth_m &&
67 this->R_Hz == parameters.R_Hz &&
68 this->numberOfPrefits == parameters.numberOfPrefits &&
69 this->EMin_log == parameters.EMin_log &&
70 this->EMax_log == parameters.EMax_log &&
71 this->TMin_ns == parameters.TMin_ns &&
72 this->TMax_ns == parameters.TMax_ns &&
73 this->ZMin_m == parameters.ZMin_m &&
74 this->resolution == parameters.resolution &&
75 this->mestimator == parameters.mestimator);
76 }

◆ ClassDef()

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

Member Data Documentation

◆ correct

JEnergyCorrection JRECONSTRUCTION::JMuonEnergy::correct

Definition at line 452 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.

Definition at line 89 of file JMuonEnergyParameters_t.hh.


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