Jpp test-rotations-old
the software that should make you happy
Loading...
Searching...
No Matches
JMuonStart.cc File Reference

Program to test JMuonStart.cc. More...

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <set>
#include <map>
#include <algorithm>
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "km3net-dataformat/offline/Head.hh"
#include "km3net-dataformat/offline/Evt.hh"
#include "km3net-dataformat/online/JDAQ.hh"
#include "km3net-dataformat/online/JDAQPMTIdentifier.hh"
#include "km3net-dataformat/online/JDAQSummaryFrame.hh"
#include "JDAQ/JDAQEventIO.hh"
#include "JDAQ/JDAQSummarysliceIO.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JModuleRouter.hh"
#include "JDetector/JPMTRouter.hh"
#include "JTrigger/JHitL0.hh"
#include "JTrigger/JBuildL0.hh"
#include "JTrigger/JTriggerToolkit.hh"
#include "JSupport/JSupport.hh"
#include "JSupport/JTriggeredFileScanner.hh"
#include "JSupport/JSummaryFileRouter.hh"
#include "JSupport/JMonteCarloFileSupportkit.hh"
#include "JAAnet/JAAnetToolkit.hh"
#include "JFit/JLine1Z.hh"
#include "JFit/JNPEHit.hh"
#include "JFit/JModel.hh"
#include "JPhysics/KM3NeT.hh"
#include "JPhysics/JGeane.hh"
#include "JReconstruction/JHitW0.hh"
#include "JReconstruction/JEvt.hh"
#include "JReconstruction/JEvtToolkit.hh"
#include "JReconstruction/JMuonStartParameters_t.hh"
#include "JReconstruction/JStart.hh"
#include "JROOT/JManager.hh"
#include "JTools/JAbstractHistogram.hh"
#include "JLang/JPredicate.hh"
#include "JLang/JComparator.hh"
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

Go to the source code of this file.

Functions

int main (int argc, char **argv)
 

Detailed Description

Program to test JMuonStart.cc.

Author
mdejong

Definition in file examples/JReconstruction/JMuonStart.cc.

Function Documentation

◆ main()

int main ( int argc,
char ** argv )

Definition at line 117 of file examples/JReconstruction/JMuonStart.cc.

118{
119 using namespace std;
120 using namespace JPP;
121 using namespace KM3NETDAQ;
122
123 typedef JTriggeredFileScanner<JEvt, JSingleFileScanner> JTriggeredFileScanner_t;
124 typedef JTriggeredFileScanner_t::multi_pointer_type multi_pointer_type;
125
127
128 JTriggeredFileScanner_t inputFile;
129 JLimit_t& numberOfEvents = inputFile.getLimit();
130 string detectorFile;
131 string outputFile;
132 JMuonStartParameters_t parameters;
133 JK40Rates rates_Hz;
134 JHistogram_t histogram;
135 int debug;
136
137 try {
138
139 parameters.numberOfPrefits = 1;
140
141 JParser<> zap("Program to test JMuonStart.");
142
143 zap['f'] = make_field(inputFile, "input file (output of JXXXMuonReconstruction.sh)");
144 zap['a'] = make_field(detectorFile);
145 zap['n'] = make_field(numberOfEvents) = JLimit::max();
146 zap['o'] = make_field(outputFile) = "start.root";
147 zap['@'] = make_field(parameters) = JPARSER::initialised();
148 zap['B'] = make_field(rates_Hz) = KM3NET::getK40Rates();
149 zap['H'] = make_field(histogram, "histogram binning") = JHistogram_t(200, -250.0, +250.0);
150 zap['d'] = make_field(debug) = 1;
151
152 zap(argc, argv);
153 }
154 catch(const exception& error) {
155 FATAL(error.what() << endl);
156 }
157
158
160
161 try {
162 load(detectorFile, detector);
163 }
164 catch(const JException& error) {
165 FATAL(error);
166 }
167
168 const JModuleRouter moduleRouter(detector);
169 const JPMTRouter pmtRouter (detector);
170
171 JSummaryFileRouter summary(inputFile);
172
173 const JTimeRange T_ns (parameters.TMin_ns, parameters.TMax_ns);
174 const JStart start(parameters.Pmin1, parameters.Pmin2, parameters.Nmax2);
175
176 typedef JHitL0 hit_type;
177
179 const JBuildL0<hit_type> buildL0;
180
181
182 JHead head;
183
184 try {
185 head = getHeader(inputFile);
186 }
187 catch(const JException& error) {
188 FATAL(error);
189 }
190
191 const Vec offset = getOffset(head);
192
193 JManager<string, TH1D> H1(new TH1D("[%].1L", NULL,
194 histogram.getNumberOfBins(), histogram.getLowerLimit(), histogram.getUpperLimit()));
195 JManager<string, TH2D> H2(new TH2D("[%].2D", NULL,
196 histogram.getNumberOfBins(), histogram.getLowerLimit(), histogram.getUpperLimit(),
197 histogram.getNumberOfBins(), histogram.getLowerLimit(), histogram.getUpperLimit()));
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
205 JDAQEvent* tev = ps;
206 JEvt* evt = ps;
207 Evt* event = ps;
208
209 summary.update(*tev);
210
211 vector<Trk>::iterator muon = event->mc_trks.end();
212
213 for (vector<Trk>::iterator track = event->mc_trks.begin(); track != event->mc_trks.end(); ++track) {
214
215 if (is_muon(*track)) {
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
231 typedef JRange<double> JRange_t;
232
233 JRange_t za(JRange_t::DEFAULT_RANGE()); // Monte Carlo
234 JRange_t zb(JRange_t::DEFAULT_RANGE()); // original result
235 JRange_t zd(JRange_t::DEFAULT_RANGE()); // new
236
237 {
238 const JRotation3D R(getDirection(*muon));
239
240 JTrack3D ta = getTrack(*muon);
241
242 ta.add(getPosition(offset));
243 ta.rotate(R);
244
245 za.setLowerLimit(ta.getZ());
246 za.setLength(fabs(muon->len));
247 }
248
249
250 buffer_type dataL0;
251
252 buildL0(*tev, moduleRouter, true, back_inserter(dataL0));
253
254
255 if (!evt->empty()) {
256
257 sort(evt->begin(), evt->end(), qualitySorter);
258
259 JEvt::const_iterator track = evt->begin();
260
261 const JRotation3D R (getDirection(*track));
262 const JLine1Z tz(getPosition (*track).rotate(R), track->getT());
263 const JFIT::JModel<JLine1Z> match(tz, parameters.roadWidth_m, T_ns);
264
265 if (track->has(JMUONSTART)) {
266 if (track->hasW(JSTART_LENGTH_METRES) && track->getW(JSTART_LENGTH_METRES) > 0.0) {
267 zb = JRange_t(tz.getZ(), tz.getZ() + track->getW(JSTART_LENGTH_METRES));
268 }
269 }
270
271 {
273
274 for (buffer_type::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {
275
276 hit_type hit(*i);
277
278 hit.rotate(R);
279
280 if (match(hit)) {
281 top[hit.getModuleID()].insert(hit.getPMTAddress());
282 }
283 }
284
285 struct JHit_t {
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
300 const JDAQSummaryFrame& frame = summary.getSummaryFrame(module->getID());
301
302 JPosition3D pos(module->getPosition());
303
304 pos.transform(R, tz.getPosition());
305
306 if (pos.getX() <= parameters.roadWidth_m) {
307
308 const double z = pos.getZ() - pos.getX() / getTanThetaC();
309 const double p = getProbability(*module, frame, rates_Hz.getMultiplesRates(), T_ns.getLength(), top[module->getID()]);
310
311 data.push_back({ z, p });
312 }
313 }
314 }
315
316 if (!data.empty()) {
317
318 sort(data.begin(), data.end(), make_comparator(&JHit_t::getZ));
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 }
339 STATUS(endl);
340
341
342 TFile out(outputFile.c_str(), "recreate");
343
344 out << H1 << H2;
345
346 out.Write();
347 out.Close();
348}
string outputFile
#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 make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2142
Monte Carlo run header.
Definition JHead.hh:1236
Detector data structure.
Definition JDetector.hh:96
Router for direct addressing of module data in detector data structure.
Router for direct addressing of PMT data in detector data structure.
Definition JPMTRouter.hh:37
Data structure for fit of straight line paralel to z-axis.
Definition JLine1Z.hh:29
JAxis3D & rotate(const JRotation3D &R)
Rotate axis.
Definition JAxis3D.hh:225
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.
Definition JVector3D.hh:115
General exception.
Definition JException.hh:24
Utility class to parse command line options.
Definition JParser.hh:1698
Auxiliary class to manage set of compatible ROOT objects (e.g. histograms) using unique keys.
Definition JManager.hh:47
File router for fast addressing of summary data.
T getLowerLimit() const
Get lower limit.
Definition JRange.hh:202
T getUpperLimit() const
Get upper limit.
Definition JRange.hh:213
Template L0 hit builder.
Definition JBuildL0.hh:38
Data structure for L0 hit.
Definition JHitL0.hh:31
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.
Definition JBillabong.cc:61
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.
Definition JGeane.hh:381
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
Definition JPerth.cc:70
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.
Definition DataQueue.cc:39
const JK40Rates & getK40Rates()
Get K40 rates.
Definition KM3NeT.hh:36
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
Definition Evt.hh:21
Type definition of range.
Definition JHead.hh:43
Detector file.
Definition JHead.hh:227
Acoustic event fit.
Acoustic hit.
Definition JBillabong.cc:70
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)...
Definition JParser.hh:68
Auxiliary class for K40 rates.
Definition JK40Rates.hh:41
const JRateL1_t & getMultiplesRates() const
Get multiples rate.
Definition JK40Rates.hh:82
Data structure for fit parameters.
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.
Definition JStart.hh:21
Auxiliary class to set-up Hit.
Definition JSirene.hh:58
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 to synchronously read DAQ events and Monte Carlo events (and optionally other events)...
Simple data structure for histogram binning.
int getNumberOfBins() const
Get number of bins.
The Vec class is a straightforward 3-d vector, which also works in pyroot.
Definition Vec.hh:13