59 int main(
int argc,
char **argv)
65 typedef JParallelFileScanner< JTypeList<JDAQEvent, JEvt> > JParallelFileScanner_t;
66 typedef JParallelFileScanner_t::multi_pointer_type multi_pointer_type;
69 JParallelFileScanner_t inputFile;
76 size_t numberOfPrefits;
83 JParser<> zap(
"Program to perform fit of muon path to data.");
88 zap[
'n'] =
make_field(numberOfEvents) = JLimit::max();
90 zap[
'R'] =
make_field(roadWidth_m) = numeric_limits<double>::max();
99 catch(
const exception& error) {
100 FATAL(error.what() << endl);
111 load(detectorFile, detector);
113 catch(
const JException& error) {
117 const JModuleRouter moduleRouter(detector);
121 const JBuildL0<JHitL0> buildL0;
124 typedef JRegressor<JLine3Z, JGandalf> JRegressor_t;
127 JRegressor_t::T_ns.setRange(-50.0, +50.0);
128 JRegressor_t::Vmax_npe = 10.0;
129 JRegressor_t::MAXIMUM_ITERATIONS = 10000;
131 JRegressor_t fit(pdfFile, TTS_ns);
133 for (
int i = 0; i != JRegressor_t::NUMBER_OF_PDFS; ++i) {
134 fit.npe[i].transform(JRegressor_t::JNPE_t::transformer_type::getDefaultTransformer());
137 fit.estimator.reset(
new JMEstimatorLorentzian());
139 fit.parameters.resize(4);
142 fit.parameters[1] = JLine3Z::pY();
143 fit.parameters[2] = JLine3Z::pDX();
144 fit.parameters[3] = JLine3Z::pDY();
146 if (fit.getRmax() < roadWidth_m) {
148 roadWidth_m = fit.getRmax();
150 WARNING(
"Set road width to [m] " << roadWidth_m << endl);
153 const double Rmax_m = 100.0;
154 const double Tmax_ns = 10.0;
161 while (inputFile.hasNext()) {
163 STATUS(
"event: " << setw(10) << inputFile.getCounter() <<
'\r');
DEBUG(endl);
165 multi_pointer_type ps = inputFile.next();
173 JEvt::iterator __end = evt->
end();
175 if (numberOfPrefits > 0) {
176 advance(__end = evt->begin(), min(numberOfPrefits, evt->size()));
179 partial_sort(evt->begin(), __end, evt->end(),
qualitySorter);
184 buildL0(*tev, moduleRouter,
true, back_inserter(dataL0));
187 if (dataL0.size() >= fit.parameters.size()) {
189 for (JEvt::const_iterator track = evt->begin(); track != __end; ++track) {
193 const JModel<JLine1Z> match(tz, roadWidth_m, JRegressor_t::T_ns);
201 for (JDataL0_t::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {
203 JHitW0 hit(*i, R_Hz);
209 top.insert(hit.getPMTIdentifier());
211 const double t1 = hit.getT() - tz.getT(hit);
213 if (tz.getDistance(hit) <= Rmax_m && t1 >= -Tmax_ns && t1 <= +Tmax_ns) {
214 Z.include(hit.getZ() - tz.getDistance(hit) /
getTanThetaC());
222 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
224 JPosition3D pos(module->getPosition());
228 if (tz.getDistance(pos) <= roadWidth_m && Z(pos.getZ() - tz.getDistance(pos) /
getTanThetaC())) {
230 for (
unsigned int i = 0; i != module->size(); ++i) {
234 JPMT pmt(module->getPMT(i));
238 buffer.push_back(JPMTW0(pmt, R_Hz, top.count(
id)));
244 const int NDF = buffer.size() - fit.parameters.size();
250 if (track->getE() > 0.1)
251 fit.E_GeV = track->getE();
255 const double chi2 = fit(JLine3Z(tz), buffer.begin(), buffer.end());
257 JTrack3D tb(fit.value);
263 out.rbegin()->setW(track->getW());
277 JSingleFileScanner<JRemove<typelist, JEvt>::typelist> io(inputFile);