基于混合蒙特卡罗方法的带权重沉积能量能谱处理

1. 问题描述

为了加速测井仪器在Geant4中的模拟，我们常需要使用减方差方法去对模拟过程进行采样。目前在核加速领域比较常见，且应用效果较好的方法是CADIS方法。
因此我们将CADIS方法移植进入了Geant4软件中。基于一个密度仪器进行了测试，由于这个方法其本质上是基于权窗进行采样，但是在实现过程中出现了以下几点问题：

1） Geant4中基于平行几何体进行权窗划分的算法中，不支持多能群的采样
2）在基于Geant4的沉积能量统计方法(G4PSEnergyDeposition)输出的Root能谱发现，探测器上的计数值，少且集中在低能区：例如

2. 问题原因

2.1 能群问题：

在存储权窗数据的方法G4WeightWindowStore中，仅对单能群进行了存储。

2.2 沉积能量能谱问题：

1） G4PSEnergyDeposition方法的权重设置
在G4PSEnergyDeposition方法中，当前输出的沉积能量=沉积能量*权重；因此如果权重特别小，会导致乘出来的沉积能量会特别特别小。这就造成了大量粒子集中在低能区。
2） Root的统计方法
传统的Root统计，基于一个事件，至多只会产生一个计数。但是正对于我们的问题关闭权重输出的能谱也显然不是正确的能谱；如下图所示：

原因是由于分裂后，粒子本身也是有权重的，代表粒子的有效贡献。因此完全关闭权重的做法也是错误的。

3. 解决方案

3.1 能群问题：

对于粒子，我们采用了n28g19分群。

采用对粒子能量单个遍历的方式，实现了代码优化：

Geant4函数 G4WeightWindowStore.cc文件

点击查看代码

#include "MyWeightWindowStore.hh"
#include "G4VPhysicalVolume.hh"
#include "G4LogicalVolume.hh"
#include "G4GeometryCellStepStream.hh"
#include "G4TransportationManager.hh"

// ***************************************************************************
// Static class variable: ptr to single instance of class
// ***************************************************************************
G4ThreadLocal MyWeightWindowStore* MyWeightWindowStore::fInstance = nullptr;

MyWeightWindowStore::
MyWeightWindowStore()
 : fWorldVolume(G4TransportationManager::GetTransportationManager()
                ->GetNavigatorForTracking()->GetWorldVolume()),
   fGeneralUpperEnergyBounds(),
   fCellToUpEnBoundLoWePairsMap(),
   fCurrentIterator(fCellToUpEnBoundLoWePairsMap.cend())
{
}

MyWeightWindowStore::
MyWeightWindowStore(const G4String& ParallelWorldName)
 : fWorldVolume(G4TransportationManager::GetTransportationManager()
                ->GetParallelWorld(ParallelWorldName)),
   fGeneralUpperEnergyBounds(),
   fCellToUpEnBoundLoWePairsMap(),
   fCurrentIterator(fCellToUpEnBoundLoWePairsMap.cend())
{
}

MyWeightWindowStore::
~MyWeightWindowStore()
{
}

G4double MyWeightWindowStore::
GetLowerWeight(const G4GeometryCell& gCell, 
                     G4double partEnergy) const
{
  SetInternalIterator(gCell);
  auto gCellIterator = fCurrentIterator;
  if (gCellIterator == fCellToUpEnBoundLoWePairsMap.cend())
  {
    Error("GetLowerWitgh() - Cell does not exist!++++++++++++");
    return 0.;
  }
  G4UpperEnergyToLowerWeightMap upEnLoWeiPairs = fCurrentIterator->second;

  G4double lowerWeight = -1;
  G4bool found = false;

  G4double *numpy = new G4double[30];
  G4int j=0;
  for (auto it = upEnLoWeiPairs.cbegin(); it != upEnLoWeiPairs.cend(); ++it)
  {
  numpy[j]=it->first;
  j++;
  }
  j=0;
  for (auto it = upEnLoWeiPairs.cbegin(); it != upEnLoWeiPairs.cend(); ++it)
  {

    if(it->first != 0.000){

    if ( numpy[j] < partEnergy && partEnergy<= it->first  )
    {
      lowerWeight = it->second;
      found = true;
    }
        j++;
  }
  }
  if (!found)
  {
    std::ostringstream err_mess;
    err_mess << "GetLowerWitgh() - Couldn't find lower weight bound." << G4endl
             << "Energy: " << partEnergy << ".";
    Error(err_mess.str());
  }
  delete numpy;
  return lowerWeight;
}

void MyWeightWindowStore::
SetInternalIterator(const G4GeometryCell& gCell) const
{
  fCurrentIterator = fCellToUpEnBoundLoWePairsMap.find(gCell);
}

G4bool MyWeightWindowStore::
IsInWorld(const G4VPhysicalVolume& aVolume) const
{
  G4bool isIn(true);
  if (!(aVolume == *fWorldVolume))
  {
    isIn = fWorldVolume->GetLogicalVolume()->IsAncestor(&aVolume);
  }
  return isIn;
}

G4bool MyWeightWindowStore::
IsKnown(const G4GeometryCell& gCell) const
{
  G4bool inWorldKnown(IsInWorld(gCell.GetPhysicalVolume()));
                      
  if ( inWorldKnown )
  {
    SetInternalIterator(gCell);
    inWorldKnown = (fCurrentIterator!=fCellToUpEnBoundLoWePairsMap.cend());
  }
  return inWorldKnown;
}

void MyWeightWindowStore::Clear()
{
  fCellToUpEnBoundLoWePairsMap.clear();
}

void MyWeightWindowStore::SetWorldVolume()
{
  G4cout << " G4IStore:: SetWorldVolume " << G4endl;
  fWorldVolume = G4TransportationManager::GetTransportationManager()
                 ->GetNavigatorForTracking()->GetWorldVolume();
  G4cout << " World volume is: " << fWorldVolume->GetName() << G4endl;
  // fGeometryCelli = new G4GeometryCellImportance;
}

void MyWeightWindowStore::SetParallelWorldVolume(const G4String& paraName)
{
  fWorldVolume = G4TransportationManager::GetTransportationManager()
                 ->GetParallelWorld(paraName);
  // fGeometryCelli = new G4GeometryCellImportance;
}

const G4VPhysicalVolume &MyWeightWindowStore::GetWorldVolume() const
{
  return *fWorldVolume;
}

const G4VPhysicalVolume* MyWeightWindowStore::
GetParallelWorldVolumePointer() const
{
  return fWorldVolume;
}

void MyWeightWindowStore::
AddLowerWeights(const G4GeometryCell& gCell,
                const std::vector<G4double>& lowerWeights)
{
  if (fGeneralUpperEnergyBounds.empty())
  {
    Error("AddLowerWeights() - No general upper energy limits set!");
  }
  if (IsKnown(gCell))
  {
    Error("AddLowerWeights() - Cell already in the store.");
  }
  if (lowerWeights.size() != fGeneralUpperEnergyBounds.size())
  {
    std::ostringstream err_mess;
    err_mess << "AddLowerWeights() - Mismatch between "
             << "number of lower weights (" << lowerWeights.size()
             << ") and energy bounds (" << fGeneralUpperEnergyBounds.size()
             << ")!";
    Error(err_mess.str());
  }
  G4UpperEnergyToLowerWeightMap map;
  G4int i = 0;

  for (auto it = fGeneralUpperEnergyBounds.cbegin(); 
            it != fGeneralUpperEnergyBounds.cend(); ++it)
  {
/*     G4cout << "    fenB =    "  <<  *it << G4endl; */
    map[*it] = lowerWeights[i];
/*     G4cout << "    lowerWeights[i] =    "  <<  lowerWeights[i] << G4endl; */
    ++i;
  }
  fCellToUpEnBoundLoWePairsMap[gCell] = map;
  //G4cout << "    fCellToUpEnBoundLoWePairsMap[gCell].size =    "  << fCellToUpEnBoundLoWePairsMap[gCell].size() << G4endl;
}

void MyWeightWindowStore::
AddUpperEboundLowerWeightPairs(const G4GeometryCell& gCell,
                               const G4UpperEnergyToLowerWeightMap& enWeMap)
{
  if (IsKnown(gCell))
  {
    Error("AddUpperEboundLowerWeightPairs() - Cell already in the store.");
  }
  if (IsKnown(gCell))
  {
    Error("AddUpperEboundLowerWeightPairs() - Cell already in the store.");
  }
  fCellToUpEnBoundLoWePairsMap[gCell] = enWeMap;

}

void MyWeightWindowStore::
SetGeneralUpperEnergyBounds(const std::set<G4double,
                            std::less<G4double> >& enBounds)
{

  if (!fGeneralUpperEnergyBounds.empty())
  {
    Error("SetGeneralUpperEnergyBounds() - Energy bounds already set.");
  }
  fGeneralUpperEnergyBounds = enBounds;
    //G4cout << "    enB =    "  <<  fGeneralUpperEnergyBounds.size() << G4endl;
}
 
void MyWeightWindowStore::Error(const G4String& msg) const
{
  G4Exception("MyWeightWindowStore::Error()",
              "GeomBias0002", FatalException, msg);
}

// ***************************************************************************
// Returns the instance of the singleton.
// Creates it in case it's called for the first time.
// ***************************************************************************
//
MyWeightWindowStore* MyWeightWindowStore::GetInstance()
{
  if (fInstance == nullptr)
  {
    fInstance = new MyWeightWindowStore();
  }
  return fInstance;    
}

// ***************************************************************************
// Returns the instance of the singleton.
// Creates it in case it's called for the first time.
// ***************************************************************************
//
MyWeightWindowStore* MyWeightWindowStore::
GetInstance(const G4String& ParallelWorldName)
{
  if (fInstance == nullptr)
  {
#ifdef G4VERBOSE
    G4cout << "G4IStore:: Creating new Parallel IStore "
           << ParallelWorldName << G4endl;
#endif
    fInstance = new MyWeightWindowStore(ParallelWorldName);
  }
  return fInstance;    
}

平行几何体设置文件，权窗差入能量及分群输入权重：

点击查看代码

G4VWeightWindowStore *OilWell_paraGeo::CreateWeightWindowStore()
{
   G4cout << " OilWell_paraGeo:: Creating Importance Store " 
         << G4endl;
 fGhostWorld = fPhysicalVolumeVector[0];
  if (!fPVolumeStore.Size())
  {
    G4Exception("OilWell_paraGeo::CreateWeightWindowStore"
               ,"OilWell_paraGeo",RunMustBeAborted
               ,"no physical volumes created yet!");
  }
 MyWeightWindowStore *wwstore = MyWeightWindowStore::GetInstance();
//创建一个能量区间进行采样
 // std::set<G4double, std::less<G4double> > enBounds;
    enBounds.insert(0.8 * MeV);
    enBounds.insert(0.6 * MeV);
    enBounds.insert(0.4 * MeV);
    enBounds.insert(0.3 * MeV);
    enBounds.insert(0.2 * MeV);
    enBounds.insert(0.045 * MeV);
    enBounds.insert(0.001 * MeV);
    enBounds.insert(0.000 * MeV);
  wwstore->SetGeneralUpperEnergyBounds(enBounds);
  //
	vector<float> ivec;
  {
//"input weights to ivec"
}

    G4double lowerWeight18 =1;
    G4double lowerWeight17 =1;
    G4double lowerWeight16 =1;
    G4double lowerWeight15 =1;
    G4double lowerWeight14 =1;
    G4double lowerWeight13 =1;

    lowerWeights.push_back(1);
    lowerWeights.push_back(1);
    lowerWeights.push_back(1);
    lowerWeights.push_back(1);
    lowerWeights.push_back(1);
    lowerWeights.push_back(1);
    lowerWeights.push_back(1);
    lowerWeights.push_back(1);

  G4GeometryCell gWorldCell(GetWorldVolumeAddress(),0);
  wwstore->AddLowerWeights(gWorldCell, lowerWeights);
G4int cell(1);
  for (cell=1; cell<=45; cell++) {
      G4GeometryCell gCell = GetGeometryCell(cell);

          G4cout << " adding cell: " << cell 
           << " replica: " << gCell.GetReplicaNumber() 
           << " name: " << gCell.GetPhysicalVolume().GetName() << G4endl;

      lowerWeight18 =ivec18[cell-1];
      lowerWeight17 =ivec17[cell-1];
      lowerWeight16 =ivec16[cell-1];
      lowerWeight15 =ivec15[cell-1];
      lowerWeight14 =ivec14[cell-1];
      lowerWeight13 =ivec13[cell-1];

      lowerWeights.clear();
      //lowerWeights.push_back(lowerWeight);
      lowerWeights.push_back(NULL);
      lowerWeights.push_back(1);
      lowerWeights.push_back(lowerWeight18);
      lowerWeights.push_back(lowerWeight17);
      lowerWeights.push_back(lowerWeight16);
      lowerWeights.push_back(lowerWeight15);
      lowerWeights.push_back(lowerWeight14);
      lowerWeights.push_back(lowerWeight13);
      wwstore->AddLowerWeights(gCell, lowerWeights);
  }
  return wwstore;
  }

3.2 沉积能量能谱问题：

为了解决这个问题，我们必须要明白Geant4中如何对粒子进行的重要性采样的；减方差的原理就是对粒子进行分裂与轮盘赌操作，以便使得更多“可能对探测器有贡献的粒子-权重大的粒子”返回到探测器，如下所示

当粒子权重高于上权窗边界后，会进行分裂，而分裂后的粒子在Geant4中与原本的粒子输出上是没有区别的，分裂粒子性质与父粒子一摸一样。其具体实现过程如下所示：

而在Geant4中，粒子的运行规律是一个Event代表粒子的生死，每个Event由多个Step构成，每个Step代表着粒子上一次发生反应到下一次反应的过程；而所有的Event（所有粒子）的模拟过程为一次Run。

明白这一点就可以利用Geant4输出的特点，得到以下数据：

由于伽马是没有办法产生伽马的，因此对于同一个事件产生的多个Gamma粒子，除了第一个，就可以知道其它都是复制产生的粒子了。因此通过单独统计这些粒子的权重和，实际上就是采样后粒子的实际贡献度。

因此现在的统计方式我们改为，对应粒子的权重和。G4PSEnergyDeposit.cc文件：

点击查看代码

G4bool MyPSEnergyDeposit::ProcessHits(G4Step* aStep,G4TouchableHistory*)
{
  G4double edep = aStep->GetTotalEnergyDeposit();
  if ( edep == 0. ) return FALSE;
  //edep *= aStep->GetPreStepPoint()->GetWeight(); // (Particle Weight)
  G4double pWeight= aStep->GetPreStepPoint()->GetWeight();
  G4int  index = GetIndex(aStep);
  G4int eventID = G4RunManager::GetRunManager()->GetCurrentEvent() -> GetEventID();
    auto endPoint = aStep->GetPostStepPoint();
    auto procName = endPoint->GetProcessDefinedStep()->GetProcessName();
    auto track        = aStep->GetTrack();
  EvtMap->add(index,edep);  
  return TRUE;
}

4. 实现效果

搭建了一个简单的测井模型，其参数如图：

3D视图：

CADIS运行效果：

能谱对比（CADIS运行30w粒子/ANALOG运行5亿个粒子）：

	Score	Rel.	FOM
CADIS	230.9	2.46%	191.6
ANALOG	36833.1	2.7%	1.69

5. 其它

在植入权窗的时候，出现了一个小插曲，就是出现了内存溢出的问题；在很少粒子的时候内存崩溃，就被程序自动杀死了；

为了解决这个问题，也是头疼很久，检查代码；在网上查了有这几种原因造成的：

在类的构造函数和析构函数中没有匹配的调用new和delete函数
两种情况下会出现这种内存泄露:

(0）在堆里创建了对象占用了内存，但是没有显示地释放对象占用的内存;

(1）在类的构造函数中动态的分配了内存，但是在析构函数中没有释放内存或者没有正确的释放内存

(2）没有正确地清除嵌套的对象指针

(3）在释放对象数组时在delete中没有使用方括号

(4）指向对象的指针数组不等同于对象数组

(5）缺少拷贝构造函数

(6）缺少重载赋值运算符

(7）nonmodifying运算符重载

(8）没有将基类的析构函数定义为虚函数
借助这个博客https://blog.csdn.net/weixin_40774605/article/details/114376434
下载了一个C++内存泄漏检查工具，排除出来原来是在加入能群的时候，new了一个数组，没有delete导致的。改了之后，程序内存问题就完美解决了。