Weihe Zeng
Tsinghua University
Dec. 18, 2017

Content

Review: Passing information from Step -> Track -> Event -> Run
Sensitive Detector
STL Containers and File I/O
- STL containers
- standard output / file stream (binary/formatted)
- Analysis Tool g4tools -> csv, ROOT, XML
ROOT basics: histogram, tree (n-tuple), function and fitting
Execute a Geant4 Program: Hard-coded v.s. Macro files.

Traditional Data Flow

120%

User Actions (passing information)

G4UserSteppingAction

// Get Current track information and "delta" data
virtual void UserSteppingAction(const G4Step*);

G4UserTrackingAction

virtual void PreUserTrackingAction(const G4Track*);  
virtual void PostUserTrackingAction(const G4Track*);

G4UserEventAction

virtual void BeginOfEventAction(const G4Event*);  
virtual void EndOfEventAction(const G4Event*);

G4UserRunAction

virtual void BeginOfRunAction(const G4Run*);  
virtual void EndOfRunAction(const G4Run*);

Sensitive Detector and Hits

Recall how we ensure that each step is within our detector:

By geometry position (based on World's coordinate);
Get volume (pointer,r name and copy#) via TouchableHandle and compare with our detector;

auto curPosi = step->GetPreStepPoint()->GetPosition();  
G4bool flagInDet = isInDetector(curPosi);  
auto volume = step->GetPreStepPoint()->GetTouchableHandle()->GetVolume();  
if ( volume == detPhysicalVolume ) {...}

Can we define some "sensitive" geometry as "detector" during Geometry Construction?

User actions: full control and straightforward, but is not a ideal "OO" way;
Sensitive Detector (SD): useful built-in utilities, clear definition of each objects, but with more code.

Components of SD

To design a sensitive detector, users should complete two abstract classes and assign SD to targeting logical volume.

Hit G4VHit: a snapshot of the physical interaction of a track in the sensitive region of a detector. In it you can store information associated with a G4Step object.
- the position and time of the step,
- the momentum and energy of the track,
- the energy deposition of the step,
- geometrical information
SD G4VSensitiveDetector: represents a detector and constructs of hit objects using information from steps along a particle track.
Defining a SD for one geometry (logical volume)
- G4VUserDetectorConstruction::ConstructSDandField
- Each SD object must have a unique name

Hits and hit collection

class B4cCalorHit : public G4VHit{  
  public:
    B4cCalorHit();
    B4cCalorHit(const B4cCalorHit&);
    virtual ~B4cCalorHit();
    // operators
    const B4cCalorHit& operator=(const B4cCalorHit&);
    G4int operator==(const B4cCalorHit&) const;
    inline void* operator new(size_t);
    inline void  operator delete(void*);
    // methods from base class
    virtual void Draw() {}
    virtual void Print();
    // methods to handle data
    void Add(G4double de, G4double dl);
    // get methods
    G4double GetEdep() const;
    G4double GetTrackLength() const;     
  private:
    G4double fEdep;        //Energy deposit in the sensitive volume
    G4double fTrackLength; //Track length in the sensitive volume
};
using B4cCalorHitsCollection = G4THitsCollection<B4cCalorHit>;

Sensitive Detector

class B4cCalorimeterSD : public G4VSensitiveDetector{  
  public:
    B4cCalorimeterSD(const G4String& name, 
                     const G4String& hitsCollectionName, 
                     G4int nofCells);
    virtual ~B4cCalorimeterSD();  

    // methods from base class
    virtual void   Initialize(
        G4HCofThisEvent* hitCollection);
    virtual G4bool ProcessHits(
        G4Step* step, G4TouchableHistory* history);
    virtual void EndOfEvent(
        G4HCofThisEvent* hitCollection);

  private:
    B4cCalorHitsCollection* fHitsCollection;
    G4int  fNofCells;
};

void B4cCalorimeterSD::Initialize(G4HCofThisEvent* hce){  
  // Create hits collection
  fHitsCollection = new B4cCalorHitsCollection(
      SensitiveDetectorName, collectionName[0]); 

  // Add this collection in hce
  auto hcID = G4SDManager::GetSDMpointer()
      ->GetCollectionID(collectionName[0]);
  hce->AddHitsCollection( hcID, fHitsCollection ); 

  // Create hits
  for (G4int i=0; i<fNofCells+1; i++ ) 
    fHitsCollection->insert(new B4cCalorHit()); 
}

G4bool B4cCalorimeterSD::ProcessHits(  
   G4Step* step, G4TouchableHistory*){  
  // energy deposit
  auto edep = step->GetTotalEnergyDeposit();  
  // step length
  G4double stepLength = step->GetStepLength();
  // Get calorimeter cell id 
  auto touchable = step->GetPreStepPoint()->GetTouchable();    
  auto layerNumber = touchable->GetReplicaNumber(1);

  // Get hit accounting data for this cell
  auto hit = (*fHitsCollection)[layerNumber];
  // Get hit for total accounting
  auto hitTotal 
    = (*fHitsCollection)[fHitsCollection->entries()-1];
  // Add values
  hit->Add(edep, stepLength);
  hitTotal->Add(edep, stepLength);       
  return true;
}

Associate to logical volume

void B4cDetectorConstruction::ConstructSDandField(){  
  // Sensitive detectors
  auto absoSD 
    = new B4cCalorimeterSD("AbsorberSD", 
        "AbsorberHitsCollection", fNofLayers);
  G4SDManager::GetSDMpointer()->AddNewDetector(absoSD);
  SetSensitiveDetector("AbsoLV",absoSD);

  auto gapSD 
    = new B4cCalorimeterSD("GapSD", 
        "GapHitsCollection", fNofLayers);
  G4SDManager::GetSDMpointer()->AddNewDetector(gapSD);
  SetSensitiveDetector("GapLV",gapSD);
}

// method of G4VUserDetectorConstruction
void G4VUserDetectorConstruction::SetSensitiveDetector  
   (const G4String& logVolName, 
    G4VSensitiveDetector* aSD, 
    G4bool multi=false)

Alternative way (NOTE: absorberLV and gapLV should be declared as private member of B4cDetectorConstruction)

// using method of G4VLogicalVolume
absorberLV->SetSensitiveDetector(absoSD);  
gapLV->SetSensitiveDetector(gapSD);

Process data in G4UserEventAction

void B4cEventAction::EndOfEventAction(const G4Event* event){  
  // Get hits collections IDs (only once)
  if ( fAbsHCID == -1 ) {
    fAbsHCID 
      = G4SDManager::GetSDMpointer()->GetCollectionID("AbsorberHitsCollection");
    fGapHCID 
      = G4SDManager::GetSDMpointer()->GetCollectionID("GapHitsCollection");  }
  // Get hits collections
  auto absoHC = GetHitsCollection(fAbsHCID, event);
  auto gapHC = GetHitsCollection(fGapHCID, event);
  // Get hit with total values
  auto absoHit = (*absoHC)[absoHC->entries()-1];
  auto gapHit = (*gapHC)[gapHC->entries()-1]; 

  // Print per event (modulo n)
  auto eventID = event->GetEventID();
  auto printModulo = G4RunManager::GetRunManager()->GetPrintProgress();
  if ( ( printModulo > 0 ) && ( eventID % printModulo == 0 ) ) {
    G4cout << "---> End of event: " << eventID << G4endl;
    PrintEventStatistics(
      absoHit->GetEdep(), absoHit->GetTrackLength(),
      gapHit->GetEdep(), gapHit->GetTrackLength());
  }  
}

Scorer (pre-defined SD) in GEANT4

example/basic/B4/B4d

void B4dDetectorConstruction::ConstructSDandField(){  
  // declare Absorber as a MultiFunctionalDetector scorer
  auto absDetector = new G4MultiFunctionalDetector("Absorber");
  G4SDManager::GetSDMpointer()->AddNewDetector(absDetector);

  G4VPrimitiveScorer* primitive;
  primitive = new G4PSEnergyDeposit("Edep");
  absDetector->RegisterPrimitive(primitive);

  primitive = new G4PSTrackLength("TrackLength");
  auto charged = new G4SDChargedFilter("chargedFilter");
  primitive ->SetFilter(charged);

  absDetector->RegisterPrimitive(primitive);  
  SetSensitiveDetector("AbsoLV",absDetector);
  ...
}

Data Container in C++

Standard Template Library (STL) in C++ provides many useful containers that store data.

Sequence containers

Sequence containers implement data structures which can be accessed sequentially.

| name | description | | ---------------------------------------- | ---------------------------------------- | | array(C++11) | static contiguous array (class template) | | vector | dynamic contiguous array (class template) | | deque | double-ended queue (class template) | | forward_list(C++11) | singly-linked list (class template) | | list | doubly-linked list (class template) |

Associative containers

Associative containers implement sorted data structures that can be quickly searched (O(log n) complexity).

| name | description | | ---------------------------------------- | ---------------------------------------- | | set | collection of unique keys, sorted by keys (class template) | | map | collection of key-value pairs, sorted by keys, keys are unique (class template) | | multiset | collection of keys, sorted by keys (class template) | | multimap | collection of key-value pairs, sorted by keys (class template) |

Unordered associative containers

Unordered associative containers implement unsorted (hashed) data structures that can be quickly searched (O(1)amortized, O(n) worst-case complexity).

| name | description | | ---------------------------------------- | ---------------------------------------- | | unorderedset(C++11) | collection of unique keys, hashed by keys | | unorderedmap(C++11) | collection of key-value pairs, hashed by keys, keys are unique (class template) | | unorderedmultiset(C++11) | collection of keys, hashed by keys (class template) | | unorderedmultimap(C++11) | collection of key-value pairs, hashed by keys (class template) |

Container adaptors

Container adaptors provide a different interface for sequential containers.

| name | description | | ---------------------------------------- | ---------------------------------------- | | stack | adapts a container to provide stack (LIFO data structure) (class template) | | queue | adapts a container to provide queue (FIFO data structure) (class template) | | priority_queue | adapts a container to provide priority queue (class template) |

Example: vector usage

#include <iostream>
#include <vector>
// Create a vector containing integers
std::vector<int> v = {7, 5, 16, 8};  
// Add two more integers to vector
v.push_back(25);  
v.push_back(13);

// Iteration Method 1
for(int n : v)  
  std::cout << n << '\n';

// Iteration Method 2
for (std::vector<int>::iterator it = v.begin();  
     it != v.end(); ++it)
  std::cout << *it << '\n';

// Iteration Method 3
for (std::vector<int>::size_type i = 0; i < v.size(); i++)  
  std::cout << v[i] << '\n'; // or v.at(i)

File I/O

150%

Example: ofstream and ifstream (binary/formatted)

#include <iostream>
#include <fstream>
#include <string>

int main() {  
  std::string filename = "Test.b";
  {
    // open mode: default std::ios::out
    std::ofstream ostrm(filename, std::ios::binary);
    double d = 3.14;
    // binary output
    ostrm.write(reinterpret_cast<char*>(&d), sizeof d); 
    // text output
    ostrm << 123 << "abc" << '\n';                      
  }
  // read back
  std::ifstream istrm(filename, std::ios::binary);
  double d;
  istrm.read(reinterpret_cast<char*>(&d), sizeof d);
  int n;
  std::string s;
  istrm >> n >> s;
  std::cout << " read back: " << d << " " << n << " " << s << '\n';
}

Analysis Tools

g4tools provides code to write and read histograms and ntuples in several formats: ROOT, XML AIDA format and CSV (comma-separated values format).

The analysis classes provide a uniform, user-friendly interface to g4tools and hide the differences according to a selected output technology from the user.

histtuple

Public interface: G4AnalysisManager

Create or get instance of G4AnalysisManager:

auto analysisManager = G4AnalysisManager::Instance();

Select output file format (in a separate header, e.g. B4Analysis.hh):

#ifndef B4Analysis_h
#define B4Analysis_h 1

//#include "g4root.hh"
//#include "g4xml.hh"
#include "g4csv.hh"
#endif

xml and csv AnalysisManager Class doesn't include SetNtupleMerging method, which should be comment out in B4RunAction.

File handling (in UserRunAction)

#include "B4Analysis.hh"

void B4RunAction::BeginOfRunAction(const G4Run* run)  
{
  // Get analysis manager
  auto analysisManager = G4AnalysisManager::Instance();

  // Open an output file
  analysisManager->OpenFile("B4");
}

void B4RunAction::EndOfRunAction(const G4Run* aRun)  
{
  // Save histograms
  auto analysisManager = G4AnalysisManager::Instance();
  analysisManager->Write();
  analysisManager->CloseFile();
}

Create and fill Histograms

Basic functions:

G4int CreateH1(const G4String& name, const G4String& title,  
               G4int nbins, G4double xmin, G4double xmax,
               const G4String& unitName = "none",
               const G4String& fcnName = "none",
               const G4String& binSchemeName = "linear");

G4int CreateH1(const G4String& name, const G4String& title,  
               const std::vector<G4double>& edges,
               const G4String& unitName = "none",
               const G4String& fcnName = "none");

G4int CreateH2(const G4String& name, const G4String& title,  
               G4int nxbins, G4double xmin, G4double xmax,
               G4int nybins, G4double ymin, G4double ymax, ...);

G4bool FillH1(G4int id, G4double value,  
              G4double weight = 1.0);
G4bool FillH2(G4int id, G4double xvalue, G4double yvalue,  
              G4double weight = 1.0);

Histograms are created in Constructor of UserRunAction:

B4RunAction::B4RunAction() : G4UserRunAction(){  
  // Create analysis manager
  // ...
  // Creating histograms (1-dimensional)
  analysisManager->CreateH1("1","Edep in absorber", 
                            100, 0., 800*MeV);
  analysisManager->CreateH1("2","Edep in gap", 
                            100, 0., 100*MeV);
}

Usage in EndOfEventAction of UserEventAction:

void B4aEventAction::EndOfEventAction(const G4Run* aRun){  
  // Fill histograms
  auto analysisManager = G4AnalysisManager::Instance();
  // FillH1(G4int id, G4double value, G4double weight = 1.0);
  analysisManager->FillH1(1, fEnergyAbs);
  analysisManager->FillH1(2, fEnergyGap);
}

Create and Fill Ntuples

G4int CreateNtuple(const G4String& name,  
                   const G4String& title);
// Create columns in the last created ntuple
G4int CreateNtupleXColumn(const G4String& name);  
void  FinishNtuple();

// Create columns in the ntuple with given id
G4int CreateNtupleXColumn(G4int ntupleId,  
                          const G4String& name);
void  FinishNtuple(G4int ntupleId);

// Methods for ntuple id = FirstNtupleId
G4bool FillNtupleXColumn(G4int id, G4Xtype value);  
G4bool AddNtupleRow();

// Methods for ntuple id > FirstNtupleId (more ntuples exist)
G4bool FillNtupleIColumn(G4int ntupleId, G4int columnId,  
                         Xtype value);
G4bool AddNtupleRow(G4int ntupleId);

[X, Xtype] can be [I, G4int], [F, G4float], [D, G4double] or [S, G4String]

Ntupless are created in Constructor of UserRunAction:

B4RunAction::B4RunAction() : G4UserRunAction(){  
  // Create analysis manager and histograms
  // ...  
  // Creating ntuple
  analysisManager->CreateNtuple("B4", "Edep and TrackL");
  analysisManager->CreateNtupleDColumn("Eabs");
  analysisManager->CreateNtupleDColumn("Labs");
  analysisManager->FinishNtuple();
}

Usage in EndOfEventAction of UserEventAction:

void B4aEventAction::EndOfEventAction(const G4Run* aRun){  
  auto analysisManager = G4AnalysisManager::Instance();
  // fill ntuple
  analysisManager->FillNtupleDColumn(0, fEdep);
  analysisManager->FillNtupleDColumn(2, fLength);
  analysisManager->AddNtupleRow();  
}

ROOT Basics

histogram

void poissonHist(){  
  auto cnt_r_h=new TH1F("count_rate",
              "Count Rate;N_{Counts};# occurencies",
              100, // Number of Bins
              -0.5, // Lower X Boundary
              15.5); // Upper X Boundary
  TRandom3 rndgen;
  // simulate the measurements
  for (int imeas=0;imeas<400;imeas++)
      cnt_r_h->Fill(rndgen.Poisson(3.6));
  // draw histogram
  auto c= new TCanvas();
  cnt_r_h->Draw();
  // Print summary
  cout << "Moments of Distribution:\n"
       << " - Mean     = " << cnt_r_h->GetMean() << " +- "
                           << cnt_r_h->GetMeanError() << "\n"
       << " - Std Dev  = " << cnt_r_h->GetStdDev() << " +- "
                           << cnt_r_h->GetStdDevError() << "\n";
}

150%

tree (ntuple)

void write_ntuple_to_file(  
  const std::string& outputFileName="conductivity_experiment.root";
  unsigned int numDataPoints=1000000);
  TFile ofile(outputFileName.c_str(),"RECREATE");
  // Initialise the TNtuple
  TTree cond_data("cond_data", "Example N-Tuple");
  // define the variables and book them for the ntuple
  float pot,cur,temp,pres;
  cond_data.Branch("Potential", &pot, "Potential/F");
  cond_data.Branch("Current", &cur, "Current/F");
  cond_data.Branch("Temperature", &temp, "Temperature/F");
  cond_data.Branch("Pressure", &pres, "Pressure/F");
  for (int i=0;i<numDataPoints;++i){
    // Fill it randomly to fake the acquired data
    pot = gRandom->Uniform(0.,10.)*gRandom->Gaus(1.,0.01);
    temp = gRandom->Uniform(250.,350.)+gRandom->Gaus(0.,0.3);
    pres = gRandom->Uniform(0.5,1.5)*gRandom->Gaus(1.,0.02);
    cur = pot/(10.+0.05*(temp-300.)-0.2*(pres-1.))*gRandom->Gaus(1.,0.01);
    cond_data.Fill();
  }
  // Save the ntuple and close the file
  cond_data.Write();
  ofile.Close();
}

function and fitting

double the_gausppar(double* vars, double* pars){  
    return pars[0]*TMath::Gaus(vars[0],pars[1],pars[2])+
        pars[3]+pars[4]*vars[0]+pars[5]*vars[0]*vars[0];}

int functionAndFitting(){  
    TF1 parabola("parabola","[0]+[1]*x+[2]*x**2",0,20);
    format_line(&parabola,kBlue,2);
    TF1 gaussian("gaussian","[0]*TMath::Gaus(x,[1],[2])",0,20);
    format_line(&gaussian,kRed,2);
    TF1 gausppar("gausppar",the_gausppar,-0,20,6);
    double a=15; double b=-1.2; double c=.03;
    double norm=4; double mean=7; double sigma=1;
    gausppar.SetParameters(norm,mean,sigma,a,b,c);

    TH1F histo("histo","Signal plus background;X vals;Y Vals",50,0,20);
    histo.SetMarkerStyle(8);
    // Fake the data
    for (int i=1;i<=5000;++i) 
        histo.Fill(gausppar.GetRandom());

    // Reset the parameters before the fit and set
    gausppar.SetParameter(0,50);
    gausppar.SetParameter(1,6);
    int npar=gausppar.GetNpar();
    for (int ipar=2;ipar<npar;++ipar) gausppar.SetParameter(ipar,1);

    // perform fit ...
    auto fitResPtr = histo.Fit(&gausppar, "S");
    // ... and retrieve fit results
    fitResPtr->Print(); // print fit results
    // get covariance Matrix an print it
    TMatrixDSym covMatrix (fitResPtr->GetCovarianceMatrix());
    covMatrix.Print();

    // Set the values of the gaussian and parabola
    for (int ipar=0;ipar<3;ipar++){
        gaussian.SetParameter(ipar,gausppar.GetParameter(ipar));
        parabola.SetParameter(ipar,gausppar.GetParameter(ipar+3));
    }
      // draw results
    histo.GetYaxis()->SetRangeUser(0,250);
    histo.DrawClone("PE");
    parabola.DrawClone("Same"); gaussian.DrawClone("Same");
    return 0;
}

histogram

Execute a GEANT4 Program

Hard-coded

int main(){  
  // Construct the default run manager
  G4RunManager* runManager = new G4RunManager;
  // Set mandatory initialization classes
  runManager->SetUserInitialization(new B1DetectorConstruction);
  runManager->SetUserInitialization(new QBBC);
  runManager->SetUserAction(new B1PrimaryGeneratorAction);
  // Set user action classes
  runManager->SetUserAction(new B1SteppingAction());
  runManager->SetUserAction(new B1EventAction());
  runManager->SetUserAction(new B1RunAction());

  // Initialize G4 kernel
  runManager->Initialize();
  // start a run
  int numberOfEvent = 1000;
  runManager->BeamOn(numberOfEvent);
  // job termination
  delete runManager;
  return 0;
}

Macro file

int main(int argc,char** argv){  
  // Construct the default run manager
  G4RunManager* runManager = new G4RunManager;
  // Set mandatory initialization classes
  ...
  // Set user action classes
  ...
  // Initialize G4 kernel
  runManager->Initialize();
  //read a macro file of commands
  G4UImanager* UI = G4UImanager::GetUIpointer();
  G4String command = "/control/execute ";
  G4String fileName = argv[1];
  UI->ApplyCommand(command+fileName);

  // job termination
  delete runManager;
  return 0;
}

In command line, run it with your macro file:

./exampleB1 run_marco.mac

Interactive Mode

int main(int argc,char** argv){  
  // Construct the default run manager
  G4RunManager* runManager = new G4RunManager;
  // Set mandatory initialization classes
  ...
  // Set user action classes
  ...
  // Initialize G4 kernel
  runManager->Initialize();
  // Define UI terminal for interactive mode
  G4UIsession * session = new G4UIterminal;
  session->SessionStart();
  delete session;
  // job termination
  delete runManager;
  return 0;
}

In the prompt G4 kernel:

Idle> /run/beamOn 1000

General Case

int main(int argc,char** argv){  
  // Construct the default run manager
  G4RunManager* runManager = new G4RunManager;
  // Set mandatory initialization classes
  ...
  // Set user action classes
  ...
  // Initialize G4 kernel
  runManager->Initialize();

#ifdef G4VIS_USE
  // Initialize visualization
  G4VisManager* visManager = new G4VisExecutive;
  visManager->Initialize();
#endif

  // Get the pointer to the User Interface manager
  G4UImanager* UImanager = G4UImanager::GetUIpointer();

  if (argc!=1) {
    // batch mode
    G4String command = "/control/execute ";
    G4String fileName = argv[1];
    UImanager->ApplyCommand(command+fileName);
  } else {
    // interactive mode : define UI session
#ifdef G4UI_USE
    G4UIExecutive* ui = new G4UIExecutive(argc, argv);
#ifdef G4VIS_USE
    UImanager->ApplyCommand("/control/execute init_vis.mac");
#else
    UImanager->ApplyCommand("/control/execute init.mac");
#endif
    ui->SessionStart();
    delete ui;
#endif
  }
#ifdef G4VIS_USE
  delete visManager;
#endif
  delete runManager;
}

Homework

Modify the original B1 example to record the following information in detector fScoringVolume:

Total energy deposition in EACH EVENT;
Total tracking length of CHARGED particles in EACH EVENT;
Energy spectrum (histogram) of this simulation run: 0-6 Mev, 1024 bins.

Run exampleB1 with this macro file:

/run/initialize
/control/verbose 2
/run/verbose 2
#
/run/printProgress 1000
/run/beamOn 10000

Hints:

step->GetTrack()->GetDefinition()->GetPDGCharge() to get charge of particle;
step->GetTotalEnergyDeposit() and step->GetStepLength() to get energy and step length.

Output: data in your favorite format (recommend using G4AnalysisManager).

References

Data Acquisition:
- Textbook: Chapter 6 & 9
- Geant4 Online document (NEW!)
C++ STL and File I/O
- C++ Containers
- File input and output
ROOT Documents and Tutorial
- ROOT primer
- Tutorial examples with Interactive Notebook
- Reference Guide
Execute Geant4 program in various ways
- Textbook: Chapter 2.10

Menu

[Tutorial] GEANT4 Introduction: Data Acquisition & ROOT Basics

Content

Traditional Data Flow

User Actions (passing information)

Sensitive Detector and Hits

Components of SD

Hits and hit collection

Sensitive Detector

Associate to logical volume

Process data in G4UserEventAction

Scorer (pre-defined SD) in GEANT4

Data Container in C++

Sequence containers

Associative containers

Unordered associative containers

Container adaptors

Example: vector usage

File I/O

Example: ofstream and ifstream (binary/formatted)

Analysis Tools

Public interface: G4AnalysisManager

File handling (in UserRunAction)

Create and fill Histograms

Create and Fill Ntuples

ROOT Basics

histogram

tree (ntuple)

function and fitting

Execute a GEANT4 Program

Hard-coded

Macro file

Interactive Mode

General Case

Homework

References