JFit_HTR.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <algorithm>

#include "evt/Head.hh"
#include "evt/Evt.hh"
#include "JDAQ/JDAQEvent.hh"
#include "JDAQ/JDAQTimeslice.hh"
#include "JDAQ/JDAQSummaryslice.hh"

#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JModuleRouter.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/JTreeScanner.hh"
#include "JSupport/JFileRecorder.hh"
#include "JSupport/JSupport.hh"

#include "JFit/JLine3Z.hh"
#include "JFit/JGandalf.hh"
#include "JFit/JFitToolkit.hh"
#include "JFit/JEvt.hh"
#include "JFit/JEvtToolkit.hh"
#include "JFit/JRegressor.hh"

#include "JTools/JConstants.hh"
#include "JTools/JFunction1D_t.hh"
#include "JTools/JFunctionalMap_t.hh"
#include "JTools/JRange.hh"

#include "JPhysics/JPDFTable.hh"
#include "JPhysics/JPDFToolkit.hh"

#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

#include "TKey.h"
#include "TH1D.h"
#include "TF1.h"

namespace {

  using namespace JPP;
  using namespace std ;

  bool selectHTRentry(const JFit& bf) {

    //double Q = bf.getQ() ;
    double dz = bf.getDZ() ;
    //int N = bf.getN() ;
    //double NDF = bf.getNDF() ;
    //double chi2 = 4.0*(Q-NDF)*NDF ; // inverse of JFIT::getQuality

    // selection criteria
    if (dz > -0.97) { return false ; }

    return true ;

  }

  double logPoison(double n, double nhat, double logP_min = -999999.9) {
    if (nhat < 0.0 || n < 0.0) { FATAL("logPoisson: n (" << n <<") or nhat (" << nhat << ") is < 0.0" << std::endl) ; }
    if (n == 0.0 && nhat == 0.0) { return 0.0 ; }
    if (n == 0.0 && nhat > 0.0) { return logP_min ; }
    if (n >= 0.0 && nhat == 0.0) { return logP_min ; }
    return TMath::Log(TMath::Poisson(n, nhat)) ;
  }

  double getLogQuality(const TH1D* data, const TH1D* model, int di, double n_min = 0.0001, double nhat_min = 0.0001) {

    double q = 0.0 ;
    const double logQ_min = -999999.9 ;
    const double logQ_empty = TMath::Max(logQ_min, logPoison(n_min, nhat_min)) ;
    int Nempty = 0 ;

    for (int i = 1 ; i <= data->GetNbinsX() ; ++i) {

      if (i+di < 1  || i+di > model->GetNbinsX()) {

        q += logQ_empty ;
        Nempty++ ;

        continue ;

      }

      double n    = TMath::Max(n_min,    data->GetBinContent(i)) ;
      double nhat = TMath::Max(nhat_min, model->GetBinContent(i+di)) ;

      q += TMath::Max(logQ_min, logPoison(n, nhat)) ;

    }

    return q ;

  }

}

/**
 * \file
 *
 * Program to fit hit time residual distributions according to a MC model.
 *
 * If the option <tt>-A</tt> is set, the time calibration will be stored in the detector file.
 * \author mkarel
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;

  string         detectorFile ;  
  string         inputFile_data;
  string         inputFile_mc;
  double         livetime_data;
  double         livetime_mc;
  string         outputFile;
  JRange<double> dtrange ;
  bool           noInterString ;
  bool           overwriteDetector;
  int            debug;

  try { 

    JParser<> zap;
    
    zap['f'] = make_field(inputFile_data);
    zap['m'] = make_field(inputFile_mc);
    zap['F'] = make_field(livetime_data);
    zap['M'] = make_field(livetime_mc); 
    zap['o'] = make_field(outputFile)          = "HTR_plots.root";
    zap['a'] = make_field(detectorFile);
    zap['A'] = make_field(overwriteDetector);
    zap['T'] = make_field(dtrange)             = JRange<double>(-5.0, 5.0);  // fit range for time change fit
    zap['I'] = make_field(noInterString) ;
    zap['d'] = make_field(debug)               = 1;
    
    zap(argc, argv);
  }
  catch(const exception& error) {
    FATAL(error.what() << endl);
  }


  JDetector detector;

  try {
    load(detectorFile, detector);
  }
  catch(const JException& error) {
    FATAL(error);
  }

  double dt_min = -100.0 ;
  double dt_max = 100.0 ;
  double dt_N =+(dt_max-dt_min)+1 ;

  cout.tie(&cerr);

  TFile infile_data(inputFile_data.c_str(),"read") ;
  TFile infile_mc(inputFile_mc.c_str(),"read") ;

  TFile outfile(outputFile.c_str(),"recreate") ;

  vector<double> timeshift(detector.size(), 0.0) ;

  for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {

    int domid = module->getID() ;

    TString DOMID(Form("%d", domid)) ;
    TTree* htr_tree_data = (TTree*)infile_data.Get(DOMID+".htr") ;
    if (htr_tree_data == NULL) { WARNING("Couldn't find model " << DOMID << ".htr in data file " << inputFile_data << " skipping DOM " << endl) ; continue ; }

    NOTICE( "processing " << DOMID << endl ) ;

    double htr_dt ;
    JFit* bestfit = new JFit() ;
    htr_tree_data->SetBranchAddress("dt",&htr_dt);
    htr_tree_data->SetBranchAddress("bestfit",&bestfit) ;

    TH1D htr_data(DOMID+".htr_data", DOMID+".htr_data", dt_N, dt_min-0.5, dt_max+0.5) ;

    Long64_t nentries = htr_tree_data->GetEntries();
    for (Long64_t i=0;i<nentries;i++) {
      htr_tree_data->GetEntry(i);

      // selection criteria
      if (selectHTRentry(*bestfit)) { htr_data.Fill(htr_dt) ; }

    }

    if (htr_data.Integral() < 1) { WARNING("No data in data file " << inputFile_data << " skipping DOM " << DOMID << endl) ; continue ; }

    // find corresponding tree in model
    TTree* htr_tree_mc = (TTree*)infile_mc.Get(htr_tree_data->GetName()) ;
    if (htr_tree_mc == NULL) { WARNING("Couldn't find model " << htr_tree_data->GetName() << " in model file " << inputFile_mc << " skipping DOM " << endl) ; continue ; }
    htr_tree_mc->SetBranchAddress("dt",&htr_dt);
    htr_tree_mc->SetBranchAddress("bestfit",&bestfit) ;

    TH1D htr_mc(DOMID+".htr_mc", DOMID+".htr_mc", dt_N, dt_min-0.5, dt_max+0.5) ;

    nentries = htr_tree_mc->GetEntries();
    for (Long64_t i=0;i<nentries;i++) {
      htr_tree_mc->GetEntry(i);

      // selection criteria
      if (selectHTRentry(*bestfit)) { htr_mc.Fill(htr_dt) ; }

    }

    if (htr_mc.Integral() < 1) { WARNING("No data in model file " << inputFile_mc << " skipping DOM " << DOMID << endl) ; continue ; } 

    htr_mc.Sumw2() ;
    htr_data.Sumw2() ;
    htr_mc.Scale(1.0/livetime_mc) ;
    htr_data.Scale(1.0/livetime_data) ;

    const double nmin = 0.01/livetime_mc ; // TODO: relatively arbitrary
    TH1D q1(DOMID+".1q", DOMID+".q1", 41, -20.5, 20.5) ;

    for (int di = 0 ; di < q1.GetNbinsX() ; di++) {
      q1.SetBinContent(di+1, getLogQuality(&htr_data, &htr_mc, q1.GetBinCenter(di+1), nmin, nmin)) ;
    }

    double dt_bestfit = q1.GetXaxis()->GetBinCenter(q1.GetMaximumBin()) ;

    // fit peak with parabola for time changes more precise than 1ns.
    TF1 f1(DOMID+".1f", "([0]*x+[1])*x+[2]", dt_bestfit+dtrange.getLowerLimit() , dt_bestfit+dtrange.getUpperLimit()) ;
    f1.SetParameter(0, 1.0) ;
    f1.SetParameter(1, 0.0) ;
    f1.SetParameter(2, 10.0) ;
    q1.Fit(&f1, "QR") ;

    if (fabs(dt_bestfit - -0.5*f1.GetParameter(1)/f1.GetParameter(0)) < 1.0) {
      dt_bestfit = -0.5*f1.GetParameter(1)/f1.GetParameter(0) ;
    }

    timeshift[distance(detector.begin(), module)] = dt_bestfit ;

    outfile.cd() ;
    TDirectory* dir = outfile.mkdir(DOMID);
    dir->cd() ;
    htr_data.Write() ;
    htr_mc.Write() ;
    q1.Write() ;

  }

  for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {

    // average shift to zero
    double ave = 0.0 ;
    int N = 0 ;
    for (JDetector::iterator mod = detector.begin(); mod != detector.end(); ++mod) {

      if (noInterString && module->getString() != mod->getString()) { continue ; }

      ave += timeshift.at(distance(detector.begin(), mod)) ;
      N++ ;

    }
    ave /= N ;

    DEBUG( "Best fit shift " << module->getID() << " : " << timeshift[distance(detector.begin(), module)]-ave << " [ns] " << endl ) ;

    if (overwriteDetector) {
      for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
        module->getPMT(pmt).addT0(timeshift[distance(detector.begin(), module)]-ave);
      }
    }

  }

  if (overwriteDetector) {

    NOTICE("Store calibration data on file " << detectorFile << endl);

    store(detectorFile, detector);

  }

  outfile.Write() ;

}