Jpp 20.0.0-rc.2
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 JEvtWeightFactorMap weightFactorMap;
78
79 string recotype;
80
81 int debug;
82
83 try {
84
85 JProperties cuts;
86
87 cuts.insert(gmake_property(energyRange));
88 cuts.insert(gmake_property(coszenithRange));
89 cuts.insert(gmake_property(containmentVolume));
90 cuts.insert(gmake_property(triggeredHitsRange));
91
92 JParser<> zap("Example program to analyse track fit results from Evt formatted data.");
93
94 zap['f'] = make_field(inputFiles);
95 zap['n'] = make_field(numberOfEvents)
96 = JLimit::max();
97 zap['a'] = make_field(detectorFile)
98 = "";
99 zap['o'] = make_field(outputFile)
100 = "histogram.root";
101 zap['R'] = make_field(recotype) =
102 track_t(),
103 shower_t();
104 zap['W'] = make_field(useWeights);
105 zap['@'] = make_field(weightFactorMap)
106 = JPARSER::initialised();
107 zap['%'] = make_field(cuts)
108 = JPARSER::initialised();
109 zap['d'] = make_field(debug)
110 = 2;
111
112 zap(argc, argv);
113 }
114 catch(const exception& error) {
115 FATAL(error.what() << endl);
116 }
117
118 if (useWeights && numberOfEvents != JLimit::max()) {
119 FATAL("Cannot apply weighting to limited number of events.");
120 }
121
122
123 // Load detector file
124
125 JDetector detector;
126
127 if (detectorFile != "") {
128 try {
129 load(detectorFile, detector);
130 }
131 catch(const JException& error) {
132 FATAL(error);
133 }
134 }
135
136 const JPMTRouter router(detector);
137
138
139 // Set reco getter function pointer
140
141 bool (*has_reconstruction)(const Evt&);
142 const Trk& (*get_best_fit)(const Evt&);
143
144 if (recotype == track_t()) {
145 has_reconstruction = &has_reconstructed_jppmuon;
146 get_best_fit = &get_best_reconstructed_jppmuon;
147 } else if (recotype == shower_t()) {
148 has_reconstruction = &has_reconstructed_jppshower;
149 get_best_fit = &get_best_reconstructed_jppshower;
150 } else {
151 FATAL("Invalid recotype: " << recotype);
152 }
153
154
155 // Histogram bins
156
157 const Int_t N_Nhits = 400;
158 const Double_t Nhits_min = 0.0 - 0.5;
159 const Double_t Nhits_max = 1000.0 - 0.5;
160
161 const Int_t N_Noverlays = 40;
162 const Double_t Noverlays_min = 0.0 - 0.5;
163 const Double_t Noverlays_max = 1000.0 - 0.5;
164
165 const Int_t N_radius = 201;
166 const Double_t radius_min = -700.0 - 0.5;
167 const Double_t radius_max = +700.0 + 0.5;
168
169 const Int_t N_height = 18;
170 const Double_t height_min = 20.0;
171 const Double_t height_max = 200.0;
172
173 const Int_t N_ToTfrac = 20;
174 const Double_t ToTfrac_min = 0.0;
175 const Double_t ToTfrac_max = 1.0;
176
177 const Int_t N_dt = 1500;
178 const Double_t dt_min = -500.0 - 0.5;
179 const Double_t dt_max = 1000.0 - 0.5;
180
181 const Int_t N_zenith = 21;
182 const Double_t zenith_min = -1.0;
183 const Double_t zenith_max = +1.0;
184
185 const Int_t N_energy = 60;
186 const Double_t energy_min = -2.0;
187 const Double_t energy_max = 4.0;
188
189 const Int_t N_LoE = 100;
190 const Double_t LoE_min = -1.0;
191 const Double_t LoE_max = 8.0;
192
193 const Int_t N_quality = 100;
194 const Double_t quality_min = -200.0;
195 const Double_t quality_max = +800.0;
196
197 const Int_t N_beta0 = 60;
198 const Double_t beta0_min = -2.0;
199 const Double_t beta0_max = +3.5;
200
201
202 // Set of histograms for data/MC-comparisons
203
204 string ordinateLabel(useWeights ? "Rate [Hz]" : "Number of events");
205
206 TFile out(outputFile.c_str(), "recreate");
207
208 TH1D hs ("hs", MAKE_CSTRING("; #snapshot hits; " << ordinateLabel),
209 N_Nhits, Nhits_min, Nhits_max);
210 TH1D hn ("hn", MAKE_CSTRING("; #triggered hits; " << ordinateLabel),
211 N_Nhits, Nhits_min, Nhits_max);
212 TH1D ho ("ho", MAKE_CSTRING("; #overlays; " << ordinateLabel),
213 N_Noverlays, Noverlays_min, Noverlays_max);
214
215 TH1D hz ("hz", MAKE_CSTRING("; cos(#theta); " << ordinateLabel),
216 N_zenith, zenith_min, zenith_max);
217 TH1D he ("he", MAKE_CSTRING("; log_{10}(E_{reco} / 1 GeV); " << ordinateLabel),
218 N_energy, energy_min, energy_max);
219 TH1D hLoE("hLoE", MAKE_CSTRING("; L/E [km/GeV]; " << ordinateLabel),
220 N_LoE, LoE_min, LoE_max);
221 TH1D ht ("ht", MAKE_CSTRING("; #Delta t [ns]; " << ordinateLabel),
222 N_dt, dt_min, dt_max);
223
224 TH1D hZ ("hZ", MAKE_CSTRING("; longitudinal ToT-barycenter [m]; " << ordinateLabel),
225 N_height, height_min, height_max);
226 TH1D hT0 ("hT0", MAKE_CSTRING("; ToT-fraction below earliest triggered hit [ns]; " << ordinateLabel),
227 N_ToTfrac, ToTfrac_min, ToTfrac_max);
228 TH1D hT1 ("hT1", MAKE_CSTRING("; ToT-fraction above earliest triggered hit [ns]; " << ordinateLabel),
229 N_ToTfrac, ToTfrac_min, ToTfrac_max);
230
231 TH1D hq ("hq", MAKE_CSTRING("; quality; " << ordinateLabel),
232 N_quality, quality_min, quality_max);
233 TH1D hb0 ("hb0", MAKE_CSTRING("; #beta_{0}; " << ordinateLabel),
234 N_beta0, beta0_min, beta0_max);
235
236 TH2D hzo ("hzo", MAKE_CSTRING("; cos(#theta); #overlays; " << ordinateLabel),
237 N_zenith, zenith_min, zenith_max,
238 N_Noverlays, Noverlays_min, Noverlays_max);
239 TH2D hze ("hze", MAKE_CSTRING("; cos(#theta); log_{10}(E_{reco} / 1 GeV); " << ordinateLabel),
240 N_zenith, zenith_min, zenith_max,
241 N_energy, energy_min, energy_max);
242 TH2D heo ("heo", MAKE_CSTRING("; log_{10}(E_{reco} / 1 GeV); #overlays; " << ordinateLabel),
243 N_energy, energy_min, energy_max,
244 N_Noverlays, Noverlays_min, Noverlays_max);
245 TH2D hzt ("hzt", MAKE_CSTRING("; cos(#theta); #Delta t [ns]; " << ordinateLabel),
246 N_zenith, zenith_min, zenith_max,
247 N_dt, dt_min, dt_max);
248 TH2D het ("het", MAKE_CSTRING("; log_{10}(E_{reco} / 1 GeV); #Delta t [ns]; " << ordinateLabel),
249 N_energy, energy_min, energy_max,
250 N_dt, dt_min, dt_max);
251 TH2D hot ("hot", MAKE_CSTRING("; #overlays; #Delta t [ns]; " << ordinateLabel),
252 N_Noverlays, Noverlays_min, Noverlays_max,
253 N_dt, dt_min, dt_max);
254
255 TH2D hzq ("hzq", MAKE_CSTRING("; cos(#theta); quality; " << ordinateLabel),
256 N_zenith, zenith_min, zenith_max,
257 N_quality, quality_min, quality_max);
258 TH2D heq ("heq", MAKE_CSTRING("; log_{10}(E_{reco} / 1 GeV); quality; " << ordinateLabel),
259 N_energy, energy_min, energy_max,
260 N_quality, quality_min, quality_max);
261 TH2D hoq ("hoq", MAKE_CSTRING("; #overlays; quality; " << ordinateLabel),
262 N_Noverlays, Noverlays_min, Noverlays_max,
263 N_quality, quality_min, quality_max);
264
265 TH2D hob0("hob0", MAKE_CSTRING("; #overlays; log_{10}(#beta_{0}); " << ordinateLabel),
266 N_Noverlays, Noverlays_min, Noverlays_max,
267 N_beta0, beta0_min, beta0_max);
268 TH2D heb0("heb0", MAKE_CSTRING("; log_{10}(E_{reco} / 1 GeV); log_{10}(#beta_{0}); " << ordinateLabel),
269 N_energy, energy_min, energy_max,
270 N_beta0, beta0_min, beta0_max);
271 TH2D hzb0("hzb0", MAKE_CSTRING("; cos(#theta); log_{10}(#beta_{0}); " << ordinateLabel),
272 N_zenith, zenith_min, zenith_max,
273 N_beta0, beta0_min, beta0_max);
274
275 TH2D hxy ("hxy", MAKE_CSTRING("; x [m]; y [m]; " << ordinateLabel),
276 N_radius, radius_min, radius_max,
277 N_radius, radius_min, radius_max);
278
279
280 // Loop over events
281
282 JEvtWeightFileScannerSet<> scanners(inputFiles);
283
284 if (scanners.setEvtWeightFactor(weightFactorMap) == 0) {
285 WARNING("No flux functions set." << endl);
286 }
287
288 STATUS(RIGHT(15) << "header ID" << RIGHT(15) << "event" << RIGHT(15) << "weight" << endl);
289
290 for (JEvtWeightFileScannerSet<>::iterator scanner = scanners.begin(); scanner != scanners.end(); ++scanner) {
291
292 for (JMultipleFileScanner<Evt> in(scanner->getFilelist(), numberOfEvents); in.hasNext(); ) {
293
294 const Evt* evt = in.next();
295
296 const double weight = (useWeights ? scanner->getWeight(*evt) : 1.0);
297
298 const size_t NsnapshotHits = evt->hits.size();
299 const size_t NtriggeredHits = count_if(evt->hits.begin(), evt->hits.end(),
300 make_predicate(&Hit::trig, (ULong64_t) 0, JComparison::gt()));
301
302 if (!triggeredHitsRange(NtriggeredHits)) { continue; }
303
304 hs.Fill(NsnapshotHits, weight);
305 hn.Fill(NtriggeredHits, weight);
306 ho.Fill(evt->overlays, weight);
307
308 if (!has_reconstruction(*evt)) { continue; }
309
310 const Trk& best_trk = (*get_best_fit)(*evt);
311 const JTrack3E& tb = getTrack(best_trk);
312
313 const double logE = log10(tb.getE());
314 const double L = weightFactorMap.getOscProb().getBaseline(-tb.getDZ());
315 const double logLoE = log10(L/tb.getE());
316
317 if (!energyRange (tb.getE()) ||
318 !coszenithRange(tb.getDZ()) ||
319 !containmentVolume.is_inside(tb.getPosition())) { // Apply cuts
320 continue;
321 }
322
323 // Extract event-weight
324
325 STATUS(RIGHT (15) << distance(scanners.begin(), scanner) <<
326 RIGHT (15) << in.getCounter() <<
327 SCIENTIFIC(15, 3) << weight << '\r'); DEBUG(endl);
328
329 // Fill histograms
330
331 hz .Fill(tb.getDZ(), weight);
332 he .Fill(logE, weight);
333 hLoE.Fill(logLoE, weight);
334 hq .Fill(best_trk.lik, weight);
335
336 hzo.Fill(tb.getDZ(), evt->overlays, weight);
337 hze.Fill(tb.getDZ(), logE, weight);
338 heo.Fill(logE, evt->overlays, weight);
339
340 hoq.Fill(evt->overlays, best_trk.lik, weight);
341 hzq.Fill(tb.getDZ(), best_trk.lik, weight);
342 heq.Fill(logE, best_trk.lik, weight);
343
344 hxy.Fill(tb.getX(), tb.getY(), weight);
345
346 if (best_trk.fitinf[JGANDALF_BETA0_RAD]) {
347
348 const double logb0 = log10(180.0 / M_PI * best_trk.fitinf[JGANDALF_BETA0_RAD]);
349
350 hb0 .Fill(logb0, weight);
351 hob0.Fill(evt->overlays, logb0, weight);
352 heb0.Fill(logE, logb0, weight);
353 hzb0.Fill(tb.getDZ(), logb0, weight);
354 }
355
356 double ToTbelow = 0.0;
357 double ToTabove = 0.0;
358 double ToTtotal = 0.0;
359 double ToTweightedZ = 0.0;
360
361 vector<Hit>::const_iterator firstHit = min_element(evt->hits.cbegin(), evt->hits.cend(),
362 make_comparator(&Hit::t));
363
364 for (vector<Hit>::const_iterator hit = evt->hits.cbegin(); hit != evt->hits.cend(); ++hit) {
365
366 if (router.hasPMT(hit->pmt_id)) {
367
368 const JPMT& pmt = router.getPMT(hit->pmt_id);
369
370 const double t0 = tb.getT(pmt);
371 const double t1 = getTime(*hit);
372
373 ht .Fill(t1 - t0, weight);
374 hot.Fill(evt->overlays, t1 - t0, weight);
375 hzt.Fill(tb.getDZ(), t1 - t0, weight);
376 het.Fill(log10(best_trk.E), t1 - t0, weight);
377
378 if (hit->trig > 0) {
379
380 ToTtotal += hit->tot;
381 ToTweightedZ += hit->tot * hit->pos.z;
382
383 if (hit->pos.z < firstHit->pos.z) {
384 ToTbelow += hit->tot;
385 } else if (hit->pos.z >= firstHit->pos.z) {
386 ToTabove += hit->tot;
387 }
388 }
389 }
390 }
391
392 hT0.Fill(ToTbelow / ToTtotal, weight);
393 hT1.Fill(ToTabove / ToTtotal, weight);
394 hZ .Fill(ToTweightedZ / ToTtotal, weight);
395 }
396 }
397
398 STATUS(endl);
399
400 out.Write();
401 out.Close();
402
403 return 0;
404}
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:72
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:1139
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::JEvtCategoryMap< JEvtWeightFactorProduct >
Template specialisation for a map between event categories and event-weight factor products.
Definition JEvtCategoryMap.hh:177
JAANET::JEvtCategoryMap< JEvtWeightFactorProduct >::getOscProb
const JOscProbHelper & getOscProb() const
Get oscillation probability calculator.
Definition JEvtCategoryMap.hh:368
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:503
JGEOMETRY3D::JCylinder3D::is_inside
bool is_inside(const JVector3D &pos) const
Check whether given point is inside cylinder.
Definition JCylinder3D.hh:208
JGEOMETRY3D::JEnergy::getE
double getE() const
Get energy.
Definition JGeometry3D/JEnergy.hh:65
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:87
JGEOMETRY3D::JTrack3E
3D track with energy.
Definition JTrack3E.hh:34
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
JPARSER::JParser
Utility class to parse command line options.
Definition JParser.hh:1698
JSUPPORT::JMultipleFileScanner
General purpose class for object reading from a list of file names.
Definition JMultipleFileScanner.hh:306
JSUPPORT::JMultipleFileScanner::hasNext
virtual bool hasNext() override
Check availability of next element.
Definition JMultipleFileScanner.hh:409
JTOOLS::JRange
Range of values.
Definition JRange.hh:42
JGANDALF_BETA0_RAD
static const int JGANDALF_BETA0_RAD
KM3NeT Data Definitions v3.6.1-2-g905a24d https://git.km3net.de/common/km3net-dataformat.
Definition fitparameters.hh:12
JAANET::getTrack
JTrack3E getTrack(const Trk &track)
Get track.
Definition JAAnetToolkit.hh:278
JDETECTOR::load
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
Definition JDetectorToolkit.hh:536
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:45
JSIRENE::getMaximumContainmentVolume
const JCylinder3D getMaximumContainmentVolume()
Forward function declarations.
Definition JVisibleEnergyToolkit.hh:48
JSYSTEM::getTime
const char * getTime()
Get current local time conform ISO-8601 standard.
Definition JDateAndTime.hh:56
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
JOSCPROB::JOscProbHelper::getBaseline
double getBaseline(const double costh) const
Get baseline for a given cosine zenith angle.
Definition JOscProbHelper.hh:180
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:58
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