46 inline const char*
const track_t() {
return "track"; }
47 inline const char*
const shower_t() {
return "shower"; }
58 int main(
int argc,
char **argv)
94 JParser<> zap(
"Example program to analyse track fit results from Evt formatted data.");
120 catch(
const exception& error) {
121 FATAL(error.what() << endl);
124 if (useWeights && numberOfEvents != JLimit::max()) {
125 FATAL(
"Cannot apply weighting to limited number of events.");
133 if (detectorFile !=
"") {
147 if (!oscProbTable.empty()) {
151 interpolator.
load(oscProbTable.c_str());
152 interpolator.set (oscParameters);
160 if (scanners.
setFlux(fluxMap) == 0) {
161 WARNING(
"No flux functions set." << endl);
167 bool (*has_reconstruction)(
const Evt&);
168 const Trk& (*get_best_fit)(
const Evt&);
173 }
else if (recotype == shower_t()) {
177 FATAL(
"Invalid recotype: " << recotype);
183 const Int_t N_Nhits = 400;
184 const Double_t Nhits_min = 0.0 - 0.5;
185 const Double_t Nhits_max = 1000.0 - 0.5;
187 const Int_t N_Noverlays = 40;
188 const Double_t Noverlays_min = 0.0 - 0.5;
189 const Double_t Noverlays_max = 1000.0 - 0.5;
191 const Int_t N_radius = 201;
192 const Double_t radius_min = -700.0 - 0.5;
193 const Double_t radius_max = +700.0 + 0.5;
195 const Int_t N_height = 18;
196 const Double_t height_min = 20.0;
197 const Double_t height_max = 200.0;
199 const Int_t N_ToTfrac = 20;
200 const Double_t ToTfrac_min = 0.0;
201 const Double_t ToTfrac_max = 1.0;
203 const Int_t N_dt = 1500;
204 const Double_t dt_min = -500.0 - 0.5;
205 const Double_t dt_max = 1000.0 - 0.5;
207 const Int_t N_zenith = 21;
208 const Double_t zenith_min = -1.0;
209 const Double_t zenith_max = +1.0;
211 const Int_t N_energy = 60;
212 const Double_t energy_min = -2.0;
213 const Double_t energy_max = 4.0;
215 const Int_t N_quality = 100;
216 const Double_t quality_min = -200.0;
217 const Double_t quality_max = +800.0;
219 const Int_t N_beta0 = 60;
220 const Double_t beta0_min = -2.0;
221 const Double_t beta0_max = +3.5;
226 string ordinateLabel(useWeights ?
"Rate [Hz]" :
"Number of events");
230 TH1D hs (
"hs",
MAKE_CSTRING(
"; #snapshot hits; " << ordinateLabel),
231 N_Nhits, Nhits_min, Nhits_max);
232 TH1D hn (
"hn",
MAKE_CSTRING(
"; #triggered hits; " << ordinateLabel),
233 N_Nhits, Nhits_min, Nhits_max);
234 TH1D ho (
"ho",
MAKE_CSTRING(
"; #overlays; " << ordinateLabel),
235 N_Noverlays, Noverlays_min, Noverlays_max);
237 TH1D hz (
"hz",
MAKE_CSTRING(
"; cos(#theta); " << ordinateLabel),
238 N_zenith, zenith_min, zenith_max);
239 TH1D he (
"he",
MAKE_CSTRING(
"; log_{10}(E_{reco} / 1 GeV); " << ordinateLabel),
240 N_energy, energy_min, energy_max);
241 TH1D ht (
"ht",
MAKE_CSTRING(
"; #Delta t [ns]; " << ordinateLabel),
242 N_dt, dt_min, dt_max);
244 TH1D hZ (
"hZ",
MAKE_CSTRING(
"; longitudinal ToT-barycenter [m]; " << ordinateLabel),
245 N_height, height_min, height_max);
246 TH1D hT0 (
"hT0",
MAKE_CSTRING(
"; ToT-fraction below earliest triggered hit [ns]; " << ordinateLabel),
247 N_ToTfrac, ToTfrac_min, ToTfrac_max);
248 TH1D hT1 (
"hT1",
MAKE_CSTRING(
"; ToT-fraction above earliest triggered hit [ns]; " << ordinateLabel),
249 N_ToTfrac, ToTfrac_min, ToTfrac_max);
251 TH1D hq (
"hq",
MAKE_CSTRING(
"; quality; " << ordinateLabel),
252 N_quality, quality_min, quality_max);
253 TH1D hb0 (
"hb0",
MAKE_CSTRING(
"; #beta_{0}; " << ordinateLabel),
254 N_beta0, beta0_min, beta0_max);
256 TH2D hzo (
"hzo",
MAKE_CSTRING(
"; cos(#theta); #overlays; " << ordinateLabel),
257 N_zenith, zenith_min, zenith_max,
258 N_Noverlays, Noverlays_min, Noverlays_max);
259 TH2D hze (
"hze",
MAKE_CSTRING(
"; cos(#theta); log_{10}(E_{reco} / 1 GeV); " << ordinateLabel),
260 N_zenith, zenith_min, zenith_max,
261 N_energy, energy_min, energy_max);
262 TH2D heo (
"heo",
MAKE_CSTRING(
"; log_{10}(E_{reco} / 1 GeV); #overlays; " << ordinateLabel),
263 N_energy, energy_min, energy_max,
264 N_Noverlays, Noverlays_min, Noverlays_max);
265 TH2D hzt (
"hzt",
MAKE_CSTRING(
"; cos(#theta); #Delta t [ns]; " << ordinateLabel),
266 N_zenith, zenith_min, zenith_max,
267 N_dt, dt_min, dt_max);
268 TH2D het (
"het",
MAKE_CSTRING(
"; log_{10}(E_{reco} / 1 GeV); #Delta t [ns]; " << ordinateLabel),
269 N_energy, energy_min, energy_max,
270 N_dt, dt_min, dt_max);
271 TH2D hot (
"hot",
MAKE_CSTRING(
"; #overlays; #Delta t [ns]; " << ordinateLabel),
272 N_Noverlays, Noverlays_min, Noverlays_max,
273 N_dt, dt_min, dt_max);
275 TH2D hzq (
"hzq",
MAKE_CSTRING(
"; cos(#theta); quality; " << ordinateLabel),
276 N_zenith, zenith_min, zenith_max,
277 N_quality, quality_min, quality_max);
278 TH2D heq (
"heq",
MAKE_CSTRING(
"; log_{10}(E_{reco} / 1 GeV); quality; " << ordinateLabel),
279 N_energy, energy_min, energy_max,
280 N_quality, quality_min, quality_max);
281 TH2D hoq (
"hoq",
MAKE_CSTRING(
"; #overlays; quality; " << ordinateLabel),
282 N_Noverlays, Noverlays_min, Noverlays_max,
283 N_quality, quality_min, quality_max);
285 TH2D hob0(
"hob0",
MAKE_CSTRING(
"; #overlays; log_{10}(#beta_{0}); " << ordinateLabel),
286 N_Noverlays, Noverlays_min, Noverlays_max,
287 N_beta0, beta0_min, beta0_max);
288 TH2D heb0(
"heb0",
MAKE_CSTRING(
"; log_{10}(E_{reco} / 1 GeV); log_{10}(#beta_{0}); " << ordinateLabel),
289 N_energy, energy_min, energy_max,
290 N_beta0, beta0_min, beta0_max);
291 TH2D hzb0(
"hzb0",
MAKE_CSTRING(
"; cos(#theta); log_{10}(#beta_{0}); " << ordinateLabel),
292 N_zenith, zenith_min, zenith_max,
293 N_beta0, beta0_min, beta0_max);
295 TH2D hxy (
"hxy",
MAKE_CSTRING(
"; x [m]; y [m]; " << ordinateLabel),
296 N_radius, radius_min, radius_max,
297 N_radius, radius_min, radius_max);
306 scanner->setLimit(numberOfEvents);
308 while (scanner->hasNext()) {
310 const Evt* evt = scanner->next();
312 const double weight = (useWeights ? scanner->getWeight(*evt) : 1.0);
314 const size_t NsnapshotHits = evt->
hits.size();
315 const size_t NtriggeredHits = count_if(evt->
hits.begin(), evt->
hits.end(),
318 if (!triggeredHitsRange(NtriggeredHits)) {
continue; }
320 hs .Fill(NsnapshotHits, weight);
321 hn .Fill(NtriggeredHits, weight);
324 if (!has_reconstruction(*evt)) {
continue; }
326 const Trk& best_trk = (*get_best_fit)(*evt);
329 const double logE = log10(tb.
getE());
331 if (!energyRange (tb.
getE()) ||
332 !coszenithRange(tb.
getDZ()) ||
340 RIGHT (15) << scanner->getCounter() <<
345 hz .Fill(tb.
getDZ(), weight);
346 he .Fill(logE, weight);
347 hq .Fill(best_trk.
lik, weight);
350 hze.Fill(tb.
getDZ(), logE, weight);
351 heo.Fill(logE, evt->
overlays, weight);
354 hzq.Fill(tb.
getDZ(), best_trk.
lik, weight);
355 heq.Fill(logE, best_trk.
lik, weight);
357 hxy.Fill(tb.
getX(), tb.
getY(), weight);
363 hb0 .Fill(logb0, weight);
364 hob0.Fill(evt->
overlays, logb0, weight);
365 heb0.Fill(logE, logb0, weight);
366 hzb0.Fill(tb.
getDZ(), logb0, weight);
369 double ToTbelow = 0.0;
370 double ToTabove = 0.0;
371 double ToTtotal = 0.0;
372 double ToTweightedZ = 0.0;
379 if (router.
hasPMT(hit->pmt_id)) {
383 const double t0 = tb.
getT(pmt);
384 const double t1 =
getTime(*hit);
386 ht .Fill(t1 - t0, weight);
387 hot.Fill(evt->
overlays, t1 - t0, weight);
388 hzt.Fill(tb.
getDZ(), t1 - t0, weight);
389 het.Fill(log10(best_trk.
E), t1 - t0, weight);
393 ToTtotal += hit->tot;
394 ToTweightedZ += hit->tot * hit->pos.z;
396 if (hit->pos.z < firstHit->pos.z) {
397 ToTbelow += hit->tot;
398 }
else if (hit->pos.z >= firstHit->pos.z) {
399 ToTabove += hit->tot;
405 hT0.Fill(ToTbelow / ToTtotal, weight);
406 hT1.Fill(ToTabove / ToTtotal, weight);
407 hZ .Fill(ToTweightedZ / ToTtotal, weight);
int main(int argc, char **argv)
Data structure for detector geometry and calibration.
Auxiliaries for defining the range of iterations of objects.
General purpose messaging.
#define DEBUG(A)
Message macros.
Scanning of objects from multiple files according a format that follows from the extension of each fi...
Direct access to PMT in detector data structure.
Utility class to parse command line options.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
#define MAKE_CSTRING(A)
Make C-string.
Utility class to parse parameter values.
#define gmake_property(A)
macros to convert (template) parameter to JPropertiesElement object
Auxiliary class to define a range between two values.
ROOT TTree parameter settings of various packages.
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
Template specialisation for a map between event categories and flux factors.
Router for direct addressing of PMT data in detector data structure.
bool hasPMT(const JObjectID &id) const
Has PMT.
const JPMT & getPMT(const JPMTAddress &address) const
Get PMT.
Data structure for PMT geometry, calibration and status.
Utility class to parse parameter values.
bool is_inside(const JVector3D &pos) const
Check whether given point is inside cylinder.
const JPosition3D & getPosition() const
Get position.
double getT(const JVector3D &pos) const
Get arrival time of Cherenkov light at given position.
double getE() const
Get energy.
double getY() const
Get y position.
double getX() const
Get x position.
double getDZ() const
Get z direction.
Data structure for single set of oscillation parameters.
Template definition of a multi-dimensional oscillation probability interpolation table.
void load(const char *const fileName)
Load oscillation probability table from file.
Utility class to parse command line options.
static const int JGANDALF_BETA0_RAD
KM3NeT Data Definitions v3.4.0-8-ge14cb17 https://git.km3net.de/common/km3net-dataformat.
JTrack3E getTrack(const Trk &track)
Get track.
double getTime(const Hit &hit)
Get true time of hit.
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
JComparator< JResult_t T::*, JComparison::lt > make_comparator(JResult_t T::*member)
Helper method to create comparator between values of data member.
JPredicate< JResult_t T::*, JComparison::eq > make_predicate(JResult_t T::*member, const JResult_t value)
Helper method to create predicate for data member.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
const JCylinder3D getMaximumContainmentVolume()
Auxiliary function to retrieve the maximum cylindrical containment volume.
KM3NeT DAQ data structures and auxiliaries.
const char *const track_t
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
std::vector< Hit > hits
list of hits
unsigned int overlays
number of overlaying triggered events
ULong64_t trig
non-zero if the hit is a trigger hit.
double t
hit time (from tdc+calibration or MC truth)
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Auxiliary class for organising Monte Carlo file scanners associated with event weighters.
size_t setFlux(const JEvtCategoryHelper &category, const JFluxHelper &flux)
Set flux function for all MC-files corresponding to a given event category.
std::vector< filescanner_type >::iterator iterator
Auxiliary class for defining the range of iterations of objects.
Auxiliary base class for list of file names.
Auxiliary data structure for floating point format specification.
The Trk class represents a Monte Carlo (MC) particle as well as a reconstructed track/shower.
std::vector< double > fitinf
place to store additional fit info, see km3net-dataformat/definitions/fitparameters....
double E
Energy [GeV] (either MC truth or reconstructed)
double lik
likelihood or lambda value (for aafit, lambda)