Jpp 19.3.0-rc.1
the software that should make you happy
Loading...
Searching...
No Matches
JMuonStart.hh
Go to the documentation of this file.
1#ifndef __JRECONSTRUCTION__JMUONSTART__
2#define __JRECONSTRUCTION__JMUONSTART__
3
4#include <string>
5#include <istream>
6#include <ostream>
7#include <iomanip>
8#include <vector>
9#include <set>
10#include <map>
11#include <algorithm>
12#include <cmath>
13
17
21
23
24#include "JTrigger/JHitL0.hh"
25#include "JTrigger/JHitR1.hh"
26#include "JTrigger/JBuildL0.hh"
27
29
30#include "JFit/JLine1Z.hh"
32#include "JFit/JFitToolkit.hh"
33#include "JFit/JModel.hh"
34#include "JFit/JNPEHit.hh"
35#include "JFit/JEnergy.hh"
36#include "JFit/JMEstimator.hh"
37
43
45#include "JPhysics/JK40Rates.hh"
46
47#include "JLang/JComparator.hh"
48#include "JLang/JPredicate.hh"
49
50
51/**
52 * \author mdejong, azegarelli, scelli
53 */
54namespace JRECONSTRUCTION {}
55namespace JPP { using namespace JRECONSTRUCTION; }
56
57namespace JRECONSTRUCTION {
58
64 using JFIT::JRegressor;
65 using JFIT::JEnergy;
68
69
70 /**
71 * Auxiliary class to determine start and end position of muon trajectory.
72 */
73 class JMuonStart :
75 public JRegressor<JEnergy>
76 {
77 public:
81
82 using JRegressor_t::operator();
83
84 /**
85 * Input data type.
86 */
87 struct input_type :
88 public JDAQEventHeader
89 {
90 /**
91 * Default constructor.
92 */
96
97
98 /**
99 * Constructor.
100 *
101 * \param header header
102 * \param in start values
103 * \param coverage coverage
104 */
106 const JEvt& in,
107 const coverage_type& coverage) :
108 JDAQEventHeader(header),
109 in(in),
111 {}
112
116 };
117
118
119 /**
120 * Constructor.
121 *
122 * \param parameters parameters
123 * \param storage storage
124 * \param rates_Hz K40 rates [Hz]
125 * \param debug debug
126 */
128 const storage_type& storage,
129 const JK40Rates& rates_Hz,
130 const int debug = 0):
131 JMuonStartParameters_t(parameters),
132 JRegressor_t(storage),
134 {
135 if (this->getRmax() < parameters.roadWidth_m) {
136 roadWidth_m = this->getRmax();
137 }
138
139 JRegressor_t::debug = debug;
140 JRegressor_t::T_ns.setRange(parameters.TMin_ns, parameters.TMax_ns);
141 }
142
143
144 /**
145 * Get input data.
146 *
147 * \param router module router
148 * \param summary summary data
149 * \param event event
150 * \param in start values
151 * \param coverage coverage
152 * \return input data
153 */
155 const JSummaryRouter& summary,
156 const JDAQEvent& event,
157 const JEvt& in,
158 const coverage_type& coverage)
159 {
160 using namespace std;
161 using namespace JTRIGGER;
162
163 input_type input(event.getDAQEventHeader(), in, coverage);
164
165 const JBuildL0 <JHitR0> buildL0;
167
168 const JDAQTimeslice timeslice(event, true);
169
170 JSuperFrame2D<JHit> buffer;
171
172 for (JDAQTimeslice::const_iterator i = timeslice.begin(); i != timeslice.end(); ++i) {
173
174 if (router.hasModule(i->getModuleID())) {
175
176 buffer(*i, router.getModule(i->getModuleID()));
177
178 buildL0(buffer, back_inserter(data[i->getModuleID()]));
179 }
180 }
181
182 for (const auto& module : router.getReference()) {
183 if (!module.empty()) {
184 input.data.push_back(module_type(module, summary.getSummaryFrame(module.getID()), data[module.getID()]));
185 }
186 }
187
188 return input;
189 }
190
191
192 /**
193 * Fit function.
194 *
195 * \param input input data
196 * \return fit results
197 */
199 {
200 using namespace std;
201 using namespace JPP;
202
203 JEvent event(JMUONSTART);
204
205 JEvt out;
206
207 // select start values
208
209 JEvt in = input.in;
210
212
213 if (!in.empty()) {
214 in.select(JHistory::is_event(in.begin()->getHistory()));
215 }
216
217 const JStart start(Pmin1, Pmin2, Nmax2);
218
219 struct JHit_t {
220
221 double getZ() const { return z; }
222 double getP() const { return p; }
223
224 double z;
225 double p;
226 double y1;
227 double yx;
228 };
229
230 for (JEvt::const_iterator track = in.begin(); track != in.end(); ++track) {
231
232 const JRotation3D R (getDirection(*track));
233 const JLine1Z tz(getPosition (*track).rotate(R), track->getT());
234
235 vector<JHit_t> data;
236
237 for (const auto& module : input.data) {
238
239 JPosition3D pos(module->getPosition());
240
241 pos.transform(R, tz.getPosition());
242
243 if (pos.getX() <= roadWidth_m) {
244
245 const double z = pos.getZ() - pos.getX() / getTanThetaC();
246 const double t = tz .getT() + (pos.getZ() + pos.getX() * getKappaC()) * getInverseSpeedOfLight();
247
248 const double p = module.getProbability(rates_Hz.getMultiplesRates(), JRegressor_t::T_ns + t);
249
250 double y1 = 0.0;
251 double yx = 0.0;
252
253 for (size_t i = 0; i != module->size(); ++i) {
254
255 if (module.getStatus(i)) {
256
257 JPMT pmt = module->getPMT(i);
258
259 pmt.transform(R, tz.getPosition());
260
261 const double npe = this->getY1(pmt);
262
263 y1 += npe;
264
265 if (count_if(module.begin(), module.end(), make_predicate(&JHitR0::getPMT, (JDAQHit::JPMT_t) i)) == 0) {
266 yx += npe;
267 }
268 }
269 }
270
271 data.push_back({ z, p, y1, yx });
272 }
273 }
274
275 double Zmin = 0.0; // start position
276 double Zmax = 0.0; // end position
277 double y1 = 0.0; // npe due to MIP between start and end position
278 double yx = 0.0; // npe due to MIP between start and end position that are not detected
279
280 if (!data.empty()) {
281
282 sort(data.begin(), data.end(), make_comparator(&JHit_t::getZ));
283
284 vector<JHit_t>::const_iterator q1 = start.find(data. begin(), data. end());
285 vector<JHit_t>::const_reverse_iterator q2 = start.find(data.rbegin(), data.rend());
286
287 if (q1 != data.end() && q2 != data.rend()) {
288
289 vector<JHit_t>::const_iterator p1 = q1; if (p1 != data. begin()) { --p1; }
290 vector<JHit_t>::const_reverse_iterator p2 = q2; if (p2 != data.rbegin()) { --p2; }
291
292 Zmin = 0.5 * (p1->getZ() + q1->getZ());
293 Zmax = 0.5 * (p2->getZ() + q2->getZ());
294
295 for (vector<JHit_t>::const_iterator i = q1; i != q2.base(); ++i) {
296 y1 += i->y1;
297 yx += i->yx;
298 }
299 }
300 }
301
302 JFit fit = *track;
303
304 fit.push_back(event());
305
306 // move track
307
308 fit.move(Zmin, getSpeedOfLight());
309
310 out.push_back(fit);
311
312 // set additional values
313
314 out.rbegin()->setW(track->getW());
315 out.rbegin()->setW(JSTART_NPE_MIP_TOTAL, y1);
316 out.rbegin()->setW(JSTART_NPE_MIP_MISSED, yx);
317 out.rbegin()->setW(JSTART_LENGTH_METRES, Zmax - Zmin);
318 out.rbegin()->setW(JPP_COVERAGE_ORIENTATION, input.coverage.orientation);
319 out.rbegin()->setW(JPP_COVERAGE_POSITION, input.coverage.position);
320 }
321
322 // apply default sorter
323
324 sort(out.begin(), out.end(), qualitySorter);
325
326 copy(input.in.begin(), input.in.end(), back_inserter(out));
327
328 return out;
329 }
330
331 private:
333 };
334}
335
336#endif
Coverage of dynamical detector calibration.
Auxiliary class to extract a subset of optical modules from a detector.
Data structure for detector geometry and calibration.
TPaveText * p1
Data regression method for JFIT::JEnergy.
Auxiliary methods to evaluate Poisson probabilities and chi2.
Basic data structure for L0 hit.
Reduced data structure for L1 hit.
Maximum likelihood estimator (M-estimators).
int debug
debug level
Definition JSirene.cc:72
Direct access to module in detector data structure.
Physics constants.
Auxiliary method to locate start and end point of muon trajectory.
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 move(const double step, const double velocity)
Move vertex along this track with given velocity.
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 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 to determine start and end position of muon trajectory.
Definition JMuonStart.hh:76
std::vector< module_type > detector_type
Definition JMuonStart.hh:80
JEvt operator()(const input_type &input)
Fit function.
JRegressor< JEnergy > JRegressor_t
Definition JMuonStart.hh:78
input_type getInput(const JModuleRouter &router, const JSummaryRouter &summary, const JDAQEvent &event, const JEvt &in, const coverage_type &coverage)
Get input data.
JMuonStart(const JMuonStartParameters_t &parameters, const storage_type &storage, const JK40Rates &rates_Hz, const int debug=0)
Constructor.
Router for fast addressing of summary data in KM3NETDAQ::JDAQSummaryslice data structure as a functio...
const JDAQSummaryFrame & getSummaryFrame() const
Get default summary frame.
JPMT_t getPMT() const
Get PMT.
Definition JHitR0.hh:60
2-dimensional frame with time calibrated data from one optical module.
const JDAQEventHeader & getDAQEventHeader() const
Get DAQ event header.
unsigned char JPMT_t
PMT channel in FPGA.
Definition JDAQHit.hh:38
Data storage class for rate measurements of all PMTs in one module.
static const int JMUONSTART
static const int JSTART_NPE_MIP_TOTAL
number of photo-electrons along the whole track from JStart.cc
static const int JPP_COVERAGE_POSITION
coverage of dynamic position calibration from any Jpp application
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 JSTART_NPE_MIP_MISSED
number of photo-electrons missed from JStart.cc
void copy(const Head &from, JHead &to)
Copy header from from to to.
Definition JHead.cc:163
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.
JComparator< JResult_t T::*, JComparison::lt > make_comparator(JResult_t T::*member)
Helper method to create comparator between values of data member.
double getKappaC()
Get average R-dependence of arrival time of Cherenkov light (a.k.a.
const double getInverseSpeedOfLight()
Get inverse speed of light.
double getTanThetaC()
Get average tangent of Cherenkov angle of water corresponding to group velocity.
const double getSpeedOfLight()
Get speed of light.
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.
Auxiliary classes and methods for triggering.
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
Template definition of a data regressor of given model.
Definition JRegressor.hh:70
Auxiliary class for K40 rates.
Definition JK40Rates.hh:41
Auxiliary class for handling module response.
Definition JModuleL0.hh:45
Data structure for fit parameters.
int Nmax2
maximal number for twofold observations
double Pmin1
minimal probability single observation
double TMin_ns
minimal time w.r.t. Cherenkov hypothesis [ns]
double Pmin2
minimal probability for twofold observations
double TMax_ns
maximal time w.r.t. Cherenkov hypothesis [ns]
input_type(const JDAQEventHeader &header, const JEvt &in, const coverage_type &coverage)
Constructor.
Auxiliary class for start or end point evaluation.
Definition JStart.hh:21
T find(T __begin, T __end) const
Get start point of muon trajectory.
Definition JStart.hh:77
Auxiliary class to set-up Hit.
Definition JSirene.hh:58