JEffectiveMass1D.cc File Reference

Example program to histogram neutrino effective mass for triggered events inside the can volume. More...

#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 "JGeometry3D/JCylinder3D.hh"
#include "JGeometry3D/JPosition3D.hh"
#include "JAAnet/JAAnetToolkit.hh"
#include "JAAnet/JHeadToolkit.hh"
#include "JAAnet/JVolume.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 mass for triggered events inside the can volume.

For genie/gSeaGen events the effective mass for triggered events inside the instrumented volume is also histogrammed.
The unit of the contents of the histogram is Mton.

Author

mdejong, bstrandberg

Definition in file JEffectiveMass1D.cc.

Function Documentation

int main	(	int	argc,
		char **	argv
	)

Definition at line 40 of file JEffectiveMass1D.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;
  double                       wall;
  bool                         logx;
  int                          numberOfBins;
  string                       detectorFile;
  int                          debug;

  
  try {

    JParser<> zap("Example program to histogram neutrino effective mass for triggered events.");

    zap['f'] = make_field(inputFile);
    zap['o'] = make_field(outputFile) = "Meff.root";
    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<Vec>(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 hM("hM", "effective mass (can) [kg]"         , numberOfBins, volume.getXmin(), volume.getXmax());
  TH1D hm("hm", "effective mass (instrumented) [kg]", numberOfBins, volume.getXmin(), volume.getXmax());
  hM.Sumw2();
  hm.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 JPosition3D& vertex   = getPosition(neutrino);
      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 hM if they are inside the can
        hM.Fill(volume.getX(E), canvol.is_inside(vertex) ? 1.0 : 0.0);

        // events are filled to hm if they are inside the instrumented volume
        hm.Fill(volume.getX(E), instvol.is_inside(vertex) ? 1.0 : 0.0);

      }
      // for other simulated events effective volume calculated from event weights, normalisation not required
      else {
        hM.Fill(volume.getX(E), volume.getW(hM.GetXaxis(), E));
      }
      
    }
    else {
      WARNING("JEffectiveMass1D: 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* hNM = (TH1D*) hM.Clone("hNM");
    hNM->Reset();

    // this histogram counts the generated events inside the instrumented volume
    TH1D* hNm = (TH1D*) hm.Clone("hNm");
    hNm->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 JPosition3D& vertex   = getPosition(neutrino);
        const double       E        = neutrino.E;
        
        hNM->Fill(volume.getX(E), canvol.is_inside(vertex) ? 1.0 : 0.0);
        hNm->Fill(volume.getX(E), instvol.is_inside(vertex) ? 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 mass
    //------------------------------------------------------------------------------
  
    hM.Divide(hNM);
    hM.Scale(canvol.getVolume()  * DENSITY_SEA_WATER * 1e-6);   // Mton
    hm.Divide(hNm);
    hm.Scale(instvol.getVolume() * DENSITY_SEA_WATER * 1e-6);   // Mton

    delete hNM;
    delete hNm;  
  }
  
  out.Write();
  out.Close();

}