63 JTriggeredFileScanner<> inputFile;
78 JParser<> zap(
"Example program to test chi2 calculations of co-variance matrix including direction uncertainty.");
83 zap[
'n'] =
make_field(numberOfEvents) = JLimit::max();
95 catch(
const exception& error) {
96 FATAL(error.what() << endl);
105 const unsigned int NUMBER_OF_HITS = 6;
106 const double STANDARD_DEVIATIONS = 3.0;
107 const double HIT_OFF = 1.0e3 * sigma_ns*sigma_ns;
113 load(detectorFile, detector);
115 catch(
const JException& error) {
119 const JModuleRouter moduleRouter(detector);
120 const JPMTRouter pmtRouter (detector);
122 const JPosition3D center = get<JPosition3D>(
getHeader(inputFile));
128 TH1D h0(
"h0", NULL, 40, 0.0, 20.0);
129 TH1D h1(
"h1", NULL, 40, 0.0, 20.0);
131 TH1D p0(
"p0", NULL, 20, 0.0, 1.0);
132 TH1D
p1(
"p1", NULL, 20, 0.0, 1.0);
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); }
144 TProfile n0(
"n0", NULL, X.size() - 1, X.data());
145 TProfile n1(
"n1", NULL, X.size() - 1, X.data());
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);
154 while (inputFile.hasNext()) {
156 STATUS(
"event: " << setw(10) << inputFile.getCounter() <<
'\r');
DEBUG(endl);
158 JTriggeredFileScanner<>::multi_pointer_type ps = inputFile.next();
161 const Evt*
event = ps;
163 const JTimeConverter converter(*event, *tev);
167 if (muon != event->mc_trks.end()) {
171 JLine1Z tz(
getPosition(*muon), converter.putTime(muon->t));
176 const double theta = alpha_deg *
PI / 180.0;
177 const double phi = gRandom->Uniform(0.0, 2*
PI);
179 const JRotation3D Rs(JAngle3D(theta,phi));
186 buildL2(*tev, moduleRouter, back_inserter(data));
188 data.erase(unique(data.begin(), data.end(), equal_to<JDAQModuleIdentifier>()), data.end());
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));
203 zbuf[address.first].push_back(JHitL0(
JDAQPMTIdentifier(
id.getModuleID(),
id.getPMTAddress()), pmt, t1));
209 sort(i->second.begin(), i->second.end(), less<JHit>());
211 data.push_back(i->second[0]);
218 data.erase(
clusterizeWeight(data.begin(), data.end(), match3B), data.end());
223 for (JData_t::iterator i = data.begin(); i != data.end(); ++i) {
227 data.erase(
clusterizeWeight(data.begin(), data.end(), match1D), data.end());
232 tz.setZ(JWeighedCenter3D(data.begin(), data.end()).getZ(),
getSpeedOfLight());
237 for (JData_t::iterator i = data.begin(); i != data.end(); ++i) {
238 i->sub(tz.getPosition());
241 tz.setPosition(JVector3D(0,0,0));
246 for (JData_t::iterator i = data.begin(); i != data.end(); ++i) {
251 if (data.size() >= NUMBER_OF_HITS + numberOfOutliers) {
253 JData_t::iterator __end = data.end();
255 for (
int n = 0;
n < numberOfOutliers &&
distance(data.begin(), __end) > NUMBER_OF_HITS; ++
n) {
258 JData_t::iterator kill = __end;
260 for (JData_t::iterator i = data.begin(); i != __end; ++i) {
262 const double y = fabs(i->getT() - tz.getT(*i)) / sigma_ns;
270 if (ymax >= STANDARD_DEVIATIONS)
271 iter_swap(kill, --__end);
276 const double chi2 =
getChi2(tz, data.begin(), __end, sigma_ns);
277 const int N =
distance(data.begin(), __end);
280 p0.Fill(TMath::Prob(chi2, N));
281 n0.Fill((
double) data.size(), (double) N / (
double) data.size());
285 if (data.size() >= NUMBER_OF_HITS + numberOfOutliers) {
290 V.set(tz, data.begin(), data.end(), alpha_deg, sigma_ns);
291 Y.set(tz, data.begin(), data.end());
295 unsigned int N = data.size();
297 for (
int n = 0; n < numberOfOutliers && N > NUMBER_OF_HITS; ++
n, --N) {
302 for (
unsigned int i = 0; i != V.size(); ++i) {
304 const double y =
getChi2(Y, V, i);
312 if (ymax >= STANDARD_DEVIATIONS * STANDARD_DEVIATIONS)
313 V.update(kill, HIT_OFF);
318 const double chi2 =
getDot(Y.begin(), Y.end(), V);
321 p1.Fill(TMath::Prob(chi2, N));
322 n1.Fill((
double) data.size(), (double) N / (
double) data.size());