Fabio Iocco / ctools

Commit 276e7423b0bf89c1cdca7f3bbf46f57af1216cc1

Authored by Jürgen Knödlseder
1 parent a71a6a96

ctbkgcube now produces a model container that replaces the original background m… 

…odel with a standard diffuse cube background model that can be used for stacked analysis. The model container can be either retrieved using the models() method, or it can be written to disk with the save() method.
Inline Side-by-side
Showing 6 changed files with 147 additions and 49 deletions
pyext/ctbkgcube.i
  View file @276e742
... ... @@ -48,6 +48,7 @@ public:
48 48 void run(void);
49 49 void save(void);
50 50 const GCTAEventCube& cube(void) const;
  51 + const GModels& models(void) const;
51 52 };
52 53  
53 54 /***********************************************************************//**
... ...
src/ctbkgcube/ctbkgcube.cpp
  View file @276e742
... ... @@ -42,6 +42,8 @@
42 42  
43 43 /* __ Coding definitions _________________________________________________ */
44 44  
  45 +/* __ Constants __________________________________________________________ */
  46 +
45 47  
46 48 /*==========================================================================
47 49 = =
... ... @@ -233,16 +235,27 @@ void ctbkgcube::run(void)
233 235 log << m_obs << std::endl;
234 236 }
235 237  
  238 + // Log input model
  239 + if (logTerse()) {
  240 + log << std::endl;
  241 + log.header1("Input model");
  242 + log << m_obs.models() << std::endl;
  243 + }
  244 +
236 245 // Write header
237 246 if (logTerse()) {
238 247 log << std::endl;
239 248 log.header1("Prepare model");
240 249 }
241 250  
242   - // Copy models from observation container
  251 + // Copy models from observation container and reset output model
  252 + // container
243 253 m_bkgmdl = m_obs.models();
  254 + m_outmdl.clear();
  255 + m_cube_model = -1;
244 256  
245   - // Remove all model that are not CTA background models from the
  257 + // Remove all models that are not CTA background models from the
  258 + // container and put all removed components in the output
246 259 // container
247 260 int num = m_bkgmdl.size();
248 261 for (int i = num-1; i >= 0; --i) {
... ... @@ -281,8 +294,10 @@ void ctbkgcube::run(void)
281 294 log << std::endl;
282 295 }
283 296  
284   - // Remove model if requested
  297 + // If removal is requested, append model to output container and
  298 + // remove it from the background model container ...
285 299 if (remove) {
  300 + m_outmdl.append(*(m_bkgmdl[i]));
286 301 m_bkgmdl.remove(i);
287 302 }
288 303  
... ... @@ -345,6 +360,27 @@ void ctbkgcube::run(void)
345 360 log << n_observations << std::endl;
346 361 }
347 362  
  363 + // Create a background model for output background cube and append
  364 + // that model to the input model in place of the original
  365 + // background model
  366 + GEnergies energies;
  367 + for (int i = 0; i < m_bkgcube.ebins(); ++i) {
  368 + energies.append(m_bkgcube.energy(i));
  369 + }
  370 + GModelSpatialDiffuseCube spatial(m_bkgcube.map(), energies, 1.0);
  371 + GModelSpectralPlaw spectral(1.0, 0.0, GEnergy(1.0, "TeV"));
  372 + GCTAModelCubeBackground model(spatial, spectral);
  373 + model.name("ctbkgcube default background model");
  374 + m_cube_model = m_outmdl.size(); // Store the slot number for save()
  375 + m_outmdl.append(model);
  376 +
  377 + // Log output model
  378 + if (logTerse()) {
  379 + log << std::endl;
  380 + log.header1("Output model");
  381 + log << m_outmdl << std::endl;
  382 + }
  383 +
348 384 // Return
349 385 return;
350 386 }
... ... @@ -369,8 +405,9 @@ void ctbkgcube::save(void)
369 405 energies.append(m_ebounds.elogmean(i));
370 406 }
371 407  
372   - // Get output filename
373   - std::string filename = (*this)["outfile"].filename();
  408 + // Get output filenames
  409 + std::string outfile = (*this)["outfile"].filename();
  410 + std::string outmodel = (*this)["outmodel"].filename();
374 411  
375 412 // Create empty FITS file
376 413 GFits fits;
... ... @@ -382,7 +419,28 @@ void ctbkgcube::save(void)
382 419 energies.write(fits);
383 420  
384 421 // Save FITS file
385   - fits.saveto(filename, clobber());
  422 + fits.saveto(outfile, clobber());
  423 +
  424 + // Write output models if filename is valid
  425 + if ((outmodel.length() > 0) && (gammalib::tolower(outmodel) != "none")) {
  426 +
  427 + // Set filename of map cube.
  428 + if (m_cube_model >= 0) {
  429 + GCTAModelCubeBackground* model =
  430 + dynamic_cast<GCTAModelCubeBackground*>(m_outmdl[m_cube_model]);
  431 + if (model != NULL) {
  432 + GModelSpatialDiffuseCube* spatial =
  433 + dynamic_cast<GModelSpatialDiffuseCube*>(model->spatial());
  434 + if (spatial != NULL) {
  435 + spatial->filename(outfile);
  436 + }
  437 + }
  438 + }
  439 +
  440 + // Save output model container
  441 + m_outmdl.save(outmodel);
  442 +
  443 + }
386 444  
387 445 // Return
388 446 return;
... ... @@ -404,8 +462,10 @@ void ctbkgcube::init_members(void)
404 462 m_obs.clear();
405 463 m_bkgcube.clear();
406 464 m_bkgmdl.clear();
  465 + m_outmdl.clear();
407 466 m_ebounds.clear();
408   - m_livetime = 0.0;
  467 + m_livetime = 0.0;
  468 + m_cube_model = -1;
409 469  
410 470 // Return
411 471 return;
... ... @@ -420,11 +480,13 @@ void ctbkgcube::init_members(void)
420 480 void ctbkgcube::copy_members(const ctbkgcube& app)
421 481 {
422 482 // Copy members
423   - m_obs = app.m_obs;
424   - m_bkgcube = app.m_bkgcube;
425   - m_bkgmdl = app.m_bkgmdl;
426   - m_ebounds = app.m_ebounds;
427   - m_livetime = app.m_livetime;
  483 + m_obs = app.m_obs;
  484 + m_bkgcube = app.m_bkgcube;
  485 + m_bkgmdl = app.m_bkgmdl;
  486 + m_outmdl = app.m_outmdl;
  487 + m_ebounds = app.m_ebounds;
  488 + m_livetime = app.m_livetime;
  489 + m_cube_model = app.m_cube_model;
428 490  
429 491 // Return
430 492 return;
... ... @@ -461,15 +523,15 @@ void ctbkgcube::get_parameters(void)
461 523 // If there are no models associated with the observations then load now
462 524 // the model definition from the XML file
463 525 if (m_obs.models().size() == 0) {
464   - if ((*this)["bkgmdl"].is_undefined()) {
465   - std::string msg = "No background model definition XML file "
466   - "specified. Please set the \"bkgmdl\" parameter "
467   - "to the background model definition XML "
468   - "filename.";
  526 + if ((*this)["inmodel"].is_undefined()) {
  527 + std::string msg = "No model definition XML file specified. "
  528 + "Please set the \"inmodel\" parameter to the "
  529 + "XML file that contains the background model "
  530 + "definition.";
469 531 throw GException::invalid_value(G_GET_PARAMETERS, msg);
470 532 }
471   - std::string bkgmdl = (*this)["bkgmdl"].filename();
472   - GModels models(bkgmdl);
  533 + std::string inmodel = (*this)["inmodel"].filename();
  534 + GModels models(inmodel);
473 535 m_obs.models(models);
474 536 }
475 537  
... ... @@ -511,9 +573,10 @@ void ctbkgcube::get_parameters(void)
511 573  
512 574 }
513 575  
514   - // Read output filename (if needed)
  576 + // Read output filenames (if needed)
515 577 if (read_ahead()) {
516   - std::string filename = (*this)["outfile"].filename();
  578 + (*this)["outfile"].filename();
  579 + (*this)["outmodel"].filename();
517 580 }
518 581  
519 582 // Return
... ...
src/ctbkgcube/ctbkgcube.hpp
  View file @276e742
... ... @@ -60,6 +60,7 @@ public:
60 60 void run(void);
61 61 void save(void);
62 62 const GCTAEventCube& cube(void) const;
  63 + const GModels& models(void) const;
63 64  
64 65 protected:
65 66 // Protected methods
... ... @@ -75,8 +76,10 @@ protected:
75 76 GObservations m_obs; //!< Observation container
76 77 GCTAEventCube m_bkgcube; //!< Background cube
77 78 GModels m_bkgmdl; //!< CTA background models
  79 + GModels m_outmdl; //!< Output models
78 80 GEbounds m_ebounds; //!< Energy boundaries
79 81 double m_livetime; //!< Total livetime for normalization
  82 + int m_cube_model; //!< Cube model index
80 83 };
81 84  
82 85  
... ... @@ -92,4 +95,17 @@ const GCTAEventCube&amp; ctbkgcube::cube(void) const
92 95 return (m_bkgcube);
93 96 }
94 97  
  98 +
  99 +/***********************************************************************//**
  100 + * @brief Return background model container
  101 + *
  102 + * @return Background model container.
  103 + ***************************************************************************/
  104 +inline
  105 +const GModels& ctbkgcube::models(void) const
  106 +{
  107 + // Return background model container
  108 + return (m_outmdl);
  109 +}
  110 +
95 111 #endif /* CTBKGCUBE_HPP */
... ...
src/ctbkgcube/ctbkgcube.par
  View file @276e742
... ... @@ -18,9 +18,10 @@
18 18 # General parameters
19 19 #===================
20 20 infile, f, a, "NONE",,, "Event list or observation definition file"
21   -bkgmdl, f, a, "NONE",,, "Background model XML file"
  21 +inmodel, f, a, "NONE",,, "Input (background) model XML file"
22 22 cntmap, f, a, "NONE",,, "Counts cube for background cube definition"
23 23 outfile, f, a, "bkgcube.fits",,, "Output background cube file"
  24 +outmodel, f, a, "NONE",,,"Output (background) model XML file"
24 25 #
25 26 # Energy binning parameters
26 27 #==========================
... ...
test/test_ctbkgcube.py
  View file @276e742
... ... @@ -66,9 +66,10 @@ class Test(gammalib.GPythonTestSuite):
66 66 # Set-up ctbkgcube
67 67 bkgcube = ctools.ctbkgcube()
68 68 bkgcube["infile"].filename(self.events_name)
69   - bkgcube["bkgmdl"].filename(self.bkg_model)
  69 + bkgcube["inmodel"].filename(self.bkg_model)
70 70 bkgcube["cntmap"].filename("NONE")
71   - bkgcube["outfile"].filename("psfcube.fits")
  71 + bkgcube["outfile"].filename("bkgcube.fits")
  72 + bkgcube["outmodel"].filename("bkgcube.xml")
72 73 bkgcube["ebinalg"].string("LOG")
73 74 bkgcube["emin"].real(0.1)
74 75 bkgcube["emax"].real(100.0)
... ...
test/test_ctools.sh
  View file @276e742
... ... @@ -85,7 +85,7 @@ if [ -s &quot;events.fits&quot; ]
85 85 then
86 86 $ECHO -n "."
87 87 else
88   - echo " events.fits file is not valid"
  88 + echo " events.fits file is not found"
89 89 exit 1
90 90 fi
91 91 #
... ... @@ -106,7 +106,7 @@ if [ -s &quot;events2.fits&quot; ]
106 106 then
107 107 $ECHO -n "."
108 108 else
109   - echo " events2.fits file is not valid"
  109 + echo " events2.fits file is not found"
110 110 exit 1
111 111 fi
112 112 $ECHO " ok"
... ... @@ -132,7 +132,7 @@ if [ -s &quot;skymap.fits&quot; ]
132 132 then
133 133 $ECHO -n "."
134 134 else
135   - $ECHO " skymap.fits file is not valid"
  135 + $ECHO " skymap.fits file is not found"
136 136 exit 1
137 137 fi
138 138 $ECHO " ok"
... ... @@ -162,7 +162,7 @@ if [ -s &quot;cntmap.fits&quot; ]
162 162 then
163 163 $ECHO -n "."
164 164 else
165   - $ECHO " cntmap.fits file is not valid"
  165 + $ECHO " cntmap.fits file is not found"
166 166 exit 1
167 167 fi
168 168 #
... ... @@ -185,7 +185,7 @@ if [ -s &quot;cntmap2.fits&quot; ]
185 185 then
186 186 $ECHO -n "."
187 187 else
188   - $ECHO " cntmap2.fits file is not valid"
  188 + $ECHO " cntmap2.fits file is not found"
189 189 exit 1
190 190 fi
191 191 $ECHO " ok"
... ... @@ -208,7 +208,7 @@ if [ -s &quot;modmap1.fits&quot; ]
208 208 then
209 209 $ECHO -n "."
210 210 else
211   - $ECHO " modmap1.fits file is not valid"
  211 + $ECHO " modmap1.fits file is not found"
212 212 exit 1
213 213 fi
214 214 #
... ... @@ -238,7 +238,7 @@ if [ -s &quot;modmap3.fits&quot; ]
238 238 then
239 239 $ECHO -n "."
240 240 else
241   - $ECHO " modmap3.fits file is not valid"
  241 + $ECHO " modmap3.fits file is not found"
242 242 exit 1
243 243 fi
244 244 #
... ... @@ -254,7 +254,7 @@ if [ -s &quot;modmap2.fits&quot; ]
254 254 then
255 255 $ECHO -n "."
256 256 else
257   - $ECHO " modmap2.fits file is not valid"
  257 + $ECHO " modmap2.fits file is not found"
258 258 exit 1
259 259 fi
260 260 $ECHO " ok"
... ... @@ -278,7 +278,7 @@ if [ -s &quot;selected_events.fits&quot; ]
278 278 then
279 279 $ECHO -n "."
280 280 else
281   - $ECHO " selected_events.fits file is not valid"
  281 + $ECHO " selected_events.fits file is not found"
282 282 exit 1
283 283 fi
284 284 $ECHO " ok"
... ... @@ -300,7 +300,7 @@ if [ -s &quot;results_binned.xml&quot; ]
300 300 then
301 301 $ECHO -n "."
302 302 else
303   - $ECHO " results_binned.xml file is not valid"
  303 + $ECHO " results_binned.xml file is not found"
304 304 exit 1
305 305 fi
306 306 #
... ... @@ -315,7 +315,7 @@ if [ -s &quot;results_unbinned.xml&quot; ]
315 315 then
316 316 $ECHO -n "."
317 317 else
318   - $ECHO " results_unbinned.xml file is not valid"
  318 + $ECHO " results_unbinned.xml file is not found"
319 319 exit 1
320 320 fi
321 321 #
... ... @@ -330,7 +330,7 @@ if [ -s &quot;results_unbinned_background_cube.xml&quot; ]
330 330 then
331 331 $ECHO -n "."
332 332 else
333   - $ECHO " results_unbinned_background_cube.xml file is not valid"
  333 + $ECHO " results_unbinned_background_cube.xml file is not found"
334 334 exit 1
335 335 fi
336 336 #
... ... @@ -345,7 +345,7 @@ if [ -s &quot;results_binned_background_cube.xml&quot; ]
345 345 then
346 346 $ECHO -n "."
347 347 else
348   - $ECHO " results_binned_background_cube.xml file is not valid"
  348 + $ECHO " results_binned_background_cube.xml file is not found"
349 349 exit 1
350 350 fi
351 351 $ECHO " ok"
... ... @@ -373,7 +373,7 @@ if [ -s &quot;tsmap.fits&quot; ]
373 373 then
374 374 $ECHO -n "."
375 375 else
376   - $ECHO " tsmap.fits file is not valid"
  376 + $ECHO " tsmap.fits file is not found"
377 377 exit 1
378 378 fi
379 379 $ECHO " ok"
... ... @@ -405,7 +405,7 @@ if [ -s &quot;expcube1.fits&quot; ]
405 405 then
406 406 $ECHO -n "."
407 407 else
408   - $ECHO " expcube1.fits file is not valid"
  408 + $ECHO " expcube1.fits file is not found"
409 409 exit 1
410 410 fi
411 411 #
... ... @@ -420,7 +420,7 @@ if [ -s &quot;expcube2.fits&quot; ]
420 420 then
421 421 $ECHO -n "."
422 422 else
423   - $ECHO " expcube2.fits file is not valid"
  423 + $ECHO " expcube2.fits file is not found"
424 424 exit 1
425 425 fi
426 426 $ECHO " ok"
... ... @@ -454,7 +454,7 @@ if [ -s &quot;psfcube1.fits&quot; ]
454 454 then
455 455 $ECHO -n "."
456 456 else
457   - $ECHO " psfcube1.fits file is not valid"
  457 + $ECHO " psfcube1.fits file is not found"
458 458 exit 1
459 459 fi
460 460 #
... ... @@ -471,7 +471,7 @@ if [ -s &quot;psfcube2.fits&quot; ]
471 471 then
472 472 $ECHO -n "."
473 473 else
474   - $ECHO " psfcube2.fits file is not valid"
  474 + $ECHO " psfcube2.fits file is not found"
475 475 exit 1
476 476 fi
477 477 $ECHO " ok"
... ... @@ -484,9 +484,10 @@ $ECHO -n &quot;Test ctbkgcube: &quot;
484 484 #
485 485 # Run 1
486 486 $ctbkgcube infile="data/crab_events.fits" \
487   - bkgmdl="data/crab.xml" \
  487 + inmodel="data/crab.xml" \
488 488 cntmap="NONE" \
489 489 outfile="bkgcube1.fits" \
  490 + outmodel="bkgcube1.xml" \
490 491 emin=0.1 \
491 492 emax=100.0 \
492 493 enumbins=20 \
... ... @@ -502,21 +503,36 @@ if [ -s &quot;bkgcube1.fits&quot; ]
502 503 then
503 504 $ECHO -n "."
504 505 else
505   - $ECHO " bkgcube1.fits file is not valid"
  506 + $ECHO " bkgcube1.fits file is not found"
  507 + exit 1
  508 +fi
  509 +if [ -s "bkgcube1.xml" ]
  510 +then
  511 + $ECHO -n "."
  512 +else
  513 + $ECHO " bkgcube1.xml file is not found"
506 514 exit 1
507 515 fi
508 516 #
509 517 # Run 2
510 518 $ctbkgcube infile="data/crab_events.fits" \
511   - bkgmdl="data/crab.xml" \
  519 + inmodel="data/crab.xml" \
512 520 cntmap="data/crab_cntmap.fits" \
513   - outfile="bkgcube2.fits"
  521 + outfile="bkgcube2.fits" \
  522 + outmodel="bkgcube2.xml"
514 523 $ECHO -n "."
515 524 if [ -s "bkgcube2.fits" ]
516 525 then
517 526 $ECHO -n "."
518 527 else
519   - $ECHO " bkgcube2.fits file is not valid"
  528 + $ECHO " bkgcube2.fits file is not found"
  529 + exit 1
  530 +fi
  531 +if [ -s "bkgcube2.fits" ]
  532 +then
  533 + $ECHO -n "."
  534 +else
  535 + $ECHO " bkgcube2.xml file is not found"
520 536 exit 1
521 537 fi
522 538 $ECHO " ok"
... ... @@ -539,7 +555,7 @@ if [ -s &quot;filtered_cube.fits&quot; ]
539 555 then
540 556 $ECHO -n "."
541 557 else
542   - $ECHO " filtered_cube.fits file is not valid"
  558 + $ECHO " filtered_cube.fits file is not found"
543 559 exit 1
544 560 fi
545 561 $ECHO " ok"
... ... @@ -562,7 +578,7 @@ if [ -s &quot;butterfly.txt&quot; ]
562 578 then
563 579 $ECHO -n "."
564 580 else
565   - $ECHO " butterfly.txt file is not valid"
  581 + $ECHO " butterfly.txt file is not found"
566 582 exit 1
567 583 fi
568 584 $ECHO " ok"
... ...