Jpp/v18.6.0-rc.1/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();

 }

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

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

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

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

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

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

main
int main(int argc, char *argv[])
Definition: Main.cc:15

JSupport.hh
ROOT TTree parameter settings of various packages.

JPDFToolkit.hh
Auxiliary methods for PDF calculations.

JFileStreamIO.hh

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

JAANET::HIT_TYPE_MUON_SCATTERED
Scattered light from muon.
Definition: JAAnetToolkit.hh:89

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

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

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

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

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

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

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

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

JAANET::HIT_TYPE_BREMSSTRAHLUNG_SCATTERED
Scattered light from Bremsstrahlung.
Definition: JAAnetToolkit.hh:93

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

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

JAxis3D.hh

JAANET::HIT_TYPE_BREMSSTRAHLUNG_DIRECT
Direct light from Bremsstrahlung.
Definition: JAAnetToolkit.hh:92

JAANET::getTime
double getTime(const Hit &hit)
Get true time of hit.
Definition: JAAnetToolkit.hh:107

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

outputFile
string outputFile
Definition: JDAQTimesliceSelector.cc:37

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

KM3NET
#define KM3NET

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

JDetector.hh
Data structure for detector geometry and calibration.

std::vector
Definition: JSTDTypes.hh:13

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

JTransformableMultiHistogram.hh

JFunctionalMap_t.hh
Various implementations of functional maps.

JVisibleEnergyToolkit.hh
Auxiliary methods for evaluating visible energies.

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

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

sin
then usage set_variable ACOUSTICS_WORKDIR $WORKDIR set_variable FORMULA sin([0]+2 *$PI *([1]+[2]*x)*x)" set_variable DY 1.0e-8 mkdir $WORKDIR for DETECTOR in $DETECTORS[*]

JPosition3D.hh

i
then rm i
Definition: JEvtReweightMupageParameterScan.sh:309

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

JFunction1D_t.hh

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

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

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

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

JGEOMETRY3D::JCylinder3D
Cylinder object.
Definition: JCylinder3D.hh:39

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

JDirection3D.hh

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

ID
Definition: raw_data_converter.cpp:43

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

JGeometry3DToolkit.hh

makedeclinationtable.x
tuple x
Definition: makedeclinationtable.py:44

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

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

JHistogram1D_t.hh

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

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

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

ANTARES::getAmbientPressure
double getAmbientPressure()
Get ambient pressure.
Definition: Antares.hh:40

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

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

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

JConstants.hh
Physics constants.

JMATH::PI
static const double PI
Mathematical constants.
Definition: JMath/JConstants.hh:20

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

JGEOMETRY3D::JVector3D::getY
double getY() const
Get y position.
Definition: JVector3D.hh:104

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

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

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

JTOOLS::JRange< double >

JMessage.hh
General purpose messaging.

JSUPPORT::JMultipleFileScanner< Head >
Template specialisation of JMultipleFileScanner for Monte Carlo header.
Definition: JMonteCarloFileSupportkit.hh:414

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

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

JMonteCarloFileSupportkit.hh

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

R
then JCookie sh JDataQuality D $DETECTOR_ID R
Definition: JDataQuality.sh:41

Hit
Definition: Hit.hh:8

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

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

JAANET::HIT_TYPE_MUON_DIRECT
Direct light from muon.
Definition: JAAnetToolkit.hh:88

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

Evt.hh

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

JParser.hh
Utility class to parse command line options.

JDetectorToolkit.hh

JGEOMETRY3D::JVector3D::getX
double getX() const
Get x position.
Definition: JVector3D.hh:94

JDispersion.hh

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

JPHYSICS::getTanThetaC
double getTanThetaC()
Get average tangent of Cherenkov angle of water corresponding to group velocity.
Definition: JPhysics/JConstants.hh:156

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

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

Z
do set_variable MODULE getModule a $WORKDIR detector_a datx L $STRING JEditDetector a $WORKDIR detector_a datx M $MODULE setz o $WORKDIR detector_a datx JEditDetector a $WORKDIR detector_b datx M $MODULE setz o $WORKDIR detector_b datx done echo Output stored at $WORKDIR detector_a datx and $WORKDIR tripod_a txt JDrawDetector2D a $WORKDIR detector_a datx a $WORKDIR detector_b datx L BL o detector $FORMAT $BATCH JDrawDetector2D T $WORKDIR tripod_a txt T $WORKDIR tripod_b txt L BL o tripod $FORMAT $BATCH JCompareDetector a $WORKDIR detector_a datx b $WORKDIR detector_b datx o $WORKDIR abc root &dev null for KEY in X Y Z
Definition: JAlignDetector.sh:81

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

JSUPPORT::JLimit::getLimit
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:84

detector
do set_variable DETECTOR_TXT $WORKDIR detector
Definition: detector-XY:fit1d.sh:38

Hit.hh

Head.hh

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

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

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

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

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

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

debug
int debug
debug level
Definition: archive-put-wiki-detectors.sh:92

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

P
then $DIR JPlotNPE PDG P
Definition: JPlotNPE-PDG.sh:62

JHead.hh

JHistogramMap_t.hh

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

getAngularAcceptance
double getAngularAcceptance(const double x)
Angular acceptence of PMT.
Definition: JDrawLED.cc:68

JPDFTransformer.hh

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