Jpp test-rotations-new
the software that should make you happy
Loading...
Searching...
No Matches
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"
13#include "JDAQ/JDAQEventIO.hh"
20#include "JSupport/JSupport.hh"
21#include "JSupport/JMeta.hh"
22#include "JPhysics/JGeanz.hh"
23
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 */
41int main(int argc, char **argv)
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}
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:72
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:2142
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
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.
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.
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
Auxiliary include file for time conversion between DAQ/trigger hit and Monte Carlo hit.