text:
# optics: general-purpose scintillator light response simulation code

## I Introduction

Scintillation detectors are today widely used in nuclear and particle physics experiments for detection of ionizing charged particles as well as photons and neutrons. These detectors are also often used in medical instrumentation, process control devices, waste management, personal protection and nuclear safeguarding .

In the simplest design, the scintillation detector consists of a scintillator assembly (a solid crystal, liquid or gaseous scintillator volume) and an optional lightguide or waveshifting guide viewed by a light-sensitive amplifying device such as a photomultiplier tube (PMT) . The fraction of energy deposited in the scintillator material by penetrating ionizing particles is absorbed and re-emitted within several nanoseconds in the form of visible or ultraviolet light, so-called fluorescent radiation. The scintillator material itself can be an organic crystal, plastics or liquid, inorganic crystal, scintillating gas or glass silicate . Scintillation light propagates through the scintillator medium and bounces off its surfaces. A fraction of the light eventually reaches the photosensitive area of the light-sensitive device where it is converted into an electrical signal. The generated current pulse is typically discriminated in a user-designed electronics circuit and integrated for a large number of events to produce the pulse-height energy spectrum.

The main considerations in particle detector design are the total light output and the detector timing response, with both factors influencing the energy resolution of the detector. The amount of scintillation light reaching the photosensitive surface should be maximized and should ideally be independent of the position and linearly dependent on the magnitude of the energy deposition. The timing response is determined by the decay constants of scintillating excitation and the size and geometry of the detector itself, as well as the characteristics of the photo-sensitive device. The light output of a scintillator primarily depends on the conversion efficiency of the deposited ionizing energy to scintillation photons. The light collected by the photosensitive surface is also dependent on the surface area, the efficiency of the light transport and the scintillator transparency to its own scintillation light. The match between spectral response of the photomultiplier, the scintillation light and the quantum efficiency of the photocathode are additional controlling factors determining photoelectron statistics.

The magnitude and linearity of the pulse-heights and timing resolution of the scintillator signal can be optimized by tuning the light transfer to the photosensitive surface. The optimization is executed by the appropriate optical treatment of the detector surfaces and by specific choices of shapes and dimensions of the scintillator/lightguide assembly and different wrapping reflectors.

The analytical solution of the light transfer and collection process can be obtained easily for only a few relatively simple geometries. If the shapes of detectors and lightguides are complex or irregular, or if the number of the detector modules is large, the Monte Carlo solution of the problem is the most practical and sometimes the only feasible calculational technique .

In order to facilitate the detector design several simple computer programs were developed in the past twenty-five years and were used to study photon transport through the elements of the scintillation detectors . These studies focused on the possibilities of improving the light collection and the uniformity of light response along the main detector axes. Refs.  compared the measurements of pulse-height spectra of the cosmic muons and radioactive sources that were moved along the detector surfaces with the simulation results. The measurements of the scintillator spatial light output nonuniformities are reported in Refs.  without comparisons to any matching simulations.

To our knowledge all published photon transport programs to date have treated the properties of rotationally symmetric or regular parallelepipedic detector shapes. The only publicly available Monte Carlo photon transport program that is flexible enough to treat a larger variety of detector shapes is the CERN Program Library GUIDE7 code .

We initially developed the optics library of subroutines to study the light transport in simple cylindrical and trapezoidal scintillator detectors . As our need to understand the realistic photon propagation in the more complicated detector geometries increased, we expanded the code to handle a wider range of detector shapes and included more realistic photon–detector boundary interactions.

In Section II we show how the optics subroutines define the detector geometry. We explain our method of specifying the physical regions, region boundaries, intersections of photon trajectories with the region interfaces and calculation of the normal vectors at the points of photon–surface interactions. Section III touches upon the problem of calculating the intersection coordinates between the directed lines and the planar, cylindrical, spherical, conical and parabolical surfaces. Section IV presents our algorithm for the photon transport through the detector volume. The overall program structure, including the required and optional input parameters and a user-controlled output format is described in Sec. V. The four progressively more complicated examples of the program applications are given in Sec. VI. An example of using the optics program in conjunction with the standard GEANT code that simulates energy depositions caused by ionizing radiation in particle detectors is demonstrated in Sec. VII. As an illustration of the physics calculation, we show how to simulate the pulse-height spectra of monoenergetic electrons and photons in a pure cesium iodide (CsI) calorimeter . Finally, the instructions for users outlining the required steps in the program installation are given in Sec. VIII.

## II Detector Geometries

### A Region Specification

We define a three-dimensional region by specifying the boundaries of the region and a point X0(3) inside the region. A convenient way of storing the region definitions is a FORTRAN record structure. While originally a record structure was a VAX FORTRAN extension to the ANSI Standard, nowadays almost all FORTRAN compilers provide these extensions . The record structure of a three-dimensional region is defined in the file REGION\_STRUCTURE.TXT:

```
C       Module REGION_STRUCTURE

C       |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
C               Record structure for 3-dimensional regions.
C       |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||

        PARAMETER       (MAX_SURF = 100)

        STRUCTURE /MATERIAL_STRUCTURE/
         REAL*4         REFRACTIVE_INDEX
         REAL*4         ATTENUATION_LENGTH
         REAL*4         SCATTERING_LENGTH
        END STRUCTURE

        STRUCTURE /REGION_STRUCTURE/

         REAL*4 X0(3)

         INTEGER*4      NPLANE
         INTEGER*4      NCONE
         INTEGER*4      NSPHERE
         INTEGER*4      NPARABOLOID

         RECORD /PLANE_STRUCTURE/       PLANE(MAX_SURF)
         RECORD /CONE_STRUCTURE/        CONE(MAX_SURF)
         RECORD /SPHERE_STRUCTURE/      SPHERE(MAX_SURF)
         RECORD /PARABOLOID_STRUCTURE/  PARABOLOID(MAX_SURF)

         RECORD /MATERIAL_STRUCTURE/    MATERIAL

        END STRUCTURE

C       |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||

C       End module REGION_STRUCTURE.

```

The maximum number of surfaces of a specific type is declared in the user parameter MAX\_SURF, while the actual number of surfaces of each type bounding a region is specified by the integer variables NPLANE, NCONE, NSPHERE or NPARABOLOID. The medium filling the region is indexed by an integer variable MATERIAL. For example, the assignments REFRACTIVE\_INDEX(1)=1.54, REFRACTIVE\_INDEX(2)=1.00 describe a scintillator medium in region 1 and air in region 2, respectively. A database of 17 different materials commonly encountered in scintillation detector design is kept in the file MATERIAL\_PROPERTIES.TXT. Presently defined media properties are the light attenuation length, scattering length, refractive index, atomic number, atomic weight, and the material volume density. A user can add additional materials and media properties or change the default values of the database entries.

Once the record structure of a region is defined, any number of variables having this structure format can be declared. Arrays of regions can be defined with the statement RECORD /REGION\_STRUCTURE/ REGION(100) that creates an array of 100 records, each with the structure format defined by the module REGION\_STRUCTURE above. The items in the structure of an individual region are referred to as “record.item”. The material of the region number 12 would therefore be referred to as “REGION(12).MATERIAL” while the $`y`$-coordinate of a point internal to the region 3 could be accessed through a variable REGION(3).X0(2).

Specific two-dimensional and three-dimensional shapes are defined as record structures in the file GEOMETRY\_STRUCTURES.TXT. This module defines the structures for a point, line, circle, plane, cylinder, cone, paraboloid, and sphere, as well as a general polygon structure specified by a set of vertices. The comments in the file explain the parameters defining each shape. A cylindrical volume is, for example, defined by the following piece of code:

```
        STRUCTURE /CYLINDER_STRUCTURE/
         REAL*4 X0 (3)  !Point on axis of cylinder.
         REAL*4 A  (3)  !Unit vector in direction of axis.
         REAL*4 RADIUS  !Radius of cylinder.
         INTEGER*4      TYPE    !Type of surface.
         RECORD /SURFACE_PROPERTIES_STRUCTURE/ PHYSICAL
        END STRUCTURE

```

Again, the user can add new geometrical shapes by defining additional structures conforming to the above format.

The program begins with a call to the subroutine DEFINE\_REGIONS where the parameters of the problem are either defined or read in from a user data file. An example input format for a cylindrical scintillator could be:

```
C Read in the geometry parameters  !Diameter of the front and back face
       READ(6,*) FRONT_D1, BACK_D1 !of region 1
       READ(6,*) Z1                !z-extent of region 1

C Define material types of regions
       REGION(1).MATERIAL=2 !Scintillator
       REGION(2).MATERIAL=1 !Air

C Define interior points of regions
       REGION(1).X0(3)=Z1/2.0
       REGION(2).X0(3)=Z1+2.0

C Define boundary surfaces of regions
       REGION(1).NPLANE=2
       REGION(1).NCONE=1

C Detector surfaces
       REGION(1).PLANE(1).TYPE=1
       REGION(1).PLANE(1).N(3)=-1.0
       REGION(1).PLANE(1).X0(3)=0.0

       REGION(1).PLANE(2).TYPE=1
       REGION(1).PLANE(2).N(3)=1.0
       REGION(1).PLANE(2).X0(3)=Z1

       REGION(1).CONE(1).TYPE=1
       REGION(1).CONE(1).A(3)=1
       REGION(1).CONE(1).SLOPE=(BACK_D1-FRONT_D1)/(2.0*Z1)

C Last volume defined contains all the others
       REGION(2).NSPHERE.=1
       REGION(2).SPHERE(1).TYPE=-2
       REGION(2).SPHERE(1).RADIUS=2.0*Z1

```

This piece of code defines a cylindrical scintillator with a z1 cm long $`z`$-axis and front and back face diameters FRONT\_D1 and FRONT\_D2, respectively. The detector is contained inside a spherical volume filled with air. The different surface types are described in Sec. II C.

### B Polygonal Geometries

The irregular polygonal volume can be described by a set of vertices stored in the POINT\_STRUCTURE record:

```
        STRUCTURE /POLYGON_STRUCTURE/
         RECORD /POINT_STRUCTURE/ VERTEX(100)   !Set of vertices.
        END STRUCTURE

```

This format gives the number of vertices of the front and back detector face in the first input line, and the $`x,y,z`$ coordinates of the front ($`z=0`$) and back detector face vertices in the second and third lines, respectively. For example, the truncated pentagonal detector geometry can be defined with a three-line entry:

```
5
2.45 3.37  0.20 3.55 0.0  0.52 7.11 0.0  0.52  8.21 3.37  0.20 5.33  5.46  0.0
0.00 6.24 22.02 2.04 0.0 22.61 8.63 0.0 22.61 10.66 6.24 22.02 5.33 10.09 21.66

```

The above input file describes a 22 cm long truncated pentagonal pyramid with the front (back) face side 3.56 (6.59) cm long.

### C Surface Specification

The recognized types of detector surfaces are distinguished by the integer variable TYPE. This parameter allows a user to specify the physical properties of the surface. The predefined surface types are perfectly transparent, perfectly reflecting, and perfectly absorbing interfaces as well as the “realistic” dielectric surface. These idealized surfaces are distinguished in the transport code by the parameter TYPE that is assigned 0, -1, -2 or positive integer values. A user should also set the optical properties of a surface through the structure SURFACE\_PROPERTIES\_STRUCTURE:

```
        STRUCTURE /SURFACE_PROPERTIES_STRUCTURE/
         REAL*4 ROUGHNESS           !Number describing roughness.
         REAL*4 REFLECTIVITY        !Fraction of incident light reflected.
         REAL*4 DIFFUSE_FRACTION    !Diffuse_reflection/total_reflection.
        END STRUCTURE

```

The REFLECTIVITY values could range from 0.0 (a perfect absorber) to 1.0 (a perfect reflector). The DIFFUSE\_FRACTION coefficient determines the ratio between the specular reflection (DIFFUSE\_FRACTION=0.0) and the diffuse scattering processes (DIFFUSE\_FRACTION=1.0, see Sec. IV B). The ROUGHNESS parameter describes the deviation of the detector sides from an ideal flat surface behavior (Sec. IV C).

It is important to note that a pair of adjacent regions that are not nested is always separated by the two surfaces. This arrangement allows for a boundary between two regions to be of a different type, depending on the region side from which a photon is approaching the interface. An example is a “two-way mirror” between two regions that could be created by making the surface of one region perfectly reflecting and the adjacent surface of the neighboring region completely transparent.

The four different detector geometries are described as examples in more detail in Sec. VI and shown in Fig. 1. The parameters of four geometries are specified in Table I. These shapes represent the real detectors whose light responses we studied both in Monte Carlo simulations and in calibration measurements using cosmic muons, monoenergetic tagged positron and photon beams, and laser pulses. The geometries are as varied as irregular truncated pyramids made of pure cesium iodide, plastic veto scintillator staves, a cylindrical plastic target, and a rectangular Plexiglas monitoring plate, respectively.

## III Intersection of Directed Lines with optics Surfaces

A straight line in three-dimensional space can be represented by the parametric equation:

$$𝑿\mathbf{=}𝑿_{\mathrm{𝟎}𝑳}\mathbf{+}𝒓𝑽\mathbf{,}$$
(1)
where $`𝑿_{\mathrm{𝟎}𝑳}`$ is a point on the line, $`𝑽`$ is a unit vector in the direction of the line, and the parameter $`𝒓`$ specifies the distance from a fixed point $`𝑿_{\mathrm{𝟎}𝑳}`$ to a point $`𝑿`$.

The parametric equations of the plane, spherical, and cylindrical surfaces are given by

$`𝒏\mathbf{}\mathbf{(}𝑿\mathbf{}𝑿_{\mathrm{𝟎}𝑷}\mathbf{)}`$ $`\mathbf{=}`$ $`\mathrm{𝟎}\mathbf{,}\text{ plane,}`$ (2)
$`\mathbf{|}𝑿\mathbf{}𝑿_{\mathrm{𝟎}𝑺}\mathbf{|}`$ $`\mathbf{=}`$ $`𝑹\mathbf{,}\text{ sphere,}`$ (3)
$`\mathbf{|}𝑿\mathbf{}𝑿_{\mathrm{𝟎}𝑪}\mathbf{}\mathbf{[}𝑨\mathbf{}\mathbf{(}𝑿\mathbf{}𝑿_{\mathrm{𝟎}𝑪}\mathbf{)}\mathbf{]}\mathbf{|}`$ $`\mathbf{=}`$ $`\mathrm{𝟎}\mathbf{,}\text{ cylinder,}`$ (4)
where $`𝒏`$ is a normal vector to the plane, $`𝑿_{\mathrm{𝟎}𝑺}`$ and $`𝑹`$ are the center and the radius of the sphere, and $`𝑿_{\mathrm{𝟎}𝑪}`$, $`𝑨`$ are the point on the cylinder axis and a unit vector in the direction of the detector $`z`$ axis, respectively.

Intersections of a straight line with these surfaces are found in, for example, Ref. . The formulas for the intersection coordinates are implemented in the optics subroutine SUBROUTINE FIND\_LINE\_REGION\_INTERSECTION(LINE,I\_REGION, X1,X2,INTERSECT,REGION,N\_REGIONS). The relevant regions should first be defined with calls to SUBROUTINE DEFINE\_REGIONS\_CRYSTAL(REGION,N\_REGIONS).

The LINE variable is specified by the three-dimensional array LINE.X0 and by the three cosines of directions LINE.V that could be calculated from the polar and azimuthal angles ($`\theta `$,$`\varphi `$):

$`\mathrm{𝙻𝙸𝙽𝙴}.𝚅(\mathrm{𝟷})`$ $`=`$ $`\mathrm{sin}\theta \mathrm{cos}\varphi ,`$ (5)
$`\mathrm{𝙻𝙸𝙽𝙴}.𝚅(\mathrm{𝟸})`$ $`=`$ $`\mathrm{sin}\theta \mathrm{sin}\varphi ,`$ (6)
$`\mathrm{𝙻𝙸𝙽𝙴}.𝚅(\mathrm{𝟹})`$ $`=`$ $`\mathrm{cos}\theta .`$ (7)
The variable I\_REGION is an index of a region that takes values from 1 to N\_REGIONS, while the REGION variable is the FORTRAN record structure described in Sec. II A.

INTERSECT is a logical variable that is .TRUE. if at least one intersection point is found and .FALSE. if there are no intersections between the line and user defined surfaces. The arguments X1 and X2 are three element arrays that contain the calculated coordinates of the intersection points.

In Fig. 2 we show the intersections of reconstructed cosmic muon trajectories with the surfaces of a cesium iodide scintillation crystal. The crystal’s geometrical shape is a truncated irregular hexagon labeled HEX–A (Fig. 1.i). Fig. 3 illustrates the cosmic muon trajectories intersecting three vertically stacked rectangular scintillator bars (Fig. 1.ii). In both cases the intersection coordinates have been calculated using the FIND\_LINE\_REGION\_INTERSECTION subroutine.

## IV Photon Tracking

### A Dielectric Reflection

The reflection and refraction of light from a perfect dielectric interface can described by the laws of geometrical optics if the scintillation photons are unpolarized. If the angles between the directions of the incident, reflected and transmitted photon with respect to the surface normal are denoted with $`\theta _i`$, $`\theta _r`$ and $`\theta _t`$, respectively, the law of reflection states that

$$\theta _r=\pi \theta _i,$$
(8)
while the law of refraction (Snell’s law) on the other hand requires that

$$n_1\mathrm{sin}\theta _i=n_2\mathrm{sin}\theta _t,$$
(9)
where $`n_1`$ and $`n_2`$ are the indices of refraction of the media in regions 1 and 2, respectively.

The user can change the program so that the full polarization of the propagating photon is taken into account. The reflection and transmission phenomena would then be described by the Fresnel formulas (see, for example, ).

### B Diffuse Reflection

A perfectly diffuse reflector is defined by the uniform scattering probability $`P(\theta )`$ into a unit solid angle $`d\mathrm{\Omega }`$, where $`\theta `$ denotes the angle between the direction of a reflected photon and the normal vector to the surface. The angular distribution of the scattered radiation is then isotropic. If the reflecting surface is a plane, the resulting distribution is often called the Lambert’s (or cosine) law :

$$I(\theta )=I_0\mathrm{cos}\theta .$$
(10)
### C Rough Surfaces

A rough surface can be modeled by introducing a set of small plane surfaces tangent to the original surface at suitable points. This description goes under the name of the facet model and is described in more detail in Ref.  where the model predictions are also compared with the experimental measurements.

In the simplest version of the facet model a rough surface is characterized by a single roughness parameter $`R`$. This parameter is zero for smooth scintillator surfaces with a perfect finish. For real detectors with uneven boundaries $`R`$ value is a positive number between 0.0 and 1.0. Polished rough surfaces with height variations of the order of a micron have the roughness coefficient on the order of 0.1.

The subroutine PERTURB\_UNIT\_VECTOR (V\_IN,ROUGHNESS,V\_OUT) calculates a random perturbation to the direction of an incident photon unit vector V\_IN following the reflection from a rough surface. The scalar product of the perturbed vector V\_OUT and the unperturbed photon direction $`\mathrm{cos}\theta `$ is equal to $`1/\sqrt{(1+\delta ^2)}`$, where $`\delta `$ has a Gaussian distribution, centered at zero, with the width given by the ROUGHNESS parameter $`R`$.

### D Absorption and Rescattering

The media properties are described by two parameters, SCATTERING\_LENGTH and ATTENUATION\_LENGTH. The code simplifies photon transport by assuming that the scattering length and the absorption length do not depend on the wavelength of the propagating light. The dispersion effects are neglected. If a user wants to consider dispersion effects, the changes should be made in the file optics.f.

## V Program Structure

### A Program Input

The program starts with the subroutine DEFINE\_REGIONS that reads a data file specifying the parameters of the detector. The call to that subroutine defines the record structures in the array REGION. A starting position of a scintillation photon is accepted next. The coordinate values can be read in from the data file PHOTONS.DAT or, alternatively, a user can select one of the predefined distributions:

1. a single point source, fixed $`x,y,z`$;
2. a through-going particle, range of $`x`$, fixed $`z`$, $`y=0`$;
3. a range of light sources on the detector axis, range of $`z`$, $`x=y=0`$;
4. a uniform light distribution throughout detector volume, range of $`x,y,z`$;
5. a uniformly distributed light sources in the center plane, range of $`x,z`$, $`y=0`$;
6. an aimed “pencil” of scintillating radiation, fixed $`x,y,z`$, and fixed directional cosines $`v_x,v_y,v_z`$;
7. a uniform distribution in the first half (first half-sector) of the detector, range of $`x,y,z`$.

All input parameters defining the detector and the scintillating light distribution, as well as the desired output format can also be entered from the tkoptics window if one is using the GUI version of the program (see Fig. 4).

### B Physics Subroutines

The initial region number is found next by calling the routine FIND\_REGION\_NUMBER. The program then steps through the following levels:

1. Find the next intersection of a photon trajectory with the predefined surface.

2. Test to see if a photon is absorbed or rescattered in transit.

2.1 If absorption or rescattering occurred, look at the current region number.

2.1.1 If a photon is in the wavelength shifter, find the point at which it is absorbed, generate a new random direction, and continue.

2.1.2 If a photon is not in the wavelength shifter, count it as absorbed and get a new particle.

2.2 If a photon is not absorbed, continue.

3. Propagate a photon to the next intersection.

4. Look into the next region to find its region number.

4.1 If the next region number is the same as the current region number, continue the photon transport and return to (1).

4.2 If we are dealing with a different region, continue to (5).

5. Look at the type of the surface a photon encounters.

5.1 If the surface is transparent, continue to the next intersection.

5.2 If the surface is reflecting, stay in the current region while changing the photon’s direction.

5.3 If the surface is absorbing, stop the photon propagation, count it as absorbed, and generate a new photon ab initio and return to step (1).

5.4 When an interface is a real surface, use the Snell’s law to determine if the photon undergoes reflection or refraction, change the direction appropriately, and step into the next region if a photon is refracted.

5.5 For the watch points, call the WATCHPOINT subroutine, then do (5.1).

5.6 For the detectors, call the DETECTOR subroutine, count a photon as detected, then do (5.3).

6. Return to (1).

This algorithm is represented by a box diagram in Fig. 5.

The following is the list of subroutines used in the photon propagation:

```
SUBROUTINE INITIALIZE ( REGION, N_REGIONS )
SUBROUTINE WRITE_STATS ( REGION, N_REGIONS )
SUBROUTINE GENERATE_PHOTON ( LINE, REGION, N_REGIONS )
SUBROUTINE LOG_PHOTON_GENERATION ( LINE )
SUBROUTINE LOG_INTERSECTION ( LINE, REGION, N_REGIONS, N_CURRENT,
                              R, X, OUTWARD_NORMAL, CLASS, ISURF, TYPE )
SUBROUTINE LOG_REFLECTION ( LINE, REGION, N_REGIONS, N_CURRENT,
                            R, X, OUTWARD_NORMAL, CLASS, ISURF, TYPE )
SUBROUTINE LOG_REFRACTION ( LINE, REGION, N_REGIONS, N_CURRENT,
                            R, X, OUTWARD_NORMAL, CLASS, ISURF, TYPE )
SUBROUTINE LOG_ABSORPTION ( LINE, REGION, N_REGIONS, N_CURRENT,
                            R, X, OUTWARD_NORMAL, CLASS, ISURF, TYPE )
SUBROUTINE LOG_ATTENUATION ( REGION, N_CURRENT, LINE, X )
SUBROUTINE DEFINE_REGION_HALF_SECTOR ( FILENAME,
SUBROUTINE DEFINE_REGIONS ( REGION, N_REGIONS )
SUBROUTINE DEFINE_REGIONS_PYRAMID ( REGION, N_REGIONS )
SUBROUTINE DEFINE_REGIONS_CONE ( REGION, N_REGIONS )
SUBROUTINE DETECTOR ( TYPE, LINE, RANGE, TIME )
SUBROUTINE WATCHPOINT ( TYPE, LINE, RANGE, TIME )
SUBROUTINE TERMINATE ( REGION, N_REGIONS )

```

The functions these subroutines perform are described in Sec. IV.

### C Program Input

The listing of the program input for the graphical user interface (GUI) version of the code is printed by default at the program termination:

```
 Enter geometry type: 1
 Ψ0 = Truncated cone
 Ψ1 = Pi-Beta CsI crystal
 Enter pyramid type: 0
 Ψ1Ψ= Pentagon
 Ψ2Ψ= Hexagon-A
 Ψ3Ψ= Hexagon-B
 Ψ4Ψ= Hexagon-C
 Ψ5Ψ= Hexagon-D
 Ψ6Ψ= Half Hexagon-D (#1)
 Ψ7Ψ= Half Hexagon-D (#2)
 Ψ0Ψ= User defined (will prompt for file name)
 Enter file name: /fs24/users/emil/tomo/tkoptics/hh1-d-mc.dat
 Enter refractive index of pyramid material: 2.10
 Enter attenuation length of pyramid material (cm): 200.0
 Enter scattering length of pyramid material (cm): 200.0
 Enter surface type for sides of pyramid: -1
 Ψ 1Ψ= normal dielectric
 Ψ-1Ψ= imperfect specular reflector
 Ψ-2Ψ= perfect absorber
 Ψ-3Ψ= imperfect diffuse reflector
 Enter reflectivity of this surface (0-1): 0.975
 Enter roughness of this surface (0-1): 0.0
 Enter diffuse fraction for this surface (0-1): 0.0
 Enter surface type for front face of pyramid: -1
 Ψ 1Ψ= normal dielectric
 Ψ-1Ψ= imperfect specular reflector
 Ψ-2Ψ= perfect absorber
 Ψ-3Ψ= imperfect diffuse reflector
 Enter reflectivity of this surface (0-1): 0.975
 Enter roughness of this surface (0-1): 0.0
 Enter diffuse fraction for this surface (0-1): 0.0
 Enter refractive index of air gap material: 1.000273
 Enter attenuation length of air gap material (cm): 1.0E+3
 Enter scattering length of air gap material (cm): 1.0E+3
 Enter surface type for pyramid wrapper: -3
 Ψ 1Ψ= normal dielectric
 Ψ-1Ψ= imperfect specular reflector
 Ψ-2Ψ= perfect absorber
 Ψ-3Ψ= imperfect diffuse reflector
 Enter reflectivity of this surface (0-1): 0.9
 Enter roughness of this surface (0-1): 0.0
 Enter diffuse fraction for this surface (0-1): 1.0
 Enter surface type for front face of wrapper: -3
 Ψ 1Ψ= normal dielectric
 Ψ-1Ψ= imperfect specular reflector
 Ψ-2Ψ= perfect absorber
 Ψ-3Ψ= imperfect diffuse reflector
 Enter reflectivity of this surface (0-1): 0.2
 Enter roughness of this surface (0-1): 0.0
 Enter diffuse fraction for this surface (0-1): 1.0
 Enter refractive index of optical joint material: 1.58
 Enter attenuation length of optical joint material (cm): 200.0
 Enter scattering length of optical joint material (cm): 200.0
 Enter refractive index of PMT window material: 1.458
 Enter attenuation length of PMT window material (cm): 100.0
 Enter scattering length of PMT window material (cm): 100.0
 ********* Xfront =   6.790000       2.865000       22.19750
 ********* AXIS =  0.0000000E+00  0.1182974      0.9929782
 ********* VECTOR =  8.4076812E-03 -5.3805402E-03 -6.0056243E-04
 ********* VECTOR =  0.0000000E+00  1.1829736E-02  9.9297822E-02
 ********* Xfront =   6.790000       2.876830       22.29680
 Enter output data type: 3
   0 = No output
   1 = Endpoint data
   2 = Tracks
   3 = Statistics
 Enter output file name: /fs24/users/emil/tomo/tkoptics/tkoptics.rz
 Enter numbers of bins: (NXBINS NYBINS NZBINS) 14 12 24
 Enter X limits of container volume: (MIN MAX) 0 14
 Enter Y limits of container volume: (MIN MAX) 0 12
 Enter Z limits of container volume: (MIN MAX) 0 24
 Enter starting distribution: 3
      0 = Single point (Fixed X,Y,Z)
      1 = Through-going particle (Range of X; fixed Z; Y=0)
      2 = Range of points along axis (Range of Z; X=Y=0)
      3 = Uniform throughout volume (Range of X,Y,Z)
      4 = Uniform in center plane (Range of X,Z; Y=0)
      5 = Aimed (Fixed X,Y,Z; Fixed VX,VY,VZ)
      6 = Uniform in first half-sector (Range of X,Y,Z)
 Enter number of photons: 5000000
 Number of photons processed =      5000000

 Number of photons attenuated =      1467024
 Number of photons absorbed =      2197182
 Number of photons overbounced =      169244

 Number of photons detected =     1166550
 Ψ(      489568 started backward )
 Ψ(      676982 started forward )

```

This output is an example of the optics simulation of a half-hexagonal CsI detector, Fig. 1.i.

### D Program Output

The user can choose three different output formats. In the ENDPOINTS format for each detected scintillation photon a record is written to an ASCII output file documenting:

1. $`x_0,y_0,z_0`$, the starting coordinates of a generated photon;
2. $`v_{x_0},v_{y_0},v_{z_0}`$, the components of the initial directional vector of the photon;
3. $`x_1,y_1,z_1`$, a point at which the photon hit the photosensitive surface;
4. $`v_{x_1},v_{y_1},v_{z_1}`$, the components of the photon’s direction vector when it hit the photosensitive surface;
5. $`t`$, the time elapsed from the photon creation to the moment of its detection;
6. $`n`$, the number of times a photon was reflected from the detector surfaces.

The program terminates by printing the photon transport statistics summary. The summary includes the number of detected photons, as well as fractions of attenuated, absorbed and “overbounced” photons (see the tail of Sec. V C listing). The limit to the number of photon bounces is set by a user. The default limit is 200 reflections. Events exceeding that limit are counted as overbounced and the propagation resumes with a new photon at the starting position.

The CERNLIB PAW macro program titled tkpawread.kumac can be used to read in the ASCII file produced by the optics code and make an HBOOK4 RZ file containing a summary Ntuple . The usage syntax is exec tkpawread infile outfile. The structure of the generated Ntuple is given in an example below:

```
          ********************************************************
          * NTUPLE ID=   10  ENTRIES= 233209   Endpoint Data     *
          ********************************************************
          *  Var numb  *   Name    *    Lower     *    Upper     *
          ********************************************************
          *      1     * xstart    * 0.180037E+00 * 0.138114E+02 *
          *      2     * ystart    * 0.146866E-03 * 0.585641E+01 *
          *      3     * zstart    * 0.561714E-02 * 0.225345E+02 *
          *      4     * vxstart   * -.999965E+00 * 0.999996E+00 *
          *      5     * vystart   * -.999988E+00 * 0.999991E+00 *
          *      6     * vzstart   * -.999990E+00 * 0.999998E+00 *
          *      7     * xstop     * 0.314337E+01 * 0.115353E+02 *
          *      8     * ystop     * 0.238158E+00 * 0.541021E+01 *
          *      9     * zstop     * 0.220957E+02 * 0.227119E+02 *
          *     10     * vxstop    * -.999665E+00 * 0.998784E+00 *
          *     11     * vystop    * -.990139E+00 * 0.999993E+00 *
          *     12     * vzstop    * -.994538E+00 * 0.999994E+00 *
          *     13     * time      * 0.976813E-02 * 0.302494E+02 *
          *     14     * bounces   * 0.000000E+00 * 0.109000E+03 *
          ********************************************************

```

The variable names correspond to the ASCII file variables. The number of Ntuple entries in this case is equal to the number of generated photons.

If a user wants to calculate the three-dimensional light collection probability function the STATISTICS output format is more appropriate. The variables written to an ASCII file and into the resulting Ntuple are:

1. x,y,z, the coordinates of the centers of the elementary cells into which the detector volume is subdivided;
2. counts, the number of detected photons originating from a given elementary cell;
3. edge, the flag indicating if an elementary cell is fully contained inside the detector volume (edge$`=`$1), or is partly outside the detector (edge$`=`$0);
4. thrown, the number of the generated photons having the origins within a given elementary cell;
5. time, the average time-of-flight of the photon (time elapsed between the generation and detection moment).

```
          ********************************************************
          * NTUPLE ID=   10  ENTRIES=   4032   Statistics        *
          ********************************************************
          *  Var numb  *   Name    *    Lower     *    Upper     *
          ********************************************************
          *      1     * x         * 0.500000E+00 * 0.135000E+02 *
          *      2     * y         * 0.500000E+00 * 0.115000E+02 *
          *      3     * z         * 0.500000E+00 * 0.235000E+02 *
          *      4     * counts    * 0.000000E+00 * 0.957900E+04 *
          *      5     * edge      * 0.000000E+00 * 0.100000E+01 *
          *      6     * thrown    * 0.000000E+00 * 0.407280E+05 *
          *      7     * time      * 0.000000E+00 * 0.463954E+01 *
          ********************************************************

```

The number of the Ntuple entries represents the number of elementary cells, and in the case above corresponds to the 14$`\times `$12$`\times `$24 cm<sup>3</sup> subdivision. The simulated detector was irregular trapezoidal pyramid shown in Fig. 1.i. The resulting histograms of the number of photon reflections as well as the distribution of the photon time-of-flight times are shown in Fig. 6. Converting the later histogram into a cumulative representation we get the intrinsic contribution of the photon arrival times to the timing lineshape of the detector system, Fig. 7.

Finally, the information on the individual photon trajectories could be captured with the TRACK output format:

```
 ******************************************************************
 * Ntuple ID = 10     Entries = 53424     Track Data
 *************************************************************************
 * Var numb * Type * Packing *    Range     *  Block   *  Name           *
 *************************************************************************
 *      1   * I*4  *         * [0,200]      * TRACK    * ntrack          *
 *      2   * I*4  *         *              * TRACK    * itrack(ntrack)  *
 *      3   * R*4  *         *              * TRACK    * xtrack(ntrack)  *
 *      4   * R*4  *         *              * TRACK    * ytrack(ntrack)  *
 *      5   * R*4  *         *              * TRACK    * ztrack(ntrack)  *
 *      6   * R*4  *         *              * TRACK    * vxtrack(ntrack) *
 *      7   * R*4  *         *              * TRACK    * vytrack(ntrack) *
 *      8   * R*4  *         *              * TRACK    * vztrack(ntrack) *
 *************************************************************************
 *  Block   *  Entries  * Unpacked * Packed *   Packing Factor           *
 *************************************************************************
 * TRACK    *  53424    * 5604     * Var.   *    Variable                *
 * Total    *    ---    * 5604     * Var.   *    Variable                *
 *************************************************************************
 * Blocks = 1            Variables = 8       Max. Columns = 1401         *
 *************************************************************************

```

The track number, starting and endpoint coordinates and the initial track direction of a photon are recorded in this file. Figure 8 illustrates the uniform distribution of the detected photons on the surface of a photocathode. The simulated detector was again the HEX-D1/2 shape shown in Fig. 1.i.

## VI Program Examples

The coordinate system is always defined with the $`z`$ axis pointing along the detector axis and the origin is placed at the front detector face. If the detector shape is polygonal, one vertex of the front detector face should always lie on the positive $`x`$ axis. All detector surfaces in examples below are normal dielectric interfaces. Some surfaces are wrapped in the specularly reflecting layer of aluminized Mylar or multiple layers of Teflon sheet. User-defined geometrical shapes are included in an encompassing sphere filled with air. The rectangular subvolume of the sphere is divided into a number of elementary cells with side lengths $`L`$, usually $``$1 cm.

### A Irregular Pentagonal, Hexagonal and Trapezoidal Scintillators

The coordinate system is as follows: the $`z`$ direction points along the axis of the crystal, with $`z=0`$ at the smaller (front) face, and $`z=L`$ at the larger (back) face. One vertex of the front face and one of the back face lie on the positive $`x`$ axis. The coordinates of the vertices for each crystal type are stored in the text data files. The nine different scintillator shapes are examined: four irregular hexagonal truncated pyramids (we label them HEX–A, HEX–B, HEX–C, and HEX–D), one regular pentagonal (PENT) and two irregular half-hexagonal truncated pyramids (HEX–D1 and HEX–D2), and two trapezohedrons (VET–1, VET–2) . The volumes of studied CsI crystals varied from 797 cm<sup>3</sup> (HEX–D1/2) to 1718 cm<sup>3</sup> (HEX–C), see Fig. 1.i. The size of the elementary cell was 1 cm<sup>3</sup>.

The geometry of the simulation is simple. The CsI scintillator crystal is wrapped in a layer of aluminized Mylar and attached via a 0.2 mm thick optical joint to a 1 mm thick photomultiplier window. Behind this window is a photocathode. The diameter of the photosensitive surface is 46 mm (HEX–D and VET) or 67 mm inches (PENT and HEX) . The entire detector is enclosed in a large absorbing sphere filled with air. The program output is written to a binary RZ file. The results of the calculation are the three-dimensional light collection probability functions illustrated on Figs. 9 and 10. The percentage probability that photons generated within two-dimensional $`x`$$`z`$ bins will generate a photoelectron on the PMT cathode is shown on top panels. The bottom panels are tomography measurements of the corresponding light responses for two representative HEX–A and HEX–D1 detectors.

### B Rectangular Plastic Veto

The cylindrical plastic hodoscope of the PIBETA detector  used in charged particle tracking consists of 20 rectangular staves, Fig 1.ii. The individual sections are optically isolated from each other by one layer of aluminized Mylar wrapping and are viewed by one inch phototubes coupled via lightguides on both ends. The simulated light response of the individual PMTs as a function of the axial position along the plastic stave is an example of the calculation of physical interest.

### C Cylindrical Plastic Target

A plastic scintillator cylindrical target is used in the center of the PIBETA detector  to stop a pion beam (Fig. 1.iii). The positrons from pion and muon decays emanating from the target center are detected subsequently in the plastic veto hodoscope and the CsI calorimeter. A typical optics calculation of the photon transport predicts the light response of the target as a function of the polar angle of the exiting positron.

### D Light-Distribution Plexiglas Plate

The lead glass calorimeter monitor for the RadPhi experiment E94-026 at JLAB  is a laser-based calorimeter monitoring system. The nitrogen laser excites a cylindrical plastic scintillator and a generated light is distributed via six fiber optics cables to the 1.2 cm thick Plexiglas sheet  facing the lead glass wall in a light-tight enclosure (Fig 1.iv). The sheet is oversized at 152 cm$`\times `$152 cm to minimize the optical nonuniformities. The two-dimensional response of the smaller prototype system predicted with the optics code, as viewed through lead glass modules is displayed on Fig. 11. The maximum predicted variation in the light output of the plate as a function of the position of an individual lead glass detector is $``$30 %.

## VII Energy Lineshape Monte Carlo

### A Interpolation of 3D Light Nonuniformity Function

The program optics can be used to calculate the light collection probability function $`f(iL,jL,kL)`$, where $`i,j,k`$ are integers and $`L`$ the step size, on a user-defined three-dimensional grid.

Starting from this three-dimensional probability function the code can calculate for an arbitrary point source of light $`(x,y,z)`$ inside the detector volume the light fraction reaching the photosensitive device. This is accomplished by means of a series of simple one-dimensional interpolations. Our method of choice is the bilinear interpolation on an elementary grid square. The accuracy will be improved if one uses not only the light output values but also the gradients and the cross derivatives of the light output probability changes. These higher order methods that result in a smoother behavior of our function go under the names of the bi-cubic interpolations or bi-cubic splines. Numerical recipes for these methods are given as FORTRAN subroutines and C functions in Ref. .

### B GEANT Simulation of Electromagnetic Showers

The three-dimensional distribution of energy deposited in the detector volumes by ionizing charged particles was calculated by the standard code GEANT . We considered monoenergetic positrons and photons of 70 MeV total energy incident on the PIBETA calorimeter made of pure CsI pyramids. In the user-written subroutine the energy deposited in every step was multiplied by the normalized value of the light collection nonuniformity function to yield the detected energy deposition. Sums are kept of both the total deposited energy and total detected light. Examples of predicted energy deposition lineshapes in the calorimeter are shown in Fig. 12. Four different simulations are plotted. The top panel shows a pulse-height spectrum for $`10^5`$ simulated photons incident on the CsI calorimeter from its center. The responses of both idealized uniform detector and the detector characterized with the optics-calculated nonuniformity functions are displayed. The bottom panel shows the difference in the response for the 70 MeV positron showers.

## VIII Installation and Testing

The computer code described in this paper is distributed in the form of uuencoded compressed tar archives named optics.uu (229 kb) and tkoptics.uu (2.4 Mb). These two files can be obtained via transfers from the URL http://pibeta.phys.virginia.edu/public\_html/optics. Same files can also be obtained directly from the authors. The UNIX utilities for the tar archive processing can be copied from the official GNU site at http://www.gnu.ai.mit.edu.

Installation of the optics software on a computer running the UNIX operating system is straightforward. Execution of the csh script uudecode taking optics or tkoptics as a parameter creates a compressed tar archive file with an extension tar.gz. The archive can be unpacked with the command gunzip. The unpacked optics library files should be placed under the /common/local/optics subdirectory. The extracted tkoptics files should be copied to the optics subdirectory on the user’s disk. If the user prefers the GUI version of the program, the installation of CERNLIB libraries, as well as of the Tk and Tcl libraries  is required. The publicly available distributions can be accessed, for example, at ftp://ftp.scriptics.com/pub/tcl/. The parameter file .tkoptics defining the input parameters echoed in the default optics window should be placed in the login directory.

The Tk and Tcl software should be installed under /common/local/tk and /common/local/tcl areas, respectively. The most recent CERN Program Library release containing subroutines used in the tkoptics program is available at the site ftp://asisftp.cern.ch.

The working area can be now named by the project’s title and should contain the required input files, like the optics.dat and tkoptics.dat. To verify that the program executes properly the user should remake the executable code by compiling and linking the source code with the make optics command . The Makefile file is included in the distribution. The program can then be started with the tkoptics command invoking a GUI interface window or it can be run in batch mode with the optics command.

## IX Acknowledgements

The authors wish to thank Micheal Sadler of the Abilene Christian University for the loan of the drift chamber tomography apparatus. Stefan Ritt and Penelope Slocum of the University of Virginia and David Lawrence of the Arizona State University have helped with the cosmic muon tomography measurements. Their help is gratefully acknowledged.

This work is supported and made possible by grants from the US National Science Foundation and the Paul Scherrer Institute.

(iii)

(iv)

FIGURE 1

FIGURE 2

FIGURE 3

FIGURE 4

FIGURE 5

FIGURE 6

FIGURE 7

FIGURE 8

(i)

(ii)

FIGURE 9

(i)

(ii)

FIGURE 10

FIGURE 11

(i)

(ii)

FIGURE 12


format: text

text:
# BeppoSAX observations of SGR 1806–20

## 1 Introduction

Soft gamma-ray repeaters (SGRs) are remarkable transient events characterized by brief ($`<`$ 1 s) and relatively soft (peak photon energy $``$20-30 keV) bursts of super-Eddington luminosity. Only four (or possibly five) SGRs are currently known (see, e.g., Hurley 2000). Soon after the discovery of the first three SGRs, the fact that they were all located within, or very close to, young supernova remnants (Vasisht et al. 1994, Kulkarni & Frail 1993), and the 8 s periodicity observed during the famous 1979 March 5 event from SGR 0525–66 (Mazets et al. 1979), suggested that SGRs involve some form of impulsive release of energy from a neutron star. This was recently confirmed with the discovery of periodicities in the 5-8 s range in the persistent X-ray counterparts of SGR 1806–20 and SGR 1900+14 (Kouveliotou et al. 1998, Hurley et al. 1999a), as well as in the bursting emission from SGR 1900+14 (Cline, Mazets & Golenetskii 1998).

The quiescent X–ray counterpart of SGR 1806–20 was discovered with the ASCA satellite (Murakami et al. 1994) within the radio nebula G10.0–0.3, most likely a supernova remnant (Kulkarni & Frail 1993). Until recently, SGR 1806–20 was thought to be associated to a luminous blue variable (LBV) star located in the non-thermal core of the radio nebula (van Kerkwijk et al. 1995). Flux variability and a change in the nebular morphology on a timescale of a few months, were also detected from the central part of the radio nebula and thought to be related to the activity from SGR 1806–20 (Frail et al. 1997). However, a more precise localization of this SGR (Hurley et al. 1999b) has shown that its position is incompatible with the LBV, which is probably responsible for the non-thermal radio emission, but may be unrelated to SGR 1806–20 .

Previous X–ray observations of SGR 1806–20 in the 1-10 keV energy range were obtained with the ASCA satellite by Sonobe et al. (1994). These authors reported a 2-10 keV luminosity of $``$1035 erg s-1 (for d=10 kpc) and the following spectral parameters: photon index $`\alpha _{ph}`$ = 2.2$`\pm `$0.2, NH = (6$`\pm `$0.2)$`\times `$1022 cm-2 and kT = 6.5$`\pm `$1.5 keV, NH = (5.2$`\pm `$0.2)$`\times `$1022 cm-2 in the case, respectively, of a power law and of a thermal bremsstrahlung fit.

## 2  Observations and Data Analysis

The BeppoSAX observations of SGR 1806–20 were performed in October 1998 and in March 1999 (see Table 1 for details). Most of our results were obtained with the Medium Energy Concentrator Spectrometers (MECS) instrument that covers the $``$ 1.6–10 keV energy range.

The MECS (Boella et al. 1997) is based on three position–sensitive gas–scintillation proportional counters placed in the focal plane of grazing incidence telescopes, providing images over a nearly circular field of view with $`28^{}`$ radius. The MECS is characterized by a total (three telescopes) effective area of $``$150 cm2 at 6 keV, a relatively good angular resolution (50% power radius of $`75^{\prime \prime }`$ at 6 keV, on-axis) and a moderate energy resolution (FWHM $``$8.5$`\sqrt{6/\mathrm{E}_{\mathrm{keV}}}`$%). Due to a failure in one of the MECS units, only two of them were available at the time of our observations.

The target was clearly detected on axis in both observations, with a net count rate in the two MECS units of $``$0.1 counts s-1.

### 2.1 Spectral Analysis

We describe the spectral results obtained in the 1.6-10 keV range by using only the MECS instrument. We verified that the inclusion of the data from the LECS instrument (Parmar et al. 1997), that extends the energy range down to $``$0.1 keV, did not significantly change the results of the fits. In fact the combination of its smaller effective area, shorter exposure times (see Table 1) and high absorption gives a small statistics in the LECS data for SGR 1806–20 . For both observations, the MECS counts for the spectral analysis were extracted from a circle of 4 centered at the source position and the background spectra were estimated from a concentric circular corona with radii $``$5 and 10. In fact the observations of standard background fields are not adequate in this case due to the presence of a substantial diffuse emission from the Galactic plane.

The spectra were equally well fit with a power law ($`\alpha _{ph}`$2) or with a thermal bremsstrahlung (kT$``$11 keV) absorbed by the interstellar medium (NH $``$6$`\times `$1022 cm-2). A blackbody spectrum gave unacceptable results (Obs A: $`\chi ^2`$ = 1.7 / 88 dof, Obs B: $`\chi ^2`$ =1.26 / 132 dof). Since the results of both observations were consistent with the same spectral shape and flux, we also performed a joint analysis of both observations to further restrict the uncertainties in the spectral parameters under the assumption of a constant source. All the derived parameters are reported in Table 2.

The observed flux for the power law best fit is $``$ 10-11 erg cm-2 s-1 in the 2-10 keV energy range. In the following we assume a distance of 15 kpc, based on the work by Corbel et al. (1997). After correcting for the effect of the interstellar absorption, this corresponds to an emitted luminosity of $``$4$`\times `$1035 erg s-1 (2-10 keV).

A comparison with the previous results reported by Sonobe et al. (1994) for the October 1993 ASCA observation seems to indicate that during the BeppoSAX observations the spectrum of SGR 1806–20 was slightly harder and the luminosity higher. Though the evidence for a spectral variation is only marginal and can probably be accounted for by systematic uncertainties related to the comparison between different instruments, the significant luminosity change (a factor $``$2) deserves a more accurate analysis owing to its important implications for the nature of SGR 1806–20 . We therefore reanalyzed the public ASCA data obtained from the HEASARC/GSFC on-line archives using the latest available response matrices.

### 2.2 Reanalysis of ASCA data

We analyzed the ASCA GIS data of the Otober 1993 observation. The data consist of two observations carried out on October 10 and 20. The spectra were accumulated from circular regions with radii 6 centered on the source position and rebinned in order to have at least 25 counts for each energy channel. As for the BeppoSAX data, we estimated the background spectra from a source free region in the same observation to properly take into account the contribution from the galactic plane X–ray emission. The power law model gave a very good fit without evidence for variations between the two ASCA observations. Combining the two data sets, we derive the following best fit values: $`\alpha _{ph}`$ = 2.25$`\pm `$0.15, NH = (6.0$`\pm `$0.5)$`\times `$1022 cm-2, unabsorbed flux (1.32$`\pm `$0.06)$`\times `$10-11 erg cm-2 s-1 in the 2-10 keV energy range. The values for the spectral index and absorption are similar to those of Sonobe et al. (1994), but we derive a slighlty greater flux, closer to the value measured with BeppoSAX.

### 2.3 Timing Analysis

The times of arrival of the MECS counts used in the spectral analysis were converted to the Solar System Barycenter and used in the timing analysis. The data were first analyzed with a standard folding technique and phase fitting was subsequently used to refine the estimate of the pulse period. The derived values are reported in Table 1. The light curve obtained in the longer observation (1999) is shown in Fig. 1 for different energy ranges. No significant variations in the shape of the folded light curve are visible as a function of time and/or energy band.

## 3 Discussion

Soft Gamma-ray Repeaters are usually interpreted in the context of the so called ”Magnetar” model (Duncan & Thompson 1992; Thompson & Duncan 1995) which is based on isolated and strongly magnetized (B$`\stackrel{>}{}`$1014 G) neutron stars. In this model the magnetic field is the main energy source, powering both the persistent X–ray (and particle) emission and the soft gamma-ray bursting activity. This involves internal heating, due to the magnetic field dissipation, and the generation of seismic activity. The latter is responsible for the soft $`\gamma `$-ray bursts, when the magnetic stresses in the neutron star crust shake the magnetosphere and accelerate particles.

If the spin-down in SGR 1806–20 is interpreted as purely due to magnetic dipole radiation losses, the observed values of P and $`\dot{P}`$ imply B$``$8 1014 G. Actually, this value is very likely an overestimate. In fact, the particle wind outflow, either continuous or in the form of strong episodic outbursts, also contributes significantly to the spin-down (Thompson & Blaes 1998). Harding et al. (2000) derived the relations to estimate magnetic field and spin-down age as a function of the particle wind duty cycle and luminosity. They showed that, for instance in SGR 1806–20 , a dipole magnetic filed of 3 1013 G is derived assuming a continuous particle wind with a luminosity of $``$1037 erg s-1.

The long term period evolution of SGR 1806–20 is shown in Fig. 2 (top left panel), where our results are plotted with all the previous measurements obtained with RXTE and ASCA (Kouveliotou et al. 1998). A linear fit to all the points gives $`\dot{P}`$=(8.27$`\pm `$0.18)$`\times `$10-11 s s-1. Woods et al. (2000) reported the results of an extensive monitoring of SGR 1806–20 with the RossiXTE satellite performed from February to August 1999. The period value we measured in March is consistent with their timing solution, which however is valid only on a limited time span. In fact our first period measurement (October 1998) is inconsistent with a backward extrapolation of the RossiXTE results. It is evident that significant variations around the average linear spin-down trend are present in SGR 1806–20 . According to Woods et al., the level of this ”timing noise” in SGR 1806–20 , is relatively larger than that expected from an extrapolation of that observed in radio pulsars.

When compared with the previous period measurements, our data show that, on a long term timescale, the period evolution of SGR 1806–20 is relatively stable, despite the bursting activity that occurred in 1997 and 1998. This is in marked contrast with the behavior of SGR 1900+14 (see Fig. 2, top right panel), that displayed much greater deviations from a constant spin-down rate. SGR 1900+14 showed two distinct periods (before and after the Summer of 1998) of nearly constant spin-down at a similar rate of $``$6.1 $`\times `$10-11 s s-1. The period increase from June to August 1998 could be due either to an enhanced spin-down rate related to the onset of bursting activity in June 1998 (Woods et al. 1999) or to a sudden discontinuous jump caused by the exceptional ”super burst” event of August 27 (Feroci et al. 1999). Though plausible explanations for both hypothesis have been advanced in the context of the magnetar model (Thompson et al. 1999) the lack of period measurements between June 9 and August 28 does not allow to discriminate between the two possibilities.

Although no Soft Gamma-ray Bursts were detected in our BeppoSAX observations, we know that SGR 1806–20 remained active during the last three years (see, e.g., Woods et al. 2000, Hurley et al. 1999b). SGR 1806–20 was also active during the 1993 ASCA observation that allowed its identification with the steady X–ray counterpart (Murakami et al. 1994). The lack of large spin-down variations in SGR 1806–20 indicates that the presence of ”normal” bursting activity is not a sufficient condition to generate a significant change in the spin-down torques. This suggests that the more likely explanation for the spin history of SGR 1900+14 is that of a discontinuous ”braking glitch” due to the 27 August event.

It is also interesting to compare SGR 1806–20 and SGR 1900+14 in terms of their luminosity and long term variability properties (see Table 3 and lower panels of Fig. 2). Our new observations of SGR 1806–20 show that this source is extremely stable in terms of spectral shape, despite the evidence of a flux increase with respect to the 1993 level. All the observations of SGR 1806–20 were obtained during periods of bursting activity. Its luminosity (for d=15 kpc) is about a factor $``$5-6 higher than that of SGR 1900+14 during its active period of September 1998. The latter source was even fainter (a factor $``$3) during observations in May 1997 (BeppoSAX, Woods et al. 1999) and May 1998 (ASCA, Hurley et al. 1999a), an extended period of quiescence in which no soft gamma-ray bursting activity was detected.

## 4 Conclusions

New observations of SGR 1806–20 with the BeppoSAX satellite and a reanalysis of archive ASCA data have shown that this source is extremely stable in terms of long term period evolution and spectral properties. We found evidence for a luminosity variation of only $``$20% between the 1993 value and the more recent BeppoSAX observations, both obtained while the source was in an active state. This stability, despite the fact that numerous bursts have been observed in the last three years from SGR 1806–20 , implies that the normal bursting activity, independent on its origin, does not significantly affect the overall (dipolar?) magnetic field configuration and/or the particle outflow stream causing the neutron star spin-down. Only exceptional events like the 27 August 1998 super outburst of SGR 1900+14 appear to cause significant jumps in the period evolution as well as a substantially increased X–ray emission.