Jpp
JShowerPositionFit.cc
Go to the documentation of this file.
1 #include <iostream>
2 #include <string>
3 
4 #include "Jeep/JParser.hh"
5 
9 #include "JDAQ/JDAQEventIO.hh"
11 
12 #include "JDetector/JDetector.hh"
15 
16 #include "JFit/JEvt.hh"
17 #include "JFit/JEvtToolkit.hh"
18 #include "JFit/JFitParameters.hh"
19 #include "JFit/JHistory.hh"
21 
22 #include "JLang/JTypeList.hh"
23 
25 #include "JSupport/JTreeScanner.hh"
27 #include "JSupport/JMeta.hh"
28 #include "JSupport/JSupport.hh"
29 
30 #include "JTrigger/JHit.hh"
31 #include "JTrigger/JAlgorithm.hh"
32 #include "JTrigger/JBind2nd.hh"
34 #include "JTrigger/JHitL0.hh"
35 #include "JTrigger/JBuildL0.hh"
36 #include "JTrigger/JBuildL2.hh"
37 #include "JTrigger/JMatch3G.hh"
38 
39 #include "JGeometry3D/JTrack3D.hh"
40 
41 #include <TFile.h>
42 #include <TH2D.h>
43 
44 int main(int argc, char **argv){
45 
46  using namespace std;
47  using namespace JSUPPORT;
48  using namespace KM3NETDAQ;
49  using namespace JFIT;
50  using namespace JAANET;
51 
52  typedef JParallelFileScanner< JTypeList<JDAQEvent, JEvt> > JParallelFileScanner_t;
53  typedef JParallelFileScanner_t::multi_pointer_type multi_pointer_type;
55 
56  JParallelFileScanner_t inputFile;
58  JLimit_t& numberOfEvents = inputFile.getLimit();
59  string detectorFile;
60  string pdfFile;
61  double Tmax_ns;
62  double ctMin;
63  double Dmax_m;
64  JRange<double> pos_grid;
65  JRange<double> time_grid;
66  double pos_step;
67  double time_step;
68  size_t numberOfPrefits;
69  int debug;
70 
71  try{
72  JParser<> zap;
73 
74  zap['f'] = make_field(inputFile);
75  zap['o'] = make_field(outputFile) = "PositionFit.root";
76  zap['a'] = make_field(detectorFile);
77  zap['F'] = make_field(pdfFile) = "PDF_ShowerPositionFit.root";
78  zap['n'] = make_field(numberOfEvents) = JLimit::max();
79  zap['N'] = make_field(Tmax_ns) = 20.0; // [ns]
80  zap['C'] = make_field(ctMin) = -0.7; // optimal value
81  zap['D'] = make_field(Dmax_m) = 80.0; // [m]
82  zap['P'] = make_field(pos_grid) = JRange<double>(-28, 28); // [m] position grid
83  zap['p'] = make_field(pos_step) = 6; // [m] density of the grid
84  zap['T'] = make_field(time_grid) = JRange<double>(-60, 60); // [ns] time grid
85  zap['t'] = make_field(time_step) = 15; // [ns] density of the grid
86  zap['R'] = make_field(numberOfPrefits) = 1;
87  zap['d'] = make_field(debug) = 2;
88 
89  if (zap.read(argc, argv) != 0) return 1;
90  }
91  catch(const exception& error){
92  FATAL(error.what() << endl);
93  }
94 
95  using namespace JTRIGGER;
96  using namespace JDETECTOR;
97  using namespace JGEOMETRY3D;
98  using namespace JTOOLS;
99 
101 
102  try {
103  load(detectorFile, detector);
104  }
105  catch(const JException& error) {
106  FATAL(error);
107  }
108 
109  TFile* pdf_file = new TFile(pdfFile.c_str());
110  TH2D* hpdf = (TH2D*)pdf_file->Get("hPDF2Dist");
111 
112  const JModuleRouter moduleRouter(detector);
113 
114  typedef vector<JHitL0> JDataL0_t;
115  typedef vector<JHitL1> JDataL1_t;
116 
117  const JBuildL0<JHitL0> buildL0;
118  const JBuildL2<JHitL1> buildL2(JL2Parameters(3, Tmax_ns, ctMin));
119  const JMatch3G<JHitL1> match3G(Dmax_m, Tmax_ns); // CAUSALITY FOR SHOWERS
120 
121  outputFile.open();
122  outputFile.put(JMeta(argc, argv));
123 
124  typedef JRegressor<JPoint4D, JSimplex> JRegressor_t;
126  JRegressor_t::MAXIMUM_ITERATIONS = 10000;
127 
128  JRegressor_t fit0(true);
129  fit0.estimator.reset(new JMEstimatorTukey(10.0));
130  JRegressor_t fit(hpdf, false);
131 
132  for (JTreeScanner<Evt> in(inputFile); inputFile.hasNext(); ) {
133 
134  STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);
135 
136  multi_pointer_type ps = inputFile.next();
137 
138  JDAQEvent* tev = ps;
139  JEvt* evt = ps;
140  JEvt out;
141 
142  JEvt::iterator __end = evt->end();
143  if (numberOfPrefits > 0) {
144  advance(__end = evt->begin(), min(numberOfPrefits, evt->size()));
145  }
146  partial_sort(evt->begin(), __end, evt->end(), qualitySorter);
147 
148  if(!evt->empty()) {
149 
150  copy(evt->begin(), __end, back_inserter(out));
151 
152  JDataL0_t dataL0;
153  JDataL0_t dataL0_selected;
154  JDataL1_t data;
155 
156  buildL0(*tev, moduleRouter, true, back_inserter(dataL0));
157  buildL2(JDAQTimeslice(*tev, false), moduleRouter, back_inserter(data));
158 
159  for (JEvt::const_iterator shower = evt->begin(); shower != __end; ++shower) {
160 
161  const JPoint4D vertex(getPosition(*shower), shower->getT());
162 
163  for (JDataL0_t::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {
164 
165  JHitL0 hit(*i);
166  const JPosition3D hit_pos(hit.getPosition());
167  const double D = hit_pos.getDistance(vertex);
168  const double t_res = hit.getT() - vertex.getT(hit_pos);
169  const JDirection3D photonDir(hit_pos - vertex.getPosition());
170  const double cosT = photonDir.getDot(hit.getDirection());
171 
172  // hit selection
173  if(data.size()>0){
174  if(D < Dmax_m && (t_res > -40 && t_res < 40) && (cosT >= -1 && cosT <= 0.1)){
175  dataL0_selected.push_back(hit);
176  }
177  } else {
178  if(D < Dmax_m && (t_res > -50 && t_res < 50) && (cosT >= -1 && cosT <= 0.1)){
179  dataL0_selected.push_back(hit);
180  }
181  }
182 
183  }
184 
185  if(dataL0_selected.size() != 0){
186 
187  double chi2 = 0.0;
188 
189  vector<double> LogChi2;
190  vector<JPoint4D> start_positions;
191 
192  for(double x = pos_grid.getLowerLimit(); x <= pos_grid.getUpperLimit();
193  x += pos_step){
194  for(double y = pos_grid.getLowerLimit(); y <= pos_grid.getUpperLimit();
195  y += pos_step){
196  for(double z = pos_grid.getLowerLimit(); z <= pos_grid.getUpperLimit();
197  z += pos_step){
198  for(double t = time_grid.getLowerLimit(); t <= time_grid.getUpperLimit();
199  t += time_step){
200 
201  chi2 = 0.0;
202 
203  JPosition3D pos_shifted(vertex.getX() + x, vertex.getY() + y, vertex.getZ() + z);
204  JPoint4D point_shifted(pos_shifted, vertex.getT() + t);
205 
206  for(JDataL0_t::const_iterator i = dataL0_selected.begin(); i != dataL0_selected.end(); ++i){
207  JHitL0 hit(*i);
208  chi2 += fit0(point_shifted, hit);
209  }
210 
211  LogChi2.push_back(chi2);
212  start_positions.push_back(point_shifted);
213 
214  }
215  }
216  }
217  }
218 
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);
225  }
226 
227  for(vector<JPoint4D>::iterator start_pos = start_positions.begin(); start_pos != start_positions.end(); ++start_pos){
228 
229  double simplex_step = 1; // [m]
230  fit0.step.resize(4);
231  fit0.step[0] = JPoint4D(JVector3D(simplex_step, 0.0, 0.0), 0.0);
232  fit0.step[1] = JPoint4D(JVector3D(0.0, simplex_step, 0.0), 0.0);
233  fit0.step[2] = JPoint4D(JVector3D(0.0, 0.0, simplex_step), 0.0);
234  fit0.step[3] = JPoint4D(JVector3D(0.0, 0.0, 0.0), 20);
235 
236  chi2 = fit0(JPoint4D(*start_pos), dataL0_selected.begin(), dataL0_selected.end());
237 
238  JVector3D fitPos0(fit0.value);
239  JPoint4D pt(fitPos0, fit0.value.getT());
240 
241  simplex_step = 0.3; // [m]
242  fit.step.resize(4);
243  fit.step[0] = JPoint4D(JVector3D(simplex_step, 0.0, 0.0), 0.0);
244  fit.step[1] = JPoint4D(JVector3D(0.0, simplex_step, 0.0), 0.0);
245  fit.step[2] = JPoint4D(JVector3D(0.0, 0.0, simplex_step), 0.0);
246  fit.step[3] = JPoint4D(JVector3D(0.0, 0.0, 0.0), 10);
247 
248  chi2 = fit(pt, dataL0_selected.begin(), dataL0_selected.end());
249 
250  int NDF = dataL0_selected.size() - fit.step.size();
251 
252  JVector3D fitPos(fit.value);
253  JTrack3D fitResult(fitPos, JVersor3Z(0,0), fit.value.getT());
254 
255  out.push_back(getFit(JHistory(shower->getHistory()).add(JSHOWERPOSITIONFIT), fitResult, getQuality(chi2, NDF), NDF));
256  out.rbegin()->setW(13, chi2);
257  out.rbegin()->setW(14, NDF);
258  }
259  } else{
260  DEBUG("Too few hits " << endl);
261  }
262  }
263  sort(out.begin(), out.end(), qualitySorter);
264  }
265  outputFile.put(out);
266  }
267  STATUS(endl);
268 
270  io >> outputFile;
271  outputFile.close();
272 
273 }
JSHOWERPOSITIONFIT
static const int JSHOWERPOSITIONFIT
Definition: reconstruction.hh:24
JMeta.hh
Head.hh
JTOOLS::JRange::getLowerLimit
T getLowerLimit() const
Get lower limit.
Definition: JRange.hh:215
KM3NETDAQ::JDAQEvent
DAQ Event.
Definition: JDAQEvent.hh:30
JHistory.hh
JTriggerParameters.hh
JTOOLS::JRange::getUpperLimit
T getUpperLimit() const
Get upper limit.
Definition: JRange.hh:226
JFileRecorder.hh
JSUPPORT::JLimit
Auxiliary class for defining the range of iterations of objects.
Definition: JLimit.hh:41
JFIT
Auxiliary classes and methods for linear and iterative data regression.
Definition: JEnergy.hh:15
JTRIGGER::JL2Parameters
Data structure for L2 parameters.
Definition: JTriggerParameters.hh:33
JTrack3D.hh
JROOT::advance
counter_type advance(counter_type &counter, const counter_type value, const counter_type limit=std::numeric_limits< counter_type >::max())
Advance counter.
Definition: JCounter.hh:35
JGEOMETRY3D::JVersor3Z
Data structure for normalised vector in positive z-direction.
Definition: JVersor3Z.hh:37
qualitySorter
bool qualitySorter(const JFIT::JFit &first, const JFIT::JFit &second)
Comparison of fit results.
Definition: quality_sorter.cc:10
JGEOMETRY3D::JDirection3D::getDirection
const JDirection3D & getDirection() const
Get direction.
Definition: JDirection3D.hh:106
JDETECTOR::load
void load(const JString &file_name, JDetector &detector)
Load detector from input file.
Definition: JDetectorToolkit.hh:476
JGEOMETRY3D::JVector3D::getDistance
double getDistance(const JVector3D &pos) const
Get distance to point.
Definition: JVector3D.hh:269
JGEOMETRY3D::JDirection3D::getDot
double getDot(const JAngle3D &angle) const
Get dot product.
Definition: JDirection3D.hh:333
std::vector< JHitL0 >
JTOOLS::JRange< double >
JEvtToolkit.hh
JGEOMETRY3D::JDirection3D
Data structure for direction in three dimensions.
Definition: JDirection3D.hh:32
Evt.hh
KM3NETDAQ::JDAQTimeslice
Data time slice.
Definition: JDAQTimeslice.hh:30
JPARSER::JParser
Utility class to parse command line options.
Definition: JParser.hh:1493
JPoint4DRegressor.hh
JFIT::JEvt
Data structure for set of track fit results.
Definition: JEvt.hh:293
JAANET::copy
void copy(const Head &from, JHead &to)
Copy header from from to to.
Definition: JHead.cc:152
JTreeScanner.hh
JTRIGGER::JHit::getT
double getT() const
Get calibrated time of hit.
Definition: JHit.hh:143
JSupport.hh
JEvt.hh
JBuildL0.hh
JDAQEventIO.hh
JGEOMETRY3D::JVector3D
Data structure for vector in three dimensions.
Definition: JVector3D.hh:33
JMatch3G.hh
JDAQTimesliceIO.hh
debug
int debug
debug level
Definition: JSirene.cc:59
JSUPPORT::JTreeScanner< Evt >
JGEOMETRY3D::JPosition3D
Data structure for position in three dimensions.
Definition: JPosition3D.hh:35
JSUPPORT::JLimit::getLimit
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:73
JHit.hh
JModuleRouter.hh
JParallelFileScanner.hh
JBind2nd.hh
JLANG::JTypeList
Type list.
Definition: JTypeList.hh:22
JTRIGGER::JMatch3G
3G match criterion.
Definition: JMatch3G.hh:29
STATUS
#define STATUS(A)
Definition: JMessage.hh:63
JBuildL2.hh
JFIT::JPoint4D
Data structure for vertex fit.
Definition: JPoint4D.hh:22
JFIT::JRegressor< JPoint4D, JSimplex >
Regressor function object for JPoint4D fit using JSimplex minimiser.
Definition: JPoint4DRegressor.hh:47
JParser.hh
JDetectorToolkit.hh
JAANET
Extensions to AAnet data format.
Definition: JAAnetToolkit.hh:36
JDETECTOR::JModuleRouter
Router for direct addressing of module data in detector data structure.
Definition: JModuleRouter.hh:34
JFIT::JMEstimatorTukey
Tukey's biweight M-estimator.
Definition: JMEstimator.hh:105
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1954
JFitParameters.hh
JSUPPORT::JParallelFileScanner
General purpose class for parallel reading of objects from a single file or multiple files.
Definition: JParallelFileScanner.hh:31
JSUPPORT::JMultipleFileScanner
General purpose class for object reading from a list of file names.
Definition: JMultipleFileScanner.hh:167
JDETECTOR::JDetector
Detector data structure.
Definition: JDetector.hh:80
JGEOMETRY3D::JPosition3D::getPosition
const JPosition3D & getPosition() const
Get position.
Definition: JPosition3D.hh:129
JAANET::getPosition
JPosition3D getPosition(const Vec &v)
Get position.
Definition: JAAnetToolkit.hh:197
JTRIGGER::JHitL0
Data structure for L0 hit.
Definition: JHitL0.hh:27
JAANET::detector
Detector file.
Definition: JHead.hh:130
JSUPPORT::JMeta
Auxiliary class for ROOT I/O of application specific meta data.
Definition: JMeta.hh:71
DEBUG
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
std
Definition: jaanetDictionary.h:36
JSUPPORT
Support classes and methods for experiment specific I/O.
Definition: JDataWriter.cc:38
KM3NETDAQ
KM3NeT DAQ data structures and auxiliaries.
Definition: DataQueue.cc:39
JTypeList.hh
JGEOMETRY3D::JTrack3D
3D track.
Definition: JTrack3D.hh:30
JDAQSummarysliceIO.hh
JGEOMETRY3D
Auxiliary classes and methods for 3D geometrical objects and operations.
Definition: JAngle3D.hh:18
JFIT::getFit
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.
Definition: JEvtToolkit.hh:116
main
int main(int argc, char **argv)
Definition: JShowerPositionFit.cc:44
JTRIGGER::JBuildL0< JHitL0 >
Template specialisation of L0 builder for JHitL0 data type.
Definition: JBuildL0.hh:101
JFIT::JHistory
Container for historical events.
Definition: JHistory.hh:95
JDetector.hh
JTRIGGER
Checksum.
Definition: JSupport/JSupport.hh:35
JTOOLS
Auxiliary classes and methods for multi-dimensional interpolations and histograms.
Definition: JAbstractCollection.hh:9
JHitL0.hh
FATAL
#define FATAL(A)
Definition: JMessage.hh:67
outputFile
string outputFile
Definition: JDAQTimesliceSelector.cc:37
JSUPPORT::JFileRecorder
Object writing to file.
Definition: JFileRecorder.hh:41
JFIT::getQuality
double getQuality(const double chi2, const int NDF)
Get quality of fit.
Definition: JEvtToolkit.hh:195
JDETECTOR
Auxiliary classes and methods for detector calibration.
Definition: JAnchor.hh:12
JLANG::JException
General exception.
Definition: JException.hh:23
JPARSER::JParser::read
int read(const int argc, const char *const argv[])
Parse the program's command line options.
Definition: JParser.hh:1791
JTRIGGER::JBuildL2
Template L2 builder.
Definition: JBuildL2.hh:45
JAlgorithm.hh