107{
111
113 typedef JParallelFileScanner_t::multi_pointer_type multi_pointer_type;
117
118 JParallelFileScanner_t inputFile;
120 string detectorFile;
121 JCalibration_t calibrationFile;
122 double Tmax_s;
123 string pdfFile;
126 int application;
129 bool batch;
131 double arrowSize = 0.003;
132 string arrowType = "|->";
133 double arrowScale = 250.0;
134 Width_t lineWidth = 2;
135 Style_t lineStyle = 1;
136 int nbinsX = 50;
137 int nbinsY = 250;
138 double T_ns = 0.0;
139 } graphics;
140 string option;
142
143
144 try {
145
147
154
156
157 JParser<> zap(
"Program to display hit probabilities.");
158
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)");
171 zap[
'O'] =
make_field(option,
"draw option") = arrow_t, histogram_t;
172 zap[
'B'] =
make_field(batch,
"batch processing");
174
175 zap(argc, argv);
176 }
177 catch(const exception& error) {
178 FATAL(error.what() << endl);
179 }
180
182 FATAL(
"Missing output file name " <<
outputFile <<
" in batch mode." << endl);
183 }
184
187 }
188
191 }
192
193
195
196 try {
198 }
201 }
202 unique_ptr<JDynamics> dynamics;
203
204 try {
205
207
208 dynamics->load(calibrationFile);
209 }
210 catch(const exception& error) {
211 if (!calibrationFile.empty()) {
213 }
214 }
215
217
219
221
223
226
228
229
230 Vec offset(0.0, 0.0, 0.0);
231
232 try {
234 } catch(const exception& error) {}
235
236
237
238
239 gROOT->SetBatch(batch);
240
241 TApplication* tp = new TApplication("user", NULL, NULL);
242 TCanvas* cv =
new TCanvas(
"display",
"", canvas.
x, canvas.
y);
243
244 unique_ptr<TStyle> gStyle(
new JStyle(
"gplot", cv->GetWw(), cv->GetWh(), graphics));
245
246 gROOT->SetStyle("gplot");
247 gROOT->ForceStyle();
248
249 const size_t NUMBER_OF_PADS = 3;
250
251 cv->SetFillStyle(4000);
252 cv->SetFillColor(kWhite);
253
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);
256
257 p1->Divide(NUMBER_OF_PADS, 1);
258
260 p2->Draw();
261
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);
271
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 };
275
276 double Xs[NUMBER_OF_PADS];
277
278 for (size_t i = 0; i != NUMBER_OF_PADS; ++i) {
279 Xs[i] = 0.003 * (Xmax[i] - Xmin[i]) * (0.5 * NUMBER_OF_PMTS);
280 }
281
282 TH2D H2[NUMBER_OF_PADS];
283 TGraph G2[NUMBER_OF_PADS];
284
285 for (size_t i = 0; i != NUMBER_OF_PADS; ++i) {
286
287 H2[i] = TH2D(
MAKE_CSTRING(
"h" << i), NULL, graphics.nbinsX, Xmin[i] - Xs[i], Xmax[i] + Xs[i], graphics.nbinsY, ymin, ymax);
288
289 H2[i].GetXaxis()->SetTitle(Xlabel[i].c_str());
290 H2[i].GetYaxis()->SetTitle("#Deltat [ns]");
291
292 H2[i].GetXaxis()->CenterTitle(true);
293 H2[i].GetYaxis()->CenterTitle(true);
294
295 H2[i].SetStats(kFALSE);
296
297 G2[i].Set(2);
298
299 G2[i].SetPoint(0, H2[i].GetXaxis()->GetXmin(), 0.0);
300 G2[i].SetPoint(1, H2[i].GetXaxis()->GetXmax(), 0.0);
301
303
304 H2[i].Draw("AXIS");
305 G2[i].Draw("SAME");
306 }
307
308
310
311 cout << "event: " << setw(8) << inputFile.getCounter() << endl;
312
313 multi_pointer_type ps = inputFile.next();
314
318
319 if (dynamics) {
320 dynamics->update(*tev);
321 }
322
323 if (mc.getEntries() != 0) {
325 }
326
328
329 if (!in->empty()) {
330
332
333 if (!event_selector(*tev, *in, event)) {
334 continue;
335 }
336
337
338 JDataL0_t dataL0;
339
340 buildL0(*tev, router, true, back_inserter(dataL0));
341
342 summary.update(*tev);
343
345
346 if (event != NULL) {
347
349
350 for (const auto& t1 : event->mc_trks) {
352 if (t1.E > muon.
getE()) {
353
355
358
359 muon =
getFit(0, ta, 0.0, 0, t1.E, 1);
360
362 }
363 }
364 }
365 }
366
367 bool monte_carlo = false;
368 size_t index = 0;
369
370 for (bool next = false; !next; ) {
371
372 for (size_t i = 0; i != NUMBER_OF_PADS; ++i) {
373 H2[i].Reset();
374 }
375
377
378 if (!monte_carlo)
379 fit = (*in)[index];
380 else
381 fit = muon;
382
386
387
388
389
390
392
393
394
396
397 for (JDataL0_t::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {
398
399 JHitW0 hit(*i, summary.getRate(i->getPMTIdentifier()));
400
401 hit.rotate(R);
402
403 if (match(hit)) {
405 }
406 }
407
408
409
410 sort(
data.begin(),
data.end(), JHitW0::compare);
411
412 JDataW0_t::iterator __end = unique(
data.begin(),
data.end(), equal_to<JDAQPMTIdentifier>());
413
414 double E_GeV = parameters.
E_GeV;
415
416
417
418
419
420
421
422
424
425 for (JDataW0_t::iterator hit =
data.begin(); hit != __end; ++hit) {
426
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);
431
433 }
434
435 const double z0 = tz.getZ();
437
439
440
441
442 ostringstream os;
445
447
448 marker[2].push_back(TMarker(z0 - tz.getZ(), 0.0, kFullCircle));
450
451 static_cast<TAttMarker&>(marker[2][0]) = TAttMarker(kRed, kFullCircle, 0.7);
452 static_cast<TAttMarker&>(marker[2][1]) = TAttMarker(kRed, kFullCircle, 0.7);
453 }
454
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);
460
461 double chi2 = 0;
462
463 for (JDataW0_t::const_iterator hit =
data.begin(); hit != __end; ++hit) {
464
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);
469
471
472 JDirection3D dir(hit->getDX(), hit->getDY(), hit->getDZ());
473
475
476 const double theta = dir.getTheta();
477 const double phi = fabs(dir.getPhi());
478
479
480 const double E = E_GeV;
481 const double dt = T_ns.constrain(hit->getT() - t1);
482
485
488 }
489
490 H1 += H0;
491
492 chi2 += H1.getChi2() - H0.getChi2();
493
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);
503
504 const double derivative = H1.getDerivativeOfChi2() - H0.getDerivativeOfChi2();
505
506 double size = derivative * graphics.arrowScale;
507
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);
512 }
513
514 const double X[NUMBER_OF_PADS] = { R, atan2(y,x), z - R/
getTanThetaC() };
515
516 const double xs = (double) (NUMBER_OF_PMTS - 2 * hit->getPMTAddress()) / (double) NUMBER_OF_PMTS;
517
518 for (size_t i = 0; i != NUMBER_OF_PADS; ++i) {
519
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());
521
522 a1.SetLineWidth(graphics.lineWidth);
523 a1.SetLineStyle(graphics.lineStyle);
524
525 arrow[i].push_back(a1);
526
527 H2[i].Fill(X[i], dt + graphics.T_ns);
528 }
529 }
530
532 os <<
" Q = " <<
FIXED(4,0) << fit.
getQ();
534 os <<
" cos(#theta) = " <<
FIXED(6,3) << fit.
getDZ();
535
536 if (monte_carlo)
537 os << " Monte Carlo";
540
541
542
543
544 TLatex title(0.05, 0.5, os.str().c_str());
545
546 title.SetTextAlign(12);
547 title.SetTextFont(42);
548 title.SetTextSize(0.6);
549
550 p2->cd();
551
552 title.Draw();
553
554 for (int i = 0; i != NUMBER_OF_PADS; ++i) {
555
557
558 if (option == arrow_t) {
559
560 for (auto& a1 : arrow[i]) {
561 a1.Draw();
562 }
563
564 for (auto& m1 : marker[i]) {
565 m1.Draw();
566 }
567 }
568
569 if (option == histogram_t) {
570 H2[i].Draw("SAME");
571 }
572 }
573
574 cv->Update();
575
576
577
578
579 if (batch) {
580
582
583 next = true;
584
585 } else {
586
587 static int count = 0;
588
589 if (count++ == 0) {
590 cout << endl << "Type '?' for possible options." << endl;
591 }
592
593 for (bool user = true; user; ) {
594
595 cout << "\n> " << flush;
596
598
599 case '?':
600 cout << endl;
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;
611 }
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;
615 break;
616
617 case 'q':
618 cout << endl;
619 return 0;
620
621 case 'u':
622 cv->Update();
623 break;
624
625 case 's':
627 break;
628
629 case '+':
630 monte_carlo = false;
631 index = (index != in->size() - 1 ? index + 1 : 0);
632 user = false;
633 break;
634
635 case '-':
636 monte_carlo = false;
637 index = (index != 0 ? index - 1 : in->size() - 1);
638 user = false;
639 break;
640
641 case 'M':
643 monte_carlo = true;
644 else
645 ERROR(endl <<
"No Monte Carlo muon available." << endl);
646 user = false;
647 break;
648
649 case 'F':
650 monte_carlo = false;
651 user = false;
652 break;
653
654 case 'L':
656 execute(
MAKE_STRING(
"make -f " <<
getPath(argv[0]) <<
"/JMakeEventSelector libs"), 3);
658 }
659 break;
660
661 case 'R':
662 tp->Run(kTRUE);
663 break;
664
665 case 'r':
666 inputFile.rewind();
667
668 default:
669 next = true;
670 user = false;
671 break;
672 }
673 }
674 }
675 }
676 }
677 }
678 cout << endl;
679}
#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.
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 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.
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.
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.