JPlotPDF_L1.cc

Go to the documentation of this file.

 
#include <string>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <vector>
 
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
 
#include "JTools/JFunction1D_t.hh"
#include "JTools/JFunctionalMap_t.hh"
#include "JTools/JAbstractHistogram.hh"
#include "JPhysics/JPDFTransformer.hh"
#include "JPhysics/JPDFTable.hh"
#include "JPhysics/JPDFTypes.hh"
#include "JPhysics/JPDFToolkit.hh"
#include "JGeometry3D/JAngle3D.hh"
#include "JGeometry3D/JRotation3D.hh"
#include "JDetector/JModule.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JDetectorSupportkit.hh"
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
 
/**
 * Auxiliary program to determine PDF of L1 hit.
 *
 * \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;
  double         R_Hz;
  JAngle3D       dir;
  double         T_ns;
  JHistogram_t   histogram;
  int            debug;
 
  try { 
 
    JParser<> zap("Auxiliary program to determine PDF of L1 hit.");
 
    zap['f'] = make_field(inputFile);
    zap['o'] = make_field(outputFile)          = "pdf.root";
    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['T'] = make_field(T_ns,      "time window [ns]")           = 10.0;          // [ns]
    zap['B'] = make_field(R_Hz,      "background rate [Hz]")       =  0.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 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_DELTARAYS,
                               SCATTERED_LIGHT_FROM_DELTARAYS,
                               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);
  }
 
  if (!histogram.is_valid()) {
 
    histogram = JHistogram_t(-50.0, +50.0);
 
    histogram.setBinWidth(0.1);
  }
  
 
  JModule module = getModule<JKM3NeT_t>(1);
 
  module.rotate(JRotation3D(dir));
 
 
  TFile out(outputFile.c_str(), "recreate");
 
  TH1D h0("h0", NULL, histogram.getNumberOfBins(), histogram.getLowerLimit(), histogram.getUpperLimit());
 
  JSplineFunction1S_t zsp;
 
  for (double t1 = histogram.getLowerLimit(); t1 <= histogram.getUpperLimit(); t1 += 0.1) {
 
    const double t2 = t1 + T_ns;
 
    JFunction1D_t::result_type y1;
    JFunction1D_t::result_type y2;
 
    for (int i = 0; i != N; ++i) {
 
      JFunction1D_t::result_type p1;
      JFunction1D_t::result_type p2;
 
      // add PMTs
 
      for (const auto& pmt : module) {
        p1 += pdf[i](R, pmt.getTheta(), pmt.getPhi(), t1);
        p2 += pdf[i](R, pmt.getTheta(), pmt.getPhi(), t2);
      }
 
      if        (is_deltarays(pdf_t[i])) {
        p1 *= getDeltaRaysFromMuon(E);
        p2 *= getDeltaRaysFromMuon(E);
      } else if (is_bremsstrahlung(pdf_t[i])) {
        p1 *= E;
        p2 *= E;
      }
 
      y1 += p1;
      y2 += p2;
    }
 
    // add background
 
    y1.f += R_Hz * 1e-9;                                         // function value
    y1.v += R_Hz * 1e-9 * (t1 - histogram.getLowerLimit());      // integral [Tmin,t1]
    y2.v += R_Hz * 1e-9 * (t2 - histogram.getLowerLimit());      // integral [Tmin,t2]
 
    zsp[t1] = y1.f * (1.0 - exp(y1.v - y2.v));                    // 1 - P(0)
  } 
 
  zsp.compile();
 
  for (int ix = 1; ix <= h0.GetNbinsX(); ++ix) {
 
    const double t1 = h0.GetBinCenter(ix);
 
    try {
 
      const JFunction1D_t::result_type result = zsp(t1);
 
      const double W = exp(-result.v) * result.f / (1.0 - exp(-result.V));
 
      h0.SetBinContent(ix, W);
    }
    catch(const exception&) {}
  }
  
  out.Write();
  out.Close();
}