Jpp 19.3.0-rc.1
the software that should make you happy
Loading...
Searching...
No Matches
JMuonPath.cc
Go to the documentation of this file.
1
2#include <string>
3#include <iostream>
4#include <iomanip>
5#include <vector>
6#include <algorithm>
7
12
13#include "JDAQ/JDAQEventIO.hh"
16
20
21#include "JTrigger/JHit.hh"
22#include "JTrigger/JFrame.hh"
24#include "JTrigger/JHitL0.hh"
25#include "JTrigger/JBuildL0.hh"
27
31#include "JSupport/JSupport.hh"
32#include "JSupport/JMeta.hh"
33
34#include "JFit/JPMTW0.hh"
35#include "JFit/JLine3Z.hh"
36#include "JFit/JGandalf.hh"
38#include "JFit/JFitToolkit.hh"
39#include "JFit/JModel.hh"
40
44
48
49#include "JPhysics/JPDFTable.hh"
51
52#include "Jeep/JParser.hh"
53#include "Jeep/JMessage.hh"
54
55
56/**
57 * \file
58 *
59 * Program to perform JFIT::JRegressor<JLine3Z,JGandalf> fit with I/O of JFIT::JEvt data.
60 * \author mdejong
61 */
62int main(int argc, char **argv)
63{
64 using namespace std;
65 using namespace JPP;
66 using namespace KM3NETDAQ;
67
68 typedef JParallelFileScanner< JTypeList<JDAQEvent, JEvt> > JParallelFileScanner_t;
69 typedef JParallelFileScanner_t::multi_pointer_type multi_pointer_type;
71
72 JParallelFileScanner_t inputFile;
74 JLimit_t& numberOfEvents = inputFile.getLimit();
75 string detectorFile;
76 string pdfFile;
77 double roadWidth_m;
78 double R_Hz;
79 size_t numberOfPrefits;
80 double TTS_ns;
81 double E_GeV;
82 int debug;
83
84 try {
85
86 JParser<> zap("Program to perform fit of muon path to data.");
87
88 zap['f'] = make_field(inputFile);
89 zap['o'] = make_field(outputFile) = "path.root";
90 zap['a'] = make_field(detectorFile);
91 zap['n'] = make_field(numberOfEvents) = JLimit::max();
92 zap['P'] = make_field(pdfFile);
93 zap['R'] = make_field(roadWidth_m) = numeric_limits<double>::max();
94 zap['B'] = make_field(R_Hz) = 6.0e3;
95 zap['N'] = make_field(numberOfPrefits) = 1;
96 zap['T'] = make_field(TTS_ns);
97 zap['E'] = make_field(E_GeV) = 1.0e3;
98 zap['d'] = make_field(debug) = 1;
99
100 zap(argc, argv);
101 }
102 catch(const exception& error) {
103 FATAL(error.what() << endl);
104 }
105
106
108
109 try {
110 load(detectorFile, detector);
111 }
112 catch(const JException& error) {
113 FATAL(error);
114 }
115
116 const JModuleRouter router(detector);
117
118
119 typedef vector<JHitL0> JDataL0_t;
120 const JBuildL0<JHitL0> buildL0;
121
122
124
125 JRegressor_t::debug = debug;
126 JRegressor_t::T_ns.setRange(-50.0, +50.0); // ns
127 JRegressor_t::Vmax_npe = 10.0; // npe
128 JRegressor_t::MAXIMUM_ITERATIONS = 10000;
129
130 JRegressor_t fit(pdfFile, TTS_ns);
131
132 fit.transform(JRegressor_t::transformer_type::getDefaultTransformer()); // get rid of R-dependent weight function
133
134 fit.estimator.reset(new JMEstimatorLorentzian());
135
136 fit.parameters.resize(4);
137
138 fit.parameters[0] = JLine3Z::pX();
139 fit.parameters[1] = JLine3Z::pY();
140 fit.parameters[2] = JLine3Z::pDX();
141 fit.parameters[3] = JLine3Z::pDY();
142
143 if (fit.getRmax() < roadWidth_m) {
144
145 roadWidth_m = fit.getRmax();
146
147 WARNING("Set road width to [m] " << roadWidth_m << endl);
148 }
149
150 const double Rmax_m = 100.0; // hit selection [m]
151 const double Tmax_ns = 10.0; // hit selection [ns]
152
153
154 outputFile.open();
155
156 outputFile.put(JMeta(argc, argv));
157
158 while (inputFile.hasNext()) {
159
160 STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);
161
162 multi_pointer_type ps = inputFile.next();
163
164 const JDAQEvent* tev = ps;
165 const JEvt* in = ps;
166
167 JEvt cp = *in;
168 JEvt out;
169
170 cp.select(numberOfPrefits, qualitySorter);
171
172
173 JDataL0_t dataL0;
174
175 buildL0(*tev, router, true, back_inserter(dataL0));
176
177
178 for (JEvt::const_iterator track = cp.begin(); track != cp.end(); ++track) {
179
180 const JRotation3D R (getDirection(*track));
181 const JLine1Z tz(getPosition (*track).rotate(R), track->getT());
182 const JModel<JLine1Z> match(tz, roadWidth_m, JRegressor_t::T_ns);
183
184 // hit selection based on start value
185
187
189
190 for (JDataL0_t::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {
191
192 JHitW0 hit(*i, R_Hz);
193
194 hit.rotate(R);
195
196 if (match(hit)) {
197
198 top.insert(hit.getPMTIdentifier());
199
200 const double t1 = hit.getT() - tz.getT(hit);
201
202 if (tz.getDistance(hit) <= Rmax_m && t1 >= -Tmax_ns && t1 <= +Tmax_ns) {
203 Z.include(hit.getZ() - tz.getDistance(hit) / getTanThetaC());
204 }
205 }
206 }
207
208
209 vector<JPMTW0> buffer;
210
211 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
212
213 JPosition3D pos(module->getPosition());
214
215 pos.rotate(R);
216
217 if (tz.getDistance(pos) <= roadWidth_m && Z(pos.getZ() - tz.getDistance(pos) / getTanThetaC())) {
218
219 for (unsigned int i = 0; i != module->size(); ++i) {
220
221 const JDAQPMTIdentifier id(module->getID(), i);
222
223 JPMT pmt(module->getPMT(i));
224
225 pmt.rotate(R);
226
227 buffer.push_back(JPMTW0(pmt, R_Hz, top.count(id)));
228 }
229 }
230 }
231
232
233 const int NDF = buffer.size() - fit.parameters.size();
234
235 if (NDF >= 0) {
236
237 // set fit parameters
238
239 if (track->getE() > 0.1)
240 fit.E_GeV = track->getE();
241 else
242 fit.E_GeV = E_GeV;
243
244 const double chi2 = fit(JLine3Z(tz), buffer.begin(), buffer.end());
245
246 JTrack3D tb(fit.value);
247
248 tb.rotate_back(R);
249
250 out.push_back(getFit(JMUONPATH, tb, getQuality(chi2), NDF));
251
252 out.rbegin()->setW(track->getW());
253 }
254 }
255
256 // apply default sorter
257
258 sort(out.begin(), out.end(), qualitySorter);
259
260 outputFile.put(out);
261 }
262 STATUS(endl);
263
265
266 io >> outputFile;
267
268 outputFile.close();
269}
string outputFile
Data structure for detector geometry and calibration.
Recording of objects on file according a format that follows from the file name extension.
Auxiliary methods to evaluate Poisson probabilities and chi2.
Various implementations of functional maps.
Basic data structure for L0 hit.
Data regression method for JFIT::JLine3Z.
General purpose messaging.
#define DEBUG(A)
Message macros.
Definition JMessage.hh:62
#define STATUS(A)
Definition JMessage.hh:63
#define FATAL(A)
Definition JMessage.hh:67
int debug
debug level
Definition JSirene.cc:72
#define WARNING(A)
Definition JMessage.hh:65
ROOT I/O of application specific meta data.
Direct access to module in detector data structure.
int main(int argc, char **argv)
Definition JMuonPath.cc:62
Auxiliary methods for PDF calculations.
Parallel scanning of objects from a single file or multiple files according a format that follows fro...
Utility class to parse command line options.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2142
Physics constants.
Scanning of objects from a single file according a format that follows from the extension of each fil...
ROOT TTree parameter settings of various packages.
Basic data structure for time and time over threshold information of hit.
Detector data structure.
Definition JDetector.hh:96
Router for direct addressing of module data in detector data structure.
Data structure for PMT geometry, calibration and status.
Definition JPMT.hh:49
Data structure for fit of straight line paralel to z-axis.
Definition JLine1Z.hh:29
static parameter_type pY()
Definition JLine1Z.hh:181
double getT(const JVector3D &pos) const
Get arrival time of Cherenkov light at given position.
Definition JLine1Z.hh:114
static parameter_type pX()
Definition JLine1Z.hh:180
double getDistance(const JVector3D &pos) const
Get distance.
Definition JLine1Z.hh:102
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
Data structure for position in three dimensions.
JPosition3D & rotate(const JRotation3D &R)
Rotate.
double getZ() const
Get z position.
Definition JVector3D.hh:115
General exception.
Definition JException.hh:24
Utility class to parse command line options.
Definition JParser.hh:1698
Auxiliary class for a hit with background rate value.
Definition JHitW0.hh:23
Object writing to file.
General purpose class for parallel reading of objects from a single file or multiple files.
Object reading from a list of files.
Range of values.
Definition JRange.hh:42
static JRange< T, JComparator_t > DEFAULT_RANGE()
Default range.
Definition JRange.hh:555
range_type & include(argument_type x)
Include given value to range.
Definition JRange.hh:397
Template L0 hit builder.
Definition JBuildL0.hh:38
double getT() const
Get calibrated time of hit.
const JDAQPMTIdentifier & getPMTIdentifier() const
Get PMT identifier.
static const int JMUONPATH
JDirection3D getDirection(const Vec &dir)
Get direction.
JPosition3D getPosition(const Vec &pos)
Get position.
double getQuality(const JEvent &evt)
Get average quality.
JFit getFit(const int id, const JMODEL::JString &string)
Get fit parameters of string.
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
bool qualitySorter(const JFit &first, const JFit &second)
Comparison of fit results.
JFIT::JRegressor< JFIT::JLine3Z, JFIT::JGandalf > JRegressor_t
Definition JPerth.cc:132
KM3NeT DAQ data structures and auxiliaries.
Definition DataQueue.cc:39
Detector file.
Definition JHead.hh:227
Acoustic event fit.
Model for fit to acoustics data.
Lorentzian M-estimator.
Auxiliary class for handling PMT geometry, rate and response.
Definition JPMTW0.hh:24
Template definition of a data regressor of given model.
Definition JRegressor.hh:70
Type list.
Definition JTypeList.hh:23
Auxiliary class for defining the range of iterations of objects.
Definition JLimit.hh:45
const JLimit & getLimit() const
Get limit.
Definition JLimit.hh:84
static counter_type max()
Get maximum counter value.
Definition JLimit.hh:128
Auxiliary class for ROOT I/O of application specific meta data.
Definition JMeta.hh:72