48 using namespace JSUPPORT;
49 using namespace KM3NETDAQ;
51 using namespace JAANET;
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);
96 using namespace JTRIGGER;
97 using namespace JDETECTOR;
98 using namespace JGEOMETRY3D;
99 using namespace JTOOLS;
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);
Utility class to parse command line options.
Data structure for direction in three dimensions.
Template specialisation of L0 builder for JHitL0 data type.
Data structure for vertex fit.
Router for direct addressing of module data in detector data structure.
General purpose class for parallel reading of objects from a single file or multiple files...
Container for historical events.
double getDistance(const JVector3D &pos) const
Get distance to point.
double getDot(const JAngle3D &angle) const
Get dot product.
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 defining the range of iterations of objects.
Regressor function object for JPoint4D fit using JSimplex minimiser.
Data structure for vector in three dimensions.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
void load(const JString &file_name, JDetector &detector)
Load detector from input file.
counter_type advance(counter_type &counter, const counter_type value, const counter_type limit=std::numeric_limits< counter_type >::max())
Advance counter.
double getQuality(const double chi2, const int NDF)
Get quality of fit.
Data structure for set of track fit results.
General purpose class for object reading from a list of file names.
Data structure for L0 hit.
bool qualitySorter(const JFIT::JFit &first, const JFIT::JFit &second)
Comparison of fit results.
void copy(const Head &from, JHead &to)
Copy header from from to to.
Data structure for position in three dimensions.
const JLimit & getLimit() const
Get limit.
Data structure for normalised vector in positive z-direction.
Tukey's biweight M-estimator.
#define DEBUG(A)
Message macros.
JPosition3D getPosition(const Vec &v)
Get position.