Jpp  17.3.2
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Functions
JAAMuonPostfit.cc File Reference

Example program to analyse track 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 "JDetector/JDetector.hh"
#include "JDetector/JPMTRouter.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JAAnet/JHeadToolkit.hh"
#include "JAAnet/JAAnetToolkit.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 track fit results from Evt formatted data.

Author
bofearraigh, bjung

Definition in file JAAMuonPostfit.cc.

Function Documentation

int main ( int  argc,
char **  argv 
)

Definition at line 54 of file JAAMuonPostfit.cc.

55 {
56  using namespace std;
57  using namespace JPP;
58  using namespace KM3NETDAQ;
59 
60  JMultipleFileScanner_t inputFiles;
61 
62  JLimit_t numberOfEvents;
63 
64  string detectorFile;
65  string outputFile;
66 
67  JRange<double> energyRange;
68  JRange<double> coszenithRange;
69  JRange<unsigned int> triggeredHitsRange;
70 
71  string weightOption;
72  int debug;
73 
74  try {
75 
76  JProperties cuts;
77 
78  cuts.insert(gmake_property(energyRange));
79  cuts.insert(gmake_property(coszenithRange));
80  cuts.insert(gmake_property(triggeredHitsRange));
81 
82  JParser<> zap("Example program to analyse track fit results from Evt formatted data.");
83 
84  zap['f'] = make_field(inputFiles);
85  zap['n'] = make_field(numberOfEvents) = JLimit::max();
86  zap['a'] = make_field(detectorFile) = "";
87  zap['o'] = make_field(outputFile) = "histogram.root";
88  zap['@'] = make_field(cuts) = JPARSER::initialised();
89  zap['W'] = make_field(weightOption) =
90  JEvtWeightOption::unweighted(),
91  JEvtWeightOption::standard(),
92  JEvtWeightOption::rescaled();
93  zap['d'] = make_field(debug) = 2;
94 
95  zap(argc, argv);
96  }
97  catch(const exception& error) {
98  FATAL(error.what() << endl);
99  }
100 
101  if (weightOption != JEvtWeightOption::unweighted() &&
102  numberOfEvents != JLimit::max()) {
103  FATAL("Cannot apply weighting to limited number of events.");
104  }
105 
106 
107  JDetector detector;
108 
109  if (detectorFile != "") {
110  try {
111  load(detectorFile, detector);
112  }
113  catch(const JException& error) {
114  FATAL(error);
115  }
116  }
117 
118 
119  const JPMTRouter router(detector);
120 
121  // Histogram bins
122 
123  const Int_t N_Nhits = 400;
124  const Double_t Nhits_min = 0.0 - 0.5;
125  const Double_t Nhits_max = 1000.0 - 0.5;
126 
127  const Int_t N_Noverlays = 40;
128  const Double_t Noverlays_min = 0.0 - 0.5;
129  const Double_t Noverlays_max = 1000.0 - 0.5;
130 
131  const Int_t N_radius = 201;
132  const Double_t radius_min = -700.0 - 0.5;
133  const Double_t radius_max = +700.0 + 0.5;
134 
135  const Int_t N_height = 18;
136  const Double_t height_min = 20.0;
137  const Double_t height_max = 200.0;
138 
139  const Int_t N_ToTfrac = 20;
140  const Double_t ToTfrac_min = 0.0;
141  const Double_t ToTfrac_max = 1.0;
142 
143  const Int_t N_dt = 1500;
144  const Double_t dt_min = -500.0 - 0.5;
145  const Double_t dt_max = 1000.0 - 0.5;
146 
147  const Int_t N_zenith = 21;
148  const Double_t zenith_min = -1.0;
149  const Double_t zenith_max = +1.0;
150 
151  const Int_t N_energy = 60;
152  const Double_t energy_min = -2.0;
153  const Double_t energy_max = 4.0;
154 
155  const Int_t N_quality = 100;
156  const Double_t quality_min = -200.0;
157  const Double_t quality_max = +800.0;
158 
159  const Int_t N_beta0 = 60;
160  const Double_t beta0_min = -2.0;
161  const Double_t beta0_max = +3.5;
162 
163 
164  // Set of histograms for data/MC-comparisons
165 
166  string ordinateLabel(weightOption == JEvtWeightOption::unweighted() ?
167  "Number of events" : "Rate [Hz]");
168 
169  TFile out(outputFile.c_str(), "recreate");
170 
171  TH1D hs ("hs", MAKE_CSTRING("; #snapshot hits; " << ordinateLabel),
172  N_Nhits, Nhits_min, Nhits_max);
173  TH1D hn ("hn", MAKE_CSTRING("; #triggered hits; " << ordinateLabel),
174  N_Nhits, Nhits_min, Nhits_max);
175  TH1D ho ("ho", MAKE_CSTRING("; #overlays; " << ordinateLabel),
176  N_Noverlays, Noverlays_min, Noverlays_max);
177 
178  TH1D hz ("hz", MAKE_CSTRING("; cos(#theta); " << ordinateLabel),
179  N_zenith, zenith_min, zenith_max);
180  TH1D he ("he", MAKE_CSTRING("; log_{10}(E_{reco} / 1 GeV); " << ordinateLabel),
181  N_energy, energy_min, energy_max);
182  TH1D ht ("ht", MAKE_CSTRING("; #Delta t [ns]; " << ordinateLabel),
183  N_dt, dt_min, dt_max);
184 
185  TH1D hZ ("hZ", MAKE_CSTRING("; longitudinal ToT-barycenter [m]; " << ordinateLabel),
186  N_height, height_min, height_max);
187  TH1D hT0 ("hT0", MAKE_CSTRING("; ToT-fraction below earliest triggered hit [ns]; " << ordinateLabel),
188  N_ToTfrac, ToTfrac_min, ToTfrac_max);
189  TH1D hT1 ("hT1", MAKE_CSTRING("; ToT-fraction above earliest triggered hit [ns]; " << ordinateLabel),
190  N_ToTfrac, ToTfrac_min, ToTfrac_max);
191 
192  TH1D hq ("hq", MAKE_CSTRING("; quality; " << ordinateLabel),
193  N_quality, quality_min, quality_max);
194  TH1D hb0 ("hb0", MAKE_CSTRING("; #beta_{0}; " << ordinateLabel),
195  N_beta0, beta0_min, beta0_max);
196 
197  TH2D hzo ("hzo", MAKE_CSTRING("; cos(#theta); #overlays; " << ordinateLabel),
198  N_zenith, zenith_min, zenith_max,
199  N_Noverlays, Noverlays_min, Noverlays_max);
200  TH2D hze ("hze", MAKE_CSTRING("; cos(#theta); log_{10}(E_{reco} / 1 GeV); " << ordinateLabel),
201  N_zenith, zenith_min, zenith_max,
202  N_energy, energy_min, energy_max);
203  TH2D heo ("heo", MAKE_CSTRING("; log_{10}(E_{reco} / 1 GeV); #overlays; " << ordinateLabel),
204  N_energy, energy_min, energy_max,
205  N_Noverlays, Noverlays_min, Noverlays_max);
206  TH2D hzt ("hzt", MAKE_CSTRING("; cos(#theta); #Delta t [ns]; " << ordinateLabel),
207  N_zenith, zenith_min, zenith_max,
208  N_dt, dt_min, dt_max);
209  TH2D het ("het", MAKE_CSTRING("; log_{10}(E_{reco} / 1 GeV); #Delta t [ns]; " << ordinateLabel),
210  N_energy, energy_min, energy_max,
211  N_dt, dt_min, dt_max);
212  TH2D hot ("hot", MAKE_CSTRING("; #overlays; #Delta t [ns]; " << ordinateLabel),
213  N_Noverlays, Noverlays_min, Noverlays_max,
214  N_dt, dt_min, dt_max);
215 
216  TH2D hzq ("hzq", MAKE_CSTRING("; cos(#theta); quality; " << ordinateLabel),
217  N_zenith, zenith_min, zenith_max,
218  N_quality, quality_min, quality_max);
219  TH2D heq ("heq", MAKE_CSTRING("; log_{10}(E_{reco} / 1 GeV); quality; " << ordinateLabel),
220  N_energy, energy_min, energy_max,
221  N_quality, quality_min, quality_max);
222  TH2D hoq ("hoq", MAKE_CSTRING("; #overlays; quality; " << ordinateLabel),
223  N_Noverlays, Noverlays_min, Noverlays_max,
224  N_quality, quality_min, quality_max);
225 
226  TH2D hob0("hob0", MAKE_CSTRING("; #overlays; log_{10}(#beta_{0}); " << ordinateLabel),
227  N_Noverlays, Noverlays_min, Noverlays_max,
228  N_beta0, beta0_min, beta0_max);
229  TH2D heb0("heb0", MAKE_CSTRING("; log_{10}(E_{reco} / 1 GeV); log_{10}(#beta_{0}); " << ordinateLabel),
230  N_energy, energy_min, energy_max,
231  N_beta0, beta0_min, beta0_max);
232  TH2D hzb0("hzb0", MAKE_CSTRING("; cos(#theta); log_{10}(#beta_{0}); " << ordinateLabel),
233  N_zenith, zenith_min, zenith_max,
234  N_beta0, beta0_min, beta0_max);
235 
236  TH2D hxy ("hxy", MAKE_CSTRING("; x [m]; y [m]; " << ordinateLabel),
237  N_radius, radius_min, radius_max,
238  N_radius, radius_min, radius_max);
239 
240 
241  JEvtWeightFileScannerSet<> scanners(inputFiles);
242 
243  STATUS(RIGHT(15) << "header ID" << RIGHT(15) << "event" << RIGHT(15) << "weight" << endl);
244 
245  for (JEvtWeightFileScannerSet<>::iterator scanner = scanners.begin(); scanner != scanners.end(); ++scanner) {
246 
247  scanner->setLimit(numberOfEvents);
248 
249  while (scanner->hasNext()) {
250 
251  const Evt* evt = scanner->next();
252 
253  if (has_reconstructed_jppmuon(*evt)) {
254 
255  DEBUG(get_best_reconstructed_jppmuon(*evt).comment << endl);
256 
257  const Trk& best_trk = get_best_reconstructed_jppmuon(*evt);
258  const JTrack3E& tb = getTrack(best_trk);
259 
260  const double logE = log10(tb.getE());
261  const size_t NsnapshotHits = evt->hits.size();
262  const size_t NtriggeredHits = count_if(evt->hits.begin(), evt->hits.end(),
263  make_predicate(&Hit::trig, (long long unsigned int) 0, JComparison::gt()));
264 
265  if (triggeredHitsRange(NtriggeredHits) &&
266  coszenithRange (tb.getDZ()) &&
267  energyRange (tb.getE())) { // Apply cuts
268 
269  // Extract event-weight
270 
271  double weight = 0.0;
272 
273  if (weightOption == JEvtWeightOption::unweighted()) {
274  weight = 1.0;
275  } else if (weightOption == JEvtWeightOption::standard()) {
276  weight = scanner->getWeight(*evt);
277  } else if (weightOption == JEvtWeightOption::rescaled()) {
278  weight = evt->w.at(WEIGHTLIST_RESCALED_EVENT_RATE);
279  }
280 
281  STATUS(RIGHT (15) << distance(scanners.begin(), scanner) <<
282  RIGHT (15) << scanner->getCounter() <<
283  SCIENTIFIC(15, 3) << weight << '\r'); DEBUG(endl);
284 
285  // Fill histograms
286 
287  hs .Fill(NsnapshotHits, weight);
288  hn .Fill(NtriggeredHits, weight);
289  ho .Fill(evt->overlays, weight);
290 
291  hz .Fill(tb.getDZ(), weight);
292  he .Fill(logE, weight);
293  hq .Fill(best_trk.lik, weight);
294 
295  hzo.Fill(tb.getDZ(), evt->overlays, weight);
296  hze.Fill(tb.getDZ(), logE, weight);
297  heo.Fill(logE, evt->overlays, weight);
298 
299  hoq.Fill(evt->overlays, best_trk.lik, weight);
300  hzq.Fill(tb.getDZ(), best_trk.lik, weight);
301  heq.Fill(logE, best_trk.lik, weight);
302 
303  hxy.Fill(tb.getX(), tb.getY(), weight);
304 
305  if (best_trk.fitinf[JGANDALF_BETA0_RAD]) {
306 
307  const double logb0 = log10(180.0 / M_PI * best_trk.fitinf[JGANDALF_BETA0_RAD]);
308 
309  hb0 .Fill(logb0, weight);
310  hob0.Fill(evt->overlays, logb0, weight);
311  heb0.Fill(logE, logb0, weight);
312  hzb0.Fill(tb.getDZ(), logb0, weight);
313  }
314 
315  double ToTbelow = 0.0;
316  double ToTabove = 0.0;
317  double ToTtotal = 0.0;
318  double ToTweightedZ = 0.0;
319 
320  vector<Hit>::const_iterator firstHit = min_element(evt->hits.cbegin(), evt->hits.cend(),
321  make_comparator(&Hit::t));
322 
323  for (vector<Hit>::const_iterator hit = evt->hits.cbegin(); hit != evt->hits.cend(); ++hit) {
324 
325  if (router.hasPMT(hit->pmt_id)) {
326 
327  const JPMT& pmt = router.getPMT(hit->pmt_id);
328 
329  const double t0 = tb.getT(pmt);
330  const double t1 = getTime(*hit);
331 
332  ht .Fill(t1 - t0, weight);
333  hot.Fill(evt->overlays, t1 - t0, weight);
334  hzt.Fill(tb.getDZ(), t1 - t0, weight);
335  het.Fill(log10(best_trk.E), t1 - t0, weight);
336 
337  if (hit->trig > 0) {
338 
339  ToTtotal += hit->tot;
340  ToTweightedZ += hit->tot * hit->pos.z;
341 
342  if (hit->pos.z < firstHit->pos.z) {
343  ToTbelow += hit->tot;
344  } else if (hit->pos.z >= firstHit->pos.z) {
345  ToTabove += hit->tot;
346  }
347  }
348  }
349  }
350 
351  hT0.Fill(ToTbelow / ToTtotal, weight);
352  hT1.Fill(ToTabove / ToTtotal, weight);
353  hZ .Fill(ToTweightedZ / ToTtotal, weight);
354  }
355  }
356  }
357  }
358 
359  STATUS(endl);
360 
361  out.Write();
362  out.Close();
363 }
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:1517
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
WEIGHTLIST_RESCALED_EVENT_RATE
static const int WEIGHTLIST_RESCALED_EVENT_RATE
Rescaled event rate [s-1].
Definition: weightlist.hh:17
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
gmake_property
#define gmake_property(A)
macro 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
Trk::comment
std::string comment
use as you like
Definition: Trk.hh:35
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
Evt::w
std::vector< double > w
MC: Weights w[0]=w1, w[1]=w2, w[2]=w3 (see e.g. Tag list or km3net-dataformat/definitions) ...
Definition: Evt.hh:42
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
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:1993
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.0.0-3-gef79250 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::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:486
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
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:484
JGEOMETRY3D::JVersor3D::getDZ
double getDZ() const
Get z direction.
Definition: JVersor3D.hh:117
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