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