JVolume1D.cc File Reference

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

Go to the source code of this file.

Functions
int	main (int argc, char **argv)

Detailed Description

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 .

Author

mdejong, bstrandberg

Definition in file JVolume1D.cc.

Function Documentation

main()

int main	(	int	argc,
		char **	argv
	)

Definition at line 108 of file JVolume1D.cc.

{
  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();
 
}