JDrawPDF.cc

Go to the documentation of this file.

 
#include <string>
#include <iostream>
#include <iomanip>
 
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
 
#include "JPhysics/JPDF.hh"
#include "JPhysics/Antares.hh"
#include "JPhysics/KM3NeT.hh"
#include "JPhysics/JPDFSupportkit.hh"
#include "JPhysics/JGeane.hh"
#include "JTools/JFunction1D_t.hh"
#include "JTools/JQuantiles.hh"
#include "JTools/JAbstractHistogram.hh"
#include "JGeometry3D/JAngle3D.hh"
#include "Jeep/JProperties.hh"
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
 
/**
 * \file
 *
 * Auxiliary program to draw PDF of Cherenkov light from muon.
 * \author mdejong, bofearraigh
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
 
  typedef JAbstractHistogram<double> JHistogram_t;
 
  string                  outputFile;
  int                     numberOfPoints;
  double                  epsilon;
  JAbsorptionLength       absorptionLength;
  JScatteringLength       scatteringLength;
  JScatteringProbability  scatteringProbability;
  double                  E;
  double                  R;
  double                  z;
  JAngle3D                dir;
  vector<int>             function;
  JHistogram_t            histogram;
  int                     debug;
 
  try { 
 
    JProperties properties;
 
    properties.insert(gmake_property(MODULE_RADIUS_M));
    properties.insert(gmake_property(absorptionLength));
    properties.insert(gmake_property(scatteringLength));
    properties.insert(gmake_property(scatteringProbability));
  
    JParser<> zap("Auxiliary program to draw PDF of Cherenkov light from muon.");
 
    zap['@'] = make_field(properties, endl
                          << "possible options absorptionLength: "      << get_keys(absorptionLength)      << endl 
                          << "possible options scatteringLength: "      << get_keys(scatteringLength)      << endl
                          << "possible options scatteringProbability: " << get_keys(scatteringProbability) << endl)  = JPARSER::initialised();
    zap['o'] = make_field(outputFile)                                                = "";
    zap['n'] = make_field(numberOfPoints,         "points for integration")          = 25;
    zap['e'] = make_field(epsilon,                "precision for integration")       = 1.0e-10;
    zap['E'] = make_field(E,                      "muon energy at vertex [GeV]")     = 1.0;
    zap['R'] = make_field(R,                      "distance of approach [m]");
    zap['z'] = make_field(z,                      "PMT z-position [m]");
    zap['D'] = make_field(dir,                    "(theta, phi) of PMT [rad]");
    zap['F'] = make_field(function,               "PDF type");
    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);
  }
  
 
  const JPDF_C
    pdf(NAMESPACE::getPhotocathodeArea(),
        NAMESPACE::getQE,
        NAMESPACE::getAngularAcceptance,
        JAbsorptionLength::getAbsorptionLength,
        JScatteringLength::getScatteringLength,
        JScatteringProbability::getScatteringProbability,
        NAMESPACE::getAmbientPressure(),
        getMinimalWavelength(),
        getMaximalWavelength(),
        numberOfPoints,
        epsilon);
 
  const double z_0  =  z  -  R/getTanThetaC();        // emission point z-position
  const double E_0  =  gWater.getE(E, z_0);           // energy of muon at emission point
 
 
  if (outputFile == "") {
 
    cout << "enter time (^C to exit) > " << flush;
 
    for (double dt; cin >> dt; ) {
 
      for (vector<int>::const_iterator F = function.begin(); F != function.end(); ++F) {
 
        cout << setw(2)         << *F             << ' '
             << SCIENTIFIC(7,1) << E              << ' '
             << FIXED(5,1)      << R              << ' '
             << FIXED(5,1)      << z              << ' '
             << FIXED(5,2)      << dir.getTheta() << ' '
             << FIXED(5,2)      << dir.getPhi()   << ' '
             << FIXED(5,1)      << dt             << ' '
             << SCIENTIFIC(9,3) << pdf.getLightFromMuon(*F, E_0, R, dir.getTheta(), dir.getPhi(), dt) << endl;
      }
    }
 
    return 0;
  }
  
 
  TFile out(outputFile.c_str(), "recreate");
 
  const double t0   =  0.0;                                   // time [ns]
 
  if (!histogram.is_valid()) {
 
    if (function.size() == 1 && function[0] == DIRECT_LIGHT_FROM_MUON) {
 
      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());
 
  JSplineFunction1S_t f1;
      
  if (E_0 >= MASS_MUON* (1.0/SIN_THETA_C_WATER)) {              // muon emission energy has to be above Cherenkov threshold
 
    if (z_0 >= 0 && z_0 <= gWater(E)) {                         // emission point is between start and end point of muon
 
      for (int i = 1; i <= h0.GetNbinsX(); ++i) {
 
        const double dt = h0.GetBinCenter(i) - t0;
 
        double value = 0.0;
 
        for (vector<int>::const_iterator F = function.begin(); F != function.end(); ++F) {
          value += pdf.getLightFromMuon(*F, E_0, R, dir.getTheta(), dir.getPhi(), dt);
        }
 
        h0.SetBinContent(i, value);
 
        f1[dt] = value;
      }
    }
  }
 
  f1.setExceptionHandler(new JSplineFunction1S_t::JDefaultResult(JMATH::zero));
  f1.compile();
 
  try {
 
    const double T_ns = 5;  // [ns]
 
    JQuantiles quantiles(f1);
 
    const double t1 = quantiles.getX();
    const double y  = f1(t1 + T_ns).v - f1(t1 - T_ns).v;
 
    DEBUG("E     "     << E                       << endl);
    DEBUG("E_0   "     << E_0                     << endl);
    DEBUG("R     "     << R                       << endl);
    DEBUG("z     "     << z                       << endl);
    DEBUG("theta "     << dir.getTheta()          << endl);
    DEBUG("phi   "     << dir.getPhi()            << endl);
    DEBUG("int   "     << quantiles.getIntegral() << endl);
    DEBUG("t1    "     << t1                      << endl);
    DEBUG("max   "     << quantiles.getY()        << endl);
    DEBUG("FWHM  "     << quantiles.getFWHM()     << endl);
    DEBUG("int[] "     << y                       << endl);
  }
  catch(const exception&) {}
 
  out.Write();
  out.Close();
}