JDrawLED.cc

Go to the documentation of this file.

 
#include <string>
#include <iostream>
#include <iomanip>
 
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
 
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
#include "JPhysics/JLED.hh"
#include "JPhysics/JPDFToolkit.hh"
#include "JPhysics/Antares.hh"
#include "JPhysics/KM3NeT.hh"
 
#include "JTools/JSpline.hh"
#include "JTools/JQuantiles.hh"
#include "JTools/JToolsToolkit.hh"
 
 
typedef std::pair<double,double> orientation;
 
std::istream& operator>>(std::istream& in, orientation& x) { return in >> x.first >> x.second; }
std::ostream& operator<<(std::ostream& out, const orientation& x) { return out << x.first << ' ' << x.second; }
 
 
/**
 * Light yield from %LED (number of p.e. per unit solid angle per unit time).
 */
class LED : public JPHYSICS::JAbstractLED {
public:
  /**
   * Constructor.
   */
  LED()
  {}
 
  /**
   * Light yield from %LED (number of p.e. per unit solid angle per unit time).
   *
   * \param  ct         zenith  angle of emission
   * \param  phi        azimuth angle of emission
   * \param  dt         time of emission [ns]
   * \return            d^2P / dOmega dt [npe/ns/sr]
   */
  double getLightFromLED(const double ct, 
                         const double phi,
                         const double dt) const
  {
    using namespace JPP;
 
    static const double sigma = 2.0;
    
    const double x = dt / sigma;
    
    return exp(-0.5*x*x) / (sigma*sqrt(2*PI)) / (4*PI);
  }
};
 
 
/**
 * Angular acceptence of PMT
 *
 * \param  x      cosine of angle of incidence
 * \return        probability
 */
inline double getAngularAcceptance(const double x)
{
  return ANTARES::getAngularAcceptance(x);
}
 
 
/**
 * \file
 *
 * Example program to draw PDF from %LED beacon.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
 
  string      outputFile;
  int         numberOfPoints;
  int         number_of_points;
  double      epsilon;
  double      D;
  double      ct;
  orientation dir;
  double      wavelength;
  double      absorptionLength;
  double      scatteringLength;
  int         debug;
 
  try { 
 
    JParser<> zap("Example program to draw PDF from LED beacon.");
 
    zap['o'] = make_field(outputFile)       = "led.root";
    zap['n'] = make_field(numberOfPoints)   = 25;
    zap['N'] = make_field(number_of_points) = 25;
    zap['e'] = make_field(epsilon)          = 1.0e-10;
    zap['A'] = make_field(absorptionLength) = 50.0;
    zap['S'] = make_field(scatteringLength) = 50.0;
    zap['R'] = make_field(D);
    zap['c'] = make_field(ct);
    zap['D'] = make_field(dir);
    zap['w'] = make_field(wavelength)       = 470;
    zap['d'] = make_field(debug)            = 2;
 
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }
 
 
  const double theta = dir.first;
  const double phi   = dir.second;
  
  using namespace JPP;
 
  // set global parameters
 
  const double P_atm = 240.0;       // ambient pressure [Atm]
  const double tmin  = -10.0;       // minimal time of emission [ns]
  const double tmax  = +10.0;       // maximal time of emission [ns]
  const double A     =  440e-4;     // ANTARES photo-cathode area [m2]
  const double R_Hz  =  0.0e3;      // singles rate [Hz]
 
  const JDispersion dispersion(P_atm);
 
  //const double ng    = dispersion.getIndexOfRefractionGroup(wavelength);
  const double lm    = scatteringLength / 0.83;
  const double lr    = scatteringLength / 0.17;
 
  const double cs    = 0.83 * 0.92;  // average cosine scattering angle
 
  const double l_abs = absorptionLength;
  const double ls    = scatteringLength;
  const double l_att = l_abs * lr / (l_abs + lr);
 
 
  LED* led = new LED();
 
  cout << "Rayleigh scattering length " << lr    << " m" << endl;
  cout << "Mie      scattering length " << lm    << " m" << endl;
  cout << "Absorption length          " << l_abs << " m" << endl;
 
  const JLED_C
    pdfMie(A,
           led,
           ANTARES::getQE,
           getAngularAcceptance,
           l_abs,
           lm,
           henyey_greenstein,
           P_atm,
           wavelength,
           tmin,
           tmax,
           JCotangent(number_of_points),
           numberOfPoints,
           epsilon);
  
 
  const JLED_C
    pdfRayleigh(A,
                led,
                ANTARES::getQE,
                getAngularAcceptance,
                l_att,
                lr,
                rayleigh,
                P_atm,
                wavelength,
                tmin,
                tmax,
                JCotangent(number_of_points),
                numberOfPoints,
                epsilon);
 
 
  TFile out(outputFile.c_str(), "recreate");
 
  TH1D h0("h0", NULL, 430, -15.0, +200.0);
  TH1D h1("h1", NULL, 430, -15.0, +200.0);
  TH1D h2("h2", NULL, 430, -15.0, +200.0);
  TH1D h3("h3", NULL, 430, -15.0, +200.0);
  TH1D ha("ha", NULL, 430, -15.0, +200.0);
 
  JSplineFunction1S_t f[4];
  JSplineFunction1S_t f1;
 
 
  for (int i = 1; i <= h1.GetNbinsX(); ++i) {
 
    const double t1 = h1.GetBinCenter(i);
 
    const double F1 = pdfMie     .getDirectLightFromLED   (D, ct, theta, phi, t1);
    const double F2 = pdfMie     .getScatteredLightFromLED(D, ct, theta, phi, t1);
    const double F3 = pdfRayleigh.getScatteredLightFromLED(D, ct, theta, phi, t1);
 
      
    h0.SetBinContent(i, F1 + F2 + F3);
    h1.SetBinContent(i, F1);
    h2.SetBinContent(i, F2);
    h3.SetBinContent(i, F3);
 
    f[0][t1] = F1;
    f[1][t1] = F2;
    f[2][t1] = F3;
    f[3][t1] = F1 + F2 + F3;
 
    f1[t1] = F1 + F2 + F3;
  }
 
  f1.compile();
 
  for (int i = 3; i != sizeof(f)/sizeof(f[0]); ++i) {
 
    f[i].compile();
 
    JQuantiles quantiles(f[i]);
 
 
   const double lz = l_abs * ls / (l_abs*(1.0-cs) + ls);
 
    NOTICE("int  " << quantiles.getIntegral() * D*D * exp(D/lz) << endl);
 
    DEBUG("D     " << D                       << endl);
    DEBUG("ct    " << ct                      << endl);
    DEBUG("theta " << theta                   << endl);
    DEBUG("phi   " << phi                     << endl);
    DEBUG("int   " << quantiles.getIntegral() << endl);
    DEBUG("t1    " << quantiles.getX()        << endl);
    DEBUG("max   " << quantiles.getY()        << endl);
    DEBUG("FWHM  " << quantiles.getFWHM()     << endl);
  }
 
 
  const double Tmin = ha.GetXaxis()->GetXmin();
  const double Tmax = ha.GetXaxis()->GetXmax();
 
  const double V = f1.rbegin()->getIntegral()  +  R_Hz * 1e-9 * (Tmax - Tmin);   // integral [Tmin,Tmax]
 
  for (int i = 1; i <= ha.GetNbinsX(); ++i) {
 
    const double t1 = ha.GetBinCenter(i);
 
    JSplineFunction1S_t::result_type p = f1(t1);
 
    double v = p.v  +  R_Hz * 1e-9 * (t1   - Tmin);
    double y = p.f  +  R_Hz * 1e-9;                   // function value
    
    const double W = exp(-v) * y / (1.0 - exp(-V));
    
    ha.SetBinContent(i,W);      
  }
 
  delete led;
 
  out.Write();
  out.Close();
}