37 #define G_APPLY_MASK "ctcubemask::apply_mask(GCTAObservation*)"
187 this->GApplication::clear();
217 log_header1(TERSE,
"Apply mask");
220 int n_observations = 0;
223 for (
int i = 0; i <
m_obs.size(); ++i) {
229 GCTAObservation*
obs =
dynamic_cast<GCTAObservation*
>(
m_obs[i]);
252 if (n_observations > 1) {
291 log_header1(TERSE, gammalib::number(
"Save counts cube",
m_obs.size()));
294 m_outcube = (*this)[
"outcube"].filename();
300 m_prefix = (*this)[
"prefix"].string();
327 log_header1(TERSE,
"Publish counts cube");
330 std::string user_name(name);
331 if (user_name.empty()) {
336 GCTAObservation*
obs =
dynamic_cast<GCTAObservation*
>(
m_obs[0]);
340 GCTAEventCube* cube =
dynamic_cast<GCTAEventCube*
>(obs->events());
344 log_value(NORMAL,
"Counts cube name", user_name);
347 cube->counts().publish(user_name);
437 (*this)[
"regfile"].query();
443 if ((*
this)[
"emin"].is_valid() && (*
this)[
"emax"].is_valid()) {
444 m_emin = (*this)[
"emin"].real();
445 m_emax = (*this)[
"emax"].real();
453 m_publish = (*this)[
"publish"].boolean();
457 (*this)[
"outcube"].query();
458 (*this)[
"prefix"].query();
481 GCTAEventCube* cube =
dynamic_cast<GCTAEventCube*
>
482 (
const_cast<GEvents*
>(obs->events()));
486 GSkyMap map = cube->counts();
487 const GEbounds& ebounds = cube->ebounds();
491 m_emin = ebounds.emin().TeV();
492 m_emax = ebounds.emax().TeV();
496 int npix = map.npix();
497 int n_ebin = ebounds.size();
502 double sum_before = 0.0;
503 for (
int i = 0; i < n_ebin; ++i) {
504 for (
int pixel = 0; pixel < npix; ++pixel) {
505 if (map(pixel, i) >= 0.0) {
506 sum_before += map(pixel, i);
512 for (
int i = 0; i < n_ebin; ++i) {
513 double emin = ebounds.emin(i).TeV() + 1.0e-6;
521 for (
int i = n_ebin-1; i >= 0; i--) {
522 double emax = ebounds.emax(i).TeV() - 1.0e-6;
530 for (
int i = 0; i < e_idx1; ++i) {
531 for (
int pixel = 0; pixel < npix; ++pixel) {
535 for (
int i = e_idx2 + 1; i < n_ebin; ++i) {
536 for (
int pixel = 0; pixel < npix; ++pixel) {
542 log_value(NORMAL,
"Selected energy band",
543 gammalib::str(ebounds.emin(e_idx1).TeV())+
" - "+
544 gammalib::str(ebounds.emax(e_idx2).TeV())+
" TeV");
548 GCTARoi roi =
get_roi(obs->pointing());
549 if (roi.is_valid()) {
550 GSkyRegionCircle circle(roi.centre().dir(), roi.radius());
551 for (
int i = e_idx1; i <= e_idx2; ++i) {
552 for (
int pixel = 0; pixel < npix; ++pixel) {
553 GSkyDir dir = map.inx2dir(pixel);
554 if (!circle.contains(dir)) {
561 log_value(NORMAL,
"Selected RoI",
562 "RA="+gammalib::str(circle.centre().ra_deg())+
" deg, "+
563 "DEC="+gammalib::str(circle.centre().dec_deg())+
" deg, "+
564 "Radius="+gammalib::str(circle.radius())+
" deg");
570 if ((*
this)[
"regfile"].is_valid()) {
571 GSkyRegions regions = GSkyRegions((*
this)[
"regfile"].filename());
572 for (
int i = e_idx1; i <= e_idx2; ++i) {
573 for (
int pixel = 0; pixel < npix; ++pixel) {
574 GSkyDir dir = map.inx2dir(pixel);
575 if (regions.contains(dir)) {
582 for (
int i = 0; i < regions.size(); ++i) {
583 log_value(NORMAL,
"Exclusion region "+gammalib::str(i+1),
588 log_value(NORMAL,
"Exclusion regions",
"None");
592 double sum_after = 0.0;
593 for (
int i = 0; i < n_ebin; ++i) {
594 for (
int pixel = 0; pixel < npix; ++pixel) {
595 if (map(pixel, i) >= 0.0) {
596 sum_after += map(pixel, i);
602 log_value(NORMAL,
"Events before masking", sum_before);
603 log_value(NORMAL,
"Events after masking", sum_after);
628 const GSkyRegionCircle* circle =
dynamic_cast<const GSkyRegionCircle*
>(®ion);
632 if (circle != NULL) {
633 rstring.append(
"RA="+gammalib::str(circle->centre().ra_deg())+
" deg, "+
634 "DEC="+gammalib::str(circle->centre().dec_deg())+
" deg, "+
635 "Radius="+gammalib::str(circle->radius())+
" deg");
655 std::vector<std::string> elements = gammalib::split(filename,
"/");
658 std::string outname =
m_prefix + elements[elements.size()-1];
674 if (
m_obs.size() > 0) {
677 GCTAObservation*
obs =
dynamic_cast<GCTAObservation*
>(
m_obs[0]);
683 log_value(NORMAL,
"Counts cube file",
m_outcube.url());
717 for (
int i = 0; i <
m_obs.size(); ++i) {
720 GCTAObservation*
obs =
dynamic_cast<GCTAObservation*
>(
m_obs[i]);
729 obs->eventfile(outfile);
732 log_value(NORMAL,
"Counts cube file", outfile);
735 obs->save(outfile, clobber());
void apply_mask(GCTAObservation *obs)
Apply mask to event cube.
void init_members(void)
Initialise class members.
bool m_select_energy
Perform energy selection.
GFilename m_outcube
Output event list or XML file.
const GObservations & obs(void) const
Return observation container.
ctcubemask & operator=(const ctcubemask &app)
Assignment operator.
void save_xml(void)
Save counts cube(s) in XML format.
void publish(const std::string &name="")
Publish counts cube.
std::string m_prefix
Prefix for multiple counts maps.
void save(void)
Save the masked event cube(s)
std::string region_string(const GSkyRegion ®ion) const
Return region string.
void copy_members(const ctcubemask &app)
Copy class members.
void get_parameters(void)
Get application parameters.
Base class for observation tools.
ctcubemask(void)
Void constructor.
void free_members(void)
Delete class members.
Cube filter tool definition.
void init_members(void)
Initialise class members.
virtual ~ctcubemask(void)
Destructor.
double m_emin
Lower energy.
ctobservation & operator=(const ctobservation &app)
Assignment operator.
std::vector< std::string > m_infiles
Input event filenames.
void process(void)
Mask data cube.
void clear(void)
Clear ctcubemask tool.
std::string set_outfile_name(const std::string &filename) const
Set output file name.
GObservations m_obs
Observation container.
bool m_publish
Publish counts cube?
void save_fits(void)
Save counts cube in FITS format.
void free_members(void)
Delete class members.
double m_emax
Upper energy.