NewProcess

Revision as of 12:24, 26 Nov 2005
217.20.209.178 (Talk | contribs)

← Go to previous diff
Current revision
217.20.214.8 (Talk | contribs)
Введение
Line 1: Line 1:
==Введение== ==Введение==
-<span style="color: darkred;"><em>Предупреждаю! Эта статья еще не дописана.</em></span>+Geant знает достаточно большое количество физических моделей. Они в него зашиты, как, просто, основная часть кода. Зайдите интереса ради в папку <tt>$G4INSTALL/source/processes</tt> и оглядитесь вокруг&nbsp;&mdash; одни только названия файлов о многом расскажут. Но, несмотря на обилие встроенных моделей, иногда может понадобиться примитивная в сущности вещь: сделать так, чтобы Geant считал по вашим сечениям. Об том, как это сделать, и написана эта статья.
-Geant знает достаточно большое количество физических моделей. Они в него зашиты, как, просто, основная часть кода. Зайдите интереса ради в папку <tt>$G4INSTALL/source/processes</tt> и оглядитесь вокруг&nbsp&mdash; одни только названия файлов о многом расскажут. Но, несмотря на обилие встроенных моделей, иногда может понадобиться примитивная в сущности вещь: сделать так, чтобы Geant считал по вашим сечениям. Об том, как это сделать, и написана эта статья.+
==Проверено на себе== ==Проверено на себе==
Line 10: Line 9:
Поглядев в файлах, я установил, что сечения были даны в виде таблицы &laquo;Энергия&mdash;Сечение&raquo;. Всего было 6 файлов для 6 реакций: (gamma, n), (gamma, 2n) и так далее. Да, и все это было для '''одного единственного''' ядра, <sup>203</sup>Tl (таллий). Поглядев в файлах, я установил, что сечения были даны в виде таблицы &laquo;Энергия&mdash;Сечение&raquo;. Всего было 6 файлов для 6 реакций: (gamma, n), (gamma, 2n) и так далее. Да, и все это было для '''одного единственного''' ядра, <sup>203</sup>Tl (таллий).
 +
 +Для задействования внешних сечений было решено написать свой процесс, поскольку тольщина таллиевого образца была достаточно солидной, чтобы гамма-кванты порядочно застревали в нем до того, как произойдет фотоядерная реакция.
 +
 +==G4DiscreteProcess==
 +
 +В моем случае можно было пренебречь вторичными фотонами, которые могут родиться в результате фотоядерной реакции, т.&nbsp;к. вероятность для такого фотона вновь провзаимодействовать с ядром ничтожна, а сами по себе фотоны меня не интересовали. Кроме того, для корректного их учета необходимо как минимум знать угловые сечения для данных реакций, а у меня было только интегральное сечение вылета одного или нескольких нейтронов&nbsp;&mdash; узнать направление вторичного фотона по этим данным мне представляется затруднительным. То же самое можно сказать и о нейтронах: я не знаю куда они полетят и с какой энергией, и, кроме того, мне это не важно. Поэтому решено было просто учитывать факт реакции, а затем просто убивать исходный фотон.
 +
 +Ядерная реакция подходит под критерий дискретного процесса (см. в документации GEANT про типы процессов), поэтому для реализации необходимо создать класс, производный от <tt>G4VDiscreteProcess</tt>.
 +
 + 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;
 + };
 +
 +Как видно, в конструкторе не делается ничего особенно важного, только загружаются сами сечения из файлов. Я не буду вдаваться в детали интерполяции значений сечений, их формата на диске и т.&nbsp;д., а остановлюсь на действительно важных деталях. Два метода в этом классе, которые и реализуют всю физику процесса это <tt>PostStepDoIt</tt> и <tt>GetMeanFreePath</tt>. Последний вызывается каждый раз при очередном шажке частицы в симуляции и определяет вероятность вызова данного конкретного процесса. <tt>PostStepDoIt</tt>&nbsp;&mdash; это то, что делает процесс, когда он все-таки происходит.
 +
 +
 + 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;
 + }
 +
 +Здесь по простой формуле вычисляется из сечения свободный пробег, который и возвращается. <tt>crossec</tt> уже переведен во внутренние единицы GEANT4 (кв. мм). Проверка с if на второй строчке это всего лишь костыль, спасающий от реальности: даны сечения только для таллия, так что ни в каких других материалах процесс не будет активизирован, т.&nbsp;к. в этом случае он вернет пробег 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;
 + }
 +
 +<tt>aParticleChange</tt>&nbsp;&mdash; это член класса '''G4VDiscreteProcess''', в котором вы должны описать все изменения, которые вызываются вашим процессом. В данном случае, как видно, я устанавливаю статус трекинга частицы в <tt>fStopAndKill</tt>, что соответствует выводу ее из процесса симуляции. Установка энергии строчкой выше, на самом деле не нужна, и здесь только ради перестраховки. Кроме того, я вывожу в '''stdout''' сообщение о срабатывании процесса.
 +
 +Подключается такой самодельный процесс точно так же, как и встроенные:
 +
 + pmanager->AddDiscreteProcess(new CuteProcess("AyAy", "/home/hatta/work/203Tl/", n));

Current revision

Введение

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

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

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

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

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

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));
Edit page