Geant4的平行世界编程

1. 头文件需要包括

 

#include "G4VUserParallelWorld.hh"

class ParallelWorld : public G4VUserParallelWorld
{ 
 public:
    ParallelWorld(G4String worldName);
    virtual ~ParallelWorld();

    virtual void ConstructSD();

 private：
    virtual void Construct();
}

 

有两个虚函数，其中一个用来构建几何，另一个构建灵敏探测器，后者不是必须的。

 

2. 构造函数的代码段，获取平行世界的名字和指针。

void ParallelWorld::Construct()

{

      G4VPhysicalVolume* GhostWorld = GetWorld();

      G4cout << "ghost World Name= " << GhostWorld->GetName() << G4endl;

fworldLogic = GhostWorld->GetLogicalVolume();

}

 

3. 构建灵敏探测器。下面的代码使用了原生多功能探测器，最后一行是设置探测器到逻辑体。

void ParallelWorld::ConstructSD()

{

       G4SDManager::GetSDMpointer()->SetVerboseLevel(1);

       G4MultiFunctionalDetector* MultiFD = new G4MultiFunctionalDetector("ghost");

       G4VPrimitiveScorer* primitiv = new G4PSDoseDeposit("dose");

       MultiFD->RegisterPrimitive(primitiv);

G4SDManager::GetSDMpointer()->AddNewDetector(MultiFD);

//SetSensitiveDetector(fScoreLogic,MultiFD);

}

 

最后一行注释了，可以在其它函数中这样设置。

  

    G4VSensitiveDetector* sd = G4SDManager::GetSDMpointer()

     ->FindSensitiveDetector("ghost");

     fScoreLogic->SetSensitiveDetector(sd);

 

3.  必须在物理列表子类的成员函数ConstructProcess()中添加AddScoringProcess()。

 

void PhysicsList::ConstructProcess()

{

    AddTransportation();

    fEmPhysicsList->ConstructProcess();

    AddScoringProcess();

}

 

 

include "G4ParallelWorldScoringProcess.hh"

void PhysicsList::AddScoringProcess()

{

    G4ParallelWorldScoringProcess* theParallelWorldScoringProcess

        = new G4ParallelWorldScoringProcess("ParaWorldScoringProc");

    theParallelWorldScoringProcess->SetParallelWorld("ParallelScoringWorld");

    theParticleIterator->reset();

while((*theParticleIterator)()){ 

    G4ParticleDefinition* particle = theParticleIterator->value();

        if(!particle->IsShortLived()){

           G4ProcessManager* pmanager = particle->GetProcessManager();

           pmanager->AddProcess(theParallelWorldScoringProcess);                         

           pmanager->SetProcessOrderingToLast(theParallelWorldScoringProcess,idxAtRest);

           pmanager->SetProcessOrdering(theParallelWorldScoringProcess,idxAlongStep,1);          

           pmanager->SetProcessOrderingToLast(theParallelWorldScoringProcess,idxPostStep);

}}}

 

4. 主程序中的代码

 

    G4VUserDetectorConstruction* detector = new DetectorConstruction;

    detector->RegisterParallelWorld(new ParallelWorld("ParallelScoringWorld"));

    runManager->SetUserInitialization(detector);

5. 可视化场景中显示平行世界，宏命令：

# Draw geometry:

/vis/drawVolume worlds

 

 Geant4内置分析工具

6 . 在运行开始时创建一个二维直方图图

 

void RunAction::BeginOfRunAction(const G4Run* aRun)

{

    G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();

    analysisManager->SetVerboseLevel(0);

    analysisManager->SetFileName(fFileName);

 

    if(!fCreateH2) analysisManager->CreateH2("h2d","dose",NX,0,NY,NZ,0,NZ);

    else analysisManager->SetH2(0,NY,0,NY,NZ,0,NZ);

    fCreateH2 = true;

    analysisManager->OpenFile();

}

 

保存文件：

 

void RunAction::save()

{

    G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();

    analysisManager->Write();

    analysisManager->CloseFile();

}

 

7. 获取平行世界对象的指针，传递变量

    G4RunManager* mgr = G4RunManager::GetRunManager();

    const G4VUserDetectorConstruction* vdet = mgr->GetUserDetectorConstruction();

    const G4VUserParallelWorld* vpar = vdet->GetParallelWorld(0);

    ParallelWorld* para=(ParallelWorld*)vpar;

 

8. 记录数据，下面的代码把数据存入二维柱状图，copy编号按行列转二维。

 

void EventAction::BeginOfEventAction(const G4Event* evt)

{

   if (fHHCID==-1) {

      fHHCID = G4SDManager::GetSDMpointer()->GetCollectionID("ghost/dose");

}}

 

void EventAction::EndOfEventAction(const G4Event* evt)

{

   G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();

   G4HCofThisEvent* HCE = evt->GetHCofThisEvent();

   if (!HCE)

    {   G4ExceptionDescription msg;

        msg << "No hits collection of this event found.\n";

        G4Exception("EventAction::EndOfEventAction()","Code001", JustWarning,           msg);

        return;

    }

    G4THitsMap* evtMap1 = static_cast*>(HCE->GetHC(fHHCID));

    G4double dose = 0.;

    std::map::iterator itr;

    for (itr = evtMap1->GetMap()->begin(); itr != evtMap1->GetMap()->end(); itr++)

    {

       G4int copyNb  = (itr->first);

       G4double d = *(itr->second);

       dose += d;

       analysisManager->FillH2(0,(G4int)(copyNb/NZ),copyNb%NZ,d/gray);

}}

 

使用ROOT画图

 

9. ROOT解释执行的C++程序，所以写一个程序 .c打开运行，或在命令行状态直接输入表达式。下面是一个简单的实例代码，简单的只要写在大括号内。

 

{

   TCanvas *c = new TCanvas("H2D","TH2D Drawing", 0, 0,1280,728);

   TFile f("camera100000.root");      //文件名

   TH2D* h2d = (TH2D*)f->Get("h2d");  // book（）的名字

   c.Divide(2,2);

   c.cd(1);

      TH1D* h1= h2d->ProjectionX(); 

      h1->GetXaxis()->CenterTitle();

      h1>SetYTitle("dose"); 

      h1->GetYaxis()->CenterTitle();

      h1->SetMinimum(0);

      h1->SetStats(kFALSE);

      h1->GetYaxis()->SetNoExponent(kTRUE);

      h1->Draw("LP*"); 

   c.cd(2);

      TH1D* hy = h2d->ProjectionY(); 

      h2->SetMinimum(0);

      h2->SetStats(kFALSE);

      h2->GetYaxis()->SetNoExponent(kTRUE);

      h2->Draw("LP*");

    c.cd(3);

      h2d->SetStats(kFALSE);

      h2d->Draw("CONT"); //"SCAT","BOX","ARR","COLZ","CONTZ","CONT1"-2-3

    c.cd(4);

      h2d->Draw("SURF2"); //"SURF"-1-2-3-4-5,"LOGO"-0-1-2,+Z      

}