1. exampleB5与exampleB4a中与root数据存储有关的代码
#ifndef B5Analysis_h
#define B5Analysis_h 1
#include "g4root.hh"
//#include "g4xml.hh"
//#include "g4csv.hh"
#endif
#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();
}
#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文件
#ifndef B4Analysis_h
#define B4Analysis_h 1
#include "g4root.hh"
//#include "g4xml.hh"
//#include "g4csv.hh"
#endif
#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();
}
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();
}
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 */
#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();
}
void B1EventAction::EndOfEventAction(const G4Event*)
{
auto analysisManager = G4AnalysisManager::Instance();
analysisManager->FillH1(0, fEdep);
}
#ifndef B1Analysis_h
#define B1Analysis_h 1
#include "g4root.hh"
#endif
exampleB1 run1.mac
/gun/particle neutron
/gun/energy 210 MeV
/analysis/setFileName B1_neutron
/run/beamOn 100
B1生成的root文件
1. 相比Histgram,怎样理解NTuple?
2. NTuple变量的神奇之处:
加事例判选后的hist与tuple对比(B4a为例,可以发现NTuple变量会强行填充数据):
图6. 加事例储存判选后Hist与NTuple的存储对比
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");
}
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");
}
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右下子图
}
}
----oneJob.sh
#PBS -N lxzjob
#PBS -l nodes=1:ppn=1
exampleB1 /data/simulation/home/liangxz/geant4_workdir/B1/run2.mac > log
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^~