50 template<
class JHit_t>
62 buffer.insert(JType_t(*i));
69 const char*
const muon_t =
"muon";
79int main(
int argc,
char **argv)
86 typedef JTriggeredFileScanner_t::multi_pointer_type multi_pointer_type;
88 JTriggeredFileScanner_t inputFile;
91 size_t numberOfPrefits;
102 JParser<> zap(
"Example program to histogram fit results.");
112 zap[
'O'] =
make_field(option) =
"LINE",
"LOGE",
"LINN",
"LOGN";
118 catch(
const exception& error) {
119 FATAL(error.what() << endl);
127 catch(
const exception& error) {
128 FATAL(error.what() << endl);
131 JVolume volume(head, option !=
"LINE");
136 cylinder.
add(offset);
138 NOTICE(
"Offset: " << offset << endl);
139 NOTICE(
"Cylinder: " << cylinder << endl);
143 TH1D job(
"job", NULL, 100, 0.5, 100.5);
145 TH1D hn(
"hn", NULL, 250, -0.5, 249.5);
146 TH1D hq(
"hq", NULL, 300, 0.0, 600.0);
147 TH1D hi(
"hi", NULL, 101, -0.5, 100.5);
148 TH1D hv(
"hv", NULL, 200, -6.0, 0.0);
149 TH1D h1(
"h1", NULL, 200, -2.0, +2.0);
150 TH1D hc(
"hc", NULL, 200, -1.0, +1.0);
151 TH1D hu(
"hu", NULL, 400, -1.0e3, 1.0e3);
153 TH1D hx(
"hx", NULL, 70, -3.0, +2.3);
154 TH1D hd(
"hd", NULL, 100, 0.0, 10.0);
155 TH1D ht(
"ht", NULL, 100, -100.0, 100.0);
157 TH1D hz0(
"hz0[start]", NULL, 400, -200.0, 200.0);
158 TH1D hz1(
"hz1[end]", NULL, 400, -200.0, 200.0);
163 TH1D E_E(
"E_E", NULL, 100, -5.0, +5.0);
164 TH2D ExE(
"ExE", NULL,
168 const Int_t ny = 28800;
169 const Double_t ymin = 0.0;
170 const Double_t ymax = 180.0;
174 if (option ==
"LINN") {
176 const double xmin = log((
double) 3);
177 const double xmax = log((
double) 300);
179 for (
double x = xmin, dx = (xmax - xmin) / 20; x <= xmax; x += dx) {
180 X.push_back((
int) exp(x));
183 }
else if (option ==
"LOGN") {
185 const double xmin = log10((
double) 3);
186 const double xmax = log10((
double) 400);
188 for (
double x = xmin, dx = (xmax - xmin) / 20; x <= xmax; x += dx) {
199 TH2D h2(
"h2", NULL, X.size() - 1, X.data(), ny, ymin, ymax);
212 while (inputFile.hasNext()) {
214 STATUS(
"event: " << setw(10) << inputFile.getCounter() <<
'\r');
DEBUG(endl);
216 multi_pointer_type ps = inputFile.next();
235 vector<Trk>::const_iterator muon =
event->mc_trks.end();
237 for (vector<Trk>::const_iterator track = event->mc_trks.begin(); track != event->mc_trks.end(); ++track) {
243 if (track->E > Emax) {
250 if (muon == event->mc_trks.end()) {
262 const double E =
event->mc_trks[0].E;
265 if (option ==
"LINN")
267 else if (option ==
"LOGN")
278 if (evt->begin() == __end) {
284 if (numberOfPrefits > 0) {
286 JEvt::iterator __q = __end;
288 advance(__end = evt->begin(), min(numberOfPrefits, (
size_t)
distance(evt->begin(), __q)));
303 }
else if (
primary == neutrino_t) {
308 JEvt::iterator best = pointing(evt->begin(), __end);
309 const Double_t beta = pointing.
getAngle(*best);
310 const double Efit = best->getE();
311 const double Eraw = best->getW(
JENERGY_ENERGY, numeric_limits<double>::min());
315 bool ok = (Efit >= Emin_GeV);
321 hn.Fill((Double_t) best->getNDF());
322 hq.Fill(best->getQ());
323 hi.Fill((Double_t)
distance(evt->begin(), best));
324 hc.Fill(best->getDZ());
330 hu.Fill(atmosphere(evt->begin(), __end));
333 hx.Fill(max(log10(beta), hx.GetXaxis()->GetXmin()));
374 E_0.Fill(volume.
getX(E,
true));
375 E_1.Fill(volume.
getX(Eraw,
true));
376 E_2.Fill(volume.
getX(Efit,
true));
377 E_E.Fill(volume.
getX(Efit) - volume.
getX(E));
378 ExE.Fill(volume.
getX(E), volume.
getX(Efit));
388 const double npe = best->getW(
JVETO_NPE);
390 const double pv = TMath::PoissonI(count, npe);
392 hv.Fill(max(log10(pv), hv.GetXaxis()->GetXmin()));
399 NOTICE(
"Number of events input " << setw(8) << right << job.GetBinContent(1) << endl);
400 NOTICE(
"Number of events with muon " << setw(8) << right << job.GetBinContent(3) << endl);
401 NOTICE(
"Number of events with fit " << setw(8) << right << job.GetBinContent(4) << endl);
402 NOTICE(
"Number of events selected " << setw(8) << right << job.GetBinContent(5) << endl);
403 NOTICE(
"Number of events with neutrino " << setw(8) << right << job.GetBinContent(6) << endl);
404 NOTICE(
"Number of events contained " << setw(8) << right << job.GetBinContent(7) << endl);
410 for (
double q : {0.50, 0.90, 0.99}) {
General purpose messaging.
#define DEBUG(A)
Message macros.
int main(int argc, char **argv)
Utility class to parse command line options.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
ROOT TTree parameter settings of various packages.
Synchronously read DAQ events and Monte Carlo events (and optionally other events).
double getRadius() const
Get radius.
Data structure for vector in two dimensions.
double getLengthSquared() const
Get length squared.
double getIntersection(const JVector3D &pos) const
Get longitudinal position along axis of position of closest approach with given position.
double getZmin() const
Get minimal z position.
bool is_inside(const JVector3D &pos) const
Check whether given point is inside cylinder.
double getZmax() const
Get maximal z position.
JCylinder3D & add(const JVector3D &pos)
Add position.
Data structure for position in three dimensions.
const JPosition3D & getPosition() const
Get position.
double getT(const JVector3D &pos) const
Get arrival time of Cherenkov light at given position.
JTime & add(const JTime &value)
Addition operator.
JTime & sub(const JTime &value)
Subtraction operator.
JVector3D & add(const JVector3D &vector)
Add vector.
Utility class to parse command line options.
Auxiliary class to evaluate atmospheric muon hypothesis.
Auxiliary class to compare fit results with respect to a reference direction (e.g....
double getAngle(const JFit &fit) const
Get angle between reference direction and fit result.
const_iterator< T > end() const
Get end of data.
const_iterator< T > begin() const
Get begin of data.
Auxiliary class to convert DAQ hit time to/from Monte Carlo hit time.
double putTime() const
Get Monte Carlo time minus DAQ/trigger time.
static const int JENERGY_ENERGY
uncorrected energy [GeV] from JMuonEnergy
static const int JSTART_NPE_MIP_TOTAL
number of photo-electrons along the whole track from JMuonStart
static const int JVETO_NPE
number of photo-electrons from JVeto.cc
static const int JSTART_LENGTH_METRES
distance between projected positions on the track of optical modules for which the response does not ...
static const int JVETO_NUMBER_OF_HITS
number of hits from JVeto.cc
Vec getOrigin(const JHead &header)
Get origin.
JDirection3D getDirection(const Vec &dir)
Get direction.
JTrack3E getTrack(const Trk &track)
Get track.
JCylinder3D getCylinder(const JHead &header)
Get cylinder corresponding to the positions of generated tracks (i.e.
bool has_neutrino(const Evt &evt)
Test whether given event has an incoming neutrino.
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.
const Trk & get_neutrino(const Evt &evt)
Get incoming neutrino.
size_t getCount(const array_type< T > &buffer, const JCompare_t &compare)
Count number of unique values.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
counter_type advance(counter_type &counter, const counter_type value, const counter_type limit=std::numeric_limits< counter_type >::max())
Advance counter.
Head getCommonHeader(const JMultipleFileScanner_t &file_list)
Get common Monte Carlo header.
KM3NeT DAQ data structures and auxiliaries.
const char *const neutrino_t
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
Auxiliary data structure for floating point format specification.
Auxiliary class for histogramming of effective volume.
Double_t getX(const Double_t E, double constrain=false) const
Get abscissa value.
Double_t getXmax() const
Get maximal abscissa value.
Double_t getXmin() const
Get minimal abscissa value.
Auxiliary class to test history.
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
General purpose sorter of fit results.
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 Trk class represents a Monte Carlo (MC) particle as well as a reconstructed track/shower.
double E
Energy [GeV] (either MC truth or reconstructed)
Reconstruction type dependent comparison of track quality.
Auxiliary include file for time conversion between DAQ/trigger hit and Monte Carlo hit.