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/JEvtWeightToolkit.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;
  JFluxMap                      fluxMaps;
  
  string                        oscProbTable;
  JOscParameters<double>        oscParameters;
  
  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(fluxMaps)
      = 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());
  }

  // 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);  
  

  // Create file scanners

  JEvtWeightFileScannerSet<> scanners(inputFiles);

  if (!oscProbTable.empty()) {

    JOscProbInterpolator<> interpolator(oscProbTable.c_str(), oscParameters);

    fluxMaps.configure(make_oscProbFunction(interpolator));
  }

  if (scanners.setFlux(fluxMaps) == 0) {
    WARNING("No flux function set." << endl);
  }
    
  STATUS(RIGHT(15) << "header ID" << RIGHT(15) << "event" << RIGHT(15) << "weight" << endl);

  // Loop over events
  
  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;
}