Jpp  17.1.0
the software that should make you happy
 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, 50, 0.0, 20.0);
130  TH1D h1("h1", NULL, 50, 0.0, 20.0);
131 
132  TH1D p0("p0", NULL, 50, 0.0, 1.0);
133  TH1D p1("p1", NULL, 50, 0.0, 1.0);
134 
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  JHit::setSlewing(!useTrue);
155 
156  while (inputFile.hasNext()) {
157 
158  STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);
159 
160  JTriggeredFileScanner<>::multi_pointer_type ps = inputFile.next();
161 
162  const JDAQEvent* tev = ps;
163  const Evt* event = ps;
164 
165  const time_converter converter(*event, *tev);
166 
167  vector<Trk>::const_iterator muon = find_if(event->mc_trks.begin(), event->mc_trks.end(), is_muon);
168 
169  if (muon != event->mc_trks.end()) {
170 
171  const JRotation3D R(getDirection(*muon));
172 
173  JLine1Z tz(getPosition(*muon), converter.putTime(muon->t));
174 
175  tz.add(center);
176  tz.rotate(R);
177 
178  const double theta = alpha_deg * PI / 180.0;
179  const double phi = gRandom->Uniform(0.0, 2*PI);
180 
181  const JRotation3D Rs(JAngle3D(theta,phi));
182 
183 
184  JData_t data;
185 
186  if (!useTrue) {
187 
188  buildL2(*tev, moduleRouter, back_inserter(data));
189 
190  data.erase(unique(data.begin(), data.end(), equal_to<JDAQModuleIdentifier>()), data.end());
191 
192  } else {
193 
195 
196  for (vector<Hit>::const_iterator i = event->mc_hits.begin(); i != event->mc_hits.end(); ++i) {
197 
198  if (!is_noise(*i)) {
199 
200  const JPMTAddress& address = pmtRouter.getAddress(i->pmt_id);
201  const JPMTIdentifier& id = pmtRouter.getIdentifier(address);
202  const JPMT& pmt = pmtRouter.getPMT(address);
203  const double t1 = converter.putTime(getTime(*i));
204 
205  zbuf[address.first].push_back(JHitL0(JDAQPMTIdentifier(id.getModuleID(), id.getPMTAddress()), pmt, t1));
206  }
207  }
208 
209  for (map<int, vector<JHitL0> >::iterator i = zbuf.begin(); i != zbuf.end(); ++i) {
210 
211  sort(i->second.begin(), i->second.end(), less<JHit>());
212 
213  data.push_back(i->second[0]);
214  }
215  }
216 
217 
218  // 3D cluster
219 
220  data.erase(clusterizeWeight(data.begin(), data.end(), match3B), data.end());
221 
222 
223  // 1D cluster
224 
225  for (JData_t::iterator i = data.begin(); i != data.end(); ++i) {
226  i->rotate(R);
227  }
228 
229  data.erase(clusterizeWeight(data.begin(), data.end(), match1D), data.end());
230 
231 
232  // move true muon to center of mass of hits
233 
234  tz.setZ(JWeighedCenter3D(data.begin(), data.end()).getZ(), getSpeedOfLight());
235 
236 
237  // set coordinate origin to true muon position
238 
239  for (JData_t::iterator i = data.begin(); i != data.end(); ++i) {
240  i->sub(tz.getPosition());
241  }
242 
243  tz.setPosition(JVector3D(0,0,0));
244 
245 
246  // apply rotational smearing to hits
247 
248  for (JData_t::iterator i = data.begin(); i != data.end(); ++i) {
249  i->rotate_back(Rs);
250  }
251 
252 
253  if (data.size() >= NUMBER_OF_HITS + numberOfOutliers) {
254 
255  JData_t::iterator __end = data.end();
256 
257  for (int n = 0; n < numberOfOutliers && distance(data.begin(), __end) > NUMBER_OF_HITS; ++n) {
258 
259  double ymax = 0.0;
260  JData_t::iterator kill = __end;
261 
262  for (JData_t::iterator i = data.begin(); i != __end; ++i) {
263 
264  const double y = fabs(i->getT() - tz.getT(*i)) / sigma_ns;
265 
266  if (y > ymax) {
267  ymax = y;
268  kill = i;
269  }
270  }
271 
272  if (ymax >= STANDARD_DEVIATIONS)
273  iter_swap(kill, --__end);
274  else
275  break;
276  }
277 
278  const double chi2 = getChi2(tz, data.begin(), __end, sigma_ns);
279  const int N = distance(data.begin(), __end);
280 
281  h0.Fill(chi2 / N);
282  p0.Fill(TMath::Prob(chi2, N));
283  n0.Fill((double) data.size(), (double) N / (double) data.size());
284  }
285 
286 
287  if (data.size() >= NUMBER_OF_HITS + numberOfOutliers) {
288 
289  JMatrixNZ V;
290  JVectorNZ Y;
291 
292  V.set(tz, data.begin(), data.end(), alpha_deg, sigma_ns);
293  Y.set(tz, data.begin(), data.end());
294 
295  V.invert();
296 
297  unsigned int N = data.size();
298 
299  for (int n = 0; n < numberOfOutliers && N > NUMBER_OF_HITS; ++n, --N) {
300 
301  double ymax = 0.0;
302  int kill = -1;
303 
304  for (unsigned int i = 0; i != V.size(); ++i) {
305 
306  const double y = getChi2(Y, V, i);
307 
308  if (y > ymax) {
309  ymax = y;
310  kill = i;
311  }
312  }
313 
314  if (ymax >= STANDARD_DEVIATIONS * STANDARD_DEVIATIONS)
315  V.update(kill, HIT_OFF);
316  else
317  break;
318  }
319 
320  const double chi2 = V.getDot(Y);
321 
322  h1.Fill(chi2 / N);
323  p1.Fill(TMath::Prob(chi2, N));
324  n1.Fill((double) data.size(), (double) N / (double) data.size());
325  }
326  }
327  }
328  STATUS(endl);
329 
330  out.Write();
331  out.Close();
332 }
Auxiliary methods to evaluate Poisson probabilities and chi2.
Router for direct addressing of PMT data in detector data structure.
Definition: JPMTRouter.hh:35
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
int main(int argc, char *argv[])
Definition: Main.cc:15
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:128
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:89
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 JShowerPostfit f $INPUT_FILE o $OUTPUT_FILE N
const JPMT & getPMT(const JPMTAddress &address) const
Get PMT.
Definition: JPMTRouter.hh:92
JVertex3D & add(const JVertex3D &value)
Addition operator.
Definition: JVertex3D.hh:87
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
static struct JTRIGGER::@82 clusterizeWeight
Anonymous struct for weighed clustering of hits.
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:457
Basic data structure for time and time over threshold information of hit.
V(JDAQEvent-JTriggerReprocessor)*1.0/(JDAQEvent+1.0e-10)
string outputFile
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
const int n
Definition: JPolint.hh:676
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:226
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, calibration and status.
Definition: JPMT.hh:43
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.
function getChi2()
int debug
debug level
Definition: JSirene.cc:67
then usage $script[distance] fi case set_variable R
Definition: JDrawLED.sh:43
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:181
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
then if[[!-f $DETECTOR]] then JDetector sh $DETECTOR fi cat $WORKDIR trigger_parameters txt<< EOFtrigger3DMuon.enabled=1;trigger3DMuon.numberOfHits=5;trigger3DMuon.gridAngle_deg=1;ctMin=0.0;TMaxLocal_ns=15.0;EOF set_variable TRIGGEREFFICIENCY_TRIGGERED_EVENTS_ONLY INPUT_FILES=() for((i=1;$i<=$NUMBER_OF_RUNS;++i));do JSirene.sh $DETECTOR $JPP_DATA/genhen.km3net_wpd_V2_0.evt.gz $WORKDIR/sirene_ ${i}.root JTriggerEfficiency.sh $DETECTOR $DETECTOR $WORKDIR/sirene_ ${i}.root $WORKDIR/trigger_efficiency_ ${i}.root $WORKDIR/trigger_parameters.txt $JPP_DATA/PMT_parameters.txt INPUT_FILES+=($WORKDIR/trigger_efficiency_ ${i}.root) done for ANGLE_DEG in $ANGLES_DEG[*];do set_variable SIGMA_NS 3.0 set_variable OUTLIERS 3 set_variable OUTPUT_FILE $WORKDIR/matrix\[${ANGLE_DEG}\deg\].root $JPP_DIR/examples/JReconstruction-f"$INPUT_FILES[*]"-o $OUTPUT_FILE-S ${SIGMA_NS}-A ${ANGLE_DEG}-O ${OUTLIERS}-d ${DEBUG}--!fiif[[$OPTION=="plot"]];then if((0));then for H1 in h0 h1;do JPlot1D-f"$WORKDIR/matrix["${^ANGLES_DEG}" deg].root:${H1}"-y"1 2e3"-Y-L TR-T""-\^"number of events [a.u.]"-> o chi2
Definition: JMatrixNZ.sh:106
Data structure for fit of straight line paralel to z-axis.
Definition: JLine1Z.hh:27
no fit printf nominal n $STRING awk v X
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.
JPosition3D & rotate(const JRotation3D &R)
Rotate.
Definition: JPosition3D.hh:186
const JPMTAddress & getAddress(const JObjectID &id) const
Get address of PMT.
Definition: JPMTRouter.hh:80
Basic data structure for L1 hit.
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
Definition: Evt.hh:20
3D match criterion with road width.
Definition: JMatch3B.hh:34
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62