Jpp  16.0.3
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/offline/Head.hh"
#include "km3net-dataformat/offline/Evt.hh"
#include "km3net-dataformat/offline/Trk.hh"
#include "km3net-dataformat/tools/reconstruction.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 51 of file JAAMuonPostfit.cc.

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