Jpp  18.0.1-rc.2
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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 "JLang/JVectorize.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 "JAAnet/JHeadToolkit.hh"
#include "JAAnet/JAAnetToolkit.hh"
#include "JAAnet/JEvtWeightToolkit.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

int main ( int  argc,
char **  argv 
)

Definition at line 56 of file JAAPostfit.cc.

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