added atomic loading function

Nathan Kelley-Hoskins · Jürgen Knödlseder · Nathan Kelley-Hoskins
1 parent cec6dca8
Showing 1 changed file with 90 additions and 0 deletions
pyext/GObservations.i
@@ -179,4 +179,94 @@ typedef GObservations::likelihood likelihood;
         for x in state[1]:
             self.append(x)
 }
+%pythoncode {
+    
+    def load_atomic_fits( self, fnames, bkg='powerlaw'):
+        """Function for loading atomic fits files, that have no outside references.
+           Each fits file should have all of the following tables:
+              EVENTS
+              TELARRAY
+              GTI
+              POINT SPREAD FUNCTION
+              EFFECTIVE AREA
+              ENERGY DISPERSION
+              BACKGROUND
+          
+          **Args:**
+            fnames : str or list of str, input fits filenames to load into self
+          
+          **Opts:**
+            bkg : str, type of spectral model to use with all backgrounds.
+                  can be one of: 
+                      constant, Const, GmodelSpectralConst,
+                      powerlaw, Plaw, GModelSpectralPlaw,
+                      BrokenPlaw, GModelSpectralBrokenPlaw,
+                      ExpPlaw, GModelSpectralExpPlaw
+        """
+        
+        import os
+        
+        if isinstance( fnames, str ) :
+          fnames = [ fnames ]
+        
+        """Loop over our input fits filenames"""
+        for ifits, fits in enumerate( fnames ) :
+          
+            ob = gammalib.GCTAObservation()
+            
+            """Convert unicode to a regular str"""
+            if isinstance(fits, bytes) :
+              fits = fits.decode('utf-8')
+            
+            ob.load( fits )
+            
+            """We need a unique identifier, so lets just use the 
+            full filename"""
+            ob.id( os.path.abspath( fits ) )
+            
+            """Load IRFs"""
+            rsp = gammalib.GCTAResponseIrf()
+            rsp.load_aeff(       fits )
+            rsp.load_psf(        fits )
+            rsp.load_edisp(      fits )
+            rsp.apply_edisp( True )
+            rsp.load_background( fits )
+            ob.response( rsp )
+            
+            """Add our GCTAObservation to the GObservations"""
+            self.append( ob )
+            
+            spec = -1
+            
+            if bkg in ['constant','Const','GModelSpectralConst'] :
+                spec = gammalib.GModelSpectralConst()
+            
+            elif bkg in [ 'powerlaw','Plaw', 'GModelSpectralPlaw'] :
+                pivot = gammalib.GEnergy( 1.0, 'TeV' )
+                spec = gammalib.GModelSpectralPlaw()
+                spec['PivotEnergy'].value( pivot.MeV() )
+                spec['PivotEnergy'].fix()
+                spec['Prefactor'  ].free()
+                spec['Prefactor'  ].value( 1.35 )
+                spec['Index'      ].free()
+                spec['Index'      ].min(  -10   )
+                spec['Index'      ].max(   10   )
+                spec['Index'      ].value( -0.05 )
+            
+            elif bkg in ['BrokenPlaw', 'GModelSpectralBrokenPlaw'] : 
+                spec = gammalib.GModelSpectralBrokenPlaw()
+            
+            elif bkg in ['ExpPlaw', 'GModelSpectralExpPlaw'] :
+                spec = gammalib.GModelSpectralExpPlaw()
+            
+            else :
+              raise ValueError( "For kwarg bkg='%s', '%s' is not a recognizable spectral model." % (bkg,bkg) )
+            
+            back = gammalib.GCTAModelIrfBackground( spec )
+            back.name( '%sbkg_%s' % ( bkg, ob.id() ) )
+            back.ids( ob.id() )
+            self.models().append( back )
+          
+}
+
 };