Jpp 20.0.0-rc.2
the software that should make you happy
Loading...
Searching...
No Matches
JShowerDirectionPrefit.hh
Go to the documentation of this file.
1#ifndef JSHOWERDIRECTIONPREFIT_INCLUDE
2#define JSHOWERDIRECTIONPREFIT_INCLUDE
3
4#include <string>
5#include <iostream>
6#include <set>
7#include <vector>
8#include <algorithm>
9#include <memory>
10
13
14#include "JTrigger/JHit.hh"
16#include "JTrigger/JHitL0.hh"
17#include "JTrigger/JHitL1.hh"
18#include "JTrigger/JHitR0.hh"
19#include "JTrigger/JHitR1.hh"
20#include "JTrigger/JBuildL0.hh"
21#include "JTrigger/JBuildL2.hh"
23#include "JTrigger/JMatch3G.hh"
24#include "JTrigger/JBind2nd.hh"
25
27
29#include "JFit/JFitToolkit.hh"
30#include "JFit/JPoint4D.hh"
31#include "JFit/JModel.hh"
32
34
39
42
44
45#include "JPhysics/KM3NeT.hh"
46
53
55
56/**
57 * \author adomi, vcarretero
58 */
59namespace JRECONSTRUCTION {}
60namespace JPP { using namespace JRECONSTRUCTION; }
61
62namespace JRECONSTRUCTION {
69 using JFIT::JRegressor;
70 using JFIT::JEnergy;
71 using JFIT::JShower3EZ;
73
74 /**
75 * class to handle the direction fit of the shower reconstruction, mainly dedicated for ORCA
76 */
79 public JRegressor<JShower3EZ, JAbstractMinimiser>
80 {
81 public:
85 using JRegressor_t::operator();
86
88
89 /**
90 * Input data type.
91 */
92 struct input_type :
93 public JDAQEventHeader
94 {
95 /**
96 * Default constructor.
97 */
99 {}
100
101
102 /**
103 * Constructor.
104 *
105 * \param header header
106 * \param in start values
107 * \param coverage coverage
108 */
110 const JEvt& in,
111 const coverage_type& coverage) :
112 JDAQEventHeader(header),
113 in(in)
114 {}
115
119 };
120
121 public:
122
123 /**
124 * Parameterized constructor
125 *
126 * \param parameters struct that holds default-optimized parameters for the reconstruction, available in $JPP_DATA.
127 * \param storage storage
128 * \param debug debug
129 */
131 const storage_type& storage,
132 const int debug = 0):
134 JRegressor_t(storage)
135 {
136 using namespace JPP;
137
138 JRegressor_t::debug = debug;
139 JRegressor_t::T_ns.setRange(parameters.TMin_ns, parameters.TMax_ns);
140 JRegressor_t::Vmax_npe = VMax_npe;
141
142 if (Emin_GeV > Emax_GeV || En <= 1) {
143 THROW(JException, "Invalid energy input " << Emin_GeV << ' ' << Emax_GeV << ' ' << En);
144 }
145
146 const double base = std::pow((Emax_GeV / Emin_GeV), 1.0 / (En - 1));
147
148 for (int i = 0; i != En; ++i) {
149 Ev.push_back(Emin_GeV * std::pow(base, i));
150 }
151 }
152
153 /**
154 * Get input data.
155 *
156 * \param router module router
157 * \param summary summary data
158 * \param event event
159 * \param in start values
160 * \param coverage coverage
161 * \return input data
162 */
164 const JSummaryRouter& summary,
165 const JDAQEvent& event,
166 const JEvt& in,
167 const coverage_type& coverage) const
168 {
169 using namespace std;
170 using namespace JTRIGGER;
171
172 input_type input(event.getDAQEventHeader(), in, coverage);
173
174 const JBuildL0 <JHitR0> buildL0;
176
177 const JDAQTimeslice timeslice(event, true);
178
179 JSuperFrame2D<JHit> buffer;
180
181 for (JDAQTimeslice::const_iterator i = timeslice.begin(); i != timeslice.end(); ++i) {
182
183 if (router.hasModule(i->getModuleID())) {
184
185 buffer(*i, router.getModule(i->getModuleID()));
186
187 buildL0(buffer, back_inserter(data[i->getModuleID()]));
188 }
189 }
190
191 for (const auto& module : router.getReference()) {
192 if (!module.empty()) {
193 input.data.push_back(module_type(module, summary.getSummaryFrame(module.getID(), R_Hz), data[module.getID()]));
194 }
195 }
196
197 return input;
198 }
199
200 /**
201 * Fit function.
202 *
203 * \param input input data
204 * \return fit results
205 */
207 {
208 using namespace std;
209 using namespace JPP;
210
212
213 JEvt out;
214
215 JEvt in = input.in;
216
218
219 if (!in.empty()) {
220 in.select(JHistory::is_event(in.begin()->getHistory()));
221 }
222
223 for (JEvt::const_iterator shower = in.begin(); shower != in.end(); ++shower) {
224
225 vector<JPMTW0> data;
226
227 const JPosition3D vertex(getPosition(*shower));
228 const double time = shower->getT();
229
230 for (const auto& module : input.data) {
231
232 JPosition3D pos(module->getPosition());
233 pos.sub(vertex);
234
235 if(pos.getLength() <= DMax_m){
236
237 const double t1 = time + pos.getLength() * getInverseSpeedOfLight() * getIndexOfRefraction();
238
239 for (size_t i = 0; i != module->size(); ++i) {
240 if (module.getStatus(i)) {
241 struct {
242
243 bool operator()(const JHitR0& hit) const
244 {
245 return (hit.getPMT() == pmt && T_ns(hit.getT()));
246 }
247
248 const JTimeRange T_ns;
249 const size_t pmt;
250
251 } match = { JRegressor_t::T_ns + t1, i };
252
253 JPMT pmt = module->getPMT(i);
254
255 pmt.sub(vertex);
256
257 data.push_back(JPMTW0(pmt, module.frame.getRate(i), count_if(module.begin(), module.end(), match)));
258
259 }
260 }
261 }
262 }
263
264 JVector3D start_dir(0,0,0);
265 for (vector<JPMTW0>::const_iterator i = data.begin(); i != data.end(); ++i) {
266 if (i->getN() != 0) start_dir.add(i->getPosition());
267 }
268
269 JOmega3D scan_directions(JDirection3D(start_dir),
271 scanAngle_deg * JMATH::PI / 180.0);
272
273 const JShower3EZ sh(JVertex3D(JVector3D(0,0,0), shower->getT()), JVersor3Z(), 1);
274
275 for (JOmega3D_t::const_iterator dir = scan_directions.begin(); dir != scan_directions.end(); ++dir) {
276
277 const JRotation3D R(*dir);
278
279 vector<double> chi2(En, 0.0);
280
281 for (vector<JPMTW0>::const_iterator i = data.begin(); i != data.end(); ++i) {
282
283 JPMTW0 pmt = *i;
284 pmt.rotate(R);
285
286 JNPE_t::result_type H1 = (*this).getH1(sh, pmt);
287 JNPE_t::result_type H0 = (*this).getH0(pmt.getR()); // getH0 = Get background hypothesis value for time integrated PDF.
288 const bool hit = pmt.getN() != 0;
289
290 const double chi2_bg = getChi2(get_value(H0), hit);
291
292 for (size_t j=0; j!=chi2.size(); ++j) {
293
294 chi2[j]+= getChi2(get_value(Ev[j] * H1 + H0), hit) - chi2_bg; // H1 is linear with E, -log-lik ratio
295
296 }
297
298 }
299
300 //Store only the lowest chi2 for a given direction
301
302 auto p = std::min_element(chi2.begin(), chi2.end());
303
304 out.push_back(getFit(JHistory(shower->getHistory(), event()),
305 JShower3D(JVertex3D(vertex,time), JDirection3D(*dir)),
306 getQuality(*p),
307 data.size(),
308 Ev[p-chi2.begin()]));
309
310 // set additional values
311 out.rbegin()->setW(JPP_COVERAGE_ORIENTATION, input.coverage.orientation);
312 out.rbegin()->setW(JPP_COVERAGE_POSITION, input.coverage.position);
313
314 }
315 }
316
317 // apply default sorter
318
319 sort(out.begin(), out.end(), qualitySorter);
320
321 copy(input.in.begin(), input.in.end(), back_inserter(out));
322
323 return out;
324
325 }
326
327 };
328}
329
330#endif
Definition of hit and track types and auxiliary methods for handling Monte Carlo data.
Algorithms for hit clustering and sorting.
Coverage of dynamical detector calibration.
Auxiliary class to extract a subset of optical modules from a detector.
Data structure for detector geometry and calibration.
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Auxiliary methods to evaluate Poisson probabilities and chi2.
Basic data structure for L0 hit.
Basic data structure for L1 hit.
Basic data structure for L0 hit.
Reduced data structure for L1 hit.
Match operator for Cherenkov light from shower in any direction.
int debug
debug level
Definition JSirene.cc:72
Data regression method for JFIT::JShower3EZ.
Basic data structure for time and time over threshold information of hit.
Properties of KM3NeT PMT and deep-sea water.
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
Abstract minimiser.
Definition JRegressor.hh:27
Data structure for fit of energy.
Data structure for set of track fit results.
void select(const JSelector_t &selector)
Select fits.
Data structure for fit of straight line in positive z-direction with energy.
Definition JShower3EZ.hh:30
JAxis3D & rotate(const JRotation3D &R)
Rotate axis.
Definition JAxis3D.hh:225
Data structure for direction in three dimensions.
Direction set covering (part of) solid angle.
Definition JOmega3D.hh:68
Data structure for position in three dimensions.
Data structure for vector in three dimensions.
Definition JVector3D.hh:36
JVector3D & add(const JVector3D &vector)
Add vector.
Definition JVector3D.hh:142
double getLength() const
Get length.
Definition JVector3D.hh:246
JVector3D & sub(const JVector3D &vector)
Subtract vector.
Definition JVector3D.hh:158
Data structure for normalised vector in positive z-direction.
Definition JVersor3Z.hh:41
General exception.
Definition JException.hh:24
const JClass_t & getReference() const
Get reference to object.
Definition JReference.hh:38
class to handle the direction fit of the shower reconstruction, mainly dedicated for ORCA
input_type getInput(const JModuleRouter &router, const JSummaryRouter &summary, const JDAQEvent &event, const JEvt &in, const coverage_type &coverage) const
Get input data.
JEvt operator()(const input_type &input)
Fit function.
JShowerDirectionPrefit(const JShowerDirectionPrefitParameters_t &parameters, const storage_type &storage, const int debug=0)
Parameterized constructor.
JRegressor< JShower3EZ, JAbstractMinimiser > JRegressor_t
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 JPP_COVERAGE_POSITION
coverage of dynamic position calibration from any Jpp application
static const int JPP_COVERAGE_ORIENTATION
coverage of dynamic orientation calibration from any Jpp application
double getChi2(const double P)
Get chi2 corresponding to given probability.
static const double PI
Mathematical constants.
double getIndexOfRefraction()
Get average index of refraction of water corresponding to group velocity.
const double getInverseSpeedOfLight()
Get inverse speed of light.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
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:455
void copy(const JFIT::JEvt::const_iterator __begin, const JFIT::JEvt::const_iterator __end, Evt &out)
Copy tracks.
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.
if((p==this->begin() &&this->getDistance(x,(p++) ->getX()) > distance_type::precision)||(p==this->end() &&this->getDistance((--p) ->getX(), x) > distance_type::precision))
Template base class for polynomial interpolations with polynomial result.
Definition JPolint.hh:775
int j
Definition JPolint.hh:801
JResultEvaluator< JResult_t >::result_type get_value(const JResult_t &value)
Helper method to recursively evaluate a to function value.
Definition JResult.hh:998
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:40
Auxiliary class to test history.
Definition JHistory.hh:157
Auxiliary class for handling PMT geometry, rate and response.
Definition JPMTW0.hh:24
int getN() const
Get number of hits.
Definition JPMTW0.hh:67
double getR() const
Get rate.
Definition JPMTW0.hh:56
Template definition of a data regressor of given model.
Definition JRegressor.hh:70
Auxiliary class for handling module response.
Definition JModuleL0.hh:45
input_type(const JDAQEventHeader &header, const JEvt &in, const coverage_type &coverage)
Constructor.