Jpp
Functions
JMatrixNZ.cc File Reference
#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <map>
#include <algorithm>
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TProfile.h"
#include "TRandom.h"
#include "evt/Head.hh"
#include "evt/Evt.hh"
#include "JDAQ/JDAQEvent.hh"
#include "JDAQ/JDAQTimeslice.hh"
#include "JAAnet/JHead.hh"
#include "JAAnet/JHeadToolkit.hh"
#include "JTools/JConstants.hh"
#include "JGeometry3D/JDirection3D.hh"
#include "JGeometry3D/JRotation3D.hh"
#include "JGeometry3D/JGeometry3DToolkit.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JModuleRouter.hh"
#include "JDetector/JPMTRouter.hh"
#include "JTrigger/JHit.hh"
#include "JTrigger/JFrame.hh"
#include "JTrigger/JTimeslice.hh"
#include "JTrigger/JHitL0.hh"
#include "JTrigger/JHitL1.hh"
#include "JTrigger/JBuildL1.hh"
#include "JTrigger/JBuildL2.hh"
#include "JTrigger/JAlgorithm.hh"
#include "JTrigger/JMatch1D.hh"
#include "JTrigger/JMatch3B.hh"
#include "JTrigger/JTimeConverter.hh"
#include "JSupport/JTriggeredFileScanner.hh"
#include "JSupport/JFileRecorder.hh"
#include "JSupport/JMonteCarloFileSupportkit.hh"
#include "JSupport/JSupport.hh"
#include "JFit/JMatrixNZ.hh"
#include "JFit/JFitToolkit.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

Go to the source code of this file.

Functions

int main (int argc, char **argv)
 

Detailed Description

Example program to test chi2 calculations based on JFIT::JMatrixNZ.

Author
mdejong

Definition in file JMatrixNZ.cc.

Function Documentation

◆ main()

int main ( int  argc,
char **  argv 
)

Definition at line 57 of file JMatrixNZ.cc.

58 {
59  using namespace std;
60  using namespace JPP;
61  using namespace KM3NETDAQ;
62 
63  JTriggeredFileScanner<> inputFile;
64  JLimit_t& numberOfEvents = inputFile.getLimit();
65  string outputFile;
66  string detectorFile;
67  double Tmax_ns;
68  double roadWidth_m;
69  double ctMin;
70  double alpha_deg;
71  double sigma_ns;
72  int numberOfOutliers;
73  bool useTrue;
74  int debug;
75 
76  try {
77 
78  JParser<> zap("Example program to test chi2 calculations of co-variance matrix including direction uncertainty.");
79 
80  zap['f'] = make_field(inputFile);
81  zap['o'] = make_field(outputFile) = "matrixnz.root";
82  zap['a'] = make_field(detectorFile);
83  zap['n'] = make_field(numberOfEvents) = JLimit::max();
84  zap['T'] = make_field(Tmax_ns) = 20.0; // [ns]
85  zap['R'] = make_field(roadWidth_m) = 150.0; // [m]
86  zap['C'] = make_field(ctMin) = 0.0; //
87  zap['A'] = make_field(alpha_deg);
88  zap['S'] = make_field(sigma_ns);
89  zap['O'] = make_field(numberOfOutliers);
90  zap['U'] = make_field(useTrue);
91  zap['d'] = make_field(debug) = 1;
92 
93  zap(argc, argv);
94  }
95  catch(const exception& error) {
96  FATAL(error.what() << endl);
97  }
98 
99 
100  using namespace KM3NETDAQ;
101 
102  cout.tie(&cerr);
103 
104 
105  const unsigned int NUMBER_OF_HITS = 6;
106  const double STANDARD_DEVIATIONS = 3.0; // [unit]
107  const double HIT_OFF = 1.0e3 * sigma_ns*sigma_ns; // [ns^2]
108 
109 
110  JDetector detector;
111 
112  try {
113  load(detectorFile, detector);
114  }
115  catch(const JException& error) {
116  FATAL(error);
117  }
118 
119  const JModuleRouter moduleRouter(detector);
120  const JPMTRouter pmtRouter (detector);
121 
122  const JPosition3D center = get<JPosition3D>(getHeader(inputFile));
123 
124 
125  TFile out(outputFile.c_str(), "recreate");
126 
127 
128  TH1D h0("h0", NULL, 40, 0.0, 20.0);
129  TH1D h1("h1", NULL, 40, 0.0, 20.0);
130 
131  TH1D p0("p0", NULL, 20, 0.0, 1.0);
132  TH1D p1("p1", NULL, 20, 0.0, 1.0);
133 
134  vector<double> X;
135 
136  {
137  double x = -0.5;
138 
139  for ( ; x < 10.0; x += 1.0) { X.push_back(x); }
140  for ( ; x < 25.0; x += 2.0) { X.push_back(x); }
141  for ( ; x < 100.0; x += 5.0) { X.push_back(x); }
142  }
143 
144  TProfile n0("n0", NULL, X.size() - 1, X.data());
145  TProfile n1("n1", NULL, X.size() - 1, X.data());
146 
147 
148  typedef vector<JHitL1> JData_t;
149  const JBuildL2<JHitL1> buildL2(JL2Parameters(2, Tmax_ns, ctMin));
150  const JMatch3B<JHitL1> match3B(roadWidth_m, Tmax_ns);
151  const JMatch1D<JHitL1> match1D(roadWidth_m, Tmax_ns);
152 
153 
154  while (inputFile.hasNext()) {
155 
156  STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);
157 
158  JTriggeredFileScanner<>::multi_pointer_type ps = inputFile.next();
159 
160  const JDAQEvent* tev = ps;
161  const Evt* event = ps;
162 
163  const JTimeConverter converter(*event, *tev);
164 
165  vector<Trk>::const_iterator muon = find_if(event->mc_trks.begin(), event->mc_trks.end(), is_muon);
166 
167  if (muon != event->mc_trks.end()) {
168 
169  const JRotation3D R(getDirection(*muon));
170 
171  JLine1Z tz(getPosition(*muon), converter.putTime(muon->t));
172 
173  tz.add(center);
174  tz.rotate(R);
175 
176  const double theta = alpha_deg * PI / 180.0;
177  const double phi = gRandom->Uniform(0.0, 2*PI);
178 
179  const JRotation3D Rs(JAngle3D(theta,phi));
180 
181 
182  JData_t data;
183 
184  if (!useTrue) {
185 
186  buildL2(*tev, moduleRouter, back_inserter(data));
187 
188  data.erase(unique(data.begin(), data.end(), equal_to<JDAQModuleIdentifier>()), data.end());
189 
190  } else {
191 
193 
194  for (vector<Hit>::const_iterator i = event->mc_hits.begin(); i != event->mc_hits.end(); ++i) {
195 
196  if (!is_noise(*i)) {
197 
198  const JPMTAddress& address = pmtRouter.getAddress(i->pmt_id);
199  const JPMTIdentifier& id = pmtRouter.getIdentifier(address);
200  const JPMT& pmt = pmtRouter.getPMT(address);
201  const double t1 = converter.putTime(getTime(*i));
202 
203  zbuf[address.first].push_back(JHitL0(JDAQPMTIdentifier(id.getModuleID(), id.getPMTAddress()), pmt, t1));
204  }
205  }
206 
207  for (map<int, vector<JHitL0> >::iterator i = zbuf.begin(); i != zbuf.end(); ++i) {
208 
209  sort(i->second.begin(), i->second.end(), less<JHit>());
210 
211  data.push_back(i->second[0]);
212  }
213  }
214 
215 
216  // 3D cluster
217 
218  data.erase(clusterizeWeight(data.begin(), data.end(), match3B), data.end());
219 
220 
221  // 1D cluster
222 
223  for (JData_t::iterator i = data.begin(); i != data.end(); ++i) {
224  i->rotate(R);
225  }
226 
227  data.erase(clusterizeWeight(data.begin(), data.end(), match1D), data.end());
228 
229 
230  // move true muon to center of mass of hits
231 
232  tz.setZ(JWeighedCenter3D(data.begin(), data.end()).getZ(), getSpeedOfLight());
233 
234 
235  // set coordinate origin to true muon position
236 
237  for (JData_t::iterator i = data.begin(); i != data.end(); ++i) {
238  i->sub(tz.getPosition());
239  }
240 
241  tz.setPosition(JVector3D(0,0,0));
242 
243 
244  // apply rotational smearing to hits
245 
246  for (JData_t::iterator i = data.begin(); i != data.end(); ++i) {
247  i->rotate_back(Rs);
248  }
249 
250 
251  if (data.size() >= NUMBER_OF_HITS + numberOfOutliers) {
252 
253  JData_t::iterator __end = data.end();
254 
255  for (int n = 0; n < numberOfOutliers && distance(data.begin(), __end) > NUMBER_OF_HITS; ++n) {
256 
257  double ymax = 0.0;
258  JData_t::iterator kill = __end;
259 
260  for (JData_t::iterator i = data.begin(); i != __end; ++i) {
261 
262  const double y = fabs(i->getT() - tz.getT(*i)) / sigma_ns;
263 
264  if (y > ymax) {
265  ymax = y;
266  kill = i;
267  }
268  }
269 
270  if (ymax >= STANDARD_DEVIATIONS)
271  iter_swap(kill, --__end);
272  else
273  break;
274  }
275 
276  const double chi2 = getChi2(tz, data.begin(), __end, sigma_ns);
277  const int N = distance(data.begin(), __end);
278 
279  h0.Fill(chi2 / N);
280  p0.Fill(TMath::Prob(chi2, N));
281  n0.Fill((double) data.size(), (double) N / (double) data.size());
282  }
283 
284 
285  if (data.size() >= NUMBER_OF_HITS + numberOfOutliers) {
286 
287  JMatrixNZ V;
288  JVectorNZ Y;
289 
290  V.set(tz, data.begin(), data.end(), alpha_deg, sigma_ns);
291  Y.set(tz, data.begin(), data.end());
292 
293  V.invert();
294 
295  unsigned int N = data.size();
296 
297  for (int n = 0; n < numberOfOutliers && N > NUMBER_OF_HITS; ++n, --N) {
298 
299  double ymax = 0.0;
300  int kill = -1;
301 
302  for (unsigned int i = 0; i != V.size(); ++i) {
303 
304  const double y = getChi2(Y, V, i);
305 
306  if (y > ymax) {
307  ymax = y;
308  kill = i;
309  }
310  }
311 
312  if (ymax >= STANDARD_DEVIATIONS * STANDARD_DEVIATIONS)
313  V.update(kill, HIT_OFF);
314  else
315  break;
316  }
317 
318  const double chi2 = getDot(Y.begin(), Y.end(), V);
319 
320  h1.Fill(chi2 / N);
321  p1.Fill(TMath::Prob(chi2, N));
322  n1.Fill((double) data.size(), (double) N / (double) data.size());
323  }
324  }
325  }
326  STATUS(endl);
327 
328  out.Write();
329  out.Close();
330 }
std::iterator
Definition: JSTDTypes.hh:18
KM3NETDAQ::JDAQEvent
DAQ Event.
Definition: JDAQEvent.hh:34
JTOOLS::getSpeedOfLight
const double getSpeedOfLight()
Number of bytes in a gigabyte.
Definition: JConstants.hh:89
JDETECTOR::load
void load(const JString &file_name, JDetector &detector)
Load detector from input file.
Definition: JDetectorToolkit.hh:456
JTOOLS::n
const int n
Definition: JPolint.hh:628
std::vector< double >
JSUPPORT::JLimit_t
JLimit JLimit_t
Type definition of limit.
Definition: JLimit.hh:215
JPARSER::JParser
Utility class to parse command line options.
Definition: JParser.hh:1493
distance
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
Definition: PhysicsEvent.hh:434
JAANET::is_muon
bool is_muon(const Trk &track)
Test whether given track is a (anti-)muon.
Definition: JAAnetToolkit.hh:367
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition: JAAnetToolkit.hh:37
JSUPPORT::getHeader
Head getHeader(const JMultipleFileScanner_t &file_list)
Get Monte Carlo header.
Definition: JMonteCarloFileSupportkit.hh:425
debug
int debug
debug level
Definition: JSirene.cc:59
JTRIGGER::clusterizeWeight
JHitIterator_t clusterizeWeight(JHitIterator_t __begin, JHitIterator_t __end, const JMatch< JHit_t > &match)
Partition data according given binary match operator.
Definition: JAlgorithm.hh:310
JSUPPORT::JLimit::getLimit
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:73
JAANET::getTime
double getTime(const Hit &hit)
Get true time of hit.
Definition: JAAnetToolkit.hh:87
std::map
Definition: JSTDTypes.hh:16
STATUS
#define STATUS(A)
Definition: JMessage.hh:63
JAANET::getDirection
JDirection3D getDirection(const Vec &v)
Get direction.
Definition: JAAnetToolkit.hh:221
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1954
KM3NETDAQ::JDAQPMTIdentifier
PMT identifier.
Definition: JDAQPMTIdentifier.hh:25
JFIT::getChi2
double getChi2(const double P)
Get chi2 corresponding to given probability.
Definition: JFitToolkit.hh:80
JAANET::getPosition
JPosition3D getPosition(const Vec &v)
Get position.
Definition: JAAnetToolkit.hh:197
DEBUG
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
std
Definition: jaanetDictionary.h:36
KM3NETDAQ
KM3NeT DAQ data structures and auxiliaries.
Definition: DataQueue.cc:39
JTOOLS::PI
static const double PI
Constants.
Definition: JConstants.hh:20
JASTRONOMY::getDot
double getDot(const JNeutrinoDirection &first, const JNeutrinoDirection &second)
Dot product.
Definition: JAstronomy.hh:409
p1
TPaveText * p1
Definition: JDrawModule3D.cc:35
FATAL
#define FATAL(A)
Definition: JMessage.hh:67
outputFile
string outputFile
Definition: JDAQTimesliceSelector.cc:37
JAANET::is_noise
bool is_noise(const Hit &hit)
Verify hit origin.
Definition: JAAnetToolkit.hh:121