粒子物理蒙特卡罗模拟库Geant4之能谱制作

$\gamma$衰变模型

---

---

抱怨！！

两日光阴浪费
在GEANT4那些物理学家写的奇怪逻辑文档中我难以找到“粒子$\gamma$衰变模型”。所以这个物理问题变成了一个文档搜寻问题。
找不到文档只能自己查找源代码，所以只能使用grep工具，结果由于我忘了使用shell通道
root@master:# grep -r gamma ./ | grep decay
以至于浪费两日时间。所以归根结底是一个Linux shell问题。：）

---

现实意义

能谱的制作在粒子物理中有着异常重要的意义。那么如何模拟粒子衰变以及收集能量沉积并进行数据分析变得非常必要。通过对被研究对象能谱的获取和分析可以直接或间接地获得物质的结构、组成元素的种类与含量等重要信息。

---

构建粒子类型、记录能谱及定义衰变链

void TrackingAction::PreUserTrackingAction(const G4Track* track)
{
  Run* run 
   = static_cast<Run*>(G4RunManager::GetRunManager()->GetNonConstCurrentRun());

  G4ParticleDefinition* particle = track->GetDefinition();
  G4String name   = particle->GetParticleName();
  fCharge = particle->GetPDGCharge();
  fMass   = particle->GetPDGMass();  

  G4double Ekin = track->GetKineticEnergy();
  G4int ID      = track->GetTrackID();

  G4bool condition = false;

  // check LifeTime
  //
  G4double meanLife = particle->GetPDGLifeTime();

  //count particles
  //
  if (ID>1) run->ParticleCount(name, Ekin, meanLife);

  //energy spectrum
  //
  G4int ih = 0;
  if (particle == G4Electron::Electron()||
      particle == G4Positron::Positron())  ih = 1;
  else if (particle == G4NeutrinoE::NeutrinoE()||
           particle == G4AntiNeutrinoE::AntiNeutrinoE()) ih = 2;
  else if (particle == G4Gamma::Gamma()) ih = 3;
  else if (particle == G4Alpha::Alpha()) ih = 4;
  else if (fCharge > 2.) ih = 5;
  if (ih) G4AnalysisManager::Instance()->FillH1(ih, Ekin);
  if (ih) G4CsvAnalysisManager::Instance()->FillH1(ih, Ekin);

  //Ion
  //
  if (fCharge > 2.) {
    //build decay chain
    if (ID == 1) fEvent->AddDecayChain(name);
      else       fEvent->AddDecayChain(" ---> " + name);
    // 
    //full chain: put at rest; if not: kill secondary      
    G4Track* tr = (G4Track*) track;
    if (fFullChain)  tr->SetTrackStatus(fStopButAlive);
      else if (ID>1) tr->SetTrackStatus(fStopAndKill);
    //
    fTime_birth = track->GetGlobalTime();
  }

  //example of saving random number seed of this fEvent, under condition
  //
  condition = (ih == 3);
  if (condition) G4RunManager::GetRunManager()->rndmSaveThisEvent();
}

严格来说、应该叫获取粒子类型，因为即将展示如何在交互命令中或者脚本中定义衰变粒子类型。

---

衰变脚本

锕系元素$Am241$的衰变脚本

/control/verbose 2
/run/verbose 1
#
/gun/particle ion
/gun/ion 95 241 #在这里定义了类型
#
/tracking/verbose 2
/run/beamOn 1
#
/control/cout/ignoreThreadsExcept 0
/tracking/verbose 0
#
/analysis/setFileName Am241 #存储数据文件为.root格式
/analysis/h1/set 1  150  0. 1500 keV    #e+ e-
/analysis/h1/set 2  150  0. 1500 keV    #neutrino
/analysis/h1/set 3  150  0. 1500 keV    #gamma
/analysis/h1/set 6  100  0. 2500 keV    #EkinTot (Q)
/analysis/h1/set 7  150  0. 15e3 keV    #P balance
/analysis/h1/set 8  100  0. 100. year   #time of life
/analysis/h1/set 9  100  1. 3. MeV      #EvisTot
#
/run/printProgress 100000
#  
/run/beamOn 1000

---

数据格式转换

.root格式是适用于欧洲原子能机构CERN的ROOT分析工具的数据格式，不适于其他数据工具使用。将其转换为.xml文件：

root@master:# sudo vi root2xml.C
{
      TFile *f = TFile::Open("Example.xml","recreate");
      TFile *mf = TFile::Open("Co60.root");
      TH1F *h = (TH1F*)mf->Get("3");
      f->cd();
      h->Write();
}  

可以得到.xml格式能谱数据：

<?xml version="1.0"?>
<root setup="2xoo" ref="null" created="2017-11-16 22:13:05" modified="2017-11-16 22:13:05" uuid="43680090-cad8-11e7-9717-0100007fbeef" version="3" file_version="53430">
  <XmlKey name="3" cycle="1" created="2017-11-16 22:13:45">
    <Object class="TH1D" ref="id0">
      <TH1D version="1">
        <TH1 version="7">
          <TNamed version="1">
            <TObject fUniqueID="0" fBits="3000008"/>
            <fName str="3"/>
            <fTitle str="energy spectrum (%): gamma"/>
          </TNamed>
          <TAttLine version="2">
            <fLineColor v="1"/>
            <fLineStyle v="1"/>
            <fLineWidth v="1"/>
          </TAttLine>
          <TAttFill version="2">
            <fFillColor v="0"/>
            <fFillStyle v="101"/>
          </TAttFill>
          <TAttMarker version="2">
            <fMarkerColor v="1"/>
            <fMarkerStyle v="1"/>
            <fMarkerSize v="1.000000e+00"/>
          </TAttMarker>
          <fNcells v="152"/>
          <fXaxis>
            <TAxis version="9">
              <TNamed version="1">
                <TObject fUniqueID="0" fBits="3000000"/>
                <fName str="xaxis"/>
                <fTitle str=" [keV]"/>
... ... ...

---

绘图

绘图脚本

{
      gROOT->Reset();

      // Draw histos filled by Geant4 simulation 
      //   
      TFile f("Am241.root");

      TCanvas* c1 = new TCanvas("c1", "  ");
      c1->SetLogy(1);
      c1->cd();
      c1->Update();

      TH1D* hist = (TH1D*)f.Get("3");
      hist->Draw("HIST");    
}  

绘图命令

root@master:# root -l Am241.root