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/offline/Head.hh"
#include "km3net-dataformat/offline/Evt.hh"
#include "km3net-dataformat/online/JDAQTriggerMask.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 "JPhysics/JConstants.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

main()

int main	(	int	argc,
		char **	argv )

Definition at line 40 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;
  JDAQTriggerMask           trigger_mask;
  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['T'] = make_field(trigger_mask, "Trigger mask")                                             = TRIGGER_MASK_ON;
    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[%]", option.c_str(), 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);
 
      const JEvtWeight& weighter = scanner->getEvtWeighter();
 
      if (!header.is_valid(&JHead::genvol)) {
        FATAL("Mising normalisation in header." << endl);
      }
 
      double Wnorm = header.genvol.numberOfEvents;
 
      if (header.is_valid(&JHead::cut_nu))
        Wnorm *= 2*PI * (header.cut_nu.cosT.getUpperLimit() - header.cut_nu.cosT.getLowerLimit());
      else
        Wnorm *= 4*PI;
 
      while (scanner->hasNext()) {
 
        STATUS("event: " << setw(10) << scanner->getCounter() << '\r'); DEBUG(endl);
 
        const Evt* event   = scanner->next();
        const Trk& primary = get_primary(*event);
 
        if (trigger_mask.hasTriggerMask(event->trigger_mask)) {
 
          double x = logx ? log10(primary.E) : primary.E;
          double y = 0.0;
          double w = (getVolume(weighter, *event)                                               /  
                      (logx ? log(10.0) * primary.E * hm->GetBinWidth(1) : hm->GetBinWidth(1))  /
                      Wnorm);
 
          if        (option == Mass_t) {
            y = w * DENSITY_SEA_WATER * 1e-6;       // Mton
          } else if (option == Volume_t) {
            y = w * 1e-9;                           // 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);
  }
}