JEvtReweightGSeaGen.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 "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"


namespace {

  /**
   * Neutrino labels for command line interface.
   */
  const map<int, string> nuLabels {
    { JAANET::TRACK_TYPE_ANTINUE,   "anue"   },
    { JAANET::TRACK_TYPE_NUE,       "nue"    },
    { JAANET::TRACK_TYPE_ANTINUMU,  "anumu"  },
    { JAANET::TRACK_TYPE_NUMU,      "numu"   },
    { JAANET::TRACK_TYPE_ANTINUTAU, "anutau" },
    { JAANET::TRACK_TYPE_NUTAU,     "nutau"  }
  };
}


/**
 * \file
 * 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\n
 * defined in `JEvtWeightFactorGSeaGen`. The indices of the physical variables are defined in\n
 * the enum `JEvtWeightFactorGSeaGen::variables`.\n
 * Parameter settings can be specified behind the TFormula separated by semicolons, e.g.:\n
 *    -\#anumu="[0]*(x[1]-[1]); p0=1; p1=2.5;"
 *
 * \author bjjung
 */
int main(int argc, char **argv)
{
  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;
}