JFitToT.cc

Go to the documentation of this file.

#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 "km3net-dataformat/online/JDAQ.hh"
#include "km3net-dataformat/online/JDAQHeader.hh"

#include "JLang/JLangToolkit.hh"
#include "JPhysics/JConstants.hh"

#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JPMTParametersMap.hh"
#include "JDetector/JPMTAnalogueSignalProcessor.hh"

#include "JROOT/JRootToolkit.hh"

#include "JTools/JRange.hh"

#include "JSupport/JMeta.hh"
#include "JSupport/JSupport.hh"
#include "JSupport/JFileRecorder.hh"
#include "JSupport/JSingleFileScanner.hh"

#include "JCalibrate/JCalibrateToT.hh"
#include "JCalibrate/JFitToT.hh"

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

namespace {

  /**
   * Wild card for histogram naming.\n
   * The wild card will replaced by the module identifier.
   */
  static const char* const WILDCARD = "%";
}


/**
 * \file
 *
 * Auxiliary program to fit time-over-threshold distributions.
 * The fitted parameters are the PMT gain and gain spread.
 * \author mkarel
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;

  typedef JRange<double>                                       JRange_t;
  typedef JTYPELIST<JMeta, JDAQHeader, JRootTypes_t>::typelist typelist;

  string                          inputFile;
  JFileRecorder<typelist>         outputFile;
  string                          detectorFile;
  string                          pmtFile;
  bool                            writeFits;
  
  double                          Wmin              =   5e3;                // Minimal histogram weight in fit range
  double                          LeftMargin        =  10.0;                // Fit range left of time-over-threshold peak 
  double                          RightMargin       =  10.0;                // Fit range right of time-over-threshold peak
  double                          gradientThreshold =  -0.005;              // Minimal right-to-left gradient value in time-over-threshold peak-search,\n
                                                                            // normalized with respect to the total number of histogram entries
  JRange_t                        gainLimits      (FITTOT_GAIN_MIN,
                                                   FITTOT_GAIN_MAX);        // gain       limits
  JRange_t                        gainSpreadLimits(FITTOT_GAINSPREAD_MIN,
                                                   FITTOT_GAINSPREAD_MAX);  // gainspread limits
  
  JRange_t                        ToTfitRange;
  string                          option;
  
  string                          regexp;
  int                             debug;

  try { 

    JProperties properties;
    
    properties.insert(gmake_property(Wmin));
    properties.insert(gmake_property(LeftMargin));
    properties.insert(gmake_property(RightMargin));
    properties.insert(gmake_property(gradientThreshold));
    properties.insert(gmake_property(gainLimits));
    properties.insert(gmake_property(gainSpreadLimits));
    
    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['O'] = make_field(option,            "ROOT fit options, see TH1::Fit")                                 = "";
    zap['@'] = make_field(properties,        "fit properties")                                                 = JPARSER::initialised();
    zap['x'] = make_field(ToTfitRange,       "ROOT fit range (time-over-threshold).")                          = JRange_t(0.0, 35.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 parametersMap;

  if (!pmtFile.empty()) {

    try {
      parametersMap.load(pmtFile.c_str()); 
    }
    catch(const JException& error) {}
  }

  // Set ROOT fit options
  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);

  outputFile.open();

  if (!outputFile.is_open()) {
    FATAL("Error opening file " << outputFile << endl);
  }
  
  outputFile.put(JMeta(argc, argv));

  //----------------------------------------------------------
  // loop over modules in detector file to fit histogram
  //----------------------------------------------------------

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

    NOTICE("Module " << setw(10) << module->getID() << ' ' << module->getLocation() << endl);

    if (module->empty()) {
      continue;
    }

    //----------------------------------------------------------
    // 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 double 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);

      const string name = MAKE_STRING(module->getID()    << '.' << FILL(2,'0') <<
                                      id.getPMTAddress() << FITTOT_SUFFIX);

      DEBUG("Processing histogram " << name << endl);

      TH1D h1(name.c_str(), NULL, ny, ymin, ymax); 
 
      h1.Sumw2();
      
      for (int iy = h2s->GetNbinsY(); iy >= 1; --iy) {
 
        const double w = h1.GetBinWidth (iy);
        const double y = h2s->GetBinContent(ix, iy);
  
        h1.SetBinContent(iy, y/w);
        h1.SetBinError  (iy, sqrt(y)/w);
      }

      JRange_t     range(ToTfitRange);
      
      const double weight = h1.Integral(h1.FindBin(range.getLowerLimit()),
                                        h1.FindBin(range.getUpperLimit()));
      
      if (weight <= Wmin) {

        WARNING("Insufficient histogram entries for PMT "  << id <<
                "(" << weight << " < " << Wmin << "); skip fit." << endl);
        
        if (writeFits) {
          outputFile.put(h1);
        }
        
        continue;
      }

      //----------------------------------------------------------
      // find the mode of the time-over-threshold distribution
      //----------------------------------------------------------

      double max     = 0.0;
      double ToTmode = JDETECTOR::TIME_OVER_THRESHOLD_NS;

      JPMTParameters& parameters = parametersMap[id];

      DEBUG(RIGHT(10) << "ToT"       << RIGHT(10) << "Counts" << RIGHT(10) << "gradient" <<
            RIGHT(10) << "threshold" << RIGHT(10) << "Mode"   << RIGHT(10) << "Max" << endl);
      
      for (int i = h1.GetBinCenter(ny-1);
               i > h1.GetBinCenter(h1.FindBin(parameters.mean_ns + parameters.sigma_ns));
             --i) {
 
         const double x = h1.GetBinCenter (i);
         const double y = h1.GetBinContent(i);
 
         const double gradient = ( (h1.GetBinContent(i-1) - h1.GetBinContent(i+1)) /
                                   (h1.GetBinCenter (i+1) - h1.GetBinCenter (i-1)) );

         DEBUG(FIXED(10,1) << x                                        << FIXED(10,1) << y       << FIXED(10,1) << gradient <<
               FIXED(10,1) << -fabs(gradientThreshold) * h1.Integral() << FIXED(10,1) << ToTmode << FIXED(10,1) << max << endl);


         if (y > max) {
 
           ToTmode = x;
           max     = y;
           
         } else if (gradient < -fabs(gradientThreshold) * h1.Integral()) {

           break;
         }
      }
      
      //----------------------------------------------------------
      // set fit range
      //----------------------------------------------------------

      if (!range.is_valid()) {
        range.setRange(ToTmode - LeftMargin,
                       ToTmode + RightMargin);
      }
      
      range.setRange(std::max(h1.GetXaxis()->GetXmin(), range.getLowerLimit()),
                     std::min(h1.GetXaxis()->GetXmax(), range.getUpperLimit()));

      DEBUG(LEFT(10) << "Fit-range: " << FIXED(20,1) << range << endl);

      //----------------------------------------------------------
      // perform minimum chi-square fit
      //----------------------------------------------------------

      JFitToT fit(parameters, range);
      
      setParLimits(fit, fit.getModelParameter(&JFitToT::JFitToTParameters::gain),
                   gainLimits.getLowerLimit(), gainLimits.getUpperLimit());
      setParLimits(fit, fit.getModelParameter(&JFitToT::JFitToTParameters::gainSpread),
                   gainSpreadLimits.getLowerLimit(), gainSpreadLimits.getUpperLimit());
      
      const double initGain       = fit.getCPU().getNPE(ToTmode);
      const double initGainSpread = initGain / 3.0;
      
      if (!gainLimits(initGain)) { // Invalid initial gain value 
        
        WARNING("Invalid initial gain for PMT " << id << "; set default values." << endl);
        
        parameters.gain       = gainLimits      .constrain(initGain);
        parameters.gainSpread = gainSpreadLimits.constrain(initGainSpread);

        if (writeFits) {
          outputFile.put(h1);
        }

        continue;
      }
      
      fit.gain       = initGain;
      fit.gainSpread = initGainSpread;
      
      const TFitResultPtr result = fit(h1, option);

      if (result.Get() == NULL) {
        FATAL("Invalid TFitResultPtr " << h1.GetName() << endl);
      }

      if (int(result) < 0 || !result->IsValid()) { // Fit failed

        WARNING("Fit failure for PMT " << id << "; set initialization values." << endl);
          
        h1.GetListOfFunctions()->Clear();

        parameters.gain       = initGain;
        parameters.gainSpread = initGainSpread;

        if (writeFits) {
          outputFile.put(h1);
        }

        continue;
      }

      if ((ymin < range.getLowerLimit()  || ymax > range.getUpperLimit()) &&
          (option.find("0") == string::npos && option.find("N") == string::npos)) {
        
        // add function with full range
        
        TF1* f1 = new TF1(MAKE_CSTRING(FITTOT_FNAME << "_fullRange"),
                          &fit,
                          &JFitToT::getValue,
                          ymin,
                          round(parameters.saturation) - 0.5,
                          JFitToT::getNumberOfModelParameters());
        f1->SetParameters(fit.GetParameters());
        f1->SetLineStyle(kDotted);
        f1->SetNpx(1000);

        h1.GetListOfFunctions()->Add(f1);
      }

      // store fit results

      const double reducedChi2 = result->Chi2() / result->Ndf();

      const JFitToTParameters values(fit.GetParameters());
      const JFitToTParameters errors(fit.GetParErrors());
      
      chi2_1d.      SetBinContent(ix, reducedChi2);
      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);
      chi2.Fill(TMath::Log10(reducedChi2));
      
      NOTICE("PMT " << id << " chi2 / NDF " << result->Chi2() << ' ' << result->Ndf() << endl);
      
      parameters.gain       = fit.gain;
      parameters.gainSpread = fit.gainSpread;
        
      if (writeFits) {
        outputFile.put(h1);
      }
    }

    if (writeFits) {
      outputFile.put(chi2_1d);
      outputFile.put(gain_1d);
      outputFile.put(gainspread_1d);
    }
  }

  {
    vector<JMeta> buffer(1, JMeta(argc, argv));

    for (JSingleFileScanner<JMeta> in(inputFile); in.hasNext(); ) {

      const JMeta* meta = in.next();

      buffer.push_back(*meta);

      outputFile.put(*meta);
    }

    for (vector<JMeta>::const_reverse_iterator i = buffer.rbegin(); i != buffer.rend(); ++i) {
      parametersMap.comment.add(*i);
    }
  }

  for (JRootFileReader<JDAQHeader> in(inputFile.c_str()); in.hasNext(); ) {
    outputFile.put(*in.next());
  }

  
  if (pmtFile != "") {
    parametersMap.store(pmtFile.c_str());
  }

  outputFile.put(gain);
  outputFile.put(chi2);

  outputFile.close();
}