Jpp - the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JReconstruction/JShowerNPE.cc
Go to the documentation of this file.
1 #include <string>
2 #include <iostream>
3 #include <iomanip>
4 #include <vector>
5 #include <algorithm>
6 
9 
10 #include "JDAQ/JDAQEventIO.hh"
12 
13 #include "JDetector/JDetector.hh"
16 
17 #include "JLang/JPredicate.hh"
18 
19 #include "JTrigger/JHitL0.hh"
20 #include "JTrigger/JBuildL0.hh"
21 
22 #include "JSupport/JSupport.hh"
25 
26 #include "JPhysics/JNPE_t.hh"
27 #include "JPhysics/JPDF_t.hh"
28 
29 #include "JFit/JPoint4D.hh"
30 #include "JFit/JModel.hh"
31 
33 #include "JReconstruction/JEvt.hh"
36 
37 #include "Jeep/JPrint.hh"
38 #include "Jeep/JParser.hh"
39 #include "Jeep/JMessage.hh"
40 
41 #include "JTools/JResult.hh"
42 
43 namespace {
44 
45  /**
46  * Auxiliary data structure for sorting of hits.
47  */
48  static const struct {
49  /**
50  * Compare hits by PMT identifier and time.
51  *
52  * \param first first hit
53  * \param second second hit
54  * \return true if first before second; else false
55  */
56  template<class T>
57  bool operator()(const T& first, const T& second) const
58  {
59  using namespace std;
60  using namespace JPP;
61  using namespace KM3NETDAQ;
62 
63  if (equal_to<JDAQPMTIdentifier>()(first, second))
64  return less<JHit>()(first, second);
65  else
66  return less<JDAQPMTIdentifier>()(first, second);
67  }
68  } compare;
69 }
70 
71 
72 /**
73  * \file
74  *
75  * Program to evaluate hit probabilities.
76  *
77  * \author mdejong and adomi
78  */
79 int main(int argc, char **argv)
80 {
81  using namespace std;
82  using namespace JPP;
83  using namespace KM3NETDAQ;
84 
85  typedef JParallelFileScanner< JTypeList<JDAQEvent, JEvt> > JParallelFileScanner_t;
86  typedef JParallelFileScanner_t::multi_pointer_type multi_pointer_type;
87 
88  JParallelFileScanner_t inputFile;
89  JLimit_t& numberOfEvents = inputFile.getLimit();
90  string detectorFile;
91  string pdfFile;
93  int debug;
94 
95  try {
96 
97  parameters.numberOfPrefits = 1;
98 
99  JParser<> zap("Program to perform detector calibration using reconstructed shower direction.");
100 
101  zap['f'] = make_field(inputFile);
102  zap['a'] = make_field(detectorFile);
103  zap['n'] = make_field(numberOfEvents) = JLimit::max();
104  zap['P'] = make_field(pdfFile);
106  zap['d'] = make_field(debug) = 1;
107 
108  zap(argc, argv);
109  }
110  catch(const exception& error) {
111  FATAL(error.what() << endl);
112  }
113 
114  //setDAQLongprint(debug >= JEEP::debug_t);
115 
116  cout.tie(&cerr);
117 
118  if (parameters.numberOfPrefits != 1) {
119  WARNING("Number of prefits " << parameters.numberOfPrefits << " != " << 1 << endl);
120  }
121 
123 
124  try {
125  load(detectorFile, detector);
126  }
127  catch(const JException& error) {
128  FATAL(error);
129  }
130 
131  const JModuleRouter router(detector);
132 
133  JSummaryFileRouter summary(inputFile, parameters.R_Hz);
134 
135  const JShowerNPE_t npe(pdfFile);
136 
137  const JTimeRange T_ns(parameters.TMin_ns, parameters.TMax_ns);
138 
139  typedef vector<JHitL0> JDataL0_t;
140  typedef vector<JHitW0> JDataW0_t;
141 
142  const JBuildL0<JHitL0> buildL0;
143 
144  while (inputFile.hasNext()) {
145 
146  STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);
147 
148  multi_pointer_type ps = inputFile.next();
149 
150  JDAQEvent* tev = ps;
151  JEvt* in = ps;
152 
153  DEBUG("event: " << *tev << endl);
154 
155  summary.update(*tev);
156 
157  in->select(parameters.numberOfPrefits, qualitySorter);
158 
159  JDataL0_t dataL0;
160 
161  buildL0(*tev, router, true, back_inserter(dataL0));
162 
163  for (JEvt::const_iterator shower = in->begin(); shower != in->end(); ++shower) {
164 
165  DEBUG("shower: " << *shower << endl);
166 
167  JRotation3D R (getDirection(*shower));
168  JPoint4D vx(getPosition(*shower).rotate(R), shower->getT());
169 
170  const JModel<JPoint4D> match(vx, parameters.roadWidth_m, T_ns);
171 
172  // hit selection based on start value
173 
174  JDataW0_t data;
175 
176  for (JDataL0_t::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {
177 
178  double rate_Hz = summary.getRate(*i);
179 
180  if (rate_Hz <= 0.0) {
181  rate_Hz = summary.getRate();
182  }
183 
184  JHitW0 hit(*i, rate_Hz);
185 
186  hit.rotate(R);
187 
188  if (match(hit)) {
189  data.push_back(hit);
190  }
191  }
192 
193  // select first hit
194 
195  sort(data.begin(), data.end(), compare);
196 
197  JDataW0_t::iterator __end = unique(data.begin(), data.end(), equal_to<JDAQPMTIdentifier>());
198 
199  DEBUG("point4d: "
200  << FIXED(12,3) << vx.getX() << ' '
201  << FIXED(12,3) << vx.getY() << ' '
202  << FIXED(12,3) << vx.getZ() << ' '
203  << FIXED(12,3) << vx.getT() << ' '
204  << FIXED( 8,3) << shower->getQ() << endl);
205 
206  double Q = 0.0;
207 
208  for (JDataW0_t::const_iterator hit = data.begin(); hit != __end; ++hit) {
209 
210  JPosition3D D(hit->getPosition());
211  D.sub(vx);
212 
213  const double x = hit->getX() - vx.getX();
214  const double y = hit->getY() - vx.getY();
215  const double z = hit->getZ() - vx.getZ();
216  const double cd = z/D.getLength(); // Delta = angle between shower direction and PMT position
217 
218  const double t1 = vx.getT() + (D.getLength() * getIndexOfRefraction() * getInverseSpeedOfLight());
219 
220  JDirection3D u(hit->getDX(), hit->getDY(), hit->getDZ()); // PMT orientation
221 
222  u.rotate(JRotation3Z(-atan2(y,x))); // rotate PMT axis to x-z plane
223 
224  const double theta = u.getTheta();
225  const double phi = fabs(u.getPhi()); // rotational symmetry of Cherenkov cone
226 
227  const double E = shower->getE();
228  const double dt = T_ns.constrain(hit->getT() - t1);
229 
230  double H1 = npe.calculate(E, D.getLength(), cd, theta, phi);
231  double H0 = hit->getR() * 1e-9 * T_ns.getLength();
232 
233  double Vmax_npe = 20.0;
234  if (H1 >= Vmax_npe) {
235  H1 *= Vmax_npe / H1;
236  }
237 
238  H1 += H0; // signal + background
239 
240  const double chi2 = getChi2(H1, true); // -log(lik)
241 
242  DEBUG("hit: "
243  << setw(10) << hit->getModuleID() << ':' << setw( 2) << setfill('0') << hit->getPMTAddress() << setfill(' ') << ' '
244  << FIXED(12,1) << E << ' '
245  << FIXED( 9,1) << R << ' '
246  << FIXED( 6,4) << theta << ' '
247  << FIXED( 6,4) << phi << ' '
248  << FIXED( 8,3) << dt << ' '
249  << FIXED(12,3) << chi2 << endl);
250 
251  Q += getQuality(chi2);
252  }
253 
254  DEBUG("quality: " << FIXED(8,3) << Q << ' ' << distance(data.begin(), __end) << endl);
255  }
256  }
257  STATUS(endl);
258 }
Utility class to parse command line options.
Definition: JParser.hh:1500
General exception.
Definition: JException.hh:23
#define WARNING(A)
Definition: JMessage.hh:65
do echo Generating $dir eval D
Definition: JDrawLED.sh:50
Q(UTCMax_s-UTCMin_s)-livetime_s
int main(int argc, char *argv[])
Definition: Main.cc:15
ROOT TTree parameter settings of various packages.
Data structure for direction in three dimensions.
Definition: JDirection3D.hh:33
double getQuality(const double chi2, const int NDF)
Get quality of fit.
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
Template specialisation of L0 builder for JHitL0 data type.
Definition: JBuildL0.hh:102
This include file containes various data structures that can be used as specific return types for the...
double getIndexOfRefraction()
Get average index of refraction of water corresponding to group velocity.
do rm f tmp H1
Data structure for vertex fit.
Definition: JPoint4D.hh:22
#define STATUS(A)
Definition: JMessage.hh:63
void update(const JDAQHeader &header)
Update router.
Detector data structure.
Definition: JDetector.hh:80
Router for direct addressing of module data in detector data structure.
Rotation matrix.
Definition: JRotation3D.hh:111
double getRate() const
Get default rate.
Auxiliary data structure for shower PDF.
Definition: JNPE_t.hh:167
Auxiliary data structure for muon PDF.
*fatal Wrong number of arguments esac JCookie sh typeset Z DETECTOR typeset Z SOURCE_RUN typeset Z TARGET_RUN set_variable PARAMETERS_FILE $WORKDIR parameters
Definition: diff-Tuna.sh:38
General purpose class for parallel reading of objects from a single file or multiple files...
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition: JParser.hh:66
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:446
Data structure for detector geometry and calibration.
then echo The file $DIR KM3NeT_00000001_00000000 root already please rename or remove it first
Basic data structure for L0 hit.
T getLength() const
Get length (difference between upper and lower limit).
Definition: JRange.hh:289
Rotation around Z-axis.
Definition: JRotation3D.hh:85
Auxiliary class for defining the range of iterations of objects.
Definition: JLimit.hh:41
I/O formatting auxiliaries.
JAxis3D & rotate(const JRotation3D &R)
Rotate axis.
Definition: JAxis3D.hh:225
Detector file.
Definition: JHead.hh:196
JDirection3D getDirection(const Vec &dir)
Get direction.
Acoustic event fit.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1961
do set_variable OUTPUT_DIRECTORY $WORKDIR T
Parallel scanning of objects from a single file or multiple files according a format that follows fro...
JPosition3D getPosition(const Vec &pos)
Get position.
T constrain(argument_type x) const
Constrain value to range.
Definition: JRange.hh:350
File router for fast addressing of summary data.
int debug
debug level
Definition: JSirene.cc:63
Auxiliary data structure for muon PDF.
General purpose messaging.
#define FATAL(A)
Definition: JMessage.hh:67
Direct access to module in detector data structure.
Auxiliary class for a hit with background rate value.
Definition: JHitW0.hh:21
then echo Variable JPP_DIR undefined exit fi source $JPP_DIR setenv sh $JPP_DIR set_variable NORTH set_variable EAST set_variable SOUTH set_variable WEST set_variable WORKDIR tmp set_variable R set_variable CT set_variable YMAX set_variable YMIN if do_usage *then usage $script[distance] fi case set_variable R
Definition: JDrawLED.sh:35
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
Data structure for fit parameters.
Utility class to parse command line options.
const double getInverseSpeedOfLight()
Get inverse speed of light.
then usage $script< input_file >< detector_file > fi set_variable OUTPUT_DIR set_variable SELECTOR JDAQTimesliceL1 set_variable DEBUG case set_variable DEBUG
Data structure for position in three dimensions.
Definition: JPosition3D.hh:36
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:73
do set_variable DETECTOR_TXT $WORKDIR detector
double u[N+1]
Definition: JPolint.hh:739
Template specialisation of class JModel to match hit with bright point.
Definition: JFit/JModel.hh:121
double getChi2(const double P)
Get chi2 corresponding to given probability.
Definition: JFitToolkit.hh:70
then fatal Wrong number of arguments fi set_variable DETECTOR $argv[1] set_variable INPUT_FILE $argv[2] eval JPrintDetector a $DETECTOR O IDENTIFIER eval JPrintDetector a $DETECTOR O SUMMARY source JAcoustics sh $DETECTOR_ID CHECK_EXIT_CODE typeset A TRIPODS get_tripods $WORKDIR tripod txt TRIPODS for EMITTER in
Definition: JCanberra.sh:38
JPosition3D & rotate(const JRotation3D &R)
Rotate.
Definition: JPosition3D.hh:186
bool qualitySorter(const JRECONSTRUCTION::JFit &first, const JRECONSTRUCTION::JFit &second)
Comparison of fit results.
double calculate(const double E, const double D, const double cd, const double theta, const double phi) const
Get PDF.
Definition: JNPE_t.hh:241
then usage $script[input file[working directory[option]]] nWhere option can be E
Definition: JMuonPostfit.sh:35