JEvtReweightMupage.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 "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

#include "JLang/JToken.hh"

#include "JAAnet/JHead.hh"
#include "JAAnet/JHeadToolkit.hh"
#include "JAAnet/JAAnetToolkit.hh"
#include "JAAnet/JEvtWeightFactorMupage.hh"
#include "JAAnet/JEvtWeightFactorFunction.hh"

#include "JSupport/JMeta.hh"
#include "JSupport/JSupport.hh"
#include "JSupport/JFileRecorder.hh"
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JEvtWeightFileScanner.hh"

#include "JGizmo/JGizmoToolkit.hh"


/**
 * \file
 * Program for reweighting mupage data according to a specifiable TFormula.
 *
 * The TFormula may take any number of parameters, but is limited to the physical variables\n
 * defined in `JEvtWeightFactorMupage`. The indices of the physical variables are defined in\n
 * the enum `JEvtWeightFactorMupage::variables`.
 *
 * \author bjung
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;

  typedef typename JTYPELIST<JAAnetTypes_t, JMetaTypes_t>::typelist   typelist;  
  typedef JToken<';'>                                                 JToken_t;
  
  
  JMultipleFileScanner_t                   inputFiles;
  JFileRecorder<typelist>                  outputFile;  
  
  string                                   formula;
  vector<JToken_t>                         parameters;
  
  int                                      debug;
  
  try {
    
    JParser<> zap;
    
    zap['f'] = make_field(inputFiles);
    zap['o'] = make_field(outputFile);
    zap['F'] = make_field(formula,        "Reweighting formula, e.g.: \"[0] + [1]*x[0]\"")  = "";
    zap['@'] = make_field(parameters,     "Parameter values, e.g.: \"p0 = 0.5; p1 = 0.2;\"") = JPARSER::initialised();
    zap['d'] = make_field(debug)                                                            = 1;

    zap(argc, argv);
  }
  catch(const exception& error) {
    FATAL(error.what() << endl);
  }

  
  // Create file scanner

  JEvtWeightFileScanner<> scanner(inputFiles);
  
  JHead& header = scanner.getHeader();
  header.createUUID();
  
  if (!is_mupage(header)) {
    FATAL("MC-header is incompatible with MUPAGE.");
  }
  
  
  // Define reweighting function
  
  JEvtWeightFactorMupage reweighter("reweighter", formula.c_str());

  if (reweighter.GetNdim() > JEvtWeightFactorMupage::NUMBER_OF_VARIABLES) {
        
    FATAL("Invalid formula " << formula << " (number of dimensions must be < " <<
          JEvtWeightFactorMupage::NUMBER_OF_VARIABLES << "; see `JAANET::JEvtWeightFactorMupage`).");
  }
  
  if (size_t(reweighter.GetNpar()) != parameters.size()) {
        
    FATAL("Not all parameters have been specified.");
  }


  // Read parameter settings

  try {
    
    for (vector<JToken_t>::const_iterator i = parameters.begin(); i != parameters.cend(); ++i) {
      
      const int index = getParameter(*i);
      
      reweighter.SetParameter(index, getValue(*i, 0));
    }
  }
  catch (JLANG::JParseError& error) {
    FATAL(error << endl);
  }
  
  
  // Set reweighting factors for mupage-files
  
  if (!scanner.setEvtWeightFactor(TRACK_TYPE_MUON,     make_weightFactor(reweighter)) ||
      !scanner.setEvtWeightFactor(TRACK_TYPE_ANTIMUON, make_weightFactor(reweighter))) {
    FATAL("Setting of reweighting factor was unsuccessful.");
  }
  
  
  // Store events with new weights

  outputFile.open();

  Head head;
  copy(header, head);
  outputFile.put(head);
  
  outputFile.put(JMeta(argc, argv));

  STATUS(RIGHT(15) << "event" << RIGHT(15) << "weight" << endl);
    
  while (scanner.hasNext()) {
      
    const Evt* event = scanner.next();

    if (event != NULL) {

      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() << "!");
      }

      Evt out = *event;

      uuid_copy(out.header_uuid, header.UUID.uuid);       

      if (out.w.size() <= WEIGHTLIST_RESCALED_EVENT_RATE) {
        out.w.resize(WEIGHTLIST_RESCALED_EVENT_RATE + 1, 0.0);
      }

      out.w.at(WEIGHTLIST_NORMALISATION) = scanner.getNormalisation(*event);
      out.w.at(WEIGHTLIST_RESCALED_EVENT_RATE) = weight;

      outputFile.put(out);
      
    } else {

      WARNING("Event " << scanner.getCounter() << " is empty; skip.");
    }
  }

  JMultipleFileScanner<JTYPELIST<JMetaTypes_t, MultiHead>::typelist> io(inputFiles);

  io >> outputFile;

  outputFile.close();
  
  return 0;
}