Jpp test-rotations-new
the software that should make you happy
Loading...
Searching...
No Matches
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"
23#include "JAAnet/JVolume.hh"
24#include "JDAQ/JDAQEventIO.hh"
26#include "JTools/JQuantile.hh"
29#include "JSupport/JSupport.hh"
32#include "JLang/JPredicate.hh"
33
34#include "Jeep/JParser.hh"
35#include "Jeep/JMessage.hh"
36
37
38namespace {
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>
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 */
79int 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) = "LINE", "LOGE", "LINN", "LOGN";
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 = getCommonHeader(inputFile);
128 }
129 catch(const exception& error) {
130 FATAL(error.what() << endl);
131 }
132
133 JVolume volume(head, option != "LINE");
134 JPosition3D offset = getPosition(getOffset(head));
135 JPosition3D origin = getPosition(getOrigin(head));
136 JCylinder3D cylinder = getCylinder(head);
137
138 cylinder.add(offset);
139
140 NOTICE("Offset: " << offset << endl);
141 NOTICE("Cylinder: " << cylinder << endl);
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 == "LINN") {
177
178 const double xmin = log((double) 3);
179 const double xmax = log((double) 300);
180
181 for (double x = xmin, dx = (xmax - xmin) / 20; x <= xmax; x += dx) {
182 X.push_back((int) exp(x));
183 }
184
185 } else if (option == "LOGN") {
186
187 const double xmin = log10((double) 3);
188 const double xmax = log10((double) 300);
189
190 for (double x = xmin, dx = (xmax - xmin) / 20; x <= xmax; x += dx) {
191 X.push_back(x);
192 }
193
194 } else {
195
196 for (double x = volume.getXmin(), dx = (volume.getXmax() - volume.getXmin()) / 20.0; x < volume.getXmax() + 0.5*dx; x += dx) {
197 X.push_back(x);
198 }
199 }
200
201 TH2D h2("h2", NULL, X.size() - 1, X.data(), ny, ymin, ymax);
202
203 TH2D Va("Va", NULL,
204 100, 0.0, cylinder.getRadius()*cylinder.getRadius(),
205 100, cylinder.getZmin() - 50.0, cylinder.getZmax() + 50.0);
206 TH2D Vb("Vb", NULL,
207 100, 0.0, cylinder.getRadius()*cylinder.getRadius(),
208 100, cylinder.getZmin() - 50.0, cylinder.getZmax() + 50.0);
209
210 JQuantile Q("Angle", true);
211 JQuantile O("Omega", true);
212
213
214 while (inputFile.hasNext()) {
215
216 STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);
217
218 multi_pointer_type ps = inputFile.next();
219
220 JDAQEvent* tev = ps;
221 JEvt* evt = ps;
222 Evt* event = ps;
223
224 const time_converter converter(*event, *tev);
225
226 job.Fill(1.0);
227
228
229 double Enu = 0.0;
230 double Emu = 0.0;
231 double Emax = 0.0;
232
233 if (has_neutrino(*event)) {
234 Enu = get_neutrino(*event).E;
235 }
236
237 vector<Trk>::const_iterator muon = event->mc_trks.end();
238
239 for (vector<Trk>::const_iterator track = event->mc_trks.begin(); track != event->mc_trks.end(); ++track) {
240
241 if (is_muon(*track)) {
242
243 Emu += track->E;
244
245 if (track->E > Emax) {
246 muon = track;
247 Emax = track->E;
248 }
249 }
250 }
251
252 if (muon == event->mc_trks.end()) {
253 continue;
254 }
255
256 job.Fill(3.0);
257
258
259 // abscissa
260
261 Double_t x = 0.0;
262
263 {
264 const double E = event->mc_trks[0].E;
265 const int N = getCount(tev->begin<JDAQTriggeredHit>(), tev->end<JDAQTriggeredHit>());
266
267 if (option == "LINN")
268 x = N;
269 else if (option == "LOGN")
270 x = log10(N);
271 else
272 x = volume.getX(E);
273 }
274
275
276 if (!evt->empty()) {
277
278 JEvt::iterator __end = partition(evt->begin(), evt->end(), JHistory::is_application(application));
279
280 if (evt->begin() == __end) {
281 continue;
282 }
283
284 job.Fill(4.0);
285
286 if (numberOfPrefits > 0) {
287
288 JEvt::iterator __q = __end;
289
290 advance(__end = evt->begin(), min(numberOfPrefits, (size_t) distance(evt->begin(), __q)));
291
292 partial_sort(evt->begin(), __end, __q, quality_sorter);
293
294 } else {
295
296 sort(evt->begin(), __end, quality_sorter);
297 }
298
299 double E = 0.0;
300 JPointing pointing;
301
302 if (primary == muon_t) {
303 E = Emu;
304 pointing = JPointing(getDirection(*muon));
305 } else if (primary == neutrino_t) {
306 E = Enu;
307 pointing = JPointing(getDirection(get_neutrino(*event)));
308 }
309
310 JEvt::iterator best = pointing(evt->begin(), __end);
311 const Double_t beta = pointing.getAngle(*best);
312 const double Efit = best->getE();
313 const double Eraw = best->getW(JENERGY_ENERGY, numeric_limits<double>::min());
314 const double mip = best->getW(JSTART_NPE_MIP, numeric_limits<double>::max());
315
316 // selection of fit result
317
318 bool ok = (Efit >= Emin_GeV && mip >= NPE);
319
320 if (ok) {
321
322 job.Fill(5.0);
323
324 hn.Fill((Double_t) best->getNDF());
325 hq.Fill(best->getQ());
326 hi.Fill((Double_t) distance(evt->begin(), best));
327 hc.Fill(best->getDZ());
328
329 // atmospheric muon hypothesis test.
330
331 if (( has_neutrino(*event) && get_neutrino(*event).dir.z >= atmosphere.dot2) || // difference in quality between best upwards and best downgoing track.
332 (!has_neutrino(*event) && best->getDZ() >= atmosphere.dot2)) { // negative (positive) values imply a down(up) going track.
333 hu.Fill(atmosphere(evt->begin(), __end));
334 }
335
336 hx.Fill(max(log10(beta), hx.GetXaxis()->GetXmin()));
337
338 Q.put(beta);
339
340 JTrack3E ta = getTrack(*muon);
341 JTrack3E tb = getTrack(*best);
342
343 ta.add(offset);
344 tb.sub(converter.putTime());
345
346 static_cast<JTrack3D&>(ta).move(ta.getIntersection(tb), getSpeedOfLight()); // move to intersection point of both tracks
347 static_cast<JTrack3D&>(tb).move(tb.getIntersection(ta), getSpeedOfLight());
348
349 hd.Fill((tb.getPosition() - ta.getPosition()).getLength());
350
351 if (has_neutrino(*event)) {
352
353 job.Fill(6.0);
354
355 const Trk neutrino = get_neutrino(*event);
356
357 JPosition3D vertex = getPosition(neutrino);
358
359 vertex.add(offset);
360
361 if (cylinder.is_inside(vertex)) {
362
363 job.Fill(7.0);
364
365 Va.Fill(JVector2D(best->getX() - origin.getX(), best->getY() - origin.getY()).getLengthSquared(), best->getZ()); //R^{2} [m^2], z [m]
366
367 JTrack3E tc = getTrack(*best);
368
369 hz0.Fill(tc.getIntersection(vertex));
370 hz1.Fill(best->getW(JSTART_LENGTH_METRES, 0.0) - gWater(muon->E));
371 }
372 }
373
374 Vb.Fill(JVector2D(best->getX() - origin.getX(), best->getY() - origin.getY()).getLengthSquared(), best->getZ()); //R^{2} [m^2], z [m]
375 ht.Fill(tb.getT() - ta.getT());
376
377 E_0.Fill(volume.getX(E, true));
378 E_1.Fill(volume.getX(Eraw, true));
379 E_2.Fill(volume.getX(Efit, true));
380 E_E.Fill(volume.getX(Efit) - volume.getX(E));
381 ExE.Fill(volume.getX(E), volume.getX(Efit));
382
383 h2.Fill(x, beta);
384
385 if (best->hasW(JSTART_NPE_MIP)) {
386 h1.Fill(log10(best->getW(JSTART_NPE_MIP)));
387 }
388
389 if (best->hasW(JVETO_NPE) && best->hasW(JVETO_NUMBER_OF_HITS)) {
390
391 const double npe = best->getW(JVETO_NPE); // number of photo-electrons
392 const int count = best->getW(JVETO_NUMBER_OF_HITS); // number of hits
393 const double pv = TMath::PoissonI(count, npe); // Poisson distribution function for (#hits,npe)
394
395 hv.Fill(max(log10(pv), hv.GetXaxis()->GetXmin()));
396 }
397 }
398 }
399 }
400 STATUS(endl);
401
402 NOTICE("Number of events input " << setw(8) << right << job.GetBinContent(1) << endl);
403 NOTICE("Number of events with muon " << setw(8) << right << job.GetBinContent(3) << endl);
404 NOTICE("Number of events with fit " << setw(8) << right << job.GetBinContent(4) << endl);
405 NOTICE("Number of events selected " << setw(8) << right << job.GetBinContent(5) << endl);
406 NOTICE("Number of events with neutrino " << setw(8) << right << job.GetBinContent(6) << endl);
407 NOTICE("Number of events contained " << setw(8) << right << job.GetBinContent(7) << endl);
408
409 if (Q.getCount() != 0) {
410 NOTICE("Median space angle [deg] " << FIXED (6,3) << Q.getQuantile(0.5) << endl);
411 }
412
413 out.Write();
414 out.Close();
415}
Definition of hit and track types and auxiliary methods for handling Monte Carlo data.
string outputFile
General purpose messaging.
#define DEBUG(A)
Message macros.
Definition JMessage.hh:62
#define STATUS(A)
Definition JMessage.hh:63
#define NOTICE(A)
Definition JMessage.hh:64
#define FATAL(A)
Definition JMessage.hh:67
int debug
debug level
Definition JSirene.cc:72
int main(int argc, char **argv)
Utility class to parse command line options.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2142
Physics constants.
ROOT TTree parameter settings of various packages.
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.
Monte Carlo run header.
Definition JHead.hh:1236
double getRadius() const
Get radius.
Definition JCircle2D.hh:144
Data structure for vector in two dimensions.
Definition JVector2D.hh:34
double getLengthSquared() const
Get length squared.
Definition JVector2D.hh:188
double getIntersection(const JVector3D &pos) const
Get longitudinal position along axis of position of closest approach with given position.
Definition JAxis3D.hh:146
double getZmin() const
Get minimal z position.
bool is_inside(const JVector3D &pos) const
Check whether given point is inside cylinder.
double getZmax() const
Get maximal z position.
JCylinder3D & add(const JVector3D &pos)
Add position.
Data structure for position in three dimensions.
const JPosition3D & getPosition() const
Get position.
double getT(const JVector3D &pos) const
Get arrival time of Cherenkov light at given position.
Definition JTrack3D.hh:126
JTime & add(const JTime &value)
Addition operator.
JTime & sub(const JTime &value)
Subtraction operator.
3D track with energy.
Definition JTrack3E.hh:32
JVector3D & add(const JVector3D &vector)
Add vector.
Definition JVector3D.hh:142
Utility class to parse command line options.
Definition JParser.hh:1698
Auxiliary class to evaluate atmospheric muon hypothesis.
Auxiliary class to compare fit results with respect to a reference direction (e.g....
double getAngle(const JFit &fit) const
Get angle between reference direction and fit result.
Template const_iterator.
Definition JDAQEvent.hh:68
const_iterator< T > end() const
Get end of data.
const_iterator< T > begin() const
Get begin of data.
Auxiliary class to convert DAQ hit time to/from Monte Carlo hit time.
double putTime() const
Get Monte Carlo time minus DAQ/trigger time.
static const int JMUONGANDALF
static const int JMUONPREFIT
static const int JMUONENERGY
static const int JMUONSIMPLEX
static const int JMUONPATH
static const int JMUONSTART
static const int JENERGY_ENERGY
uncorrected energy [GeV] from JEnergy.cc
static const int JVETO_NPE
number of photo-electrons from JVeto.cc
static const int JSTART_LENGTH_METRES
distance between projected positions on the track of optical modules for which the response does not ...
static const int JVETO_NUMBER_OF_HITS
number of hits from JVeto.cc
static const int JSTART_NPE_MIP
number of photo-electrons up to the barycentre from JStart.cc
Vec getOrigin(const JHead &header)
Get origin.
JDirection3D getDirection(const Vec &dir)
Get direction.
JTrack3E getTrack(const Trk &track)
Get track.
JCylinder3D getCylinder(const JHead &header)
Get cylinder corresponding to the positions of generated tracks (i.e.
bool has_neutrino(const Evt &evt)
Test whether given event has an incoming neutrino.
JPosition3D getPosition(const Vec &pos)
Get position.
bool is_muon(const Trk &track)
Test whether given track is a (anti-)muon.
Vec getOffset(const JHead &header)
Get offset.
const Trk & get_neutrino(const Evt &evt)
Get incoming neutrino.
size_t getCount(const array_type< T > &buffer, const JCompare_t &compare)
Count number of unique values.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
counter_type advance(counter_type &counter, const counter_type value, const counter_type limit=std::numeric_limits< counter_type >::max())
Advance counter.
Head getCommonHeader(const JMultipleFileScanner_t &file_list)
Get common Monte Carlo header.
KM3NeT DAQ data structures and auxiliaries.
Definition DataQueue.cc:39
const char *const neutrino_t
Definition io_ascii.hh:29
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
Definition Evt.hh:21
Auxiliary data structure for floating point format specification.
Definition JManip.hh:448
Auxiliary class for histogramming of effective volume.
Definition JVolume.hh:29
Double_t getX(const Double_t E, double constrain=false) const
Get abscissa value.
Definition JVolume.hh:115
Double_t getXmax() const
Get maximal abscissa value.
Definition JVolume.hh:102
Double_t getXmin() const
Get minimal abscissa value.
Definition JVolume.hh:91
Primary particle.
Definition JHead.hh:1174
Acoustic event fit.
Auxiliary class to test history.
Definition JHistory.hh:167
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition JParser.hh:68
General purpose sorter of fit results.
Auxiliary class for defining the range of iterations of objects.
Definition JLimit.hh:45
const JLimit & getLimit() const
Get limit.
Definition JLimit.hh:84
static counter_type max()
Get maximum counter value.
Definition JLimit.hh:128
Auxiliary class to synchronously read DAQ events and Monte Carlo events (and optionally other events)...
Auxiliary data structure for running average, standard deviation and quantiles.
Definition JQuantile.hh:46
double getQuantile(const double Q, const Quantile_t option=forward_t) const
Get quantile.
Definition JQuantile.hh:349
void put(const double x, const double w=1.0)
Put value.
Definition JQuantile.hh:133
long long int getCount() const
Get total count.
Definition JQuantile.hh:186
The Trk class represents a Monte Carlo (MC) particle as well as a reconstructed track/shower.
Definition Trk.hh:15
Vec dir
track direction
Definition Trk.hh:18
double E
Energy [GeV] (either MC truth or reconstructed)
Definition Trk.hh:20
double z
Definition Vec.hh:14
Reconstruction type dependent comparison of track quality.
Auxiliary include file for time conversion between DAQ/trigger hit and Monte Carlo hit.