JEvtReweight.cc

Go to the documentation of this file.

#include <uuid/uuid.h>
#include <iostream>
#include <string>
#include <vector>
 
#include "km3net-dataformat/definitions/weightlist.hh"
 
#include "km3net-dataformat/offline/MultiHead.hh"
#include "km3net-dataformat/offline/Head.hh"
#include "km3net-dataformat/offline/Evt.hh"
 
#include "JDAQ/JDAQEventIO.hh"
#include "JDAQ/JDAQTimesliceIO.hh"
#include "JDAQ/JDAQSummarysliceIO.hh"
 
#include "JTrigger/JTriggerParameters.hh"
 
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
#include "JOscProb/JOscParameters.hh"
#include "JOscProb/JOscProbInterpolator.hh"
#include "JOscProb/JOscProbHelper.hh"
 
#include "JAAnet/JHead.hh"
#include "JAAnet/JMultiHead.hh"
#include "JAAnet/JHeadToolkit.hh"
 
#include "JAAnet/JAAnetToolkit.hh"
#include "JAAnet/JParticleTypes.hh"
 
#include "JAAnet/JEvtCategoryMap.hh"
 
#include "JReconstruction/JEvt.hh"
 
#include "JSupport/JMeta.hh"
#include "JSupport/JSupport.hh"
#include "JSupport/JFileRecorder.hh"
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JEvtWeightFileScannerSet.hh"
#include "JSupport/JMonteCarloFileSupportkit.hh"
 
#include "JGizmo/JGizmoToolkit.hh"
 
 
/**
 * \file
 * Program for reweighting Monte-Carlo simulations.
 *
 * \author bjjung
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;
 
  typedef JTYPELIST<JRemove<JAllTypes_t, JDAQTimesliceTypes_t>::typelist, JFIT::JEvt>::typelist   typelist;
  
  
  JMultipleFileScanner_t  inputFiles;
  JFileRecorder<typelist> outputFile;
 
  string                  oscProbTable;
  JOscParameters<double>  oscParameters;
 
  JEvtWeightFactorMap     factors{JOscProbInterpolator<>()};
  
  int                     debug;
  
  try {
    
    JParser<> zap;
    
    zap['f'] = make_field(inputFiles);
    zap['o'] = make_field(outputFile);
    zap['P'] = make_field(oscProbTable)  = JPARSER::initialised();
    zap['@'] = make_field(oscParameters) = JPARSER::initialised();
    zap['%'] = make_field(factors)       = JPARSER::initialised();
    zap['d'] = make_field(debug)         = 1;
 
    zap(argc, argv);
  }
  catch(const exception& error) {
    FATAL(error.what() << endl);
  }  
  
  
  // Load oscillation probability table
 
  if (!oscProbTable.empty()) {
  
    JOscProbInterpolator<>& interpolator = static_cast<JOscProbInterpolator<>&>(factors.oscProb.getOscProbInterface());
  
    interpolator.load(oscProbTable.c_str());
    interpolator.set (oscParameters);
  }
  
  
  // Reweight events
 
  JEvtWeightFileScannerSet<> scanners(inputFiles);
  
  if (scanners.setEvtWeightFactor(factors) == 0) {
    FATAL("No weight factors could be set." << endl);
  }
  
  JMultiHead eventHeaders = getMultiHeader(inputFiles); 
  
  JHead commonHeader = scanners.begin()->getHeader();
 
  outputFile.open();
 
  outputFile.put(JMeta(argc, argv));
 
  STATUS(RIGHT(15) << "event" << RIGHT(15) << "weight" << endl);
  
  for (JEvtWeightFileScannerSet<>::iterator scanner = scanners.begin(); scanner != scanners.end(); ++scanner) {
 
    JHead header = scanner->getHeader();
    header.createUUID(); // Ensure UUID is set
 
    commonHeader = commonHeader.getMatch(header);
 
    eventHeaders.merge(header);
 
    while (scanner->hasNext()) {
      
      Evt* event = scanner->next();
 
      if (event == NULL) { 
        WARNING("Event " << scanner->getCounter() << " is empty; skip.");
        continue; 
      }
            
      const double weight = scanner->getWeight(*event);
        
      STATUS(RIGHT     (15)    << scanner->getCounter() <<
             SCIENTIFIC(15, 3) << weight << '\r'); DEBUG(endl);
       
      if (!isfinite(weight)) {
 
        WARNING("Non-finite reweighting factor " << weight << 
                " for event " << scanner->getCounter() << "!");
      }
 
      if (event->w.size() <= WEIGHTLIST_RESCALED_EVENT_RATE) {
        event->w.resize(WEIGHTLIST_RESCALED_EVENT_RATE + 1, 0.0);
      }
 
      const JUUID headerUUID = (uuid_is_null(event->header_uuid) == 0 ? // is valid
                                eventHeaders.getHeaderUUID(*event) : header.UUID);
 
      event->w[WEIGHTLIST_NORMALISATION] = eventHeaders.getNormalisation(headerUUID);
      event->w.at(WEIGHTLIST_RESCALED_EVENT_RATE) = weight;
      
      uuid_copy(event->header_uuid, headerUUID.uuid);
      
      outputFile.put(*event);
    }
  }
 
  Head newHead;
  copy(commonHeader, newHead);
  
  outputFile.put(newHead);
 
  outputFile.put(static_cast<const MultiHead&>(eventHeaders));
  
  JMultipleFileScanner<JRemove<typelist, JAAnetTypes_t>::typelist> io(inputFiles);
 
  io >> outputFile;
  
  outputFile.close();
  
  return 0;
}