NewProcess

Введение

Geant знает достаточно большое количество физических моделей. Они в него зашиты, как, просто, основная часть кода. Зайдите интереса ради в папку $G4INSTALL/source/processes и оглядитесь вокруг — одни только названия файлов о многом расскажут. Но, несмотря на обилие встроенных моделей, иногда может понадобиться примитивная в сущности вещь: сделать так, чтобы Geant считал по вашим сечениям. Об том, как это сделать, и написана эта статья.

Проверено на себе

Лично мне пришлось заняться этим на первый взгляд бессмысленным делом, когда руководство разбраковало в пух и прах встроенные в Geant адронные модели и дало мне несколько файлов с фразой «Вот тебе надежные сечения». Не объяснив при этом даже формат хранения данных… но это к делу не относится.

Поглядев в файлах, я установил, что сечения были даны в виде таблицы «Энергия—Сечение». Всего было 6 файлов для 6 реакций: (gamma, n), (gamma, 2n) и так далее. Да, и все это было для одного единственного ядра, ²⁰³Tl (таллий).

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

G4DiscreteProcess

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

Ядерная реакция подходит под критерий дискретного процесса (см. в документации GEANT про типы процессов), поэтому для реализации необходимо создать класс, производный от G4VDiscreteProcess.

class CuteProcess: public G4VDiscreteProcess
{
public:
  /// Конструктор.
  /// @param dirname: директория с сечениями;
  /// @param nn: число выделяемых нейтронов;
  CuteProcess(const G4String& name,
              const G4String& dirname,
              int nn): 
    G4VDiscreteProcess(name, fUserDefined),
    nneutrons(nn)
  {
    loadFromDir(dirname, nn);
  }
  virtual G4VParticleChange *PostStepDoIt(const G4Track& track, const G4Step& step);
protected:
  virtual G4double GetMeanFreePath(const G4Track& aTrack,
                                   G4double previousStepSize,
                                   G4ForceCondition* condition);
  void loadFromDir(const G4String& dirname, int nn);
private:
  Maths::Interpolation::Linear* xs;
  int nneutrons;
};

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

G4double CuteProcess::GetMeanFreePath(const G4Track& aTrack,
                                      G4double previousStepSize,
                                      G4ForceCondition* condition)
{
  double energy = aTrack.GetKineticEnergy();
  if (aTrack.GetMaterial()->GetName() != "Thallium" 
  || aTrack.GetDefinition() != G4Gamma::GammaDefinition())
    return 1e100;
  double crosssec = xs->getValue(energy); // Сечение из файла.
  double l = 1/(6.022e23/mole*11.85*g/cm3/(204.3833*g/mole)*crosssec);
  if(l<0) l = 1e100;
  return l;
}

Здесь по простой формуле вычисляется из сечения свободный пробег, который и возвращается. crossec уже переведен во внутренние единицы GEANT4 (кв. мм). Проверка с if на второй строчке это всего лишь костыль, спасающий от реальности: даны сечения только для таллия, так что ни в каких других материалах процесс не будет активизирован, т. к. в этом случае он вернет пробег 1e100.

Вот, что мы делаем, когда процесс-таки вызывается:

G4VParticleChange* CuteProcess::PostStepDoIt(const G4Track& track, const G4Step& step)
{
  aParticleChange.Initialize(track);
  aParticleChange.ProposeEnergy(0.);
  aParticleChange.ProposeTrackStatus(fStopAndKill);
  std::cout << "gen_neutron: " << nneutrons << " more" << std::endl;
  return &aParticleChange;
}

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

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

pmanager->AddDiscreteProcess(new CuteProcess("AyAy", "/home/hatta/work/203Tl/", n));