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

Program to histogram event-by-event data of shower 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/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 "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"
#include "JPhysics/JPDFTransformer.hh"
#include "JPhysics/JPDF.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 shower light for making PDFs.

Author
lquinn

Definition in file JHistHDG.cc.

Function Documentation

int main ( int  argc,
char **  argv 
)

Definition at line 52 of file JHistHDG.cc.

53 {
54  using namespace std;
55  using namespace JPP;
56 
57  JMultipleFileScanner<Evt> inputFile;
58  JLimit_t& numberOfEvents = inputFile.getLimit();
59  string outputFile;
60  string detectorFile;
61  bool hadronicMode;
62  int debug;
63 
64  try {
65 
66  JParser<> zap("Program to histogram event-by-event data of shower light for making PDFs.");
67 
68  zap['f'] = make_field(inputFile);
69  zap['o'] = make_field(outputFile);
70  zap['n'] = make_field(numberOfEvents) = JLimit::max();
71  zap['a'] = make_field(detectorFile);
72  zap['H'] = make_field(hadronicMode);
73  zap['d'] = make_field(debug) = 1;
74 
75  zap(argc, argv);
76  }
77  catch(const exception &error) {
78  FATAL(error.what() << endl);
79  }
80 
81  if(hadronicMode) {
82  NOTICE("Plotting hadronic showers." << endl);
83  } else {
84  NOTICE("Plotting EM showers." << endl);
85  }
86 
87 
88  JDetector detector;
89 
90  try {
91  load(detectorFile, detector);
92  }
93  catch(const JException& error) {
94  FATAL(error);
95  }
96 
97  const JPosition3D center = get<JPosition3D>(getHeader(inputFile));
98 
99  NOTICE("Apply detector offset " << center << endl);
100 
101  detector -= center;
102 
103  const JPMTRouter pmtRouter(detector);
104 
105  const double P_atm = NAMESPACE::getAmbientPressure();
106  const double wmin = getMinimalWavelength();
107  const double wmax = getMaximalWavelength();
108 
109  const JDispersion dispersion(P_atm);
110 
111  const double ng[] = { dispersion.getIndexOfRefractionGroup(wmax),
112  dispersion.getIndexOfRefractionGroup(wmin) };
113 
114  typedef JHistogram1D_t::abscissa_type abscissa_type;
115 
116  typedef JTransformableMultiHistogram<JHistogram1D_t,
117  JMAPLIST<JHistogramMap_t,
118  JHistogramMap_t,
119  JHistogramGridMap_t,
120  JHistogramGridMap_t>::maplist> JMultiHistogram_t;
121 
122  typedef JPDFTransformer<4, abscissa_type> JFunction4DTransformer_t;
123 
124  JMultiHistogram_t h0; // occurrence rate of PMT (used for normalisation)
125  JMultiHistogram_t h1; // light from cascade
126 
127  h1.transformer.reset(new JFunction4DTransformer_t(21.5, 2, ng[0], 0.0, JGeant(JGeanx(1.00, -2.2)), 1e-2, NAMESPACE::getAngularAcceptance, 0.05));
128 
129  set<double> C; // cosine emission angle
130 
131  JQuadrature qeant(-1.0, 1.0, 20, JGeanx(1.00, -2.2));
132 
133  for (JQuadrature::const_iterator i = qeant.begin(); i != qeant.end(); ++i) {
134  C.insert(i->getX());
135  }
136 
137  C.insert(-1.01);
138  C.insert(-1.00);
139  C.insert( 1.00);
140  C.insert( 1.01);
141 
142  const double R[] = {0.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, 60.0, 70.0, 80.0};
143  const double Dmax_m = 80.0;
144 
145  for (int i = 0; i != sizeof(R)/sizeof(R[0]); ++i) {
146 
147  const double R_m = R[i];
148 
149  for (set<double>::const_iterator c = C.begin(); c != C.end(); ++c) {
150 
151  const double cd = *c;
152 
153  const double grid = 10.0 + 0.0 * R_m/100.0; // [deg]
154  const double alpha = 2.0 * sqrt(1.0 - cos(grid * PI / 180.0)); // azimuth angle unit step size
155 
156  const int number_of_theta_points = max(2, (int) (180.0/(1.4 * grid)));
157  const double theta_step = PI / (number_of_theta_points + 1);
158 
159  for (double theta = -0.5*theta_step; theta < PI + theta_step; theta += theta_step) {
160 
161  const int number_of_phi_points = max(2, (int) (PI * sin(theta) / alpha));
162  const double phi_step = PI / (number_of_phi_points + 1);
163 
164  for (double phi = -0.5*phi_step; phi < PI + phi_step; phi += phi_step) {
165 
166  for (JMultiHistogram_t* buffer[] = { &h0, &h1, NULL }, **histogram = buffer; *histogram != NULL; ++histogram) {
167  (**histogram)[R_m][cd][theta][phi];
168  }
169  }
170  }
171  }
172  }
173 
174  double buffer[] = { 0.0, 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 };
175  for (JMultiHistogram_t::super_iterator
176  i1 = h1.super_begin(); i1 != h1.super_end(); ++i1) {
177 
178  for (int j = 0; j != sizeof(buffer)/sizeof(buffer[0]); ++j) {
179  i1.getValue()[buffer[j]];
180  }
181  }
182 
183  while (inputFile.hasNext()) {
184 
185  NOTICE("event: " << setw(10) << inputFile.getCounter() << '\r'); STATUS(endl);
186 
187  const Evt* event = inputFile.next();
188  Trk track;
189 
190  if(count_electrons(*event) && !hadronicMode) {
191  track = get_electron(*event);
192  }
193  else if (count_hadrons(*event) && hadronicMode){
194  track = get_hadronic_cascade(*event);
195  }
196  else {
197  WARNING("No electron/hadrons found. Perhaps you meant to run this program in hadronic mode." << endl);
198  continue;
199  }
200 
201  double E = track.E;
202  JPosition3D pos = getPosition (track);
203  JDirection3D dir = getDirection(track);
204  double t0 = track.t;
205 
206  const JRotation3D R(dir);
207 
208  pos.rotate(R);
209 
210  for (vector<Hit>::const_iterator hit = event->mc_hits.begin(); hit != event->mc_hits.end(); ++hit) {
211 
212  try {
213 
214  JAxis3D axis = pmtRouter.getPMT(hit->pmt_id);
215  axis.transform(R, pos);
216 
217  if ((!hadronicMode && from_electron(*hit)) || (hadronicMode && from_hadron(*hit))) {
218 
219  const double t1 = t0 + axis.getLength() * getInverseSpeedOfLight() * getIndexOfRefraction();
220 
221  const double D_m = axis.getLength();
222  const double cd = axis.getZ()/D_m;
223  const double theta = axis.getTheta();
224  const double phi = fabs(axis.getPhi());
225  const double dt = getTime(*hit) - t1;
226  const double npe = getNPE(*hit);
227 
228  if (D_m < Dmax_m) {
229  h1.fill(D_m, cd, theta, phi, dt, npe/E);
230  }
231  }
232  }
233  catch(const exception& error) {
234  FATAL(error.what() << endl);
235  }
236  }
237 
238  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
239 
240  JPosition3D P = module->getPosition();
241 
242  P.rotate(R);
243 
244  P.sub(pos);
245 
246  if (P.getLength() < h0.getXmax()) {
247 
248  for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
249 
250  JAxis3D axis = *pmt;
251  axis.transform(R, pos);
252 
253  h0.fill(axis.getLength(), axis.getZ()/axis.getLength(), axis.getTheta(), fabs(axis.getPhi()), 0.0, 1.0);
254  }
255  }
256  }
257  }
258 
259  double integral = 0;
260  for (JMultiHistogram_t::super_iterator i = h0.super_begin(); i != h0.super_end(); ++i) {
261  integral+=i.getValue().getIntegral();
262  }
263  DEBUG("Integral:\t" << integral << endl);
264 
265  // output
266 
267  JFileStreamWriter out(outputFile.c_str());
268 
269  NOTICE("Storing, " << flush);
270 
271  for (const JMultiHistogram_t* buffer[] = { &h0, &h1, NULL }, **i = buffer; *i != NULL; ++i) {
272  out.store(**i);
273  }
274 
275  out.close();
276  NOTICE("done." << endl);
277 }
int count_electrons(const Evt &evt)
Count the number of electrons in a given event.
Utility class to parse command line options.
Definition: JParser.hh:1410
#define WARNING(A)
Definition: JMessage.hh:63
double getIndexOfRefraction()
Get average index of refraction of water.
Definition: JConstants.hh:111
#define STATUS(A)
Definition: JMessage.hh:61
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 getMinimalWavelength()
Get minimal wavelength for PDF evaluations.
Definition: JPDFToolkit.hh:38
double getTime(const Hit &hit)
Get true time of hit.
string outputFile
Head getHeader(const JMultipleFileScanner_t &file_list)
Get Monte Carlo header.
JLimit JLimit_t
Type definition of limit.
Definition: JLimit.hh:214
const double getInverseSpeedOfLight()
Get inverse speed of light.
Definition: JConstants.hh:100
double getMaximalWavelength()
Get maximal wavelength for PDF evaluations.
Definition: JPDFToolkit.hh:49
Trk get_hadronic_cascade(const Evt &evt)
Auxiliary function to get average direction and total energy of a hadronic cascade.
const Trk & get_electron(const Evt &evt)
Get first electron from the event tracklist.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1836
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.
void load(const JString &file_name, JDetector &detector)
Load detector from input file.
int debug
debug level
Definition: JSirene.cc:59
bool from_hadron(const Hit &hit)
Test whether given hit was produced by a hadronic shower.
#define FATAL(A)
Definition: JMessage.hh:65
collection_type::abscissa_type abscissa_type
JDirection3D getDirection(const Vec &v)
Get direction.
double getNPE(const Hit &hit)
Get true charge of hit.
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:72
static const double C
Speed of light in vacuum [m/ns].
Definition: JConstants.hh:22
bool from_electron(const Hit &hit)
Test whether given hit was produced by an electron.
int count_hadrons(const Evt &evt)
Count the number of hadrons in a given event (not including neutral pions)
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:60
JPosition3D getPosition(const Vec &v)
Get position.
double getAngularAcceptance(const double x)
Angular acceptence of PMT.
Definition: JDrawLED.cc:84