JEvtReweightMupage.cc

Go to the documentation of this file.

#include <iostream>
#include <iomanip>
#include <string>
#include <vector>

#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/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<Head, Evt, JMeta>::typelist   typelist;
  typedef JToken<';'>                                      JToken_t;
  
  
  string                                   inputFile;
  string                                   outputFile;
  
  string                                   formula;
  vector<JToken_t>                         parameters;
  
  int                                      debug;
  
  try {
    
    JParser<> zap;
    
    zap['f'] = make_field(inputFile);
    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(inputFile);
  
  const JHead& header = scanner.getHeader();
  
  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_ANTIMUON, make_weightFactor(reweighter))) {
    FATAL("Setting of reweighting factor was unsuccessful.");
  }
  
  
  // Store events with new weights
    
  JFileRecorder<typelist> out(outputFile.c_str());

  out.open();

  Head head;
  copy(header, head);
  out.put(head);
    
  out.put(JMeta(argc, argv));
  
  for (JMultipleFileScanner<JMeta> in(scanner.getFilelist()); in.hasNext(); ) {
    out.put(*in.next());
  }

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

    if (event != NULL) {

      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() < RESCALED_WEIGHT_INDEX + 1) {
        event->w.resize(RESCALED_WEIGHT_INDEX + 1);
      }
      
      event->w.at(RESCALED_WEIGHT_INDEX) = weight;

      out.put(*event);
      
    } else {

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

  out.close();
  
  return 0;
}