Jpp master_rocky-44-g75b7c4f75
the software that should make you happy
Loading...
Searching...
No Matches
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 */
55inline 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
76int 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.
JFunction1D_t.hh
JFunctionalMap_t.hh
Various implementations of functional maps.
main
int main(int argc, char **argv)
Definition JHistHDE.cc:76
getHitRemapping
std::vector< Int_t > getHitRemapping(const std::vector< Trk > *tracklist)
Definition JHistHDE.cc:55
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:2142
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::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:1698
JPARSER::JParser::read
int read(const int argc, const char *const argv[])
Parse the program's command line options.
Definition JParser.hh:1992
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
General purpose class for object reading from a list of file names.
Definition JMultipleFileScanner.hh:306
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::JCollection< JElement2D_t >::const_iterator
container_type::const_iterator const_iterator
Definition JCollection.hh:91
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
JAANET::getDirection
JDirection3D getDirection(const Vec &dir)
Get direction.
Definition JAAnetToolkit.hh:214
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:505
JLANG::to_string
std::string to_string(const T &value)
Convert value to string.
Definition JLangToolkit.hh:207
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
JSYSTEM::getTime
const char * getTime()
Get current local time conform ISO-8601 standard.
Definition JDateAndTime.hh:55
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.12.0