Jpp  18.2.1-ARCA-DF-PATCH
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JFitK40.hh
Go to the documentation of this file.
1 #ifndef __JCALIBRATE_JGANDALFK40__
2 #define __JCALIBRATE_JGANDALFK40__
3 
4 #include <vector>
5 #include <map>
6 #include <memory>
7 #include <limits>
8 #include <cmath>
9 
11 
12 #include "JLang/JException.hh"
13 #include "JLang/JManip.hh"
14 
15 #include "JTools/JRange.hh"
16 
17 #include "JDetector/JModule.hh"
18 
19 #include "JMath/JVectorND.hh"
20 #include "JMath/JMatrixNS.hh"
21 #include "JMath/JMath.hh"
22 #include "JMath/JGauss.hh"
23 #include "JMath/JMathToolkit.hh"
24 
25 #include "JFit/JMEstimator.hh"
26 
28 #include "JCalibrate/JTDC_t.hh"
29 
30 
31 /**
32  * \author mdejong
33  */
34 
35 namespace JCALIBRATE {}
36 namespace JPP { using namespace JCALIBRATE; }
37 
38 namespace JCALIBRATE {
39 
42  using JDETECTOR::JModule;
43  using JMATH::JMath;
44  using JFIT::JMEstimator;
45 
46 
47  /**
48  * Fit options.
49  */
50  enum JOption_t {
51  FIT_PMTS_t = 1, //!< fit parameters of PMTs
52  FIT_PMTS_AND_ANGULAR_DEPENDENCE_t = 2, //!< fit parameters of PMTs and angular dependence of K40 rate
53  FIT_PMTS_AND_BACKGROUND_t = 3, //!< fit parameters of PMTs and background
54  FIT_PMTS_QE_FIXED_t = 4, //!< fit parameters of PMTs with QE fixed
55  FIT_MODEL_t = 5 //!< fit parameters of K40 rate and TTSs of PMTs
56  };
57 
58  static const int INVALID_INDEX = -1; //!< invalid index
59 
60 
61  /**
62  * Data structure for measured coincidence rate of pair of PMTs.
63  */
64  struct rate_type {
65  /**
66  * Default constructor.
67  */
69  dt_ns(0.0),
70  value(0.0),
71  error(0.0)
72  {}
73 
74 
75  /**
76  * Constructor.
77  *
78  * \param dt_ns time difference [ns]
79  * \param value value of rate [Hz/ns]
80  * \param error error of rate [Hz/ns]
81  */
82  rate_type(double dt_ns,
83  double value,
84  double error) :
85  dt_ns(dt_ns),
86  value(value),
87  error(error)
88  {}
89 
90  double dt_ns; //!< time difference [ns]
91  double value; //!< value of rate [Hz/ns]
92  double error; //!< error of rate [Hz/ns]
93  };
94 
95 
96  /**
97  * Data structure for measured coincidence rates of all pairs of PMTs in optical module.
98  */
99  struct data_type :
100  public std::map<pair_type, std::vector<rate_type> >
101  {};
102 
103 
104  /**
105  * Auxiliary class for fit parameter with optional limits.
106  */
107  class JParameter_t :
108  public JMath<JParameter_t>
109  {
110  public:
111  /**
112  * Fit options.
113  */
114  enum FIT_t {
115  FREE_t = 0, //!< free
116  FIXED_t //!< fixed
117  };
118 
119 
120  /**
121  * Type definition for range of parameter values.
122  */
124 
125 
126  /**
127  * Default constructor.
128  */
130  {
131  set(0.0);
132  }
133 
134 
135  /**
136  * Constructor.
137  *
138  * \param value value
139  * \param range range
140  */
141  JParameter_t(const double value,
143  range(range)
144  {
145  set(value);
146  }
147 
148 
149  /**
150  * Negate parameter.
151  *
152  * \return this parameter
153  */
155  {
156  set(-get());
157 
158  return *this;
159  }
160 
161 
162  /**
163  * Add parameter.
164  *
165  * \param parameter parameter
166  * \return this parameter
167  */
168  JParameter_t& add(const JParameter_t& parameter)
169  {
170  set(get() + parameter.get());
171 
172  return *this;
173  }
174 
175 
176  /**
177  * Subtract parameter.
178  *
179  * \param parameter parameter
180  * \return this parameter
181  */
182  JParameter_t& sub(const JParameter_t& parameter)
183  {
184  set(get() - parameter.get());
185 
186  return *this;
187  }
188 
189 
190  /**
191  * Scale parameter.
192  *
193  * \param factor multiplication factor
194  * \return this parameter
195  */
196  JParameter_t& mul(const double factor)
197  {
198  set(get() * factor);
199 
200  return *this;
201  }
202 
203 
204  /**
205  * Scale parameter.
206  *
207  * \param factor division factor
208  * \return this parameter
209  */
210  JParameter_t& div(const double factor)
211  {
212  set(get() / factor);
213 
214  return *this;
215  }
216 
217 
218  /**
219  * Scale parameter.
220  *
221  * \param first first parameter
222  * \param second second parameter
223  * \return this parameter
224  */
226  {
227  set(first.get() * second.get());
228 
229  return *this;
230  }
231 
232 
233  /**
234  * Check if parameter is free.
235  *
236  * \return true if free; else false
237  */
238  bool isFree() const
239  {
240  return option == FREE_t;
241  }
242 
243 
244  /**
245  * Check if parameter is fixed.
246  *
247  * \return true if fixed; else false
248  */
249  bool isFixed() const
250  {
251  return option == FIXED_t;
252  }
253 
254 
255  /**
256  * Check if parameter is bound.
257  *
258  * \return true if bound; else false
259  */
260  bool isBound() const
261  {
262  return range.is_valid();
263  }
264 
265 
266  /**
267  * Set current value.
268  */
269  void set()
270  {
271  option = FREE_t;
272  }
273 
274 
275  /**
276  * Fix current value.
277  */
278  void fix()
279  {
280  option = FIXED_t;
281  }
282 
283 
284  /**
285  * Get value.
286  *
287  * \return value
288  */
289  double get() const
290  {
291  if (isBound())
292  return range.getLowerLimit() + 0.5 * range.getLength() * (sin(value) + 1.0);
293  else
294  return value;
295  }
296 
297 
298  /**
299  * Set value.
300  *
301  * \param value value
302  */
303  void set(const double value)
304  {
305  if (isBound())
306  this->value = asin(2.0 * (range.constrain(value) - range.getLowerLimit()) / range.getLength() - 1.0);
307  else
308  this->value = value;
309 
310  set();
311  }
312 
313 
314  /**
315  * Fix value.
316  *
317  * \param value value
318  */
319  void fix(const double value)
320  {
321  set(value);
322 
323  fix();
324  }
325 
326 
327  /**
328  * Get derivative of value.
329  *
330  * \return derivative of value
331  */
332  double getDerivative() const
333  {
334  if (isBound())
335  return 0.5 * range.getLength() * cos(value);
336  else
337  return 1.0;
338  }
339 
340 
341  /**
342  * Type conversion operator.
343  *
344  * \return value
345  */
346  double operator()() const
347  {
348  return get();
349  }
350 
351 
352  /**
353  * Type conversion operator.
354  *
355  * \return value
356  */
357  operator double() const
358  {
359  return get();
360  }
361 
362 
363  /**
364  * Assignment operator.
365  *
366  * \param value value
367  * \return this parameter
368  */
369  JParameter_t& operator=(double value)
370  {
371  set(value);
372 
373  return *this;
374  }
375 
376 
377  /**
378  * Read parameter from input stream.
379  *
380  * \param in input stream
381  * \param object parameter
382  * \return input stream
383  */
384  friend inline std::istream& operator>>(std::istream& in, JParameter_t& object)
385  {
386  return in >> object.value;
387  }
388 
389 
390  /**
391  * Write parameter to output stream.
392  *
393  * \param out output stream
394  * \param object parameter
395  * \return output stream
396  */
397  friend inline std::ostream& operator<<(std::ostream& out, const JParameter_t& object)
398  {
399  using namespace std;
400 
401  out << FIXED(12,6) << object.get() << ' '
402  << setw(5) << (object.isFixed() ? "fixed" : " ") << ' ';
403 
404  if (object.isBound()) {
405  out << FIXED(12,6) << object.value << ' ';
406  out << FIXED(12,6) << object.range.getLowerLimit() << ' '
407  << FIXED(12,6) << object.range.getUpperLimit();
408  }
409 
410  return out;
411  }
412 
413 
414  double value = 0.0;
417  };
418 
419 
420  /**
421  * Fit parameters for single PMT.
422  */
424 
425  static constexpr double QE_MIN = 0.0; //!< minimal QE
426  static constexpr double QE_MAX = 2.0; //!< maximal QE
427  static constexpr double TTS_NS = 2.0; //!< start value transition-time spread [ns]
428 
429  /**
430  * Default constructor.
431  */
433  {
434  reset();
435  }
436 
437 
438  /**
439  * Reset.
440  */
441  void reset()
442  {
443  status = true;
444 
445  QE .set(0.0);
446  TTS.set(0.0);
447  t0 .set(0.0);
448  bg .set(0.0);
449  }
450 
451 
452  /**
453  * Get default values.
454  *
455  * \return parameters
456  */
458  {
460 
461  parameters.QE.range = JParameter_t::range_type(QE_MIN, QE_MAX);
462 
463  parameters.status = true;
464 
465  parameters.QE .set(1.0);
466  parameters.TTS.set(TTS_NS);
467  parameters.t0 .set(0.0);
468  parameters.bg .set(0.0);
469 
470  return parameters;
471  }
472 
473 
474  /**
475  * Get number of fit parameters.
476  *
477  * \return number of parameters
478  */
479  inline size_t getN() const
480  {
481  return ((QE. isFree() ? 1 : 0) +
482  (TTS.isFree() ? 1 : 0) +
483  (t0 .isFree() ? 1 : 0) +
484  (bg .isFree() ? 1 : 0));
485  }
486 
487 
488  /**
489  * Get index of parameter.
490  *
491  * \param p pointer to data member
492  * \return index
493  */
495  {
496  if (!(this->*p).isFree()) {
497  return INVALID_INDEX;
498  }
499 
500  int N = 0;
501 
502  if (p == &JPMTParameters_t::QE ) { return N; }; if (QE .isFree()) { ++N; }
503  if (p == &JPMTParameters_t::TTS) { return N; }; if (TTS.isFree()) { ++N; }
504  if (p == &JPMTParameters_t::t0 ) { return N; }; if (t0 .isFree()) { ++N; }
505  if (p == &JPMTParameters_t::bg ) { return N; }; if (bg .isFree()) { ++N; }
506 
507  return INVALID_INDEX;
508  }
509 
510 
511  /**
512  * Disable PMT.
513  */
514  void disable()
515  {
516  status = false;
517 
518  QE .fix(0.0);
519  TTS.fix(TTS);
520  t0 .fix(0.0);
521  bg .fix(0.0);
522  }
523 
524 
525  /**
526  * Enable PMT.
527  */
528  void enable()
529  {
530  status = true;
531 
532  QE .set();
533  TTS.set();
534  t0 .set();
535  bg .set();
536  }
537 
538 
539  /**
540  * Write PMT parameters to output stream.
541  *
542  * \param out output stream
543  * \param object PMT parameters
544  * \return output stream
545  */
546  friend inline std::ostream& operator<<(std::ostream& out, const JPMTParameters_t& object)
547  {
548  using namespace std;
549 
550  out << "QE " << FIXED(7,3) << object.QE << endl;
551  out << "TTS " << FIXED(7,3) << object.TTS << endl;
552  out << "t0 " << FIXED(7,3) << object.t0 << endl;
553  out << "bg " << FIXED(7,3) << object.bg << endl;
554 
555  return out;
556  }
557 
558 
559  bool status; //!< status
560  JParameter_t QE; //!< relative quantum efficiency [unit]
561  JParameter_t TTS; //!< transition-time spread [ns]
562  JParameter_t t0; //!< time offset [ns]
563  JParameter_t bg; //!< background [Hz/ns]
564  };
565 
566 
567  /**
568  * Fit parameters for two-fold coincidence rate due to K40.
569  */
571  /**
572  * Default constructor.
573  */
575  {
576  reset();
577  }
578 
579 
580  /**
581  * Reset.
582  */
583  void reset()
584  {
585  R .set(0.0);
586  p1.set(0.0);
587  p2.set(0.0);
588  p3.set(0.0);
589  p4.set(0.0);
590  cc.set(0.0);
591  }
592 
593 
594  JParameter_t R; //!< maximal coincidence rate [Hz]
595  JParameter_t p1; //!< 1st order angle dependence coincidence rate
596  JParameter_t p2; //!< 2nd order angle dependence coincidence rate
597  JParameter_t p3; //!< 3rd order angle dependence coincidence rate
598  JParameter_t p4; //!< 4th order angle dependence coincidence rate
599  JParameter_t cc; //!< fraction of signal correlated background
600  };
601 
602 
603  /**
604  * Fit parameters for two-fold coincidence rate due to K40.
605  */
606  struct JK40Parameters :
608  {
609  /**
610  * Default constructor.
611  */
613  {}
614 
615 
616  /**
617  * Get default values.
618  *
619  * The default parameter values are set to those obtained from a designated simulation
620  * of K40 decays (see http://wiki.km3net.de/index.php/OMGsim_simulations_for_K40_fit).\n
621  * If you change these values, you may also want to change the corresponding values in JK40DefaultSimulator.hh.
622  *
623  * \return parameters
624  */
625  static const JK40Parameters& getInstance()
626  {
627  static JK40Parameters parameters;
628 
629  parameters.R .set(18.460546);
630  parameters.p1.set( 3.0767);
631  parameters.p2.set(-1.2078);
632  parameters.p3.set( 0.9905);
633  parameters.p4.set( 0.9379);
634  parameters.cc.set( 0.0);
635 
636  return parameters;
637  }
638 
639 
640  /**
641  * Get K40 parameters.
642  *
643  * \return K40 parameters
644  */
646  {
647  return static_cast<const JK40Parameters&>(*this);
648  }
649 
650 
651  /**
652  * Set K40 parameters.
653  *
654  * \param parameters K40 parameters
655  */
657  {
658  static_cast<JK40Parameters_t&>(*this) = parameters;
659  }
660 
661 
662  /**
663  * Get number of fit parameters.
664  *
665  * \return number of parameters
666  */
667  inline size_t getN() const
668  {
669  return ((R .isFree() ? 1 : 0) +
670  (p1.isFree() ? 1 : 0) +
671  (p2.isFree() ? 1 : 0) +
672  (p3.isFree() ? 1 : 0) +
673  (p4.isFree() ? 1 : 0) +
674  (cc.isFree() ? 1 : 0));
675  }
676 
677 
678  /**
679  * Get index of parameter.
680  *
681  * \param p pointer to data member
682  * \return index
683  */
685  {
686  if (!(this->*p).isFree()) {
687  return INVALID_INDEX;
688  }
689 
690  int N = 0;
691 
692  if (p == &JK40Parameters::R) { return N; } if (R .isFree()) { ++N; }
693  if (p == &JK40Parameters::p1) { return N; } if (p1.isFree()) { ++N; }
694  if (p == &JK40Parameters::p2) { return N; } if (p2.isFree()) { ++N; }
695  if (p == &JK40Parameters::p3) { return N; } if (p3.isFree()) { ++N; }
696  if (p == &JK40Parameters::p4) { return N; } if (p4.isFree()) { ++N; }
697  if (p == &JK40Parameters::cc) { return N; } if (cc.isFree()) { ++N; }
698 
699  return INVALID_INDEX;
700  }
701 
702 
703  /**
704  * Get K40 coincidence rate as a function of cosine angle between PMT axes.
705  *
706  * \param ct cosine angle between PMT axes
707  * \return rate [Hz]
708  */
709  double getValue(const double ct) const
710  {
711  return R * exp(ct*(p1+ct*(p2+ct*(p3+ct*p4))) - (p1+p2+p3+p4));
712  }
713 
714 
715  /**
716  * Get gradient.
717  *
718  * \param ct cosine angle between PMT axes
719  * \return gradient
720  */
721  const JK40Parameters_t& getGradient(const double ct) const
722  {
723  gradient.reset();
724 
725  const double rate = getValue(ct);
726  const double ct2 = ct * ct;
727 
728  if (R .isFree()) { gradient.R = rate / R; }
729  if (p1.isFree()) { gradient.p1 = rate * ct - rate; }
730  if (p2.isFree()) { gradient.p2 = rate * ct2 - rate; }
731  if (p3.isFree()) { gradient.p3 = rate * ct2 * ct - rate; }
732  if (p4.isFree()) { gradient.p4 = rate * ct2 * ct2 - rate; }
733  if (cc.isFree()) { gradient.cc = rate; }
734 
735  return gradient;
736  }
737 
738  private:
740  };
741 
742 
743  /**
744  * Fit model.
745  *
746  * In the absence of TDC constraints, the average time offset is fixed to zero.
747  */
748  struct JModel :
749  public JK40Parameters,
750  public JModule,
751  public JCombinatorics_t
752  {
756 
757 
758  /**
759  * Auxiliary data structure for derived quantities of a given PMT pair.
760  */
761  struct real_type {
762  double ct; //!< cosine angle between PMT axes
763  double t0; //!< time offset [ns]
764  double sigma; //!< total width [ns]
765  double signal; //!< combined signal
766  double background; //!< combined background
767  };
768 
769 
770  /**
771  * Default constructor.
772  */
774  {}
775 
776 
777  /**
778  * Constructor.
779  *
780  * \param module detector module
781  * \param parameters K40 parameters
782  * \param TDC TDC constraints
783  * \param option option
784  */
785  JModel(const JModule& module,
786  const JK40Parameters& parameters,
787  const JTDC_t::range_type& TDC,
788  const int option) :
789  JModule (module),
790  JCombinatorics_t(module)
791  {
792  setK40Parameters(parameters);
793 
794  for (int i = 0; i != NUMBER_OF_PMTS; ++i) {
795  this->parameters[i] = JPMTParameters_t::getInstance();
796  }
797 
798  for (JTDC_t::const_iterator i = TDC.first; i != TDC.second; ++i) {
799  this->parameters[i->second].t0.fix();
800  }
801 
802  this->index = (TDC.first == TDC.second ? 0 : INVALID_INDEX);
803 
804  setOption(option);
805  }
806 
807 
808  /**
809  * Constructor.
810  *
811  * \param module detector module
812  * \param parameters K40 parameters
813  */
814  JModel(const JModule& module,
815  const JK40Parameters& parameters) :
816  JModule (module),
817  JCombinatorics_t(module)
818  {
819  setK40Parameters(parameters);
820 
821  for (int i = 0; i != NUMBER_OF_PMTS; ++i) {
822  this->parameters[i] = JPMTParameters_t::getInstance();
823  }
824  }
825 
826 
827  /**
828  * Get fit option.
829  *
830  * \return option
831  */
833  {
834  return option;
835  }
836 
837 
838  /**
839  * Set fit option.
840  *
841  * \param option option
842  */
843  inline void setOption(const int option)
844  {
845  switch (option) {
846 
847  case FIT_PMTS_t:
848 
849  R .fix();
850  p1.fix();
851  p2.fix();
852  p3.fix();
853  p4.fix();
854 
855  for (int i = 0; i != NUMBER_OF_PMTS; ++i) {
856  parameters[i].bg.fix();
857  }
858 
859  parameters[index].t0.fix();
860 
861  break;
862 
864 
865  R .fix();
866 
867  for (int i = 0; i != NUMBER_OF_PMTS; ++i) {
868  parameters[i].bg.fix();
869  }
870 
871  parameters[index].t0.fix();
872 
873  break;
874 
876 
877  R .fix();
878  p1.fix();
879  p2.fix();
880  p3.fix();
881  p4.fix();
882 
883  parameters[index].t0.fix();
884 
885  break;
886 
887  case FIT_PMTS_QE_FIXED_t:
888 
889  R .fix();
890  p1.fix();
891  p2.fix();
892  p3.fix();
893  p4.fix();
894 
895  for (int i = 0; i != NUMBER_OF_PMTS; ++i) {
896  parameters[i].QE.fix();
897  parameters[i].bg.fix();
898  }
899 
900  parameters[index].t0.fix();
901 
902  break;
903 
904  case FIT_MODEL_t:
905 
906  cc.fix();
907 
908  for (int i = 0; i != NUMBER_OF_PMTS; ++i) {
909  parameters[i].QE.fix();
910  parameters[i].t0.fix();
911  parameters[i].bg.fix();
912  }
913 
914  break;
915 
916  default:
917 
918  THROW(JValueOutOfRange, "Invalid option " << option);
919  }
920 
921  this->option = static_cast<JOption_t>(option);
922  }
923 
924 
925  /**
926  * Check if time offset is fixed.
927  *
928  * \return true if time offset fixed; else false
929  */
930  bool hasFixedTimeOffset() const
931  {
932  return index != INVALID_INDEX;
933  }
934 
935 
936  /**
937  * Get time offset.
938  *
939  * \return time offset
940  */
941  double getFixedTimeOffset() const
942  {
943  double t0 = 0.0;
944  size_t N = 0;
945 
946  for (int i = 0; i != NUMBER_OF_PMTS; ++i) {
947  if (parameters[i].t0.isFree()) {
948  t0 += parameters[i].t0;
949  N += 1;
950  }
951  }
952 
953  return t0 /= N;
954  }
955 
956 
957  /**
958  * Get index of PMT used for fixed time offset
959  *
960  * \return index
961  */
962  int getIndex() const
963  {
964  return index;
965  }
966 
967 
968  /**
969  * Set index of PMT used for fixed time offset
970  */
971  void setIndex()
972  {
973  if (index != INVALID_INDEX) {
974 
975  for (int i = 0; i != NUMBER_OF_PMTS; ++i) {
976 
977  if (parameters[i].status) {
978 
979  index = i;
980 
981  parameters[index].t0.fix();
982 
983  return;
984  }
985  }
986 
987  THROW(JValueOutOfRange, "No valid index.");
988  }
989  }
990 
991 
992  /**
993  * Get number of fit parameters.
994  *
995  * \return number of parameters
996  */
997  inline size_t getN() const
998  {
999  size_t N = JK40Parameters::getN();
1000 
1001  for (int i = 0; i != NUMBER_OF_PMTS; ++i) {
1002  N += parameters[i].getN();
1003  }
1004 
1005  return N;
1006  }
1007 
1008 
1009  /**
1010  * Get index of parameter.
1011  *
1012  * \param pmt PMT
1013  * \return index
1014  */
1015  int getIndex(int pmt) const
1016  {
1017  int N = JK40Parameters::getN();
1018 
1019  for (int i = 0; i != pmt; ++i) {
1020  N += parameters[i].getN();
1021  }
1022 
1023  return N;
1024  }
1025 
1026 
1027  /**
1028  * Get index of parameter.
1029  *
1030  * \param pmt PMT
1031  * \param p pointer to data member
1032  * \return index
1033  */
1034  int getIndex(int pmt, JParameter_t JPMTParameters_t::*p) const
1035  {
1036  return getIndex(pmt) + parameters[pmt].getIndex(p);
1037  }
1038 
1039 
1040  /**
1041  * Get intrinsic K40 arrival time spread.
1042  *
1043  * \return sigma [ns]
1044  */
1045  double getSigmaK40() const
1046  {
1047  return this->sigmaK40_ns;
1048  }
1049 
1050 
1051  /**
1052  * Set intrinsic K40 arrival time spread.
1053  *
1054  * \param sigma sigma [ns]
1055  */
1056  void setSigmaK40(const double sigma)
1057  {
1058  this->sigmaK40_ns = sigma;
1059  }
1060 
1061 
1062  /**
1063  * Get derived parameters.
1064  *
1065  * \param pair PMT pair
1066  * \return parameters
1067  */
1068  const real_type& getReal(const pair_type& pair) const
1069  {
1070  real.ct = JPP::getDot((*this)[pair.first].getDirection(), (*this)[pair.second].getDirection());
1071 
1072  real.t0 = (pair.first == this->index ? -this->parameters[pair.second].t0() :
1073  pair.second == this->index ? +this->parameters[pair.first ].t0() :
1074  this->parameters[pair.first].t0() - this->parameters[pair.second].t0());
1075 
1076  real.sigma = sqrt(this->parameters[pair.first ].TTS() * this->parameters[pair.first ].TTS() +
1077  this->parameters[pair.second].TTS() * this->parameters[pair.second].TTS() +
1078  this->getSigmaK40() * this->getSigmaK40());
1079 
1080  real.signal = this->parameters[pair.first].QE() * this->parameters[pair.second].QE();
1081 
1082  real.background = this->parameters[pair.first].bg() + this->parameters[pair.second].bg();
1083 
1084  return real;
1085  }
1086 
1087 
1088  /**
1089  * Get K40 coincidence rate.
1090  *
1091  * \param pair PMT pair
1092  * \param dt_ns time difference [ns]
1093  * \return rate [Hz/ns]
1094  */
1095  double getValue(const pair_type& pair, const double dt_ns) const
1096  {
1097  using namespace std;
1098  using namespace JPP;
1099 
1100  const real_type& real = getReal(pair);
1101 
1102  const JGauss gauss(real.t0, real.sigma, real.signal);
1103 
1104  const double R1 = this->getValue(real.ct);
1105  const double R2 = gauss.getValue(dt_ns);
1106 
1107  return real.background + R1 * (cc() + R2);
1108  }
1109 
1110 
1111  /**
1112  * Get K40 coincidence rate.
1113  *
1114  * \param pair PMT pair
1115  * \return rate [Hz]
1116  */
1117  double getValue(const pair_type& pair) const
1118  {
1119  using namespace std;
1120  using namespace JPP;
1121 
1122  const real_type& real = getReal(pair);
1123 
1124  const double R1 = this->getValue(real.ct);
1125  const double R2 = real.signal;
1126 
1127  return real.background + R1 * (cc() + R2);
1128  }
1129 
1130 
1131  /**
1132  * Write model parameters to output stream.
1133  *
1134  * \param out output stream
1135  * \param object model parameters
1136  * \return output stream
1137  */
1138  friend inline std::ostream& operator<<(std::ostream& out, const JModel& object)
1139  {
1140  using namespace std;
1141 
1142  out << "Module " << setw(10) << object.getID() << endl;
1143  out << "option " << object.option << endl;
1144  out << "index " << object.index << endl;
1145  out << "Rate [Hz] " << FIXED(12,6) << object.R << endl;
1146  out << "p1 " << FIXED(12,6) << object.p1 << endl;
1147  out << "p2 " << FIXED(12,6) << object.p2 << endl;
1148  out << "p3 " << FIXED(12,6) << object.p3 << endl;
1149  out << "p4 " << FIXED(12,6) << object.p4 << endl;
1150  out << "cc " << FIXED(12,6) << object.cc << endl;
1151 
1152  for (int i = 0; i != NUMBER_OF_PMTS; ++i) {
1153  out << "PMT[" << FILL(2,'0') << i << FILL() << "]." << object.parameters[i].status << endl << object.parameters[i];
1154  }
1155 
1156  return out;
1157  }
1158 
1160 
1161  private:
1162  int index; //!< index of PMT used for fixed time offset
1163  double sigmaK40_ns = 0.54; //!< intrinsic K40 arrival time spread [ns]
1165  mutable real_type real;
1166  };
1167 
1168 
1169  /**
1170  * Fit.
1171  */
1172  class JFit
1173  {
1174  public:
1175  /**
1176  * Result type.
1177  */
1178  struct result_type {
1179  double chi2;
1180  int ndf;
1181  };
1182 
1183  typedef std::shared_ptr<JMEstimator> estimator_type;
1184 
1185 
1186  /**
1187  * Constructor
1188  *
1189  * \param option M-estimator
1190  * \param debug debug
1191  */
1192  JFit(const int option, const int debug) :
1193  debug(debug)
1194  {
1195  using namespace JPP;
1196 
1197  estimator.reset(getMEstimator(option));
1198  }
1199 
1200 
1201  /**
1202  * Fit.
1203  *
1204  * \param data data
1205  * \return chi2, NDF
1206  */
1208  {
1209  using namespace std;
1210  using namespace JPP;
1211 
1212 
1213  value.setIndex();
1214 
1215  const size_t N = value.getN();
1216 
1217  V.resize(N);
1218  Y.resize(N);
1219  h.resize(N);
1220 
1221  int ndf = 0;
1222 
1223  for (data_type::const_iterator ix = data.begin(); ix != data.end(); ++ix) {
1224 
1225  const pair_type& pair = ix->first;
1226 
1227  if (value.parameters[pair.first ].status &&
1228  value.parameters[pair.second].status) {
1229 
1230  ndf += ix->second.size();
1231  }
1232  }
1233 
1234  ndf -= value.getN();
1235 
1236 
1237  lambda = LAMBDA_MIN;
1238 
1239  double precessor = numeric_limits<double>::max();
1240 
1241  for (numberOfIterations = 0; numberOfIterations != MAXIMUM_ITERATIONS; ++numberOfIterations) {
1242 
1243  DEBUG("step: " << numberOfIterations << endl);
1244 
1245  evaluate(data);
1246 
1247  DEBUG("lambda: " << FIXED(12,5) << lambda << endl);
1248  DEBUG("chi2: " << FIXED(12,3) << successor << endl);
1249 
1250  if (successor < precessor) {
1251 
1252  if (numberOfIterations != 0) {
1253 
1254  if (fabs(precessor - successor) < EPSILON*fabs(precessor)) {
1255  return { successor / estimator->getRho(1.0), ndf };
1256  }
1257 
1258  if (lambda > LAMBDA_MIN) {
1259  lambda /= LAMBDA_DOWN;
1260  }
1261  }
1262 
1263  precessor = successor;
1264  previous = value;
1265 
1266  } else {
1267 
1268  value = previous;
1269  lambda *= LAMBDA_UP;
1270 
1271  if (lambda > LAMBDA_MAX) {
1272  return { precessor / estimator->getRho(1.0), ndf }; // no improvement found
1273  }
1274 
1275  evaluate(data);
1276  }
1277 
1278  if (debug >= debug_t) {
1279 
1280  size_t row = 0;
1281 
1282  if (value.R .isFree()) { cout << "R " << FIXED(12,5) << Y[row] << endl; ++row; }
1283  if (value.p1.isFree()) { cout << "p1 " << FIXED(12,5) << Y[row] << endl; ++row; }
1284  if (value.p2.isFree()) { cout << "p2 " << FIXED(12,5) << Y[row] << endl; ++row; }
1285  if (value.p3.isFree()) { cout << "p3 " << FIXED(12,5) << Y[row] << endl; ++row; }
1286  if (value.p4.isFree()) { cout << "p4 " << FIXED(12,5) << Y[row] << endl; ++row; }
1287  if (value.cc.isFree()) { cout << "cc " << FIXED(12,3) << Y[row] << endl; ++row; }
1288 
1289  for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
1290  if (value.parameters[pmt].QE .isFree()) { cout << "PMT[" << setw(2) << pmt << "].QE " << FIXED(12,5) << Y[row] << endl; ++row; }
1291  if (value.parameters[pmt].TTS.isFree()) { cout << "PMT[" << setw(2) << pmt << "].TTS " << FIXED(12,5) << Y[row] << endl; ++row; }
1292  if (value.parameters[pmt].t0 .isFree()) { cout << "PMT[" << setw(2) << pmt << "].t0 " << FIXED(12,5) << Y[row] << endl; ++row; }
1293  if (value.parameters[pmt].bg .isFree()) { cout << "PMT[" << setw(2) << pmt << "].bg " << FIXED(12,5) << Y[row] << endl; ++row; }
1294  }
1295  }
1296 
1297  // force definite positiveness
1298 
1299  for (size_t i = 0; i != N; ++i) {
1300 
1301  if (V(i,i) < PIVOT) {
1302  V(i,i) = PIVOT;
1303  }
1304 
1305  h[i] = 1.0 / sqrt(V(i,i));
1306  }
1307 
1308  // normalisation
1309 
1310  for (size_t i = 0; i != N; ++i) {
1311  for (size_t j = 0; j != i; ++j) {
1312  V(j,i) *= h[i] * h[j];
1313  V(i,j) = V(j,i);
1314  }
1315  }
1316 
1317  for (size_t i = 0; i != N; ++i) {
1318  V(i,i) = 1.0 + lambda;
1319  }
1320 
1321  // solve A x = b
1322 
1323  for (size_t col = 0; col != N; ++col) {
1324  Y[col] *= h[col];
1325  }
1326 
1327  try {
1328  V.solve(Y);
1329  }
1330  catch (const exception& error) {
1331 
1332  ERROR("JGandalf: " << error.what() << endl << V << endl);
1333 
1334  break;
1335  }
1336 
1337  // update value
1338 
1339  const double factor = 2.0;
1340 
1341  size_t row = 0;
1342 
1343  if (value.R .isFree()) { value.R -= factor * h[row] * Y[row]; ++row; }
1344  if (value.p1.isFree()) { value.p1 -= factor * h[row] * Y[row]; ++row; }
1345  if (value.p2.isFree()) { value.p2 -= factor * h[row] * Y[row]; ++row; }
1346  if (value.p3.isFree()) { value.p3 -= factor * h[row] * Y[row]; ++row; }
1347  if (value.p4.isFree()) { value.p4 -= factor * h[row] * Y[row]; ++row; }
1348  if (value.cc.isFree()) { value.cc -= factor * h[row] * Y[row]; ++row; }
1349 
1350  for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
1351  if (value.parameters[pmt].QE .isFree()) { value.parameters[pmt].QE -= factor * h[row] * Y[row]; ++row; }
1352  if (value.parameters[pmt].TTS.isFree()) { value.parameters[pmt].TTS -= factor * h[row] * Y[row]; ++row; }
1353  if (value.parameters[pmt].t0 .isFree()) { value.parameters[pmt].t0 -= factor * h[row] * Y[row]; ++row; }
1354  if (value.parameters[pmt].bg .isFree()) { value.parameters[pmt].bg -= factor * h[row] * Y[row]; ++row; }
1355  }
1356  }
1357 
1358  return { precessor / estimator->getRho(1.0), ndf };
1359  }
1360 
1361 
1362  static constexpr int MAXIMUM_ITERATIONS = 10000; //!< maximal number of iterations.
1363  static constexpr double EPSILON = 1.0e-4; //!< maximal distance to minimum.
1364  static constexpr double LAMBDA_MIN = 0.01; //!< minimal value control parameter
1365  static constexpr double LAMBDA_MAX = 100.0; //!< maximal value control parameter
1366  static constexpr double LAMBDA_UP = 10.0; //!< multiplication factor control parameter
1367  static constexpr double LAMBDA_DOWN = 10.0; //!< multiplication factor control parameter
1368  static constexpr double PIVOT = std::numeric_limits<double>::epsilon(); //!< minimal value diagonal element of matrix
1369 
1370  int debug;
1371  estimator_type estimator; //!< M-Estimator function
1372 
1373  double lambda;
1377 
1378  private:
1379  /**
1380  * Evaluation of fit.
1381  *
1382  * \param data data
1383  */
1384  void evaluate(const data_type& data)
1385  {
1386  using namespace std;
1387  using namespace JPP;
1388 
1389  typedef JModel::real_type real_type;
1390 
1391 
1392  successor = 0.0;
1393 
1394  V.reset();
1395  Y.reset();
1396 
1397 
1398  // model parameter indices
1399 
1400  const struct M_t {
1401  M_t(const JModel& model)
1402  {
1403  R = model.getIndex(&JK40Parameters_t::R);
1404  p1 = model.getIndex(&JK40Parameters_t::p1);
1405  p2 = model.getIndex(&JK40Parameters_t::p2);
1406  p3 = model.getIndex(&JK40Parameters_t::p3);
1407  p4 = model.getIndex(&JK40Parameters_t::p4);
1408  cc = model.getIndex(&JK40Parameters_t::cc);
1409  }
1410 
1411  int R;
1412  int p1;
1413  int p2;
1414  int p3;
1415  int p4;
1416  int cc;
1417 
1418  } M(value);
1419 
1420 
1421  // PMT parameter indices
1422 
1423  struct I_t {
1424  I_t(const JModel& model, const int pmt) :
1425  QE (INVALID_INDEX),
1426  TTS(INVALID_INDEX),
1427  t0 (INVALID_INDEX),
1428  bg (INVALID_INDEX)
1429  {
1430  const int index = model.getIndex(pmt);
1431 
1432  int N = 0;
1433 
1434  if (model.parameters[pmt].QE .isFree()) { QE = index + N; ++N; }
1435  if (model.parameters[pmt].TTS.isFree()) { TTS = index + N; ++N; }
1436  if (model.parameters[pmt].t0 .isFree()) { t0 = index + N; ++N; }
1437  if (model.parameters[pmt].bg .isFree()) { bg = index + N; ++N; }
1438  }
1439 
1440  int QE;
1441  int TTS;
1442  int t0;
1443  int bg;
1444  };
1445 
1446 
1448 
1449  buffer_type buffer;
1450 
1451  for (data_type::const_iterator ix = data.begin(); ix != data.end(); ++ix) {
1452 
1453  const pair_type& pair = ix->first;
1454 
1455  if (value.parameters[pair.first ].status &&
1456  value.parameters[pair.second].status) {
1457 
1458  const real_type& real = value.getReal(pair);
1459 
1460  const JGauss gauss(real.t0, real.sigma, real.signal);
1461 
1462  const double R1 = value.getValue (real.ct);
1463  const JK40Parameters_t& R1p = value.getGradient(real.ct);
1464 
1465  const std::pair<I_t, I_t> PMT(I_t(value, pair.first),
1466  I_t(value, pair.second));
1467 
1468  for (const rate_type& iy : ix->second) {
1469 
1470  const double R2 = gauss.getValue (iy.dt_ns);
1471  const JGauss& R2p = gauss.getGradient(iy.dt_ns);
1472 
1473  const double R = real.background + R1 * (value.cc() + R2);
1474  const double u = (iy.value - R) / iy.error;
1475  const double w = -estimator->getPsi(u) / iy.error;
1476 
1477  successor += estimator->getRho(u);
1478 
1479  buffer.clear();
1480 
1481  if (M.R != INVALID_INDEX) { buffer.push_back({M.R, w * (value.cc() + R2) * R1p.R () * value.R .getDerivative()}); }
1482  if (M.p1 != INVALID_INDEX) { buffer.push_back({M.p1, w * (value.cc() + R2) * R1p.p1() * value.p1.getDerivative()}); }
1483  if (M.p2 != INVALID_INDEX) { buffer.push_back({M.p2, w * (value.cc() + R2) * R1p.p2() * value.p2.getDerivative()}); }
1484  if (M.p3 != INVALID_INDEX) { buffer.push_back({M.p3, w * (value.cc() + R2) * R1p.p3() * value.p3.getDerivative()}); }
1485  if (M.p4 != INVALID_INDEX) { buffer.push_back({M.p4, w * (value.cc() + R2) * R1p.p4() * value.p4.getDerivative()}); }
1486  if (M.cc != INVALID_INDEX) { buffer.push_back({M.cc, w * R1 * R1p.cc() * value.cc.getDerivative()}); }
1487 
1488  if (PMT.first .QE != INVALID_INDEX) { buffer.push_back({PMT.first .QE , w * R1 * R2p.signal * value.parameters[pair.second].QE () * value.parameters[pair.first ].QE .getDerivative()}); }
1489  if (PMT.second.QE != INVALID_INDEX) { buffer.push_back({PMT.second.QE , w * R1 * R2p.signal * value.parameters[pair.first ].QE () * value.parameters[pair.second].QE .getDerivative()}); }
1490  if (PMT.first .TTS != INVALID_INDEX) { buffer.push_back({PMT.first .TTS, w * R1 * R2p.sigma * value.parameters[pair.first ].TTS() * value.parameters[pair.first ].TTS.getDerivative() / real.sigma}); }
1491  if (PMT.second.TTS != INVALID_INDEX) { buffer.push_back({PMT.second.TTS, w * R1 * R2p.sigma * value.parameters[pair.second].TTS() * value.parameters[pair.second].TTS.getDerivative() / real.sigma}); }
1492  if (PMT.first .t0 != INVALID_INDEX) { buffer.push_back({PMT.first .t0, w * R1 * R2p.mean * value.parameters[pair.first ].t0 .getDerivative() * +1.0}); }
1493  if (PMT.second.t0 != INVALID_INDEX) { buffer.push_back({PMT.second.t0, w * R1 * R2p.mean * value.parameters[pair.second].t0 .getDerivative() * -1.0}); }
1494  if (PMT.first .bg != INVALID_INDEX) { buffer.push_back({PMT.first .bg, w * value.parameters[pair.first ].bg .getDerivative()}); }
1495  if (PMT.second.bg != INVALID_INDEX) { buffer.push_back({PMT.second.bg, w * value.parameters[pair.second].bg .getDerivative()}); }
1496 
1497  for (buffer_type::const_iterator row = buffer.begin(); row != buffer.end(); ++row) {
1498 
1499  Y[row->first] += row->second;
1500 
1501  V[row->first][row->first] += row->second * row->second;
1502 
1503  for (buffer_type::const_iterator col = buffer.begin(); col != row; ++col) {
1504  V[row->first][col->first] += row->second * col->second;
1505  V[col->first][row->first] = V[row->first][col->first];
1506  }
1507  }
1508  }
1509  }
1510  }
1511  }
1512 
1513  JMATH::JVectorND Y; // gradient
1514  double successor;
1516  std::vector<double> h; // normalisation vector
1517  };
1518 }
1519 
1520 #endif
1521 
1522 
Data structure for measured coincidence rate of pair of PMTs.
Definition: JFitK40.hh:64
std::vector< double > h
Definition: JFitK40.hh:1516
result_type operator()(const data_type &data)
Fit.
Definition: JFitK40.hh:1207
JOption_t
Fit options.
Definition: JFitK40.hh:50
JTOOLS::JRange< double > range_type
Type definition for range of parameter values.
Definition: JFitK40.hh:123
size_t getN() const
Get number of fit parameters.
Definition: JFitK40.hh:479
data_type w[N+1][M+1]
Definition: JPolint.hh:867
Exceptions.
void set()
Set current value.
Definition: JFitK40.hh:269
double getValue(const JScale_t scale)
Get numerical value corresponding to scale.
Definition: JScale.hh:47
int getIndex(JParameter_t JPMTParameters_t::*p) const
Get index of parameter.
Definition: JFitK40.hh:494
const real_type & getReal(const pair_type &pair) const
Get derived parameters.
Definition: JFitK40.hh:1068
debug
Definition: JMessage.hh:29
T getLowerLimit() const
Get lower limit.
Definition: JRange.hh:202
int getIndex(JParameter_t JK40Parameters::*p) const
Get index of parameter.
Definition: JFitK40.hh:684
JParameter_t t0
time offset [ns]
Definition: JFitK40.hh:562
Auxiliary base class for aritmetic operations of derived class types.
Definition: JMath.hh:109
JK40Parameters()
Default constructor.
Definition: JFitK40.hh:612
TPaveText * p1
Auxiliary methods for geometrical methods.
Data structure for a composite optical module.
Definition: JModule.hh:68
void setSigmaK40(const double sigma)
Set intrinsic K40 arrival time spread.
Definition: JFitK40.hh:1056
JParameter_t R
maximal coincidence rate [Hz]
Definition: JFitK40.hh:594
bool hasFixedTimeOffset() const
Check if time offset is fixed.
Definition: JFitK40.hh:930
static const int INVALID_INDEX
invalid index
Definition: JFitK40.hh:58
size_t getN() const
Get number of fit parameters.
Definition: JFitK40.hh:667
Interface for maximum likelihood estimator (M-estimator).
Definition: JMEstimator.hh:20
JParameter_t & add(const JParameter_t &parameter)
Add parameter.
Definition: JFitK40.hh:168
JParameter_t & sub(const JParameter_t &parameter)
Subtract parameter.
Definition: JFitK40.hh:182
JK40Parameters_t()
Default constructor.
Definition: JFitK40.hh:574
JParameter_t & negate()
Negate parameter.
Definition: JFitK40.hh:154
double getDot(const JNeutrinoDirection &first, const JNeutrinoDirection &second)
Dot product.
Definition: JAstronomy.hh:674
const JK40Parameters_t & getGradient(const double ct) const
Get gradient.
Definition: JFitK40.hh:721
static const JPBS_t PMT(3, 4, 2, 3)
PBS of photo-multiplier tube (PMT)
void setIndex()
Set index of PMT used for fixed time offset.
Definition: JFitK40.hh:971
void setK40Parameters(const JK40Parameters_t &parameters)
Set K40 parameters.
Definition: JFitK40.hh:656
int index
index of PMT used for fixed time offset
Definition: JFitK40.hh:1162
void set(const double value)
Set value.
Definition: JFitK40.hh:303
JOption_t option
Definition: JFitK40.hh:1164
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Definition: JException.hh:712
then set_variable singlesRate set_variable doublesRate set_variable numberOfSlices echo Generating random background echo Singles rate
double getDerivative() const
Get derivative of value.
Definition: JFitK40.hh:332
double getValue(const pair_type &pair) const
Get K40 coincidence rate.
Definition: JFitK40.hh:1117
Data structure for a pair of indices.
#define R1(x)
*fatal Wrong number of arguments esac JCookie sh typeset Z DETECTOR typeset Z SOURCE_RUN typeset Z TARGET_RUN set_variable PARAMETERS_FILE $WORKDIR parameters
Definition: diff-Tuna.sh:38
int getIndex() const
Get index of PMT used for fixed time offset.
Definition: JFitK40.hh:962
JK40Parameters_t gradient
Definition: JFitK40.hh:739
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
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:446
V(JDAQEvent-JTriggerReprocessor)*1.0/(JDAQEvent+1.0e-10)
JModel()
Default constructor.
Definition: JFitK40.hh:773
double background
combined background
Definition: JFitK40.hh:766
JPMTParameters_t()
Default constructor.
Definition: JFitK40.hh:432
std::vector< JHitW0 > buffer_type
hits
Definition: JPerth.cc:67
then echo The file $DIR KM3NeT_00000001_00000000 root already please rename or remove it first
double getValue(const double ct) const
Get K40 coincidence rate as a function of cosine angle between PMT axes.
Definition: JFitK40.hh:709
friend std::ostream & operator<<(std::ostream &out, const JModel &object)
Write model parameters to output stream.
Definition: JFitK40.hh:1138
const double sigma[]
Definition: JQuadrature.cc:74
double dt_ns
time difference [ns]
Definition: JFitK40.hh:90
real_type real
Definition: JFitK40.hh:1165
T getLength() const
Get length (difference between upper and lower limit).
Definition: JRange.hh:289
then usage set_variable ACOUSTICS_WORKDIR $WORKDIR set_variable FORMULA sin([0]+2 *$PI *([1]+[2]*x)*x)" set_variable DY 1.0e-8 mkdir $WORKDIR for DETECTOR in $DETECTORS[*]
JFit(const int option, const int debug)
Constructor.
Definition: JFitK40.hh:1192
void setOption(const int option)
Set fit option.
Definition: JFitK40.hh:843
JParameter_t bg
background [Hz/ns]
Definition: JFitK40.hh:563
fit parameters of PMTs
Definition: JFitK40.hh:51
JMATH::JMatrixNS V
Definition: JFitK40.hh:1376
void fix()
Fix current value.
Definition: JFitK40.hh:278
JParameter_t & operator=(double value)
Assignment operator.
Definition: JFitK40.hh:369
Fit parameters for two-fold coincidence rate due to K40.
Definition: JFitK40.hh:606
JParameter_t QE
relative quantum efficiency [unit]
Definition: JFitK40.hh:560
JParameter_t p3
3rd order angle dependence coincidence rate
Definition: JFitK40.hh:597
int numberOfIterations
Definition: JFitK40.hh:1375
int getIndex(const int first, const int second) const
Get index of pair of indices.
void evaluate(const data_type &data)
Evaluation of fit.
Definition: JFitK40.hh:1384
double ct
cosine angle between PMT axes
Definition: JFitK40.hh:762
double error
error of rate [Hz/ns]
Definition: JFitK40.hh:92
fit parameters of PMTs and background
Definition: JFitK40.hh:53
estimator_type estimator
M-Estimator function.
Definition: JFitK40.hh:1371
JParameter_t()
Default constructor.
Definition: JFitK40.hh:129
rate_type(double dt_ns, double value, double error)
Constructor.
Definition: JFitK40.hh:82
fit parameters of PMTs and angular dependence of K40 rate
Definition: JFitK40.hh:52
#define ERROR(A)
Definition: JMessage.hh:66
double t0
time offset [ns]
Definition: JFitK40.hh:763
const JK40Parameters & getK40Parameters() const
Get K40 parameters.
Definition: JFitK40.hh:645
friend std::ostream & operator<<(std::ostream &out, const JParameter_t &object)
Write parameter to output stream.
Definition: JFitK40.hh:397
T constrain(argument_type x) const
Constrain value to range.
Definition: JRange.hh:350
int getIndex()
Get index for user I/O manipulation.
Definition: JManip.hh:26
void fix(const double value)
Fix value.
Definition: JFitK40.hh:319
double sigma
total width [ns]
Definition: JFitK40.hh:764
void enable()
Enable PMT.
Definition: JFitK40.hh:528
JParameter_t TTS
transition-time spread [ns]
Definition: JFitK40.hh:561
N x N symmetric matrix.
Definition: JMatrixNS.hh:28
static const JPMTParameters_t & getInstance()
Get default values.
Definition: JFitK40.hh:457
FIT_t
Fit options.
Definition: JFitK40.hh:114
Auxiliary data structure for sequence of same character.
Definition: JManip.hh:328
Fit parameters for single PMT.
Definition: JFitK40.hh:423
I/O manipulators.
JParameter_t p4
4th order angle dependence coincidence rate
Definition: JFitK40.hh:598
JParameter_t & mul(const JParameter_t &first, const JParameter_t &second)
Scale parameter.
Definition: JFitK40.hh:225
fit parameters of PMTs with QE fixed
Definition: JFitK40.hh:54
p2
Definition: module-Z:fit.sh:74
JOption_t getOption() const
Get fit option.
Definition: JFitK40.hh:832
JMATH::JVectorND Y
Definition: JFitK40.hh:1513
void reset(T &value)
Reset value.
then JCookie sh JDataQuality D $DETECTOR_ID R
Definition: JDataQuality.sh:41
bool is_valid() const
Check validity of range.
Definition: JRange.hh:311
size_t getN() const
Get number of fit parameters.
Definition: JFitK40.hh:997
friend std::istream & operator>>(std::istream &in, JParameter_t &object)
Read parameter from input stream.
Definition: JFitK40.hh:384
then usage $script< input file >[option[primary[working directory]]] nWhere option can be N
Definition: JMuonPostfit.sh:40
bool isFree() const
Check if parameter is free.
Definition: JFitK40.hh:238
Auxiliary class to define a range between two values.
JParameter_t p2
2nd order angle dependence coincidence rate
Definition: JFitK40.hh:596
Auxiliary data structure for derived quantities of a given PMT pair.
Definition: JFitK40.hh:761
fit parameters of K40 rate and TTSs of PMTs
Definition: JFitK40.hh:55
then set_variable NUMBER_OF_TESTS else set_variable NUMBER_OF_TESTS fi function gauss()
double getFixedTimeOffset() const
Get time offset.
Definition: JFitK40.hh:941
JParameter_t(const double value, const range_type &range=range_type::DEFAULT_RANGE())
Constructor.
Definition: JFitK40.hh:141
double signal
combined signal
Definition: JFitK40.hh:765
PMT combinatorics for optical module.
JParameter_t & mul(const double factor)
Scale parameter.
Definition: JFitK40.hh:196
double getValue(const pair_type &pair, const double dt_ns) const
Get K40 coincidence rate.
Definition: JFitK40.hh:1095
Data structure for measured coincidence rates of all pairs of PMTs in optical module.
Definition: JFitK40.hh:99
static const JK40Parameters & getInstance()
Get default values.
Definition: JFitK40.hh:625
set_variable TDC
Definition: JPrintTDC.sh:20
JModel(const JModule &module, const JK40Parameters &parameters, const JTDC_t::range_type &TDC, const int option)
Constructor.
Definition: JFitK40.hh:785
JModel(const JModule &module, const JK40Parameters &parameters)
Constructor.
Definition: JFitK40.hh:814
Base class for data structures with artithmetic capabilities.
int getIndex(int pmt) const
Get index of parameter.
Definition: JFitK40.hh:1015
double getSigmaK40() const
Get intrinsic K40 arrival time spread.
Definition: JFitK40.hh:1045
int j
Definition: JPolint.hh:792
Exception for accessing a value in a collection that is outside of its range.
Definition: JException.hh:178
int getIndex(int pmt, JParameter_t JPMTParameters_t::*p) const
Get index of parameter.
Definition: JFitK40.hh:1034
Fit parameters for two-fold coincidence rate due to K40.
Definition: JFitK40.hh:570
JMEstimator * getMEstimator(const int type)
Get M-Estimator.
Definition: JMEstimator.hh:203
static JRange< double, std::less< double > > DEFAULT_RANGE()
Default range.
Definition: JRange.hh:555
KM3NeT DAQ constants, bit handling, etc.
static const int NUMBER_OF_PMTS
Total number of PMTs in module.
Definition: JDAQ.hh:26
double u[N+1]
Definition: JPolint.hh:865
bool isBound() const
Check if parameter is bound.
Definition: JFitK40.hh:260
bool isFixed() const
Check if parameter is fixed.
Definition: JFitK40.hh:249
std::shared_ptr< JMEstimator > estimator_type
Definition: JFitK40.hh:1183
Fit model.
Definition: JFitK40.hh:748
then fatal Wrong number of arguments fi set_variable DETECTOR $argv[1] set_variable INPUT_FILE $argv[2] eval JPrintDetector a $DETECTOR O IDENTIFIER eval JPrintDetector a $DETECTOR O SUMMARY JAcoustics sh $DETECTOR_ID source JAcousticsToolkit sh CHECK_EXIT_CODE typeset A EMITTERS get_tripods $WORKDIR tripod txt EMITTERS get_transmitters $WORKDIR transmitter txt EMITTERS for EMITTER in
Definition: JCanberra.sh:46
rate_type()
Default constructor.
Definition: JFitK40.hh:68
double successor
Definition: JFitK40.hh:1514
void disable()
Disable PMT.
Definition: JFitK40.hh:514
JPMTParameters_t parameters[NUMBER_OF_PMTS]
Definition: JFitK40.hh:1159
JParameter_t & div(const double factor)
Scale parameter.
Definition: JFitK40.hh:210
then fatal Wrong number of arguments fi set_variable DETECTOR $argv[1] set_variable STRING $argv[2] set_array QUANTILES set_variable FORMULA *[0] exp(-0.5 *(x-[1])*(x-[1])/([2]*[2]))" set_variable MODULE `getModule -a $DETECTOR -L "$STRING 0"` source JAcousticsToolkit.sh typeset -A TRIPODS get_tripods $WORKDIR/tripod.txt TRIPODS XMEAN
double get() const
Get value.
Definition: JFitK40.hh:289
friend std::ostream & operator<<(std::ostream &out, const JPMTParameters_t &object)
Write PMT parameters to output stream.
Definition: JFitK40.hh:546
p3
Definition: module-Z:fit.sh:74
Maximum likelihood estimator (M-estimators).
double operator()() const
Type conversion operator.
Definition: JFitK40.hh:346
double value
value of rate [Hz/ns]
Definition: JFitK40.hh:91
const double epsilon
Definition: JQuadrature.cc:21
int debug
debug level
Nx1 matrix.
Definition: JVectorND.hh:21
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
JParameter_t p1
1st order angle dependence coincidence rate
Definition: JFitK40.hh:595
JParameter_t cc
fraction of signal correlated background
Definition: JFitK40.hh:599
Data structure for optical module.
Auxiliary class for fit parameter with optional limits.
Definition: JFitK40.hh:107