Jpp  18.2.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 
15 
16 #include "JLang/JException.hh"
17 #include "JLang/JPredicate.hh"
18 
21 #include "JGeometry3D/JAxis3D.hh"
22 #include "JGeometry3D/JTrack3D.hh"
23 #include "JGeometry3D/JTrack3E.hh"
24 #include "JGeometry3D/JVertex3D.hh"
26 
27 #include "JPhysics/JConstants.hh"
28 #include "JPhysics/JTimeRange.hh"
29 #include "JPhysics/JGeane.hh"
30 
31 #include "JAAnet/JParticleTypes.hh"
32 #include "JAAnet/JPDB.hh"
33 
34 
35 /**
36  * \file
37  *
38  * Definition of hit and track types and auxiliary methods for handling Monte Carlo data.
39  * \author mdejong
40  */
41 namespace JAANET {}
42 namespace JPP { using namespace JAANET; }
43 
44 /**
45  * Extensions to Evt data format.
46  */
47 namespace JAANET {
48 
52 
60 
64 
65  /**
66  * AAnet shower fit reconstruction type.
67  */
68  static const int AASHOWER_RECONSTRUCTION_TYPE = 101;
69 
70  /**
71  * Enumeration of hit types based on km3 codes.
72  */
73  enum JHitType_t {
74  HIT_TYPE_MUON_DIRECT = +5, //!< Direct light from muon
75  HIT_TYPE_MUON_SCATTERED = -5, //!< Scattered light from muon
76  HIT_TYPE_DELTARAYS_DIRECT = +4, //!< Direct light from delta-rays
77  HIT_TYPE_DELTARAYS_SCATTERED = -4, //!< Scattered light from delta-rays
78  HIT_TYPE_BREMSSTRAHLUNG_DIRECT = +3, //!< Direct light from Bremsstrahlung
79  HIT_TYPE_BREMSSTRAHLUNG_SCATTERED = -3, //!< Scattered light from Bremsstrahlung
80  HIT_TYPE_SHOWER_DIRECT = +99, //!< Direct light from primary shower
81  HIT_TYPE_SHOWER_SCATTERED = -99, //!< Scattered light from primary shower
82  HIT_TYPE_NOISE = -1, //!< Random noise
83  HIT_TYPE_UNKNOWN = 0 //!< Unknown source
84  };
85 
86 
87  /**
88  * Get true time of hit.
89  *
90  * \param hit hit
91  * \return true time of hit
92  */
93  inline double getTime(const Hit& hit)
94  {
95  return hit.t;
96  }
97 
98 
99  /**
100  * Get true charge of hit.
101  *
102  * \param hit hit
103  * \return true charge of hit
104  */
105  inline double getNPE(const Hit& hit)
106  {
107  return hit.a;
108  }
109 
110 
111  /**
112  * Verify hit origin.
113  *
114  * \param hit hit
115  * \return true if noise; else false
116  */
117  inline bool is_noise(const Hit& hit)
118  {
119  return hit.type == HIT_TYPE_NOISE;
120  }
121 
122 
123  /**
124  * Get time range (i.e.\ time between earliest and latest hit) of Monte Carlo event.\n
125  * Note that the global event time in not taken into account.
126  *
127  * \param event event
128  * \return time range
129  */
130  inline JTimeRange getTimeRange(const Evt& event)
131  {
133 
134  for (std::vector<Hit>::const_iterator hit = event.mc_hits.begin(); hit != event.mc_hits.end(); ++hit) {
135  if (!is_noise(*hit)) {
136  time_range.include(getTime(*hit));
137  }
138  }
139 
140  return time_range;
141  }
142 
143 
144  /**
145  * Get time range (i.e.\ time between earliest and latest hit) of Monte Carlo event.\n
146  * The resulting time range is larger than or equal to the given time window.\n
147  * Note that the global event time in not taken into account.
148  *
149  * \param event event
150  * \param T_ns time window [ns]
151  * \return time range
152  */
153  inline JTimeRange getTimeRange(const Evt& event, const JTimeRange& T_ns)
154  {
155  JTimeRange time_range = getTimeRange(event);
156 
157  if (!time_range.is_valid()) {
158  time_range = T_ns;
159  }
160 
161  if (time_range.getLength() < T_ns.getLength()) {
162  time_range.setLength(T_ns.getLength());
163  }
164 
165  return time_range;
166  }
167 
168 
169  /**
170  * Get position.
171  *
172  * \param pos position
173  * \return position
174  */
175  inline JPosition3D getPosition(const Vec& pos)
176  {
177  return JPosition3D(pos.x, pos.y, pos.z);
178  }
179 
180 
181  /**
182  * Get position.
183  *
184  * \param pos position
185  * \return position
186  */
187  inline Vec getPosition(const JPosition3D& pos)
188  {
189  return Vec(pos.getX(), pos.getY(), pos.getZ());
190  }
191 
192 
193  /**
194  * Get position.
195  *
196  * \param track track
197  * \return position
198  */
199  inline JPosition3D getPosition(const Trk& track)
200  {
201  return getPosition(track.pos);
202  }
203 
204 
205  /**
206  * Get direction.
207  *
208  * \param dir direction
209  * \return direction
210  */
211  inline JDirection3D getDirection(const Vec& dir)
212  {
213  return JDirection3D(dir.x, dir.y, dir.z);
214  }
215 
216 
217  /**
218  * Get direction.
219  *
220  * \param dir direction
221  * \return direction
222  */
223  inline Vec getDirection(const JDirection3D& dir)
224  {
225  return Vec(dir.getDX(), dir.getDY(), dir.getDZ());
226  }
227 
228 
229  /**
230  * Get direction.
231  *
232  * \param track track
233  * \return direction
234  */
235  inline JDirection3D getDirection(const Trk& track)
236  {
237  return getDirection(track.dir);
238  }
239 
240 
241  /**
242  * Get axis.
243  *
244  * \param track track
245  * \return axis
246  */
247  inline JAxis3D getAxis(const Trk& track)
248  {
249  return JAxis3D(getPosition(track), getDirection(track));
250  }
251 
252 
253  /**
254  * Get track.
255  *
256  * \param track track
257  * \return track
258  */
259  inline JTrack3E getTrack(const Trk& track)
260  {
261  return JTrack3E(JTrack3D(getAxis(track), track.t), track.E);
262  }
263 
264 
265  /**
266  * Get transformation.
267  *
268  * \param track track
269  * \return transformation
270  */
272  {
273  return JTransformation3D(getPosition(track), getDirection(track));
274  }
275 
276 
277  /**
278  * Get track information.
279  *
280  * \param track track
281  * \param index index
282  * \param value default value
283  * \return actual value
284  */
285  inline double getW(const Trk& track, const size_t index, const double value)
286  {
287  if (index < track.fitinf.size())
288  return track.fitinf[index];
289  else
290  return value;
291  }
292 
293 
294  /**
295  * Test whether given track is a photon or neutral pion.
296  *
297  * \param track track
298  * \return true if photon or neutral pion; else false
299  */
300  inline bool is_photon (const Trk& track) { return ( track.type == TRACK_TYPE_PHOTON ||
301  abs(track.type) == TRACK_TYPE_NEUTRAL_PION); }
302 
303  /**
304  * Test whether given track is a neutrino.
305  *
306  * \param track track
307  * \return true if neutrino; else false
308  */
309  inline bool is_neutrino(const Trk& track) { return (abs(track.type) == TRACK_TYPE_NUE ||
310  abs(track.type) == TRACK_TYPE_NUMU ||
311  abs(track.type) == TRACK_TYPE_NUTAU); }
312 
313  /**
314  * Test whether given track is a (anti-)electron.
315  *
316  * \param track track
317  * \return true if (anti-)electron; else false
318  */
319  inline bool is_electron(const Trk& track) { return (abs(track.type) == TRACK_TYPE_ELECTRON); }
320 
321  /**
322  * Test whether given track is a (anti-)muon.
323  *
324  * \param track track
325  * \return true if (anti-)muon; else false
326  */
327  inline bool is_muon (const Trk& track) { return (abs(track.type) == TRACK_TYPE_MUON); }
328 
329  /**
330  * Test whether given track is a (anti-)tau.
331  *
332  * \param track track
333  * \return true if (anti-)tau; else false
334  */
335  inline bool is_tau (const Trk& track) { return (abs(track.type) == TRACK_TYPE_TAU); }
336 
337  /**
338  * Test whether given track is a charged pion.
339  *
340  * \param track track
341  * \return true if charged pion; else false
342  */
343  inline bool is_pion (const Trk& track) { return (abs(track.type) == TRACK_TYPE_CHARGED_PION_PLUS); }
344 
345  /**
346  * Test whether given track is a (anti-)proton.
347  *
348  * \param track track
349  * \return true if (anti-)proton; else false
350  */
351  inline bool is_proton (const Trk& track) { return (abs(track.type) == TRACK_TYPE_PROTON); }
352 
353  /**
354  * Test whether given track is a (anti-)neutron.
355  *
356  * \param track track
357  * \return true if (anti-)neutron; else false
358  */
359  inline bool is_neutron (const Trk& track) { return (abs(track.type) == TRACK_TYPE_NEUTRON); }
360 
361  /**
362  * Test whether given track is a lepton
363  *
364  * \param track track
365  * \return true if lepton; else fails
366  */
367  inline bool is_lepton (const Trk& track) { return (is_neutrino(track) ||
368  is_electron(track) ||
369  is_muon (track) ||
370  is_tau (track)); }
371 
372  /**
373  * Test whether given track is a charged lepton
374  *
375  * \param track track
376  * \return true if lepton; else fails
377  */
378  inline bool is_charged_lepton (const Trk& track) { return (is_electron(track) ||
379  is_muon (track) ||
380  is_tau (track)); }
381 
382  /**
383  * Test whether given track is a hadron
384  *
385  * \param track track
386  * \return true if hadron; else fails
387  */
388  inline bool is_hadron (const Trk& track) { return (abs(track.type) != TRACK_TYPE_PHOTON &&
389  !is_lepton(track)); }
390 
391  /**
392  * Test whether given track is a meson (is hadron and third digit of PDG code is zero)
393  *
394  * \param track track
395  * \return true if hadron; else fails
396  */
397  inline bool is_meson (const Trk& track) { return (is_hadron(track) &&
398  ((int)(track.type / 1000)) % 10 == 0); }
399 
400  /**
401  * Test whether given track is a baryon (is hadron and third digit of PDG code is not zero)
402  *
403  * \param track track
404  * \return true if hadron; else fails
405  */
406  inline bool is_baryon (const Trk& track) { return (is_hadron(track) &&
407  ((int)(track.type / 1000)) % 10 != 0); }
408 
409  /**
410  * Test whether given track corresponds to given particle type.
411  *
412  * \param track track
413  * \param type particle type
414  * \return true if matches the corresponding particle; else false
415  */
416  inline bool has_particleID (const Trk& track, const int type) { return (track.type == type); }
417 
418  /**
419  * Test whether given track corresponds to an initial state particle.
420  *
421  * \param track track
422  * \return true if track corresponds to an initial state particle; else false
423  */
424  inline bool is_initialstate(const Trk& track)
425  {
426  if (Evt::ROOT_IO_VERSION >= 14) {
427 
428  return (track.status == TRK_ST_PRIMARYNEUTRINO ||
429  track.status == TRK_ST_PRIMARYCOSMIC ||
430  track.status == TRK_ST_MUONBUNDLE ||
431  track.status == TRK_ST_ININUCLEI);
432 
433  } else if (Evt::ROOT_IO_VERSION >= 11) {
434 
435  return (track.mother_id == TRK_MOTHER_UNDEFINED ||
436  track.mother_id == TRK_MOTHER_NONE);
437 
438  } else {
439 
440  return is_neutrino(track) && track.id == 0;
441  }
442  }
443 
444  /**
445  * Test whether given track corresponds to a final state particle.
446  *
447  * \param track track
448  * \return true if track corresponds to a final state particle; else false
449  */
450  inline bool is_finalstate(const Trk& track)
451  {
452  if (Evt::ROOT_IO_VERSION >= 16) {
453  return track.status == TRK_ST_FINALSTATE;
454  } else {
455  return !is_initialstate(track);
456  }
457  }
458 
459  /**
460  * Test whether given event has an incoming neutrino.
461  *
462  * \param evt event
463  * \return true if neutrino; else false
464  */
465  inline bool has_neutrino(const Evt& evt)
466  {
467  if (Evt::ROOT_IO_VERSION >= 14) {
468 
469  std::vector<Trk>::const_iterator i = find_if(evt.mc_trks.cbegin(), evt.mc_trks.cend(),
471  return i != evt.mc_trks.cend();
472 
473  } else {
474 
475  return !evt.mc_trks.empty() && is_neutrino(evt.mc_trks[0]);
476  }
477  }
478 
479  /**
480  * Get incoming neutrino.
481  *
482  * \param evt event
483  * \return neutrino
484  */
485  inline const Trk& get_neutrino(const Evt& evt)
486  {
487  if (Evt::ROOT_IO_VERSION >= 14) {
488 
489  std::vector<Trk>::const_iterator i = find_if(evt.mc_trks.cbegin(), evt.mc_trks.cend(),
491 
492  if (i != evt.mc_trks.cend()) { return *i; }
493 
494  } else if (!evt.mc_trks.empty() && is_neutrino(evt.mc_trks[0])) {
495 
496  return evt.mc_trks[0];
497  }
498 
499  THROW(JIndexOutOfRange, "JAANET::get_neutrino(): Event " << evt.id << " has no neutrino.");
500  }
501 
502  /**
503  * Get vertex.
504  *
505  * \param track track
506  * \return vertex
507  */
508  inline JVertex3D getVertex(const Trk& track)
509  {
510  return JVertex3D(getPosition(track), track.t);
511  }
512 
513  /**
514  * Get event vertex.
515  *
516  * \param event event
517  * \return event vertex
518  */
519  inline JVertex3D getVertex(const Evt& event)
520  {
521  using namespace std;
522  using namespace JPP;
523 
524  if (has_neutrino(event)) {
525 
526  const Trk& neutrino = get_neutrino(event);
527 
528  return getVertex(neutrino);
529 
530  } else if (!event.mc_trks.empty()) {
531 
532  vector<Trk>::const_iterator i = find_if(event.mc_trks.begin(), event.mc_trks.end(), &is_initialstate);
533 
534  if (i != event.mc_trks.end()) {
535  return getVertex(*i);
536  } else {
537  THROW(JValueOutOfRange, "getVertex(): No initial state track found.");
538  }
539 
540  } else if (!event.trks.empty()) { // For reconstructed DAQ events
541 
542  const Trk& track = get_best_reconstructed_track<JPP_RECONSTRUCTION_TYPE>(event);
543 
544  return getVertex(track);
545 
546  } else {
547 
548  THROW(JValueOutOfRange, "getVertex(): Could not find a valid event vertex.");
549  }
550  }
551 
552  /**
553  * Count the number of electrons in a given event
554  *
555  * \param evt event
556  * \return number of electrons
557  */
558  inline int count_electrons(const Evt& evt)
559  {
560  return std::count_if(evt.mc_trks.begin(), evt.mc_trks.end(), is_electron);
561  }
562 
563  /**
564  * Test whether given event has an electron.
565  *
566  * \param evt event
567  * \return true if event has electron; else false
568  */
569  inline bool has_electron(const Evt& evt)
570  {
571  return count_electrons(evt) != 0;
572  }
573 
574  /**
575  * Get first electron from the event tracklist
576  *
577  * \param evt event
578  * \return electron
579  */
580  inline const Trk& get_electron(const Evt& evt)
581  {
582  if (count_electrons(evt) > 0)
583  return *std::find_if(evt.mc_trks.begin(), evt.mc_trks.end(), is_electron);
584  else
585  THROW(JIndexOutOfRange, "This event has no electron.");
586  }
587 
588  /**
589  * Test whether given hit was produced by an electron
590  *
591  * Warning: This method will only work with the output of KM3Sim.
592  * For JSirene or KM3, you should check if <tt>track.id == hit.origin</tt> instead.
593  *
594  * \param hit hit
595  * \return true if hit was produced by an electron; else false
596  */
597  inline bool from_electron(const Hit& hit)
598  {
599  return hit.type == GEANT4_TYPE_ELECTRON || hit.type == GEANT4_TYPE_ANTIELECTRON;
600  }
601 
602  /**
603  * Count the number of muons in a given event
604  *
605  * \param evt event
606  * \return number of muons
607  */
608  inline int count_muons(const Evt& evt)
609  {
610  return std::count_if(evt.mc_trks.begin(), evt.mc_trks.end(), is_muon);
611  }
612 
613  /**
614  * Test whether given event has a muon.
615  *
616  * \param evt event
617  * \return true if event has muon; else false
618  */
619  inline bool has_muon(const Evt& evt)
620  {
621  return count_muons(evt) != 0;
622  }
623 
624  /**
625  * Get first muon from the event tracklist
626  *
627  * \param evt event
628  * \return muon
629  */
630  inline const Trk& get_muon(const Evt& evt)
631  {
632  if (count_muons(evt) > 0)
633  return *std::find_if(evt.mc_trks.begin(), evt.mc_trks.end(), is_muon);
634  else
635  THROW(JIndexOutOfRange, "This event has no muon.");
636  }
637 
638  /**
639  * Test whether given hit was produced by a muon
640  *
641  * Warning: This method will only work with the output of KM3Sim.
642  * For JSirene or KM3, you should check if <tt>track.id == hit.origin</tt> instead.
643  *
644  * \param hit hit
645  * \return true if hit was produced by a muon; else false
646  */
647  inline bool from_muon(const Hit& hit)
648  {
649  return hit.type == GEANT4_TYPE_MUON || hit.type == GEANT4_TYPE_ANTIMUON;
650  }
651 
652  /**
653  * Count the number of taus in a given event
654  *
655  * \param evt event
656  * \return number of taus
657  */
658  inline int count_taus(const Evt& evt)
659  {
660  return std::count_if(evt.mc_trks.begin(), evt.mc_trks.end(), is_tau);
661  }
662 
663  /**
664  * Test whether given event has a tau.
665  *
666  * \param evt event
667  * \return true if event has tau; else false
668  */
669  inline bool has_tau(const Evt& evt)
670  {
671  return count_taus(evt) != 0;
672  }
673 
674  /**
675  * Get first tau from the event tracklist
676  *
677  * \param evt event
678  * \return tau
679  */
680  inline const Trk& get_tau(const Evt& evt)
681  {
682  if (count_taus(evt) > 0)
683  return *std::find_if(evt.mc_trks.begin(), evt.mc_trks.end(), is_tau);
684  else
685  THROW(JIndexOutOfRange, "This event has no tau.");
686  }
687 
688  /**
689  * Test whether given hit was produced by a tau
690  *
691  * Warning: This method will only work with the output of KM3Sim.
692  * For JSirene or KM3, you should check if <tt>track.id == hit.origin</tt> instead.
693  *
694  * \param hit hit
695  * \return true if hit was produced by a tau; else false
696  */
697  inline bool from_tau(const Hit& hit)
698  {
699  return hit.type == GEANT4_TYPE_TAU || hit.type == GEANT4_TYPE_ANTITAU;
700  }
701 
702  /**
703  * Count the number of hadrons in a given event (not including neutral pions)
704  *
705  * \param evt event
706  * \return number of hadrons
707  */
708  inline int count_hadrons(const Evt& evt)
709  {
710  return std::count_if(evt.mc_trks.begin(), evt.mc_trks.end(), is_hadron);
711  }
712 
713  /**
714  * Test whether given hit was produced by a hadronic shower
715  *
716  * Warning: This method will only work with the output of KM3Sim.
717  * For JSirene or KM3, you should check if <tt>track.id == hit.origin</tt> instead.
718  *
719  * \param hit hit
720  * \return true if hit was produced by a hadron; else false
721  */
722  inline bool from_hadron(const Hit& hit)
723  {
724  return (!from_electron(hit) && !from_muon(hit) && !from_tau(hit));
725  }
726 
727 
728  /**
729  * Add time offset to mc event;
730  * according to current definition, the absolute time of the event is defined by the track "t" attribute;
731  * this could change in the future if the global attribute mc_t is assigned to this purpose.
732  *
733  * \param evt event
734  * \param tOff time offset [ns]
735  */
736  inline void add_time_offset(Evt& evt, const double tOff)
737  {
738  for (std::vector<Trk>::iterator p = evt.mc_trks.begin(); p != evt.mc_trks.end(); p++) {
739  p->t += tOff;
740  }
741  }
742 
743 
744  /**
745  * Get kinetic energy of particle with given mass.
746  *
747  * \param E energy [GeV]
748  * \param m mass [GeV]
749  * \return energy [GeV]
750  */
751  inline double getKineticEnergy(const double E, const double m)
752  {
753  if (E > m) {
754  return sqrt((E - m) * (E + m));
755  } else {
756  return 0.0;
757  }
758  }
759 
760 
761  /**
762  * Get track kinetic energy.
763  *
764  * \param trk track
765  * \return kinetic energy [GeV]
766  */
767  inline double getKineticEnergy(const Trk& trk)
768  {
769  const double energy = trk.E;
770  const double mass = JPDB::getInstance().getPDG(trk.type).mass;
771 
772  return getKineticEnergy(energy, mass);
773  }
774 
775 
776  /**
777  * Get initial state energy of a neutrino interaction.\n
778  * This method includes baryon number conservation.
779  *
780  * \param evt event
781  * \return energy [GeV]
782  */
783  inline double getE0(const Evt& evt)
784  {
785  using namespace std;
786 
787  const Trk& neutrino = get_neutrino(evt);
788 
789  double E0 = neutrino.E;
790 
791  for (vector<Trk>::const_iterator track = evt.mc_trks.cbegin(); track != evt.mc_trks.end(); ++track) {
792 
793  if (!is_finalstate(*track)) { continue; }
794 
795  if (track->type == TRACK_TYPE_NEUTRON ||
796  track->type == TRACK_TYPE_SIGMA_MINUS ||
797  track->type == TRACK_TYPE_NEUTRAL_SIGMA) {
798 
799  E0 += MASS_NEUTRON;
800 
801  } else if (track->type == TRACK_TYPE_PROTON ||
802  track->type == TRACK_TYPE_LAMBDA ||
803  track->type == TRACK_TYPE_SIGMA_PLUS) {
804 
805  E0 += MASS_PROTON;
806 
807  } else if (track->type == TRACK_TYPE_ANTINEUTRON) {
808 
809  E0 -= MASS_NEUTRON;
810 
811  } else if (track->type == TRACK_TYPE_ANTIPROTON ||
812  track->type == TRACK_TYPE_ANTILAMBDA) {
813 
814  E0 -= MASS_PROTON;
815  }
816  }
817 
818  return E0;
819  }
820 
821 
822  /**
823  * Get final state energy of a neutrino interaction.\n
824  * This method includes muon energy loss.
825  *
826  * \param evt event
827  * \return energy [GeV]
828  */
829  inline double getE1(const Evt& evt)
830  {
831  using namespace std;
832  using namespace JPP;
833 
834  double E1 = 0.0;
835 
836  const Trk& neutrino = get_neutrino(evt);
837 
838  for (vector<Trk>::const_iterator track = evt.mc_trks.cbegin(); track != evt.mc_trks.cend(); ++track) {
839 
840  if (!is_finalstate(*track)) { continue; }
841 
842  if (is_muon(*track)) {
843 
844  const Trk& mother = ( track->mother_id > TRK_MOTHER_UNDEFINED ?
845  evt.mc_trks[track->mother_id] : neutrino );
846 
847  const double distance = (track->pos - mother.pos).len();
848 
849  E1 += gWater(track->E, -distance);
850 
851  } else {
852 
853  E1 += track->E;
854  }
855  }
856 
857  return E1;
858  }
859 
860 
861  /**
862  * Get momentum vector of the initial state of a neutrino interaction.\n
863  * This method assumes neutrino DIS on a stationary nucleus
864  *
865  * \param evt event
866  * \return energy [GeV]
867  */
868  inline Vec getP0(const Evt& evt)
869  {
870  const Trk& neutrino = get_neutrino(evt);
871 
872  return neutrino.E * neutrino.dir;
873  }
874 
875 
876  /**
877  * Get momentum vector of the final state of a neutrino interaction.\n
878  * This method includes muon energy losses.
879  *
880  * \param evt event
881  * \return final state momentum vector [GeV]
882  */
883  inline Vec getP1(const Evt& evt)
884  {
885  using namespace std;
886  using namespace JPP;
887 
888  Vec P1(0,0,0);
889 
890  const Trk& neutrino = get_neutrino(evt);
891 
892  for (vector<Trk>::const_iterator track = evt.mc_trks.cbegin(); track != evt.mc_trks.cend(); ++track) {
893 
894  if (!is_finalstate(*track)) { continue; }
895 
896  double kineticEnergy = 0.0;
897 
898  if (is_muon(*track)) {
899 
900  const Trk& mother = ( track->mother_id > TRK_MOTHER_UNDEFINED ?
901  evt.mc_trks[track->mother_id] : neutrino );
902 
903  const double distance = (track->pos - mother.pos).len();
904  const double energy = gWater(track->E, -distance);
905 
906  kineticEnergy = getKineticEnergy(energy, MASS_MUON);
907 
908  } else {
909 
910  kineticEnergy = getKineticEnergy(*track);
911  }
912 
913  P1 += kineticEnergy * track->dir;
914  }
915 
916  return P1;
917  }
918 
919 
920  /**
921  * Get final state invariant mass.
922  *
923  * \param event event
924  * \return invariant mass [GeV]
925  */
926  inline double getInvariantMass(const Evt& event)
927  {
928  using namespace std;
929  using namespace JPP;
930 
931  double Minv = 0.0;
932 
933  for (vector<Trk>::const_iterator track = event.mc_trks.cbegin(); track != event.mc_trks.cend(); ++track) {
934 
935  if (!is_finalstate(*track)) { continue; }
936 
937  const double mass = JPDB::getInstance().getPDG(track->type).mass;
938 
939  Minv += mass;
940  }
941 
942  return Minv;
943  }
944 }
945 
946 #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.
JVertex3D getVertex(const Trk &track)
Get vertex.
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.
double getInvariantMass(const Evt &event)
Get final state invariant mass.
JHitType_t
Enumeration of hit types based on km3 codes.
Exceptions.
Scattered light from primary shower.
JPredicate< JResult_t T::*, JComparison::eq > make_predicate(JResult_t T::*member, const JResult_t value)
Helper method to create predicate for data member.
Definition: JPredicate.hh:128
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 $script< input file >[option[primary[working directory]]] nWhere option can be E
Definition: JMuonPostfit.sh:40
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.
static const int TRK_MOTHER_NONE
mother id of a particle if it has no parent
Definition: trkmembers.hh:13
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:712
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.
bool is_finalstate(const Trk &track)
Test whether given track corresponds to a final state particle.
int mother_id
MC id of the parent particle.
Definition: Trk.hh:29
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, see km3net-dataformat/definitions/fitparameters.csv
Definition: Trk.hh:32
Physics constants.
static const int TRK_ST_PRIMARYNEUTRINO
initial state neutrino (&#39;neutrino&#39; tag in evt files from gseagen and genhen).
Definition: trkmembers.hh:16
static const int TRK_ST_MUONBUNDLE
initial state muon bundle (mupage)
Definition: trkmembers.hh:18
static const int TRK_ST_ININUCLEI
Initial state nuclei (gseagen)
Definition: trkmembers.hh:19
static const int TRK_ST_PRIMARYCOSMIC
initial state cosmic ray (&#39;track_primary&#39; tag in evt files from corant).
Definition: trkmembers.hh:17
double getY() const
Get y position.
Definition: JVector3D.hh:104
bool is_initialstate(const Trk &track)
Test whether given track corresponds to an initial state particle.
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.
static const int TRK_MOTHER_UNDEFINED
KM3NeT Data Definitions v3.1.0-13-g917945e https://git.km3net.de/common/km3net-dataformat.
Definition: trkmembers.hh:12
int type
MC: particle type in PDG encoding.
Definition: Trk.hh:24
int id
track identifier
Definition: Trk.hh:16
int status
MC status code, see km3net-dataformat/definitions/trkmembers.csv for values.
Definition: Trk.hh:28
then for APP in event gandalf start energy
Definition: JMuonMCEvt.sh:44
static const int TRK_ST_FINALSTATE
for MC: the particle must be processed by detector simulation (&#39;track_in&#39; tag in evt files)...
Definition: trkmembers.hh:15
bool is_pion(const Trk &track)
Test whether given track is a charged pion.
Vec pos
postion [m] of the track at time t
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
std::vector< Trk > trks
list of reconstructed tracks (can be several because of prefits,showers, etc).
Definition: Evt.hh:39
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:48
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:22
double getNPE(const Hit &hit)
Get true charge of hit.
Data structure for position in three dimensions.
Definition: JPosition3D.hh:36
static int ROOT_IO_VERSION
Streamer version as obtained from ROOT file.
Definition: Evt.hh:274
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:178
static JRange< double, std::less< double > > DEFAULT_RANGE()
Default range.
Definition: JRange.hh:555
Exception for accessing an index in a collection that is outside of its range.
Definition: JException.hh:106
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:49
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:20
JTOOLS::JRange< double > JTimeRange
Type definition for time range (unit [s]).
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.