9#include "TApplication.h"
104 return (i >= 0 && i < (
int) trk.
fitinf.size());
115 inline double getW(
const Trk& trk,
const int i)
129 inline double getW(
const Trk& trk,
const int i,
const double value)
145 void setW(
Trk& trk,
const int i,
const double value)
147 if (i >= (
int) trk.
fitinf.size()) {
148 trk.
fitinf.resize(i + 1, 0.0);
168 inline void execute(
const std::string& command,
int debug)
175 istream in(process.getInputStreamBuffer());
177 for (
string buffer;
getline(in, buffer); ) {
178 DEBUG(buffer << endl);
182 const char*
const histogram_t =
"histogram";
183 const char*
const arrow_t =
"arrow";
201 JLimit_t& numberOfEvents = inputFile.getLimit();
210 double arrowSize = 0.003;
211 string arrowType =
"|->";
212 double arrowScale = 250.0;
213 Width_t lineWidth = 2;
214 Style_t lineStyle = 1;
236 JParser<> zap(
"Program to display hit probabilities.");
238 zap[
'w'] =
make_field(canvas,
"size of canvas <nx>x<ny> [pixels]") =
JCanvas(1200, 600);
239 zap[
'f'] =
make_field(inputFile,
"input file (output of JXXXMuonReconstruction.sh)");
247 zap[
'O'] =
make_field(option,
"draw option") = arrow_t, histogram_t;
248 zap[
'B'] =
make_field(batch,
"batch processing");
253 catch(
const exception& error) {
254 FATAL(error.what() << endl);
258 FATAL(
"Missing output file name " <<
outputFile <<
" in batch mode." << endl);
283 gROOT->SetBatch(batch);
285 TApplication* tp =
new TApplication(
"user", NULL, NULL);
286 TCanvas* cv =
new TCanvas(
"display",
"", canvas.
x, canvas.
y);
288 unique_ptr<TStyle> gStyle(
new JStyle(
"gplot", cv->GetWw(), cv->GetWh(), graphics));
290 gROOT->SetStyle(
"gplot");
293 const size_t NUMBER_OF_PADS = 3;
295 cv->SetFillStyle(4000);
296 cv->SetFillColor(kWhite);
298 TPad*
p1 =
new TPad(
"p1", NULL, 0.0, 0.00, 1.0, 0.95);
299 TPad* p2 =
new TPad(
"p2", NULL, 0.0, 0.95, 1.0, 1.00);
301 p1->Divide(NUMBER_OF_PADS, 1);
306 const double Dmax = 1000.0;
307 const double Rmin = 0.0;
308 const double Rmax = min(parameters.
roadWidth_m, 0.4 * Dmax);
309 const double Tmin = min(parameters.
TMin_ns, -10.0);
310 const double Tmax = max(parameters.
TMax_ns, +100.0);
311 const double Amin = 0.002 * (Tmax - Tmin);
312 const double Amax = 0.8 * (Tmax - Tmin);
313 const double ymin = Tmin - (option == arrow_t ? 0.2 * Amax : 0.0);
314 const double ymax = Tmax + (option == arrow_t ? 0.5 * Amax : 0.0);
316 const string Xlabel[NUMBER_OF_PADS] = {
"R [m]",
"#phi [rad]",
"z [m]" };
317 const double Xmin [NUMBER_OF_PADS] = { Rmin, -
PI, -0.3 * Dmax };
318 const double Xmax [NUMBER_OF_PADS] = { Rmax, +
PI, +0.3 * Dmax };
320 double Xs[NUMBER_OF_PADS];
322 for (
size_t i = 0; i != NUMBER_OF_PADS; ++i) {
323 Xs[i] = 0.003 * (Xmax[i] - Xmin[i]) * (0.5 * NUMBER_OF_PMTS);
326 TH2D H2[NUMBER_OF_PADS];
327 TGraph G2[NUMBER_OF_PADS];
329 for (
size_t i = 0; i != NUMBER_OF_PADS; ++i) {
331 H2[i] = TH2D(
MAKE_CSTRING(
"h" << i), NULL, graphics.nbinsX, Xmin[i] - Xs[i], Xmax[i] + Xs[i], graphics.nbinsY, ymin, ymax);
333 H2[i].GetXaxis()->SetTitle(Xlabel[i].c_str());
334 H2[i].GetYaxis()->SetTitle(
"#Deltat [ns]");
336 H2[i].GetXaxis()->CenterTitle(
true);
337 H2[i].GetYaxis()->CenterTitle(
true);
339 H2[i].SetStats(kFALSE);
343 G2[i].SetPoint(0, H2[i].GetXaxis()->GetXmin(), 0.0);
344 G2[i].SetPoint(1, H2[i].GetXaxis()->GetXmax(), 0.0);
355 cout <<
"event: " << setw(8) << inputFile.
getCounter() << endl;
357 const Evt* evt = inputFile.
next();
363 if (!event_selector(fit, *evt)) {
369 for (
const Hit& hit : evt->
hits) {
386 for (
const auto& trk : evt->
mc_trks) {
388 if (trk.E > muon.
E) {
400 bool monte_carlo =
false;
402 for (
bool next =
false; !next; ) {
404 for (
size_t i = 0; i != NUMBER_OF_PADS; ++i) {
429 for (JDataL0_t::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {
442 sort(data.begin(), data.end(), JHitW0::compare);
444 JDataW0_t::iterator __end = unique(data.begin(), data.end(), equal_to<JDAQPMTIdentifier>());
446 double E_GeV = parameters.
E_GeV;
457 for (JDataW0_t::iterator hit = data.begin(); hit != __end; ++hit) {
459 const double x = hit->getX() - tz.
getX();
460 const double y = hit->getY() - tz.
getY();
461 const double z = hit->getZ();
462 const double R = sqrt(x*x + y*y);
467 const double z0 = tz.
getZ();
481 marker[2].push_back(TMarker(z0 - tz.
getZ(), 0.0, kFullCircle));
484 static_cast<TAttMarker&
>(marker[2][0]) = TAttMarker(kRed, kFullCircle, 0.7);
485 static_cast<TAttMarker&
>(marker[2][1]) = TAttMarker(kRed, kFullCircle, 0.7);
496 for (JDataW0_t::const_iterator hit = data.begin(); hit != __end; ++hit) {
498 const double x = hit->getX() - tz.
getX();
499 const double y = hit->getY() - tz.
getY();
500 const double z = hit->getZ() - tz.
getZ();
501 const double R = sqrt(x*x + y*y);
505 JDirection3D dir(hit->getDX(), hit->getDY(), hit->getDZ());
509 const double theta = dir.
getTheta();
510 const double phi = fabs(dir.
getPhi());
513 const double E = E_GeV;
514 const double dt = T_ns.
constrain(hit->getT() - t1);
525 chi2 += H1.getChi2() - H0.getChi2();
528 << setw(8) << hit->getModuleID() <<
'.' <<
FILL(2,
'0') << (
int) hit->getPMTAddress() <<
FILL() <<
' '
530 <<
FIXED(7,2) << R <<
' '
531 <<
FIXED(7,4) << theta <<
' '
532 <<
FIXED(7,4) << phi <<
' '
533 <<
FIXED(7,3) << dt <<
' '
534 <<
FIXED(7,3) << H1.getChi2() <<
' '
535 <<
FIXED(7,3) << H0.getChi2() << endl);
537 const double derivative = H1.getDerivativeOfChi2() - H0.getDerivativeOfChi2();
539 double size = derivative * graphics.arrowScale;
541 if (fabs(size) < Amin) {
542 size = (size > 0.0 ? +Amin : -Amin);
543 }
else if (fabs(size) > Amax) {
544 size = (size > 0.0 ? +Amax : -Amax);
547 const double X[NUMBER_OF_PADS] = { R, atan2(y,x), z - R/
getTanThetaC() };
549 const double xs = (double) (NUMBER_OF_PMTS - 2 * hit->getPMTAddress()) / (double) NUMBER_OF_PMTS;
551 for (
size_t i = 0; i != NUMBER_OF_PADS; ++i) {
553 TArrow a1(X[i] + xs*Xs[i], dt + graphics.T_ns, X[i] + xs*Xs[i], dt + graphics.T_ns + size, graphics.arrowSize, graphics.arrowType.c_str());
555 a1.SetLineWidth(graphics.lineWidth);
556 a1.SetLineStyle(graphics.lineStyle);
558 arrow[i].push_back(a1);
560 H2[i].Fill(X[i], dt + graphics.T_ns);
565 os <<
" Q = " <<
FIXED(4,0) << trk.
lik
566 <<
'/' <<
FIXED(4,0) << -chi2;
567 os <<
" E = " <<
SCIENTIFIC(7,1) << trk.
E <<
" [GeV]";
568 os <<
" cos(#theta) = " <<
FIXED(6,3) << trk.
dir.
z;
571 os <<
" Monte Carlo";
578 TLatex title(0.05, 0.5, os.str().c_str());
580 title.SetTextAlign(12);
581 title.SetTextFont(42);
582 title.SetTextSize(0.6);
588 for (
int i = 0; i != NUMBER_OF_PADS; ++i) {
592 if (option == arrow_t) {
594 for (
auto& a1 : arrow[i]) {
598 for (
auto& m1 : marker[i]) {
603 if (option == histogram_t) {
621 static int count = 0;
624 cout << endl <<
"Type '?' for possible options." << endl;
627 for (
bool user =
true; user; ) {
629 cout <<
"\n> " << flush;
635 cout <<
"possible options: " << endl;
636 cout <<
'q' <<
" -> " <<
"exit application" << endl;
637 cout <<
'u' <<
" -> " <<
"update canvas" << endl;
638 cout <<
's' <<
" -> " <<
"save graphics to file" << endl;
639 cout <<
'M' <<
" -> " <<
"Monte Carlo true muon information" << endl;
640 cout <<
'F' <<
" -> " <<
"fit information" << endl;
642 cout <<
'L' <<
" -> " <<
"reload event selector" << endl;
644 cout <<
'r' <<
" -> " <<
"rewind input file" << endl;
645 cout <<
'R' <<
" -> " <<
"switch to ROOT mode (quit ROOT to continue)" << endl;
646 cout <<
'p' <<
" -> " <<
"print event information" << endl;
647 cout <<
' ' <<
" -> " <<
"next event (as well as any other key)" << endl;
666 ERROR(endl <<
"No Monte Carlo muon available." << endl);
677 execute(
MAKE_STRING(
"make -f " <<
getPath(argv[0]) <<
"/JMakeEventSelector libs"), 3);
690 for (
const auto& trk : evt->
mc_trks) {
691 cout <<
"MC "; trk.
print(cout); cout << endl;
693 for (
const auto& trk : evt->
trks) {
694 cout <<
"fit "; trk.
print(cout); cout << endl;
696 for (
const auto& hit : evt->
hits) {
697 cout <<
"hit "; hit.
print(cout); cout << endl;
int main(int argc, char **argv)
KM3NeT DAQ constants, bit handling, etc.
Basic data structure for L0 hit.
Keyboard settings for unbuffered input.
General purpose messaging.
#define DEBUG(A)
Message macros.
Auxiliary data structure for muon PDF.
Utility class to parse command line options.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
I/O formatting auxiliaries.
#define MAKE_CSTRING(A)
Make C-string.
#define MAKE_STRING(A)
Make string.
Utility class to parse parameter values.
#define gmake_property(A)
macros to convert (template) parameter to JPropertiesElement object
Scanning of objects from a single file according a format that follows from the extension of each fil...
ROOT TTree parameter settings of various packages.
Utility class to parse parameter values.
Data structure for fit of straight line paralel to z-axis.
double getT(const JVector3D &pos) const
Get arrival time of Cherenkov light at given position.
double getZ(const JPosition3D &pos) const
Get point of emission of Cherenkov light along muon path.
void setZ(const double z, const double velocity)
Set z-position of vertex.
JAxis3D & rotate(const JRotation3D &R)
Rotate axis.
Data structure for direction in three dimensions.
JDirection3D & rotate(const JRotation3D &R)
Rotate.
double getY() const
Get y position.
double getX() const
Get x position.
double getTheta() const
Get theta angle.
double getPhi() const
Get phi angle.
Utility class to parse command line options.
Auxiliary class for a hit with background rate value.
Data structure for size of TCanvas.
int y
number of pixels in Y
int x
number of pixels in X
Wrapper class around ROOT TStyle.
Object reading from a list of files.
virtual void rewind() override
Rewind.
virtual bool hasNext() override
Check availability of next element.
counter_type getCounter() const
Get counter.
virtual const pointer_type & next() override
Get next element.
File router for fast addressing of summary data.
void update(const JDAQHeader &header)
Update router.
double getRate(const JDAQPMTIdentifier &id, const double rate_Hz) const
Get rate.
Enable unbuffered terminal input.
Streaming of input and output from Linux command.
Data structure for L0 hit.
static void setSlewing(const bool slewing)
Set slewing option.
Data structure for UTC time.
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 JMUONSTART
static const int JSTART_LENGTH_METRES
distance between projected positions on the track of optical modules for which the response does not ...
Extensions to Evt data format.
JDirection3D getDirection(const Vec &dir)
Get direction.
bool hasW(const Trk &trk, const int i)
Check availability of value.
bool has_muon(const Evt &evt)
Test whether given event has a muon.
void setW(Trk &trk, const int i, const double value)
Set associated value.
double getW(const Trk &track, const size_t index, const double value)
Get track information.
JPosition3D getPosition(const Vec &pos)
Get position.
bool is_muon(const Trk &track)
Test whether given track is a (anti-)muon.
std::string getPath(const std::string &file_name)
Get path, i.e. part before last JEEP::PATHNAME_SEPARATOR if any.
double getAngle(const JQuaternion3D &first, const JQuaternion3D &second)
Get space angle between quanternions.
std::istream & getline(std::istream &in, JString &object)
Read string from input stream until end of line.
std::string replace(const std::string &input, const std::string &target, const std::string &replacement)
Replace tokens in string.
static const double PI
Mathematical constants.
const double getInverseSpeedOfLight()
Get inverse speed of light.
double getTanThetaC()
Get average tangent of Cherenkov angle of water corresponding to group velocity.
const double getSpeedOfLight()
Get speed of light.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
const JFit & get_best_reconstructed_track(const JEvt &evt, JTrackSelector_t selector, JQualitySorter_t comparator)
Get best reconstructed track.
bool has_reconstructed_track(const JEvt &evt, JTrackSelector_t selector)
Test whether given event has a track according selection.
KM3NeT DAQ data structures and auxiliaries.
static const char WILDCARD
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
int frame_index
from the raw data
int run_id
DAQ run identifier.
std::vector< Hit > hits
list of hits
std::vector< Trk > mc_trks
MC: list of MC truth tracks.
int det_id
detector identifier from DAQ
ULong64_t trigger_counter
trigger counter
void print(std::ostream &out=std::cout) const
Print event.
std::vector< Trk > trks
list of reconstructed tracks (can be several because of prefits,showers, etc).
TTimeStamp t
UTC time of the timeslice, or the event_time for MC. (default: 01 Jan 1970 00:00:00)
Auxiliary data structure for sequence of same character.
Auxiliary data structure for floating point format specification.
int dom_id
module identifier from the data (unique in the detector).
void print(std::ostream &out=std::cout) const
Print hit.
Vec dir
hit direction; i.e. direction of the PMT
unsigned int channel_id
PMT channel id {0,1, .., 30} local to moduke.
unsigned int tot
tot value as stored in raw data (int for pyroot)
double t
hit time (from tdc+calibration or MC truth)
JEventSelector()
Default constructor.
static bool select(const Trk &trk, const Evt &evt)
Default event selection.
Model for fit to acoustics data.
void(*)(Args...) function
bool is_valid() const
Check validity of function.
void reload()
Reload function from shared library.
Auxiliary data structure for muon PDF.
JFunction1D_t::result_type result_type
result_type calculate(const double E, const double R, const double theta, const double phi, const double t1) const
Get PDF.
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Data structure for fit parameters.
double TTS_ns
transition-time spread [ns]
double TMin_ns
minimal time w.r.t. Cherenkov hypothesis [ns]
double roadWidth_m
road width [m]
double TMax_ns
maximal time w.r.t. Cherenkov hypothesis [ns]
double VMax_npe
maximum number of of photo-electrons
double R_Hz
default rate [Hz]
size_t numberOfPrefits
number of prefits
Auxiliary class for defining the range of iterations of objects.
static counter_type max()
Get maximum counter value.
Auxiliary data structure for floating point format specification.
The Trk class represents a Monte Carlo (MC) particle as well as a reconstructed track/shower.
void print(std::ostream &out=std::cout) const
Print track.
std::vector< double > fitinf
place to store additional fit info, see km3net-dataformat/definitions/fitparameters....
double E
Energy [GeV] (either MC truth or reconstructed)
double t
track time [ns] (when the particle is at pos )
double len
length, if applicable [m]
double lik
likelihood or lambda value (for aafit, lambda)
Range of reconstruction stages.
Auxiliary include file for time conversion between DAQ/trigger hit and Monte Carlo hit.