JPlotCDF.cc

Go to the documentation of this file.

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

#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TProfile.h"
#include "TRandom3.h"

#include "JTools/JFunction1D_t.hh"
#include "JTools/JFunctionalMap_t.hh"
#include "JTools/JAbstractHistogram.hh"
#include "JPhysics/JPDFTransformer.hh"
#include "JPhysics/JCDFTable.hh"
#include "JPhysics/JPDFTypes.hh"
#include "JGeometry3D/JAngle3D.hh"
#include "Jeep/JTimer.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"


/**
 * \file
 *
 * Program to verify generation of arrival times of Cherenkov photons from a muon using tabulated CDF.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;

  typedef JAbstractHistogram<double> JHistogram_t;

  string        inputFile;
  string        outputFile;
  double        E;
  double        R;
  JAngle3D      dir;
  int           numberOfEvents;
  JHistogram_t  histogram;
  int           debug;

  try { 

    JParser<> zap("Program to verify generation of arrival times of Cherenkov photons from a muon using tabulated CDF.");

    zap['f'] = make_field(inputFile);
    zap['o'] = make_field(outputFile)                              = "";
    zap['E'] = make_field(E,         "muon energy [GeV]")          = 1.0;
    zap['R'] = make_field(R,         "distance of approach [m]");
    zap['D'] = make_field(dir,       "(theta, phi) of PMT [rad]");
    zap['N'] = make_field(numberOfEvents)                          = 0;
    zap['H'] = make_field(histogram, "histogram binning")          = JHistogram_t();
    zap['d'] = make_field(debug)                                   = 0;

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


  typedef JSplineFunction1D_t                                     JFunction1D_t;
  //typedef JPolint1Function1D_t                                    JFunction1D_t;
  typedef JMAPLIST<JPolint1FunctionalMap,
                   JPolint1FunctionalGridMap,
                   JPolint1FunctionalGridMap>::maplist            JMapList_t;
  typedef JCDFTable<JFunction1D_t, JMapList_t>                    JCDF_t;

  JCDF_t cdf;


  try {

    NOTICE("loading input from file " << inputFile << "... " << flush);

    cdf.load(inputFile.c_str());

    NOTICE("OK" << endl);
  }
  catch(const JException& error) {
    FATAL(error.what() << endl);
  }


  if (outputFile == "") {

    for ( ; ; ) {

      double  x;

      cout << "> " << flush;
      cin  >> x;

      try {
        
        const double npe = cdf.getNPE (R, dir.getTheta(), dir.getPhi()) * E;        
        const double t   = cdf.getTime(R, dir.getTheta(), dir.getPhi(), x);

        cout << " --> " << t << ' ' << npe;
      }
      catch(const exception& error) {
        ERROR(error.what() << endl);
      }

      cout << endl;
    }

  } else {
  
   
    TFile out(outputFile.c_str(), "recreate");
  
    int function = 1;

    if (inputFile.find(getLabel(DIRECT_LIGHT_FROM_MUON)) == string::npos) {
      function = 0;
    }

    //const double z0   = -100.0;                                 // [m]
    //const double t0   = (-z0 +  R * getTanThetaC()) / C;        // time [ns]
    const double t0   =  0.0;                                   // time [ns]

    if (!histogram.is_valid()) {

      if (function == 0) {

        histogram = JHistogram_t(t0 - 20.0, t0 + 50.0);

        histogram.setBinWidth(0.1);

      } else {

        histogram = JHistogram_t(t0 - 20.0, t0 + 500.0);

        histogram.setBinWidth(0.5);
      }
    }

    TH1D h0("h0", NULL, histogram.getNumberOfBins(), histogram.getLowerLimit(), histogram.getUpperLimit());
    TH1D h1("h1", NULL, 100000, 0.0, +1.0);

    h0.Sumw2();
  
    for (int j = 1; j <= h1.GetNbinsX(); ++j) {
    
      try {
      
        const double x = h1.GetBinCenter(j);
        const double t = cdf.getTime(R, dir.getTheta(), dir.getPhi(), x);
      
        h1.SetBinContent(j, t);
      }
      catch(const exception& error) {
        ERROR(error.what() << endl);
      }
    }
    
    if (numberOfEvents > 0) {
    
  
      JTimer timer;
    
      timer.reset();
    
      for (int i = 0; i != numberOfEvents; ++i) {
      
        if (i%1000== 0) { NOTICE(i << '\r'); }
      
        timer.start();
      
        try {
        
          const double npe = cdf.getNPE(R, dir.getTheta(), dir.getPhi()) * E;
          
          for (int j = gRandom->Poisson(npe); j != 0; --j) {
          
            const double x = gRandom->Rndm();
            const double t = cdf.getTime(R, dir.getTheta(), dir.getPhi(), x);
          
            h0.Fill(t + t0);
          }
        }
        catch(const exception& error) {
          NOTICE(error.what() << endl);
        }
      
        timer.stop();
      }
      NOTICE(endl);
  
      const double w = 1.0 / (double) numberOfEvents;
    
      if (debug > JEEP::notice_t) {
        timer.print(cout, w);
      }

      // normalise histogram (dP/dt)
    
      for (int j = 1; j <= h0.GetNbinsX(); ++j) {
      
        const double y  = h0.GetBinContent(j);
        const double z  = h0.GetBinError  (j);
        const double dx = h0.GetBinWidth  (j);
      
        h0.SetBinContent(j, y / (dx*numberOfEvents*E));
        h0.SetBinError  (j, z / (dx*numberOfEvents*E));
      }
    }

    out.Write();
    out.Close();
  }
}