JEvtReweightGSeaGen.cc File Reference

Program for reweighting gSeaGen 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 "Jeep/JProperties.hh"
#include "JLang/JToken.hh"
#include "JAAnet/JHead.hh"
#include "JAAnet/JHeadToolkit.hh"
#include "JAAnet/JAAnetToolkit.hh"
#include "JAAnet/JParticleTypes.hh"
#include "JAAnet/JEvtWeightFactorGSeaGen.hh"
#include "JAAnet/JEvtWeightFactorFunction.hh"
#include "JAAnet/JAtmosphericNeutrinoFlux.hh"
#include "JAAnet/JFlux.hh"
#include "JOscProb/JOscParameters.hh"
#include "JOscProb/JOscProbInterpolator.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 gSeaGen data according to a specifiable TFormula.

One TFormula can be specified for each neutrino type. The TFormula may take any number of parameters, but is limited to the physical variables
defined in JEvtWeightFactorGSeaGen. The indices of the physical variables are defined in
the enum JEvtWeightFactorGSeaGen::variables.
Parameter settings can be specified behind the TFormula separated by semicolons, e.g.:
-#anumu="[0]*(x[1]-[1]); p0=1; p1=2.5;"

Author

bjjung

Definition in file JEvtReweightGSeaGen.cc.

Function Documentation

int main	(	int	argc,
		char **	argv
	)

Definition at line 69 of file JEvtReweightGSeaGen.cc.

{
  using namespace std;
  using namespace JPP;
  
  typedef typename JTYPELIST<JAAnetTypes_t, JMetaTypes_t>::typelist   typelist;
  typedef JEvtWeightFactorFunction<JAtmosphericNeutrinoFlux, JFlux>   JFluxFunction_t;
  
  
  JMultipleFileScanner_t            inputFiles;
  JFileRecorder<typelist>           outputFile;

  string                            oscProbTable;
  JOscParameters                    oscParameters;

  map<int, JEvtWeightFactorGSeaGen> formulae;    
  
  int                               debug;
  
  try {
    
    JProperties formulaMap(JEquationParameters("=","\n\r","./","#"));

    for (map<int, string>::const_iterator i = nuLabels.cbegin(); i != nuLabels.cend(); ++i) {
      formulaMap[i->second] = formulae[i->first];
    }
    
    JParser<> zap;
    
    zap['f'] = make_field(inputFiles);
    zap['o'] = make_field(outputFile);
    zap['P'] = make_field(oscProbTable);
    zap['@'] = make_field(oscParameters) = JPARSER::initialised();
    zap['#'] = make_field(formulaMap)    = JPARSER::initialised();
    zap['d'] = make_field(debug)         = 1;

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


  // Create atmospheric neutrino flux function
  
  const JFluxFunction_t fluxFunction = make_atmosphericNeutrinoFluxFunction(oscProbTable.c_str(), oscParameters);

  
  // Create file scanner

  JEvtWeightFileScanner<> scanner(inputFiles);
  
  JHead& header = scanner.getHeader();
  header.createUUID();
  
  if (!is_gseagen(header)) {
    FATAL("MC-header is incompatible with gSeaGen.");
  }
  
  
  // Define reweighting functions

  for (map<int, JEvtWeightFactorGSeaGen>::iterator i = formulae.begin(); i != formulae.end(); ++i) {

    const TString expression = i->second.GetExpFormula();
    
    if        (!expression.IsNull() && i->second.IsValid() && i->second.GetNdim() <= JEvtWeightFactorGSeaGen::NUMBER_OF_VARIABLES) {

      i->second.configure(fluxFunction);

      scanner.setEvtWeightFactor(i->first, make_weightFactor(i->second));
      
    } else if (expression.IsNull()) {

      continue;
      
    } else {
      
      FATAL("Invalid formula " << expression << " (number of dimensions must be < " <<
            JEvtWeightFactorGSeaGen::NUMBER_OF_VARIABLES << "; see `JAANET::JEvtWeightFactorGSeaGen`).");
    }
  }
  
  
  // 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;
}