Starting with Fermi/LAT

From WikiVirgo

Jump to: navigation, search

Руководство по обработке данных Fermi/LAT

(Автор: Zho-zho)

В этом руководстве я расскажу о базовых приёмах обработки данных инструмента LAT, установленного на космическом телескопе Fermi Предполагается, что программное обеспечение Fermi Science Tools уже установлено и инициализировано

Цель данного руководства состоит не в том, чтобы дать исчерпывающее описание всех процедур, связанных с обработкой данных инструмента Fermi/LAT, а в том чтобы осветить основные моменты, связанные с ней, поделится опытом и предложить решения основных проблем, возникающих перед человеком, который использовать данные Fermi/LAT в своей работе.

Contents

О Fermi/LAT

Космический телескоп Fermi был запущен 11-го июня 2008-го года. Изначально он назывался GLAST, но позднее был переименован в Fermi. Два главных инструмента на его борту — LAT (Large Area Telescope) и GBM (GLAST Burst Monitor). LAT используется для обычных наблюдений источников, в то время как GBM — в основном для наблюдений гамма-всплесков. В этом руководстве я расскажу, как обрабатывать данные LAT (часто пишут Fermi/LAT).

Итак. LAT наблюдает небо в диапазоне энергий от 30 МэВ и до ~ 1 ТэВ, хотя «нормальным» для него считается диапазон 100 МэВ - 300 ГэВ. Ниже 100 МэВ возрастающий фон крайне затрудняет обработку данных, а выше 300 ГэВ возможны проблемы с калибровкой. Fermi был специально запущен таким образом, чтобы постоянно сканировать всё небо: полный оборот вокруг своей оси он делает за ~1.6 часа, а обзор всего неба - за два оборота, т.е. 3.2 часа. Телескоп иногда перестаёт вращаться и наводится только на один объект, но такое случается только в экстраординарных случаях (например, вспышка «стандартной свечи» - пульсара в Крабовидной туманности в сентябре 2010-го года).

Как известно, любой телескоп строит изображение точечного источника не в виде точки, а немного его размазывает. Профиль этого размазывания носит название PSF (от английского Point-Spread Function), и его ширина определяет разрешающую способность телескопа. В случае Fermi/LAT PSF очень сильно зависит от энергии — θPSF = 2*(E / 1 ГэВ) − 0.8 градуса (95% доверительный интервал), так что если на 100 МэВ-ах он составляет 10 градусов, то на 100 ГэВ-ах — только 0.2 градуса.

Эти особенности спутника следует учитывать при анализе полученных с него данных.


Загрузка данных

Данные Fermi/LAT можно скачать по следующей ссылке:

http://fermi.gsfc.nasa.gov/cgi-bin/ssc/LAT/LATDataQuery.cgi

Страничка загрузки достаточно простая — я думаю, ни у кого с ней сложностей не возникнет. ;-) Для анализа необходимо скачать файл с зарегистрированными событиями (photon file) и файл, содержащий историю положения и ориентации спутника в пространстве (spacecraft file). Следует иметь ввиду, что если содержимое файла событий зависит от участка неба, который вы выбрали, то содержимое файла спутника зависит только от времени, так что если нужно проанализировать наблюдения двух объектов за один и тот же промежуток времени, то файл истории спутника достаточно скачать лишь один раз.

Давайте, для примера, скачаем данные для небезызвестного пульсара в Крабовидной туманности, например с 1-го августа по 24-е сентября 2010-го года. Для этого в поле Object «Name Or Coordinates» указываем Crab, а в поле «Observations Dates» пишем 2010-08-01 00:00:00, 2010-09-24 00:00:00. Скачаем данные для интервала энергий 1-300 ГэВ, для чего в поле «Energy Range» укажем 1000,300000. Внизу странички укажем, что именно мы хотим скачать — ставим галочки возле «Photon data» и «Spacecraft data».

Готово. Нажимаем Start search. На следующей страничке сервер порадует нас сообщением, что запрос успешно принят, сообщит краткую информацию о нашем запросе и предложит ссылку, по которой можно будет скачать данные. Эту краткую информацию недурственно где-нибудь сохранить, так как она может позже пригодиться — как минимум напомнить, что именно мы скачивали. Я обычно сохраняю её в файл Query.txt, который в нашем случае будет содержать следующее:

Search Center (RA,Dec) 	=	(83.6331,22.0145)
Radius 	=	15 degrees
Start Time (MET)	=	302313600 seconds (2010-08-01T00:00:00)
Stop Time (MET)	=	306979200 seconds (2010-09-24T00:00:00)
Minimum Energy 	=	1000 MeV
Maximum Energy 	=	300000 MeV

Теперь можно проследовать по предложенной ссылке, дождаться пока сервер закончит формировать заказанные нами файлы и скачать их.

Инициализация Science Tools (только в VIRGO классе Физ.фак.-а)

Fermi Science Tools на компьютерах VIRGO физического факультета настроено для работы как с оболочкой Bash (скрипт login_fermi.sh), так и с оболочкой Csh (скрипт login_fermi.csh). В дальнейшем все примеры показаны для оболочки Bash. Для перехода в Bash и инициализации достаточно набрать следущее:

[yevgen@Hyperion Fermi]$ bash
[yevgen@Hyperion Fermi]$ . /virgo/scripts/login_fermi.sh

Применяем фильтры

Сервер Fermi называет файлы совершенного ужасающими именами вроде L100924055816E0D2F37E91_PH00.fits и L100924055816E0D2F37E91_SC00.fits. Настоятельно рекомендую сменить их на что-то менее мозгодробильное, вроде ph.fits и sc.fits. Кстати, если данных много, сервер Fermi любит разбивать их на небольшие файлы (после нашего переименования пускай это будут ph0.fits, ph1.fits, ph2.fits, ph3.fits и т. д.). Это может быть удобно для скачивания, но здорово мешает при обработке, потому лично я их всегда склеиваю в один файл ph.fits следующим образом:

mv *PH00* ph0.fits
mv *PH01* ph1.fits
mv *PH02* ph2.fits
mv *PH03* ph3.fits
mv *SC* sc.fits
echo -n "Processing EVENTS extension... "
ftmerge 'ph0.fits[events],ph1.fits[events],ph2.fits[events],ph3.fits[events]' \
 events.fits clobber=yes lastkey="TStop"
fdelhdu events.fits[2] confirm=no proceed=yes
echo "Done!"
echo -n "Processing GTI extension... "
ftmerge 'ph0.fits[gti],ph1.fits[gti],ph2.fits[gti],ph3.fits[gti]' gti.fits \
 clobber=yes
fdelhdu gti.fits[1] confirm=no proceed=yes
echo "Done!"
echo "Combining files..."
fappend gti.fits[gti] events.fits
echo "Removing temporary files..."
mv -f events.fits ph.fits
rm -f gti.fits
echo "Done!"

Вы можете облачить это в какой нибудь скрипт (я так и сделал) и запускать по мере необходимости.

Теперь можно собственно и отфильтровать наши данные. Зачем это необходимо? В первую очередь это нужно сделать потому что не все события, зарегистрированные спутником, соответствуют фотонам от источника — есть ещё и космические лучи и прочий (возможный) «мусор», который туда попадает. Fermi сам оценивает, какие именно события чему соответствуют, но поскольку он не знает, что нам понадобится, то в фотонный файл попадает всё. Вообще рекомендации касательно фильтров событий можно найти здесь:

http://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/Cicerone/Cicerone_Data_Exploration/Data_preparation.html

Фотонам соответствует класс событий 2, именно его мы и отфильтруем.

[yevgen@Hyperion Fermi]$ gtselect evclsmin=2 evclsmax=2
Input FT1 file[] ph.fits
Output FT1 file[] ph_f.fits
RA for new search center (degrees) (0:360) [] 83.6331
Dec for new search center (degrees) (-90:90) [] 22.0145
radius of new search region (degrees) (0:180) [] 10
start time (MET in s) (0:) [] 302313600
end time (MET in s) (0:) [] 306979200
lower energy limit (MeV) (0:) [] 1000
upper energy limit (MeV) (0:) [] 300000
maximum zenith angle value (degrees) (0:180) [] 105
Done.

Значения координат, начального и конечного моментов времени соответствуют нашему запросу (вот зачем был нужен Query.txt).

Теперь нужно учесть, что спутник мог работать не всё время. Возможные перебои в штатной работе инструмента решаются путём использования GTI (Good Time Interval). Чтобы дописать GTI к нашим данным нужно использовать команду gtmktime:

[yevgen@Hyperion Fermi]$ gtmktime
Spacecraft data file[] sc.fits
Filter expression[] (IN_SAA!=T) && (DATA_QUAL==1)
Apply ROI-based zenith angle cut[] yes
Event data file[] ph_f.fits
Output event file name[] ph_f_g.fits

Готово! Теперь можно приступать собственно к анализу данных.

Снимки неба

Что ж, давайте посмотрим, как же выглядит наш пульсар в гамма-диапазоне. Для этого воспользуемся командой gtbin, которой «скормим» отфильтрованные данные:

[yevgen@Hyperion Fermi]$ gtbin
This is gtbin version ScienceTools-v9r17p0-fssc-20100906
Type of output file (CCUBE|CMAP|LC|PHA1|PHA2) [] CMAP
Event data file name[] ph_f_g.fits
Output file name[] ph_cMap.fits
Spacecraft data file name[] sc.fits
Size of the X axis in pixels[] 220
Size of the Y axis in pixels[] 220
Image scale (in degrees/pixel)[] 0.1
Coordinate system (CEL - celestial, GAL -galactic) (CEL|GAL) [] CEL
First coordinate of image center in degrees (RA or galactic l)[] 83.6331
Second coordinate of image center in degrees (DEC or galactic b)[] 22.0145
Rotation angle of image axis, in degrees[] 0
Projection method e.g. AIT|ARC|CAR|GLS|MER|NCP|SIN|STG|TAN:[] CAR


Получившуюся картинку можно посмотреть программой ds9.

Image:Crab_image.png

Разумеется, её можно ещё немного причесать. На энергии в 1 ГэВ 68% фотонов, благодаря PSF, размываются в кружок 0.8 градуса, так что картинку стоит смазать гаусианой (из доступного в ds9 она ближе всего к истинной форме PSF) с ядром размером в 8 пикслей — ведь мы строили картинку в масштабе 0.1 градуса/пиксель. Кроме того, давайте для понятности добавим координатную сетку, скажем в галактических координатах.

Image:Crab_image_smoothed.png

Полученная в итоге картинка заметно более радует глаз, на ней даже видно плоскость галактики.

Строим кривую блеска

Для того, чтобы построить кривую блеска, нужно выбрать только те фотоны, которые пришли от источника (хотя бы в своём большинстве). Для тех, кто уже строил кривые блеска с XMM-Newton замечу, что процедура подобная до слёз. Итак, будем строить кривую блеска в интервале энергий 1-300 ГэВ. На энергии в 1 ГэВ PSF Fermi/LAT составляет 2 градуса (по 95-ти процентам фотонов) и с повышением энергии только падает, так что достаточно будет собрать фотоны с кружёчка радиусом в 2 градуса вокруг нашего пульсара.

[yevgen@Hyperion Fermi]$ gtselect
Input FT1 file[ph.fits] 
Output FT1 file[ph_f.fits] 
RA for new search center (degrees) (0:360) [83.6331] 
Dec for new search center (degrees) (-90:90) [22.0145] 
radius of new search region (degrees) (0:180) [10] 2
start time (MET in s) (0:) [302313600] 
end time (MET in s) (0:) [306979200] 
lower energy limit (MeV) (0:) [1000] 
upper energy limit (MeV) (0:) [300000] 
maximum zenith angle value (degrees) (0:180) [105] 
Done.

Как видите, параметры указанные при предыдущем вызове команды запоминаются, так что можно явно указать только то, что в нужно изменить, в данном случае — радиус кружёчка. Теперь применим фильтры.

[yevgen@Hyperion Fermi]$ gtmktime
Spacecraft data file[sc.fits]
Filter expression[(IN_SAA!=T) && (DATA_QUAL==1)]
Apply ROI-based zenith angle cut[yes]
Event data file[ph_f.fits]
Output event file name[ph_f_g.fits]

И наконец построим кривую блеска, скажем, с шагом в 5 дней (432000 секунд). Для этого воспользуемся командой gtbin с параметром LC (что обозначает Light Curve — кривая блеска):

[yevgen@Hyperion Fermi]$ gtbin
This is gtbin version ScienceTools-v9r17p0-fssc-20100906
Type of output file (CCUBE|CMAP|LC|PHA1|PHA2) []  LC
Event data file name[]  ph_f_g.fits
Output file name[]  ph_lc.fits
Spacecraft data file name[]  sc.fits
Algorithm for defining time bins (FILE|LIN|SNR) []  LIN
Start value for first time bin in MET[]  302313600
Stop value for last time bin in MET[]  306979200
Width of linearly uniform time bins in seconds[] 432000

Теперь кривую блеска можно нарисовать. Вот что у нас получилось:

Image:Crab_lc.png

Как видим, файл с кривой блеска — ph_lc.fits — содержит в себе набор колонок:

Time — время каждого бина

TimeDel — ширина бина

Counts — число фотонов в бине

Error — ошибка числа фотонов в бине.

Чего же не хватает этой кривой блеска? Как видим, единицы измерения там — отсчёты (counts). Это величина сугубо инструментальная, и для того чтобы привязать её к чему-то её нужно поделить на эффективную площадь телескопа и время, в течении которого телескоп наблюдал объект. Всё было бы просто, если бы в каждом бине спутник был ориентирован по отношению к объекту одинаковым образом и наблюдал бы его одинаковое время, но, к сожалению, это, вообще говоря, не так. Выйти из этого положения нам поможет программа gtexposure, которая допишет в наш файл с кривой блеска экспозицию — произведение эффективной площади на время наблюдения для каждого бина.

Для правильного вычисления экспозиции, вообще говоря, необходимо знать спектр источника. Gtexposure считает что спектр имеет степенной вид, и потому спросит нас о спектральном индексе. Средний спектральный индекс нашего пульсара составляет 2.2, потому так и укажем:

[yevgen@Hyperion Fermi]$ gtexposure
Light curve file[] ph_lc.fits
Spacecraft file[] sc.fits 
Response functions[] P6_V3_DIFFUSE
Source model XML file[] none
Photon index for spectral weighting[] -2.2

Теперь можно нарисовать кривую блеска в «правильных» единицах — фотонах/(см2*сек).

Image:Crab_lc_exp_corr.png

Вот и всё, кривая блеска готова. Как видим, в конце наблюдения пульсар резко вспыхнул. Кстати говоря, эта вспышка (22-го сентября 2010-го года) вызвала немало шума, так как этот пульсар считался своеобразной «стандартной свечёй» и по нему калибровали многие телескопы.

Тут следует заметить, что по-хорошему необходимо построить такую же кривую блеска для фона — соседнего участка неба свободного от источников — чтобы убедиться, что наблюдаемые изменения потока от источника реальны и не обусловлены колебаниями фона. Но в данном конкретном случае, смею вас заверить, что вспышка настоящая. 8-]

Personal tools