Jpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JMuonGandalf.hh
Go to the documentation of this file.
1 #ifndef JMUONGANDALF_INCLUDE
2 #define JMUONGANDALF_INCLUDE
3 
4 #include <string>
5 #include <iostream>
6 #include <iomanip>
7 #include <vector>
8 #include <algorithm>
9 #include <memory>
10 
11 #include "JDAQ/JDAQEvent.hh"
12 #include "JDAQ/JDAQTimeslice.hh"
13 #include "JDAQ/JDAQSummaryslice.hh"// <- not sure why this dependency is needed
14 
15 #include "JTrigger/JHit.hh"
16 #include "JTrigger/JTimeslice.hh"
17 #include "JTrigger/JHitL0.hh"
18 #include "JTrigger/JBuildL0.hh"
19 #include "JTrigger/JTriggerParameters.hh"// <- not sure why this dependency is needed
20 
21 #include "JFit/JLine1Z.hh"
22 #include "JFit/JLine1ZEstimator.hh"
23 #include "JFit/JFitToolkit.hh"
24 #include "JFit/JEvt.hh"
25 #include "JFit/JEvtToolkit.hh"
26 #include "JFit/JFitParameters.hh"
27 #include "JFit/JFitStatus.hh"
29 #include "JFit/JLine3Z.hh"
30 #include "JFit/JModel.hh"
31 #include "JFit/JSimplex.hh"
32 #include "JFit/JLine3ZRegressor.hh"
33 #include "JFit/JHitW0.hh"
34 
36 
37 /**
38  * \author mdejong, gmaggi
39  */
40 
41 namespace JFIT
42 {
43  /**
44  * class to handle Muon Gandalf reconstruction
45  * first, two other angular reconstructionS should be ran : JMuonPrefit and JMuonSimplex
46  */
48  {
49  private:
51 
54 
56 
57  public:
58 
59  /**
60  * Default constructor
61  */
63  {}
64 
65  /**
66  * Parameterized constructor
67  *
68  * \param router JSharedPointer of JModuleRouter, this is built via detector file.
69  * \param parameters struct that holds default-optimized parameters for the reconstruction, for ARCA and ORCA configuration. These default parameters could be found in $JPP_DATA.
70  * \param pdfFile PDF file descriptor
71  */
73  const JFIT::JMuonGandalfParameters_t &parameters,
74  std::string pdfFile):
75  router_(router),
76  parameters_(parameters)
77  {
78  JRegressor_t::debug = 1;//see what to do with this
79  JRegressor_t::T_ns.setRange(-50.0, +450.0);// ns
80  JRegressor_t::Vmax_npe = 10.0;
81  JRegressor_t::MAXIMUM_ITERATIONS = 10000;
82 
83  fitRegresor_ = {pdfFile, parameters.TTS_ns};
84 
85  fitRegresor_.parameters.resize(5);
86 
87  fitRegresor_.parameters[0] = JFIT::JLine3Z::pX();
88  fitRegresor_.parameters[1] = JFIT::JLine3Z::pY();
89  fitRegresor_.parameters[2] = JFIT::JLine3Z::pT();
90  fitRegresor_.parameters[3] = JFIT::JLine3Z::pDX();
91  fitRegresor_.parameters[4] = JFIT::JLine3Z::pDY();
92  }
93 
95 
96  /**
97  * Declaration of Member function that actually performs the reconstruction
98  *
99  * \param timeSlice which is contructed via a JDAQEvent
100  * \param InPreFits input JEvt which is a container of prefits, normally JMuonSimplex
101  * \param OutFits non const (output) JEvt.
102  */
103  void
104  getJEvt(const KM3NETDAQ::JDAQTimeslice &timeSlice,
105  JFIT::JEvt &InPreFits,
106  JFIT::JEvt &OutFits) const;
107 
108  /**
109  * Compare hits by PMT identifier and time.
110  *
111  * \param first first hit
112  * \param second second hit
113  * \return true if first before second; else false
114  */
115  static bool
116  compare(const JTRIGGER::JHitL0& first,
117  const JTRIGGER::JHitL0& second)
118  {
119  if (std::equal_to<KM3NETDAQ::JDAQPMTIdentifier>()(first, second))
120  return std::less<JTRIGGER::JHit>()(first, second);
121  else
122  return std::less<KM3NETDAQ::JDAQPMTIdentifier>()(first, second);
123  }
124 
125 
126  };
127 }
128 
129 /**
130  * Member function definition
131  */
132 void
134  JFIT::JEvt &InPreFits,
135  JFIT::JEvt &OutFits) const
136 {
137  if ( InPreFits.empty() ) return;
138 
139  const bool reprocess = parameters_.reprocess;
140  double roadWidth_m = parameters_.roadWidth_m;
141  const double R_Hz = parameters_.R_Hz;
142  const size_t numberOfPrefits = parameters_.numberOfPrefits;
143  const double E_GeV = parameters_.E_GeV;
144 
145  typedef std::vector<JTRIGGER::JHitL0> JDataL0_t;
146  typedef std::vector<JHitW0> JDataW0_t;
147 
149 
150  if (fitRegresor_.getRmax() < roadWidth_m) {
151  roadWidth_m = fitRegresor_.getRmax();
152  }
153 
154  JEvt::iterator __end = InPreFits.end();
155 
156  if (reprocess) {
157  __end = std::partition( InPreFits.begin(), __end,
159  );
160  }
161 
162  if (InPreFits.begin() != __end) {
163 
164  std::copy(InPreFits.begin(), __end, std::back_inserter(OutFits));
165 
166  if (numberOfPrefits > 0) {
167  std::advance(__end = InPreFits.begin(), std::min(numberOfPrefits, OutFits.size()));
168  }
169 
170  std::partial_sort(InPreFits.begin(), __end, InPreFits.end(), qualitySorter);
171 
172  __end = std::partition(InPreFits.begin(), __end,
173  JFIT::JHistory::is_event(InPreFits.begin()->getHistory())
174  );
175 
176  JDataL0_t dataL0;
177 
178  buildL0(timeSlice, *router_, std::back_inserter(dataL0));
179 
180  if (dataL0.size() >= fitRegresor_.parameters.size()) {
181 
182  for (JEvt::const_iterator track = InPreFits.begin(); track != __end; ++track) {
183 
185  const JFIT::JLine1Z tz(getPosition (*track).rotate(R), track->getT());
186  const JFIT::JModel<JFIT::JLine1Z> match(tz, roadWidth_m,
188 
189  //hit selection based on start value
190 
191  JDataW0_t data;
192 
193  for (JDataL0_t::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {
194 
195  JFIT::JHitW0 hit(*i, R_Hz);
196 
197  hit.rotate(R);
198 
199  if (match(hit)) {
200  data.push_back(hit);
201  }
202  }
203 
204  //select first hit
205 
206  std::sort(data.begin(), data.end(), compare);
207 
208  JDataW0_t::iterator __end = std::unique( data.begin(), data.end(),
209  std::equal_to<KM3NETDAQ::JDAQPMTIdentifier>()
210  );
211 
212  const int NDF = std::distance(data.begin(), __end) - fitRegresor_.parameters.size();
213 
214  if (NDF >= 0) {
215 
216  //set fit parameters
217 
218  if (track->getE() > 0.1)
219  fitRegresor_.E_GeV = track->getE();
220  else
221  fitRegresor_.E_GeV = E_GeV;
222 
223  const double chi2 = fitRegresor_(JLine3Z(tz), data.begin(), __end);
224 
226 
227  tb.rotate_back(R);
228 
229  const double energy(0);
230  OutFits.push_back(JFIT::getFit(JFIT::JHistory(track->getHistory()).add(JFIT::JFitApplication_t::JMUONGANDALF), tb,
231  JFIT::getQuality(chi2), NDF,energy,JFIT::JFitStatus_t::OKAY));
232 
233  //set additional weights
234 
235  OutFits.rbegin()->setW(JGANDALF_BETA0_RAD,sqrt(fitRegresor_.error.getDX() * fitRegresor_.error.getDX() +
236  fitRegresor_.error.getDY() * fitRegresor_.error.getDY()));
237  OutFits.rbegin()->setW(JGANDALF_BETA1_RAD,sqrt(fitRegresor_.error.getDX() * fitRegresor_.error.getDY()));
238  OutFits.rbegin()->setW(JGANDALF_CHI2, chi2);
239  OutFits.rbegin()->setW(JGANDALF_NUMBER_OF_HITS, distance(data.begin(), __end));
240  OutFits.rbegin()->setW(JGANDALF_LAMBDA, fitRegresor_.lambda);
241  OutFits.rbegin()->setW(JGANDALF_NUMBER_OF_ITERATIONS, fitRegresor_.numberOfIterations);
242  }
243  }
244  }
245  //apply default sorter
246 
247  std::sort(OutFits.begin(), OutFits.end(), qualitySorter);
248  }
249 
250 }
251 
252 #endif
253 
Auxiliary methods to evaluate Poisson probabilities and chi2.
Data regression method for JFIT::JLine3Z.
class to handle Muon Gandalf reconstruction first, two other angular reconstructionS should be ran : ...
Definition: JMuonGandalf.hh:47
JRegressor_t fitRegresor_
Definition: JMuonGandalf.hh:55
number of iterations from JGandalf.cc
Definition for fit results, A good fit should hold: OKAY.
Rotation matrix.
Definition: JRotation3D.hh:108
void getJEvt(const KM3NETDAQ::JDAQTimeslice &timeSlice, JFIT::JEvt &InPreFits, JFIT::JEvt &OutFits) const
Declaration of Member function that actually performs the reconstruction.
Auxiliary class to test history.
Definition: JHistory.hh:149
Data structure for fit of straight line in positive z-direction.
Definition: JLine3Z.hh:35
Container for historical events.
Definition: JHistory.hh:95
static parameter_type pT()
Definition: JLine1Z.hh:182
chi2 from JGandalf.cc
Linear fit of JFIT::JLine1Z.
JMuonGandalf(const JLANG::JSharedPointer< const JDETECTOR::JModuleRouter > &router, const JFIT::JMuonGandalfParameters_t &parameters, std::string pdfFile)
Parameterized constructor.
Definition: JMuonGandalf.hh:72
JLANG::JSharedPointer< const JDETECTOR::JModuleRouter > router_
Definition: JMuonGandalf.hh:52
angular resolution [rad] from JGandalf.cc
JFit getFit(const JHistory &history, const JTrack3D &track, const double Q, const int NDF, const double energy=0.0, const int status=0)
Get fit.
Definition: JEvtToolkit.hh:116
Basic data structure for L0 hit.
Definition of fit parameters from various applications.
Basic data structure for time and time over threshold information of hit.
static parameter_type pDX()
Definition: JLine3Z.hh:319
JAxis3D & rotate(const JRotation3D &R)
Rotate axis.
Definition: JAxis3D.hh:235
JFIT::JRegressor< JFIT::JLine3Z, JFIT::JGandalf > JRegressor_t
Definition: JMuonGandalf.hh:50
control parameter from JGandalf.cc
Auxiliary class to test history.
Definition: JHistory.hh:101
Auxiliary class for a hit with background rate value.
Definition: JHitW0.hh:21
JDirection3D getDirection(const JFit &fit)
Get direction.
Definition: JEvtToolkit.hh:75
Data time slice.
Auxiliary class to match data points with given model.
Definition: JModel.hh:24
angular resolution [rad] from JGandalf.cc
int debug
debug level
Definition: JSirene.cc:59
static parameter_type pY()
Definition: JLine1Z.hh:181
static parameter_type pDY()
Definition: JLine3Z.hh:320
counter_type advance(counter_type &counter, const counter_type value, const counter_type limit=std::numeric_limits< counter_type >::max())
Advance counter.
Definition: JCounter.hh:35
number of hits from JGandalf.cc
static parameter_type pX()
Definition: JLine1Z.hh:180
JPosition3D getPosition(const JFit &fit)
Get position.
Definition: JEvtToolkit.hh:63
JFIT::JMuonGandalfParameters_t parameters_
Definition: JMuonGandalf.hh:53
double getQuality(const double chi2, const int NDF)
Get quality of fit.
Definition: JEvtToolkit.hh:195
Data structure for set of track fit results.
Definition: JEvt.hh:312
Data structure for L0 hit.
Definition: JHitL0.hh:27
Data structure for fit of straight line paralel to z-axis.
Definition: JLine1Z.hh:27
void copy(const Head &from, JHead &to)
Copy header from from to to.
Definition: JHead.cc:40
JMuonGandalf()
Default constructor.
Definition: JMuonGandalf.hh:62
bool qualitySorter(const JFit &first, const JFit &second)
Comparison of fit results.
Definition: JEvtToolkit.hh:223
static bool compare(const JTRIGGER::JHitL0 &first, const JTRIGGER::JHitL0 &second)
Compare hits by PMT identifier and time.
Template L0 hit builder.
Definition: JBuildL0.hh:35
JPosition3D & rotate(const JRotation3D &R)
Rotate.
Definition: JPosition3D.hh:185