58int main(
int argc,
char **argv)
92 JParser<> zap(
"Example program to analyse track fit results from Evt formatted data.");
114 catch(
const exception& error) {
115 FATAL(error.what() << endl);
118 if (useWeights && numberOfEvents !=
JLimit::max()) {
119 FATAL(
"Cannot apply weighting to limited number of events.");
127 if (detectorFile !=
"") {
141 bool (*has_reconstruction)(
const Evt&);
142 const Trk& (*get_best_fit)(
const Evt&);
144 if (recotype == track_t()) {
147 }
else if (recotype == shower_t()) {
151 FATAL(
"Invalid recotype: " << recotype);
157 const Int_t N_Nhits = 400;
158 const Double_t Nhits_min = 0.0 - 0.5;
159 const Double_t Nhits_max = 1000.0 - 0.5;
161 const Int_t N_Noverlays = 40;
162 const Double_t Noverlays_min = 0.0 - 0.5;
163 const Double_t Noverlays_max = 1000.0 - 0.5;
165 const Int_t N_radius = 201;
166 const Double_t radius_min = -700.0 - 0.5;
167 const Double_t radius_max = +700.0 + 0.5;
169 const Int_t N_height = 18;
170 const Double_t height_min = 20.0;
171 const Double_t height_max = 200.0;
173 const Int_t N_ToTfrac = 20;
174 const Double_t ToTfrac_min = 0.0;
175 const Double_t ToTfrac_max = 1.0;
177 const Int_t N_dt = 1500;
178 const Double_t dt_min = -500.0 - 0.5;
179 const Double_t dt_max = 1000.0 - 0.5;
181 const Int_t N_zenith = 21;
182 const Double_t zenith_min = -1.0;
183 const Double_t zenith_max = +1.0;
185 const Int_t N_energy = 60;
186 const Double_t energy_min = -2.0;
187 const Double_t energy_max = 4.0;
189 const Int_t N_LoE = 100;
190 const Double_t LoE_min = -1.0;
191 const Double_t LoE_max = 8.0;
193 const Int_t N_quality = 100;
194 const Double_t quality_min = -200.0;
195 const Double_t quality_max = +800.0;
197 const Int_t N_beta0 = 60;
198 const Double_t beta0_min = -2.0;
199 const Double_t beta0_max = +3.5;
204 string ordinateLabel(useWeights ?
"Rate [Hz]" :
"Number of events");
208 TH1D hs (
"hs",
MAKE_CSTRING(
"; #snapshot hits; " << ordinateLabel),
209 N_Nhits, Nhits_min, Nhits_max);
210 TH1D hn (
"hn",
MAKE_CSTRING(
"; #triggered hits; " << ordinateLabel),
211 N_Nhits, Nhits_min, Nhits_max);
212 TH1D ho (
"ho",
MAKE_CSTRING(
"; #overlays; " << ordinateLabel),
213 N_Noverlays, Noverlays_min, Noverlays_max);
215 TH1D hz (
"hz",
MAKE_CSTRING(
"; cos(#theta); " << ordinateLabel),
216 N_zenith, zenith_min, zenith_max);
217 TH1D he (
"he",
MAKE_CSTRING(
"; log_{10}(E_{reco} / 1 GeV); " << ordinateLabel),
218 N_energy, energy_min, energy_max);
219 TH1D hLoE(
"hLoE",
MAKE_CSTRING(
"; L/E [km/GeV]; " << ordinateLabel),
220 N_LoE, LoE_min, LoE_max);
221 TH1D ht (
"ht",
MAKE_CSTRING(
"; #Delta t [ns]; " << ordinateLabel),
222 N_dt, dt_min, dt_max);
224 TH1D hZ (
"hZ",
MAKE_CSTRING(
"; longitudinal ToT-barycenter [m]; " << ordinateLabel),
225 N_height, height_min, height_max);
226 TH1D hT0 (
"hT0",
MAKE_CSTRING(
"; ToT-fraction below earliest triggered hit [ns]; " << ordinateLabel),
227 N_ToTfrac, ToTfrac_min, ToTfrac_max);
228 TH1D hT1 (
"hT1",
MAKE_CSTRING(
"; ToT-fraction above earliest triggered hit [ns]; " << ordinateLabel),
229 N_ToTfrac, ToTfrac_min, ToTfrac_max);
231 TH1D hq (
"hq",
MAKE_CSTRING(
"; quality; " << ordinateLabel),
232 N_quality, quality_min, quality_max);
233 TH1D hb0 (
"hb0",
MAKE_CSTRING(
"; #beta_{0}; " << ordinateLabel),
234 N_beta0, beta0_min, beta0_max);
236 TH2D hzo (
"hzo",
MAKE_CSTRING(
"; cos(#theta); #overlays; " << ordinateLabel),
237 N_zenith, zenith_min, zenith_max,
238 N_Noverlays, Noverlays_min, Noverlays_max);
239 TH2D hze (
"hze",
MAKE_CSTRING(
"; cos(#theta); log_{10}(E_{reco} / 1 GeV); " << ordinateLabel),
240 N_zenith, zenith_min, zenith_max,
241 N_energy, energy_min, energy_max);
242 TH2D heo (
"heo",
MAKE_CSTRING(
"; log_{10}(E_{reco} / 1 GeV); #overlays; " << ordinateLabel),
243 N_energy, energy_min, energy_max,
244 N_Noverlays, Noverlays_min, Noverlays_max);
245 TH2D hzt (
"hzt",
MAKE_CSTRING(
"; cos(#theta); #Delta t [ns]; " << ordinateLabel),
246 N_zenith, zenith_min, zenith_max,
247 N_dt, dt_min, dt_max);
248 TH2D het (
"het",
MAKE_CSTRING(
"; log_{10}(E_{reco} / 1 GeV); #Delta t [ns]; " << ordinateLabel),
249 N_energy, energy_min, energy_max,
250 N_dt, dt_min, dt_max);
251 TH2D hot (
"hot",
MAKE_CSTRING(
"; #overlays; #Delta t [ns]; " << ordinateLabel),
252 N_Noverlays, Noverlays_min, Noverlays_max,
253 N_dt, dt_min, dt_max);
255 TH2D hzq (
"hzq",
MAKE_CSTRING(
"; cos(#theta); quality; " << ordinateLabel),
256 N_zenith, zenith_min, zenith_max,
257 N_quality, quality_min, quality_max);
258 TH2D heq (
"heq",
MAKE_CSTRING(
"; log_{10}(E_{reco} / 1 GeV); quality; " << ordinateLabel),
259 N_energy, energy_min, energy_max,
260 N_quality, quality_min, quality_max);
261 TH2D hoq (
"hoq",
MAKE_CSTRING(
"; #overlays; quality; " << ordinateLabel),
262 N_Noverlays, Noverlays_min, Noverlays_max,
263 N_quality, quality_min, quality_max);
265 TH2D hob0(
"hob0",
MAKE_CSTRING(
"; #overlays; log_{10}(#beta_{0}); " << ordinateLabel),
266 N_Noverlays, Noverlays_min, Noverlays_max,
267 N_beta0, beta0_min, beta0_max);
268 TH2D heb0(
"heb0",
MAKE_CSTRING(
"; log_{10}(E_{reco} / 1 GeV); log_{10}(#beta_{0}); " << ordinateLabel),
269 N_energy, energy_min, energy_max,
270 N_beta0, beta0_min, beta0_max);
271 TH2D hzb0(
"hzb0",
MAKE_CSTRING(
"; cos(#theta); log_{10}(#beta_{0}); " << ordinateLabel),
272 N_zenith, zenith_min, zenith_max,
273 N_beta0, beta0_min, beta0_max);
275 TH2D hxy (
"hxy",
MAKE_CSTRING(
"; x [m]; y [m]; " << ordinateLabel),
276 N_radius, radius_min, radius_max,
277 N_radius, radius_min, radius_max);
284 if (scanners.
setFlux(fluxMap) == 0) {
285 WARNING(
"No flux functions set." << endl);
292 scanner->setLimit(numberOfEvents);
294 while (scanner->hasNext()) {
296 const Evt* evt = scanner->next();
298 const double weight = (useWeights ? scanner->getWeight(*evt) : 1.0);
300 const size_t NsnapshotHits = evt->
hits.size();
301 const size_t NtriggeredHits = count_if(evt->
hits.begin(), evt->
hits.end(),
304 if (!triggeredHitsRange(NtriggeredHits)) {
continue; }
306 hs.Fill(NsnapshotHits, weight);
307 hn.Fill(NtriggeredHits, weight);
310 if (!has_reconstruction(*evt)) {
continue; }
312 const Trk& best_trk = (*get_best_fit)(*evt);
315 const double logE = log10(tb.
getE());
317 const double logLoE = log10(L/tb.
getE());
319 if (!energyRange (tb.
getE()) ||
320 !coszenithRange(tb.
getDZ()) ||
328 RIGHT (15) << scanner->getCounter() <<
333 hz .Fill(tb.
getDZ(), weight);
334 he .Fill(logE, weight);
335 hLoE.Fill(logLoE, weight);
336 hq .Fill(best_trk.
lik, weight);
339 hze.Fill(tb.
getDZ(), logE, weight);
340 heo.Fill(logE, evt->
overlays, weight);
343 hzq.Fill(tb.
getDZ(), best_trk.
lik, weight);
344 heq.Fill(logE, best_trk.
lik, weight);
346 hxy.Fill(tb.
getX(), tb.
getY(), weight);
352 hb0 .Fill(logb0, weight);
353 hob0.Fill(evt->
overlays, logb0, weight);
354 heb0.Fill(logE, logb0, weight);
355 hzb0.Fill(tb.
getDZ(), logb0, weight);
358 double ToTbelow = 0.0;
359 double ToTabove = 0.0;
360 double ToTtotal = 0.0;
361 double ToTweightedZ = 0.0;
363 vector<Hit>::const_iterator firstHit = min_element(evt->
hits.cbegin(), evt->
hits.cend(),
366 for (vector<Hit>::const_iterator hit = evt->
hits.cbegin(); hit != evt->
hits.cend(); ++hit) {
368 if (router.
hasPMT(hit->pmt_id)) {
372 const double t0 = tb.
getT(pmt);
373 const double t1 =
getTime(*hit);
375 ht .Fill(t1 - t0, weight);
376 hot.Fill(evt->
overlays, t1 - t0, weight);
377 hzt.Fill(tb.
getDZ(), t1 - t0, weight);
378 het.Fill(log10(best_trk.
E), t1 - t0, weight);
382 ToTtotal += hit->tot;
383 ToTweightedZ += hit->tot * hit->pos.z;
385 if (hit->pos.z < firstHit->pos.z) {
386 ToTbelow += hit->tot;
387 }
else if (hit->pos.z >= firstHit->pos.z) {
388 ToTabove += hit->tot;
394 hT0.Fill(ToTbelow / ToTtotal, weight);
395 hT1.Fill(ToTabove / ToTtotal, weight);
396 hZ .Fill(ToTweightedZ / ToTtotal, weight);