Jpp/public/JHistHDE_8cc_source.html

#include <string>

#include <iostream>

#include <fstream>

#include <iomanip>

#include <stdexcept>


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

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

#include "km3net-dataformat/offline/Hit.hh"

#include "km3net-dataformat/offline/Trk.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/JDispersion.hh"

#include "JPhysics/Antares.hh"

#include "JPhysics/KM3NeT.hh"

#include "JIO/JFileStreamIO.hh"


#include "JGeometry3D/JGeometry3DToolkit.hh"

#include "JGeometry3D/JPosition3D.hh"

#include "JGeometry3D/JDirection3D.hh"

#include "JGeometry3D/JAxis3D.hh"


#include "JPhysics/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"


#include "JPhysics/JPDFTransformer.hh"

#include "JPhysics/JPDF.hh"


#include "JAAnet/JParticleTypes.hh"


/**

 * \file

 *

 * Program to histogram event-by-event and track-by-track for making PDFs of lights of individual particles.

 * \author jseneca

 */


inline std::vector<Int_t> getHitRemapping(const std::vector<Trk>* tracklist)

{

  using namespace std;

  using namespace JPP;


  vector<Int_t> hit_remapping;


  for (vector<Trk>::const_iterator i = tracklist->begin(); i != tracklist->end(); ++i) {


    //Skip pathological neutrino that is not the primary incoming neutrino


    if(i->type == TRACK_TYPE_NUMU && i != tracklist->begin()) {

      continue;

    }


    hit_remapping.push_back(i->id);

  }


  return hit_remapping;

}


int main(int argc, char **argv)

{

  using namespace std;

  using namespace JPP;


  JMultipleFileScanner<Evt> inputFile;

  //JLimit_t&   numberOfEvents = inputFile.getLimit();

  string      outputFile = "J%s.HDE"; //'s' for single (particle), or simulation

  string      detectorFile;

  vector<int> particles;

  uint        fix_bug;

  int         debug;


  try {


    JParser<> zap("Program to histogram event-by-event data of shower 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['p'] = make_field(particles);

    zap['b'] = make_field(fix_bug);

    zap['d'] = make_field(debug)                  = 1;


    if (zap.read(argc, argv) != 0)

      return 1;

  }

  catch(const exception &error) {

    FATAL(error.what() << endl);

  }


  JDetector detector;


  if(outputFile.find('%') != string::npos){

    sort(particles.begin(), particles.end());

    string prefix_t = "";

    for(vector<int>::iterator ptype = particles.begin(); ptype != particles.end(); ptype++){

      prefix_t += to_string((long long int)*ptype) + "_";

    }

    prefix_t.erase(prefix_t.size() - 1);

    string::size_type pos_1 = outputFile.find('%');

    outputFile.replace(pos_1, 1, prefix_t);

  }


  try {

    load(detectorFile, detector);

  }

  catch(const JException& error) {

    FATAL(error);

  }


  JHead head(JMultipleFileScanner<Head>(inputFile).getHeader());  // Monte Carlo header


  const Vec  offset = getOffset(head);


  NOTICE("Apply detector offset " << offset << endl);


  detector -= getPosition(offset);


  const JPMTRouter pmtRouter(detector);


  const double P_atm = NAMESPACE::getAmbientPressure();

  const double wmin  = getMinimalWavelength();

  const double wmax  = getMaximalWavelength();


  const JDispersion dispersion(P_atm);


  const double ng[] = { dispersion.getIndexOfRefractionGroup(wmax),

                        dispersion.getIndexOfRefractionGroup(wmin) };


  typedef JHistogram1D_t::abscissa_type                            abscissa_type;


  typedef JTransformableMultiHistogram<JHistogram1D_t,

                                       JMAPLIST<JHistogramMap_t,

                                                JHistogramMap_t,

                                                JHistogramMap_t,

                                                JHistogramGridMap_t,

                                                JHistogramGridMap_t>::maplist> JMultiHistogram_t;


  typedef JPDFTransformer<5, abscissa_type>  JFunction5DTransformer_t;


  JMultiHistogram_t h0;   // occurrence rate of PMT (used for normalisation)

  JMultiHistogram_t h1;   // light from track cascade


  h1.transformer.reset(new JFunction5DTransformer_t(21.5, 2, ng[0], 0.0, JGeant(JGeanx(1.00, -2.2)), 1e-2, NAMESPACE::getAngularAcceptance, 0.05));


  set<double> C;  // cosine emission angle


  JQuadrature qeant(-1.0, 1.0, 30, JGeanx(1.00, -2.2));


  for (JQuadrature::const_iterator i = qeant.begin(); i != qeant.end(); ++i) {

    C.insert(i->getX());

  }


  C.insert(-1.01);

  C.insert(-1.00);

  C.insert( 1.00);

  C.insert( 1.01);


  const double E_b[] = {0.01, 0.2, 0.4, 0.6, 0.8, 1.0, 1.5, 2.0, 3.0, 4.0, 5.0, 6.0, 8.0, 10.0, 12.0, 14.0, 16.0, 18.0, 20.0, 25.0, 30.0, 40.0, 55.0, 70.0, 85.0, 100.0};


  const double R[]   = {0.0, 1.0, 2.5, 5.0, 7.5, 10.0, 12.5, 15.0, 17.5, 20.0, 25.0, 30.0, 35.0, 40.0, 50.0, 65.0, 80.0};

  const double Dmax_m = R[sizeof(R)/sizeof(R[0]) - 1]; //Last element of R


  for (int iE = 0; iE != sizeof(E_b)/sizeof(E_b[0]); ++iE) {


    const double E_m = E_b[iE];


    for (int i = 0; i != sizeof(R)/sizeof(R[0]); ++i) {


      const double R_m = R[i];


      for (set<double>::const_iterator c = C.begin(); c != C.end(); ++c) {


        const double cd = *c;


        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* p : { &h0, &h1 }) {

              (*p)[E_m][R_m][cd][theta][phi];

            }

          }

        }

      }

    }

  }


  double buffer[] = {-15.0,-10.0,-7.5,-5,-4,-3,-2,-1.5,-1.0,-0.5,-0.1, 0.0,0.1, 0.5, 1.0, 1.5, 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 };

  for (JMultiHistogram_t::super_iterator

         i1 = h1.super_begin(); i1 != h1.super_end(); ++i1) {

    for (int j = 0; j != sizeof(buffer)/sizeof(buffer[0]); ++j) {

      i1.getValue()[buffer[j]];

    }

  }//Buffer defines the binning of the inner histograms (dt)


  long long int h0_fillcount = 0;

  long long int h1_fillcount = 0;


  while (inputFile.hasNext()) {


    //NOTICE("event: " << setw(10) << inputFile.getCounter() << '\r'); STATUS(endl);

    //I suspect the printing above slows down the program considerably


    const Evt* event = inputFile.next();

    const vector<Trk>* mc_tracks = &(event->mc_trks);


    vector<int> hit_remap;


    //*************************************************

    //Fixes 2016 atmospheric neutrino simulation hit-making neutrino bug

    //Remove after fix.

    if(fix_bug == 0){

      for (vector<Trk>::const_iterator i = mc_tracks->begin(); i != mc_tracks->end(); ++i) {

        hit_remap.push_back(i->id);

      }

    } else if(fix_bug == 1){

      hit_remap = getHitRemapping(mc_tracks);

    } else if(fix_bug == 2){

      bool skipevent = false;

      for (vector<Trk>::const_iterator i = mc_tracks->begin(); i != mc_tracks->end(); ++i) {

        if(i->type == 14){

          skipevent = true;

          break;

        }

        hit_remap.push_back(i->id);

      }

      if(skipevent) continue;

    }

    //*************************************************


    double t0 = (*mc_tracks)[1].t;


    for (vector<Hit>::const_iterator hit = event->mc_hits.begin(); hit != event->mc_hits.end(); ++hit)  {


      try {


        if(hit->origin >= (int)(*mc_tracks).size()){

          std::out_of_range err("Hit origin not in tracklist. Avoided segfault");

          throw err;

        }


        Int_t corr_origin  = hit_remap[hit->origin];

        Trk curr_track     = (*mc_tracks)[corr_origin]; //Get track from corrected origin

        JDirection3D dir = getDirection(curr_track);    const JRotation3D R(dir);

        JPosition3D pos = getPosition(curr_track);                 pos.rotate(R);

        JAxis3D axis = pmtRouter.getPMT(hit->pmt_id);

        axis.transform(R, pos);


        const double E     = curr_track.E;

        const double t1    = t0 + axis.getLength() * getInverseSpeedOfLight() * getIndexOfRefraction();

        const double D_m   = axis.getLength();

        const double cd    = axis.getZ()/D_m;

        const double theta = axis.getTheta();

        const double phi   = fabs(axis.getPhi());

        const double dt    = getTime(*hit) - t1;

        const double npe   = getNPE (*hit);


        if(D_m < Dmax_m){

          h1.fill(E, D_m, cd, theta, phi, dt, npe);

          h1_fillcount += 1;

        }

      }

      catch(const exception& error) {

        std::cout << "Fatal error for event: " << inputFile.getCounter() << std::endl;

        FATAL(error.what() << endl);

      }

    }


    //                               ignore primary neutrino ->    ^

    for (vector<Trk>::const_iterator tr = event->mc_trks.begin() + 1; tr != event->mc_trks.end(); ++tr) {


      if(find(particles.begin(), particles.end(), tr->type) == particles.end()){

        continue;

      }


      //cout << "Particle accepted: " << tr->type << endl;


      JDirection3D dir = getDirection(*tr);    const JRotation3D R(dir);

      JPosition3D pos = getPosition(*tr);       pos.rotate(R);


      for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {


        JPosition3D P = module->getPosition();


        P.rotate(R);


        P.sub(pos);


        if (P.getLength() < Dmax_m) {


          for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {


            JAxis3D axis = *pmt;

            axis.transform(R, pos);


            h0.fill(tr->E, axis.getLength(), axis.getZ()/axis.getLength(), axis.getTheta(), fabs(axis.getPhi()), 0.0, 1.0);

            h0_fillcount += 1;

          }

        }

      }

    }


    if(h1_fillcount > h0_fillcount){

      std::cout << "WARNING: recorded hits in normalization histogram that were not recorded in normalization histogram. This should not happen." << std::endl;

      std::cout << "h1_fillcount: " << h1_fillcount << ", h0_fillcount: " << h0_fillcount << std::endl;

    }

  };


  double integral = 0;

  for (JMultiHistogram_t::super_iterator i = h0.super_begin(); i != h0.super_end(); ++i) {

    integral+=i.getValue().getIntegral();

  }

  DEBUG("Integral:\t" << integral << endl);


  // output


  JFileStreamWriter out(outputFile.c_str());


  NOTICE("Storing " << outputFile << ", " << flush);


  for (const JMultiHistogram_t* p : { &h0, &h1 }) {

    out.store(*p);

  }


  out.close();

  NOTICE("JHistHDE done. " << endl);

  return 1;

}


Antares.hh
Properties of Antares PMT and deep-sea water.

Evt.hh

Head.hh

Hit.hh

JAAnetToolkit.hh
Definition of hit and track types and auxiliary methods for handling Monte Carlo data.

JAxis3D.hh

outputFile
string outputFile
Definition JDAQTimesliceSelector.cc:37

JDetectorToolkit.hh

JDetector.hh
Data structure for detector geometry and calibration.

JDirection3D.hh

JDispersion.hh

JFileStreamIO.hh

JFunction1D_t.hh

JFunctionalMap_t.hh
Various implementations of functional maps.

JGeometry3DToolkit.hh

JHead.hh

main
int main(int argc, char **argv)
Definition JHistHDE.cc:76

getHitRemapping
std::vector< Int_t > getHitRemapping(const std::vector< Trk > *tracklist)
Definition JHistHDE.cc:55

JHistogram1D_t.hh

JHistogramMap_t.hh

JMessage.hh
General purpose messaging.

DEBUG
#define DEBUG(A)
Message macros.
Definition JMessage.hh:62

NOTICE
#define NOTICE(A)
Definition JMessage.hh:64

FATAL
#define FATAL(A)
Definition JMessage.hh:67

debug
int debug
debug level
Definition JSirene.cc:72

JMonteCarloFileSupportkit.hh

JMultipleFileScanner.hh
Scanning of objects from multiple files according a format that follows from the extension of each fi...

JPDFToolkit.hh
Auxiliary methods for PDF calculations.

JPDFTransformer.hh

JPDF.hh

JPMTRouter.hh
Direct access to PMT in detector data structure.

JParser.hh
Utility class to parse command line options.

make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2142

JParticleTypes.hh
Definition of particle types.

JConstants.hh
Physics constants.

JPosition3D.hh

JSupport.hh
ROOT TTree parameter settings of various packages.

JTransformableMultiHistogram.hh

KM3NeT.hh
Properties of KM3NeT PMT and deep-sea water.

Trk.hh

JAANET::JHead
Monte Carlo run header.
Definition JHead.hh:1236

JDETECTOR::JDetector
Detector data structure.
Definition JDetector.hh:96

JDETECTOR::JPMTRouter
Router for direct addressing of PMT data in detector data structure.
Definition JPMTRouter.hh:37

JGEOMETRY3D::JAxis3D
Axis object.
Definition JAxis3D.hh:41

JGEOMETRY3D::JAxis3D::transform
void transform(const JAxis3D &axis)
Transform axis to reference frame of given axis.
Definition JAxis3D.hh:359

JGEOMETRY3D::JDirection3D
Data structure for direction in three dimensions.
Definition JDirection3D.hh:35

JGEOMETRY3D::JPosition3D
Data structure for position in three dimensions.
Definition JPosition3D.hh:38

JGEOMETRY3D::JPosition3D::rotate
JPosition3D & rotate(const JRotation3D &R)
Rotate.
Definition JPosition3D.hh:186

JGEOMETRY3D::JRotation3D
Rotation matrix.
Definition JRotation3D.hh:114

JGEOMETRY3D::JVector3D::getLength
double getLength() const
Get length.
Definition JVector3D.hh:246

JGEOMETRY3D::JVector3D::getZ
double getZ() const
Get z position.
Definition JVector3D.hh:115

JGEOMETRY3D::JVector3D::sub
JVector3D & sub(const JVector3D &vector)
Subtract vector.
Definition JVector3D.hh:158

JGEOMETRY3D::JVersor3D::getTheta
double getTheta() const
Get theta angle.
Definition JVersor3D.hh:128

JGEOMETRY3D::JVersor3D::getPhi
double getPhi() const
Get phi angle.
Definition JVersor3D.hh:144

JIO::JFileStreamWriter
Binary buffered file output.
Definition JFileStreamIO.hh:96

JIO::JFileStreamWriter::close
void close()
Close file.
Definition JFileStreamIO.hh:142

JIO::JWriter::store
JWriter & store(const JSerialisable &object)
Write object.
Definition JSerialisable.hh:176

JLANG::JException
General exception.
Definition JException.hh:24

JOSCPROB::JOscProbInterpolator
Template definition of a multi-dimensional oscillation probability interpolation table.
Definition JOscProbInterpolator.hh:71

JOSCPROB::JOscProbInterpolator::read
JReader & read(JReader &in) override final
Read from input.
Definition JOscProbInterpolator.hh:258

JPHYSICS::JDispersion
Implementation of dispersion for water in deep sea.
Definition JDispersion.hh:28

JPHYSICS::JGeant
Function object for the probability density function of photon emission from EM-shower as a function ...
Definition JGeant.hh:32

JPHYSICS::JGeanx
Probability density function of photon emission from EM-shower as a function of cosine of the emissio...
Definition JGeanx.hh:32

JTOOLS::JCollection< JElement2D_t >::const_iterator
container_type::const_iterator const_iterator
Definition JCollection.hh:91

JTOOLS::JHistogram1D::abscissa_type
collection_type::abscissa_type abscissa_type
Definition JHistogram1D.hh:123

JTOOLS::JQuadrature
Type definition for numerical integration.
Definition JQuadrature.hh:39

JTOOLS::JTransformableMultiHistogram
Transformable multidimensional histogram.
Definition JTransformableMultiHistogram.hh:39

std::set
Definition JSTDTypes.hh:16

std::vector
Definition JSTDTypes.hh:15

JAANET::getDirection
JDirection3D getDirection(const Vec &dir)
Get direction.
Definition JAAnetToolkit.hh:218

JAANET::getPosition
JPosition3D getPosition(const Vec &pos)
Get position.
Definition JAAnetToolkit.hh:182

JAANET::getNPE
double getNPE(const Hit &hit)
Get true charge of hit.
Definition JAAnetToolkit.hh:112

JAANET::getOffset
Vec getOffset(const JHead &header)
Get offset.
Definition JHeadToolkit.hh:291

JAANET::TRACK_TYPE_NUMU
@ TRACK_TYPE_NUMU
Definition JParticleTypes.hh:99

JDETECTOR::load
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
Definition JDetectorToolkit.hh:520

JLANG::to_string
std::string to_string(const T &value)
Convert value to string.
Definition JLangToolkit.hh:207

JPHYSICS::getIndexOfRefraction
double getIndexOfRefraction()
Get average index of refraction of water corresponding to group velocity.
Definition JPhysics/JConstants.hh:134

JPHYSICS::getMinimalWavelength
double getMinimalWavelength()
Get minimal wavelength for PDF evaluations.
Definition JPDFToolkit.hh:27

JPHYSICS::getInverseSpeedOfLight
const double getInverseSpeedOfLight()
Get inverse speed of light.
Definition JPhysics/JConstants.hh:123

JPHYSICS::getMaximalWavelength
double getMaximalWavelength()
Get maximal wavelength for PDF evaluations.
Definition JPDFToolkit.hh:38

JPHYSICS::C
static const double C
Physics constants.
Definition JPhysics/JConstants.hh:25

JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition JAAnetToolkit.hh:45

JSUPPORT::getHeader
Head getHeader(const JMultipleFileScanner_t &file_list)
Get Monte Carlo header.
Definition JMonteCarloFileSupportkit.hh:547

JSYSTEM::getTime
const char * getTime()
Get current local time conform ISO-8601 standard.
Definition JDateAndTime.hh:56

JTOOLS::JHistogram1D_t
JHistogram1D< JElement2D< double, double >, JCollection > JHistogram1D_t
Type definition of a 1 dimensional histogram based on a JCollection.
Definition JHistogram1D_t.hh:22

JTOOLS::j
int j
Definition JPolint.hh:801

std
Definition JSTDTypes.hh:14

Evt
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
Definition Evt.hh:21

JAANET::detector
Detector file.
Definition JHead.hh:227

JTOOLS::JHistogramGridMap_t
Type definition of a JHistogramMap based on JGridMap implementation.
Definition JHistogramMap_t.hh:36

JTOOLS::JHistogramMap_t
Type definition of a JHistogramMap based on JMap implementation.
Definition JHistogramMap_t.hh:27

JTOOLS::JMAPLIST
Auxiliary class for recursive map list generation.
Definition JMapList.hh:109

Trk
The Trk class represents a Monte Carlo (MC) particle as well as a reconstructed track/shower.
Definition Trk.hh:15

Vec
The Vec class is a straightforward 3-d vector, which also works in pyroot.
Definition Vec.hh:13