Jpp 20.0.0-27-g39925593c-D
the software that should make you happy
Loading...
Searching...
No Matches
Functions
JMuonPostfit.cc File Reference

Example program to histogram fit results for the muon reconstruction chain. More...

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <set>
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TProfile.h"
#include "TMath.h"
#include "km3net-dataformat/offline/Head.hh"
#include "km3net-dataformat/offline/Evt.hh"
#include "km3net-dataformat/online/JDAQModuleIdentifier.hh"
#include "km3net-dataformat/definitions/reconstruction.hh"
#include "km3net-dataformat/definitions/fitparameters.hh"
#include "km3net-dataformat/tools/time_converter.hh"
#include "JAAnet/JHead.hh"
#include "JAAnet/JHeadToolkit.hh"
#include "JAAnet/JAAnetToolkit.hh"
#include "JAAnet/JVolume.hh"
#include "JDAQ/JDAQEventIO.hh"
#include "JPhysics/JConstants.hh"
#include "JTools/JQuantile.hh"
#include "JSupport/JTriggeredFileScanner.hh"
#include "JSupport/JMonteCarloFileSupportkit.hh"
#include "JSupport/JSupport.hh"
#include "JReconstruction/JEvt.hh"
#include "JReconstruction/JEvtToolkit.hh"
#include "JLang/JPredicate.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

Example program to histogram fit results for the muon reconstruction chain.

Author
mdejong, bofearraigh.

Definition in file JMuonPostfit.cc.

Function Documentation

◆ main()

int main ( int argc,
char ** argv )

Definition at line 79 of file JMuonPostfit.cc.

80{
81 using namespace std;
82 using namespace JPP;
83 using namespace KM3NETDAQ;
84
85 typedef JTriggeredFileScanner<JEvt> JTriggeredFileScanner_t;
86 typedef JTriggeredFileScanner_t::multi_pointer_type multi_pointer_type;
87
88 JTriggeredFileScanner_t inputFile;
89 JLimit_t& numberOfEvents = inputFile.getLimit();
90 string outputFile;
91 size_t numberOfPrefits;
92 JQualitySorter quality_sorter;
93 JAtmosphericMuon atmosphere;
94 double Emin_GeV;
95 int application;
96 string option;
97 string primary;
98 int debug;
99
100 try {
101
102 JParser<> zap("Example program to histogram fit results.");
103
104 zap['f'] = make_field(inputFile);
105 zap['o'] = make_field(outputFile) = "postfit.root";
106 zap['n'] = make_field(numberOfEvents) = JLimit::max();
107 zap['N'] = make_field(numberOfPrefits) = 1;
108 zap['L'] = make_field(quality_sorter) = JPARSER::initialised();
109 zap['E'] = make_field(Emin_GeV) = 0.0;
110 zap['A'] = make_field(application) = JMUONPREFIT, JMUONSIMPLEX, JMUONGANDALF, JMUONENERGY, JMUONSTART, JMUONPATH;
111 zap['a'] = make_field(atmosphere) = JAtmosphericMuon(90.0, 90.0);
112 zap['O'] = make_field(option) = "LINE", "LOGE", "LINN", "LOGN";
113 zap['p'] = make_field(primary) = muon_t, neutrino_t;
114 zap['d'] = make_field(debug) = 2;
115
116 zap(argc, argv);
117 }
118 catch(const exception& error) {
119 FATAL(error.what() << endl);
120 }
121
122 JHead head;
123
124 try {
125 head = getCommonHeader(inputFile);
126 }
127 catch(const exception& error) {
128 FATAL(error.what() << endl);
129 }
130
131 JVolume volume(head, option != "LINE");
132 JPosition3D offset = getPosition(getOffset(head));
133 JPosition3D origin = getPosition(getOrigin(head));
134 JCylinder3D cylinder = getCylinder(head);
135
136 cylinder.add(offset);
137
138 NOTICE("Offset: " << offset << endl);
139 NOTICE("Cylinder: " << cylinder << endl);
140
141 TFile out(outputFile.c_str(), "recreate");
142
143 TH1D job("job", NULL, 100, 0.5, 100.5);
144
145 TH1D hn("hn", NULL, 250, -0.5, 249.5); // NDF
146 TH1D hq("hq", NULL, 300, 0.0, 600.0); // quality parameter
147 TH1D hi("hi", NULL, 101, -0.5, 100.5); // index of the best event in order of descending angular resolution
148 TH1D hc("hc", NULL, 200, -1.0, +1.0); // dz (z-slope)
149 TH1D hu("hu", NULL, 400, -1.0e3, 1.0e3); // quality difference between best up-going and down-going track (check of atmospheric muon)
150
151 TH1D hx("hx", NULL, 70, -3.0, +2.3); // angular deviation [log(deg)]
152 TH1D hd("hd", NULL, 100, 0.0, 10.0); // distance between true and best event at intersection point [m]
153 TH1D ht("ht", NULL, 100, -100.0, 100.0); // time between true and best track [ns]
154
155 TH1D hz0("hz0[start]", NULL, 400, -200.0, 200.0); // start point [m]
156 TH1D hz1("hz1[end]", NULL, 400, -200.0, 200.0); // end point [m]
157
158 TH1D E_0("E_0", NULL, 100, volume.getXmin(), volume.getXmax()); // true muon energy
159 TH1D E_1("E_1", NULL, 100, volume.getXmin(), volume.getXmax()); // uncorrected energy of best track
160 TH1D E_2("E_2", NULL, 100, volume.getXmin(), volume.getXmax()); // corrected energy of best track
161 TH1D E_E("E_E", NULL, 100, -5.0, +5.0); // ratio of reconstructed energy and true energy [log(E/E)]
162 TH2D ExE("ExE", NULL,
163 40, volume.getXmin(), volume.getXmax(),
164 40, volume.getXmin(), volume.getXmax()); // (Etrue,Ereco)
165
166 const Int_t ny = 28800;
167 const Double_t ymin = 0.0; // [deg]
168 const Double_t ymax = 180.0; // [deg]
169
170 vector<double> X;
171
172 if (option == "LINN") {
173
174 const double xmin = log((double) 3);
175 const double xmax = log((double) 300);
176
177 for (double x = xmin, dx = (xmax - xmin) / 20; x <= xmax; x += dx) {
178 X.push_back((int) exp(x));
179 }
180
181 } else if (option == "LOGN") {
182
183 const double xmin = log10((double) 3);
184 const double xmax = log10((double) 400);
185
186 for (double x = xmin, dx = (xmax - xmin) / 20; x <= xmax; x += dx) {
187 X.push_back(x);
188 }
189
190 } else {
191
192 for (double x = volume.getXmin(), dx = (volume.getXmax() - volume.getXmin()) / 20.0; x < volume.getXmax() + 0.5*dx; x += dx) {
193 X.push_back(x);
194 }
195 }
196
197 TH2D h2("h2", NULL, X.size() - 1, X.data(), ny, ymin, ymax);
198
199 TH2D Va("Va", NULL,
200 100, 0.0, cylinder.getRadius()*cylinder.getRadius(),
201 100, cylinder.getZmin() - 50.0, cylinder.getZmax() + 50.0);
202 TH2D Vb("Vb", NULL,
203 100, 0.0, cylinder.getRadius()*cylinder.getRadius(),
204 100, cylinder.getZmin() - 50.0, cylinder.getZmax() + 50.0);
205
206 JQuantile Q("Angle");
207 JQuantile O("Omega");
208
209
210 while (inputFile.hasNext()) {
211
212 STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);
213
214 multi_pointer_type ps = inputFile.next();
215
216 JDAQEvent* tev = ps;
217 JEvt* evt = ps;
218 Evt* event = ps;
219
220 const time_converter converter(*event, *tev);
221
222 job.Fill(1.0);
223
224
225 double Enu = 0.0;
226 double Emu = 0.0;
227 double Emax = 0.0;
228
229 if (has_neutrino(*event)) {
230 Enu = get_neutrino(*event).E;
231 }
232
233 vector<Trk>::const_iterator muon = event->mc_trks.end();
234
235 for (vector<Trk>::const_iterator track = event->mc_trks.begin(); track != event->mc_trks.end(); ++track) {
236
237 if (is_muon(*track)) {
238
239 Emu += track->E;
240
241 if (track->E > Emax) {
242 muon = track;
243 Emax = track->E;
244 }
245 }
246 }
247
248 if (muon == event->mc_trks.end()) {
249 continue;
250 }
251
252 job.Fill(3.0);
253
254
255 // abscissa
256
257 Double_t x = 0.0;
258
259 {
260 const double E = event->mc_trks[0].E;
261 const int N = getCount(tev->begin<JDAQTriggeredHit>(), tev->end<JDAQTriggeredHit>());
262
263 if (option == "LINN")
264 x = N;
265 else if (option == "LOGN")
266 x = log10(N);
267 else
268 x = volume.getX(E);
269 }
270
271
272 if (!evt->empty()) {
273
274 JEvt::iterator __end = partition(evt->begin(), evt->end(), JHistory::is_application(application));
275
276 if (evt->begin() == __end) {
277 continue;
278 }
279
280 job.Fill(4.0);
281
282 if (numberOfPrefits > 0) {
283
284 JEvt::iterator __q = __end;
285
286 advance(__end = evt->begin(), min(numberOfPrefits, (size_t) distance(evt->begin(), __q)));
287
288 partial_sort(evt->begin(), __end, __q, quality_sorter);
289
290 } else {
291
292 sort(evt->begin(), __end, quality_sorter);
293 }
294
295 double E = 0.0;
296 JPointing pointing;
297
298 if (primary == muon_t) {
299 E = Emu;
300 pointing = JPointing(getDirection(*muon));
301 } else if (primary == neutrino_t) {
302 E = Enu;
303 pointing = JPointing(getDirection(get_neutrino(*event)));
304 }
305
306 JEvt::iterator best = pointing(evt->begin(), __end);
307 const Double_t beta = pointing.getAngle(*best);
308 const double Efit = best->getE();
309 const double Eraw = best->getW(JENERGY_ENERGY, numeric_limits<double>::min());
310
311 // selection of fit result
312
313 bool ok = (Efit >= Emin_GeV);
314
315 if (ok) {
316
317 job.Fill(5.0);
318
319 hn.Fill((Double_t) best->getNDF());
320 hq.Fill(best->getQ());
321 hi.Fill((Double_t) distance(evt->begin(), best));
322 hc.Fill(best->getDZ());
323
324 // atmospheric muon hypothesis test.
325
326 if (( has_neutrino(*event) && get_neutrino(*event).dir.z >= atmosphere.dot2) || // difference in quality between best upwards and best downgoing track.
327 (!has_neutrino(*event) && best->getDZ() >= atmosphere.dot2)) { // negative (positive) values imply a down(up) going track.
328 hu.Fill(atmosphere(evt->begin(), __end));
329 }
330
331 hx.Fill(max(log10(beta), hx.GetXaxis()->GetXmin()));
332
333 Q.put(beta);
334
335 JTrack3E ta = getTrack(*muon);
336 JTrack3E tb = getTrack(*best);
337
338 ta.add(offset);
339 tb.sub(converter.putTime());
340
341 static_cast<JTrack3D&>(ta).move(ta.getIntersection(tb), getSpeedOfLight()); // move to intersection point of both tracks
342 static_cast<JTrack3D&>(tb).move(tb.getIntersection(ta), getSpeedOfLight());
343
344 hd.Fill((tb.getPosition() - ta.getPosition()).getLength());
345
346 if (has_neutrino(*event)) {
347
348 job.Fill(6.0);
349
350 const Trk neutrino = get_neutrino(*event);
351
352 JPosition3D vertex = getPosition(neutrino);
353
354 vertex.add(offset);
355
356 if (cylinder.is_inside(vertex)) {
357
358 job.Fill(7.0);
359
360 Va.Fill(JVector2D(best->getX() - origin.getX(), best->getY() - origin.getY()).getLengthSquared(), best->getZ()); //R^{2} [m^2], z [m]
361
362 JTrack3E tc = getTrack(*best);
363
364 hz0.Fill(tc.getIntersection(vertex));
365 hz1.Fill(best->getW(JSTART_LENGTH_METRES, 0.0) - gWater(muon->E));
366 }
367 }
368
369 Vb.Fill(JVector2D(best->getX() - origin.getX(), best->getY() - origin.getY()).getLengthSquared(), best->getZ()); //R^{2} [m^2], z [m]
370 ht.Fill(tb.getT() - ta.getT());
371
372 E_0.Fill(volume.getX(E, true));
373 E_1.Fill(volume.getX(Eraw, true));
374 E_2.Fill(volume.getX(Efit, true));
375 E_E.Fill(volume.getX(Efit) - volume.getX(E));
376 ExE.Fill(volume.getX(E), volume.getX(Efit));
377
378 h2.Fill(x, beta);
379 }
380 }
381 }
382 STATUS(endl);
383
384 NOTICE("Number of events input " << setw(8) << right << job.GetBinContent(1) << endl);
385 NOTICE("Number of events with muon " << setw(8) << right << job.GetBinContent(3) << endl);
386 NOTICE("Number of events with fit " << setw(8) << right << job.GetBinContent(4) << endl);
387 NOTICE("Number of events selected " << setw(8) << right << job.GetBinContent(5) << endl);
388 NOTICE("Number of events with neutrino " << setw(8) << right << job.GetBinContent(6) << endl);
389 NOTICE("Number of events contained " << setw(8) << right << job.GetBinContent(7) << endl);
390
391 if (Q.getCount() != 0) {
392
393 NOTICE("Median space angle [deg] " << FIXED (6,3) << Q.getQuantile(0.5) << endl);
394
395 for (double q : {0.50, 0.90, 0.99}) {
396 NOTICE("Space angle " << FIXED(5,1) << (100.0 * q) << "% quantile [deg] " << FIXED(6,3) << Q.getQuantile(q) << endl);
397 }
398 }
399
400 out.Write();
401 out.Close();
402}
outputFile
string outputFile
Definition JDAQTimesliceSelector.cc:37
DEBUG
#define DEBUG(A)
Message macros.
Definition JMessage.hh:62
STATUS
#define STATUS(A)
Definition JMessage.hh:63
NOTICE
#define NOTICE(A)
Definition JMessage.hh:64
FATAL
#define FATAL(A)
Definition JMessage.hh:67
debug
int debug
debug level
Definition JSirene.cc:72
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2142
distance
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
Definition build/antares-dataformat/PhysicsEvent.hh:434
JAANET::JHead
Monte Carlo run header.
Definition JHead.hh:1236
JGEOMETRY2D::JVector2D
Data structure for vector in two dimensions.
Definition JVector2D.hh:34
JGEOMETRY2D::JVector2D::getLengthSquared
double getLengthSquared() const
Get length squared.
Definition JVector2D.hh:188
JGEOMETRY3D::JPosition3D
Data structure for position in three dimensions.
Definition JPosition3D.hh:38
JGEOMETRY3D::JTrack3E
3D track with energy.
Definition JTrack3E.hh:34
JOSCPROB::JOscProbInterpolator
Template definition of a multi-dimensional oscillation probability interpolation table.
Definition JOscProbInterpolator.hh:71
JRECONSTRUCTION::JAtmosphericMuon
Auxiliary class to evaluate atmospheric muon hypothesis.
Definition JReconstruction/JEvtToolkit.hh:929
JRECONSTRUCTION::JPointing
Auxiliary class to compare fit results with respect to a reference direction (e.g....
Definition JReconstruction/JEvtToolkit.hh:688
time_converter
Auxiliary class to convert DAQ hit time to/from Monte Carlo hit time.
Definition time_converter.hh:20
JENERGY_ENERGY
static const int JENERGY_ENERGY
uncorrected energy [GeV] see JRECONSTRUCTION::JMuonEnergy
Definition fitparameters.hh:16
JSTART_LENGTH_METRES
static const int JSTART_LENGTH_METRES
distance between projected positions on the track of optical modules for which the response does not ...
Definition fitparameters.hh:27
JAANET::getOrigin
Vec getOrigin(const JHead &header)
Get origin.
Definition JHeadToolkit.hh:325
JAANET::getDirection
JDirection3D getDirection(const Vec &dir)
Get direction.
Definition JAAnetToolkit.hh:218
JAANET::getTrack
JTrack3E getTrack(const Trk &track)
Get track.
Definition JAAnetToolkit.hh:278
JAANET::getCylinder
JCylinder3D getCylinder(const JHead &header)
Get cylinder corresponding to the positions of generated tracks (i.e.
Definition JHeadToolkit.hh:347
JAANET::has_neutrino
bool has_neutrino(const Evt &evt)
Test whether given event has an incoming neutrino.
Definition JAAnetToolkit.hh:520
JAANET::getPosition
JPosition3D getPosition(const Vec &pos)
Get position.
Definition JAAnetToolkit.hh:182
JAANET::is_muon
bool is_muon(const Trk &track)
Test whether given track is a (anti-)muon.
Definition JAAnetToolkit.hh:366
JAANET::getOffset
Vec getOffset(const JHead &header)
Get offset.
Definition JHeadToolkit.hh:291
JAANET::get_neutrino
const Trk & get_neutrino(const Evt &evt)
Get incoming neutrino.
Definition JAAnetToolkit.hh:542
JLANG::getCount
size_t getCount(const array_type< T > &buffer, const JCompare_t &compare)
Count number of unique values.
Definition JVectorize.hh:262
JPHYSICS::getSpeedOfLight
const double getSpeedOfLight()
Get speed of light.
Definition JPhysics/JConstants.hh:112
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition JAAnetToolkit.hh:45
JROOT::advance
counter_type advance(counter_type &counter, const counter_type value, const counter_type limit=std::numeric_limits< counter_type >::max())
Advance counter.
Definition JROOT/JCounter.hh:35
JSUPPORT::getCommonHeader
Head getCommonHeader(const JMultipleFileScanner_t &file_list)
Get common Monte Carlo header.
Definition JMonteCarloFileSupportkit.hh:559
KM3NETDAQ
KM3NeT DAQ data structures and auxiliaries.
Definition DataQueue.cc:39
create_update_issues.origin
str origin
Definition create_update_issues.py:4
mc_keys::neutrino_t
const char *const neutrino_t
Definition io_ascii.hh:29
Evt
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
Definition Evt.hh:21
FIXED
Auxiliary data structure for floating point format specification.
Definition JManip.hh:448
JAANET::JVolume
Auxiliary class for histogramming of effective volume.
Definition JVolume.hh:29
JAANET::primary
Primary particle.
Definition JHead.hh:1174
JACOUSTICS::JEvt
Acoustic event fit.
Definition JAcoustics/JEvt.hh:310
JFIT::JHistory::is_application
Auxiliary class to test history.
Definition JHistory.hh:188
JPARSER::initialised
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition JParser.hh:68
JRECONSTRUCTION::JQualitySorter
General purpose sorter of fit results.
Definition JReconstruction/JEvtToolkit.hh:426
JSUPPORT::JLimit
Auxiliary class for defining the range of iterations of objects.
Definition JLimit.hh:45
JSUPPORT::JLimit::max
static counter_type max()
Get maximum counter value.
Definition JLimit.hh:128
JTOOLS::JQuantile
Auxiliary data structure for running average, standard deviation and quantiles.
Definition JQuantile.hh:41
Trk
The Trk class represents a Monte Carlo (MC) particle as well as a reconstructed track/shower.
Definition Trk.hh:15
quality_sorter
Reconstruction type dependent comparison of track quality.
Definition tools/reconstruction.hh:85
Generated by doxygen 1.10.0