JEffectiveMassOffline1D.cc File Reference

Example program to histogram neutrino effective mass for triggered events for offline reconstruction files, using the event weights only. More...

#include <string>
#include <iostream>
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "km3net-dataformat/definitions/weightlist.hh"
#include "km3net-dataformat/definitions/w2list_gseagen.hh"
#include "km3net-dataformat/offline/Head.hh"
#include "km3net-dataformat/offline/Evt.hh"
#include "km3net-dataformat/offline/Vec.hh"
#include "JAAnet/JVolume.hh"
#include "JAAnet/JHeadToolkit.hh"
#include "JAAnet/JAAnetToolkit.hh"
#include "JSupport/JSupport.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 for offline reconstruction files, using the event weights only.

The unit of the contents of the histogram is or .

Author

bjung, mdejong

Definition in file JEffectiveMassOffline1D.cc.

Function Documentation

int main	(	int	argc,
		char **	argv
	)

Definition at line 41 of file JEffectiveMassOffline1D.cc.

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

  JMultipleFileScanner_t    inputFiles;
  string                    outputFile;
  bool                      logx;
  int                       numberOfBins;
  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['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['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 JHead header = scanner->getHeader();

      if (!is_gseagen(header)) {
        ERROR("Not a gSeaGen file; you may want to use JVolume1D" << endl);
      }
    
      const JVolume volume(header, logx);

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

    JManager<int, TH1D> hm(new TH1D("hm[%]", option.c_str(), numberOfBins, Xmin, Xmax));

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

      const JHead header = scanner->getHeader();

      const double Wnorm = header.genvol.numberOfEvents * log(10) * hm->GetBinWidth(1) * 4*M_PI;

      NOTICE("Scanning file type " << distance(scanners.begin(), scanner) << endl);
      DEBUG (header << 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 double weight  = (event->w[WEIGHTLIST_DIFFERENTIAL_EVENT_RATE] *
                                event->w2list[W2LIST_GSEAGEN_WATER_INT_LEN]  /
                                event->w2list[W2LIST_GSEAGEN_P_EARTH]        /
                                primary.E                                    /
                                Wnorm);                // cm^3

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

        double y = 0.0;

        if        (option == Mass_t) {
          y = weight * DENSITY_SEA_WATER * 1e-6;       // Mton
        } else if (option == Volume_t) {
          y = weight * 1e-15;                          // km^3
        }

        hm[primary.type]->Fill(x, y);
      }
      STATUS(endl);
    }
  
    TFile out(outputFile.c_str(), "recreate");

    out << hm;
  
    out.Close();
  }
  catch(const JException& error) {
    FATAL(error << endl << "You may want to use \"offline\" formatted file.");
  }
}