Jpp 19.3.0-rc.1
the software that should make you happy
Loading...
Searching...
No Matches
JMuonEnergy.hh
Go to the documentation of this file.
1#ifndef __JRECONSTRUCTION__JMUONENERGY__
2#define __JRECONSTRUCTION__JMUONENERGY__
3
4#include <string>
5#include <istream>
6#include <ostream>
7#include <iomanip>
8#include <set>
9#include <vector>
10#include <algorithm>
11#include <cmath>
12
16
20
22
23#include "JTrigger/JHitR0.hh"
24#include "JTrigger/JBuildL0.hh"
25
27
28#include "JFit/JLine1Z.hh"
30#include "JFit/JFitToolkit.hh"
31#include "JFit/JModel.hh"
32#include "JFit/JNPEHit.hh"
33#include "JFit/JEnergy.hh"
34
41
42#include "JPhysics/JGeane.hh"
43
44
45/**
46 * \author mdejong, azegarelli, scelli
47 */
48namespace JRECONSTRUCTION {}
49namespace JPP { using namespace JRECONSTRUCTION; }
50
51namespace JRECONSTRUCTION {
52
58 using JFIT::JRegressor;
59 using JFIT::JEnergy;
61
62
63 /**
64 * Auxiliary class to to determine muon energy.
65 */
66 class JMuonEnergy :
68 public JRegressor<JEnergy>
69 {
70 public:
74
75 using JRegressor_t::operator();
76
77 /**
78 * Input data type.
79 */
80 struct input_type :
81 public JDAQEventHeader
82 {
83 /**
84 * Default constructor.
85 */
87 {}
88
89
90 /**
91 * Constructor.
92 *
93 * \param header header
94 * \param in start values
95 * \param coverage coverage
96 */
98 const JEvt& in,
99 const coverage_type& coverage) :
100 JDAQEventHeader(header),
101 in(in)
102 {}
103
107 };
108
109
110 /**
111 * Constructor.
112 *
113 * \param parameters parameters
114 * \param storage storage
115 * \param correct energy correction
116 * \param debug debug
117 */
119 const storage_type& storage,
121 const int debug = 0):
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 }
137
138
139 /**
140 * Get input data.
141 *
142 * \param router module router
143 * \param summary summary data
144 * \param event event
145 * \param in start values
146 * \param coverage coverage
147 * \return input data
148 */
150 const JSummaryRouter& summary,
151 const JDAQEvent& event,
152 const JEvt& in,
153 const coverage_type& coverage) const
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 }
185
186
187 /**
188 * Auxiliary class for energy estimation.
189 */
190 struct JResult {
191 /**
192 * Constructor.
193 *
194 * \param x Energy [log(E/GeV)]
195 * \param chi2 chi2
196 */
197 JResult(const JEnergy& x = 0.0,
198 const double chi2 = std::numeric_limits<double>::max())
199 {
200 this->x = x;
201 this->chi2 = chi2;
202 }
203
204 /**
205 * Type conversion.
206 *
207 * \return true if valid chi2; else false
208 */
209 operator bool() const
210 {
211 return chi2 != std::numeric_limits<double>::max();
212 }
213
214 JEnergy x; //!< Energy
215 double chi2; //!< chi2
216 };
217
218
219 /**
220 * Fit function.
221 *
222 * \param input input data
223 * \return fit results
224 */
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
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
255 vector<JNPEHit> data;
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 }
452
454 };
455}
456
457#endif
458
Coverage of dynamical detector calibration.
Auxiliary class to extract a subset of optical modules from a detector.
Data structure for detector geometry and calibration.
Data regression method for JFIT::JEnergy.
Auxiliary methods to evaluate Poisson probabilities and chi2.
Energy loss of muon.
Basic data structure for L0 hit.
#define DEBUG(A)
Message macros.
Definition JMessage.hh:62
int debug
debug level
Definition JSirene.cc:72
Direct access to module in detector data structure.
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
Router for direct addressing of module data in detector data structure.
bool hasModule(const JObjectID &id) const
Has module.
const JModule & getModule(const JObjectID &id) const
Get module parameters.
Data structure for a composite optical module.
Definition JModule.hh:75
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 set of track fit results.
void select(const JSelector_t &selector)
Select fits.
Data structure for track fit results with history and optional associated values.
void setE(const double E)
Set energy.
Data structure for fit of straight line paralel to z-axis.
Definition JLine1Z.hh:29
double getT(const JVector3D &pos) const
Get arrival time of Cherenkov light at given position.
Definition JLine1Z.hh:114
void transform(const JAxis3D &axis)
Transform axis to reference frame of given axis.
Definition JAxis3D.hh:359
Data structure for position in three dimensions.
void transform(const JRotation3D &R, const JVector3D &pos)
Transform position.
JPosition3D & rotate(const JRotation3D &R)
Rotate.
const JPosition3D & getPosition() const
Get position.
double getY() const
Get y position.
Definition JVector3D.hh:104
double getZ() const
Get z position.
Definition JVector3D.hh:115
double getX() const
Get x position.
Definition JVector3D.hh:94
const JClass_t & getReference() const
Get reference to object.
Definition JReference.hh:38
Auxiliary class for correction of energy determined by JEnergy.cc.
Auxiliary class to to determine muon energy.
JMuonEnergy(const JMuonEnergyParameters_t &parameters, const storage_type &storage, const JEnergyCorrection &correct, const int debug=0)
Constructor.
JEvt operator()(const input_type &input)
Fit function.
JRegressor< JEnergy > JRegressor_t
std::vector< module_type > detector_type
input_type getInput(const JModuleRouter &router, const JSummaryRouter &summary, const JDAQEvent &event, const JEvt &in, const coverage_type &coverage) const
Get input data.
Router for fast addressing of summary data in KM3NETDAQ::JDAQSummaryslice data structure as a functio...
const JDAQSummaryFrame & getSummaryFrame() const
Get default summary frame.
Reduced data structure for L0 hit.
Definition JHitR0.hh:27
JPMT_t getPMT() const
Get PMT.
Definition JHitR0.hh:60
double getT() const
Get calibrated time of hit.
2-dimensional frame with time calibrated data from one optical module.
const JDAQEventHeader & getDAQEventHeader() const
Get DAQ event header.
Data storage class for rate measurements of all PMTs in one module.
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.
JMEstimator * getMEstimator(const int type)
Get M-Estimator.
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.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Model fits to data.
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 classes and methods for triggering.
Auxiliary data structure for sequence of same character.
Definition JManip.hh:330
Auxiliary data structure for floating point format specification.
Definition JManip.hh:448
Data structure for coverage of detector by dynamical calibrations.
Definition JCoverage.hh:19
double position
coverage of detector by available position calibration [0,1]
Definition JCoverage.hh:42
double orientation
coverage of detector by available orientation calibration [0,1]
Definition JCoverage.hh:41
Auxiliary class for historical event.
Definition JHistory.hh:36
Auxiliary class to test history.
Definition JHistory.hh:136
double getY0() const
Get light yield due to random background.
Definition JK40.hh:49
Auxiliary class for simultaneously handling light yields and response of PMT.
Definition JNPEHit.hh:21
size_t getN() const
Get number of hits.
Definition JNPEHit.hh:52
double getYB() const
Get light yield due to bremsstrahlung.
Definition JNPE.hh:94
double getY1() const
Get light yield due to minimum ionizing particle.
Definition JNPE.hh:72
double getYA() const
Get light yield due to delta-rays.
Definition JNPE.hh:83
Template definition of a data regressor of given model.
Definition JRegressor.hh:70
Auxiliary class for handling module response.
Definition JModuleL0.hh:45
double TMin_ns
minimal time w.r.t. Cherenkov hypothesis [ns]
double resolution
energy resolution [log10(GeV)]
double EMin_log
minimal energy [log10(GeV)]
double TMax_ns
maximal time w.r.t. Cherenkov hypothesis [ns]
double EMax_log
maximal energy [log10(GeV)]
Auxiliary class for energy estimation.
JResult(const JEnergy &x=0.0, const double chi2=std::numeric_limits< double >::max())
Constructor.
input_type(const JDAQEventHeader &header, const JEvt &in, const coverage_type &coverage)
Constructor.
Auxiliary data structure for floating point format specification.
Definition JManip.hh:488