jean perrot / ctools

Browse Code »

Commit 3d85a1fb504a47f96da12814f0eeb2a405197fda

Authored by ldivenere 2017-04-05 11:57:11 +0200

1 parent 1c2f5a5a

added ctprob files

Inline Side-by-side

Showing 5 changed files with 1144 additions and 0 deletions

src/ctprob/Makefile.am 0 → 100644

View file @3d85a1f

	1	+#############################################################################
	2	+# Makefile.am for ctmodel tool #
	3	+# ------------------------------------------------------------------------- #
	4	+# Copyright (C) 2016 Juergen Knoedlseder #
	5	+# ------------------------------------------------------------------------- #
	6	+# #
	7	+# This program is free software: you can redistribute it and/or modify #
	8	+# it under the terms of the GNU General Public License as published by #
	9	+# the Free Software Foundation, either version 3 of the License, or #
	10	+# (at your option) any later version. #
	11	+# #
	12	+# This program is distributed in the hope that it will be useful, #
	13	+# but WITHOUT ANY WARRANTY; without even the implied warranty of #
	14	+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the #
	15	+# GNU General Public License for more details. #
	16	+# #
	17	+# You should have received a copy of the GNU General Public License #
	18	+# along with this program. If not, see <http://www.gnu.org/licenses/>. #
	19	+# #
	20	+#############################################################################
	21	+
	22	+# Set include directories
	23	+INCLUDES = -I$(top_srcdir)/src/support
	24	+
	25	+# Program name
	26	+bin_PROGRAMS = ctmodel
	27	+
	28	+# Sources
	29	+sources = ctmodel.cpp ctmodel.hpp
	30	+
	31	+# Libtool library
	32	+noinst_LTLIBRARIES = libctmodel.la
	33	+libctmodel_la_SOURCES = $(sources)
	34	+libctmodel_la_LDFLAGS = @OPENMP_LIB_LDFLAGS@
	35	+
	36	+# Program definition
	37	+ctmodel_SOURCES = main.cpp
	38	+ctmodel_LDFLAGS = $(LIB_LINK_OPTION) @LDFLAGS@
	39	+ctmodel_CPPFLAGS = @CPPFLAGS@
	40	+ctmodel_LDADD = libctmodel.la $(top_builddir)/src/support/libsupport.la
	41	+
	42	+# Parameter files
	43	+pfilesdir = $(prefix)/syspfiles
	44	+pfiles_DATA = $(srcdir)/ctmodel.par
	45	+
	46	+# Extra files to be distributed in the tar file
	47	+EXTRA_DIST = $(pfiles_DATA)
...	...

src/ctprob/ctprob.cpp 0 → 100644

View file @3d85a1f

	1	+/***************************************************************************
	2	+ * ctmodel - Model cube generation tool *
	3	+ * ----------------------------------------------------------------------- *
	4	+ * copyright (C) 2012-2016 by Juergen Knoedlseder *
	5	+ * ----------------------------------------------------------------------- *
	6	+ * *
	7	+ * This program is free software: you can redistribute it and/or modify *
	8	+ * it under the terms of the GNU General Public License as published by *
	9	+ * the Free Software Foundation, either version 3 of the License, or *
	10	+ * (at your option) any later version. *
	11	+ * *
	12	+ * This program is distributed in the hope that it will be useful, *
	13	+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
	14	+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
	15	+ * GNU General Public License for more details. *
	16	+ * *
	17	+ * You should have received a copy of the GNU General Public License *
	18	+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
	19	+ * *
	20	+ ***************************************************************************/
	21	+/**
	22	+ * @file ctmodel.cpp
	23	+ * @brief Model cube generation tool implementation
	24	+ * @author Juergen Knoedlseder
	25	+ */
	26	+
	27	+/* __ Includes ___________________________________________________________ */
	28	+#ifdef HAVE_CONFIG_H
	29	+#include <config.h>
	30	+#endif
	31	+#include <cstdio>
	32	+#include "ctmodel.hpp"
	33	+#include "GTools.hpp"
	34	+
	35	+/* __ Method name definitions ____________________________________________ */
	36	+#define G_GET_PARAMETERS "ctmodel::get_parameters()"
	37	+#define G_FILL_CUBE "ctmodel::fill_cube(GCTAObservation*)"
	38	+
	39	+/* __ Debug definitions __________________________________________________ */
	40	+
	41	+/* __ Coding definitions _________________________________________________ */
	42	+
	43	+/* __ Constants __________________________________________________________ */
	44	+const GEnergy g_energy_margin(1.0e-12, "TeV");
	45	+
	46	+
	47	+/*==========================================================================
	48	+ = =
	49	+ = Constructors/destructors =
	50	+ = =
	51	+ ==========================================================================*/
	52	+
	53	+/*********************************************************************//
	54	+ * @brief Void constructor
	55	+ ***************************************************************************/
	56	+ctmodel::ctmodel(void) : ctobservation(CTMODEL_NAME, CTMODEL_VERSION)
	57	+{
	58	+ // Initialise members
	59	+ init_members();
	60	+
	61	+ // Return
	62	+ return;
	63	+}
	64	+
	65	+
	66	+/*********************************************************************//
	67	+ * @brief Observations constructor
	68	+ *
	69	+ * param[in] obs Observation container.
	70	+ *
	71	+ * This method creates an instance of the class by copying an existing
	72	+ * observations container.
	73	+ ***************************************************************************/
	74	+ctmodel::ctmodel(const GObservations& obs) :
	75	+ ctobservation(CTMODEL_NAME, CTMODEL_VERSION, obs)
	76	+{
	77	+ // Initialise members
	78	+ init_members();
	79	+
	80	+ // Return
	81	+ return;
	82	+}
	83	+
	84	+
	85	+
	86	+/*********************************************************************//
	87	+ * @brief Command line constructor
	88	+ *
	89	+ * @param[in] argc Number of arguments in command line.
	90	+ * @param[in] argv Array of command line arguments.
	91	+ ***************************************************************************/
	92	+ctmodel::ctmodel(int argc, char *argv[]) :
	93	+ ctobservation(CTMODEL_NAME, CTMODEL_VERSION, argc, argv)
	94	+{
	95	+ // Initialise members
	96	+ init_members();
	97	+
	98	+ // Return
	99	+ return;
	100	+}
	101	+
	102	+
	103	+/*********************************************************************//
	104	+ * @brief Copy constructor
	105	+ *
	106	+ * @param[in] app Application.
	107	+ ***************************************************************************/
	108	+ctmodel::ctmodel(const ctmodel& app) : ctobservation(app)
	109	+{
	110	+ // Initialise members
	111	+ init_members();
	112	+
	113	+ // Copy members
	114	+ copy_members(app);
	115	+
	116	+ // Return
	117	+ return;
	118	+}
	119	+
	120	+
	121	+/*********************************************************************//
	122	+ * @brief Destructor
	123	+ ***************************************************************************/
	124	+ctmodel::~ctmodel(void)
	125	+{
	126	+ // Free members
	127	+ free_members();
	128	+
	129	+ // Return
	130	+ return;
	131	+}
	132	+
	133	+
	134	+/*==========================================================================
	135	+ = =
	136	+ = Operators =
	137	+ = =
	138	+ ==========================================================================*/
	139	+
	140	+/*********************************************************************//
	141	+ * @brief Assignment operator
	142	+ *
	143	+ * @param[in] app Application.
	144	+ * @return Application.
	145	+ ***************************************************************************/
	146	+ctmodel& ctmodel::operator=(const ctmodel& app)
	147	+{
	148	+ // Execute only if object is not identical
	149	+ if (this != &app) {
	150	+
	151	+ // Copy base class members
	152	+ this->ctobservation::operator=(app);
	153	+
	154	+ // Free members
	155	+ free_members();
	156	+
	157	+ // Initialise members
	158	+ init_members();
	159	+
	160	+ // Copy members
	161	+ copy_members(app);
	162	+
	163	+ } // endif: object was not identical
	164	+
	165	+ // Return this object
	166	+ return *this;
	167	+}
	168	+
	169	+
	170	+/*==========================================================================
	171	+ = =
	172	+ = Public methods =
	173	+ = =
	174	+ ==========================================================================*/
	175	+
	176	+/*********************************************************************//
	177	+ * @brief Clear ctmodel tool
	178	+ *
	179	+ * Clears ctmodel tool.
	180	+ ***************************************************************************/
	181	+void ctmodel::clear(void)
	182	+{
	183	+ // Free members
	184	+ free_members();
	185	+ this->ctobservation::free_members();
	186	+ this->ctool::free_members();
	187	+
	188	+ // Clear base class (needed to conserve tool name and version)
	189	+ this->GApplication::clear();
	190	+
	191	+ // Initialise members
	192	+ this->ctool::init_members();
	193	+ this->ctobservation::init_members();
	194	+ init_members();
	195	+
	196	+ // Write header into logger
	197	+ log_header();
	198	+
	199	+ // Return
	200	+ return;
	201	+}
	202	+
	203	+
	204	+/*********************************************************************//
	205	+ * @brief Generate the model map(s)
	206	+ *
	207	+ * This method reads the task parameters from the parfile, sets up the
	208	+ * observation container, loops over all CTA observations in the container
	209	+ * and generates a model map for each CTA observation.
	210	+ ***************************************************************************/
	211	+void ctmodel::run(void)
	212	+{
	213	+ // If we're in debug mode then all output is also dumped on the screen
	214	+ if (logDebug()) {
	215	+ log.cout(true);
	216	+ }
	217	+
	218	+ // Get task parameters
	219	+ get_parameters();
	220	+
	221	+ // Set energy dispersion flags of all CTA observations and save old
	222	+ // values in save_edisp vector
	223	+ std::vector<bool> save_edisp = set_edisp(m_obs, m_apply_edisp);
	224	+
	225	+ // Write input observation container into logger
	226	+ log_observations(NORMAL, m_obs, "Input observation");
	227	+
	228	+ // Write input model container into logger
	229	+ log_models(NORMAL, m_obs.models(), "Input model");
	230	+
	231	+ // Write header
	232	+ log_header1(TERSE, "Generate model cube");
	233	+
	234	+ // Loop over all observations in the container
	235	+ for (int i = 0; i < m_obs.size(); ++i) {
	236	+
	237	+ // Write header for the current observation
	238	+ log_header3(TERSE, get_obs_header(m_obs[i]));
	239	+
	240	+ // Get CTA observation
	241	+ GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
	242	+
	243	+ // Skip observation if it's not CTA
	244	+ if (obs == NULL) {
	245	+ std::string msg = " Skipping "+m_obs[i]->instrument()+
	246	+ " observation";
	247	+ log_string(NORMAL, msg);
	248	+ continue;
	249	+ }
	250	+
	251	+ // Fill cube and leave loop if we are binned mode (meaning we
	252	+ // only have one binned observation)
	253	+ if (m_binned) {
	254	+ fill_cube(obs);
	255	+ break;
	256	+ }
	257	+
	258	+ // Skip observation if we have a binned observation
	259	+ if (obs->eventtype() == "CountsCube") {
	260	+ std::string msg = " Skipping binned "+m_obs[i]->instrument()+
	261	+ " observation";
	262	+ log_string(NORMAL, msg);
	263	+ continue;
	264	+ }
	265	+
	266	+ // Fill the cube
	267	+ fill_cube(obs);
	268	+
	269	+ // Dispose events to free memory if event file exists on disk
	270	+ if (obs->eventfile().length() > 0 && obs->eventfile().exists()) {
	271	+ obs->dispose_events();
	272	+ }
	273	+
	274	+ } // endfor: looped over observations
	275	+
	276	+ // Write model cube into header
	277	+ log_header1(NORMAL, "Model cube");
	278	+ log_string(NORMAL, m_cube.print());
	279	+
	280	+ // Restore energy dispersion flags of all CTA observations
	281	+ restore_edisp(m_obs, save_edisp);
	282	+
	283	+ // Optionally publish model cube
	284	+ if (m_publish) {
	285	+ publish();
	286	+ }
	287	+
	288	+ // Return
	289	+ return;
	290	+}
	291	+
	292	+
	293	+/*********************************************************************//
	294	+ * @brief Save model cube
	295	+ *
	296	+ * Saves the model cube into a FITS file specified using the "outfile"
	297	+ * task parameter.
	298	+ ***************************************************************************/
	299	+void ctmodel::save(void)
	300	+{
	301	+ // Write header
	302	+ log_header1(TERSE, "Save model cube");
	303	+
	304	+ // Get model cube filename
	305	+ m_outcube = (*this)["outcube"].filename();
	306	+
	307	+ // Save only if filename is non-empty and cube has some size
	308	+ if (!m_outcube.is_empty() && m_cube.size() > 0) {
	309	+ m_cube.save(m_outcube, clobber());
	310	+ }
	311	+
	312	+ // Write into logger what has been done
	313	+ std::string fname = (m_outcube.is_empty()) ? "NONE" : m_outcube.url();
	314	+ if (m_cube.size() == 0) {
	315	+ fname.append(" (cube is empty, no file created)");
	316	+ }
	317	+ log_value(NORMAL, "Model cube file", fname);
	318	+
	319	+ // Return
	320	+ return;
	321	+}
	322	+
	323	+
	324	+/*********************************************************************//
	325	+ * @brief Publish model cube
	326	+ *
	327	+ * @param[in] name Model cube name.
	328	+ ***************************************************************************/
	329	+void ctmodel::publish(const std::string& name)
	330	+{
	331	+ // Write header into logger
	332	+ log_header1(TERSE, "Publish model cube");
	333	+
	334	+ // Set default name is user name is empty
	335	+ std::string user_name(name);
	336	+ if (user_name.empty()) {
	337	+ user_name = CTMODEL_NAME;
	338	+ }
	339	+
	340	+ // Log filename
	341	+ log_value(NORMAL, "Model cube name", user_name);
	342	+
	343	+ // Publish model cube
	344	+ m_cube.counts().publish(user_name);
	345	+
	346	+ // Return
	347	+ return;
	348	+}
	349	+
	350	+
	351	+/*********************************************************************//
	352	+ * @brief Set model cube
	353	+ *
	354	+ * @param[in] cube Model cube.
	355	+ *
	356	+ * Set model cube and set all cube bins to zero.
	357	+ ***************************************************************************/
	358	+void ctmodel::cube(const GCTAEventCube& cube)
	359	+{
	360	+ // Set cube
	361	+ m_cube = cube;
	362	+
	363	+ // Set all cube bins to zero
	364	+ for (int i = 0; i < m_cube.size(); ++i) {
	365	+ m_cube[i]->counts(0.0);
	366	+ }
	367	+
	368	+ // Signal that cube has been set
	369	+ m_has_cube = true;
	370	+
	371	+ // Return
	372	+ return;
	373	+}
	374	+
	375	+
	376	+/*==========================================================================
	377	+ = =
	378	+ = Private methods =
	379	+ = =
	380	+ ==========================================================================*/
	381	+
	382	+/*********************************************************************//
	383	+ * @brief Initialise class members
	384	+ ***************************************************************************/
	385	+void ctmodel::init_members(void)
	386	+{
	387	+ // Initialise members
	388	+ m_outcube.clear();
	389	+ m_apply_edisp = false;
	390	+ m_publish = false;
	391	+ m_chatter = static_cast<GChatter>(2);
	392	+
	393	+ // Initialise protected members
	394	+ m_cube.clear();
	395	+ m_gti.clear();
	396	+ m_has_cube = false;
	397	+ m_append_cube = false;
	398	+ m_binned = false;
	399	+
	400	+ // Return
	401	+ return;
	402	+}
	403	+
	404	+
	405	+/*********************************************************************//
	406	+ * @brief Copy class members
	407	+ *
	408	+ * @param[in] app Application.
	409	+ ***************************************************************************/
	410	+void ctmodel::copy_members(const ctmodel& app)
	411	+{
	412	+ // Copy attributes
	413	+ m_outcube = app.m_outcube;
	414	+ m_apply_edisp = app.m_apply_edisp;
	415	+ m_publish = app.m_publish;
	416	+ m_chatter = app.m_chatter;
	417	+
	418	+ // Copy protected members
	419	+ m_cube = app.m_cube;
	420	+ m_gti = app.m_gti;
	421	+ m_has_cube = app.m_has_cube;
	422	+ m_append_cube = app.m_append_cube;
	423	+ m_binned = app.m_binned;
	424	+
	425	+ // Return
	426	+ return;
	427	+}
	428	+
	429	+
	430	+/*********************************************************************//
	431	+ * @brief Delete class members
	432	+ ***************************************************************************/
	433	+void ctmodel::free_members(void)
	434	+{
	435	+ // Return
	436	+ return;
	437	+}
	438	+
	439	+
	440	+/*********************************************************************//
	441	+ * @brief Get application parameters
	442	+ *
	443	+ * Get all task parameters from parameter file or (if required) by querying
	444	+ * the user. The parameters are read in the correct order.
	445	+ ***************************************************************************/
	446	+void ctmodel::get_parameters(void)
	447	+{
	448	+ // Reset cube append flag
	449	+ m_append_cube = false;
	450	+
	451	+ // If there are no observations in container then load them via user
	452	+ // parameters.
	453	+ if (m_obs.size() == 0) {
	454	+ get_obs();
	455	+ }
	456	+
	457	+ // ... otherwise add response information and energy boundaries in case
	458	+ // that they are missing
	459	+ else {
	460	+ setup_observations(m_obs);
	461	+ }
	462	+
	463	+ // If we have now excactly one CTA observation (but no cube has yet been
	464	+ // appended to the observation) then check whether this observation
	465	+ // is a binned observation, and if yes, extract the counts cube for
	466	+ // model generation
	467	+ if ((m_obs.size() == 1) && (m_append_cube == false)) {
	468	+
	469	+ // Get CTA observation
	470	+ GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[0]);
	471	+
	472	+ // Continue only if observation is a CTA observation
	473	+ if (obs != NULL) {
	474	+
	475	+ // Check for binned observation
	476	+ if (obs->eventtype() == "CountsCube") {
	477	+
	478	+ // Set cube from binned observation
	479	+ GCTAEventCube* evtcube = dynamic_cast<GCTAEventCube*>
	480	+ (const_cast<GEvents*>(obs->events()));
	481	+
	482	+ cube(*evtcube);
	483	+
	484	+ // Signal that cube has been set
	485	+ m_has_cube = true;
	486	+
	487	+ // Signal that we are in binned mode
	488	+ m_binned = true;
	489	+
	490	+ } // endif: observation was binned
	491	+
	492	+ } // endif: observation was CTA
	493	+
	494	+ } // endif: had exactly one observation
	495	+
	496	+ // Read model definition file if required
	497	+ if (m_obs.models().size() == 0) {
	498	+
	499	+ // Get model filename
	500	+ std::string inmodel = (*this)["inmodel"].filename();
	501	+
	502	+ // Load models from file
	503	+ m_obs.models(inmodel);
	504	+
	505	+ } // endif: there were no models
	506	+
	507	+ // Get energy dispersion flag parameters
	508	+ m_apply_edisp = (*this)["edisp"].boolean();
	509	+
	510	+ // If we do not have yet a counts cube for model computation then check
	511	+ // whether we should read it from the "incube" parameter or whether we
	512	+ // should create it from scratch using the task parameters
	513	+ if (!m_has_cube) {
	514	+
	515	+ // Read cube definition file
	516	+ std::string incube = (*this)["incube"].filename();
	517	+
	518	+ // If the cube filename is valid the load the cube and set all cube
	519	+ // bins to zero
	520	+ if (is_valid_filename(incube)) {
	521	+
	522	+ // Load cube from given file
	523	+ m_cube.load(incube);
	524	+
	525	+ // Set all cube bins to zero
	526	+ for (int i = 0; i < m_cube.size(); ++i) {
	527	+ m_cube[i]->counts(0.0);
	528	+ }
	529	+
	530	+ } // endif: cube filename was valid
	531	+
	532	+ // ... otherwise create a cube from the user parameters
	533	+ else {
	534	+ m_cube = create_cube(m_obs);
	535	+ }
	536	+
	537	+ // Signal that cube has been set
	538	+ m_has_cube = true;
	539	+
	540	+ } // endif: we had no cube yet
	541	+
	542	+ // Get remaining parameters
	543	+ m_publish = (*this)["publish"].boolean();
	544	+ m_chatter = static_cast<GChatter>((*this)["chatter"].integer());
	545	+
	546	+ // Read optionally output cube filenames
	547	+ if (read_ahead()) {
	548	+ m_outcube = (*this)["outcube"].filename();
	549	+ }
	550	+
	551	+ // If cube should be appended to first observation then do that now.
	552	+ // This is a kluge that makes sure that the cube is passed as part
	553	+ // of the observation in case that a cube response is used. The kluge
	554	+ // is needed because the GCTACubeSourceDiffuse::set method needs to
	555	+ // get the full event cube from the observation. It is also at this
	556	+ // step that the GTI, which may just be a dummy GTI when create_cube()
	557	+ // has been used, will be set.
	558	+ if (m_append_cube) {
	559	+
	560	+ //TODO: Check that energy boundaries are compatible
	561	+
	562	+ // Attach GTI of observations to model cube
	563	+ m_cube.gti(m_obs[0]->events()->gti());
	564	+
	565	+ // Attach model cube to observations
	566	+ m_obs[0]->events(m_cube);
	567	+
	568	+ } // endif: cube was scheduled for appending
	569	+
	570	+ // Write parameters into logger
	571	+ log_parameters(TERSE);
	572	+
	573	+ // Return
	574	+ return;
	575	+}
	576	+
	577	+
	578	+/*********************************************************************//
	579	+ * @brief Get observation container
	580	+ *
	581	+ * Get an observation container according to the user parameters. The method
	582	+ * supports loading of a individual FITS file or an observation definition
	583	+ * file in XML format.
	584	+ *
	585	+ * If the input filename is empty, the method checks for the existence of the
	586	+ * "expcube", "psfcube" and "bkgcube" parameters. If file names have been
	587	+ * specified, the method loads the files and creates a dummy events cube that
	588	+ * is appended to the observation container.
	589	+ *
	590	+ * If no file names are specified for the "expcube", "psfcube" or "bkgcube"
	591	+ * parameters, the method reads the necessary parameters to build a CTA
	592	+ * observation from scratch.
	593	+ *
	594	+ * The method sets m_append_cube = true and m_binned = true in case that
	595	+ * a stacked observation is requested (as detected by the presence of the
	596	+ * "expcube", "psfcube", and "bkgcube" parameters). In that case, it appended
	597	+ * a dummy event cube to the observation.
	598	+ *
	599	+ * @todo Support stacked energy dispersion
	600	+ ***************************************************************************/
	601	+void ctmodel::get_obs(void)
	602	+{
	603	+ // Get the filename from the input parameters
	604	+ std::string filename = (*this)["inobs"].filename();
	605	+
	606	+ // If no observation definition file has been specified then read all
	607	+ // parameters that are necessary to create an observation from scratch
	608	+ if (!is_valid_filename(filename)) {
	609	+
	610	+ // Get response cube filenames
	611	+ std::string expcube = (*this)["expcube"].filename();
	612	+ std::string psfcube = (*this)["psfcube"].filename();
	613	+ std::string bkgcube = (*this)["bkgcube"].filename();
	614	+
	615	+ // If the filenames are valid then build an observation from cube
	616	+ // response information
	617	+ if (is_valid_filename(expcube) && is_valid_filename(psfcube) &&
	618	+ is_valid_filename(bkgcube)) {
	619	+
	620	+ // Get exposure, PSF and background cubes
	621	+ GCTACubeExposure exposure(expcube);
	622	+ GCTACubePsf psf(psfcube);
	623	+ GCTACubeBackground background(bkgcube);
	624	+
	625	+ // Create energy boundaries
	626	+ GEbounds ebounds = create_ebounds();
	627	+
	628	+ // Create dummy sky map cube
	629	+ GSkyMap map("CAR","GAL",0.0,0.0,1.0,1.0,1,1,ebounds.size());
	630	+
	631	+ // Create event cube
	632	+ GCTAEventCube cube(map, ebounds, exposure.gti());
	633	+
	634	+ // Create CTA observation
	635	+ GCTAObservation cta;
	636	+ cta.events(cube);
	637	+ cta.response(exposure, psf, background);
	638	+
	639	+ // Append observation to container
	640	+ m_obs.append(cta);
	641	+
	642	+ // Signal that we are in binned mode
	643	+ m_binned = true;
	644	+
	645	+ // Signal that we appended a cube
	646	+ m_append_cube = true;
	647	+
	648	+ } // endif: cube response information was available
	649	+
	650	+ // ... otherwise build an observation from IRF response information
	651	+ else {
	652	+
	653	+ // Create CTA observation
	654	+ GCTAObservation cta = create_cta_obs();
	655	+
	656	+ // Set response
	657	+ set_obs_response(&cta);
	658	+
	659	+ // Append observation to container
	660	+ m_obs.append(cta);
	661	+
	662	+ }
	663	+
	664	+ } // endif: filename was not valid
	665	+
	666	+ // ... otherwise we have a file name
	667	+ else {
	668	+
	669	+ // If file is a FITS file then create an empty CTA observation
	670	+ // and load file into observation
	671	+ if (GFilename(filename).is_fits()) {
	672	+
	673	+ // Allocate empty CTA observation
	674	+ GCTAObservation cta;
	675	+
	676	+ // Load data
	677	+ cta.load(filename);
	678	+
	679	+ // Set response
	680	+ set_obs_response(&cta);
	681	+
	682	+ // Append observation to container
	683	+ m_obs.append(cta);
	684	+
	685	+ // Signal that no XML file should be used for storage
	686	+ m_use_xml = false;
	687	+
	688	+ }
	689	+
	690	+ // ... otherwise load file into observation container
	691	+ else {
	692	+
	693	+ // Load observations from XML file
	694	+ m_obs.load(filename);
	695	+
	696	+ // For all observations that have no response, set the response
	697	+ // from the task parameters
	698	+ set_response(m_obs);
	699	+
	700	+ // Set observation boundary parameters (emin, emax, rad)
	701	+ set_obs_bounds(m_obs);
	702	+
	703	+ // Signal that XML file should be used for storage
	704	+ m_use_xml = true;
	705	+
	706	+ } // endelse: file was an XML file
	707	+
	708	+ }
	709	+
	710	+ // Return
	711	+ return;
	712	+
	713	+}
	714	+
	715	+
	716	+/*********************************************************************//
	717	+ * @brief Fill model into model cube
	718	+ *
	719	+ * @param[in] obs CTA observation.
	720	+ *
	721	+ * Adds the expected number of events for a given observation to the events
	722	+ * that are already found in the model cube. The method also updates the
	723	+ * GTI of the model cube so that cube GTI is a list of the GTIs of all
	724	+ * observations that were used to generate the model cube.
	725	+ ***************************************************************************/
	726	+void ctmodel::fill_cube(const GCTAObservation* obs)
	727	+{
	728	+ // Get references to GTI and energy boundaries for the event list
	729	+ const GGti& gti = obs->events()->gti();
	730	+ const GEbounds& obs_ebounds = obs->ebounds();
	731	+
	732	+ // Get cube energy boundaries
	733	+ const GEbounds& cube_ebounds = m_cube.ebounds();
	734	+
	735	+ // Get counts cube usage flags
	736	+ std::vector<bool> usage = cube_layer_usage(cube_ebounds, obs_ebounds);
	737	+
	738	+ // Initialise empty, invalid RoI
	739	+ GCTARoi roi;
	740	+
	741	+ // Retrieve RoI in case we have an unbinned observation
	742	+ if (obs->eventtype() == "EventList") {
	743	+ roi = obs->roi();
	744	+ }
	745	+
	746	+ // Initialise statistics
	747	+ double sum = 0.0;
	748	+ int num_outside_ebds = 0;
	749	+ int num_outside_roi = 0;
	750	+
	751	+ // Setup cube GTIs for this observation
	752	+ m_cube.gti(obs->events()->gti());
	753	+
	754	+ // Loop over all cube bins
	755	+ for (int i = 0; i < m_cube.size(); ++i) {
	756	+
	757	+ // Get cube bin
	758	+ GCTAEventBin* bin = m_cube[i];
	759	+
	760	+ // Determine counts cube energy bin
	761	+ int iebin = cube_ebounds.index(bin->energy());
	762	+
	763	+ // Skip bin if the corresponding counts cube energy bin is not fully
	764	+ // contained in the event list energy range. This avoids having
	765	+ // partially filled bins.
	766	+ if (!usage[iebin]) {
	767	+ num_outside_ebds++;
	768	+ continue;
	769	+ }
	770	+
	771	+ // If RoI is valid then skip bin if it is outside the RoI of the
	772	+ // observation
	773	+ if (roi.is_valid() && !roi.contains(*bin)) {
	774	+ num_outside_roi++;
	775	+ continue;
	776	+ }
	777	+
	778	+ // Get actual bin value
	779	+ double value = bin->counts();
	780	+
	781	+ // Compute model value for cube bin
	782	+ double model = m_obs.models().eval(bin, obs) * bin->size();
	783	+
	784	+ // Add model to actual value
	785	+ value += model;
	786	+ sum += model;
	787	+
	788	+ // Store value
	789	+ bin->counts(value);
	790	+
	791	+ } // endfor: looped over all cube bins
	792	+
	793	+ // Append GTIs of observation to list of GTIs
	794	+ m_gti.extend(gti);
	795	+
	796	+ // Update GTIs
	797	+ m_cube.gti(m_gti);
	798	+
	799	+ // Log filling results
	800	+ log_value(NORMAL, "Model events in cube", sum);
	801	+ log_value(NORMAL, "Bins outside energy range", num_outside_ebds);
	802	+ log_value(NORMAL, "Bins outside RoI", num_outside_roi);
	803	+
	804	+ // Write model cube into header
	805	+ log_header2(EXPLICIT, "Model cube");
	806	+ log_string(EXPLICIT, m_cube.print(m_chatter));
	807	+
	808	+ // Return
	809	+ return;
	810	+}
	811	+
...	...

src/ctprob/ctprob.hpp 0 → 100644

View file @3d85a1f

	1	+/***************************************************************************
	2	+ * ctmodel - Model cube generation tool *
	3	+ * ----------------------------------------------------------------------- *
	4	+ * copyright (C) 2012-2016 by Juergen Knoedlseder *
	5	+ * ----------------------------------------------------------------------- *
	6	+ * *
	7	+ * This program is free software: you can redistribute it and/or modify *
	8	+ * it under the terms of the GNU General Public License as published by *
	9	+ * the Free Software Foundation, either version 3 of the License, or *
	10	+ * (at your option) any later version. *
	11	+ * *
	12	+ * This program is distributed in the hope that it will be useful, *
	13	+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
	14	+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
	15	+ * GNU General Public License for more details. *
	16	+ * *
	17	+ * You should have received a copy of the GNU General Public License *
	18	+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
	19	+ * *
	20	+ ***************************************************************************/
	21	+/**
	22	+ * @file ctmodel.hpp
	23	+ * @brief Model cube generation tool definition
	24	+ * @author Juergen Knoedlseder
	25	+ */
	26	+
	27	+#ifndef CTMODEL_HPP
	28	+#define CTMODEL_HPP
	29	+
	30	+/* __ Includes ___________________________________________________________ */
	31	+#include "GammaLib.hpp"
	32	+#include "GCTALib.hpp"
	33	+#include "ctobservation.hpp"
	34	+
	35	+/* __Definitions _________________________________________________________ */
	36	+#define CTMODEL_NAME "ctmodel"
	37	+#define CTMODEL_VERSION "1.1.0"
	38	+
	39	+
	40	+/*********************************************************************//
	41	+ * @class ctmodel
	42	+ *
	43	+ * @brief Model cube generation tool
	44	+ *
	45	+ * This class creates counts model map(s). The definition of the counts model
	46	+ * can be taken from a predefined observation container, from the counts maps
	47	+ * found in an observation definition XML file, from an individual counts
	48	+ * map, or from the parameter definition.
	49	+ *
	50	+ * Results are stored in an observation container that can be written to disk
	51	+ * in form of FITS files (model maps) and an updated observation definition
	52	+ * XML file.
	53	+ ***************************************************************************/
	54	+class ctmodel : public ctobservation {
	55	+
	56	+public:
	57	+ // Constructors and destructors
	58	+ ctmodel(void);
	59	+ explicit ctmodel(const GObservations& obs);
	60	+ ctmodel(int argc, char *argv[]);
	61	+ ctmodel(const ctmodel& app);
	62	+ virtual ~ctmodel(void);
	63	+
	64	+ // Operators
	65	+ ctmodel& operator=(const ctmodel& app);
	66	+
	67	+ // Methods
	68	+ void clear(void);
	69	+ void run(void);
	70	+ void save(void);
	71	+ void publish(const std::string& name = "");
	72	+ const GCTAEventCube& cube(void) const;
	73	+ void cube(const GCTAEventCube& cube);
	74	+ void models(const GModels& models);
	75	+
	76	+protected:
	77	+ // Protected methods
	78	+ void init_members(void);
	79	+ void copy_members(const ctmodel& app);
	80	+ void free_members(void);
	81	+ void get_parameters(void);
	82	+ void get_obs(void);
	83	+ void fill_cube(const GCTAObservation* obs);
	84	+ bool has_cube(void) const;
	85	+
	86	+ // User parameters
	87	+ GFilename m_outcube; //!< Output model cube
	88	+ bool m_apply_edisp; //!< Apply energy dispersion?
	89	+ bool m_publish; //!< Publish model cube?
	90	+ GChatter m_chatter; //!< Chattiness
	91	+
	92	+ // Protected members
	93	+ GCTAEventCube m_cube; //!< Model cube
	94	+ GGti m_gti; //!< Model cube GTIs
	95	+ bool m_has_cube; //!< Signal if cube has been set or loaded
	96	+ bool m_append_cube; //!< Signal that cube should be appended
	97	+ bool m_binned; //!< Signals that we are in binned mode
	98	+};
	99	+
	100	+
	101	+/*********************************************************************//
	102	+ * @brief Return model cube
	103	+ *
	104	+ * @return Reference to model cube.
	105	+ ***************************************************************************/
	106	+inline
	107	+const GCTAEventCube& ctmodel::cube(void) const
	108	+{
	109	+ return m_cube;
	110	+}
	111	+
	112	+
	113	+/*********************************************************************//
	114	+ * @brief Signal if cube has been set or loaded
	115	+ *
	116	+ * @return True if cube has been set or loaded.
	117	+ ***************************************************************************/
	118	+inline
	119	+bool ctmodel::has_cube(void) const
	120	+{
	121	+ return m_has_cube;
	122	+}
	123	+
	124	+
	125	+/*********************************************************************//
	126	+ * @brief Set models
	127	+ *
	128	+ * @param[in] models Model container.
	129	+ *
	130	+ * Set model container that should be used for model generation.
	131	+ ***************************************************************************/
	132	+inline
	133	+void ctmodel::models(const GModels& models)
	134	+{
	135	+ m_obs.models(models);
	136	+ return;
	137	+}
	138	+
	139	+#endif /* CTMODEL_HPP */
...	...

src/ctprob/ctprob.par 0 → 100644

View file @3d85a1f

	1	+#############################################################################
	2	+# #
	3	+# CTA Science Analysis Tools #
	4	+# #
	5	+# ------------------------------------------------------------------------- #
	6	+# #
	7	+# File: ctmodel.par #
	8	+# #
	9	+# Author: Juergen Knoedlseder #
	10	+# Jurgen.Knodlseder@irap.omp.eu #
	11	+# IRAP #
	12	+# #
	13	+# Purpose: Parameter file for the ctmodel application #
	14	+# #
	15	+#############################################################################
	16	+#
	17	+# General parameters
	18	+#===================
	19	+inobs, f, a, "NONE",,, "Input event list, counts cube or observation definition XML file"
	20	+incube, f, a, "NONE",,, "Input counts cube file to extract model cube definition"
	21	+inmodel, f, a, "$CTOOLS/share/models/crab.xml",,, "Input model definition XML file"
	22	+expcube, f, a, "NONE",,, "Input exposure cube file (only needed for stacked analysis)"
	23	+psfcube, f, a, "NONE",,, "Input PSF cube file (only needed for stacked analysis)"
	24	+edispcube,f, a, "NONE",,, "Input energy dispersion cube file (only needed for stacked analysis)"
	25	+bkgcube, f, a, "NONE",,, "Input background cube file (only needed for stacked analysis)"
	26	+caldb, s, a, "prod2",,, "Calibration database"
	27	+irf, s, a, "South_0.5h",,, "Instrument response function"
	28	+edisp, b, h, no,,, "Apply energy dispersion?"
	29	+outcube, f, a, "modcube.fits",,, "Output model cube file"
	30	+
	31	+#
	32	+# Observation parameters
	33	+#=======================
	34	+ra, r, a, 83.63,0,360, "RA of pointing (degrees)"
	35	+dec, r, a, 22.01,-90,90, "Dec of pointing (degrees)"
	36	+rad, r, a, 5.0,0,180, "Radius of FOV (degrees)"
	37	+tmin, r, a, 0.0,,, "Start time (MET in s)"
	38	+tmax, r, a, 1800.0,,, "End time (MET in s)"
	39	+deadc, r, h, 0.98,0,1, "Average deadtime correction factor"
	40	+
	41	+#
	42	+# Energy binning parameters
	43	+#==========================
	44	+ebinalg, s, h, "LOG",FILE\|LIN\|LOG,, "Algorithm for defining energy bins"
	45	+emin, r, a, 0.1,,, "Start value for first energy bin in TeV"
	46	+emax, r, a, 100.0,,, "Stop value for last energy bin in TeV"
	47	+enumbins, i, a, 20,,, "Number of energy bins"
	48	+ebinfile, f, h, "NONE",,, "Name of the file containing the energy bin definition"
	49	+
	50	+#
	51	+# Spatial binning parameters
	52	+#===========================
	53	+usepnt, b, h, no,,, "Use pointing instead of xref/yref parameters?"
	54	+nxpix, i, a, 200,,, "Size of the X axis in pixels"
	55	+nypix, i, a, 200,,, "Size of the Y axis in pixels"
	56	+binsz, r, a, 0.02,,, "Image scale (in degrees/pixel)"
	57	+coordsys, s, a, "CEL", CEL\|GAL,,"Coordinate system (CEL - celestial, GAL - galactic)"
	58	+xref, r, a, 83.63,0,360, "First coordinate of image center in degrees (RA or galactic l)"
	59	+yref, r, a, 22.01,-90,90, "Second coordinate of image center in degrees (DEC or galactic b)"
	60	+proj, s, a, "CAR", AIT\|AZP\|CAR\|MER\|MOL\|STG\|TAN,, "Projection method"
	61	+
	62	+#
	63	+# Standard parameters
	64	+#====================
	65	+publish, b, h, no,,, "Publish model cube on VO Hub?"
	66	+chatter, i, h, 2,0,4, "Chattiness of output"
	67	+clobber, b, h, yes,,, "Overwrite existing output files with new output files?"
	68	+debug, b, h, no,,, "Debugging mode activated"
	69	+mode, s, h, "ql",,, "Mode of automatic parameters"
	70	+logfile, f, h, "ctmodel.log",,, "Log filename"
...	...

src/ctprob/main.cpp 0 → 100644

View file @3d85a1f

	1	+/***************************************************************************
	2	+ * ctmodel - CTA counts model tool *
	3	+ * ----------------------------------------------------------------------- *
	4	+ * copyright (C) 2012-2016 by Juergen Knoedlseder *
	5	+ * ----------------------------------------------------------------------- *
	6	+ * *
	7	+ * This program is free software: you can redistribute it and/or modify *
	8	+ * it under the terms of the GNU General Public License as published by *
	9	+ * the Free Software Foundation, either version 3 of the License, or *
	10	+ * (at your option) any later version. *
	11	+ * *
	12	+ * This program is distributed in the hope that it will be useful, *
	13	+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
	14	+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
	15	+ * GNU General Public License for more details. *
	16	+ * *
	17	+ * You should have received a copy of the GNU General Public License *
	18	+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
	19	+ * *
	20	+ ***************************************************************************/
	21	+/**
	22	+ * @file ctmodel/main.cpp
	23	+ * @brief CTA counts model tool main code
	24	+ * @author Juergen Knoedlseder
	25	+ */
	26	+
	27	+/* __ Includes ___________________________________________________________ */
	28	+#ifdef HAVE_CONFIG_H
	29	+#include <config.h>
	30	+#endif
	31	+#include "support.hpp"
	32	+#include "ctmodel.hpp"
	33	+
	34	+
	35	+/*********************************************************************//
	36	+ * @brief Main entry point of application
	37	+ *
	38	+ * @param[in] argc Number of command line arguments.
	39	+ * @param[in] argv Command line arguments.
	40	+ *
	41	+ * This is the main entry point of the ctmodel application. It allocates a
	42	+ * ctmodel object and executes the application. Any exceptions that occur
	43	+ * will be catched and corresponding error messages written in the
	44	+ * application logger and into the standard output.
	45	+ ***************************************************************************/
	46	+int main (int argc, char *argv[])
	47	+{
	48	+ // Initialise return code
	49	+ int rc = 1;
	50	+
	51	+ // Initialise pointer on application
	52	+ ctmodel* application = NULL;
	53	+
	54	+ // Execute application
	55	+ try {
	56	+
	57	+ // Create instance of application
	58	+ application = new ctmodel(argc, argv);
	59	+
	60	+ // Execute ctool
	61	+ rc = execute_ctool(application);
	62	+
	63	+ }
	64	+
	65	+ catch (std::exception &e) {
	66	+
	67	+ // Report exception
	68	+ report_ctool_failure("ctmodel", e.what());
	69	+
	70	+ }
	71	+
	72	+ // Delete application
	73	+ delete application;
	74	+
	75	+ // Return
	76	+ return rc;
	77	+}
...	...