45 int main(
int argc,
char **argv){
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());
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);
double getT() const
Get calibrated time of hit.
Utility class to parse command line options.
Data structure for direction in three dimensions.
Algorithms for hit clustering and sorting.
Template specialisation of L0 builder for JHitL0 data type.
Data structure for vertex fit.
const JDirection3D & getDirection() const
Get direction.
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.
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.
Data structure for detector geometry and calibration.
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.
Basic data structure for L0 hit.
Definition of fit parameters from various applications.
Basic data structure for time and time over threshold information of hit.
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
Parallel scanning of objects from a single file or multiple files according a format that follows fro...
void load(const JString &file_name, JDetector &detector)
Load detector from input file.
const JPosition3D & getPosition() const
Get position.
counter_type advance(counter_type &counter, const counter_type value, const counter_type limit=std::numeric_limits< counter_type >::max())
Advance counter.
Direct access to module in detector data structure.
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.
Utility class to parse command line options.
Data structure for L0 hit.
ROOT TTree parameter settings.
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.
Match operator for Cherenkov light from shower in any direction.
#define DEBUG(A)
Message macros.
JPosition3D getPosition(const Vec &v)
Get position.
int main(int argc, char *argv[])