Jpp test-rotations-old
the software that should make you happy
Loading...
Searching...
No Matches
JRECONSTRUCTION::JMuonPrefit Struct Reference

Wrapper class to make pre-fit of muon trajectory. More...

#include <JMuonPrefit.hh>

Inheritance diagram for JRECONSTRUCTION::JMuonPrefit:
JRECONSTRUCTION::JMuonPrefitParameters_t JFIT::JEstimator< JModel_t > TObject

Classes

struct  cmz
 Auxiliary data structure for sorting of hits. More...
 
struct  input_type
 Input data type. More...
 

Public Types

typedef JEstimator< JLine1ZJEstimator_t
 
typedef JTRIGGER::JHitR1 hit_type
 
typedef std::vector< hit_typebuffer_type
 

Public Member Functions

 JMuonPrefit (const JMuonPrefitParameters_t &parameters, const int debug=0)
 Constructor.
 
 JMuonPrefit (const JMuonPrefitParameters_t &parameters, const JOmega3D &omega, const int debug=0)
 Constructor.
 
input_type getInput (const JModuleRouter &router, const JDAQEvent &event, const coverage_type &coverage) const
 Get input data.
 
JEvt operator() (const input_type &input)
 Fit function.
 
void reset ()
 Reset fit parameters.
 
bool equals (const JMuonPrefitParameters_t &parameters) const
 Equality.
 
 ClassDef (JMuonPrefitParameters_t, 1)
 

Public Attributes

JOmega3D omega
 
int debug
 
int factoryLimit
 factory limit for combinatorics
 
int NMaxHits
 maximal number of hits
 
double sigma_ns
 time resolution [ns]
 
double gridAngle_deg
 grid angle for directions [deg]
 
bool useL0
 option for L0 hit use
 
int numberOfOutliers
 maximum number of outliers
 
size_t numberOfPrefits
 number of prefits
 
double DZMax
 maximal slope for downward pointing solutions
 
size_t numberOfDZMax
 additional number of downward pointing solutions
 
double TMaxLocal_ns
 time window for local coincidences [ns]
 
double ctMin
 minimal cosine space angle between PMT axes
 
double roadWidth_m
 road width [m]
 
double Qmin
 minimal quality step
 

Static Public Attributes

static const struct JRECONSTRUCTION::JMuonPrefit::cmz cmz
 

Private Attributes

buffer_type data
 
JMatrixNZ V
 
JVectorNZ Y
 

Detailed Description

Wrapper class to make pre-fit of muon trajectory.

The JMuonPrefit fit is used to generate start values for subsequent fits (usually JMuonSimplex and JMuonGandalf).
To this end, a scan of directions is made and the time and transverse positions of the track are fitted for each direction (JFIT::JEstimator<JLine1Z>).
The directions are spaced by the parameters JMuonPrefitParameters_t::gridAngle_deg.
This angle corresponds to the envisaged angular accuracy of the result.
The probability that one of the results is less than this angle away from the correct value, multiple start values should be considered (JMuonPrefitParameters_t::numberOfPrefits).
Note that the CPU time scales with the inverse of the square of this angle.
The chi-squared is based on the time residuals.

Definition at line 79 of file JMuonPrefit.hh.

Member Typedef Documentation

◆ JEstimator_t

◆ hit_type

◆ buffer_type

Constructor & Destructor Documentation

◆ JMuonPrefit() [1/2]

JRECONSTRUCTION::JMuonPrefit::JMuonPrefit ( const JMuonPrefitParameters_t & parameters,
const int debug = 0 )
inline

Constructor.

Parameters
parametersparameters
debugdebug

Definition at line 126 of file JMuonPrefit.hh.

127 :
128 JMuonPrefitParameters_t(parameters),
129 omega (parameters.gridAngle_deg * JMATH::PI/180.0),
130 debug (debug)
131 {}
static const double PI
Mathematical constants.
double gridAngle_deg
grid angle for directions [deg]

◆ JMuonPrefit() [2/2]

JRECONSTRUCTION::JMuonPrefit::JMuonPrefit ( const JMuonPrefitParameters_t & parameters,
const JOmega3D & omega,
const int debug = 0 )
inline

Constructor.

Parameters
parametersparameters
omegadirections
debugdebug

Definition at line 141 of file JMuonPrefit.hh.

143 :
144 JMuonPrefitParameters_t(parameters),
145 omega (omega),
146 debug (debug)
147 {}

Member Function Documentation

◆ getInput()

input_type JRECONSTRUCTION::JMuonPrefit::getInput ( const JModuleRouter & router,
const JDAQEvent & event,
const coverage_type & coverage ) const
inline

Get input data.

Parameters
routermodule router
eventevent
coveragecoverage
Returns
input data

Definition at line 158 of file JMuonPrefit.hh.

162 {
163 using namespace std;
164 using namespace JTRIGGER;
165 using namespace KM3NETDAQ;
166
167 const JBuildL0<hit_type> buildL0;
169
170 input_type input(event.getDAQEventHeader(), coverage);
171
172 buffer_type& dataL0 = input.dataL0;
173 buffer_type& dataL1 = input.dataL1;
174
175 buildL2(event, router, !useL0, back_inserter(dataL1));
176
177 // 3D cluster of unique optical modules
178
180
181 sort(dataL1.begin(), dataL1.end(), hit_type::compare);
182
183 buffer_type::iterator __end = dataL1.end();
184
185 __end = unique(dataL1.begin(), __end, equal_to<JDAQModuleIdentifier>());
186
187 __end = clusterizeWeight(dataL1.begin(), __end, match3B);
188
189 dataL1.erase(__end, dataL1.end());
190
191
192 if (useL0) {
193
194 buildL0(event, router, true, back_inserter(dataL0));
195
196 __end = dataL0.end();
197
198 for (buffer_type::iterator i = dataL0.begin(); i != __end; ) {
199
200 if (match3B.count(*i, dataL1.begin(), dataL1.end()) != 0)
201 ++i;
202 else
203 swap(*i, *--__end);
204 }
205
206 dataL0.erase(__end, dataL0.end());
207 }
208
209 return input;
210 }
Template L0 hit builder.
Definition JBuildL0.hh:38
Template L2 builder.
Definition JBuildL2.hh:49
static const struct JTRIGGER::JHitR1::compare compare
3D match criterion with road width.
Definition JMatch3B.hh:36
const JDAQEventHeader & getDAQEventHeader() const
Get DAQ event header.
Auxiliary classes and methods for triggering.
JHitIterator_t clusterizeWeight(JHitIterator_t __begin, JHitIterator_t __end, const JMatch_t &match)
Partition data according given binary match operator.
KM3NeT DAQ data structures and auxiliaries.
Definition DataQueue.cc:39
double ctMin
minimal cosine space angle between PMT axes
double TMaxLocal_ns
time window for local coincidences [ns]
Data structure for L2 parameters.

◆ operator()()

JEvt JRECONSTRUCTION::JMuonPrefit::operator() ( const input_type & input)
inline

Fit function.

Parameters
inputinput data
Returns
fit results

Definition at line 219 of file JMuonPrefit.hh.

220 {
221 using namespace std;
222 using namespace JPP;
223
224 const double STANDARD_DEVIATIONS = 3.0; // [unit]
225 const double HIT_OFF = 1.0e3 * sigma_ns * sigma_ns; // [ns^2]
226
228
229 const buffer_type& dataL0 = input.dataL0;
230 const buffer_type& dataL1 = input.dataL1;
231
232 data.reserve(dataL0.size() +
233 dataL1.size());
234
235 JEvent event(JMUONPREFIT);
236
237 JEvt out;
238
239 for (JOmega3D_t::const_iterator dir = omega.begin(); dir != omega.end(); ++dir) {
240
241 const JRotation3D R(*dir);
242
243
244 buffer_type::iterator __end = copy(dataL1.begin(), dataL1.end(), data.begin());
245
246 for (buffer_type::iterator i = data.begin(); i != __end; ++i) {
247 i->rotate(R);
248 }
249
250
251 // reduce data
252
253 if (distance(data.begin(), __end) > NMaxHits) {
254
255 buffer_type::iterator __p = data.begin();
256
257 advance(__p, NMaxHits);
258
259 partial_sort(data.begin(), __p, __end, cmz);
260
261 __end = __p;
262 }
263
264
265 // 1D cluster
266
267 __end = clusterizeWeight(data.begin(), __end, match1D);
268
269 if (useL0) {
270
271 buffer_type::iterator p = __end; // begin L0 data
272 buffer_type::iterator q = copy(dataL0.begin(), dataL0.end(), p); // end L0 data
273
274 for (buffer_type::iterator i = p; i != q; ++i) {
275
276 if (find_if(data.begin(), __end, make_predicate(&hit_type::getModuleID, i->getModuleID())) == __end) {
277
278 i->rotate(R);
279
280 if (match1D.count(*i, data.begin(), __end) != 0) {
281 *p = *i;
282 ++p;
283 }
284 }
285 }
286
287 __end = clusterize(__end, p, match1D);
288 }
289
290
291 if (distance(data.begin(), __end) <= NUMBER_OF_PARAMETERS) {
292 continue;
293 }
294
295
296 // 1D fit
297
298 JLine1Z tz;
299 double chi2 = numeric_limits<double>::max();
300 int NDF = distance(data.begin(), __end) - NUMBER_OF_PARAMETERS;
301 int N = getCount(data.begin(), __end);
302
303
304 if (distance(data.begin(), __end) <= factoryLimit) {
305
306 int number_of_outliers = numberOfOutliers;
307
308 if (number_of_outliers > NDF - 1) {
309 number_of_outliers = NDF - 1;
310 }
311
312 double ymin = numeric_limits<double>::max();
313
314 buffer_type::iterator __end1 = __end;
315
316 for (int n = 0; n <= number_of_outliers; ++n, --__end1) {
317
318 sort(data.begin(), __end, hit_type::compare);
319
320 do {
321 /*
322 if (getNumberOfStrings(router, data.begin(), __end1) < 2) {
323 continue;
324 }
325 */
326 try {
327
328 (*this)(data.begin(), __end1);
329
330 V.set(*this, data.begin(), __end1, gridAngle_deg, sigma_ns);
331 Y.set(*this, data.begin(), __end1);
332
333 V.invert();
334
335 double y = getChi2(Y, V);
336
337 if (y <= -(STANDARD_DEVIATIONS * STANDARD_DEVIATIONS)) {
338
339 WARNING(endl << "chi2(1) " << y << endl);
340
341 } else {
342
343 if (y < 0.0) {
344 y = 0.0;
345 }
346
347 if (y < ymin) {
348 ymin = y;
349 tz = *this;
350 chi2 = ymin;
351 NDF = distance(data.begin(), __end1) - NUMBER_OF_PARAMETERS;
352 N = getCount(data.begin(), __end1);
353 }
354 }
355 }
356 catch(const exception& error) {}
357
358 } while (next_permutation(data.begin(), __end1, __end, hit_type::compare));
359
360 ymin -= STANDARD_DEVIATIONS * STANDARD_DEVIATIONS;
361 }
362
363 } else {
364
365 const int number_of_outliers = NDF - 1;
366
367 try {
368
369 (*this)(data.begin(), __end);
370
371 V.set(*this, data.begin(), __end, gridAngle_deg, sigma_ns);
372 Y.set(*this, data.begin(), __end);
373
374 V.invert();
375
376 for (int n = 0; n <= number_of_outliers; ++n) {
377
378 double ymax = 0.0;
379 int k = -1;
380
381 for (size_t i = 0; i != Y.size(); ++i) {
382
383 double y = getChi2(Y, V, i);
384
385 if (y > ymax) {
386 ymax = y;
387 k = i;
388 }
389 }
390
391 if (ymax < STANDARD_DEVIATIONS * STANDARD_DEVIATIONS) {
392 break;
393 }
394
395 try {
396
397 V.update(k, +HIT_OFF);
398
399 this->update(data.begin(), __end, V);
400
401 Y.set(*this, data.begin(), __end);
402
403 tz = *this;
404 NDF -= 1;
405 N -= getCount(data[k]);
406 }
407 catch(const exception& error) {
408
409 V.update(k, -HIT_OFF);
410
411 static_cast<JLine1Z&>(*this) = tz;
412
413 Y.set(*this, data.begin(), __end);
414
415 break;
416 }
417 }
418
419 chi2 = getChi2(Y, V);
420 tz = *this;
421 }
422 catch(const exception& error) {}
423 }
424
425 if (chi2 != numeric_limits<double>::max()) {
426
427 tz.rotate_back(R);
428
429 out.push_back(getFit(event(), tz, *dir, getQuality(chi2, N, NDF), NDF));
430
431 // set additional values
432
433 out.rbegin()->setW(JPP_COVERAGE_ORIENTATION, input.coverage.orientation);
434 out.rbegin()->setW(JPP_COVERAGE_POSITION, input.coverage.position);
435 }
436 }
437
438
439 if (numberOfPrefits > 0) {
440
441 JEvt::iterator __end = out.end();
442
443 if (Qmin > 0) { // sort distinct maxima
444
445 __end = out.begin(); // output
446
447 for (JEvt::iterator p1 = out.begin(); p1 != out.end() && (size_t) distance(out.begin(), __end) < numberOfPrefits; ) {
448
449 JEvt::iterator p2 = p1;
450
451 for (JEvt::iterator i = p1; i != out.end(); ++i) {
452 if (i->getQ() > p2->getQ()) {
453 p2 = i;
454 }
455 }
456
457 swap(*p2, *p1);
458
459 p2 = p1++;
460
461 sort(p1, out.end(), JPointing(*p2));
462
463 for (double Q = p2->getQ(); p1 != out.end() && (p1->getQ() >= p2->getQ() - Qmin || p1->getQ() <= Q); Q = (p1++)->getQ()) {}
464
465 swap(*(__end++), *p2);
466 }
467
468 } else if (numberOfPrefits < out.size()) { // sort subset
469
470 advance(__end = out.begin(), numberOfPrefits);
471
472 partial_sort(out.begin(), __end, out.end(), qualitySorter);
473
474 } else { // sort all
475
476 sort(out.begin(), __end, qualitySorter);
477 }
478
479 // add downward pointing solutions if available but not yet sufficient
480
481 int nz = numberOfDZMax - count_if(out.begin(), __end, make_predicate(&JFit::getDZ, DZMax, JComparison::le()));
482
483 if (nz > 0) {
484
485 JEvt::iterator __p = __end;
486 JEvt::iterator __q = __end = partition(__p, out.end(), make_predicate(&JFit::getDZ, DZMax, JComparison::le()));
487
488 if (nz < distance(__p, __q)) {
489
490 advance(__end = __p, nz);
491
492 partial_sort(__p, __end, __q, qualitySorter);
493
494 } else {
495
496 sort(__p, __end, qualitySorter);
497 }
498 }
499
500 out.erase(__end, out.end());
501
502 } else {
503
504 sort(out.begin(), out.end(), qualitySorter);
505 }
506
507 return out;
508 }
TPaveText * p1
#define WARNING(A)
Definition JMessage.hh:65
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
double getDZ() const
Get Z-slope.
Data structure for fit of straight line paralel to z-axis.
Definition JLine1Z.hh:29
void set(const JVector3D &pos, T __begin, T __end, const double alpha, const double sigma)
Set co-variance matrix.
Definition JMatrixNZ.hh:85
void set(const JLine1Z &track, T __begin, T __end)
Set time residual vector.
Definition JVectorNZ.hh:68
JPosition3D & rotate_back(const JRotation3D &R)
Rotate back.
Auxiliary class to compare fit results with respect to a reference direction (e.g....
1D match criterion.
Definition JMatch1D.hh:33
int getModuleID() const
Get module identifier.
static const int JMUONPREFIT
static const int JPP_COVERAGE_POSITION
coverage of dynamic position calibration from any Jpp application
static const int JPP_COVERAGE_ORIENTATION
coverage of dynamic orientation calibration from any Jpp application
void copy(const Head &from, JHead &to)
Copy header from from to to.
Definition JHead.cc:163
double getChi2(const double P)
Get chi2 corresponding to given probability.
size_t getCount(const array_type< T > &buffer, const JCompare_t &compare)
Count number of unique values.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
double getQuality(const double chi2, const int N, const int NDF)
Get quality of fit.
bool qualitySorter(const JFit &first, const JFit &second)
Comparison of fit results.
JFit getFit(const JHistory &history, const JTrack3D &track, const double Q, const int NDF, const double energy=0.0, const int status=SINGLE_STAGE)
Get fit.
counter_type advance(counter_type &counter, const counter_type value, const counter_type limit=std::numeric_limits< counter_type >::max())
Advance counter.
const int n
Definition JPolint.hh:791
bool next_permutation(T __begin, T __last, T __end, JComparator_t compare, std::bidirectional_iterator_tag)
Implementation of method next_permutation for bidirectional iterators.
JHitIterator_t clusterize(JHitIterator_t __begin, JHitIterator_t __end, const JMatch_t &match, const int Nmin=1)
Partition data according given binary match operator.
Definition JAlgorithm.hh:38
Acoustic event fit.
void update(const size_t k, const double value)
Update inverted matrix at given diagonal element.
Definition JMatrixNS.hh:446
void invert()
Invert matrix according LDU decomposition.
Definition JMatrixNS.hh:75
double DZMax
maximal slope for downward pointing solutions
int factoryLimit
factory limit for combinatorics
size_t numberOfDZMax
additional number of downward pointing solutions
static const struct JRECONSTRUCTION::JMuonPrefit::cmz cmz

◆ reset()

void JRECONSTRUCTION::JMuonPrefitParameters_t::reset ( )
inlineinherited

Reset fit parameters.

Definition at line 39 of file JMuonPrefitParameters_t.hh.

40 {
41 factoryLimit = 8;
42 NMaxHits = 50;
43 sigma_ns = 5;
44 gridAngle_deg = 1;
45 useL0 = false;
47 numberOfPrefits = 12;
48 DZMax = 0.0;
49 numberOfDZMax = 1;
50 TMaxLocal_ns = 18.0;
51 ctMin = 0;
52 roadWidth_m = 200;
53 Qmin = 0.0;
54 }

◆ equals()

bool JRECONSTRUCTION::JMuonPrefitParameters_t::equals ( const JMuonPrefitParameters_t & parameters) const
inlineinherited

Equality.

Parameters
parametersfit parameters
Returns
true if equals; else false

Definition at line 62 of file JMuonPrefitParameters_t.hh.

63 {
64 return (this->factoryLimit == parameters.factoryLimit &&
65 this->NMaxHits == parameters.NMaxHits &&
66 this->sigma_ns == parameters.sigma_ns &&
67 this->gridAngle_deg == parameters.gridAngle_deg &&
68 this->useL0 == parameters.useL0 &&
69 this->numberOfOutliers == parameters.numberOfOutliers &&
70 this->numberOfPrefits == parameters.numberOfPrefits &&
71 this->DZMax == parameters.DZMax &&
72 this->numberOfDZMax == parameters.numberOfDZMax &&
73 this->TMaxLocal_ns == parameters.TMaxLocal_ns &&
74 this->roadWidth_m == parameters.roadWidth_m &&
75 this->Qmin == parameters.Qmin);
76 }

◆ ClassDef()

JRECONSTRUCTION::JMuonPrefitParameters_t::ClassDef ( JMuonPrefitParameters_t ,
1  )
inherited

Member Data Documentation

◆ cmz

const struct JRECONSTRUCTION::JMuonPrefit::cmz JRECONSTRUCTION::JMuonPrefit::cmz
static

◆ omega

JOmega3D JRECONSTRUCTION::JMuonPrefit::omega

Definition at line 533 of file JMuonPrefit.hh.

◆ debug

int JRECONSTRUCTION::JMuonPrefit::debug

Definition at line 534 of file JMuonPrefit.hh.

◆ data

buffer_type JRECONSTRUCTION::JMuonPrefit::data
private

Definition at line 537 of file JMuonPrefit.hh.

◆ V

JMatrixNZ JRECONSTRUCTION::JMuonPrefit::V
private

Definition at line 538 of file JMuonPrefit.hh.

◆ Y

JVectorNZ JRECONSTRUCTION::JMuonPrefit::Y
private

Definition at line 539 of file JMuonPrefit.hh.

◆ factoryLimit

int JRECONSTRUCTION::JMuonPrefitParameters_t::factoryLimit
inherited

factory limit for combinatorics

Definition at line 80 of file JMuonPrefitParameters_t.hh.

◆ NMaxHits

int JRECONSTRUCTION::JMuonPrefitParameters_t::NMaxHits
inherited

maximal number of hits

Definition at line 81 of file JMuonPrefitParameters_t.hh.

◆ sigma_ns

double JRECONSTRUCTION::JMuonPrefitParameters_t::sigma_ns
inherited

time resolution [ns]

Definition at line 82 of file JMuonPrefitParameters_t.hh.

◆ gridAngle_deg

double JRECONSTRUCTION::JMuonPrefitParameters_t::gridAngle_deg
inherited

grid angle for directions [deg]

Definition at line 83 of file JMuonPrefitParameters_t.hh.

◆ useL0

bool JRECONSTRUCTION::JMuonPrefitParameters_t::useL0
inherited

option for L0 hit use

Definition at line 84 of file JMuonPrefitParameters_t.hh.

◆ numberOfOutliers

int JRECONSTRUCTION::JMuonPrefitParameters_t::numberOfOutliers
inherited

maximum number of outliers

Definition at line 85 of file JMuonPrefitParameters_t.hh.

◆ numberOfPrefits

size_t JRECONSTRUCTION::JMuonPrefitParameters_t::numberOfPrefits
inherited

number of prefits

Definition at line 86 of file JMuonPrefitParameters_t.hh.

◆ DZMax

double JRECONSTRUCTION::JMuonPrefitParameters_t::DZMax
inherited

maximal slope for downward pointing solutions

Definition at line 87 of file JMuonPrefitParameters_t.hh.

◆ numberOfDZMax

size_t JRECONSTRUCTION::JMuonPrefitParameters_t::numberOfDZMax
inherited

additional number of downward pointing solutions

Definition at line 88 of file JMuonPrefitParameters_t.hh.

◆ TMaxLocal_ns

double JRECONSTRUCTION::JMuonPrefitParameters_t::TMaxLocal_ns
inherited

time window for local coincidences [ns]

Definition at line 89 of file JMuonPrefitParameters_t.hh.

◆ ctMin

double JRECONSTRUCTION::JMuonPrefitParameters_t::ctMin
inherited

minimal cosine space angle between PMT axes

Definition at line 90 of file JMuonPrefitParameters_t.hh.

◆ roadWidth_m

double JRECONSTRUCTION::JMuonPrefitParameters_t::roadWidth_m
inherited

road width [m]

Definition at line 91 of file JMuonPrefitParameters_t.hh.

◆ Qmin

double JRECONSTRUCTION::JMuonPrefitParameters_t::Qmin
inherited

minimal quality step

Definition at line 92 of file JMuonPrefitParameters_t.hh.


The documentation for this struct was generated from the following file: