Jpp test-rotations-old
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();
87
88 /**
89 * Input data type.
90 */
91 struct input_type :
92 public JDAQEventHeader
93 {
94 /**
95 * Default constructor.
96 */
98 {}
99
100
101 /**
102 * Constructor.
103 *
104 * \param header header
105 * \param in start values
106 * \param coverage coverage
107 */
109 const JEvt& in,
110 const coverage_type& coverage) :
111 JDAQEventHeader(header),
112 in(in)
113 {}
114
118 };
119
120 public:
121
122 /**
123 * Parameterized constructor
124 *
125 * \param parameters struct that holds default-optimized parameters for the reconstruction, available in $JPP_DATA.
126 * \param storage storage
127 * \param debug debug
128 */
130 const storage_type& storage,
131 const int debug = 0):
133 JRegressor_t(storage)
134 {
135 using namespace JPP;
136
137 JRegressor_t::debug = debug;
138 JRegressor_t::T_ns.setRange(parameters.TMin_ns, parameters.TMax_ns);
139 JRegressor_t::Vmax_npe = VMax_npe;
140
141 if (Emin_GeV > Emax_GeV || En <= 1) {
142 THROW(JException, "Invalid energy input " << Emin_GeV << ' ' << Emax_GeV << ' ' << En);
143 }
144
145 const double base = std::pow((Emax_GeV / Emin_GeV), 1.0 / (En - 1));
146
147 for (int i = 0; i != En; ++i) {
148 Ev.push_back(Emin_GeV * std::pow(base, i));
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 input_type input(event.getDAQEventHeader(), in, coverage);
172
173 const JBuildL0 <JHitR0> buildL0;
175
176 const JDAQTimeslice timeslice(event, true);
177
178 JSuperFrame2D<JHit> buffer;
179
180 for (JDAQTimeslice::const_iterator i = timeslice.begin(); i != timeslice.end(); ++i) {
181
182 if (router.hasModule(i->getModuleID())) {
183
184 buffer(*i, router.getModule(i->getModuleID()));
185
186 buildL0(buffer, back_inserter(data[i->getModuleID()]));
187 }
188 }
189
190 for (const auto& module : router.getReference()) {
191 if (!module.empty()) {
192 input.data.push_back(module_type(module, summary.getSummaryFrame(module.getID(), R_Hz), data[module.getID()]));
193 }
194 }
195
196 return input;
197 }
198
199 /**
200 * Fit function.
201 *
202 * \param input input data
203 * \return fit results
204 */
206 {
207 using namespace std;
208 using namespace JPP;
209
211
212 JEvt out;
213
214 JEvt in = input.in;
215
217
218 if (!in.empty()) {
219 in.select(JHistory::is_event(in.begin()->getHistory()));
220 }
221
222 for (JEvt::const_iterator shower = in.begin(); shower != in.end(); ++shower) {
223
224 vector<JPMTW0> data;
225
226 const JPosition3D vertex(getPosition(*shower));
227 const double time = shower->getT();
228
229 for (const auto& module : input.data) {
230
231 JPosition3D pos(module->getPosition());
232 pos.sub(vertex);
233
234 if(pos.getLength() <= DMax_m){
235
236 const double t1 = time + pos.getLength() * getInverseSpeedOfLight() * getIndexOfRefraction();
237
238 for (size_t i = 0; i != module->size(); ++i) {
239 if (module.getStatus(i)) {
240 struct {
241
242 bool operator()(const JHitR0& hit) const
243 {
244 return (hit.getPMT() == pmt && T_ns(hit.getT()));
245 }
246
247 const JTimeRange T_ns;
248 const size_t pmt;
249
250 } match = { JRegressor_t::T_ns + t1, i };
251
252 JPMT pmt = module->getPMT(i);
253
254 pmt.sub(vertex);
255
256 data.push_back(JPMTW0(pmt, module.frame.getRate(i), count_if(module.begin(), module.end(), match)));
257
258 }
259 }
260 }
261 }
262
263 JVector3D start_dir(0,0,0);
264 for (vector<JPMTW0>::const_iterator i = data.begin(); i != data.end(); ++i) {
265 if (i->getN() != 0) start_dir.add(i->getPosition());
266 }
267
268 JOmega3D scan_directions(JDirection3D(start_dir),
270 scanAngle_deg * JMATH::PI / 180.0);
271
272 const JShower3EZ sh(JVertex3D(JVector3D(0,0,0), shower->getT()), JVersor3Z(), 1);
273
274 for (JOmega3D_t::const_iterator dir = scan_directions.begin(); dir != scan_directions.end(); ++dir) {
275
276 const JRotation3D R(*dir);
277
278 vector<double> chi2(En, 0.0);
279
280 for (vector<JPMTW0>::const_iterator i = data.begin(); i != data.end(); ++i) {
281
282 JPMTW0 pmt = *i;
283 pmt.rotate(R);
284
285 JNPE_t::result_type H1 = (*this).getH1(sh, pmt);
286 JNPE_t::result_type H0 = (*this).getH0(pmt.getR()); // getH0 = Get background hypothesis value for time integrated PDF.
287 const bool hit = pmt.getN() != 0;
288
289 const double chi2_bg = getChi2(get_value(H0), hit);
290
291 for (size_t j=0; j!=chi2.size(); ++j) {
292
293 chi2[j]+= getChi2(get_value(Ev[j] * H1 + H0), hit) - chi2_bg; // H1 is linear with E, -log-lik ratio
294
295 }
296
297 }
298
299 //Store only the lowest chi2 for a given direction
300
301 auto p = std::min_element(chi2.begin(), chi2.end());
302
303 out.push_back(getFit(JHistory(shower->getHistory(), event()),
304 JShower3E(JVertex3D(vertex,time), JDirection3D(*dir), Ev[p-chi2.begin()]),
305 getQuality(*p),
306 data.size(),
307 Ev[p-chi2.begin()]));
308
309 // set additional values
310 out.rbegin()->setW(JPP_COVERAGE_ORIENTATION, input.coverage.orientation);
311 out.rbegin()->setW(JPP_COVERAGE_POSITION, input.coverage.position);
312
313 }
314 }
315
316 // apply default sorter
317
318 sort(out.begin(), out.end(), qualitySorter);
319
320 copy(input.in.begin(), input.in.end(), back_inserter(out));
321
322 return out;
323
324 }
325
326 };
327}
328
329#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.
Definition JEnergy.hh:31
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 JSHOWERDIRECTIONPREFIT
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
void copy(const Head &from, JHead &to)
Copy header from from to to.
Definition JHead.cc:163
double getChi2(const double P)
Get chi2 corresponding to given probability.
JTrack3E JShower3E
Type definition of 3D shower with energy.
Definition JShower3E.hh:19
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).
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.
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:36
Auxiliary class to test history.
Definition JHistory.hh:136
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.