Jpp test-rotations-old
the software that should make you happy
Loading...
Searching...
No Matches
JAAPostfit.cc
Go to the documentation of this file.
1
2#include <string>
3#include <iostream>
4#include <iomanip>
5#include <cmath>
6
7#include "TROOT.h"
8#include "TFile.h"
9#include "TH1D.h"
10
12
14
18
19#include "Jeep/JParser.hh"
20#include "Jeep/JMessage.hh"
21#include "Jeep/JProperties.hh"
22
26
28
30
34
35#include "JSupport/JLimit.hh"
36#include "JSupport/JSupport.hh"
40
41#include "JTools/JRange.hh"
42
43
44namespace {
45
46 inline const char* const track_t() { return "track"; }
47 inline const char* const shower_t() { return "shower"; }
48}
49
50
51/**
52 * \file
53 *
54 * Example program to analyse fit results from Evt formatted data.
55 *
56 * \author bofearraigh, bjjung
57 */
58int main(int argc, char **argv)
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 JFluxMap fluxMap;
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(fluxMap)
107 zap['%'] = make_field(cuts)
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
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.setFlux(fluxMap) == 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 scanner->setLimit(numberOfEvents);
293
294 while (scanner->hasNext()) {
295
296 const Evt* evt = scanner->next();
297
298 const double weight = (useWeights ? scanner->getWeight(*evt) : 1.0);
299
300 const size_t NsnapshotHits = evt->hits.size();
301 const size_t NtriggeredHits = count_if(evt->hits.begin(), evt->hits.end(),
302 make_predicate(&Hit::trig, (long long unsigned int) 0, JComparison::gt()));
303
304 if (!triggeredHitsRange(NtriggeredHits)) { continue; }
305
306 hs.Fill(NsnapshotHits, weight);
307 hn.Fill(NtriggeredHits, weight);
308 ho.Fill(evt->overlays, weight);
309
310 if (!has_reconstruction(*evt)) { continue; }
311
312 const Trk& best_trk = (*get_best_fit)(*evt);
313 const JTrack3E& tb = getTrack(best_trk);
314
315 const double logE = log10(tb.getE());
316 const double L = fluxMap.getOscProb().getBaseline(-tb.getDZ());
317 const double logLoE = log10(L/tb.getE());
318
319 if (!energyRange (tb.getE()) ||
320 !coszenithRange(tb.getDZ()) ||
321 !containmentVolume.is_inside(tb.getPosition())) { // Apply cuts
322 continue;
323 }
324
325 // Extract event-weight
326
327 STATUS(RIGHT (15) << distance(scanners.begin(), scanner) <<
328 RIGHT (15) << scanner->getCounter() <<
329 SCIENTIFIC(15, 3) << weight << '\r'); DEBUG(endl);
330
331 // Fill histograms
332
333 hz .Fill(tb.getDZ(), weight);
334 he .Fill(logE, weight);
335 hLoE.Fill(logLoE, weight);
336 hq .Fill(best_trk.lik, weight);
337
338 hzo.Fill(tb.getDZ(), evt->overlays, weight);
339 hze.Fill(tb.getDZ(), logE, weight);
340 heo.Fill(logE, evt->overlays, weight);
341
342 hoq.Fill(evt->overlays, best_trk.lik, weight);
343 hzq.Fill(tb.getDZ(), best_trk.lik, weight);
344 heq.Fill(logE, best_trk.lik, weight);
345
346 hxy.Fill(tb.getX(), tb.getY(), weight);
347
348 if (best_trk.fitinf[JGANDALF_BETA0_RAD]) {
349
350 const double logb0 = log10(180.0 / M_PI * best_trk.fitinf[JGANDALF_BETA0_RAD]);
351
352 hb0 .Fill(logb0, weight);
353 hob0.Fill(evt->overlays, logb0, weight);
354 heb0.Fill(logE, logb0, weight);
355 hzb0.Fill(tb.getDZ(), logb0, weight);
356 }
357
358 double ToTbelow = 0.0;
359 double ToTabove = 0.0;
360 double ToTtotal = 0.0;
361 double ToTweightedZ = 0.0;
362
363 vector<Hit>::const_iterator firstHit = min_element(evt->hits.cbegin(), evt->hits.cend(),
365
366 for (vector<Hit>::const_iterator hit = evt->hits.cbegin(); hit != evt->hits.cend(); ++hit) {
367
368 if (router.hasPMT(hit->pmt_id)) {
369
370 const JPMT& pmt = router.getPMT(hit->pmt_id);
371
372 const double t0 = tb.getT(pmt);
373 const double t1 = getTime(*hit);
374
375 ht .Fill(t1 - t0, weight);
376 hot.Fill(evt->overlays, t1 - t0, weight);
377 hzt.Fill(tb.getDZ(), t1 - t0, weight);
378 het.Fill(log10(best_trk.E), t1 - t0, weight);
379
380 if (hit->trig > 0) {
381
382 ToTtotal += hit->tot;
383 ToTweightedZ += hit->tot * hit->pos.z;
384
385 if (hit->pos.z < firstHit->pos.z) {
386 ToTbelow += hit->tot;
387 } else if (hit->pos.z >= firstHit->pos.z) {
388 ToTabove += hit->tot;
389 }
390 }
391 }
392 }
393
394 hT0.Fill(ToTbelow / ToTtotal, weight);
395 hT1.Fill(ToTabove / ToTtotal, weight);
396 hZ .Fill(ToTweightedZ / ToTtotal, weight);
397 }
398 }
399
400 STATUS(endl);
401
402 out.Write();
403 out.Close();
404
405 return 0;
406}
int main(int argc, char **argv)
Definition JAAPostfit.cc:58
Definition of hit and track types and auxiliary methods for handling Monte Carlo data.
string outputFile
Data structure for detector geometry and calibration.
Auxiliaries for defining the range of iterations of objects.
General purpose messaging.
#define DEBUG(A)
Message macros.
Definition JMessage.hh:62
#define STATUS(A)
Definition JMessage.hh:63
#define FATAL(A)
Definition JMessage.hh:67
int debug
debug level
Definition JSirene.cc:72
#define WARNING(A)
Definition JMessage.hh:65
Scanning of objects from multiple files according a format that follows from the extension of each fi...
Direct access to PMT in detector data structure.
Utility class to parse command line options.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2142
Physics constants.
#define MAKE_CSTRING(A)
Make C-string.
Definition JPrint.hh:72
Utility class to parse parameter values.
#define gmake_property(A)
macros to convert (template) parameter to JPropertiesElement object
Auxiliary class to define a range between two values.
ROOT TTree parameter settings of various packages.
Auxiliary methods for evaluating visible energies.
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
Template specialisation for a map between event categories and flux factors.
const JOscProbHelper & getOscProb() const
Get oscillation probability calculator.
Detector data structure.
Definition JDetector.hh:96
Router for direct addressing of PMT data in detector data structure.
Definition JPMTRouter.hh:37
bool hasPMT(const JObjectID &id) const
Has PMT.
const JPMT & getPMT(const JPMTAddress &address) const
Get PMT.
Definition JPMTRouter.hh:92
Data structure for PMT geometry, calibration and status.
Definition JPMT.hh:49
Utility class to parse parameter values.
bool is_inside(const JVector3D &pos) const
Check whether given point is inside cylinder.
const JPosition3D & getPosition() const
Get position.
double getT(const JVector3D &pos) const
Get arrival time of Cherenkov light at given position.
Definition JTrack3D.hh:126
3D track with energy.
Definition JTrack3E.hh:32
double getE() const
Get energy.
Definition JTrack3E.hh:93
double getY() const
Get y position.
Definition JVector3D.hh:104
double getX() const
Get x position.
Definition JVector3D.hh:94
double getDZ() const
Get z direction.
Definition JVersor3D.hh:117
General exception.
Definition JException.hh:24
Utility class to parse command line options.
Definition JParser.hh:1698
Range of values.
Definition JRange.hh:42
static const int JGANDALF_BETA0_RAD
KM3NeT Data Definitions v3.6.0 https://git.km3net.de/common/km3net-dataformat.
JTrack3E getTrack(const Trk &track)
Get track.
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
JComparator< JResult_t T::*, JComparison::lt > make_comparator(JResult_t T::*member)
Helper method to create comparator between values of data member.
@ RIGHT
Definition JTwosome.hh:18
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
const JCylinder3D getMaximumContainmentVolume()
Forward function declarations.
const char * getTime()
Get current local time conform ISO-8601 standard.
KM3NeT DAQ data structures and auxiliaries.
Definition DataQueue.cc:39
const char *const track_t
Definition io_ascii.hh:27
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
Definition Evt.hh:21
std::vector< Hit > hits
list of hits
Definition Evt.hh:38
unsigned int overlays
number of overlaying triggered events
Definition Evt.hh:32
ULong64_t trig
non-zero if the hit is a trigger hit.
Definition Hit.hh:18
double t
hit time (from tdc+calibration or MC truth)
Definition Hit.hh:23
Detector file.
Definition JHead.hh:227
double getBaseline(const double costh) const
Get baseline for a given cosine zenith angle.
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition JParser.hh:68
Auxiliary class for organising Monte Carlo file scanners associated with event weighters.
size_t setFlux(const JEvtCategoryHelper &category, const JFluxHelper &flux)
Set flux function for all MC-files corresponding to a given event category.
std::vector< filescanner_type >::iterator iterator
Auxiliary class for defining the range of iterations of objects.
Definition JLimit.hh:45
static counter_type max()
Get maximum counter value.
Definition JLimit.hh:128
Auxiliary base class for list of file names.
Auxiliary data structure for alignment of data.
Definition JManip.hh:298
Auxiliary data structure for floating point format specification.
Definition JManip.hh:488
The Trk class represents a Monte Carlo (MC) particle as well as a reconstructed track/shower.
Definition Trk.hh:15
std::vector< double > fitinf
place to store additional fit info, see km3net-dataformat/definitions/fitparameters....
Definition Trk.hh:32
double E
Energy [GeV] (either MC truth or reconstructed)
Definition Trk.hh:20
double lik
likelihood or lambda value (for aafit, lambda)
Definition Trk.hh:23
Auxiliary methods for selection of reconstructed tracks.
const Trk & get_best_reconstructed_jppshower(const Evt &evt)
Get best reconstructed shower.
const Trk & get_best_reconstructed_jppmuon(const Evt &evt)
Get best reconstructed muon.
bool has_reconstructed_jppshower(const Evt &evt)
Test whether given event has a track with shower reconstruction.
bool has_reconstructed_jppmuon(const Evt &evt)
Test whether given event has a track with muon reconstruction.