Jpp master_rocky-44-g75b7c4f75
the software that should make you happy
Loading...
Searching...
No Matches
Functions
JAAPostfit.cc File Reference

Example program to analyse fit results from Evt formatted data. More...

#include <string>
#include <iostream>
#include <iomanip>
#include <cmath>
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "km3net-dataformat/definitions/weightlist.hh"
#include "km3net-dataformat/tools/reconstruction.hh"
#include "km3net-dataformat/offline/Head.hh"
#include "km3net-dataformat/offline/Evt.hh"
#include "km3net-dataformat/offline/Trk.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
#include "Jeep/JProperties.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JPMTRouter.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JPhysics/JConstants.hh"
#include "JSirene/JVisibleEnergyToolkit.hh"
#include "JAAnet/JHeadToolkit.hh"
#include "JAAnet/JAAnetToolkit.hh"
#include "JAAnet/JEvtCategoryMap.hh"
#include "JSupport/JLimit.hh"
#include "JSupport/JSupport.hh"
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JEvtWeightFileScannerSet.hh"
#include "JSupport/JMonteCarloFileSupportkit.hh"
#include "JTools/JRange.hh"

Go to the source code of this file.

Functions

int main (int argc, char **argv)
 

Detailed Description

Example program to analyse fit results from Evt formatted data.

Author
bofearraigh, bjjung

Definition in file JAAPostfit.cc.

Function Documentation

◆ main()

int main ( int argc,
char ** argv )

Definition at line 58 of file JAAPostfit.cc.

59{
60 using namespace std;
61 using namespace JPP;
62 using namespace KM3NETDAQ;
63
64 JMultipleFileScanner_t inputFiles;
65
66 JLimit_t numberOfEvents;
67
68 string detectorFile;
69 string outputFile;
70
71 JRange<double> energyRange;
72 JRange<double> coszenithRange;
73 JCylinder3D containmentVolume = getMaximumContainmentVolume();
74 JRange<unsigned int> triggeredHitsRange;
75
76 bool useWeights;
77 string oscProbTable;
78 JOscParameters<double> oscParameters;
79 JFluxMap fluxMap{JOscProbInterpolator<>()};
80
81 string recotype;
82
83 int debug;
84
85 try {
86
87 JProperties cuts;
88
89 cuts.insert(gmake_property(energyRange));
90 cuts.insert(gmake_property(coszenithRange));
91 cuts.insert(gmake_property(containmentVolume));
92 cuts.insert(gmake_property(triggeredHitsRange));
93
94 JParser<> zap("Example program to analyse track fit results from Evt formatted data.");
95
96 zap['f'] = make_field(inputFiles);
97 zap['n'] = make_field(numberOfEvents)
98 = JLimit::max();
99 zap['a'] = make_field(detectorFile)
100 = "";
101 zap['o'] = make_field(outputFile)
102 = "histogram.root";
103 zap['R'] = make_field(recotype) =
104 track_t(),
105 shower_t();
106 zap['W'] = make_field(useWeights);
107 zap['@'] = make_field(fluxMap)
108 = JPARSER::initialised();
109 zap['P'] = make_field(oscProbTable)
110 = JPARSER::initialised();
111 zap['#'] = make_field(oscParameters)
112 = JPARSER::initialised();
113 zap['%'] = make_field(cuts)
114 = JPARSER::initialised();
115 zap['d'] = make_field(debug)
116 = 2;
117
118 zap(argc, argv);
119 }
120 catch(const exception& error) {
121 FATAL(error.what() << endl);
122 }
123
124 if (useWeights && numberOfEvents != JLimit::max()) {
125 FATAL("Cannot apply weighting to limited number of events.");
126 }
127
128
129 // Load detector file
130
131 JDetector detector;
132
133 if (detectorFile != "") {
134 try {
135 load(detectorFile, detector);
136 }
137 catch(const JException& error) {
138 FATAL(error);
139 }
140 }
141
142 const JPMTRouter router(detector);
143
144
145 // Load oscillation probability table
146
147 if (!oscProbTable.empty()) {
148
149 JOscProbInterpolator<>& interpolator = static_cast<JOscProbInterpolator<>&>(fluxMap.oscProb.getOscProbInterface());
150
151 interpolator.load(oscProbTable.c_str());
152 interpolator.set (oscParameters);
153 }
154
155
156 // Create file scanners
157
158 JEvtWeightFileScannerSet<> scanners(inputFiles);
159
160 if (scanners.setFlux(fluxMap) == 0) {
161 WARNING("No flux functions set." << endl);
162 }
163
164
165 // Set reco getter function pointer
166
167 bool (*has_reconstruction)(const Evt&);
168 const Trk& (*get_best_fit)(const Evt&);
169
170 if (recotype == track_t()) {
171 has_reconstruction = &has_reconstructed_jppmuon;
172 get_best_fit = &get_best_reconstructed_jppmuon;
173 } else if (recotype == shower_t()) {
174 has_reconstruction = &has_reconstructed_jppshower;
175 get_best_fit = &get_best_reconstructed_jppshower;
176 } else {
177 FATAL("Invalid recotype: " << recotype);
178 }
179
180
181 // Histogram bins
182
183 const Int_t N_Nhits = 400;
184 const Double_t Nhits_min = 0.0 - 0.5;
185 const Double_t Nhits_max = 1000.0 - 0.5;
186
187 const Int_t N_Noverlays = 40;
188 const Double_t Noverlays_min = 0.0 - 0.5;
189 const Double_t Noverlays_max = 1000.0 - 0.5;
190
191 const Int_t N_radius = 201;
192 const Double_t radius_min = -700.0 - 0.5;
193 const Double_t radius_max = +700.0 + 0.5;
194
195 const Int_t N_height = 18;
196 const Double_t height_min = 20.0;
197 const Double_t height_max = 200.0;
198
199 const Int_t N_ToTfrac = 20;
200 const Double_t ToTfrac_min = 0.0;
201 const Double_t ToTfrac_max = 1.0;
202
203 const Int_t N_dt = 1500;
204 const Double_t dt_min = -500.0 - 0.5;
205 const Double_t dt_max = 1000.0 - 0.5;
206
207 const Int_t N_zenith = 21;
208 const Double_t zenith_min = -1.0;
209 const Double_t zenith_max = +1.0;
210
211 const Int_t N_energy = 60;
212 const Double_t energy_min = -2.0;
213 const Double_t energy_max = 4.0;
214
215 const Int_t N_quality = 100;
216 const Double_t quality_min = -200.0;
217 const Double_t quality_max = +800.0;
218
219 const Int_t N_beta0 = 60;
220 const Double_t beta0_min = -2.0;
221 const Double_t beta0_max = +3.5;
222
223
224 // Set of histograms for data/MC-comparisons
225
226 string ordinateLabel(useWeights ? "Rate [Hz]" : "Number of events");
227
228 TFile out(outputFile.c_str(), "recreate");
229
230 TH1D hs ("hs", MAKE_CSTRING("; #snapshot hits; " << ordinateLabel),
231 N_Nhits, Nhits_min, Nhits_max);
232 TH1D hn ("hn", MAKE_CSTRING("; #triggered hits; " << ordinateLabel),
233 N_Nhits, Nhits_min, Nhits_max);
234 TH1D ho ("ho", MAKE_CSTRING("; #overlays; " << ordinateLabel),
235 N_Noverlays, Noverlays_min, Noverlays_max);
236
237 TH1D hz ("hz", MAKE_CSTRING("; cos(#theta); " << ordinateLabel),
238 N_zenith, zenith_min, zenith_max);
239 TH1D he ("he", MAKE_CSTRING("; log_{10}(E_{reco} / 1 GeV); " << ordinateLabel),
240 N_energy, energy_min, energy_max);
241 TH1D ht ("ht", MAKE_CSTRING("; #Delta t [ns]; " << ordinateLabel),
242 N_dt, dt_min, dt_max);
243
244 TH1D hZ ("hZ", MAKE_CSTRING("; longitudinal ToT-barycenter [m]; " << ordinateLabel),
245 N_height, height_min, height_max);
246 TH1D hT0 ("hT0", MAKE_CSTRING("; ToT-fraction below earliest triggered hit [ns]; " << ordinateLabel),
247 N_ToTfrac, ToTfrac_min, ToTfrac_max);
248 TH1D hT1 ("hT1", MAKE_CSTRING("; ToT-fraction above earliest triggered hit [ns]; " << ordinateLabel),
249 N_ToTfrac, ToTfrac_min, ToTfrac_max);
250
251 TH1D hq ("hq", MAKE_CSTRING("; quality; " << ordinateLabel),
252 N_quality, quality_min, quality_max);
253 TH1D hb0 ("hb0", MAKE_CSTRING("; #beta_{0}; " << ordinateLabel),
254 N_beta0, beta0_min, beta0_max);
255
256 TH2D hzo ("hzo", MAKE_CSTRING("; cos(#theta); #overlays; " << ordinateLabel),
257 N_zenith, zenith_min, zenith_max,
258 N_Noverlays, Noverlays_min, Noverlays_max);
259 TH2D hze ("hze", MAKE_CSTRING("; cos(#theta); log_{10}(E_{reco} / 1 GeV); " << ordinateLabel),
260 N_zenith, zenith_min, zenith_max,
261 N_energy, energy_min, energy_max);
262 TH2D heo ("heo", MAKE_CSTRING("; log_{10}(E_{reco} / 1 GeV); #overlays; " << ordinateLabel),
263 N_energy, energy_min, energy_max,
264 N_Noverlays, Noverlays_min, Noverlays_max);
265 TH2D hzt ("hzt", MAKE_CSTRING("; cos(#theta); #Delta t [ns]; " << ordinateLabel),
266 N_zenith, zenith_min, zenith_max,
267 N_dt, dt_min, dt_max);
268 TH2D het ("het", MAKE_CSTRING("; log_{10}(E_{reco} / 1 GeV); #Delta t [ns]; " << ordinateLabel),
269 N_energy, energy_min, energy_max,
270 N_dt, dt_min, dt_max);
271 TH2D hot ("hot", MAKE_CSTRING("; #overlays; #Delta t [ns]; " << ordinateLabel),
272 N_Noverlays, Noverlays_min, Noverlays_max,
273 N_dt, dt_min, dt_max);
274
275 TH2D hzq ("hzq", MAKE_CSTRING("; cos(#theta); quality; " << ordinateLabel),
276 N_zenith, zenith_min, zenith_max,
277 N_quality, quality_min, quality_max);
278 TH2D heq ("heq", MAKE_CSTRING("; log_{10}(E_{reco} / 1 GeV); quality; " << ordinateLabel),
279 N_energy, energy_min, energy_max,
280 N_quality, quality_min, quality_max);
281 TH2D hoq ("hoq", MAKE_CSTRING("; #overlays; quality; " << ordinateLabel),
282 N_Noverlays, Noverlays_min, Noverlays_max,
283 N_quality, quality_min, quality_max);
284
285 TH2D hob0("hob0", MAKE_CSTRING("; #overlays; log_{10}(#beta_{0}); " << ordinateLabel),
286 N_Noverlays, Noverlays_min, Noverlays_max,
287 N_beta0, beta0_min, beta0_max);
288 TH2D heb0("heb0", MAKE_CSTRING("; log_{10}(E_{reco} / 1 GeV); log_{10}(#beta_{0}); " << ordinateLabel),
289 N_energy, energy_min, energy_max,
290 N_beta0, beta0_min, beta0_max);
291 TH2D hzb0("hzb0", MAKE_CSTRING("; cos(#theta); log_{10}(#beta_{0}); " << ordinateLabel),
292 N_zenith, zenith_min, zenith_max,
293 N_beta0, beta0_min, beta0_max);
294
295 TH2D hxy ("hxy", MAKE_CSTRING("; x [m]; y [m]; " << ordinateLabel),
296 N_radius, radius_min, radius_max,
297 N_radius, radius_min, radius_max);
298
299
300 // Loop over events
301
302 STATUS(RIGHT(15) << "header ID" << RIGHT(15) << "event" << RIGHT(15) << "weight" << endl);
303
304 for (JEvtWeightFileScannerSet<>::iterator scanner = scanners.begin(); scanner != scanners.end(); ++scanner) {
305
306 scanner->setLimit(numberOfEvents);
307
308 while (scanner->hasNext()) {
309
310 const Evt* evt = scanner->next();
311
312 const double weight = (useWeights ? scanner->getWeight(*evt) : 1.0);
313
314 const size_t NsnapshotHits = evt->hits.size();
315 const size_t NtriggeredHits = count_if(evt->hits.begin(), evt->hits.end(),
316 make_predicate(&Hit::trig, (long long unsigned int) 0, JComparison::gt()));
317
318 if (!triggeredHitsRange(NtriggeredHits)) { continue; }
319
320 hs .Fill(NsnapshotHits, weight);
321 hn .Fill(NtriggeredHits, weight);
322 ho .Fill(evt->overlays, weight);
323
324 if (!has_reconstruction(*evt)) { continue; }
325
326 const Trk& best_trk = (*get_best_fit)(*evt);
327 const JTrack3E& tb = getTrack(best_trk);
328
329 const double logE = log10(tb.getE());
330
331 if (!energyRange (tb.getE()) ||
332 !coszenithRange(tb.getDZ()) ||
333 !containmentVolume.is_inside(tb.getPosition())) { // Apply cuts
334 continue;
335 }
336
337 // Extract event-weight
338
339 STATUS(RIGHT (15) << distance(scanners.begin(), scanner) <<
340 RIGHT (15) << scanner->getCounter() <<
341 SCIENTIFIC(15, 3) << weight << '\r'); DEBUG(endl);
342
343 // Fill histograms
344
345 hz .Fill(tb.getDZ(), weight);
346 he .Fill(logE, weight);
347 hq .Fill(best_trk.lik, weight);
348
349 hzo.Fill(tb.getDZ(), evt->overlays, weight);
350 hze.Fill(tb.getDZ(), logE, weight);
351 heo.Fill(logE, evt->overlays, weight);
352
353 hoq.Fill(evt->overlays, best_trk.lik, weight);
354 hzq.Fill(tb.getDZ(), best_trk.lik, weight);
355 heq.Fill(logE, best_trk.lik, weight);
356
357 hxy.Fill(tb.getX(), tb.getY(), weight);
358
359 if (best_trk.fitinf[JGANDALF_BETA0_RAD]) {
360
361 const double logb0 = log10(180.0 / M_PI * best_trk.fitinf[JGANDALF_BETA0_RAD]);
362
363 hb0 .Fill(logb0, weight);
364 hob0.Fill(evt->overlays, logb0, weight);
365 heb0.Fill(logE, logb0, weight);
366 hzb0.Fill(tb.getDZ(), logb0, weight);
367 }
368
369 double ToTbelow = 0.0;
370 double ToTabove = 0.0;
371 double ToTtotal = 0.0;
372 double ToTweightedZ = 0.0;
373
374 vector<Hit>::const_iterator firstHit = min_element(evt->hits.cbegin(), evt->hits.cend(),
375 make_comparator(&Hit::t));
376
377 for (vector<Hit>::const_iterator hit = evt->hits.cbegin(); hit != evt->hits.cend(); ++hit) {
378
379 if (router.hasPMT(hit->pmt_id)) {
380
381 const JPMT& pmt = router.getPMT(hit->pmt_id);
382
383 const double t0 = tb.getT(pmt);
384 const double t1 = getTime(*hit);
385
386 ht .Fill(t1 - t0, weight);
387 hot.Fill(evt->overlays, t1 - t0, weight);
388 hzt.Fill(tb.getDZ(), t1 - t0, weight);
389 het.Fill(log10(best_trk.E), t1 - t0, weight);
390
391 if (hit->trig > 0) {
392
393 ToTtotal += hit->tot;
394 ToTweightedZ += hit->tot * hit->pos.z;
395
396 if (hit->pos.z < firstHit->pos.z) {
397 ToTbelow += hit->tot;
398 } else if (hit->pos.z >= firstHit->pos.z) {
399 ToTabove += hit->tot;
400 }
401 }
402 }
403 }
404
405 hT0.Fill(ToTbelow / ToTtotal, weight);
406 hT1.Fill(ToTabove / ToTtotal, weight);
407 hZ .Fill(ToTweightedZ / ToTtotal, weight);
408 }
409 }
410
411 STATUS(endl);
412
413 out.Write();
414 out.Close();
415
416 return 0;
417}
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
FATAL
#define FATAL(A)
Definition JMessage.hh:67
debug
int debug
debug level
Definition JSirene.cc:69
WARNING
#define WARNING(A)
Definition JMessage.hh:65
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2142
MAKE_CSTRING
#define MAKE_CSTRING(A)
Make C-string.
Definition JPrint.hh:72
gmake_property
#define gmake_property(A)
macros to convert (template) parameter to JPropertiesElement object
Definition JProperties.hh:1293
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 PhysicsEvent.hh:434
JAANET::JEvtCategoryMap< JFluxHelper >
Template specialisation for a map between event categories and flux factors.
Definition JEvtCategoryMap.hh:515
JDETECTOR::JDetector
Detector data structure.
Definition JDetector.hh:96
JDETECTOR::JPMTRouter
Router for direct addressing of PMT data in detector data structure.
Definition JPMTRouter.hh:37
JDETECTOR::JPMT
Data structure for PMT geometry, calibration and status.
Definition JPMT.hh:49
JEEP::JProperties
Utility class to parse parameter values.
Definition JProperties.hh:501
JGEOMETRY3D::JCylinder3D::is_inside
bool is_inside(const JVector3D &pos) const
Check whether given point is inside cylinder.
Definition JCylinder3D.hh:208
JGEOMETRY3D::JPosition3D::getPosition
const JPosition3D & getPosition() const
Get position.
Definition JPosition3D.hh:130
JGEOMETRY3D::JTrack3D::getT
double getT(const JVector3D &pos) const
Get arrival time of Cherenkov light at given position.
Definition JTrack3D.hh:126
JGEOMETRY3D::JTrack3E
3D track with energy.
Definition JTrack3E.hh:32
JGEOMETRY3D::JTrack3E::getE
double getE() const
Get energy.
Definition JTrack3E.hh:93
JGEOMETRY3D::JVector3D::getY
double getY() const
Get y position.
Definition JVector3D.hh:104
JGEOMETRY3D::JVector3D::getX
double getX() const
Get x position.
Definition JVector3D.hh:94
JGEOMETRY3D::JVersor3D::getDZ
double getDZ() const
Get z direction.
Definition JVersor3D.hh:117
JLANG::JException
General exception.
Definition JException.hh:24
JOSCPROB::JOscParameters
Data structure for single set of oscillation parameters.
Definition JOscParameters.hh:52
JOSCPROB::JOscProbInterpolator
Template definition of a multi-dimensional oscillation probability interpolation table.
Definition JOscProbInterpolator.hh:65
JOSCPROB::JOscProbInterpolator::load
void load(const char *const fileName)
Load oscillation probability table from file.
Definition JOscProbInterpolator.hh:147
JPARSER::JParser
Utility class to parse command line options.
Definition JParser.hh:1698
JTOOLS::JRange
Range of values.
Definition JRange.hh:42
JGANDALF_BETA0_RAD
static const int JGANDALF_BETA0_RAD
KM3NeT Data Definitions v3.4.0-8-ge14cb17 https://git.km3net.de/common/km3net-dataformat.
Definition fitparameters.hh:12
JAANET::getTrack
JTrack3E getTrack(const Trk &track)
Get track.
Definition JAAnetToolkit.hh:262
JDETECTOR::load
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
Definition JDetectorToolkit.hh:505
JLANG::make_comparator
JComparator< JResult_t T::*, JComparison::lt > make_comparator(JResult_t T::*member)
Helper method to create comparator between values of data member.
Definition JLang/JComparator.hh:185
JLANG::RIGHT
@ RIGHT
Definition JTwosome.hh:18
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition JAAnetToolkit.hh:43
JSIRENE::getMaximumContainmentVolume
const JCylinder3D getMaximumContainmentVolume()
Auxiliary function to retrieve the maximum cylindrical containment volume.
Definition JVisibleEnergyToolkit.hh:46
JSYSTEM::getTime
const char * getTime()
Get current local time conform ISO-8601 standard.
Definition JDateAndTime.hh:55
KM3NETDAQ
KM3NeT DAQ data structures and auxiliaries.
Definition DataQueue.cc:39
mc_keys::track_t
const char *const track_t
Definition io_ascii.hh:27
Evt
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
Definition Evt.hh:21
Evt::hits
std::vector< Hit > hits
list of hits
Definition Evt.hh:38
Evt::overlays
unsigned int overlays
number of overlaying triggered events
Definition Evt.hh:32
Hit::trig
ULong64_t trig
non-zero if the hit is a trigger hit.
Definition Hit.hh:18
Hit::t
double t
hit time (from tdc+calibration or MC truth)
Definition Hit.hh:23
JAANET::detector
Detector file.
Definition JHead.hh:227
JPARSER::initialised
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition JParser.hh:68
JSUPPORT::JEvtWeightFileScannerSet
Auxiliary class for organising Monte Carlo file scanners associated with event weighters.
Definition JEvtWeightFileScannerSet.hh:61
JSUPPORT::JEvtWeightFileScannerSet::iterator
std::vector< filescanner_type >::iterator iterator
Definition JEvtWeightFileScannerSet.hh:72
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
JSUPPORT::JMultipleFileScanner_t
Auxiliary base class for list of file names.
Definition JMultipleFileScanner.hh:51
RIGHT
Auxiliary data structure for alignment of data.
Definition JManip.hh:298
SCIENTIFIC
Auxiliary data structure for floating point format specification.
Definition JManip.hh:488
Trk
The Trk class represents a Monte Carlo (MC) particle as well as a reconstructed track/shower.
Definition Trk.hh:15
Trk::fitinf
std::vector< double > fitinf
place to store additional fit info, see km3net-dataformat/definitions/fitparameters....
Definition Trk.hh:32
Trk::E
double E
Energy [GeV] (either MC truth or reconstructed)
Definition Trk.hh:20
Trk::lik
double lik
likelihood or lambda value (for aafit, lambda)
Definition Trk.hh:23
get_best_reconstructed_jppshower
const Trk & get_best_reconstructed_jppshower(const Evt &evt)
Get best reconstructed shower.
Definition tools/reconstruction.hh:399
get_best_reconstructed_jppmuon
const Trk & get_best_reconstructed_jppmuon(const Evt &evt)
Get best reconstructed muon.
Definition tools/reconstruction.hh:387
has_reconstructed_jppshower
bool has_reconstructed_jppshower(const Evt &evt)
Test whether given event has a track with shower reconstruction.
Definition tools/reconstruction.hh:317
has_reconstructed_jppmuon
bool has_reconstructed_jppmuon(const Evt &evt)
Test whether given event has a track with muon reconstruction.
Definition tools/reconstruction.hh:305
Generated by doxygen 1.12.0