Geant4 学习心得—AIDA的使用I
AIDA(Abstract Interfaces for Data Analysis)主要是提供Hbook,ROOT,或者XML格式输出。在学习的过程中我们看到好几种不同的调用方式,在这里我们逐个总结出来。
extended\electromagnetic\TestEM4例子的AIDA调用方式:
(1)在https://www.doczj.com/doc/6d14515248.html,中定义hbook 一维谱(或者root,XML格式的文件)
#ifdef G4ANALYSIS_USE
#include "AIDA/AIDA.h" //需要的头文件
#endif
RunAction::RunAction():af(0), tree(0)
{
histo[0] = 0;
#ifdef G4ANALYSIS_USE
// Creating the analysis factory
af = AIDA_createAnalysisFactory();
if (af) {
// Creating the tree factory
AIDA::ITreeFactory* tf = af->createTreeFactory();
// Creating a tree mapped to an hbook file.
G4bool readOnly = false;
G4bool createNew = true;
G4String options = "--noErrors uncompress";
tree = tf->create("TestEm4.hbook","hbook",readOnly,createNew, options);
//tree = tf->create("TestEm4.root","root",readOnly,createNew, options);
//tree = tf->create("TestEm4.XML" ,"XML" ,readOnly,createNew, options);
delete tf;
if (tree) {
// Creating a histogram factory
AIDA::IHistogramFactory* hf = af->createHistogramFactory(*tree);
// Creating the histogram
histo[0]=hf->createHistogram1D // 头文件中定义了histo[1]
("1","total energy deposit in HpGe Detector(keV)",3000,1.,3000);
// energy ragin 1-5000keV
delete hf;
}
}
#endif
}
RunAction::~RunAction()
#ifdef G4ANALYSIS_USE
tree->commit(); // Writing the histograms to the file
tree->close(); // and closing the tree (and the file)
delete tree;
delete af;
#endif
}
(2)在https://www.doczj.com/doc/6d14515248.html,中累积谱
#ifdef G4ANALYSIS_USE
#include "AIDA/IHistogram1D.h"
#endif
void EventAction::BeginOfEventAction( const G4Event* evt)
{
G4int evtNb = evt->GetEventID();
//additional initializations
TotalEnergyDeposit = 0.; // 初始化能量变量
}
void EventAction::EndOfEventAction( const G4Event* evt)
{
#ifdef G4ANALYSIS_USE
Run->GetHisto(0)->fill(TotalEnergyDeposit/MeV);
#endif
}
总结:这是最简单的定义一维Histogram的办法。
Geant4 学习心得—AIDA的使用II
AIDA(Abstract Interfaces for Data Analysis)主要是提供Hbook,ROOT,或者XML格式输出。在学习的过程中有几种不同的调用方式。
extended\hadronic\Hadr01例子的AIDA调用方式:
(1)在https://www.doczj.com/doc/6d14515248.html,中调用HistoManager管理器
void RunAction::EndOfRunAction(const G4Run*)
{
HistoManager::GetPointer()->EndOfRun(); //获得Histo的指针}
(2)在https://www.doczj.com/doc/6d14515248.html,中定义HisyoManger的子函数等HistoManager* HistoManager::fManager = 0;
HistoManager* HistoManager::GetPointer()
{
if(!fManager) {
static HistoManager manager;
fManager = &manager;
}
return fManager;
}
(3)在https://www.doczj.com/doc/6d14515248.html,中定义Hbook文件等
#include "Histo.hh"
#ifdef G4ANALYSIS_USE
#include
#include "HistoMessenger.hh"
#endif
Histo::Histo(G4int ver)
{
verbose = ver;
histName = "histo";
histType = "hbook";
option = "--noErrors uncompress";
nHisto = 0;
defaultAct = 1;
tupleName = "tuple.hbook";
tupleId = "100";
tupleList = "";
ntup = 0;
messenger = 0;
#ifdef G4ANALYSIS_USE
messenger = new HistoMessenger(this); // messange 指针tree = 0;
af = 0;
#endif
}
Histo::~Histo()
{
#ifdef G4ANALYSIS_USE
delete messenger;
delete af;
#endif
}
void Histo::book()
{
#ifdef G4ANALYSIS_USE
// Creating the analysis factory
if(!af) af = AIDA_createAnalysisFactory();
// Creating the tree factory
AIDA::ITreeFactory* tf = af->createTreeFactory();
// Creating a tree mapped to a new hbook file.
G4String name = histDir + histName + histExt + "." + histType;
tree = tf->create(name, histType, false, true, option);
delete tf;
// Creating a histogram factory, whose histograms will be handled by the tree AIDA::IHistogramFactory* hf = af->createHistogramFactory( *tree );
// Creating an 1-dimensional histograms in the root directory of the tree
for(G4int i=0; i if(active) { G4String idd; if(histType == "root") idd = "h" + ids; else idd = ids; histo = hf->createHistogram1D(idd, tittles, bins, xmin, xmax); } else { histo = 0; } } delete hf; // Creating a tuple factory, whose tuples will be handled by the tree if(tupleList != "") { AIDA::ITupleFactory* tpf = af->createTupleFactory( *tree ); ntup = tpf->create(tupleId, tupleName, tupleList); delete tpf; } #endif } void Histo::save() { #ifdef G4ANALYSIS_USE // Write histogram file tree->commit(); tree->close(); delete tree; tree = 0; #endif } void Histo::reset() { #ifdef G4ANALYSIS_USE delete tree; tree = 0; #endif } void Histo::setFileType(const G4String& nam) { if(nam == "hbook" || nam == "root" || nam == "aida") histType = nam; else if(nam == "XML" || nam == "xml") histType = "aida"; } void Histo::add1D(const G4String& id, const G4String& name, G4int nb, G4double x1, G4double x2, G4double u) { if(nHisto > 0) { for(G4int i=0; i if(ids == id) return; } } nHisto++; x1 /= u; x2 /= u; active.push_back(defaultAct); bins.push_back(nb); xmin.push_back(x1); xmax.push_back(x2); unit.push_back(u); ids.push_back(id); tittles.push_back(name); histo.push_back(0); } void Histo::setHisto1D(G4int i, G4int nb, G4double x1, G4double x2, G4double u) { if(i>=0 && i bins = nb; xmin = x1; xmax = x2; unit = u; } } void Histo::fill(G4int i, G4double x, G4double w) { #ifdef G4ANALYSIS_USE if(i>=0 && i histo->fill(x/unit, w); } #endif } void Histo::scale(G4int i, G4double x) { #ifdef G4ANALYSIS_USE if(i>=0 && i histo->scale(x); } #endif } void Histo::addTuple(const G4String& w1, const G4String& w2, const G4String& w3) { tupleId = w1; tupleName = w2; tupleList = w3; } void Histo::fillTuple(const G4String& parname, G4double x) { #ifdef G4ANALYSIS_USE if(ntup) ntup->fill(ntup->findColumn(parname), (float)x); #endif } void Histo::addRow() { #ifdef G4ANALYSIS_USE if(ntup) ntup->addRow(); #endif } void Histo::print(G4int i) { #ifdef G4ANALYSIS_USE if(i>=0 && i G4double step = (xmax - xmin)/G4double( bins); G4double x = xmin - step*0.5; G4double y, maxX=0, maxY=0; G4int maxJ=0; for(G4int j=0; j x += step; y = histo->binHeight(j); if(maxY < y) {maxY = y; maxX = x; maxJ = j;} } } #endif } (4)在RUNACTION中开始和结束累谱。 void RunAction::BeginOfRunAction(const G4Run* aRun) { G4int id = aRun->GetRunID(); (HistoManager::GetPointer())->BeginOfRun(); void RunAction::EndOfRunAction(const G4Run*) { HistoManager::GetPointer()->EndOfRun(); //开始累谱 } (5)在灵敏探测器的定义中 TargetSD::TargetSD(const G4String& name) :G4VSensitiveDetector(name) { theHisto = HistoManager::GetPointer(); } G4bool TargetSD::ProcessHits(G4Step* aStep, G4TouchableHistory*) { theHisto->AddTargetStep(aStep); return true; } CheckVolumeSD::CheckVolumeSD(const G4String& name) :G4VSensitiveDetector(name) { theHisto = HistoManager::GetPointer(); } G4bool CheckVolumeSD::ProcessHits(G4Step* aStep, G4TouchableHistory*) { const G4Track* track = aStep->GetTrack(); if(track->GetTrackID() > 1) theHisto->AddLeakingParticle(track); return true; Geant4 学习心得—AIDA的使用III AIDA(Abstract Interfaces for Data Analysis)主要是提供Hbook,ROOT,或者XML格式输出。在学习的过程中有几种不同的调用方式。 extended\biasing\B02例子的AIDA调用方式: (1)在https://www.doczj.com/doc/6d14515248.html,中调用AIDA // AIDA stuff #include "AIDA/IAnalysisFactory.h" #include "AIDA/ITree.h" #include "AIDA/ITreeFactory.h" #include "AIDA/IHistogram1D.h" #include "AIDA/IHistogramFactory.h" int main(int , char **) { // create a histogram for a special scorer for the last cell AIDA::IAnalysisFactory *af = AIDA_createAnalysisFactory(); AIDA::ITreeFactory *tf = af->createTreeFactory(); AIDA::ITree *tree = tf->create("b02.hbook", "hbook",false,true); AIDA::IHistogramFactory *hf = af->createHistogramFactory( *tree ); AIDA::IHistogram1D *h = hf->createHistogram1D("10","w*sl vs. e", 30, 0., 20*MeV); tree->commit(); tree->close(); delete af; delete tf; delete tree; delete runManager; return 0; } (2)在Score中调用Histo,在hh头文件中定义了fHisto B02CellScorer::B02CellScorer(AIDA::IHistogram1D *h):fHisto(h){} void B02CellScorer::ScoreAnExitingStep(const G4Step &aStep, const G4GeometryCell &pre_gCell) { fG4CellScorer.ScoreAnExitingStep(aStep, pre_gCell); FillHisto(aStep); } void B02CellScorer::FillHisto(const G4Step &aStep) { G4StepPoint *p = 0; p = aStep.GetPreStepPoint(); G4double sl = aStep.GetStepLength(); G4double slw = sl * p->GetWeight(); fHisto->fill(p->GetKineticEnergy(), slw); } Geant4 学习心得—AIDA的使用IV AIDA(Abstract Interfaces for Data Analysis)主要是提供Hbook,ROOT,或者XML格式输出。在学习的过程中有几种不同的调用方式。 extended\Analysis\AnaEx01例子的AIDA调用方式:采用https://www.doczj.com/doc/6d14515248.html, (1)在https://www.doczj.com/doc/6d14515248.html,中定义所有的Histogram,Tuple 以及Plotter等 #ifdef G4ANALYSIS_USE #include "G4ios.hh" #include "G4SDManager.hh" #include "G4Run.hh" #include "G4Event.hh" #include "G4HCofThisEvent.hh" #include #include #include #include #include #include #include #include "AnaEx01CalorHit.hh" #include "AnaEx01AnalysisManager.hh" AnaEx01AnalysisManager::AnaEx01AnalysisManager(AIDA::IAnalysisFactory* aAIDA) :fCalorimeterCollID(-1),fAIDA(aAIDA),fTree(0),fEAbs(0),fLAbs(0),fEGap(0), fLGap(0),fTuple(0) { // Could fail if no AIDA implementation found : if(!fAIDA) { G4cout << "AIDA analysis factory not found." << G4endl; return; } AIDA::ITreeFactory* treeFactory = fAIDA->createTreeFactory(); if(!treeFactory) return; // Create a tree-like container to handle histograms. // This tree is associated to a AnaEx01. std::string h_name_EAbs("EAbs"); std::string h_name_LAbs("LAbs"); std::string h_name_EGap("EGap"); std::string h_name_LGap("LGap"); std::string t_name("AnaEx01"); // File format : std::string format("xml"); //std::string format("hbook"); //std::string format("root"); std::string file("AnaEx01"); std::string ext; ext = "."+format; std::string opts; if(format=="hbook") { opts = "compress=no"; h_name_EAbs = "1"; h_name_LAbs = "2"; h_name_EGap = "3"; h_name_LGap = "4"; t_name = "101"; } else if(format=="root") { opts = "compress=yes"; } else if(format=="xml") { ext = ".aida"; opts = "compress=no"; } else { G4cout << "storage format \"" << format << "\"" << " not handled in this example." << G4endl; return; } file += ext; fTree = treeFactory->create(file,format,false,true,opts); // Factories are not "managed" by an AIDA analysis system. // They must be deleted by the AIDA user code. delete treeFactory; if(!fTree) { G4cout << "can't create tree associated to file \"" << file << "\"." << G4endl; return; } fTree->mkdir("histograms"); fTree->cd("histograms"); // Create an histo factory that will create histo in the tree : AIDA::IHistogramFactory* histoFactory = fAIDA->createHistogramFactory(*fTree); if(histoFactory) { fEAbs = histoFactory->createHistogram1D(h_name_EAbs,"EAbs",100,0,100); if(!fEAbs) G4cout << "can't create histo EAbs." << G4endl; fLAbs = histoFactory->createHistogram1D(h_name_LAbs,"LAbs",100,0,100); if(!fLAbs) G4cout << "can't create histo LAbs." << G4endl; fEGap = histoFactory->createHistogram1D(h_name_EGap,"EGap",100,0,10); if(!fEGap) G4cout << "can't create histo EGap." << G4endl; fLGap = histoFactory->createHistogram1D(h_name_LGap,"LGap",100,0,100); if(!fLGap) G4cout << "can't create histo LGap." << G4endl; delete histoFactory; } fTree->cd(".."); fTree->mkdir("tuples"); fTree->cd("tuples"); // Get a tuple factory : AIDA::ITupleFactory* tupleFactory = fAIDA->createTupleFactory(*fTree); if(tupleFactory) { // Create a tuple : fTuple = tupleFactory->create(t_name,"AnaEx01", "double EAbs,double LAbs,double EGap,double LGap"); if(!fTuple) G4cout << "can't create tuple." << G4endl; delete tupleFactory; } fTree->cd(".."); } AnaEx01AnalysisManager::~AnaEx01AnalysisManager() { } void AnaEx01AnalysisManager::BeginOfRun(const G4Run* aRun){ G4cout << "### Run " << aRun->GetRunID() << " start." << G4endl; } void AnaEx01AnalysisManager::EndOfRun(const G4Run*){ if(fTree) fTree->commit(); if(fEAbs) { G4cout << "Histo : EAbs : mean " << fEAbs->mean() << " rms : " << fEAbs->rms() << G4endl; G4cout << "Histo : LAbs : mean " << fLAbs->mean() << " rms : " << fLAbs->rms() << G4endl; G4cout << "Histo : EGap : mean " << fEGap->mean() << " rms : " << fEGap->rms() << G4endl; G4cout << "Histo : LGap : mean " << fLGap->mean() << " rms : " << fLGap->rms() << G4endl; } } void AnaEx01AnalysisManager::BeginOfEvent(const G4Event*){ if(fCalorimeterCollID==-1) { G4SDManager* SDman = G4SDManager::GetSDMpointer(); fCalorimeterCollID = SDman->GetCollectionID("CalCollection"); } } void AnaEx01AnalysisManager::EndOfEvent(const G4Event* aEvent){ if(!fEAbs) return; // No histo booked ! if(!fTuple) return; // No tuple booked ! //G4int evtNb = aEvent->GetEventID(); G4HCofThisEvent* HCE = aEvent->GetHCofThisEvent(); AnaEx01CalorHitsCollection* CHC = HCE ? (AnaEx01CalorHitsCollection*)(HCE->GetHC(fCalorimeterCollID)) : 0; if(CHC) { G4int n_hit = CHC->entries(); for (G4int i=0;i G4double EAbs = (*CHC)->GetEdepAbs(); G4double LAbs = (*CHC)->GetTrakAbs(); G4double EGap = (*CHC)->GetEdepGap(); G4double LGap = (*CHC)->GetTrakGap(); fEAbs->fill(EAbs); fLAbs->fill(LAbs); fEGap->fill(EGap); fLGap->fill(LGap); fTuple->fill(0,EAbs); fTuple->fill(1,LAbs); fTuple->fill(2,EGap); fTuple->fill(3,LGap); fTuple->addRow(); } } } void AnaEx01AnalysisManager::Step(const G4Step*){} #endif (2)头文件定义为: #ifndef AnaEx01AnalysisManager_h #define AnaEx01AnalysisManager_h 1 #ifdef G4ANALYSIS_USE class G4Run; class G4Event; class G4Step; namespace AIDA { class IAnalysisFactory; class ITree; class IHistogram1D; class ITuple; } class AnaEx01AnalysisManager { public: AnaEx01AnalysisManager(AIDA::IAnalysisFactory*); virtual ~AnaEx01AnalysisManager(); public: virtual void BeginOfRun(const G4Run*); virtual void EndOfRun(const G4Run*); virtual void BeginOfEvent(const G4Event*); virtual void EndOfEvent(const G4Event*); virtual void Step(const G4Step*); private: int fCalorimeterCollID; AIDA::IAnalysisFactory* fAIDA; AIDA::ITree* fTree; AIDA::IHistogram1D* fEAbs; AIDA::IHistogram1D* fLAbs; AIDA::IHistogram1D* fEGap; AIDA::IHistogram1D* fLGap; AIDA::ITuple* fTuple; }; #else class AnaEx01AnalysisManager; #endif #endif (3)在RunAction中调用,填谱 #ifdef G4ANALYSIS_USE #include "AnaEx01AnalysisManager.hh" #endif #include "AnaEx01RunAction.hh" AnaEx01RunAction::AnaEx01RunAction( AnaEx01AnalysisManager* aAnalysisManager ):fAnalysisManager(aAnalysisManager){} AnaEx01RunAction::~AnaEx01RunAction(){} void AnaEx01RunAction::BeginOfRunAction(const G4Run* aRun) { #ifdef G4ANALYSIS_USE if(fAnalysisManager) fAnalysisManager->BeginOfRun(aRun); #endif } void AnaEx01RunAction::EndOfRunAction(const G4Run* aRun){ #ifdef G4ANALYSIS_USE if(fAnalysisManager) fAnalysisManager->EndOfRun(aRun); #endif } Geant4 学习心得—如何转动实体 我们通过A01例子来学习如何设置实体的转动 (1) 在DetectorConstructtion.hh中定义RotationMatrix: #include "G4RotationMatrix.hh" class A01DetectorConstruction : public G4VUserDetectorConstruction { private: G4RotationMatrix* armRotation; } (2)在https://www.doczj.com/doc/6d14515248.html,中定义转动角度,在G4VPlacement中的第一项设置转动。armRotation = new G4RotationMatrix(); armRotation->rotateY(armAngle); // second arm secondArmPhys = new G4PVPlacement(armRotation,G4ThreeVector(x,0.,z),secondArmLogical, "secondArmPhys",worldLogical,0,0); (3)A01中https://www.doczj.com/doc/6d14515248.html,的设置磁场转动 // placement of Tube G4RotationMatrix* fieldRot = new G4RotationMatrix(); fieldRot->rotateX(90.*deg); new G4PVPlacement(fieldRot,G4ThreeVector(),magneticLogical, "magneticPhysical",worldLogical,0,0); (4) 在ExN04的https://www.doczj.com/doc/6d14515248.html,中设置的转动 G4RotationMatrix rm; rm.rotateZ(phi); muoncounter_phys = new G4PVPlacement(G4Transform3D(rm,G4ThreeVector(x,y,z)), muoncounter_log, "muoncounter_P", experimentalHall_log,false,i);