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 "JTools/JElement.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"


/**
 * Scaling of absorption and scattering length.
 */
double  absorptionLengthFactor;
double  scatteringLengthFactor;

inline double getAbsorptionLength(const double lambda)
{
  return absorptionLengthFactor * NAMESPACE::getAbsorptionLength(lambda);
}


inline double getScatteringLength(const double lambda)
{
  return scatteringLengthFactor * NAMESPACE::getScatteringLength(lambda);
}


/**
 * \file
 *
 * Auxiliary program to draw PDF of Cherenkov light from muon.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;

  typedef JAbstractHistogram<double> JHistogram_t;

  string        outputFile;
  int           numberOfPoints;
  double        epsilon;
  double        E;
  double        R;
  JAngle3D      dir;
  vector<int>   function;
  JHistogram_t  histogram;
  int           debug;

  try { 

    JProperties properties;

    properties.insert(gmake_property(JPHYSICS::MODULE_RADIUS_M));
  
    JParser<> zap("Auxiliary program to draw PDF of Cherenkov light from muon.");

    zap['@'] = make_field(properties)             = 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['A'] = make_field(absorptionLengthFactor, "scaling factor")             = 1.0;
    zap['S'] = make_field(scatteringLengthFactor, "scaling factor")             = 1.0;
    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['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,
        getAbsorptionLength,
        getScatteringLength,
        NAMESPACE::getScatteringProbability,
        NAMESPACE::getAmbientPressure(),
        getMinimalWavelength(),
        getMaximalWavelength(),
        numberOfPoints,
        epsilon);


  if (outputFile == "") {

    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,2)      << dir.getTheta() << ' '
             << FIXED(5,2)      << dir.getPhi()   << ' '
             << FIXED(5,1)      << dt             << ' '
             << SCIENTIFIC(9,3) << pdf.getLightFromMuon(*F, E, R, dir.getTheta(), dir.getPhi(), dt) * E << endl;
      }
    }

    return 0;
  }
  
 
  TFile out(outputFile.c_str(), "recreate");

  //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.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;
 
  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, R, dir.getTheta(), dir.getPhi(), dt);
    }

    h0.SetBinContent(i, value);

    f1[dt] = value;
  }

  f1.setExceptionHandler(new JSplineFunction1S_t::JDefaultResult(JMATH::zero));
  f1.compile();

  const double T = 5;  // [ns]

  JQuantiles quantiles(f1);

  const double t1 = quantiles.getX();
  const double y  = f1(t1 + T).v - f1(t1 - T).v;

  DEBUG("E     " << E                       << endl);
  DEBUG("R     " << R                       << 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);  

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