JEvtReweightMupage.cc File Reference

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

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

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 43 of file JEvtReweightMupage.cc.

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