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 "JGeometry3D/JCylinder3D.hh"
#include "JGeometry3D/JPosition3D.hh"

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

#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"

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


int   debug;

namespace {

  using namespace JPP;

  /**
   * Test event from genie/gSeaGen simulation.
   *
   * The function checks that the vertex of the neutrino event is inside the can volume.
   *
   * \param  cylinder     cylinder
   * \param  event        event
   * \return              true if accepted; else false
   */
  inline bool test_event(const JCylinder3D& cylinder, const Evt& event)
  {
    using namespace std;
    using namespace JPP;

    Bool_t accepted = true;
    
    if ( has_neutrino(event) ) {      
      accepted = accepted && cylinder.is_inside( getPosition( get_neutrino(event) ) );
    }
    else {

      // future/further logic for atm_muon events here, if any
      
      accepted = false;
    }

    return accepted;
  }

  /**
   * Function to infer the coordinate origin from the header.
   *
   * This bulky logic exists because historically gSeaGen and Jpp have used different coordinate systems. For gSeaGen files the coordinate origin can be syncronised to Jpp, which is performed here.
   *
   * \param header     Monte Carlo simulation file header
   * \return           Coordinate origin in the simulation
   */
  inline Vec get_coord_origin(const JHead &header)
  {
    using namespace std;
    
    Vec coord_origin(0, 0, 0);

    if        (header.is_valid(&JHead::coord_origin)) {
      
      coord_origin = header.coord_origin;

      NOTICE("JVolume1D::get_coord_origin() Coordinate origin (x, y, z) " << coord_origin << " read from header" << endl);

    } else if (header.is_valid(&JHead::can) ) {

      if (header.can.zmin < 0) coord_origin.z = -header.can.zmin;
    
      NOTICE("JVolume1D::get_coord_origin() Coordinate origin (x, y, z) " << coord_origin << " calculated from can dimensions (zmin, zmax, r) " << header.can.zmin << " " << header.can.zmax <<  " " << header.can.r << endl);

    } else {

      FATAL("JVolume1D::get_coord_origin() Error in determining the coordinate origin.");
    }

    return coord_origin;    
  }
}

/**
 * \file
 *
 * Example program to histogram neutrino effective volume for triggered events.  For genie/gSeaGen events a histogram depicting the fraction of events that triggered the detector inside the instrumented volume is also shown.
 * 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;

  //------------------------------------------------------------------------------
  // parse command line arguments, check that header can be found in the input file
  //------------------------------------------------------------------------------
  
  JTriggeredFileScanner<>      inputFile;
  string                       outputFile;
  JLimit&                      numberOfEvents = inputFile.getLimit();
  double                       wall;
  bool                         logx;
  int                          numberOfBins;
  string                       detectorFile;

  
  try {

    JParser<> zap("Example program to histogram neutrino effective volume for triggered events. For genie/gSeaGen events a histogram depicting the fraction of events that triggered the detector inside the instrumented volume is also shown.");

    zap['f'] = make_field(inputFile);
    zap['o'] = make_field(outputFile) = "volume.root";
    zap['n'] = make_field(numberOfEvents) = JLimit::max();
    zap['R'] = make_field(wall, "Addition margin around the volume in which the considered events must reside") = 0.0;
    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['a'] = make_field(detectorFile, "Detector file: if not provided, trigger fraction is not calculated") = "";
    
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }

  //-------------------------------------------------------------------------
  // deal with the header and init the detector
  //-------------------------------------------------------------------------

  Head head;

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

  if (detectorFile != "") {

    try {
      load(detectorFile, detector);
    }
    catch(const JException& error) {
      FATAL(error);
    }
  }
    
  //------------------------------------------------------------------------------
  // create the can volume from the dimensions stored in the header and instrumented volume from detector file
  //------------------------------------------------------------------------------
  
  const JHead buffer(head);

  const bool genie = is_gseagen(buffer);
  
  // the can volume comes from the generator header and the coordinate system is consistent by construction
  const JVolume volume(head, logx);
  JCylinder3D canvol = get<JCylinder3D>(buffer);
  canvol.addMargin(wall);

  // as the detector is constructed from Jpp detector file, the coordinate system needs to be synced
  // to the one used in the event generator
  detector -= JAANET::getPosition( get_coord_origin(buffer) );
  JCylinder3D instvol( detector.begin(), detector.end() );
  instvol.addMargin(wall);

  NOTICE("JVolume1D: Instrumented volume dimensions (zmin, zmax, r): " << instvol.getZmin() << " " << instvol.getZmax() << " " << instvol.getRadius() << endl );
  
  //------------------------------------------------------------------------------
  // initialise output
  //------------------------------------------------------------------------------

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

  TH1D hV("hV", "effective volume in km^3"           , numberOfBins, volume.getXmin(), volume.getXmax());
  TH1D hF("hF", "n_trig/n_gen in instrumented volume", numberOfBins, volume.getXmin(), volume.getXmax());
  hV.Sumw2();
  hF.Sumw2();

  //------------------------------------------------------------------------------
  // loop over triggered events in the input file
  //------------------------------------------------------------------------------

  NOTICE("JVolume1D: Scanning triggered events." << endl);
  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;

      // genie/gSeaGen events filled with weight 1, histograms normalised after the second loop below
      if ( genie ) {
      
        // events are filled to hV if they are inside the can
        hV.Fill(volume.getX(E), test_event(canvol, *event) ? 1.0 : 0.0);

        // events are filled to hF if they are inside the instrumented volume
        hF.Fill(volume.getX(E), test_event(instvol, *event) ? 1.0 : 0.0);

      }
      // for other simulated events effective volume calculated from event weights, normalisation not required
      else {
        hV.Fill(volume.getX(E), volume.getW(hV.GetXaxis(), E));
      }
      
    }
    else {
      WARNING("JVolume1D: cannot find neutrino in triggered event " << inputFile.getCounter() );
    }
    
  } // end loop over triggered events
  STATUS(endl);
  
  //------------------------------------------------------------------------------
  // loop over generated events in the input file; FATAL if generator events have not been stored 
  //------------------------------------------------------------------------------

  if ( genie ) {
  
    // this histogram counts the generated events inside the can volume
    TH1D* hNV = (TH1D*) hV.Clone("hNV");
    hNV->Reset();

    // this histogram counts the generated events inside the instrumented volume
    TH1D* hNF = (TH1D*) hF.Clone("hNF");
    hNF->Reset();

    NOTICE("JVolume1D: Scanning generated events." << endl);

    JMultipleFileScanner<Evt> in(inputFile);

    while ( in.hasNext() ) {

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

      const Evt* event = in.next();
      
      if (has_neutrino(*event)) {

        const Trk&   neutrino = get_neutrino(*event);
        const double E        = neutrino.E;
        
        hNV->Fill(volume.getX(E), test_event(canvol, *event) ? 1.0 : 0.0);
        hNF->Fill(volume.getX(E), test_event(instvol, *event) ? 1.0 : 0.0);

      }
      else {
        WARNING("JVolume1D: cannot find neutrino in generated event " << inputFile.getCounter() );
      }
      
    }
    STATUS(endl);

    if ( in.getCounter() == 0 ) {
      FATAL("JVolume1D: generated events not stored in the input file, JTriggerEfficiency should be run without option -O");
    }

    //------------------------------------------------------------------------------
    // convert 'generated' and 'triggered' histograms to effective volume and trigger fraction
    //------------------------------------------------------------------------------
  
    hV.Divide(hNV);
    hV.Scale(canvol.getVolume()*1e-9);    // km^3
    hF.Divide(hNF);

    delete hNV;
    delete hNF;  
  
  }
  
  out.Write();
  out.Close();

}