54 typedef JParallelFileScanner_t::multi_pointer_type multi_pointer_type;
57 JParallelFileScanner_t inputFile;
65 JRange<double> pos_grid;
66 JRange<double> time_grid;
69 size_t numberOfPrefits;
78 zap[
'F'] =
make_field(pdfFile) =
"PDF_ShowerPositionFit.root";
79 zap[
'n'] =
make_field(numberOfEvents) = JLimit::max();
83 zap[
'P'] =
make_field(pos_grid) = JRange<double>(-28, 28);
85 zap[
'T'] =
make_field(time_grid) = JRange<double>(-60, 60);
90 if (zap.
read(argc, argv) != 0)
return 1;
92 catch(
const exception& error){
93 FATAL(error.what() << endl);
106 catch(
const JException& error) {
110 TFile* pdf_file =
new TFile(pdfFile.c_str());
111 TH2D* hpdf = (TH2D*)pdf_file->Get(
"hPDF2Dist");
127 JRegressor_t::MAXIMUM_ITERATIONS = 10000;
129 JRegressor_t fit0(
true);
131 JRegressor_t fit(hpdf,
false);
135 STATUS(
"event: " << setw(10) << inputFile.getCounter() <<
'\r');
DEBUG(endl);
137 multi_pointer_type ps = inputFile.next();
143 JEvt::iterator __end = evt->end();
144 if (numberOfPrefits > 0) {
145 advance(__end = evt->begin(), min(numberOfPrefits, evt->size()));
147 partial_sort(evt->begin(), __end, evt->end(),
qualitySorter);
151 copy(evt->begin(), __end, back_inserter(out));
154 JDataL0_t dataL0_selected;
157 buildL0(*tev, moduleRouter,
true, back_inserter(dataL0));
158 buildL2(
JDAQTimeslice(*tev,
false), moduleRouter, back_inserter(data));
160 for (JEvt::const_iterator shower = evt->begin(); shower != __end; ++shower) {
164 for (JDataL0_t::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {
169 const double t_res = hit.getT() - vertex.getT(hit_pos);
170 const JDirection3D photonDir(hit_pos - vertex.getPosition());
171 const double cosT = photonDir.
getDot(hit.getDirection());
175 if(D < Dmax_m && (t_res > -40 && t_res < 40) && (cosT >= -1 && cosT <= 0.1)){
176 dataL0_selected.push_back(hit);
179 if(D < Dmax_m && (t_res > -50 && t_res < 50) && (cosT >= -1 && cosT <= 0.1)){
180 dataL0_selected.push_back(hit);
186 if(dataL0_selected.size() != 0){
193 for(
double x = pos_grid.getLowerLimit(); x <= pos_grid.getUpperLimit();
195 for(
double y = pos_grid.getLowerLimit(); y <= pos_grid.getUpperLimit();
197 for(
double z = pos_grid.getLowerLimit(); z <= pos_grid.getUpperLimit();
199 for(
double t = time_grid.getLowerLimit(); t <= time_grid.getUpperLimit();
204 JPosition3D pos_shifted(vertex.getX() + x, vertex.getY() + y, vertex.getZ() + z);
205 JPoint4D point_shifted(pos_shifted, vertex.getT() + t);
207 for(JDataL0_t::const_iterator i = dataL0_selected.begin(); i != dataL0_selected.end(); ++i){
209 chi2 += fit0(point_shifted, hit);
212 LogChi2.push_back(chi2);
213 start_positions.push_back(point_shifted);
220 while(LogChi2.size() > 35){
222 it = max_element(LogChi2.begin(), LogChi2.end());
223 int pos = it - LogChi2.begin();
224 LogChi2.erase(LogChi2.begin() + pos);
225 start_positions.erase(start_positions.begin() + pos);
230 double simplex_step = 1;
237 chi2 = fit0(
JPoint4D(*start_pos), dataL0_selected.begin(), dataL0_selected.end());
240 JPoint4D pt(fitPos0, fit0.value.getT());
249 chi2 = fit(pt, dataL0_selected.begin(), dataL0_selected.end());
251 int NDF = dataL0_selected.size() - fit.step.size();
257 out.rbegin()->setW(13, chi2);
258 out.rbegin()->setW(14, NDF);
261 DEBUG(
"Too few hits " << endl);