JAddHDF.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <fstream>
#include <iomanip>
 
#include "JPhysics/JPDFTable.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/JVersor3D.hh"
#include "JPhysics/JConstants.hh"
 
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
#include "JPhysics/JPDFTransformer.hh"
 
 
/**
 * \file
 *
 * Program to add multi-dimensional histograms of muon light
 * \author lquinn
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
 
  vector<string>      inputFile;
  string              outputFile;
  int                 debug;
 
  try {
 
    JParser<> zap("Program to add multi-dimensional histograms of muon light");
 
    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,
                                                JHistogramGridMap_t,
                                                JHistogramGridMap_t>::maplist> JMultiHistogram_t;
 
  typedef JPDFTransformer<3, abscissa_type>                         JFunctionTransformer_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 JFunctionTransformer_t());
  h2.transformer.reset(new JFunctionTransformer_t());
 
  bool add = false;
 
  for(vector<string>::const_iterator file_name = inputFile.begin(); file_name != inputFile.end(); ++file_name) {
 
    NOTICE("loading input from file " << *file_name << "... " << flush);
 
    try {
 
      JFileStreamReader in(file_name->c_str());
 
      if (add) {
 
        JMultiHistogram_t p0;
        JMultiHistogram_t p1;
        JMultiHistogram_t p2;
 
        p1.transformer.reset(new JFunctionTransformer_t());
        p2.transformer.reset(new JFunctionTransformer_t());
 
        for(JMultiHistogram_t* p : { &p0, &p1, &p2 }) {
          in.load(*p);
        }
 
        in.close();
 
        NOTICE("done." << endl);
 
        // Add new histogram to first histogram
 
        for(JMultiHistogram_t::super_iterator
              i0 = h0.super_begin(), i1 = h1.super_begin(), i2 = h2.super_begin(),
              j0 = p0.super_begin(), j1 = p1.super_begin(), j2 = p2.super_begin();
            i0 != h0.super_end(); ++i0, ++i1, ++i2, ++j0, ++j1, ++j2)
          {
            i0.getValue().add(j0.getValue());
            i1.getValue().add(j1.getValue());
            i2.getValue().add(j2.getValue());
          }
      }
      else {
 
        for(JMultiHistogram_t* p : { &h0, &h1, &h2 }) {
          in.load(*p);
        }
 
        in.close();
 
        NOTICE("done." << endl);
 
        add = true;
      }
    }
    catch(const JException& error) {
      FATAL(error.what() << endl);
    }
  }
 
  JFileStreamWriter out(outputFile.c_str());
 
  NOTICE("Storing, " << flush);
 
  for (const JMultiHistogram_t* p : { &h0, &h1, &h2 }) {
    out.store(*p);
  }
 
  out.close();
  NOTICE("done." << endl);
}