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

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

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.

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";
}

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;
}

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