JVolume1D.cc

Go to the documentation of this file.

#include <string>
#include <iostream>

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

#include "km3net-dataformat/offline/Head.hh"
#include "km3net-dataformat/offline/Evt.hh"
#include "km3net-dataformat/offline/Vec.hh"
#include "JDAQ/JDAQEventIO.hh"
#include "JAAnet/JAAnetToolkit.hh"
#include "JAAnet/JVolume.hh"

#include "JSupport/JTriggeredFileScanner.hh"
#include "JSupport/JMonteCarloFileSupportkit.hh"
#include "JSupport/JSupport.hh"

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

/**
 * \file
 *
 * Example program to histogram neutrino effective volume for triggered events.\n
 * For genie/gSeaGen events, use application JEffectiveMassOnline1D or JEffectiveMassOffline1D.\n
 * The unit of the contents of the histogram is \f$km^{3}\f$.
 * \author mdejong, bstrandberg
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;
  
  JTriggeredFileScanner<> inputFile;
  string   outputFile;
  JLimit&  numberOfEvents = inputFile.getLimit();
  bool     logx;
  int      numberOfBins;
  int      debug;
  
  try {

    JParser<> zap("Example program to histogram neutrino effective volume for triggered events."\
                  "\nFor genie/gSeaGen events, use application JEffectiveMass1D.");

    zap['f'] = make_field(inputFile);
    zap['o'] = make_field(outputFile)     = "volume.root";
    zap['n'] = make_field(numberOfEvents) = JLimit::max();
    zap['X'] = make_field(logx,         "Use logarithm of energy");
    zap['N'] = make_field(numberOfBins, "Number of bins in the energy range of the MC simulation") = 10;
    zap['d'] = make_field(debug)          = 1;
    
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }

  Head head;

  try {
    head = getHeader(inputFile);
  }
  catch(const JException& error) {
    FATAL(error);
  }
  
  const JHead buffer(head);

  if (!is_genhen(buffer)) {
    FATAL("Use JEffectiveMass1D" << endl);
  }
  
  const JVolume volume(head, logx);

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

  TH1D hV("hV", "effective volume in km^3", numberOfBins, volume.getXmin(), volume.getXmax());

  hV.Sumw2();

  while (inputFile.hasNext()) {

    STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);

    JTriggeredFileScanner<>::multi_pointer_type ps = inputFile.next();

    const Evt* event = ps;
    
    if (has_neutrino(*event)) {

      const Trk&   neutrino = get_neutrino(*event);
      const double E        = neutrino.E;

      hV.Fill(volume.getX(E), volume.getW(hV.GetXaxis(), E));
    }
    else {
      WARNING("JVolume1D: cannot find neutrino in triggered event " << inputFile.getCounter() );
    }
  }
  STATUS(endl);
  
  out.Write();
  out.Close();
}