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 "JDAQ/JDAQEventIO.hh"
#include "JDAQ/JDAQTimesliceIO.hh"
#include "JDAQ/JDAQSummarysliceIO.hh"
#include "JTrigger/JTriggerParameters.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/JMultiHead.hh"
#include "JAAnet/JHeadToolkit.hh"
#include "JAAnet/JAAnetToolkit.hh"
#include "JAAnet/JParticleTypes.hh"
#include "JAAnet/JEvtWeightFactorGSeaGen.hh"
#include "JAAnet/JEvtWeightFactorFunction.hh"
#include "JAAnet/JEvtWeightFactorMultiParticle.hh"
#include "JAAnet/JAtmosphericNeutrinoFlux.hh"
#include "JAAnet/JFlux.hh"
#include "JOscProb/JOscParameters.hh"
#include "JOscProb/JOscProbInterpolator.hh"
#include "JReconstruction/JEvt.hh"
#include "JSupport/JMeta.hh"
#include "JSupport/JSupport.hh"
#include "JSupport/JFileRecorder.hh"
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JEvtWeightFileScannerSet.hh"
#include "JSupport/JMonteCarloFileSupportkit.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 80 of file JEvtReweightGSeaGen.cc.

{
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;

  typedef JTYPELIST<JRemove<JAllTypes_t, JDAQTimesliceTypes_t>::typelist, JFIT::JEvt>::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 flux function

  JFluxMultiParticle fluxFunction;
  
  const JFluxFunction_t atmFlux = make_atmosphericNeutrinoFluxFunction(oscProbTable.c_str(), oscParameters);

  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(atmFlux);

      fluxFunction.insert(i->first, make_fluxFunction(i->second));
      
    } else if (expression.IsNull()) {

      fluxFunction.insert(i->first, atmFlux);
      
    } else {
      
      FATAL("Invalid formula " << expression << " (number of dimensions must be < " <<
            JEvtWeightFactorGSeaGen::NUMBER_OF_VARIABLES << "; see `JAANET::JEvtWeightFactorGSeaGen`).");
    }
  }
  
  
  // Reweight events

  outputFile.open();

  outputFile.put(JMeta(argc, argv));

  JEvtWeightFileScannerSet<> scanners(inputFiles);
  scanners.setFlux(fluxFunction);

  JMultiHead eventHeaders = getMultiHeader(inputFiles); 
  
  JHead commonHeader = scanners.begin()->getHeader();

  STATUS(RIGHT(15) << "event" << RIGHT(15) << "weight" << endl);
  
  for (JEvtWeightFileScannerSet<>::iterator scanner = scanners.begin(); scanner != scanners.end(); ++scanner) {

    JHead header = scanner->getHeader();
    header.createUUID(); // Ensure UUID is set

    commonHeader = commonHeader.getMatch(header);

    eventHeaders.merge(header);

    while (scanner->hasNext()) {
      
      Evt* event = scanner->next();

      if (event == NULL) { 
        WARNING("Event " << scanner->getCounter() << " is empty; skip.");
        continue; 
      }
            
      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() << "!");
      }

      if (event->w.size() <= WEIGHTLIST_RESCALED_EVENT_RATE) {
        event->w.resize(WEIGHTLIST_RESCALED_EVENT_RATE + 1, 0.0);
      }

      const JUUID headerUUID = (uuid_is_null(event->header_uuid) == 0 ? // is valid
                                eventHeaders.getHeaderUUID(*event) : header.UUID);

      event->w[WEIGHTLIST_NORMALISATION] = eventHeaders.getNormalisation(headerUUID);
      event->w.at(WEIGHTLIST_RESCALED_EVENT_RATE) = weight;
      
      uuid_copy(event->header_uuid, headerUUID.uuid);
      
      outputFile.put(*event);
    }
  }

  Head newHead;
  copy(commonHeader, newHead);
  
  outputFile.put(newHead);

  outputFile.put(static_cast<const MultiHead&>(eventHeaders));
  
  JMultipleFileScanner<JRemove<typelist, JAAnetTypes_t>::typelist> io(inputFiles);

  io >> outputFile;
  
  outputFile.close();
  
  return 0;
}