JCalibrateMuon.cc File Reference

Program to perform detector calibration using reconstructed muon trajectories. More...

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <set>
#include <algorithm>
#include <memory>
#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 "JDetector/JPMTParametersMap.hh"
#include "JDynamics/JDynamics.hh"
#include "JTrigger/JHitL0.hh"
#include "JTrigger/JBuildL0.hh"
#include "JSupport/JSupport.hh"
#include "JSupport/JParallelFileScanner.hh"
#include "JSupport/JSummaryFileRouter.hh"
#include "JSupport/JMeta.hh"
#include "JReconstruction/JEvt.hh"
#include "JReconstruction/JMuonGandalf.hh"
#include "JReconstruction/JMuonGandalfParameters_t.hh"
#include "JTools/JAbstractHistogram.hh"
#include "JROOT/JRootToolkit.hh"
#include "JROOT/JManager.hh"
#include "JCalibrate/JCalibrateMuon.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 detector calibration using reconstructed muon trajectories.

It is recommended to perform the detector calibration after JMuonGandalf.cc and to accordingly set option JMuonGandalfParameters_t::numberOfPrefits = 1, so that 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 JCalibrateMuon.cc.

Function Documentation

main()

int main	(	int	argc,
		char **	argv )

Definition at line 62 of file JCalibrateMuon.cc.

{
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;
  
  typedef JParallelFileScanner< JTypeList<JDAQEvent, JFIT::JEvt> > JParallelFileScanner_t;
  typedef JParallelFileScanner_t::multi_pointer_type               multi_pointer_type;
  typedef JMultipleFileScanner<calibration_types>                  JCalibration_t;
  typedef JAbstractHistogram<double>                               histogram_type;
 
  JParallelFileScanner_t   inputFile;
  string                   outputFile;
  JLimit_t&                numberOfEvents = inputFile.getLimit();
  string                   detectorFile;
  JCalibration_t           calibrationFile;
  double                   Tmax_s;
  string                   pdfFile;
  JMuonGandalfParameters_t parameters;
  JPMTParametersMap        pmtParameters;
  histogram_type           calibrate;
  string                   target;
  int                      debug;
 
  try { 
 
    parameters.numberOfPrefits = 1;
 
    JParser<> zap("Program to perform detector calibration using reconstructed muon trajectories.");
    
    zap['f'] = make_field(inputFile);
    zap['o'] = make_field(outputFile)          = "gandalf.root";
    zap['a'] = make_field(detectorFile);
    zap['+'] = make_field(calibrationFile)     = JPARSER::initialised();
    zap['T'] = make_field(Tmax_s)              = 100.0;
    zap['n'] = make_field(numberOfEvents)      = JLimit::max();
    zap['F'] = make_field(pdfFile);
    zap['c'] = make_field(calibrate, "histogram for time calibration.");
    zap['P'] = make_field(pmtParameters,       "PMT simulation data (or corresponding file name)")                 = JPARSER::initialised();
    zap['@'] = make_field(parameters)          = JPARSER::initialised();
    zap['R'] = make_field(target,    "option for time calibration.")  =  module_t, string_t;
    zap['d'] = make_field(debug)               = 1;
    
    zap(argc, argv);
  }
  catch(const exception& error) {
    FATAL(error.what() << endl);
  }
 
 
  if (!calibrate.is_valid()) {
    FATAL("Invalid calibration data " << calibrate << endl);
  }
 
  if (parameters.numberOfPrefits != 1) {
    WARNING("Number of prefits " << parameters.numberOfPrefits << " != " << 1 << endl);
  }
 
  JDetector detector;
 
  try {
    load(detectorFile, detector);
  }
  catch(const JException& error) {
    FATAL(error);
  }
 
  unique_ptr<JDynamics> dynamics;
 
  if (!calibrationFile.empty()) {
 
    try {
 
      dynamics.reset(new JDynamics(detector, Tmax_s));
 
      dynamics->load(calibrationFile);
    }
    catch(const exception& error) {
      FATAL(error.what());
    }
  }
 
  const JRouter_t* rpm = NULL;
 
  if      (target == module_t)
    rpm = new JModuleRouter_t(detector);
  else if (target == string_t)
    rpm = new JStringRouter_t(detector);
  else
    FATAL("No valid target; check option -R" << endl);
 
  const JModuleRouter router(dynamics ? dynamics->getDetector() : detector);
 
  JSummaryFileRouter  summary(inputFile);
 
  const JMuonGandalf::storage_type storage(pdfFile, JTimeRange(parameters.TMin_ns, parameters.TMax_ns), parameters.TTS_ns);
 
  JMuonGandalf fit(parameters, storage, pmtParameters, debug);
 
  typedef vector<JHitL0>  JDataL0_t;
  typedef vector<JHitW0>  JDataW0_t;
 
  const JBuildL0<JHitL0>  buildL0;
 
 
  typedef JManager<string, TProfile>  JManager_t;
 
  TH2D       ha("ha", NULL, 256, -0.5, 255.5, 
                calibrate.getNumberOfBins(), calibrate.getLowerLimit(), calibrate.getUpperLimit());
 
  TProfile   hb("hb", NULL, 256, -0.5, 255.5);
 
  TProfile   hr("hr", NULL,  60,  0.0, 150.0);
 
  TH2D       h2("h2", NULL, rpm->getN(), -0.5, rpm->getN() - 0.5,
                calibrate.getNumberOfBins(), calibrate.getLowerLimit(), calibrate.getUpperLimit());
 
  JManager_t g1(new TProfile("%", NULL, rpm->getN(), -0.5, rpm->getN() - 0.5, -1.0, +1.0));
 
 
  while (inputFile.hasNext()) {
 
    STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);
 
    multi_pointer_type ps = inputFile.next();
 
    JDAQEvent*  tev = ps;
    JFIT::JEvt* in  = ps;
 
    summary.update(*tev);
 
    if (dynamics) {
      dynamics->update(*tev);
    }
 
    in->select(parameters.numberOfPrefits, qualitySorter);
 
    JDataL0_t dataL0;
      
    buildL0(*tev, router, true, back_inserter(dataL0));
      
    for (JFIT::JEvt::const_iterator track = in->begin(); track != in->end(); ++track) {
 
      const JRotation3D  R (getDirection(*track));
      const JLine1Z      tz(getPosition (*track).rotate(R), track->getT());
      JRange<double>     Z_m;
 
      if (track->hasW(JSTART_LENGTH_METRES)) {
        Z_m = JZRange(0.0, track->getW(JSTART_LENGTH_METRES))  +  JRange<double>(parameters.ZMin_m, parameters.ZMax_m);
      }
 
      const JFIT::JModel<JLine1Z> match(tz, parameters.roadWidth_m, JRange<double>(parameters.TMin_ns, parameters.TMax_ns), Z_m);
 
      // hit selection based on start value
 
      JDataW0_t data;
          
      for (JDataL0_t::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {
 
        const JPMTIdentifier id(i->getModuleID(), i->getPMTAddress());
 
        const JPMTParameters& wip = pmtParameters.getPMTParameters(id);
 
        const int    type = wip.getType();
        const double QE   = wip.QE;
        const double R_Hz = summary.getRate(i->getPMTIdentifier(), parameters.R_Hz);
 
        JHitW0 hit(*i, type, QE, R_Hz);
            
        hit.rotate(R);
 
        if (match(hit)) {
          data.push_back(hit);
        }
      }
          
      // select first hit
          
      sort(data.begin(), data.end(), JHitL0::compare);
 
      JDataW0_t::iterator __end = unique(data.begin(), data.end(), equal_to<JDAQPMTIdentifier>());
          
 
      // set fit parameters
 
      if (track->getE() > 0.1)
        fit.JRegressor_t::E_GeV   = track->getE();
      else
        fit.JRegressor_t::E_GeV   = parameters.E_GeV;
 
      set<int> buffer;
 
      for (JDataW0_t::const_iterator hit = data.begin(); hit != __end; ++hit) {
        buffer.insert(rpm->getIndex(hit->getModuleID()));
      }
 
      const JLine3Z ta(tz);
 
      for (const int index : buffer) {
 
        JDataW0_t::iterator q = partition(data.begin(), __end, [rpm, index](const JHitW0& hit) { return rpm->getIndex(hit.getModuleID()) != index; });
 
        if (distance(data.begin(), q) - fit.parameters.size() > 0) {
 
          fit(ta, data.begin(), q);                   // re-fit
 
          for (JDataW0_t::const_iterator hit = q; hit != __end; ++hit) {
 
            const double t1       = fit.value.getT(hit->getPosition());
            JTrack3D     gradient = fit(fit.value, *hit).gradient;
 
            gradient.rotate_back(R);
 
            ha.Fill(hit->getToT(), hit->getT1() - t1);
            hb.Fill(hit->getToT(), gradient.getT());
 
            hr.Fill(fit.value.getDistance(*hit), gradient.getT());
 
            h2.Fill(index, hit->getT() - t1);
 
            g1["T"]->Fill(index, gradient.getT());
            g1["X"]->Fill(index, gradient.getX());
            g1["Y"]->Fill(index, gradient.getY());
            g1["Z"]->Fill(index, gradient.getZ());
          }
        }
      }
    }
  }
  STATUS(endl);
 
  TFile out(outputFile.c_str(), "recreate");
 
  putObject(out, JMeta(argc, argv));
 
  out << ha;
  out << hb;
  out << hr;
  out << h2;
  out << g1;
 
  out.Write();
  out.Close();
 
  delete rpm;
}