JEventShapeVariables.cc

Go to the documentation of this file.

#include <iostream>
#include <string>
#include <vector>
 
#include "km3net-dataformat/offline/Evt.hh"
#include "km3net-dataformat/offline/Trk.hh"
 
#include "JAAnet/JAAnetToolkit.hh"
#include "JAAnet/JLorentzBoost.hh"
#include "JAAnet/JEvtCategoryMap.hh"
 
#include "JOscProb/JOscParameters.hh"
#include "JOscProb/JOscProbInterpolator.hh"
 
#include "JSirene/JEventShapeVariables.hh"
 
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
 
#include "JSupport/JSupport.hh"
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JEvtWeightFileScannerSet.hh"
#include "JSupport/JMonteCarloFileSupportkit.hh"
 
#include "Jeep/JeepToolkit.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
#include "TFile.h"
#include "TH1D.h"
 
 
/**
 * \file
 *
 * Auxiliary program to plot event shape variables.
 * \author bjjung
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
 
  JMultipleFileScanner_t        inputFiles;  
  string                        outputFile;
  string                        detectorFile;
  
  JLimit_t                      numberOfEvents;
 
  bool                          useBoost;
 
  bool                          useWeights;
  string                        oscProbTable;
  JOscParameters<double>        oscParameters;
  JFluxMap                      fluxMap{JOscProbInterpolator<>()};  
  
  int                           debug;
 
  try { 
 
    JParser<> zap("Auxiliary program to test event shape variables.");
 
    zap['f'] = make_field(inputFiles);
    zap['o'] = make_field(outputFile);
    zap['a'] = make_field(detectorFile)        = "";
    zap['n'] = make_field(numberOfEvents)      = JLimit::max();
    zap['B'] = make_field(useBoost);
    zap['W'] = make_field(useWeights);
    zap['@'] = make_field(fluxMap)
      = JPARSER::initialised();
    zap['P'] = make_field(oscProbTable)
      = JPARSER::initialised();
    zap['#'] = make_field(oscParameters)
      = JPARSER::initialised();
    zap['d'] = make_field(debug)               = 3;
    
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }
 
  if (useWeights && numberOfEvents != JLimit::max()) {
    FATAL("Cannot apply weighting to limited number of events.");
  }
  
 
  // Load detector file
  
  JDetector detector;
 
  if (!detectorFile.empty()) {
    try {
      load(detectorFile, detector);
    }
    catch(const JException& error) {
      FATAL(error);
    }
  }
 
  JCylinder3D fiducialVolume;
 
  if (!detector.empty()) {
    fiducialVolume.set(detector.begin(), detector.end());
  }
 
  
  // Load oscillation probability table
  
  if (!oscProbTable.empty()) {
  
    JOscProbInterpolator<>& interpolator = static_cast<JOscProbInterpolator<>&>(fluxMap.oscProb.getOscProbInterface());
  
    interpolator.load(oscProbTable.c_str());
    interpolator.set (oscParameters);
  }
 
 
  // Create file scanners
 
  JEvtWeightFileScannerSet<> scanners(inputFiles);
  
  if (scanners.setFlux(fluxMap) == 0) {
    WARNING("No flux function set." << endl);
  }
  
 
  // Define histograms
 
  TFile out(outputFile.c_str(), "recreate");
 
  TH1D hT("hT", "thrust",       50, 0.5, 1.0);
  
  TH1D hA("hA", "aplanarity",   50, 0.0, 0.5);
  TH1D hS("hS", "sphericity",  100, 0.0, 1.0);
  TH1D hc("hc", "circularity", 100, 0.0, 1.0); 
  TH1D hp("hp", "planar flow", 100, 0.0, 1.0);
  
  TH1D hC("hC", "C-variable",  100, 0.0, 1.0);
  TH1D hD("hD", "D-variable",  100, 0.0, 1.0);
  
  TH1D hH10("hH10", "H10",     100, 0.0, 1.0);
  TH1D hH20("hH20", "H20",     100, 0.0, 1.0);
  TH1D hH30("hH30", "H30",     100, 0.0, 1.0);
  TH1D hH40("hH40", "H40",     100, 0.0, 1.0);
 
  TH1D hq("hq", "Logarithmic eigenvalue ratio hit inertia tensor", 200, -20.0, 0.0);  
    
    
  // Loop over events
 
  STATUS(RIGHT(15) << "header ID" << RIGHT(15) << "event" << RIGHT(15) << "weight" << endl);
  
  for (JEvtWeightFileScannerSet<>::iterator scanner = scanners.begin(); scanner != scanners.end(); ++scanner) {
    
    scanner->setLimit(numberOfEvents);
 
    const JHead& header = scanner->getHeader();
 
    if (detector.empty()) { // Use can volume if detector file has not been specified
      fiducialVolume = getCylinder(header);
    }
    
    while (scanner->hasNext()) {
      
      Evt* event = scanner->next();
      
      const double weight = (useWeights ? scanner->getWeight(*event) : 1.0);
      
      STATUS(RIGHT     (15)    << distance(scanners.begin(), scanner) <<
             RIGHT     (15)    << scanner->getCounter()               <<
             SCIENTIFIC(15, 3) << weight << '\r'); DEBUG(endl);
 
      if (useBoost) { boostToCOM(*event); }
 
      const JVertex3D vertex = getVertex(*event);
 
      const JSphericityTensor S1(*event, 1.0);
      const JSphericityTensor S2(*event, 2.0);
 
      const JFoxWolframMoments H(*event, fiducialVolume, 4);
 
      const JHitInertiaTensor Q(*event, vertex);
 
      const double thrust = getThrust(*event);
 
      const double aplanarity  = S2.getAplanarity();
      const double sphericity  = S2.getSphericity();
      const double circularity = S2.getCircularity();
      const double planarflow  = S2.getPlanarFlow();
 
      const double C = S1.getC();
      const double D = S1.getD();
 
      const double q = Q.getEigenvalueRatio();
 
      hT.Fill(thrust,        weight);
 
      hA.Fill(aplanarity,    weight);
      hS.Fill(sphericity,    weight);
      hc.Fill(circularity,   weight);
      hp.Fill(planarflow,    weight);
 
      hC.Fill(C,             weight);
      hD.Fill(D,             weight);
 
      hH10.Fill(H[1]/H[0],   weight);
      hH20.Fill(H[2]/H[0],   weight);
      hH30.Fill(H[3]/H[0],   weight);
      hH40.Fill(H[4]/H[0],   weight);      
 
      hq.Fill(log10(q), weight);
    }
  }
  
  STATUS(endl); 
 
  out.Write();
  out.Close();
 
  return 0;
}