Jpp  19.1.0-rc.1
the software that should make you happy
JShowerMCEvtSmear.cc
Go to the documentation of this file.
1 #include <string>
2 #include <iostream>
3 #include <iomanip>
4 
9 
10 #include "JAAnet/JHead.hh"
11 #include "JAAnet/JHeadToolkit.hh"
12 #include "JAAnet/JAAnetToolkit.hh"
13 #include "JDAQ/JDAQEventIO.hh"
14 #include "JDAQ/JDAQTimesliceIO.hh"
20 #include "JSupport/JSupport.hh"
21 #include "JSupport/JMeta.hh"
22 #include "JPhysics/JGeanz.hh"
23 
24 #include "JReconstruction/JEvt.hh"
26 
27 #include "Jeep/JParser.hh"
28 #include "Jeep/JMessage.hh"
29 
30 #include <TRandom3.h>
31 /**
32  * \file
33  *
34  * Auxiliary program to store Monte Carlo information from a neutrino or the primary electron in JFIT::JEvt format and
35  * smear the values of the position and time around a sphere of 4D - length. They are equally spaced following
36  * https://marcalexa.github.io/superfibonacci/
37  * By default the neutrino is considered, a bool variable allows to select the primary electron instead.
38  *
39  * \author vcarretero
40  */
41 int main(int argc, char **argv)
42 {
43  using namespace std;
44  using namespace JPP;
45  using namespace KM3NETDAQ;
46 
48 
49  JTriggeredFileScanner<> inputFile;
51  JLimit_t& numberOfEvents = inputFile.getLimit();
52  int debug;
53  bool take_electron;
54  double length;
55  int numberOfPoints;
56  try {
57 
58  JParser<> zap("Auxiliary program to store Monte Carlo true shower in format for subsequent fits.");
59 
60  zap['f'] = make_field(inputFile, "output of JTriggerEfficiency[.sh]");
61  zap['o'] = make_field(outputFile) = "mcevt.root";
62  zap['n'] = make_field(numberOfEvents) = JLimit::max();
63  zap['d'] = make_field(debug) = 2;
64  zap['e'] = make_field(take_electron);
65  zap['r'] = make_field(length) = 1;
66  zap['N'] = make_field(numberOfPoints) = 100;
67  zap(argc, argv);
68  }
69  catch(const exception& error) {
70  FATAL(error.what() << endl);
71  }
72 
73 
74  const Vec offset = getOffset(getHeader(inputFile));
75 
76  outputFile.open();
77 
78  outputFile.put(JMeta(argc, argv));
79 
80  while (inputFile.hasNext()) {
81  STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);
82 
84 
85  const JDAQEvent* tev = ps;
86  const Evt* event = ps;
87  JEvt out;
88  const time_converter converter(*event, *tev);
89 
91 
92  if(!take_electron) fermion = find_if(event->mc_trks.begin(), event->mc_trks.end(), is_neutrino);
93  else fermion = find_if(event->mc_trks.begin(), event->mc_trks.end(), is_electron);
94 
95  if (fermion != event->mc_trks.end()) {
96 
97  // if fermion = neutrino => fermion position = neutrino interaction vertex
98  double time = converter.putTime();
99  JVector3D position = getPosition(*fermion);
100 
101  // if fermion = electron => fermion position and time = EM shower brightest point
102  if(take_electron){
103  JVector3D electron_dir = getDirection(*fermion);
104  double shower_elongation = geanz.getMaximum(fermion->E);
105  electron_dir *= shower_elongation;
106  position.add(electron_dir);
107  time += (shower_elongation * getInverseSpeedOfLight());
108  }
109  TRandom3 rand;
110 
111  double dn = 1.0 / (double)numberOfPoints;
112  double mc0 = 1.0 / sqrt(2.0);
113  double mc1 = 1.0 / 1.533751168755204288118041;
114 
115  for (int i = 0; i < numberOfPoints; i++)
116  {
117  double s = (double)i+0.5;
118  double ab = 2.0 * M_PI * s;
119  double alpha = ab * mc0;
120  double beta = ab * mc1;
121  s *= dn;
122  double r = length * sqrt(s);
123  double R = length * sqrt(1.0-s);
124  double x = r*sin(alpha);
125  double y = r*cos(alpha);
126  double z = R*sin(beta);
127  double ct = R*cos(beta);
128 
129  JVector3D positionS(position.getX() + x,position.getY() + y,position.getZ() + z);
130  JTrack3D td(positionS, getDirection(*fermion), time + ct * getInverseSpeedOfLight() * getIndexOfRefraction());
131 
132  JTrack3E ta(td, fermion->E);
133 
134  ta.add(getPosition(offset));
135  out.push_back(getFit(JMCEVT, ta, 0, 0.0, ta.getE()));
136  }
137  }
138 
139  outputFile.put(out);
140  }
141 
142  STATUS(endl);
143 
145 
146  io >> outputFile;
147 
148  outputFile.close();
149 }
Definition of hit and track types and auxiliary methods for handling Monte Carlo data.
string outputFile
Recording of objects on file according a format that follows from the file name extension.
Longitudinal emission profile EM-shower.
General purpose messaging.
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
#define STATUS(A)
Definition: JMessage.hh:63
#define FATAL(A)
Definition: JMessage.hh:67
int debug
debug level
Definition: JSirene.cc:69
ROOT I/O of application specific meta data.
Utility class to parse command line options.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:2158
int numberOfPoints
Definition: JResultPDF.cc:22
int main(int argc, char **argv)
ROOT TTree parameter settings of various packages.
Synchronously read DAQ events and Monte Carlo events (and optionally other events).
JTime & add(const JTime &value)
Addition operator.
3D track with energy.
Definition: JTrack3E.hh:32
double getE() const
Get energy.
Definition: JTrack3E.hh:93
Data structure for vector in three dimensions.
Definition: JVector3D.hh:36
JVector3D & add(const JVector3D &vector)
Add vector.
Definition: JVector3D.hh:142
double getY() const
Get y position.
Definition: JVector3D.hh:104
double getZ() const
Get z position.
Definition: JVector3D.hh:115
double getX() const
Get x position.
Definition: JVector3D.hh:94
Utility class to parse command line options.
Definition: JParser.hh:1714
double getMaximum(const double E) const
Get depth of shower maximum.
Definition: JGeanz.hh:162
Object writing to file.
General purpose class for object reading from a list of file names.
counter_type getCounter() const
Get counter.
Auxiliary class to convert DAQ hit time to/from Monte Carlo hit time.
double putTime() const
Get Monte Carlo time minus DAQ/trigger time.
static const int JMCEVT
JDirection3D getDirection(const Vec &dir)
Get direction.
bool is_electron(const Trk &track)
Test whether given track is a (anti-)electron.
bool is_neutrino(const Trk &track)
Test whether given track is a neutrino.
JPosition3D getPosition(const Vec &pos)
Get position.
Vec getOffset(const JHead &header)
Get offset.
JFit getFit(const int id, const JMODEL::JString &string)
Get fit parameters of string.
static const JGeanz geanz(1.85, 0.62, 0.54)
Function object for longitudinal EM-shower profile.
double getIndexOfRefraction()
Get average index of refraction of water corresponding to group velocity.
const double getInverseSpeedOfLight()
Get inverse speed of light.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Head getHeader(const JMultipleFileScanner_t &file_list)
Get Monte Carlo header.
data_type r[M+1]
Definition: JPolint.hh:868
KM3NeT DAQ data structures and auxiliaries.
Definition: DataQueue.cc:39
Definition: JSTDTypes.hh:14
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
Definition: Evt.hh:21
Acoustic event fit.
General purpose class for multiple pointers.
Type list.
Definition: JTypeList.hh:23
Auxiliary class for defining the range of iterations of objects.
Definition: JLimit.hh:45
Auxiliary class for ROOT I/O of application specific meta data.
Definition: JMeta.hh:72
Auxiliary class to synchronously read DAQ events and Monte Carlo events (and optionally other events)...
virtual bool hasNext() override
Check availability of next element.
virtual const multi_pointer_type & next() override
Get next element.
The Vec class is a straightforward 3-d vector, which also works in pyroot.
Definition: Vec.hh:13
Auxiliary include file for time conversion between DAQ/trigger hit and Monte Carlo hit.