Jpp  18.2.1-ARCA-DF-PATCH
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 "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 60 of file JAAPostfit.cc.

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