Jpp/f0c5bc59/JHistHDG_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/JDispersion.hh"

 #include "JPhysics/JConstants.hh"

 #include "JPhysics/Antares.hh"

 #include "JPhysics/KM3NeT.hh"


 #include "JIO/JFileStreamIO.hh"


 #include "JGeometry3D/JGeometry3DToolkit.hh"

 #include "JGeometry3D/JDirection3D.hh"

 #include "JGeometry3D/JCylinder3D.hh"

 #include "JGeometry3D/JPosition3D.hh"

 #include "JGeometry3D/JAxis3D.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 "JSirene/JVisibleEnergyToolkit.hh"


 /**

  * \file

  *

  * Program to histogram event-by-event data of shower light for making PDFs.

  * \author lquinn

  */

 int main(int argc, char **argv)

 {

   using namespace std;

   using namespace JPP;


   JMultipleFileScanner<Evt> inputFile;

   JLimit_t&   numberOfEvents = inputFile.getLimit();

   string      outputFile;

   string      detectorFile;

   bool        hadronicMode;

   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['H'] = make_field(hadronicMode);

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


     zap(argc, argv);

   }

   catch(const exception &error) {

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

   }


   if(hadronicMode) {

     NOTICE("Plotting hadronic showers." << endl);

   } else {

     NOTICE("Plotting EM showers." << endl);

   }


   JDetector detector;


   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 JCylinder3D can(detector.begin(), detector.end());


   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,

                                                 JHistogramGridMap_t,

                                                 JHistogramGridMap_t>::maplist> JMultiHistogram_t;


   typedef JPDFTransformer<4, abscissa_type>                         JFunction4DTransformer_t;


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

   JMultiHistogram_t h1;   // light from cascade


   h1.transformer.reset(new JFunction4DTransformer_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, 20, 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 R[] = {0.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, 60.0, 70.0, 80.0};

   const double Dmax_m = 80.0;


   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)[R_m][cd][theta][phi];

           }

         }

       }

     }

   }


   double buffer[] = { 0.0, 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]];

     }

   }


   while (inputFile.hasNext()) {


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


     const Evt* event = inputFile.next();


     if (!has_neutrino(*event)) {

       WARNING("Event " << inputFile.getCounter() << " does not correspond to a neutrino interaction; skip.");

       continue;

     }


     if (!has_electron(*event) || count_hadrons(*event) == 0) {

       WARNING("No electron/hadrons found; skip.");

       continue;

     }


     const Trk& neutrino = get_neutrino(*event);


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

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


     JVertex3D vertex = getVertex(neutrino);


     const JRotation3D R(getDirection(EvisVector));


     vertex.rotate(R);


     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.getPosition());


         if ((!hadronicMode && from_electron(*hit)) || (hadronicMode && from_hadron(*hit))) {


           const double t1    = vertex.getT() + 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(D_m, cd, theta, phi, dt, npe/Evis);

           }

         }

       }

       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 (P.getLength() < h0.getXmax()) {


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


           JAxis3D axis = *pmt;


           axis.transform(R, vertex);


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

         }

       }

     }

   }


   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, " << flush);


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

     out.store(*p);

   }


   out.close();

   NOTICE("done." << endl);

 }

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

JCylinder3D.hh

outputFile
string outputFile
Definition: JDAQTimesliceSelector.cc:37

JDetectorToolkit.hh

JDetector.hh
Data structure for detector geometry and calibration.

JDirection3D.hh

JDispersion.hh

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

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: JHistHDG.cc:56

JHistogram1D_t.hh

JHistogramMap_t.hh

JMessage.hh
General purpose messaging.

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

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

WARNING
#define WARNING(A)
Definition: JMessage.hh:65

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

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::JPosition3D::getPosition
const JPosition3D & getPosition() const
Get position.
Definition: JPosition3D.hh:130

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

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::getIndexOfRefractionGroup
virtual double getIndexOfRefractionGroup(const double lambda) const
Index of refraction for group velocity.
Definition: JDispersionInterface.hh:52

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

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

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

JSUPPORT::JMultipleFileScanner< Evt >

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::JQuadrature
Type definition for numerical integration.
Definition: JQuadrature.hh:39

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

std::set
Definition: JSTDTypes.hh:16

std::vector
Definition: JSTDTypes.hh:15

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

JAANET::get_neutrino
const Trk & get_neutrino(const Evt &evt)
Get incoming neutrino.
Definition: JAAnetToolkit.hh:513

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

JAANET::from_hadron
bool from_hadron(const Hit &hit)
Test whether given hit was produced by a hadronic shower.
Definition: JAAnetToolkit.hh:934

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

JAANET::has_neutrino
bool has_neutrino(const Evt &evt)
Test whether given event has an incoming neutrino.
Definition: JAAnetToolkit.hh:491

JAANET::has_electron
bool has_electron(const Evt &evt)
Test whether given event has an electron.
Definition: JAAnetToolkit.hh:694

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

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

JAANET::from_electron
bool from_electron(const Hit &hit)
Test whether given hit was produced by an electron.
Definition: JAAnetToolkit.hh:722

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

JAANET::count_hadrons
int count_hadrons(const Evt &evt)
Count the number of hadrons in a given event (not including neutral pions)
Definition: JAAnetToolkit.hh:920

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

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

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

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

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

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

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

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

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

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

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