Jpp test-rotations-old
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
134 this->estimator.reset(getMEstimator(mestimator));
135 }
136
137
138 /**
139 * Get input data.
140 *
141 * \param router module router
142 * \param summary summary data
143 * \param event event
144 * \param in start values
145 * \param coverage coverage
146 * \return input data
147 */
149 const JSummaryRouter& summary,
150 const JDAQEvent& event,
151 const JEvt& in,
152 const coverage_type& coverage) const
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 }
184
185
186 /**
187 * Auxiliary class for energy estimation.
188 */
189 struct JResult {
190 /**
191 * Constructor.
192 *
193 * \param x Energy [log(E/GeV)]
194 * \param chi2 chi2
195 */
196 JResult(const JEnergy& x = 0.0,
197 const double chi2 = std::numeric_limits<double>::max())
198 {
199 this->x = x;
200 this->chi2 = chi2;
201 }
202
203 /**
204 * Type conversion.
205 *
206 * \return true if valid chi2; else false
207 */
208 operator bool() const
209 {
210 return chi2 != std::numeric_limits<double>::max();
211 }
212
213 JEnergy x; //!< Energy
214 double chi2; //!< chi2
215 };
216
217
218 /**
219 * Fit function.
220 *
221 * \param input input data
222 * \return fit results
223 */
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
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
254 vector<JNPEHit> data;
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 }
451
453 };
454}
455
456#endif
457
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 JDAQModuleIdentifier &module) const
Get 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 resolution
energy resolution [log10(GeV)]
double EMin_log
minimal energy [log10(GeV)]
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