JFitToT.cc File Reference

#include <string>
#include <iostream>
#include <iomanip>
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TF1.h"
#include "TFitResult.h"
#include "TMath.h"
#include "JLang/JLangToolkit.hh"
#include "JTools/JConstants.hh"
#include "JDAQ/JDAQ.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JPMTParametersMap.hh"
#include "JROOT/JRootToolkit.hh"
#include "JTools/JRange.hh"
#include "JSupport/JMeta.hh"
#include "JCalibrate/JCalibrateToT.hh"
#include "JCalibrate/JFitToT.hh"
#include "Jeep/JPrint.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

Auxiliary program to fit time-over-threshold distributions. The fitted parameters are the PMT gain and gain spread.

Author

mkarel

Definition in file JFitToT.cc.

Function Documentation

main()

int main	(	int	argc,
		char **	argv
	)

Definition at line 48 of file JFitToT.cc.

{
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;
 
  typedef JRange<double>      JRange_t;
 
  string    inputFile;
  string    outputFile;
  string    detectorFile;
  string    pmtFile;
  bool      writeFits;
  JRange_t  X;
  string    option;
  double    WMin;
  string    regexp;
  int       debug;
 
  try { 
 
    JParser<> zap("Auxiliary program to fit time-over-threshold distributions.");
    
    zap['f'] = make_field(inputFile,         "input file (output from JCalibrateToT).");
    zap['o'] = make_field(outputFile,        "output file.")                                                   = "fit.root";
    zap['a'] = make_field(detectorFile,      "detector file.");
    zap['P'] = make_field(pmtFile,           "specify PMT file name that can be given to JTriggerEfficiency.") = "";
    zap['w'] = make_field(writeFits,         "write fit results.");
    zap['x'] = make_field(X,                 "ROOT fit range (time-over-threshold).")                          = JRange_t(10.0, 35.0);
    zap['O'] = make_field(option,            "ROOT fit option, see TH1::Fit.")                                 = "";
    zap['W'] = make_field(WMin,              "minimal histogram contents.")                                    = 1000.0;
    zap['R'] = make_field(regexp,            "regular expression for histogram name.")                         = MAKE_CSTRING(WILDCARD << _2SToT);
    zap['d'] = make_field(debug,             "debug.")                                                         = 1;
 
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }
 
 
  //----------------------------------------------------------
  // load detector file and PMT parameters
  //----------------------------------------------------------
 
  JDetector detector;
 
  try { 
    load(detectorFile, detector);
  }
  catch(const JException& error) {
    FATAL(error);
  }
 
 
  JPMTParametersMap parameters;
 
  if (pmtFile != "") {
    try { 
      parameters.load(pmtFile.c_str());
    }
    catch(const JException& error) {}
  }
 
  parameters.comment.add(JMeta(argc, argv));
 
     
  if (option.find('R') == string::npos) { option += 'R'; }
  if (option.find('S') == string::npos) { option += 'S'; }
  if (option.find('Q') == string::npos && debug < JEEP::debug_t) { option += 'Q'; }
 
 
  TFile* in = TFile::Open(inputFile.c_str(), "exist");
 
  if (in == NULL || !in->IsOpen()) {
    FATAL("File: " << inputFile << " not opened." << endl);
  }
 
 
  TH2D gain("gain", NULL,
            100,  0.0, 2.0,
            100,  0.0, 1.0);
  TH1D chi2("chi2", NULL, 100,  0.0, 5.0);
 
 
  TFile out(outputFile.c_str(), "recreate");
 
  //----------------------------------------------------------
  // loop over modules in detector file to fit histogram
  //----------------------------------------------------------
 
  for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {
 
    DEBUG("Module " << module->getID() << endl);
 
    //----------------------------------------------------------
    // get histogram for this module
    //----------------------------------------------------------
 
    TH2D* h2s = (TH2D*) in->Get(replace(regexp, WILDCARD, MAKE_STRING(module->getID())).c_str());
 
    if (h2s == NULL) {
 
      WARNING("No histogram " << module->getID() << "; skip fit." << endl);
 
      continue;
    }
 
    const int    ny   = h2s->GetYaxis()->GetNbins();
    const double ymin = h2s->GetYaxis()->GetXmin();
    const double ymax = h2s->GetYaxis()->GetXmax();
 
    TH1D chi2_1d      (MAKE_CSTRING(module->getID() << ".1chi2"),       NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS-0.5);
    TH1D gain_1d      (MAKE_CSTRING(module->getID() << ".1gain"),       NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS-0.5);
    TH1D gainspread_1d(MAKE_CSTRING(module->getID() << ".1gainspread"), NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS-0.5);
 
    for (int ix = 1; ix <= h2s->GetXaxis()->GetNbins(); ++ix) {
 
      const JPMTIdentifier id(module->getID(), ix-1);
 
      TH1D h1(MAKE_CSTRING(module->getID() << "." << id.getPMTAddress() << ".1ToT"), NULL, ny, ymin, ymax);
 
      h1.Sumw2();
 
      Double_t W = 0.0;
      
      for (int iy = 1; iy <= h2s->GetNbinsY(); ++iy) {
 
        const double x = h1.GetBinCenter(iy);
        const double z = h1.GetBinWidth (iy);
        const double y = h2s->GetBinContent(ix, iy);
 
        if (X(x)) {
          W += y;
        }
 
        h1.SetBinContent(iy, y/z);
        h1.SetBinError  (iy, sqrt(y)/z);
      }
 
      if (W <= WMin) {
 
        WARNING("Not enough entries for PMT " << id << ' ' << W << "; skip fit." << endl);
 
        continue;
      } 
 
      h1.Scale(1.0/h1.Integral());
 
 
      //----------------------------------------------------------
      // set-up fitting procedure
      //----------------------------------------------------------
      
      JFitToT fit(parameters.getPMTParameters(id),
                  max(ymin, X.first),
                  min(ymax, X.second));
 
 
      const TFitResultPtr result = fit(h1, option.c_str());
 
      chi2.Fill(TMath::Log10(result->Chi2() / result->Ndf()));
 
      if (debug >= JEEP::notice_t) {
        cout << "PMT " << id << " chi2 / NDF " << result->Chi2() << ' ' << result->Ndf() << ' ' << (result->IsValid() ? "" : "failed") << endl;
      }
 
      if (!result->IsValid()) {
 
        WARNING("Fit failure for PMT " << id << "; skip." << endl);
 
        continue;
      }
 
      if (ymin < X.getLowerLimit() || ymax > X.getUpperLimit()) {
 
        // add function with full range
        
        TF1* p = new TF1("f1&",
                         JFitToT::getValue,
                         ymin,
                         ymax,
                         JFitToT::getNumberOfModelParameters());
 
        p->SetParameters(fit.GetParameters());
        p->SetNpx(1000);
 
        h1.GetListOfFunctions()->Add(p);
      }
 
 
      // store fit results
 
      const JFitToTParameters values(fit.GetParameters());
      const JFitToTParameters errors(fit.GetParErrors());
 
      chi2_1d.      SetBinContent(ix, result->Chi2() / result->Ndf());
      gain_1d.      SetBinContent(ix, values.gain);
      gain_1d.      SetBinError  (ix, errors.gain);
      gainspread_1d.SetBinContent(ix, values.gainSpread);
      gainspread_1d.SetBinError  (ix, errors.gainSpread);
 
      gain.Fill(fit.gain, fit.gainSpread);
 
      if (fit.gain > 0.0) {
 
        parameters[id].gain       = fit.gain;
        parameters[id].gainSpread = fit.gainSpread;
 
      } else {
 
        WARNING("Error fit of PMT " << id << "; set default values." << endl);
 
        parameters[id].gain       = 0.0;
        parameters[id].gainSpread = parameters.getDefaultPMTParameters().gainSpread;
      }
        
      if (writeFits) {
        out << h1;
      }
    }
 
    if (writeFits) {
      out << chi2_1d << gain_1d << gainspread_1d;
    }
  }
    
  if (pmtFile != "") {
    parameters.store(pmtFile.c_str());
  }
 
  out << gain << chi2;
 
  out.Write();
  out.Close();
}