JGandalf.cc File Reference

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <set>
#include <algorithm>
#include "TROOT.h"
#include "TFile.h"
#include "TProfile.h"
#include "TH2D.h"
#include "km3net-dataformat/offline/Head.hh"
#include "km3net-dataformat/offline/Evt.hh"
#include "JDAQ/JDAQEventIO.hh"
#include "JDAQ/JDAQTimesliceIO.hh"
#include "JDAQ/JDAQSummarysliceIO.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JModuleRouter.hh"
#include "JLang/JPredicate.hh"
#include "JTrigger/JHit.hh"
#include "JTrigger/JFrame.hh"
#include "JTrigger/JTimeslice.hh"
#include "JTrigger/JHitL0.hh"
#include "JTrigger/JBuildL0.hh"
#include "JTrigger/JTriggerParameters.hh"
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JParallelFileScanner.hh"
#include "JSupport/JFileRecorder.hh"
#include "JSupport/JSupport.hh"
#include "JSupport/JMeta.hh"
#include "JSupport/JSummaryFileRouter.hh"
#include "JFit/JHitW0.hh"
#include "JFit/JPMTW0.hh"
#include "JFit/JLine3Z.hh"
#include "JFit/JGandalf.hh"
#include "JFit/JLine3ZRegressor.hh"
#include "JFit/JFitToolkit.hh"
#include "JFit/JEvt.hh"
#include "JFit/JEvtToolkit.hh"
#include "JFit/JRegressor.hh"
#include "JFit/JFitParameters.hh"
#include "JFit/JModel.hh"
#include "JTools/JConstants.hh"
#include "JTools/JFunction1D_t.hh"
#include "JTools/JFunctionalMap_t.hh"
#include "JTools/JAbstractHistogram.hh"
#include "JROOT/JRootToolkit.hh"
#include "JGizmo/JManager.hh"
#include "JPhysics/JPDFTable.hh"
#include "JPhysics/JPDFToolkit.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

Go to the source code of this file.

Functions
int	main (int argc, char **argv)

Detailed Description

Program to perform JFIT::JRegressor<JLine3Z,JGandalf> fit with I/O of JFIT::JEvt data.

This program can also be used for detector calibrations. To this end, the output from a previous run of this program can be used in conjunction with option -c. It is recommended to accordingly set option -N 1 so that only the best track from the previous output is used. Several histograms will be produced which can subsequently be processed with JHobbit.cc for the determination of the time calibration of optical modules and JFrodo.cc for the time-slewing correction.

Author

mdejong

Definition in file JGandalf.cc.

Function Documentation

main()

int main	(	int	argc,
		char **	argv
	)

Definition at line 98 of file JGandalf.cc.

{
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;
  
  typedef JParallelFileScanner< JTypeList<JDAQEvent, JEvt> >       JParallelFileScanner_t;
  typedef JParallelFileScanner_t::multi_pointer_type               multi_pointer_type;
  typedef JTYPELIST<JAllTypes_t, JEvt, TObject>::typelist          typelist;
  typedef JAbstractHistogram<double>                               histogram_type;
 
  JParallelFileScanner_t  inputFile;
  JFileRecorder<typelist> outputFile;
  JLimit_t&      numberOfEvents = inputFile.getLimit();
  string         detectorFile;
  string         pdfFile;
  double         roadWidth_m;
  double         R_Hz;
  size_t         numberOfPrefits;
  double         TTS_ns;
  JTimeRange     T_compress;
  double         E_GeV;
  JZRange        z;
  bool           reprocess;
  histogram_type calibrate;
  int            debug;
 
  try { 
 
    JParser<> zap("Program to perform PDF fit of muon trajectory to data.");
    
    zap['f'] = make_field(inputFile);
    zap['o'] = make_field(outputFile)          = "gandalf.root";
    zap['a'] = make_field(detectorFile);
    zap['n'] = make_field(numberOfEvents)      = JLimit::max();
    zap['P'] = make_field(pdfFile);
    zap['R'] = make_field(roadWidth_m)         = numeric_limits<double>::max();
    zap['B'] = make_field(R_Hz,      "Default background rate if no summary data are available.")  = 6.0e3;
    zap['N'] = make_field(numberOfPrefits)     = 0;
    zap['T'] = make_field(TTS_ns,    "Time smearing of PDF.");
    zap['C'] = make_field(T_compress)          = JTimeRange();
    zap['E'] = make_field(E_GeV,     "Default energy if no energy estimate is available.")         = 1.0e3;
    zap['z'] = make_field(z,         "Extension of longitudinal range of muon track.")             = JZRange(0.0, 0.0);
    zap['r'] = make_field(reprocess);
    zap['c'] = make_field(calibrate, "Histogram for time calibration per optical module.")         = histogram_type();
    zap['d'] = make_field(debug)               = 1;
    
    zap(argc, argv);
  }
  catch(const exception& error) {
    FATAL(error.what() << endl);
  }
 
  cout.tie(&cerr);
 
  JDetector detector;
 
  try {
    load(detectorFile, detector);
  }
  catch(const JException& error) {
    FATAL(error);
  }
 
  const JModuleRouter moduleRouter (detector);
  JSummaryFileRouter  summaryRouter(inputFile, R_Hz);
 
 
  typedef vector<JHitL0>  JDataL0_t;
  typedef vector<JHitW0>  JDataW0_t;
  const JBuildL0<JHitL0>  buildL0;
 
 
  typedef JRegressor<JLine3Z, JGandalf> JRegressor_t;
 
  JRegressor_t::debug              = debug;
  JRegressor_t::T_ns.setRange(-50.0, +450.0);   // ns
  JRegressor_t::Vmax_npe           =   10.0;    // npe
  JRegressor_t::MAXIMUM_ITERATIONS = 1000;
 
  JRegressor_t fit(pdfFile, TTS_ns);
  
  if (T_compress != JTimeRange()) {
    fit.compress(T_compress);
  }
 
  fit.parameters.resize(5);
 
  fit.parameters[0] = JLine3Z::pX();
  fit.parameters[1] = JLine3Z::pY();
  fit.parameters[2] = JLine3Z::pT();
  fit.parameters[3] = JLine3Z::pDX();
  fit.parameters[4] = JLine3Z::pDY();
 
  if (fit.getRmax() < roadWidth_m) {
    roadWidth_m = fit.getRmax();
  }
 
  TH2D*     ha = NULL;
  TProfile* hb = NULL;
  TH2D*     h2 = NULL;
 
  typedef JManager<string, TProfile>  JManager_t;
 
  JManager_t gd;
 
  if (calibrate.is_valid()) {
 
    if (numberOfPrefits != 1) {
      WARNING("Running in calibration mode but number of prefits " << numberOfPrefits << " != " << 1 << endl);
    }
 
    ha = new TH2D    ("ha", NULL, 256, -0.5, 255.5, 
                      calibrate.getNumberOfBins(), calibrate.getLowerLimit(), calibrate.getUpperLimit());
    hb = new TProfile("hb", NULL, 256, -0.5, 255.5);
 
    h2 = new TH2D    ("h2", NULL, detector.size(), -0.5, detector.size() - 0.5,
                      calibrate.getNumberOfBins(), calibrate.getLowerLimit(), calibrate.getUpperLimit());
 
    gd = JManager_t(new TProfile("%", NULL, detector.size(), -0.5, detector.size() - 0.5));
  }
 
  outputFile.open();
 
  outputFile.put(JMeta(argc, argv));
 
  while (inputFile.hasNext()) {
 
    STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);
 
    multi_pointer_type ps = inputFile.next();
 
    JDAQEvent* tev = ps;
    JEvt*      evt = ps;
    JEvt       out;
 
    summaryRouter.update(*tev);
 
    JEvt::iterator __end = evt->end();
 
    if (reprocess) {
      __end = partition(evt->begin(), __end, JHistory::is_not_event(JMUONGANDALF));
    }
    
    if (evt->begin() != __end) {
 
      copy(evt->begin(), __end, back_inserter(out));
    
      if (numberOfPrefits > 0) {
        advance(__end = evt->begin(), min(numberOfPrefits, out.size()));
      }
 
      partial_sort(evt->begin(), __end, evt->end(), qualitySorter);
      
      __end = partition(evt->begin(), __end, JHistory::is_event(evt->begin()->getHistory()));
 
 
      JDataL0_t dataL0;
      
      buildL0(*tev, moduleRouter, true, back_inserter(dataL0));
      
 
      if (dataL0.size() >= fit.parameters.size()) {
 
        for (JEvt::const_iterator track = evt->begin(); track != __end; ++track) {
 
          const JRotation3D  R (getDirection(*track));
          const JLine1Z      tz(getPosition (*track).rotate(R), track->getT());
          JZRange            Z;                      // default infinite track length for hit selection
      
          if (track->hasW(JSTART_LENGTH_METRES)) {   // use output of JStart if available
            Z = JZRange(0.0, track->getW(JSTART_LENGTH_METRES))  +  z; 
          }
 
          const JModel<JLine1Z> match(tz, roadWidth_m, JRegressor_t::T_ns, Z);
 
          // hit selection based on start value
 
          JDataW0_t data;
          
          for (JDataL0_t::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {
            
            double rate_Hz = summaryRouter.getRate(*i);
 
            if (rate_Hz <= 0.0) {
              rate_Hz = summaryRouter.getRate();
            }
 
            JHitW0 hit(*i, rate_Hz);
            
            hit.rotate(R);
 
            if (match(hit)) {
              data.push_back(hit);
            }
          }
          
          // select first hit
          
          sort(data.begin(), data.end(), compare);
          
          JDataW0_t::iterator __end = unique(data.begin(), data.end(), equal_to<JDAQPMTIdentifier>());
          
 
          const int NDF = distance(data.begin(), __end) - fit.parameters.size();
 
          if (NDF >= 0) {
 
            // set fit parameters
 
            if (track->getE() > 0.1)
              fit.E_GeV   = track->getE();
            else
              fit.E_GeV   = E_GeV;
 
            const double chi2 = fit(JLine3Z(tz), data.begin(), __end);
            
            JTrack3D tb(fit.value);
            
            tb.rotate_back(R);
 
            out.push_back(getFit(JHistory(track->getHistory()).add(JMUONGANDALF), tb, getQuality(chi2), NDF));
 
            // set additional values
 
            out.rbegin()->setW(track->getW());
            out.rbegin()->setW(JGANDALF_BETA0_RAD,            sqrt(fit.error.getDX() * fit.error.getDX()  +
                                                                   fit.error.getDY() * fit.error.getDY()));
            out.rbegin()->setW(JGANDALF_BETA1_RAD,            sqrt(fit.error.getDX() * fit.error.getDY()));
            out.rbegin()->setW(JGANDALF_CHI2,                 chi2);
            out.rbegin()->setW(JGANDALF_NUMBER_OF_HITS,       distance(data.begin(), __end));
            out.rbegin()->setW(JGANDALF_LAMBDA,               fit.lambda);
            out.rbegin()->setW(JGANDALF_NUMBER_OF_ITERATIONS, fit.numberOfIterations);
 
            if (calibrate.is_valid()) {
 
              set<int> buffer;
 
              for (JDataW0_t::const_iterator hit = data.begin(); hit != __end; ++hit) {
                buffer.insert(hit->getModuleID());
              }
 
              const JLine3Z ta = fit.value;
 
              for (set<int>::const_iterator id = buffer.begin(); id != buffer.end(); ++id) {
 
                JDataW0_t::iterator q = partition(data.begin(), __end, make_predicate(&JHitW0::getModuleID, *id, JComparison::ne()));
 
                if (distance(data.begin(), q) - fit.parameters.size() > 0) {
 
                  fit(ta, data.begin(), q);     // re-fit
 
                  sort(q, __end, compare);      // select first hit in module
 
                  const int    index    = moduleRouter.getIndex(*id);
                  const double t1       = fit.value.getT(q->getPosition());
                  JTrack3D     gradient = fit(fit.value, *q).gradient;
 
                  gradient.rotate_back(R);
 
                  ha->Fill(q->getToT(), q->getT1() - t1);
                  hb->Fill(q->getToT(), gradient.getT());
 
                  h2->Fill(index, q->getT() - t1);
 
                  gd["T"]->Fill(index, gradient.getT());
                  gd["X"]->Fill(index, gradient.getX());
                  gd["Y"]->Fill(index, gradient.getY());
                  gd["Z"]->Fill(index, gradient.getZ());
                }
              }
            }
          }
        }
      }
 
      // apply default sorter
 
      sort(out.begin(), out.end(), qualitySorter);
    }
 
    outputFile.put(out);
  }
  STATUS(endl);
 
  if (calibrate.is_valid()) {
 
    outputFile.put(*ha);
    outputFile.put(*hb);
    outputFile.put(*h2);
 
    for (JManager_t::const_iterator i = gd.begin(); i != gd.end(); ++i) {
      outputFile.put(*(i->second));
    }
  }
 
  JSingleFileScanner<JRemove<typelist, JEvt>::typelist> io(inputFile);
 
  io >> outputFile;
 
  outputFile.close();
}