57 JMultipleFileScanner<Evt> inputFile;
66 JParser<> zap(
"Program to histogram event-by-event data of shower light for making PDFs.");
70 zap[
'n'] =
make_field(numberOfEvents) = JLimit::max();
77 catch(
const exception &error) {
78 FATAL(error.what() << endl);
82 NOTICE(
"Plotting hadronic showers." << endl);
84 NOTICE(
"Plotting EM showers." << endl);
91 load(detectorFile, detector);
93 catch(
const JException& error) {
97 const JPosition3D center = get<JPosition3D>(
getHeader(inputFile));
99 NOTICE(
"Apply detector offset " << center << endl);
103 const JPMTRouter pmtRouter(detector);
109 const JDispersion dispersion(P_atm);
111 const double ng[] = { dispersion.getIndexOfRefractionGroup(wmax),
112 dispersion.getIndexOfRefractionGroup(wmin) };
117 JMAPLIST<JHistogramMap_t,
120 JHistogramGridMap_t>::maplist> JMultiHistogram_t;
122 typedef JPDFTransformer<4, abscissa_type> JFunction4DTransformer_t;
124 JMultiHistogram_t h0;
125 JMultiHistogram_t h1;
131 JQuadrature qeant(-1.0, 1.0, 20, JGeanx(1.00, -2.2));
133 for (JQuadrature::const_iterator i = qeant.begin(); i != qeant.end(); ++i) {
142 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};
143 const double Dmax_m = 80.0;
145 for (
int i = 0; i !=
sizeof(R)/
sizeof(R[0]); ++i) {
147 const double R_m = R[i];
151 const double cd = *c;
153 const double grid = 10.0 + 0.0 * R_m/100.0;
154 const double alpha = 2.0 * sqrt(1.0 - cos(grid *
PI / 180.0));
156 const int number_of_theta_points = max(2, (
int) (180.0/(1.4 * grid)));
157 const double theta_step =
PI / (number_of_theta_points + 1);
159 for (
double theta = -0.5*theta_step; theta <
PI + theta_step; theta += theta_step) {
161 const int number_of_phi_points = max(2, (
int) (
PI * sin(theta) / alpha));
162 const double phi_step =
PI / (number_of_phi_points + 1);
164 for (
double phi = -0.5*phi_step; phi <
PI + phi_step; phi += phi_step) {
166 for (JMultiHistogram_t* buffer[] = { &h0, &h1, NULL }, **histogram = buffer; *histogram != NULL; ++histogram) {
167 (**histogram)[R_m][cd][theta][phi];
174 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 };
175 for (JMultiHistogram_t::super_iterator
176 i1 = h1.super_begin(); i1 != h1.super_end(); ++i1) {
178 for (
int j = 0; j !=
sizeof(buffer)/
sizeof(buffer[0]); ++j) {
179 i1.getValue()[buffer[j]];
183 while (inputFile.hasNext()) {
185 NOTICE(
"event: " << setw(10) << inputFile.getCounter() <<
'\r');
STATUS(endl);
187 const Evt*
event = inputFile.next();
197 WARNING(
"No electron/hadrons found. Perhaps you meant to run this program in hadronic mode." << endl);
206 const JRotation3D R(dir);
214 JAxis3D axis = pmtRouter.getPMT(hit->pmt_id);
215 axis.transform(R, pos);
221 const double D_m = axis.getLength();
222 const double cd = axis.getZ()/D_m;
223 const double theta = axis.getTheta();
224 const double phi = fabs(axis.getPhi());
225 const double dt =
getTime(*hit) - t1;
226 const double npe =
getNPE(*hit);
229 h1.fill(D_m, cd, theta, phi, dt, npe/E);
233 catch(
const exception& error) {
234 FATAL(error.what() << endl);
238 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
240 JPosition3D P = module->getPosition();
246 if (P.getLength() < h0.getXmax()) {
248 for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
251 axis.transform(R, pos);
253 h0.fill(axis.getLength(), axis.getZ()/axis.getLength(), axis.getTheta(), fabs(axis.getPhi()), 0.0, 1.0);
260 for (JMultiHistogram_t::super_iterator i = h0.super_begin(); i != h0.super_end(); ++i) {
261 integral+=i.getValue().getIntegral();
263 DEBUG(
"Integral:\t" << integral << endl);
269 NOTICE(
"Storing, " << flush);
271 for (
const JMultiHistogram_t* buffer[] = { &h0, &h1, NULL }, **i = buffer; *i != NULL; ++i) {
int count_electrons(const Evt &evt)
Count the number of electrons in a given event.
Utility class to parse command line options.
Structure to store the ToT mean and standard deviation of the hits produced by a nanobeacon in a sour...
double getMinimalWavelength()
Get minimal wavelength for PDF evaluations.
double getTime(const Hit &hit)
Get true time of hit.
Head getHeader(const JMultipleFileScanner_t &file_list)
Get Monte Carlo header.
JLimit JLimit_t
Type definition of limit.
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.
const Trk & get_electron(const Evt &evt)
Get first electron from the event tracklist.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
double getAmbientPressure()
Get ambient pressure.
void load(const JString &file_name, JDetector &detector)
Load detector from input file.
bool from_hadron(const Hit &hit)
Test whether given hit was produced by a hadronic shower.
JDirection3D getDirection(const Vec &v)
Get direction.
double getNPE(const Hit &hit)
Get true charge of hit.
const JLimit & getLimit() const
Get limit.
bool from_electron(const Hit &hit)
Test whether given hit was produced by an electron.
int count_hadrons(const Evt &evt)
Count the number of hadrons in a given event (not including neutral pions)
#define DEBUG(A)
Message macros.
JPosition3D getPosition(const Vec &v)
Get position.
double getAngularAcceptance(const double x)
Angular acceptence of PMT.