56 int main(
int argc,
char **argv)
70 JParser<> zap(
"Program to histogram event-by-event data of shower light for making PDFs.");
74 zap[
'n'] =
make_field(numberOfEvents) = JLimit::max();
81 catch(
const exception &error) {
82 FATAL(error.what() << endl);
86 NOTICE(
"Plotting hadronic showers." << endl);
88 NOTICE(
"Plotting EM showers." << endl);
103 NOTICE(
"Apply detector offset " << center << endl);
131 JMultiHistogram_t h0;
132 JMultiHistogram_t
h1;
140 for (JQuadrature::const_iterator i = qeant.begin(); i != qeant.end(); ++i) {
149 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};
150 const double Dmax_m = 80.0;
152 for (
int i = 0; i !=
sizeof(
R)/
sizeof(
R[0]); ++i) {
154 const double R_m =
R[i];
158 const double cd = *
c;
160 const double grid = 10.0 + 0.0 * R_m/100.0;
161 const double alpha = 2.0 * sqrt(1.0 - cos(grid *
PI / 180.0));
163 const int number_of_theta_points = max(2, (
int) (180.0/(1.4 * grid)));
164 const double theta_step =
PI / (number_of_theta_points + 1);
166 for (
double theta = -0.5*theta_step; theta <
PI + theta_step; theta += theta_step) {
168 const int number_of_phi_points = max(2, (
int) (
PI * sin(theta) / alpha));
169 const double phi_step =
PI / (number_of_phi_points + 1);
171 for (
double phi = -0.5*phi_step; phi <
PI + phi_step; phi += phi_step) {
173 for (JMultiHistogram_t* p : { &h0, &
h1 }) {
174 (*p)[R_m][cd][theta][phi];
181 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 };
183 for (JMultiHistogram_t::super_iterator i1 =
h1.super_begin(); i1 !=
h1.super_end(); ++i1) {
184 for (
int j = 0;
j !=
sizeof(buffer)/
sizeof(buffer[0]); ++
j) {
185 i1.getValue()[buffer[
j]];
189 while (inputFile.hasNext()) {
191 NOTICE(
"event: " << setw(10) << inputFile.getCounter() <<
'\r');
STATUS(endl);
193 const Evt*
event = inputFile.next();
196 WARNING(
"Event " << inputFile.getCounter() <<
" does not correspond to a neutrino interaction; skip.");
201 WARNING(
"No electron/hadrons found; skip.");
218 JAxis3D axis = pmtRouter.getPMT(hit->pmt_id);
226 const double cd = axis.
getZ() / D_m;
227 const double theta = axis.
getTheta();
228 const double phi = fabs(axis.
getPhi());
229 const double dt =
getTime(*hit) - t1;
230 const double npe =
getNPE (*hit);
233 h1.fill(D_m, cd, theta, phi, dt, npe/Evis.
len());
237 catch(
const exception& error) {
238 FATAL(error.what() << endl);
242 for (JDetector::const_iterator module =
detector.begin(); module !=
detector.end(); ++module) {
252 for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
264 for (JMultiHistogram_t::super_iterator i = h0.super_begin(); i != h0.super_end(); ++i) {
265 integral+=i.getValue().getIntegral();
267 DEBUG(
"Integral:\t" << integral << endl);
273 NOTICE(
"Storing, " << flush);
275 for (
const JMultiHistogram_t* p : { &h0, &
h1 }) {
Router for direct addressing of PMT data in detector data structure.
JVertex3D getVertex(const Trk &track)
Get vertex.
Utility class to parse command line options.
int main(int argc, char *argv[])
ROOT TTree parameter settings of various packages.
double getIndexOfRefraction()
Get average index of refraction of water corresponding to group velocity.
bool has_neutrino(const Evt &evt)
Test whether given event has an incoming neutrino.
Implementation of dispersion for water in deep sea.
bool has_electron(const Evt &evt)
Test whether given event has an electron.
Properties of Antares PMT and deep-sea water.
then for HISTOGRAM in h0 h1
double getMinimalWavelength()
Get minimal wavelength for PDF evaluations.
double getTime(const Hit &hit)
Get true time of hit.
double getPhi() const
Get phi angle.
Data structure for detector geometry and calibration.
void transform(const JAxis3D &axis)
Transform axis to reference frame of given axis.
Various implementations of functional maps.
double len() const
Get length.
Head getHeader(const JMultipleFileScanner_t &file_list)
Get Monte Carlo header.
static const double C
Physics constants.
Auxiliary class for defining the range of iterations of objects.
JVector3D & sub(const JVector3D &vector)
Subtract vector.
Properties of KM3NeT PMT and deep-sea water.
The Vec class is a straightforward 3-d vector, which also works in pyroot.
double getMaximalWavelength()
Get maximal wavelength for PDF evaluations.
JDirection3D getDirection(const Vec &dir)
Get direction.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
double getTheta() const
Get theta angle.
double getAmbientPressure()
Get ambient pressure.
static const double PI
Mathematical constants.
Direct access to PMT in detector data structure.
const JPosition3D & getPosition() const
Get position.
double getLength() const
Get length.
then usage $script[distance] fi case set_variable R
General purpose messaging.
bool from_hadron(const Hit &hit)
Test whether given hit was produced by a hadronic shower.
Scanning of objects from multiple files according a format that follows from the extension of each fi...
Probability density function of photon emission from EM-shower as a function of cosine of the emissio...
virtual double getIndexOfRefractionGroup(const double lambda) const
Index of refraction for group velocity.
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
General purpose class for object reading from a list of file names.
Utility class to parse command line options.
const double getInverseSpeedOfLight()
Get inverse speed of light.
double getNPE(const Hit &hit)
Get true charge of hit.
Data structure for position in three dimensions.
const JLimit & getLimit() const
Get limit.
$WORKDIR ev_configure_domsimulator txt echo process $DOM_SIMULATOR $i $SOURCE_HOST[$index] csh c(setenv ROOTSYS $ROOTSYS &&source $JPP_DIR/setenv.csh $JPP_DIR &&($DOM_SIMULATOR\-u\$NAME\$\-H\$SERVER\$\-M\$LOGGER\$\-d $DEBUG</dev/null > &/dev/null &))'
Function object for the probability density function of photon emission from EM-shower as a function ...
const Trk & get_neutrino(const Evt &evt)
Get incoming neutrino.
do set_variable DETECTOR_TXT $WORKDIR detector
Vec getVisibleEnergy(const Trk &track, const JCylinder3D &can)
Get the visible energy of a track.
bool from_electron(const Hit &hit)
Test whether given hit was produced by an electron.
double getT(const JVector3D &pos) const
Get arrival time of Cherenkov light at given position.
The Trk class represents a Monte Carlo (MC) particle as well as a reconstructed track/shower.
JPosition3D & rotate(const JRotation3D &R)
Rotate.
Binary buffered file output.
double getZ() const
Get z position.
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
int count_hadrons(const Evt &evt)
Count the number of hadrons in a given event (not including neutral pions)
#define DEBUG(A)
Message macros.
The cylinder used for photon tracking.
double getAngularAcceptance(const double x)
Angular acceptence of PMT.