118{
122
124 typedef JTriggeredFileScanner_t::multi_pointer_type multi_pointer_type;
125
127
128 JTriggeredFileScanner_t inputFile;
130 string detectorFile;
136
137 try {
138
140
141 JParser<> zap(
"Program to test JMuonStart.");
142
143 zap[
'f'] =
make_field(inputFile,
"input file (output of JXXXMuonReconstruction.sh)");
151
152 zap(argc, argv);
153 }
154 catch(const exception& error) {
155 FATAL(error.what() << endl);
156 }
157
158
160
161 try {
163 }
166 }
167
170
172
175
177
180
181
183
184 try {
186 }
189 }
190
192
198
199 while (inputFile.hasNext()) {
200
201 STATUS(
"event: " << setw(10) << inputFile.getCounter() <<
'\r');
DEBUG(endl);
202
203 multi_pointer_type ps = inputFile.next();
204
208
209 summary.update(*tev);
210
212
214
216 if (muon == event->mc_trks.end() || track->E > muon->E) {
217 muon = track;
218 }
219 }
220 }
221
222 if (muon == event->mc_trks.end()) {
223 continue;
224 }
225
226 if (muon->len == 0.0) {
227 muon->len =
gWater(muon->E);
228 }
229
230
232
233 JRange_t za(JRange_t::DEFAULT_RANGE());
234 JRange_t zb(JRange_t::DEFAULT_RANGE());
235 JRange_t zd(JRange_t::DEFAULT_RANGE());
236
237 {
239
241
244
245 za.setLowerLimit(ta.
getZ());
246 za.setLength(fabs(muon->len));
247 }
248
249
251
252 buildL0(*tev, moduleRouter, true, back_inserter(dataL0));
253
254
255 if (!evt->empty()) {
256
258
259 JEvt::const_iterator track = evt->begin();
260
264
268 }
269 }
270
271 {
273
274 for (buffer_type::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {
275
277
278 hit.rotate(R);
279
280 if (match(hit)) {
281 top[hit.getModuleID()].insert(hit.getPMTAddress());
282 }
283 }
284
286
287 double getZ() const { return z; }
288 double getP()
const {
return p; }
289
290 double z;
291 double p;
292 };
293
295
296 for (JDetector::const_iterator module =
detector.begin(); module !=
detector.end(); ++module) {
297
298 if (summary.hasSummaryFrame(module->getID())) {
299
301
303
304 pos.transform(R, tz.getPosition());
305
307
308 const double z = pos.getZ() - pos.getX() /
getTanThetaC();
310
311 data.push_back({ z, p });
312 }
313 }
314 }
315
317
319
320 vector<JHit_t>::const_iterator track_start = start.find(
data. begin(),
data. end());
321 vector<JHit_t>::const_reverse_iterator track_end = start.find(
data.rbegin(),
data.rend());
322
323 if (track_start ==
data. end()) { track_start =
data. begin(); }
324 if (track_end ==
data.rend()) { track_end =
data.rbegin(); }
325
326 zd =
JRange_t(track_start->getZ(), track_end->getZ());
327
328 zd.add(tz.getZ());
329 }
330 }
331 }
332
333 H1["B"]->Fill(zb.getLength() - za.getLength());
334 H1["D"]->Fill(zd.getLength() - za.getLength());
335
336 H2["B"]->Fill(zb.getLowerLimit() - za.getLowerLimit(), zb.getUpperLimit() - za.getUpperLimit());
337 H2["D"]->Fill(zd.getLowerLimit() - za.getLowerLimit(), zd.getUpperLimit() - za.getUpperLimit());
338 }
340
341
343
344 out << H1 << H2;
345
346 out.Write();
347 out.Close();
348}
#define DEBUG(A)
Message macros.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Router for direct addressing of module data in detector data structure.
Router for direct addressing of PMT data in detector data structure.
Data structure for fit of straight line paralel to z-axis.
JAxis3D & rotate(const JRotation3D &R)
Rotate axis.
Data structure for position in three dimensions.
JPosition3D & rotate(const JRotation3D &R)
Rotate.
JTime & add(const JTime &value)
Addition operator.
double getZ() const
Get z position.
Utility class to parse command line options.
Auxiliary class to manage set of compatible ROOT objects (e.g. histograms) using unique keys.
File router for fast addressing of summary data.
Data structure for L0 hit.
Data storage class for rate measurements of all PMTs in one module.
static const int JMUONSTART
static const int JSTART_LENGTH_METRES
distance between projected positions on the track of optical modules for which the response does not ...
JDirection3D getDirection(const Vec &dir)
Get direction.
JTrack3E getTrack(const Trk &track)
Get track.
JPosition3D getPosition(const Vec &pos)
Get position.
bool is_muon(const Trk &track)
Test whether given track is a (anti-)muon.
Vec getOffset(const JHead &header)
Get offset.
JAbstractHistogram< double > JHistogram_t
Type definition for scan along axis.
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
double getProbability(const size_t N, const size_t M, const JK40Rates &R_Hz, const double T_ns)
Get probability due to random background.
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.
static const JGeaneWater gWater
Function object for energy loss of muon in sea water.
double getTanThetaC()
Get average tangent of Cherenkov angle of water corresponding to group velocity.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
std::vector< JHitW0 > buffer_type
hits
bool qualitySorter(const JFit &first, const JFit &second)
Comparison of fit results.
Head getHeader(const JMultipleFileScanner_t &file_list)
Get Monte Carlo header.
KM3NeT DAQ data structures and auxiliaries.
const JK40Rates & getK40Rates()
Get K40 rates.
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
Type definition of range.
Auxiliary class to match data points with given model.
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Auxiliary class for K40 rates.
const JRateL1_t & getMultiplesRates() const
Get multiples rate.
Data structure for fit parameters.
double roadWidth_m
road width [m]
size_t numberOfPrefits
number of prefits
int Nmax2
maximal number for twofold observations
double Pmin1
minimal probability single observation
double TMin_ns
minimal time w.r.t. Cherenkov hypothesis [ns]
double Pmin2
minimal probability for twofold observations
double TMax_ns
maximal time w.r.t. Cherenkov hypothesis [ns]
Auxiliary class for start or end point evaluation.
Auxiliary class to set-up Hit.
Auxiliary class for defining the range of iterations of objects.
const JLimit & getLimit() const
Get limit.
static counter_type max()
Get maximum counter value.
Auxiliary class to synchronously read DAQ events and Monte Carlo events (and optionally other events)...
The Vec class is a straightforward 3-d vector, which also works in pyroot.