JEffectiveMassOffline.cc

Go to the documentation of this file.

#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"
 
 
namespace {
}
 
 
namespace {
 
  const char* const Mass_t    =  "Mass";   //!< Effective mass   option
  const char* const Volume_t  =  "Volume"; //!< Effecitve volume option
 
 
  /**
   * Retrieve weight normalisation from given header.
   *
   * \param  header             header
   */
  double getWnorm(const JAANET::JHead& header)
  {
    using namespace JPP;
    
    if (!header.is_valid(&JHead::genvol)) {
      THROW(JNoValue, "getWnorm(): Mising normalisation in header.");
    }
    
    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;
  
    return Wnorm;
  }
}
 
/**
 * \file
 *
 * Example program to histogram neutrino effective mass for triggered events for offline reconstruction files, using the event weights only.\n
 * The unit of the contents of the histogram is \f$Mton\f$ or \f$km^{3}\f$. 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
 */
int main(int argc, char **argv)
{
  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();
}