JEffectiveMass1D.cc

Go to the documentation of this file.

#include <string>
#include <iostream>

#include "km3net-dataformat/offline/Head.hh"
#include "km3net-dataformat/offline/Evt.hh"
#include "km3net-dataformat/offline/Vec.hh"

#include "JDAQ/JDAQEventIO.hh"

#include "JGeometry3D/JCylinder3D.hh"
#include "JGeometry3D/JPosition3D.hh"

#include "JAAnet/JAAnetToolkit.hh"
#include "JAAnet/JHeadToolkit.hh"
#include "JAAnet/JVolume.hh"

#include "JSupport/JTriggeredFileScanner.hh"
#include "JSupport/JEvtWeightFileScannerSet.hh"
#include "JSupport/JMonteCarloFileSupportkit.hh"
#include "JSupport/JSupport.hh"

#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"

#include "JPhysics/JConstants.hh"
#include "JPhysics/JGeane.hh"

#include "JReconstruction/JEvt.hh"
#include "JReconstruction/JEvtToolkit.hh"

#include "JROOT/JManager.hh"

#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"


namespace {

  /**
   * Check if given event has a muon intersecting the given cylinder.
   *
   * \param  event         event
   * \param  cylinder      cylinder
   * \return               true if event has intersecting muon; else false.
   */
  bool hasIntersectingMuon(const Evt&                      event,
                           const JGEOMETRY3D::JCylinder3D& cylinder)
  {
    using namespace std;
    using namespace JPP;

    for (vector<Trk>::const_iterator i = event.mc_trks.cbegin(); i != event.mc_trks.cend(); ++i) {
      
      if (is_finalstate(*i) && is_muon(*i)) {

        const JCylinder3D::intersection_type& intersection = cylinder.getIntersection(getAxis(*i));
        
        const double Lmuon = gWater.getX(i->E, MASS_MUON / getSinThetaC());
        const double Leff  = (intersection.first < 0.0        ? 
                              min(Lmuon, intersection.second) : 
                              min(Lmuon, intersection.second) - intersection.first);

        if (Leff > 0.0) { return true; }
      }
    }

    return false;
  }
}


/**
 * \file
 *
 * Example program to histogram neutrino effective mass for triggered events inside the can volume.\n
 * For genie/gSeaGen events the effective mass for triggered events inside the instrumented volume is also histogrammed.\n
 * The unit of the contents of the histogram is Mton.
 * \author mdejong, bstrandberg
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;


  //------------------------------------------------------------------------------
  // parse command line arguments, check that header can be found in the input file
  //------------------------------------------------------------------------------
  
  JMultipleFileScanner_t       inputFiles;
  string                       outputFile;
  double                       wall;
  bool                         logx;
  int                          numberOfBins;
  string                       detectorFile;
  int                          debug;

  
  try {

    JParser<> zap("Example program to histogram neutrino effective mass for triggered events.");

    zap['f'] = make_field(inputFiles);
    zap['o'] = make_field(outputFile) = "Meff.root";
    zap['R'] = make_field(wall, "Addition margin around the volume in which the considered events must reside") = 0.0;
    zap['X'] = make_field(logx, "Use logarithm of energy");
    zap['N'] = make_field(numberOfBins, "Number of bins in the energy range of the MC simulation") = 10;
    zap['d'] = make_field(debug) = 1;
    zap['a'] = make_field(detectorFile, "Detector file: if not provided, trigger fraction is not calculated") = "";
    
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }

  //-------------------------------------------------------------------------
  // Initialise the detector
  //-------------------------------------------------------------------------
  
  JDetector detector;

  if (detectorFile != "") {

    try {
      load(detectorFile, detector);
    }
    catch(const JException& error) {
      FATAL(error);
    }
  }


  //-------------------------------------------------------------------------
  // Retrieve auxiliary header info
  //-------------------------------------------------------------------------

  JCylinder3D canvol;
  JCylinder3D instvol;

  vector<JVolume> volumes;

  double Xmin = numeric_limits<double>::max();
  double Xmax = numeric_limits<double>::lowest();

  JEvtWeightFileScannerSet<> scanners(inputFiles);

  for (JEvtWeightFileScannerSet<>::const_iterator scanner = scanners.cbegin(); scanner != scanners.cend(); ++scanner) {
    
    const JHead& header = scanner->getHeader();
    
    const JVolume volume(header, logx);

    if (volume.getXmin() < Xmin) { Xmin = volume.getXmin(); }
    if (volume.getXmax() > Xmax) { Xmax = volume.getXmax(); }

    volumes.push_back(volume);

    // the can volume comes from the generator header and the coordinate system is consistent by construction
    JCylinder3D can = get<JCylinder3D>(header);
    canvol.addMargin(wall);

    if (can.getVolume() > canvol.getVolume()) { canvol = can; }

    // as the detector is constructed from Jpp detector file, the coordinate system needs to be synced
    // to the one used in the event generator
    JPosition3D center = getPosition(get<Vec>(header));
    detector -= center;

    JCylinder3D inst( detector.begin(), detector.end() );
    instvol.addMargin(wall);

    if (inst.getVolume() > instvol.getVolume()) { instvol = inst; }

    detector += center;
  }

  
  //------------------------------------------------------------------------------
  // initialise output
  //------------------------------------------------------------------------------

  TFile out(outputFile.c_str(), "recreate");

  TH1::SetDefaultSumw2();

  JManager<int, TH1D> hM(new TH1D("hM[%]", "effective mass (can) [kg]"         , numberOfBins, Xmin, Xmax));
  JManager<int, TH1D> hm(new TH1D("hm[%]", "effective mass (instrumented) [kg]", numberOfBins, Xmin, Xmax));
  JManager<int, TH1D> hNM(new TH1D("hNM[%]", "normalization effective mass (can) [kg]"         , numberOfBins, Xmin, Xmax));
  JManager<int, TH1D> hNm(new TH1D("hNm[%]", "normalization effective mass (instrumented) [kg]", numberOfBins, Xmin, Xmax));


  for (JEvtWeightFileScannerSet<>::iterator scanner = scanners.begin(); scanner != scanners.end(); ++scanner) {

    const size_t scannerIndex = distance(scanners.begin(), scanner);

    const JVolume& volume = volumes[scannerIndex];

    NOTICE("JEffectiveMass1D: Instrumented volume dimensions (zmin, zmax, r): " << instvol.getZmin() << " " << instvol.getZmax() << " " << instvol.getRadius() << endl );

    //------------------------------------------------------------------------------
    // loop over generated events in the input files
    //------------------------------------------------------------------------------

    NOTICE("Scanning generated events for file type " << scannerIndex << endl);
    
    while (scanner->hasNext()) {

      STATUS("event: " << setw(10) << scanner->getCounter() << '\r'); DEBUG(endl);

      const Evt* event = scanner->next();

      const Trk& primary = get_primary(*event);

      const JVertex3D& vertex  = getVertex(*event);

      const bool isValid1 = (canvol.is_inside(vertex)  || hasIntersectingMuon(*event, canvol));
      const bool isValid2 = (instvol.is_inside(vertex) || hasIntersectingMuon(*event, instvol));

      // events are filled to hM if they are inside the can
      hNM[primary.type]->Fill(volume.getX(primary.E), isValid1 ? scanner->getWeight(*event) : 0.0);

      // events are filled to hm if they are inside the instrumented volume
      hNm[primary.type]->Fill(volume.getX(primary.E), isValid2 ? scanner->getWeight(*event) : 0.0);

    } // end loop over generated events
    STATUS(endl);

    if (scanner->getCounter() == 0) {
      FATAL("JEffectiveMass1D: generated events not stored in the input file, JTriggerEfficiency should be run without option -O");
    }

    //------------------------------------------------------------------------------
    // loop over triggered events in the input files
    //------------------------------------------------------------------------------

    NOTICE("Scanning triggered events for file type " << scannerIndex << endl);
    
    JEvtWeightFileScanner<JTriggeredFileScanner<> > in(scanner->getFilelist());
    
    while (in.hasNext()) {

      STATUS("event: " << setw(10) << in.getCounter() << '\r'); DEBUG(endl);

      JTriggeredFileScanner<>::multi_pointer_type mp = in.next();

      const Evt* event = mp;

      const Trk& primary = get_primary(*event);

      const JVertex3D& vertex  = getVertex(*event);

      const bool isValid1 = (canvol.is_inside(vertex)  || hasIntersectingMuon(*event, canvol));
      const bool isValid2 = (instvol.is_inside(vertex) || hasIntersectingMuon(*event, instvol));
        
      hM[primary.type]->Fill(volume.getX(primary.E), isValid1 ? in.getWeight(*event) : 0.0);
      hm[primary.type]->Fill(volume.getX(primary.E), isValid2 ? in.getWeight(*event) : 0.0);

    } // end loop over generated events
    STATUS(endl);
  }
  

  //------------------------------------------------------------------------------
  // convert 'generated' and 'triggered' histograms to effective mass
  //------------------------------------------------------------------------------

  vector<int> pdgIDs = get_keys(hM);

  for (vector<int>::const_iterator i = pdgIDs.begin(); i != pdgIDs.end(); ++i) {

    const int pdgID = *i;

    hM[pdgID]->Divide(hNM[pdgID]);
    hM[pdgID]->Scale(canvol.getVolume()  * DENSITY_SEA_WATER * 1e-6); // Mton
    hm[pdgID]->Divide(hNm[pdgID]);
    hm[pdgID]->Scale(instvol.getVolume() * DENSITY_SEA_WATER * 1e-6); // Mton
  }

  out << hM << hm;
  
  out.Close();
}