JDemoPDF.cc

Go to the documentation of this file.

 
#include <string>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <assert.h>
 
#include "TROOT.h"
#include "TFile.h"
#include "TKey.h"
#include "TH1D.h"
#include "TProfile.h"
 
#include "JTools/JFunction1D_t.hh"
#include "JTools/JFunctionalMap_t.hh"
#include "JTools/JQuantiles.hh"
#include "JPhysics/JPDFTable.hh"
#include "JPhysics/JPDFTypes.hh"
#include "JGeometry3D/JAngle3D.hh"
#include "Jeep/JTimer.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
 
/**
 * \file
 * Demonstration program to plot RMS of arrival time of first hit as a function of
 * the minimal distance of approach of a muon to the PMT.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
 
  string      inputFile;
  string      outputFile;  
  double      E_GeV;
  double      R_Hz;
  JAngle3D    dir;
  int         debug;
 
  try { 
 
    JParser<> zap("Demonstration program to plot RMS of arrival time of first hit as a function of the minimal distance of approach of a muon to the PMT.");
 
    zap['f'] = make_field(inputFile);
    zap['o'] = make_field(outputFile)       = "demo.root";
    zap['E'] = make_field(E_GeV,     "muon energy [GeV]");
    zap['R'] = make_field(R_Hz,      "background rate [Hz]");
    zap['D'] = make_field(dir,       "(theta, phi) of PMT [rad]");
    zap['d'] = make_field(debug)            = 0;
 
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }
 
 
  typedef JSplineFunction1S_t                                     JFunction1D_t;
  typedef JMAPLIST<JPolint1FunctionalMap,
                   JPolint1FunctionalGridMap,
                   JPolint1FunctionalGridMap>::maplist            JMapList_t;
  typedef JPDFTable<JFunction1D_t, JMapList_t>                    JPDF_t;
 
  JFunction1D_t::JSupervisor supervisor(new JFunction1D_t::JDefaultResult(zero));
 
  const JPDFType_t pdf_t[] = { DIRECT_LIGHT_FROM_MUON,
                               SCATTERED_LIGHT_FROM_MUON,
                               DIRECT_LIGHT_FROM_EMSHOWERS,
                               SCATTERED_LIGHT_FROM_EMSHOWERS };
    
  const int  N = sizeof(pdf_t) / sizeof(pdf_t[0]);
 
  JPDF_t pdf[N];
 
  for (int i = 0; i != N; ++i) {
 
    try {
 
      const string file_name = getFilename(inputFile, pdf_t[i]);
      
      NOTICE("loading input from file " << file_name << "... " << flush);
      
      pdf[i].load(file_name.c_str());
      
      NOTICE("OK" << endl);
    }
    catch(const JException& error) {
      FATAL(error.what() << endl);
    }
    
    pdf[i].setExceptionHandler(supervisor);
  }
 
 
  vector<double>           X;        // transverse distance [m]
 
  if (X.empty()) {
    X.push_back(  5.0);
    X.push_back( 15.0);
    X.push_back( 25.0);
    X.push_back( 35.0);
    X.push_back( 45.0);
    X.push_back( 55.0);
    X.push_back( 65.0);
    X.push_back( 75.0);
    X.push_back( 85.0);
    X.push_back( 95.0);
  }
  
 
  TFile out(outputFile.c_str(), "recreate");
 
  vector<TH1D*> buffer;
 
  for (vector<double>::const_iterator i = X.begin(); i != X.end(); ++i) {
 
    ostringstream os;
 
    os << "t[" << *i << "]";
    
    buffer.push_back(new TH1D(os.str().c_str(), NULL, 5000, -250.0, +250.0));
  }
 
  vector<double> x;
 
  {
    vector<double>::const_iterator j = X.begin();
    vector<double>::const_iterator i = j++;
 
    x.push_back(*i - 0.5*(*j-*i));
 
    for (; j != X.end(); ++i, ++j) 
      x.push_back(0.5*(*i+*j));
 
    --i;
    --j;
 
    x.push_back(*j + 0.5*(*j-*i));
  }
 
  TH1D rms("rms", NULL, X.size(), &x[0]);
 
 
  JFunction1D_t::result_type p[4];
 
  JTimer timer;
 
  timer.reset();
 
  const double Tmin = -250.0;     // [ns]
  const double Tmax = +250.0;     // [ns]
 
  int n = 0;
 
  for (size_t i = 0; i != X.size(); ++i) {
 
    TH1D*        h = buffer[i];
    const double R = X[i];
 
    JSplineFunction1S_t zsp;
 
    for (int j = 1; j <= h->GetNbinsX(); ++j) {
 
      const double t1 = h->GetBinCenter(j);
 
      timer.start();
 
      ++n;
 
      for (int k = 0; k != 4; ++k) {
        p[k] = pdf[k](R, dir.getTheta(), dir.getPhi(), t1);
      }
 
      // integral [Tmin,t1]
      const double v = 
        p[0].v         +
        p[1].v         +
        p[2].v * E_GeV +
        p[3].v * E_GeV +
        R_Hz * 1e-9 * (t1   - Tmin);
 
      // integral [Tmin,Tmax]
      const double V = 
        p[0].V         +
        p[1].V         +
        p[2].V * E_GeV +
        p[3].V * E_GeV +
        R_Hz * 1e-9 * (Tmax - Tmin);
 
      // function value
      const double y = 
        p[0].f         +
        p[1].f         +
        p[2].f * E_GeV +
        p[3].f * E_GeV + 
        R_Hz * 1e-9;
 
      const double W = exp(-v) * y / (1.0 - exp(-V));
 
      timer.stop();
 
      zsp[t1] = W;
 
      h->SetBinContent(j,W);
    } 
 
    zsp.compile();
 
    try {
 
      JQuantiles result(zsp);
 
      const double sig = result.getFWHM() * 0.5 / sqrt(2.0*log(2.0));
 
      rms.Fill(R, sig);
 
      NOTICE("integral " << R << ' ' << result.getIntegral() << endl);
    }
    catch(const exception&) {}
  }
 
  if (n != 0) {
 
    const float w = 1.0 / (float) n;
 
    if (debug > notice_t) {
      timer.print(cout, w);
    }
  }
 
  out.Write();
  out.Close();
}