Jpp  15.0.1-rc.2-highQE
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JMuonPostfit.cc
Go to the documentation of this file.
1 #include <string>
2 #include <iostream>
3 #include <iomanip>
4 #include <vector>
5 #include <set>
6 
7 #include "TROOT.h"
8 #include "TFile.h"
9 #include "TH1D.h"
10 #include "TH2D.h"
11 #include "TProfile.h"
12 #include "TMath.h"
13 
20 #include "JAAnet/JHead.hh"
21 #include "JAAnet/JHeadToolkit.hh"
22 #include "JAAnet/JAAnetToolkit.hh"
23 #include "JAAnet/JVolume.hh"
24 #include "JDAQ/JDAQEventIO.hh"
25 #include "JPhysics/JConstants.hh"
26 #include "JTools/JQuantile.hh"
29 #include "JSupport/JSupport.hh"
30 #include "JReconstruction/JEvt.hh"
32 #include "JLang/JPredicate.hh"
33 
34 #include "Jeep/JParser.hh"
35 #include "Jeep/JMessage.hh"
36 
37 
38 namespace {
39 
41 
42 
43  /**
44  * Count number of unique identifiers in event.
45  *
46  * \param __begin begin of data
47  * \param __end end of data
48  * \return number of unique identifiers
49  */
50  template<class JHit_t>
51  inline int getCount(JDAQEvent::const_iterator<JHit_t> __begin,
53  {
54  using namespace std;
55  using namespace KM3NETDAQ;
56 
57  typedef JDAQModuleIdentifier JType_t;
58 
59  set<JType_t> buffer;
60 
61  for (JDAQEvent::const_iterator<JHit_t> i = __begin; i != __end; ++i) {
62  buffer.insert(JType_t(*i));
63  }
64 
65  return buffer.size();
66  }
67 
68 
69  const char* const muon_t = "muon";
70  const char* const neutrino_t = "neutrino";
71 }
72 
73 /**
74  * \file
75  *
76  * Example program to histogram fit results for the muon reconstruction chain.
77  * \author mdejong, bofearraigh.
78  */
79 int main(int argc, char **argv)
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;
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;
110  zap['E'] = make_field(Emin_GeV) = 0.0;
111  zap['M'] = make_field(NPE) = 0.0;
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 
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 }
static const int JMUONSTART
Utility class to parse command line options.
Definition: JParser.hh:1500
General exception.
Definition: JException.hh:23
Double_t getXmin() const
Get minimal abscissa value.
Definition: JVolume.hh:89
Q(UTCMax_s-UTCMin_s)-livetime_s
int main(int argc, char *argv[])
Definition: Main.cc:15
ROOT TTree parameter settings of various packages.
static const int JMUONPATH
static const int JVETO_NPE
number of photo-electrons from JVeto.cc
Auxiliary class to compare fit results with respect to a reference direction (e.g. true muon).
JTrack3E getTrack(const Trk &track)
Get track.
static const int JENERGY_ENERGY
uncorrected energy [GeV] from JEnergy.cc
double z
Definition: Vec.hh:14
Synchronously read DAQ events and Monte Carlo events (and optionally other events).
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
JTime & sub(const JTime &value)
Subtraction operator.
const char *const neutrino_t
Definition: io_ascii.hh:29
Auxiliary data structure for running average, standard deviation and quantiles.
Definition: JQuantile.hh:43
#define STATUS(A)
Definition: JMessage.hh:63
bool has_neutrino(const Evt &evt)
Test whether given event has an incoming neutrino.
Template const_iterator.
Definition: JDAQEvent.hh:68
double getIntersection(const JVector3D &pos) const
Get longitudinal position along axis of position of closest approach with given position.
Definition: JAxis3D.hh:146
bool is_muon(const Trk &track)
Test whether given track is a (anti-)muon.
Vec dir
track direction
Definition: Trk.hh:18
General purpose sorter of fit results.
Auxiliary class to synchronously read DAQ events and Monte Carlo events (and optionally other events)...
Double_t getXmax() const
Get maximal abscissa value.
Definition: JVolume.hh:100
then for HISTOGRAM in h0 h1
Definition: JMatrixNZ.sh:71
Double_t getX(const Double_t E, double constrain=false) const
Get abscissa value.
Definition: JVolume.hh:113
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition: JParser.hh:66
static const JGeaneWater gWater
Function object for energy loss of muon in sea water.
Definition: JGeane.hh:323
Auxiliary class to convert DAQ hit time to/from Monte Carlo hit time.
JCylinder3D & add(const JVector3D &pos)
Add position.
Definition: JCylinder3D.hh:156
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:446
string outputFile
double E
Energy [GeV] (either MC truth or reconstructed)
Definition: Trk.hh:20
3D track with energy.
Definition: JTrack3E.hh:30
then usage E
Definition: JMuonPostfit.sh:35
double getQuantile(const double Q, const bool reverse=false) const
Get quantile.
Definition: JQuantile.hh:309
static const int JSTART_NPE_MIP
number of photo-electrons up to the barycentre from JStart.cc
JTime & add(const JTime &value)
Addition operator.
Head getHeader(const JMultipleFileScanner_t &file_list)
Get Monte Carlo header.
Auxiliary class for histogramming of effective volume.
Definition: JVolume.hh:26
Auxiliary class for defining the range of iterations of objects.
Definition: JLimit.hh:41
const_iterator< T > end() const
Get end of data.
static const int JMUONPREFIT
Cylinder object.
Definition: JCylinder3D.hh:39
double getAngle(const JFit &fit) const
Get angle between reference direction and fit result.
JDirection3D getDirection(const Vec &dir)
Get direction.
Acoustic event fit.
const_iterator< T > begin() const
Get begin of data.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1961
Auxiliary class to test history.
Definition: JHistory.hh:104
static const int JMUONGANDALF
JPosition3D getPosition(const Vec &pos)
Get position.
Definition of hit and track types and auxiliary methods for handling Monte Carlo data.
#define NOTICE(A)
Definition: JMessage.hh:64
then break fi done getCenter read X Y Z let X
void put(const double x, const double w=1.0)
Put value.
Definition: JQuantile.hh:133
Physics constants.
double putTime() const
Get Monte Carlo time minus DAQ/trigger time.
static const int JVETO_NUMBER_OF_HITS
number of hits from JVeto.cc
int debug
debug level
Definition: JSirene.cc:63
const JPosition3D & getPosition() const
Get position.
Definition: JPosition3D.hh:130
long long int getCount() const
Get total count.
Definition: JQuantile.hh:186
Reconstruction type dependent comparison of track quality.
Auxiliary class to evaluate atmospheric muon hypothesis.
double getT(const JVector3D &pos) const
Get arrival time of Cherenkov light at given position.
Definition: JTrack3D.hh:126
General purpose messaging.
counter_type advance(counter_type &counter, const counter_type value, const counter_type limit=std::numeric_limits< counter_type >::max())
Advance counter.
Auxiliary include file for time conversion between DAQ/trigger hit and Monte Carlo hit...
Monte Carlo run header.
Definition: JHead.hh:1113
#define FATAL(A)
Definition: JMessage.hh:67
const double getSpeedOfLight()
Get speed of light.
static const int JSTART_LENGTH_METRES
distance between first and last hits in metres from JStart.cc
std::vector< int > count
Definition: JAlgorithm.hh:180
static const int JMUONSIMPLEX
Utility class to parse command line options.
then usage $script< input_file >< detector_file > fi set_variable OUTPUT_DIR set_variable SELECTOR JDAQTimesliceL1 set_variable DEBUG case set_variable DEBUG
int getCount(const T &hit)
Get hit count.
Primary particle.
Definition: JHead.hh:1053
Data structure for position in three dimensions.
Definition: JPosition3D.hh:36
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:73
const Trk & get_neutrino(const Evt &evt)
Get incoming neutrino.
The Trk class represents a Monte Carlo (MC) particle as well as a reconstructed track/shower.
Definition: Trk.hh:14
JVector3D & add(const JVector3D &vector)
Add vector.
Definition: JVector3D.hh:142
static const int JMUONENERGY
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
Definition: Evt.hh:19