Jpp test-rotations-new
the software that should make you happy
Loading...
Searching...
No Matches
JShowerMCEvtSmear.cc File Reference

Auxiliary program to store Monte Carlo information from a neutrino or the primary electron in JFIT::JEvt format and smear the values of the position and time around a sphere of 4D - length. More...

Go to the source code of this file.

Functions

int main (int argc, char **argv)
 

Detailed Description

Auxiliary program to store Monte Carlo information from a neutrino or the primary electron in JFIT::JEvt format and smear the values of the position and time around a sphere of 4D - length.

They are equally spaced following https://marcalexa.github.io/superfibonacci/ By default the neutrino is considered, a bool variable allows to select the primary electron instead.

Author
vcarretero

Definition in file JShowerMCEvtSmear.cc.

Function Documentation

◆ main()

int main ( int argc,
char ** argv )

Definition at line 41 of file JShowerMCEvtSmear.cc.

42{
43 using namespace std;
44 using namespace JPP;
45 using namespace KM3NETDAQ;
46
48
51 JLimit_t& numberOfEvents = inputFile.getLimit();
52 int debug;
53 bool take_electron;
54 double length;
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
90 vector<Trk>::const_iterator fermion;
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}
string outputFile
#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:72
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2142
int numberOfPoints
Definition JResultPDF.cc:22
3D track with energy.
Definition JTrack3E.hh:32
Data structure for vector in three dimensions.
Definition JVector3D.hh:36
double getY() const
Get y position.
Definition JVector3D.hh:104
JVector3D & add(const JVector3D &vector)
Add vector.
Definition JVector3D.hh:142
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:1698
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.
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.
KM3NeT DAQ data structures and auxiliaries.
Definition DataQueue.cc:39
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
static counter_type max()
Get maximum counter value.
Definition JLimit.hh:128
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