Jpp
JMuonPostfit.cc
Go to the documentation of this file.
1 #include <string>
2 #include <iostream>
3 #include <iomanip>
4 #include <vector>
5 
6 #include "TROOT.h"
7 #include "TFile.h"
8 #include "TH1D.h"
9 #include "TH2D.h"
10 #include "TProfile.h"
11 #include "TMath.h"
12 
13 #include "evt/Head.hh"
14 #include "evt/Evt.hh"
15 #include "JAAnet/JHead.hh"
16 #include "JAAnet/JHeadToolkit.hh"
17 #include "JDAQ/JDAQEvent.hh"
18 #include "JTools/JConstants.hh"
19 #include "JTools/JQuantile.hh"
20 #include "JTrigger/JTimeConverter.hh"
21 #include "JSupport/JTriggeredFileScanner.hh"
22 #include "JSupport/JMonteCarloFileSupportkit.hh"
23 #include "JSupport/JSupport.hh"
24 #include "JFit/JEvt.hh"
25 #include "JFit/JEvtToolkit.hh"
26 #include "JFit/JFitParameters.hh"
27 #include "JGizmo/JVolume.hh"
28 #include "JLang/JPredicate.hh"
29 
30 #include "Jeep/JParser.hh"
31 #include "Jeep/JMessage.hh"
32 
33 
34 namespace {
35 
36  using KM3NETDAQ::JDAQEvent;
37 
38 
39  /**
40  * Count number of unique identifiers in event.
41  *
42  * \param __begin begin of data
43  * \param __end end of data
44  * \return number of unique identifiers
45  */
46  template<class JHit_t>
47  inline int getCount(JDAQEvent::const_iterator<JHit_t> __begin,
48  JDAQEvent::const_iterator<JHit_t> __end)
49  {
50  using namespace std;
51  using namespace KM3NETDAQ;
52 
53  typedef JDAQModuleIdentifier JType_t;
54 
55  set<JType_t> buffer;
56 
57  for (JDAQEvent::const_iterator<JHit_t> i = __begin; i != __end; ++i) {
58  buffer.insert(JType_t(*i));
59  }
60 
61  return buffer.size();
62  }
63 
64 }
65 /**
66  * \file
67  *
68  * Example program to histogram fit results for the muon reconstruction chain.
69  * \author mdejong, bofearraigh.
70  */
71 int main(int argc, char **argv)
72 {
73  using namespace std;
74  using namespace JPP;
75  using namespace KM3NETDAQ;
76 
77  typedef JTriggeredFileScanner<JEvt> JTriggeredFileScanner_t;
78  typedef JTriggeredFileScanner_t::multi_pointer_type multi_pointer_type;
79 
80  JTriggeredFileScanner_t inputFile;
81  JLimit_t& numberOfEvents = inputFile.getLimit();
82  string outputFile;
83  size_t numberOfPrefits;
84  JQualitySorter quality_sorter;
85  JAtmosphericMuon atmosphere;
86  double Emin_GeV;
87  double NPE;
88  int application;
89  string option;
90  int debug;
91 
92  try {
93 
94  JParser<> zap("Example program to histogram fit results.");
95 
96  zap['f'] = make_field(inputFile);
97  zap['o'] = make_field(outputFile) = "postfit.root";
98  zap['n'] = make_field(numberOfEvents) = JLimit::max();
99  zap['N'] = make_field(numberOfPrefits) = 1;
100  zap['L'] = make_field(quality_sorter) = JPARSER::initialised();
101  zap['E'] = make_field(Emin_GeV) = 0.0;
102  zap['M'] = make_field(NPE) = 0.0;
103  zap['A'] = make_field(application) = JMUONPREFIT, JMUONSIMPLEX, JMUONGANDALF, JMUONENERGY, JMUONSTART, JMUONPATH;
104  zap['a'] = make_field(atmosphere) = JAtmosphericMuon(90.0, 90.0);
105  zap['O'] = make_field(option) = "E", "N", "LINE", "LOGE";
106  zap['d'] = make_field(debug) = 2;
107 
108  zap(argc, argv);
109  }
110  catch(const exception& error) {
111  FATAL(error.what() << endl);
112  }
113 
114  JHead head;
115 
116  try {
117  head = getHeader(inputFile);
118  }
119  catch(const JException& error) {
120  FATAL(error);
121  }
122 
123  const JVolume volume(head, option != "LINE");
124  const JPosition3D center = get<JPosition3D>(head);
125  JCylinder3D cylinder = get<JCylinder3D>(head);
126 
127  cylinder.add(center);
128 
129  NOTICE("Reposition can [m]: " << cylinder << endl);
130 
131  const double EMIN_GEV = 10; // set to low value for ORCA..
132  const double E_nu_min = head.cut_nu.Emin; //neutrino minimum energy
133  const double E_nu_max = head.cut_nu.Emax; //neutrino maximum energy
134 
135  TFile out(outputFile.c_str(), "recreate");
136 
137  TH1D job("job", NULL, 100, 0.5, 100.5);
138 
139  TH1D hn("hn", NULL, 250, -0.5, 249.5); //NDF
140  TH1D hq("hq", NULL, 300, 0.0, 600.0); //quality parameter
141  TH1D hi("hi", NULL, 101, -0.5, 100.5); //index of the best event in order of descending angular resolution
142  TH1D hv("hv", NULL, 200, -6.0, 0.0); //log of Poisson distribution of (#hits,number photo-electrons)
143  TH1D h1("h1", NULL, 200, -2.0, +2.0); //number of photo-electrons along the whole track
144  TH1D hc("hc", NULL, 200, -1.0, +1.0); //dz (z-slope)
145  TH1D hu("hu", NULL, 400, -1.0e3, 1.0e3); //quality difference between best up-going and down-going track (check of atmospheric muon)
146 
147  TH1D hx("hx", NULL, 70, -3.0, +2.3); // angular deviation [log(deg)]
148  TH1D hd("hd", NULL, 100, 0.0, 10.0); // distance between true and best event at intersection point [m]
149  TH1D hz("hz", NULL, 100, -200.0, 200.0); // intersection of track with neutrino vertex [m]
150  TH1D ht("ht", NULL, 100, -100.0, 100.0); // time between true and best track [ns]
151 
152  TH1D e0("e0", NULL, 100, volume.getXmin(), volume.getXmax()); // true track energy [log(E)]
153  TH1D e1("e1", NULL, 100, volume.getXmin(), volume.getXmax()); // uncorrected energy [log(E)]
154  TH1D e2("e2", NULL, 100, volume.getXmin(), volume.getXmax()); // energy of best track[log(E)]
155  TH1D er("er", NULL, 100, -5.0, +5.0); // ratio of reconstructed energy and true energy [log(E/E)]
156  TH2D ee("ee", NULL,
157  40, volume.getXmin(), volume.getXmax(),
158  40, volume.getXmin(), volume.getXmax()); //(Etrue,Ereco)
159 
160  TH2D muon_angle_Emu("muon_angle_Emu", NULL, 30, E_nu_min, E_nu_max, 100, 0,90); // energy vs angular deviation between true muon and reco muon direction
161  TH2D muon_angle_Enu("muon_angle_Enu", NULL, 30, E_nu_min, E_nu_max, 100, 0,90); // energy vs angular deviation between true muon and reco muon direction
162  TH2D neutrino_angle_Enu("neutrino_angle_Enu", NULL, 30, E_nu_min, E_nu_max, 100, 0,90); // energy vs angular deviation between neutrino and reco muon direction
163  TH2D MC_angle_Enu("MC_angle_Enu", NULL, 30, E_nu_min, E_nu_max, 100, 0,90); // energy vs angular deviation between neutrino and true muon direction
164 
165  const Int_t ny = 28800;
166  const Double_t ymin = 0.0; // [deg]
167  const Double_t ymax = 180.0; // [deg]
168 
169  vector<double> X;
170 
171  if (option.find('E') != string::npos) {
172 
173  for (double x = volume.getXmin(); x <= volume.getXmax(); x += (volume.getXmax() - volume.getXmin()) / 20) {
174  X.push_back(x);
175  }
176 
177  } else {
178 
179  double x = -0.5;
180 
181  for ( ; x <= 15.5; x += 1.0) { X.push_back(x); }
182  for ( ; x <= 25.5; x += 2.0) { X.push_back(x); }
183  for ( ; x <= 50.5; x += 5.0) { X.push_back(x); }
184  for ( ; x <= 100.5; x += 10.0) { X.push_back(x); }
185  for ( ; x <= 250.5; x += 20.0) { X.push_back(x); }
186  }
187 
188  TH2D h2("h2", NULL, X.size() - 1, X.data(), ny, ymin, ymax);
189  TProfile he("he", NULL, X.size() - 1, X.data());
190 
191  TH2D ha("ha", NULL, 100, 0.0, 0.5e6, 100, -100.0, +900.0);
192  TH2D hb("hb", NULL, 100, 0.0, 0.5e6, 100, -100.0, +900.0);
193 
194  JQuantile Q("Angle", true);
195  JQuantile O("Omega", true);
196 
197 
198  while (inputFile.hasNext()) {
199 
200  STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);
201 
202  multi_pointer_type ps = inputFile.next();
203 
204  JDAQEvent* tev = ps;
205  JEvt* evt = ps;
206  Evt* event = ps;
207 
208  const JTimeConverter converter(*event, *tev);
209 
210  job.Fill(1.0);
211 
212 
213  double Emu = 0.0;
214  double Emax = 0.0;
215 
216  vector<Trk>::const_iterator muon = event->mc_trks.end();
217 
218  for (vector<Trk>::const_iterator track = event->mc_trks.begin(); track != event->mc_trks.end(); ++track) {
219 
220  if (is_muon(*track)) {
221 
222  Emu += track->E;
223 
224  if (track->E > Emax) {
225  muon = track;
226  Emax = track->E;
227  }
228  }
229  }
230 
231  if (muon == event->mc_trks.end()) {
232  continue;
233  }
234 
235  job.Fill(3.0);
236 
237 
238  // abscissa
239 
240  Double_t x = 0.0;
241 
242  if (option.find('E') != string::npos)
243  x = volume.getX(event->mc_trks[0].E);
244  else
245  x = getCount(tev->begin<JDAQTriggeredHit>(), tev->end<JDAQTriggeredHit>());
246 
247 
248  // weight for efficiency determination
249 
250  Double_t W = 0.0;
251 
252  if (!evt->empty()) {
253 
254  JEvt::iterator __end = partition(evt->begin(), evt->end(), JHistory::is_event(application));
255 
256  if (evt->begin() == __end) {
257  continue;
258  }
259 
260  job.Fill(4.0);
261 
262  if (numberOfPrefits > 0) {
263 
264  JEvt::iterator __q = __end;
265 
266  advance(__end = evt->begin(), min(numberOfPrefits, (size_t) distance(evt->begin(), __q)));
267 
268  partial_sort(evt->begin(), __end, __q, quality_sorter);
269 
270  } else {
271 
272  sort(evt->begin(), __end, quality_sorter);
273  }
274 
275  const JPointing pointing(getDirection(*muon));
276 
277  JEvt::iterator best = pointing(evt->begin(), __end);
278  const Double_t beta = pointing.getAngle(*best);
279  const double Efit = best->getE();
280  const double Eraw = best->getW(JENERGY_ENERGY, numeric_limits<double>::min());
281  const double mip = best->getW(JSTART_NPE_MIP, numeric_limits<double>::max());
282 
283  // selection of fit result
284 
285  bool ok = (Efit >= Emin_GeV &&
286  mip >= NPE);
287 
288  if (ok) {
289 
290  W = 1.0;
291 
292  job.Fill(5.0);
293 
294  hn.Fill((Double_t) best->getNDF());
295  hq.Fill(best->getQ());
296  hi.Fill((Double_t) distance(evt->begin(), best));
297  hc.Fill(best->getDZ());
298 
299  if (( has_neutrino(*event) && get_neutrino(*event).dir.z >= atmosphere.dot2) || //Atmospheric muon hypothesis test. Difference in quality between best upwards and best downgoing track.
300  (!has_neutrino(*event) && best->getDZ() >= atmosphere.dot2)) { //negative values imply a down going track, positive an up-going track.
301  hu.Fill(atmosphere(evt->begin(), __end));
302  }
303 
304  hx.Fill(max(log10(beta), hx.GetXaxis()->GetXmin()));
305 
306  Q.put(beta);
307 
308  const Trk neutrino = get_neutrino(*event);
309  JTrack3E ta = getTrack(*muon);
310  JTrack3E tb = getTrack(*best);
311 
312  ta.add(center);
313  tb.sub(converter.putTime());
314 
315  static_cast<JTrack3D&>(ta).move(ta.getIntersection(tb), getSpeedOfLight()); //move to intersection point of both tracks
316  static_cast<JTrack3D&>(tb).move(tb.getIntersection(ta), getSpeedOfLight());
317 
318  hd.Fill((tb.getPosition() - ta.getPosition()).getLength());
319 
320  if (has_neutrino(*event)) {
321 
322  job.Fill(6.0);
323 
324  JPosition3D vertex = getPosition(neutrino);
325 
326  vertex.add(center);
327 
328  if (cylinder.is_inside(vertex)) {
329 
330  job.Fill(7.0);
331 
332  JTrack3E tc = getTrack(*best);
333 
334  ha.Fill(best->getX()*best->getX() + best->getY()*best->getY(), best->getZ()); //R^{2} [m^2],
335  hz.Fill(tc.getIntersection(vertex));
336  }
337  }
338 
339  if (best->getE() >= EMIN_GEV) {
340  hb.Fill(best->getX()*best->getX() + best->getY()*best->getY(), best->getZ()); //R^{2} [m^2], z [m]
341  }
342  ht.Fill(tb.getT() - ta.getT());
343 
344  e0.Fill(volume.getX(Emu, true));
345  e1.Fill(volume.getX(Eraw, true));
346  e2.Fill(volume.getX(Efit, true));
347  er.Fill(volume.getX(Efit) - volume.getX(Emu));
348  ee.Fill(volume.getX(Emu), volume.getX(Efit));
349  h2.Fill(x, beta);
350 
351  const double Enu = neutrino.E; //true neutrino E
352 
353  neutrino_angle_Enu.Fill(Enu, (180/M_PI)* acos(getDirection(*best).getDot(getDirection(neutrino))));
354  muon_angle_Enu.Fill(Enu,(180/M_PI)* acos(getDirection(*best).getDot(getDirection(*muon))));
355  muon_angle_Emu.Fill(Emu,(180/M_PI)* acos(getDirection(*best).getDot(getDirection(*muon))));
356  MC_angle_Enu.Fill(Enu, (180/M_PI)* acos(getDirection(*muon).getDot(getDirection(neutrino))));
357 
358 
359  if (best->hasW(JSTART_NPE_MIP)) {
360  h1.Fill(log10(best->getW(JSTART_NPE_MIP)));
361  }
362 
363  if (best->hasW(JVETO_NPE) && best->hasW(JVETO_NUMBER_OF_HITS)) {
364 
365  const double npe = best->getW(JVETO_NPE); //number of photoelectrons
366  const int count = best->getW(JVETO_NUMBER_OF_HITS); //number of hits
367  const double pv = TMath::PoissonI(count, npe); //Poisson distribution function for (#hits,npe)
368 
369  hv.Fill(max(log10(pv), hv.GetXaxis()->GetXmin()));
370  }
371  }
372  }
373 
374  he.Fill(x, W); //efficiency
375  }
376  STATUS(endl);
377 
378  NOTICE("Number of events input " << setw(8) << right << job.GetBinContent(1) << endl);
379  NOTICE("Number of events with muon " << setw(8) << right << job.GetBinContent(3) << endl);
380  NOTICE("Number of events with fit " << setw(8) << right << job.GetBinContent(4) << endl);
381  NOTICE("Number of events selected " << setw(8) << right << job.GetBinContent(5) << endl);
382  NOTICE("Number of events with neutrino " << setw(8) << right << job.GetBinContent(6) << endl);
383  NOTICE("Number of events contained " << setw(8) << right << job.GetBinContent(7) << endl);
384 
385  if (Q.getCount() != 0) {
386  NOTICE("Median space angle [deg] " << FIXED (6,3) << Q.getQuantile(0.5) << endl);
387  }
388 
389  out.Write();
390  out.Close();
391 }
JFIT::JMUONSIMPLEX
JSimplex.cc.
Definition: JFitApplications.hh:24
JFIT::JMUONENERGY
JEnergy.cc.
Definition: JFitApplications.hh:26
JPredicate.hh
KM3NETDAQ::JDAQEvent
DAQ Event.
Definition: JDAQEvent.hh:34
FIXED
Auxiliary data structure for floating point format specification.
Definition: JPrint.hh:481
JTOOLS::getSpeedOfLight
const double getSpeedOfLight()
Number of bytes in a gigabyte.
Definition: JConstants.hh:89
KM3NETDAQ::JDAQEvent::begin
const_iterator< T > begin() const
Get begin of data.
JMessage.hh
main
int main(int argc, char **argv)
Definition: JMuonPostfit.cc:71
JFIT::JMUONGANDALF
JGandalf.cc.
Definition: JFitApplications.hh:25
JHead.hh
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: JCounter.hh:35
JPARSER::initialised
Empty structure for specification of parser element that is initialised (i.e.
Definition: JParser.hh:63
JAANET::JHead::cut_nu
JAANET::cut_nu cut_nu
Definition: JHead.hh:1081
JFIT::JMUONSTART
JStart.cc.
Definition: JFitApplications.hh:27
JGEOMETRY3D::JDirection3D::getDot
double getDot(const JAngle3D &angle) const
Get dot product.
Definition: JDirection3D.hh:333
std::vector< double >
JEvtToolkit.hh
JSUPPORT::JLimit_t
JLimit JLimit_t
Type definition of limit.
Definition: JLimit.hh:215
JAANET::getTrack
JTrack3E getTrack(const Trk &track)
Get track.
Definition: JAAnetToolkit.hh:256
JTimeConverter.hh
JPARSER::JParser
Utility class to parse command line options.
Definition: JParser.hh:1493
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
JTriggeredFileScanner.hh
NOTICE
#define NOTICE(A)
Definition: JMessage.hh:64
JQuantile.hh
JFIT::JVETO_NUMBER_OF_HITS
number of hits from JVeto.cc
Definition: JFitParameters.hh:29
JAANET::is_muon
bool is_muon(const Trk &track)
Test whether given track is a (anti-)muon.
Definition: JAAnetToolkit.hh:367
std::set
Definition: JSTDTypes.hh:13
JAANET::JHead
Monte Carlo run header.
Definition: JHead.hh:839
KM3NETDAQ::JDAQEvent::end
const_iterator< T > end() const
Get end of data.
JHeadToolkit.hh
JSupport.hh
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition: JAAnetToolkit.hh:37
KM3NETDAQ::JDAQTriggeredHit
DAQ triggered hit.
Definition: JDAQTriggeredHit.hh:25
JEvt.hh
JSUPPORT::getHeader
Head getHeader(const JMultipleFileScanner_t &file_list)
Get Monte Carlo header.
Definition: JMonteCarloFileSupportkit.hh:425
JFIT::JMUONPREFIT
JPrefit.cc.
Definition: JFitApplications.hh:23
debug
int debug
debug level
Definition: JSirene.cc:59
JConstants.hh
JAANET::get_neutrino
const Trk & get_neutrino(const Evt &evt)
Get incoming neutrino.
Definition: JAAnetToolkit.hh:438
JSUPPORT::JLimit::getLimit
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:73
JAANET::has_neutrino
bool has_neutrino(const Evt &evt)
Test whether given event has an incoming neutrino.
Definition: JAAnetToolkit.hh:427
KM3NETDAQ::JDAQModuleIdentifier
Module identifier.
Definition: JDAQModuleIdentifier.hh:24
JFIT::JENERGY_ENERGY
uncorrected energy [GeV] from JEnergy.cc
Definition: JFitParameters.hh:21
JFIT::JMUONPATH
JPath.cc.
Definition: JFitApplications.hh:41
STATUS
#define STATUS(A)
Definition: JMessage.hh:63
JParser.hh
JFIT::JSTART_NPE_MIP
number of photo-electrons up to the barycentre from JStart.cc
Definition: JFitParameters.hh:25
JAANET::getDirection
JDirection3D getDirection(const Vec &v)
Get direction.
Definition: JAAnetToolkit.hh:221
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1954
JFitParameters.hh
JAANET::getPosition
JPosition3D getPosition(const Vec &v)
Get position.
Definition: JAAnetToolkit.hh:197
DEBUG
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
std
Definition: jaanetDictionary.h:36
KM3NETDAQ
KM3NeT DAQ data structures and auxiliaries.
Definition: DataQueue.cc:39
JVolume.hh
JDAQEvent.hh
JAANET::cut::Emax
double Emax
Maximal energy [GeV].
Definition: JHead.hh:211
JAANET::cut::Emin
double Emin
Minimal energy [GeV].
Definition: JHead.hh:210
JMonteCarloFileSupportkit.hh
FATAL
#define FATAL(A)
Definition: JMessage.hh:67
outputFile
string outputFile
Definition: JDAQTimesliceSelector.cc:37
JFIT::getCount
int getCount(const JHitL0 &hit)
Get hit count.
Definition: JEvtToolkit.hh:728
JFIT::JVETO_NPE
number of photo-electrons from JVeto.cc
Definition: JFitParameters.hh:28
Generated by   doxygen 1.8.16