Jpp 19.3.0-rc.2
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 JRegressor_t::T_ns.setRange(TMin_ns, TMax_ns);
134
135 this->estimator.reset(getMEstimator(mestimator));
136 }
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).
double TMin_ns
minimal time w.r.t. Cherenkov hypothesis [ns]
double TMax_ns
maximal time w.r.t. Cherenkov hypothesis [ns]

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 149 of file JMuonEnergy.hh.

154 {
155 using namespace std;
156 using namespace JTRIGGER;
157
158 input_type input(event.getDAQEventHeader(), in, coverage);
159
160 const JBuildL0 <JHitR0> buildL0;
162
163 const JDAQTimeslice timeslice(event, true);
164
165 JSuperFrame2D<JHit> buffer;
166
167 for (JDAQTimeslice::const_iterator i = timeslice.begin(); i != timeslice.end(); ++i) {
168
169 if (router.hasModule(i->getModuleID())) {
170
171 buffer(*i, router.getModule(i->getModuleID()));
172
173 buildL0(buffer, back_inserter(data[i->getModuleID()]));
174 }
175 }
176
177 for (const auto& module : router.getReference()) {
178 if (!module.empty()) {
179 input.data.push_back(module_type(module, summary.getSummaryFrame(module.getID()), data[module.getID()]));
180 }
181 }
182
183 return input;
184 }
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
Get default 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 225 of file JMuonEnergy.hh.

226 {
227 using namespace std;
228 using namespace JPP;
229
230 JEvent event(JMUONENERGY);
231
232 JEvt out;
233
234 // select start values
235
236 JEvt in = input.in;
237
238 in.select(numberOfPrefits, qualitySorter);
239
240 if (!in.empty()) {
241 in.select(JHistory::is_event(in.begin()->getHistory()));
242 }
243
244 for (JEvt::const_iterator track = in.begin(); track != in.end(); ++track) {
245
246 const JRotation3D R (getDirection(*track));
247 const JLine1Z tz(getPosition (*track).rotate(R), track->getT());
248
249 double zmin = numeric_limits<double>::lowest();
250
251 if (track->hasW(JSTART_LENGTH_METRES) && track->getW(JSTART_LENGTH_METRES) > 0.0) {
252 zmin = this->ZMin_m;
253 }
254
256
257 for (const auto& module : input.data) {
258
259 JPosition3D pos(module->getPosition());
260
261 pos.transform(R, tz.getPosition());
262
263 if (pos.getX() <= roadWidth_m) {
264
265 const double z1 = pos.getZ() - pos.getX() / getTanThetaC();
266 const double t1 = tz .getT() + (pos.getZ() + pos.getX() * getKappaC()) * getInverseSpeedOfLight();
267
268 if (z1 >= zmin) {
269
270 for (size_t i = 0; i != module->size(); ++i) {
271
272 if (module.getStatus(i)) {
273
274 const struct {
275
276 bool operator()(const JHitR0& hit) const
277 {
278 return (hit.getPMT() == pmt && T_ns(hit.getT()));
279 }
280
281 const JTimeRange T_ns;
282 const size_t pmt;
283
284 } match = { JRegressor_t::T_ns + t1, i };
285
286 JPMT pmt = module->getPMT(i);
287
288 pmt.transform(R, tz.getPosition());
289
290 const JNPEHit hit(this->getNPE(pmt, module.frame.getRate(i)), count_if(module.begin(), module.end(), match));
291
292 DEBUG("hit: " << setw(8) << module->getID() << '.' << FILL(2,'0') << i << ' '
293 << FIXED(7,3) << hypot(pmt.getX(), pmt.getY()) << ' '
294 << FIXED(7,3) << module.frame.getRate(i) * 1.0e-3 << ' '
295 << SCIENTIFIC(9,3) << hit.getY0() << ' '
296 << SCIENTIFIC(9,3) << hit.getY1() << ' '
297 << SCIENTIFIC(9,3) << hit.getYA() << ' '
298 << SCIENTIFIC(9,3) << hit.getYB() << ' '
299 << setw(2) << hit.getN() << endl);
300
301 data.push_back(hit);
302 }
303 }
304 }
305 }
306 }
307
308 const int NDF = distance(data.begin(), data.end()) - 1;
309
310 if (NDF >= 0) {
311
312 // 5-point search between given limits
313
314 const int N = 5;
315
316 JResult result[N];
317
318 for (int i = 0; i != N; ++i) {
319 result[i].x = EMin_log + i * (EMax_log - EMin_log) / (N-1);
320 }
321
322 map<double, double> buffer;
323
324 do {
325
326 int j = 0;
327
328 for (int i = 0; i != N; ++i) {
329
330 if (!result[i]) {
331
332 const JEnergy x = result[i].x;
333 const double chi2 = (*this)(x, data.begin(), data.end());
334
335 result[i].chi2 = chi2;
336 buffer[chi2] = x.getE();
337 }
338
339 if (result[i].chi2 < result[j].chi2) {
340 j = i;
341 }
342 }
343
344
345 for (int i = 0; i != N; ++i) {
346 DEBUG(' ' << FIXED(5,2) << result[i].x << ' ' << FIXED(9,3) << result[i].chi2);
347 }
348 DEBUG(endl);
349
350 // squeeze range
351
352 switch (j) {
353
354 case 0:
355 result[4] = result[1];
356 result[2] = result[0];
357 result[0] = JResult(2*result[2].x - result[4].x);
358 break;
359
360 case 1:
361 result[4] = result[2];
362 result[2] = result[1];
363 break;
364
365 case 2:
366 result[0] = result[1];
367 result[4] = result[3];
368 break;
369
370 case 3:
371 result[0] = result[2];
372 result[2] = result[3];
373 break;
374
375 case 4:
376 result[0] = result[3];
377 result[2] = result[4];
378 result[4] = JResult(2*result[2].x - result[0].x);
379 break;
380 }
381
382 result[1] = JResult(0.5 * (result[0].x + result[2].x));
383 result[3] = JResult(0.5 * (result[2].x + result[4].x));
384
385 } while (result[4].x - result[0].x > resolution);
386
387
388 if (result[1].chi2 != result[3].chi2) {
389
390 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);
391 result[2].chi2 = (*this)(result[2].x, data.begin(), data.end());
392
393 }
394
395 const double chi2 = result[2].chi2;
396 const double E = result[2].x.getE();
397
398 // calculate additional variables
399
400 double Emin = numeric_limits<double>::max();
401 double Emax = numeric_limits<double>::lowest();
402
403 for (map<double, double>::const_iterator i = buffer.begin(); i != buffer.end() && i->first <= chi2 + 1.0; ++i) {
404 if (i->second < Emin) { Emin = i->second; }
405 if (i->second > Emax) { Emax = i->second; }
406 }
407
408 const double mu_range = gWater(E); // range of a muon with the reconstructed energy
409
410 double noise_likelihood = 0.0; // log-likelihood of every hit being from K40
411 int number_of_hits = 0; // number of hits selected for JEnergy
412
413 for (vector<JNPEHit>::const_iterator i = data.begin(); i != data.end(); ++i) {
414 noise_likelihood += log10(getP(i->getY0(), i->getN())); // probability of getting the observed multiplicity with K40
415 number_of_hits += i->getN(); // hit multiplicity
416 }
417
418 JFit fit = *track;
419
420 fit.push_back(event());
421
422 // set corrected energy
423
424 fit.setE(correct(E));
425
426 out.push_back(fit);
427
428 // set additional values
429
430 out.rbegin()->setW(track->getW());
431 out.rbegin()->setW(JENERGY_ENERGY, E);
432 out.rbegin()->setW(JENERGY_CHI2, chi2);
433 out.rbegin()->setW(JENERGY_MUON_RANGE_METRES, mu_range);
434 out.rbegin()->setW(JENERGY_NOISE_LIKELIHOOD, noise_likelihood);
435 out.rbegin()->setW(JENERGY_NDF, NDF);
436 out.rbegin()->setW(JENERGY_NUMBER_OF_HITS, number_of_hits);
437 out.rbegin()->setW(JENERGY_MINIMAL_ENERGY, Emin);
438 out.rbegin()->setW(JENERGY_MAXIMAL_ENERGY, Emax);
439 out.rbegin()->setW(JPP_COVERAGE_ORIENTATION, input.coverage.orientation);
440 out.rbegin()->setW(JPP_COVERAGE_POSITION, input.coverage.position);
441 }
442 }
443
444 // apply default sorter
445
446 sort(out.begin(), out.end(), qualitySorter);
447
448 copy(input.in.begin(), input.in.end(), back_inserter(out));
449
450 return out;
451 }
#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

◆ 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 453 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: