JEMShower.cc File Reference

Auxiliary program to draw npe as a function of EM-energy. More...

#include <string>
#include <iostream>
#include <iomanip>
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "JPhysics/JPDF.hh"
#include "JPhysics/KM3NeT.hh"
#include "JPhysics/JConstants.hh"
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
#include <cmath>
#include "JTools/JFunction1D_t.hh"
#include "JPhysics/JDispersion.hh"
#include "JPhysics/JPDFToolkit.hh"

Go to the source code of this file.

Functions
int	main (int argc, char **argv)

Detailed Description

Auxiliary program to draw npe as a function of EM-energy.

Author

mdejong

Definition in file JEMShower.cc.

Function Documentation

int main	(	int	argc,
		char **	argv
	)

Definition at line 147 of file JEMShower.cc.

{
  using namespace std;

  string      outputFile;
  int         numberOfPoints;
  double      epsilon;
  int         debug;

  try { 

    JParser<> zap("Auxiliary program to draw npe as a function of EM-energy.");

    zap['o'] = make_field(outputFile)      = "geanc.root";
    zap['n'] = make_field(numberOfPoints)  = 25;
    zap['e'] = make_field(epsilon)         = 1.0e-10;
    zap['d'] = make_field(debug)           = 2;

    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }
  
  
  using namespace JPP;

  const JPDF_C
    pdf(NAMESPACE::getPhotocathodeArea(),
        NAMESPACE::getQE,
        NAMESPACE::getAngularAcceptance,
        NAMESPACE::getAbsorptionLength,
        NAMESPACE::getScatteringLength,
        NAMESPACE::getScatteringProbability,
        NAMESPACE::getAmbientPressure(),
        getMinimalWavelength(),
        getMaximalWavelength(),
        numberOfPoints,
        epsilon);


  double xmin = -4.0;
  double xmax =  0.0;

  const double precision = 1.0e-10;

  while (fabs(xmax - xmin) > precision) {

    const double x = 0.5 * (xmin + xmax);
    const double E = pow(10.0, x);

    const double y = getEMShowerCorrection(E);

    if (y < 0.5) 
      xmin = x;
    else
      xmax = x;
  }

  const double EMIN =  pow(10.0, 0.5*(xmax + xmin));

  NOTICE("Threshold kinetic energy [GeV] " << sqrt((EMIN + MASS_ELECTRON) * (EMIN - MASS_ELECTRON)) << endl);

 
  TFile out(outputFile.c_str(), "recreate");

  TH1D h0("h0", NULL, 10000, -4.0, +4.0);
  TH1D h1("h1", NULL, 10000, -4.0, +4.0);
 
  for (int i = 1; i <= h0.GetNbinsX(); ++i) {

    const Double_t x = h0.GetBinCenter(i);
    const Double_t E = pow(10.0, x);

    h0.SetBinContent(i, E * geanc() * pdf.getNumberOfPhotons());
    h1.SetBinContent(i, getEMShowerCorrection(E));
  }

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