JEvtReweightMupage.cc File Reference

Program for reweighting mupage data according to a specifiable TFormula. More...

#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"

Go to the source code of this file.

Functions
int	main (int argc, char **argv)

Detailed Description

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
defined in JEvtWeightFactorMupage. The indices of the physical variables are defined in
the enum JEvtWeightFactorMupage::variables.

Author

bjung

Definition in file JEvtReweightMupage.cc.

Function Documentation

int main	(	int	argc,
		char **	argv
	)

Definition at line 39 of file JEvtReweightMupage.cc.

{
  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;
}