80{
84
86 typedef JTriggeredFileScanner_t::multi_pointer_type multi_pointer_type;
87
88 JTriggeredFileScanner_t inputFile;
91 size_t numberOfPrefits;
94 double Emin_GeV;
95 int application;
96 string option;
99
100 try {
101
102 JParser<> zap(
"Example program to histogram fit results.");
103
112 zap[
'O'] =
make_field(option) =
"LINE",
"LOGE",
"LINN",
"LOGN";
115
116 zap(argc, argv);
117 }
118 catch(const exception& error) {
119 FATAL(error.what() << endl);
120 }
121
123
124 try {
126 }
127 catch(const exception& error) {
128 FATAL(error.what() << endl);
129 }
130
131 JVolume volume(head, option !=
"LINE");
135
136 cylinder.
add(offset);
137
138 NOTICE(
"Offset: " << offset << endl);
139 NOTICE(
"Cylinder: " << cylinder << endl);
140
142
143 TH1D job("job", NULL, 100, 0.5, 100.5);
144
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);
152
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);
156
157 TH1D hz0("hz0[start]", NULL, 400, -200.0, 200.0);
158 TH1D hz1("hz1[end]", NULL, 400, -200.0, 200.0);
159
160 TH1D E_0("E_0", NULL, 100, volume.getXmin(), volume.getXmax());
161 TH1D E_1("E_1", NULL, 100, volume.getXmin(), volume.getXmax());
162 TH1D E_2("E_2", NULL, 100, volume.getXmin(), volume.getXmax());
163 TH1D E_E("E_E", NULL, 100, -5.0, +5.0);
164 TH2D ExE("ExE", NULL,
165 40, volume.getXmin(), volume.getXmax(),
166 40, volume.getXmin(), volume.getXmax());
167
168 const Int_t ny = 28800;
169 const Double_t ymin = 0.0;
170 const Double_t ymax = 180.0;
171
173
174 if (option == "LINN") {
175
176 const double xmin = log((double) 3);
177 const double xmax = log((double) 300);
178
179 for (
double x = xmin, dx = (xmax - xmin) / 20;
x <= xmax;
x += dx) {
180 X.push_back((int) exp(x));
181 }
182
183 } else if (option == "LOGN") {
184
185 const double xmin = log10((double) 3);
186 const double xmax = log10((double) 400);
187
188 for (
double x = xmin, dx = (xmax - xmin) / 20;
x <= xmax;
x += dx) {
189 X.push_back(x);
190 }
191
192 } else {
193
194 for (double x = volume.getXmin(), dx = (volume.getXmax() - volume.getXmin()) / 20.0; x < volume.getXmax() + 0.5*dx; x += dx) {
195 X.push_back(x);
196 }
197 }
198
199 TH2D h2("h2", NULL, X.size() - 1, X.data(), ny, ymin, ymax);
200
201 TH2D Va("Va", NULL,
204 TH2D Vb("Vb", NULL,
207
210
211
212 while (inputFile.hasNext()) {
213
214 STATUS(
"event: " << setw(10) << inputFile.getCounter() <<
'\r');
DEBUG(endl);
215
216 multi_pointer_type ps = inputFile.next();
217
221
223
224 job.Fill(1.0);
225
226
227 double Enu = 0.0;
228 double Emu = 0.0;
229 double Emax = 0.0;
230
233 }
234
235 vector<Trk>::const_iterator muon = event->mc_trks.end();
236
237 for (vector<Trk>::const_iterator track = event->mc_trks.begin(); track != event->mc_trks.end(); ++track) {
238
240
241 Emu += track->E;
242
243 if (track->E > Emax) {
244 muon = track;
245 Emax = track->E;
246 }
247 }
248 }
249
250 if (muon == event->mc_trks.end()) {
251 continue;
252 }
253
254 job.Fill(3.0);
255
256
257
258
260
261 {
262 const double E = event->mc_trks[0].E;
264
265 if (option == "LINN")
267 else if (option == "LOGN")
269 else
271 }
272
273
274 if (!evt->empty()) {
275
277
278 if (evt->begin() == __end) {
279 continue;
280 }
281
282 job.Fill(4.0);
283
284 if (numberOfPrefits > 0) {
285
286 JEvt::iterator __q = __end;
287
288 advance(__end = evt->begin(), min(numberOfPrefits, (
size_t)
distance(evt->begin(), __q)));
289
291
292 } else {
293
295 }
296
297 double E = 0.0;
299
301 E = Emu;
303 }
else if (
primary == neutrino_t) {
304 E = Enu;
306 }
307
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());
312
313
314
315 bool ok = (Efit >= Emin_GeV);
316
317 if (ok) {
318
319 job.Fill(5.0);
320
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());
325
326
327
330 hu.Fill(atmosphere(evt->begin(), __end));
331 }
332
333 hx.Fill(max(log10(beta), hx.GetXaxis()->GetXmin()));
334
335 Q.put(beta);
336
339
341 tb.
sub(converter.putTime());
342
345
347
349
350 job.Fill(6.0);
351
353
355
357
359
360 job.Fill(7.0);
361
363
365
368 }
369 }
370
373
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));
379
380 h2.Fill(x, beta);
381
384 }
385
387
388 const double npe = best->getW(
JVETO_NPE);
390 const double pv = TMath::PoissonI(count, npe);
391
392 hv.Fill(max(log10(pv), hv.GetXaxis()->GetXmin()));
393 }
394 }
395 }
396 }
398
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);
405
406 if (Q.getCount() != 0) {
407
408 NOTICE(
"Median space angle [deg] " <<
FIXED (6,3) << Q.getQuantile(0.5) << endl);
409
410 for (double q : {0.50, 0.90, 0.99}) {
411 NOTICE(
"Space angle " <<
FIXED(5,1) << (100.0 * q) <<
"% quantile [deg] " <<
FIXED(6,3) << Q.getQuantile(q) << endl);
412 }
413 }
414
415 out.Write();
416 out.Close();
417}
#define DEBUG(A)
Message macros.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
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.
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.
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.