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