Jpp test-rotations-new
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 }
141
142
143 /**
144 * Get input data.
145 *
146 * \param router module router
147 * \param summary summary data
148 * \param event event
149 * \param in start values
150 * \param coverage coverage
151 * \return input data
152 */
154 const JSummaryRouter& summary,
155 const JDAQEvent& event,
156 const JEvt& in,
157 const coverage_type& coverage) const
158 {
159 using namespace std;
160 using namespace JTRIGGER;
161
162 input_type input(event.getDAQEventHeader(), in, coverage);
163
164 const JBuildL0 <JHitR0> buildL0;
166
167 const JDAQTimeslice timeslice(event, true);
168
169 JSuperFrame2D<JHit> buffer;
170
171 for (JDAQTimeslice::const_iterator i = timeslice.begin(); i != timeslice.end(); ++i) {
172
173 if (router.hasModule(i->getModuleID())) {
174
175 buffer(*i, router.getModule(i->getModuleID()));
176
177 buildL0(buffer, back_inserter(data[i->getModuleID()]));
178 }
179 }
180
181 for (const auto& module : router.getReference()) {
182 if (!module.empty()) {
183 input.data.push_back(module_type(module, summary.getSummaryFrame(module.getID(), R_Hz), data[module.getID()]));
184 }
185 }
186
187 return input;
188 }
189
190
191 /**
192 * Fit function.
193 *
194 * \param input input data
195 * \return fit results
196 */
198 {
199 using namespace std;
200 using namespace JPP;
201
202 JEvent event(JMUONSTART);
203
204 JEvt out;
205
206 // select start values
207
208 JEvt in = input.in;
209
211
212 if (!in.empty()) {
213 in.select(JHistory::is_event(in.begin()->getHistory()));
214 }
215
216 const JStart start(Pmin1, Pmin2, Nmax2);
217
218 struct JHit_t {
219
220 double getZ() const { return z; }
221 double getP() const { return p; }
222
223 double z;
224 double p;
225 double y1;
226 double yx;
227 };
228
229 for (JEvt::const_iterator track = in.begin(); track != in.end(); ++track) {
230
231 const JRotation3D R (getDirection(*track));
232 const JLine1Z tz(getPosition (*track).rotate(R), track->getT());
233
234 vector<JHit_t> data;
235
236 for (const auto& module : input.data) {
237
238 JPosition3D pos(module->getPosition());
239
240 pos.transform(R, tz.getPosition());
241
242 if (pos.getX() <= roadWidth_m) {
243
244 const double z = pos.getZ() - pos.getX() / getTanThetaC();
245 const double t = tz .getT() + (pos.getZ() + pos.getX() * getKappaC()) * getInverseSpeedOfLight();
246
247 const double p = module.getProbability(rates_Hz.getMultiplesRates(), T_ns + t);
248
249 double y1 = 0.0;
250 double yx = 0.0;
251
252 for (size_t i = 0; i != module->size(); ++i) {
253
254 if (module.getStatus(i)) {
255
256 JPMT pmt = module->getPMT(i);
257
258 pmt.transform(R, tz.getPosition());
259
260 const double npe = this->getY1(pmt);
261
262 y1 += npe;
263
264 if (count_if(module.begin(), module.end(), make_predicate(&JHitR0::getPMT, (JDAQHit::JPMT_t) i)) == 0) {
265 yx += npe;
266 }
267 }
268 }
269
270 data.push_back({ z, p, y1, yx });
271 }
272 }
273
274 double Zmin = 0.0; // start position
275 double Zmax = 0.0; // end position
276 double y1 = 0.0; // npe due to MIP between start and end position
277 double yx = 0.0; // npe due to MIP between start and end position that are not detected
278
279 if (!data.empty()) {
280
281 sort(data.begin(), data.end(), make_comparator(&JHit_t::getZ));
282
283 vector<JHit_t>::const_iterator q1 = start.find(data. begin(), data. end());
284 vector<JHit_t>::const_reverse_iterator q2 = start.find(data.rbegin(), data.rend());
285
286 if (q1 != data.end() && q2 != data.rend()) {
287
288 vector<JHit_t>::const_iterator p1 = q1; if (p1 != data. begin()) { --p1; }
289 vector<JHit_t>::const_reverse_iterator p2 = q2; if (p2 != data.rbegin()) { --p2; }
290
291 Zmin = 0.5 * (p1->getZ() + q1->getZ());
292 Zmax = 0.5 * (p2->getZ() + q2->getZ());
293
294 for (vector<JHit_t>::const_iterator i = q1; i != q2.base(); ++i) {
295 y1 += i->y1;
296 yx += i->yx;
297 }
298 }
299 }
300
301 JFit fit = *track;
302
303 fit.push_back(event());
304
305 // move track
306
307 fit.move(Zmin, getSpeedOfLight());
308
309 out.push_back(fit);
310
311 // set additional values
312
313 out.rbegin()->setW(track->getW());
314 out.rbegin()->setW(JSTART_NPE_MIP_TOTAL, y1);
315 out.rbegin()->setW(JSTART_NPE_MIP_MISSED, yx);
316 out.rbegin()->setW(JSTART_LENGTH_METRES, Zmax - Zmin);
317 out.rbegin()->setW(JPP_COVERAGE_ORIENTATION, input.coverage.orientation);
318 out.rbegin()->setW(JPP_COVERAGE_POSITION, input.coverage.position);
319 }
320
321 // apply default sorter
322
323 sort(out.begin(), out.end(), qualitySorter);
324
325 copy(input.in.begin(), input.in.end(), back_inserter(out));
326
327 return out;
328 }
329
330 private:
332 };
333}
334
335#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
input_type getInput(const JModuleRouter &router, const JSummaryRouter &summary, const JDAQEvent &event, const JEvt &in, const coverage_type &coverage) const
Get input data.
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
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 JDAQModuleIdentifier &module) const
Get 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 Pmin2
minimal probability for twofold observations
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