JFrodo.cc

Go to the documentation of this file.

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

#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TF1.h"
#include "TMath.h"
#include "TFitResult.h"

#include "JTools/JRange.hh"
#include "JDetector/JCalibration.hh"
#include "JDetector/JPMTAnalogueSignalProcessor.hh"

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

namespace {

  /**
   * Name of histogram with fit results from JCalibrateMuon.
   */
  const char* h2_t           =  "ha";

  const char* gauss_t        =  "Gauss";           //!< Gaussian
  const char* emg_t          =  "EMG";             //!< Exponentially Modified Gaussian


  using namespace JPP;

  JPMTAnalogueSignalProcessor cpu;

  /**
   * Fit function.
   *
   * \param  x            abscissa  values
   * \param  parameters   parameter values
   * \return              ordinate  value
   */
  Double_t getRisetime(Double_t* x, Double_t* parameters)
  {
    const double t0  = parameters[0];
    const double p1  = parameters[1];
    const double t1  = parameters[2];

    const double npe = cpu.getNPE(x[0]);

    double ts = t0  +  cpu.getRiseTime(npe)  +  t1 * (1.0 - exp(-npe*p1));

    return ts;
  }
}


/**
 * \file
 * Fit time-slewing histogram in output of JCalibrateMuon.cc in calibration mode.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;

  typedef JRange<double>  JRange_t;

  vector<string>  inputFile;
  string          outputFile;
  string          function;
  JRange_t        T_ns;
  bool            writeFits;
  JPMTParameters  parameters;
  JRange_t        X;
  string          option;
  int             debug;

  try {

    JProperties properties = parameters.getProperties();

    JParser<> zap("Program to fit time-slewing histogram in output of JCalibrateMuon.");

    zap['f'] = make_field(inputFile,         "input files (output from JCalibrateMuon).");
    zap['o'] = make_field(outputFile,        "output file.")                                    = "frodo.root";
    zap['F'] = make_field(function,          "fit function")                                    = gauss_t, emg_t;
    zap['T'] = make_field(T_ns,              "ROOT fit range around maximum [ns].")             = JRange_t(-7.5,+7.5);
    zap['w'] = make_field(writeFits,         "write fit results.");
    zap['P'] = make_field(properties)       = JPARSER::initialised();
    zap['x'] = make_field(X,                 "ROOT fit range for time-over threshold")          = JRange_t(0.0, 256.0);
    zap['O'] = make_field(option,            "ROOT fit option, see TH1::Fit.")                  = "";
    zap['d'] = make_field(debug)            = 1;

    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }


  if (!T_ns.is_valid()) {
    FATAL("Invalid fit range [ns] " << T_ns << endl);
  }

  cpu.setPMTParameters(parameters);

  if (option.find('S') == string::npos) { option += 'S'; }
  if (option.find('Q') == string::npos && debug < JEEP::debug_t) { option += 'Q'; }

  // rebin time-over-threshold

  const Double_t x[] = {
    -0.5,  3.5,  6.5,  9.5, 12.5, 15.5, 18.5, 
    20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5,
    30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5,
    40.5, 42.5, 44.5, 46.5, 48.5,
    50.5, 52.5, 54.5, 56.5, 58.5,
    60.5, 65.5,
    70.5, 75.5,
    80.5, 85.5,
    90.5, 95.5,
    100.5, 120.5, 140.5, 160.5, 180.5, 200.5, 250.5
  };

  const int N = sizeof(x)/sizeof(x[0]);

  TH2D* h2 = NULL;

  for (vector<string>::const_iterator i = inputFile.begin(); i != inputFile.end(); ++i) {

    NOTICE("Processing " << *i << endl);

    TFile in(i->c_str(), "exist");

    if (!in.IsOpen()) {
      FATAL("File " << *i << " not opened." << endl);
    }

    TH2D* p = dynamic_cast<TH2D*>(in.Get(h2_t));
      
    if (p == NULL) {
      FATAL("File " << *i << " has no histogram <" << h2_t << ">." << endl);
    }

    if (h2 == NULL)
      h2 = (TH2D*) p->Clone();
    else
      h2->Add(p);

    h2->SetDirectory(0);

    in.Close();
  }

  if (h2 == NULL) {
    FATAL("No histogram <" << h2_t << ">." << endl);
  }

  // histogram fit results

  TH1D h0("h0", NULL, N-1, x);

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

  for (int i = 1; i <= h0.GetXaxis()->GetNbins(); ++i) {

    const Int_t imin = h2->GetXaxis()->FindBin(h0.GetXaxis()->GetBinLowEdge(i));
    const Int_t imax = h2->GetXaxis()->FindBin(h0.GetXaxis()->GetBinUpEdge (i));

    TH1D h1(MAKE_CSTRING(i << ".1D"), NULL, h2->GetYaxis()->GetNbins(), h2->GetYaxis()->GetXmin(), h2->GetYaxis()->GetXmax());

    // start values
    
    Double_t ymax   =  0.0;
    Double_t x0     =  0.0;       // [ns]
    Double_t sigma  =  3.5;       // [ns]
    
    for (int iy = 1; iy <= h2->GetYaxis()->GetNbins(); ++iy) {

      Double_t x = h1.GetBinCenter (iy);
      Double_t y = 0.0;

      for (int ix = imin; ix != imax; ++ix) {
        y += h2->GetBinContent(ix,iy);
      }

      h1.SetBinContent(iy, y);

      if (y > ymax) {
        ymax = y;
        x0   = x;
      }
    }


    TF1* f1 = NULL;

    if        (function == gauss_t) {

      f1 = new TF1(function.c_str(), "[0]*TMath::Gaus(x,[1],[2]) + [3]");

      f1->SetParameter(0, ymax);
      f1->SetParameter(1, x0);
      f1->SetParameter(2, sigma);
      f1->SetParameter(3, 0.0);

      f1->SetParLimits(3, 0.0, ymax);
    
    } else if (function == emg_t) {

      f1 = new TF1(function.c_str(), "[0]*TMath::Exp(0.5*[3]*(2.0*[1]+[3]*[2]*[2]-2.0*x))*TMath::Erfc(([1]+[3]*[2]*[2]-x)/(TMath::Sqrt(2.0)*[2])) +[4]");

      f1->SetParameter(0, ymax);
      f1->SetParameter(1, x0 - 0.5*sigma);
      f1->SetParameter(2, sigma);
      f1->SetParameter(3, 0.04);
      f1->SetParameter(4, 0.0);

      f1->SetParLimits(4, 0.0, ymax);

    } else {

      FATAL("Unknown fit function " << function << endl);
    }

    // fit

    Double_t xmin = x0 + T_ns.getLowerLimit();
    Double_t xmax = x0 + T_ns.getUpperLimit();

    if (xmin < h1.GetXaxis()->GetXmin()) { xmin = h1.GetXaxis()->GetXmin(); }
    if (xmax > h1.GetXaxis()->GetXmax()) { xmax = h1.GetXaxis()->GetXmax(); }
      
    TFitResultPtr result = h1.Fit(f1, option.c_str(), "same", xmin, xmax);

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

    if (debug >= notice_t || !result->IsValid()) {
      cout << "Bin: "
           << setw(4)    << i                     << ' '
           << FIXED(7,3) << f1->GetParameter(1)   << " +/- "
           << FIXED(7,3) << f1->GetParError (1)   << ' '
           << FIXED(7,3) << result->Chi2()        << '/'
           << result->Ndf()                       << ' '
           << (result->IsValid() ? "" : "failed") << endl;
    }

    if (result->IsValid()) {
      h0.SetBinContent(i, f1->GetParameter(1));
      h0.SetBinError  (i, f1->GetParError (1));
    }

    if (writeFits) {
      h1.Write();
    }

    delete f1;
  }

  // Fit

  TF1 f1("f1", getRisetime, x[0], x[N-1], 3);

  f1.SetParameter(0,   0.0);
  f1.SetParameter(1,   2.0e-2);
  f1.FixParameter(2,  -7.0);

  TFitResultPtr result = h0.Fit(&f1, "S", "same", X.constrain(x[0]), X.constrain(x[N-1]));

  if (debug >= notice_t || !result->IsValid()) {
    cout << "Time-slewing correction"           << endl;
    cout << FIXED(7,3) << result->Chi2()        << '/'
         << result->Ndf()                       << ' '
         << (result->IsValid() ? "" : "failed") << endl;

    for (int i = 0; i != f1.GetNpar(); ++i) {
      cout << "parameter[" << i << "] = " << FIXED(8,4) << f1.GetParameter(i) << " +/- " << FIXED(8,4) << f1.GetParError (i) << endl;
    }
  }

  cout << "// Produced by JFrodo.cc; to be included on JGetRiseTime.hh." << endl;
 
  const double t0 = f1.Eval(TIME_OVER_THRESHOLD_NS);

  for (int i = 0; i != 256; ++i) {
    cout << "this->push_back(" << FIXED(6,2) << f1.Eval((Double_t) i) - t0 << ");" << endl;
  }
 
  h0.Write();

  out.Close();
}