Jpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JMuonPDF.cc
Go to the documentation of this file.
1 #include <string>
2 #include <iostream>
3 #include <iomanip>
4 #include <vector>
5 #include <algorithm>
6 
9 
10 #include "JDAQ/JDAQEventIO.hh"
12 
13 #include "JDetector/JDetector.hh"
16 
17 #include "JLang/JPredicate.hh"
18 
19 #include "JTrigger/JHitL0.hh"
20 #include "JTrigger/JBuildL0.hh"
21 
22 #include "JSupport/JSupport.hh"
25 
26 #include "JPhysics/JPDF_t.hh"
27 
28 #include "JFit/JLine1Z.hh"
29 #include "JFit/JHitW0.hh"
30 #include "JFit/JModel.hh"
31 
32 #include "JReconstruction/JEvt.hh"
35 
36 #include "Jeep/JPrint.hh"
37 #include "Jeep/JParser.hh"
38 #include "Jeep/JMessage.hh"
39 
40 
41 namespace {
42 
43  /**
44  * Auxiliary data structure for sorting of hits.
45  */
46  static const struct {
47  /**
48  * Compare hits by PMT identifier and time.
49  *
50  * \param first first hit
51  * \param second second hit
52  * \return true if first before second; else false
53  */
54  template<class T>
55  bool operator()(const T& first, const T& second) const
56  {
57  using namespace std;
58  using namespace JPP;
59  using namespace KM3NETDAQ;
60 
61  if (equal_to<JDAQPMTIdentifier>()(first, second))
62  return less<JHit>()(first, second);
63  else
64  return less<JDAQPMTIdentifier>()(first, second);
65  }
66  } compare;
67 }
68 
69 
70 /**
71  * \file
72  *
73  * Program to evaluate hit probabilities.
74  *
75  * \author mdejong
76  */
77 int main(int argc, char **argv)
78 {
79  using namespace std;
80  using namespace JPP;
81  using namespace KM3NETDAQ;
82 
83  typedef JParallelFileScanner< JTypeList<JDAQEvent, JEvt> > JParallelFileScanner_t;
84  typedef JParallelFileScanner_t::multi_pointer_type multi_pointer_type;
85 
86  JParallelFileScanner_t inputFile;
87  JLimit_t& numberOfEvents = inputFile.getLimit();
88  string detectorFile;
89  string pdfFile;
91  int debug;
92 
93  try {
94 
95  parameters.numberOfPrefits = 1;
96 
97  JParser<> zap("Program to perform detector calibration using reconstructed muon trajectories.");
98 
99  zap['f'] = make_field(inputFile);
100  zap['a'] = make_field(detectorFile);
101  zap['n'] = make_field(numberOfEvents) = JLimit::max();
102  zap['P'] = make_field(pdfFile);
104  zap['d'] = make_field(debug) = 1;
105 
106  zap(argc, argv);
107  }
108  catch(const exception& error) {
109  FATAL(error.what() << endl);
110  }
111 
112  //setDAQLongprint(debug >= JEEP::debug_t);
113 
114  cout.tie(&cerr);
115 
116  if (parameters.numberOfPrefits != 1) {
117  WARNING("Number of prefits " << parameters.numberOfPrefits << " != " << 1 << endl);
118  }
119 
121 
122  try {
123  load(detectorFile, detector);
124  }
125  catch(const JException& error) {
126  FATAL(error);
127  }
128 
129  const JModuleRouter router(detector);
130 
131  JSummaryFileRouter summary(inputFile, parameters.R_Hz);
132 
133  const JMuonPDF_t pdf(pdfFile, parameters.TTS_ns);
134 
135  const JTimeRange T_ns(parameters.TMin_ns, parameters.TMax_ns);
136 
137  typedef vector<JHitL0> JDataL0_t;
138  typedef vector<JHitW0> JDataW0_t;
139 
140  const JBuildL0<JHitL0> buildL0;
141 
142  while (inputFile.hasNext()) {
143 
144  STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);
145 
146  multi_pointer_type ps = inputFile.next();
147 
148  JDAQEvent* tev = ps;
149  JEvt* in = ps;
150 
151  DEBUG("event: " << *tev << endl);
152 
153  summary.update(*tev);
154 
155  in->select(parameters.numberOfPrefits, qualitySorter);
156 
157  JDataL0_t dataL0;
158 
159  buildL0(*tev, router, true, back_inserter(dataL0));
160 
161  for (JEvt::const_iterator track = in->begin(); track != in->end(); ++track) {
162 
163  DEBUG("track: " << *track << endl);
164 
165  JRotation3D R (getDirection(*track));
166  JLine1Z tz(getPosition (*track).rotate(R), track->getT());
167  JRange<double> Z_m;
168 
169  if (track->hasW(JSTART_LENGTH_METRES)) {
170  Z_m = JRange<double>(0.0, track->getW(JSTART_LENGTH_METRES)) + JRange<double>(parameters.ZMin_m, parameters.ZMax_m);
171  }
172 
173  const JModel<JLine1Z> match(tz, parameters.roadWidth_m, T_ns, Z_m);
174 
175  // hit selection based on start value
176 
177  JDataW0_t data;
178 
179  for (JDataL0_t::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {
180 
181  double rate_Hz = summary.getRate(*i);
182 
183  if (rate_Hz <= 0.0) {
184  rate_Hz = summary.getRate();
185  }
186 
187  JHitW0 hit(*i, rate_Hz);
188 
189  hit.rotate(R);
190 
191  if (match(hit)) {
192  data.push_back(hit);
193  }
194  }
195 
196  // select first hit
197 
198  sort(data.begin(), data.end(), compare);
199 
200  JDataW0_t::iterator __end = unique(data.begin(), data.end(), equal_to<JDAQPMTIdentifier>());
201 
202 
203  double E_GeV = parameters.E_GeV;
204  /*
205  if (track->getE() > 0.1) {
206  E_GeV = track->getE();
207  }
208  */
209 
210  DEBUG("line1z: "
211  << FIXED(12,3) << tz.getX() << ' '
212  << FIXED(12,3) << tz.getY() << ' '
213  << FIXED(12,3) << tz.getZ() << ' '
214  << FIXED(12,3) << tz.getT() << ' '
215  << FIXED(12,1) << E_GeV << ' '
216  << FIXED( 8,3) << track->getQ() << endl);
217 
218  // move track to bary-center of hits
219  /*
220  int N = 0;
221  double z = 0.0;
222 
223  for (JDataW0_t::const_iterator hit = data.begin(); hit != __end; ++hit) {
224  N += 1;
225  z += hit->getZ();
226  }
227 
228  if (N != 0) {
229  tz.setZ(z/N, getSpeedOfLight());
230  }
231  */
232 
233  double Q = 0.0;
234 
235  for (JDataW0_t::const_iterator hit = data.begin(); hit != __end; ++hit) {
236 
237  const double x = hit->getX() - tz.getX();
238  const double y = hit->getY() - tz.getY();
239  const double z = hit->getZ() - tz.getZ();
240  const double R = sqrt(x*x + y*y);
241 
242  const double t1 = tz.getT() + (z + R * getTanThetaC()) * getInverseSpeedOfLight();
243 
244  JDirection3D u(hit->getDX(), hit->getDY(), hit->getDZ()); // PMT orientation
245 
246  u.rotate(JRotation3Z(-atan2(y,x))); // rotate PMT axis to x-z plane
247 
248  const double theta = u.getTheta();
249  const double phi = fabs(u.getPhi()); // rotational symmetry of Cherenkov cone
250 
251  //const double E = gWater.getE(E_GeV, z); // correct for energy loss
252  const double E = E_GeV;
253  const double dt = T_ns.constrain(hit->getT() - t1);
254 
255  JMuonPDF_t::result_type H1 = pdf.calculate(E, R, theta, phi, dt);
256  JMuonPDF_t::result_type H0(hit->getR() * 1e-9, 0.0, T_ns);
257 
258  H1 += H0; // signal + background
259 
260  const double chi2 = H1.getChi2() - H0.getChi2(); // likelihood ratio
261 
262  DEBUG("hit: "
263  << setw(10) << hit->getModuleID() << ':' << setw( 2) << setfill('0') << hit->getPMTAddress() << setfill(' ') << ' '
264  << FIXED(12,1) << E << ' '
265  << FIXED( 9,1) << R << ' '
266  << FIXED( 6,4) << theta << ' '
267  << FIXED( 6,4) << phi << ' '
268  << FIXED( 8,3) << dt << ' '
269  << FIXED(12,3) << chi2 << endl);
270 
271  Q += getQuality(chi2);
272  }
273 
274  DEBUG("quality: " << FIXED(8,3) << Q << ' ' << distance(data.begin(), __end) << endl);
275  }
276  }
277  STATUS(endl);
278 }
Utility class to parse command line options.
Definition: JParser.hh:1493
General exception.
Definition: JException.hh:23
#define WARNING(A)
Definition: JMessage.hh:65
ROOT TTree parameter settings.
Data structure for direction in three dimensions.
Definition: JDirection3D.hh:32
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
Template specialisation of L0 builder for JHitL0 data type.
Definition: JBuildL0.hh:102
#define STATUS(A)
Definition: JMessage.hh:63
void update(const JDAQHeader &header)
Update router.
Detector data structure.
Definition: JDetector.hh:80
Router for direct addressing of module data in detector data structure.
Rotation matrix.
Definition: JRotation3D.hh:111
esac print_variable DETECTOR INPUT_FILE OUTPUT_FILE CDF for TYPE in
Definition: JSirene.sh:45
double getRate() const
Get default rate.
*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
General purpose class for parallel reading of objects from a single file or multiple files...
Template specialisation of class JModel to match hit with muon trajectory along z-axis.
Definition: JFit/JModel.hh:34
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition: JParser.hh:63
Auxiliary data structure for floating point format specification.
Definition: JPrint.hh:481
Data structure for detector geometry and calibration.
double getQuality(const double chi2, const int N, const int NDF)
Get quality of fit.
Definition: JEvtToolkit.hh:206
then echo The file $DIR KM3NeT_00000001_00000000 root already please rename or remove it first
Basic data structure for L0 hit.
Rotation around Z-axis.
Definition: JRotation3D.hh:85
result_type calculate(const double E, const double R, const double theta, const double phi, const double t1) const
Get PDF.
Definition: JPDF_t.hh:117
const double getInverseSpeedOfLight()
Get inverse speed of light.
Definition: JConstants.hh:100
Auxiliary class for defining the range of iterations of objects.
Definition: JLimit.hh:41
I/O formatting auxiliaries.
JAxis3D & rotate(const JRotation3D &R)
Rotate axis.
Definition: JAxis3D.hh:236
Detector file.
Definition: JHead.hh:130
JFunction1D_t::result_type result_type
Definition: JPDF_t.hh:29
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1954
do set_variable OUTPUT_DIRECTORY $WORKDIR T
Parallel scanning of objects from a single file or multiple files according a format that follows fro...
Auxiliary class for a hit with background rate value.
Definition: JHitW0.hh:21
T constrain(argument_type x) const
Constrain value to range.
Definition: JRange.hh:351
File router for fast addressing of summary data.
void load(const JString &file_name, JDetector &detector)
Load detector from input file.
int debug
debug level
Definition: JSirene.cc:61
then usage $script[distance] fi case set_variable R
Definition: JDrawLED.sh:40
General purpose messaging.
Auxiliary data structure for muon PDF.
Definition: JPDF_t.hh:19
#define FATAL(A)
Definition: JMessage.hh:67
Direct access to module in detector data structure.
static const int JSTART_LENGTH_METRES
distance between first and last hits in metres from JStart.cc
Data structure for set of track fit results.
Definition: JEvt.hh:294
Utility class to parse command line options.
Data structure for fit of straight line paralel to z-axis.
Definition: JLine1Z.hh:27
JDirection3D getDirection(const Vec &v)
Get direction.
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:73
double u[N+1]
Definition: JPolint.hh:706
double getTanThetaC()
Get average tangent of Cherenkov angle of water.
Definition: JConstants.hh:133
JPosition3D & rotate(const JRotation3D &R)
Rotate.
Definition: JPosition3D.hh:185
void select(const JSelector_t &selector)
Select fits.
Definition: JEvt.hh:314
bool qualitySorter(const JRECONSTRUCTION::JFit &first, const JRECONSTRUCTION::JFit &second)
Comparison of fit results.
then usage $script[input file[working directory[option]]] nWhere option can be E
Definition: JMuonPostfit.sh:37
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
JPosition3D getPosition(const Vec &v)
Get position.
int main(int argc, char *argv[])
Definition: Main.cpp:15