JEffectiveMass1D.cc File Reference

Example program to histogram neutrino effective mass for triggered events by counting 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 "JDAQ/JDAQEventIO.hh"
#include "JGeometry3D/JCylinder3D.hh"
#include "JGeometry3D/JPosition3D.hh"
#include "JAAnet/JVolume.hh"
#include "JAAnet/JHeadToolkit.hh"
#include "JAAnet/JAAnetToolkit.hh"
#include "JSupport/JSupport.hh"
#include "JSupport/JTriggeredFileScanner.hh"
#include "JSupport/JEvtWeightFileScannerSet.hh"
#include "JSupport/JMonteCarloFileSupportkit.hh"
#include "JROOT/JManager.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 by counting events inside the can volume.

All generated events should therefore be available in the input file(s).
This implies that option JSirene -N 0 should be used and option JTriggerEfficiency -O should not be used.
Alternatively, the option -C can be used to specify a cylinder as the can.
The unit of the contents of the histogram is or .

Author

mdejong, bstrandberg, bjung

Definition in file JEffectiveMass1D.cc.

Function Documentation

int main	(	int	argc,
		char **	argv
	)

Definition at line 47 of file JEffectiveMass1D.cc.

{
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;

  JMultipleFileScanner_t    inputFiles;
  string                    outputFile;
  double                    margin;
  bool                      logx;
  int                       numberOfBins;
  JCylinder3D               cylinder;
  std::string               option;
  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(margin,       "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['C'] = make_field(cylinder,     "User cylinder")                                                        = JPARSER::initialised();
    zap['O'] = make_field(option,       "Result option")                                                        = Mass_t, Volume_t;
    zap['d'] = make_field(debug)           = 2;
    
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }

  try{

    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 JVolume volume(scanner->getHeader(), logx);

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

    JManager<int, TH1D> hM (new TH1D("hM[%]",  NULL, numberOfBins, Xmin, Xmax));
    JManager<int, TH1D> hNM(new TH1D("hNM[%]", NULL, numberOfBins, Xmin, Xmax));

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

      const JHead header = scanner->getHeader();

      NOTICE("Scanning file type " << scanner->getName() << endl);
      DEBUG (header << endl);

      JCylinder3D can = cylinder;

      if (can.getVolume() == 0.0) {
        can = getCylinder(header);
      }

      can.addMargin(margin);

      NOTICE("Cylinder, including margin: " << can  << endl);

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

      size_t N = 0;
    
      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);

        if (event->mc_hits.empty()) {
          N += 1;
        }

        const double x = logx ? log10(primary.E) : primary.E;

        if (can.is_inside(vertex)) {
          hNM[primary.type]->Fill(x);
        }
      }
      STATUS(endl);

      if (N == 0) {
        ERROR("No events with zero hits -> use application JEffectiveMassOffline1D." << endl);
      }

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

      for (JTriggeredFileScanner<> in(scanner->getFilelist()); in.hasNext(); ) {

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

        const Evt*   event   = in.next();
        const Trk&   primary = get_primary(*event);

        double x = logx ? log10(primary.E) : primary.E;
        double y = 0.0;

        if        (option == Mass_t) {
          y = can .getVolume() * DENSITY_SEA_WATER * 1e-6;    // Mton
        } else if (option == Volume_t) {
          y = can .getVolume() * 1.0e-9;                      // km^3
        }

        hM[primary.type]->Fill(x, y);    
      }
      STATUS(endl);
    }

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

    for (int pdg : get_keys(hM)) {
      hM[pdg]->Divide(hNM[pdg]);
    }

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

    out << hM;
  
    out.Close();
  }
  catch(const JException& error) {
    FATAL(error << endl);
  }
}