GEANT4中AIDA的调用

鸿飞冥冥

Geant4 学习心得—AIDA的使用 I

AIDA(Abstract Interfaces for Data Analysis)主要是提供Hbook，ROOT，或者XML格式输出。在学习的过程中我们看到好几种不同的调用方式，在这里我们逐个总结出来。

extended\electromagnetic\TestEM4例子的AIDA调用方式：
（1）在RUNACTION.CC中定义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）在EVENTACTION.CC中累积谱

#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的办法。

转自www.52mc.net

鸿飞冥冥

Geant4 学习心得—AIDA的使用 II

AIDA(Abstract Interfaces for Data Analysis)主要是提供Hbook，ROOT，或者XML格式输出。在学习的过程中有几种不同的调用方式。

extended\hadronic\Hadr01例子的AIDA调用方式：

（1）在RUNACTION.CC中调用HistoManager管理器
void RunAction::EndOfRunAction(const G4Run*)
{
  HistoManager::GetPointer()->EndOfRun();   //获得Histo的指针
}

（2）在HistoManger.CC中定义HisyoManger的子函数等
HistoManager* HistoManager::fManager = 0;

HistoManager* HistoManager::GetPointer()
{
  if(!fManager) {
    static HistoManager manager;
    fManager = &manager;
  }
  return fManager;
}

（3）在Histo.CC中定义Hbook文件等
#include "Histo.hh"
#ifdef G4ANALYSIS_USE
#include <AIDA/AIDA.h>
#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<nHisto; 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<nHisto; 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<nHisto) {
    bins = nb;
    xmin = x1;
    xmax = x2;
    unit = u;
  }
}

void Histo::fill(G4int i, G4double x, G4double w)
{
#ifdef G4ANALYSIS_USE   
  if(i>=0 && i<nHisto) {
    histo->fill(x/unit, w);
  }  
#endif
}


void Histo::scale(G4int i, G4double x)
{
#ifdef G4ANALYSIS_USE   
  if(i>=0 && i<nHisto) {
    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<nHisto) {
    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<bins; 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:rocessHits(G4Step* aStep, G4TouchableHistory*)
{
  theHisto->AddTargetStep(aStep);
  return true;
}

CheckVolumeSD::CheckVolumeSD(const G4String& name)
:G4VSensitiveDetector(name)
{
  theHisto = HistoManager::GetPointer();
}

G4bool CheckVolumeSD:rocessHits(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）在MAIN.CC中调用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调用方式：采用AnalysisManager.cc

（1）在AnalysisManager.cc中定义所有的Histogram，Tuple 以及 Plotter等
#ifdef G4ANALYSIS_USE
#include "G4ios.hh"
#include "G4SDManager.hh"
#include "G4Run.hh"
#include "G4Event.hh"
#include "G4HCofThisEvent.hh"
#include <AIDA/IAnalysisFactory.h>
#include <AIDA/ITreeFactory.h>
#include <AIDA/ITupleFactory.h>
#include <AIDA/IHistogramFactory.h>
#include <AIDA/ITree.h>
#include <AIDA/IHistogram1D.h>
#include <AIDA/ITuple.h>
#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.<format> file.   

  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<n_hit;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）在DetectorConstruction.cc中定义转动角度，在G4VPlacement中的第一项设置转动。
  armRotation = new G4RotationMatrix();
  armRotation->rotateY(armAngle);

// second arm
  secondArmPhys
    = new G4PVPlacement(armRotation,G4ThreeVector(x,0.,z),secondArmLogical,
                        "secondArmPhys",worldLogical,0,0);

（3）A01中DetectorConstruction.cc的设置磁场转动
// placement of Tube  
  G4RotationMatrix* fieldRot = new G4RotationMatrix();
  fieldRot->rotateX(90.*deg);
  new G4PVPlacement(fieldRot,G4ThreeVector(),magneticLogical,
                    "magneticPhysical",worldLogical,0,0);

(4) 在ExN04的DetectorConstruction.cc中设置的转动

  G4RotationMatrix rm;
    rm.rotateZ(phi);
    muoncounter_phys = new G4PVPlacement(G4Transform3D(rm,G4ThreeVector(x,y,z)),
                          muoncounter_log, "muoncounter_P",
                          experimentalHall_log,false,i);

鸿飞冥冥

好的，贴完了，

根据网上CNSCOTT网友的博文格式写的，在此表示感谢，呵呵