JEffectiveMass1D.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 "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"


/**
 * \file
 *
 * Example program to histogram neutrino effective mass for triggered events inside the can volume.\n
 * For genie/gSeaGen events the effective mass for triggered events inside the instrumented volume is also histogrammed.\n
 * The unit of the contents of the histogram is Mton.
 * \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;
  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();

}