Jpp  19.1.0
the software that should make you happy
JMuonPDF.cc
Go to the documentation of this file.
1 #include <string>
2 #include <iostream>
3 #include <iomanip>
4 #include <vector>
5 #include <algorithm>
6 
7 #include "JDAQ/JDAQEventIO.hh"
9 
10 #include "JDetector/JDetector.hh"
13 
14 #include "JTrigger/JHitL0.hh"
15 #include "JTrigger/JBuildL0.hh"
16 
17 #include "JSupport/JSupport.hh"
20 
21 #include "JPhysics/JPDF_t.hh"
22 #include "JPhysics/JNPE_t.hh"
23 
24 #include "JFit/JLine1Z.hh"
25 #include "JFit/JModel.hh"
26 
28 #include "JReconstruction/JEvt.hh"
31 
32 #include "Jeep/JPrint.hh"
33 #include "Jeep/JParser.hh"
34 #include "Jeep/JMessage.hh"
35 
36 
37 namespace {
38 
39  /**
40  * Auxiliary data structure for sorting of hits.
41  */
42  static const struct {
43  /**
44  * Compare hits by PMT identifier and time.
45  *
46  * \param first first hit
47  * \param second second hit
48  * \return true if first before second; else false
49  */
50  template<class T>
51  bool operator()(const T& first, const T& second) const
52  {
53  if (first.getPMTIdentifier() == second.getPMTIdentifier())
54  return first.getT() < second.getT();
55  else
56  return first.getPMTIdentifier() < second.getPMTIdentifier();
57  }
58  } compare;
59 }
60 
61 
62 /**
63  * \file
64  *
65  * Program to evaluate hit probabilities.
66  *
67  * \author mdejong
68  */
69 int main(int argc, char **argv)
70 {
71  using namespace std;
72  using namespace JPP;
73  using namespace KM3NETDAQ;
74 
75  typedef JParallelFileScanner< JTypeList<JDAQEvent, JEvt> > JParallelFileScanner_t;
76  typedef JParallelFileScanner_t::multi_pointer_type multi_pointer_type;
77 
78  JParallelFileScanner_t inputFile;
79  JLimit_t& numberOfEvents = inputFile.getLimit();
80  string detectorFile;
81  string pdfFile;
82  JMuonGandalfParameters_t parameters;
83  int debug;
84 
85  try {
86 
87  parameters.numberOfPrefits = 1;
88 
89  JParser<> zap("Program to evaluate hit probabilities.");
90 
91  zap['f'] = make_field(inputFile);
92  zap['a'] = make_field(detectorFile);
93  zap['n'] = make_field(numberOfEvents) = JLimit::max();
94  zap['P'] = make_field(pdfFile);
95  zap['@'] = make_field(parameters) = JPARSER::initialised();
96  zap['d'] = make_field(debug) = 1;
97 
98  zap(argc, argv);
99  }
100  catch(const exception& error) {
101  FATAL(error.what() << endl);
102  }
103 
104  //setDAQLongprint(debug >= JEEP::debug_t);
105 
106 
107  if (parameters.numberOfPrefits != 1) {
108  WARNING("Number of prefits " << parameters.numberOfPrefits << " != " << 1 << endl);
109  }
110 
112 
113  try {
114  load(detectorFile, detector);
115  }
116  catch(const JException& error) {
117  FATAL(error);
118  }
119 
120  const JModuleRouter router(detector);
121 
122  JSummaryFileRouter summary(inputFile, parameters.R_Hz);
123 
124  const JMuonPDF_t pdf(pdfFile, parameters.TTS_ns);
125  const JMuonNPE_t npe(pdfFile);
126 
127  const JTimeRange T_ns(parameters.TMin_ns, parameters.TMax_ns);
128 
129  typedef vector<JHitL0> JDataL0_t;
130  typedef vector<JHitW0> JDataW0_t;
131 
132  const JBuildL0<JHitL0> buildL0;
133 
134  while (inputFile.hasNext()) {
135 
136  STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);
137 
138  multi_pointer_type ps = inputFile.next();
139 
140  JDAQEvent* tev = ps;
141  JEvt* in = ps;
142 
143  DEBUG("event: " << *tev << endl);
144 
145  summary.update(*tev);
146 
147  in->select(parameters.numberOfPrefits, qualitySorter);
148 
149  JDataL0_t dataL0;
150 
151  buildL0(*tev, router, true, back_inserter(dataL0));
152 
153  for (JEvt::const_iterator track = in->begin(); track != in->end(); ++track) {
154 
155  DEBUG("track: " << *track << endl);
156 
157  JRotation3D R (getDirection(*track));
158  JLine1Z tz(getPosition (*track).rotate(R), track->getT());
159  JRange<double> Z_m;
160 
161  if (track->hasW(JSTART_LENGTH_METRES)) {
162  Z_m = JRange<double>(0.0, track->getW(JSTART_LENGTH_METRES)) + JRange<double>(parameters.ZMin_m, parameters.ZMax_m);
163  }
164 
165  const JModel<JLine1Z> match(tz, parameters.roadWidth_m, T_ns, Z_m);
166 
167  // hit selection based on start value
168 
169  JDataW0_t data;
170 
171  for (JDataL0_t::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {
172 
173  JHitW0 hit(*i, summary.getRate(*i));
174 
175  hit.rotate(R);
176 
177  if (match(hit)) {
178  data.push_back(hit);
179  }
180  }
181 
182  // select first hit in PMT
183 
184  sort(data.begin(), data.end(), compare);
185 
186  JDataW0_t::iterator __end = unique(data.begin(), data.end(), equal_to<JDAQPMTIdentifier>());
187 
188 
189  double E_GeV = parameters.E_GeV;
190  /*
191  if (track->getE() > 0.1) {
192  E_GeV = track->getE();
193  }
194  */
195 
196  DEBUG("line1z: "
197  << FIXED(12,3) << tz.getX() << ' '
198  << FIXED(12,3) << tz.getY() << ' '
199  << FIXED(12,3) << tz.getZ() << ' '
200  << FIXED(12,3) << tz.getT() << ' '
201  << FIXED(12,1) << E_GeV << ' '
202  << FIXED( 8,3) << track->getQ() << endl);
203 
204  // move track to bary-center of hits
205  /*
206  int N = 0;
207  double z = 0.0;
208 
209  for (JDataW0_t::const_iterator hit = data.begin(); hit != __end; ++hit) {
210  N += 1;
211  z += hit->getZ();
212  }
213 
214  if (N != 0) {
215  tz.setZ(z/N, getSpeedOfLight());
216  }
217  */
218 
219  double Q = 0.0;
220 
221  for (JDataW0_t::const_iterator hit = data.begin(); hit != __end; ++hit) {
222 
223  const double x = hit->getX() - tz.getX();
224  const double y = hit->getY() - tz.getY();
225  const double z = hit->getZ() - tz.getZ();
226  const double R = sqrt(x*x + y*y);
227 
228  const double t1 = tz.getT() + (z + R * getTanThetaC()) * getInverseSpeedOfLight();
229 
230  JDirection3D u(hit->getDX(), hit->getDY(), hit->getDZ()); // PMT orientation
231 
232  u.rotate(JRotation3Z(-atan2(y,x))); // rotate PMT axis to x-z plane
233 
234  const double theta = u.getTheta();
235  const double phi = fabs(u.getPhi()); // rotational symmetry of Cherenkov cone
236 
237  //const double E = gWater.getE(E_GeV, z); // correct for energy loss
238  const double E = E_GeV;
239  const double dt = T_ns.constrain(hit->getT() - t1);
240 
241  JMuonPDF_t::result_type H1 = pdf.calculate(E, R, theta, phi, dt);
242  JMuonPDF_t::result_type H0(hit->getR() * 1e-9, 0.0, T_ns);
243 
244  H1 += H0; // signal + background
245 
246  const double chi2 = H1.getChi2() - H0.getChi2(); // likelihood ratio
247 
248  DEBUG("hit: "
249  << setw(10) << hit->getModuleID() << ':' << setw( 2) << setfill('0') << hit->getPMTAddress() << setfill(' ') << ' '
250  << FIXED(12,1) << E << ' '
251  << FIXED( 9,1) << R << ' '
252  << FIXED( 6,4) << theta << ' '
253  << FIXED( 6,4) << phi << ' '
254  << FIXED( 8,3) << dt << ' '
255  << FIXED(12,3) << chi2 << endl);
256 
257  Q += getQuality(chi2);
258  }
259 
260  DEBUG("quality: " << FIXED(8,3) << Q << ' ' << distance(data.begin(), __end) << endl);
261 
262 
263  double Y = 0.0;
264 
265  for (JDetector::const_iterator i = detector.begin(); i != detector.end(); ++i) {
266 
267  JPosition3D pos(i->getPosition());
268 
269  pos.transform(R, tz.getPosition());
270 
271  if (pos.getX() <= parameters.roadWidth_m) {
272 
273  JModule module = *i;
274 
275  module.transform(R, tz.getPosition());
276 
277  for (JModule::iterator pmt = module.begin(); pmt != module.end(); ++pmt) {
278 
279  const double R = sqrt(pmt->getX()*pmt->getX() + pmt->getY()*pmt->getY());
280  const double theta = pmt->getTheta();
281  const double phi = fabs(pmt->getPhi());
282  const double y = npe.calculate(1.0, R, theta, phi); // MIP
283 
284  DEBUG("PMT: "
285  << setw(10) << module.getID() << ':' << setw( 2) << setfill('0') << distance(module.begin(),pmt) << setfill(' ') << ' '
286  << FIXED(9,1) << R << ' '
287  << FIXED(6,4) << theta << ' '
288  << FIXED(6,4) << phi << ' '
289  << SCIENTIFIC(7,2) << y << endl);
290 
291  Y += y;
292  }
293  }
294  }
295 
296  DEBUG("npe: " << FIXED(7,1) << Y << endl);
297  }
298  }
299  STATUS(endl);
300 }
Data structure for detector geometry and calibration.
Basic data structure for L0 hit.
General purpose messaging.
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
#define STATUS(A)
Definition: JMessage.hh:63
#define FATAL(A)
Definition: JMessage.hh:67
int debug
debug level
Definition: JSirene.cc:69
#define WARNING(A)
Definition: JMessage.hh:65
Direct access to module in detector data structure.
int main(int argc, char **argv)
Definition: JMuonPDF.cc:69
Auxiliary data structure for muon PDF.
Auxiliary data structure for muon PDF.
Parallel scanning of objects from a single file or multiple files according a format that follows fro...
Utility class to parse command line options.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:2142
I/O formatting auxiliaries.
ROOT TTree parameter settings of various packages.
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
Detector data structure.
Definition: JDetector.hh:96
Router for direct addressing of module data in detector data structure.
Data structure for a composite optical module.
Definition: JModule.hh:75
void transform(const JRotation3D &R, const JVector3D &pos)
Transformation of geometry (see method JGEOMETRY3D::JPosition3D::transform(const JRotation3D&,...
Definition: JModule.hh:345
Data structure for fit of straight line paralel to z-axis.
Definition: JLine1Z.hh:29
double getT(const JVector3D &pos) const
Get arrival time of Cherenkov light at given position.
Definition: JLine1Z.hh:114
double getZ(const JPosition3D &pos) const
Get point of emission of Cherenkov light along muon path.
Definition: JLine1Z.hh:134
JAxis3D & rotate(const JRotation3D &R)
Rotate axis.
Definition: JAxis3D.hh:225
Data structure for direction in three dimensions.
Definition: JDirection3D.hh:35
Data structure for position in three dimensions.
Definition: JPosition3D.hh:38
JPosition3D & rotate(const JRotation3D &R)
Rotate.
Definition: JPosition3D.hh:186
const JPosition3D & getPosition() const
Get position.
Definition: JPosition3D.hh:130
void transform(const JRotation3D &R, const JVector3D &pos)
Transform position.
Definition: JPosition3D.hh:331
Rotation matrix.
Definition: JRotation3D.hh:114
Rotation around Z-axis.
Definition: JRotation3D.hh:87
double getY() const
Get y position.
Definition: JVector3D.hh:104
double getX() const
Get x position.
Definition: JVector3D.hh:94
General exception.
Definition: JException.hh:24
int getID() const
Get identifier.
Definition: JObjectID.hh:50
Utility class to parse command line options.
Definition: JParser.hh:1698
Auxiliary class for a hit with background rate value.
Definition: JHitW0.hh:23
General purpose class for parallel reading of objects from a single file or multiple files.
File router for fast addressing of summary data.
void update(const JDAQHeader &header)
Update router.
double getRate() const
Get default rate.
T constrain(argument_type x) const
Constrain value to range.
Definition: JRange.hh:350
Template specialisation of L0 builder for JHitL0 data type.
Definition: JBuildL0.hh:105
static const int JSTART_LENGTH_METRES
distance between first and last hits in metres from JStart.cc
JDirection3D getDirection(const Vec &dir)
Get direction.
JPosition3D getPosition(const Vec &pos)
Get position.
double getQuality(const double chi2, const int NDF)
Get quality of fit.
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
const double getInverseSpeedOfLight()
Get inverse speed of light.
double getTanThetaC()
Get average tangent of Cherenkov angle of water corresponding to group velocity.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
bool qualitySorter(const JFit &first, const JFit &second)
Comparison of fit results.
double u[N+1]
Definition: JPolint.hh:865
KM3NeT DAQ data structures and auxiliaries.
Definition: DataQueue.cc:39
Definition: JSTDTypes.hh:14
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:448
Detector file.
Definition: JHead.hh:227
Acoustic event fit.
Template specialisation of class JModel to match hit with muon trajectory along z-axis.
Definition: JFit/JModel.hh:36
double calculate(const double E, const double R, const double theta, const double phi) const
Get PDF.
Definition: JNPE_t.hh:108
Auxiliary data structure for muon PDF.
Definition: JPDF_t.hh:135
JFunction1D_t::result_type result_type
Definition: JPDF_t.hh:145
result_type calculate(const double E, const double R, const double theta, const double phi, const double t1) const
Get PDF.
Definition: JPDF_t.hh:233
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition: JParser.hh:68
double TMin_ns
minimal time w.r.t. Cherenkov hypothesis [ns]
double TMax_ns
maximal time w.r.t. Cherenkov hypothesis [ns]
Auxiliary class for defining the range of iterations of objects.
Definition: JLimit.hh:45
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:84
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:488