JHistHDF.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <fstream>
#include <iomanip>
 
#include "km3net-dataformat/offline/Head.hh"
#include "km3net-dataformat/offline/Evt.hh"
#include "km3net-dataformat/offline/Hit.hh"
 
#include "JAAnet/JHead.hh"
#include "JAAnet/JAAnetToolkit.hh"
 
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JMonteCarloFileSupportkit.hh"
#include "JSupport/JSupport.hh"
 
#include "JPhysics/JPDFToolkit.hh"
#include "JPhysics/JPDFTransformer.hh"
#include "JPhysics/JDispersion.hh"
#include "JPhysics/Antares.hh"
#include "JPhysics/KM3NeT.hh"
#include "JIO/JFileStreamIO.hh"
 
#include "JGeometry3D/JPosition3D.hh"
#include "JGeometry3D/JDirection3D.hh"
#include "JGeometry3D/JAxis3D.hh"
#include "JGeometry3D/JGeometry3DToolkit.hh"
 
#include "JTools/JConstants.hh"
#include "JTools/JHistogram1D_t.hh"
#include "JTools/JHistogramMap_t.hh"
#include "JTools/JTransformableMultiHistogram.hh"
#include "JTools/JFunction1D_t.hh"
#include "JTools/JFunctionalMap_t.hh"
 
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JPMTRouter.hh"
 
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
 
namespace {
 
  using namespace JPP;
 
  /**
   * Hit classifier class
   */
  template<class JMultiHistogram_t>
  class JHitClassifier
  {
  public:
    /**
     *  Constructor
     *  \param head         MC file header
     *  \param detector     detector
     */
    JHitClassifier(const JHead& head, JDetector& detector) :
      KM3Sim(is_km3sim(head))
    {
      if (KM3Sim) {
        detector -= JCenter3D(detector.begin(), detector.end());   // For KM3Sim, you get it from the detector file
      }
      else {
        detector -= get<JPosition3D>(head);                        // For KM3, the MC origin is stored in the header
      }
    }
 
 
    /**
     *  Function to classify hits
     * \param  track          muon track
     * \param  hit            hit
     * \return                1 for muon; 2 for shower; -1 for anything else
     */
    int operator()(const Trk& track, const Hit& hit) const
    {
      if (KM3Sim)
        return classify_km3sim(track, hit);
      else
        return classify_km3(track, hit);
    }
 
 
    /**
     * Auxillary function to calculate the shower energy.
     * For KM3Sim, (ORCA neutrinos < 100GeV) this will be the hadronic shower energy.
     * For KM3 and JSirene, (ARCA or atmospheric muons) this will be the bremsstrahlung showers along the track.
     *
     * \param event           Monte Carlo event
     * \param track           muon track
     * \param axis            PMT axis
     *
     * \return                energy [GeV]
     */
    inline double getEnergy(const Evt& event, const Trk& track, const JAxis3D& axis) const
    {
      if (KM3Sim)
        return get_hadronic_cascade(event).E;
      else
        return gWater(track.E, axis.getZ());
    }
 
 
  private:
    const bool KM3Sim;        //!< true if KM3Sim was used (ORCA neutrinos < 100 GeV); else false
 
 
    /**
     * Function for classifying hits generated with KM3
     * 
     * \param  track          muon track
     * \param  hit            hit
     * \return                1 for muon; 2 for shower; -1 for anything else
     */
    int classify_km3(const Trk& track, const Hit& hit) const
    {
      if (track.id == hit.origin) {
        switch(hit.type) {
        case HIT_TYPE_MUON_DIRECT:
          return 1;
        case HIT_TYPE_MUON_SCATTERED:
          return 1;
        case HIT_TYPE_BREMSSTRAHLUNG_DIRECT:
          return 2;
        case HIT_TYPE_BREMSSTRAHLUNG_SCATTERED:
          return 2;
        default:
          return -1;
        }
      }
      return -1;
    }
 
 
    /**
     * Function for classifying hits generated with KM3Sim
     *
     * \param  track          muon track
     * \param  hit            hit
     * 
     * \return                1 for muon; 2 for shower
     */
    int classify_km3sim(const Trk& track, const Hit& hit) const
    {
      if (from_muon(hit)) 
        return 1;
      else
        return 2;
    }
  };
}
 
 
/**
 * \file
 *
 * Program to histogram event-by-event data of muon light for making PDFs.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
 
  JMultipleFileScanner<Evt> inputFile;
  JLimit_t&   numberOfEvents = inputFile.getLimit();
  string      outputFile;
  string      detectorFile;
  int         debug;
 
  try {
 
    JParser<> zap("Program to histogram event-by-event data of muon light for making PDFs.");
 
    zap['f'] = make_field(inputFile);
    zap['o'] = make_field(outputFile);
    zap['n'] = make_field(numberOfEvents)      = JLimit::max();
    zap['a'] = make_field(detectorFile);
    zap['d'] = make_field(debug)               = 1;
 
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }
 
 
  using namespace KM3NET;
 
 
  JDetector detector;
 
  try {
    load(detectorFile, detector);
  }
  catch(const JException& error) {
    FATAL(error);
  }
 
  JHead head(JMultipleFileScanner<Head>(inputFile).getHeader());      // Monte Carlo header
 
  const JPMTRouter pmtRouter(detector);
 
  const double P_atm = NAMESPACE::getAmbientPressure();
  const double wmax  = getMaximalWavelength();
    
  const JDispersion dispersion(P_atm);
 
  typedef JHistogram1D_t::abscissa_type                             abscissa_type;
 
  typedef JTransformableMultiHistogram<JHistogram1D_t, 
                                       JMAPLIST<JHistogramMap_t,
                                                JHistogramGridMap_t,
                                                JHistogramGridMap_t>::maplist> JMultiHistogram_t;
 
  typedef JPDFTransformer<3, abscissa_type>                         JFunction3DTransformer_t;
 
  JMultiHistogram_t h0;   // occurrence rate of PMT (used for normalisation)
  JMultiHistogram_t h1;   // light from muon
  JMultiHistogram_t h2;   // light from showers
 
  const JHitClassifier<JMultiHistogram_t> classifier(head, detector);
 
  const double cmin = dispersion.getKmin (wmax);
 
  h1.transformer.reset(new JFunction3DTransformer_t(35.0, 2, cmin, 0.0,  NAMESPACE::getAngularAcceptance, 0.10));
  h2.transformer.reset(new JFunction3DTransformer_t(35.0, 2, cmin, 0.0,  NAMESPACE::getAngularAcceptance, 0.10));
 
 
  // configure bins
 
  const double R[] = { 0.0, 5.0, 10.0, 20.0, 30.0, 40.0, 50.0, 70.0, 80.0, 90.0, 100.0, 170.0, 250.0 };
 
  for (int i = 0; i != sizeof(R)/sizeof(R[0]); ++i) {
 
    const double R_m = R[i];
    
    const double grid  = 10.0 +  0.0 * R_m/100.0;                         // [deg]
    const double alpha = 2.0 * sqrt(1.0 - cos(grid * PI / 180.0));        // azimuth angle unit step size
    
    const int    number_of_theta_points = max(2, (int) (180.0/(1.4 * grid)));
    const double theta_step             = PI / (number_of_theta_points + 1);
    
    for (double theta = -0.5*theta_step; theta < PI + theta_step; theta += theta_step) {
      
      const int    number_of_phi_points = max(2, (int) (PI * sin(theta) / alpha));
      const double phi_step             = PI / (number_of_phi_points + 1);
      
      for (double phi = -0.5*phi_step; phi < PI + phi_step; phi += phi_step) {
 
        for (JMultiHistogram_t* buffer[] = { &h0, &h1, &h2, NULL }, **histogram = buffer; *histogram != NULL; ++histogram) {
          (**histogram)[R_m][theta][phi];
        }
      }
    }
  }
 
  for (JMultiHistogram_t::super_iterator 
         i1 = h1.super_begin(),
         i2 = h2.super_begin(); i1 != h1.super_end(); ++i1, ++i2) {
 
    for (double buffer[] = { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 7.5, 10.0, 15.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 85.0, 100.0, 125.0, 150.0, 200.0, 500.0, -1.0 }, *x = buffer; *x >= 0.0; ++x) {
      i1.getValue()[*x];
      i2.getValue()[*x];
    }
  }
 
 
  while (inputFile.hasNext()) {
 
    NOTICE("event: " << setw(10) << inputFile.getCounter() << '\r'); STATUS(endl);
 
    const Evt* event = inputFile.next();
 
    for (vector<Trk>::const_iterator track = event->mc_trks.begin(); track != event->mc_trks.end(); ++track) {
 
      if (is_muon(*track)) {
 
        // track parameters
 
        double       E   = track->E;
        JPosition3D  pos = getPosition (*track);
        JDirection3D dir = getDirection(*track);
        double       t0  = track->t;
 
        const JRotation3D R(dir);
 
        pos.rotate(R);
 
        JRange<double> Z(0.0, gWater(E));   // range of muon
 
        for (vector<Hit>::const_iterator hit = event->mc_hits.begin(); hit != event->mc_hits.end(); ++hit) {
        
          try {
        
            JAxis3D axis = pmtRouter.getPMT(hit->pmt_id);
        
            axis.transform(R, pos);
        
            const double E1 = classifier.getEnergy(*event, *track, axis);
            const double t1 = t0 + (axis.getZ() + axis.getX()*getTanThetaC()) / getSpeedOfLight();
        
            const double R_m   = axis.getX();
            const double theta = axis.getTheta();
            const double phi   = fabs(axis.getPhi());
            const double dt    = getTime(*hit) - t1;
            const double npe   = getNPE(*hit);
 
            switch (classifier(*track, *hit)) {
 
            case 1:
              h1.fill(R_m, theta, phi, dt, npe);              // light from muon
              break;
 
            case 2:
              h2.fill(R_m, theta, phi, dt, npe/max(1.0, E1)); // light from shower
              break;
            }
          }
          catch(const exception& error) {
            FATAL(error.what() << endl);
          }
        }
 
 
        for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
          
          JPosition3D P = module->getPosition();
 
          P.rotate(R);
 
          P.sub(pos);
 
          if (sqrt(P.getX()*P.getX() + P.getY()*P.getY()) <= h0.getXmax()) {
 
            for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
          
              JAxis3D axis = *pmt;
 
              axis.transform(R, pos);
 
              if (Z(axis.getZ() - axis.getX()/getTanThetaC())) {
                h0.fill(axis.getX(), axis.getTheta(), fabs(axis.getPhi()), 0.0, 1.0);
              }
            }
          }
        }
      }
    }
  }
  NOTICE(endl);
 
 
  // output
 
  JFileStreamWriter out(outputFile.c_str());
 
  for (const JMultiHistogram_t* buffer[] = { &h0, &h1, &h2, NULL }, **i = buffer; *i != NULL; ++i) {
    out.store(**i);
  }
 
  out.close();
}