109{
113
115 typedef JParallelFileScanner_t::multi_pointer_type multi_pointer_type;
119
120 JParallelFileScanner_t inputFile;
122 string detectorFile;
123 JCalibration_t calibrationFile;
124 double Tmax_s;
125 string pdfFile;
128 int application;
131 bool batch;
133 double arrowSize = 0.003;
134 string arrowType = "|->";
135 double arrowScale = 250.0;
136 Width_t lineWidth = 2;
137 Style_t lineStyle = 1;
138 int nbinsX = 50;
139 int nbinsY = 250;
140 double T_ns = 0.0;
141 } graphics;
142 string option;
144
145
146 try {
147
149
156
158
159 JParser<> zap(
"Program to display hit probabilities.");
160
161 zap[
'w'] =
make_field(canvas,
"size of canvas <nx>x<ny> [pixels]") =
JCanvas(1200, 600);
162 zap[
'f'] =
make_field(inputFile,
"input file (output of JXXXMuonReconstruction.sh)");
173 zap[
'O'] =
make_field(option,
"draw option") = arrow_t, histogram_t;
174 zap[
'B'] =
make_field(batch,
"batch processing");
176
177 zap(argc, argv);
178 }
179 catch(const exception& error) {
180 FATAL(error.what() << endl);
181 }
182
184 FATAL(
"Missing output file name " <<
outputFile <<
" in batch mode." << endl);
185 }
186
189 }
190
193 }
194
195
197
198 try {
200 }
203 }
204 unique_ptr<JDynamics> dynamics;
205
206 try {
207
209
210 dynamics->load(calibrationFile);
211 }
212 catch(const exception& error) {
213 if (!calibrationFile.empty()) {
215 }
216 }
217
218 const double Zbed = 0.0;
219
221
223
224 if (cylinder.getZmin() < Zbed) {
225 cylinder.setZmin(Zbed);
226 }
227
229
231
233
235
238
240
241
242 Vec offset(0.0, 0.0, 0.0);
243
244 try {
246 } catch(const exception& error) {}
247
248
249
250
251 gROOT->SetBatch(batch);
252
253 TApplication* tp = new TApplication("user", NULL, NULL);
254 TCanvas* cv =
new TCanvas(
"display",
"", canvas.
x, canvas.
y);
255
256 if (!batch) {
257 ((TRootCanvas *) cv->GetCanvasImp())->Connect("CloseWindow()", "TApplication", tp, "Terminate()");
258 }
259
260 unique_ptr<TStyle> gStyle(
new JStyle(
"gplot", cv->GetWw(), cv->GetWh(), graphics));
261
262 gROOT->SetStyle("gplot");
263 gROOT->ForceStyle();
264
265 const size_t NUMBER_OF_PADS = 3;
266
267 cv->SetFillStyle(4000);
268 cv->SetFillColor(kWhite);
269
270 TPad*
p1 =
new TPad(
"p1", NULL, 0.0, 0.00, 1.0, 0.95);
271 TPad* p2 = new TPad("p2", NULL, 0.0, 0.95, 1.0, 1.00);
272
273 p1->Divide(NUMBER_OF_PADS, 1);
274
276 p2->Draw();
277
279 const double Rmin = 0.0;
280 const double Rmax = min(parameters.
roadWidth_m, 0.4 * Dmax);
281 const double Tmin = min(parameters.
TMin_ns, -10.0);
282 const double Tmax = max(parameters.
TMax_ns, +100.0);
283 const double Amin = 0.002 * (Tmax - Tmin);
284 const double Amax = 0.8 * (Tmax - Tmin);
285 const double ymin = Tmin - (option == arrow_t ? 0.2 * Amax : 0.0);
286 const double ymax = Tmax + (option == arrow_t ? 0.5 * Amax : 0.0);
287
288 const string Xlabel[NUMBER_OF_PADS] = { "R [m]", "#phi [rad]", "z [m]" };
289 const double Xmin [NUMBER_OF_PADS] = { Rmin, -
PI, -0.4 * Dmax };
290 const double Xmax [NUMBER_OF_PADS] = { Rmax, +
PI, +0.4 * Dmax };
291
292 double Xs[NUMBER_OF_PADS];
293
294 for (size_t i = 0; i != NUMBER_OF_PADS; ++i) {
295 Xs[i] = 0.003 * (Xmax[i] - Xmin[i]) * (0.5 * NUMBER_OF_PMTS);
296 }
297
298 TH2D H2[NUMBER_OF_PADS];
299 TGraph G2[NUMBER_OF_PADS];
300
301 for (size_t i = 0; i != NUMBER_OF_PADS; ++i) {
302
303 H2[i] = TH2D(
MAKE_CSTRING(
"h" << i), NULL, graphics.nbinsX, Xmin[i] - Xs[i], Xmax[i] + Xs[i], graphics.nbinsY, ymin, ymax);
304
305 H2[i].GetXaxis()->SetTitle(Xlabel[i].c_str());
306 H2[i].GetYaxis()->SetTitle("#Deltat [ns]");
307
308 H2[i].GetXaxis()->CenterTitle(true);
309 H2[i].GetYaxis()->CenterTitle(true);
310
311 H2[i].SetStats(kFALSE);
312
313 G2[i].Set(2);
314
315 G2[i].SetPoint(0, H2[i].GetXaxis()->GetXmin(), 0.0);
316 G2[i].SetPoint(1, H2[i].GetXaxis()->GetXmax(), 0.0);
317
319
320 H2[i].Draw("AXIS");
321 G2[i].Draw("SAME");
322 }
323
324
326
327 cout << "event: " << setw(8) << inputFile.getCounter() << endl;
328
329 multi_pointer_type ps = inputFile.next();
330
334
335 if (dynamics) {
336 dynamics->update(*tev);
337 }
338
339 if (mc.getEntries() != 0) {
341 }
342
344
345 if (!in->empty()) {
346
348
349 if (!event_selector(*tev, *in, event)) {
350 continue;
351 }
352
353
354 JDataL0_t dataL0;
355
356 buildL0(*tev, router, true, back_inserter(dataL0));
357
358 summary.update(*tev);
359
361
362 if (event != NULL) {
363
365
366 for (const auto& t1 : event->mc_trks) {
368 if (t1.E > muon.
getE()) {
369
371
374
375 muon =
getFit(0, ta, 0.0, 0, t1.E, 1);
376
378 }
379 }
380 }
381 }
382
383 bool monte_carlo = false;
384 size_t index = 0;
385
386 for (bool next = false; !next; ) {
387
388 for (size_t i = 0; i != NUMBER_OF_PADS; ++i) {
389 H2[i].Reset();
390 }
391
393
394 if (!monte_carlo)
395 fit = (*in)[index];
396 else
397 fit = muon;
398
402
403
404
405
406
408
409
410
412
413 for (JDataL0_t::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {
414
415 JHitW0 hit(*i, summary.getRate(i->getPMTIdentifier(), parameters.
R_Hz));
416
417 hit.rotate(R);
418
419 if (match(hit)) {
421 }
422 }
423
424
425
426 sort(
data.begin(),
data.end(), JHitW0::compare);
427
428 JDataW0_t::iterator __end = unique(
data.begin(),
data.end(), equal_to<JDAQPMTIdentifier>());
429
430 double E_GeV = parameters.
E_GeV;
431
432
433
434
435
436
437
438
440
441 for (JDataW0_t::iterator hit =
data.begin(); hit != __end; ++hit) {
442
443 const double x = hit->getX() - tz.getX();
444 const double y = hit->getY() - tz.getY();
445 const double z = hit->getZ();
446 const double R = sqrt(x*x + y*y);
447
449 }
450
451 const double z0 = tz.getZ();
453
455
456
457
458 ostringstream os;
461
463
464 marker[2].push_back(TMarker(z0 - tz.getZ(), 0.0, kFullCircle));
466
467 static_cast<TAttMarker&>(marker[2][0]) = TAttMarker(kRed, kFullCircle, 0.7);
468 static_cast<TAttMarker&>(marker[2][1]) = TAttMarker(kRed, kFullCircle, 0.7);
469 }
470
472 <<
FIXED(7,2) << tz.getX() <<
' '
473 <<
FIXED(7,2) << tz.getY() <<
' '
474 <<
FIXED(7,2) << tz.getZ() <<
' '
475 <<
FIXED(12,2) << tz.getT() << endl);
476
477 double chi2 = 0;
478
479 for (JDataW0_t::const_iterator hit =
data.begin(); hit != __end; ++hit) {
480
481 const double x = hit->getX() - tz.getX();
482 const double y = hit->getY() - tz.getY();
483 const double z = hit->getZ() - tz.getZ();
484 const double R = sqrt(x*x + y*y);
485
487
488 JDirection3D dir(hit->getDX(), hit->getDY(), hit->getDZ());
489
491
492 const double theta = dir.getTheta();
493 const double phi = fabs(dir.getPhi());
494
495
496 const double E = E_GeV;
497 const double dt = T_ns.constrain(hit->getT() - t1);
498
501
504 }
505
506 H1 += H0;
507
508 chi2 += H1.getChi2() - H0.getChi2();
509
511 << setw(8) << hit->getModuleID() <<
'.' <<
FILL(2,
'0') << (
int) hit->getPMTAddress() <<
FILL() <<
' '
513 <<
FIXED(7,2) << R <<
' '
514 <<
FIXED(7,4) << theta <<
' '
515 <<
FIXED(7,4) << phi <<
' '
516 <<
FIXED(7,3) << dt <<
' '
517 <<
FIXED(7,3) << H1.getChi2() <<
' '
518 <<
FIXED(7,3) << H0.getChi2() << endl);
519
520 const double derivative = H1.getDerivativeOfChi2() - H0.getDerivativeOfChi2();
521
522 double size = derivative * graphics.arrowScale;
523
524 if (fabs(size) < Amin) {
525 size = (size > 0.0 ? +Amin : -Amin);
526 } else if (fabs(size) > Amax) {
527 size = (size > 0.0 ? +Amax : -Amax);
528 }
529
530 const double X[NUMBER_OF_PADS] = { R, atan2(y,x), z - R/
getTanThetaC() };
531
532 const double xs = (double) (NUMBER_OF_PMTS - 2 * hit->getPMTAddress()) / (double) NUMBER_OF_PMTS;
533
534 for (size_t i = 0; i != NUMBER_OF_PADS; ++i) {
535
536 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());
537
538 a1.SetLineWidth(graphics.lineWidth);
539 a1.SetLineStyle(graphics.lineStyle);
540
541 arrow[i].push_back(a1);
542
543 H2[i].Fill(X[i], dt + graphics.T_ns);
544 }
545 }
546
547 double zmax = 0.0;
548
549 for (size_t i = 0; i != NUMBER_OF_PADS; ++i) {
550 if (H2[i].GetMaximum() > zmax) {
551 zmax = H2[i].GetMaximum();
552 }
553 }
554
555 zmax *= 1.2;
556
557 for (size_t i = 0; i != NUMBER_OF_PADS; ++i) {
558 H2[i].SetMaximum(zmax);
559 }
560
562 os <<
" Q = " <<
FIXED(4,0) << fit.
getQ()
563 <<
'/' <<
FIXED(4,0) << -chi2;
565 os <<
" cos(#theta) = " <<
FIXED(6,3) << fit.
getDZ();
567
568 if (monte_carlo)
569 os << " Monte Carlo";
572
573
574
575
576 TLatex title(0.05, 0.5, os.str().c_str());
577
578 title.SetTextAlign(12);
579 title.SetTextFont(42);
580 title.SetTextSize(0.6);
581
582 p2->cd();
583
584 title.Draw();
585
586 for (int i = 0; i != NUMBER_OF_PADS; ++i) {
587
589
590 if (option == arrow_t) {
591
592 for (auto& a1 : arrow[i]) {
593 a1.Draw();
594 }
595
596 for (auto& m1 : marker[i]) {
597 m1.Draw();
598 }
599 }
600
601 if (option == histogram_t) {
602 H2[i].Draw("SAME");
603 }
604 }
605
606 cv->Update();
607
608
609
610
611 if (batch) {
612
614
615 next = true;
616
617 } else {
618
619 static int count = 0;
620
621 if (count++ == 0) {
622 cout << endl << "Type '?' for possible options." << endl;
623 }
624
625 for (bool user = true; user; ) {
626
628
630
631 ts.move(intersection.first);
632
633 cout << "\n> " << flush;
634
636
637 case '?':
638 cout << endl;
639 cout << "possible options: " << endl;
640 cout << 'p' << " -> " << "print information" << endl;
641 cout << 'q' << " -> " << "exit application" << endl;
642 cout << 'u' << " -> " << "update canvas" << endl;
643 cout << 's' << " -> " << "save graphics to file" << endl;
644 cout << '+' << " -> " << "next fit" << endl;
645 cout << '-' << " -> " << "previous fit" << endl;
646 cout << 'M' << " -> " << "Monte Carlo true muon information" << endl;
647 cout << 'F' << " -> " << "fit information" << endl;
649 cout << 'L' << " -> " << "reload event selector" << endl;
650 }
651 cout << 'r' << " -> " << "rewind input file" << endl;
652 cout << 'R' << " -> " << "switch to ROOT mode (quit ROOT to continue)" << endl;
653 cout << ' ' << " -> " << "next event (as well as any other key)" << endl;
654 break;
655
656 case 'p':
657
658 cout << endl;
659 cout <<
"intersection: " <<
FIXED(6,1) << intersection.first <<
' '<<
FIXED(6,1) << intersection.second << endl;
660 cout << "entry point: "
661 <<
FIXED(6,1) << ts.getX() - cylinder.getX() <<
' '
662 <<
FIXED(6,1) << ts.getY() - cylinder.getY() <<
' '
663 <<
FIXED(6,1) << ts.getZ() << endl;
666 }
667 break;
668
669 case 'q':
670 cout << endl;
671 return 0;
672
673 case 'u':
674 cv->Update();
675 break;
676
677 case 's':
679 break;
680
681 case '+':
682 monte_carlo = false;
683 index = (index != in->size() - 1 ? index + 1 : 0);
684 user = false;
685 break;
686
687 case '-':
688 monte_carlo = false;
689 index = (index != 0 ? index - 1 : in->size() - 1);
690 user = false;
691 break;
692
693 case 'M':
695 monte_carlo = true;
696 else
697 ERROR(endl <<
"No Monte Carlo muon available." << endl);
698 user = false;
699 break;
700
701 case 'F':
702 monte_carlo = false;
703 user = false;
704 break;
705
706 case 'L':
708 execute(
MAKE_STRING(
"make -f " <<
getPath(argv[0]) <<
"/JMakeEventSelector libs"), 3);
710 }
711 break;
712
713 case 'R':
714 tp->Run(kTRUE);
715 break;
716
717 case 'r':
718 inputFile.rewind();
719
720 default:
721 next = true;
722 user = false;
723 break;
724 }
725 }
726 }
727 }
728 }
729 }
730 cout << endl;
731}
#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
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.
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.
const std::vector< double > & getW() const
Get associated values.
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.
Object reading from a list of files.
File router for fast addressing of summary data.
Template definition for direct access of elements in ROOT TChain.
Enable unbuffered terminal input.
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 projected positions on the track of optical modules for which the response does not ...
JAxis3D getAxis(const Trk &track)
Get axis.
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.
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).
bool qualitySorter(const JFit &first, const JFit &second)
Comparison of fit results.
JRECONSTRUCTION::JWeight getWeight
Head getHeader(const JMultipleFileScanner_t &file_list)
Get Monte Carlo header.
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.
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.
Auxiliary class to match data points with given model.
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.
static counter_type max()
Get maximum counter value.
Auxiliary data structure for floating point format specification.
The Vec class is a straightforward 3-d vector, which also works in pyroot.