JEffectiveMass1D.cc File Reference

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

#include <string>
#include <iostream>
#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/JEvtWeightFileScannerSet.hh"
#include "JSupport/JMonteCarloFileSupportkit.hh"
#include "JSupport/JSupport.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JPhysics/JConstants.hh"
#include "JPhysics/JGeane.hh"
#include "JReconstruction/JEvt.hh"
#include "JReconstruction/JEvtToolkit.hh"
#include "JROOT/JManager.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"

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 84 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
  //------------------------------------------------------------------------------
  
  JMultipleFileScanner_t       inputFiles;
  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(inputFiles);
    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);
  }

  //-------------------------------------------------------------------------
  // Initialise the detector
  //-------------------------------------------------------------------------
  
  JDetector detector;

  if (detectorFile != "") {

    try {
      load(detectorFile, detector);
    }
    catch(const JException& error) {
      FATAL(error);
    }
  }


  //-------------------------------------------------------------------------
  // Retrieve auxiliary header info
  //-------------------------------------------------------------------------

  JCylinder3D canvol;
  JCylinder3D instvol;

  vector<JVolume> volumes;

  double Xmin = numeric_limits<double>::max();
  double Xmax = numeric_limits<double>::lowest();

  JEvtWeightFileScannerSet<> scanners(inputFiles);

  for (JEvtWeightFileScannerSet<>::const_iterator scanner = scanners.cbegin(); scanner != scanners.cend(); ++scanner) {
    
    const JHead& header = scanner->getHeader();
    
    const JVolume volume(header, logx);

    if (volume.getXmin() < Xmin) { Xmin = volume.getXmin(); }
    if (volume.getXmax() > Xmax) { Xmax = volume.getXmax(); }

    volumes.push_back(volume);

    // the can volume comes from the generator header and the coordinate system is consistent by construction
    JCylinder3D can = get<JCylinder3D>(header);
    canvol.addMargin(wall);

    if (can.getVolume() > canvol.getVolume()) { canvol = can; }

    // as the detector is constructed from Jpp detector file, the coordinate system needs to be synced
    // to the one used in the event generator
    JPosition3D center = getPosition(get<Vec>(header));
    detector -= center;

    JCylinder3D inst( detector.begin(), detector.end() );
    instvol.addMargin(wall);

    if (inst.getVolume() > instvol.getVolume()) { instvol = inst; }

    detector += center;
  }

  
  //------------------------------------------------------------------------------
  // initialise output
  //------------------------------------------------------------------------------

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

  TH1::SetDefaultSumw2();

  JManager<int, TH1D> hM(new TH1D("hM[%]", "effective mass (can) [kg]"         , numberOfBins, Xmin, Xmax));
  JManager<int, TH1D> hm(new TH1D("hm[%]", "effective mass (instrumented) [kg]", numberOfBins, Xmin, Xmax));
  JManager<int, TH1D> hNM(new TH1D("hNM[%]", "normalization effective mass (can) [kg]"         , numberOfBins, Xmin, Xmax));
  JManager<int, TH1D> hNm(new TH1D("hNm[%]", "normalization effective mass (instrumented) [kg]", numberOfBins, Xmin, Xmax));


  for (JEvtWeightFileScannerSet<>::iterator scanner = scanners.begin(); scanner != scanners.end(); ++scanner) {

    const size_t scannerIndex = distance(scanners.begin(), scanner);

    const JVolume& volume = volumes[scannerIndex];

    NOTICE("JEffectiveMass1D: Instrumented volume dimensions (zmin, zmax, r): " << instvol.getZmin() << " " << instvol.getZmax() << " " << instvol.getRadius() << endl );

    //------------------------------------------------------------------------------
    // loop over generated events in the input files
    //------------------------------------------------------------------------------

    NOTICE("Scanning generated events for file type " << scannerIndex << endl);
    
    while (scanner->hasNext()) {

      STATUS("event: " << setw(10) << scanner->getCounter() << '\r'); DEBUG(endl);

      const Evt* event = scanner->next();

      const Trk& primary = get_primary(*event);

      const JVertex3D& vertex  = getVertex(*event);

      const bool isValid1 = (canvol.is_inside(vertex)  || hasIntersectingMuon(*event, canvol));
      const bool isValid2 = (instvol.is_inside(vertex) || hasIntersectingMuon(*event, instvol));

      // events are filled to hM if they are inside the can
      hNM[primary.type]->Fill(volume.getX(primary.E), isValid1 ? scanner->getWeight(*event) : 0.0);

      // events are filled to hm if they are inside the instrumented volume
      hNm[primary.type]->Fill(volume.getX(primary.E), isValid2 ? scanner->getWeight(*event) : 0.0);

    } // end loop over generated events
    STATUS(endl);

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

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

    NOTICE("Scanning triggered events for file type " << scannerIndex << endl);
    
    JEvtWeightFileScanner<JTriggeredFileScanner<> > in(scanner->getFilelist());
    
    while (in.hasNext()) {

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

      JTriggeredFileScanner<>::multi_pointer_type mp = in.next();

      const Evt* event = mp;

      const Trk& primary = get_primary(*event);

      const JVertex3D& vertex  = getVertex(*event);

      const bool isValid1 = (canvol.is_inside(vertex)  || hasIntersectingMuon(*event, canvol));
      const bool isValid2 = (instvol.is_inside(vertex) || hasIntersectingMuon(*event, instvol));
        
      hM[primary.type]->Fill(volume.getX(primary.E), isValid1 ? in.getWeight(*event) : 0.0);
      hm[primary.type]->Fill(volume.getX(primary.E), isValid2 ? in.getWeight(*event) : 0.0);

    } // end loop over generated events
    STATUS(endl);
  }
  

  //------------------------------------------------------------------------------
  // convert 'generated' and 'triggered' histograms to effective mass
  //------------------------------------------------------------------------------

  vector<int> pdgIDs = get_keys(hM);

  for (vector<int>::const_iterator i = pdgIDs.begin(); i != pdgIDs.end(); ++i) {

    const int pdgID = *i;

    hM[pdgID]->Divide(hNM[pdgID]);
    hM[pdgID]->Scale(canvol.getVolume()  * DENSITY_SEA_WATER * 1e-6); // Mton
    hm[pdgID]->Divide(hNm[pdgID]);
    hm[pdgID]->Scale(instvol.getVolume() * DENSITY_SEA_WATER * 1e-6); // Mton
  }

  out << hM << hm;
  
  out.Close();
}