Geant4的root文件读取

梁佐佐

发布于 2020-09-04 17:29:39

2.4K0

发布于 2020-09-04 17:29:39

文章被收录于专栏：人芳觅人芳觅

来不及解释了，超长超稳教程送给大家~。~

exampleB5

B5Analysis.hh

#ifndef B5Analysis_h
#define B5Analysis_h 1
#include "g4root.hh"
//#include "g4xml.hh"
//#include "g4csv.hh"
#endif

B5RunAction.cc

#include "B5Analysis.hh"

B5RunAction::B5RunAction(B5EventAction*eventAction): G4UserRunAction(),fEventAction(eventAction)

{
auto analysisManager = G4AnalysisManager::Instance();
G4cout << "Using " << analysisManager->GetType() << G4endl;
analysisManager->SetNtupleMerging(true);
analysisManager->SetVerboseLevel(1);
analysisManager->SetFileName("B5");
// Book histograms, ntuple // Creating 1D histograms
analysisManager->CreateH1("Chamber1","Drift Chamber 1 # Hits", 50, 0., 50); // h1 Id = 0
analysisManager->CreateH1("Chamber2","Drift Chamber 2 # Hits", 50, 0., 50); // h1 Id = 1
// Creating 2D histograms
analysisManager->CreateH2("Chamber1 XY","Drift Chamber 1 X vs Y",50, -1000., 1000, 50, -300., 300.);// h2 Id = 0
analysisManager->CreateH2("Chamber2 XY","Drift Chamber 2 X vs Y",50, -1500., 1500, 50, -300., 300.);// h2 Id = 1
// Creating ntuple
if ( fEventAction ) {
analysisManager->CreateNtuple("B5", "Hits");
analysisManager->CreateNtupleIColumn("Dc1Hits"); // column Id = 0
analysisManager->CreateNtupleIColumn("Dc2Hits"); // column Id = 1
analysisManager->CreateNtupleDColumn("ECEnergy");// column Id = 2
analysisManager->CreateNtupleDColumn("HCEnergy");// column Id = 3
analysisManager->CreateNtupleDColumn("Time1"); // column Id = 4
analysisManager->CreateNtupleDColumn("Time2"); // column Id = 5
analysisManager->CreateNtupleDColumn("ECEnergyVector", fEventAction->GetEmCalEdep());// column Id = 6
problem // std::vector <<>> fEmCalEdep;
analysisManager->CreateNtupleDColumn("HCEnergyVector", fEventAction->GetHadCalEdep());// column Id = 7
analysisManager->FinishNtuple();
}
}

B5RunAction::~B5RunAction()

{
delete G4AnalysisManager::Instance();
}

void B5RunAction::BeginOfRunAction(const G4Run* /run/)

{
// Get analysis manager
auto analysisManager = G4AnalysisManager::Instance();
// Open an output file, it can be overwritten in a macro
analysisManager->OpenFile();
}

void B5RunAction::EndOfRunAction(const G4Run* /run/)

{
// save histograms & ntuple
auto analysisManager = G4AnalysisManager::Instance();
analysisManager->Write();
analysisManager->CloseFile();
}

B5EventAction.cc

#include "B5Analysis.hh"

void B5EventAction::EndOfEventAction(const G4Event* event)

{
// Fill histograms & ntuple
// Get analysis manage
auto analysisManager = G4AnalysisManager::Instance();
// Fill histograms
analysisManager->FillH1(0, nhit );
analysisManager->FillH1(1, nhit );
analysisManager->FillH2(0, localPos.x(), localPos.y());
analysisManager->FillH2(1, localPos.x(), localPos.y());
analysisManager->FillNtupleIColumn(0, dHC1->entries());
analysisManager->FillNtupleIColumn(1, dHC1->entries());
analysisManager->FillNtupleDColumn(2, totalEmE);
analysisManager->FillNtupleDColumn(3, totalHadE);
analysisManager->FillNtupleDColumn(4,(hHC1)[i]->GetTime());
analysisManager->FillNtupleDColumn(5,(hHC2)[i]->GetTime());
analysisManager->AddNtupleRow();
}

B5生成的root文件

exampleB4a

B4Analysis.hh

#ifndef B4Analysis_h
#define B4Analysis_h 1
#include "g4root.hh"
//#include "g4xml.hh"
//#include "g4csv.hh"
#endif

B4RunAction.cc

#include "B4Analysis.hh"

B4RunAction::B4RunAction(): G4UserRunAction()

{
auto analysisManager = G4AnalysisManager::Instance();
G4cout << "Using " << analysisManager->GetType() << G4endl;
analysisManager->SetVerboseLevel(1);
analysisManager->SetNtupleMerging(true);
analysisManager->CreateH1("Eabs","Edep in absorber", 100, 0., 800MeV);
analysisManager->CreateH1("Egap","Edep in gap", 100, 0., 100MeV);
analysisManager->CreateH1("Labs","trackL in absorber", 100, 0., 1m);
analysisManager->CreateH1("Lgap","trackL in gap", 100, 0., 50cm);
analysisManager->CreateNtuple("B4", "Edep and TrackL");
analysisManager->CreateNtupleDColumn("Eabs");
analysisManager->CreateNtupleDColumn("Egap");
analysisManager->CreateNtupleDColumn("Labs");
analysisManager->CreateNtupleDColumn("Lgap");
analysisManager->FinishNtuple();
}

B4RunAction::~B4RunAction()

{
delete G4AnalysisManager::Instance();
}

void B4RunAction::BeginOfRunAction(const G4Run* /run/)

{
auto analysisManager = G4AnalysisManager::Instance();
G4String fileName = "B4";
analysisManager->OpenFile(fileName);
}

void B4RunAction::EndOfRunAction(const G4Run* /run/)

{
auto analysisManager = G4AnalysisManager::Instance();
G4cout << " EAbs : mean = "
<< G4BestUnit(analysisManager->GetH1(0)->mean(), "Energy")
<< " rms = "
<< G4BestUnit(analysisManager->GetH1(0)->rms(), "Energy") << G4endl;
analysisManager->Write();
analysisManager->CloseFile();
}

B4aEventAction.cc

void B4aEventAction::EndOfEventAction(const G4Event* event)

{
// fill histograms
analysisManager->FillH1(0, fEnergyAbs);
analysisManager->FillH1(1, fEnergyGap);
analysisManager->FillH1(2, fTrackLAbs);
analysisManager->FillH1(3, fTrackLGap);
// fill ntuple
analysisManager->FillNtupleDColumn(0, fEnergyAbs);
analysisManager->FillNtupleDColumn(1, fEnergyGap);
analysisManager->FillNtupleDColumn(2, fTrackLAbs);
analysisManager->FillNtupleDColumn(3, fTrackLGap);
analysisManager->AddNtupleRow();
}

Test run2.mac

exampleB4a -m run2.mac

## 50 MeV e-
/run/beamOn 1000

B4a生成的root文件

1. Ntuple

2. Histgram

以B4a的运行run2.mac得到的B4.root为例，怎样出图？

在../B4a/目录下 root B4.root [ ] B4->Draw("Eabs") 便可得到：

/* 这里的单位重新采用了keV */

2. exampleB1的root读取测试

B1RunAction.cc

#include "B1Analysis.hh" B1RunAction::B1RunAction(): G4UserRunAction()

{
auto analysisManager = G4AnalysisManager::Instance();
G4cout << "Using " << analysisManager->GetType() << G4endl;
analysisManager->SetNtupleMerging(true);
analysisManager->SetVerboseLevel(1); analysisManager->SetFileName("lxzb1");
analysisManager->CreateH1("Eabs","Edep in absorber", 100, 0., 200*MeV);
}

B1RunAction::~B1RunAction()

{
delete G4AnalysisManager::Instance();
}

void B1RunAction::BeginOfRunAction(const G4Run*)

{
auto analysisManager = G4AnalysisManager::Instance();
analysisManager->OpenFile();
}

void B1RunAction::EndOfRunAction(const G4Run* run)

{
auto analysisManager = G4AnalysisManager::Instance();
analysisManager->Write();
analysisManager->CloseFile();
}

B1EventAction.cc

void B1EventAction::EndOfEventAction(const G4Event*)

{
auto analysisManager = G4AnalysisManager::Instance();
analysisManager->FillH1(0, fEdep);
}

B1Analysis.hh

#ifndef B1Analysis_h
#define B1Analysis_h 1
#include "g4root.hh"
#endif

Test run1.mac

exampleB1 run1.mac

/gun/particle neutron
/gun/energy 210 MeV
/analysis/setFileName B1_neutron
/run/beamOn 100

B1生成的root文件

3. 分析Histgram与NTuple的存储特点

1. 相比Histgram，怎样理解NTuple？

先来看一下root手册中对“多元组”NTuple的用法介绍，看小节 7.2 N-tuples in ROOT " http://web.mit.edu/root_v6.12/ROOT-Primer.pdf " ;
NTuple好处在于保留存储数据的真实值，不会像Hist一样丢掉数据点的具体数值（只能反应出Bin数值，即数据失真）；
一个NTuple变量可以同时存储多个变量，然而为了清晰，将一个NTuple变量劈成（Branch）了很多Columns，于是有了CreateNtupleDColumn("var")，其中的D表示存储数值为double型（故I表示int型）；
当以一维直方图显示单个NTuple变量时，hist的bin数默认为100个。

2. NTuple变量的神奇之处：

在EventAction中填充数据时，每个Event下均会填充，if ==0 时填充数值0；
NTuple中存储的变量可以按照EventID的顺序还原存入数值，它是个适应于多线程中保存所有数据的强大容器。

加事例判选后的hist与tuple对比（B4a为例，可以发现NTuple变量会强行填充数据）：

图6. 加事例储存判选后Hist与NTuple的存储对比

**4. 给出后续*.root文件读取的代码示范**

1.直接读取直方图

root plothist.C

----plothist.C

{
using namespace std;
gROOT->Reset();
TFile* infile = TFile::Open("B4.root");
TCanvas* c1 = new TCanvas("c1", " ");
TH1D* hist1 = (TH1D*)infile->Get("Eabs");
hist1->Draw("HIST");
}

2. 从直方图目录中读取直方图

root plothist.C

----plothist.C

{
using namespace std;
gROOT->Reset();
TFile f = TFile("AnaEx01.root");
TCanvas* c1 = new TCanvas("c1", " ");
TDirectory* dir = f.Get("histo");
TH1D* hist1 = (TH1D*)dir->Get("1");
hist1->Draw("HIST");
}

3. 从多元组NTuple变量中读取不同Column中的数据

root plothist.C

----plothist.C

{
using namespace std;
gROOT->Reset();
TChain in_chain("B4"); /* 其中"B4"为创建的NTuple变量名 */
in_chain.Add("B4.root");
double Edep;
in_chain.SetBranchAddress("Eabs", &Edep); /* 将Column"Eabs"传递给Edep变量 */
for (int irow=0;irow<in_chain.GetEntries();++irow)
{
in_chain.GetEntry(irow);
cout<<Edep<<endl;
//将所有Event下保存的Edep输出，不满足存储条件的事例将输出0，见图6右下子图
}
}

5. 服务器后台qsub提交作业

1.提交单次作业：qsub oneJob.sh

----oneJob.sh

#PBS -N lxzjob
#PBS -l nodes=1:ppn=1
exampleB1 /data/simulation/home/liangxz/geant4_workdir/B1/run2.mac > log

2.提交多个作业：./multijob.sh (预先 chmod +x multijob.sh)

1) ----oneJob.sh

#PBS -N lxzjob  
#PBS -l nodes=1:ppn=1  
exampleB1  /data/simulation/home/liangxz/geant4_workdir/B1/${NAME} > log_${NAME}  
#NAME表示选择的mac文件

2) ----multihob.sh

#!/bin/bash  
# i 代表要申请多少核心  
i=1  
while(( \${i}<=5 ))
do   
   qsub -v NAME="run${i}.mac" oneJob.sh  
   let i++  
done

喜欢的话，分享一下吧~^o^~

本文参与腾讯云自媒体分享计划，分享自微信公众号。

原始发表：2020-06-07，如有侵权请联系 cloudcommunity@tencent.com 删除

bash

本文分享自人芳觅微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体分享计划，欢迎热爱写作的你一起参与！

登录后参与评论

0 条评论

热度