Jpp 19.3.0-rc.2
the software that should make you happy
Loading...
Searching...
No Matches
JMuonGandalf.hh
Go to the documentation of this file.
1#ifndef __JRECONSTRUCTION__JMUONGANDALF__
2#define __JRECONSTRUCTION__JMUONGANDALF__
3
4#include <string>
5#include <iostream>
6#include <iomanip>
7#include <vector>
8#include <algorithm>
9
13
14#include "JTrigger/JHitL0.hh"
15#include "JTrigger/JBuildL0.hh"
16
17#include "JFit/JFitToolkit.hh"
18#include "JFit/JLine1Z.hh"
19#include "JFit/JLine3Z.hh"
20#include "JFit/JModel.hh"
21#include "JFit/JGandalf.hh"
23
28
29#include "JPhysics/JPDFTable.hh"
31#include "JPhysics/JGeane.hh"
32
34
36
38
40
41#include "JTools/JRange.hh"
42
43#include "Jeep/JMessage.hh"
44
45
46/**
47 * \author mdejong, gmaggi
48 */
49
50namespace JRECONSTRUCTION {}
51namespace JPP { using namespace JRECONSTRUCTION; }
52
53namespace JRECONSTRUCTION {
54
58 using JFIT::JRegressor;
59 using JFIT::JLine3Z;
60 using JFIT::JGandalf;
62
63
64 /**
65 * Wrapper class to make final fit of muon trajectory.
66 *
67 * The JMuonGandalf fit uses one or more start values (usually taken from the output of JMuonSimplex).\n
68 * All hits of which the PMT position lies within a set road width (JMuonGandalfParameters_t::roadWidth_m) and
69 * time is within a set window (JMuonGandalfParameters_t::TMin_ns, JMuonGandalfParameters_t::TMax_ns) around the Cherenkov hypothesis are taken.\n
70 * In case there are multiple hits from the same PMT is the specified window,
71 * the first hit is taken and the other hits are discarded.\n
72 * The PDF is accordingly evaluated, i.e.\ the normalised probability for a first hit at the given time of the hit is taken.
73 * The normalisation is consistently based on the specified time window.\n
74 * Note that this hit selection is unbiased with respect to the PDF of a single PMT.
75 */
76 struct JMuonGandalf :
78 public JRegressor<JLine3Z, JGandalf>
79 {
83
84 using JRegressor_t::operator();
85
86 /**
87 * Input data type.
88 */
89 struct input_type :
90 public JDAQEventHeader
91 {
92 /**
93 * Default constructor.
94 */
96 {}
97
98
99 /**
100 * Constructor.
101 *
102 * \param header header
103 * \param in start values
104 * \param coverage coverage
105 */
107 const JEvt& in,
108 const coverage_type& coverage) :
109 JDAQEventHeader(header),
110 in(in),
112 {}
113
117 };
118
119
120 /**
121 * Constructor
122 *
123 * \param parameters parameters
124 * \param storage storage
125 * \param debug debug
126 */
128 const storage_type& storage,
129 const int debug = 0) :
130 JMuonGandalfParameters_t(parameters),
131 JRegressor_t(storage)
132 {
133 if (this->getRmax() < roadWidth_m) {
134 roadWidth_m = this->getRmax();
135 }
136
137 JRegressor_t::debug = debug;
138 JRegressor_t::T_ns.setRange(TMin_ns, TMax_ns);
139 JRegressor_t::Vmax_npe = VMax_npe;
140 JRegressor_t::MAXIMUM_ITERATIONS = NMax;
141
142 this->parameters.resize(5);
143
144 this->parameters[0] = JLine3Z::pX();
145 this->parameters[1] = JLine3Z::pY();
146 this->parameters[2] = JLine3Z::pT();
147 this->parameters[3] = JLine3Z::pDX();
148 this->parameters[4] = JLine3Z::pDY();
149 }
150
151
152 /**
153 * Get input data.
154 *
155 * \param router module router
156 * \param summary summary data
157 * \param event event
158 * \param in start values
159 * \param coverage coverage
160 * \return input data
161 */
163 const JSummaryRouter& summary,
164 const JDAQEvent& event,
165 const JEvt& in,
166 const coverage_type& coverage) const
167 {
168 using namespace std;
169 using namespace JTRIGGER;
170
171 const JBuildL0<JHitL0> buildL0;
172
173 input_type input(event.getDAQEventHeader(), in, coverage);
174
175 vector<JHitL0> data;
176
177 buildL0(event, router, true, back_inserter(data));
178
179 for (const auto& hit : data) {
180 input.data.push_back(hit_type(hit, summary.getRate(hit.getPMTIdentifier())));
181 }
182
183 return input;
184 }
185
186
187 /**
188 * Fit function.
189 *
190 * \param input input data
191 * \return fit results
192 */
194 {
195 using namespace std;
196 using namespace JFIT;
197 using namespace JGEOMETRY3D;
198
199 JEvent event(JMUONGANDALF);
200
201 JEvt out;
202
203 const buffer_type& data = input.data;
204
205 // select start values
206
207 JEvt in = input.in;
208
210
211 if (!in.empty()) {
212 in.select(JHistory::is_event(in.begin()->getHistory()));
213 }
214
215 for (JEvt::const_iterator track = in.begin(); track != in.end(); ++track) {
216
217 const JRotation3D R (getDirection(*track));
218 const JLine1Z tz(getPosition (*track).rotate(R), track->getT());
219 JRange<double> Z_m;
220
221 if (track->hasW(JSTART_LENGTH_METRES) &&
222 track->getW(JSTART_LENGTH_METRES) > 0.0) {
223 Z_m = JZRange(ZMin_m, ZMax_m + track->getW(JSTART_LENGTH_METRES));
224 }
225
226 const JModel<JLine1Z> match(tz, roadWidth_m, JRegressor_t::T_ns, Z_m);
227
228 // hit selection based on start value
229
230 buffer_type buffer;
231
232 for (buffer_type::const_iterator i = data.begin(); i != data.end(); ++i) {
233
234 hit_type hit(*i);
235
236 hit.rotate(R);
237
238 if (match(hit)) {
239 buffer.push_back(hit);
240 }
241 }
242
243 // select first hit
244
245 sort(buffer.begin(), buffer.end(), JHitL0::compare);
246
247 buffer_type::iterator __end = unique(buffer.begin(), buffer.end(), equal_to<JDAQPMTIdentifier>());
248
249
250 const int NDF = distance(buffer.begin(), __end) - this->parameters.size();
251
252 if (NDF > 0) {
253
254 // set fit parameters
255
256 if (track->getE() > 0.1)
257 JRegressor_t::E_GeV = track->getE();
258 else
259 JRegressor_t::E_GeV = this->JMuonGandalfParameters_t::E_GeV;
260
261 const double chi2 = (*this)(JLine3Z(tz), buffer.begin(), __end);
262
263 // check error matrix
264
265 bool status = true;
266
267 for (size_t i = 0; i != this->V.size(); ++i) {
268 if (std::isnan(this->V(i,i)) || this->V(i,i) < 0.0) {
269 status = false;
270 }
271 }
272
273 if (status) {
274
275 JTrack3D tb(this->value);
276
277 tb.rotate_back(R);
278
279 out.push_back(getFit(JHistory(track->getHistory(), event()), tb, getQuality(chi2), NDF));
280
281 // set additional values
282
283 out.rbegin()->setV(this->V.size(), this->V);
284
285 out.rbegin()->setW(track->getW());
286 out.rbegin()->setW(JGANDALF_BETA0_RAD, sqrt(this->error.getDX() * this->error.getDX() +
287 this->error.getDY() * this->error.getDY()));
288 out.rbegin()->setW(JGANDALF_BETA1_RAD, sqrt(this->error.getDX() * this->error.getDY()));
289 out.rbegin()->setW(JGANDALF_CHI2, chi2);
290 out.rbegin()->setW(JGANDALF_NUMBER_OF_HITS, distance(buffer.begin(), __end));
291 out.rbegin()->setW(JGANDALF_LAMBDA, this->lambda);
292 out.rbegin()->setW(JGANDALF_NUMBER_OF_ITERATIONS, this->numberOfIterations);
293 out.rbegin()->setW(JPP_COVERAGE_ORIENTATION, input.coverage.orientation);
294 out.rbegin()->setW(JPP_COVERAGE_POSITION, input.coverage.position);
295 }
296 }
297 }
298
299 // apply default sorter
300
301 sort(out.begin(), out.end(), qualitySorter);
302
303 copy(input.in.begin(), input.in.end(), back_inserter(out));
304
305 return out;
306 }
307 };
308}
309
310#endif
311
Coverage of dynamical detector calibration.
Auxiliary methods to evaluate Poisson probabilities and chi2.
Energy loss of muon.
Basic data structure for L0 hit.
Data regression method for JFIT::JLine3Z.
General purpose messaging.
int debug
debug level
Definition JSirene.cc:72
Direct access to module in detector data structure.
Auxiliary methods for PDF calculations.
Auxiliary class to define a range between two values.
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.
Data structure for set of track fit results.
void select(const JSelector_t &selector)
Select fits.
Fit method based on the Levenberg-Marquardt method.
Definition JGandalf.hh:87
Data structure for fit of straight line paralel to z-axis.
Definition JLine1Z.hh:29
static parameter_type pY()
Definition JLine1Z.hh:181
static parameter_type pX()
Definition JLine1Z.hh:180
static parameter_type pT()
Definition JLine1Z.hh:182
Data structure for fit of straight line in positive z-direction.
Definition JLine3Z.hh:40
static parameter_type pDY()
Definition JLine3Z.hh:320
static parameter_type pDX()
Definition JLine3Z.hh:319
JAxis3D & rotate_back(const JRotation3D &R)
Rotate back axis.
Definition JAxis3D.hh:240
JAxis3D & rotate(const JRotation3D &R)
Rotate axis.
Definition JAxis3D.hh:225
JPosition3D & rotate(const JRotation3D &R)
Rotate.
Auxiliary class for a hit with background rate value.
Definition JHitW0.hh:23
Router for fast addressing of summary data in KM3NETDAQ::JDAQSummaryslice data structure as a functio...
double getRate() const
Get default rate.
Range of values.
Definition JRange.hh:42
Template L0 hit builder.
Definition JBuildL0.hh:38
const JDAQEventHeader & getDAQEventHeader() const
Get DAQ event header.
static const int JMUONGANDALF
static const int JGANDALF_LAMBDA
control parameter from JGandalf.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 JGANDALF_BETA0_RAD
KM3NeT Data Definitions v3.5.1-1-gd514d72 https://git.km3net.de/common/km3net-dataformat.
static const int JGANDALF_NUMBER_OF_ITERATIONS
number of iterations from JGandalf.cc
static const int JPP_COVERAGE_ORIENTATION
coverage of dynamic orientation calibration from any Jpp application
static const int JGANDALF_BETA1_RAD
angular resolution [rad] from JGandalf.cc
static const int JGANDALF_NUMBER_OF_HITS
number of hits from JGandalf.cc
static const int JGANDALF_CHI2
chi2 from JGandalf.cc
void copy(const Head &from, JHead &to)
Copy header from from to to.
Definition JHead.cc:163
Auxiliary classes and methods for linear and iterative data regression.
Definition JEnergy.hh:15
JTOOLS::JRange< double > JZRange
Auxiliary classes and methods for 3D geometrical objects and operations.
Definition JAngle3D.hh:19
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Model fits to data.
double getQuality(const double chi2, const int N, const int NDF)
Get quality of fit.
JPosition3D getPosition(const JFit &fit)
Get position.
JFIT::JHistory JHistory
Definition JHistory.hh:405
bool qualitySorter(const JFit &first, const JFit &second)
Comparison of fit results.
JFit getFit(const JHistory &history, const JTrack3D &track, const double Q, const int NDF, const double energy=0.0, const int status=SINGLE_STAGE)
Get fit.
JDirection3D getDirection(const JFit &fit)
Get direction.
Auxiliary classes and methods for triggering.
Model for fit to acoustics data.
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
double TMin_ns
minimal time w.r.t. Cherenkov hypothesis [ns]
double TMax_ns
maximal time w.r.t. Cherenkov hypothesis [ns]
double VMax_npe
maximum number of of photo-electrons
input_type(const JDAQEventHeader &header, const JEvt &in, const coverage_type &coverage)
Constructor.
Wrapper class to make final fit of muon trajectory.
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< hit_type > buffer_type
JMuonGandalf(const JMuonGandalfParameters_t &parameters, const storage_type &storage, const int debug=0)
Constructor.
JRegressor< JLine3Z, JGandalf > JRegressor_t
JEvt operator()(const input_type &input)
Fit function.
Auxiliary data structure for sorting of hits.
Definition JHitL0.hh:85