Jpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JMatrixNZ.cc
Go to the documentation of this file.
1 #include <string>
2 #include <iostream>
3 #include <iomanip>
4 #include <vector>
5 #include <map>
6 #include <algorithm>
7 
8 #include "TROOT.h"
9 #include "TFile.h"
10 #include "TH1D.h"
11 #include "TProfile.h"
12 #include "TRandom.h"
13 
17 #include "JDAQ/JDAQEventIO.hh"
18 #include "JDAQ/JDAQTimesliceIO.hh"
19 #include "JAAnet/JHead.hh"
20 #include "JAAnet/JHeadToolkit.hh"
21 #include "JAAnet/JAAnetToolkit.hh"
22 #include "JPhysics/JConstants.hh"
26 #include "JDetector/JDetector.hh"
29 #include "JDetector/JPMTRouter.hh"
30 #include "JTrigger/JHit.hh"
31 #include "JTrigger/JFrame.hh"
32 #include "JTrigger/JTimeslice.hh"
33 #include "JTrigger/JHitL0.hh"
34 #include "JTrigger/JHitL1.hh"
35 #include "JTrigger/JBuildL1.hh"
36 #include "JTrigger/JBuildL2.hh"
37 #include "JTrigger/JAlgorithm.hh"
38 #include "JTrigger/JMatch1D.hh"
39 #include "JTrigger/JMatch3B.hh"
43 #include "JSupport/JSupport.hh"
44 
45 #include "JFit/JMatrixNZ.hh"
46 #include "JFit/JFitToolkit.hh"
47 
48 #include "Jeep/JParser.hh"
49 #include "Jeep/JMessage.hh"
50 
51 
52 /**
53  * \file
54  *
55  * Example program to test chi2 calculations based on JFIT::JMatrixNZ.
56  * \author mdejong
57  */
58 int main(int argc, char **argv)
59 {
60  using namespace std;
61  using namespace JPP;
62  using namespace KM3NETDAQ;
63 
64  JTriggeredFileScanner<> inputFile;
65  JLimit_t& numberOfEvents = inputFile.getLimit();
66  string outputFile;
67  string detectorFile;
68  double Tmax_ns;
69  double roadWidth_m;
70  double ctMin;
71  double alpha_deg;
72  double sigma_ns;
73  int numberOfOutliers;
74  bool useTrue;
75  int debug;
76 
77  try {
78 
79  JParser<> zap("Example program to test chi2 calculations of co-variance matrix including direction uncertainty.");
80 
81  zap['f'] = make_field(inputFile);
82  zap['o'] = make_field(outputFile) = "matrixnz.root";
83  zap['a'] = make_field(detectorFile);
84  zap['n'] = make_field(numberOfEvents) = JLimit::max();
85  zap['T'] = make_field(Tmax_ns) = 20.0; // [ns]
86  zap['R'] = make_field(roadWidth_m) = 150.0; // [m]
87  zap['C'] = make_field(ctMin) = 0.0; //
88  zap['A'] = make_field(alpha_deg);
89  zap['S'] = make_field(sigma_ns);
90  zap['O'] = make_field(numberOfOutliers);
91  zap['U'] = make_field(useTrue);
92  zap['d'] = make_field(debug) = 1;
93 
94  zap(argc, argv);
95  }
96  catch(const exception& error) {
97  FATAL(error.what() << endl);
98  }
99 
100 
101  using namespace KM3NETDAQ;
102 
103  cout.tie(&cerr);
104 
105 
106  const unsigned int NUMBER_OF_HITS = 6;
107  const double STANDARD_DEVIATIONS = 3.0; // [unit]
108  const double HIT_OFF = 1.0e3 * sigma_ns*sigma_ns; // [ns^2]
109 
110 
112 
113  try {
114  load(detectorFile, detector);
115  }
116  catch(const JException& error) {
117  FATAL(error);
118  }
119 
120  const JModuleRouter moduleRouter(detector);
121  const JPMTRouter pmtRouter (detector);
122 
123  const JPosition3D center = get<JPosition3D>(getHeader(inputFile));
124 
125 
126  TFile out(outputFile.c_str(), "recreate");
127 
128 
129  TH1D h0("h0", NULL, 40, 0.0, 20.0);
130  TH1D h1("h1", NULL, 40, 0.0, 20.0);
131 
132  TH1D p0("p0", NULL, 20, 0.0, 1.0);
133  TH1D p1("p1", NULL, 20, 0.0, 1.0);
134 
135  vector<double> X;
136 
137  {
138  double x = -0.5;
139 
140  for ( ; x < 10.0; x += 1.0) { X.push_back(x); }
141  for ( ; x < 25.0; x += 2.0) { X.push_back(x); }
142  for ( ; x < 100.0; x += 5.0) { X.push_back(x); }
143  }
144 
145  TProfile n0("n0", NULL, X.size() - 1, X.data());
146  TProfile n1("n1", NULL, X.size() - 1, X.data());
147 
148 
149  typedef vector<JHitL1> JData_t;
150  const JBuildL2<JHitL1> buildL2(JL2Parameters(2, Tmax_ns, ctMin));
151  const JMatch3B<JHitL1> match3B(roadWidth_m, Tmax_ns);
152  const JMatch1D<JHitL1> match1D(roadWidth_m, Tmax_ns);
153 
154 
155  while (inputFile.hasNext()) {
156 
157  STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);
158 
159  JTriggeredFileScanner<>::multi_pointer_type ps = inputFile.next();
160 
161  const JDAQEvent* tev = ps;
162  const Evt* event = ps;
163 
164  const time_converter converter(*event, *tev);
165 
166  vector<Trk>::const_iterator muon = find_if(event->mc_trks.begin(), event->mc_trks.end(), is_muon);
167 
168  if (muon != event->mc_trks.end()) {
169 
170  const JRotation3D R(getDirection(*muon));
171 
172  JLine1Z tz(getPosition(*muon), converter.putTime(muon->t));
173 
174  tz.add(center);
175  tz.rotate(R);
176 
177  const double theta = alpha_deg * PI / 180.0;
178  const double phi = gRandom->Uniform(0.0, 2*PI);
179 
180  const JRotation3D Rs(JAngle3D(theta,phi));
181 
182 
183  JData_t data;
184 
185  if (!useTrue) {
186 
187  buildL2(*tev, moduleRouter, back_inserter(data));
188 
189  data.erase(unique(data.begin(), data.end(), equal_to<JDAQModuleIdentifier>()), data.end());
190 
191  } else {
192 
194 
195  for (vector<Hit>::const_iterator i = event->mc_hits.begin(); i != event->mc_hits.end(); ++i) {
196 
197  if (!is_noise(*i)) {
198 
199  const JPMTAddress& address = pmtRouter.getAddress(i->pmt_id);
200  const JPMTIdentifier& id = pmtRouter.getIdentifier(address);
201  const JPMT& pmt = pmtRouter.getPMT(address);
202  const double t1 = converter.putTime(getTime(*i));
203 
204  zbuf[address.first].push_back(JHitL0(JDAQPMTIdentifier(id.getModuleID(), id.getPMTAddress()), pmt, t1));
205  }
206  }
207 
208  for (map<int, vector<JHitL0> >::iterator i = zbuf.begin(); i != zbuf.end(); ++i) {
209 
210  sort(i->second.begin(), i->second.end(), less<JHit>());
211 
212  data.push_back(i->second[0]);
213  }
214  }
215 
216 
217  // 3D cluster
218 
219  data.erase(clusterizeWeight(data.begin(), data.end(), match3B), data.end());
220 
221 
222  // 1D cluster
223 
224  for (JData_t::iterator i = data.begin(); i != data.end(); ++i) {
225  i->rotate(R);
226  }
227 
228  data.erase(clusterizeWeight(data.begin(), data.end(), match1D), data.end());
229 
230 
231  // move true muon to center of mass of hits
232 
233  tz.setZ(JWeighedCenter3D(data.begin(), data.end()).getZ(), getSpeedOfLight());
234 
235 
236  // set coordinate origin to true muon position
237 
238  for (JData_t::iterator i = data.begin(); i != data.end(); ++i) {
239  i->sub(tz.getPosition());
240  }
241 
242  tz.setPosition(JVector3D(0,0,0));
243 
244 
245  // apply rotational smearing to hits
246 
247  for (JData_t::iterator i = data.begin(); i != data.end(); ++i) {
248  i->rotate_back(Rs);
249  }
250 
251 
252  if (data.size() >= NUMBER_OF_HITS + numberOfOutliers) {
253 
254  JData_t::iterator __end = data.end();
255 
256  for (int n = 0; n < numberOfOutliers && distance(data.begin(), __end) > NUMBER_OF_HITS; ++n) {
257 
258  double ymax = 0.0;
259  JData_t::iterator kill = __end;
260 
261  for (JData_t::iterator i = data.begin(); i != __end; ++i) {
262 
263  const double y = fabs(i->getT() - tz.getT(*i)) / sigma_ns;
264 
265  if (y > ymax) {
266  ymax = y;
267  kill = i;
268  }
269  }
270 
271  if (ymax >= STANDARD_DEVIATIONS)
272  iter_swap(kill, --__end);
273  else
274  break;
275  }
276 
277  const double chi2 = getChi2(tz, data.begin(), __end, sigma_ns);
278  const int N = distance(data.begin(), __end);
279 
280  h0.Fill(chi2 / N);
281  p0.Fill(TMath::Prob(chi2, N));
282  n0.Fill((double) data.size(), (double) N / (double) data.size());
283  }
284 
285 
286  if (data.size() >= NUMBER_OF_HITS + numberOfOutliers) {
287 
288  JMatrixNZ V;
289  JVectorNZ Y;
290 
291  V.set(tz, data.begin(), data.end(), alpha_deg, sigma_ns);
292  Y.set(tz, data.begin(), data.end());
293 
294  V.invert();
295 
296  unsigned int N = data.size();
297 
298  for (int n = 0; n < numberOfOutliers && N > NUMBER_OF_HITS; ++n, --N) {
299 
300  double ymax = 0.0;
301  int kill = -1;
302 
303  for (unsigned int i = 0; i != V.size(); ++i) {
304 
305  const double y = getChi2(Y, V, i);
306 
307  if (y > ymax) {
308  ymax = y;
309  kill = i;
310  }
311  }
312 
313  if (ymax >= STANDARD_DEVIATIONS * STANDARD_DEVIATIONS)
314  V.update(kill, HIT_OFF);
315  else
316  break;
317  }
318 
319  const double chi2 = V.getDot(Y);
320 
321  h1.Fill(chi2 / N);
322  p1.Fill(TMath::Prob(chi2, N));
323  n1.Fill((double) data.size(), (double) N / (double) data.size());
324  }
325  }
326  }
327  STATUS(endl);
328 
329  out.Write();
330  out.Close();
331 }
Auxiliary methods to evaluate Poisson probabilities and chi2.
Router for direct addressing of PMT data in detector data structure.
Definition: JPMTRouter.hh:33
Data structure for angles in three dimensions.
Definition: JAngle3D.hh:33
Utility class to parse command line options.
Definition: JParser.hh:1500
General exception.
Definition: JException.hh:23
1D match criterion.
Definition: JMatch1D.hh:31
ROOT TTree parameter settings of various packages.
TPaveText * p1
Match operator for Cherenkov light from muon with given direction.
void set(const JVector3D &pos, T __begin, T __end, const double alpha, const double sigma)
Set co-variance matrix.
Definition: JMatrixNZ.hh:85
Algorithms for hit clustering and sorting.
Synchronously read DAQ events and Monte Carlo events (and optionally other events).
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
JPMTIdentifier getIdentifier(const JPMTAddress &address) const
Get identifier of PMT.
Definition: JPMTRouter.hh:126
static double getDot(const variance &first, const variance &second)
Get dot product.
Definition: JMatrixNZ.hh:196
#define STATUS(A)
Definition: JMessage.hh:63
Detector data structure.
Definition: JDetector.hh:80
bool is_muon(const Trk &track)
Test whether given track is a (anti-)muon.
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.
Rotation matrix.
Definition: JRotation3D.hh:111
Auxiliary class to synchronously read DAQ events and Monte Carlo events (and optionally other events)...
then for HISTOGRAM in h0 h1
Definition: JMatrixNZ.sh:69
const JPMT & getPMT(const JPMTAddress &address) const
Get PMT.
Definition: JPMTRouter.hh:90
JVertex3D & add(const JVertex3D &value)
Addition operator.
Definition: JVertex3D.hh:85
Auxiliary class to convert DAQ hit time to/from Monte Carlo hit time.
then fatal Wrong number of arguments fi set_variable STRING $argv[1] set_variable DETECTORXY_TXT $WORKDIR $DETECTORXY_TXT tail read X Y CHI2 RMS printf optimum n $X $Y $CHI2 $RMS awk v Y
double getTime(const Hit &hit)
Get true time of hit.
void update(const size_t k, const double value)
Update inverted matrix at given diagonal element.
Definition: JMatrixNS.hh:456
Basic data structure for time and time over threshold information of hit.
string outputFile
static struct JTRIGGER::@76 clusterizeWeight
Anonymous struct for weighed clustering of hits.
Data structure for detector geometry and calibration.
bool is_noise(const Hit &hit)
Verify hit origin.
Head getHeader(const JMultipleFileScanner_t &file_list)
Get Monte Carlo header.
Basic data structure for L0 hit.
int first
index of module in detector data structure
Auxiliary class for defining the range of iterations of objects.
Definition: JLimit.hh:41
Template L2 builder.
Definition: JBuildL2.hh:45
Determination of the time residual vector of hits for a track along z-axis (JFIT::JLine1Z).
Definition: JVectorNZ.hh:21
Detector file.
Definition: JHead.hh:196
Match operator for Cherenkov light from muon in any direction.
JDirection3D getDirection(const Vec &dir)
Get direction.
Data structure for vector in three dimensions.
Definition: JVector3D.hh:34
Determination of the co-variance matrix of hits for a track along z-axis (JFIT::JLine1Z).
Definition: JMatrixNZ.hh:28
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1961
Data structure for PMT geometry and calibration.
Definition: JPMT.hh:47
JPosition3D getPosition(const Vec &pos)
Get position.
Definition of hit and track types and auxiliary methods for handling Monte Carlo data.
Physics constants.
double putTime() const
Get Monte Carlo time minus DAQ/trigger time.
static const double PI
Mathematical constants.
void set(const JLine1Z &track, T __begin, T __end)
Set time residual vector.
Definition: JVectorNZ.hh:68
Direct access to PMT in detector data structure.
int debug
debug level
Definition: JSirene.cc:63
then usage $script[distance] fi case set_variable R
Definition: JDrawLED.sh:40
General purpose messaging.
Auxiliary include file for time conversion between DAQ/trigger hit and Monte Carlo hit...
#define FATAL(A)
Definition: JMessage.hh:67
void invert()
Invert matrix according LDU decomposition.
Definition: JMatrixNS.hh:80
Direct access to module in detector data structure.
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:157
Data structure for L2 parameters.
const double getSpeedOfLight()
Get speed of light.
Address of PMT in detector data structure.
Definition: JPMTAddress.hh:32
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
Utility class to parse command line options.
Data structure for L0 hit.
Definition: JHitL0.hh:27
alias put_queue eval echo n
Definition: qlib.csh:19
Data structure for fit of straight line paralel to z-axis.
Definition: JLine1Z.hh:27
Data structure for position in three dimensions.
Definition: JPosition3D.hh:36
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:73
do set_variable DETECTOR_TXT $WORKDIR detector
General purpose class for multiple pointers.
double getChi2(const double P)
Get chi2 corresponding to given probability.
Definition: JFitToolkit.hh:79
JPosition3D & rotate(const JRotation3D &R)
Rotate.
Definition: JPosition3D.hh:186
then usage $script[input file[working directory[option]]] nWhere option can be N
Definition: JMuonPostfit.sh:37
const JPMTAddress & getAddress(const JObjectID &id) const
Get address of PMT.
Definition: JPMTRouter.hh:78
Basic data structure for L1 hit.
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
Definition: Evt.hh:19
3D match criterion with road width.
Definition: JMatch3B.hh:34
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
int main(int argc, char *argv[])
Definition: Main.cpp:15