JLegolas.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
#include <fstream>
#include <vector>
#include <cmath>
#include <limits>
 
#include "TROOT.h"
#include "TFile.h"
#include "TGraph.h"
#include "TGraphErrors.h"
#include "TGraphSmooth.h"
#include "TF1.h"
#include "TH1D.h"
#include "TSpline.h"
 
#include "JROOT/JMinimizer.hh"
 
#include "JTools/JRange.hh"
 
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
namespace {
 
  const char* const h0_1 = "h0";
}
 
/**
 * \file
 *
 * Auxiliary program to model transition time distribution generator from data (e.g.\ output of JPulsar.cc).
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
 
  typedef   JRange<double>  JRange_t;
 
  string    inputFile;
  string    outputFile;
  Double_t  bass;
  Double_t  span;
  string    pdf;
  string    cdf;
  JRange_t  T_ns;
  int       debug;
 
  try {
 
    JParser<> zap("Auxiliary program to model transition time distribution generator from data.");
    
    zap['f'] = make_field(inputFile,    "input file (output from JPulsar)");
    zap['o'] = make_field(outputFile)   = "";
    zap['B'] = make_field(bass,         "see TGraphSmooth")                                                = 0.00;
    zap['S'] = make_field(span,         "see TGraphSmooth")                                                = 0.01;
    zap['P'] = make_field(pdf,          "PDF output file; for inclusion in JPMTTransitTimeProbability.hh") = "";
    zap['C'] = make_field(cdf,          "CDF output file; for inclusion in JPMTTransitTimeGenerator.hh")   = "";
    zap['T'] = make_field(T_ns,         "time range for PDF and CDF")                                      = JRange_t(-20.0, +100.0);
    zap['d'] = make_field(debug)        = 1;
 
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }
 
 
  TH1D* h1 = NULL;
 
  TFile in(inputFile.c_str(), "exist");
 
  if (!in.IsOpen()) {
    FATAL("File " << inputFile << " not opened." << endl);
  }
 
  h1 = dynamic_cast<TH1D*>(in.Get(h0_1));
 
  if (h1 == NULL) {
    FATAL("No histogram <" << h0_1 << ">." << endl);
  }
 
 
  // Fit function
 
  TF1 f1("f1", "[0]*exp(-0.5*(x-[1])*(x-[1])/([2]*[2])) + [3]");
  
  // start values
    
  Double_t ymax  = 0.0;
  Double_t x0    = 0.0;
  Double_t sigma = 2.0;
  Double_t ymin  = 0.0;
  
  for (Int_t i = 1; i <= h1->GetNbinsX(); ++i) {
 
    const Double_t x = h1->GetBinCenter (i);
    const Double_t y = h1->GetBinContent(i);
 
    if (y > ymax) {
      ymax = y;
      x0   = x;
    }
  } 
  
  f1.SetParameter(0, ymax);
  f1.SetParameter(1, x0);
  f1.SetParameter(2, sigma);
  f1.SetParameter(3, ymin);  
  
  // fit
 
  h1->Fit(&f1, "LL", "same", x0 - 5 * sigma, x0 + 5 * sigma);
  
  // copy results
  
  ymax  = f1.GetParameter(0);
  x0    = f1.GetParameter(1);
  sigma = f1.GetParameter(2);
  ymin  = f1.GetParameter(3);
 
  
  // count pre- and delayed pulses and determine background
  
  Double_t yy = 0.0;
  Double_t yl = 0.0;
  Double_t yr = 0.0;
 
  Int_t    n0 = 0.0;  
  Double_t y0 = 0.0;
 
  for (Int_t i = 1; i <= h1->GetNbinsX(); ++i) {
 
    const Double_t x = h1->GetBinCenter (i);
    const Double_t y = h1->GetBinContent(i);
 
    if (T_ns(x - x0)) {
    
      yy += y;
    
      if      (x < x0 - 5.0 * sigma)
        yl += y;
      else if (x > x0 + 5.0 * sigma)
        yr += y;
 
    } else {
 
      n0 += 1;
      y0 += y;
    }
  }
 
  if (n0 != 0) {
    y0 /= n0;
  }
 
  // subtract average background
 
  yy -= y0;
  yl -= y0;
  yr -= y0;
    
  NOTICE("Number of pulses: " << yy      << endl);
  NOTICE("Pre-pulses:       " << yl / yy << endl);
  NOTICE("Delayed pulses:   " << yr / yy << endl);
 
 
  // extract data for smoothing
  
  typedef vector<Double_t>  JVector_t;
 
  JVector_t X;
  JVector_t Y;
  JVector_t EX;
  JVector_t EY;
 
  for (Int_t i = 1; i <= h1->GetNbinsX(); ++i) {
 
    const Double_t x = h1->GetBinCenter (i)  -  x0;   // zero at position of Gauss
    const Double_t y = h1->GetBinContent(i)  -  y0;   // subtract background
 
    if (T_ns(x)) {
      X .push_back(x);
      Y .push_back(y);
      EX.push_back(0.0);
      EY.push_back(sqrt(y));
    }
  }
  
  const int N = X.size();
  
  if (N <= 25) {
    FATAL("Not enough points." << endl);
  }
 
 
  TGraphErrors g0(N, X.data(), Y.data(), EX.data(), EY.data());
 
  g0.SetName("g0");
 
  TGraphErrors g1(g0);  // copy
 
  g1.SetName("g1");
 
 
  // divide TGraph data by TF1 fit function
 
  for (int i = 0; i != N; ++i) {
      
    const Double_t x = g1.GetX()[i];
    const Double_t f = f1.Eval(x + x0);
      
    g1.GetY() [i] /= f;
    g1.GetEY()[i] /= f;
  }
 
 
  // smooth TGraph
 
  TGraphSmooth gs("gs");
    
  TGraph* gout = gs.SmoothSuper(&g1, "", bass, span, false);
 
    
  // copy back results
    
  for (int i = 0; i != N; ++i) {
    g1.GetY()[i] = gout->GetY()[i];
  }
 
 
  // multiply TGraph data by TF1 fit function
 
  for (int i = 0; i != N; ++i) {
      
    const Double_t x = g1.GetX()[i];
    const Double_t f = f1.Eval(x + x0);
      
    g1.GetY() [i] *= f;
    g1.GetEY()[i] *= f;
  }
 
 
  if (pdf != "") {
 
    // normalise TGraph
 
    TGraph g2(g1);  // copy
 
    Double_t W = 0.0;
    
    for (int i = 0; i != N; ++i) {
      W += g2.GetY()[i];
    }
    
    for (int i = 0; i != N; ++i) {
      g2.GetY()[i] /= W;
    }
 
    ofstream out(pdf.c_str());
 
    out << "\t// produced by JLegolas.cc" << endl;
 
    const Double_t xmin = T_ns.getLowerLimit();
    const Double_t xmax = T_ns.getUpperLimit();
    const Double_t dx   = 0.25;
 
    for (Double_t x = xmin; x <= xmax + 0.5*dx; x += dx) {
      
      Double_t y = g2.Eval(x);
 
      if (y < 0.0) {
        y = 0.0;
      }
 
      out << "\t(*this)" 
          << "[" 
          << showpos   << FIXED(7,2) << x 
          << "] = " 
          << noshowpos << FIXED(8,6) << y
          << ";" << endl; 
    }
 
    out.close();
  }
 
  if (cdf != "") {
 
    // cumulate TGraph
    
    TGraph g2(g1);  // copy
 
    for (int i = 0; i+1 != N; ++i) {
      g2.GetY()[i+1] += g2.GetY()[i];
    }
 
    {    
      const Double_t xmin = g2.GetX()[ 0 ];
      const Double_t xmax = g2.GetX()[N-1];
      const Double_t dx   = (xmax - xmin) / N;
 
      const Double_t ymin = g2.GetY()[ 0 ];
      const Double_t ymax = g2.GetY()[N-1];
    
      for (int i = 0; i != N; ++i) {
      
        const Double_t x = g2.GetX()[i];
        const Double_t y = g2.GetY()[i];
      
        g2.GetX()[i] = (y - ymin) / ymax;
        g2.GetY()[i] =  x + 0.5 * dx;
      }
    }
    
    ofstream out(cdf.c_str());
 
    out << "\t// produced by JLegolas.cc" << endl;
    
    const Double_t xmin = 0.0;
    const Double_t xmax = 1.0;
    const Double_t dx   = 0.001;
 
    for (Double_t x = xmin; x <= xmax + 0.5*dx; x += dx) {
      
      out << "\t(*this)" 
          << "[" 
          << noshowpos << FIXED(6,4) << x 
          << "] = " 
          << showpos   << FIXED(9,5) << g2.Eval(x) 
          << ";" << endl; 
    }
 
    out.close();
  }
 
 
  if (outputFile != "") {
 
    // copy of original histogram data for overlaying with TGraph
    
    const Double_t xmin = X[ 0 ];
    const Double_t xmax = X[N-1];
    const Double_t dx   = (xmax - xmin) / N;
 
    TH1D h2("pmt", NULL, N, xmin - 0.5*dx, xmax + 0.5*dx);
    
    h2.Sumw2();
 
    // normalise
 
    Double_t W = 0.0;
 
    for (int i = 0; i != N; ++i) {
      W += Y[i];
    }
 
    W *= dx;
 
    for (int i = 0; i != N; ++i) {
 
      h2.SetBinContent(i+1, Y [i]/W);
      h2.SetBinError  (i+1, EY[i]/W);
 
      g1.GetY() [i] /= W;
      g1.GetEY()[i] /= W;
    }
    
    TFile out(outputFile.c_str(), "recreate");
    
    h1->Write();
    h2.Write();
    g0.Write();
    g1.Write();
    
    out.Write();
    out.Close();
  }
}