Jpp  18.5.2
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Functions
JMuonPDF.cc File Reference

Program to evaluate hit probabilities. More...

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <algorithm>
#include "JDAQ/JDAQEventIO.hh"
#include "JDAQ/JDAQSummarysliceIO.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JModuleRouter.hh"
#include "JTrigger/JHitL0.hh"
#include "JTrigger/JBuildL0.hh"
#include "JSupport/JSupport.hh"
#include "JSupport/JParallelFileScanner.hh"
#include "JSupport/JSummaryFileRouter.hh"
#include "JPhysics/JPDF_t.hh"
#include "JPhysics/JNPE_t.hh"
#include "JFit/JLine1Z.hh"
#include "JFit/JModel.hh"
#include "JReconstruction/JHitW0.hh"
#include "JReconstruction/JEvt.hh"
#include "JReconstruction/JEvtToolkit.hh"
#include "JReconstruction/JMuonGandalfParameters_t.hh"
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

Go to the source code of this file.

Functions

int main (int argc, char **argv)
 

Detailed Description

Program to evaluate hit probabilities.

Author
mdejong

Definition in file JMuonPDF.cc.

Function Documentation

int main ( int  argc,
char **  argv 
)

Definition at line 69 of file JMuonPDF.cc.

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;
82  JMuonGandalfParameters_t parameters;
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);
95  zap['@'] = make_field(parameters) = JPARSER::initialised();
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 
111  JDetector detector;
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 }
JPARSER::JParser
Utility class to parse command line options.
Definition: JParser.hh:1514
JLANG::JException
General exception.
Definition: JException.hh:24
WARNING
#define WARNING(A)
Definition: JMessage.hh:65
Q
Q(UTCMax_s-UTCMin_s)-livetime_s
JGEOMETRY3D::JDirection3D
Data structure for direction in three dimensions.
Definition: JDirection3D.hh:33
JDETECTOR::JModule
Data structure for a composite optical module.
Definition: JModule.hh:67
E
then usage $script< input file >[option[primary[working directory]]] nWhere option can be E
Definition: JMuonPostfit.sh:40
JACOUSTICS::getQuality
double getQuality(const double chi2, const int NDF)
Get quality of fit.
Definition: JAcoustics/JEvtToolkit.hh:29
distance
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
Definition: PhysicsEvent.hh:434
JTRIGGER::JBuildL0< JHitL0 >
Template specialisation of L0 builder for JHitL0 data type.
Definition: JBuildL0.hh:102
STATUS
#define STATUS(A)
Definition: JMessage.hh:63
JDETECTOR::JDetector
Detector data structure.
Definition: JDetector.hh:89
JDETECTOR::JModuleRouter
Router for direct addressing of module data in detector data structure.
Definition: JModuleRouter.hh:36
JGEOMETRY3D::JRotation3D
Rotation matrix.
Definition: JRotation3D.hh:111
parameters
*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
JSUPPORT::JParallelFileScanner
General purpose class for parallel reading of objects from a single file or multiple files...
Definition: JParallelFileScanner.hh:33
JFIT::JModel< JLine1Z >
Template specialisation of class JModel to match hit with muon trajectory along z-axis.
Definition: JFit/JModel.hh:34
JPARSER::initialised
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition: JParser.hh:83
Y
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
Definition: detector-XY:fit1d.sh:33
FIXED
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:446
JGEOMETRY3D::JPosition3D::transform
void transform(const JRotation3D &R, const JVector3D &pos)
Transform position.
Definition: JPosition3D.hh:331
JRECONSTRUCTION::qualitySorter
bool qualitySorter(const JFit &first, const JFit &second)
Comparison of fit results.
Definition: JReconstruction/JEvtToolkit.hh:244
JGEOMETRY3D::JRotation3Z
Rotation around Z-axis.
Definition: JRotation3D.hh:85
JSUPPORT::JLimit
Auxiliary class for defining the range of iterations of objects.
Definition: JLimit.hh:41
JGEOMETRY3D::JAxis3D::rotate
JAxis3D & rotate(const JRotation3D &R)
Rotate axis.
Definition: JAxis3D.hh:225
JAANET::detector
Detector file.
Definition: JHead.hh:226
JAANET::getDirection
JDirection3D getDirection(const Vec &dir)
Get direction.
Definition: JAAnetToolkit.hh:211
JMuonPDF_t::result_type
JFunction1D_t::result_type result_type
Definition: JPDF_t.hh:145
JACOUSTICS::JEvt
Acoustic event fit.
Definition: JAcoustics/JEvt.hh:299
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1989
JLANG::JObjectID::getID
int getID() const
Get identifier.
Definition: JObjectID.hh:50
JAANET::getPosition
JPosition3D getPosition(const Vec &pos)
Get position.
Definition: JAAnetToolkit.hh:175
JSUPPORT::JSummaryFileRouter
File router for fast addressing of summary data.
Definition: JSummaryFileRouter.hh:36
JRECONSTRUCTION::JMuonGandalfParameters_t
Data structure for fit parameters.
Definition: JMuonGandalfParameters_t.hh:21
JMuonPDF_t
Auxiliary data structure for muon PDF.
Definition: JPDF_t.hh:135
FATAL
#define FATAL(A)
Definition: JMessage.hh:67
R
then JCookie sh JDataQuality D $DETECTOR_ID R
Definition: JDataQuality.sh:41
JRECONSTRUCTION::JHitW0
Auxiliary class for a hit with background rate value.
Definition: JHitW0.hh:21
JSTART_LENGTH_METRES
static const int JSTART_LENGTH_METRES
distance between first and last hits in metres from JStart.cc
Definition: fitparameters.hh:25
JDETECTOR::load
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
Definition: JDetectorToolkit.hh:483
JMuonNPE_t
Definition: JNPE_t.hh:24
in
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
JDETECTOR::JModule::transform
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
JPHYSICS::getInverseSpeedOfLight
const double getInverseSpeedOfLight()
Get inverse speed of light.
Definition: JPhysics/JConstants.hh:123
chi2
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
JFIT::JLine1Z
Data structure for fit of straight line paralel to z-axis.
Definition: JLine1Z.hh:27
JPHYSICS::getTanThetaC
double getTanThetaC()
Get average tangent of Cherenkov angle of water corresponding to group velocity.
Definition: JPhysics/JConstants.hh:156
JGEOMETRY3D::JPosition3D
Data structure for position in three dimensions.
Definition: JPosition3D.hh:36
JSUPPORT::JLimit::getLimit
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:84
detector
do set_variable DETECTOR_TXT $WORKDIR detector
Definition: detector-XY:fit1d.sh:38
JTOOLS::u
double u[N+1]
Definition: JPolint.hh:865
JGEOMETRY3D::JPosition3D::rotate
JPosition3D & rotate(const JRotation3D &R)
Rotate.
Definition: JPosition3D.hh:186
SCIENTIFIC
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:486
debug
int debug
debug level
Definition: archive-put-wiki-detectors.sh:92
DEBUG
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
Generated by   doxygen 1.8.5