Jpp  18.5.0
the software that should make you happy
 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 
7 #include "JDAQ/JDAQEventIO.hh"
9 
10 #include "JDetector/JDetector.hh"
13 
14 #include "JTrigger/JHitL0.hh"
15 #include "JTrigger/JBuildL0.hh"
16 
17 #include "JSupport/JSupport.hh"
20 
21 #include "JPhysics/JPDF_t.hh"
22 #include "JPhysics/JNPE_t.hh"
23 
24 #include "JFit/JLine1Z.hh"
25 #include "JFit/JModel.hh"
26 
28 #include "JReconstruction/JEvt.hh"
31 
32 #include "Jeep/JPrint.hh"
33 #include "Jeep/JParser.hh"
34 #include "Jeep/JMessage.hh"
35 
36 
37 namespace {
38 
39  /**
40  * Auxiliary data structure for sorting of hits.
41  */
42  static const struct {
43  /**
44  * Compare hits by PMT identifier and time.
45  *
46  * \param first first hit
47  * \param second second hit
48  * \return true if first before second; else false
49  */
50  template<class T>
51  bool operator()(const T& first, const T& second) const
52  {
53  if (first.getPMTIdentifier() == second.getPMTIdentifier())
54  return first.getT() < second.getT();
55  else
56  return first.getPMTIdentifier() < second.getPMTIdentifier();
57  }
58  } compare;
59 }
60 
61 
62 /**
63  * \file
64  *
65  * Program to evaluate hit probabilities.
66  *
67  * \author mdejong
68  */
69 int main(int argc, char **argv)
70 {
71  using namespace std;
72  using namespace JPP;
73  using namespace KM3NETDAQ;
74 
75  typedef JParallelFileScanner< JTypeList<JDAQEvent, JEvt> > JParallelFileScanner_t;
76  typedef JParallelFileScanner_t::multi_pointer_type multi_pointer_type;
77 
78  JParallelFileScanner_t inputFile;
79  JLimit_t& numberOfEvents = inputFile.getLimit();
80  string detectorFile;
81  string pdfFile;
83  int debug;
84 
85  try {
86 
87  parameters.numberOfPrefits = 1;
88 
89  JParser<> zap("Program to evaluate hit probabilities.");
90 
91  zap['f'] = make_field(inputFile);
92  zap['a'] = make_field(detectorFile);
93  zap['n'] = make_field(numberOfEvents) = JLimit::max();
94  zap['P'] = make_field(pdfFile);
96  zap['d'] = make_field(debug) = 1;
97 
98  zap(argc, argv);
99  }
100  catch(const exception& error) {
101  FATAL(error.what() << endl);
102  }
103 
104  //setDAQLongprint(debug >= JEEP::debug_t);
105 
106 
107  if (parameters.numberOfPrefits != 1) {
108  WARNING("Number of prefits " << parameters.numberOfPrefits << " != " << 1 << endl);
109  }
110 
112 
113  try {
114  load(detectorFile, detector);
115  }
116  catch(const JException& error) {
117  FATAL(error);
118  }
119 
120  const JModuleRouter router(detector);
121 
122  JSummaryFileRouter summary(inputFile, parameters.R_Hz);
123 
124  const JMuonPDF_t pdf(pdfFile, parameters.TTS_ns);
125  const JMuonNPE_t npe(pdfFile);
126 
127  const JTimeRange T_ns(parameters.TMin_ns, parameters.TMax_ns);
128 
129  typedef vector<JHitL0> JDataL0_t;
130  typedef vector<JHitW0> JDataW0_t;
131 
132  const JBuildL0<JHitL0> buildL0;
133 
134  while (inputFile.hasNext()) {
135 
136  STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);
137 
138  multi_pointer_type ps = inputFile.next();
139 
140  JDAQEvent* tev = ps;
141  JEvt* in = ps;
142 
143  DEBUG("event: " << *tev << endl);
144 
145  summary.update(*tev);
146 
147  in->select(parameters.numberOfPrefits, qualitySorter);
148 
149  JDataL0_t dataL0;
150 
151  buildL0(*tev, router, true, back_inserter(dataL0));
152 
153  for (JEvt::const_iterator track = in->begin(); track != in->end(); ++track) {
154 
155  DEBUG("track: " << *track << endl);
156 
157  JRotation3D R (getDirection(*track));
158  JLine1Z tz(getPosition (*track).rotate(R), track->getT());
159  JRange<double> Z_m;
160 
161  if (track->hasW(JSTART_LENGTH_METRES)) {
162  Z_m = JRange<double>(0.0, track->getW(JSTART_LENGTH_METRES)) + JRange<double>(parameters.ZMin_m, parameters.ZMax_m);
163  }
164 
165  const JModel<JLine1Z> match(tz, parameters.roadWidth_m, T_ns, Z_m);
166 
167  // hit selection based on start value
168 
169  JDataW0_t data;
170 
171  for (JDataL0_t::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {
172 
173  JHitW0 hit(*i, summary.getRate(*i));
174 
175  hit.rotate(R);
176 
177  if (match(hit)) {
178  data.push_back(hit);
179  }
180  }
181 
182  // select first hit in PMT
183 
184  sort(data.begin(), data.end(), compare);
185 
186  JDataW0_t::iterator __end = unique(data.begin(), data.end(), equal_to<JDAQPMTIdentifier>());
187 
188 
189  double E_GeV = parameters.E_GeV;
190  /*
191  if (track->getE() > 0.1) {
192  E_GeV = track->getE();
193  }
194  */
195 
196  DEBUG("line1z: "
197  << FIXED(12,3) << tz.getX() << ' '
198  << FIXED(12,3) << tz.getY() << ' '
199  << FIXED(12,3) << tz.getZ() << ' '
200  << FIXED(12,3) << tz.getT() << ' '
201  << FIXED(12,1) << E_GeV << ' '
202  << FIXED( 8,3) << track->getQ() << endl);
203 
204  // move track to bary-center of hits
205  /*
206  int N = 0;
207  double z = 0.0;
208 
209  for (JDataW0_t::const_iterator hit = data.begin(); hit != __end; ++hit) {
210  N += 1;
211  z += hit->getZ();
212  }
213 
214  if (N != 0) {
215  tz.setZ(z/N, getSpeedOfLight());
216  }
217  */
218 
219  double Q = 0.0;
220 
221  for (JDataW0_t::const_iterator hit = data.begin(); hit != __end; ++hit) {
222 
223  const double x = hit->getX() - tz.getX();
224  const double y = hit->getY() - tz.getY();
225  const double z = hit->getZ() - tz.getZ();
226  const double R = sqrt(x*x + y*y);
227 
228  const double t1 = tz.getT() + (z + R * getTanThetaC()) * getInverseSpeedOfLight();
229 
230  JDirection3D u(hit->getDX(), hit->getDY(), hit->getDZ()); // PMT orientation
231 
232  u.rotate(JRotation3Z(-atan2(y,x))); // rotate PMT axis to x-z plane
233 
234  const double theta = u.getTheta();
235  const double phi = fabs(u.getPhi()); // rotational symmetry of Cherenkov cone
236 
237  //const double E = gWater.getE(E_GeV, z); // correct for energy loss
238  const double E = E_GeV;
239  const double dt = T_ns.constrain(hit->getT() - t1);
240 
241  JMuonPDF_t::result_type H1 = pdf.calculate(E, R, theta, phi, dt);
242  JMuonPDF_t::result_type H0(hit->getR() * 1e-9, 0.0, T_ns);
243 
244  H1 += H0; // signal + background
245 
246  const double chi2 = H1.getChi2() - H0.getChi2(); // likelihood ratio
247 
248  DEBUG("hit: "
249  << setw(10) << hit->getModuleID() << ':' << setw( 2) << setfill('0') << hit->getPMTAddress() << setfill(' ') << ' '
250  << FIXED(12,1) << E << ' '
251  << FIXED( 9,1) << R << ' '
252  << FIXED( 6,4) << theta << ' '
253  << FIXED( 6,4) << phi << ' '
254  << FIXED( 8,3) << dt << ' '
255  << FIXED(12,3) << chi2 << endl);
256 
257  Q += getQuality(chi2);
258  }
259 
260  DEBUG("quality: " << FIXED(8,3) << Q << ' ' << distance(data.begin(), __end) << endl);
261 
262 
263  double Y = 0.0;
264 
265  for (JDetector::const_iterator i = detector.begin(); i != detector.end(); ++i) {
266 
267  JPosition3D pos(i->getPosition());
268 
269  pos.transform(R, tz.getPosition());
270 
271  if (pos.getX() <= parameters.roadWidth_m) {
272 
273  JModule module = *i;
274 
275  module.transform(R, tz.getPosition());
276 
277  for (JModule::iterator pmt = module.begin(); pmt != module.end(); ++pmt) {
278 
279  const double R = sqrt(pmt->getX()*pmt->getX() + pmt->getY()*pmt->getY());
280  const double theta = pmt->getTheta();
281  const double phi = fabs(pmt->getPhi());
282  const double y = npe.calculate(1.0, R, theta, phi); // MIP
283 
284  DEBUG("PMT: "
285  << setw(10) << module.getID() << ':' << setw( 2) << setfill('0') << distance(module.begin(),pmt) << setfill(' ') << ' '
286  << FIXED(9,1) << R << ' '
287  << FIXED(6,4) << theta << ' '
288  << FIXED(6,4) << phi << ' '
289  << SCIENTIFIC(7,2) << y << endl);
290 
291  Y += y;
292  }
293  }
294  }
295 
296  DEBUG("npe: " << FIXED(7,1) << Y << endl);
297  }
298  }
299  STATUS(endl);
300 }
Utility class to parse command line options.
Definition: JParser.hh:1514
General exception.
Definition: JException.hh:24
#define WARNING(A)
Definition: JMessage.hh:65
double calculate(const double E, const double R, const double theta, const double phi) const
Get PDF.
Definition: JNPE_t.hh:108
Q(UTCMax_s-UTCMin_s)-livetime_s
int main(int argc, char *argv[])
Definition: Main.cc:15
ROOT TTree parameter settings of various packages.
Data structure for direction in three dimensions.
Definition: JDirection3D.hh:33
Data structure for a composite optical module.
Definition: JModule.hh:67
then usage $script< input file >[option[primary[working directory]]] nWhere option can be E
Definition: JMuonPostfit.sh:40
double getQuality(const double chi2, const int NDF)
Get quality of fit.
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:89
Router for direct addressing of module data in detector data structure.
Rotation matrix.
Definition: JRotation3D.hh:111
double getRate() const
Get default rate.
Auxiliary data structure for muon PDF.
*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:83
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
Data structure for detector geometry and calibration.
void transform(const JRotation3D &R, const JVector3D &pos)
Transform position.
Definition: JPosition3D.hh:331
bool qualitySorter(const JFit &first, const JFit &second)
Comparison of fit results.
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:233
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:225
Detector file.
Definition: JHead.hh:226
JDirection3D getDirection(const Vec &dir)
Get direction.
JFunction1D_t::result_type result_type
Definition: JPDF_t.hh:145
Acoustic event fit.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1989
do set_variable OUTPUT_DIRECTORY $WORKDIR T
int getID() const
Get identifier.
Definition: JObjectID.hh:50
Parallel scanning of objects from a single file or multiple files according a format that follows fro...
JPosition3D getPosition(const Vec &pos)
Get position.
T constrain(argument_type x) const
Constrain value to range.
Definition: JRange.hh:350
File router for fast addressing of summary data.
Auxiliary data structure for muon PDF.
General purpose messaging.
Auxiliary data structure for muon PDF.
Definition: JPDF_t.hh:135
#define FATAL(A)
Definition: JMessage.hh:67
Direct access to module in detector data structure.
then JCookie sh JDataQuality D $DETECTOR_ID R
Definition: JDataQuality.sh:41
Auxiliary class for a hit with background rate value.
Definition: JHitW0.hh:21
static const int JSTART_LENGTH_METRES
distance between first and last hits in metres from JStart.cc
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
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:48
Utility class to parse command line options.
void transform(const JRotation3D &R, const JVector3D &pos)
Transformation of geometry (see method JGEOMETRY3D::JPosition3D::transform(const JRotation3D&amp;, const JVector3D&amp;)).
Definition: JModule.hh:345
const double getInverseSpeedOfLight()
Get inverse speed of light.
then if[[!-f $DETECTOR]] then JDetector sh $DETECTOR fi cat $WORKDIR trigger_parameters txt<< EOFtrigger3DMuon.enabled=1;trigger3DMuon.numberOfHits=5;trigger3DMuon.gridAngle_deg=1;ctMin=0.0;TMaxLocal_ns=15.0;EOF set_variable TRIGGEREFFICIENCY_TRIGGERED_EVENTS_ONLY INPUT_FILES=() for((i=1;$i<=$NUMBER_OF_RUNS;++i));do JSirene.sh $DETECTOR $JPP_DATA/genhen.km3net_wpd_V2_0.evt.gz $WORKDIR/sirene_ ${i}.root JTriggerEfficiency.sh $DETECTOR $DETECTOR $WORKDIR/sirene_ ${i}.root $WORKDIR/trigger_efficiency_ ${i}.root $WORKDIR/trigger_parameters.txt $JPP_DATA/PMT_parameters.txt INPUT_FILES+=($WORKDIR/trigger_efficiency_ ${i}.root) done for ANGLE_DEG in $ANGLES_DEG[*];do set_variable SIGMA_NS 3.0 set_variable OUTLIERS 3 set_variable OUTPUT_FILE $WORKDIR/matrix\[${ANGLE_DEG}\deg\].root $JPP_DIR/examples/JReconstruction-f"$INPUT_FILES[*]"-o $OUTPUT_FILE-S ${SIGMA_NS}-A ${ANGLE_DEG}-O ${OUTLIERS}-d ${DEBUG}--!fiif[[$OPTION=="plot"]];then if((0));then for H1 in h0 h1;do JPlot1D-f"$WORKDIR/matrix["${^ANGLES_DEG}" deg].root:${H1}"-y"1 2e3"-Y-L TR-T""-\^"number of events [a.u.]"-> o chi2
Definition: JMatrixNZ.sh:106
Data structure for fit of straight line paralel to z-axis.
Definition: JLine1Z.hh:27
double getTanThetaC()
Get average tangent of Cherenkov angle of water corresponding to group velocity.
Data structure for position in three dimensions.
Definition: JPosition3D.hh:36
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:84
do set_variable DETECTOR_TXT $WORKDIR detector
double u[N+1]
Definition: JPolint.hh:865
JPosition3D & rotate(const JRotation3D &R)
Rotate.
Definition: JPosition3D.hh:186
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:486
int debug
debug level
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62