JEffectiveMassOffline.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 "JSirene/JVisibleEnergyToolkit.hh"
#include "JROOT/JManager.hh"
#include "JTools/JAbstractHistogram.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 . It produces histograms based on true neutrino energy (hm0), total visible energy (hm1), leptonic visible energy (hm2) and hadronic visible energy (hm3).

Author

bjung, mdejong

Definition in file JEffectiveMassOffline.cc.

Function Documentation

main()

int main	(	int	argc,
		char **	argv )

Definition at line 74 of file JEffectiveMassOffline.cc.

{
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;
 
  typedef JAbstractHistogram<double>   JAbstractHistogram_t;
  
 
  JMultipleFileScanner_t    inputFiles;
  string                    outputFile;
  
  JAbstractHistogram_t      X;
  bool                      logx;
  double                    margin;
  double                    epsilon;
  JCylinder3D               cylinder;
 
  std::string               option;
  
  JDAQTriggerMask           trigger_mask;
  
  int                       debug;
 
  
  try {
 
    JParser<> zap("Example program to histogram neutrino effective mass for triggered events using offline files.");
 
    zap['f'] = make_field(inputFiles);
    zap['o'] = make_field(outputFile)      = "Meff.root";
    zap['x'] = make_field(X,            "Histogram binning")
      = JAbstractHistogram_t(100, 1.0, 100.0);
    zap['X'] = make_field(logx,         "Use logarithm of energy");
    zap['C'] = make_field(cylinder,     "Cylinder to compute visible energy depositions")
      = JPARSER::initialised();
    zap['m'] = make_field(margin,       "Extra margin to consider visible energy depositions")
      = 0.0;
    zap['e'] = make_field(epsilon,      "Minimal amount of visible energy assigned to events") = 1e-8;
    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);
  }
 
 
  const JEvtWeightFileScannerSet<> scanners(inputFiles);  
  
  const JMultiHead multiheader = getMultiHeader(inputFiles);
 
  const double Xmin    = X.getLowerLimit();
  const double Xmax    = X.getUpperLimit();
  const int    Xnbins  = X.getNumberOfBins();
 
  JManager<int, TH1D> hm0(new TH1D("hm0[%]", option.c_str(), Xnbins, Xmin, Xmax));
  JManager<int, TH1D> hm1(new TH1D("hm1[%]", option.c_str(), Xnbins, Xmin, Xmax));
  JManager<int, TH1D> hm2(new TH1D("hm2[%]", option.c_str(), Xnbins, Xmin, Xmax));
  JManager<int, TH1D> hm3(new TH1D("hm3[%]", option.c_str(), Xnbins, Xmin, Xmax));
 
  for (JEvtWeightFileScannerSet<>::const_iterator scanner = scanners.begin(); scanner != scanners.end(); ++scanner) {
 
    JHead header = scanner->getHeader();
 
    JEvtWeightHelper weighter(scanner->getEvtWeighter());
 
    NOTICE("Scanning file type " << scanner->getName() << endl);
    DEBUG (header << endl);
      
    for (JMultipleFileScanner<Evt> in(scanner->getFilelist()); in.hasNext(); ) {
 
      STATUS("event: " << setw(10) << in.getCounter() << '\r'); DEBUG(endl);
 
      Evt* event   = in.next();
      const Trk& primary = get_primary(*event);
 
 
      if (trigger_mask.hasTriggerMask(event->trigger_mask)) {
        for (vector<Trk>::iterator track = event->mc_trks.begin(); track != event->mc_trks.end(); ++track) {
          track->clearusr(); // Prevent visible energy from being calculated based on `mc_usr_keys::energy_lost_in_can`   
        }
        double Wnorm = 0.0;
          
        if        (header.is_valid(&JHead::genvol)) {
            
          Wnorm = getWnorm(header);
            
        } else if (uuid_is_null(event->header_uuid) == 0) {
 
          const JUUID mergedUUID = multiheader.getHeaderUUID(*event);
          
          MultiHead::const_iterator i = multiheader.find(mergedUUID.uuid, true);
 
          if (i == multiheader.cend()) {
            FATAL("Could not find header corresponding to UUID: " << JUUID(event->header_uuid) << endl);
          }
 
          const JHead header1(*i);
          DEBUG(header1 << endl);
 
          Wnorm = getWnorm(header1);
 
          weighter.configure(getEventWeighter(header1));
            
        } else {
            
          FATAL("Mising normalisation in event/header." << endl);
        }
 
        JCylinder3D can = cylinder;
 
        if (can.getVolume() == 0.0) {
          can = getCylinder(header);
        }
 
        can.addMargin(margin);
        
        double Evis    = getVisibleEnergy(*event, can);
        double EvisLL  = (has_leading_lepton(*event) ?
                          getVisibleEnergyLeadingLepton(*event, can) :
                          0.0);
        double EvisHad = Evis - EvisLL;
 
        if (Evis    == 0) Evis    = epsilon;
        if (EvisLL  == 0) EvisLL  = epsilon;
        if (EvisHad == 0) EvisHad = epsilon;
        
        const double x0 = logx ? log10(primary.E) : primary.E;
        const double x1 = logx ? log10(Evis)      : Evis;
        const double x2 = logx ? log10(EvisLL)    : EvisLL;
        const double x3 = logx ? log10(EvisHad)   : EvisHad;
 
        double y = 0.0;
        double w = (getVolume(weighter.getEvtWeighter(), *event)                              /
                    (logx ? log(10.0) * primary.E * hm0->GetBinWidth(1) : hm0->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
        }
 
        const int interactionID0 = getInteractionID(*event, false, false);
        const int interactionID1 = getInteractionID(*event, true,  false);
 
        hm0[interactionID0]->Fill(x0, y);
        hm1[interactionID0]->Fill(x1, y);
        hm2[interactionID0]->Fill(x2, y);
        hm3[interactionID0]->Fill(x3, y);
 
        hm0[interactionID1]->Fill(x0, y);
        hm2[interactionID1]->Fill(x2, y);
        hm1[interactionID1]->Fill(x1, y);
        hm3[interactionID1]->Fill(x3, y);
 
        if (abs(primary.type) == TRACK_TYPE_NUTAU) {
 
          const int interactionID2 = getInteractionID(*event, false, true);
          const int interactionID3 = getInteractionID(*event, true,  true);
 
          hm0[interactionID2]->Fill(x0, y);
          hm1[interactionID2]->Fill(x1, y);
          hm2[interactionID2]->Fill(x2, y);
          hm3[interactionID2]->Fill(x3, y);
 
          hm0[interactionID3]->Fill(x0, y);
          hm2[interactionID3]->Fill(x2, y);
          hm1[interactionID3]->Fill(x1, y);
          hm3[interactionID3]->Fill(x3, y);
        }
      }
    }
    STATUS(endl);
  }
  
  TFile out(outputFile.c_str(), "recreate");
 
  out << hm0 << hm1 << hm2 << hm3;
 
  out.Close();
}