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(), parameters.
R_Hz));
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 <<
'/' <<
FIXED(4,0) << -chi2;
567 os <<
" E = " <<
SCIENTIFIC(7,1) << trk.
E <<
" [GeV]";
568 os <<
" cos(#theta) = " <<
FIXED(6,3) << trk.
dir.
z;
569
570 if (monte_carlo)
571 os << " Monte Carlo";
574
575
576
577
578 TLatex title(0.05, 0.5, os.str().c_str());
579
580 title.SetTextAlign(12);
581 title.SetTextFont(42);
582 title.SetTextSize(0.6);
583
584 p2->cd();
585
586 title.Draw();
587
588 for (int i = 0; i != NUMBER_OF_PADS; ++i) {
589
591
592 if (option == arrow_t) {
593
594 for (auto& a1 : arrow[i]) {
595 a1.Draw();
596 }
597
598 for (auto& m1 : marker[i]) {
599 m1.Draw();
600 }
601 }
602
603 if (option == histogram_t) {
604 H2[i].Draw("SAME");
605 }
606 }
607
608 cv->Update();
609
610
611
612
613 if (batch) {
614
616
617 next = true;
618
619 } else {
620
621 static int count = 0;
622
623 if (count++ == 0) {
624 cout << endl << "Type '?' for possible options." << endl;
625 }
626
627 for (bool user = true; user; ) {
628
629 cout << "\n> " << flush;
630
632
633 case '?':
634 cout << endl;
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;
643 }
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;
648 break;
649
650 case 'q':
651 cout << endl;
652 return 0;
653
654 case 'u':
655 cv->Update();
656 break;
657
658 case 's':
660 break;
661
662 case 'M':
664 monte_carlo = true;
665 else
666 ERROR(endl <<
"No Monte Carlo muon available." << endl);
667 user = false;
668 break;
669
670 case 'F':
671 monte_carlo = false;
672 user = false;
673 break;
674
675 case 'L':
677 execute(
MAKE_STRING(
"make -f " <<
getPath(argv[0]) <<
"/JMakeEventSelector libs"), 3);
679 }
680 break;
681
682 case 'R':
683 tp->Run(kTRUE);
684 break;
685
686 case 'p':
687 cout << endl;
690 for (
const auto& trk : evt->
mc_trks) {
691 cout <<
"MC "; trk.
print(cout); cout << endl;
692 }
693 for (
const auto& trk : evt->
trks) {
694 cout <<
"fit "; trk.
print(cout); cout << endl;
695 }
696 for (
const auto& hit : evt->
hits) {
697 cout <<
"hit "; hit.
print(cout); cout << endl;
698 }
699 }
700 break;
701
702 case 'r':
704
705 default:
706 next = true;
707 user = false;
708 break;
709 }
710 }
711 }
712 }
713 }
714 }
715 cout << endl;
716}
#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 projected positions on the track of optical modules for which the response does not ...
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.