Jpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Functions
JHistHDF.cc File Reference

Program to histogram event-by-event data of muon light for making PDFs. More...

#include <string>
#include <iostream>
#include <fstream>
#include <iomanip>
#include "evt/Head.hh"
#include "evt/Evt.hh"
#include "evt/Hit.hh"
#include "JAAnet/JHead.hh"
#include "JAAnet/JHeadToolkit.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/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/JConstants.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"

Go to the source code of this file.

Functions

int main (int argc, char **argv)
 

Detailed Description

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

Author
mdejong

Definition in file JHistHDF.cc.

Function Documentation

int main ( int  argc,
char **  argv 
)

Definition at line 164 of file JHistHDF.cc.

165 {
166  using namespace std;
167  using namespace JPP;
168 
169  JMultipleFileScanner<Evt> inputFile;
170  JLimit_t& numberOfEvents = inputFile.getLimit();
171  string outputFile;
172  string detectorFile;
173  int debug;
174 
175  try {
176 
177  JParser<> zap("Program to histogram event-by-event data of muon light for making PDFs.");
178 
179  zap['f'] = make_field(inputFile);
180  zap['o'] = make_field(outputFile);
181  zap['n'] = make_field(numberOfEvents) = JLimit::max();
182  zap['a'] = make_field(detectorFile);
183  zap['d'] = make_field(debug) = 1;
184 
185  zap(argc, argv);
186  }
187  catch(const exception &error) {
188  FATAL(error.what() << endl);
189  }
190 
191 
192  using namespace KM3NET;
193 
194 
196 
197  try {
198  load(detectorFile, detector);
199  }
200  catch(const JException& error) {
201  FATAL(error);
202  }
203 
204  JHead head(JMultipleFileScanner<Head>(inputFile).getHeader()); // Monte Carlo header
205 
206  const JPMTRouter pmtRouter(detector);
207 
208  const double P_atm = NAMESPACE::getAmbientPressure();
209  const double wmax = getMaximalWavelength();
210 
211  const JDispersion dispersion(P_atm);
212 
213  typedef JHistogram1D_t::abscissa_type abscissa_type;
214 
218  JHistogramGridMap_t>::maplist> JMultiHistogram_t;
219 
220  typedef JPDFTransformer<3, abscissa_type> JFunction3DTransformer_t;
221 
222  JMultiHistogram_t h0; // occurrence rate of PMT (used for normalisation)
223  JMultiHistogram_t h1; // light from muon
224  JMultiHistogram_t h2; // light from showers
225 
226  const JHitClassifier<JMultiHistogram_t> classifier(head, detector);
227 
228  const double cmin = dispersion.getKmin (wmax);
229 
230  h1.transformer.reset(new JFunction3DTransformer_t(35.0, 2, cmin, 0.0, NAMESPACE::getAngularAcceptance, 0.10));
231  h2.transformer.reset(new JFunction3DTransformer_t(35.0, 2, cmin, 0.0, NAMESPACE::getAngularAcceptance, 0.10));
232 
233 
234  // configure bins
235 
236  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 
238  for (int i = 0; i != sizeof(R)/sizeof(R[0]); ++i) {
239 
240  const double R_m = R[i];
241 
242  const double grid = 10.0 + 0.0 * R_m/100.0; // [deg]
243  const double alpha = 2.0 * sqrt(1.0 - cos(grid * PI / 180.0)); // azimuth angle unit step size
244 
245  const int number_of_theta_points = max(2, (int) (180.0/(1.4 * grid)));
246  const double theta_step = PI / (number_of_theta_points + 1);
247 
248  for (double theta = -0.5*theta_step; theta < PI + theta_step; theta += theta_step) {
249 
250  const int number_of_phi_points = max(2, (int) (PI * sin(theta) / alpha));
251  const double phi_step = PI / (number_of_phi_points + 1);
252 
253  for (double phi = -0.5*phi_step; phi < PI + phi_step; phi += phi_step) {
254 
255  for (JMultiHistogram_t* buffer[] = { &h0, &h1, &h2, NULL }, **histogram = buffer; *histogram != NULL; ++histogram) {
256  (**histogram)[R_m][theta][phi];
257  }
258  }
259  }
260  }
261 
262  for (JMultiHistogram_t::super_iterator
263  i1 = h1.super_begin(),
264  i2 = h2.super_begin(); i1 != h1.super_end(); ++i1, ++i2) {
265 
266  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) {
267  i1.getValue()[*x];
268  i2.getValue()[*x];
269  }
270  }
271 
272 
273  while (inputFile.hasNext()) {
274 
275  NOTICE("event: " << setw(10) << inputFile.getCounter() << '\r'); STATUS(endl);
276 
277  const Evt* event = inputFile.next();
278 
279  for (vector<Trk>::const_iterator track = event->mc_trks.begin(); track != event->mc_trks.end(); ++track) {
280 
281  if (is_muon(*track)) {
282 
283  // track parameters
284 
285  double E = track->E;
286  JPosition3D pos = getPosition (*track);
287  JDirection3D dir = getDirection(*track);
288  double t0 = track->t;
289 
290  const JRotation3D R(dir);
291 
292  pos.rotate(R);
293 
294  JRange<double> Z(0.0, gWater(E)); // range of muon
295 
296  for (vector<Hit>::const_iterator hit = event->mc_hits.begin(); hit != event->mc_hits.end(); ++hit) {
297 
298  try {
299 
300  JAxis3D axis = pmtRouter.getPMT(hit->pmt_id);
301 
302  axis.transform(R, pos);
303 
304  const double E1 = classifier.getEnergy(*event, *track, axis);
305  const double t1 = t0 + (axis.getZ() + axis.getX()*getTanThetaC()) / getSpeedOfLight();
306 
307  const double R_m = axis.getX();
308  const double theta = axis.getTheta();
309  const double phi = fabs(axis.getPhi());
310  const double dt = getTime(*hit) - t1;
311  const double npe = getNPE(*hit);
312 
313  switch (classifier(*track, *hit)) {
314 
315  case 1:
316  h1.fill(R_m, theta, phi, dt, npe); // light from muon
317  break;
318 
319  case 2:
320  h2.fill(R_m, theta, phi, dt, npe/max(1.0, E1)); // light from shower
321  break;
322  }
323  }
324  catch(const exception& error) {
325  FATAL(error.what() << endl);
326  }
327  }
328 
329 
330  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
331 
332  JPosition3D P = module->getPosition();
333 
334  P.rotate(R);
335 
336  P.sub(pos);
337 
338  if (sqrt(P.getX()*P.getX() + P.getY()*P.getY()) <= h0.getXmax()) {
339 
340  for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
341 
342  JAxis3D axis = *pmt;
343 
344  axis.transform(R, pos);
345 
346  if (Z(axis.getZ() - axis.getX()/getTanThetaC())) {
347  h0.fill(axis.getX(), axis.getTheta(), fabs(axis.getPhi()), 0.0, 1.0);
348  }
349  }
350  }
351  }
352  }
353  }
354  }
355  NOTICE(endl);
356 
357 
358  // output
359 
360  JFileStreamWriter out(outputFile.c_str());
361 
362  for (const JMultiHistogram_t* buffer[] = { &h0, &h1, &h2, NULL }, **i = buffer; *i != NULL; ++i) {
363  out.store(**i);
364  }
365 
366  out.close();
367 }
Router for direct addressing of PMT data in detector data structure.
Definition: JPMTRouter.hh:33
Utility class to parse command line options.
Definition: JParser.hh:1410
General exception.
Definition: JException.hh:40
#define KM3NET
Definition: JDAQ.hh:20
JWriter & store(const JSerialisable &object)
Write object.
static const JGeane gWater(2.67e-1 *JTOOLS::DENSITY_SEA_WATER, 3.4e-4 *JTOOLS::DENSITY_SEA_WATER)
Function object for Energy loss of muon in sea water.
Data structure for direction in three dimensions.
Definition: JDirection3D.hh:32
const double getSpeedOfLight()
Number of bytes in a gigabyte.
Definition: JConstants.hh:89
#define STATUS(A)
Definition: JMessage.hh:61
Detector data structure.
Definition: JDetector.hh:77
Implementation of dispersion for water in deep sea.
Definition: JDispersion.hh:26
bool is_muon(const Trk &track)
Test whether given track is a (anti-)muon.
Template definition of transformer of the Probability Density Functions of the time response of a PMT...
Rotation matrix.
Definition: JRotation3D.hh:108
static const double PI
Constants.
Definition: JConstants.hh:20
Structure to store the ToT mean and standard deviation of the hits produced by a nanobeacon in a sour...
double getTime(const Hit &hit)
Get true time of hit.
double getPhi() const
Get phi angle.
Definition: JVersor3D.hh:139
string outputFile
void transform(const JAxis3D &axis)
Transform axis to reference frame of given axis.
Definition: JAxis3D.hh:354
Head getHeader(const JMultipleFileScanner_t &file_list)
Get Monte Carlo header.
Axis object.
Definition: JAxis3D.hh:37
Auxiliary class for defining the range of iterations of objects.
Definition: JLimit.hh:41
JVector3D & sub(const JVector3D &vector)
Subtract vector.
Definition: JVector3D.hh:156
double getMaximalWavelength()
Get maximal wavelength for PDF evaluations.
Definition: JPDFToolkit.hh:49
Detector file.
Definition: JHead.hh:126
Transformable multidimensional histogram.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1836
double getTheta() const
Get theta angle.
Definition: JVersor3D.hh:123
Type definition of a JHistogramMap based on JGridMap implementation.
double getAmbientPressure()
Get ambient pressure.
Definition: Antares.hh:30
#define NOTICE(A)
Definition: JMessage.hh:62
JHistogram1D< JElement2D< double, double >, JCollection > JHistogram1D_t
Type definition of a 1 dimensional histogram based on a JCollection.
double getY() const
Get y position.
Definition: JVector3D.hh:102
void load(const JString &file_name, JDetector &detector)
Load detector from input file.
int debug
debug level
Definition: JSirene.cc:59
const JPosition3D & getPosition() const
Get position.
Definition: JPosition3D.hh:129
Auxiliary class for recursive map list generation.
Definition: JMapList.hh:108
Type definition of a JHistogramMap based on JMap implementation.
Monte Carlo run header.
Definition: JHead.hh:727
Template specialisation of JMultipleFileScanner for Monte Carlo header.
#define FATAL(A)
Definition: JMessage.hh:65
General purpose class for object reading from a list of file names.
double getX() const
Get x position.
Definition: JVector3D.hh:92
collection_type::abscissa_type abscissa_type
JDirection3D getDirection(const Vec &v)
Get direction.
double getNPE(const Hit &hit)
Get true charge of hit.
Data structure for position in three dimensions.
Definition: JPosition3D.hh:35
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:72
double getTanThetaC()
Get average tangent of Cherenkov angle of water.
Definition: JConstants.hh:122
JPosition3D & rotate(const JRotation3D &R)
Rotate.
Definition: JPosition3D.hh:185
Binary buffered file output.
double getZ() const
Get z position.
Definition: JVector3D.hh:113
JPosition3D getPosition(const Vec &v)
Get position.
double getAngularAcceptance(const double x)
Angular acceptence of PMT.
Definition: JDrawLED.cc:84