Jpp  18.6.0-rc.1
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JHistHDF.cc
Go to the documentation of this file.
1 #include <string>
2 #include <iostream>
3 #include <fstream>
4 #include <iomanip>
5 
9 
10 #include "JAAnet/JHead.hh"
11 #include "JAAnet/JAAnetToolkit.hh"
12 
15 #include "JSupport/JSupport.hh"
16 
17 #include "JPhysics/JPDFToolkit.hh"
19 #include "JPhysics/JDispersion.hh"
20 #include "JPhysics/JConstants.hh"
21 #include "JPhysics/Antares.hh"
22 #include "JPhysics/KM3NeT.hh"
23 
24 #include "JIO/JFileStreamIO.hh"
25 
28 #include "JGeometry3D/JAxis3D.hh"
30 
31 #include "JTools/JHistogram1D_t.hh"
34 #include "JTools/JFunction1D_t.hh"
36 
37 #include "JDetector/JDetector.hh"
39 #include "JDetector/JPMTRouter.hh"
40 
42 
43 #include "Jeep/JParser.hh"
44 #include "Jeep/JMessage.hh"
45 
46 
47 namespace {
48 
49  using namespace JPP;
50 
51 
52  /**
53  * Function for classifying hits generated with KM3 or JSirene.
54  *
55  * \param track muon track
56  * \param hit hit
57  * \return 1 for muon; 2 for shower; -1 for anything else
58  */
59  inline int classify(const Trk& track, const Hit& hit)
60  {
61  if (track.id == hit.origin) {
62  switch(hit.type) {
64  return 1;
66  return 1;
68  return 2;
70  return 2;
71  default:
72  return -1;
73  }
74  }
75  return -1;
76  }
77 
78 
79  /**
80  * Function for classifying hits generated with KM3Sim.
81  *
82  * \param track muon track
83  * \param hit hit
84  * \return 1 for muon; 2 for shower
85  */
86  inline int classify_km3sim(const Trk& track, const Hit& hit)
87  {
88  if (from_muon(hit))
89  return 1;
90  else
91  return 2;
92  }
93 }
94 
95 
96 /**
97  * \file
98  *
99  * Program to histogram event-by-event data of muon light for making PDFs.
100  * \author mdejong
101  */
102 int main(int argc, char **argv)
103 {
104  using namespace std;
105  using namespace JPP;
106  using namespace KM3NET;
107 
108  JMultipleFileScanner<Evt> inputFile;
109  JLimit_t& numberOfEvents = inputFile.getLimit();
110  string outputFile;
111  string detectorFile;
112  int debug;
113 
114  try {
115 
116  JParser<> zap("Program to histogram event-by-event data of muon light for making PDFs.");
117 
118  zap['f'] = make_field(inputFile);
119  zap['o'] = make_field(outputFile);
120  zap['n'] = make_field(numberOfEvents) = JLimit::max();
121  zap['a'] = make_field(detectorFile);
122  zap['d'] = make_field(debug) = 1;
123 
124  zap(argc, argv);
125  }
126  catch(const exception &error) {
127  FATAL(error.what() << endl);
128  }
129 
130 
132 
133  try {
134  load(detectorFile, detector);
135  }
136  catch(const JException& error) {
137  FATAL(error);
138  }
139 
140  JHead head(JMultipleFileScanner<Head>(inputFile).getHeader()); // Monte Carlo header
141 
142  const bool KM3Sim = is_km3sim(head);
143  const Vec offset = getOffset(head);
144 
145  NOTICE("Apply detector offset " << offset << endl);
146 
147  detector -= getPosition(offset);
148 
149  const JCylinder3D can(detector.begin(), detector.end());
150 
151  const JPMTRouter pmtRouter(detector);
152 
153  const double P_atm = NAMESPACE::getAmbientPressure();
154  const double wmax = getMaximalWavelength();
155 
156  const JDispersion dispersion(P_atm);
157 
158  typedef JHistogram1D_t::abscissa_type abscissa_type;
159 
163  JHistogramGridMap_t>::maplist> JMultiHistogram_t;
164 
165  typedef JPDFTransformer<3, abscissa_type> JFunction3DTransformer_t;
166 
167  JMultiHistogram_t h0; // occurrence rate of PMT (used for normalisation)
168  JMultiHistogram_t h1; // light from muon
169  JMultiHistogram_t h2; // light from showers
170 
171  const double cmin = dispersion.getKmin (wmax);
172 
173  h1.transformer.reset(new JFunction3DTransformer_t(35.0, 2, cmin, 0.0, NAMESPACE::getAngularAcceptance, 0.10));
174  h2.transformer.reset(new JFunction3DTransformer_t(35.0, 2, cmin, 0.0, NAMESPACE::getAngularAcceptance, 0.10));
175 
176 
177  // configure bins
178 
179  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 };
180 
181  for (int i = 0; i != sizeof(R)/sizeof(R[0]); ++i) {
182 
183  const double R_m = R[i];
184 
185  const double grid = 10.0 + 0.0 * R_m/100.0; // [deg]
186  const double alpha = 2.0 * sqrt(1.0 - cos(grid * PI / 180.0)); // azimuth angle unit step size
187 
188  const int number_of_theta_points = max(2, (int) (180.0/(1.4 * grid)));
189  const double theta_step = PI / (number_of_theta_points + 1);
190 
191  for (double theta = -0.5*theta_step; theta < PI + theta_step; theta += theta_step) {
192 
193  const int number_of_phi_points = max(2, (int) (PI * sin(theta) / alpha));
194  const double phi_step = PI / (number_of_phi_points + 1);
195 
196  for (double phi = -0.5*phi_step; phi < PI + phi_step; phi += phi_step) {
197 
198  for (JMultiHistogram_t* p : { &h0, &h1, &h2 }) {
199  (*p)[R_m][theta][phi];
200  }
201  }
202  }
203  }
204 
205  for (JMultiHistogram_t::super_iterator
206  i1 = h1.super_begin(),
207  i2 = h2.super_begin(); i1 != h1.super_end(); ++i1, ++i2) {
208 
209  for (double x : { 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} ) {
210  i1.getValue()[x];
211  i2.getValue()[x];
212  }
213  }
214 
215 
216  while (inputFile.hasNext()) {
217 
218  NOTICE("event: " << setw(10) << inputFile.getCounter() << '\r'); STATUS(endl);
219 
220  const Evt* event = inputFile.next();
221 
222  for (vector<Trk>::const_iterator track = event->mc_trks.begin(); track != event->mc_trks.end(); ++track) {
223 
224  if (is_muon(*track)) {
225 
226  // track parameters
227 
228  const double Evis = getVisibleEnergy (*event, can);
229  const Vec EvisVector = getVisibleEnergyVector(*event, can);
230 
231  JVertex3D vertex = getVertex(*track);
232 
233  const JRotation3D R(getDirection(EvisVector));
234 
235  vertex.rotate(R);
236 
237  JRange<double> Z(0.0, gWater(track->E)); // range of muon
238 
239  for (vector<Hit>::const_iterator hit = event->mc_hits.begin(); hit != event->mc_hits.end(); ++hit) {
240 
241  try {
242 
243  JAxis3D axis = pmtRouter.getPMT(hit->pmt_id);
244 
245  axis.transform(R, vertex);
246 
247  const double t1 = vertex.getT() + (axis.getZ() + axis.getX()*getTanThetaC()) / getSpeedOfLight();
248 
249  const double R_m = axis.getX();
250  const double theta = axis.getTheta();
251  const double phi = fabs(axis.getPhi());
252  const double dt = getTime(*hit) - t1;
253  const double npe = getNPE (*hit);
254 
255  const int ID = (KM3Sim ? classify_km3sim(*track, *hit) : classify(*track, *hit));
256 
257  switch (ID) {
258 
259  case 1:
260  h1.fill(R_m, theta, phi, dt, npe); // light from muon
261  break;
262 
263  case 2:
264  h2.fill(R_m, theta, phi, dt, npe / max(1.0, Evis)); // light from shower
265  break;
266  }
267  }
268  catch(const exception& error) {
269  FATAL(error.what() << endl);
270  }
271  }
272 
273 
274  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
275 
276  JPosition3D P = module->getPosition();
277 
278  P.rotate(R);
279 
280  P.sub(vertex);
281 
282  if (sqrt(P.getX()*P.getX() + P.getY()*P.getY()) <= h0.getXmax()) {
283 
284  for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
285 
286  JAxis3D axis = *pmt;
287 
288  axis.transform(R, vertex);
289 
290  if (Z(axis.getZ() - axis.getX()/getTanThetaC())) {
291  h0.fill(axis.getX(), axis.getTheta(), fabs(axis.getPhi()), 0.0, 1.0);
292  }
293  }
294  }
295  }
296  }
297  }
298  }
299  NOTICE(endl);
300 
301 
302  // output
303 
304  JFileStreamWriter out(outputFile.c_str());
305 
306  for (const JMultiHistogram_t* p : { &h0, &h1, &h2 }) {
307  out.store(*p);
308  }
309 
310  out.close();
311 }
Router for direct addressing of PMT data in detector data structure.
Definition: JPMTRouter.hh:35
JVertex3D getVertex(const Trk &track)
Get vertex.
Utility class to parse command line options.
Definition: JParser.hh:1711
General exception.
Definition: JException.hh:24
double getKmin(const double lambda) const
Get smallest index of refraction for Bremsstrahlung light (i.e. point at which dt/dz = 0)...
JWriter & store(const JSerialisable &object)
Write object.
int main(int argc, char *argv[])
Definition: Main.cc:15
ROOT TTree parameter settings of various packages.
Auxiliary methods for PDF calculations.
Vec getOffset(const JHead &header)
Get offset.
Scattered light from muon.
#define STATUS(A)
Definition: JMessage.hh:63
Detector data structure.
Definition: JDetector.hh:89
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 function (PDF) of the time response of ...
Rotation matrix.
Definition: JRotation3D.hh:111
Properties of Antares PMT and deep-sea water.
int pmt_id
global PMT identifier as found in evt files
Definition: Hit.hh:20
Scattered light from Bremsstrahlung.
Vec getVisibleEnergyVector(const Trk &track, const JCylinder3D &can=getMaximumContainmentVolume())
Get the visible energy vector of a track.
static const JGeaneWater gWater
Function object for energy loss of muon in sea water.
Definition: JGeane.hh:381
Direct light from Bremsstrahlung.
double getTime(const Hit &hit)
Get true time of hit.
double getPhi() const
Get phi angle.
Definition: JVersor3D.hh:144
string outputFile
double E
Energy [GeV] (either MC truth or reconstructed)
Definition: Trk.hh:20
#define KM3NET
int origin
track id of the track that created this hit (mc only)
Definition: Hit.hh:29
Data structure for detector geometry and calibration.
void transform(const JAxis3D &axis)
Transform axis to reference frame of given axis.
Definition: JAxis3D.hh:359
Various implementations of functional maps.
Auxiliary methods for evaluating visible energies.
Head getHeader(const JMultipleFileScanner_t &file_list)
Get Monte Carlo header.
Axis object.
Definition: JAxis3D.hh:38
then usage set_variable ACOUSTICS_WORKDIR $WORKDIR set_variable FORMULA sin([0]+2 *$PI *([1]+[2]*x)*x)" set_variable DY 1.0e-8 mkdir $WORKDIR for DETECTOR in $DETECTORS[*]
Auxiliary class for defining the range of iterations of objects.
Definition: JLimit.hh:41
JVector3D & sub(const JVector3D &vector)
Subtract vector.
Definition: JVector3D.hh:158
Properties of KM3NeT PMT and deep-sea water.
The Vec class is a straightforward 3-d vector, which also works in pyroot.
Definition: Vec.hh:12
double getMaximalWavelength()
Get maximal wavelength for PDF evaluations.
Definition: JPDFToolkit.hh:46
Cylinder object.
Definition: JCylinder3D.hh:39
Detector file.
Definition: JHead.hh:226
JDirection3D getDirection(const Vec &dir)
Get direction.
Monte Carlo run header.
Definition: JHead.hh:1234
Transformable multidimensional histogram.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:2158
double getTheta() const
Get theta angle.
Definition: JVersor3D.hh:128
Type definition of a JHistogramMap based on JGridMap implementation.
JPosition3D getPosition(const Vec &pos)
Get position.
double getAmbientPressure()
Get ambient pressure.
Definition: Antares.hh:40
Definition of hit and track types and auxiliary methods for handling Monte Carlo data.
#define NOTICE(A)
Definition: JMessage.hh:64
double getVisibleEnergy(const Trk &track, const JCylinder3D &can=getMaximumContainmentVolume())
Get the visible energy of a track.
Physics constants.
static const double PI
Mathematical constants.
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:104
Direct access to PMT in detector data structure.
Auxiliary class for recursive map list generation.
Definition: JMapList.hh:108
Type definition of a JHistogramMap based on JMap implementation.
General purpose messaging.
Template specialisation of JMultipleFileScanner for Monte Carlo header.
#define FATAL(A)
Definition: JMessage.hh:67
Scanning of objects from multiple files according a format that follows from the extension of each fi...
int id
track identifier
Definition: Trk.hh:16
then JCookie sh JDataQuality D $DETECTOR_ID R
Definition: JDataQuality.sh:41
Definition: Hit.hh:8
const double getSpeedOfLight()
Get speed of light.
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
Direct light from muon.
General purpose class for object reading from a list of file names.
Utility class to parse command line options.
double getX() const
Get x position.
Definition: JVector3D.hh:94
collection_type::abscissa_type abscissa_type
double getTanThetaC()
Get average tangent of Cherenkov angle of water corresponding to group velocity.
bool is_km3sim(const JHead &header)
Check for detector simulation.
double getNPE(const Hit &hit)
Get true charge of hit.
do set_variable MODULE getModule a $WORKDIR detector_a datx L $STRING JEditDetector a $WORKDIR detector_a datx M $MODULE setz o $WORKDIR detector_a datx JEditDetector a $WORKDIR detector_b datx M $MODULE setz o $WORKDIR detector_b datx done echo Output stored at $WORKDIR detector_a datx and $WORKDIR tripod_a txt JDrawDetector2D a $WORKDIR detector_a datx a $WORKDIR detector_b datx L BL o detector $FORMAT $BATCH JDrawDetector2D T $WORKDIR tripod_a txt T $WORKDIR tripod_b txt L BL o tripod $FORMAT $BATCH JCompareDetector a $WORKDIR detector_a datx b $WORKDIR detector_b datx o $WORKDIR abc root &dev null for KEY in X Y Z
Data structure for position in three dimensions.
Definition: JPosition3D.hh:36
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:84
do set_variable DETECTOR_TXT $WORKDIR detector
double getT(const JVector3D &pos) const
Get arrival time of Cherenkov light at given position.
Definition: JVertex3D.hh:147
The Trk class represents a Monte Carlo (MC) particle as well as a reconstructed track/shower.
Definition: Trk.hh:14
JPosition3D & rotate(const JRotation3D &R)
Rotate.
Definition: JPosition3D.hh:186
Binary buffered file output.
int type
particle type or parametrisation used for hit (mc only)
Definition: Hit.hh:28
double getZ() const
Get z position.
Definition: JVector3D.hh:115
int debug
debug level
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
Definition: Evt.hh:20
then $DIR JPlotNPE PDG P
Definition: JPlotNPE-PDG.sh:62
The cylinder used for photon tracking.
Definition: JHead.hh:575
double getAngularAcceptance(const double x)
Angular acceptence of PMT.
Definition: JDrawLED.cc:68
bool from_muon(const Hit &hit)
Test whether given hit was produced by a muon.