Jpp  18.3.1
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 "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 "JPhysics/JConstants.hh"
#include "JGeometry3D/JGeometry3DToolkit.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 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  JFluxMap fluxMaps;
78 
79  string oscProbTable;
80  JOscParameters<> oscParameters;
81 
82  string recotype;
83 
84  int debug;
85 
86  try {
87 
88  JProperties cuts;
89 
90  cuts.insert(gmake_property(energyRange));
91  cuts.insert(gmake_property(coszenithRange));
92  cuts.insert(gmake_property(containmentVolume));
93  cuts.insert(gmake_property(triggeredHitsRange));
94 
95  JParser<> zap("Example program to analyse track fit results from Evt formatted data.");
96 
97  zap['f'] = make_field(inputFiles);
98  zap['n'] = make_field(numberOfEvents)
99  = JLimit::max();
100  zap['a'] = make_field(detectorFile)
101  = "";
102  zap['o'] = make_field(outputFile)
103  = "histogram.root";
104  zap['R'] = make_field(recotype) =
105  track_t(),
106  shower_t();
107  zap['W'] = make_field(useWeights);
108  zap['@'] = make_field(fluxMaps)
109  = JPARSER::initialised();
110  zap['P'] = make_field(oscProbTable)
111  = JPARSER::initialised();
112  zap['#'] = make_field(oscParameters)
113  = JPARSER::initialised();
114  zap['%'] = make_field(cuts)
115  = JPARSER::initialised();
116  zap['d'] = make_field(debug)
117  = 2;
118 
119  zap(argc, argv);
120  }
121  catch(const exception& error) {
122  FATAL(error.what() << endl);
123  }
124 
125  if (useWeights && numberOfEvents != JLimit::max()) {
126  FATAL("Cannot apply weighting to limited number of events.");
127  }
128 
129 
130  // Load detector file
131 
132  JDetector detector;
133 
134  if (detectorFile != "") {
135  try {
136  load(detectorFile, detector);
137  }
138  catch(const JException& error) {
139  FATAL(error);
140  }
141  }
142 
143  const JPMTRouter router(detector);
144 
145 
146  // Set reco getter function pointer
147 
148  bool (*has_reconstruction)(const Evt&);
149  const Trk& (*get_best_fit)(const Evt&);
150 
151  if (recotype == track_t()) {
152  has_reconstruction = &has_reconstructed_jppmuon;
153  get_best_fit = &get_best_reconstructed_jppmuon;
154  } else if (recotype == shower_t()) {
155  has_reconstruction = &has_reconstructed_jppshower;
156  get_best_fit = &get_best_reconstructed_jppshower;
157  } else {
158  FATAL("Invalid recotype: " << recotype);
159  }
160 
161 
162  // Histogram bins
163 
164  const Int_t N_Nhits = 400;
165  const Double_t Nhits_min = 0.0 - 0.5;
166  const Double_t Nhits_max = 1000.0 - 0.5;
167 
168  const Int_t N_Noverlays = 40;
169  const Double_t Noverlays_min = 0.0 - 0.5;
170  const Double_t Noverlays_max = 1000.0 - 0.5;
171 
172  const Int_t N_radius = 201;
173  const Double_t radius_min = -700.0 - 0.5;
174  const Double_t radius_max = +700.0 + 0.5;
175 
176  const Int_t N_height = 18;
177  const Double_t height_min = 20.0;
178  const Double_t height_max = 200.0;
179 
180  const Int_t N_ToTfrac = 20;
181  const Double_t ToTfrac_min = 0.0;
182  const Double_t ToTfrac_max = 1.0;
183 
184  const Int_t N_dt = 1500;
185  const Double_t dt_min = -500.0 - 0.5;
186  const Double_t dt_max = 1000.0 - 0.5;
187 
188  const Int_t N_zenith = 21;
189  const Double_t zenith_min = -1.0;
190  const Double_t zenith_max = +1.0;
191 
192  const Int_t N_energy = 60;
193  const Double_t energy_min = -2.0;
194  const Double_t energy_max = 4.0;
195 
196  const Int_t N_quality = 100;
197  const Double_t quality_min = -200.0;
198  const Double_t quality_max = +800.0;
199 
200  const Int_t N_beta0 = 60;
201  const Double_t beta0_min = -2.0;
202  const Double_t beta0_max = +3.5;
203 
204 
205  // Set of histograms for data/MC-comparisons
206 
207  string ordinateLabel(useWeights ? "Rate [Hz]" : "Number of events");
208 
209  TFile out(outputFile.c_str(), "recreate");
210 
211  TH1D hs ("hs", MAKE_CSTRING("; #snapshot hits; " << ordinateLabel),
212  N_Nhits, Nhits_min, Nhits_max);
213  TH1D hn ("hn", MAKE_CSTRING("; #triggered hits; " << ordinateLabel),
214  N_Nhits, Nhits_min, Nhits_max);
215  TH1D ho ("ho", MAKE_CSTRING("; #overlays; " << ordinateLabel),
216  N_Noverlays, Noverlays_min, Noverlays_max);
217 
218  TH1D hz ("hz", MAKE_CSTRING("; cos(#theta); " << ordinateLabel),
219  N_zenith, zenith_min, zenith_max);
220  TH1D he ("he", MAKE_CSTRING("; log_{10}(E_{reco} / 1 GeV); " << ordinateLabel),
221  N_energy, energy_min, energy_max);
222  TH1D ht ("ht", MAKE_CSTRING("; #Delta t [ns]; " << ordinateLabel),
223  N_dt, dt_min, dt_max);
224 
225  TH1D hZ ("hZ", MAKE_CSTRING("; longitudinal ToT-barycenter [m]; " << ordinateLabel),
226  N_height, height_min, height_max);
227  TH1D hT0 ("hT0", MAKE_CSTRING("; ToT-fraction below earliest triggered hit [ns]; " << ordinateLabel),
228  N_ToTfrac, ToTfrac_min, ToTfrac_max);
229  TH1D hT1 ("hT1", MAKE_CSTRING("; ToT-fraction above earliest triggered hit [ns]; " << ordinateLabel),
230  N_ToTfrac, ToTfrac_min, ToTfrac_max);
231 
232  TH1D hq ("hq", MAKE_CSTRING("; quality; " << ordinateLabel),
233  N_quality, quality_min, quality_max);
234  TH1D hb0 ("hb0", MAKE_CSTRING("; #beta_{0}; " << ordinateLabel),
235  N_beta0, beta0_min, beta0_max);
236 
237  TH2D hzo ("hzo", MAKE_CSTRING("; cos(#theta); #overlays; " << ordinateLabel),
238  N_zenith, zenith_min, zenith_max,
239  N_Noverlays, Noverlays_min, Noverlays_max);
240  TH2D hze ("hze", MAKE_CSTRING("; cos(#theta); log_{10}(E_{reco} / 1 GeV); " << ordinateLabel),
241  N_zenith, zenith_min, zenith_max,
242  N_energy, energy_min, energy_max);
243  TH2D heo ("heo", MAKE_CSTRING("; log_{10}(E_{reco} / 1 GeV); #overlays; " << ordinateLabel),
244  N_energy, energy_min, energy_max,
245  N_Noverlays, Noverlays_min, Noverlays_max);
246  TH2D hzt ("hzt", MAKE_CSTRING("; cos(#theta); #Delta t [ns]; " << ordinateLabel),
247  N_zenith, zenith_min, zenith_max,
248  N_dt, dt_min, dt_max);
249  TH2D het ("het", MAKE_CSTRING("; log_{10}(E_{reco} / 1 GeV); #Delta t [ns]; " << ordinateLabel),
250  N_energy, energy_min, energy_max,
251  N_dt, dt_min, dt_max);
252  TH2D hot ("hot", MAKE_CSTRING("; #overlays; #Delta t [ns]; " << ordinateLabel),
253  N_Noverlays, Noverlays_min, Noverlays_max,
254  N_dt, dt_min, dt_max);
255 
256  TH2D hzq ("hzq", MAKE_CSTRING("; cos(#theta); quality; " << ordinateLabel),
257  N_zenith, zenith_min, zenith_max,
258  N_quality, quality_min, quality_max);
259  TH2D heq ("heq", MAKE_CSTRING("; log_{10}(E_{reco} / 1 GeV); quality; " << ordinateLabel),
260  N_energy, energy_min, energy_max,
261  N_quality, quality_min, quality_max);
262  TH2D hoq ("hoq", MAKE_CSTRING("; #overlays; quality; " << ordinateLabel),
263  N_Noverlays, Noverlays_min, Noverlays_max,
264  N_quality, quality_min, quality_max);
265 
266  TH2D hob0("hob0", MAKE_CSTRING("; #overlays; log_{10}(#beta_{0}); " << ordinateLabel),
267  N_Noverlays, Noverlays_min, Noverlays_max,
268  N_beta0, beta0_min, beta0_max);
269  TH2D heb0("heb0", MAKE_CSTRING("; log_{10}(E_{reco} / 1 GeV); log_{10}(#beta_{0}); " << ordinateLabel),
270  N_energy, energy_min, energy_max,
271  N_beta0, beta0_min, beta0_max);
272  TH2D hzb0("hzb0", MAKE_CSTRING("; cos(#theta); log_{10}(#beta_{0}); " << ordinateLabel),
273  N_zenith, zenith_min, zenith_max,
274  N_beta0, beta0_min, beta0_max);
275 
276  TH2D hxy ("hxy", MAKE_CSTRING("; x [m]; y [m]; " << ordinateLabel),
277  N_radius, radius_min, radius_max,
278  N_radius, radius_min, radius_max);
279 
280 
281  JEvtWeightFileScannerSet<> scanners(inputFiles);
282 
283  if (!oscProbTable.empty()) {
284 
285  JOscProbInterpolator<> interpolator(oscProbTable.c_str(), oscParameters);
286 
287  fluxMaps.configure(make_oscProbFunction(interpolator));
288  }
289 
290  if (scanners.setFlux(fluxMaps) == 0) {
291  WARNING("No flux function set." << endl);
292  }
293 
294  STATUS(RIGHT(15) << "header ID" << RIGHT(15) << "event" << RIGHT(15) << "weight" << endl);
295 
296  for (JEvtWeightFileScannerSet<>::iterator scanner = scanners.begin(); scanner != scanners.end(); ++scanner) {
297 
298  scanner->setLimit(numberOfEvents);
299 
300  while (scanner->hasNext()) {
301 
302  const Evt* evt = scanner->next();
303 
304  const double weight = (useWeights ? scanner->getWeight(*evt) : 1.0);
305 
306  const size_t NsnapshotHits = evt->hits.size();
307  const size_t NtriggeredHits = count_if(evt->hits.begin(), evt->hits.end(),
308  make_predicate(&Hit::trig, (long long unsigned int) 0, JComparison::gt()));
309 
310  if (!triggeredHitsRange(NtriggeredHits)) { continue; }
311 
312  hs .Fill(NsnapshotHits, weight);
313  hn .Fill(NtriggeredHits, weight);
314  ho .Fill(evt->overlays, weight);
315 
316  if (!has_reconstruction(*evt)) { continue; }
317 
318  const Trk& best_trk = (*get_best_fit)(*evt);
319  const JTrack3E& tb = getTrack(best_trk);
320 
321  const double logE = log10(tb.getE());
322 
323  if (!energyRange (tb.getE()) ||
324  !coszenithRange(tb.getDZ()) ||
325  !containmentVolume.is_inside(tb.getPosition())) { // Apply cuts
326  continue;
327  }
328 
329  // Extract event-weight
330 
331  STATUS(RIGHT (15) << distance(scanners.begin(), scanner) <<
332  RIGHT (15) << scanner->getCounter() <<
333  SCIENTIFIC(15, 3) << weight << '\r'); DEBUG(endl);
334 
335  // Fill histograms
336 
337  hz .Fill(tb.getDZ(), weight);
338  he .Fill(logE, weight);
339  hq .Fill(best_trk.lik, weight);
340 
341  hzo.Fill(tb.getDZ(), evt->overlays, weight);
342  hze.Fill(tb.getDZ(), logE, weight);
343  heo.Fill(logE, evt->overlays, weight);
344 
345  hoq.Fill(evt->overlays, best_trk.lik, weight);
346  hzq.Fill(tb.getDZ(), best_trk.lik, weight);
347  heq.Fill(logE, best_trk.lik, weight);
348 
349  hxy.Fill(tb.getX(), tb.getY(), weight);
350 
351  if (best_trk.fitinf[JGANDALF_BETA0_RAD]) {
352 
353  const double logb0 = log10(180.0 / M_PI * best_trk.fitinf[JGANDALF_BETA0_RAD]);
354 
355  hb0 .Fill(logb0, weight);
356  hob0.Fill(evt->overlays, logb0, weight);
357  heb0.Fill(logE, logb0, weight);
358  hzb0.Fill(tb.getDZ(), logb0, weight);
359  }
360 
361  double ToTbelow = 0.0;
362  double ToTabove = 0.0;
363  double ToTtotal = 0.0;
364  double ToTweightedZ = 0.0;
365 
366  vector<Hit>::const_iterator firstHit = min_element(evt->hits.cbegin(), evt->hits.cend(),
367  make_comparator(&Hit::t));
368 
369  for (vector<Hit>::const_iterator hit = evt->hits.cbegin(); hit != evt->hits.cend(); ++hit) {
370 
371  if (router.hasPMT(hit->pmt_id)) {
372 
373  const JPMT& pmt = router.getPMT(hit->pmt_id);
374 
375  const double t0 = tb.getT(pmt);
376  const double t1 = getTime(*hit);
377 
378  ht .Fill(t1 - t0, weight);
379  hot.Fill(evt->overlays, t1 - t0, weight);
380  hzt.Fill(tb.getDZ(), t1 - t0, weight);
381  het.Fill(log10(best_trk.E), t1 - t0, weight);
382 
383  if (hit->trig > 0) {
384 
385  ToTtotal += hit->tot;
386  ToTweightedZ += hit->tot * hit->pos.z;
387 
388  if (hit->pos.z < firstHit->pos.z) {
389  ToTbelow += hit->tot;
390  } else if (hit->pos.z >= firstHit->pos.z) {
391  ToTabove += hit->tot;
392  }
393  }
394  }
395  }
396 
397  hT0.Fill(ToTbelow / ToTtotal, weight);
398  hT1.Fill(ToTabove / ToTtotal, weight);
399  hZ .Fill(ToTweightedZ / ToTtotal, weight);
400  }
401  }
402 
403  STATUS(endl);
404 
405  out.Write();
406  out.Close();
407 
408  return 0;
409 }
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:24
WARNING
#define WARNING(A)
Definition: JMessage.hh:65
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:1293
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:497
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
Evt::overlays
unsigned int overlays
number of overlaying triggered events
Definition: Evt.hh:32
JGEOMETRY3D::JCylinder3D
Cylinder object.
Definition: JCylinder3D.hh:39
JAANET::detector
Detector file.
Definition: JHead.hh:226
JGEOMETRY3D::JTrack3E::getE
double getE() const
Get energy.
Definition: JTrack3E.hh:93
JGEOMETRY3D::getMaximumContainmentVolume
const JCylinder3D getMaximumContainmentVolume()
Auxiliary function to retrieve the maximum cylindrical containment volume.
Definition: JGeometry3DToolkit.hh:250
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})"
JOSCPROB::make_oscProbFunction
JOscProbFunction< JFunction_t > make_oscProbFunction(const JFunction_t &function)
Auxiliary method for creating an interface to an oscillation probability function.
Definition: JOscProbFunction.hh:109
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.2.0-6-ge33aafe https://git.km3net.de/common/km3net-dataformat.
Definition: fitparameters.hh:12
JGEOMETRY3D::JVector3D::getY
double getY() const
Get y position.
Definition: JVector3D.hh:104
JGEOMETRY3D::JPosition3D::getPosition
const JPosition3D & getPosition() const
Get position.
Definition: JPosition3D.hh:130
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
JAANET::JFluxMap
Auxiliary class for parsing multiparticle fluxes.
Definition: JEvtWeightToolkit.hh:203
FATAL
#define FATAL(A)
Definition: JMessage.hh:67
JSUPPORT::JMultipleFileScanner_t
Auxiliary base class for list of file names.
Definition: JMultipleFileScanner.hh:49
JOSCPROB::JOscParameters
Data structure for single set of oscillation parameters.
Definition: JOscParameters.hh:50
JDETECTOR::load
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
Definition: JDetectorToolkit.hh:483
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:61
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:48
Evt::hits
std::vector< Hit > hits
list of hits
Definition: Evt.hh:38
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
SCIENTIFIC
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:486
JGEOMETRY3D::JVersor3D::getDZ
double getDZ() const
Get z direction.
Definition: JVersor3D.hh:117
JOSCPROB::JOscProbInterpolator
Template definition of a multi-dimensional oscillation probability interpolation table.
Definition: JOscProbInterpolator.hh:61
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