8 #include "TApplication.h"
102 return (i >= 0 && i < (
int) trk.
fitinf.size());
127 inline double getW(
const Trk& trk,
const int i,
const double value)
143 void setW(
Trk& trk,
const int i,
const double value)
145 if (i >= (
int) trk.
fitinf.size()) {
146 trk.
fitinf.resize(i + 1, 0.0);
158 const char WILDCARD =
'%';
173 istream
in(process.getInputStreamBuffer());
175 for (
string buffer;
getline(
in, buffer); ) {
176 DEBUG(buffer << endl);
189 int main(
int argc,
char **argv)
193 using namespace KM3NETDAQ;
214 JParser<> zap(
"Program to display hit probabilities.");
216 zap[
'w'] =
make_field(canvas,
"size of canvas <nx>x<ny> [pixels]") =
JCanvas(1200, 600);
217 zap[
'f'] =
make_field(inputFile,
"input file (output of JXXXMuonReconstruction.sh)");
218 zap[
'n'] =
make_field(numberOfEvents) = JLimit::max();
227 zap[
'B'] =
make_field(batch,
"batch processing");
232 catch(
const exception& error) {
233 FATAL(error.what() << endl);
237 FATAL(
"Missing output file name " <<
outputFile <<
" in batch mode." << endl);
244 if (
outputFile.find(WILDCARD) == string::npos) {
245 FATAL(
"Output file name " <<
outputFile <<
" has no wild card '" << WILDCARD <<
"'" << endl);
261 gROOT->SetBatch(batch);
263 TApplication* tp =
new TApplication(
"user", NULL, NULL);
264 TCanvas* cv =
new TCanvas(
"display",
"", canvas.x, canvas.y);
268 gROOT->SetStyle(
"gplot");
271 const size_t NUMBER_OF_PADS = 3;
273 cv->SetFillStyle(4000);
274 cv->SetFillColor(kWhite);
275 cv->Divide(NUMBER_OF_PADS, 1);
277 const double Dmax = 1000.0;
278 const double Rmin = 0.0;
279 const double Rmax = min(
parameters.roadWidth_m, 0.4 * Dmax);
280 const double Tmin = min(
parameters.TMin_ns, -10.0);
281 const double Tmax = max(
parameters.TMax_ns, +100.0);
282 const double Amin = 0.002 * (Tmax - Tmin);
283 const double Amax = 0.8 * (Tmax - Tmin);
284 const double ymin = min(-Amax, Tmin - 0.3 * Amax);
285 const double ymax = max(+Amax, Tmax + 0.5 * Amax);
287 const string Xlabel[NUMBER_OF_PADS] = {
"R [m]",
"#phi [rad]",
"z [m]" };
288 const double Xmin [NUMBER_OF_PADS] = { Rmin, -
PI, -0.5 * Dmax };
289 const double Xmax [NUMBER_OF_PADS] = { Rmax, +
PI, +0.5 * Dmax };
291 double Xs[NUMBER_OF_PADS];
293 for (
size_t i = 0;
i != NUMBER_OF_PADS; ++
i) {
297 TH2D H2[NUMBER_OF_PADS];
298 TGraph G2[NUMBER_OF_PADS];
300 for (
size_t i = 0;
i != NUMBER_OF_PADS; ++
i) {
302 H2[
i] = TH2D(
MAKE_CSTRING(
"h" <<
i), NULL, 100, Xmin[
i] - Xs[
i], Xmax[i] + Xs[i], 100, ymin, ymax);
304 H2[
i].GetXaxis()->SetTitle(Xlabel[i].c_str());
305 H2[
i].GetYaxis()->SetTitle(
"#Deltat [ns]");
307 H2[
i].GetXaxis()->CenterTitle(
true);
308 H2[
i].GetYaxis()->CenterTitle(
true);
310 H2[
i].SetStats(kFALSE);
314 G2[
i].SetPoint(0, H2[i].GetXaxis()->GetXmin(), 0.0);
315 G2[
i].SetPoint(1, H2[i].GetXaxis()->GetXmax(), 0.0);
324 while (inputFile.hasNext()) {
326 cout <<
"\revent: " << setw(8) << inputFile.getCounter() << flush;
328 Evt* evt = inputFile.next();
330 if (has_reconstructed_track<JPP_RECONSTRUCTION_TYPE>(*evt,
rec_stages_range(application))) {
332 Trk fit = get_best_reconstructed_track<JPP_RECONSTRUCTION_TYPE>(*evt,
rec_stages_range(application));
334 if (!event_selector(fit, *evt)) {
341 for (
const Hit& hit : evt->hits) {
358 for (
const auto& trk : evt->mc_trks) {
360 if (trk.E > muon.
E) {
372 bool next_event =
false;
373 bool monte_carlo =
false;
375 while (!next_event) {
398 for (JDataL0_t::const_iterator
i = dataL0.begin();
i != dataL0.end(); ++
i) {
411 sort(data.begin(), data.end(), JHitW0::compare);
413 JDataW0_t::iterator __end = unique(data.begin(), data.end(), equal_to<JDAQPMTIdentifier>());
426 const double z0 = tz.
getZ();
434 TLatex latex [NUMBER_OF_PADS];
438 for (
int i = 0;
i != NUMBER_OF_PADS; ++
i) {
440 latex[
i].SetTextAlign(12);
441 latex[
i].SetTextFont(42);
442 latex[
i].SetTextSize(0.06);
444 latex[
i].SetX(H2[
i].GetXaxis()->GetXmin() + 0.05 * (H2[
i].GetXaxis()->GetXmax() - H2[
i].GetXaxis()->GetXmin()));
445 latex[
i].SetY(H2[
i].GetYaxis()->GetXmax() - 0.05 * (H2[
i].GetYaxis()->GetXmax() - H2[
i].GetYaxis()->GetXmin()));
450 marker[2].push_back(TMarker(z0 - tz.
getZ(), 0.0, kFullCircle));
453 static_cast<TAttMarker&
>(marker[2][0]) = TAttMarker(kRed, kFullCircle, 0.7);
454 static_cast<TAttMarker&
>(marker[2][1]) = TAttMarker(kRed, kFullCircle, 0.7);
459 for (JDataW0_t::const_iterator hit = data.begin(); hit != __end; ++hit) {
461 const double x = hit->getX() - tz.
getX();
462 const double y = hit->getY() - tz.
getY();
463 const double z = hit->getZ() - tz.
getZ();
464 const double R = sqrt(x*x + y*y);
468 JDirection3D dir(hit->getDX(), hit->getDY(), hit->getDZ());
472 const double theta = dir.
getTheta();
473 const double phi = fabs(dir.getPhi());
476 const double E = E_GeV;
477 const double dt = T_ns.
constrain(hit->getT() - t1);
488 chi2 += H1.getChi2() - H0.getChi2();
490 const double derivative = H1.getDerivativeOfChi2() - H0.getDerivativeOfChi2();
492 double size = derivative * arrowScale;
494 if (fabs(size) < Amin) {
495 size = (size > 0.0 ? +Amin : -Amin);
496 }
else if (fabs(size) > Amax) {
497 size = (size > 0.0 ? +Amax : -Amax);
500 const double X[NUMBER_OF_PADS] = {
R, atan2(y,x), z - R/
getTanThetaC() };
504 for (
size_t i = 0;
i != NUMBER_OF_PADS; ++
i) {
505 arrow[
i].push_back(TArrow(X[
i] + xs*Xs[
i], dt, X[i] + xs*Xs[i], dt + size, arrowSize, arrowType.c_str()));
510 ":" << evt->frame_index <<
511 " " << evt->trigger_counter));
526 latex[2].SetTitle(
"Monte Carlo");
532 for (
int i = 0;
i != NUMBER_OF_PADS; ++
i) {
538 for (
auto& a1 : arrow[
i]) {
542 for (
auto& m1 : marker[i]) {
560 for (
int c = 0;
c == 0; ) {
564 static bool first =
true;
567 cout << endl <<
"Type '?' for possible options." << endl;
573 cout <<
"\n> " << flush;
579 cout <<
"possible options: " << endl;
580 cout <<
'q' <<
" -> " <<
"exit application" << endl;
581 cout <<
'u' <<
" -> " <<
"update canvas" << endl;
582 cout <<
's' <<
" -> " <<
"save graphics to file" << endl;
583 cout <<
'M' <<
" -> " <<
"Monte Carlo true muon information" << endl;
584 cout <<
'F' <<
" -> " <<
"fit information" << endl;
585 if (event_selector.is_valid()) {
586 cout <<
'L' <<
" -> " <<
"reload event selector" << endl;
588 cout <<
'r' <<
" -> " <<
"rewind input file" << endl;
589 cout <<
'R' <<
" -> " <<
"switch to ROOT mode (quit ROOT to continue)" << endl;
590 cout <<
'p' <<
" -> " <<
"print event information" << endl;
591 cout <<
' ' <<
" -> " <<
"next event (as well as any other key)" << endl;
613 ERROR(endl <<
"No Monte Carlo muon available." << endl);
626 if (event_selector.is_valid()) {
627 execute(
MAKE_STRING(
"make -f " <<
getPath(argv[0]) <<
"/JMakeEventSelector libs"), 3);
628 event_selector.reload();
641 for (
const auto& trk : evt->mc_trks) {
642 cout <<
"MC "; trk.
print(cout); cout << endl;
644 for (
const auto& trk : evt->trks) {
645 cout <<
"fit "; trk.
print(cout); cout << endl;
647 for (
const auto& hit : evt->hits) {
648 cout <<
"hit "; hit.
print(cout); cout << endl;
static const int JMUONSTART
Utility class to parse command line options.
double getAngle(const JQuaternion3D &first, const JQuaternion3D &second)
Get space angle between quanternions.
int main(int argc, char *argv[])
ROOT TTree parameter settings of various packages.
TString replace(const TString &target, const TRegexp ®exp, const T &replacement)
Replace regular expression in input by given replacement.
Data structure for direction in three dimensions.
then usage $script< input file >[option[primary[working directory]]] nWhere option can be E
double t
track time [ns] (when the particle is at pos )
JEventSelector()
Default constructor.
std::string getPath(const std::string &file_name)
Get path, i.e. part before last JEEP::PATHNAME_SEPARATOR if any.
void update(const JDAQHeader &header)
Update router.
bool is_muon(const Trk &track)
Test whether given track is a (anti-)muon.
void setW(Trk &trk, const int i, const double value)
Set associated value.
double getRate() const
Get default rate.
*fatal Wrong number of arguments esac JCookie sh typeset Z DETECTOR typeset Z SOURCE_RUN typeset Z TARGET_RUN set_variable PARAMETERS_FILE $WORKDIR parameters
#define MAKE_CSTRING(A)
Make C-string.
Template specialisation of class JModel to match hit with muon trajectory along z-axis.
double getZ(const JPosition3D &pos) const
Get point of emission of Cherenkov light along muon path.
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Auxiliary class to convert DAQ hit time to/from Monte Carlo hit time.
Auxiliary data structure for floating point format specification.
Data structure for UTC time.
double E
Energy [GeV] (either MC truth or reconstructed)
Range of reconstruction stages.
then echo The file $DIR KM3NeT_00000001_00000000 root already please rename or remove it first
#define MAKE_STRING(A)
Make string.
void print(std::ostream &out=std::cout) const
Print hit.
Basic data structure for L0 hit.
result_type calculate(const double E, const double R, const double theta, const double phi, const double t1) const
Get PDF.
The template JSinglePointer class can be used to hold a pointer to an object.
Scanning of objects from a single file according a format that follows from the extension of each fil...
double getW(const Trk &track, const size_t index, const double value)
Get track information.
Auxiliary class for defining the range of iterations of objects.
static const int JMUONPREFIT
I/O formatting auxiliaries.
JAxis3D & rotate(const JRotation3D &R)
Rotate axis.
JDirection3D getDirection(const Vec &dir)
Get direction.
JFunction1D_t::result_type result_type
JDirection3D & rotate(const JRotation3D &R)
Rotate.
Keyboard settings for unbuffered input.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Enable unbuffered terminal input.
const array_type< JValue_t > & make_array(const JValue_t(&array)[N])
Method to create array of values.
static const int JMUONGANDALF
double getTheta() const
Get theta angle.
JPosition3D getPosition(const Vec &pos)
Get position.
std::vector< double > fitinf
place to store additional fit info, see km3net-dataformat/definitions/fitparameters.csv
std::istream & getline(std::istream &in, JString &object)
Read string from input stream until end of line.
double len
length, if applicable [m]
double putTime() const
Get Monte Carlo time minus DAQ/trigger time.
static const double PI
Mathematical constants.
File router for fast addressing of summary data.
double getY() const
Get y position.
Auxiliary data structure for muon PDF.
Vec dir
hit direction; i.e. direction of the PMT
Data structure for fit parameters.
General purpose messaging.
Auxiliary data structure for sequence of same character.
Auxiliary include file for time conversion between DAQ/trigger hit and Monte Carlo hit...
Auxiliary data structure for muon PDF.
Streaming of input and output from Linux command.
$WORKDIR ev_configure_dqsimulator txt echo process $DQ_SIMULATOR $i $SOURCE_HOST[$index] csh c(setenv ROOTSYS $ROOTSYS &&source $JPP_DIR/setenv.csh $JPP_DIR &&($DQ_SIMULATOR\-u\$NAME\$\-H\$SERVER\$\-M\$LOGGER\$\-d $DEBUG</dev/null > &/dev/null &))'
then JCookie sh JDataQuality D $DETECTOR_ID R
Auxiliary class for a hit with background rate value.
const double getSpeedOfLight()
Get speed of light.
static const int JSTART_LENGTH_METRES
distance between first and last hits in metres from JStart.cc
bool has_muon(const Evt &evt)
Test whether given event has a muon.
double getT(const JVector3D &pos) const
Get arrival time of Cherenkov light at given position.
static const int JMUONSIMPLEX
Utility class to parse command line options.
double t
hit time (from tdc+calibration or MC truth)
Data structure for L0 hit.
const double getInverseSpeedOfLight()
Get inverse speed of light.
int dom_id
module identifier from the data (unique in the detector).
void setZ(const double z, const double velocity)
Set z-position of vertex.
then if[[!-f $DETECTOR]] then JDetector sh $DETECTOR fi cat $WORKDIR trigger_parameters txt<< EOFtrigger3DMuon.enabled=1;trigger3DMuon.numberOfHits=5;trigger3DMuon.gridAngle_deg=1;ctMin=0.0;TMaxLocal_ns=15.0;EOF set_variable TRIGGEREFFICIENCY_TRIGGERED_EVENTS_ONLY INPUT_FILES=() for((i=1;$i<=$NUMBER_OF_RUNS;++i));do JSirene.sh $DETECTOR $JPP_DATA/genhen.km3net_wpd_V2_0.evt.gz $WORKDIR/sirene_ ${i}.root JTriggerEfficiency.sh $DETECTOR $DETECTOR $WORKDIR/sirene_ ${i}.root $WORKDIR/trigger_efficiency_ ${i}.root $WORKDIR/trigger_parameters.txt $JPP_DATA/PMT_parameters.txt INPUT_FILES+=($WORKDIR/trigger_efficiency_ ${i}.root) done for ANGLE_DEG in $ANGLES_DEG[*];do set_variable SIGMA_NS 3.0 set_variable OUTLIERS 3 set_variable OUTPUT_FILE $WORKDIR/matrix\[${ANGLE_DEG}\deg\].root $JPP_DIR/examples/JReconstruction-f"$INPUT_FILES[*]"-o $OUTPUT_FILE-S ${SIGMA_NS}-A ${ANGLE_DEG}-O ${OUTLIERS}-d ${DEBUG}--!fiif[[$OPTION=="plot"]];then if((0));then for H1 in h0 h1;do JPlot1D-f"$WORKDIR/matrix["${^ANGLES_DEG}" deg].root:${H1}"-y"1 2e3"-Y-L TR-T""-\^"number of events [a.u.]"-> o chi2
Data structure for fit of straight line paralel to z-axis.
unsigned int tot
tot value as stored in raw data (int for pyroot)
double getX() const
Get x position.
no fit printf nominal n $STRING awk v X
double getTanThetaC()
Get average tangent of Cherenkov angle of water corresponding to group velocity.
unsigned int channel_id
PMT channel id {0,1, .., 30} local to moduke.
Object reading from a list of files.
const JLimit & getLimit() const
Get limit.
void print(std::ostream &out=std::cout) const
Print track.
KM3NeT DAQ constants, bit handling, etc.
static const int NUMBER_OF_PMTS
Total number of PMTs in module.
then fatal Wrong number of arguments fi set_variable DETECTOR $argv[1] set_variable INPUT_FILE $argv[2] eval JPrintDetector a $DETECTOR O IDENTIFIER eval JPrintDetector a $DETECTOR O SUMMARY JAcoustics sh $DETECTOR_ID source JAcousticsToolkit sh CHECK_EXIT_CODE typeset A EMITTERS get_tripods $WORKDIR tripod txt EMITTERS get_transmitters $WORKDIR transmitter txt EMITTERS for EMITTER in
The Trk class represents a Monte Carlo (MC) particle as well as a reconstructed track/shower.
Auxiliary data structure for floating point format specification.
Wrapper class around ROOT TStyle.
static bool select(const Trk &trk, const Evt &evt)
Default event selection.
static const int JMUONENERGY
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
Data structure for size of TCanvas.
#define DEBUG(A)
Message macros.
bool hasW(const Trk &trk, const int i)
Check availability of value.