Jpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Functions
JShowerPositionFit.cc File Reference
#include <iostream>
#include <string>
#include "Jeep/JParser.hh"
#include "evt/Head.hh"
#include "evt/Det.hh"
#include "JDAQ/JDAQTimeslice.hh"
#include "JDAQ/JDAQEvent.hh"
#include "JDAQ/JDAQSummaryslice.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JModuleRouter.hh"
#include "JFit/JEvt.hh"
#include "JFit/JEvtToolkit.hh"
#include "JFit/JFitParameters.hh"
#include "JFit/JHistory.hh"
#include "JFit/JPoint4DRegressor.hh"
#include "JLang/JTypeList.hh"
#include "JSupport/JParallelFileScanner.hh"
#include "JSupport/JTreeScanner.hh"
#include "JSupport/JFileRecorder.hh"
#include "JSupport/JMeta.hh"
#include "JSupport/JSupport.hh"
#include "JTrigger/JHit.hh"
#include "JTrigger/JAlgorithm.hh"
#include "JTrigger/JBind2nd.hh"
#include "JTrigger/JTriggerParameters.hh"
#include "JTrigger/JHitL0.hh"
#include "JTrigger/JBuildL0.hh"
#include "JTrigger/JBuildL2.hh"
#include "JTrigger/JMatch3G.hh"
#include "JGeometry3D/JTrack3D.hh"
#include <TFile.h>
#include <TH2D.h>

Go to the source code of this file.

Functions

int main (int argc, char **argv)
 

Function Documentation

int main ( int  argc,
char **  argv 
)

Definition at line 45 of file JShowerPositionFit.cc.

45  {
46 
47  using namespace std;
48  using namespace JSUPPORT;
49  using namespace KM3NETDAQ;
50  using namespace JFIT;
51  using namespace JAANET;
52 
53  typedef JParallelFileScanner< JTypeList<JDAQEvent, JEvt> > JParallelFileScanner_t;
54  typedef JParallelFileScanner_t::multi_pointer_type multi_pointer_type;
56 
57  JParallelFileScanner_t inputFile;
59  JLimit_t& numberOfEvents = inputFile.getLimit();
60  string detectorFile;
61  string pdfFile;
62  double Tmax_ns;
63  double ctMin;
64  double Dmax_m;
65  JRange<double> pos_grid;
66  JRange<double> time_grid;
67  double pos_step;
68  double time_step;
69  size_t numberOfPrefits;
70  int debug;
71 
72  try{
73  JParser<> zap;
74 
75  zap['f'] = make_field(inputFile);
76  zap['o'] = make_field(outputFile) = "PositionFit.root";
77  zap['a'] = make_field(detectorFile);
78  zap['F'] = make_field(pdfFile) = "PDF_ShowerPositionFit.root";
79  zap['n'] = make_field(numberOfEvents) = JLimit::max();
80  zap['N'] = make_field(Tmax_ns) = 20.0; // [ns]
81  zap['C'] = make_field(ctMin) = -0.7; // optimal value
82  zap['D'] = make_field(Dmax_m) = 80.0; // [m]
83  zap['P'] = make_field(pos_grid) = JRange<double>(-28, 28); // [m] position grid
84  zap['p'] = make_field(pos_step) = 6; // [m] density of the grid
85  zap['T'] = make_field(time_grid) = JRange<double>(-60, 60); // [ns] time grid
86  zap['t'] = make_field(time_step) = 15; // [ns] density of the grid
87  zap['R'] = make_field(numberOfPrefits) = 1;
88  zap['d'] = make_field(debug) = 2;
89 
90  if (zap.read(argc, argv) != 0) return 1;
91  }
92  catch(const exception& error){
93  FATAL(error.what() << endl);
94  }
95 
96  using namespace JTRIGGER;
97  using namespace JDETECTOR;
98  using namespace JGEOMETRY3D;
99  using namespace JTOOLS;
100 
102 
103  try {
104  load(detectorFile, detector);
105  }
106  catch(const JException& error) {
107  FATAL(error);
108  }
109 
110  TFile* pdf_file = new TFile(pdfFile.c_str());
111  TH2D* hpdf = (TH2D*)pdf_file->Get("hPDF2Dist");
112 
113  const JModuleRouter moduleRouter(detector);
114 
115  typedef vector<JHitL0> JDataL0_t;
116  typedef vector<JHitL1> JDataL1_t;
117 
118  const JBuildL0<JHitL0> buildL0;
119  const JBuildL2<JHitL1> buildL2(3, Tmax_ns, ctMin);
120  const JMatch3G<JHitL1> match3G(Dmax_m, Tmax_ns); // CAUSALITY FOR SHOWERS
121 
122  outputFile.open();
123  outputFile.put(JMeta(argc, argv));
124 
125  typedef JRegressor<JPoint4D, JSimplex> JRegressor_t;
127  JRegressor_t::MAXIMUM_ITERATIONS = 10000;
128 
129  JRegressor_t fit0(true);
130  fit0.estimator.reset(new JMEstimatorTukey(10.0));
131  JRegressor_t fit(hpdf, false);
132 
133  for (JTreeScanner<Evt> in(inputFile); inputFile.hasNext(); ) {
134 
135  STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);
136 
137  multi_pointer_type ps = inputFile.next();
138 
139  JDAQEvent* tev = ps;
140  JEvt* evt = ps;
141  JEvt out;
142 
143  JEvt::iterator __end = evt->end();
144  if (numberOfPrefits > 0) {
145  advance(__end = evt->begin(), min(numberOfPrefits, evt->size()));
146  }
147  partial_sort(evt->begin(), __end, evt->end(), qualitySorter);
148 
149  if(!evt->empty()) {
150 
151  copy(evt->begin(), __end, back_inserter(out));
152 
153  JDataL0_t dataL0;
154  JDataL0_t dataL0_selected;
155  JDataL1_t data;
156 
157  buildL0(*tev, moduleRouter, true, back_inserter(dataL0));
158  buildL2(JDAQTimeslice(*tev, false), moduleRouter, back_inserter(data));
159 
160  for (JEvt::const_iterator shower = evt->begin(); shower != __end; ++shower) {
161 
162  const JPoint4D vertex(getPosition(*shower), shower->getT());
163 
164  for (JDataL0_t::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {
165 
166  JHitL0 hit(*i);
167  const JPosition3D hit_pos(hit.getPosition());
168  const double D = hit_pos.getDistance(vertex);
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());
172 
173  // hit selection
174  if(data.size()>0){
175  if(D < Dmax_m && (t_res > -40 && t_res < 40) && (cosT >= -1 && cosT <= 0.1)){
176  dataL0_selected.push_back(hit);
177  }
178  } else {
179  if(D < Dmax_m && (t_res > -50 && t_res < 50) && (cosT >= -1 && cosT <= 0.1)){
180  dataL0_selected.push_back(hit);
181  }
182  }
183 
184  }
185 
186  if(dataL0_selected.size() != 0){
187 
188  double chi2 = 0.0;
189 
190  vector<double> LogChi2;
191  vector<JPoint4D> start_positions;
192 
193  for(double x = pos_grid.getLowerLimit(); x <= pos_grid.getUpperLimit();
194  x += pos_step){
195  for(double y = pos_grid.getLowerLimit(); y <= pos_grid.getUpperLimit();
196  y += pos_step){
197  for(double z = pos_grid.getLowerLimit(); z <= pos_grid.getUpperLimit();
198  z += pos_step){
199  for(double t = time_grid.getLowerLimit(); t <= time_grid.getUpperLimit();
200  t += time_step){
201 
202  chi2 = 0.0;
203 
204  JPosition3D pos_shifted(vertex.getX() + x, vertex.getY() + y, vertex.getZ() + z);
205  JPoint4D point_shifted(pos_shifted, vertex.getT() + t);
206 
207  for(JDataL0_t::const_iterator i = dataL0_selected.begin(); i != dataL0_selected.end(); ++i){
208  JHitL0 hit(*i);
209  chi2 += fit0(point_shifted, hit);
210  }
211 
212  LogChi2.push_back(chi2);
213  start_positions.push_back(point_shifted);
214 
215  }
216  }
217  }
218  }
219 
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);
226  }
227 
228  for(vector<JPoint4D>::iterator start_pos = start_positions.begin(); start_pos != start_positions.end(); ++start_pos){
229 
230  double simplex_step = 1; // [m]
231  fit0.step.resize(4);
232  fit0.step[0] = JPoint4D(JVector3D(simplex_step, 0.0, 0.0), 0.0);
233  fit0.step[1] = JPoint4D(JVector3D(0.0, simplex_step, 0.0), 0.0);
234  fit0.step[2] = JPoint4D(JVector3D(0.0, 0.0, simplex_step), 0.0);
235  fit0.step[3] = JPoint4D(JVector3D(0.0, 0.0, 0.0), 20);
236 
237  chi2 = fit0(JPoint4D(*start_pos), dataL0_selected.begin(), dataL0_selected.end());
238 
239  JVector3D fitPos0(fit0.value);
240  JPoint4D pt(fitPos0, fit0.value.getT());
241 
242  simplex_step = 0.3; // [m]
243  fit.step.resize(4);
244  fit.step[0] = JPoint4D(JVector3D(simplex_step, 0.0, 0.0), 0.0);
245  fit.step[1] = JPoint4D(JVector3D(0.0, simplex_step, 0.0), 0.0);
246  fit.step[2] = JPoint4D(JVector3D(0.0, 0.0, simplex_step), 0.0);
247  fit.step[3] = JPoint4D(JVector3D(0.0, 0.0, 0.0), 10);
248 
249  chi2 = fit(pt, dataL0_selected.begin(), dataL0_selected.end());
250 
251  int NDF = dataL0_selected.size() - fit.step.size();
252 
253  JVector3D fitPos(fit.value);
254  JTrack3D fitResult(fitPos, JVersor3Z(0,0), fit.value.getT());
255 
256  out.push_back(getFit(JHistory(shower->getHistory()).add(JSHOWERPOSITIONFIT), fitResult, getQuality(chi2, NDF), NDF));
257  out.rbegin()->setW(13, chi2);
258  out.rbegin()->setW(14, NDF);
259  }
260  } else{
261  DEBUG("Too few hits " << endl);
262  }
263  }
264  sort(out.begin(), out.end(), qualitySorter);
265  }
266  outputFile.put(out);
267  }
268  STATUS(endl);
269 
271  io >> outputFile;
272  outputFile.close();
273 
274 }
Auxiliary class for ROOT I/O of application specific meta data.
Definition: JMeta.hh:71
Object writing to file.
Utility class to parse command line options.
Definition: JParser.hh:1410
3G match criterion.
Definition: JMatch3G.hh:29
Data structure for direction in three dimensions.
Definition: JDirection3D.hh:32
Template specialisation of L0 builder for JHitL0 data type.
Definition: JBuildL0.hh:102
Data structure for vertex fit.
Definition: JPoint4D.hh:22
#define STATUS(A)
Definition: JMessage.hh:61
JShowerPositionFit.cc.
Detector data structure.
Definition: JDetector.hh:77
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.
Definition: JHistory.hh:95
string outputFile
double getDistance(const JVector3D &pos) const
Get distance to point.
Definition: JVector3D.hh:266
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.
Definition: JEvtToolkit.hh:116
Type list.
Definition: JTypeList.hh:22
Auxiliary class for defining the range of iterations of objects.
Definition: JLimit.hh:41
Template L2 builder.
Definition: JBuildL2.hh:47
Detector file.
Definition: JHead.hh:126
Regressor function object for JPoint4D fit using JSimplex minimiser.
Data structure for vector in three dimensions.
Definition: JVector3D.hh:32
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1836
Data time slice.
void load(const JString &file_name, JDetector &detector)
Load detector from input file.
int debug
debug level
Definition: JSirene.cc:59
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
#define FATAL(A)
Definition: JMessage.hh:65
double getQuality(const double chi2, const int NDF)
Get quality of fit.
Definition: JEvtToolkit.hh:195
Data structure for set of track fit results.
Definition: JEvt.hh:312
General purpose class for object reading from a list of file names.
Data structure for L0 hit.
Definition: JHitL0.hh:27
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.
Definition: JHead.cc:40
Data structure for position in three dimensions.
Definition: JPosition3D.hh:35
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:72
Data structure for normalised vector in positive z-direction.
Definition: JVersor3Z.hh:36
Tukey&#39;s biweight M-estimator.
Definition: JMEstimator.hh:105
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:60
JPosition3D getPosition(const Vec &v)
Get position.