Jpp  19.1.0-rc.1
the software that should make you happy
JHistHDE.cc
Go to the documentation of this file.
1 #include <string>
2 #include <iostream>
3 #include <fstream>
4 #include <iomanip>
5 #include <stdexcept>
6 
7 #include "km3net-dataformat/offline/Head.hh"
8 #include "km3net-dataformat/offline/Evt.hh"
9 #include "km3net-dataformat/offline/Hit.hh"
10 #include "km3net-dataformat/offline/Trk.hh"
11 
12 #include "JAAnet/JHead.hh"
13 #include "JAAnet/JAAnetToolkit.hh"
14 
15 #include "JSupport/JMultipleFileScanner.hh"
16 #include "JSupport/JMonteCarloFileSupportkit.hh"
17 #include "JSupport/JSupport.hh"
18 
19 #include "JPhysics/JPDFToolkit.hh"
20 #include "JPhysics/JDispersion.hh"
21 #include "JPhysics/Antares.hh"
22 #include "JPhysics/KM3NeT.hh"
23 #include "JIO/JFileStreamIO.hh"
24 
25 #include "JGeometry3D/JGeometry3DToolkit.hh"
26 #include "JGeometry3D/JPosition3D.hh"
27 #include "JGeometry3D/JDirection3D.hh"
28 #include "JGeometry3D/JAxis3D.hh"
29 
30 #include "JPhysics/JConstants.hh"
31 #include "JTools/JHistogram1D_t.hh"
32 #include "JTools/JHistogramMap_t.hh"
33 #include "JTools/JTransformableMultiHistogram.hh"
34 #include "JTools/JFunction1D_t.hh"
35 #include "JTools/JFunctionalMap_t.hh"
36 
37 #include "JDetector/JDetector.hh"
38 #include "JDetector/JDetectorToolkit.hh"
39 #include "JDetector/JPMTRouter.hh"
40 
41 #include "Jeep/JParser.hh"
42 #include "Jeep/JMessage.hh"
43 
44 #include "JPhysics/JPDFTransformer.hh"
45 #include "JPhysics/JPDF.hh"
46 
47 #include "JAAnet/JParticleTypes.hh"
48 
49 /**
50  * \file
51  *
52  * Program to histogram event-by-event and track-by-track for making PDFs of lights of individual particles.
53  * \author jseneca
54  */
55 inline std::vector<Int_t> getHitRemapping(const std::vector<Trk>* tracklist)
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 }
75 
76 int main(int argc, char **argv)
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 }
Antares.hh
Properties of Antares PMT and deep-sea water.
JAAnetToolkit.hh
Definition of hit and track types and auxiliary methods for handling Monte Carlo data.
outputFile
string outputFile
Definition: JDAQTimesliceSelector.cc:37
JDetector.hh
Data structure for detector geometry and calibration.
getAngularAcceptance
double getAngularAcceptance(const double x)
Angular acceptence of PMT.
Definition: JDrawLED.cc:68
JFunction1D_t.hh
JFunctionalMap_t.hh
Various implementations of functional maps.
getHitRemapping
std::vector< Int_t > getHitRemapping(const std::vector< Trk > *tracklist)
Definition: JHistHDE.cc:55
main
int main(int argc, char **argv)
Definition: JHistHDE.cc:76
JHistogram1D_t.hh
JHistogramMap_t.hh
JMessage.hh
General purpose messaging.
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
JMultipleFileScanner.hh
Scanning of objects from multiple files according a format that follows from the extension of each fi...
JPDFToolkit.hh
Auxiliary methods for PDF calculations.
JPMTRouter.hh
Direct access to PMT in detector data structure.
JParser.hh
Utility class to parse command line options.
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:2158
JParticleTypes.hh
Definition of particle types.
JConstants.hh
Physics constants.
JSupport.hh
ROOT TTree parameter settings of various packages.
KM3NeT.hh
Properties of KM3NeT PMT and deep-sea water.
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
JDETECTOR::JPMTRouter::getPMT
const JPMT & getPMT(const JPMTAddress &address) const
Get PMT.
Definition: JPMTRouter.hh:92
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
JIO::JFileStreamWriter::close
void close()
Close file.
Definition: JFileStreamIO.hh:142
JIO::JWriter::store
JWriter & store(const JSerialisable &object)
Write object.
Definition: JSerialisable.hh:173
JLANG::JException
General exception.
Definition: JException.hh:24
JPARSER::JParser
Utility class to parse command line options.
Definition: JParser.hh:1714
JPARSER::JParser::read
int read(const int argc, const char *const argv[])
Parse the program's command line options.
Definition: JParser.hh:2008
JPHYSICS::JDispersionInterface::getIndexOfRefractionGroup
virtual double getIndexOfRefractionGroup(const double lambda) const
Index of refraction for group velocity.
Definition: JDispersionInterface.hh:52
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
JAANET::TRACK_TYPE_NUMU
@ TRACK_TYPE_NUMU
Definition: JParticleTypes.hh:99
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
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition: JAAnetToolkit.hh:43
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
std
Definition: JSTDTypes.hh:14
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