38 #if defined(__unix__) || defined(__unix) || defined(unix) || (defined(__APPLE__) && defined(__MACH__))
40 #include <sys/resource.h>
41 #if defined(__APPLE__) && defined(__MACH__)
42 #include <mach/mach.h>
43 #elif (defined(_AIX) || defined(__TOS__AIX__)) || (defined(__sun__) || defined(__sun) || defined(sun) && (defined(__SVR4) || defined(__svr4__)))
46 #elif defined(__linux__) || defined(__linux) || defined(linux) || defined(__gnu_linux__)
53 #define G_SETUP_OBSERVATION "ctool::setup_observations(GObservations&)"
54 #define G_SETUP_MODELS "ctool::setup_models(GObservations&, std::string&)"
55 #define G_GET_MEAN_POINTING "ctool::get_mean_pointing(GObservations&)"
56 #define G_CREATE_EBOUNDS "ctool::create_ebounds()"
57 #define G_CREATE_ENERGIES "ctool::create_energies()"
58 #define G_RESTORE_EDISP "ctool::restore_edisp(GObservations&, "\
60 #define G_SET_OBS_RESPONSE "ctool::set_obs_response(GCTAObservation*)"
61 #define G_PROVIDE_HELP "ctool::provide_help()"
91 const std::string& version) :
GApplication(name, version)
117 const std::string& version,
118 const GApplicationPars& pars) :
GApplication(name, version, pars)
150 const std::string& version,
227 this->GApplication::operator=(app);
331 std::setlocale(LC_ALL,
"en_US.UTF-8");
375 char* ptr = std::getenv(
"CTOOLS");
377 std::string path = std::string(ptr) +
"/syspfiles";
378 std::string filename = path +
"/" + par_filename();
379 if (access(filename.c_str(), R_OK) == 0) {
381 m_pars.syspfiles(path);
382 m_pars.load(par_filename(), m_args);
432 const bool& response,
437 if (obs.size() == 0) {
441 if ((*
this)[
"inobs"].is_undefined()) {
442 std::string msg =
"The \"inobs\" parameter \""+
443 (*this)[
"inobs"].current_value()+
444 "\" is not a valid filename. Please specify "
445 "a valid event list, counts cube or "
446 "observation definition XML file.";
451 GFilename filename = (*this)[
"inobs"].filename();
454 if (!filename.exists()) {
455 std::string msg =
"The \"inobs\" parameter \""+filename.url()+
456 "\" specifies a file that does not "
457 "exist. Please specify a valid event "
458 "list, counts cube or observation "
459 "definition XML file.";
465 else if (filename.is_fits()) {
468 GCTAObservation cta(filename);
472 if (!list && cta.eventtype() ==
"EventList") {
473 std::string msg =
"The \"inobs\" parameter \""+filename.url()+
474 "\" specifies an event list, but this "
475 "tool does not accept event lists.";
481 if (!cube && cta.eventtype() ==
"CountsCube") {
482 std::string msg =
"The \"inobs\" parameter \""+filename.url()+
483 "\" specifies a counts cube, but this "
484 "tool does not accept counts cubes.";
546 const std::string& name)
551 if (obs.models().size() == 0) {
555 if ((*
this)[
"inmodel"].is_undefined()) {
556 std::string msg =
"The \"inmodel\" parameter \""+
557 (*this)[
"inmodel"].current_value()+
558 "\" is not a valid filename. Please specify "
559 "a valid model definition XML file.";
564 GFilename filename = (*this)[
"inmodel"].filename();
567 obs.models(GModels(filename.url()));
577 if (!obs.models().contains(name)) {
578 std::string msg =
"Model \""+name+
"\" not found in model "
579 "container. Please add a model with that name "
580 "or check for possible typos.";
618 std::string ebinalg = (*this)[
"ebinalg"].string();
622 if (ebinalg ==
"FILE") {
625 GFilename ebinfile = (*this)[
"ebinfile"].filename();
633 #pragma omp critical(ctool_create_ebounds)
639 if (!ebinfile.has_extname()) {
646 GFits file(ebinfile.url());
647 if (file.contains(gammalib::extname_ebounds)) {
649 ebounds.load(ebinfile);
651 else if (file.contains(
"ENERGYBINS")) {
653 ebinfile = ebinfile.url() +
"[ENERGYBINS]";
654 ebounds.load(ebinfile);
656 else if (file.contains(gammalib::extname_energies)) {
658 ebounds.set(GEnergies(ebinfile));
662 msg =
"No extension with name \""+
663 gammalib::extname_ebounds+
"\", "
664 "\"ENERGYBINS\" or \""+
665 gammalib::extname_energies+
"\" found in FITS "
666 "file \""+ebinfile+
"\". Please specify a "
667 "valid energy binning file.";
677 GFits file(ebinfile.url());
682 if (!file.contains(ebinfile.extname())) {
683 msg =
"No extension \""+ebinfile.extname()+
"\" "
684 "found in energy binning file \""+
685 ebinfile.url()+
"\". Please provide a valid "
693 GFitsTable& table = *file.table(ebinfile.extname());
697 if (table.ncols() == 1) {
699 ebounds.set(GEnergies(ebinfile));
703 ebounds.load(ebinfile);
714 catch (
const std::exception& e) {
733 double emin = (*this)[
"emin"].real();
734 double emax = (*this)[
"emax"].real();
735 int enumbins = (*this)[
"enumbins"].integer();
736 double ebingamma = 1.0;
739 if (ebinalg ==
"POW") {
740 ebingamma = (*this)[
"ebingamma"].real();
744 ebounds = GEbounds(enumbins, GEnergy(emin,
"TeV"), GEnergy(emax,
"TeV"),
782 std::string ebinalg = (*this)[
"ebinalg"].string();
786 if (ebinalg ==
"FILE") {
789 GFilename ebinfile = (*this)[
"ebinfile"].filename();
797 #pragma omp critical(ctool_create_energies)
803 if (!ebinfile.has_extname()) {
810 GFits file(ebinfile.url());
811 if (file.contains(gammalib::extname_energies)) {
813 energies.load(ebinfile);
815 else if (file.contains(gammalib::extname_ebounds)) {
817 energies.set(GEbounds(ebinfile));
819 else if (file.contains(
"ENERGYBINS")) {
821 ebinfile = ebinfile.url() +
"[ENERGYBINS]";
822 energies.set(GEbounds(ebinfile));
826 msg =
"No extension with name \""+
827 gammalib::extname_energies+
"\", "
828 "\"ENERGYBINS\" or \""+
829 gammalib::extname_ebounds+
"\" found in FITS "
830 "file \""+ebinfile+
"\". Please specify a "
831 "valid energy binning file.";
841 GFits file(ebinfile.url());
846 if (!file.contains(ebinfile.extname())) {
847 msg =
"No extension \""+ebinfile.extname()+
"\" "
848 "found in energy node file \""+
849 ebinfile.url()+
"\". Please provide a valid "
857 GFitsTable& table = *file.table(ebinfile.extname());
861 if (table.ncols() == 1) {
863 energies.load(ebinfile);
867 energies.set(GEbounds(ebinfile));
878 catch (
const std::exception& e) {
897 double emin = (*this)[
"emin"].real();
898 double emax = (*this)[
"emax"].real();
899 int enumbins = (*this)[
"enumbins"].integer();
900 double ebingamma = 1.0;
903 if (ebinalg ==
"POW") {
904 ebingamma = (*this)[
"ebingamma"].real();
908 energies = GEnergies(enumbins, GEnergy(emin,
"TeV"), GEnergy(emax,
"TeV"),
940 std::string coordsys = (*this)[
"coordsys"].string();
941 std::string proj = (*this)[
"proj"].string();
946 bool usepnt = (*this)[
"usepnt"].boolean();
948 xref = (*this)[
"xref"].real();
949 yref = (*this)[
"yref"].real();
953 double binsz = (*this)[
"binsz"].real();
954 int nxpix = (*this)[
"nxpix"].integer();
955 int nypix = (*this)[
"nypix"].integer();
964 if (gammalib::toupper(coordsys) ==
"GAL") {
970 yref = pnt.dec_deg();
976 GSkyMap map = GSkyMap(proj, coordsys, xref, yref, -binsz, binsz,
1016 gti.append(GTime(0.0), GTime(0.1234));
1019 map.nmaps(ebounds.size());
1022 GCTAEventCube cube = GCTAEventCube(map, ebounds, gti);
1059 double deadc = (*this)[
"deadc"].real();
1062 GCTAObservation obs;
1072 GTimeReference ref = (has_par(
"mjdref"))
1073 ? GTimeReference((*
this)[
"mjdref"].real(),
"s",
"TT",
"LOCAL")
1074 : GTimeReference(G_CTA_MJDREF,
"s",
"TT",
"LOCAL");
1077 std::string instrument = (has_par(
"instrument"))
1078 ? (*
this)[
"instrument"].string()
1090 list.ebounds(ebounds);
1096 obs.instrument(instrument);
1098 obs.ontime(gti.ontime());
1099 obs.livetime(gti.ontime()*deadc);
1122 if ((*
this)[
"inobs"].is_undefined()) {
1123 std::string msg =
"A valid file needs to be specified for the "
1124 "\"inobs\" parameter, yet \""+
1125 (*this)[
"inobs"].current_value()+
1127 " Specify a valid observation definition or "
1128 "FITS file to proceed.";
1129 throw GException::invalid_value(method, msg);
1133 GFilename filename = (*this)[
"inobs"].filename();
1152 if ((*
this)[
"inobs"].is_valid()) {
1155 GFilename filename = (*this)[
"inobs"].filename();
1158 if (filename.is_fits()) {
1161 bool is_cube =
false;
1168 GCTAEventCube cube(filename);
1182 std::string msg =
"A counts cube has been specified for the "
1183 "\"inobs\" parameter, yet no counts cube "
1184 "can be specified as input observation."
1185 " Instead, specify an event list or an "
1186 "observation definition file.";
1187 throw GException::invalid_value(method, msg);
1211 this->GApplication::log_parameters(chatter);
1215 if (has_par(
"edisp") && !(*
this)[
"edisp"].
boolean()) {
1218 GChatter chattiness =
static_cast<GChatter
>((&m_pars[
"chatter"])->integer());
1222 if (chattiness >= chatter) {
1226 log <<
"WARNING: Energy dispersion will *NOT* be considered for the "
1227 "computation. To consider" << std::endl;
1228 log <<
" energy dispersion, please set the \"edisp\" "
1229 "parameter to \"yes\". Be aware that" << std::endl;
1230 log <<
" using energy dispersion will considerably slow "
1231 "down the computations." << std::endl;
1252 const GObservations& obs,
1253 const std::string& what)
1256 GChatter chattiness =
static_cast<GChatter
>((*this)[
"chatter"].integer());
1260 if (chattiness >= chatter) {
1264 log.header1(gammalib::number(what, obs.size()));
1267 log << obs.print(chattiness) << std::endl;
1286 const GModels& models,
1287 const std::string& what)
1290 GChatter chattiness =
static_cast<GChatter
>((*this)[
"chatter"].integer());
1294 if (chattiness >= chatter) {
1298 log.header1(gammalib::number(what, models.size()));
1301 log << models.print(chattiness) << std::endl;
1332 bool usepnt = (has_par(
"usepnt") && (*this)[
"usepnt"].boolean());
1337 roi.centre(GCTAInstDir(pnt.dir()));
1342 else if ((*
this)[
"ra"].is_valid() && (*this)[
"dec"].is_valid()) {
1348 GCTAInstDir instdir(skydir);
1351 roi.centre(instdir);
1360 if (usepnt && (*
this)[
"rad"].is_valid()) {
1361 roi.radius((*
this)[
"rad"].real());
1384 if ((*
this)[
"emin"].is_valid() && (*
this)[
"emax"].is_valid()) {
1387 double e_min = (*this)[
"emin"].real();
1388 double e_max = (*this)[
"emax"].real();
1395 ebounds.append(emin, emax);
1421 if ((*
this)[
"tmin"].is_valid() && (*
this)[
"tmax"].is_valid()) {
1425 GTime tmin = (*this)[
"tmin"].time(ref);
1426 GTime tmax = (*this)[
"tmax"].time(ref);
1429 gti.append(tmin, tmax);
1451 GCTAPointing pnt(skydir);
1468 double ra = (*this)[
"ra"].real();
1469 double dec = (*this)[
"dec"].real();
1473 skydir.radec_deg(ra, dec);
1494 std::string prefix = (*this)[
"prefix"].string();
1497 GFilename fname(filename);
1500 std::vector<std::string> elements = gammalib::split(fname.url(),
"/");
1503 std::string outname = (elements.size() > 0) ?
1504 prefix + elements[elements.size()-1] :
"";
1507 std::string strip =
".gz";
1508 std::string::size_type i_start = outname.find(strip);
1509 if (i_start != std::string::npos) {
1510 outname.erase(i_start, strip.length());
1534 (*this)[
"emin"].real();
1535 (*this)[
"emax"].real();
1538 bool stacked = ((*this)[
"enumbins"].integer() > 0) ?
true :
false;
1543 (*this)[
"coordsys"].string();
1544 (*this)[
"proj"].string();
1545 (*this)[
"xref"].real();
1546 (*this)[
"yref"].real();
1547 (*this)[
"nxpix"].integer();
1548 (*this)[
"nypix"].integer();
1549 (*this)[
"binsz"].real();
1574 bool onoff = ((*this)[
"method"].string() ==
"ONOFF");
1580 (*this)[
"inexclusion"].query();
1581 (*this)[
"emin"].real();
1582 (*this)[
"emax"].real();
1583 (*this)[
"enumbins"].integer();
1584 (*this)[
"coordsys"].string();
1585 (*this)[
"xref"].real();
1586 (*this)[
"yref"].real();
1587 if ((*
this)[
"srcshape"].
string() ==
"CIRCLE"){
1588 (*this)[
"rad"].real();
1590 if ((*
this)[
"bkgmethod"].
string() ==
"REFLECTED"){
1591 (*this)[
"bkgregmin"].integer();
1593 (*this)[
"maxoffset"].real();
1594 (*this)[
"etruemin"].real();
1595 (*this)[
"etruemax"].real();
1596 (*this)[
"etruebins"].integer();
1617 for (
int i = 0; i < obs.size(); ++i) {
1620 GCTAObservation* cta =
dynamic_cast<GCTAObservation*
>(obs[i]);
1626 if (!cta->has_response()) {
1653 std::vector<bool> old_edisp(obs.size(),
false);
1656 for (
int i = 0; i < obs.size(); ++i) {
1659 GCTAObservation* cta =
dynamic_cast<GCTAObservation*
>(obs[i]);
1665 old_edisp[i] = cta->response()->apply_edisp();
1668 cta->response()->apply_edisp(edisp);
1693 if (edisp.size() != obs.size()) {
1694 std::string msg =
"Size of energy dispersion flag vector ("+
1695 gammalib::str(edisp.size())+
" elements) is "
1696 "incompatible with number of observations in the "
1697 "observation container ("+gammalib::str(obs.size())+
1703 for (
int i = 0; i < obs.size(); ++i) {
1706 GCTAObservation* cta =
dynamic_cast<GCTAObservation*
>(obs[i]);
1710 cta->response()->apply_edisp(edisp[i]);
1755 bool has_response =
false;
1760 if (dynamic_cast<const GCTAEventCube*>(obs->events()) != NULL) {
1764 if (has_par(
"expcube") && has_par(
"psfcube") &&
1765 has_par(
"bkgcube") && has_par(
"edispcube")) {
1770 bool query_edisp =
true;
1771 if (has_par(
"edisp")) {
1772 query_edisp = (*this)[
"edisp"].boolean();
1776 if ((*
this)[
"expcube"].is_valid() &&
1777 (*
this)[
"psfcube"].is_valid() &&
1778 (*
this)[
"bkgcube"].is_valid()) {
1782 GCTACubeExposure exposure((*
this)[
"expcube"].filename());
1783 GCTACubePsf psf((*
this)[
"psfcube"].filename());
1784 if (query_edisp && (*
this)[
"edispcube"].is_valid()) {
1785 (*this)[
"edispcube"].filename();
1787 GCTACubeBackground background((*
this)[
"bkgcube"].filename());
1793 if ((*
this)[
"edispcube"].is_valid()) {
1794 GCTACubeEdisp edisp((*
this)[
"edispcube"].filename());
1795 obs->response(exposure, psf, edisp, background);
1803 if (has_par(
"edisp") && (*
this)[
"edisp"].
boolean()) {
1804 std::string msg =
"Energy dispersion requested but "
1805 "no energy dispersion cube was "
1812 obs->response(exposure, psf, background);
1820 obs->response(exposure, psf, background);
1824 has_response =
true;
1833 if (!has_response) {
1836 std::string database = (*this)[
"caldb"].string();
1837 std::string irf = (*this)[
"irf"].string();
1840 std::string instrument = (has_par(
"instrument"))
1841 ? (*
this)[
"instrument"].string()
1848 std::string parameter =
"parameter name=\"Calibration\""
1849 " database=\""+database+
"\""
1850 " response=\""+irf+
"\"";
1856 xml.attribute(
"instrument", instrument);
1859 xml.append(parameter);
1862 GCTAResponseIrf response(xml);
1865 obs->response(response);
1892 if ((*
this)[
"inobs"].is_undefined()) {
1906 obs.append(cta_obs);
1914 std::string filename = (*this)[
"inobs"].filename();
1918 if (GFilename(filename).is_fits()) {
1921 GCTAObservation cta_obs(filename);
1932 obs.append(cta_obs);
1987 for (
int i = 0; i < obs.size(); ++i) {
1990 const GCTAObservation* cta_obs =
dynamic_cast<const GCTAObservation*
>(obs[i]);
1991 if (cta_obs != NULL) {
1994 const GCTAPointing& pnt = cta_obs->pointing();
1997 xref += pnt.dir().ra_deg();
1998 yref += pnt.dir().dec_deg();
2007 std::string msg =
"No valid CTA observation has been found in "
2008 "observation list, hence no pointing information "
2009 "could be extracted. Use the \"usepnt=no\" "
2010 "option and specify the pointing explicitly.";
2015 double ra = xref / double(n_pnt);
2016 double dec = yref / double(n_pnt);
2019 GSkyDir dir = GSkyDir();
2020 dir.radec_deg(ra,dec);
2042 #if defined(__APPLE__) && defined(__MACH__)
2043 #ifdef MACH_TASK_BASIC_INFO
2044 struct mach_task_basic_info info;
2045 mach_msg_type_number_t infoCount = MACH_TASK_BASIC_INFO_COUNT;
2046 if (task_info(mach_task_self( ), MACH_TASK_BASIC_INFO,
2047 (task_info_t)&info, &infoCount) == KERN_SUCCESS) {
2048 rss = (size_t)info.resident_size;
2051 struct task_basic_info info;
2052 mach_msg_type_number_t info_count = TASK_BASIC_INFO_COUNT;
2053 if (task_info(mach_task_self(), TASK_BASIC_INFO,
2054 (task_info_t)&info, &info_count) == KERN_SUCCESS) {
2055 rss = (size_t)info.resident_size;
2059 #elif defined(__linux__) || defined(__linux) || defined(linux) || defined(__gnu_linux__)
2061 if ((fp = fopen(
"/proc/self/statm",
"r" )) != NULL) {
2062 if (fscanf( fp,
"%*s%ld", &rss ) == 1) {
2063 rss *= (size_t)sysconf(_SC_PAGESIZE);
2093 std::string header = obs->instrument() +
" observation";
2096 if (!obs->name().empty()) {
2097 header +=
" \"" + obs->name() +
"\"";
2101 if (!obs->id().empty()) {
2102 header +=
" (id=" + obs->id() +
")";
2120 const GCTAObservation& obs)
2123 GEnergies engs = energies;
2126 GEbounds ebounds = obs.ebounds();
2129 for (
int iebin = 0; iebin <= ebounds.size(); ++iebin) {
2133 if (iebin < ebounds.size()) {
2134 energy = ebounds.emin(iebin);
2137 energy = ebounds.emax();
2143 for (
int k = 0; k < engs.size(); ++k) {
2147 if (std::abs(engs[k].MeV()-energy.MeV()) < 1.0) {
2153 if (engs[k] > energy) {
2165 engs.insert(index, energy);
2168 engs.append(energy);
2172 log_value(NORMAL,
"Insert energy", energy.print());
2195 const GEbounds& list_ebounds)
const
2198 const GEnergy energy_margin(0.01,
"GeV");
2201 std::vector<bool> usage(cube_ebounds.size(),
true);
2204 for (
int iebin = 0; iebin < cube_ebounds.size(); ++iebin) {
2214 usage[iebin] = list_ebounds.contains(cube_ebounds.emin(iebin) +
2216 cube_ebounds.emax(iebin) -
2267 const std::string& infile,
2268 const std::string& evtname,
2269 const std::string& gtiname,
2270 const std::string& outfile)
const
2273 if (obs->eventtype() ==
"EventList") {
2276 GFilename outname(outfile);
2277 std::string outevt = evtname;
2278 std::string outgti = gtiname;
2279 if (outname.has_extname()) {
2280 std::vector<std::string> extnames =
2281 gammalib::split(outname.extname(),
";");
2282 if (extnames.size() > 0) {
2283 std::string extname = gammalib::strip_whitespace(extnames[0]);
2284 if (!extname.empty()) {
2288 if (extnames.size() > 1) {
2289 std::string extname = gammalib::strip_whitespace(extnames[1]);
2290 if (!extname.empty()) {
2297 #pragma omp critical(ctool_save_event_list)
2304 obs->write(outfits, outevt, outgti);
2311 GFits infits(infile);
2312 for (
int extno = 1; extno < infits.size(); ++extno) {
2313 GFitsHDU* hdu = infits.at(extno);
2314 if (hdu->extname() != evtname &&
2315 hdu->extname() != gtiname &&
2316 hdu->extname() != outevt &&
2317 hdu->extname() != outgti) {
2318 outfits.append(*hdu);
2329 outfits.saveto(outname.url(), clobber());
2353 const std::string& evtname)
const
2356 std::string gtiname = gammalib::extname_gti;
2359 if (!filename.empty()) {
2362 GCTAEventList events(filename+
"["+evtname+
"]");
2365 gtiname = events.gtiname();
2389 std::string warning;
2392 int n = energies.size();
2397 double logEmin = std::log10(energies[0].TeV());
2398 double logEmax = std::log10(energies[n-1].TeV());
2399 nrequired = int((logEmax - logEmin) * 25.0);
2403 if (n < nrequired) {
2404 warning.append(
"\nWARNING: Only "+gammalib::str(n-1)+
" energy bins "
2405 "have been requested. This may be too few energy "
2407 "\n At least 25 bins per decade in energy are "
2408 "recommended, which for the given"
2409 "\n energy range would be "+
2410 gammalib::str(nrequired)+
" bins. Consider increasing "
2411 "the number of energy bins.");
2430 std::string warning;
2433 std::string fname = std::string(filename);
2434 std::string suffix = gammalib::tolower(fname.substr(fname.length()-4,4));
2437 if (suffix !=
".xml") {
2438 warning.append(
"WARNING: Name of observation definition output file "
2439 "\""+filename+
"\" does not terminate with \".xml\"."
2440 "\n This is not an error, but may be misleading. "
2441 "It is recommended to use the suffix"
2442 "\n \".xml\" for observation definition files.");
2473 std::string helpfile = name()+
".txt";
2476 char* ptr = std::getenv(
"CTOOLS");
2478 std::string msg =
"CTOOLS environment variable not set, cannot "
2479 "display help file. Please set the CTOOLS "
2480 "environment variable.";
2486 std::string fname = std::string(ptr) +
"/share/help/" + helpfile;
2487 FILE* fptr = fopen(fname.c_str(),
"r");
2490 while (fgets(line, n, fptr) != NULL) {
2491 std::cout << std::string(line);
2496 std::cout <<
"No help available for "+name()+
"." << std::endl;