44 int main(
int argc,
char **argv){
53 typedef JParallelFileScanner_t::multi_pointer_type multi_pointer_type;
56 JParallelFileScanner_t inputFile;
68 size_t numberOfPrefits;
77 zap[
'F'] =
make_field(pdfFile) =
"PDF_ShowerPositionFit.root";
78 zap[
'n'] =
make_field(numberOfEvents) = JLimit::max();
89 if (zap.
read(argc, argv) != 0)
return 1;
91 catch(
const exception& error){
92 FATAL(error.what() << endl);
109 TFile* pdf_file =
new TFile(pdfFile.c_str());
110 TH2D* hpdf = (TH2D*)pdf_file->Get(
"hPDF2Dist");
126 JRegressor_t::MAXIMUM_ITERATIONS = 10000;
128 JRegressor_t fit0(
true);
130 JRegressor_t fit(hpdf,
false);
134 STATUS(
"event: " << setw(10) << inputFile.getCounter() <<
'\r');
DEBUG(endl);
136 multi_pointer_type ps = inputFile.next();
142 JEvt::iterator __end = evt->end();
143 if (numberOfPrefits > 0) {
144 advance(__end = evt->begin(), min(numberOfPrefits, evt->size()));
146 partial_sort(evt->begin(), __end, evt->end(),
qualitySorter);
150 copy(evt->begin(), __end, back_inserter(out));
153 JDataL0_t dataL0_selected;
156 buildL0(*tev, moduleRouter,
true, back_inserter(dataL0));
157 buildL2(
JDAQTimeslice(*tev,
false), moduleRouter, back_inserter(data));
159 for (JEvt::const_iterator shower = evt->begin(); shower != __end; ++shower) {
163 for (JDataL0_t::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {
168 const double t_res = hit.
getT() - vertex.getT(hit_pos);
169 const JDirection3D photonDir(hit_pos - vertex.getPosition());
174 if(D < Dmax_m && (t_res > -40 && t_res < 40) && (cosT >= -1 && cosT <= 0.1)){
175 dataL0_selected.push_back(hit);
178 if(D < Dmax_m && (t_res > -50 && t_res < 50) && (cosT >= -1 && cosT <= 0.1)){
179 dataL0_selected.push_back(hit);
185 if(dataL0_selected.size() != 0){
203 JPosition3D pos_shifted(vertex.getX() + x, vertex.getY() + y, vertex.getZ() + z);
204 JPoint4D point_shifted(pos_shifted, vertex.getT() + t);
206 for(JDataL0_t::const_iterator i = dataL0_selected.begin(); i != dataL0_selected.end(); ++i){
208 chi2 += fit0(point_shifted, hit);
211 LogChi2.push_back(chi2);
212 start_positions.push_back(point_shifted);
219 while(LogChi2.size() > 35){
221 it = max_element(LogChi2.begin(), LogChi2.end());
222 int pos = it - LogChi2.begin();
223 LogChi2.erase(LogChi2.begin() + pos);
224 start_positions.erase(start_positions.begin() + pos);
229 double simplex_step = 1;
236 chi2 = fit0(
JPoint4D(*start_pos), dataL0_selected.begin(), dataL0_selected.end());
239 JPoint4D pt(fitPos0, fit0.value.getT());
248 chi2 = fit(pt, dataL0_selected.begin(), dataL0_selected.end());
250 int NDF = dataL0_selected.size() - fit.step.size();
256 out.rbegin()->setW(13, chi2);
257 out.rbegin()->setW(14, NDF);
260 DEBUG(
"Too few hits " << endl);