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


namespace {

  /**
   * Neutrino labels for command line interface.
   */
  const std::map<int, std::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;
  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;
}