80{
84
86 typedef JTriggeredFileScanner_t::multi_pointer_type multi_pointer_type;
87
89 JLimit_t& numberOfEvents = inputFile.getLimit();
91 size_t numberOfPrefits;
94 double Emin_GeV;
95 double NPE;
97 string option;
100 double radius;
103
104 try {
105
106 JParser<> zap(
"Example program to histogram fit results.");
107
123
125 }
126 catch(const exception& error) {
128 }
129
131
133
134 try {
136 }
137 catch(const exception& error) {
139 }
140
145
147
149
151
156
158
160
162
163 TH1D hn(
"hn",
"NDF (Number of Degrees of Freedom)", 250, -0.5, 249.5);
164 TH1D hq(
"hq",
"Fit Quality", 300, 0.0, 600.0);
165 TH1D hi(
"hi",
"Fit Index vs Best Resolution", 101, -0.5, 100.5);
166 TH1D hd(
"hd",
"Distance between Reco and True position", 2000, 0.0, 100.0);
167 TH1D hz(
"hz",
"Longitudinal Distance between Reco and MC lepton position", 100, -200.0, 200.0);
168 TH1D ht(
"ht",
"Time difference between Reco and MC lepton", 100, 0, 100.0);
169 TH1D h4D(
"h4D",
"4D-distance between reco and true vertex", 2000, 0.0, 100.0);
170
171 TH1D ha(
"ha",
"Angle between reco direction and true direction", 1000, 0.0, 180.0);
172
173 TH1D e0(
"e0",
"True lepton energy", 100, volume.getXmin(), volume.getXmax());
174 TH1D e2(
"e2",
"Reconstructed energy", 100, volume.getXmin(), volume.getXmax());
175 TH1D er(
"er",
"Ratio of reconstructed energy and true energy", 100, -5.0, +5.0);
176 TH1D ea(
"ea",
"Distribution of Energy error for all the events; E_{reco}-E_{true}", 100, -30, +30);
177 TH1D ed3_5GeV(
"ed3_5GeV",
"Distribution of Energy error for events with MC energy #in [3,5] GeV; E_{reco}-E_{true}", 100, -30, +30);
178 TH1D ed8_11GeV(
"ed8_11GeV",
"Distribution of Energy error for events with MC energy #in [8,11] GeV; E_{reco}-E_{true}", 100, -30, +30);
179 TH1D ed15_21GeV(
"ed15_21GeV",
"Distribution of Energy error for events with MC energy #in [15,21] GeV; E_{reco}-E_{true}", 100, -30, +30);
180 TH2D ee(
"ee",
"; E_{true} [GeV]; E_{reco} [GeV]",
181 100, volume.getXmin(), volume.getXmax(),
182 100, volume.getXmin(), volume.getXmax());
186
188
189 if (option.find('E') != string::npos) {
190
191 for (
double x = volume.getXmin();
x <= volume.getXmax();
x += (volume.getXmax() - volume.getXmin()) / 20) {
192 X.push_back(x);
193 }
194
195 } else {
196
198
199 for ( ;
x <= 15.5;
x += 1.0) { X.push_back(x); }
200 for ( ;
x <= 25.5;
x += 2.0) { X.push_back(x); }
201 for ( ;
x <= 50.5;
x += 5.0) { X.push_back(x); }
202 for ( ;
x <= 100.5;
x += 10.0) { X.push_back(x); }
203 for ( ;
x <= 250.5;
x += 20.0) { X.push_back(x); }
204 }
205
207 TH2D h3(
"h3",
NULL, X.size() - 1, X.data(), 2000, 0.0, 100.0);
209 TProfile hfracE(
"hfracE",
"Fractional Energy Error", X.size() - 1, X.data());
210 TProfile he(
"he",
"Reconstruction Efficiency", X.size() - 1, X.data());
211 TProfile htheta_nu_lep(
"htheta_nu_lep",
"Angle between neutrino and primary lepton", X.size() - 1, X.data());
212 TH1D hln1d(
"hln1d",
"Angle between neutrino and lepton",40,0,4);
213 TH2D hb(
"hb",
NULL, 100, 0.0, 0.5e6, 100, -100.0, +900.0);
214 TH2S hmca(
"hmca",
NULL, X.size() - 1, X.data(), 1000, 0, 180);
216 TH2D hY(
"hY",
"Reco Bjorken Y vs MC Bjorken Y", 40, 0, 1, 40, 0, 1);
217 TH2D hby3_5GeV(
"hby3_5GeV",
"Reco Bjorken Y vs MC Bjorken Y for events with Reco E #in [3, 5] GeV", 40, 0, 1, 40, 0, 1);
218 TH2D hby8_10GeV(
"hby8_10GeV",
"Reco Bjorken Y vs MC Bjorken Y for events with Reco E #in [8, 10] GeV", 40, 0, 1, 40, 0, 1);
219
226 while (inputFile.hasNext()) {
227
229
230 multi_pointer_type ps = inputFile.next();
231
235
237
239
243
245
247
250
251 }
252
254
256
258
260
264 }
265
266
268
273 }
274 }
275 }
276
277 if (
lepton == event->mc_trks.end()) {
278 continue;
279 }
280
282
284
285
286
288
289 if (option.find('E') != string::npos){
290
293
294 } else{
296 }
297
298
300
301 if (!evt->empty()) {
302
304
305 if (evt->begin() == __end) {
306 continue;
307 }
308
310
311
312 if (numberOfPrefits > 0 ) {
313
314 JEvt::iterator
__q = __end;
315
316 advance(__end = evt->begin(), min(numberOfPrefits, (
size_t)
distance(evt->begin(),
__q)));
317
319
320 } else {
321
323 }
324
325 JEvt::iterator
best = evt->begin();
326
330
332
339 }
340
342
347 } else {
349 }
351
354
356
357
358 bool ok = (
Efit >= Emin_GeV);
359
361
362 W = 1.0;
363
365
370
372
375 } else {
377 }
378
380
382
383 double zenith =
tb.getDirection().getTheta() * 180 /
PI;
384
385
391
404 } else {
415 }
416
418
420
422
424
426
428
429 hz.Fill(
tb.getIntersection(
mc_vx));
430 }
431 }
432
436
438 }
440
442
445
446 }
448
449 e0.Fill(volume.getX(
Enu,
true));
450 er.Fill(volume.getX(
Efit) - volume.getX(
Enu));
451 ee.Fill(volume.getX(
Enu), volume.getX(
Efit));
456 hfracE.Fill(x, fabs(volume.getX(
Efit) - volume.getX(
Enu))/volume.getX(
Enu));
457
461
462 } else {
463
464 e0.Fill(volume.getX(
Elep,
true));
465 er.Fill(volume.getX(
Efit) - volume.getX(
Elep));
466 ee.Fill(volume.getX(
Elep), volume.getX(
Efit));
472
476
477 }
478 e2.Fill(volume.getX(
Efit,
true));
480 }
481 }
482 }
484
486
487 NOTICE(
"Number of events input " <<
setw(8) << right <<
job.GetBinContent(1) <<
endl);
489 NOTICE(
"Number of events with electron " <<
setw(8) << right <<
job.GetBinContent(3) <<
endl);
490 } else{
491 NOTICE(
"Number of events with muon " <<
setw(8) << right <<
job.GetBinContent(3) <<
endl);
492 }
493 NOTICE(
"Number of events with fit " <<
setw(8) << right <<
job.GetBinContent(4) <<
endl);
494 NOTICE(
"Number of events selected " <<
setw(8) << right <<
job.GetBinContent(5) <<
endl);
495 NOTICE(
"Number of events with neutrino " <<
setw(8) << right <<
job.GetBinContent(6) <<
endl);
496 NOTICE(
"Number of events contained " <<
setw(8) << right <<
job.GetBinContent(7) <<
endl);
497
498 if (Q.getCount() != 0) {
499 NOTICE(
"Median, 70, 80 and 90% quantiles of space angle [deg] " <<
FIXED (6,3) << Q.getQuantile(0.5) <<
" "<< Q.getQuantile(0.7) <<
" "<< Q.getQuantile(0.8)<<
" "<< Q.getQuantile(0.9)<<
endl);
500 }
501 if (QE.getCount() != 0) {
502 NOTICE(
"Median, 70%, 80% and 90% quantiles of energy ratio log(Ereco/Etrue) " <<
FIXED (6,3) << QE.getQuantile(0.5) <<
" "<< QE.getQuantile(0.7) <<
" "<< QE.getQuantile(0.8)<<
" "<< QE.getQuantile(0.9) <<
endl);
503 }
504 if(
Qp.getCount() != 0){
505 NOTICE(
"Median distance to vertex: " <<
setw(8) <<
Qp.getQuantile(0.5) <<
" "<<
Qp.getQuantile(0.7)<<
" "<<
Qp.getQuantile(0.8)<<
" "<<
Qp.getQuantile(0.9)<<
endl);
506 NOTICE(
"Median time to vertex: " <<
setw(8) << Qt.getQuantile(0.5) <<
" "<< Qt.getQuantile(0.7)<<
" "<< Qt.getQuantile(0.8)<<
" "<<Qt.getQuantile(0.9)<<
endl);
507 NOTICE(
"Median 4D-distance to vertex: " <<
setw(8) <<
Q4d.getQuantile(0.5) <<
" "<<
Q4d.getQuantile(0.7)<<
" "<<
Q4d.getQuantile(0.8)<<
" "<<
Q4d.getQuantile(0.9)<<
endl);
508 }
509 out.Write();
510 out.Close();
511}
#define DEBUG(A)
Message macros.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Data structure for circle in two dimensions.
Data structure for position in two dimensions.
Data structure for vector in two dimensions.
double getLengthSquared() const
Get length squared.
Data structure for position in three dimensions.
Data structure for vector in three dimensions.
Template definition of a multi-dimensional oscillation probability interpolation table.
Auxiliary class to evaluate atmospheric muon hypothesis.
Auxiliary class to compare fit results with respect to a reference direction (e.g....
Auxiliary class to compare fit results with respect to a reference position.
Auxiliary class to compare fit results with respect to a reference energy.
Auxiliary class to convert DAQ hit time to/from Monte Carlo hit time.
Vec getOrigin(const JHead &header)
Get origin.
JDirection3D getDirection(const Vec &dir)
Get direction.
bool is_electron(const Trk &track)
Test whether given track is a (anti-)electron.
JShower3E getShower(const Trk &shower)
Get shower.
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.
double getAngle(const JQuaternion3D &first, const JQuaternion3D &second)
Get space angle between quanternions.
size_t getCount(const array_type< T > &buffer, const JCompare_t &compare)
Count number of unique values.
T pow(const T &x, const double y)
Power .
static const double PI
Mathematical constants.
static const JGeanz geanz(1.85, 0.62, 0.54)
Function object for longitudinal EM-shower profile.
double getIndexOfRefraction()
Get average index of refraction of water corresponding to group velocity.
const double getInverseSpeedOfLight()
Get inverse speed of light.
const double getSpeedOfLight()
Get speed of light.
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.
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.
static counter_type max()
Get maximum counter value.
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.