JFit_HTR.cc File Reference

Program to fit hit time residual distributions according to a MC model. More...

#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"

Go to the source code of this file.

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

Detailed Description

Program to fit hit time residual distributions according to a MC model.

If the option -A is set, the time calibration will be stored in the detector file.

Author

mkarel

Definition in file JFit_HTR.cc.

Function Documentation

int main	(	int	argc,
		char **	argv
	)

Definition at line 120 of file JFit_HTR.cc.

{
  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() ;

}