JDrawLED.cc File Reference

Example program to draw PDF from LED beacon. More...

#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"

Go to the source code of this file.

Classes
class	LED
	Light yield from LED (number of p.e. More...

Typedefs
typedef std::pair< double, double >	orientation

Functions
std::istream &	operator>> (std::istream &in, orientation &x)
std::ostream &	operator<< (std::ostream &out, const orientation &x)
double	getAngularAcceptance (const double x)
	Angular acceptence of PMT.
int	main (int argc, char **argv)

Detailed Description

Example program to draw PDF from LED beacon.

Author

mdejong

Definition in file JDrawLED.cc.

Typedef Documentation

orientation

typedef std::pair<double,double> orientation

Definition at line 22 of file JDrawLED.cc.

Function Documentation

operator>>()

std::istream & operator>>	(	std::istream &	in,
		orientation &	x )

Definition at line 24 of file JDrawLED.cc.

{ return in >> x.first >> x.second; }

operator<<()

std::ostream & operator<<	(	std::ostream &	out,
		const orientation &	x )

Definition at line 25 of file JDrawLED.cc.

{ return out << x.first << ' ' << x.second; }

getAngularAcceptance()

double getAngularAcceptance ( const double x )

inline

Angular acceptence of PMT.

Parameters

x	cosine of angle of incidence

Returns

probability

Definition at line 68 of file JDrawLED.cc.

{
  return ANTARES::getAngularAcceptance(x);
}

main()

int main	(	int	argc,
		char **	argv )

Definition at line 80 of file JDrawLED.cc.

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