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.