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. More...
double	scatteringLengthFactor

Detailed Description

Auxiliary program to draw PDF of Cherenkov light from muon.

Author

mdejong, bofearraigh

Definition in file JDrawPDF.cc.

Function Documentation

double getAbsorptionLength ( const double  lambda )

inline

Definition at line 31 of file JDrawPDF.cc.

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

double getScatteringLength ( const double  lambda )

inline

Definition at line 37 of file JDrawPDF.cc.

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

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();

  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("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);

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

Variable Documentation

double absorptionLengthFactor

Scaling of absorption and scattering length.

Definition at line 28 of file JDrawPDF.cc.

double scatteringLengthFactor

Definition at line 29 of file JDrawPDF.cc.