Jpp  debug
the software that should make you happy
Functions
JHistHDE.cc File Reference

Program to histogram event-by-event and track-by-track for making PDFs of lights of individual particles. More...

#include <string>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <stdexcept>
#include "km3net-dataformat/offline/Head.hh"
#include "km3net-dataformat/offline/Evt.hh"
#include "km3net-dataformat/offline/Hit.hh"
#include "km3net-dataformat/offline/Trk.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/Antares.hh"
#include "JPhysics/KM3NeT.hh"
#include "JIO/JFileStreamIO.hh"
#include "JGeometry3D/JGeometry3DToolkit.hh"
#include "JGeometry3D/JPosition3D.hh"
#include "JGeometry3D/JDirection3D.hh"
#include "JGeometry3D/JAxis3D.hh"
#include "JPhysics/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"
#include "JPhysics/JPDFTransformer.hh"
#include "JPhysics/JPDF.hh"
#include "JAAnet/JParticleTypes.hh"

Go to the source code of this file.

Functions

std::vector< Int_t > getHitRemapping (const std::vector< Trk > *tracklist)
 
int main (int argc, char **argv)
 

Detailed Description

Program to histogram event-by-event and track-by-track for making PDFs of lights of individual particles.

Author
jseneca

Definition in file JHistHDE.cc.

Function Documentation

◆ getHitRemapping()

std::vector<Int_t> getHitRemapping ( const std::vector< Trk > *  tracklist)
inline

Definition at line 55 of file JHistHDE.cc.

56 {
57  using namespace std;
58  using namespace JPP;
59 
60  vector<Int_t> hit_remapping;
61 
62  for (vector<Trk>::const_iterator i = tracklist->begin(); i != tracklist->end(); ++i) {
63 
64  //Skip pathological neutrino that is not the primary incoming neutrino
65 
66  if(i->type == TRACK_TYPE_NUMU && i != tracklist->begin()) {
67  continue;
68  }
69 
70  hit_remapping.push_back(i->id);
71  }
72 
73  return hit_remapping;
74 }
JAANET::TRACK_TYPE_NUMU
@ TRACK_TYPE_NUMU
Definition: JParticleTypes.hh:99
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition: JAAnetToolkit.hh:43
std
Definition: JSTDTypes.hh:14

◆ main()

int main ( int  argc,
char **  argv 
)

Definition at line 76 of file JHistHDE.cc.

77 {
78  using namespace std;
79  using namespace JPP;
80 
81  JMultipleFileScanner<Evt> inputFile;
82  //JLimit_t& numberOfEvents = inputFile.getLimit();
83  string outputFile = "J%s.HDE"; //'s' for single (particle), or simulation
84  string detectorFile;
85  vector<int> particles;
86  uint fix_bug;
87  int debug;
88 
89  try {
90 
91  JParser<> zap("Program to histogram event-by-event data of shower light for making PDFs.");
92 
93  zap['f'] = make_field(inputFile);
94  zap['o'] = make_field(outputFile);
95  //zap['n'] = make_field(numberOfEvents) = JLimit::max();
96  zap['a'] = make_field(detectorFile);
97  zap['p'] = make_field(particles);
98  zap['b'] = make_field(fix_bug);
99  zap['d'] = make_field(debug) = 1;
100 
101  if (zap.read(argc, argv) != 0)
102  return 1;
103  }
104  catch(const exception &error) {
105  FATAL(error.what() << endl);
106  }
107 
108  JDetector detector;
109 
110  if(outputFile.find('%') != string::npos){
111  sort(particles.begin(), particles.end());
112  string prefix_t = "";
113  for(vector<int>::iterator ptype = particles.begin(); ptype != particles.end(); ptype++){
114  prefix_t += to_string((long long int)*ptype) + "_";
115  }
116  prefix_t.erase(prefix_t.size() - 1);
117  string::size_type pos_1 = outputFile.find('%');
118  outputFile.replace(pos_1, 1, prefix_t);
119  }
120 
121  try {
122  load(detectorFile, detector);
123  }
124  catch(const JException& error) {
125  FATAL(error);
126  }
127 
128  JHead head(JMultipleFileScanner<Head>(inputFile).getHeader()); // Monte Carlo header
129 
130  const Vec offset = getOffset(head);
131 
132  NOTICE("Apply detector offset " << offset << endl);
133 
134  detector -= getPosition(offset);
135 
136  const JPMTRouter pmtRouter(detector);
137 
138  const double P_atm = NAMESPACE::getAmbientPressure();
139  const double wmin = getMinimalWavelength();
140  const double wmax = getMaximalWavelength();
141 
142  const JDispersion dispersion(P_atm);
143 
144  const double ng[] = { dispersion.getIndexOfRefractionGroup(wmax),
145  dispersion.getIndexOfRefractionGroup(wmin) };
146 
147  typedef JHistogram1D_t::abscissa_type abscissa_type;
148 
149  typedef JTransformableMultiHistogram<JHistogram1D_t,
150  JMAPLIST<JHistogramMap_t,
151  JHistogramMap_t,
152  JHistogramMap_t,
153  JHistogramGridMap_t,
154  JHistogramGridMap_t>::maplist> JMultiHistogram_t;
155 
156  typedef JPDFTransformer<5, abscissa_type> JFunction5DTransformer_t;
157 
158  JMultiHistogram_t h0; // occurrence rate of PMT (used for normalisation)
159  JMultiHistogram_t h1; // light from track cascade
160 
161  h1.transformer.reset(new JFunction5DTransformer_t(21.5, 2, ng[0], 0.0, JGeant(JGeanx(1.00, -2.2)), 1e-2, NAMESPACE::getAngularAcceptance, 0.05));
162 
163  set<double> C; // cosine emission angle
164 
165  JQuadrature qeant(-1.0, 1.0, 30, JGeanx(1.00, -2.2));
166 
167  for (JQuadrature::const_iterator i = qeant.begin(); i != qeant.end(); ++i) {
168  C.insert(i->getX());
169  }
170 
171  C.insert(-1.01);
172  C.insert(-1.00);
173  C.insert( 1.00);
174  C.insert( 1.01);
175 
176  const double E_b[] = {0.01, 0.2, 0.4, 0.6, 0.8, 1.0, 1.5, 2.0, 3.0, 4.0, 5.0, 6.0, 8.0, 10.0, 12.0, 14.0, 16.0, 18.0, 20.0, 25.0, 30.0, 40.0, 55.0, 70.0, 85.0, 100.0};
177 
178  const double R[] = {0.0, 1.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, 65.0, 80.0};
179  const double Dmax_m = R[sizeof(R)/sizeof(R[0]) - 1]; //Last element of R
180 
181  for (int iE = 0; iE != sizeof(E_b)/sizeof(E_b[0]); ++iE) {
182 
183  const double E_m = E_b[iE];
184 
185  for (int i = 0; i != sizeof(R)/sizeof(R[0]); ++i) {
186 
187  const double R_m = R[i];
188 
189  for (set<double>::const_iterator c = C.begin(); c != C.end(); ++c) {
190 
191  const double cd = *c;
192 
193  const double grid = 10.0 + 0.0 * R_m/100.0; // [deg]
194  const double alpha = 2.0 * sqrt(1.0 - cos(grid * PI / 180.0)); // azimuth angle unit step size
195 
196  const int number_of_theta_points = (max(2, (int) (180.0/(1.4 * grid))));
197  const double theta_step = PI / (number_of_theta_points + 1);
198 
199  for (double theta = -0.5*theta_step; theta < PI + theta_step; theta += theta_step) {
200  const int number_of_phi_points = (max(2, (int) (PI * sin(theta) / alpha)));
201  const double phi_step = PI / (number_of_phi_points + 1);
202 
203  for (double phi = -0.5*phi_step; phi < PI + phi_step; phi += phi_step) {
204 
205  for (JMultiHistogram_t* p : { &h0, &h1 }) {
206  (*p)[E_m][R_m][cd][theta][phi];
207  }
208  }
209  }
210  }
211  }
212  }
213 
214  double buffer[] = {-15.0,-10.0,-7.5,-5,-4,-3,-2,-1.5,-1.0,-0.5,-0.1, 0.0,0.1, 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 };
215  for (JMultiHistogram_t::super_iterator
216  i1 = h1.super_begin(); i1 != h1.super_end(); ++i1) {
217  for (int j = 0; j != sizeof(buffer)/sizeof(buffer[0]); ++j) {
218  i1.getValue()[buffer[j]];
219  }
220  }//Buffer defines the binning of the inner histograms (dt)
221 
222 
223  long long int h0_fillcount = 0;
224  long long int h1_fillcount = 0;
225 
226  while (inputFile.hasNext()) {
227 
228  //NOTICE("event: " << setw(10) << inputFile.getCounter() << '\r'); STATUS(endl);
229  //I suspect the printing above slows down the program considerably
230 
231  const Evt* event = inputFile.next();
232  const vector<Trk>* mc_tracks = &(event->mc_trks);
233 
234  vector<int> hit_remap;
235 
236  //*************************************************
237  //Fixes 2016 atmospheric neutrino simulation hit-making neutrino bug
238  //Remove after fix.
239  if(fix_bug == 0){
240  for (vector<Trk>::const_iterator i = mc_tracks->begin(); i != mc_tracks->end(); ++i) {
241  hit_remap.push_back(i->id);
242  }
243  } else if(fix_bug == 1){
244  hit_remap = getHitRemapping(mc_tracks);
245  } else if(fix_bug == 2){
246  bool skipevent = false;
247  for (vector<Trk>::const_iterator i = mc_tracks->begin(); i != mc_tracks->end(); ++i) {
248  if(i->type == 14){
249  skipevent = true;
250  break;
251  }
252  hit_remap.push_back(i->id);
253  }
254  if(skipevent) continue;
255  }
256  //*************************************************
257 
258  double t0 = (*mc_tracks)[1].t;
259 
260  for (vector<Hit>::const_iterator hit = event->mc_hits.begin(); hit != event->mc_hits.end(); ++hit) {
261 
262  try {
263 
264  if(hit->origin >= (int)(*mc_tracks).size()){
265  std::out_of_range err("Hit origin not in tracklist. Avoided segfault");
266  throw err;
267  }
268 
269  Int_t corr_origin = hit_remap[hit->origin];
270  Trk curr_track = (*mc_tracks)[corr_origin]; //Get track from corrected origin
271  JDirection3D dir = getDirection(curr_track); const JRotation3D R(dir);
272  JPosition3D pos = getPosition(curr_track); pos.rotate(R);
273  JAxis3D axis = pmtRouter.getPMT(hit->pmt_id);
274  axis.transform(R, pos);
275 
276  const double E = curr_track.E;
277  const double t1 = t0 + axis.getLength() * getInverseSpeedOfLight() * getIndexOfRefraction();
278  const double D_m = axis.getLength();
279  const double cd = axis.getZ()/D_m;
280  const double theta = axis.getTheta();
281  const double phi = fabs(axis.getPhi());
282  const double dt = getTime(*hit) - t1;
283  const double npe = getNPE (*hit);
284 
285  if(D_m < Dmax_m){
286  h1.fill(E, D_m, cd, theta, phi, dt, npe);
287  h1_fillcount += 1;
288  }
289  }
290  catch(const exception& error) {
291  std::cout << "Fatal error for event: " << inputFile.getCounter() << std::endl;
292  FATAL(error.what() << endl);
293  }
294  }
295 
296  // ignore primary neutrino -> ^
297  for (vector<Trk>::const_iterator tr = event->mc_trks.begin() + 1; tr != event->mc_trks.end(); ++tr) {
298 
299  if(find(particles.begin(), particles.end(), tr->type) == particles.end()){
300  continue;
301  }
302 
303  //cout << "Particle accepted: " << tr->type << endl;
304 
305 
306  JDirection3D dir = getDirection(*tr); const JRotation3D R(dir);
307  JPosition3D pos = getPosition(*tr); pos.rotate(R);
308 
309  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
310 
311  JPosition3D P = module->getPosition();
312 
313  P.rotate(R);
314 
315  P.sub(pos);
316 
317  if (P.getLength() < Dmax_m) {
318 
319  for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
320 
321  JAxis3D axis = *pmt;
322  axis.transform(R, pos);
323 
324  h0.fill(tr->E, axis.getLength(), axis.getZ()/axis.getLength(), axis.getTheta(), fabs(axis.getPhi()), 0.0, 1.0);
325  h0_fillcount += 1;
326  }
327  }
328  }
329  }
330 
331  if(h1_fillcount > h0_fillcount){
332  std::cout << "WARNING: recorded hits in normalization histogram that were not recorded in normalization histogram. This should not happen." << std::endl;
333  std::cout << "h1_fillcount: " << h1_fillcount << ", h0_fillcount: " << h0_fillcount << std::endl;
334  }
335  };
336 
337  double integral = 0;
338  for (JMultiHistogram_t::super_iterator i = h0.super_begin(); i != h0.super_end(); ++i) {
339  integral+=i.getValue().getIntegral();
340  }
341  DEBUG("Integral:\t" << integral << endl);
342 
343  // output
344 
345  JFileStreamWriter out(outputFile.c_str());
346 
347  NOTICE("Storing " << outputFile << ", " << flush);
348 
349  for (const JMultiHistogram_t* p : { &h0, &h1 }) {
350  out.store(*p);
351  }
352 
353  out.close();
354  NOTICE("JHistHDE done. " << endl);
355  return 1;
356 }
outputFile
string outputFile
Definition: JDAQTimesliceSelector.cc:37
getAngularAcceptance
double getAngularAcceptance(const double x)
Angular acceptence of PMT.
Definition: JDrawLED.cc:68
getHitRemapping
std::vector< Int_t > getHitRemapping(const std::vector< Trk > *tracklist)
Definition: JHistHDE.cc:55
DEBUG
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
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
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:2158
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
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::JDirection3D
Data structure for direction in three dimensions.
Definition: JDirection3D.hh:35
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
JIO::JFileStreamWriter
Binary buffered file output.
Definition: JFileStreamIO.hh:96
JLANG::JException
General exception.
Definition: JException.hh:24
JPARSER::JParser
Utility class to parse command line options.
Definition: JParser.hh:1714
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::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
ANTARES::getAmbientPressure
double getAmbientPressure()
Get ambient pressure.
Definition: Antares.hh:40
JAANET::getDirection
JDirection3D getDirection(const Vec &dir)
Get direction.
Definition: JAAnetToolkit.hh:214
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::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
JDETECTOR::load
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
Definition: JDetectorToolkit.hh:484
JLANG::to_string
std::string to_string(const T &value)
Convert value to string.
Definition: JLangToolkit.hh:207
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
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
Evt
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
Definition: Evt.hh:21
JAANET::detector
Detector file.
Definition: JHead.hh:227
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
Trk::E
double E
Energy [GeV] (either MC truth or reconstructed)
Definition: Trk.hh:20
Vec
The Vec class is a straightforward 3-d vector, which also works in pyroot.
Definition: Vec.hh:13
Generated by doxygen 1.9.1