48 inline const char*
const track_t() {
return "track"; }
49 inline const char*
const shower_t() {
return "shower"; }
60 int main(
int argc,
char **argv)
64 using namespace KM3NETDAQ;
97 JParser<> zap(
"Example program to analyse track fit results from Evt formatted data.");
123 catch(
const exception& error) {
124 FATAL(error.what() << endl);
127 if (useWeights && numberOfEvents != JLimit::max()) {
128 FATAL(
"Cannot apply weighting to limited number of events.");
136 if (detectorFile !=
"") {
150 bool (*has_reconstruction)(
const Evt&);
151 const Trk& (*get_best_fit)(
const Evt&);
156 }
else if (recotype == shower_t()) {
160 FATAL(
"Invalid recotype: " << recotype);
166 const Int_t N_Nhits = 400;
167 const Double_t Nhits_min = 0.0 - 0.5;
168 const Double_t Nhits_max = 1000.0 - 0.5;
170 const Int_t N_Noverlays = 40;
171 const Double_t Noverlays_min = 0.0 - 0.5;
172 const Double_t Noverlays_max = 1000.0 - 0.5;
174 const Int_t N_radius = 201;
175 const Double_t radius_min = -700.0 - 0.5;
176 const Double_t radius_max = +700.0 + 0.5;
178 const Int_t N_height = 18;
179 const Double_t height_min = 20.0;
180 const Double_t height_max = 200.0;
182 const Int_t N_ToTfrac = 20;
183 const Double_t ToTfrac_min = 0.0;
184 const Double_t ToTfrac_max = 1.0;
186 const Int_t N_dt = 1500;
187 const Double_t dt_min = -500.0 - 0.5;
188 const Double_t dt_max = 1000.0 - 0.5;
190 const Int_t N_zenith = 21;
191 const Double_t zenith_min = -1.0;
192 const Double_t zenith_max = +1.0;
194 const Int_t N_energy = 60;
195 const Double_t energy_min = -2.0;
196 const Double_t energy_max = 4.0;
198 const Int_t N_quality = 100;
199 const Double_t quality_min = -200.0;
200 const Double_t quality_max = +800.0;
202 const Int_t N_beta0 = 60;
203 const Double_t beta0_min = -2.0;
204 const Double_t beta0_max = +3.5;
209 string ordinateLabel(useWeights ?
"Rate [Hz]" :
"Number of events");
213 TH1D hs (
"hs",
MAKE_CSTRING(
"; #snapshot hits; " << ordinateLabel),
214 N_Nhits, Nhits_min, Nhits_max);
215 TH1D hn (
"hn",
MAKE_CSTRING(
"; #triggered hits; " << ordinateLabel),
216 N_Nhits, Nhits_min, Nhits_max);
217 TH1D ho (
"ho",
MAKE_CSTRING(
"; #overlays; " << ordinateLabel),
218 N_Noverlays, Noverlays_min, Noverlays_max);
220 TH1D hz (
"hz",
MAKE_CSTRING(
"; cos(#theta); " << ordinateLabel),
221 N_zenith, zenith_min, zenith_max);
222 TH1D he (
"he",
MAKE_CSTRING(
"; log_{10}(E_{reco} / 1 GeV); " << ordinateLabel),
223 N_energy, energy_min, energy_max);
224 TH1D ht (
"ht",
MAKE_CSTRING(
"; #Delta t [ns]; " << ordinateLabel),
225 N_dt, dt_min, dt_max);
227 TH1D hZ (
"hZ",
MAKE_CSTRING(
"; longitudinal ToT-barycenter [m]; " << ordinateLabel),
228 N_height, height_min, height_max);
229 TH1D hT0 (
"hT0",
MAKE_CSTRING(
"; ToT-fraction below earliest triggered hit [ns]; " << ordinateLabel),
230 N_ToTfrac, ToTfrac_min, ToTfrac_max);
231 TH1D hT1 (
"hT1",
MAKE_CSTRING(
"; ToT-fraction above earliest triggered hit [ns]; " << ordinateLabel),
232 N_ToTfrac, ToTfrac_min, ToTfrac_max);
234 TH1D hq (
"hq",
MAKE_CSTRING(
"; quality; " << ordinateLabel),
235 N_quality, quality_min, quality_max);
236 TH1D hb0 (
"hb0",
MAKE_CSTRING(
"; #beta_{0}; " << ordinateLabel),
237 N_beta0, beta0_min, beta0_max);
239 TH2D hzo (
"hzo",
MAKE_CSTRING(
"; cos(#theta); #overlays; " << ordinateLabel),
240 N_zenith, zenith_min, zenith_max,
241 N_Noverlays, Noverlays_min, Noverlays_max);
242 TH2D hze (
"hze",
MAKE_CSTRING(
"; cos(#theta); log_{10}(E_{reco} / 1 GeV); " << ordinateLabel),
243 N_zenith, zenith_min, zenith_max,
244 N_energy, energy_min, energy_max);
245 TH2D heo (
"heo",
MAKE_CSTRING(
"; log_{10}(E_{reco} / 1 GeV); #overlays; " << ordinateLabel),
246 N_energy, energy_min, energy_max,
247 N_Noverlays, Noverlays_min, Noverlays_max);
248 TH2D hzt (
"hzt",
MAKE_CSTRING(
"; cos(#theta); #Delta t [ns]; " << ordinateLabel),
249 N_zenith, zenith_min, zenith_max,
250 N_dt, dt_min, dt_max);
251 TH2D het (
"het",
MAKE_CSTRING(
"; log_{10}(E_{reco} / 1 GeV); #Delta t [ns]; " << ordinateLabel),
252 N_energy, energy_min, energy_max,
253 N_dt, dt_min, dt_max);
254 TH2D hot (
"hot",
MAKE_CSTRING(
"; #overlays; #Delta t [ns]; " << ordinateLabel),
255 N_Noverlays, Noverlays_min, Noverlays_max,
256 N_dt, dt_min, dt_max);
258 TH2D hzq (
"hzq",
MAKE_CSTRING(
"; cos(#theta); quality; " << ordinateLabel),
259 N_zenith, zenith_min, zenith_max,
260 N_quality, quality_min, quality_max);
261 TH2D heq (
"heq",
MAKE_CSTRING(
"; log_{10}(E_{reco} / 1 GeV); quality; " << ordinateLabel),
262 N_energy, energy_min, energy_max,
263 N_quality, quality_min, quality_max);
264 TH2D hoq (
"hoq",
MAKE_CSTRING(
"; #overlays; quality; " << ordinateLabel),
265 N_Noverlays, Noverlays_min, Noverlays_max,
266 N_quality, quality_min, quality_max);
268 TH2D hob0(
"hob0",
MAKE_CSTRING(
"; #overlays; log_{10}(#beta_{0}); " << ordinateLabel),
269 N_Noverlays, Noverlays_min, Noverlays_max,
270 N_beta0, beta0_min, beta0_max);
271 TH2D heb0(
"heb0",
MAKE_CSTRING(
"; log_{10}(E_{reco} / 1 GeV); log_{10}(#beta_{0}); " << ordinateLabel),
272 N_energy, energy_min, energy_max,
273 N_beta0, beta0_min, beta0_max);
274 TH2D hzb0(
"hzb0",
MAKE_CSTRING(
"; cos(#theta); log_{10}(#beta_{0}); " << ordinateLabel),
275 N_zenith, zenith_min, zenith_max,
276 N_beta0, beta0_min, beta0_max);
278 TH2D hxy (
"hxy",
MAKE_CSTRING(
"; x [m]; y [m]; " << ordinateLabel),
279 N_radius, radius_min, radius_max,
280 N_radius, radius_min, radius_max);
285 if (!oscProbTable.empty()) {
292 if (scanners.
setFlux(fluxMaps) == 0) {
293 WARNING(
"No flux function set." << endl);
300 scanner->setLimit(numberOfEvents);
302 while (scanner->hasNext()) {
304 const Evt* evt = scanner->next();
306 const double weight = (useWeights ? scanner->getWeight(*evt) : 1.0);
308 const size_t NsnapshotHits = evt->
hits.size();
309 const size_t NtriggeredHits = count_if(evt->
hits.begin(), evt->
hits.end(),
312 if (!triggeredHitsRange(NtriggeredHits)) {
continue; }
314 hs .Fill(NsnapshotHits, weight);
315 hn .Fill(NtriggeredHits, weight);
318 if (!has_reconstruction(*evt)) {
continue; }
320 const Trk& best_trk = (*get_best_fit)(*evt);
325 if (!energyRange (tb.
getE()) ||
326 !coszenithRange(tb.
getDZ()) ||
334 RIGHT (15) << scanner->getCounter() <<
339 hz .Fill(tb.
getDZ(), weight);
340 he .Fill(logE, weight);
341 hq .Fill(best_trk.
lik, weight);
343 hzo.Fill(tb.
getDZ(), evt->overlays, weight);
344 hze.Fill(tb.
getDZ(), logE, weight);
345 heo.Fill(logE, evt->overlays, weight);
347 hoq.Fill(evt->overlays, best_trk.
lik, weight);
348 hzq.Fill(tb.
getDZ(), best_trk.
lik, weight);
349 heq.Fill(logE, best_trk.
lik, weight);
351 hxy.Fill(tb.
getX(), tb.
getY(), weight);
357 hb0 .Fill(logb0, weight);
358 hob0.Fill(evt->overlays, logb0, weight);
359 heb0.Fill(logE, logb0, weight);
360 hzb0.Fill(tb.
getDZ(), logb0, weight);
363 double ToTbelow = 0.0;
364 double ToTabove = 0.0;
365 double ToTtotal = 0.0;
366 double ToTweightedZ = 0.0;
373 if (router.
hasPMT(hit->pmt_id)) {
377 const double t0 = tb.
getT(pmt);
378 const double t1 =
getTime(*hit);
380 ht .Fill(t1 - t0, weight);
381 hot.Fill(evt->overlays, t1 - t0, weight);
382 hzt.Fill(tb.
getDZ(), t1 - t0, weight);
383 het.Fill(
log10(best_trk.
E), t1 - t0, weight);
387 ToTtotal += hit->tot;
388 ToTweightedZ += hit->tot * hit->pos.z;
390 if (hit->pos.z < firstHit->pos.z) {
391 ToTbelow += hit->tot;
392 }
else if (hit->pos.z >= firstHit->pos.z) {
393 ToTabove += hit->tot;
399 hT0.Fill(ToTbelow / ToTtotal, weight);
400 hT1.Fill(ToTabove / ToTtotal, weight);
401 hZ .Fill(ToTweightedZ / ToTtotal, weight);
Router for direct addressing of PMT data in detector data structure.
Utility class to parse command line options.
JPredicate< JResult_t T::*, JComparison::eq > make_predicate(JResult_t T::*member, const JResult_t value)
Helper method to create predicate for data member.
int main(int argc, char *argv[])
ROOT TTree parameter settings of various packages.
JComparator< JResult_t T::*, JComparison::lt > make_comparator(JResult_t T::*member)
Helper method to create comparator between values of data member.
JTrack3E getTrack(const Trk &track)
Get track.
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
#define gmake_property(A)
macros to convert (template) parameter to JPropertiesElement object
Utility class to parse parameter values.
Data structure for single set of oscillation parameters.
#define MAKE_CSTRING(A)
Make C-string.
bool hasPMT(const JObjectID &id) const
Has PMT.
const JPMT & getPMT(const JPMTAddress &address) const
Get PMT.
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
double getTime(const Hit &hit)
Get true time of hit.
double E
Energy [GeV] (either MC truth or reconstructed)
Data structure for detector geometry and calibration.
Utility class to parse parameter values.
Auxiliary class for defining the range of iterations of objects.
unsigned int overlays
number of overlaying triggered events
double getE() const
Get energy.
const JCylinder3D getMaximumContainmentVolume()
Auxiliary function to retrieve the maximum cylindrical containment volume.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
set_variable E_E log10(E_{fit}/E_{#mu})"
JOscProbFunction< JFunction_t > make_oscProbFunction(const JFunction_t &function)
Auxiliary method for creating an interface to an oscillation probability function.
Data structure for PMT geometry, calibration and status.
std::vector< double > fitinf
place to store additional fit info, see km3net-dataformat/definitions/fitparameters.csv
static const int JGANDALF_BETA0_RAD
KM3NeT Data Definitions v3.1.0-13-g917945e https://git.km3net.de/common/km3net-dataformat.
Auxiliary class for parsing multiparticle fluxes.
double getY() const
Get y position.
Direct access to PMT in detector data structure.
const JPosition3D & getPosition() const
Get position.
Auxiliary methods to convert data members or return values of member methods of a set of objects to a...
double getT(const JVector3D &pos) const
Get arrival time of Cherenkov light at given position.
General purpose messaging.
Scanning of objects from multiple files according a format that follows from the extension of each fi...
Auxiliary base class for list of file names.
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
const char *const track_t
std::vector< filescanner_type >::iterator iterator
Auxiliary class to define a range between two values.
ULong64_t trig
non-zero if the hit is a trigger hit.
double lik
likelihood or lambda value (for aafit, lambda)
Utility class to parse command line options.
double t
hit time (from tdc+calibration or MC truth)
double getX() const
Get x position.
size_t setFlux(const int type, const JFlux &flux)
Set flux function for all MC-files corresponding to a given PDG code.
Auxiliary class for organising Monte Carlo file scanners associated with event weighters.
std::vector< Hit > hits
list of hits
do set_variable DETECTOR_TXT $WORKDIR detector
Auxiliaries for defining the range of iterations of objects.
The Trk class represents a Monte Carlo (MC) particle as well as a reconstructed track/shower.
Auxiliary data structure for floating point format specification.
double getDZ() const
Get z direction.
const array_type< JKey_t > & get_keys(const std::map< JKey_t, JValue_t, JComparator_t, JAllocator_t > &data)
Method to create array of keys of map.
Template definition of a multi-dimensional oscillation probability interpolation table.
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
#define DEBUG(A)
Message macros.