Jpp  master_rocky-43-ge265d140c
the software that should make you happy
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
Cylinder object.
Definition: JCylinder3D.hh:41
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
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
JAANET::getTime
double getTime(const Hit &hit)
Get true time of hit.
Definition: JAAnetToolkit.hh:96
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
JLANG::make_predicate
JPredicate< JResult_t T::*, JComparison::eq > make_predicate(JResult_t T::*member, const JResult_t value)
Helper method to create predicate for data member.
Definition: JPredicate.hh:128
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
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
std
Definition: JSTDTypes.hh:14
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::JMultipleFileScanner_t
Auxiliary base class for list of file names.
Definition: JMultipleFileScanner.hh:51
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
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
get_best_reconstructed_jppmuon
const Trk & get_best_reconstructed_jppmuon(const Evt &evt)
Get best reconstructed muon.
Definition: tools/reconstruction.hh:387
Generated by doxygen 1.9.1