Jpp/test-rotations-old/JHistHDF_8cc_source.html

#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/JConstants.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/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 "JSirene/JVisibleEnergyToolkit.hh"


#include "Jeep/JParser.hh"

#include "Jeep/JMessage.hh"


namespace {


  using namespace JPP;


  /**

   * Function for classifying hits generated with KM3 or JSirene.

   *

   * \param  track          muon track

   * \param  hit            hit

   * \return                1 for muon; 2 for shower; -1 for anything else

   */

  inline int classify(const Trk& track, const Hit& hit)

  {

    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

   */

  inline int classify_km3sim(const Trk& track, const Hit& hit)

  {

    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;

  using namespace KM3NET;


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

  }


  JDetector detector;


  try {

    load(detectorFile, detector);

  }

  catch(const JException& error) {

    FATAL(error);

  }


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


  const bool KM3Sim = is_km3sim(head);

  const Vec  offset = getOffset(head);


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


  detector -= getPosition(offset);


  const JCylinder3D can(detector.begin(), detector.end());


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

          (*p)[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 x : { 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} ) {

      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


        const double Evis       = getVisibleEnergy      (*event, can);

        const Vec    EvisVector = getVisibleEnergyVector(*event, can);


        JVertex3D vertex = getVertex(*track);


        const JRotation3D R(getDirection(EvisVector));


        vertex.rotate(R);


        JRange<double> Z(0.0, gWater(track->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, vertex);


            const double t1    = vertex.getT() + (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);


            const int    ID    = (KM3Sim ? classify_km3sim(*track, *hit) : classify(*track, *hit));


            switch (ID) {


            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, Evis)); // 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(vertex);


          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, vertex);


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

    out.store(*p);

  }


  out.close();

}


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 JHistHDF.cc:102

JHistogram1D_t.hh

JHistogramMap_t.hh

JMessage.hh
General purpose messaging.

STATUS
#define STATUS(A)
Definition JMessage.hh:63

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

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

JConstants.hh
Physics constants.

JPosition3D.hh

JSupport.hh
ROOT TTree parameter settings of various packages.

JTransformableMultiHistogram.hh

JVisibleEnergyToolkit.hh
Auxiliary methods for evaluating visible energies.

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

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

JDETECTOR::JPMTRouter::getPMT
const JPMT & getPMT(const JPMTAddress &address) const
Get PMT.
Definition JPMTRouter.hh:92

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::JCylinder3D
Cylinder object.
Definition JCylinder3D.hh:41

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::getY
double getY() const
Get y position.
Definition JVector3D.hh:104

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::JVector3D::getX
double getX() const
Get x position.
Definition JVector3D.hh:94

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

JGEOMETRY3D::JVertex3D
3D vertex.
Definition JVertex3D.hh:36

JGEOMETRY3D::JVertex3D::getT
double getT(const JVector3D &pos) const
Get arrival time of Cherenkov light at given position.
Definition JVertex3D.hh:147

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:173

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

JPARSER::JParser
Utility class to parse command line options.
Definition JParser.hh:1698

JPHYSICS::JDispersionInterface::getKmin
double getKmin(const double lambda) const
Get smallest index of refraction for Bremsstrahlung light (i.e. point at which dt/dz = 0).
Definition JDispersionInterface.hh:91

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

JPHYSICS::JPDFTransformer
Template definition of transformer of the probability density function (PDF) of the time response of ...
Definition JPDFTransformer.hh:620

JSUPPORT::JMultipleFileScanner
General purpose class for object reading from a list of file names.
Definition JMultipleFileScanner.hh:306

JSUPPORT::JMultipleFileScanner::hasNext
virtual bool hasNext() override
Check availability of next element.
Definition JMultipleFileScanner.hh:409

JSUPPORT::JMultipleFileScanner::getCounter
counter_type getCounter() const
Get counter.
Definition JMultipleFileScanner.hh:382

JSUPPORT::JMultipleFileScanner::next
virtual const pointer_type & next() override
Get next element.
Definition JMultipleFileScanner.hh:457

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

JTOOLS::JRange
Range of values.
Definition JRange.hh:42

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

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

JAANET::is_km3sim
bool is_km3sim(const JHead &header)
Check for detector simulation.
Definition JHeadToolkit.hh:157

JAANET::from_muon
bool from_muon(const Hit &hit)
Test whether given hit was produced by a muon.
Definition JAAnetToolkit.hh:793

JAANET::getVertex
JVertex3D getVertex(const Trk &track)
Get vertex.
Definition JAAnetToolkit.hh:654

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

JAANET::is_muon
bool is_muon(const Trk &track)
Test whether given track is a (anti-)muon.
Definition JAAnetToolkit.hh:338

JAANET::HIT_TYPE_MUON_DIRECT
@ HIT_TYPE_MUON_DIRECT
Direct light from muon.
Definition JAAnetToolkit.hh:77

JAANET::HIT_TYPE_MUON_SCATTERED
@ HIT_TYPE_MUON_SCATTERED
Scattered light from muon.
Definition JAAnetToolkit.hh:78

JAANET::HIT_TYPE_BREMSSTRAHLUNG_SCATTERED
@ HIT_TYPE_BREMSSTRAHLUNG_SCATTERED
Scattered light from Bremsstrahlung.
Definition JAAnetToolkit.hh:82

JAANET::HIT_TYPE_BREMSSTRAHLUNG_DIRECT
@ HIT_TYPE_BREMSSTRAHLUNG_DIRECT
Direct light from Bremsstrahlung.
Definition JAAnetToolkit.hh:81

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

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

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

JPHYSICS::gWater
static const JGeaneWater gWater
Function object for energy loss of muon in sea water.
Definition JGeane.hh:381

JPHYSICS::getSpeedOfLight
const double getSpeedOfLight()
Get speed of light.
Definition JPhysics/JConstants.hh:112

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

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

JSIRENE::getVisibleEnergy
double getVisibleEnergy(const Trk &, const JCylinder3D &)
Get the visible energy of a track.
Definition JVisibleEnergyToolkit.hh:73

JSIRENE::getVisibleEnergyVector
Vec getVisibleEnergyVector(const Trk &track, const JCylinder3D &can=getMaximumContainmentVolume())
Get the visible energy vector of a track.
Definition JVisibleEnergyToolkit.hh:128

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:55

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

KM3NET
Name space for KM3NeT.
Definition Jpp.hh:16

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

Hit
Definition Hit.hh:10

Hit::pmt_id
int pmt_id
global PMT identifier as found in evt files
Definition Hit.hh:20

Hit::origin
int origin
track id of the track that created this hit (mc only)
Definition Hit.hh:29

Hit::type
int type
particle type or parametrisation used for hit (mc only)
Definition Hit.hh:28

ID
Definition raw_data_converter.cpp:43

JAANET::can
The cylinder used for photon tracking.
Definition JHead.hh:575

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

JSUPPORT::JLimit
Auxiliary class for defining the range of iterations of objects.
Definition JLimit.hh:45

JSUPPORT::JLimit::max
static counter_type max()
Get maximum counter value.
Definition JLimit.hh:128

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

Trk::id
int id
track identifier
Definition Trk.hh:16

Trk::E
double E
Energy [GeV] (either MC truth or reconstructed)
Definition Trk.hh:20

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