JPMT.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>

#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"

#include "JLang/Jpp.hh"
#include "JPhysics/Antares.hh"
#include "JPhysics/KM3NeT.hh"
#include "JPhysics/KM3NeT2D.hh"
#include "JROOT/JManager.hh"
#include "JROOT/JRootToolkit.hh"
#include "JGizmo/JGizmoToolkit.hh"

#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"


/**
 * \file
 *
 * Auxiliary program to plot QE and angular acceptance of PMT.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;

  string   outputFile;
  bool     pdf;
  int      debug;

  try {

    JParser<> zap("Auxiliary program to plot QE and angular acceptance of PMT.");

    zap['o'] = make_field(outputFile)  = "pmt.root";
    zap['P'] = make_field(pdf);
    zap['d'] = make_field(debug)       = 2;

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

  using namespace NAMESPACE;
  using namespace JPP;


  TH1D h0(MAKE_CSTRING("QE["  << getNamespace() << "]"),  NULL, 340, 280.0, 900.0);
  TH1D h1(MAKE_CSTRING("PMT[" << getNamespace() << "]"),  NULL, 1000, -1.0, +1.0);

  typedef JManager<double, TH1D>   JManager_t;
  
  JManager_t qe ((TH1D*) h0.Clone("QE[%]"),  '%', JFormat_t(3, 1, ios::fixed));
  JManager_t pmt((TH1D*) h1.Clone("PMT[%]"), '%', JFormat_t(4, 0, ios::fixed));

  for (int i = 1; i <= h0.GetNbinsX(); ++i) {

    const double y = h0.GetBinCenter (i);

    h0.SetBinContent(i, getQE(y));

    for (double x = -1.0; x <= +1.0; x += 0.1) {
      if (KM3NET::getAngularAcceptance(x) > 0.0) {
        qe [floor(x*10.0)/10.0]->SetBinContent(i, KM3NET2D::getPhotocathodeArea2D(x, y) / KM3NET::getPhotocathodeArea() / KM3NET::getAngularAcceptance(x));
      }
    }
  }


  double W = 0.0;

  for (int i = 1; i <= h1.GetNbinsX(); ++i) {

    const double x  = h1.GetBinCenter (i);
    const double dx = h1.GetBinWidth  (i);

    h1.SetBinContent(i, getPhotocathodeArea() * getAngularAcceptance(x) * 1.0e4);

    W += getPhotocathodeArea() * getAngularAcceptance(x) * dx * 1.0e4;

    for (double y = 340; y <= 640; y += 20.0) {
      if (KM3NET::getQE(y) > 0.0) {
        pmt[y]->SetBinContent(i, KM3NET2D::getPhotocathodeArea2D(x, y) * 1.0e4 / KM3NET::getQE(y));
      }
    }
  }

  NOTICE(getNamespace() << " PMT average photo-cathode area " << FIXED(5,2) << W << " [cm^2]" << endl);

  if (pdf) {

    convertToPDF(h0, "NW");
    convertToPDF(h1, "NW");
    
    for (JManager_t::iterator i = qe .begin(); i != qe .end(); ++i) { convertToPDF(*i->second, "NW"); } 
    for (JManager_t::iterator i = pmt.begin(); i != pmt.end(); ++i) { convertToPDF(*i->second, "NW"); } 
  }

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

  out << h0 << h1;
  out << qe << pmt;

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