diff --git a/ChangeLog b/ChangeLog index a170c31..2ff6f52 100644 --- a/ChangeLog +++ b/ChangeLog @@ -11,6 +11,7 @@ Add ctulimit tool (#1363) Allow ctmodel usage for stacked analysis from command line Correct rate check in ctobssim for map cubes (#1452) + Handle background cubes for stacked analysis (#1371) 2015-01-08 Juergen Knoedlseder diff --git a/NEWS b/NEWS index 94bcb5b..d133223 100644 --- a/NEWS +++ b/NEWS @@ -2,7 +2,7 @@ New Features and Important Changes in ctools 1.0.0 Author(s) : Juergen Knoedlseder -23 March 2015 +28 March 2015 1. Introduction @@ -50,12 +50,14 @@ None 8. ctlike - Maximum likelihood fitting -------------------------------------- -None +The bkgcube parameter has been added to allow specification of a background +cube for stacked analysis. 9. cttsmap - Generation of Test Statistics map ---------------------------------------------- -None +The bkgcube parameter has been added to allow specification of a background +cube for stacked analysis. 10. ctskymap - CTA sky mapping @@ -71,6 +73,9 @@ the original IRF information. This makes generation of model cubes substantially faster. The old method of accessing the IRF response for model computation is still possible. +The bkgcube parameter has been added to allow specification of a background +cube for stacked analysis. + 12. ctexpcube - CTA exposure cube generation -------------------------------------------- @@ -94,13 +99,17 @@ None 16. ctbutterfly - Generate butterfly diagrams --------------------------------------------- -None +The bkgcube parameter has been added to allow specification of a background +cube for stacked analysis. 17. ctulimit - Compute upper limit ---------------------------------- Added an initial version of the upper limit computation tool. +The bkgcube parameter has been added to allow specification of a background +cube for stacked analysis. + 18. Examples ------------ diff --git a/doc/source/reference_manual/cspull.rst b/doc/source/reference_manual/cspull.rst index 7d0d69e..8e96bf2 100644 --- a/doc/source/reference_manual/cspull.rst +++ b/doc/source/reference_manual/cspull.rst @@ -26,9 +26,15 @@ General parameters ``outfile [file]`` ASCII file containing the individual pull values. -``ntrials [integer]`` - Number of samples for generating the pull distribution. - +``expcube [file]`` + Exposure cube file (only needed for stacked analysis). + +``psfcube [file]`` + PSF cube file (only needed for stacked analysis). + +``bkgcube [file]`` + Background cube file (only needed for stacked analysis). + ``caldb [string]`` Calibration database. @@ -38,6 +44,9 @@ General parameters ``(edisp = no) [boolean]`` Apply energy dispersion to response computation. +``ntrials [integer]`` + Number of samples for generating the pull distribution. + ``ra [real]`` Right Ascension of CTA pointing (J2000, in degrees). diff --git a/doc/source/reference_manual/ctbutterfly.rst b/doc/source/reference_manual/ctbutterfly.rst index 2fc3abf..fd3a334 100644 --- a/doc/source/reference_manual/ctbutterfly.rst +++ b/doc/source/reference_manual/ctbutterfly.rst @@ -25,12 +25,15 @@ General parameters ``srcname [string]`` Name of model component for which butterfly should be computed. -``expcube = NONE [file]`` +``expcube [file]`` Exposure cube file (only needed for stacked analysis). -``psfcube = NONE [file]`` +``psfcube [file]`` PSF cube file (only needed for stacked analysis). +``bkgcube [file]`` + Background cube file (only needed for stacked analysis). + ``caldb [string]`` Calibration database. diff --git a/doc/source/reference_manual/ctlike.rst b/doc/source/reference_manual/ctlike.rst index 3205405..348cad2 100644 --- a/doc/source/reference_manual/ctlike.rst +++ b/doc/source/reference_manual/ctlike.rst @@ -34,6 +34,9 @@ General parameters ``psfcube [file]`` PSF cube file (only needed for stacked analysis). +``bkgcube [file]`` + Background cube file (only needed for stacked analysis). + ``caldb [string]`` Calibration database. diff --git a/doc/source/reference_manual/ctmodel.rst b/doc/source/reference_manual/ctmodel.rst index 4d5da11..b5eb7ef 100644 --- a/doc/source/reference_manual/ctmodel.rst +++ b/doc/source/reference_manual/ctmodel.rst @@ -43,6 +43,9 @@ General parameters ``psfcube [file]`` PSF cube file (only needed for stacked analysis). +``bkgcube [file]`` + Background cube file (only needed for stacked analysis). + ``caldb [string]`` Calibration database. diff --git a/doc/source/reference_manual/cttsmap.rst b/doc/source/reference_manual/cttsmap.rst index 09fecf8..d79ff70 100644 --- a/doc/source/reference_manual/cttsmap.rst +++ b/doc/source/reference_manual/cttsmap.rst @@ -31,6 +31,9 @@ General parameters ``psfcube [file]`` PSF cube file (only needed for stacked analysis). +``bkgcube [file]`` + Background cube file (only needed for stacked analysis). + ``caldb [string]`` Calibration database. diff --git a/doc/source/reference_manual/ctulimit.rst b/doc/source/reference_manual/ctulimit.rst index 55d0d62..df6df5b 100644 --- a/doc/source/reference_manual/ctulimit.rst +++ b/doc/source/reference_manual/ctulimit.rst @@ -24,12 +24,15 @@ General parameters ``srcname [string]`` Name of model component for which upper limit should be computed. -``expcube = NONE [file]`` +``expcube [file]`` Exposure cube file (only needed for stacked analysis). -``psfcube = NONE [file]`` +``psfcube [file]`` PSF cube file (only needed for stacked analysis). +``bkgcube [file]`` + Background cube file (only needed for stacked analysis). + ``caldb [string]`` Calibration database. diff --git a/pyext/ctbkgcube.i b/pyext/ctbkgcube.i index 33e1f11..c8631b0 100644 --- a/pyext/ctbkgcube.i +++ b/pyext/ctbkgcube.i @@ -1,7 +1,7 @@ /*************************************************************************** * ctbkgcube - Background cube generation tool * * ----------------------------------------------------------------------- * - * copyright (C) 2014 by Chia-Chun Lu * + * copyright (C) 2014-2015 by Chia-Chun Lu * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -44,11 +44,11 @@ public: virtual ~ctbkgcube(void); // Methods - void clear(void); - void run(void); - void save(void); - const GCTAEventCube& cube(void) const; - const GModels& models(void) const; + void clear(void); + void run(void); + void save(void); + const GCTACubeBackground& cube(void) const; + const GModels& models(void) const; }; /***********************************************************************//** diff --git a/scripts/cspull.py b/scripts/cspull.py index f2972b6..9abf3ae 100755 --- a/scripts/cspull.py +++ b/scripts/cspull.py @@ -98,8 +98,11 @@ class cspull(ctools.cscript): pars.append(gammalib.GApplicationPar("outfile","f","a","pull.dat","","","Output file name")) pars.append(gammalib.GApplicationPar("expcube","s","a","NONE","","","Exposure cube file (only needed for stacked analysis)")) pars.append(gammalib.GApplicationPar("psfcube","s","a","NONE","","","PSF cube file (only needed for stacked analysis)")) + pars.append(gammalib.GApplicationPar("bkgcube","s","a","NONE","","","Background cube file (only needed for stacked analysis)")) pars.append(gammalib.GApplicationPar("caldb","s","a","dummy","","","Calibration database")) pars.append(gammalib.GApplicationPar("irf","s","a","cta_dummy_irf","","","Instrument response function")) + pars.append(gammalib.GApplicationPar("edisp","b","h","no","","","Apply energy dispersion?")) + pars.append(gammalib.GApplicationPar("ntrials","i","a","10","","","Number of trials")) pars.append(gammalib.GApplicationPar("ra","r","a","83.6331","0","360","RA of pointing (deg)")) pars.append(gammalib.GApplicationPar("dec","r","a","22.0145","-90","90","Dec of pointing (deg)")) pars.append(gammalib.GApplicationPar("emin","r","a","0.1","","","Lower energy limit (TeV)")) @@ -107,14 +110,12 @@ class cspull(ctools.cscript): pars.append(gammalib.GApplicationPar("tmin","r","h","0.0","","","Start time (MET in s)")) pars.append(gammalib.GApplicationPar("tmax","r","a","1800.0","","","Duration (in s)")) pars.append(gammalib.GApplicationPar("enumbins","i","a","0","","","Number of energy bins (0=unbinned)")) - pars.append(gammalib.GApplicationPar("npix","i","h","200","","","Number of pixels for binned")) - pars.append(gammalib.GApplicationPar("binsz","r","h","0.05","","","Pixel size for binned (deg/pixel)")) + pars.append(gammalib.GApplicationPar("npix","i","a","200","","","Number of pixels for binned")) + pars.append(gammalib.GApplicationPar("binsz","r","a","0.05","","","Pixel size for binned (deg/pixel)")) pars.append(gammalib.GApplicationPar("deadc","r","h","0.95","0","1","Deadtime correction factor")) pars.append(gammalib.GApplicationPar("rad","r","h","5.0","0","180","Radius of ROI (deg)")) pars.append(gammalib.GApplicationPar("pattern","s","h","single","","","Observation pattern (single/four)")) pars.append(gammalib.GApplicationPar("offset","r","h","1.5","","","Observation pattern offset (deg)")) - pars.append(gammalib.GApplicationPar("ntrials","i","a","10","","","Number of trials")) - pars.append(gammalib.GApplicationPar("edisp","b","h","no","","","Apply energy dispersion?")) pars.append_standard() pars.append(gammalib.GApplicationPar("logfile","f","h","cspull.log","","","Log filename")) pars.save(parfile) @@ -140,7 +141,7 @@ class cspull(ctools.cscript): # Get number of energy bins self.m_enumbins = self["enumbins"].integer() - + # Read parameters for binned if requested if not self.m_enumbins == 0: self.m_npix = self["npix"].integer() diff --git a/src/ctbkgcube/ctbkgcube.cpp b/src/ctbkgcube/ctbkgcube.cpp index da6a057..acdbb86 100644 --- a/src/ctbkgcube/ctbkgcube.cpp +++ b/src/ctbkgcube/ctbkgcube.cpp @@ -34,7 +34,6 @@ /* __ Method name definitions ____________________________________________ */ #define G_GET_PARAMETERS "ctbkgcube::get_parameters()" -#define G_FILL_CUBE "ctbkgcube::fill_cube(GCTAObservation*)" /* __ Debug definitions __________________________________________________ */ @@ -209,9 +208,6 @@ void ctbkgcube::run(void) log.cout(true); } - // Initialise run parameters - m_livetime = 0.0; - // Get task parameters get_parameters(); @@ -248,15 +244,14 @@ void ctbkgcube::run(void) // Copy models from observation container and reset output model // container - m_bkgmdl = m_obs.models(); + GModels models_orig = m_obs.models(); + m_bkgmdl = m_obs.models(); m_outmdl.clear(); - m_cube_model = -1; // Remove all models that are not CTA background models from the // container and put all removed components in the output // container int num = m_bkgmdl.size(); - for (int i = num-1; i >= 0; --i) { // Flag removal @@ -309,75 +304,28 @@ void ctbkgcube::run(void) log.header1("Generate background cube"); } - // Initialise observation counter - int n_observations = 0; - - // Loop over all observations in the container - for (int i = 0; i < m_obs.size(); ++i) { - - // Get CTA observation - GCTAObservation* obs = dynamic_cast(m_obs[i]); - - // Continue only if observation is a CTA observation - if (obs != NULL) { - - // Write header for observation - if (logTerse()) { - if (obs->name().length() > 1) { - log.header3("Observation "+obs->name()); - } - else { - log.header3("Observation "+gammalib::str(i)); - } - } - - // Increment number of observations - n_observations++; - - // Fill the cube - fill_cube(obs); - - } // endif: CTA observation found + // Assign background models to container + m_obs.models(m_bkgmdl); - } // endfor: looped over observations + // Fill background cube from observations + m_background.fill(m_obs); - // Log results - if (logTerse()) { - log << gammalib::parformat("Accumulated livetime"); - log << m_livetime << " sec" << std::endl; - } + // Create a background model for the output background cube and append + // that model to the input model in place of the original + // background models + // TODO: We might think of creating the spectral model via user parameter + GModelSpectralPlaw spectral(1.0, 0.0, GEnergy(1.0, "TeV")); + GCTAModelCubeBackground model(spectral); - // Re-normalize cube to get units of counts/MeV/s/sr - if (m_livetime > 0.0) { - - // Loop over all bins in background cube and divide the content - // by the total livetime. - for (int i = 0; i < m_bkgcube.size(); ++i) { - GCTAEventBin* bin = m_bkgcube[i]; - double rate = bin->counts() / m_livetime; - bin->counts(rate); - } + // Set model name + model.name("BackgroundModel"); - } // endif: livetime was positive + // Set model instrument + // TODO: Account for possibility to have observations from different + // IACTs in the same container + model.instruments("CTA,HESS,MAGIC,VERITAS"); - // Log results - if (logTerse()) { - log << gammalib::parformat("Number of observations"); - log << n_observations << std::endl; - } - - // Create a background model for output background cube and append - // that model to the input model in place of the original - // background model - GEnergies energies; - for (int i = 0; i < m_bkgcube.ebins(); ++i) { - energies.append(m_bkgcube.energy(i)); - } - GModelSpatialDiffuseCube spatial(m_bkgcube.map(), energies, 1.0); - GModelSpectralPlaw spectral(1.0, 0.0, GEnergy(1.0, "TeV")); - GCTAModelCubeBackground model(spatial, spectral); - model.name("ctbkgcube default background model"); - m_cube_model = m_outmdl.size(); // Store the slot number for save() + // Append model to output container m_outmdl.append(model); // Log output model @@ -387,6 +335,9 @@ void ctbkgcube::run(void) log << m_outmdl << std::endl; } + // Recover original models + m_obs.models(models_orig); + // Return return; } @@ -405,47 +356,19 @@ void ctbkgcube::save(void) log.header1("Save background cube"); } - // Create energies container from energy boundaries - GEnergies energies; - for (int i = 0; i < m_ebounds.size(); ++i) { - energies.append(m_ebounds.elogmean(i)); - } - // Get output filenames std::string outfile = (*this)["outcube"].filename(); std::string outmodel = (*this)["outmodel"].filename(); - // Create empty FITS file - GFits fits; - - // Write background cube - m_bkgcube.map().write(fits); - - // Write energies - energies.write(fits); - - // Save FITS file - fits.saveto(outfile, clobber()); + // Save background cube + m_background.save(outfile, clobber()); // Write output models if filename is valid if ((outmodel.length() > 0) && (gammalib::tolower(outmodel) != "none")) { - // Set filename of map cube. - if (m_cube_model >= 0) { - GCTAModelCubeBackground* model = - dynamic_cast(m_outmdl[m_cube_model]); - if (model != NULL) { - GModelSpatialDiffuseCube* spatial = - dynamic_cast(model->spatial()); - if (spatial != NULL) { - spatial->filename(outfile); - } - } - } - - // Save output model container + // Save output model for binned analyses m_outmdl.save(outmodel); - + } // Return @@ -468,12 +391,9 @@ void ctbkgcube::init_members(void) m_outcube.clear(); m_outmodel.clear(); m_obs.clear(); - m_bkgcube.clear(); + m_background.clear(); m_bkgmdl.clear(); m_outmdl.clear(); - m_ebounds.clear(); - m_livetime = 0.0; - m_cube_model = -1; // Return return; @@ -491,12 +411,9 @@ void ctbkgcube::copy_members(const ctbkgcube& app) m_outmodel = app.m_outmodel; m_outcube = app.m_outcube; m_obs = app.m_obs; - m_bkgcube = app.m_bkgcube; + m_background = app.m_background; m_bkgmdl = app.m_bkgmdl; m_outmdl = app.m_outmdl; - m_ebounds = app.m_ebounds; - m_livetime = app.m_livetime; - m_cube_model = app.m_cube_model; // Return return; @@ -536,22 +453,6 @@ void ctbkgcube::get_parameters(void) } // endif: there was no observation in the container - - // If there are no models associated with the observations then load now - // the model definition from the XML file - if (m_obs.models().size() == 0) { - if ((*this)["inmodel"].is_undefined()) { - std::string msg = "No model definition XML file specified. " - "Please set the \"inmodel\" parameter to the " - "XML file that contains the background model " - "definition."; - throw GException::invalid_value(G_GET_PARAMETERS, msg); - } - std::string inmodel = (*this)["inmodel"].filename(); - GModels models(inmodel); - m_obs.models(models); - } - // Get the incube filename std::string incube = (*this)["incube"].filename(); @@ -560,28 +461,35 @@ void ctbkgcube::get_parameters(void) (gammalib::strip_whitespace(incube) == "")) { // Create an event cube based on task parameters - m_bkgcube = create_cube(m_obs); + GCTAEventCube cube = create_cube(m_obs); + + // Define background cube + m_background = GCTACubeBackground(cube); } // endif: filename was not valid // ... otherwise setup the background cube from the counts map else { - - // Load event cube from filename - m_bkgcube.load(incube); - // Loop over all bins in background cube - for (int i = 0; i < m_bkgcube.size(); ++i) { - // Get event bin - GCTAEventBin* bin = m_bkgcube[i]; - bin->counts(0.0); - - } + // Define background cube from file + m_background = GCTACubeBackground(incube); } // endelse: cube was loaded from file - // Get energy definition - m_ebounds = m_bkgcube.ebounds(); + // If there are no models associated with the observations then load now + // the model definition from the XML file + if (m_obs.models().size() == 0) { + if ((*this)["inmodel"].is_undefined()) { + std::string msg = "No model definition XML file specified. " + "Please set the \"inmodel\" parameter to the " + "XML file that contains the background model " + "definition."; + throw GException::invalid_value(G_GET_PARAMETERS, msg); + } + std::string inmodel = (*this)["inmodel"].filename(); + GModels models(inmodel); + m_obs.models(models); + } // Read output filenames (if needed) if (read_ahead()) { @@ -592,92 +500,3 @@ void ctbkgcube::get_parameters(void) // Return return; } - -/***********************************************************************//** - * @brief Generate background cube - * - * @param[in] obs Pointer to CTA observation. - * - * Fills the background cube with the model value for a CTA observation. - * This method will not need to load the event file as it gets the event - * ROI and GTI directly from the observation. - ***************************************************************************/ -void ctbkgcube::fill_cube(GCTAObservation* obs) -{ - - // Continue only if observation pointer is valid - if (obs != NULL) { - - // Initialise statistics - double sum = 0.0; - int n_bins_inside = 0; - int n_bins_outside = 0; - - // Extract region of interest from CTA observation - GCTARoi roi = obs->roi(); - - // Set GTI of actual observations as the GTI of the event cube - m_bkgcube.gti(obs->gti()); - - // Get observation livetime - double livetime = obs->livetime(); - - // Loop over all bins in background cube - for (int i = 0; i < m_bkgcube.size(); ++i) { - - // Get event bin - GCTAEventBin* bin = m_bkgcube[i]; - - // Continue only if binned in contained in ROI - if (roi.contains(*bin)) { - - // Compute model value for event bin. The model value is - // given in counts/MeV/s/sr. - double model = m_bkgmdl.eval(*bin, *obs); - - // Compute number events - sum += model * bin->size(); - - // Multiply by livetime to get the correct weighting for - // each observation. We divide by the total livetime later - // to get the background model in units of counts/MeV/s/sr. - model *= livetime; - - // Add existing number of counts - model += bin->counts(); - - // Store cumulated value (units: counts/MeV/sr) - bin->counts(model); - - // Increment contained bin counter - n_bins_inside++; - - } // endif: bin was contained in RoI - - // ... otherwise count bin as an outsider - else { - n_bins_outside++; - } - - } // endfor: looped over all bins - - // Accumulate livetime - m_livetime += livetime; - - // Log results - if (logTerse()) { - log << gammalib::parformat("Bins within RoI"); - log << n_bins_inside << std::endl; - log << gammalib::parformat("Bins outside RoI"); - log << n_bins_outside << std::endl; - log << gammalib::parformat("Background events in cube"); - log << sum << std::endl; - log << gammalib::parformat("Cube livetime"); - log << livetime << " sec" << std::endl; - } - - } // endif: observation pointer was not valid - - // Return - return; -} diff --git a/src/ctbkgcube/ctbkgcube.hpp b/src/ctbkgcube/ctbkgcube.hpp index dae4042..9d64ef1 100644 --- a/src/ctbkgcube/ctbkgcube.hpp +++ b/src/ctbkgcube/ctbkgcube.hpp @@ -56,11 +56,11 @@ public: ctbkgcube& operator=(const ctbkgcube& app); // Methods - void clear(void); - void run(void); - void save(void); - const GCTAEventCube& cube(void) const; - const GModels& models(void) const; + void clear(void); + void run(void); + void save(void); + const GCTACubeBackground& cube(void) const; + const GModels& models(void) const; protected: // Protected methods @@ -71,31 +71,29 @@ protected: void fill_cube(GCTAObservation* obs); // Protected members - std::string m_outcube; //!< Filename of output cube - std::string m_outmodel; //!< Filename of output XML model - GObservations m_obs; //!< Observation container - GCTAEventCube m_bkgcube; //!< Background cube - GModels m_bkgmdl; //!< CTA background models - GModels m_outmdl; //!< Output models - GEbounds m_ebounds; //!< Energy boundaries - double m_livetime; //!< Total livetime for normalization - int m_cube_model; //!< Cube model index + std::string m_outcube; //!< Filename of output cube + std::string m_outmodel; //!< Filename of output XML model + GObservations m_obs; //!< Observation container + GCTACubeBackground m_background; //!< Background cube response + GModels m_bkgmdl; //!< CTA background models + GModels m_outmdl; //!< Output models }; /***********************************************************************//** - * @brief Return event cube containing background rate + * @brief Return background response cube containing background rate * - * @return Background event cube. + * @return Background response cube. ***************************************************************************/ inline -const GCTAEventCube& ctbkgcube::cube(void) const +const GCTACubeBackground& ctbkgcube::cube(void) const { // Return background cube - return (m_bkgcube); + return (m_background); } + /***********************************************************************//** * @brief Return background model container * diff --git a/src/ctbkgcube/ctbkgcube.par b/src/ctbkgcube/ctbkgcube.par index 272d150..400642d 100644 --- a/src/ctbkgcube/ctbkgcube.par +++ b/src/ctbkgcube/ctbkgcube.par @@ -18,10 +18,10 @@ # Observational parameters #=================== inobs, f, a, "NONE",,, "Input event list or observation definition file" -inmodel, f, a, "NONE",,, "Input (background) model XML file" +inmodel, f, a, "NONE",,, "Input model XML file" incube, f, a, "NONE",,, "Counts cube for background cube definition" outcube, f, a, "bkgcube.fits",,, "Output background cube file" -outmodel, f, a, "NONE",,, "Output (background) model XML file" +outmodel, f, a, "NONE",,, "Output model XML file" caldb, s, a, "dummy",,, "Calibration database" irf, s, a, "cta_dummy_irf",,, "Instrument response function" # diff --git a/src/ctbutterfly/ctbutterfly.par b/src/ctbutterfly/ctbutterfly.par index b823d03..2d5e224 100644 --- a/src/ctbutterfly/ctbutterfly.par +++ b/src/ctbutterfly/ctbutterfly.par @@ -22,6 +22,7 @@ inmodel, f, a, "$CTOOLS/share/models/crab.xml",,, "Source model" srcname, s, a, "Crab",,,"Source of interest" expcube, f, a, "NONE",,, "Exposure cube file (only needed for stacked analysis)" psfcube, f, a, "NONE",,, "PSF cube file (only needed for stacked analysis)" +bkgcube, f, a, "NONE",,, "Background cube file (only needed for stacked analysis)" caldb, s, a, "dummy",,, "Calibration database" irf, s, a, "cta_dummy_irf",,, "Instrument response function" outfile, f, a, "butterfly.txt",,, "Output ascii file" diff --git a/src/ctlike/ctlike.par b/src/ctlike/ctlike.par index 19dc39f..9cf0136 100644 --- a/src/ctlike/ctlike.par +++ b/src/ctlike/ctlike.par @@ -21,6 +21,7 @@ inobs, f, a, "events.fits",,, "Event list, counts cube or observation definit inmodel, f, a, "$CTOOLS/share/models/crab.xml",,, "Source model" expcube, f, a, "NONE",,, "Exposure cube file (only needed for stacked analysis)" psfcube, f, a, "NONE",,, "PSF cube file (only needed for stacked analysis)" +bkgcube, f, a, "NONE",,, "Background cube file (only needed for stacked analysis)" caldb, s, a, "dummy",,, "Calibration database" irf, s, a, "cta_dummy_irf",,, "Instrument response function" outmodel, f, a, "crab_results.xml",,, "Source model output file" diff --git a/src/ctmodel/ctmodel.par b/src/ctmodel/ctmodel.par index 9cd20a9..9506397 100644 --- a/src/ctmodel/ctmodel.par +++ b/src/ctmodel/ctmodel.par @@ -27,6 +27,7 @@ edisp, b, h, no,,, "Apply energy dispersion?" inobs, f, a, "NONE",,, "Event list, counts cube or observation definition file" expcube, f, a, "NONE",,, "Exposure cube file (only needed for stacked analysis)" psfcube, f, a, "NONE",,, "PSF cube file (only needed for stacked analysis)" +bkgcube, f, a, "NONE",,, "Background cube file (only needed for stacked analysis)" caldb, s, a, "dummy",,, "Calibration database" irf, s, a, "cta_dummy_irf",,, "Instrument response function" ra, r, a, 83.63,0,360, "RA of pointing (degrees)" diff --git a/src/cttsmap/cttsmap.par b/src/cttsmap/cttsmap.par index 862c497..ffe2699 100644 --- a/src/cttsmap/cttsmap.par +++ b/src/cttsmap/cttsmap.par @@ -22,6 +22,7 @@ inmodel, f, a, "$CTOOLS/share/models/crab.xml",,, "Source model" srcname, s, a, "Crab",,, "Test source" expcube, f, a, "NONE",,, "Exposure cube file (only needed for stacked analysis)" psfcube, f, a, "NONE",,, "PSF cube file (only needed for stacked analysis)" +bkgcube, f, a, "NONE",,, "Background cube file (only needed for stacked analysis)" caldb, s, a, "dummy",,, "Calibration database" irf, s, a, "cta_dummy_irf",,, "Instrument response function" outmap, f, a, "tsmap.fits",,, "Output Test Statistic map" diff --git a/src/ctulimit/ctulimit.par b/src/ctulimit/ctulimit.par index d3cf453..a577450 100644 --- a/src/ctulimit/ctulimit.par +++ b/src/ctulimit/ctulimit.par @@ -22,6 +22,7 @@ inmodel, f, a, "$CTOOLS/share/models/crab.xml",,, "Source model" srcname, s, a, "Crab",,, "Source of interest" expcube, f, a, "NONE",,, "Exposure cube file (only needed for stacked analysis)" psfcube, f, a, "NONE",,, "PSF cube file (only needed for stacked analysis)" +bkgcube, f, a, "NONE",,, "Background cube file (only needed for stacked analysis)" caldb, s, a, "dummy",,, "Calibration database" irf, s, a, "cta_dummy_irf",,, "Instrument response function" diff --git a/src/support/ctool.cpp b/src/support/ctool.cpp index 391bfc1..c169bc5 100644 --- a/src/support/ctool.cpp +++ b/src/support/ctool.cpp @@ -732,23 +732,26 @@ void ctool::set_obs_response(GCTAObservation* obs) // If the observation contains a counts cube, then first check whether // the expcube and psfcube parameters exist and are not NONE if (dynamic_cast(obs->events()) != NULL) { - if (has_par("expcube") && has_par("psfcube")) { + if (has_par("expcube") && has_par("psfcube") && has_par("bkgcube")) { // Get filenames std::string expcube = (*this)["expcube"].filename(); std::string psfcube = (*this)["psfcube"].filename(); + std::string bkgcube = (*this)["bkgcube"].filename(); // Extract response information if available - if ((expcube != "NONE") && (psfcube != "NONE") && + if ((expcube != "NONE") && (psfcube != "NONE") && (bkgcube != "NONE") && (gammalib::strip_whitespace(expcube) != "") && - (gammalib::strip_whitespace(psfcube) != "")) { + (gammalib::strip_whitespace(psfcube) != "") && + (gammalib::strip_whitespace(bkgcube) != "")) { // Get exposure and PSF cubes - GCTACubeExposure exposure(expcube); - GCTACubePsf psf(psfcube); + GCTACubeExposure exposure(expcube); + GCTACubePsf psf(psfcube); + GCTACubeBackground background(bkgcube); // Set reponse - obs->response(exposure, psf); + obs->response(exposure, psf, background); // Signal response availability has_response = true; diff --git a/test/data/model_background_cube_const.xml b/test/data/model_background_cube_const.xml deleted file mode 100644 index b337a8f..0000000 --- a/test/data/model_background_cube_const.xml +++ /dev/null @@ -1,11 +0,0 @@ - - - - - - - - - - - diff --git a/test/data/model_background_cube.xml b/test/data/model_cube_background.xml index 6491db6..2cea6fe 100644 --- a/test/data/model_background_cube.xml +++ b/test/data/model_cube_background.xml @@ -6,8 +6,5 @@ - - - diff --git a/test/data/test_cube.fits b/test/data/test_cube.fits deleted file mode 100644 index 0d7988d..0000000 --- a/test/data/test_cube.fits +++ /dev/null diff --git a/test/test_ctools.sh b/test/test_ctools.sh index 6303fa4..d3ef41c 100755 --- a/test/test_ctools.sh +++ b/test/test_ctools.sh @@ -91,27 +91,6 @@ else echo " events.fits file is not found" exit 1 fi -# -# Run 2 -$ctobssim inmodel="data/model_background_cube.xml" \ - outevents="events2.fits" \ - caldb="irf" \ - irf="cta_dummy_irf" \ - ra=84.17263 \ - dec=22.01444 \ - rad=1.0 \ - tmin=0.0 \ - tmax=1.0 \ - emin=0.3 \ - emax=8.0 chatter=4 -$ECHO -n "." -if [ -s "events2.fits" ] -then - $ECHO -n "." -else - echo " events2.fits file is not found" - exit 1 -fi $ECHO " ok" @@ -148,7 +127,7 @@ $ECHO -n "Test ctbin: " # # Run 1 $ctbin inobs="events.fits" \ - outcube="cntmap.fits" \ + outcube="cntmap1.fits" \ emin=0.1 \ emax=100.0 \ enumbins=20 \ @@ -161,27 +140,27 @@ $ctbin inobs="events.fits" \ yref=22.01 \ proj="CAR" $ECHO -n "." -if [ -s "cntmap.fits" ] +if [ -s "cntmap1.fits" ] then $ECHO -n "." else - $ECHO " cntmap.fits file is not found" + $ECHO " cntmap1.fits file is not found" exit 1 fi # # Run 2 -$ctbin inobs="events2.fits" \ +$ctbin inobs="events.fits" \ outcube="cntmap2.fits" \ emin=0.3 \ emax=8.0 \ - enumbins=10 \ + enumbins=7 \ ebinalg="LOG" \ - nxpix=50 \ - nypix=50 \ - binsz=0.02 \ + nxpix=30 \ + nypix=30 \ + binsz=0.05 \ coordsys="CEL" \ - xref=84.17263 \ - yref=22.01444 \ + xref=86.63 \ + yref=22.01 \ proj="CAR" $ECHO -n "." if [ -s "cntmap2.fits" ] @@ -193,6 +172,168 @@ else fi $ECHO " ok" +# +# Test ctexpcube +# ============== +$ECHO -n "Test ctexpcube: " +# +# Run 1 +$ctexpcube inobs="data/crab_events.fits" \ + incube="NONE" \ + outcube="expcube1.fits" \ + caldb="irf" \ + irf="cta_dummy_irf" \ + emin=0.1 \ + emax=100.0 \ + enumbins=20 \ + nxpix=200 \ + nypix=200 \ + binsz=0.02 \ + coordsys="CEL" \ + xref=83.63 \ + yref=22.01 \ + proj="CAR" +$ECHO -n "." +if [ -s "expcube1.fits" ] +then +$ECHO -n "." +else +$ECHO " expcube1.fits file is not found" +exit 1 +fi +# +# Run 2 +$ctexpcube inobs="events.fits" \ + incube="cntmap2.fits" \ + outcube="expcube2.fits" \ + caldb="irf" \ + irf="cta_dummy_irf" +$ECHO -n "." +if [ -s "expcube2.fits" ] +then +$ECHO -n "." +else +$ECHO " expcube2.fits file is not found" +exit 1 +fi +$ECHO " ok" + + +# +# Test ctpsfcube +# ============== +$ECHO -n "Test ctpsfcube: " +# +# Run 1 +$ctpsfcube inobs="data/crab_events.fits" \ + incube="NONE" \ + outcube="psfcube1.fits" \ + caldb="irf" \ + irf="cta_dummy_irf" \ + emin=0.1 \ + emax=100.0 \ + enumbins=20 \ + nxpix=10 \ + nypix=10 \ + binsz=0.4 \ + coordsys="CEL" \ + xref=83.63 \ + yref=22.01 \ + proj="CAR" \ + amax=0.3 \ + anumbins=10 +$ECHO -n "." +if [ -s "psfcube1.fits" ] +then +$ECHO -n "." +else +$ECHO " psfcube1.fits file is not found" +exit 1 +fi +# +# Run 2 +$ctpsfcube inobs="events.fits" \ + incube="cntmap2.fits" \ + outcube="psfcube2.fits" \ + caldb="irf" \ + irf="cta_dummy_irf" \ + amax=0.3 \ + anumbins=2 +$ECHO -n "." +if [ -s "psfcube2.fits" ] +then +$ECHO -n "." +else +$ECHO " psfcube2.fits file is not found" +exit 1 +fi +$ECHO " ok" + + +# +# Test ctbkgcube +# ============== +$ECHO -n "Test ctbkgcube: " +# +# Run 1 +$ctbkgcube inobs="data/crab_events.fits" \ + inmodel="data/crab.xml" \ + incube="NONE" \ + outcube="bkgcube1.fits" \ + outmodel="bkgcube1.xml" \ + caldb="irf" \ + irf="cta_dummy_irf" \ + emin=0.1 \ + emax=100.0 \ + enumbins=20 \ + nxpix=10 \ + nypix=10 \ + binsz=0.4 \ + coordsys="CEL" \ + xref=83.63 \ + yref=22.01 \ + proj="CAR" +$ECHO -n "." +if [ -s "bkgcube1.fits" ] +then +$ECHO -n "." +else +$ECHO " bkgcube1.fits file is not found" +exit 1 +fi +if [ -s "bkgcube1.xml" ] +then +$ECHO -n "." +else +$ECHO " bkgcube1.xml file is not found" +exit 1 +fi +# +# Run 2 +$ctbkgcube inobs="events.fits" \ + inmodel="data/crab.xml" \ + incube="cntmap2.fits" \ + caldb="irf" \ + irf="cta_dummy_irf" \ + outcube="bkgcube2.fits" \ + outmodel="bkgcube2.xml" +$ECHO -n "." +if [ -s "bkgcube2.fits" ] +then +$ECHO -n "." +else +$ECHO " bkgcube2.fits file is not found" +exit 1 +fi +if [ -s "bkgcube2.fits" ] +then +$ECHO -n "." +else +$ECHO " bkgcube2.xml file is not found" +exit 1 +fi +$ECHO " ok" + # @@ -201,11 +342,12 @@ $ECHO " ok" $ECHO -n "Test ctmodel: " # # Run 1 -$ctmodel inobs="cntmap.fits" \ - incube="cntmap.fits" \ +$ctmodel inobs="cntmap1.fits" \ + incube="cntmap1.fits" \ outcube="modmap1.fits" \ expcube="NONE" \ psfcube="NONE" \ + bkgcube="NONE" \ caldb="irf" \ irf="cta_dummy_irf" \ inmodel="data/crab.xml" @@ -217,13 +359,35 @@ else $ECHO " modmap1.fits file is not found" exit 1 fi + # # Run 2 +$ctmodel inobs="cntmap2.fits" \ + incube="cntmap2.fits" \ + outcube="modmap2.fits" \ + expcube="expcube2.fits" \ + psfcube="psfcube2.fits" \ + bkgcube="bkgcube2.fits" \ + caldb="irf" \ + irf="cta_dummy_irf" \ + inmodel="bkgcube2.xml" +$ECHO -n "." +if [ -s "modmap2.fits" ] +then +$ECHO -n "." +else +$ECHO " modmap2.fits file is not found" +exit 1 +fi + +# +# Run 3 $ctmodel inobs="NONE" \ incube="NONE" \ outcube="modmap3.fits" \ expcube="NONE" \ psfcube="NONE" \ + bkgcube="NONE" \ caldb="irf" \ irf="cta_dummy_irf" \ inmodel="data/crab.xml" \ @@ -250,24 +414,6 @@ else $ECHO " modmap3.fits file is not found" exit 1 fi -# -# Run 3 -$ctmodel inobs="cntmap2.fits" \ - incube="cntmap2.fits" \ - outcube="modmap2.fits" \ - expcube="NONE" \ - psfcube="NONE" \ - caldb="irf" \ - irf="cta_dummy_irf" \ - inmodel="data/model_background_cube.xml" -$ECHO -n "." -if [ -s "modmap2.fits" ] -then - $ECHO -n "." -else - $ECHO " modmap2.fits file is not found" - exit 1 -fi $ECHO " ok" @@ -301,11 +447,12 @@ $ECHO " ok" $ECHO -n "Test ctlike: " # # Run 1 -$ctlike inobs="cntmap.fits" \ +$ctlike inobs="cntmap1.fits" \ inmodel="data/crab.xml" \ outmodel="results_binned.xml" \ expcube="NONE" \ psfcube="NONE" \ + bkgcube="NONE" \ caldb="irf" \ irf="cta_dummy_irf" $ECHO -n "." @@ -331,36 +478,23 @@ else $ECHO " results_unbinned.xml file is not found" exit 1 fi + # # Run 3 -$ctlike inobs="events2.fits" \ - inmodel="data/model_background_cube.xml" \ - outmodel="results_unbinned_background_cube.xml" \ - caldb="irf" \ - irf="cta_dummy_irf" -$ECHO -n "." -if [ -s "results_unbinned_background_cube.xml" ] -then - $ECHO -n "." -else - $ECHO " results_unbinned_background_cube.xml file is not found" - exit 1 -fi -# -# Run 4 $ctlike inobs="cntmap2.fits" \ - inmodel="data/model_background_cube.xml" \ - outmodel="results_binned_background_cube.xml" \ - expcube="NONE" \ - psfcube="NONE" \ + inmodel="bkgcube2.xml" \ + outmodel="results_binned_cube_background.xml" \ + expcube="expcube2.fits" \ + psfcube="psfcube2.fits" \ + bkgcube="bkgcube2.fits" \ caldb="irf" \ irf="cta_dummy_irf" $ECHO -n "." -if [ -s "results_binned_background_cube.xml" ] +if [ -s "results_binned_cube_background.xml" ] then $ECHO -n "." else - $ECHO " results_binned_background_cube.xml file is not found" + $ECHO " results_binned_cube_background.xml file is not found" exit 1 fi $ECHO " ok" @@ -394,167 +528,6 @@ fi $ECHO " ok" -# -# Test ctexpcube -# ============== -$ECHO -n "Test ctexpcube: " -# -# Run 1 -$ctexpcube inobs="data/crab_events.fits" \ - incube="NONE" \ - outcube="expcube1.fits" \ - caldb="irf" \ - irf="cta_dummy_irf" \ - emin=0.1 \ - emax=100.0 \ - enumbins=20 \ - nxpix=200 \ - nypix=200 \ - binsz=0.02 \ - coordsys="CEL" \ - xref=83.63 \ - yref=22.01 \ - proj="CAR" -$ECHO -n "." -if [ -s "expcube1.fits" ] -then - $ECHO -n "." -else - $ECHO " expcube1.fits file is not found" - exit 1 -fi -# -# Run 2 -$ctexpcube inobs="data/crab_events.fits" \ - incube="data/crab_cntmap.fits" \ - outcube="expcube2.fits" \ - caldb="irf" \ - irf="cta_dummy_irf" -$ECHO -n "." -if [ -s "expcube2.fits" ] -then - $ECHO -n "." -else - $ECHO " expcube2.fits file is not found" - exit 1 -fi -$ECHO " ok" - - -# -# Test ctpsfcube -# ============== -$ECHO -n "Test ctpsfcube: " -# -# Run 1 -$ctpsfcube inobs="data/crab_events.fits" \ - incube="NONE" \ - outcube="psfcube1.fits" \ - caldb="irf" \ - irf="cta_dummy_irf" \ - emin=0.1 \ - emax=100.0 \ - enumbins=20 \ - nxpix=10 \ - nypix=10 \ - binsz=0.4 \ - coordsys="CEL" \ - xref=83.63 \ - yref=22.01 \ - proj="CAR" \ - amax=0.3 \ - anumbins=10 -$ECHO -n "." -if [ -s "psfcube1.fits" ] -then - $ECHO -n "." -else - $ECHO " psfcube1.fits file is not found" - exit 1 -fi -# -# Run 2 -$ctpsfcube inobs="data/crab_events.fits" \ - incube="data/crab_cntmap.fits" \ - outcube="psfcube2.fits" \ - caldb="irf" \ - irf="cta_dummy_irf" \ - amax=0.3 \ - anumbins=2 -$ECHO -n "." -if [ -s "psfcube2.fits" ] -then - $ECHO -n "." -else - $ECHO " psfcube2.fits file is not found" - exit 1 -fi -$ECHO " ok" - - -# -# Test ctbkgcube -# ============== -$ECHO -n "Test ctbkgcube: " -# -# Run 1 -$ctbkgcube inobs="data/crab_events.fits" \ - inmodel="data/crab.xml" \ - incube="NONE" \ - outcube="bkgcube1.fits" \ - outmodel="bkgcube1.xml" \ - caldb="irf" \ - irf="cta_dummy_irf" \ - emin=0.1 \ - emax=100.0 \ - enumbins=20 \ - nxpix=10 \ - nypix=10 \ - binsz=0.4 \ - coordsys="CEL" \ - xref=83.63 \ - yref=22.01 \ - proj="CAR" -$ECHO -n "." -if [ -s "bkgcube1.fits" ] -then - $ECHO -n "." -else - $ECHO " bkgcube1.fits file is not found" - exit 1 -fi -if [ -s "bkgcube1.xml" ] -then - $ECHO -n "." -else - $ECHO " bkgcube1.xml file is not found" - exit 1 -fi -# -# Run 2 -$ctbkgcube inobs="data/crab_events.fits" \ - inmodel="data/crab.xml" \ - incube="data/crab_cntmap.fits" \ - caldb="irf" \ - irf="cta_dummy_irf" \ - outcube="bkgcube2.fits" \ - outmodel="bkgcube2.xml" -$ECHO -n "." -if [ -s "bkgcube2.fits" ] -then - $ECHO -n "." -else - $ECHO " bkgcube2.fits file is not found" - exit 1 -fi -if [ -s "bkgcube2.fits" ] -then - $ECHO -n "." -else - $ECHO " bkgcube2.xml file is not found" - exit 1 -fi -$ECHO " ok" # @@ -608,10 +581,10 @@ $ECHO " ok" # =============== $ECHO -n "Test ctulimit: " $ctulimit inobs="data/crab_events.fits.gz" \ - inmodel="data/crab.xml" \ - srcname="Crab" \ - caldb="irf" \ - irf="cta_dummy_irf" + inmodel="data/crab.xml" \ + srcname="Crab" \ + caldb="irf" \ + irf="cta_dummy_irf" $ECHO -n "." if [ -s "ctulimit.log" ] then