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 "JGizmo/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"
 
 
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   - MC 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 ( is_valid(header.coord_origin) ) {
      
      coord_origin = header.coord_origin;
 
      // using NOTICE here gives a linker error...
      cout << "JVolume1D::get_coord_origin() Coordinate origin (x, y, z) " << coord_origin << " read from header" << endl;
 
    }
    
    else if ( is_valid(header.can) ) {
 
      if (header.can.zmin < 0) coord_origin.z = -header.can.zmin;
    
      cout << "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;
  int                          debug;
  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_genie(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();
 
}