51 template<
class JMultiHistogram_t>
64 detector -=
JCenter3D(detector.begin(), detector.end());
67 detector -= get<JPosition3D>(head);
78 int operator()(
const Trk& track,
const Hit& hit)
const
81 return classify_km3sim(track, hit);
83 return classify_km3(track, hit);
98 inline double getEnergy(
const Evt& event,
const Trk& track,
const JAxis3D& axis)
const
118 int classify_km3(
const Trk& track,
const Hit& hit)
const
146 int classify_km3sim(
const Trk& track,
const Hit& hit)
const
163 int main(
int argc,
char **argv)
176 JParser<> zap(
"Program to histogram event-by-event data of muon light for making PDFs.");
180 zap[
'n'] =
make_field(numberOfEvents) = JLimit::max();
186 catch(
const exception &error) {
187 FATAL(error.what() << endl);
197 load(detectorFile, detector);
221 JMultiHistogram_t h0;
222 JMultiHistogram_t
h1;
223 JMultiHistogram_t h2;
225 const JHitClassifier<JMultiHistogram_t> classifier(head, detector);
227 const double cmin = dispersion.
getKmin (wmax);
235 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 };
237 for (
int i = 0; i !=
sizeof(
R)/
sizeof(
R[0]); ++i) {
239 const double R_m =
R[i];
241 const double grid = 10.0 + 0.0 * R_m/100.0;
242 const double alpha = 2.0 * sqrt(1.0 - cos(grid *
PI / 180.0));
244 const int number_of_theta_points = max(2, (
int) (180.0/(1.4 * grid)));
245 const double theta_step =
PI / (number_of_theta_points + 1);
247 for (
double theta = -0.5*theta_step; theta <
PI + theta_step; theta += theta_step) {
249 const int number_of_phi_points = max(2, (
int) (
PI * sin(theta) / alpha));
250 const double phi_step =
PI / (number_of_phi_points + 1);
252 for (
double phi = -0.5*phi_step; phi <
PI + phi_step; phi += phi_step) {
254 for (JMultiHistogram_t* buffer[] = { &h0, &
h1, &h2, NULL }, **histogram = buffer; *histogram != NULL; ++histogram) {
255 (**histogram)[R_m][theta][phi];
261 for (JMultiHistogram_t::super_iterator
262 i1 =
h1.super_begin(),
263 i2 = h2.super_begin(); i1 !=
h1.super_end(); ++i1, ++i2) {
265 for (
double buffer[] = { 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, -1.0 }, *x = buffer; *x >= 0.0; ++x) {
272 while (inputFile.hasNext()) {
274 NOTICE(
"event: " << setw(10) << inputFile.getCounter() <<
'\r');
STATUS(endl);
276 const Evt*
event = inputFile.next();
287 double t0 = track->
t;
303 const double E1 = classifier.getEnergy(*event, *track, axis);
306 const double R_m = axis.
getX();
307 const double theta = axis.
getTheta();
308 const double phi = fabs(axis.
getPhi());
309 const double dt =
getTime(*hit) - t1;
310 const double npe =
getNPE(*hit);
312 switch (classifier(*track, *hit)) {
315 h1.fill(R_m, theta, phi, dt, npe);
319 h2.fill(R_m, theta, phi, dt, npe/max(1.0, E1));
323 catch(
const exception& error) {
324 FATAL(error.what() << endl);
329 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
339 for (JModule::const_iterator
pmt = module->begin();
pmt != module->end(); ++
pmt) {
361 for (
const JMultiHistogram_t* buffer[] = { &h0, &
h1, &h2, NULL }, **i = buffer; *i != NULL; ++i) {
Router for direct addressing of PMT data in detector data structure.
Utility class to parse command line options.
double getKmin(const double lambda) const
Get smallest index of refraction for Bremsstrahlung light (i.e. point at which dt/dz = 0)...
JWriter & store(const JSerialisable &object)
Write object.
ROOT TTree parameter settings.
Data structure for direction in three dimensions.
double t
track time [ns] (when the particle is at pos )
Scattered light from muon.
Implementation of dispersion for water in deep sea.
bool is_muon(const Trk &track)
Test whether given track is a (anti-)muon.
Optical properties of Antares deep-sea site.
int pmt_id
global PMT identifier as found in evt files
then for HISTOGRAM in h0 h1
Scattered light from Bremsstrahlung.
const JPMT & getPMT(const JPMTAddress &address) const
Get PMT.
static const JGeaneWater gWater
Function object for energy loss of muon in sea water.
Direct light from Bremsstrahlung.
double getTime(const Hit &hit)
Get true time of hit.
double getPhi() const
Get phi angle.
double E
Energy [GeV] (either MC truth or reconstructed)
int origin
track id of the track that created this hit
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.
esac $JPP_DIR examples JDetector JTransitTime o $OUTPUT_FILE n N $NPE T $TTS_NS d $DEBUG for HISTOGRAM in tts tt2 pmt
Head getHeader(const JMultipleFileScanner_t &file_list)
Get Monte Carlo header.
Auxiliary class for defining the range of iterations of objects.
JVector3D & sub(const JVector3D &vector)
Subtract vector.
Optical properties of KM3NeT deep-sea site.
double getMaximalWavelength()
Get maximal wavelength for PDF evaluations.
Trk get_hadronic_cascade(const Evt &evt)
Auxiliary function to get average direction and total energy of a hadronic cascade.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
double getTheta() const
Get theta angle.
double getAmbientPressure()
Get ambient pressure.
double getY() const
Get y position.
void load(const JString &file_name, JDetector &detector)
Load detector from input file.
Direct access to PMT in detector data structure.
then usage $script[distance] fi case set_variable R
General purpose messaging.
Template specialisation of JMultipleFileScanner for Monte Carlo header.
Scanning of objects from multiple files according a format that follows from the extension of each fi...
General purpose class for object reading from a list of file names.
Utility class to parse command line options.
std::vector< Hit > mc_hits
MC: list of MC truth hits.
double getX() const
Get x position.
JDirection3D getDirection(const Vec &v)
Get direction.
bool is_km3sim(const JHead &header)
Check for detector simulation.
double getNPE(const Hit &hit)
Get true charge of hit.
Data structure for position in three dimensions.
const JLimit & getLimit() const
Get limit.
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.
int type
particle type or parametrisation used for hit (mc only)
double getZ() const
Get z position.
std::vector< Trk > mc_trks
MC: list of MC truth tracks.
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
then usage $script[input file[working directory[option]]] nWhere option can be E
JPosition3D getPosition(const Vec &v)
Get position.
double getAngularAcceptance(const double x)
Angular acceptence of PMT.
int main(int argc, char *argv[])
bool from_muon(const Hit &hit)
Test whether given hit was produced by a muon.