46{
49
51
53 JLimit_t& numberOfEvents = inputFile.getLimit();
55 string detectorFile;
58
59 try {
60
61 JParser<> zap(
"Example program to verify Monte Carlo data.");
62
69
70 zap(argc, argv);
71 }
72 catch(const exception &error) {
73 FATAL(error.what() << endl);
74 }
75
76
78
79 if (detectorFile != "") {
80
81 try {
83 }
86 }
87 }
88
90
93
94 NOTICE(
"Apply detector offset " << offset << endl);
95
97
99
101
102 for ( ;
x < -10.0;
x += 5.0) { X.push_back(x); }
103 for ( ;
x < +20.0;
x += 1.0) { X.push_back(x); }
104 for ( ;
x < +50.0;
x += 2.0) { X.push_back(x); }
105 for ( ;
x < +100.0;
x += 5.0) { X.push_back(x); }
106 for ( ;
x < +250.0;
x += 10.0) { X.push_back(x); }
107 for ( ;
x < +500.0;
x += 25.0) { X.push_back(x); }
108 for ( ;
x < +900.0;
x += 50.0) { X.push_back(x); }
109
111
113
115
117
118 const Evt*
event = inputFile.
next();
119
120 size_t number_of_muons = 0;
121
122 for (vector<Trk>::const_iterator track = event->mc_trks.begin(); track != event->mc_trks.end(); ++track) {
124 number_of_muons += 1;
125 }
126 }
127
128 if (!M(number_of_muons)) {
129 continue;
130 }
131
132 for (vector<Hit>::const_iterator hit = event->mc_hits.begin(); hit != event->mc_hits.end(); ++hit) {
133
134 npe[hit->type] += hit->a;
135
136 const double t1 =
getTime(*hit);
137
138 vector<Trk>::const_iterator trk = find_if(event->mc_trks.begin(),
139 event->mc_trks.end(),
140 make_predicate(&
Trk::id, hit->origin));
141
142 if (trk != event->mc_trks.end()) {
143
144 if (router.hasPMT(hit->pmt_id)) {
145
146 const JPMT& pmt = router.getPMT(hit->pmt_id);
147 double t0 = numeric_limits<double>::max();
148
151
153
156
158
159 } else {
160
161 continue;
162 }
163
164 H1[pmt.
getID()]->Fill(t1 - t0, hit->a);
165
166 H1->Fill(t1 - t0, hit->a);
167 }
168 }
169 }
170 }
172
174
175 const double W = 1.0 / (double) inputFile.
getCounter();
176
178
179 for (auto& i : H1) {
181 }
182 }
183
185
186 cout << "photon-electron statistics" << endl;
187
188 for (const auto& i : npe) {
189 cout << setw(3) << showpos << i.first << ' ' << setw(8) << noshowpos << i.second << endl;
190 }
191 }
192
194
195 out << *H1 << H1;
196
197 out.Write();
198 out.Close();
199}
#define DEBUG(A)
Message macros.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Router for direct addressing of PMT data in detector data structure.
Data structure for PMT geometry, calibration and status.
double getT(const JVector3D &pos) const
Get arrival time of Cherenkov light at given position.
double getT(const JVector3D &pos) const
Get arrival time of Cherenkov light at given position.
int getID() const
Get identifier.
Utility class to parse command line options.
Auxiliary class to manage set of compatible ROOT objects (e.g. histograms) using unique keys.
General purpose class for object reading from a list of file names.
virtual bool hasNext() override
Check availability of next element.
counter_type getCounter() const
Get counter.
virtual const pointer_type & next() override
Get next element.
JTrack3E getTrack(const Trk &track)
Get track.
bool is_electron(const Trk &track)
Test whether given track is a (anti-)electron.
JShower3E getShower(const Trk &shower)
Get shower.
JPosition3D getPosition(const Vec &pos)
Get position.
bool is_muon(const Trk &track)
Test whether given track is a (anti-)muon.
Vec getOffset(const JHead &header)
Get offset.
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.
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
void convertToPDF(TH1 &h1, const std::string &option="NW", const double factor=1.0)
Convert 1D histogram to PDF.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Head getHeader(const JMultipleFileScanner_t &file_list)
Get Monte Carlo header.
const char * getTime()
Get current local time conform ISO-8601 standard.
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
Type definition of range.
Auxiliary class for defining the range of iterations of objects.
static counter_type max()
Get maximum counter value.
The Vec class is a straightforward 3-d vector, which also works in pyroot.