76 using namespace KM3NETDAQ;
78 typedef JTriggeredFileScanner<JEvt> JTriggeredFileScanner_t;
79 typedef JTriggeredFileScanner_t::multi_pointer_type multi_pointer_type;
81 JTriggeredFileScanner_t inputFile;
84 size_t numberOfPrefits;
85 JQualitySorter quality_sorter;
86 JAtmosphericMuon atmosphere;
95 JParser<> zap(
"Example program to histogram fit results.");
99 zap[
'n'] =
make_field(numberOfEvents) = JLimit::max();
105 zap[
'a'] =
make_field(atmosphere) = JAtmosphericMuon(90.0, 90.0);
106 zap[
'O'] =
make_field(option) =
"E",
"N",
"LINE",
"LOGE";
111 catch(
const exception& error) {
112 FATAL(error.what() << endl);
123 catch(
const JException& error) {
127 const JVolume volume(head, option !=
"LINE");
128 const JPosition3D center = get<JPosition3D>(head);
129 JCylinder3D cylinder = get<JCylinder3D>(head);
131 cylinder.add(center);
133 NOTICE(
"Reposition can [m]: " << cylinder << endl);
135 const double EMIN_GEV = 10e3;
140 TH1D job(
"job", NULL, 100, 0.5, 100.5);
142 TH1D hn(
"hn", NULL, 250, -0.5, 249.5);
143 TH1D hq(
"hq", NULL, 300, 0.0, 600.0);
144 TH1D hi(
"hi", NULL, 101, -0.5, 100.5);
145 TH1D hv(
"hv", NULL, 200, -6.0, 0.0);
146 TH1D h1(
"h1", NULL, 200, -2.0, +2.0);
147 TH1D hc(
"hc", NULL, 200, -1.0, +1.0);
148 TH1D hu(
"hu", NULL, 400, -1.0e3, 1.0e3);
150 TH1D hx(
"hx", NULL, 100, -3.0, +2.3);
151 TH1D hd(
"hd", NULL, 100, 0.0, 10.0);
152 TH1D hz(
"hz", NULL, 100, -200.0, 200.0);
153 TH1D ht(
"ht", NULL, 100, -100.0, 100.0);
155 TH1D e0(
"e0", NULL, 100, volume.getXmin(), volume.getXmax());
156 TH1D e1(
"e1", NULL, 100, volume.getXmin(), volume.getXmax());
157 TH1D e2(
"e2", NULL, 100, volume.getXmin(), volume.getXmax());
158 TH1D er(
"er", NULL, 100, -5.0, +5.0);
160 40, volume.getXmin(), volume.getXmax(),
161 40, volume.getXmin(), volume.getXmax());
163 const Int_t ny = 28800;
164 const Double_t ymin = 0.0;
165 const Double_t ymax = 180.0;
169 if (option.find(
'E') != string::npos) {
171 for (
double x = volume.getXmin(); x <= volume.getXmax(); x += (volume.getXmax() - volume.getXmin()) / 20) {
179 for ( ; x <= 15.5; x += 1.0) { X.push_back(x); }
180 for ( ; x <= 25.5; x += 2.0) { X.push_back(x); }
181 for ( ; x <= 50.5; x += 5.0) { X.push_back(x); }
182 for ( ; x <= 100.5; x += 10.0) { X.push_back(x); }
183 for ( ; x <= 250.5; x += 20.0) { X.push_back(x); }
186 TH2D h2(
"h2", NULL, X.size() - 1, X.data(), ny, ymin, ymax);
187 TProfile he(
"he", NULL, X.size() - 1, X.data());
189 TH2D ha(
"ha", NULL, 100, 0.0, 0.5e6, 100, -100.0, +900.0);
190 TH2D hb(
"hb", NULL, 100, 0.0, 0.5e6, 100, -100.0, +900.0);
192 JQuantile Q(
"Angle",
true);
193 JQuantile O(
"Omega",
true);
196 while (inputFile.hasNext()) {
198 STATUS(
"event: " << setw(10) << inputFile.getCounter() <<
'\r');
DEBUG(endl);
200 multi_pointer_type ps = inputFile.next();
206 const JTimeConverter converter(*
event, *tev);
222 if (track->E > Emax) {
229 if (muon ==
event->mc_trks.end()) {
240 if (option.find(
'E') != string::npos)
241 x = volume.getX(
event->mc_trks[0].E);
252 JEvt::iterator __end = partition(evt->begin(), evt->end(), JHistory::is_event(application));
254 if (evt->begin() == __end) {
260 if (numberOfPrefits > 0) {
262 JEvt::iterator __q = __end;
264 advance(__end = evt->begin(), min(numberOfPrefits, (
size_t) distance(evt->begin(), __q)));
266 partial_sort(evt->begin(), __end, __q, quality_sorter);
270 sort(evt->begin(), __end, quality_sorter);
275 JEvt::iterator best = pointing(evt->begin(), __end);
276 const Double_t beta = pointing.getAngle(*best);
277 const double Efit = best->getE();
278 const double Eraw = best->getW(
JENERGY_ENERGY, numeric_limits<double>::min());
279 const double mip = best->getW(
JSTART_NPE_MIP, numeric_limits<double>::max());
285 bool ok = (Efit >= Emin_GeV &&
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());
301 hu.Fill(atmosphere(evt->begin(), __end));
304 hx.Fill(max(log10(beta), hx.GetXaxis()->GetXmin()));
312 tb.sub(converter.putTime());
316 static_cast<JTrack3D&
>(ta).move(ta.getIntersection(tb),
getSpeedOfLight());
317 static_cast<JTrack3D&
>(tb).move(tb.getIntersection(ta),
getSpeedOfLight());
319 hd.Fill((tb.getPosition() - ta.getPosition()).getLength());
329 if (cylinder.is_inside(vertex)) {
335 ha.Fill(best->getX()*best->getX() + best->getY()*best->getY(), best->getZ());
336 hz.Fill(tc.getIntersection(vertex));
340 if (best->getE() >= EMIN_GEV) {
341 hb.Fill(best->getX()*best->getX() + best->getY()*best->getY(), best->getZ());
344 ht.Fill(tb.getT() - ta.getT());
346 e0.Fill(volume.getX(Emu,
true));
347 e1.Fill(volume.getX(Eraw,
true));
348 e2.Fill(volume.getX(Efit,
true));
349 er.Fill(volume.getX(Efit) - volume.getX(Emu));
350 ee.Fill(volume.getX(Emu), volume.getX(Efit));
360 const double npe = best->getW(
JVETO_NPE);
362 const double pv = TMath::PoissonI(count, npe);
364 hv.Fill(max(log10(pv), hv.GetXaxis()->GetXmin()));
373 NOTICE(
"Number of events input " << setw(8) << right << job.GetBinContent(1) << endl);
374 NOTICE(
"Number of events with muon " << setw(8) << right << job.GetBinContent(3) << endl);
375 NOTICE(
"Number of events with fit " << setw(8) << right << job.GetBinContent(4) << endl);
376 NOTICE(
"Number of events selected " << setw(8) << right << job.GetBinContent(5) << endl);
377 NOTICE(
"Number of events with neutrino " << setw(8) << right << job.GetBinContent(6) << endl);
378 NOTICE(
"Number of events contained " << setw(8) << right << job.GetBinContent(7) << endl);
380 if (Q.getCount() != 0) {
381 NOTICE(
"Median space angle [deg] " <<
FIXED (6,3) << Q.getQuantile(0.5) << endl);
Utility class to parse command line options.
number of photo-electrons from JVeto.cc
JTrack3E getTrack(const Trk &track)
Get track.
double getCount(TH1D *hptr, int muon_threshold)
bool has_neutrino(const Evt &evt)
Test whether given event has an incoming neutrino.
bool is_muon(const Trk &track)
Test whether given track is a (anti-)muon.
Structure to store the ToT mean and standard deviation of the hits produced by a nanobeacon in a sour...
Empty structure for specification of parser element that is initialised (i.e.
Auxiliary data structure for floating point format specification.
number of hits from JVeto.cc
Head getHeader(const JMultipleFileScanner_t &file_list)
Get Monte Carlo header.
JLimit JLimit_t
Type definition of limit.
const_iterator< T > end() const
Get end of data.
const_iterator< T > begin() const
Get begin of data.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
counter_type advance(counter_type &counter, const counter_type value, const counter_type limit=std::numeric_limits< counter_type >::max())
Advance counter.
uncorrected energy [GeV] from JEnergy.cc
number of photo-electrons up to the barycentre from JStart.cc
JDirection3D getDirection(const Vec &v)
Get direction.
const JLimit & getLimit() const
Get limit.
const Trk & get_neutrino(const Evt &evt)
Get incoming neutrino.
#define DEBUG(A)
Message macros.
JPosition3D getPosition(const Vec &v)
Get position.