JMakeHDF.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <fstream>
#include <iomanip>
 
#include "JPhysics/JPDFTransformer.hh"
#include "JPhysics/JPDFTable.hh"
#include "JPhysics/JPDFTypes.hh"
 
#include "JIO/JFileStreamIO.hh"
 
#include "JTools/JHistogram1D_t.hh"
#include "JTools/JHistogramMap_t.hh"
#include "JTools/JTransformableMultiHistogram.hh"
#include "JTools/JFunction1D_t.hh"
#include "JTools/JFunctionalMap_t.hh"
 
#include "JGeometry3D/JAngle3D.hh"
#include "JGeometry3D/JDirection3D.hh"
#include "JPhysics/JConstants.hh"
 
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
 
/**
 * \file
 *
 * Program to convert multi-dimensional histograms of muon light to multi-dimensional PDFs.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
 
  string      inputFile;
  string      fileDescriptor;
  int         debug;
 
  try {
 
    JParser<> zap("Program to convert multi-dimensional histograms of muon light to multi-dimensional PDFs.");
 
    zap['f'] = make_field(inputFile);
    zap['o'] = make_field(fileDescriptor)      = "J%p.dat";
    zap['d'] = make_field(debug)               = 1;
 
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }
 
 
  if (!hasWildCard(fileDescriptor)) {
    FATAL("error file descriptor " << fileDescriptor << endl);
  }
 
  typedef JHistogram1D_t::abscissa_type                             abscissa_type;
 
  typedef JTransformableMultiHistogram<JHistogram1D_t, 
                                       JMAPLIST<JHistogramMap_t,
                                                JHistogramGridMap_t,
                                                JHistogramGridMap_t>::maplist> JMultiHistogram_t;
 
  typedef JPDFTransformer<3, abscissa_type>                         JFunction3DTransformer_t;
 
  JMultiHistogram_t h0;   // occurrence rate of PMT (used for normalisation)
  JMultiHistogram_t h1;   // light from muon
  JMultiHistogram_t h2;   // light from EM showers
 
  h1.transformer.reset(new JFunction3DTransformer_t());
  h2.transformer.reset(new JFunction3DTransformer_t());
 
 
  // input
 
  JFileStreamReader in(inputFile.c_str());
 
  for (JMultiHistogram_t* p : { &h0, &h1, &h2 }) {
    in.load(*p);
  }
 
  in.close();
 
 
  // rebin histograms according angle of incidence on PMT and normalise contents to PMT occurrence rate
 
  for (JMultiHistogram_t::super_iterator
         i0 = h0.super_begin(),
         i1 = h1.super_begin(),
         i2 = h2.super_begin(); i0 != h0.super_end(); ++i0, ++i1, ++i2) {
 
    const double W = i0.getValue().getIntegral();
 
    if (W != 0.0) {
 
      const JDirection3D u(JAngle3D(i0->second->first, i0->second->second->first));
      const JDirection3D v(getSinThetaC(), 0.0, getCosThetaC());
 
      int number_of_bins = (int) (2 + u.getDot(v));
 
      for (JHistogram1D_t* p : { &i1.getValue(), &i2.getValue() }) {
        p->rebin(JHistogram1D_t::JRebin(number_of_bins));
        p->div(W);
      }
    }
  }
 
 
  // output
 
  struct tuple {
    const JPDFType_t         type;
    const JMultiHistogram_t& value;
  };
 
  tuple ntuple[] = {                       // map histogram to PDF type
    { SCATTERED_LIGHT_FROM_MUON,      h1 },
    { SCATTERED_LIGHT_FROM_EMSHOWERS, h2 }
  };
 
 
  typedef JPDFTable<JPolint1Function1D_t, 
                    JMAPLIST<JPolint1FunctionalMap,
                             JPolint1FunctionalGridMap,
                             JPolint1FunctionalGridMap>::maplist> JPDF_t;
  
  for (int i = 0; i != sizeof(ntuple)/sizeof(ntuple[0]); ++i) {
 
    const string file_name = getFilename(fileDescriptor, ntuple[i].type);
 
    try {
 
      NOTICE("storing output to file " << file_name << "... " << flush);
 
      const JPDF_t pdf(ntuple[i].value);
 
      pdf.store(file_name.c_str());
 
      NOTICE("OK" << endl);
    }
    catch(const JException& error) {
      FATAL(error.what() << endl);
    }
  }
}