113 typedef JParallelFileScanner_t::multi_pointer_type multi_pointer_type;
118 JParallelFileScanner_t inputFile;
121 JCalibration_t calibrationFile;
131 double arrowSize = 0.003;
132 string arrowType =
"|->";
133 double arrowScale = 250.0;
134 Width_t lineWidth = 2;
135 Style_t lineStyle = 1;
157 JParser<> zap(
"Program to display hit probabilities.");
159 zap[
'w'] =
make_field(canvas,
"size of canvas <nx>x<ny> [pixels]") =
JCanvas(1200, 600);
160 zap[
'f'] =
make_field(inputFile,
"input file (output of JXXXMuonReconstruction.sh)");
164 zap[
'n'] =
make_field(numberOfEvents) = JLimit::max();
171 zap[
'O'] =
make_field(option,
"draw option") = arrow_t, histogram_t;
172 zap[
'B'] =
make_field(batch,
"batch processing");
177 catch(
const exception& error) {
178 FATAL(error.what() << endl);
182 FATAL(
"Missing output file name " <<
outputFile <<
" in batch mode." << endl);
202 unique_ptr<JDynamics> dynamics;
208 dynamics->load(calibrationFile);
210 catch(
const exception& error) {
211 if (!calibrationFile.empty()) {
230 Vec offset(0.0, 0.0, 0.0);
234 }
catch(
const exception& error) {}
239 gROOT->SetBatch(batch);
241 TApplication* tp =
new TApplication(
"user", NULL, NULL);
242 TCanvas* cv =
new TCanvas(
"display",
"", canvas.
x, canvas.
y);
244 unique_ptr<TStyle> gStyle(
new JStyle(
"gplot", cv->GetWw(), cv->GetWh(), graphics));
246 gROOT->SetStyle(
"gplot");
249 const size_t NUMBER_OF_PADS = 3;
251 cv->SetFillStyle(4000);
252 cv->SetFillColor(kWhite);
254 TPad*
p1 =
new TPad(
"p1", NULL, 0.0, 0.00, 1.0, 0.95);
255 TPad* p2 =
new TPad(
"p2", NULL, 0.0, 0.95, 1.0, 1.00);
257 p1->Divide(NUMBER_OF_PADS, 1);
263 const double Rmin = 0.0;
264 const double Rmax = min(parameters.
roadWidth_m, 0.4 * Dmax);
265 const double Tmin = min(parameters.
TMin_ns, -10.0);
266 const double Tmax = max(parameters.
TMax_ns, +100.0);
267 const double Amin = 0.002 * (Tmax - Tmin);
268 const double Amax = 0.8 * (Tmax - Tmin);
269 const double ymin = Tmin - (option == arrow_t ? 0.2 * Amax : 0.0);
270 const double ymax = Tmax + (option == arrow_t ? 0.5 * Amax : 0.0);
272 const string Xlabel[NUMBER_OF_PADS] = {
"R [m]",
"#phi [rad]",
"z [m]" };
273 const double Xmin [NUMBER_OF_PADS] = { Rmin, -
PI, -0.3 * Dmax };
274 const double Xmax [NUMBER_OF_PADS] = { Rmax, +
PI, +0.3 * Dmax };
276 double Xs[NUMBER_OF_PADS];
278 for (
size_t i = 0; i != NUMBER_OF_PADS; ++i) {
282 TH2D H2[NUMBER_OF_PADS];
283 TGraph G2[NUMBER_OF_PADS];
285 for (
size_t i = 0; i != NUMBER_OF_PADS; ++i) {
287 H2[i] = TH2D(
MAKE_CSTRING(
"h" << i), NULL, graphics.nbinsX, Xmin[i] - Xs[i], Xmax[i] + Xs[i], graphics.nbinsY, ymin, ymax);
289 H2[i].GetXaxis()->SetTitle(Xlabel[i].c_str());
290 H2[i].GetYaxis()->SetTitle(
"#Deltat [ns]");
292 H2[i].GetXaxis()->CenterTitle(
true);
293 H2[i].GetYaxis()->CenterTitle(
true);
295 H2[i].SetStats(kFALSE);
299 G2[i].SetPoint(0, H2[i].GetXaxis()->GetXmin(), 0.0);
300 G2[i].SetPoint(1, H2[i].GetXaxis()->GetXmax(), 0.0);
311 cout <<
"event: " << setw(8) << inputFile.getCounter() << endl;
313 multi_pointer_type ps = inputFile.next();
320 dynamics->update(*tev);
323 if (mc.getEntries() != 0) {
333 if (!event_selector(*tev, *in, event)) {
340 buildL0(*tev, router,
true, back_inserter(dataL0));
342 summary.update(*tev);
350 for (
const auto& t1 : event->mc_trks) {
352 if (t1.E > muon.
getE()) {
359 muon =
getFit(0, ta, 0.0, 0, t1.E, 1);
367 bool monte_carlo =
false;
370 for (
bool next =
false; !next; ) {
372 for (
size_t i = 0; i != NUMBER_OF_PADS; ++i) {
397 for (JDataL0_t::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {
399 JHitW0 hit(*i, summary.getRate(i->getPMTIdentifier()));
410 sort(
data.begin(),
data.end(), JHitW0::compare);
412 JDataW0_t::iterator __end = unique(
data.begin(),
data.end(), equal_to<JDAQPMTIdentifier>());
414 double E_GeV = parameters.
E_GeV;
425 for (JDataW0_t::iterator hit =
data.begin(); hit != __end; ++hit) {
427 const double x = hit->getX() - tz.getX();
428 const double y = hit->getY() - tz.getY();
429 const double z = hit->getZ();
430 const double R = sqrt(
x*
x +
y*
y);
435 const double z0 = tz.getZ();
448 marker[2].push_back(TMarker(z0 - tz.getZ(), 0.0, kFullCircle));
451 static_cast<TAttMarker&
>(marker[2][0]) = TAttMarker(kRed, kFullCircle, 0.7);
452 static_cast<TAttMarker&
>(marker[2][1]) = TAttMarker(kRed, kFullCircle, 0.7);
456 <<
FIXED(7,2) << tz.getX() <<
' '
457 <<
FIXED(7,2) << tz.getY() <<
' '
458 <<
FIXED(7,2) << tz.getZ() <<
' '
459 <<
FIXED(12,2) << tz.getT() << endl);
463 for (JDataW0_t::const_iterator hit =
data.begin(); hit != __end; ++hit) {
465 const double x = hit->getX() - tz.getX();
466 const double y = hit->getY() - tz.getY();
467 const double z = hit->getZ() - tz.getZ();
468 const double R = sqrt(
x*
x +
y*
y);
472 JDirection3D dir(hit->getDX(), hit->getDY(), hit->getDZ());
476 const double theta = dir.getTheta();
477 const double phi = fabs(dir.getPhi());
480 const double E = E_GeV;
481 const double dt = T_ns.constrain(hit->getT() - t1);
492 chi2 += H1.getChi2() - H0.getChi2();
495 << setw(8) << hit->getModuleID() <<
'.' <<
FILL(2,
'0') << (
int) hit->getPMTAddress() <<
FILL() <<
' '
497 <<
FIXED(7,2) << R <<
' '
498 <<
FIXED(7,4) << theta <<
' '
499 <<
FIXED(7,4) << phi <<
' '
500 <<
FIXED(7,3) << dt <<
' '
501 <<
FIXED(7,3) << H1.getChi2() <<
' '
502 <<
FIXED(7,3) << H0.getChi2() << endl);
504 const double derivative = H1.getDerivativeOfChi2() - H0.getDerivativeOfChi2();
506 double size = derivative * graphics.arrowScale;
508 if (fabs(size) < Amin) {
509 size = (size > 0.0 ? +Amin : -Amin);
510 }
else if (fabs(size) > Amax) {
511 size = (size > 0.0 ? +Amax : -Amax);
514 const double X[NUMBER_OF_PADS] = { R, atan2(
y,
x), z - R/
getTanThetaC() };
518 for (
size_t i = 0; i != NUMBER_OF_PADS; ++i) {
520 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());
522 a1.SetLineWidth(graphics.lineWidth);
523 a1.SetLineStyle(graphics.lineStyle);
525 arrow[i].push_back(a1);
527 H2[i].Fill(X[i], dt + graphics.T_ns);
532 os <<
" Q = " <<
FIXED(4,0) << fit.
getQ();
534 os <<
" cos(#theta) = " <<
FIXED(6,3) << fit.
getDZ();
537 os <<
" Monte Carlo";
544 TLatex title(0.05, 0.5, os.str().c_str());
546 title.SetTextAlign(12);
547 title.SetTextFont(42);
548 title.SetTextSize(0.6);
554 for (
int i = 0; i != NUMBER_OF_PADS; ++i) {
558 if (option == arrow_t) {
560 for (
auto& a1 : arrow[i]) {
564 for (
auto& m1 : marker[i]) {
569 if (option == histogram_t) {
587 static int count = 0;
590 cout << endl <<
"Type '?' for possible options." << endl;
593 for (
bool user =
true; user; ) {
595 cout <<
"\n> " << flush;
601 cout <<
"possible options: " << endl;
602 cout <<
'q' <<
" -> " <<
"exit application" << endl;
603 cout <<
'u' <<
" -> " <<
"update canvas" << endl;
604 cout <<
's' <<
" -> " <<
"save graphics to file" << endl;
605 cout <<
'+' <<
" -> " <<
"next fit" << endl;
606 cout <<
'-' <<
" -> " <<
"previous fit" << endl;
607 cout <<
'M' <<
" -> " <<
"Monte Carlo true muon information" << endl;
608 cout <<
'F' <<
" -> " <<
"fit information" << endl;
610 cout <<
'L' <<
" -> " <<
"reload event selector" << endl;
612 cout <<
'r' <<
" -> " <<
"rewind input file" << endl;
613 cout <<
'R' <<
" -> " <<
"switch to ROOT mode (quit ROOT to continue)" << endl;
614 cout <<
' ' <<
" -> " <<
"next event (as well as any other key)" << endl;
631 index = (index != in->size() - 1 ? index + 1 : 0);
637 index = (index != 0 ? index - 1 : in->size() - 1);
645 ERROR(endl <<
"No Monte Carlo muon available." << endl);
656 execute(
MAKE_STRING(
"make -f " <<
getPath(argv[0]) <<
"/JMakeEventSelector libs"), 3);
#define DEBUG(A)
Message macros.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
TString replace(const TString &target, const TRegexp ®exp, const T &replacement)
Replace regular expression in input by given replacement.
#define MAKE_CSTRING(A)
Make C-string.
#define MAKE_STRING(A)
Make string.
#define gmake_property(A)
macros to convert (template) parameter to JPropertiesElement object
Router for direct addressing of module data in detector data structure.
Utility class to parse parameter values.
Data structure for set of track fit results.
void select(const JSelector_t &selector)
Select fits.
Data structure for track fit results with history and optional associated values.
void setW(const std::vector< double > &W)
Set associated values.
const std::vector< double > & getW() const
Get associated values.
double getDZ() const
Get Z-slope.
double getE() const
Get energy.
int getStatus() const
Get status of the fit; negative values should refer to a bad fit.
double getQ() const
Get quality.
double getT() const
Get time.
bool hasW(const int i) const
Check availability of value.
Data structure for fit of straight line paralel to z-axis.
Data structure for direction in three dimensions.
JTime & add(const JTime &value)
Addition operator.
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.
General purpose class for object reading from a list of file names.
General purpose class for parallel reading of objects from a single file or multiple files.
File router for fast addressing of summary data.
Enable unbuffered terminal input.
Template specialisation of L0 builder for JHitL0 data type.
int getRunNumber() const
Get run number.
int getFrameIndex() const
Get frame index.
JTriggerCounter_t getCounter() const
Get trigger counter.
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 JMUONENERGY
static const int JMUONSTART
static const int JSTART_LENGTH_METRES
distance between first and last hits in metres from JStart.cc
JDirection3D getDirection(const Vec &dir)
Get direction.
JTrack3E getTrack(const Trk &track)
Get track.
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.
JFit getFit(const int id, const JMODEL::JString &string)
Get fit parameters of string.
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
double getMaximalDistance(const JDetector &detector, const bool option=false)
Get maximal distance between modules in detector.
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.
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).
bool qualitySorter(const JFit &first, const JFit &second)
Comparison of fit results.
Head getHeader(const JMultipleFileScanner_t &file_list)
Get Monte Carlo header.
KM3NeT DAQ data structures and auxiliaries.
static const char WILDCARD
static const int NUMBER_OF_PMTS
Total number of PMTs in module.
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
Auxiliary data structure for sequence of same character.
Auxiliary data structure for floating point format specification.
Dynamic detector calibration.
bool is_valid() const
Check validity of function.
void reload()
Reload function from shared library.
Auxiliary class to test history.
Template specialisation of class JModel to match hit with muon trajectory along z-axis.
Auxiliary class for recursive type list generation.
Auxiliary data structure for muon PDF.
JFunction1D_t::result_type result_type
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.
const JLimit & getLimit() const
Get limit.
Auxiliary data structure for floating point format specification.
The Vec class is a straightforward 3-d vector, which also works in pyroot.