«Цветков Егор Александрович МЕТОДЫ РАСЧЕТА НЕАДДИТИВНЫХ ФУНКЦИОНАЛОВ В ПРИКЛАДНЫХ ЗАДАЧАХ РАДИАЦИОННОЙ ФИЗИКИ 05.13.18 — математическое моделирование, численные методы и комплексы программ Диссертация на соискание ученой ...»

Министерство образования и наук

и Российской Федерации

Московский физико-технический институт

(государственный университет)

На правах рукописи

УДК 519.245

Цветков Егор Александрович

МЕТОДЫ РАСЧЕТА НЕАДДИТИВНЫХ ФУНКЦИОНАЛОВ В

ПРИКЛАДНЫХ ЗАДАЧАХ РАДИАЦИОННОЙ ФИЗИКИ

05.13.18 — математическое моделирование, численные методы и комплексы программ Диссертация на соискание ученой степени кандидата физико–математических наук

Научный руководитель — к.т.н. Шаховский Валентин Владимирович Автор:

Москва – Оглавление Введение 1 Математические модели многодетекторных приборов 1.1 Свойство аддитивности..................... 1.1.1 Аддитивность по столкновениям.............. 1.1.2 Аддитивность по траекториям............... 1.1.3 Неаддитивные функционалы................ 1.1.4 Аддитивность по времени.................. 1.1.4.1 Определение...................... 1.1.4.2 Свойства......................... 1.2 Требования к разрабатываемым математическим моделям. 1.3 Моделирование сцинтилляционного детектора........ 1.3.1 Выбор модели сцинтилляционного детектора....... 1.3.2 Вычисление потерянной энергии.............. 1.3.3 Входные и выходные данные моделей детектора..... 1.3.4 Формулировка используемого алгоритма......... 1.3.4.1 Выделение точек пересечения частицами границ детектора........................ 1.3.4.2 Случай = ().................... 1.3.4.3 Случай = (1, 2,..., )............. 1.4 Моделирование спектрометра................. 1.5 Моделирование электронной схемы совпадений....... 1.6 Моделирование интроскопа, работающего по методу меченых нейтронов.......................... 1.6.1 Метод меченых нейтронов................. 1.6.2 Моделирование прибора................... 1.7 Моделирование комптоновского гамма-спектрометра.... 1.7.1 Приложения и принцип работы комптоновских гаммаспектрометров........................ 1.7.2 Модель комптоновского гамма-спектрометра....... 1.7.3 Определение типа совпадений............... 1.8 Результаты............................ 2 Весовой алгоритм вычисления среднего неаддитивных функционалов 2.1 Существующие способы вычисления среднего небольцмановских функционалов..................... 2.2 Концепция супертреков..................... 2.3 Формулировка предлагаемого алгоритма........... 2.4 Фазовое пространсво...................... 2.5 Вероятностное пространство.................. 2.5.1 Способы описания ветвящихся траекторий........ 2.5.2 Математическое описание супертреков.......... 2.5.2.1 Дерево супертрека................... 2.5.2.2 Кодирование структуры супертреков......... 2.5.2.3 Нумерация вершин супертреков........... 2.5.2.4 Математическое описание супертреков........ 2.5.3 Плотность вероятности................... 2.5.3.1 Физические вероятности................ 2.5.3.2 Плотность вероятности реализации супертрека... 2.5.4 Построение вероятностного пространства......... 2.5.4.1 Обзор литературы................... 2.5.4.2 Непосредственное построение вероятностного пространства........................ 2.5.4.3 Вероятностное пространство (S,, )...... 2.7.2 Разыгрывание со смещённой плотностью вероятности.. 2.7.5.2 Обоснование несмещённости оценки * при использовании метода DXTRAN в одной вершине..... 2.7.5.3 Обоснование несмещённости оценки * при использовании метода DXTRAN в нескольких вершинах. 2.8 Случай функционала, не аддитивного по траекториям... 2.10 Способ вычисления весовой оценки для схемы совпадений. 2.10.1 Представление отлика схемы совпадений в виде суммы. 2.10.3 Вычисление среднего значения методом Монте-Карло.. 3.2 Вычисление отклика детектора (1, 2,..., ) по данным 3.2.1 Условия, при которых возможно восстановление супертреков............................. 3.2.4 Восстановление необходимых для вычисления значения случайной величины данных о супертреке в случае 3.2.5 Вспомогательные поверхности для сбора информации о 4.5 Реализация алгоритма для вычисления неаддитивных 5.4 Вычисление оставленной в объёме детектора энергии с 5.6 Оптимизация геометрической конфигурации комптоновского гамма-спектрометра..................... 5.6.2 Оптимизация геометрической конфигурации....... А Построение вероятностного пространства........... Б.3 Перемена порядка интегрирования в повторном интеграле Г Обоснование метода Монте-Карло............... Г.1 Вероятностное пространство и случайная величина... Г.2 Классический метод Монте-Карло............. Список условных обозначений • x – пространственные координаты, задающие положение частицы;

• – единичный вектор направления скорости;

• – энергия частицы;

• – момент времени;

• = (x,,, ) – фазовые координаты частицы;

• – количество статистических испытаний;

• – максимально возможное количество частиц, получаемое в результате столкновения;

• – количество аргументов у функционала (1, 2,..., );

• – количество взаимоисключающих реализаций супертрека ;

• S – множество всех супертреков;

• S – подмножество множества всех супертреков S, в которое входят все супертреки с одинаковой структурой;

• – количество вершин у супертреков из множества S без учёта точки рождения первичной частицы;

• – время, в течение которого интегрируется сигнал ФЭУ;

• – мёртвое время сцинтилляционного детектора;

• – разрешающее время схемы совпадений;

• – аналоговая оценка среднего значения функционала;

• * – весовая оценка среднего значения функционала;

• () – отклик детектора на соударение в точке ;

• () – отклик детектора на супертрек ;

• () – вероятность поглощения в точке ;

• ( 1, 2,..., ) – плотность вероятности того, что после столкновения в точке, из реакции выйдет частиц, первая из которых испытает следующее столкновение в 1, вторая – в 2,..., • () – вероятность того, что при столкновении в точке из реакции выйдет частиц;

• () – плотность вероятности реализации супертрека ;

• ( |0, 1,..., 1 ) – условная плотность вероятности того, что при разыгрывании -ой вершины супертрека из S она будет иметь координаты, при условии, что все предыдущие вершины имели координаты 0, 1,..., 1 ;

• – статистический вес;

• + – события, соответствующие вхождению частиц в область детектора;

• – события, соответствующие выхождению частиц из области детектора;

• – энергия, оставленная в детекторе при единичном акте столкновения;

• ( ) – длина оптического пути от точки до точки ;

• I() – индикатор события.

Введение Актуальность В настоящее время интерес к разработке приборов, регистрирующих ядерное излучение, постоянно возрастает. Областями человеческой деятельности, в которых используются схемы совпадений, являются гаммаастрономия, обеспечение безопасности [21, 22, 42, 47] (досмотр грузов и контейнеров, экологическая безопасность), медицина [35, 44, 63, 71, 79] и некоторые другие. Успехи в области создания новых сцинтилляционных материалов позволяют улучшить такие характеристики приборов, работающих по схемам совпадений или антисовпадений, как отношение сигнал-шум и пространственное разрешение[43].

В сфере обеспечения безопасности важным направлением является разработка технических средств, позволяющих производить досмотр грузов, не вскрывая контейнеры. Наиболее перспективным в смысле высокой информативности и компактности является метод меченых нейтронов (ММН, в зарубежной терминологии – API, associated particle imaging). Эксперименты в лабораторных условиях показали возможность создания компактной системы измерения спектров гаммаизлучения от неупругого взаимодействия нейтронов с высоким пространственным разрешением на базе портативных нейтронных (, ) генераторов со встроенным многоэлементным детектором альфа-частиц.

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

В настоящее время в ряде отечественных и зарубежных лабораторий разработаны или планируются к разработке макетные образцы систем обнаружения и идентификации ВВ. Так в Радиевом институте им. В.Г.

Хлопина (г. С.-Петербург) по заказу дирекции программы НАТО «Наука для мира» для демонстрации возможностей метода меченых нейтронов создана серия приборов-интроскопов под общим названием SENNA.

Внешний вид приборов SENNA иллюстрирует рис. 1.

Вес прибора около 40 кг, однако все результаты по обнаружению получены с применением дополнительной коллимации детекторов железными конусами общим весом 10 кг. Для обнаружения ВВ весом около 400—500 грамм требуется время 5—10 минут, а весом 3 кг около минуты.

Дальность управления 10—15 метров ограничена длиной сигнального кабеля; питание от сети 220 В. Сокращение предела обнаружения (вес ВВ) и времени обнаружения невозможно из-за перегрузки каналов оцифровки сигналов, хотя при замене нейтронного генератора можно было бы повысить мощность нейтронного потока и, соответственно, улучшить статистику отсчетов. Жёсткая компоновка общей конструкции, фиксирующая положение генератора, модуля гамма-детектора и защиты и ограничивает подход к обследуемому объекту. В НПЦ «Аспект» г. Дубна разработана возимая установка для досмотра автомобильной техники на основе метода меченых нейтронов. Вес установки (без автомобиля) кг; порог обнаружения ВВ 10 кг за 10 минут. Внешний вид установки представлен на рис. 2.

Проект Евросоюза, названный EURITRACK, предназначается для обследования большегрузных контейнеров и находится в настоящее время в процессе проектирования. Предполагаемый объем – несколько кубометров. Внешний вид проектируемого терминального комплекса EURITRACK представлен на рис. 3.

Рис. 2. Схематичный вид прибора Aspect.

Рис. 3. Схематичный вид установки EURITRACK.

Комптоновский гамма-спектрометр позволяет определить направление прилёта гамма-квантов по измерениям потерянной ими энергии при последовательных рассеяниях. Ярким примером такого прибора является гамма-телескоп COMPTEL, работающий по схеме двойных рассеяний.

Комптоновский гамма-спектрометр состоит из нескольких слоёв детекторов, параметры которых подбираются таким образом, чтобы во всех слоях кроме последнего гамма-кванты испытывали преимущественно одиночное комптоновское рассеяние, а в последнем – полное поглощение.

Зная энергии, потерянные гамма-квантом в каждом слое и примерные координаты точек столкновения, можно определить конус возможных направлений прилёта гамма-кванта. Приборы на основе комптоновских гамма-спектрометров позволяют получить трёхмерную картину распределения источников гамма-излучения в некотором объёме, что делает его применение интересным в медицинской томографии [75] и в сфере безопасности [66].

Разработка новых приборов, как правило, включает в себя этап расчётно-теоретического обоснования их параметров и характеристик.

Для этого наиболее часто используется метод Монте-Карло, в котором проводится большое количество статистических испытаний для определения среднего значения некоторого функционала, заданного на множестве случайных траекторий.

Если при моделировании траекторий частиц используется физическая плотность вероятности, то такой метод Монте-Карло называется аналоговым. В некоторых случаях для достижения приемлемой статистической погрешности расчётов, проводимых аналоговыми методами Монте-Карло, требуется огромное количество испытаний. Например, для набора базовых спектров, позволяющих идентифицировать вещество при использовании ММН, в работе [33] потребовалось по 4, 85 · 1010 статистических испытаний для каждого вещества. Это послужило причиной разработки весовых методов Монте-Карло, в которых статистические испытания проводятся с использованием модельной вероятности, отличной от физической. Модельная вероятность выбирается так, чтобы увеличить количество попаданий частиц в область детекторов. Как при аналоговом, так и при весовом моделировании, зависимость относительной статистической погрешности от количества испытаний имеет вид Для весовых методов коэффициент пропорциональности существенно меньше, чем для аналоговых. За счёт этого при использовании весовых методов требуемый уровень статистической погрешности достигается при меньшем количестве испытаний, чем при использовании аналоговых методов.

В настоящее время весовые методы хорошо проработаны для вычисления значения функционалов определённого класса, названных в работе [37] больцмановскими. Под больцмановским функционалом понимается функционал, который может быть вычислен при известной одночастичной плотности распределения, являющейся решением уравнения Больцмана. Примером больцмановских функционалов являются среднее количество пересечений частицами некоторой поверхности или среднее количество столкновений внутри некоторого объёма. Если физическая величина зависит от совместного влияния нескольких частиц, то она не может быть представлена в виде больцмановского функционала[37].

Например, на отклик приборов, работающих по схемам совпадений, оказывают совместное влияние частицы, находящиеся в разных областях пространства. Функционал, описывающий отклик таких приборов, должен зависеть от совместной плотности распределения частиц. Большинство имитационных моделей приборов, регистрирующих ядерное излучение, не удаётся адекватно описать функционалом от плотности столкновений.

Современной тенденцией является построение сложных имитационных моделей приборов, учитывающих различные физические явления.

Например, в работах [36, 43, 82, 86] используется модель сцинтилляционного детектора, имитирующая оптическое распространение фотонов внутри сцинтилляционного вещества. Испускание сцинтилляционных фотонов в точках столкновения гамма-квантов с элементами вещества ведётся по закону Пуассона с переменной интенсивностью [77, 87], что ещё больше усложняет модель. Для учёта более тонких явлений таких, как инертность ФЭУ, производится моделирование испускания электронов фотокатодом и прохождения их между электродами ФЭУ.

В работе [76] приведена аналитическая модель, имитирующая поведение ФЭУ. Разработка универсальной программы, моделирующей движение сцинтилляционных фотонов, работу ФЭУ и учитывающей неидеальность электронной схемы, регистрирующей совпадения, проведена в [88].

В перечисленных работах модели детекторов не могут быть описаны больцмановскими функционалами, моделирование ведётся аналоговыми методами, а учёт различных физических явлений существенно снижает скорость расчётов. В программах MCNP-DSP [90] и MVP [70], предназначенных для решения широкого круга задач, и в программе KENONR [50], в случаях, когда важно учитывать совместное распределение частиц, также используются аналоговые методы расчёта траекторий частиц.

Проблема ускорения расчётов, проводимых по методу Монте-Карло, стоит очень остро и во многих других работах.

Современное состояние проблемы В литературе [3, 7, 11, 15, 16, 24, 25] весовые методы предлагаются для вычисления аддитивных по столкновениям функционалов, определённых на множестве реализаций случайного процесса переноса излучения. Значение такого функционала () равно сумме вкладов от отдельных столкновений, то есть где — случайная траектория частицы, — точки соударения частицы, движущейся по этой траектории, = 1, 2,.... Примером такого функционала является количество столкновений в выделенном объёме.

Очевидно, что аддитивные по столкновениям функционалы являются частным случаем больцмановских функционалов.

Ещё одним важным частным случаем больцмановских функционалов являются функционалы, представимые в виде Примером такого функционала является количество пересечений частицами выбранной поверхности.

Наиболее известными весовыми методами, применяемыми при вычислении среднего значения аддитивных по столкновениям функционалов, являются разыгрывание со смещённой плотностью вероятности, расщепление, русская рулетка и DXTRAN. При разыгрывании со смещённой плотностью вероятности искусственно повышается вероятность попадания частиц в область сосредоточения детекторов. При использовании метода расщепления частица расщепляется на несколько новых частиц, что даёт возможность точнее оценить вклад в показания детекторов исходной частицы. Метод русской рулетки позволяет прервать моделирование траектории частицы, если известно, что ожидаемый вклад частицы в показания детектора слишком мал. В методе DXTRAN при каждом столкновении частицы происходит оценка потока излучения внутрь выбранной сферы, окружающей область сосредоточения детекторов [28].

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

В большей части литературы построение универсальных методов для вычисления небольцмановских функционалов ограничивается аналоговой схемой, а весовые методы строятся для различных частных случаев. В классических работах С. М. Ермакова и Л. Яноши рассмотрены весовые методы вычисления первых двух моментов от аддитивных по столкновениям функционалов (дисперсия аддитивного функционала вычисляется как среднее квадрата, то есть неаддитивного функционала).

Также в работах С. М. Ермакова особое внимание уделяется случаю ветвящихся траекторий, так как он всегда представляет отдельную сложность при проведении доказательств несмещённости весовых оценок.

В работах В. В. Учайкина и А. В. Лаппы рассмотрены функционалы «столкновительно-трекового» класса, обобщающего класс аддитивных по столкновениям функционалов. В работах А. В. Лаппы построены неимитационные методы, вычисляющие любой момент от аддитивного по траекториям функционала. В [1, 41] построение весового алгоритма вычисления небольцмановских функционалов ведётся на конкретном примере вычисления энергии, оставленной частицами в детекторе, а среда считается неразмножающей. В [25] рассмотрен способ вычисления произвольного функционала при помощи поливариантного разложения, но на практике этот способ применим, когда все члены разложения, начиная с некоторого номера, близки к нулю. Также известны способы для вычисления любого момента аддитивного функционала [12, 13].

Основополагающими работами, предлагающими универсальный подход к вычислению небольцмановских функционалов весовыми методами, являются работы [37, 38]. В них предлагается приписывать статистический вес всей ветвящейся траектории целиком, а не каждой частице в отдельности, как это делается в схеме Неймана-Улама при вычислении среднего аддитивных по столкновениям функционалов.

Ветвящаяся траектория в этих работах названа супертреком, так как она рассматривается как неделимая коллекция отдельных неветвящихся траекторий частиц. В перечисленных работах предлагаются два способа получения супертреков.

В первом способе, названном способом развёртывания (deconvolution approach), результатом моделирования истории одной первичной частицы является ветвящаяся траектория, содержащая вершины двух типов. В вершинах первого типа происходит разветвление траектории естесвенным образом, то есть в результате моделирования реакции, происходящей при столкновении частицы. В вершинах второго типа разветвление происходит в результате применения какого-либо весового метода. Достигнув такой вершины, частица пошла бы дальше только по одной из ветвей. Выбирая в каждой вершине второго типа по одной ветви, получаем ветвящуюся траекторию, состоящую только из вершин первого типа, что и должно быть для физически осуществимой траектории. Вес этой траектории равен произведению весовых множителей выбранных ветвей. В работах [37, 38] описанный подход применяется для вычисления потерянной в чувствительном объёме детектора энергии частиц. В работах [83, 84] описанный подход предложен для вычисления отклика приборов, работающих по схемам совпадений. Стоит отметить, что в работах [83, 84] рассматриваются только совпадения, вызванные частицами, принадлежащими одной и той же ветвящейся траектории.

Совпадения, вызванные частицами из разных ветвящихся траекторий, не учитываются.

Во втором способе, названным способом разыгрывания супертреков (supertrack approach) результатом статистического испытания является вся ветвящаяся траектория целиком. Весовые методы применяются не к отдельной частице, а ко всей траектории целиком. Например, при применении метода расщепления в некоторой вершине из одной траектории получаются несколько новых, совпадающих друг с другом до этой вершины.

Многие весовые методы, разработанные для вычисления аддитивных по столкновениям функционалов, очень легко переносятся на случай развёртывания или разыгрывания супертреков. В работах [37, 38] в рамках концепции супертреков сформулированы весовые методы разыгрывания со смещённой плотностью вероятности, расщепления, русской рулетки и DXTRAN.

Метод развёртывания нашел воплощение во многих работах, например, в программах MCNP5 [40, 68, 81, 83], MCNPX [59] и MCBEND [78]. Метод разыгрывания супертреков кажется простым, но он требует разработки нового или глубокой переработки существующего программного обеспечения для его реализации. В работах [39, 48] выполнена модификация программы MCNP для поддержки концепции супертреков.

В работе [84] используется метод, очень похожий на метод разыгрывания супертреков.

В перечисленных работах, использующих концепцию супертреков, ведётся вычисление значений небольцмановских функционалов определённого класса, в которых частицы, относящиеся к разным ветвящимся траекториям, делают независимые вклады в значение функционала.

Другими словами, отклик детектора на траектории 1, 2,..., можно представить в виде Такие функционалы называются аддитивными по траекториям. Не все функционалы обладают этим свойством, например, иногда важно учитывать совпадения, вызванные фоновым излучением или отражёнными от удаленных предметов частицами. Такие совпадения могут быть учтены в математических моделях детекторов, не делающих различия между частицами, принадлежащими разным супертрекам, и частицами, принадлежащими одному и тому же супертреку.

В настоящей работе предложенная в [38] концепция супертреков обоснована и обобщена на случай вычисления неаддитивных по траекториям функционалов. Математическое обоснование несмещённости оценки среднего значения функционала ведётся традиционным способом через усреднение по всем возможным ветвящимся траекториям, в то время как в работе [37] акцент делается на обоснование предложенных способов через разбор большого количества примеров.

Цель Целью настоящей работы является разработка весовых методов Монте-Карло для вычисления отклика приборов, описываемых неаддитивными функционалами.

Методы исследований Методами исследований, использованными в настоящей работе, являются методы вычислительной математики, математического анализа, теории вероятностей и математической статистики, теории весовых методов Монте-Карло и теории переноса излучения. Программный код выполнен на языках Fortran-90, C++ и bash. Для расчётов использовался многопроцессорный вычислительный кластер с поддержкой HP-MPI и многопроцессорная рабочая станция.

Научная новизна 1. Концепция супертреков обобщена на случай вычисления функционалов, учитывающих совпадения между частицами, принадлежащими разным ветвящимся траекториям.

2. Впервые получено доказательство несмещённости весовых оценок в рамках концепции супертреков посредством их осреднения на множестве всех ветвящихся траекторий.

3. Построено представление функций отклика физических приборов, работающих по схемам совпадений, в виде неаддитивных по траекториям функционалов. Введено новое свойство аддитивности по времени, которым обладают функционалы (1, 2,..., ), описывающие работу физических приборов. Впервые доказано, что для оценки среднего значения таких функционалов достаточно единственной реализации кортежа (1, 2,..., ) из ветвящихся траекторий.

4. Впервые для учёта случайных совпадений применены весовые методы, разработанные автором и реализованные в виде комплекса программ.

Практическая значимость работы Основной практической ценностью полученных результатов является многократное ускорение расчётов в прикладных задачах радиационной физики.

Разработан универсальный программный комплекс, позволяющий оптимизировать геометрические параметры приборов, работающих по схемам совпадений.

Показано, что наиболее эффективной геометрической конфигурацией комптоновского гамма-спектрометра с точки зрения отношения количества истинных совпадений к полному количеству зарегистрированных совпадений является конфигурация при меньших расстояниях между слоями детекторов.

Проведен расчёт базовых спектров углерода, азота и кислорода в методе меченых нейтронов.

Разработанные программы использовались для теоретической оценки количества совпадений в работах [9, 62, 69, 72].

В ходе исследований эмпирически проверены статистические свойства некоторых популярных в настоящее время программных генераторов псевдослучайных чисел при помощи теста на равномерное заполнение единичного гиперкуба [30]. Тесты позволили обнаружить ошибку в реализации одного из генераторов библиотеки CLHEP. Для остальных генераторов найдены участки генерируемых ими последовательностей с хорошими и плохими статистическими свойствами.

Достоверность Достоверность полученных теоретических результатов обеспечивается их проверкой при помощи численного решения тестовых задач.

Многократные сравнения результатов весового моделирования с результатами аналогового моделирования подтверждают несмещённость предлагаемой в настоящей работе оценки среднего значения функционала.

Проведены проверки разработанных программных модулей на задачах с известными теоретическими решениями.

Разработанная модель сцинтилляционного детектора получена обобщением известной и хорошо себя зарекомендовавшей в работах многих авторов модели детектора на случай неаддитивных по траекториям функционалов. Результаты оптимизации геометрической конфигурации комптоновского гамма-спектрометра имеют наглядный физический смысл. Для модели интроскопа, работающего по методу меченых нейтронов, проведено сравнение расчётного спектра углерода с экспериментальным.

Положения, выносимые на защиту 1. Обобщение концепции супертреков на случай неаддитивных функционалов.

2. Обоснование концепции супертреков через доказательство несмещённости весовой оценки посредством усреднения по множеству ветвящихся траекторий.

3. Представление откликов приборов, работающих по схемам совпадений, в виде неаддитивных функционалов. В число приборов входит комптоновский гамма-спектрометр и интроскоп, работающий по методу меченых нейтронов.

4. Математическая модель комптоновского гамма-спектрометра, различающая истинные и ложные совпадения разных типов.

5. Комплекс программ, позволяющий вычислять среднее значение неаддитивных по траекториям функционалов.

Личный вклад автора Личный вклад соискателя заключается в обосновании и обобщении концепции супертреков, разработке способов применения весовых методов Монте-Карло для вычисления неаддитивных функционалов, разработке математических моделей комптоновского гамма-спектрометра и интроскопа, работающего по методу меченых нейтронов, и их программной реализации, тестировании разработанного программного обеспечения. Все результаты, представленные в диссертации, получены соискателем самостоятельно или при непосредственном его участии.

Апробация работы Результаты, изложенные в диссертации, докладывались и обсуждались на следующих конференциях и семинарах • 7-ая международная конференция «Люминесцентные детекторы и преобразователи ионизирующего излучения» LUMDETR- (Краков, Польша, 2009 г.);

• Пятая Всероссийская конференция «Проблемы обеспечения взрывобезопасности и противодействия терроризму» (Санкт– Петербург, 2010 г.);

• семинар кафедры компьютерного моделирования ФАЛТ МФТИ;

• семинар 11 отдела Института прикладной математики им. Келдыша «Вычислительные методы и математическое моделирование»;

• семинар кафедры статистического моделирования математикомеханического факультета СПбГУ;

• научные конференции МФТИ «Современные проблемы фундаментальных и прикладных наук», в секциях кафедр автоматизированных биотехнических систем и высшей математики;

Разработанные программы были использованы для оценки количества случайных совпадений в работах [9, 62, 69, 72].

По теме проводимых исследований соискателем опубликовано работ, из них 3 — в изданиях из списка ВАК. Семь работ выполнены без соавторов.

Получено свидетельство о государственной регистрации программы для ЭВМ № 2012613937 «Программный комплекс для вычисления неаддитивных функционалов весовыми методами Монте-Карло».

Структура работы Первая глава работы посвящена построению математических моделей приборов, работающих по схемам совпадений. Отклик прибора описывается функционалом (1, 2,..., ), где 1, 2,..., – реализации ветвящихся траекторий. Показано, что в зависимости от того, какие физические эффекты необходимо учитывать, функционал может быть аддитивным по столкновениям, по траекториям или не быть аддитивным вообще. Модели приборов, работающих по схемам совпадений, описываются неаддитивными функционалами.

Введено свойство аддитивности по времени. Пусть траектории разыгрываются в порядке увеличения момента времени рождения первичной частицы. Пусть в интервале времени (1, 2 ) в системе не существует частиц. Пусть траектории, существующие до момента времени 1, имеют номера от 1 до 1 1, а траектории, существующие после 2 - от 1 до. Пусть 2 1 >, где – некоторый параметр. Тогда, если то будем говорить, что функционал (1, 2,..., ) удовлетворяет свойству аддитивности по времени с параметром.

Проведя разбиение на слагаемые по всем таким интервалам времени (1, 2 ), получаем, что (1, 2,..., ) равен сумме случайного количества слагаемых. Доказано, что все слагаемые, быть может, кроме последнего независимы в совокупности и одинаково распределены.

Выполнен обзор статей, посвящённых моделированию сцинтилляционных детекторов и приборов на их основе. Обоснован выбор простейшей модели сцинтилляционного детектора, в которой за отклик детектора принимается оставленная частицами в чувствительном объёме энергия.

Описаны математические модели электронного дифференциального спектрометра и электронной схемы совпадений, обрабатывающих сигналы сцинтилляционных детекторов.

На основе описанных моделей разработана модель интроскопа, работающего по методу меченых нейтронов, и комптоновского гаммаспектрометра. Модель интроскопа работает в приближении аддитивности по траекториям. Модель комптоновского гамма-спектрометра проводит раздельных подсчёт истинных событий и ложных событий разного типа. Предложен критерий для оптимизации геометрической конфигурации комптоновского гамма-спектрометра.

Модели приведены к одному формату входных данных и могут быть описаны функционалом (1, 2,..., ). Это позволит абстрагироваться от особенностей конкретной модели и построить общие универсальные весовые методы.

Вторая глава посвящена построению универсальных численных методов для вычисления среднего значения неаддитивных функционалов (1, 2,..., ), с учётом ветвлений траекторий. Обзор литературы, проведённый в этой главе, показал, что наиболее универсальным имитационным способом вычисления неаддитивных функционалов является концепция супертреков.

Выписана предлагаемая в настоящей работе оценка среднего значения неаддитивного функционала. Обоснование несмещённости предлагаемой оценки проведено через доказательство равенства математического ожидания весовой оценки математическому ожиданию величины, полученной аналоговым моделированием. Вычисление математического ожидания в обоих случаях проведено путём усреднения по всем возможным ветвящимся траекториям.

Третья глава посвящена техническим аспектам построения вычислительных моделей приборов, работающих по схемам совпадений. В этой главе обоснован выбор кода программы MCNP в качестве транспортного кода. Особое внимание уделено подключению кода MCNP и сбору необходимой для развёртывания супертреков информации.

Показано, что задача восстановления необходимой информации о ветвящихся траекториях по данным кода MCNP может быть решена сравнительно просто в двух случаях – в схемах совпадений, и в задачах, в которых несущественна корреляция между супертреками. В соответствии с описанным в третьей главе способом подключения кода MCNP разработано программное обеспечение для генерации супертреков, которые являются входными данными подпрограмм, разработанных во второй главе.

Четвертая глава посвящена тестированию разработанного программного обеспечения. Тестирование проведено сначала для каждого модуля по отдельности, а затем для всего программного обеспечения в целом. Эмпирическая проверка несмещённости предлагаемой в настоящей работе оценки случайной величины проведена на задаче о бросании нескольких точек на единичный отрезок. Проверка правильности подключения кода MCNP проводилась путём сравнения получаемых от этого кода данных с данными трассировки траекторий частиц. Проверка правильности работы функций, имитирующих функционирование электронной схемы совпадений, проводилась на простой задаче, для которой известна аналитическая оценка количества совпадений в зависимости от количества испущенных источником частиц. Для комптоновского гамма-спектрометра предложен критерий, при помощи которого можно отделить смоделированные истинные совпадения от ложных. Показано, как предложенный критерий может быть применен для оптимизации геометрической конфигурации прибора, работающего по схеме двойных рассеяний. Расчёты показали, что использование предлагаемых весовых методов существенно ускоряет скорость расчётов по сравнению с аналоговыми методами.

Заключение содержит основные результаты работы.

Приложение А посвящено построению вероятностного пространства на множестве всех ветвящихся траекторий.

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

Приложение В содержит описание способа подключения транспортного кода программы MCNP.

Приложение Г для справки содержит строгое математическое обоснование метода Монте-Карло.

Приложение ?? посвящено расчёту базовых спектров углерода, азота и кислорода при помощи разработанной модели интроскопа, работающего по методу меченых нейтронов.

Глава Математические модели многодетекторных приборов 1.1 Свойство аддитивности 1.1.1 Аддитивность по столкновениям Пусть траектория частицы состоит из точек фазового пространства 0, 1,...,, соответствующих рождению и последующим столкновениям частицы. Тогда для функционала (), описывающего отклик аддитивного по столкновениям детектора, на который влияют только столкновения частиц, справедливо представление Именно для оценки среднего значения функционалов, обладающих таким свойством аддитивности, разработано достаточное количество весовых методов Монте-Карло. Такими функционалами можно представить различные характеристики поля излучения, например, количество столкновений частиц в заданном элементе объёма. Как правило, результатом расчётов являются средние значения этих величин, отнесённые к количеству частиц, испущенных источником. Для устранения возможной путаницы будем называть рассмотренное свойство аддитивности (1.1) свойством аддитивности по столкновениям.

Свойство аддитивности по столкновениям иногда удобно обобщить следующим образом что даёт возможность отнести к аддитивным функционалам ток или поток частиц через заданную поверхность. Например, для вычисления тока частиц через поверхность величина ( +1 ) принимается равной 1, если участок траектории от до +1 пересекает поверхность, иначе она равна 0.

Ещё более общий класс аддитивных функционалов рассматривается в работе [14], в которой вводятся функционалы, представимые в виде суммы столкновительной части, трековой, мгновенной и поверхностой.

Столкновительная часть представляется в виде суммы вкладов от всех столкновений, трековая - в виде суммы от всех прямолинейных участков, мгновенная часть определяется положением всех частиц в заранее определённые моменты времени, а поверхностая - суммой по всем пересечениям частицами выбранной поверхности.

1.1.2 Аддитивность по траекториям Если моделируется работа спектрометра на базе сцинтилляционного детектора, то за отклик такого прибора обычно принимается количество попаданий в выбранный канал спектрометра. Пусть — реализация ветвящейся траектории. Функционал, описывающий количество попаданий в -ый канал спектрометра, имеет вид где () — суммарная энергия, потерянная частицей в чувствительном объёме детектора, () — функция, равная 1, если энергия попадает в интервал энергий -ого канала спектрометра, и равная 0, если не попадает.

Рассмотрим пример (см. рис. 4). Пусть в детектор попадает некоторая частица, которая в точках 1 и 2 теряет энергии 1 и 2. Обозначим номера каналов спектрометра, в которые попадают эти энергии, через 1 и 2. Если бы произошло столкновение только в точке, то отклик детектора ( ) на столкновение был бы равен 1, если номер рассматриваемого канала совпадает с номером, иначе он равен 0.

То есть ( ) = I( = ), где I() — индикатор события. Номер канала спектрометра, в который попадает сумма энергий = 1 + обозначим через 3. Отклик детектора на оба столкновения равен 1, если = 3, иначе он равен 0, то есть (1, 2 ) = I( = 3 ). Очевидно, что (1, 2 ) = (1 ) + (2 ), то есть рассмотренный функционал уже не является аддитивным по столкновениям.

Рис. 4. Вычисление оставленной в детекторе энергии.

Отклик детектора на траекторий 1, 2,..., по-прежнему можно представить в виде суммы то есть отклики на каждую траекторию расчитываются отдельно друг от друга. Функционалы, обладающие свойством (1.2), будем называть аддитивными по траекториям.

Для вычисления средних значений неаддитивных по столкновениям, но аддитивных по траекториям функционалов весовые методы разработаны недостаточно. В виде таких функционалов можно представить количество обнаруженных совпадений, если достаточно учитывать совпадения, вызванные одной и той же частицей. Так сделано в работах [57], [74], [80], посвященных моделированию схем совпадений.

1.1.3 Неаддитивные функционалы Для учёта совпадений, вызванных частицами, относящимися к разным ветвящимся траекториям, функционалы вида (1.2) уже не годятся.

В этом случае функционал, описывающий количество зарегистрированных прибором совпадений, можно записать в самом общем виде где 1, 2,..., — реализации случайных траекторий, — количество частиц, испущенных источником, при этом ни одно из рассмотренных выше свойств аддитивности не предполагается.

Приведём примеры неаддитивных функционалов. Пусть в детектор сначала попадает первая частица(см. рис. 5), траектория которой обозначена через 1. В точке 1 частица теряет энергию 1. Пусть этой энергии соответствует канал спектрометра 1, то есть (1 ) = I( = 1 ). Для второй частицы, траекторию которой обозначим через 2, потерянную в детекторе в точке 2 энергию обозначим через 2. Пусть она попадает в канал спектрометра с номером 2, (2 ) = I( = 2 ). Сумма откликов детектора на эти две траектории (1 ) + (2 ) равна 2, если номер рассматриваемого канала спектрометра совпадает с номерами 1 и одновременно ( = 1 = 2 ), равна 1, если номер совпадает либо с 1, либо с 2, иначе она равна нулю.

Рис. 5. Вычисление энергии сцинтилляционной вспышки в случае одновременного попадания двух частиц.

В случае, когда в детектор попадают две эти частицы одновременно, энергия сцинтилляционной вспышки равна 1 + 2. Номер канала спектрометра, в который попадает эта энергия, обозначим через 3, (1, 2 ) = I( = 3 ). Очевидно, что (1, 2 ) = (1 ) + (2 ).

Рассмотрим схему совпадений (cм. рис. 6) в качестве второго примера. Прибор срабатывает тогда, когда срабатывают оба детектора и 2 одновременно. Пусть в результате соударения частицы в точке образовалось две вторичных частицы. Одна из них испытала соударение в точке 1 внутри детектора 1, вторая — в 2 внутри детектора 2. В этом случае регистрируется совпадение, то есть отклик прибора (1, 2 ) = 1. Если соударение произошло только в одном детекторе, то совпадение не регистрируется, то есть (1 ) = (2 ) = 0. Следовательно, (1, 2 ) = (1 ) + (2 ), то есть отклик прибора в этой схеме совпадений не описывается аддитивным по столкновениям функционалом.

Рис. 6. Регистрация совпадения, вызванного одновременным Можно показать, что отклик прибора, работающего по схеме совпадений, не всегда описывается аддитивным по траекториям функционалом. На рис. 7 показан случай, когда вторая частица, родившаяся в точке не попадает в детектор 2, но в него случайно попадает частица, относящаяся к траектории 2 (в таком случае говорят о случайном совпадении). Очевидно, что (1, 2 ) = (1 ) + (2 ). Следовательно, в задачах, в которых важно учесть случайные совпадения, отклик прибора не может быть описан аддитивным по траекториям функционалом.

Аддитивный по столкновениям функционал является частным случаем аддитивного по траекториям функционала. Если функционал аддитивен по траекториям, то он будет частным случаем функционала вида (1.3), неаддитивного ни в одном смысле.

1.1.4 Аддитивность по времени 1.1.4.1 Определение Модели приборов, рассматриваемые в настоящей работе, не являются аддитивными ни по столкновениям, ни по траекториям. Однако модели физических приборов всё же могут быть описаны аддитивными в некотором смысле функционалами. Например, электронная схема, обрабатывающая сигнал ФЭУ в сцинтилляционном детекторе, при превышении сигналом некоторого порога, переключается в режим интегрирования сигнала. Спустя некоторое время полученное значение энергии передаётся спектрометру, который прибавляет единицу к счётчику числа попаданий в соответствующий интервал энергий. После этого прибор возвращается в состояние ожидания следующей частицы. Рассмотренный алгоритм работы наводит на мысль о том, что функционал (1, 2,..., ) в некоторых случаях всё же можно разбить на сумму слагаемых.

Будем предполагать, что супертреки 1, 2,..., упорядочены по моменту времени рождения частицы. Упорядочить их всегда можно, если (1, 2,..., ) обладает свойством симметрии, то есть не меняет своего значения при перестановке любых двух аргументов.

Обозначим через 0 интервал времени, необходимый прибору, чтобы вернуться в исходное состояние. Обозначим через, = 1, 2,..., 1, интервалы времени, в течение которых в системе не существует частиц вообще, причем будем рассматривать только интервалы длительнее 0, то есть Рассмотрим один такой интервал = (1, 2 ) с номером. Все траектории можно разбить на две группы. К первой группе относятся траектории, для которых все моменты времени существования частиц внутри поверхности меньше 1. Ко второй группе относятся траектории, для которых эти моменты времени больше 2. Пусть индексы траекторий, относящихся к первой группе, меняются от 1 до, а ко второй — от + до. Таким образом, построены номера, причём 1 1 < 2 <... < Если выполняется то будем говорить, что функция детектора обладает свойством аддитивности по времени с параметром 0.

Если интервал времени между испусканием частиц источником бесконечно большой, то моменты существования частиц из любых двух траекторий разделены бесконечным интервалом времени. Тогда в выражении (1.5) разбиение на слагаемые происходит при каждом значении индекса 1, и выражение (1.5) переходит в (1.2).

1.1.4.2 Свойства Пусть разыграна бесконечная неубывающая последовательность { } моментов времени рождения первичных частиц. Также имеется бесконечная последовательность супертреков { }, причём -ый супертрек порождается -ой первичной частицей. Пусть на множестве всех пар последовательностей { } и { } задан функционал (1, 2,... ; 1, 2,...) (здесь моменты времени рождения первичной частицы специально отделяются от описания супертрека).

Обозначим { } последовательность интервалов времени между рождениями частиц, = +1. Элементы последовательности распределены показательно с плотностью вероятности и независимы. Будем считать функционал (1, 2,... ; 1, 2,...) однородным по времени, то есть представимым в виде (1, 2,... ; 1, 2,...) = (1, 2,... ; 1, 2,...). Другими словами, показания детектора не зависят от выбора начала отсчёта времени, а также от момента включения прибора, лишь бы прибор был включен до рождения самой первой частицы.

Пусть после поглощения -ой частицы выполняется условие (1.4).

Рассмотрим две случайных величины 1 = (1, 2,... ; 1, 2,...) и 2 = (+1, +2,... ; +1, +2,...). Элементы хвоста {+ } последовательности { } распределены так же, как и элементы всей последовательности целиком, не зависимо от того, является ли случайной величиной или нет. То же самое можно сказать и про последовательности {+ } и { }.

Следовательно, случайные величины 1 и 2 одинаково распределены.

Для доказательства одинаковой распределённости слагаемых в формуле (1.5) осталось взять в качестве (1, 2,... ; 1, 2,...) значение первого слагаемого в этой формуле (отклик детектора на начальные члены последовательностей { } и { } до момента обрыва их условием (1.4)), а в качестве взять 1. При таком выборе получаем, что является первым слагаемым формулы (1.5), а 2 – вторым слагаемым.

Так как все члены последовательностей { } и { } имеют одинаковое распределение, то 1 и 2 то же имеют одинаковое распределение.

При этом первое слагаемое зависит от членов последовательностей с номерами от 1 до 1 1, а второе – от 1 до 2. Так как члены последовательностей { } и { } разыгрываются независимо друг от друга и независимо от номера, то независимыми будут первые два слагаемых в (1.5) 1.

Таким же способом легко показать, что любые первые слагаемых независимы в совокупности и одинаково распределены.

Возвращаясь к обозначениям в формуле (1.5) получаем, что все слагаемые в этой формуле, быть может, кроме последнего независимы в совокупности и имеют одинаковое распределение. Последнее (-ое слагаемое) может иметь другое распределение, так как величины 1, 2,..., 1 выбираются из условия (1.4), а величина = задаётся константой.

Требование однородности функционала = (1, 2,... ; 1, 2,...) можно заменить на = (, 1, 2,... ; 1, 2,...), где – интервал времени до рождения первой частицы с момента времени включения прибора (или момента времени, когда прибор вышел из мёртвого состояния и готов регистрировать новые события). Но в этом случае принципиально, что бы интервалы времени были распределены показательно, в то время как ранее это свойство не использовалось. Только для показательного распределения где () – вероятность того, что время ожидания рождения следующей частицы превысит, () = (свойство показательного распределения, известное, как парадокс ожидания). Интервал времени до рождения следующей частицы не зависит от того, в какой момент времени прибор включился, и тоже имеет показательное распределение с плотностью.

1 = 1 ({ }, { }), то есть зависит от всех членов обоих последовательностей, как и 1 = 1.2 Требования к разрабатываемым математическим моделям Как было уже показано, модель детектора может описываться неаддитивным функционалом или функционалом, обладающим одним из свойств аддитивности. Выбор, каким свойством аддитивности должен обладать функционал, зависит от того, события каких типов должны учитываться в разрабатываемой модели прибора. Наиболее общим случаем является неаддитивный функционал. Дальнейшие рассуждения строятся именно для случая неаддитивного функционала, если это не оговорено особо.

Разрабатываемые модели приборов являются имитационными моделями, так как имитационные методы Монте-Карло обладают большей универсальностью и потенциальной возможностью расширения их функциональности. При имитационном моделировании разыгрываются ветвящиеся траектории частиц, которые являются входными данными для моделей детекторов. Ветвящаяся траектория описывается деревом (теория графов) и фазовыми координатами частиц в вершинах дерева.

Фазовые координаты состоят из пространственных координат, направления скорости, энергии и времени. Никаких других данных о траекториях представленные модели не используют.

Вопрос о том, как разыгрывать ветвящиеся траектории частиц, в настоящей главе не рассматривается.

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

Дополнительным требованием является раздельный подсчёт истинных и различного типа ложных совпадений. Это требование наиболее актуально в модели комптоновского гамма-спектрометра.

Моделируемыми приборами являются спектрометр на основе сцинтилляционного детектора, комптоновский гамма-спектрометр, интроскоп, работающий по методу меченых нейтронов.

1.3 Моделирование сцинтилляционного детектора 1.3.1 Выбор модели сцинтилляционного детектора Моделированию сцинтилляционных детекторов посвящены многие работы, например, [43, 77, 86, 87]. В них под сцинтилляционным детектором понимается сборка, состоящая из чувствительного сцинтилляционного вещества, ФЭУ и электронной схемы, обрабатывающей сигналы ФЭУ. В разных моделях можно выделить части, отвечающие за моделирование одного или нескольких из перечисленных элементов. Учет влияния других элементов детектора на выходной сигнал, например, влияния усилителя сигнала ФЭУ, как правило не проводится. Перечисленные три элемента сцинтилляционного детектора отличаются друг от друга с точки зрения возможностей точного моделирования их работы. Наиболее сложными для моделирования и вносящими наибольшую погрешность являются процессы, проходящие в сцинтилляционном веществе, в то время как преобразование аналогового сигнала ФЭУ в цифровой и преобразование цифрового сигнала внутри электронной схемы может быть смоделировано сравнительно просто. Учет искажений, вносимых ФЭУ (например, учёт инерционности ФЭУ), в большинстве случаев не является необходимым и проводится не во всех работах.

В простейшем случае принимается, что отклик сцинтилляционного детектора пропорционален энергии, оставленной частицей в чувствительном объеме детектора. В частности, принимается, что отклик детектора не зависит от того, в какой части чувствительного объема произошла сцинтилляционная вспышка. Для детекторов протяженной формы справедливость этого предположения проверялась в работе [55].

В этой работе показано, что для протяженного стержня LSO длиной 3 см энергия, соответствующая середине фотопика, может меняться в два раза в зависимости от того, какая область стержня облучается.

Сечение стержня — квадрат со стороной 2—3 мм. Для кристаллов другой формы положение фотопика может меняться еще значительнее. Также известно, что энергия сцинтилляционной вспышки не всегда пропорциональна потерянной частицей энергии. В работе [73] показано, что для кристалла LSO пропорциональность нарушается при энергии падающих квантов в диапазоне от 5 до 60 кэВ. Несмотря на перечисленные факты, предположение о пропорциональности отклика потерянной энергии и независимости его от того, в какой части сцинтиллятора произошла вспышка, вполне приемлемо для большинства задач (в том числе и для LSO в рассматриваемом нами диапазоне энергий в несколько МэВ).

Расчет величины потерянной энергии прост в реализации и присутствует во многих программах, моделирующих траектории частиц методом Монте-Карло.

В настоящей работе за измеренную энергию сцинтилляционной вспышки принимается суммарная энергия, потерянная частицами в чувствительном объёме детектора.

Выбор момента пересечения частицей границы детектора за момент времени начала сцинтилляционной вспышки может быть обоснован тем, что время интегрирования (порядка 1 мкс) существенно больше характерного времени нахождения частицы внутри чувствительной области детектора (порядка / 1 10 нс, 10 см — размер детектора, — скорость частицы, сравнима со скоростью света) и больше характерного времени сцинтилляционной вспышки (до нескольких сотен наносекунд). В обоих случаях за время интегрирования будут учтены все взаимодействия частицы внутри сцинтилляционного вещества. Но таким выбором вносится систематическая погрешность в определение момента времени начала вспышки (порядка /), которая в настоящей работе не существенна.

Для более точного определения момента, когда фототок превысит пороговое значение, в работах [43, 88] моделируется форма сигнала ФЭУ.

В них производится моделирование распространения сцинтилляционных фотонов с учётом известных законов нарастания и спада интенсивности сцинтилляционной вспышки и с учётом качества обработки поверхностей кристаллов. В этих же работах для моделирования работы ФЭУ используется функция отклика ФЭУ на попадание фотона на фотокатод. Помимо этого в работе [43] учитывается качество обработки поверхностей кристаллов. Возможность моделирования сцинтилляционных фотонов присутствует в GEANT4 [54].

В работе [87] приводится формула, отражающая закон спада во времени интенсивности испускания фотонов в точке, где произошла сцинтилляционная вспышка. Зависимость с двумя экспонентами, отражающая нарастание и спад сигнала, приведена, например, в работе [60].

Более детальное описание обоих зависимостей приведено в [77]. Моделирование движения электронов внутри ФЭУ при расчете схем совпадений встречается редко. В работах [43] используется аналитическая формула, приведённая в [65], отражающая отклик ФЭУ на попадание на фотокатод фотона. Фундаментальные ограничения на точность определения момента попадания фотона на фотокатод, следующие из статистического распределения моментов появления фотоэлектронов, приведены в [76].

Используемая в настоящей работе модель детектора может быть с успехом заменена на одну из моделей, рассмотренных выше. Несмотря на то, что некоторые из этих моделей являются сугубо аналоговыми, они вписываются в общий случай неаддитивного функционала и могут быть сопряжены с разрабатываемыми весовыми методами.

1.3.2 Вычисление потерянной энергии Потерянную энергию можно вычислить несколькими разными способами. Один из них — суммирование потерь энергии частицей в каждой точке взаимодействия внутри чувствительного объёма детектора. Для этого в некоторых работах даже вводится дополнительная фазовая координата, равная потерянной частицей энергии [1]. Альтернативным методом, используемым в настоящей работе, является вычисление разности энергии частицы на входе в детектор и на выходе из него. Выбор в пользу этого метода был сделан после анализа возможностей реализации обоих способов с использованием кода программы MCNP.

Потерянная энергия может вычисляться двумя способами: как для отдельной частицы, так и для нескольких частиц, если они попадают в детектор одновременно.

В первом способе при вычислении отклика детектора учитываются только частицы, порожденные одной и той же первичной частицей, включая саму первичную частицу. Другими словами, для каждой ветвящейся траектории вычисляется отклик (), а отклик детектора на несколько траекторий равен сумме откликов на каждую траекторию в отдельности (выполняется аддитивность по траекториям). Этот способ соответствует счетчику F8 программы MCNP. На рис. 8, иллюстрирующем рассматриваемый способ, энергия сцинтилляционной вспышки Во втором способе учитывается, что во время интегрирования сигнала ФЭУ может появиться еще одна частица, относящаяся к другой траектории. Новая частица вызывает сцинтилляционные вспышки, накладывающиеся на вспышки от первой частицы. Для учета этого производится имитация работы электронного блока, интегрирующего сигнал ФЭУ. За момент времени начала сцинтилляционной вспышки принимается момент пересечения частицей границы детектора. От этого момента времени отсчитывается интервал времени, в течение которого при расчете потерянной энергии учитываются все частицы, попадающие внутрь детектора. На рис. 9 энергия сцинтилляционной вспышки попрежнему равна 1 2 3 + 4 5, но эта величина складывается из слагаемых, относящихся к разным траекториям 1 и 2.

Очевидно, что необходимым условием влияния частицы на отклик детектора является пересечение частицей поверхности детектора. СущеРис. 8. Вычисление энергии сцинтилляционной вспышки. Первый Рис. 9. Вычисление энергии сцинтилляционной вспышки. Второй ственными для вычисления отклика детектора параметрами траектории являются энергия частицы при входе в детектор, энергия частицы при выходе из детектора, моменты времени пересечения частицей границ детектора.

За входные данные всех рассматриваемых ниже моделей принимаются точки пересечения частицами поверхности детектора. В нестационарном случае в описание точек пересечения входит время.

Выходными данными является последовательность импульсов, снимаемых с интегрирующего сигнал ФЭУ блока электроники. Каждый импульс характеризуется энергией и временем прихода фронта. В рамках выбранной математической модели сцинтилляционного детектора энергия импульса равна энергии вспышки и равна энергии, потерянной частицами в чувствительном веществе детектора.

1.3.4 Формулировка используемого алгоритма 1.3.4.1 Выделение точек пересечения частицами границ детектора Сформулируем алгоритм, по которому в настоящей работе вычисляется потерянная энергия. Примем за отчёт детектора энергию частицы, оставленную внутри некоторой пространственной области (в этой области находится чувствительное вещество детектора). Обозначим поверхность, ограничивающую эту область, через.

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

Для каждой точки пересечения на основании указанных данных о траектории можно вычислить моменты времени, в которые частица, двигаясь вдоль ребра, попадает в каждую из точек пересечения. Множество всех точек пересечения прямолинейных участков траектории с поверхностью будем называть пересечением траектории с поверхностью. Пусть — момент времени пересечения частицей поверхности, — энергия частицы, — направления движения частицы, равное 1 при движении частицы внутрь области, и 1 — наружу. Совокупность = {,, } будем называть событием. Событие соответствует пересечению частицей границы детектора.

Множество всех событий распадается на два подмножества + и. В подмножество + входят те события, для которых = 1, в подмножество — события, для которых = 1. О любых двух событиях 1 + и 2, образованных пересечением одной и той же траектории с поверхностью, будем говорить, что событие следует за событием 1, если путь в дереве, описывающем ветвящуюся траекторию, от 1 до 2 полностью лежит внутри области. В случае неразмножающей среды физически это означает, что после входа в область в момент времени 1 частица вышла из нее в момент времени 2. В случае размножающей среды за одним событием 1 может следовать несколько событий 2, 3 и т. д., причём все дочерние частицы, родившиеся внутри, первый раз покидают область в моменты времени соответственно 2, 3 и т. д.

На рис. 10 приведён пример разбиения пересечений ветвящейся траектории с границей области на подмножества + и. В точке в результате столкновения родилось две вторичных частицы. Попаданию первой частицы в область соответствует событие +. Внутри области из этой частицы образовалось опять две вторичных частицы, при вылете которых из области происходят события и. Вторая частица, родившаяся в точке, попадает в область два раза. Первое посещение соответствует событиям + и, второе — + и. За событием 1 следуют события 1,1 и 1,2, за событием 2 —, за событием 3 — 3,1.

Рис. 10. Пересечение ветвящейся траектории и границы области.

Пусть + = {+, +, + } — множество всех событий, полученных пересечением разыгранных методом Монте-Карло траекторий с областью и соответствующих вхождению частиц в область, + = соответствующие вылету частиц из области, = 1, = 1, 2,...,.

Потерянная при -ом попадании в область энергия равна эта величина имеет смысл вклада в энергию сцинтилляционной вспышки, сделанного событием +. Определим теперь, какие события + внесут вклад в энергию сцинтилляционной вспышки. Возможны два способа отбора событий.

1.3.4.2 Случай = () В первом способе отклик детектора рассчитывается для каждой траектории в отдельности ((1, 2,..., ) аддитивен по траекториям).

При этом энергия вспышки, вызванная частицами траектории рассчитывается по формуле где суммирование ведётся по всем номерам, при которых события + соответствуют точкам вхождения в область частиц выбранной траектории. Физически этот случай означает, что время между испусканием частиц источником несоизмеримо больше характерной временной длительности вспышки.

Для удобства перечислим основные шаги алгоритма в случае = () вместе:

1. Для каждой траектории найти все пересечения с границей области и определить множества событий {+ } и { }.

2. Для каждого события + вычислить величину по формуле (1.6).

3. Для каждой траектории определить энергию сцинтилляционной вспышки по формуле (1.7).

Во втором способе производится выбор событий +, попадающих в окно интегрирования. Будем считать, что множество {+ } упорядочено по моменту времени пересечения области, то есть +1 +. Среди моментов времени + выберем первые моментов таким образом, чтобы где — длительность временного окна, в течение которого производится интегрирование сигнала ФЭУ. Если количество элементов во множестве { } мало и не позволяет удовлетворить второму неравенству, то за номер принимается последний номер в упорядоченном множестве { }. Все сцинтилляционные вспышки с энергиями 1, 2,..., попадают внутрь окна, в течение которого производится интегрирование сигнала ФЭУ. Измеренная детектором энергия равна причём за время открытия окна интегрирования принимается время прихода фронта первой вспышки, попавшей в окно времени интегрирования, которое равно 1. Определим теперь индекс таким образом, что где — мёртвое время детектора, в течение которого все вспышки игнорируются. Аналогично выбору, если второе неравенство не может быть удовлетворено, то в качестве номера + надо выбрать последний номер в упорядоченном множестве { }. Из множества {+ } уберем первые + событий, после чего алгоритм возвращается в исходное состояние, но с уже меньшим количеством событий. Повторяя перечисленные шаги до тех пор, пока множество {+ } не пусто, получаем последовательно энергии второй, третьей и т. д. сцинтилляционных вспышек 2, 3,..., и моменты времени их возникновения 2, 3,..., Для удобства перечислим основные шаги изложенного алгоритма в 1. Для каждой траектории найти все пересечения с границей области и определить множества событий {+ } и { }.

2. Для каждого события + вычислить величину по формуле (1.6).

3. Множество пар чисел {(, )} отсортировать по возрастанию времени пересечения границы области.

4. Определить последовательность измеренных детектором энергий (1.8), (1.9), (1.10) интервал 0, за который модель сцинтилляционного детектора успевает вернуться в исходное состояние, равен сумме +.

1.4 Моделирование спектрометра Объектом моделирования является многоканальный спектрометр, функциями которого является определение, в какой канал попадает зарегистрированный импульс, и увеличение счётчика этого канала.

За входные данные модели спектрометра принимается совокупность выходных данных моделей детекторов. Это последовательность измеренных детектором энергий и моменты времени для каждого детектора.

Пусть — энергия сцинтилляционной вспышки, пусть — энергетические границы каналов, = 0, 1,...,. Модель спектрометра определяет номер канала такой, что Затем количество попаданий в этот канал увеличивается на единицу.

Выходными данными модели является вектор, компоненты которого показывают количество попаданий в каждый из каналов спектрометра. Все выводимые в настоящей работе формулы доказываются для скалярных случайных величин, но они легко могут быть перенесены на случай векторных случайных величин.

Временной интервал 0, за который спектрометр успевает вернуться в исходное состояние (см. п. 1.1.4), может быть выбран равным временному интервалу 0 модели сцинтилляционного детектора. Неучитывание параметров электронной схемы оправдывается тем, что характерное время реакции электронной схемы существенно меньше времени интегрирования сигнала ФЭУ в сцинтилляционном детекторе. Последний параметр как раз и определяет 1.5 Моделирование электронной схемы совпадений Объектом моделирования является электронная схема, обрабатывающая сигналы ФЭУ, в составе комптоновского гамма-спектрометра.

Схемы совпадений характеризуются следующими основными параметрами: разрешающим временем (максимальный временной сдвиг между входными сигналами, при котором они регистрируются как одновременные, от англ. permission time), чувствительностью (минимальный уровень входных сигналов, поступающих одновременно на все входы схемы совпадений), мёртвым временем (от англ. paralysis time).

В настоящей работе учитываются только разрешающее время и мёртвое время схемы совпадений, в остальном схема совпадений считается идеальной. Мёртвое время обычно намного меньше разрешающего времени.

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

Модель электронной схемы совпадений имеет три состояния (см.

рис. 11): основное состояние, состояние ожидания сигнала от второго детектора и состояние отработки мёртвого времени. Выделим условия, при которых модель переходит из одного состояния в другое.

Рис. 11. Диаграмма состояний модели электронной схемы совпадений.

В начальный момент времени модель находится в основном состоянии и ждёт появления сигнала от одного из сцинтилляционных детекторов. При поступлении сигнала модель переходит в состояние ожидания появления сигнала от второго детектора. Здесь возможны три альтернативных пути развития событий. Первый путь реализуется, если пришёл сигнал от того же самого детектора. Тогда модель ведёт себя так же, как если бы она находилась в основном состоянии на момент прихода этого сигнала. Второй путь реализуется, если пришёл сигнал от другого детектора. В этом случае модель регистрирует совпадение и переходит в состояние отработки мёртвого времени. Третий путь — если в течение времени не пришёл ни один сигнал. Тогда модель опять переходит в состояние отработки мёртвого времени. В состоянии отработки мёртвого времени схема находится в течение времени и ни на что не реагирует, затем переходит в основное состояние.

Временной интервал 0 для электронной схемы совпадений, за который она успевает вернуться в исходное состояние при отсутствии сигналов от сцинтилляционных детекторов, составляет +. Для модели, состоящей из моделей сцинтилляционных детекторов и модели схемы совпадений, в качестве интервала 0 выбирается сумма интервалов 0 для обоих моделей, равная + + +.

1.6 Моделирование интроскопа, работающего по методу меченых нейтронов 1.6.1 Метод меченых нейтронов В настоящее время возможности применения метода меченых нейтронов (ММН) для идентификации скрытых веществ исследуется во многих лабораториях мира [17, 26, 42, 47, 49, 61]. К практическим применениям ММН относятся обнаружение запрещённых к транспортировке веществ в закрытых контейнерах, поиск взрывчатых веществ, каротаж скважин и другие.

Суть метода меченых нейтронов (ММН) заключается в следующем.

При бомбардировке тритиевой мишени потоком дейтронов в ходе реакции образуются быстрые моноэнергетические нейтроны с энергией 14, 1 МэВ, и альфа-частицы с энергией 3, 5 МэВ, которые разлетаются под углом примерно 180 градусов в системе центра масс. Нейтроны образуют поток зондирующего излучения.

Направление и момент времени вылета альфа-частицы измеряется специальным многосекционным детектором. По этим данным вычисляется направление и момент времени вылета сопутствующего ей нейтрона.

Нейтроны, проникая в инспектируемый объект, взаимодействуют с его веществом, и в результате протекания реакции образуют гамма-кванты с характерным для инспектируемого вещества энергетическим спектром.

Разница во времени между регистрацией гамма-кванта и альфа частицы позволяет судить о месте взаимодействия нейтрона с веществом и, таким образом, позволяет получить информацию о пространственном положении изотопа, испустившего гамма-квант. Измерение спектра индуцированных гамма-квантов позволяет идентифицировать вещество.

Во многих приложениях ММН остро стоит вопрос об уменьшении времени, необходимом для измерения спектра гамма-излучения, но достаточного для проведения идентификации. К решению этого вопроса подходят с двух сторон. Первый путь – оптимизация параметров измерительного прибора. Второй – разработка новых алгоритмов распознавания измеренных спектров. В обоих случаях используются методы математического моделирования, так как они позволяют существенно сэкономить время на разработку прибора и ресурсы.

1.6.2 Моделирование прибора Приборы, работающие по методу ММН, состоят из генератора меченых нейтронов, одного или нескольких детекторов гамма-излучения и электронной схемы, обрабатывающей сигналы детекторов. Алгоритм обработки сингалов отдельных детекторов сложен и для системы EURITRACK описан, например, в [53]. Возможно построение упрощенной модели прибора, работающего по методу меченых нейтронов, на основе уже разработанных моделей сцинтилляционного детектора и электронной схемы совпадений.

При моделировании траекторий нейтронов методом Монте-Карло время и направление вылета нейтрона известно, для их определения разыгрывать траекторию альфа-частицы нет необходимости. Вместо этого каждый нейтрон, покинувший источник, считается меченым. В настоящей работе траектории альфа-частиц не разыгрываются.

Предполагается, что функционал, описывающий отклик прибора, аддитивен по траекториям. При этом выбрасываются из рассмотрения случайные совпадения, вызванные возвращением многократно рассеянных гамма-квантов. Другими словами, не рассматриваются случайные совпадения, вызванные неправильным сопоставлением попавшего в детектор гамма-кванта и испущенного нейтрона. Это предположение оправдывается тем, что эксплуатация прибора в условиях, когда указанные случайные совпадения будут играть существенную роль, не предполагается.

Построенная в этой главе модель прибора характеризуется двумя параметрами – интервал времени от момента вылета нейтрона до момента открытия временного окна и длительность временного окна. Оба параметра выбираются так, чтобы во временное окно попадали гаммакванты, приходящие из инспектируемого объекта.

Модель интроскопа (см. рис. 12), работающего по методу меченых нейтронов, построена из модели сцинтилляционного детектора (п. 1.3.4) и модели спектрометра (п. 1.4). К модели спектрометра добавлена проверка условия, позволяющего осуществить стробирование, то есть выделить только те сигналы, которые попадают во временной интервал.

Входными данными всей модели являются события + и, соответстсвующие пересечению траекторий 1, 2,..., с границами чувствительного вещества детектора. Выходными данными является массив, отражающий количество попаданий в каждый из каналов спектрометра.

Рис. 12. Модель прибора, работающего по методу меченых нейтронов.

1.7 Моделирование комптоновского гаммаспектрометра 1.7.1 Приложения и принцип работы комптоновских Областью приложения гамма-телескопов помимо гамма-астрономии является дистанционное обнаружение ядерного оружия, утечек в местах захоронения радиоактивных отходов и на территориях действующих атомных электростанций, контроль за перемещением делящихся материалов, медицинская томография и другие.

Одной из разновидностей приборов, позволяющих определить направление прилёта гамма-квантов, являются комптоновские гаммаспектрометры. Ярким примером комптоновского гамма-спектрометра является гамма-телескоп COMPTEL.

На рис. 13 приведена геометрическая конфигурация комптоновского гамма-спектрометра. Рассмотрим принцип работы этого прибора.

Прибор состоит из двух слоев сцинтилляционных детекторов. Пусть гамма-квант испытал одиночное комптоновское рассеяние в одном из детекторов первого слоя и полное поглощение в детекторе второго слоя.

Тогда можно определить исходную энергию гамма-кванта и энергию, потерянную им при комптоновском рассеянии.

Потерянная при рассеянии энергия измеряется первым детектором, потерянная при полном поглощении энергия 1 – вторым. Исходная энергия 0 определяется как сумма 0 = 1 +.

Зная исходную энергию и потерянную энергию можно вычислить угол, на который отклонился гамма-квант при комптоновском рассеянии, из соотношения где 2 – энергия покоя электрона. С учётом пространственного расположения детекторов можно определить конус возможных направлений на источник. В идеальном случае конусы, соответствующие таким совпадениям, пересекаются в пространстве по лучу, направленному на точечный источник.

В качестве разрешающего времени электронной схемы совпадений выбирается время, достаточное, чтобы рассеянный в первом слое гамма-квант достиг второго слоя.

Рис. 13. Схема комптоновского гамма-спектрометра 1.7.2 Модель комптоновского гамма-спектрометра Структура модели комптоновского гамма-спектрометра повторяет структуру физического прибора (рис. 14). Один и тот же алгоритм, рассмотренный при построении модели сцинтилляционного детектора в п.

1.3.4, используется для моделирования работы всех сцинтилляционных детекторов в приборе. При моделировании работы электронной схемы совпадений используется алгоритм, рассмотренный в п. 1.5.

Временной интервал 0, за который модель комптоновского гаммаспектрометра успевает вернуться в исходное состояние, выбирается равным + + +, как это указано в п. 1.5.

Разработанная модель комптоновского гамма-спектрометра имеет две отличительные черты. Первая из них – это учёт совпадений, вызванных разными частицами (функционал, описывающий модель, не является аддитивным по траекториям). Для учёта таких совпадений необходимо обращать внимание только на момент времени обнаружения частицы каким-либо из детекторов, не учитывая номер разыгрываемой ветвящейся траектории. В большинстве программ, предназначенных для подсчёта совпадений, программа возвращается в исходное состояние Рис. 14. Модель комптоновского гамма-спектрометра.

после окончания моделирования каждой ветвящейся траектории, то есть вычисляется функционал (1 ) + (2 ) +... + ( ). Так сделано, например, в работах [46, 75, 84].

Прогонка модели выполняется в два этапа. Сначала производится моделирование работы сцинтилляционных детекторов. Выходными данными на этом шаге являются хронологические последовательности откликов (, ), полученные для каждого сцинтилляционного детектора. Затем из откликов всех детекторов образуется единая последовательность, которая в хронологическом порядке подаётся на вход модели электронной схемы совпадений.

Вторая особенность разрабатываемой модели комптоновского гамма-спектрометра состоит в учёте типа реакций, произошедщих в детекторе первого слоя и детекторе второго слоя. Учёт типа реакций позволяет определить, является совпадение истинным или ложным.

После раздельного подсчёта количества истинных и ложных совпадений можно вычислить отношение где – полное количество зарегистрированных совпадений, – количество зарегистрированных истинных совпадений. Коэффициент имеет смысл избирательности прибора. Геометрические параметры прибора должны быть выбраны так, чтобы максимизировать при определённых условиях работы прибора. Как в рамках модели различать истинные совпадения и ложные совпадения разных типов описано в п.

1.7.3.

Для комптоновского спектрометра с тремя рассеяниями похожие исследования проводились в работе [64], целью которых была оптимизация геометрических параметров прибора, при которых регистрируется максимальное количество истинных совпадений. Однако в этой работе не проводится различие между ложными совпадениями различных типов, и не учитываются ложные совпадения, вызванные разными частицами (прибор описывается аддитивным по траекториям функционалом).

1.7.3 Определение типа совпадений В случае, когда событие в первом детекторе вызвано одиночным комптоновским рассеянием, а событие во втором детекторе – поглощением гамма-кванта, то совпадение является истинным. Другие совпадения являются ложными. На практике отличить истинные совпадения от ложных невозможно. Поэтому при разработке приборов, работающих на эффекте двойных рассеяний, важно минимизировать вероятность ложных совпадений.

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

Для обеспечения последующих этапов разрабатываемого алгоритма информацией, из которой можно установить истинность или ложность совпадения, каждый отчёт сцинтилляционного детектора дополняется следующей информацией: число вошедших в детектор частиц, число вышедших частиц и число соударений внутри детектора за время интегрирования сигнала ФЭУ. С учётом сделанных ограничений на моделирование физических процессов можно утверждать, что если все три числа равны единице, то внутри детектора произошло одиночное комптоновское рассеяние. Если число вышедших частиц равно нулю, а число вошедших равно единице, то в детекторе произошло поглощение гамма-кванта.

Итак, совпадение будет истинным при выполнении следующих шести условий:

1. количество частиц, вошедших в детектор первого слоя, равно 1;

2. количество частиц, вышедших из детектора первого слоя, равно 1;

3. количество соударений внутри детектора первого слоя равно 1;

4. количество частиц, вошедших в детектор второго слоя, равно 1;

5. количество частиц, вышедших из детектора второго слоя, равно 0.

6. частицы, вошедшие в детектор первого слоя и в детектор второго слоя, принадлежат одной и той же ветвящейся траектории.

Проверка этих шести условий позволяет отнести каждое зарегистрированное двойное совпадение к одной из 26 групп, причём только одна из этих групп будет соответствовать истинным совпадениям. Подсчитав число попаданий в эту группу можно вычислить коэффициент. При помощи анализа распределения количества совпадений по этим группам можно указать способы повышения коэффициента.

1.8 Результаты Функционал, описывающий отклик прибора, регистрирующего ядерное излучение, может быть аддитивным по столкновениями, аддитивным по траекториям или неаддитивным. Выполнение того или иного свойства аддитивности связано с тем, события каких типов должны быть учтены в математической модели прибора. Неаддитивные функционалы являются наиболее общим случаем. Показано, что неаддитивные функционалы, описывающие отклик физических приборов, всё же могут быть разбиты на сумму слагаемых. Условием разбиения является отсутствие частиц в течение длительного интервала времени. Такое свойство аддитивности названо аддитивностью по времени. Доказана независимость и одинаковая распределённость слагаемых, на которые распадается аддитивный по времени функционал.

В качестве отклика сцинтилляционного детектора выбрана модель, согласно которой за отклик детектора принимается энергия, оставленная всеми частицами за время интегрирования сигнала ФЭУ. Выбранная модель обобщена на случай неаддитивного по траекториям функционала.

Модель сцинтилляционного детектора дополнена возможностью учёта типа реакции, что позволяет отделить в рамках построенной на её основе модели комптоновского гамма-спектрометра истинные совпадения от ложных. Предложен критерий для оптимизации геометрических параметров комптоновского гамма-спектрометра.

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

Все изложенные модели физических приборов являются имитационными. Модели приборов, работающих по схемам совпадений, строятся из моделей сцинтилляционного детектора и электронной схемы совпадений.

Модели приведены к одному формату входных данных и могут быть описаны функционалом (1, 2,..., ). Это позволит абстрагироваться от особенностей конкретной модели и построить общие универсальные весовые методы.

Глава Весовой алгоритм вычисления среднего неаддитивных функционалов среднего небольцмановских функционалов В классических работах, в частности в книгах В. В. Учайкина, А. А. Лагутина [25], С. М. Ермакова [7], Г. А. Михайлова [16], разработка весовых методов ведётся для вычисления функционалов определённых классов. Наиболее изученным классом является класс больцмановских функционалов, представимых в виде где — фазовые координаты частицы, () — плотность столкновений, () — функция детектора, интегрирование ведётся по всему фазовому пространству. Также разработаны методы для вычисления квадрата аддитивного функционала или любой его натуральной степени [12, 13].

Построение общего численного метода Монте-Карло для произвольных функционалов в этих и многих других работах ограничивается аналоговым методом.

В работе [25] рассматриваются детекторы, отклик которых представим в виде главной части, являющейся аддитивной, и малой неаддитивной добавки. Для вычисления показаний детекторов такого типа развит подход, называемый поливариантным разложением. Теоретически такой подход годится для вычисления среднего любых неаддитивных функционалов, но на практике он применим в тех случаях, когда все члены разложения, начиная с некоторого номера, становятся равными нулю.

В некоторых работах учитывается только одна исходящая частица даже в тех случаях, когда из реакции может выходить несколько частиц.

Несмещённость оценки при этом достигается корректировкой статистического веса. Так сделано, например, в работах [1, 41], в которой возможность учёта только одной выходящей частицы позволила разработать сопряжённые методы вычисления небольцмановских функционалов. Возможность учёта только одной частицы при моделировании реакции физически означает то, что вероятность одновременного попадания в детектор нескольких ветвей траектории мала. Это условие выполняется, например, в задачах, где чувствительная область детектора мала.

Предложенный в работах [1, 41] общий подход для расчёта небольцмановских функционалов заключается в расширении фазового пространства дополнительной переменной. Например, в случае расчёта спектра потерянной энергии — это поглощённая энергия частицы, в случае моделирования схем совпадений — это индикаторы внесения вклада в показания различных детекторов.

Обобщая сказанное выше, можно утверждать, что в большей части литературы, посвящённой весовому моделированию, рассматриваются различные специальные классы небольцмановских функционалов. Совершенно очевидно, что модели детекторов во многих задачах могут не попадать ни в один из рассмотренных классов. Также можно заключить, что учёт ветвлений в траекториях представляет собой отдельную сложность.

От указанных недостатков свободен общий теоретиковероятностный подход, согласно которому вычисляется среднее значение случайной величины, заданной на множестве всех ветвящихся траекторий, по некоторой вероятностной мере. Этот подход, намеченный в классических учебниках [23, 58], адаптирован для численных расчётов в работах [7, 18]. В работах [3, 8, 16] на достаточно строгом уровне изложены весовые методы Монте-Карло для вычисления среднего случайной величины, заданной на абстрактном вероятностном пространстве.

В настоящей работе концепция супертреков обосновывается в рамках общего теоретико-вероятностного подхода применением весовых методов Монте-Карло для вычисления интеграла где интегрирование ведётся по всем возможным ветвящимся траекториям, — вероятностная мера, () — усредняемый функционал. Исходом случайного испытания в методе Монте-Карло является ветвящаяся траектория. При вычислении этого интеграла можно ветвящиеся траектории разыгрывать со смещенной плотностью вероятности (), а в качестве весового множителя использовать () = ()/ (), где () — несмещённая (физическая) плотность вероятности реализации траектории.

Концепция супертреков привлекает отсутствием существенных ограничений на функционал (). Ценность концепции супертреков состоит в её универсальности, которая проявляется в возможности построения весового алгоритма моделирования во всех случаях, когда возможно построение аналогового алгоритма. Однако в работе [37] строгому обоснованию концепции супертреков уделено недостаточно внимания, но приводится разбор многочисленного множества примеров.

2.2 Концепция супертреков В работе [37] предлагается «запоминать» при моделировании физически совместно осуществимые траектории первичной и всех вторичных частиц, и при вычислении отклика детектора рассматривать такие траектории совместно. Совокупность траектории первичной частицы и всех вторичных частиц в [37] названа супертреком. Супертрек является исходом статистического испытания при вычислении интеграла (2.1) методом Монте-Карло.

Супертрек — это минимальная единица, к которой могут быть применены весовые методы Монте-Карло. Статистический вес приписывается всему супертреку, а не каждому столкновению в отдельности.

Новый супертрек создается вместе с рождением новой частицы в некоррелированном источнике. Этот супертрек имеет вес = и содержит траекторию только вновь рожденной частицы. Если в результате какого-нибудь физического процесса рождаются вторичные частицы, то их траектории добавляются к тому же супертреку, которому принадлежит траектория первичной частицы.

Рассмотрим весовые методы разыгрывания супертреков. При разыгрывании траектории частицы со смещённой плотностью вероятности без использования супертреков её вес умножается на корректирующий множитель. При использовании концепции супертреков на этот множитель умножается вес всего супертрека.

При расщеплении траектории частицы образуются несколько взаимоисключающих траекторий той же самой частицы. Они являются взаимоисключающими, так как физическая частица пошла бы только по одной из них. Если расщепляется одна из траекторий частиц, входящих в супертрек, то расщеплению подвергается весь супертрек. В каждый из новых супертреков входит по одной из получившихся ветвей траектории, подвергшейся расщеплению, и все остальные траектории из исходного супертрека. При этом каждый из новых супертреков по-прежнему содержит только физически осуществимые вместе траектории, и не одна из траекторий не потеряна. Но новые супертреки являются взаимоисключающими, так как они содержат взаимоисключающие траектории одной и той же частицы. Новые супертреки имеют суммарный вес, равный весу исходного супертрека.

При игре в русскую рулетку обрывается весь супертрек (его вес становится равным нулю). Было бы неправильно обрывать отдельные траектории частиц, входящих в супертрек. Аналогичные рассуждения можно привести и для других весовых методов. В частности, при применении метода DXTRAN (см., например, [28, 90]) супертрек расщепляется на два новых супертрека. В один из них входит DXTRANветвь траектории частицы, в другой — обычная траектория. После расщепления супертрек, содержащий обычную частицу (не-DXTRAN), играет в русскую рулетку. Более подробно суть метода DXTRAN описана в п. 2.7.5.1 настоящей главы.

Вместо вычисления величины где S — супертрек, в концепции супертреков предлагается использовать весовые методы при моделировании траекторий частиц и вычислять величину где — вес супертрека. В рамках концепции супертреков предложен способ вычисления веса в процессе разыгрывания супертрека при применении метода расщепления частицы, метода русской рулетки, поглощения уменьшением веса и разыгрывания траекторий со смещенной плотностью вероятности рассеяния. При использовании этого способа математическое ожидание величины * будет равно математическому ожиданию величины. Особое внимание уделено возможности реализации этого подхода с использованием программы MCNP. Однако автор концепции ограничился рассмотрением функционалов, аддитивных в смысле (1.2), где каждое слагаемое () зависит только от одного супертрека. В настоящей работе приводится обобщение концепции супертреков на случай, когда 2.3 Формулировка предлагаемого алгоритма Для того, чтобы подчеркнуть использование предложенных в [37] весовых методов будем использовать термин «супертрек» для обозначения ветвящейся траектории.

Пусть в результате весового моделирования супертреков -ый супертрек расщепился на супертреков с весами, = 1.., = 1... Предлагаемый алгоритм сводится к вычислению среднего значения случайной величины вместо вычисления среднего значения случайной величины Целью этой главы является доказательство несмещённости оценки, то есть доказательство равенства Практический смысл этого равенства состоит в том, что вместо проведения большого количества статистических испытаний для вычисления среднего величины, можно вычислять среднее величины *, проведя гораздо меньшее количество испытаний с использованием весовых методов. К весовым методам, применимость которых к вычислению * доказана в настоящей работе, относятся: моделирование со смещенной плотностью вероятности, расщепление, русская рулетка и DXTRAN.

2.4 Фазовое пространсво Будем описывать состояние частицы тремя величинами где x — пространственные координаты частицы, — единичный вектор направления скорости частицы, — «обобщенная» энергия частицы.

Она совпадает с энергией частицы в случае, если моделируются частицы только одного типа. Точка является элементом фазового пространства В фазовом пространстве выделим полуоткрытый прямоугольник = { : <, = 1, 2,..., 7} (расчётная область), где обозначает -ую координату точки из фазового пространства, < для всех = 1, 2,..., 7. Будем предполагать, что.

В случае моделирования блужданий частиц разных типов в переменной учитывается не только энергия, но и тип частицы, используя концепцию обобщенной частицы. Согласно этой концепции тип частицы определяется дополнительной переменной, принимающей дискретные значения. Пусть — физическая энергия частицы, а — тип этой частицы, = 1, 2,...,, где — количество типов частиц в задаче.

Пусть множество возможных значений энергий частиц ограничено, то есть 0 <. Тогда под величиной понимается Приведённая формула устанавливает взаимнооднозначное соответствие между физической энергией частицы и её типом с одной стороны и используемой в дальнейшем величиной с другой стороны (по этой причине величину будем называть в дальнейшем просто энергией).

Благодаря этому при определении математических ожиданий вычисляемых величин можно проводить усреднение только по энергии, не проводя усреднение по типам частиц.

2.5 Вероятностное пространство 2.5.1 Способы описания ветвящихся траекторий В классических работах [23, 58] рассмотрены разные случаи ветвящихся процессов, из которых нас будут интересовать только процессы с дискретным временем (время – номер поколения). Эти процессы можно поделить на две группы. В первой группе относятся процессы, в которых в каждый момент времени важно знать только фазовые координаты частиц из текущей популяции. Ко второй группе – процессы, в которых кроме фазовых координат частиц важно знать еще и родословную частиц (какое место занимает частица в дереве поколений).