195{
199
201 JLimit_t& numberOfEvents = inputFile.getLimit();
202 string pdfFile;
205 int application;
208 bool batch;
210 double arrowSize = 0.003;
211 string arrowType = "|->";
212 double arrowScale = 250.0;
213 Width_t lineWidth = 2;
214 Style_t lineStyle = 1;
215 int nbinsX = 50;
216 int nbinsY = 250;
217 double T_ns = 0.0;
218 } graphics;
219 string option;
221
222
223 try {
224
226
233
235
236 JParser<> zap(
"Program to display hit probabilities.");
237
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");
250
251 zap(argc, argv);
252 }
253 catch(const exception& error) {
254 FATAL(error.what() << endl);
255 }
256
258 FATAL(
"Missing output file name " <<
outputFile <<
" in batch mode." << endl);
259 }
260
263 }
264
267 }
268
270
272
274
276
279
280
281
282
283 gROOT->SetBatch(batch);
284
285 TApplication* tp = new TApplication("user", NULL, NULL);
286 TCanvas* cv =
new TCanvas(
"display",
"", canvas.
x, canvas.
y);
287
288 unique_ptr<TStyle> gStyle(
new JStyle(
"gplot", cv->GetWw(), cv->GetWh(), graphics));
289
290 gROOT->SetStyle("gplot");
291 gROOT->ForceStyle();
292
293 const size_t NUMBER_OF_PADS = 3;
294
295 cv->SetFillStyle(4000);
296 cv->SetFillColor(kWhite);
297
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);
300
301 p1->Divide(NUMBER_OF_PADS, 1);
302
304 p2->Draw();
305
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);
315
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 };
319
320 double Xs[NUMBER_OF_PADS];
321
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);
324 }
325
326 TH2D H2[NUMBER_OF_PADS];
327 TGraph G2[NUMBER_OF_PADS];
328
329 for (size_t i = 0; i != NUMBER_OF_PADS; ++i) {
330
331 H2[i] = TH2D(
MAKE_CSTRING(
"h" << i), NULL, graphics.nbinsX, Xmin[i] - Xs[i], Xmax[i] + Xs[i], graphics.nbinsY, ymin, ymax);
332
333 H2[i].GetXaxis()->SetTitle(Xlabel[i].c_str());
334 H2[i].GetYaxis()->SetTitle("#Deltat [ns]");
335
336 H2[i].GetXaxis()->CenterTitle(true);
337 H2[i].GetYaxis()->CenterTitle(true);
338
339 H2[i].SetStats(kFALSE);
340
341 G2[i].Set(2);
342
343 G2[i].SetPoint(0, H2[i].GetXaxis()->GetXmin(), 0.0);
344 G2[i].SetPoint(1, H2[i].GetXaxis()->GetXmax(), 0.0);
345
347
348 H2[i].Draw("AXIS");
349 G2[i].Draw("SAME");
350 }
351
352
354
355 cout <<
"event: " << setw(8) << inputFile.
getCounter() << endl;
356
357 const Evt* evt = inputFile.
next();
358
360
362
363 if (!event_selector(fit, *evt)) {
364 continue;
365 }
366
367 JDataL0_t dataL0;
368
369 for (
const Hit& hit : evt->
hits) {
373 }
374
379
381
383
385
386 for (
const auto& trk : evt->
mc_trks) {
388 if (trk.E > muon.
E) {
389
390 muon = trk;
392
394 }
395 }
396 }
397 }
398
399
400 bool monte_carlo = false;
401
402 for (bool next = false; !next; ) {
403
404 for (size_t i = 0; i != NUMBER_OF_PADS; ++i) {
405 H2[i].Reset();
406 }
407
409
410 if (!monte_carlo)
411 trk = fit;
412 else
413 trk = muon;
414
418
419
420
421
422
424
425
426
428
429 for (JDataL0_t::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {
430
431 JHitW0 hit(*i, summary.getRate(i->getPMTIdentifier()));
432
433 hit.rotate(R);
434
435 if (match(hit)) {
437 }
438 }
439
440
441
442 sort(
data.begin(),
data.end(), JHitW0::compare);
443
444 JDataW0_t::iterator __end = unique(
data.begin(),
data.end(), equal_to<JDAQPMTIdentifier>());
445
446 double E_GeV = parameters.
E_GeV;
447
448
449
450
451
452
453
454
456
457 for (JDataW0_t::iterator hit =
data.begin(); hit != __end; ++hit) {
458
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);
463
465 }
466
467 const double z0 = tz.getZ();
469
471
472
473
474
475 ostringstream os;
478
480
481 marker[2].push_back(TMarker(z0 - tz.getZ(), 0.0, kFullCircle));
483
484 static_cast<TAttMarker&>(marker[2][0]) = TAttMarker(kRed, kFullCircle, 0.7);
485 static_cast<TAttMarker&>(marker[2][1]) = TAttMarker(kRed, kFullCircle, 0.7);
486 }
487
489 <<
FIXED(7,2) << tz.getX() <<
' '
490 <<
FIXED(7,2) << tz.getY() <<
' '
491 <<
FIXED(7,2) << tz.getZ() <<
' '
492 <<
FIXED(12,2) << tz.getT() << endl);
493
494 double chi2 = 0;
495
496 for (JDataW0_t::const_iterator hit =
data.begin(); hit != __end; ++hit) {
497
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);
502
504
505 JDirection3D dir(hit->getDX(), hit->getDY(), hit->getDZ());
506
508
509 const double theta = dir.getTheta();
510 const double phi = fabs(dir.getPhi());
511
512
513 const double E = E_GeV;
514 const double dt = T_ns.constrain(hit->getT() - t1);
515
518
521 }
522
523 H1 += H0;
524
525 chi2 += H1.getChi2() - H0.getChi2();
526
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);
536
537 const double derivative = H1.getDerivativeOfChi2() - H0.getDerivativeOfChi2();
538
539 double size = derivative * graphics.arrowScale;
540
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);
545 }
546
547 const double X[NUMBER_OF_PADS] = { R, atan2(y,x), z - R/
getTanThetaC() };
548
549 const double xs = (double) (NUMBER_OF_PMTS - 2 * hit->getPMTAddress()) / (double) NUMBER_OF_PMTS;
550
551 for (size_t i = 0; i != NUMBER_OF_PADS; ++i) {
552
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());
554
555 a1.SetLineWidth(graphics.lineWidth);
556 a1.SetLineStyle(graphics.lineStyle);
557
558 arrow[i].push_back(a1);
559
560 H2[i].Fill(X[i], dt + graphics.T_ns);
561 }
562 }
563
565 os <<
" Q = " <<
FIXED(4,0) << trk.
lik;
566 os <<
" E = " <<
SCIENTIFIC(7,1) << trk.
E <<
" [GeV]";
567 os <<
" cos(#theta) = " <<
FIXED(6,3) << trk.
dir.
z;
568
569 if (monte_carlo)
570 os << " Monte Carlo";
573
574
575
576
577 TLatex title(0.05, 0.5, os.str().c_str());
578
579 title.SetTextAlign(12);
580 title.SetTextFont(42);
581 title.SetTextSize(0.6);
582
583 p2->cd();
584
585 title.Draw();
586
587 for (int i = 0; i != NUMBER_OF_PADS; ++i) {
588
590
591 if (option == arrow_t) {
592
593 for (auto& a1 : arrow[i]) {
594 a1.Draw();
595 }
596
597 for (auto& m1 : marker[i]) {
598 m1.Draw();
599 }
600 }
601
602 if (option == histogram_t) {
603 H2[i].Draw("SAME");
604 }
605 }
606
607 cv->Update();
608
609
610
611
612 if (batch) {
613
615
616 next = true;
617
618 } else {
619
620 static int count = 0;
621
622 if (count++ == 0) {
623 cout << endl << "Type '?' for possible options." << endl;
624 }
625
626 for (bool user = true; user; ) {
627
628 cout << "\n> " << flush;
629
631
632 case '?':
633 cout << endl;
634 cout << "possible options: " << endl;
635 cout << 'q' << " -> " << "exit application" << endl;
636 cout << 'u' << " -> " << "update canvas" << endl;
637 cout << 's' << " -> " << "save graphics to file" << endl;
638 cout << 'M' << " -> " << "Monte Carlo true muon information" << endl;
639 cout << 'F' << " -> " << "fit information" << endl;
641 cout << 'L' << " -> " << "reload event selector" << endl;
642 }
643 cout << 'r' << " -> " << "rewind input file" << endl;
644 cout << 'R' << " -> " << "switch to ROOT mode (quit ROOT to continue)" << endl;
645 cout << 'p' << " -> " << "print event information" << endl;
646 cout << ' ' << " -> " << "next event (as well as any other key)" << endl;
647 break;
648
649 case 'q':
650 cout << endl;
651 return 0;
652
653 case 'u':
654 cv->Update();
655 break;
656
657 case 's':
659 break;
660
661 case 'M':
663 monte_carlo = true;
664 else
665 ERROR(endl <<
"No Monte Carlo muon available." << endl);
666 user = false;
667 break;
668
669 case 'F':
670 monte_carlo = false;
671 user = false;
672 break;
673
674 case 'L':
676 execute(
MAKE_STRING(
"make -f " <<
getPath(argv[0]) <<
"/JMakeEventSelector libs"), 3);
678 }
679 break;
680
681 case 'R':
682 tp->Run(kTRUE);
683 break;
684
685 case 'p':
686 cout << endl;
689 for (
const auto& trk : evt->
mc_trks) {
690 cout <<
"MC "; trk.
print(cout); cout << endl;
691 }
692 for (
const auto& trk : evt->
trks) {
693 cout <<
"fit "; trk.
print(cout); cout << endl;
694 }
695 for (
const auto& hit : evt->
hits) {
696 cout <<
"hit "; hit.
print(cout); cout << endl;
697 }
698 }
699 break;
700
701 case 'r':
703
704 default:
705 next = true;
706 user = false;
707 break;
708 }
709 }
710 }
711 }
712 }
713 }
714 cout << endl;
715}
#define DEBUG(A)
Message macros.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
#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
Utility class to parse parameter values.
Data structure for fit of straight line paralel to z-axis.
Data structure for direction in three dimensions.
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.
Enable unbuffered terminal input.
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 first and last hits in metres from JStart.cc
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::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)
Model for fit to acoustics data.
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
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.
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.