JDrawPDF.cc File Reference

Auxiliary program to draw PDF of Cherenkov light from muon. More...

#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"
#include "JPhysics/JGeane.hh"

Go to the source code of this file.

Functions
double	getAbsorptionLength (const double lambda)
double	getScatteringLength (const double lambda)
int	main (int argc, char **argv)

Variables
double	absorptionLengthFactor
	Scaling of absorption and scattering length.
double	scatteringLengthFactor

Detailed Description

Auxiliary program to draw PDF of Cherenkov light from muon.

Author

mdejong, bofearraigh

Definition in file JDrawPDF.cc.

Function Documentation

getAbsorptionLength()

double getAbsorptionLength ( const double lambda )

inline

Definition at line 31 of file JDrawPDF.cc.

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

getScatteringLength()

double getScatteringLength ( const double lambda )

inline

Definition at line 37 of file JDrawPDF.cc.

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

main()

int main	(	int	argc,
		char **	argv )

Definition at line 49 of file JDrawPDF.cc.

{
  using namespace std;
  using namespace JPP;
 
  typedef JAbstractHistogram<double> JHistogram_t;
 
  string        outputFile;
  int           numberOfPoints;
  double        epsilon;
  double        E;
  double        R;
  double        z;
  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 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,
        getAbsorptionLength,
        getScatteringLength,
        NAMESPACE::getScatteringProbability,
        NAMESPACE::getAmbientPressure(),
        getMinimalWavelength(),
        getMaximalWavelength(),
        numberOfPoints,
        epsilon);
 
  const double z_emission = z - R/getTanThetaC(); // emission point z-position
  const double E_emission  = gWater.getE(E, z_emission ); // Get energy of muon at emission point
 
  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,1)      << z              << ' '
             << FIXED(5,2)      << dir.getTheta() << ' '
             << FIXED(5,2)      << dir.getPhi()   << ' '
             << FIXED(5,1)      << dt             << ' '
             << SCIENTIFIC(9,3) << pdf.getLightFromMuon(*F, E_emission, R, dir.getTheta(), dir.getPhi(), dt) * E_emission << 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_emission > MASS_MUON* (1/COS_THETA_C_WATER) ) { // muon emission energy has to be above Cherenkov threshold
 
    if (z_emission >= 0 && z_emission <= 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_emission, 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_emission" << E_emission              << 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();
}

Variable Documentation

absorptionLengthFactor

double absorptionLengthFactor

Scaling of absorption and scattering length.

Definition at line 28 of file JDrawPDF.cc.

scatteringLengthFactor

double scatteringLengthFactor

Definition at line 29 of file JDrawPDF.cc.