Johan Bregeon / ctools

Commit b2026ebc69f7753e7c6857527966e2d700f4f54b

Authored by Jürgen Knödlseder
2 parents d7f924ca 412c7728

Merge branch '1566-csspec-update-for-binned-mode' into integration

Inline Side-by-side
Showing 10 changed files with 179 additions and 225 deletions
cscripts/csspec.py
  View file @b2026eb
... ... @@ -43,7 +43,7 @@ class csspec(ctools.cscript):
43 43 self.version = "1.0.0"
44 44  
45 45 # Initialise some members
46   - self.obs = None
  46 + self.obs = gammalib.GObservations()
47 47 self.fits = None
48 48  
49 49 # Initialise some members
... ... @@ -106,7 +106,7 @@ class csspec(ctools.cscript):
106 106 pars.append(gammalib.GApplicationPar("outfile","f","a","spectrum.fits","","","Output file name"))
107 107 pars.append(gammalib.GApplicationPar("expcube","f","a","NONE","","","Exposure cube file (only needed for stacked analysis)"))
108 108 pars.append(gammalib.GApplicationPar("psfcube","f","a","NONE","","","PSF cube file (only needed for stacked analysis)"))
109   - pars.append(gammalib.GApplicationPar("bkgcube","s","a","NONE","","","Background cube file (only needed for stacked analysis)"))
  109 + pars.append(gammalib.GApplicationPar("bkgcube","f","a","NONE","","","Background cube file (only needed for stacked analysis)"))
110 110 pars.append(gammalib.GApplicationPar("caldb","s","a","prod2","","","Calibration database"))
111 111 pars.append(gammalib.GApplicationPar("irf","s","a","South_50h","","","Instrument response function"))
112 112 pars.append(gammalib.GApplicationPar("edisp","b","h","no","","","Apply energy dispersion?"))
... ... @@ -114,16 +114,6 @@ class csspec(ctools.cscript):
114 114 pars.append(gammalib.GApplicationPar("emax","r","h","100.0","","","Upper energy limit for spectral points(TeV)"))
115 115 pars.append(gammalib.GApplicationPar("enumbins","i","a","20","","","Number of spectral points"))
116 116 pars.append(gammalib.GApplicationPar("ebinalg","s", "h","LOG","FILE|LIN|LOG","", "Algorithm for defining energy bins"))
117   - pars.append(gammalib.GApplicationPar("binned","b","a","no","yes|no","","Use binned analysis in each energy bin"))
118   - pars.append(gammalib.GApplicationPar("nebins","i","h","5","","","Number of energy bins per spectral point"))
119   - pars.append(gammalib.GApplicationPar("coordsys","s","a","CEL","CEL|GAL","","Coordinate System"))
120   - pars.append(gammalib.GApplicationPar("proj","s","a","CAR","AIT|AZP|CAR|MER|MOL|STG|TAN","","Projection method"))
121   - pars.append(gammalib.GApplicationPar("xref","r","a","83.63","0","360","First coordinate of image center in degrees (RA or galactic l)"))
122   - pars.append(gammalib.GApplicationPar("yref","r","a","22.01","-90","90","Second coordinate of image center in degrees (DEC or galactic b)"))
123   - pars.append(gammalib.GApplicationPar("nxpix","i","a","200","","","Size of the X axis in pixels"))
124   - pars.append(gammalib.GApplicationPar("nypix","i","a","200","","","Size of the Y axis in pixels"))
125   - pars.append(gammalib.GApplicationPar("binsz","r","a","0.02","","","Pixel size (deg/pixel)"))
126   - pars.append(gammalib.GApplicationPar("anumbins","i","h","200","","","Angular separation bins for PSF cube in binned mode"))
127 117 pars.append(gammalib.GApplicationPar("calc_ts","b","h","yes","yes|no","","Compute TS value in each bin"))
128 118 pars.append(gammalib.GApplicationPar("calc_ulim","b","h","yes","yes|no","","Compute upper limit in each bin"))
129 119 pars.append(gammalib.GApplicationPar("fix_srcs","b","h","yes","yes|no","","Fix other skymodel parameters"))
... ... @@ -139,10 +129,18 @@ class csspec(ctools.cscript):
139 129 """
140 130 Get parameters from parfile and setup the observation.
141 131 """
  132 +
142 133 # Set observation if not done before
143 134 if self.obs == None or self.obs.size() == 0:
144 135 self.require_inobs("csspec::get_parameters()")
145 136 self.obs = self.get_observations()
  137 +
  138 + # Check if we have one binned cta observation, i.e. if we are in binned mode
  139 + self.m_binned_mode = False
  140 + if self.obs.size() == 1:
  141 + if self.obs[0].classname() == "GCTAObservation":
  142 + if self.obs[0].eventtype() == "CountsCube":
  143 + self.m_binned_mode = True
146 144  
147 145 # Set models if we have none
148 146 if self.obs.models().size() == 0:
... ... @@ -151,25 +149,44 @@ class csspec(ctools.cscript):
151 149 # Get source name
152 150 self.m_srcname = self["srcname"].string()
153 151  
154   - # Get ebounds
155   - self.m_ebounds = self.create_ebounds()
156   -
157   - # Get binning flag
158   - self.m_binned = self["binned"].boolean()
159   -
  152 + # Get ebounds
  153 + if self.m_binned_mode:
  154 +
  155 + # Extract cube ebounds
  156 + cube_ebounds = self.obs[0].events().ebounds()
  157 +
  158 + # Read user parameters
  159 + self.m_enumbins = self["enumbins"].integer()
  160 + self.m_emin = self["emin"].real() - 1e-6 #rounding tolerance
  161 + self.m_emax = self["emax"].real() + 1e-6 #rounding tolerance
  162 +
  163 + # Compute number of energy layers used for one spectral bin
  164 + n_layers = cube_ebounds.size() / self.m_enumbins
  165 +
  166 + # Don't allow zero layers
  167 + if n_layers == 0:
  168 + n_layers = 1
  169 +
  170 + # Create new ebounds object and fill it
  171 + self.m_ebounds = gammalib.GEbounds()
  172 + for i in range(cube_ebounds.size())[::n_layers]:
  173 + if i + n_layers > cube_ebounds.size():
  174 + break
  175 +
  176 + # Get new emin and emax of bin
  177 + emin = cube_ebounds.emin(i)
  178 + emax = cube_ebounds.emax(i + n_layers - 1)
  179 +
  180 + # Append energy bin if within user range
  181 + if self.m_emin <= emin.TeV() and self.m_emax >= emax.TeV():
  182 + self.m_ebounds.append(emin, emax)
  183 +
  184 + # Unbinned mode
  185 + else:
  186 + self.m_ebounds = self.create_ebounds()
  187 +
160 188 # Get edisp flag
161 189 self.m_edisp = self["edisp"].boolean()
162   -
163   - if self.m_binned:
164   - self.m_xref = self["xref"].real()
165   - self.m_yref = self["yref"].real()
166   - self.m_nxpix = self["nxpix"].integer()
167   - self.m_nypix = self["nypix"].integer()
168   - self.m_binsz = self["binsz"].real()
169   - self.m_coordsys = self["coordsys"].string()
170   - self.m_proj = self["proj"].string()
171   - self.m_ebins = self["nebins"].integer()
172   - self.m_anumbins = self["anumbins"].integer()
173 190  
174 191 # Read other parameters
175 192 self.m_outfile = self["outfile"].filename()
... ... @@ -286,7 +303,8 @@ class csspec(ctools.cscript):
286 303 # Write header
287 304 if self.logTerse():
288 305 self.log("\n")
289   - self.log.header1("Generate spectrum")
  306 + self.log.header1("Generate spectrum")
  307 + self.log(str(self.m_ebounds))
290 308  
291 309 # Initialise FITS Table with extension "SPECTRUM"
292 310 table = gammalib.GFitsBinTable(self.m_ebounds.size())
... ... @@ -298,19 +316,19 @@ class csspec(ctools.cscript):
298 316  
299 317 # Create FITS table columns
300 318 energy = gammalib.GFitsTableDoubleCol("Energy", self.m_ebounds.size())
301   - energy.unit("TeV")
302 319 energy_low = gammalib.GFitsTableDoubleCol("ed_Energy", self.m_ebounds.size())
303   - energy_low.unit("TeV")
304 320 energy_high = gammalib.GFitsTableDoubleCol("eu_Energy", self.m_ebounds.size())
305   - energy_high.unit("TeV")
306 321 flux = gammalib.GFitsTableDoubleCol("Flux", self.m_ebounds.size())
307   - flux.unit("erg/cm2/s")
308 322 flux_err = gammalib.GFitsTableDoubleCol("e_Flux", self.m_ebounds.size())
309   - flux_err.unit("erg/cm2/s")
310 323 TSvalues = gammalib.GFitsTableDoubleCol("TS", self.m_ebounds.size())
311 324 ulim_values = gammalib.GFitsTableDoubleCol("UpperLimit", self.m_ebounds.size())
312   - ulim_values.unit("erg/cm2/s")
313 325 Npred_values = gammalib.GFitsTableDoubleCol("Npred", self.m_ebounds.size())
  326 + energy.unit("TeV")
  327 + energy_low.unit("TeV")
  328 + energy_high.unit("TeV")
  329 + flux.unit("erg/cm2/s")
  330 + flux_err.unit("erg/cm2/s")
  331 + ulim_values.unit("erg/cm2/s")
314 332  
315 333  
316 334 # Loop over energy bins
... ... @@ -334,8 +352,8 @@ class csspec(ctools.cscript):
334 352 energy_low[i] = (elogmean - emin).TeV()
335 353 energy_high[i] = (emax - elogmean).TeV()
336 354  
337   - # unbinned analysis
338   - if not self.m_binned:
  355 + # use ctselect for unbinned analysis
  356 + if not self.m_binned_mode:
339 357  
340 358 # Log information
341 359 if self.logExplicit():
... ... @@ -355,121 +373,33 @@ class csspec(ctools.cscript):
355 373 # Retrieve observation
356 374 obs = select.obs()
357 375  
358   - # Binned analysis
  376 + # use ctcubemask for binned analysis
359 377 else:
360 378  
361 379 # Header
362 380 if self.logTerse():
363   - self.log.header3("Binning events")
364   -
365   - # Bin events
366   - bin = ctools.ctbin(self.obs)
367   - bin["usepnt"] = False
368   - bin["ebinalg"] = "LOG"
369   - bin["xref"] = self.m_xref
370   - bin["yref"] = self.m_yref
371   - bin["binsz"] = self.m_binsz
372   - bin["nxpix"] = self.m_nxpix
373   - bin["nypix"] = self.m_nypix
374   - bin["enumbins"] = self.m_ebins
375   - bin["emin"] = emin.TeV()
376   - bin["emax"] = emax.TeV()
377   - bin["coordsys"] = self.m_coordsys
378   - bin["proj"] = self.m_proj
379   - bin.run()
380   -
381   - # Header
382   - if self.logTerse():
383   - self.log.header3("Creating exposure cube")
384   -
385   - # Create exposure cube
386   - expcube = ctools.ctexpcube(self.obs)
387   - expcube["incube"] = "NONE"
388   - expcube["usepnt"] = False
389   - expcube["ebinalg"] = "LOG"
390   - expcube["edisp"] = self.m_edisp
391   - expcube["xref"] = self.m_xref
392   - expcube["yref"] = self.m_yref
393   - expcube["binsz"] = self.m_binsz
394   - expcube["nxpix"] = self.m_nxpix
395   - expcube["nypix"] = self.m_nypix
396   - expcube["enumbins"] = self.m_ebins
397   - expcube["emin"] = emin.TeV()
398   - expcube["emax"] = emax.TeV()
399   - expcube["coordsys"] = self.m_coordsys
400   - expcube["proj"] = self.m_proj
401   - expcube.run()
402   -
403   - # Header
404   - if self.logTerse():
405   - self.log.header3("Creating PSF cube")
406   -
407   - # Create psf cube
408   - psfcube = ctools.ctpsfcube(self.obs)
409   - psfcube["incube"] = "NONE"
410   - psfcube["usepnt"] = False
411   - psfcube["ebinalg"] = "LOG"
412   - psfcube["edisp"] = self.m_edisp
413   - psfcube["xref"] = self.m_xref
414   - psfcube["yref"] = self.m_yref
415   - psfcube["binsz"] = self.m_binsz*20.0 # Slow Psf variation
416   - psfcube["nxpix"] = self.m_nxpix/20+1
417   - psfcube["nypix"] = self.m_nypix/20+1
418   - psfcube["enumbins"] = self.m_ebins
419   - psfcube["anumbins"] = self.m_anumbins
420   - psfcube["emin"] = emin.TeV()
421   - psfcube["emax"] = emax.TeV()
422   - psfcube["coordsys"] = self.m_coordsys
423   - psfcube["proj"] = self.m_proj
424   - psfcube.run()
425   -
426   - # Header
427   - if self.logTerse():
428   - self.log.header3("Creating background cube")
429   -
430   - # Create background cube
431   - bkgcube = ctools.ctbkgcube(self.obs)
432   - bkgcube["incube"] = "NONE"
433   - bkgcube["usepnt"] = False
434   - bkgcube["ebinalg"] = "LOG"
435   - bkgcube["xref"] = self.m_xref
436   - bkgcube["yref"] = self.m_yref
437   - bkgcube["binsz"] = self.m_binsz
438   - bkgcube["nxpix"] = self.m_nxpix
439   - bkgcube["nypix"] = self.m_nypix
440   - bkgcube["enumbins"] = self.m_ebins
441   - bkgcube["emin"] = emin.TeV()
442   - bkgcube["emax"] = emax.TeV()
443   - bkgcube["coordsys"] = self.m_coordsys
444   - bkgcube["proj"] = self.m_proj
445   - bkgcube.run()
  381 + self.log.header3("Selecting events")
446 382  
  383 + # Select layers
  384 + cubemask = ctools.ctcubemask(self.obs)
  385 + cubemask["regfile"] = "NONE"
  386 + cubemask["ra"] = "UNDEFINED"
  387 + cubemask["dec"] = "UNDEFINED"
  388 + cubemask["rad"] = "UNDEFINED"
  389 + cubemask["emin"] = emin.TeV()
  390 + cubemask["emax"] = emax.TeV()
  391 + cubemask.run()
  392 +
447 393 # Set new binned observation
448   - obs = bin.obs()
449   -
450   - # Set precomputed binned response
451   - obs[0].response(expcube.expcube(), psfcube.psfcube(), bkgcube.bkgcube())
452   -
453   - # Get new models
454   - models = bkgcube.models()
455   -
456   - # Fix background models if required
457   - if self.m_fix_bkg:
458   - for model in models:
459   - if not model.classname() == "GModelSky":
460   - for par in model:
461   - par.fix()
462   -
463   - # Set new models to binned observation
464   - obs.models(models)
  394 + obs = cubemask.obs()
465 395  
466 396 # Header
467 397 if self.logTerse():
468 398 self.log.header3("Performing fit")
469 399  
470 400 # Likelihood
471   - like = ctools.ctlike(obs)
472   - like["edisp"].boolean(self.m_edisp)
  401 + like = ctools.ctlike(obs)
  402 + like["edisp"] = self.m_edisp
473 403 like.run()
474 404  
475 405 # Skip bin if no event was present
... ... @@ -518,9 +448,9 @@ class csspec(ctools.cscript):
518 448 if self.m_calc_ts:
519 449 TS = source.ts()
520 450  
521   - # Compute Npred value
  451 + # Compute Npred value (only works for unbinned analysis)
522 452 Npred = 0.0
523   - if not self.m_binned:
  453 + if not self.m_binned_mode:
524 454 for observation in like.obs():
525 455 Npred += observation.npred(source)
526 456  
... ... @@ -529,8 +459,8 @@ class csspec(ctools.cscript):
529 459  
530 460 # Compute flux error
531 461 parvalue = source.spectral()[0].value()
532   - rel_error = source.spectral()[0].error()/parvalue
533   - e_flux = fitted_flux*rel_error
  462 + rel_error = source.spectral()[0].error() / parvalue
  463 + e_flux = fitted_flux * rel_error
534 464  
535 465 # Set values for storage
536 466 TSvalues[i] = TS
... ...
doc/source/reference_manual/csspec.rst
  View file @b2026eb
... ... @@ -19,8 +19,6 @@ of the source detection. Optionally, it also computes an upper flux limit
19 19 that is particularily useful in case that the source is not significantly
20 20 detected within an energy bin (hidden parameter ``calc_ulim``).
21 21  
22   -The script works on event lists only.
23   -
24 22 On output, the script will provide a FITS file with the fitted source
25 23 spectrum.
26 24  
... ... @@ -56,6 +54,9 @@ General parameters
56 54 ``irf [string]``
57 55 Instrumental response function.
58 56  
  57 +``(edisp = no) [boolean]``
  58 + Apply energy dispersion to response computation?
  59 +
59 60 ``emin [real]``
60 61 Lower energy limit of events (in TeV).
61 62  
... ... @@ -68,36 +69,6 @@ General parameters
68 69 ``(ebinalg = LOG) <FILE|LIN|LOG> [string]``
69 70 Algorithm for defining energy bins.
70 71  
71   -``binned [boolean]``
72   - Specifies whether a binned computation should be used.
73   -
74   -``(nebins = 5) [integer]``
75   - Number of bins per spectral point for binned analysis.
76   -
77   -``coordsys <CEL|GAL> [string]``
78   - Coordinate system (CEL - celestial, GAL - galactic).
79   -
80   -``proj <AIT|AZP|CAR|MER|MOL|STG|TAN> [string]``
81   - Projection method.
82   -
83   -``xref [real]``
84   - Right Ascension / Galactic longitude of image centre (J2000, in degrees).
85   -
86   -``yref [real]``
87   - Declination / Galactic latitude of image centre (J2000, in degrees).
88   -
89   -``nxpix [integer]``
90   - Size of the Right Ascension / Galactic longitude axis (in pixels).
91   -
92   -``nypix [integer]``
93   - Size of the Declination / Galactic latitude axis (in pixels).
94   -
95   -``binsz [real]``
96   - Pixel size (in degrees/pixel).
97   -
98   -``(anumbins = 200) [integer]``
99   - Number of angular separation bins.
100   -
101 72 ``(calc_ts = yes) [boolean]``
102 73 Compute TS for each spectral point?
103 74  
... ...
src/ctcubemask/ctcubemask.cpp
  View file @b2026eb
... ... @@ -349,6 +349,10 @@ void ctcubemask::save(void)
349 349 ***************************************************************************/
350 350 void ctcubemask::get_parameters(void)
351 351 {
  352 + // Initialise selection flags
  353 + m_select_energy = true;
  354 + m_select_roi = true;
  355 +
352 356 // If there are no observations in container then load them via user
353 357 // parameters
354 358 if (m_obs.size() == 0) {
... ... @@ -359,7 +363,7 @@ void ctcubemask::get_parameters(void)
359 363 std::string msg = "A valid file needs to be specified for the "
360 364 "\"inobs\" parameter, yet \""+
361 365 (*this)["inobs"].filename()+"\" was given."
362   - " Specify a vaild observation definition or "
  366 + " Specify a valid observation definition or "
363 367 "event list FITS file to proceed.";
364 368 throw GException::invalid_value(G_GET_PARAMETERS, msg);
365 369 }
... ... @@ -373,13 +377,36 @@ void ctcubemask::get_parameters(void)
373 377 m_regfile = (*this)["regfile"].filename();
374 378 m_usepnt = (*this)["usepnt"].boolean();
375 379 if (!m_usepnt) {
376   - m_ra = (*this)["ra"].real();
377   - m_dec = (*this)["dec"].real();
  380 +
  381 + // Check RA/DEC parameters for validity to read
  382 + if ((*this)["ra"].is_valid() && (*this)["dec"].is_valid()) {
  383 + m_ra = (*this)["ra"].real();
  384 + m_dec = (*this)["dec"].real();
  385 + m_select_roi = true;
  386 + }
  387 + else {
  388 + m_select_roi = false;
  389 + }
  390 + }
  391 +
  392 + // Check if radius is vaild for a RoI selection
  393 + if (m_select_roi && (*this)["rad"].is_valid()) {
  394 + m_rad = (*this)["rad"].real();
378 395 }
379   - m_rad = (*this)["rad"].real();
380   - m_emin = (*this)["emin"].real();
381   - m_emax = (*this)["emax"].real();
382   -
  396 + else {
  397 + m_select_roi = false;
  398 + }
  399 +
  400 + // Check for sanity of energy selection parameters
  401 + if ((*this)["emin"].is_valid() && (*this)["emax"].is_valid()) {
  402 + m_emin = (*this)["emin"].real();
  403 + m_emax = (*this)["emax"].real();
  404 + m_select_energy = true;
  405 + }
  406 + else {
  407 + m_select_energy = false;
  408 + }
  409 +
383 410 // Optionally read ahead parameters so that they get correctly
384 411 // dumped into the log file
385 412 if (read_ahead()) {
... ... @@ -410,6 +437,12 @@ void ctcubemask::apply_mask(GCTAObservation* obs)
410 437 GSkyMap map = cube->map();
411 438 const GEbounds& ebounds = cube->ebounds();
412 439  
  440 + // If no energy selection is required set energy boundaries to cube boundaries
  441 + if (!m_select_energy) {
  442 + m_emin = ebounds.emin().TeV();
  443 + m_emax = ebounds.emax().TeV();
  444 + }
  445 +
413 446 // Initialise energy selection
414 447 int n_ebin = ebounds.size();
415 448 int e_idx1 = 0;
... ... @@ -455,24 +488,26 @@ void ctcubemask::apply_mask(GCTAObservation* obs)
455 488 }
456 489  
457 490 // Set all pixels inside selected energy bands but outside ROI
458   - // to -1.0
459   - GSkyRegionCircle roi(m_ra, m_dec, m_rad);
460   - for (int i = e_idx1; i <= e_idx2; ++i) {
461   - for (int pixel = 0; pixel < npix; ++pixel) {
462   - GSkyDir dir = map.inx2dir(pixel);
463   - if (!roi.contains(dir)) {
464   - map(pixel,i) = -1.0;
  491 + // to -1.0 if requested
  492 + if (m_select_roi) {
  493 + GSkyRegionCircle roi(m_ra, m_dec, m_rad);
  494 + for (int i = e_idx1; i <= e_idx2; ++i) {
  495 + for (int pixel = 0; pixel < npix; ++pixel) {
  496 + GSkyDir dir = map.inx2dir(pixel);
  497 + if (!roi.contains(dir)) {
  498 + map(pixel,i) = -1.0;
  499 + }
465 500 }
466 501 }
467   - }
468 502  
469   - // Log selected energy band
470   - if (logTerse()) {
471   - log << gammalib::parformat("Selected ROI");
472   - log << "RA=" << m_ra << " deg, ";
473   - log << "DEC=" << m_dec << " deg, ";
474   - log << "Radius=" << m_rad << " deg";
475   - log << std::endl;
  503 + // Log selected energy band
  504 + if (logTerse()) {
  505 + log << gammalib::parformat("Selected ROI");
  506 + log << "RA=" << m_ra << " deg, ";
  507 + log << "DEC=" << m_dec << " deg, ";
  508 + log << "Radius=" << m_rad << " deg";
  509 + log << std::endl;
  510 + }
476 511 }
477 512  
478 513 // Set all pixels inside selected energy bands and inside exclusion
... ... @@ -538,6 +573,8 @@ void ctcubemask::init_members(void)
538 573 // Initialise protected members
539 574 m_obs.clear();
540 575 m_infiles.clear();
  576 + m_select_energy = true;
  577 + m_select_roi = true;
541 578  
542 579 // Return
543 580 return;
... ... @@ -563,8 +600,10 @@ void ctcubemask::copy_members(const ctcubemask&amp; app)
563 600 m_emax = app.m_emax;
564 601  
565 602 // Copy protected members
566   - m_obs = app.m_obs;
567   - m_infiles = app.m_infiles;
  603 + m_obs = app.m_obs;
  604 + m_infiles = app.m_infiles;
  605 + m_select_energy = app.m_select_energy;
  606 + m_select_roi = app.m_select_roi;
568 607  
569 608 // Return
570 609 return;
... ...
src/ctcubemask/ctcubemask.hpp
  View file @b2026eb
... ... @@ -90,8 +90,10 @@ protected:
90 90 double m_emax; //!< Upper energy
91 91  
92 92 // Protected members
93   - GObservations m_obs; //!< Observations container
94   - std::vector<std::string> m_infiles; //!< Input event filenames
  93 + GObservations m_obs; //!< Observations container
  94 + std::vector<std::string> m_infiles; //!< Input event filenames
  95 + bool m_select_energy; //!< Perform energy selection
  96 + bool m_select_roi; //!< Perform ROI selection
95 97 };
96 98  
97 99  
... ...
test/data/crab_bkgcube.fits.gz 0 → 100644
View file @b2026eb
No preview for this file type
test/data/crab_bkgcube.xml 0 → 100644
  View file @b2026eb
  1 +<?xml version="1.0" encoding="UTF-8" standalone="no"?>
  2 +<source_library title="source library">
  3 + <source name="Crab" type="PointSource">
  4 + <spectrum type="PowerLaw">
  5 + <parameter name="Prefactor" value="5.7" error="0" scale="1e-16" min="1e-07" max="1000" free="1" />
  6 + <parameter name="Index" value="2.48" error="0" scale="-1" min="0" max="5" free="1" />
  7 + <parameter name="Scale" value="0.3" scale="1e+06" min="0.01" max="1000" free="0" />
  8 + </spectrum>
  9 + <spatialModel type="SkyDirFunction">
  10 + <parameter name="RA" value="83.6331" scale="1" min="-360" max="360" free="0" />
  11 + <parameter name="DEC" value="22.0145" scale="1" min="-90" max="90" free="0" />
  12 + </spatialModel>
  13 + </source>
  14 + <source name="BackgroundModel" type="CTACubeBackground" instrument="CTA,HESS,MAGIC,VERITAS">
  15 + <spectrum type="PowerLaw">
  16 + <parameter name="Prefactor" value="1" error="0" scale="1" min="0.01" max="100" free="1" />
  17 + <parameter name="Index" value="0" error="0" scale="1" min="-5" max="5" free="1" />
  18 + <parameter name="Scale" value="1" scale="1e+06" free="0" />
  19 + </spectrum>
  20 + </source>
  21 +</source_library>
... ...
test/data/crab_cntmap_small.fits.gz 0 → 100644
View file @b2026eb
No preview for this file type
test/data/crab_expcube.fits.gz 0 → 100644
View file @b2026eb
No preview for this file type
test/data/crab_psfcube.fits.gz 0 → 100644
View file @b2026eb
No preview for this file type
test/test_cscripts.sh.in
  View file @b2026eb
... ... @@ -209,7 +209,6 @@ chmod +x @abs_top_srcdir@/cscripts/csspec.py
209 209 srcname="Crab" \
210 210 emin=0.1 \
211 211 emax=100.0 \
212   - binned=no \
213 212 outfile="spectrum.fits" > csspec.log
214 213 $ECHO -n "."
215 214 if [ -s "spectrum.fits" ]
... ... @@ -228,30 +227,22 @@ $ECHO &quot; ok&quot;
228 227 $ECHO -n "Test csspec (binned): "
229 228 chmod +x @abs_top_srcdir@/cscripts/csspec.py
230 229 @abs_top_srcdir@/cscripts/csspec.py \
231   - inobs="data/crab_events.fits" \
232   - inmodel="data/crab.xml" \
233   - caldb="irf" \
234   - irf="cta_dummy_irf" \
235   - enumbins=3 \
236   - srcname="Crab" \
  230 + inobs="data/crab_cntmap_small.fits.gz" \
  231 + expcube="data/crab_expcube.fits.gz" \
  232 + psfcube="data/crab_psfcube.fits.gz" \
  233 + bkgcube="data/crab_bkgcube.fits.gz" \
  234 + inmodel="data/crab_bkgcube.xml" \
  235 + enumbins=2 \
237 236 emin=0.1 \
238 237 emax=100.0 \
239   - binned=yes \
240   - nxpix=10 \
241   - nypix=10 \
242   - binsz=0.1 \
243   - nebins=2 \
244   - coordsys="CEL" \
245   - xref=83.63 \
246   - yref=22.01 \
247   - proj="CAR" \
248   - outfile="spectrum.fits" > csspec.log
  238 + srcname="Crab" \
  239 + outfile="spectrum_binned.fits" > csspec.log
249 240 $ECHO -n "."
250   -if [ -s "spectrum.fits" ]
  241 +if [ -s "spectrum_binned.fits" ]
251 242 then
252 243 $ECHO -n "."
253 244 else
254   - $ECHO " spectrum.fits file is not valid"
  245 + $ECHO " spectrum_binned.fits file is not valid"
255 246 exit 1
256 247 fi
257 248 $ECHO " ok"
... ... @@ -309,13 +300,13 @@ chmod +x @abs_top_srcdir@/cscripts/cslightcrv.py
309 300 nxpix=10 \
310 301 nypix=10 \
311 302 binsz=0.1 \
312   - outfile="lightcurve.fits" > cslightcrv.log
  303 + outfile="lightcurve_binned.fits" > cslightcrv.log
313 304 $ECHO -n "."
314   -if [ -s "lightcurve.fits" ]
  305 +if [ -s "lightcurve_binned.fits" ]
315 306 then
316 307 $ECHO -n "."
317 308 else
318   - $ECHO " lightcurve.fits file is not valid"
  309 + $ECHO " lightcurve_binned.fits file is not valid"
319 310 exit 1
320 311 fi
321 312 $ECHO " ok"
... ...