Jpp  16.0.0-rc.1
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Functions
JMuonPostfit.cc File Reference

Example program to histogram fit results for the muon reconstruction chain. More...

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <set>
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TProfile.h"
#include "TMath.h"
#include "km3net-dataformat/offline/Head.hh"
#include "km3net-dataformat/offline/Evt.hh"
#include "km3net-dataformat/online/JDAQModuleIdentifier.hh"
#include "km3net-dataformat/definitions/reconstruction.hh"
#include "km3net-dataformat/definitions/fitparameters.hh"
#include "km3net-dataformat/tools/time_converter.hh"
#include "JAAnet/JHead.hh"
#include "JAAnet/JHeadToolkit.hh"
#include "JAAnet/JAAnetToolkit.hh"
#include "JAAnet/JVolume.hh"
#include "JDAQ/JDAQEventIO.hh"
#include "JPhysics/JConstants.hh"
#include "JTools/JQuantile.hh"
#include "JSupport/JTriggeredFileScanner.hh"
#include "JSupport/JMonteCarloFileSupportkit.hh"
#include "JSupport/JSupport.hh"
#include "JReconstruction/JEvt.hh"
#include "JReconstruction/JEvtToolkit.hh"
#include "JLang/JPredicate.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

Go to the source code of this file.

Functions

int main (int argc, char **argv)
 

Detailed Description

Example program to histogram fit results for the muon reconstruction chain.

Author
mdejong, bofearraigh.

Definition in file JMuonPostfit.cc.

Function Documentation

int main ( int  argc,
char **  argv 
)

Definition at line 79 of file JMuonPostfit.cc.

80 {
81  using namespace std;
82  using namespace JPP;
83  using namespace KM3NETDAQ;
84 
85  typedef JTriggeredFileScanner<JEvt> JTriggeredFileScanner_t;
86  typedef JTriggeredFileScanner_t::multi_pointer_type multi_pointer_type;
87 
88  JTriggeredFileScanner_t inputFile;
89  JLimit_t& numberOfEvents = inputFile.getLimit();
90  string outputFile;
91  size_t numberOfPrefits;
92  JQualitySorter quality_sorter;
93  JAtmosphericMuon atmosphere;
94  double Emin_GeV;
95  double NPE;
96  int application;
97  string option;
98  string primary;
99  int debug;
100 
101  try {
102 
103  JParser<> zap("Example program to histogram fit results.");
104 
105  zap['f'] = make_field(inputFile);
106  zap['o'] = make_field(outputFile) = "postfit.root";
107  zap['n'] = make_field(numberOfEvents) = JLimit::max();
108  zap['N'] = make_field(numberOfPrefits) = 1;
109  zap['L'] = make_field(quality_sorter) = JPARSER::initialised();
110  zap['E'] = make_field(Emin_GeV) = 0.0;
111  zap['M'] = make_field(NPE) = 0.0;
112  zap['A'] = make_field(application) = JMUONPREFIT, JMUONSIMPLEX, JMUONGANDALF, JMUONENERGY, JMUONSTART, JMUONPATH;
113  zap['a'] = make_field(atmosphere) = JAtmosphericMuon(90.0, 90.0);
114  zap['O'] = make_field(option) = "E", "N", "LINE", "LOGE";
115  zap['p'] = make_field(primary) = muon_t, neutrino_t;
116  zap['d'] = make_field(debug) = 2;
117 
118  zap(argc, argv);
119  }
120  catch(const exception& error) {
121  FATAL(error.what() << endl);
122  }
123 
124  JHead head;
125 
126  try {
127  head = getHeader(inputFile);
128  }
129  catch(const JException& error) {
130  FATAL(error);
131  }
132 
133  const JVolume volume(head, option != "LINE");
134  const JPosition3D center = get<JPosition3D>(head);
135  JCylinder3D cylinder = get<JCylinder3D>(head);
136 
137  cylinder.add(center);
138 
139  NOTICE("Center [m]: " << center << endl);
140  NOTICE("Reposition can [m]: " << cylinder << endl);
141 
142 
143  TFile out(outputFile.c_str(), "recreate");
144 
145  TH1D job("job", NULL, 100, 0.5, 100.5);
146 
147  TH1D hn("hn", NULL, 250, -0.5, 249.5); // NDF
148  TH1D hq("hq", NULL, 300, 0.0, 600.0); // quality parameter
149  TH1D hi("hi", NULL, 101, -0.5, 100.5); // index of the best event in order of descending angular resolution
150  TH1D hv("hv", NULL, 200, -6.0, 0.0); // log of Poisson distribution of (#hits,number photo-electrons)
151  TH1D h1("h1", NULL, 200, -2.0, +2.0); // number of photo-electrons along the whole track
152  TH1D hc("hc", NULL, 200, -1.0, +1.0); // dz (z-slope)
153  TH1D hu("hu", NULL, 400, -1.0e3, 1.0e3); // quality difference between best up-going and down-going track (check of atmospheric muon)
154 
155  TH1D hx("hx", NULL, 70, -3.0, +2.3); // angular deviation [log(deg)]
156  TH1D hd("hd", NULL, 100, 0.0, 10.0); // distance between true and best event at intersection point [m]
157  TH1D ht("ht", NULL, 100, -100.0, 100.0); // time between true and best track [ns]
158 
159  TH1D hz0("hz0[start]", NULL, 400, -200.0, 200.0); // start point [m]
160  TH1D hz1("hz1[end]", NULL, 400, -200.0, 200.0); // end point [m]
161 
162  TH1D E_0("E_0", NULL, 100, volume.getXmin(), volume.getXmax()); // true muon energy
163  TH1D E_1("E_1", NULL, 100, volume.getXmin(), volume.getXmax()); // uncorrected energy of best track
164  TH1D E_2("E_2", NULL, 100, volume.getXmin(), volume.getXmax()); // corrected energy of best track
165  TH1D E_E("E_E", NULL, 100, -5.0, +5.0); // ratio of reconstructed energy and true energy [log(E/E)]
166  TH2D ExE("ExE", NULL,
167  40, volume.getXmin(), volume.getXmax(),
168  40, volume.getXmin(), volume.getXmax()); // (Etrue,Ereco)
169 
170  const Int_t ny = 28800;
171  const Double_t ymin = 0.0; // [deg]
172  const Double_t ymax = 180.0; // [deg]
173 
174  vector<double> X;
175 
176  if (option.find('E') != string::npos) {
177 
178  for (double x = volume.getXmin(); x <= volume.getXmax(); x += (volume.getXmax() - volume.getXmin()) / 20) {
179  X.push_back(x);
180  }
181 
182  } else {
183 
184  double x = -0.5;
185 
186  for ( ; x <= 15.5; x += 1.0) { X.push_back(x); }
187  for ( ; x <= 25.5; x += 2.0) { X.push_back(x); }
188  for ( ; x <= 50.5; x += 5.0) { X.push_back(x); }
189  for ( ; x <= 100.5; x += 10.0) { X.push_back(x); }
190  for ( ; x <= 250.5; x += 20.0) { X.push_back(x); }
191  }
192 
193  TH2D h2("h2", NULL, X.size() - 1, X.data(), ny, ymin, ymax);
194 
195  TH2D Va("Va", NULL,
196  100, 0.0, cylinder.getRadius()*cylinder.getRadius(),
197  100, cylinder.getZmin() - 50.0, cylinder.getZmax() + 50.0);
198  TH2D Vb("Vb", NULL,
199  100, 0.0, cylinder.getRadius()*cylinder.getRadius(),
200  100, cylinder.getZmin() - 50.0, cylinder.getZmax() + 50.0);
201 
202  JQuantile Q("Angle", true);
203  JQuantile O("Omega", true);
204 
205 
206  while (inputFile.hasNext()) {
207 
208  STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);
209 
210  multi_pointer_type ps = inputFile.next();
211 
212  JDAQEvent* tev = ps;
213  JEvt* evt = ps;
214  Evt* event = ps;
215 
216  const time_converter converter(*event, *tev);
217 
218  job.Fill(1.0);
219 
220 
221  double Enu = 0.0;
222  double Emu = 0.0;
223  double Emax = 0.0;
224 
225  if (has_neutrino(*event)) {
226  Enu = get_neutrino(*event).E;
227  }
228 
229  vector<Trk>::const_iterator muon = event->mc_trks.end();
230 
231  for (vector<Trk>::const_iterator track = event->mc_trks.begin(); track != event->mc_trks.end(); ++track) {
232 
233  if (is_muon(*track)) {
234 
235  Emu += track->E;
236 
237  if (track->E > Emax) {
238  muon = track;
239  Emax = track->E;
240  }
241  }
242  }
243 
244  if (muon == event->mc_trks.end()) {
245  continue;
246  }
247 
248  job.Fill(3.0);
249 
250 
251  // abscissa
252 
253  Double_t x = 0.0;
254 
255  if (option.find('E') != string::npos)
256  x = volume.getX(event->mc_trks[0].E);
257  else
258  x = getCount(tev->begin<JDAQTriggeredHit>(), tev->end<JDAQTriggeredHit>());
259 
260  if (!evt->empty()) {
261 
262  JEvt::iterator __end = partition(evt->begin(), evt->end(), JHistory::is_event(application));
263 
264  if (evt->begin() == __end) {
265  continue;
266  }
267 
268  job.Fill(4.0);
269 
270  if (numberOfPrefits > 0) {
271 
272  JEvt::iterator __q = __end;
273 
274  advance(__end = evt->begin(), min(numberOfPrefits, (size_t) distance(evt->begin(), __q)));
275 
276  partial_sort(evt->begin(), __end, __q, quality_sorter);
277 
278  } else {
279 
280  sort(evt->begin(), __end, quality_sorter);
281  }
282 
283  double E = 0.0;
284  JPointing pointing;
285 
286  if (primary == muon_t) {
287  E = Emu;
288  pointing = JPointing(getDirection(*muon));
289  } else if (primary == neutrino_t) {
290  E = Enu;
291  pointing = JPointing(getDirection(get_neutrino(*event)));
292  }
293 
294  JEvt::iterator best = pointing(evt->begin(), __end);
295  const Double_t beta = pointing.getAngle(*best);
296  const double Efit = best->getE();
297  const double Eraw = best->getW(JENERGY_ENERGY, numeric_limits<double>::min());
298  const double mip = best->getW(JSTART_NPE_MIP, numeric_limits<double>::max());
299 
300  // selection of fit result
301 
302  bool ok = (Efit >= Emin_GeV && mip >= NPE);
303 
304  if (ok) {
305 
306  job.Fill(5.0);
307 
308  hn.Fill((Double_t) best->getNDF());
309  hq.Fill(best->getQ());
310  hi.Fill((Double_t) distance(evt->begin(), best));
311  hc.Fill(best->getDZ());
312 
313  // atmospheric muon hypothesis test.
314 
315  if (( has_neutrino(*event) && get_neutrino(*event).dir.z >= atmosphere.dot2) || // difference in quality between best upwards and best downgoing track.
316  (!has_neutrino(*event) && best->getDZ() >= atmosphere.dot2)) { // negative (positive) values imply a down(up) going track.
317  hu.Fill(atmosphere(evt->begin(), __end));
318  }
319 
320  hx.Fill(max(log10(beta), hx.GetXaxis()->GetXmin()));
321 
322  Q.put(beta);
323 
324  JTrack3E ta = getTrack(*muon);
325  JTrack3E tb = getTrack(*best);
326 
327  ta.add(center);
328  tb.sub(converter.putTime());
329 
330  static_cast<JTrack3D&>(ta).move(ta.getIntersection(tb), getSpeedOfLight()); // move to intersection point of both tracks
331  static_cast<JTrack3D&>(tb).move(tb.getIntersection(ta), getSpeedOfLight());
332 
333  hd.Fill((tb.getPosition() - ta.getPosition()).getLength());
334 
335  if (has_neutrino(*event)) {
336 
337  job.Fill(6.0);
338 
339  const Trk neutrino = get_neutrino(*event);
340 
341  JPosition3D vertex = getPosition(neutrino);
342 
343  vertex.add(center);
344 
345  if (cylinder.is_inside(vertex)) {
346 
347  job.Fill(7.0);
348 
349  Va.Fill(best->getX()*best->getX() + best->getY()*best->getY(), best->getZ()); //R^{2} [m^2], z [m]
350 
351  JTrack3E tc = getTrack(*best);
352 
353  hz0.Fill(tc.getIntersection(vertex));
354  hz1.Fill(best->getW(JSTART_LENGTH_METRES, 0.0) - gWater(muon->E));
355  }
356  }
357 
358  Vb.Fill(best->getX()*best->getX() + best->getY()*best->getY(), best->getZ()); //R^{2} [m^2], z [m]
359  ht.Fill(tb.getT() - ta.getT());
360 
361  E_0.Fill(volume.getX(E, true));
362  E_1.Fill(volume.getX(Eraw, true));
363  E_2.Fill(volume.getX(Efit, true));
364  E_E.Fill(volume.getX(Efit) - volume.getX(E));
365  ExE.Fill(volume.getX(E), volume.getX(Efit));
366 
367  h2.Fill(x, beta);
368 
369  if (best->hasW(JSTART_NPE_MIP)) {
370  h1.Fill(log10(best->getW(JSTART_NPE_MIP)));
371  }
372 
373  if (best->hasW(JVETO_NPE) && best->hasW(JVETO_NUMBER_OF_HITS)) {
374 
375  const double npe = best->getW(JVETO_NPE); // number of photo-electrons
376  const int count = best->getW(JVETO_NUMBER_OF_HITS); // number of hits
377  const double pv = TMath::PoissonI(count, npe); // Poisson distribution function for (#hits,npe)
378 
379  hv.Fill(max(log10(pv), hv.GetXaxis()->GetXmin()));
380  }
381  }
382  }
383  }
384  STATUS(endl);
385 
386  NOTICE("Number of events input " << setw(8) << right << job.GetBinContent(1) << endl);
387  NOTICE("Number of events with muon " << setw(8) << right << job.GetBinContent(3) << endl);
388  NOTICE("Number of events with fit " << setw(8) << right << job.GetBinContent(4) << endl);
389  NOTICE("Number of events selected " << setw(8) << right << job.GetBinContent(5) << endl);
390  NOTICE("Number of events with neutrino " << setw(8) << right << job.GetBinContent(6) << endl);
391  NOTICE("Number of events contained " << setw(8) << right << job.GetBinContent(7) << endl);
392 
393  if (Q.getCount() != 0) {
394  NOTICE("Median space angle [deg] " << FIXED (6,3) << Q.getQuantile(0.5) << endl);
395  }
396 
397  out.Write();
398  out.Close();
399 }
JMUONSTART
static const int JMUONSTART
Definition: definitions/reconstruction.hh:18
JPARSER::JParser
Utility class to parse command line options.
Definition: JParser.hh:1500
JLANG::JException
General exception.
Definition: JException.hh:23
Q
Q(UTCMax_s-UTCMin_s)-livetime_s
JMUONPATH
static const int JMUONPATH
Definition: definitions/reconstruction.hh:44
JVETO_NPE
static const int JVETO_NPE
number of photo-electrons from JVeto.cc
Definition: fitparameters.hh:23
JRECONSTRUCTION::JPointing
Auxiliary class to compare fit results with respect to a reference direction (e.g. true muon).
Definition: JReconstruction/JEvtToolkit.hh:523
E
then usage E
Definition: JMuonPostfit.sh:35
JAANET::getTrack
JTrack3E getTrack(const Trk &track)
Get track.
Definition: JAAnetToolkit.hh:250
JENERGY_ENERGY
static const int JENERGY_ENERGY
uncorrected energy [GeV] from JEnergy.cc
Definition: fitparameters.hh:16
Vec::z
double z
Definition: Vec.hh:14
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
JGEOMETRY3D::JTime::sub
JTime & sub(const JTime &value)
Subtraction operator.
Definition: JGeometry3D/JTime.hh:81
mc_keys::neutrino_t
const char *const neutrino_t
Definition: io_ascii.hh:29
JTOOLS::JQuantile
Auxiliary data structure for running average, standard deviation and quantiles.
Definition: JQuantile.hh:43
STATUS
#define STATUS(A)
Definition: JMessage.hh:63
JAANET::has_neutrino
bool has_neutrino(const Evt &evt)
Test whether given event has an incoming neutrino.
Definition: JAAnetToolkit.hh:427
JGEOMETRY3D::JAxis3D::getIntersection
double getIntersection(const JVector3D &pos) const
Get longitudinal position along axis of position of closest approach with given position.
Definition: JAxis3D.hh:146
JAANET::is_muon
bool is_muon(const Trk &track)
Test whether given track is a (anti-)muon.
Definition: JAAnetToolkit.hh:330
Trk::dir
Vec dir
track direction
Definition: Trk.hh:18
JRECONSTRUCTION::JQualitySorter
General purpose sorter of fit results.
Definition: JReconstruction/JEvtToolkit.hh:259
JSUPPORT::JTriggeredFileScanner
Auxiliary class to synchronously read DAQ events and Monte Carlo events (and optionally other events)...
Definition: JTriggeredFileScanner.hh:39
h1
then for HISTOGRAM in h0 h1
Definition: JMatrixNZ.sh:71
JPARSER::initialised
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition: JParser.hh:66
JPHYSICS::gWater
static const JGeaneWater gWater
Function object for energy loss of muon in sea water.
Definition: JGeane.hh:323
time_converter
Auxiliary class to convert DAQ hit time to/from Monte Carlo hit time.
Definition: time_converter.hh:19
JGEOMETRY3D::JCylinder3D::add
JCylinder3D & add(const JVector3D &pos)
Add position.
Definition: JCylinder3D.hh:156
FIXED
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:446
outputFile
string outputFile
Definition: JDAQTimesliceSelector.cc:37
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
JSTART_NPE_MIP
static const int JSTART_NPE_MIP
number of photo-electrons up to the barycentre from JStart.cc
Definition: fitparameters.hh:20
JGEOMETRY3D::JTime::add
JTime & add(const JTime &value)
Addition operator.
Definition: JGeometry3D/JTime.hh:67
JSUPPORT::getHeader
Head getHeader(const JMultipleFileScanner_t &file_list)
Get Monte Carlo header.
Definition: JMonteCarloFileSupportkit.hh:471
JAANET::JVolume
Auxiliary class for histogramming of effective volume.
Definition: JVolume.hh:29
JSUPPORT::JLimit
Auxiliary class for defining the range of iterations of objects.
Definition: JLimit.hh:41
KM3NETDAQ::JDAQEvent::end
const_iterator< T > end() const
Get end of data.
JMUONPREFIT
static const int JMUONPREFIT
Definition: definitions/reconstruction.hh:14
JGEOMETRY3D::JCylinder3D
Cylinder object.
Definition: JCylinder3D.hh:39
JRECONSTRUCTION::JPointing::getAngle
double getAngle(const JFit &fit) const
Get angle between reference direction and fit result.
Definition: JReconstruction/JEvtToolkit.hh:572
JAANET::getDirection
JDirection3D getDirection(const Vec &dir)
Get direction.
Definition: JAAnetToolkit.hh:202
JACOUSTICS::JEvt
Acoustic event fit.
Definition: JAcoustics/JEvt.hh:184
KM3NETDAQ::JDAQEvent::begin
const_iterator< T > begin() const
Get begin of data.
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})"
JFIT::JHistory::is_event
Auxiliary class to test history.
Definition: JHistory.hh:104
JMUONGANDALF
static const int JMUONGANDALF
Definition: definitions/reconstruction.hh:16
JAANET::getPosition
JPosition3D getPosition(const Vec &pos)
Get position.
Definition: JAAnetToolkit.hh:166
NOTICE
#define NOTICE(A)
Definition: JMessage.hh:64
X
then break fi done getCenter read X Y Z let X
Definition: acoustics-fit.sh:523
JVETO_NUMBER_OF_HITS
static const int JVETO_NUMBER_OF_HITS
number of hits from JVeto.cc
Definition: fitparameters.hh:24
debug
int debug
debug level
Definition: JSirene.cc:63
JGEOMETRY3D::JPosition3D::getPosition
const JPosition3D & getPosition() const
Get position.
Definition: JPosition3D.hh:130
quality_sorter
Reconstruction type dependent comparison of track quality.
Definition: tools/reconstruction.hh:85
JRECONSTRUCTION::JAtmosphericMuon
Auxiliary class to evaluate atmospheric muon hypothesis.
Definition: JReconstruction/JEvtToolkit.hh:709
JGEOMETRY3D::JTrack3D::getT
double getT(const JVector3D &pos) const
Get arrival time of Cherenkov light at given position.
Definition: JTrack3D.hh:126
JROOT::advance
counter_type advance(counter_type &counter, const counter_type value, const counter_type limit=std::numeric_limits< counter_type >::max())
Advance counter.
Definition: JROOT/JCounter.hh:35
JAANET::JHead
Monte Carlo run header.
Definition: JHead.hh:1164
FATAL
#define FATAL(A)
Definition: JMessage.hh:67
JPHYSICS::getSpeedOfLight
const double getSpeedOfLight()
Get speed of light.
Definition: JPhysics/JConstants.hh:105
JSTART_LENGTH_METRES
static const int JSTART_LENGTH_METRES
distance between first and last hits in metres from JStart.cc
Definition: fitparameters.hh:22
JTRIGGER::count
std::vector< int > count
Definition: JAlgorithm.hh:180
JMUONSIMPLEX
static const int JMUONSIMPLEX
Definition: definitions/reconstruction.hh:15
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
JTRIGGER::getCount
int getCount(const T &hit)
Get hit count.
Definition: JTrigger/JHit.hh:277
JAANET::primary
Primary particle.
Definition: JHead.hh:1104
JGEOMETRY3D::JPosition3D
Data structure for position in three dimensions.
Definition: JPosition3D.hh:36
JSUPPORT::JLimit::getLimit
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:73
JAANET::get_neutrino
const Trk & get_neutrino(const Evt &evt)
Get incoming neutrino.
Definition: JAAnetToolkit.hh:438
Trk
The Trk class represents a Monte Carlo (MC) particle as well as a reconstructed track/shower.
Definition: Trk.hh:14
JGEOMETRY3D::JVector3D::add
JVector3D & add(const JVector3D &vector)
Add vector.
Definition: JVector3D.hh:142
JMUONENERGY
static const int JMUONENERGY
Definition: definitions/reconstruction.hh:17
Evt
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
Definition: Evt.hh:19
Generated by   doxygen 1.8.5