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; }


inline double henyey_greenstein(const double ct)
{
  static const double a = 0.924;

  return ANTARES::henyey_greenstein(a, ct);
}


inline double rayleigh(const double ct)
{
  static const double a = 0.853;

  return ANTARES::rayleigh(a, ct);
}


/**
 * 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();
}