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 "JDynamics/JDynamics.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/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 "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

int main	(	int	argc,
		char **	argv
	)

Definition at line 66 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 JTYPELIST<JCOMPASS::JOrientation,
                    JACOUSTICS::JEvt>::typelist                    calibration_types;
  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;
  histogram_type           calibrate;
  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['P'] = make_field(pdfFile);
    zap['c'] = make_field(calibrate, "Histogram for time calibration per optical module.");
    zap['@'] = make_field(parameters)          = JPARSER::initialised();
    zap['d'] = make_field(debug)               = 1;
    
    zap(argc, argv);
  }
  catch(const exception& error) {
    FATAL(error.what() << endl);
  }

  cout.tie(&cerr);

  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);
  }


  getMechanics.load(detector.getID());

  unique_ptr<JDynamics> dynamics;

  try {

    dynamics.reset(new JDynamics(detector, Tmax_s));

    dynamics->load(calibrationFile);
  }
  catch(const exception& error) {
    if (!calibrationFile.empty()) {
      FATAL(error.what());
    }
  }

  const JModuleRouter router(dynamics ? dynamics->getDetector() : detector);

  JSummaryFileRouter  summary(inputFile, parameters.R_Hz);

  JMuonGandalf fit(parameters, router, summary, pdfFile, 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, detector.size(), -0.5, detector.size() - 0.5,
                calibrate.getNumberOfBins(), calibrate.getLowerLimit(), calibrate.getUpperLimit());

  JManager_t g1(new TProfile("%", NULL, detector.size(), -0.5, detector.size() - 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) {
            
        double rate_Hz = summary.getRate(*i);

        if (rate_Hz <= 0.0) {
          rate_Hz = summary.getRate();
        }

        JHitW0 hit(*i, rate_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(hit->getModuleID());
      }

      const JLine3Z ta(tz);

      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

          for (JDataW0_t::const_iterator hit = q; hit != __end; ++hit) {

            const int    index    = router.getIndex(*id);
            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();
}