Jpp  pmt_effective_area_update
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JMuonPath.cc
Go to the documentation of this file.
1 
2 #include <string>
3 #include <iostream>
4 #include <iomanip>
5 #include <vector>
6 #include <algorithm>
7 
11 #include "JDAQ/JDAQEventIO.hh"
12 #include "JDAQ/JDAQTimesliceIO.hh"
14 
15 #include "JDetector/JDetector.hh"
18 
19 #include "JTrigger/JHit.hh"
20 #include "JTrigger/JFrame.hh"
21 #include "JTrigger/JTimeslice.hh"
22 #include "JTrigger/JHitL0.hh"
23 #include "JTrigger/JBuildL0.hh"
25 
29 #include "JSupport/JSupport.hh"
30 #include "JSupport/JMeta.hh"
31 
32 #include "JFit/JPMTW0.hh"
33 #include "JFit/JLine3Z.hh"
34 #include "JFit/JGandalf.hh"
35 #include "JFit/JLine3ZRegressor.hh"
36 #include "JFit/JFitToolkit.hh"
37 #include "JFit/JModel.hh"
38 
40 #include "JReconstruction/JEvt.hh"
42 
43 #include "JPhysics/JConstants.hh"
44 #include "JTools/JFunction1D_t.hh"
46 
47 #include "JPhysics/JPDFTable.hh"
48 #include "JPhysics/JPDFToolkit.hh"
49 
50 #include "Jeep/JParser.hh"
51 #include "Jeep/JMessage.hh"
52 
53 
54 /**
55  * \file
56  *
57  * Program to perform JFIT::JRegressor<JLine3Z,JGandalf> fit with I/O of JFIT::JEvt data.
58  * \author mdejong
59  */
60 int main(int argc, char **argv)
61 {
62  using namespace std;
63  using namespace JPP;
64  using namespace KM3NETDAQ;
65 
66  typedef JParallelFileScanner< JTypeList<JDAQEvent, JEvt> > JParallelFileScanner_t;
67  typedef JParallelFileScanner_t::multi_pointer_type multi_pointer_type;
69 
70  JParallelFileScanner_t inputFile;
72  JLimit_t& numberOfEvents = inputFile.getLimit();
73  string detectorFile;
74  string pdfFile;
75  double roadWidth_m;
76  double R_Hz;
77  size_t numberOfPrefits;
78  double TTS_ns;
79  double E_GeV;
80  int debug;
81 
82  try {
83 
84  JParser<> zap("Program to perform fit of muon path to data.");
85 
86  zap['f'] = make_field(inputFile);
87  zap['o'] = make_field(outputFile) = "path.root";
88  zap['a'] = make_field(detectorFile);
89  zap['n'] = make_field(numberOfEvents) = JLimit::max();
90  zap['P'] = make_field(pdfFile);
91  zap['R'] = make_field(roadWidth_m) = numeric_limits<double>::max();
92  zap['B'] = make_field(R_Hz) = 6.0e3;
93  zap['N'] = make_field(numberOfPrefits) = 1;
94  zap['T'] = make_field(TTS_ns);
95  zap['E'] = make_field(E_GeV) = 1.0e3;
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 
105  cout.tie(&cerr);
106 
107 
109 
110  try {
111  load(detectorFile, detector);
112  }
113  catch(const JException& error) {
114  FATAL(error);
115  }
116 
117  const JModuleRouter router(detector);
118 
119 
120  typedef vector<JHitL0> JDataL0_t;
121  const JBuildL0<JHitL0> buildL0;
122 
123 
124  typedef JRegressor<JLine3Z, JGandalf> JRegressor_t;
125 
127  JRegressor_t::T_ns.setRange(-50.0, +50.0); // ns
128  JRegressor_t::Vmax_npe = 10.0; // npe
129  JRegressor_t::MAXIMUM_ITERATIONS = 10000;
130 
131  JRegressor_t fit(pdfFile, TTS_ns);
132 
133  for (int i = 0; i != JRegressor_t::NUMBER_OF_PDFS; ++i) {
134  fit.npe[i].transform(JRegressor_t::JNPE_t::transformer_type::getDefaultTransformer()); // get rid of R-dependent weight function
135  }
136 
137  fit.estimator.reset(new JMEstimatorLorentzian());
138 
139  fit.parameters.resize(4);
140 
141  fit.parameters[0] = JLine3Z::pX();
142  fit.parameters[1] = JLine3Z::pY();
143  fit.parameters[2] = JLine3Z::pDX();
144  fit.parameters[3] = JLine3Z::pDY();
145 
146  if (fit.getRmax() < roadWidth_m) {
147 
148  roadWidth_m = fit.getRmax();
149 
150  WARNING("Set road width to [m] " << roadWidth_m << endl);
151  }
152 
153  const double Rmax_m = 100.0; //!< hit selection [m]
154  const double Tmax_ns = 10.0; //!< hit selection [ns]
155 
156 
157  outputFile.open();
158 
159  outputFile.put(JMeta(argc, argv));
160 
161  while (inputFile.hasNext()) {
162 
163  STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);
164 
165  multi_pointer_type ps = inputFile.next();
166 
167  const JDAQEvent* tev = ps;
168  const JEvt* in = ps;
169 
170  JEvt cp = *in;
171  JEvt out;
172 
173  cp.select(numberOfPrefits, qualitySorter);
174 
175 
176  JDataL0_t dataL0;
177 
178  buildL0(*tev, router, true, back_inserter(dataL0));
179 
180 
181  for (JEvt::const_iterator track = cp.begin(); track != cp.end(); ++track) {
182 
183  const JRotation3D R (getDirection(*track));
184  const JLine1Z tz(getPosition (*track).rotate(R), track->getT());
185  const JModel<JLine1Z> match(tz, roadWidth_m, JRegressor_t::T_ns);
186 
187  // hit selection based on start value
188 
190 
192 
193  for (JDataL0_t::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {
194 
195  JHitW0 hit(*i, R_Hz);
196 
197  hit.rotate(R);
198 
199  if (match(hit)) {
200 
201  top.insert(hit.getPMTIdentifier());
202 
203  const double t1 = hit.getT() - tz.getT(hit);
204 
205  if (tz.getDistance(hit) <= Rmax_m && t1 >= -Tmax_ns && t1 <= +Tmax_ns) {
206  Z.include(hit.getZ() - tz.getDistance(hit) / getTanThetaC());
207  }
208  }
209  }
210 
211 
212  vector<JPMTW0> buffer;
213 
214  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
215 
216  JPosition3D pos(module->getPosition());
217 
218  pos.rotate(R);
219 
220  if (tz.getDistance(pos) <= roadWidth_m && Z(pos.getZ() - tz.getDistance(pos) / getTanThetaC())) {
221 
222  for (unsigned int i = 0; i != module->size(); ++i) {
223 
224  const JDAQPMTIdentifier id(module->getID(), i);
225 
226  JPMT pmt(module->getPMT(i));
227 
228  pmt.rotate(R);
229 
230  buffer.push_back(JPMTW0(pmt, R_Hz, top.count(id)));
231  }
232  }
233  }
234 
235 
236  const int NDF = buffer.size() - fit.parameters.size();
237 
238  if (NDF >= 0) {
239 
240  // set fit parameters
241 
242  if (track->getE() > 0.1)
243  fit.E_GeV = track->getE();
244  else
245  fit.E_GeV = E_GeV;
246 
247  const double chi2 = fit(JLine3Z(tz), buffer.begin(), buffer.end());
248 
249  JTrack3D tb(fit.value);
250 
251  tb.rotate_back(R);
252 
253  out.push_back(getFit(JMUONPATH, tb, getQuality(chi2), NDF));
254 
255  out.rbegin()->setW(track->getW());
256  }
257  }
258 
259  // apply default sorter
260 
261  sort(out.begin(), out.end(), qualitySorter);
262 
263  outputFile.put(out);
264  }
265  STATUS(endl);
266 
268 
269  io >> outputFile;
270 
271  outputFile.close();
272 }
Auxiliary methods to evaluate Poisson probabilities and chi2.
Data regression method for JFIT::JLine3Z.
Auxiliary class for ROOT I/O of application specific meta data.
Definition: JMeta.hh:70
double getT() const
Get calibrated time of hit.
Object writing to file.
Utility class to parse command line options.
Definition: JParser.hh:1500
General exception.
Definition: JException.hh:23
#define WARNING(A)
Definition: JMessage.hh:65
int main(int argc, char *argv[])
Definition: Main.cc:15
ROOT TTree parameter settings of various packages.
Auxiliary methods for PDF calculations.
static const int JMUONPATH
range_type include(argument_type x)
Include given value to range.
Definition: JRange.hh:398
double getQuality(const double chi2, const int NDF)
Get quality of fit.
Template specialisation of L0 builder for JHitL0 data type.
Definition: JBuildL0.hh:102
const JDAQPMTIdentifier & getPMTIdentifier() const
Get PMT identifier.
Auxiliary class for handling PMT geometry, rate and response.
Definition: JPMTW0.hh:22
#define STATUS(A)
Definition: JMessage.hh:63
Detector data structure.
Definition: JDetector.hh:81
Recording of objects on file according a format that follows from the file name extension.
Router for direct addressing of module data in detector data structure.
Rotation matrix.
Definition: JRotation3D.hh:111
General purpose class for parallel reading of objects from a single file or multiple files...
Template specialisation of class JModel to match hit with muon trajectory along z-axis.
Definition: JFit/JModel.hh:34
Data structure for fit of straight line in positive z-direction.
Definition: JLine3Z.hh:36
Basic data structure for time and time over threshold information of hit.
string outputFile
Data structure for detector geometry and calibration.
Various implementations of functional maps.
Basic data structure for L0 hit.
Type list.
Definition: JTypeList.hh:22
Scanning of objects from a single file according a format that follows from the extension of each fil...
Auxiliary class for defining the range of iterations of objects.
Definition: JLimit.hh:41
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
Parallel scanning of objects from a single file or multiple files according a format that follows fro...
Data structure for PMT geometry, calibration and status.
Definition: JPMT.hh:42
ROOT I/O of application specific meta data.
JPosition3D getPosition(const Vec &pos)
Get position.
Physics constants.
int debug
debug level
Definition: JSirene.cc:63
then usage $script[distance] fi case set_variable R
Definition: JDrawLED.sh:43
General purpose messaging.
#define FATAL(A)
Definition: JMessage.hh:67
Direct access to module in detector data structure.
JFit getFit(const JHistory &history, const JTrack3D &track, const double Q, const int NDF, const double energy=0.0, const int status=0)
Get fit.
Auxiliary class for a hit with background rate value.
Definition: JHitW0.hh:21
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
Regressor function object for JLine3Z fit using JGandalf minimiser.
Utility class to parse command line options.
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 fit of straight line paralel to z-axis.
Definition: JLine1Z.hh:27
then cp
double getTanThetaC()
Get average tangent of Cherenkov angle of water corresponding to group velocity.
Object reading from a list of files.
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
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:40
JPosition3D & rotate(const JRotation3D &R)
Rotate.
Definition: JPosition3D.hh:186
double getZ() const
Get z position.
Definition: JVector3D.hh:115
bool qualitySorter(const JRECONSTRUCTION::JFit &first, const JRECONSTRUCTION::JFit &second)
Comparison of fit results.
Lorentzian M-estimator.
Definition: JMEstimator.hh:79