Jpp
Functions
JShowerPostfit.cc File Reference
#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TProfile.h"
#include "TMath.h"
#include "evt/Head.hh"
#include "evt/Evt.hh"
#include "JAAnet/JHead.hh"
#include "JAAnet/JHeadToolkit.hh"
#include "JDAQ/JDAQEvent.hh"
#include "JTools/JConstants.hh"
#include "JTools/JQuantile.hh"
#include "JTrigger/JTimeConverter.hh"
#include "JSupport/JTriggeredFileScanner.hh"
#include "JSupport/JMonteCarloFileSupportkit.hh"
#include "JSupport/JSupport.hh"
#include "JFit/JEvt.hh"
#include "JFit/JEvtToolkit.hh"
#include "JFit/JFitParameters.hh"
#include "JGizmo/JVolume.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.

Author
mdejong and adomi -> This code is an adaptation of JPostfit to JShowerFit reconstruction

Definition in file JShowerPostfit.cc.

Function Documentation

◆ main()

int main ( int  argc,
char **  argv 
)

Definition at line 73 of file JShowerPostfit.cc.

74 {
75  using namespace std;
76  using namespace JPP;
77  using namespace KM3NETDAQ;
78 
79  typedef JTriggeredFileScanner<JEvt> JTriggeredFileScanner_t;
80  typedef JTriggeredFileScanner_t::multi_pointer_type multi_pointer_type;
81 
82  JTriggeredFileScanner_t inputFile;
83  JLimit_t& numberOfEvents = inputFile.getLimit();
84  string outputFile;
85  size_t numberOfPrefits;
86  JQualitySorter quality_sorter;
87  JAtmosphericMuon atmosphere;
88  double Emin_GeV;
89  double NPE;
90  int application;
91  string option;
92  bool isMuon;
93  int debug;
94 
95  try {
96 
97  JParser<> zap("Example program to histogram fit results.");
98 
99  zap['f'] = make_field(inputFile);
100  zap['o'] = make_field(outputFile) = "postfit.root";
101  zap['n'] = make_field(numberOfEvents) = JLimit::max();
102  zap['N'] = make_field(numberOfPrefits) = 1;
103  zap['L'] = make_field(quality_sorter) = JPARSER::initialised();
104  zap['E'] = make_field(Emin_GeV) = 0.0;
105  zap['M'] = make_field(NPE) = 0.0;
106  zap['A'] = make_field(application) = JSHOWERPREFIT, JSHOWERPOSITIONFIT, JSHOWERCOMPLETEFIT;
107  zap['a'] = make_field(atmosphere) = JAtmosphericMuon(90.0, 90.0);
108  zap['O'] = make_field(option) = "E", "N", "LINE", "LOGE";
109  zap['I'] = make_field(isMuon);
110  zap['d'] = make_field(debug) = 2;
111 
112  zap(argc, argv);
113  }
114  catch(const exception& error) {
115  FATAL(error.what() << endl);
116  }
117 
118  cout << "APPLICATION " << application << endl;
119 
120  JHead head;
121 
122  try {
123  head = getHeader(inputFile);
124  }
125  catch(const JException& error) {
126  FATAL(error);
127  }
128 
129  const JVolume volume(head, option != "LINE");
130  const JPosition3D center = get<JPosition3D>(head);
131  JCylinder3D cylinder = get<JCylinder3D>(head);
132 
133  cylinder.add(center);
134 
135  NOTICE("Reposition can [m]: " << cylinder << endl);
136 
137  const double EMIN_GEV = 10e3; // MUPAGE
138 
139  TFile out(outputFile.c_str(), "recreate");
140 
141  TH1D job("job", NULL, 100, 0.5, 100.5);
142 
143  TH1D hn("hn", NULL, 250, -0.5, 249.5);
144  TH1D hq("hq", NULL, 300, 0.0, 600.0);
145  TH1D hi("hi", NULL, 101, -0.5, 100.5);
146  TH1D hv("hv", NULL, 200, -6.0, 0.0);
147  TH1D h1("h1", NULL, 200, -2.0, +2.0);
148  TH1D hc("hc", NULL, 200, -1.0, +1.0);
149  TH1D hu("hu", NULL, 400, -1.0e3, 1.0e3);
150 
151  TH1D hx("hx", NULL, 100, -3.0, +2.3); // [log(deg)]
152  TH1D hd("hd", NULL, 100, 0.0, 10.0); // [m]
153  TH1D hz("hz", NULL, 100, -200.0, 200.0); // [m]
154  TH1D ht("ht", NULL, 100, -100.0, 100.0); // [ns]
155 
156  TH1D e0("e0", NULL, 100, volume.getXmin(), volume.getXmax()); // [log(E)]
157  TH1D e1("e1", NULL, 100, volume.getXmin(), volume.getXmax()); // [log(E)]
158  TH1D e2("e2", NULL, 100, volume.getXmin(), volume.getXmax()); // [log(E)]
159  TH1D er("er", NULL, 100, -5.0, +5.0); // [log(E/E)]
160  TH2D ee("ee", NULL,
161  40, volume.getXmin(), volume.getXmax(),
162  40, volume.getXmin(), volume.getXmax());
163 
164  TH2D hAngle ("hAngle", ";E_{TRUE};angle [#circ]", 100, volume.getXmin(), volume.getXmax(), 360, 0, 180.0);
165 
166  const Int_t ny = 28800;
167  const Double_t ymin = 0.0; // [deg]
168  const Double_t ymax = 180.0; // [deg]
169 
170  vector<double> X;
171 
172  if (option.find('E') != string::npos) {
173 
174  for (double x = volume.getXmin(); x <= volume.getXmax(); x += (volume.getXmax() - volume.getXmin()) / 20) {
175  X.push_back(x);
176  }
177 
178  } else {
179 
180  double x = -0.5;
181 
182  for ( ; x <= 15.5; x += 1.0) { X.push_back(x); }
183  for ( ; x <= 25.5; x += 2.0) { X.push_back(x); }
184  for ( ; x <= 50.5; x += 5.0) { X.push_back(x); }
185  for ( ; x <= 100.5; x += 10.0) { X.push_back(x); }
186  for ( ; x <= 250.5; x += 20.0) { X.push_back(x); }
187  }
188 
189  TH2D h2("h2", NULL, X.size() - 1, X.data(), ny, ymin, ymax);
190  TProfile he("he", NULL, X.size() - 1, X.data());
191 
192  TH2D ha("ha", NULL, 100, 0.0, 0.5e6, 100, -100.0, +900.0);
193  TH2D hb("hb", NULL, 100, 0.0, 0.5e6, 100, -100.0, +900.0);
194 
195  JQuantile Q("Angle", true);
196  JQuantile O("Omega", true);
197 
198 
199  while (inputFile.hasNext()) {
200 
201  STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);
202 
203  multi_pointer_type ps = inputFile.next();
204 
205  JDAQEvent* tev = ps;
206  JEvt* evt = ps;
207  Evt* event = ps;
208 
209  const JTimeConverter converter(*event, *tev);
210 
211  job.Fill(1.0);
212 
213  double Emu = 0.0;
214  double Emax = 0.0;
215 
216  vector<Trk>::const_iterator lepton = event->mc_trks.end(); // can be electron or muon
217  // vector<Trk>::const_iterator electron = event->mc_trks.end();
218 
219  for (vector<Trk>::const_iterator shower = event->mc_trks.begin(); shower != event->mc_trks.end(); ++shower) {
220 
221  if (!isMuon && is_electron(*shower)) {
222 
223  Emu += shower->E;
224 
225  if (shower->E > Emax) {
226  lepton = shower;
227  Emax = shower->E;
228  }
229  } else if(isMuon && is_muon(*shower)){
230 
231  Emu += shower->E;
232  if (shower->E > Emax) {
233  lepton = shower;
234  Emax = shower->E;
235  }
236  }
237  }
238 
239  if (lepton == event->mc_trks.end()) {
240  continue;
241  }
242 
243  job.Fill(3.0);
244 
245 
246  // abscissa
247 
248  Double_t x = 0.0;
249 
250  if (option.find('E') != string::npos)
251  x = volume.getX(event->mc_trks[0].E);
252  else
253  x = getCount(tev->begin<JDAQTriggeredHit>(), tev->end<JDAQTriggeredHit>());
254 
255 
256  // weight for efficiency determination
257 
258  Double_t W = 0.0;
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  const JPointing pointing(getDirection(*lepton));
284 
285  JEvt::iterator best = pointing(evt->begin(), __end);
286  const Double_t beta = pointing.getAngle(*best);
287  const double Efit = best->getE();
288  const double Eraw = best->getW(JENERGY_ENERGY, numeric_limits<double>::min());
289  const double mip = best->getW(JSTART_NPE_MIP, numeric_limits<double>::max());
290  //const double npe = best->getW(JVETO_NPE, 0.0);
291  //const int count = best->getW(JVETO_NUMBER_OF_HITS, 0.0);
292 
293  // selection of fit result
294 
295  bool ok = (Efit >= Emin_GeV &&
296  mip >= NPE);
297 
298  if (ok) {
299 
300  W = 1.0;
301 
302  job.Fill(5.0);
303 
304  hn.Fill((Double_t) best->getNDF());
305  hq.Fill(best->getQ());
306  hi.Fill((Double_t) distance(evt->begin(), best));
307  hc.Fill(best->getDZ());
308 
309  if (( has_neutrino(*event) && get_neutrino(*event).dir.z >= atmosphere.dot2) ||
310  (!has_neutrino(*event) && best->getDZ() >= atmosphere.dot2)) {
311  hu.Fill(atmosphere(evt->begin(), __end));
312  }
313 
314  hx.Fill(max(log10(beta), hx.GetXaxis()->GetXmin()));
315 
316  Q.put(beta);
317 
318  JTrack3E ta = getTrack(*lepton);
319  JTrack3E tb = getTrack(*best);
320 
321  ta.add(center);
322  tb.sub(converter.putTime());
323 
324  //ta.move(ta.getIntersection(tb), getSpeedOfLight(), gWater);
325  //tb.move(tb.getIntersection(ta), getSpeedOfLight(), gWater);
326  static_cast<JTrack3D&>(ta).move(ta.getIntersection(tb), getSpeedOfLight());
327  static_cast<JTrack3D&>(tb).move(tb.getIntersection(ta), getSpeedOfLight());
328 
329  hd.Fill((tb.getPosition() - ta.getPosition()).getLength());
330 
331  if (has_neutrino(*event)) {
332 
333  job.Fill(6.0);
334 
335  JPosition3D vertex = getPosition(get_neutrino(*event));
336 
337  vertex.add(center);
338 
339  if (cylinder.is_inside(vertex)) {
340 
341  job.Fill(7.0);
342 
343  JTrack3E tc = getTrack(*best);
344 
345  ha.Fill(best->getX()*best->getX() + best->getY()*best->getY(), best->getZ());
346  hz.Fill(tc.getIntersection(vertex));
347  }
348  }
349 
350  if (best->getE() >= EMIN_GEV) {
351  hb.Fill(best->getX()*best->getX() + best->getY()*best->getY(), best->getZ());
352  }
353 
354  ht.Fill(tb.getT() - ta.getT());
355 
356  e0.Fill(volume.getX(Emu, true));
357  e1.Fill(volume.getX(Eraw, true));
358  e2.Fill(volume.getX(Efit, true));
359  er.Fill(volume.getX(Efit) - volume.getX(Emu));
360  ee.Fill(volume.getX(Emu), volume.getX(Efit));
361 
362  hAngle.Fill(Efit, beta);
363  h2.Fill(x, beta);
364 
365  if (best->hasW(JSTART_NPE_MIP)) {
366  h1.Fill(log10(best->getW(JSTART_NPE_MIP)));
367  }
368 
369  if (best->hasW(JVETO_NPE) && best->hasW(JVETO_NUMBER_OF_HITS)) {
370 
371  const double npe = best->getW(JVETO_NPE);
372  const int count = best->getW(JVETO_NUMBER_OF_HITS);
373  const double pv = TMath::PoissonI(count, npe);
374 
375  hv.Fill(max(log10(pv), hv.GetXaxis()->GetXmin()));
376  }
377  }
378  }
379 
380  he.Fill(x, W);
381  }
382  STATUS(endl);
383 
384  NOTICE("Number of events input " << setw(8) << right << job.GetBinContent(1) << endl);
385  if(!isMuon){
386  NOTICE("Number of events with electron " << setw(8) << right << job.GetBinContent(3) << endl);
387  } else{
388  NOTICE("Number of events with muon " << setw(8) << right << job.GetBinContent(3) << endl);
389  }
390  NOTICE("Number of events with fit " << setw(8) << right << job.GetBinContent(4) << endl);
391  NOTICE("Number of events selected " << setw(8) << right << job.GetBinContent(5) << endl);
392  NOTICE("Number of events with neutrino " << setw(8) << right << job.GetBinContent(6) << endl);
393  NOTICE("Number of events contained " << setw(8) << right << job.GetBinContent(7) << endl);
394 
395  if (Q.getCount() != 0) {
396  NOTICE("Median space angle [deg] " << FIXED (6,3) << Q.getQuantile(0.5) << endl);
397  }
398 
399  out.Write();
400  out.Close();
401 }
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.
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
JFIT::JSHOWERPOSITIONFIT
JShowerPositionFit.cc.
Definition: JFitApplications.hh:33
std::vector< double >
JFIT::JSHOWERPREFIT
JShowerPrefit.cc.
Definition: JFitApplications.hh:32
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
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
NOTICE
#define NOTICE(A)
Definition: JMessage.hh:64
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
JAANET::JHead
Monte Carlo run header.
Definition: JHead.hh:839
KM3NETDAQ::JDAQEvent::end
const_iterator< T > end() const
Get end of data.
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
JSUPPORT::getHeader
Head getHeader(const JMultipleFileScanner_t &file_list)
Get Monte Carlo header.
Definition: JMonteCarloFileSupportkit.hh:425
debug
int debug
debug level
Definition: JSirene.cc:59
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
JFIT::JENERGY_ENERGY
uncorrected energy [GeV] from JEnergy.cc
Definition: JFitParameters.hh:21
STATUS
#define STATUS(A)
Definition: JMessage.hh:63
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
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
JAANET::is_electron
bool is_electron(const Trk &track)
Test whether given track is a (anti-)electron.
Definition: JAAnetToolkit.hh:359
KM3NETDAQ
KM3NeT DAQ data structures and auxiliaries.
Definition: DataQueue.cc:39
JFIT::JSHOWERCOMPLETEFIT
JShowerFit.cc.
Definition: JFitApplications.hh:34
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