54 inline double getMuonWeight(
const double E,
71 inline double getMuonWeight(
const double E,
79 return R *
exp(+R/l_att) * getMuonWeight(E, npe);
90 inline double getNeutrinoWeight(
const double E,
105 inline double getNeutrinoWeight(
const double E,
109 static const double l_att = 45.0;
111 return D*D *
exp(+D/l_att) * getNeutrinoWeight(E, npe);
126 return first.second > second.second;
137 inline TH1D* makeSortedH1(TH1D*
in,
const TString&
name)
146 for (Int_t ix = 1; ix <= in->GetXaxis()->GetNbins(); ++ix) {
148 const Int_t
x = (Int_t) in->GetXaxis()->GetBinCenter(ix);
149 const Double_t y = in->GetBinContent(ix);
152 data.push_back(make_pair(x, y));
158 sort(data.begin(), data.end(), compare);
162 TH1D* out =
new TH1D(name, NULL, data.size(), -0.5, data.size() - 0.5);
164 for (Int_t i = 0; i < (Int_t) data.size(); i++) {
166 const Int_t ix = i + 1;
168 out->GetXaxis()->SetBinLabel(ix,
MAKE_CSTRING(fixed << setprecision(0) << data[i].
first));
169 out->SetBinContent(ix, data[i].second);
183 int main(
int argc,
char **argv)
201 JParser<> zap(
"Example program to verify Monte Carlo data.");
204 zap[
'n'] =
make_field(numberOfEvents) = JLimit::max();
212 catch(
const exception &error) {
213 FATAL(error.what() << endl);
225 if (detectorFile !=
"") {
246 NOTICE(
"Apply detector offset " << center << endl);
257 NOTICE(
"JDomino: Instrumented volume dimensions (x, y, r, zmin, zmax): " << inst_vol << endl);
271 if (header.is_valid(&JHead::cut_nu)) {
272 xmin =
log10(header.cut_nu.E.getLowerLimit());
273 xmax =
log10(header.cut_nu.E.getUpperLimit());
274 }
else if (header.is_valid(&JHead::cut_in)) {
275 xmin =
log10(header.cut_in.E.getLowerLimit());
276 xmax =
log10(header.cut_in.E.getUpperLimit());
278 const int nx = (int) ((xmax - xmin) / 0.1);
280 const Double_t
T[] = { -50.0, -20.0, -10.0, -5.0, -2.0, 0.0, +2.0, +5.0, +10.0, +15.0, +20.0, +30.0, +40.0, +50.0,
281 +75.0, +100.0, +150.0, +200.0, +250.0, +300.0, +400.0, +500.0, 1000.0 };
287 TH1D hits(
"hits", NULL, 100, 0.0, 8.0);
288 TH1D trks(
"trks", NULL, 10001, -5000.5, 5000.5);
289 TH1D job (
"job" , NULL, 10001, -5000.5, 5000.5);
291 TProfile hits_per_E_in (
"hits_per_E_in",
"average number of hits per E_{#nu} inside instrumented volume", nx, xmin, xmax);
292 TProfile hits_per_E_out(
"hits_per_E_out",
"average number of hits per E_{#nu} outside instrumented volume", nx, xmin, xmax);
296 JManager_t npe_per_type(
new TH1D(
"npe[%]", NULL, 5000, 0.5, 5000.5),
'%', ios::fmtflags(ios::showpos));
298 TH1D npe_per_pmt (
"pmt", NULL, 100001, -0.5, 100000.5);
305 TH2D nuExD(
"nuExD", NULL, nx, xmin, xmax, 100, 0.0, 1000.0);
306 TH2D* nuExD_tmp = (TH2D*) nuExD.Clone(
"nuExD.tmp");
308 TH2D nuExc(
"nuExc", NULL, nx, xmin, xmax, 100, -1.0, +1.0);
309 TH2D* nuExc_tmp = (TH2D*) nuExc.Clone(
"nuExc.tmp");
311 TH2D nuDxc(
"nuDxc", NULL, 50, 0.0, 1000.0, 100, -1.0, +1.0);
312 TH2D* nuDxc_tmp = (TH2D*) nuDxc.Clone(
"nuDxc.tmp");
314 TH2D nuDxU(
"nuDxU", NULL, 50, 0.0, 1000.0, 100, -1.0, +1.0);
315 TH2D* nuDxU_tmp = (TH2D*) nuDxU.Clone(
"nuDxU.tmp");
317 TH2D nuDxT(
"nuDxT", NULL, 50, 0.0, 1000.0,
getSize(
T) - 1,
T);
318 TH2D* nuDxT_tmp = (TH2D*) nuDxT.Clone(
"nuDxT.tmp");
331 TH2D muExR(
"muExR", NULL, nx, xmin, xmax, 30, 0.0, 300.0);
332 TH2D* muExR_tmp = (TH2D*) muExR.Clone(
"muExR.tmp");
334 TH2D muRxU(
"muRxU", NULL, 15, 0.0, 300.0, 100, -1.0, +1.0);
335 TH2D* muRxU_tmp = (TH2D*) muRxU.Clone(
"muRxU.tmp");
337 TH2D muRxT(
"muRxT", NULL, 15, 0.0, 300.0,
getSize(
T) - 1,
T);
338 TH2D* muRxT_tmp = (TH2D*) muRxT.Clone(
"muRxT.tmp");
349 while (inputFile.hasNext()) {
351 STATUS(
"event: " << setw(10) << inputFile.getCounter() <<
'\r');
DEBUG(endl);
353 const Evt*
event = inputFile.next();
361 hits.Fill(
log10((Double_t) event->mc_hits.size()));
364 trks.Fill(track->type);
385 const int type = hit->type;
386 const double t1 =
getTime(*hit);
387 const double npe =
getNPE (*hit);
389 npe_pmt[hit->pmt_id].put(npe);
391 npe_type[0] .put(npe);
392 npe_type[type].put(npe);
394 if (hit_types.empty() || hit_types.count(type) != 0) {
397 event->mc_trks.end(),
400 if (track == event->mc_trks.end()) {
401 ERROR(
"Hit " << *hit <<
" has no origin." << endl);
405 if (count_if(event->mc_trks.begin(),
406 event->mc_trks.end(),
408 ERROR(
"Hit " << *hit <<
" has ambiguous origin." << endl);
412 job.Fill((
double) track->type, npe);
414 if (router.
hasPMT(hit->pmt_id)) {
422 const double E = muon.
getE();
423 const double x =
log10(E);
425 const double t0 = muon.
getT (pmt);
426 const double ct1 = muon.
getDot (pmt);
428 muExR.Fill(x, R, getMuonWeight(E, npe));
429 muRxU.Fill(R, ct1, getMuonWeight(E, R, npe));
430 muRxT.Fill(R, t1 - t0, getMuonWeight(E, R, npe));
434 const double E = neutrino.
getE();
435 const double x =
log10(E);
437 const double t0 = vertex.
getT(pmt);
439 const double ct1 = vertex.
getDot(pmt);
441 nuExD.Fill(x, D, getNeutrinoWeight(E, npe));
442 nuExc.Fill(x, ct0, getNeutrinoWeight(E, npe));
443 nuDxc.Fill(D, ct0, getNeutrinoWeight(E, D, npe));
444 nuDxU.Fill(D, ct1, getNeutrinoWeight(E, D, npe));
445 nuDxT.Fill(D, t1 - t0, getNeutrinoWeight(E, D, npe));
456 for (map_type::const_iterator i = npe_pmt.begin(); i != npe_pmt.end(); ++i) {
457 npe_per_pmt.Fill(i->second.getTotal());
460 for (map_type::const_iterator i = npe_type.begin(); i != npe_type.end(); ++i) {
461 npe_per_type[i->first]->Fill(i->second.getTotal());
475 const double E = muon.
getE();
476 const double x =
log10(E);
478 for (JDetector::const_iterator module =
detector.begin(); module !=
detector.end(); ++module) {
482 muExR_tmp->Fill(x, R, module->size());
484 if (R < muRxU.GetXaxis()->GetXmax()) {
485 for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
486 muRxU_tmp->Fill(R, muon.
getDot(*pmt));
490 if (R < muRxT.GetXaxis()->GetXmax()) {
491 for (Int_t iy = 1; iy <= muRxT_tmp->GetYaxis()->GetNbins(); ++iy) {
492 muRxT_tmp->Fill(R, muRxT_tmp->GetYaxis()->GetBinCenter(iy), muRxT_tmp->GetYaxis()->GetBinWidth(iy));
506 const double E = neutrino.
getE();
507 const double x =
log10(E);
509 if (inst_vol.is_inside(vertex))
510 hits_per_E_in .Fill(x, NPE);
512 hits_per_E_out.Fill(x, NPE);
514 for (JDetector::const_iterator module =
detector.begin(); module !=
detector.end(); ++module) {
517 const double ct0 = neutrino.
getDot(module->getPosition());
519 nuExD_tmp->Fill(x, D, module->size());
520 nuExc_tmp->Fill(x, ct0, module->size());
521 nuDxc_tmp->Fill(D, ct0, module->size());
523 if (D < nuDxU.GetXaxis()->GetXmax()) {
524 for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
525 nuDxU_tmp->Fill(D, neutrino.
getDot(*pmt));
529 if (D < nuDxT.GetXaxis()->GetXmax()) {
530 for (Int_t iy = 1; iy <= nuDxT_tmp->GetYaxis()->GetNbins(); ++iy) {
531 nuDxT_tmp->Fill(D, nuDxT_tmp->GetYaxis()->GetBinCenter(iy), nuDxT_tmp->GetYaxis()->GetBinWidth(iy));
545 TH1D *job_sorted = makeSortedH1(&job,
"hits_by_pdg");
546 TH1D *trks_sorted = makeSortedH1(&trks,
"trks_sorted");
548 if (inputFile.getCounter() != 0) {
550 const Double_t W = 1.0 / ((Double_t) inputFile.getCounter());
552 for (TH1D* p : { &hits, &npe_per_pmt, &job, &trks, job_sorted, trks_sorted }) {
556 for (JManager_t::iterator i = npe_per_type.begin(); i != npe_per_type.end(); ++i) {
560 nuExD.Divide(nuExD_tmp);
561 nuExc.Divide(nuExc_tmp);
562 nuDxc.Divide(nuDxc_tmp);
563 nuDxU.Divide(nuDxU_tmp);
564 nuDxT.Divide(nuDxT_tmp);
566 muExR.Divide(muExR_tmp);
567 muRxU.Divide(muRxU_tmp);
568 muRxT.Divide(muRxT_tmp);
578 out << hits << job << *job_sorted << trks << *trks_sorted << hits_per_E_in << hits_per_E_out << npe_per_type << npe_per_pmt;
580 out << nuExD << nuExc << nuDxc << nuDxU << nuDxT;
581 out << muExR << muRxU << muRxT;
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.
JTrack3E getTrack(const Trk &track)
Get track.
double geanc()
Equivalent muon track length per unit shower energy.
bool has_neutrino(const Evt &evt)
Test whether given event has an incoming neutrino.
bool is_muon(const Trk &track)
Test whether given track is a (anti-)muon.
double getDot(const JAxis3D &axis) const
Get cosine angle of impact of Cherenkov light on PMT.
then echo Enter input within $TIMEOUT_S seconds echo n User name
#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)...
static const JGeaneWater gWater
Function object for energy loss of muon in sea water.
Dynamic ROOT object management.
double getTime(const Hit &hit)
Get true time of hit.
double getDistance(const JVector3D &pos) const
Get distance to point.
size_t getSize(T(&array)[N])
Get size of c-array.
Data structure for detector geometry and calibration.
then fatal Wrong number of arguments fi set_variable DETECTOR $argv[1] set_variable STRING $argv[2] set_array QUANTILES set_variable FORMULA *[0] exp(-0.5 *(x-[1])*(x-[1])/([2]*[2]))" set_variable MODULE `getModule -a $DETECTOR -L "$STRING 0"` source JAcoustics.sh -- typeset -A TRIPODS get_tripods $WORKDIR/tripod.txt TRIPODS XMEAN
then echo The file $DIR KM3NeT_00000001_00000000 root already please rename or remove it first
Head getHeader(const JMultipleFileScanner_t &file_list)
Get Monte Carlo header.
Auxiliary class for defining the range of iterations of objects.
Auxiliary class to manage set of compatible ROOT objects (e.g. histograms) using unique keys...
I/O formatting auxiliaries.
double getE() const
Get energy.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
set_variable E_E log10(E_{fit}/E_{#mu})"
do set_variable OUTPUT_DIRECTORY $WORKDIR T
Data structure for PMT geometry, calibration and status.
Direct access to PMT in detector data structure.
const JPosition3D & getPosition() const
Get position.
then usage $script[distance] fi case set_variable R
double getT(const JVector3D &pos) const
Get arrival time of Cherenkov light at given position.
General purpose messaging.
virtual double getB() const override
Get energy loss constant.
Scanning of objects from multiple files according a format that follows from the extension of each fi...
JVertex3D getVertex() const
Get vertex of this track.
double getDistance(const JVector3D &pos) const
Get distance.
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
General purpose class for object reading from a list of file names.
Utility class to parse command line options.
then usage $script< input_file >< detector_file > fi set_variable OUTPUT_DIR set_variable SELECTOR JDAQTimesliceL1 set_variable DEBUG case set_variable DEBUG
double getNPE(const Hit &hit)
Get true charge of hit.
Data structure for position in three dimensions.
const JLimit & getLimit() const
Get limit.
const Trk & get_neutrino(const Evt &evt)
Get incoming neutrino.
do set_variable DETECTOR_TXT $WORKDIR detector
double getT(const JVector3D &pos) const
Get arrival time of Cherenkov light at given position.
then fatal Wrong number of arguments fi set_variable DETECTOR $argv[1] set_variable INPUT_FILE $argv[2] eval JPrintDetector a $DETECTOR O IDENTIFIER eval JPrintDetector a $DETECTOR O SUMMARY source JAcoustics sh $DETECTOR_ID CHECK_EXIT_CODE typeset A TRIPODS get_tripods $WORKDIR tripod txt TRIPODS for EMITTER in
do echo Generating $dir eval D
double getDot(const JAxis3D &axis) const
Get cosine angle of impact of Cherenkov light on PMT.
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
double getSinThetaC()
Get average sine of Cherenkov angle of water corresponding to group velocity.