Jpp  17.0.0-rc.1
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JAAnetToolkit.hh
Go to the documentation of this file.
1 #ifndef __JAANET__JAANETTOOLKIT__
2 #define __JAANET__JAANETTOOLKIT__
3 
4 #include <cstdlib>
5 #include <algorithm>
6 
14 
15 #include "JLang/JException.hh"
16 
19 #include "JGeometry3D/JAxis3D.hh"
20 #include "JGeometry3D/JTrack3D.hh"
21 #include "JGeometry3D/JTrack3E.hh"
22 #include "JGeometry3D/JVertex3D.hh"
24 
25 #include "JPhysics/JConstants.hh"
26 #include "JPhysics/JTimeRange.hh"
27 #include "JPhysics/JGeane.hh"
28 
29 #include "JAAnet/JParticleTypes.hh"
30 #include "JAAnet/JPDB.hh"
31 
32 
33 /**
34  * \file
35  *
36  * Definition of hit and track types and auxiliary methods for handling Monte Carlo data.
37  * \author mdejong
38  */
39 namespace JAANET {}
40 namespace JPP { using namespace JAANET; }
41 
42 /**
43  * Extensions to Evt data format.
44  */
45 namespace JAANET {
46 
49 
57 
61 
62  /**
63  * Index of rescaled weight.
64  */
65  static const int RESCALED_WEIGHT_INDEX = 3;
66 
67  /**
68  * AAnet shower fit reconstruction type.
69  */
70  static const int AASHOWER_RECONSTRUCTION_TYPE = 101;
71 
72  /**
73  * Enumeration of hit types based on km3 codes.
74  */
75  enum JHitType_t {
76  HIT_TYPE_MUON_DIRECT = +5, //!< Direct light from muon
77  HIT_TYPE_MUON_SCATTERED = -5, //!< Scattered light from muon
78  HIT_TYPE_DELTARAYS_DIRECT = +4, //!< Direct light from delta-rays
79  HIT_TYPE_DELTARAYS_SCATTERED = -4, //!< Scattered light from delta-rays
80  HIT_TYPE_BREMSSTRAHLUNG_DIRECT = +3, //!< Direct light from Bremsstrahlung
81  HIT_TYPE_BREMSSTRAHLUNG_SCATTERED = -3, //!< Scattered light from Bremsstrahlung
82  HIT_TYPE_SHOWER_DIRECT = +99, //!< Direct light from primary shower
83  HIT_TYPE_SHOWER_SCATTERED = -99, //!< Scattered light from primary shower
84  HIT_TYPE_NOISE = -1, //!< Random noise
85  HIT_TYPE_UNKNOWN = 0 //!< Unknown source
86  };
87 
88 
89  /**
90  * Get true time of hit.
91  *
92  * \param hit hit
93  * \return true time of hit
94  */
95  inline double getTime(const Hit& hit)
96  {
97  return hit.t;
98  }
99 
100 
101  /**
102  * Get true charge of hit.
103  *
104  * \param hit hit
105  * \return true charge of hit
106  */
107  inline double getNPE(const Hit& hit)
108  {
109  return hit.a;
110  }
111 
112 
113  /**
114  * Verify hit origin.
115  *
116  * \param hit hit
117  * \return true if noise; else false
118  */
119  inline bool is_noise(const Hit& hit)
120  {
121  return hit.type == HIT_TYPE_NOISE;
122  }
123 
124 
125  /**
126  * Get time range (i.e.\ time between earliest and latest hit) of Monte Carlo event.\n
127  * Note that the global event time in not taken into account.
128  *
129  * \param event event
130  * \return time range
131  */
132  inline JTimeRange getTimeRange(const Evt& event)
133  {
135 
136  for (std::vector<Hit>::const_iterator hit = event.mc_hits.begin(); hit != event.mc_hits.end(); ++hit) {
137  if (!is_noise(*hit)) {
138  time_range.include(getTime(*hit));
139  }
140  }
141 
142  return time_range;
143  }
144 
145 
146  /**
147  * Get time range (i.e.\ time between earliest and latest hit) of Monte Carlo event.\n
148  * The resulting time range is larger than or equal to the given time window.\n
149  * Note that the global event time in not taken into account.
150  *
151  * \param event event
152  * \param T_ns time window [ns]
153  * \return time range
154  */
155  inline JTimeRange getTimeRange(const Evt& event, const JTimeRange& T_ns)
156  {
157  JTimeRange time_range = getTimeRange(event);
158 
159  if (!time_range.is_valid()) {
160  time_range = T_ns;
161  }
162 
163  if (time_range.getLength() < T_ns.getLength()) {
164  time_range.setLength(T_ns.getLength());
165  }
166 
167  return time_range;
168  }
169 
170 
171  /**
172  * Get position.
173  *
174  * \param pos position
175  * \return position
176  */
177  inline JPosition3D getPosition(const Vec& pos)
178  {
179  return JPosition3D(pos.x, pos.y, pos.z);
180  }
181 
182 
183  /**
184  * Get position.
185  *
186  * \param pos position
187  * \return position
188  */
189  inline Vec getPosition(const JPosition3D& pos)
190  {
191  return Vec(pos.getX(), pos.getY(), pos.getZ());
192  }
193 
194 
195  /**
196  * Get position.
197  *
198  * \param track track
199  * \return position
200  */
201  inline JPosition3D getPosition(const Trk& track)
202  {
203  return getPosition(track.pos);
204  }
205 
206 
207  /**
208  * Get direction.
209  *
210  * \param dir direction
211  * \return direction
212  */
213  inline JDirection3D getDirection(const Vec& dir)
214  {
215  return JDirection3D(dir.x, dir.y, dir.z);
216  }
217 
218 
219  /**
220  * Get direction.
221  *
222  * \param dir direction
223  * \return direction
224  */
225  inline Vec getDirection(const JDirection3D& dir)
226  {
227  return Vec(dir.getDX(), dir.getDY(), dir.getDZ());
228  }
229 
230 
231  /**
232  * Get direction.
233  *
234  * \param track track
235  * \return direction
236  */
237  inline JDirection3D getDirection(const Trk& track)
238  {
239  return getDirection(track.dir);
240  }
241 
242 
243  /**
244  * Get axis.
245  *
246  * \param track track
247  * \return axis
248  */
249  inline JAxis3D getAxis(const Trk& track)
250  {
251  return JAxis3D(getPosition(track), getDirection(track));
252  }
253 
254 
255  /**
256  * Get track.
257  *
258  * \param track track
259  * \return track
260  */
261  inline JTrack3E getTrack(const Trk& track)
262  {
263  return JTrack3E(JTrack3D(getAxis(track), track.t), track.E);
264  }
265 
266 
267  /**
268  * Get vertex.
269  *
270  * \param track track
271  * \return vertex
272  */
273  inline JVertex3D getVertex(const Trk& track)
274  {
275  return JVertex3D(getPosition(track), track.t);
276  }
277 
278 
279  /**
280  * Get transformation.
281  *
282  * \param track track
283  * \return transformation
284  */
286  {
287  return JTransformation3D(getPosition(track), getDirection(track));
288  }
289 
290 
291  /**
292  * Get track information.
293  *
294  * \param track track
295  * \param index index
296  * \param value default value
297  * \return actual value
298  */
299  inline double getW(const Trk& track, const size_t index, const double value)
300  {
301  if (index < track.fitinf.size())
302  return track.fitinf[index];
303  else
304  return value;
305  }
306 
307 
308  /**
309  * Test whether given track is a photon or neutral pion.
310  *
311  * \param track track
312  * \return true if photon or neutral pion; else false
313  */
314  inline bool is_photon (const Trk& track) { return ( track.type == TRACK_TYPE_PHOTON ||
315  abs(track.type) == TRACK_TYPE_NEUTRAL_PION); }
316 
317  /**
318  * Test whether given track is a neutrino.
319  *
320  * \param track track
321  * \return true if neutrino; else false
322  */
323  inline bool is_neutrino(const Trk& track) { return (abs(track.type) == TRACK_TYPE_NUE ||
324  abs(track.type) == TRACK_TYPE_NUMU ||
325  abs(track.type) == TRACK_TYPE_NUTAU); }
326 
327  /**
328  * Test whether given track is a (anti-)electron.
329  *
330  * \param track track
331  * \return true if (anti-)electron; else false
332  */
333  inline bool is_electron(const Trk& track) { return (abs(track.type) == TRACK_TYPE_ELECTRON); }
334 
335  /**
336  * Test whether given track is a (anti-)muon.
337  *
338  * \param track track
339  * \return true if (anti-)muon; else false
340  */
341  inline bool is_muon (const Trk& track) { return (abs(track.type) == TRACK_TYPE_MUON); }
342 
343  /**
344  * Test whether given track is a (anti-)tau.
345  *
346  * \param track track
347  * \return true if (anti-)tau; else false
348  */
349  inline bool is_tau (const Trk& track) { return (abs(track.type) == TRACK_TYPE_TAU); }
350 
351  /**
352  * Test whether given track is a charged pion.
353  *
354  * \param track track
355  * \return true if charged pion; else false
356  */
357  inline bool is_pion (const Trk& track) { return (abs(track.type) == TRACK_TYPE_CHARGED_PION_PLUS); }
358 
359  /**
360  * Test whether given track is a (anti-)proton.
361  *
362  * \param track track
363  * \return true if (anti-)proton; else false
364  */
365  inline bool is_proton (const Trk& track) { return (abs(track.type) == TRACK_TYPE_PROTON); }
366 
367  /**
368  * Test whether given track is a (anti-)neutron.
369  *
370  * \param track track
371  * \return true if (anti-)neutron; else false
372  */
373  inline bool is_neutron (const Trk& track) { return (abs(track.type) == TRACK_TYPE_NEUTRON); }
374 
375  /**
376  * Test whether given track is a lepton
377  *
378  * \param track track
379  * \return true if lepton; else fails
380  */
381  inline bool is_lepton (const Trk& track) { return (is_neutrino(track) ||
382  is_electron(track) ||
383  is_muon (track) ||
384  is_tau (track)); }
385 
386  /**
387  * Test whether given track is a charged lepton
388  *
389  * \param track track
390  * \return true if lepton; else fails
391  */
392  inline bool is_charged_lepton (const Trk& track) { return (is_electron(track) ||
393  is_muon (track) ||
394  is_tau (track)); }
395 
396  /**
397  * Test whether given track is a hadron
398  *
399  * \param track track
400  * \return true if hadron; else fails
401  */
402  inline bool is_hadron (const Trk& track) { return (abs(track.type) != TRACK_TYPE_PHOTON &&
403  !is_lepton(track)); }
404 
405  /**
406  * Test whether given track is a meson (is hadron and third digit of PDG code is zero)
407  *
408  * \param track track
409  * \return true if hadron; else fails
410  */
411  inline bool is_meson (const Trk& track) { return (is_hadron(track) &&
412  ((int)(track.type / 1000)) % 10 == 0); }
413 
414  /**
415  * Test whether given track is a baryon (is hadron and third digit of PDG code is not zero)
416  *
417  * \param track track
418  * \return true if hadron; else fails
419  */
420  inline bool is_baryon (const Trk& track) { return (is_hadron(track) &&
421  ((int)(track.type / 1000)) % 10 != 0); }
422 
423  /**
424  * Test whether given track corresponds to given particle type.
425  *
426  * \param track track
427  * \param type particle type
428  * \return true if matches the corresponding particle; else false
429  */
430  inline bool has_particleID (const Trk& track, const int type) { return (track.type == type); }
431 
432  /**
433  * Test whether given event has an incoming neutrino.
434  *
435  * \param evt event
436  * \return true if neutrino; else false
437  */
438  inline bool has_neutrino(const Evt& evt)
439  {
440  return !evt.mc_trks.empty() && is_neutrino(evt.mc_trks[0]);
441  }
442 
443  /**
444  * Get incoming neutrino.
445  *
446  * \param evt event
447  * \return neutrino
448  */
449  inline const Trk& get_neutrino(const Evt& evt)
450  {
451  if (has_neutrino(evt))
452  return evt.mc_trks[0];
453  else
454  THROW(JIndexOutOfRange, "JAANET::get_neutrino(): Event " << evt.id << " has no neutrino.");
455  }
456 
457  /**
458  * Count the number of electrons in a given event
459  *
460  * \param evt event
461  * \return number of electrons
462  */
463  inline int count_electrons(const Evt& evt)
464  {
465  return std::count_if(evt.mc_trks.begin(), evt.mc_trks.end(), is_electron);
466  }
467 
468  /**
469  * Test whether given event has an electron.
470  *
471  * \param evt event
472  * \return true if event has electron; else false
473  */
474  inline bool has_electron(const Evt& evt)
475  {
476  return count_electrons(evt) != 0;
477  }
478 
479  /**
480  * Get first electron from the event tracklist
481  *
482  * \param evt event
483  * \return electron
484  */
485  inline const Trk& get_electron(const Evt& evt)
486  {
487  if (count_electrons(evt) > 0)
488  return *std::find_if(evt.mc_trks.begin(), evt.mc_trks.end(), is_electron);
489  else
490  THROW(JIndexOutOfRange, "This event has no electron.");
491  }
492 
493  /**
494  * Test whether given hit was produced by an electron
495  *
496  * Warning: This method will only work with the output of KM3Sim.
497  * For JSirene or KM3, you should check if <tt>track.id == hit.origin</tt> instead.
498  *
499  * \param hit hit
500  * \return true if hit was produced by an electron; else false
501  */
502  inline bool from_electron(const Hit& hit)
503  {
504  return hit.type == GEANT4_TYPE_ELECTRON || hit.type == GEANT4_TYPE_ANTIELECTRON;
505  }
506 
507  /**
508  * Count the number of muons in a given event
509  *
510  * \param evt event
511  * \return number of muons
512  */
513  inline int count_muons(const Evt& evt)
514  {
515  return std::count_if(evt.mc_trks.begin(), evt.mc_trks.end(), is_muon);
516  }
517 
518  /**
519  * Test whether given event has a muon.
520  *
521  * \param evt event
522  * \return true if event has muon; else false
523  */
524  inline bool has_muon(const Evt& evt)
525  {
526  return count_muons(evt) != 0;
527  }
528 
529  /**
530  * Get first muon from the event tracklist
531  *
532  * \param evt event
533  * \return muon
534  */
535  inline const Trk& get_muon(const Evt& evt)
536  {
537  if (count_muons(evt) > 0)
538  return *std::find_if(evt.mc_trks.begin(), evt.mc_trks.end(), is_muon);
539  else
540  THROW(JIndexOutOfRange, "This event has no muon.");
541  }
542 
543  /**
544  * Test whether given hit was produced by a muon
545  *
546  * Warning: This method will only work with the output of KM3Sim.
547  * For JSirene or KM3, you should check if <tt>track.id == hit.origin</tt> instead.
548  *
549  * \param hit hit
550  * \return true if hit was produced by a muon; else false
551  */
552  inline bool from_muon(const Hit& hit)
553  {
554  return hit.type == GEANT4_TYPE_MUON || hit.type == GEANT4_TYPE_ANTIMUON;
555  }
556 
557  /**
558  * Count the number of taus in a given event
559  *
560  * \param evt event
561  * \return number of taus
562  */
563  inline int count_taus(const Evt& evt)
564  {
565  return std::count_if(evt.mc_trks.begin(), evt.mc_trks.end(), is_tau);
566  }
567 
568  /**
569  * Test whether given event has a tau.
570  *
571  * \param evt event
572  * \return true if event has tau; else false
573  */
574  inline bool has_tau(const Evt& evt)
575  {
576  return count_taus(evt) != 0;
577  }
578 
579  /**
580  * Get first tau from the event tracklist
581  *
582  * \param evt event
583  * \return tau
584  */
585  inline const Trk& get_tau(const Evt& evt)
586  {
587  if (count_taus(evt) > 0)
588  return *std::find_if(evt.mc_trks.begin(), evt.mc_trks.end(), is_tau);
589  else
590  THROW(JIndexOutOfRange, "This event has no tau.");
591  }
592 
593  /**
594  * Test whether given hit was produced by a tau
595  *
596  * Warning: This method will only work with the output of KM3Sim.
597  * For JSirene or KM3, you should check if track.id == hit.origin instead.
598  *
599  * \param hit hit
600  * \return true if hit was produced by a tau; else false
601  */
602  inline bool from_tau(const Hit& hit)
603  {
604  return hit.type == GEANT4_TYPE_TAU || hit.type == GEANT4_TYPE_ANTITAU;
605  }
606 
607  /**
608  * Count the number of hadrons in a given event (not including neutral pions)
609  *
610  * \param evt event
611  * \return number of hadrons
612  */
613  inline int count_hadrons(const Evt& evt)
614  {
615  return std::count_if(evt.mc_trks.begin(), evt.mc_trks.end(), is_hadron);
616  }
617 
618  /**
619  * Test whether given hit was produced by a hadronic shower
620  *
621  * Warning: This method will only work with the output of KM3Sim.
622  * For JSirene or KM3, you should check if track.id == hit.origin instead.
623  *
624  * \param hit hit
625  * \return true if hit was produced by a hadron; else false
626  */
627  inline bool from_hadron(const Hit& hit)
628  {
629  return (!from_electron(hit) && !from_muon(hit) && !from_tau(hit));
630  }
631 
632 
633  /**
634  * Add time offset to mc event;
635  * according to current definition, the absolute time of the event is defined by the track "t" attribute;
636  * this could change in the future if the global attribute mc_t is assigned to this purpose.
637  *
638  * \param evt event
639  * \param tOff time offset [ns]
640  */
641  inline void add_time_offset(Evt& evt, const double tOff)
642  {
643  for (std::vector<Trk>::iterator p = evt.mc_trks.begin(); p != evt.mc_trks.end(); p++) {
644  p->t += tOff;
645  }
646  }
647 
648 
649  /**
650  * Get kinetic energy of particle with given mass.
651  *
652  * \param E energy [GeV]
653  * \param m mass [GeV]
654  * \return energy [GeV]
655  */
656  inline double getKineticEnergy(const double E, const double m)
657  {
658  if (E > m) {
659  return sqrt((E - m) * (E + m));
660  } else {
661  return 0.0;
662  }
663  }
664 
665 
666  /**
667  * Get track kinetic energy.
668  *
669  * \param trk track
670  * \return kinetic energy [GeV]
671  */
672  inline double getKineticEnergy(const Trk& trk)
673  {
674  const double energy = trk.E;
675  const double mass = JPDB::getInstance().getPDG(trk.type).mass;
676 
677  return getKineticEnergy(energy, mass);
678  }
679 
680 
681  /**
682  * Get initial state energy of a neutrino interaction.\n
683  * This method includes baryon number conservation.
684  *
685  * Note: If the given event does not correspond to a neutrino interaction\n
686  * this method returns 0.0 GeV.
687  *
688  * \param evt event
689  * \return energy [GeV]
690  */
691  inline double getE0(const Evt& evt)
692  {
693  double E0 = 0.0;
694 
695  if (!evt.mc_trks.empty() && is_neutrino(evt.mc_trks[0])) {
696 
697  const Trk& neutrino = evt.mc_trks[0];
698 
699  E0 += neutrino.E;
700 
701  for (size_t i = 1; i != evt.mc_trks.size(); ++i) {
702 
703  const Trk& track = evt.mc_trks[i];
704 
705  if (track.type == TRACK_TYPE_NEUTRON ||
706  track.type == TRACK_TYPE_SIGMA_MINUS) {
707  E0 += MASS_NEUTRON;
708  }
709 
710  if (track.type == TRACK_TYPE_PROTON ||
711  track.type == TRACK_TYPE_LAMBDA ||
712  track.type == TRACK_TYPE_SIGMA_PLUS) {
713  E0 += MASS_PROTON;
714  }
715 
716  if (track.type == TRACK_TYPE_ANTINEUTRON) {
717  E0 -= MASS_NEUTRON;
718  }
719 
720  if (track.type == TRACK_TYPE_ANTIPROTON ||
721  track.type == TRACK_TYPE_ANTILAMBDA) {
722  E0 -= MASS_PROTON;
723  }
724  }
725  }
726 
727  return E0;
728  }
729 
730 
731  /**
732  * Get final state energy of a neutrino interaction.\n
733  * This method includes muon energy loss.
734  *
735  * \param evt event
736  * \return energy [GeV]
737  */
738  inline double getE1(const Evt& evt)
739  {
740  double E1 = 0.0;
741 
742  if (!evt.mc_trks.empty() && is_neutrino(evt.mc_trks[0])) {
743 
744  const Trk& neutrino = evt.mc_trks[0];
745 
746  for (size_t i = 1; i != evt.mc_trks.size(); ++i) {
747 
748  const Trk& track = evt.mc_trks[i];
749 
750  if (is_muon(track)) {
751 
752  const double distance = (track.pos - neutrino.pos).len();
753 
754  E1 += JPHYSICS::gWater(track.E, -distance);
755 
756  } else {
757 
758  E1 += track.E;
759  }
760  }
761  }
762 
763  return E1;
764  }
765 
766 
767  /**
768  * Get momentum vector of the initial state of a neutrino interaction.\n
769  * This method assumes neutrino DIS on a stationary nucleus
770  *
771  * \param evt event
772  * \return energy [GeV]
773  */
774  inline Vec getP0(const Evt& evt)
775  {
776  Vec P0(0,0,0);
777 
778  if (!evt.mc_trks.empty() && is_neutrino(evt.mc_trks[0])) {
779 
780  const Trk& neutrino = evt.mc_trks[0];
781 
782  P0 = neutrino.E * neutrino.dir;
783  }
784 
785  return P0;
786  }
787 
788 
789  /**
790  * Get momentum vector of the final state of a neutrino interaction.\n
791  * This method includes muon energy losses.
792  *
793  * Note: If the given event does not correspond to a neutrino interaction\n
794  * a null vector will be returned.
795  *
796  * \param evt event
797  * \return final state momentum vector [GeV]
798  */
799  inline Vec getP1(const Evt& evt)
800  {
801  Vec P1(0,0,0);
802 
803  if (!evt.mc_trks.empty() && is_neutrino(evt.mc_trks[0])) {
804 
805  const Trk& neutrino = evt.mc_trks[0];
806 
807  for (size_t i = 1; i != evt.mc_trks.size(); ++i) {
808 
809  const Trk& track = evt.mc_trks[i];
810 
811  double kineticEnergy = 0.0;
812 
813  if (is_muon(track)) {
814 
815  const double distance = (track.pos - neutrino.pos).len();
816  const double energy = JPHYSICS::gWater(track.E, -distance);
817 
818  kineticEnergy = getKineticEnergy(energy, JPHYSICS::MASS_MUON);
819 
820  } else {
821 
822  kineticEnergy = getKineticEnergy(track);
823  }
824 
825  P1 += kineticEnergy * track.dir;
826  }
827  }
828 
829  return P1;
830  }
831 }
832 
833 #endif
bool from_tau(const Hit &hit)
Test whether given hit was produced by a tau.
Vec getP1(const Evt &evt)
Get momentum vector of the final state of a neutrino interaction.
static const int RESCALED_WEIGHT_INDEX
Index of rescaled weight.
JVertex3D getVertex(const Trk &track)
Get vertex.
static const JRange< double, std::less< double > > DEFAULT_RANGE
Default range.
Definition: JRange.hh:555
bool is_electron(const Trk &track)
Test whether given track is a (anti-)electron.
int count_electrons(const Evt &evt)
Count the number of electrons in a given event.
JHitType_t
Enumeration of hit types based on km3 codes.
Exceptions.
Scattered light from primary shower.
bool is_proton(const Trk &track)
Test whether given track is a (anti-)proton.
const Trk & get_muon(const Evt &evt)
Get first muon from the event tracklist.
Data structure for direction in three dimensions.
Definition: JDirection3D.hh:33
then usage E
Definition: JMuonPostfit.sh:35
JTrack3E getTrack(const Trk &track)
Get track.
Direct light from delta-rays.
double t
track time [ns] (when the particle is at pos )
Definition: Trk.hh:19
double z
Definition: Vec.hh:14
Scattered light from muon.
JTransformation3D getTransformation(const Trk &track)
Get transformation.
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
double getE0(const Evt &evt)
Get initial state energy of a neutrino interaction.
Energy loss of muon.
void add_time_offset(Evt &evt, const double tOff)
Add time offset to mc event; according to current definition, the absolute time of the event is defin...
bool has_neutrino(const Evt &evt)
Test whether given event has an incoming neutrino.
bool is_muon(const Trk &track)
Test whether given track is a (anti-)muon.
bool has_electron(const Evt &evt)
Test whether given event has an electron.
JTOOLS::JRange< double > JTimeRange
Type definition for time range (unit [ns]).
int count_muons(const Evt &evt)
Count the number of muons in a given event.
Auxiliary methods for selection of reconstructed tracks.
Vec dir
track direction
Definition: Trk.hh:18
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Definition: JException.hh:696
double getE1(const Evt &evt)
Get final state energy of a neutrino interaction.
static const double MASS_MUON
muon mass [GeV]
bool has_particleID(const Trk &track, const int type)
Test whether given track corresponds to given particle type.
bool is_charged_lepton(const Trk &track)
Test whether given track is a charged lepton.
Scattered light from Bremsstrahlung.
static const JGeaneWater gWater
Function object for energy loss of muon in sea water.
Definition: JGeane.hh:381
const JParticle & getPDG(const int pdg) const
Get particle corresponding to given PDG code.
Definition: JPDB.hh:234
Direct light from Bremsstrahlung.
double getTime(const Hit &hit)
Get true time of hit.
bool is_baryon(const Trk &track)
Test whether given track is a baryon (is hadron and third digit of PDG code is not zero) ...
double E
Energy [GeV] (either MC truth or reconstructed)
Definition: Trk.hh:20
double y
Definition: Vec.hh:14
3D track with energy.
Definition: JTrack3E.hh:30
Direct light from primary shower.
static const int AASHOWER_RECONSTRUCTION_TYPE
AAnet shower fit reconstruction type.
static const JPDB & getInstance()
Get particle data book.
Definition: JPDB.hh:125
bool is_noise(const Hit &hit)
Verify hit origin.
JTimeRange getTimeRange(const Evt &event)
Get time range (i.e. time between earliest and latest hit) of Monte Carlo event.
static const double MASS_NEUTRON
neutron mass [GeV]
bool is_neutrino(const Trk &track)
Test whether given track is a neutrino.
Axis object.
Definition: JAxis3D.hh:38
bool is_lepton(const Trk &track)
Test whether given track is a lepton.
double getW(const Trk &track, const size_t index, const double value)
Get track information.
double a
hit amplitude (in p.e.)
Definition: Hit.hh:24
double x
Definition: Vec.hh:14
The Vec class is a straightforward 3-d vector, which also works in pyroot.
Definition: Vec.hh:12
double getKineticEnergy(const double E, const double m)
Get kinetic energy of particle with given mass.
JDirection3D getDirection(const Vec &dir)
Get direction.
const Trk & get_electron(const Evt &evt)
Get first electron from the event tracklist.
double getDY() const
Get y direction.
Definition: JVersor3D.hh:106
double getDX() const
Get x direction.
Definition: JVersor3D.hh:95
JAxis3D getAxis(const Trk &track)
Get axis.
JPosition3D getPosition(const Vec &pos)
Get position.
std::vector< double > fitinf
place to store additional fit info, for jgandalf, see JFitParameters.hh
Definition: Trk.hh:32
Physics constants.
double getY() const
Get y position.
Definition: JVector3D.hh:104
bool is_photon(const Trk &track)
Test whether given track is a photon or neutral pion.
bool from_hadron(const Hit &hit)
Test whether given hit was produced by a hadronic shower.
int type
MC: particle type in PDG encoding.
Definition: Trk.hh:24
then for APP in event gandalf start energy
Definition: JMuonMCEvt.sh:44
bool is_pion(const Trk &track)
Test whether given track is a charged pion.
Vec pos
postion of the track at time t [m]
Definition: Trk.hh:17
double mass
mass of particle [GeV]
Definition: JPDB.hh:92
int count_taus(const Evt &evt)
Count the number of taus in a given event.
Definition: Hit.hh:8
const Trk & get_tau(const Evt &evt)
Get first tau from the event tracklist.
Direct light from muon.
then P1
bool has_muon(const Evt &evt)
Test whether given event has a muon.
std::vector< Hit > mc_hits
MC: list of MC truth hits.
Definition: Evt.hh:45
double t
hit time (from tdc+calibration or MC truth)
Definition: Hit.hh:23
Scattered light from delta-rays.
static const double MASS_PROTON
proton mass [GeV]
Definition of particle types.
double getX() const
Get x position.
Definition: JVector3D.hh:94
bool has_tau(const Evt &evt)
Test whether given event has a tau.
int id
offline event identifier
Definition: Evt.hh:21
double getNPE(const Hit &hit)
Get true charge of hit.
Data structure for position in three dimensions.
Definition: JPosition3D.hh:36
bool is_neutron(const Trk &track)
Test whether given track is a (anti-)neutron.
Exception for accessing a value in a collection that is outside of its range.
Definition: JException.hh:162
Exception for accessing an index in a collection that is outside of its range.
Definition: JException.hh:90
Vec getP0(const Evt &evt)
Get momentum vector of the initial state of a neutrino interaction.
const Trk & get_neutrino(const Evt &evt)
Get incoming neutrino.
bool is_hadron(const Trk &track)
Test whether given track is a hadron.
bool is_tau(const Trk &track)
Test whether given track is a (anti-)tau.
bool from_electron(const Hit &hit)
Test whether given hit was produced by an electron.
The Trk class represents a Monte Carlo (MC) particle as well as a reconstructed track/shower.
Definition: Trk.hh:14
int type
particle type or parametrisation used for hit (mc only)
Definition: Hit.hh:28
double getZ() const
Get z position.
Definition: JVector3D.hh:115
double getDZ() const
Get z direction.
Definition: JVersor3D.hh:117
std::vector< Trk > mc_trks
MC: list of MC truth tracks.
Definition: Evt.hh:46
bool is_meson(const Trk &track)
Test whether given track is a meson (is hadron and third digit of PDG code is zero) ...
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
Definition: Evt.hh:19
int count_hadrons(const Evt &evt)
Count the number of hadrons in a given event (not including neutral pions)
bool from_muon(const Hit &hit)
Test whether given hit was produced by a muon.