JMakeHDE.cc

Go to the documentation of this file.

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

#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"

#include "JPhysics/JPDFTransformer.hh"


/**
 * \file
 *
 * Program to convert 5-dimensional histograms of shower light to multi-dimensional PDFs.
 * \author jseneca
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;

  string         inputFile;
  string         outputFile;
  int            debug;
  try {

    JParser<> zap("Program to convert multi-dimensional histograms of shower light to multi-dimensional PDFs.");
    
    zap['f'] = make_field(inputFile);
    zap['o'] = make_field(outputFile);
    zap['d'] = make_field(debug)               = 1;
    
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }


  typedef JHistogram1D_t::abscissa_type                             abscissa_type;

  typedef JTransformableMultiHistogram<JHistogram1D_t, 
                              JMAPLIST<JHistogramMap_t,
                                       JHistogramMap_t,
                                       JHistogramMap_t,
                                   JHistogramGridMap_t,
                                   JHistogramGridMap_t>::maplist>   JMultiHistogram_t;

  typedef JPDFTransformer<5, abscissa_type>                         JFunction5DTransformer_t;

  JMultiHistogram_t h0;   // occurrence rate of PMT (used for normalisation)
  JMultiHistogram_t h1;   // light from cascade

  h1.transformer.reset(new JFunction5DTransformer_t());

  // input

  NOTICE("Loading in input, " << flush);

  JFileStreamReader in(inputFile.c_str());

  for (JMultiHistogram_t* p : { &h0, &h1 }) {
    in.load(*p);
  }

  in.close();

  NOTICE("done." << endl);

  // 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(); i1 != h1.super_end(); ++i0, ++i1) {

    const double W = i0.getValue().getIntegral();

    if (W != 0.0) {
      
      const JDirection3D u(JAngle3D(i0->second->second->first, i0->second->second->second->first));
      const JDirection3D v(getSinThetaC(), 0.0, getCosThetaC());
      
      int number_of_bins = (int) (2 + u.getDot(v));
      
      i1.getValue().rebin(JHistogram1D_t::JRebin(number_of_bins));
      i1.getValue().div(W);
    }    
  }

  // output

  struct tuple {
    const JPDFType_t         type;
    const JMultiHistogram_t& value;
  };
  
  typedef JPDFTable<JPolint1Function1D_t, 
                    JMAPLIST<JPolint1FunctionalMap,
                             JPolint1FunctionalMap,
                             JPolint1FunctionalMap,
                             JPolint1FunctionalGridMap,
                             JPolint1FunctionalGridMap>::maplist> JPDF_t;

  try {
    NOTICE("Storing output to file " << outputFile << "... " << flush);

    const JPDF_t pdf(h1);

    pdf.store(outputFile.c_str());

    NOTICE("OK" << endl);
  }
  catch(const JException& error) {
    FATAL(error.what() << endl);
  }
  NOTICE("JMakeHDE done." << endl);
}