JTTS.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 "TString.h"
#include "TRegexp.h"
#include "TGraph.h"
#include "TGraphErrors.h"
#include "TGraphSmooth.h"
#include "TF1.h"
#include "TH1D.h"
#include "TSpline.h"

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


/**
 * \file
 *
 * Auxiliary program to model transition time distribution generator from data.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;

  string    inputFile;
  string    outputFile;
  Double_t  bass;
  Double_t  span;
  Double_t  xbin;
  int       numberOfBins;
  string    pdf;
  string    cdf;
  int       debug;

  try {

    JParser<> zap("Auxiliary program to model transition time distribution generator from data.");
    
    zap['f'] = make_field(inputFile);
    zap['o'] = make_field(outputFile)   = "";
    zap['B'] = make_field(bass)         = 0.0;
    zap['S'] = make_field(span)         = 0.0;
    zap['x'] = make_field(xbin)         = 1.0;              // [ns]
    zap['n'] = make_field(numberOfBins) = 1000;
    zap['P'] = make_field(pdf)          = "";
    zap['C'] = make_field(cdf)          = "";
    zap['d'] = make_field(debug)        = 1;

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


  const TRegexp regexp("V[0-9]+");

  TString buffer(inputFile.c_str());
  
  Ssiz_t len = 0;
  Ssiz_t pos = buffer.Index(regexp, &len);

  int version = 0;
  
  if (pos != -1) {
    buffer  = buffer(pos+1, len-1);
    version = buffer.Atoi();
  }
  
  NOTICE("File version " << version << endl);
  

  typedef vector<Double_t>  JVector_t;

  JVector_t X;
  JVector_t Y;
  JVector_t EX;
  JVector_t EY;

  ifstream in(inputFile.c_str());

  if (version == 0) {
    
    for (Double_t x = 0.0, y; in >> y; x += xbin) {
      X .push_back(x);
      Y .push_back(y);
      EX.push_back(0.0);
      EY.push_back(sqrt(y));
    }
    
  } else {
  
    in.ignore(numeric_limits<streamsize>::max(), '\n');
    
    for (Double_t x, y, dy, rms; in >> x >> y >> dy >> rms; ) {
      X .push_back(x);
      Y .push_back(y);
      EX.push_back(0.0);
      EY.push_back(dy);
    }
  }
  
  in.close();


  int N = X.size();
  
  if (N < 25) {
    FATAL("Not enough points." << endl);
  }

  if (version != 0) {

    xbin = (X[N-1] - X[0]) / (N - 1);
    
    NOTICE("Bin width [ns] " << xbin << endl);
  }

  // Histogram for log-likelihood fit

  TH1D h1("h1", NULL, N, X[0] - 0.5*xbin, X[N-1] + 0.5*xbin);
  
  h1.Sumw2();

  for (int i = 0; i != N; ++i) {
    h1.SetBinContent(i+1, Y [i]);
    h1.SetBinError  (i+1, EY[i]);
  }


  // 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 i = 0; i != N; ++i) {
    if (Y[i] > ymax) {
      ymax = Y[i];
      x0   = X[i];
    }
  } 
  
  f1.SetParameter(0, ymax);
  f1.SetParameter(1, x0);
  f1.SetParameter(2, sigma);
  f1.SetParameter(3, ymin);  
  

  h1.Fit(&f1, "LL", "same", x0 - 5 * sigma, x0 + 5 * sigma);
  
  
  // center data at position of Gauss
  
  ymax  = f1.GetParameter(0);
  x0    = f1.GetParameter(1);
  sigma = f1.GetParameter(2);
  ymin  = f1.GetParameter(3);

  if (version == 0) {

    for (JVector_t::iterator x = X.begin(); x != X.end(); ++x) {
      *x -= x0;
    }
    
    f1.SetParameter(1, x0 = 0);    
    
    if (ymin <= 0.0) {
      f1.SetParameter(3, 1e-4 * ymax);    // avoid division by zero
    }
  }
  
  // count pre- and delayed pulses
  
  Double_t yy = 0.0;
  Double_t yl = 0.0;
  Double_t yr = 0.0;
  
  for (int i = 0; i != N; ++i) {
    
    const Double_t x = X[i];
    const Double_t y = Y[i];      
    
    yy += y;
    
    if      (x < x0 - 5.0 * sigma)
      yl += y;
    else if (x > x0 + 5.0 * sigma)
      yr += y;
  }
    
  NOTICE("Number of pulses: " << yy      << endl);
  NOTICE("Pre-pulses:       " << yl / yy << endl);
  NOTICE("Delayed pulses:   " << yr / yy << endl);

  if (version == 0) {

    // skip leading and trailing zeros for smoothing operation
  
    int n1 =   0;
    int n2 = N - 1;
    
    while (n1 != N - 1 && Y[n1] == 0) ++n1;
    while (n2 !=   0   && Y[n2] == 0) --n2;
    
    //if (n1 !=   0  ) --n1;
    //if (n2 != N - 1) ++n2;
    
    if (n2 <= n1) {
      FATAL("No non-zero data." << endl);
    }
    
    X  = JVector_t(X .data() + n1, X .data() + n2);
    Y  = JVector_t(Y .data() + n1, Y .data() + n2);
    EX = JVector_t(EX.data() + n1, EX.data() + n2);
    EY = JVector_t(EY.data() + n1, EY.data() + n2);
    
    N  = X.size();
  }

  TGraphErrors g0(N, X.data(), Y.data(), EX.data(), EY.data());

  g0.SetName("g0");

  TGraphErrors g1(g0);  // copy

  g1.SetName("g1");


  if (version == 0) {

    // 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);
      
      g1.GetY() [i] /= f;
      g1.GetEY()[i] /= f;
    }

    // smooth TGraph
    
    TGraphSmooth   gs("gs");
    
    TGraph* gout = gs.SmoothSuper(&g1, "", bass, span, false, g1.GetEY());
    
    
    // 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);
      
      g1.GetY() [i] *= f;
      g1.GetEY()[i] *= f;
    }
  }

  if (pdf != "") {

    TGraph g2(g1);

    Double_t W = 0.0;
    
    for (int i = 0; i != N; ++i) {
      W += g2.GetY()[i];
    }
    
    W *= xbin;
    
    for (int i = 0; i != N; ++i) {
      g2.GetY()[i] /= W;
    }

    ofstream out(pdf.c_str());

    const Double_t xmin = g2.GetX()[ 0 ];
    const Double_t xmax = g2.GetX()[N-1];

    for (Double_t x = xmin, dx = (xmax - xmin) / numberOfBins; x <= xmax + 0.5*dx; x += dx) {
      
      out << "\t(*this)" 
          << "[" 
          << showpos   << FIXED(7,2) << x 
          << "] = " 
          << noshowpos << FIXED(6,4) << g2.Eval(x) 
          << ";" << endl; 
    }

    out.close();
  }

  if (cdf != "") {
    
    TGraph g2(g1);

    for (int i = 0, j = 1; j != N; ++i, ++j) {
      g2.GetY()[j] += g2.GetY()[i];
    }
    
    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 * xbin;
    }
    
    ofstream out(cdf.c_str());
    
    const Double_t xmin = 0.0;
    const Double_t xmax = 1.0;

    for (Double_t x = xmin, dx = (xmax - xmin) / numberOfBins; 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 histogram centered at peak position

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

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