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/JCDFTable.hh"
#include "JPhysics/JPDFTypes.hh"
#include "JPhysics/JDeltaRays.hh"
#include "JGeometry3D/JAngle3D.hh"
#include "JROOT/JManager.hh"
#include "JROOT/JRootToolkit.hh"
#include "JGizmo/JGizmoToolkit.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;
 
  vector<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)                                   = 2;
 
    zap(argc, argv);
  }
  catch(const exception& error) {
    FATAL(error.what() << endl);
  }
 
 
  typedef JHermiteSplineFunction1D_t                              JFunction1D_t;
  //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;
 
  const int N = inputFile.size();
 
  int    type[N];
  JCDF_t cdf [N];
 
  try {
 
    for (int i = 0; i != N; ++i) {
 
      NOTICE("loading input from file " << inputFile[i] << "... " << flush);
 
      type[i] = getPDFType(inputFile[i]);
 
      cdf [i].load(inputFile[i].c_str());
 
      NOTICE("OK" << endl);
    }
  }
  catch(const JException& error) {
    FATAL(error.what() << endl);
  }
 
 
  if (outputFile == "") {
 
    for ( ; ; ) {
 
      double  x;
 
      cout << "> " << flush;
      cin  >> x;
 
      cout << " --> ";
 
      for (int i = 0; i != N; ++i) {
      
        try {
        
          double npe = cdf[i].getNPE (R, dir.getTheta(), dir.getPhi());
          double t   = cdf[i].getTime(R, dir.getTheta(), dir.getPhi(), x);
 
          if        (is_bremsstrahlung(type[i])) {
            npe *= E;
          } else if (is_deltarays(type[i])) {
            npe *= JDeltaRays::getEnergyLossFromMuon(E);
          }
 
          cout << ' ' << FIXED(6,2) << t << ' ' << FIXED(5,2) << npe;
        }
        catch(const exception& error) {
          ERROR(error.what() << endl);
        }
      }
 
      cout << endl;
    }
 
    return 0;
  }
   
 
  int function = 0;
 
  if (inputFile.size() == 1 && getPDFType(inputFile[0]) == DIRECT_LIGHT_FROM_MUON) {
    function = 1;
  }
 
  TH1D* h0 = NULL;
 
  //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()) {
 
    vector<Double_t> X;
 
    double x = -100.0;
 
    {
      for ( ; x <  -10.0; x +=  5.0) { X.push_back(t0 + x); }
      for ( ; x <  +20.0; x +=  1.0) { X.push_back(t0 + x); }
      for ( ; x <  +50.0; x +=  2.0) { X.push_back(t0 + x); }
    }
    if (function != 1) {
      for ( ; x < +100.0; x +=  5.0) { X.push_back(t0 + x); }
      for ( ; x < +250.0; x += 10.0) { X.push_back(t0 + x); }
      for ( ; x < +500.0; x += 25.0) { X.push_back(t0 + x); }
      for ( ; x < +900.0; x += 50.0) { X.push_back(t0 + x); }
    }
 
    h0 = new TH1D("h0", NULL,  X.size() - 1, X.data());
 
  } else {
 
    h0 = new TH1D("h0", NULL, histogram.getNumberOfBins(), histogram.getLowerLimit(), histogram.getUpperLimit());
  }
 
  h0->Sumw2();
  
  JManager<int, TH1D> H1(new TH1D("h1[%]", NULL, 100000, 0.0, +1.0));
  
  for (int i = 0; i != N; ++i) {
    
    for (int j = 1; j <= H1->GetNbinsX(); ++j) {
    
      try {
      
        const double x = H1->GetBinCenter(j);
        const double t = cdf[i].getTime(R, dir.getTheta(), dir.getPhi(), x);
      
        H1[i]->SetBinContent(j, t);
      }
      catch(const exception& error) {
        ERROR(error.what() << endl);
      }
    }
  }
    
  if (numberOfEvents > 0) {
  
    JTimer timer;
    
    timer.reset();
    
    for (int counter = 0; counter != numberOfEvents; ++counter) {
      
      if (counter%1000== 0) { 
        STATUS(setw(10) << counter << '\r'); DEBUG(endl);
      }
      
      timer.start();
      
      for (int i = 0; i != N; ++i) {
 
        try {
        
          double npe = cdf[i].getNPE(R, dir.getTheta(), dir.getPhi());
 
          if        (is_bremsstrahlung(type[i])) {
            npe *= E;
          } else if (is_deltarays(type[i])) {
            npe *= JDeltaRays::getEnergyLossFromMuon(E);
          }
 
          for (int j = gRandom->Poisson(npe); j != 0; --j) {
          
            const double x = gRandom->Rndm();
            const double t = cdf[i].getTime(R, dir.getTheta(), dir.getPhi(), x);
          
            h0->Fill(t + t0);
          }
        }
        catch(const exception& error) {
          NOTICE(error.what() << endl);
        }
      }
      
      timer.stop();
    }
    STATUS(endl);
  
    const double W = 1.0 / (double) numberOfEvents;
    
    if (debug > JEEP::status_t) {
      timer.print(cout, W);
    }
 
    convertToPDF(*h0, "WE", W);
  }
 
  TFile out(outputFile.c_str(), "recreate");
 
  out << *h0 << H1;
  
  out.Write();
  out.Close();
}