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	henyey_greenstein (const double ct)
double	rayleigh (const double ct)
double	getAngularAcceptance (const double x)
	Angular acceptence of PMT. More...
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

typedef std::pair<double,double> orientation

Definition at line 22 of file JDrawLED.cc.

Function Documentation

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

Definition at line 24 of file JDrawLED.cc.

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

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

Definition at line 25 of file JDrawLED.cc.

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

double henyey_greenstein ( const double  ct )

inline

Definition at line 28 of file JDrawLED.cc.

{
  static const double a = 0.924;

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

double rayleigh ( const double  ct )

inline

Definition at line 36 of file JDrawLED.cc.

{
  static const double a = 0.853;

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

double getAngularAcceptance ( const double  x )

inline

Angular acceptence of PMT.

Parameters

x	cosine of angle of incidence

Returns

probability

Definition at line 84 of file JDrawLED.cc.

{
  return ANTARES::getAngularAcceptance(x);
}

int main	(	int	argc,
		char **	argv
	)

Definition at line 96 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();
}