gears.cc

Go to the documentation of this file.

#include <vector>
using namespace std;
#include <G4SteppingManager.hh>
#include <G4SteppingVerbose.hh>
class Output : public G4SteppingVerbose {
protected:
  void Record(); 
public:
  Output(); 
  void TrackingStarted() {
    G4SteppingVerbose::TrackingStarted();
    Record();
  } 
  void StepInfo() {
    G4SteppingVerbose::StepInfo();
    Record();
  } 
  void Reset() {
    trk.clear();
    stp.clear();
    vlm.clear();
    pro.clear();
    pdg.clear();
    pid.clear();
    xx.clear();
    yy.clear();
    zz.clear();
    dt.clear();
    de.clear();
    dl.clear();
    l.clear();
    x.clear();
    y.clear();
    z.clear();
    t.clear();
    k.clear();
    p.clear();
    px.clear();
    py.clear();
    pz.clear();
    q.clear();
    et.clear();
  }
  void SetSteppingVerbose(int level) { fManager->SetVerboseLevel(level); }
  int GetSteppingVerbose() { return fManager->GetverboseLevel(); }

  vector<int> trk;   
  vector<int> stp;   
  vector<int> vlm;   
  vector<int> pro;   
  vector<int> pdg;   
  vector<int> pid;   
  vector<double> xx; 
  vector<double> yy; 
  vector<double> zz; 
  vector<double> dt; 
  vector<double> de; 
  vector<double> dl; 
  vector<double> l;  
  vector<double> x;  
  vector<double> y;  
  vector<double> z;  
  vector<double> t;  
  vector<double> k;  
  vector<double> p;  
  vector<double> px; 
  vector<double> py; 
  vector<double> pz; 
  vector<double> q;  
  vector<double> et; 
};
//______________________________________________________________________________
//
#include <G4AnalysisManager.hh>
Output::Output() : G4SteppingVerbose() {
  auto manager = G4AnalysisManager::Instance();
  manager->CreateNtuple("t", "Geant4 step points");
  manager->CreateNtupleIColumn("n"); // total number of recorded hits
  manager->CreateNtupleIColumn("m"); // max copy number of sensitive volume
  manager->CreateNtupleIColumn("trk", trk);
  manager->CreateNtupleIColumn("stp", stp);
  manager->CreateNtupleIColumn("vlm", vlm);
  manager->CreateNtupleIColumn("pro", pro);
  manager->CreateNtupleIColumn("pdg", pdg);
  manager->CreateNtupleIColumn("pid", pid);
  manager->CreateNtupleDColumn("xx", xx);
  manager->CreateNtupleDColumn("yy", yy);
  manager->CreateNtupleDColumn("zz", zz);
  manager->CreateNtupleDColumn("dt", dt);
  manager->CreateNtupleDColumn("de", de);
  manager->CreateNtupleDColumn("dl", dl);
  manager->CreateNtupleDColumn("l", l);
  manager->CreateNtupleDColumn("x", x);
  manager->CreateNtupleDColumn("y", y);
  manager->CreateNtupleDColumn("z", z);
  manager->CreateNtupleDColumn("t", t);
  manager->CreateNtupleDColumn("k", k);
  manager->CreateNtupleDColumn("p", p);
  manager->CreateNtupleDColumn("px", px);
  manager->CreateNtupleDColumn("py", py);
  manager->CreateNtupleDColumn("pz", pz);
  manager->CreateNtupleDColumn("q", q);
  manager->CreateNtupleDColumn("et", et);
  manager->FinishNtuple();
}
//______________________________________________________________________________
//
#include <G4NavigationHistory.hh>
void Output::Record() {
  if (GetSilent() == 1) // CopyState() won't be called in G4SteppingVerbose
    CopyState();        // point fTrack, fStep, etc. to right places

  G4TouchableHandle handle = fStep->GetPreStepPoint()->GetTouchableHandle();
  int copyNo = handle->GetReplicaNumber();
  if (copyNo <= 0)
    return; // skip uninteresting volumes (copy No. of world == 0)
  if (trk.size() >= 10000) {
    G4cout << "GEARS: # of step points >=10000. Recording stopped." << G4endl;
    fTrack->SetTrackStatus(fKillTrackAndSecondaries);
    return;
  }

  trk.push_back(fTrack->GetTrackID());
  stp.push_back(fTrack->GetCurrentStepNumber());
  vlm.push_back(copyNo);
  pdg.push_back(fTrack->GetDefinition()->GetPDGEncoding());
  pid.push_back(fTrack->GetParentID());
  if (stp.back() == 0) { // step zero
    if (pid.back() != 0) // not primary particle
      pro.push_back(fTrack->GetCreatorProcess()->GetProcessType() * 1000 +
                    fTrack->GetCreatorProcess()->GetProcessSubType());
    else
      pro.push_back(1000); // primary particle
  } else {
    const G4VProcess *pr = fStep->GetPostStepPoint()->GetProcessDefinedStep();
    if (pr)
      pro.push_back(pr->GetProcessType() * 1000 + pr->GetProcessSubType());
    else
      pro.push_back(900); // not sure why pr can be zero
  }

  k.push_back(fTrack->GetKineticEnergy() / CLHEP::keV);
  p.push_back(fTrack->GetMomentum().mag() / CLHEP::keV);
  q.push_back(fStep->GetPostStepPoint()->GetCharge());
  l.push_back(fTrack->GetTrackLength() / CLHEP::mm);

  px.push_back(fTrack->GetMomentumDirection().x());
  py.push_back(fTrack->GetMomentumDirection().y());
  pz.push_back(fTrack->GetMomentumDirection().z());
  de.push_back(fStep->GetTotalEnergyDeposit() / CLHEP::keV);
  dl.push_back(fTrack->GetStepLength() / CLHEP::mm);

  t.push_back(fTrack->GetGlobalTime() / CLHEP::ns);
  x.push_back(fTrack->GetPosition().x() / CLHEP::mm);
  y.push_back(fTrack->GetPosition().y() / CLHEP::mm);
  z.push_back(fTrack->GetPosition().z() / CLHEP::mm);

  G4ThreeVector pos = handle->GetHistory()->GetTopTransform().TransformPoint(
      fStep->GetPostStepPoint()->GetPosition());
  xx.push_back(pos.x() / CLHEP::mm);
  yy.push_back(pos.y() / CLHEP::mm);
  zz.push_back(pos.z() / CLHEP::mm);
  dt.push_back(fTrack->GetLocalTime() / CLHEP::ns);

  if (de.back() > 0 &&
      G4StrUtil::contains(handle->GetVolume()->GetName(), "(S)")) {
    if (et.size() < (unsigned int)copyNo + 1)
      et.resize((unsigned int)copyNo + 1);
    et[copyNo] += de.back();
    et[0] += de.back();
  }
}
//______________________________________________________________________________
//
#include <G4OpticalSurface.hh>
struct BorderSurface {
  G4String name;           
  G4String v1;             
  G4String v2;             
  G4OpticalSurface *optic; 
  BorderSurface *next;     
};
//______________________________________________________________________________
//
#include <G4tgrLineProcessor.hh>
class LineProcessor : public G4tgrLineProcessor {
private:
  G4MaterialPropertiesTable *
  CreateMaterialPropertiesTable(const vector<G4String> &words,
                                size_t idxOfWords);

public:
  LineProcessor() : G4tgrLineProcessor(), Border(0) {};
  ~LineProcessor() {
    while (Border) { // deleting G4OpticalSurface is done in Geant4
      BorderSurface *next = Border->next;
      delete Border;
      Border = next;
    }
  }
  G4bool ProcessLine(const vector<G4String> &words);
  BorderSurface *Border; 
};
//______________________________________________________________________________
//
#include <G4NistManager.hh>
#include <G4UImessenger.hh>
#include <G4tgbMaterialMgr.hh>
G4bool LineProcessor::ProcessLine(const vector<G4String> &words) {
  // process default text geometry tags
  G4bool processed = G4tgrLineProcessor::ProcessLine(words);
  if (processed)
    return true; // no new tags involved

  // process added tags: prop & surf
  G4String tag = words[0];
  G4StrUtil::to_lower(tag);          // to lower cases
  if (tag.substr(0, 5) == ":prop") { // set optical properties of a material
    G4NistManager *mgr = G4NistManager::Instance();
    mgr->SetVerbose(2);
    G4Material *m = mgr->FindOrBuildMaterial(words[1]);
    if (m == NULL) // if not in NIST, then build in tgb
      m = G4tgbMaterialMgr::GetInstance()->FindOrBuildG4Material(words[1]);
    G4cout << "GEARS: Set optical properties of " << words[1] << ":" << G4endl;
    m->SetMaterialPropertiesTable(CreateMaterialPropertiesTable(words, 2));
    return true;
  } else if (tag.substr(0, 5) == ":surf") { // define an optical surface
    auto bdr = new BorderSurface;
    bdr->next = Border; // save current border pointer
    Border = bdr;       // overwrite current border pointer
    bdr->name = words[1];
    bdr->v1 = words[2];
    bdr->v2 = words[3];
    bdr->optic = new G4OpticalSurface(words[1]);

    G4cout << "GEARS: Set optical properties of " << bdr->name << ":"<< G4endl;
    size_t i = 4;
    while (i < words.size()) { // loop over optical surface setup lines
      G4String setting = words[i], value = words[i + 1];
      G4StrUtil::to_lower(setting); G4StrUtil::to_lower(value);
      if (setting == "property") {
        i++;
        break;
      } else if (setting == "type") {
        if (value == "dielectric_metal")
          bdr->optic->SetType(dielectric_metal);
        else if (value == "dielectric_dielectric")
          bdr->optic->SetType(dielectric_dielectric);
        else if (value == "dielectric_lut")
          bdr->optic->SetType(dielectric_LUT);
        else if (value == "dielectric_lutdavis")
          bdr->optic->SetType(dielectric_LUTDAVIS);
        else if (value == "dielectric_dichroic")
          bdr->optic->SetType(dielectric_dichroic);
        else if (value == "firsov")
          bdr->optic->SetType(firsov);
        else if (value == "x_ray")
          bdr->optic->SetType(x_ray);
        else if (value == "coated")
          bdr->optic->SetType(coated);
        else
          G4cout << "GEARS: Ignore unknown surface type " << value << G4endl;
        G4cout << "GEARS: type = " << value << G4endl;
      } else if (setting == "model") {
        if (value == "glisur")
          bdr->optic->SetModel(glisur);
        else if (value == "unified")
          bdr->optic->SetModel(unified);
        else if (value == "lut")
          bdr->optic->SetModel(LUT);
        else if (value == "davis")
          bdr->optic->SetModel(DAVIS);
        else if (value == "dichroic")
          bdr->optic->SetModel(dichroic);
        else
          G4cout << "GEARS: Ignore unknown surface model " << value << G4endl;
        G4cout << "GEARS: model = " << value << G4endl;
      } else if (setting == "finish") {
        if (value == "polished") bdr->optic->SetFinish(polished);
        else if (value == "polishedfrontpainted") bdr->optic->SetFinish(polishedfrontpainted);
        else if (value == "polishedbackpainted") bdr->optic->SetFinish(polishedbackpainted);
        else if (value == "ground") bdr->optic->SetFinish(ground);
        else if (value == "groundfrontpainted") bdr->optic->SetFinish(groundfrontpainted);
        else if (value == "groundbackpainted") bdr->optic->SetFinish(groundbackpainted);
        // --- LBNL LUT Model Finishes ---
        else if (value == "polishedlumirrorair") bdr->optic->SetFinish(polishedlumirrorair);
        else if (value == "polishedlumirrorglue") bdr->optic->SetFinish(polishedlumirrorglue);
        else if (value == "polishedair") bdr->optic->SetFinish(polishedair);
        else if (value == "polishedteflonair") bdr->optic->SetFinish(polishedteflonair);
        else if (value == "polishedtioair") bdr->optic->SetFinish(polishedtioair);
        else if (value == "polishedtyvekair") bdr->optic->SetFinish(polishedtyvekair);
        else if (value == "polishedvm2000air") bdr->optic->SetFinish(polishedvm2000air);
        else if (value == "polishedvm2000glue") bdr->optic->SetFinish(polishedvm2000glue);

        else if (value == "etchedlumirrorair") bdr->optic->SetFinish(etchedlumirrorair);
        else if (value == "etchedlumirrorglue") bdr->optic->SetFinish(etchedlumirrorglue);
        else if (value == "etchedair") bdr->optic->SetFinish(etchedair);
        else if (value == "etchedteflonair") bdr->optic->SetFinish(etchedteflonair);
        else if (value == "etchedtioair") bdr->optic->SetFinish(etchedtioair);
        else if (value == "etchedtyvekair") bdr->optic->SetFinish(etchedtyvekair);
        else if (value == "etchedvm2000air") bdr->optic->SetFinish(etchedvm2000air);
        else if (value == "etchedvm2000glue") bdr->optic->SetFinish(etchedvm2000glue);

        else if (value == "groundlumirrorair") bdr->optic->SetFinish(groundlumirrorair);
        else if (value == "groundlumirrorglue") bdr->optic->SetFinish(groundlumirrorglue);
        else if (value == "groundair") bdr->optic->SetFinish(groundair);
        else if (value == "groundteflonair") bdr->optic->SetFinish(groundteflonair);
        else if (value == "groundtioair") bdr->optic->SetFinish(groundtioair);
        else if (value == "groundtyvekair") bdr->optic->SetFinish(groundtyvekair);
        else if (value == "groundvm2000air") bdr->optic->SetFinish(groundvm2000air);
        else if (value == "groundvm2000glue") bdr->optic->SetFinish(groundvm2000glue);
        // --- DAVIS Model Finishes (Lowercase Conversion) ---
        else if (value == "rough_lut") bdr->optic->SetFinish(Rough_LUT);
        else if (value == "roughteflon_lut") bdr->optic->SetFinish(RoughTeflon_LUT);
        else if (value == "roughesr_lut") bdr->optic->SetFinish(RoughESR_LUT);
        else if (value == "roughesrgrease_lut") bdr->optic->SetFinish(RoughESRGrease_LUT);
        else if (value == "polished_lut") bdr->optic->SetFinish(Polished_LUT);
        else if (value == "polishedteflon_lut") bdr->optic->SetFinish(PolishedTeflon_LUT);
        else if (value == "polishedesr_lut") bdr->optic->SetFinish(PolishedESR_LUT);
        else if (value == "polishedesrgrease_lut") bdr->optic->SetFinish(PolishedESRGrease_LUT);
        else if (value == "detector_lut") bdr->optic->SetFinish(Detector_LUT);

        else {
          G4cout << "GEARS: Unknown surface finish: " << value << "!" << G4endl;
          abort();
        }
        G4cout << "GEARS: finish = " << value << G4endl;
      } else if (setting == "sigma_alpha") {
        bdr->optic->SetSigmaAlpha(G4UIcommand::ConvertToDouble(value));
        G4cout<<"GEARS: sigma_alpha = "<<bdr->optic->GetSigmaAlpha()<<G4endl;
      } else
        G4cout << "GEARS: Ignore unknown surface setting " << value << G4endl;
      i += 2;
    }
    if (i < words.size()) // break while loop because of "property"
      bdr->optic->SetMaterialPropertiesTable(
          CreateMaterialPropertiesTable(words, i));
    return true;
  } else
    return false;
}
//______________________________________________________________________________
//
#include <G4tgrUtils.hh>
G4MaterialPropertiesTable *
LineProcessor::CreateMaterialPropertiesTable(const vector<G4String> &words,
                                             size_t idxOfWords) {
  bool photonEnergyUnDefined = true;
  int cnt = 0;             // number of photon energy values
  double *energies = NULL; // photon energy values
  auto table = new G4MaterialPropertiesTable();
  for (size_t i = idxOfWords; i < words.size(); i++) {
    G4String property = words[i];
    G4StrUtil::to_upper(property);
    if (G4StrUtil::contains(property, "TIMECONSTANT") ||
        G4StrUtil::contains(property, "SCINTILLATIONYIELD") ||
        property == "RESOLUTIONSCALE" || property == "YIELDRATIO") {
      table->AddConstProperty(property, G4tgrUtils::GetDouble(words[i + 1]));
      G4cout << "GEARS: " << property << " = " << words[i + 1] << G4endl;
      i++; // property value has been used
    } else if (property.substr(0, 12) == "PHOTON_ENERG") {
      photonEnergyUnDefined = false;
      cnt = G4UIcommand::ConvertToInt(words[i + 1]); // get array size
      energies = new double[cnt];                    // create energy array
      for (int j = 0; j < cnt; j++)
        energies[j] = G4tgrUtils::GetDouble(words[i + 2 + j]);
      i = i + 1 + cnt; // array has been used
    } else {           // wavelength-dependent properties
      if (photonEnergyUnDefined) {
        G4cout << "GEARS: photon energies undefined, "
               << "ignore all wavelength-dependent properties!" << G4endl;
        break;
      }
      double *values = new double[cnt];
      for (int j = 0; j < cnt; j++)
        values[j] = G4tgrUtils::GetDouble(words[i + 1 + j]);
      G4cout << "GEARS: " << property << " = " << values[0] << ", " << values[1]
             << ", ..." << G4endl;
      table->AddProperty(property, energies, values, cnt);
      delete[] values;
      i = i + cnt;
    }
  }
  delete[] energies;
  return table;
}
//______________________________________________________________________________
//
#include <G4tgbDetectorBuilder.hh>
class TextDetectorBuilder : public G4tgbDetectorBuilder {
public:
  TextDetectorBuilder() : G4tgbDetectorBuilder() {
    fLineProcessor = new LineProcessor();
  }
  ~TextDetectorBuilder() { delete fLineProcessor; }
  const G4tgrVolume *ReadDetector(); 

  G4VPhysicalVolume *ConstructDetector(const G4tgrVolume *topVol);

private:
  LineProcessor *fLineProcessor; 
};
//______________________________________________________________________________
//
#include <G4tgbVolumeMgr.hh>
#include <G4tgrFileReader.hh>
#include <G4tgrVolumeMgr.hh>
const G4tgrVolume *TextDetectorBuilder::ReadDetector() {
  G4tgrFileReader *reader = G4tgrFileReader::GetInstance();
  reader->SetLineProcessor(fLineProcessor);
  reader->ReadFiles();
  G4tgrVolumeMgr *mgr = G4tgrVolumeMgr::GetInstance();
  const G4tgrVolume *world = mgr->GetTopVolume();
  return world;
}
//______________________________________________________________________________
//
#include <G4LogicalBorderSurface.hh>
G4VPhysicalVolume *
TextDetectorBuilder::ConstructDetector(const G4tgrVolume *topVol) {
  G4VPhysicalVolume *world = G4tgbDetectorBuilder::ConstructDetector(topVol);

  G4tgbVolumeMgr *tgbVolmgr = G4tgbVolumeMgr::GetInstance();
  BorderSurface *border = fLineProcessor->Border;
  while (border) {
    G4String physV1 = border->v1.substr(0, border->v1.find(":"));
    G4String physV2 = border->v2.substr(0, border->v2.find(":"));
    int copyNo1 = atoi(border->v1.substr(border->v1.find(":") + 1).data());
    int copyNo2 = atoi(border->v2.substr(border->v2.find(":") + 1).data());
    G4LogicalVolume *m1 = tgbVolmgr->FindG4PhysVol(physV1)->GetMotherLogical();
    G4LogicalVolume *m2 = tgbVolmgr->FindG4PhysVol(physV2)->GetMotherLogical();
    // search for physics volumes on the sides of the border
    G4VPhysicalVolume *v1 = 0, *v2 = 0;
    for (int i = 0; i < (int)m1->GetNoDaughters(); i++) {
      v1 = m1->GetDaughter(i);
      if (v1->GetCopyNo() == copyNo1)
        break;
    }
    for (int i = 0; i < (int)m2->GetNoDaughters(); i++) {
      v2 = m2->GetDaughter(i);
      if (v2->GetCopyNo() == copyNo2)
        break;
    }
    if (v1 && v2) {
      new G4LogicalBorderSurface(border->name, v1, v2, border->optic);
      G4cout << "Border surface " << border->name << " in between " << physV1
             << ":" << copyNo1 << " and " << physV2 << ":" << copyNo2
             << " added" << G4endl;
    }
    border = border->next;
  }
  return world;
}
//______________________________________________________________________________
//
#include <G4UIcmdWith3VectorAndUnit.hh>
#include <G4UIcmdWithAString.hh>
#include <G4UIdirectory.hh>
#include <G4UImessenger.hh>
#include <G4VUserDetectorConstruction.hh>
class Detector : public G4VUserDetectorConstruction, public G4UImessenger {
public:
  Detector();
  ~Detector() {
    delete fCmdSetB;
    delete fCmdSrc;
    delete fCmdOut;
  }
  G4VPhysicalVolume *Construct(); 
  void SetNewValue(G4UIcommand *cmd, G4String value); 

private:
  G4UIcmdWith3VectorAndUnit *fCmdSetB; 
  G4UIcmdWithAString *fCmdSrc;         
  G4UIcmdWithAString *fCmdOut;         
  G4UIcommand *fCmdLmt;                
  G4VPhysicalVolume *fWorld;
};
//______________________________________________________________________________
//
Detector::Detector()
    : G4VUserDetectorConstruction(), G4UImessenger(), fWorld(0) {
#ifdef hasGDML
  fCmdOut = new G4UIcmdWithAString("/geometry/export", this);
  fCmdOut->SetGuidance("Export geometry gdml file name");
  fCmdOut->SetParameterName("gdml geometry output", false);
  fCmdOut->AvailableForStates(G4State_Idle);
#else
  fCmdOut = 0;
#endif

  fCmdSrc = new G4UIcmdWithAString("/geometry/source", this);
  fCmdSrc->SetGuidance("Set geometry source file name");
  fCmdSrc->SetParameterName("text geometry input", false);
  fCmdSrc->AvailableForStates(G4State_PreInit);

  fCmdLmt = new G4UIcommand("/tracking/setStepLimit", this);
  fCmdLmt->SetGuidance("set max step length for a logical volume");
  fCmdLmt->AvailableForStates(G4State_Idle);
  G4UIparameter *p0 = new G4UIparameter("logical volume name", 's', false);
  G4UIparameter *p1 = new G4UIparameter("step length in mm", 'd', false);
  fCmdLmt->SetParameter(p0);
  fCmdLmt->SetParameter(p1);

  fCmdSetB = new G4UIcmdWith3VectorAndUnit("/geometry/SetB", this);
  fCmdSetB->SetGuidance("Set uniform magnetic field value.");
  fCmdSetB->SetParameterName("Bx", "By", "Bz", false);
  fCmdSetB->SetUnitCategory("Magnetic flux density");
}
//______________________________________________________________________________
//
#include <G4FieldManager.hh>
#include <G4LogicalVolumeStore.hh>
#include <G4TransportationManager.hh>
#include <G4UniformMagField.hh>
#include <G4UserLimits.hh>
#ifdef hasGDML
#include "G4GDMLParser.hh"
#endif
void Detector::SetNewValue(G4UIcommand *cmd, G4String value) {
  if (cmd == fCmdSetB) {
    auto field = new G4UniformMagField(0, 0, 0);
    field->SetFieldValue(fCmdSetB->GetNew3VectorValue(value));
    G4FieldManager *mgr =
        G4TransportationManager::GetTransportationManager()->GetFieldManager();
    mgr->SetDetectorField(field);
    mgr->CreateChordFinder(field);
    G4cout << "GEARS: Magnetic field is set to " << value << G4endl;
#ifdef hasGDML
  } else if (cmd == fCmdOut) {
    G4GDMLParser paser;
    paser.Write(value, fWorld);
#endif
  } else if (cmd == fCmdLmt) {
    istringstream iss(value);
    G4String vlm, limit;
    iss >> vlm >> limit;
    auto v = G4LogicalVolumeStore::GetInstance()->GetVolume(vlm);
    if (v) {
      v->SetUserLimits(new G4UserLimits(stof(limit)));
      G4cout << "GEARS: max step length in " << vlm << ": " << limit << " mm"
             << G4endl;
    }
  } else {                                            // cmd==fCmdSrc
    if (value.substr(value.length() - 4) != "gdml") { // text geometry input
      G4tgbVolumeMgr *mgr = G4tgbVolumeMgr::GetInstance();
      mgr->AddTextFile(value);
      auto tgb = new TextDetectorBuilder;
      mgr->SetDetectorBuilder(tgb);
      fWorld = mgr->ReadAndConstructDetector();
#ifdef hasGDML
    } else { // GDML input
      G4GDMLParser parser;
      parser.Read(value);
      fWorld = parser.GetWorldVolume();
#endif
    }
  }
}
//______________________________________________________________________________
//
#include "G4Box.hh"
#include "G4PVPlacement.hh"
G4VPhysicalVolume *Detector::Construct() {
  if (fWorld == NULL) {
    G4cout << "GEARS: no detector specified, set to a 10x10x10 m^3 box."
           << G4endl;
    auto box = new G4Box("hall", 5 * CLHEP::m, 5 * CLHEP::m, 5 * CLHEP::m);
    G4NistManager *nist = G4NistManager::Instance();
    G4Material *vacuum = nist->FindOrBuildMaterial("G4_Galactic");
    auto v = new G4LogicalVolume(box, vacuum, "hall");
    fWorld = new G4PVPlacement(0, G4ThreeVector(), v, "hall", 0, 0, 0);
  }
  return fWorld;
}
//______________________________________________________________________________
//
#include <G4GeneralParticleSource.hh>
#include <G4VUserPrimaryGeneratorAction.hh>
class Generator : public G4VUserPrimaryGeneratorAction {
private:
  G4GeneralParticleSource *fSource;

public:
  Generator() : G4VUserPrimaryGeneratorAction(), fSource(0) {
    fSource = new G4GeneralParticleSource;
  }
  virtual ~Generator() { delete fSource; }
  virtual void GeneratePrimaries(G4Event *evt) {
    fSource->GeneratePrimaryVertex(evt);
  } 
};
//______________________________________________________________________________
//
#include <G4Run.hh>
#include <G4RunManagerFactory.hh>
#include <G4UserRunAction.hh>
class RunAction : public G4UserRunAction {
public:
  void BeginOfRunAction(const G4Run *) {
    auto a = G4AnalysisManager::Instance();
    if (a->GetFileName() == "")
      return;
    a->OpenFile();
    Output *o = ((Output *)G4VSteppingVerbose::GetInstance());
    if (o->GetSteppingVerbose() == 0) { // in case of /tracking/verbose 0
      o->SetSilent(1);                  // avoid screen dump
      o->SetSteppingVerbose(1);         // enable calling StepInfo() in
                                // G4SteppingManager
    }
  } 
  void EndOfRunAction(const G4Run *) {
    auto a = G4AnalysisManager::Instance();
    if (a->GetFileName() != "") {
      a->Write();
      a->CloseFile();
    }
  } 
};
//______________________________________________________________________________
//
void SaveAndResetEvent() {
  auto a = G4AnalysisManager::Instance();
  Output *o = ((Output *)G4VSteppingVerbose::GetInstance());
  if (a->GetFileName() != "" && o->stp.size() != 0) {
    a->FillNtupleIColumn(0, o->stp.size());
    a->FillNtupleIColumn(1, o->et.size() - 1);
    a->AddNtupleRow();
  } // save n-tuple if it is not empty and output file name is specified
  o->Reset(); // reset Output member variables for new record
} 
//______________________________________________________________________________
//
#include <G4UserEventAction.hh>
class EventAction : public G4UserEventAction {
public:
  void EndOfEventAction(const G4Event *) { SaveAndResetEvent(); }
};
//______________________________________________________________________________
//
#include <G4UIcmdWithADoubleAndUnit.hh>
#include <G4UserStackingAction.hh>
class StackingAction : public G4UserStackingAction, public G4UImessenger {
private:
  double fT0;         
  double fTimeWindow; 
  G4UIcmdWithADoubleAndUnit *fCmdT; 
public:
  StackingAction()
      : G4UserStackingAction(), G4UImessenger(), fT0(0), fTimeWindow(0),
        fCmdT(0) {
    fCmdT = new G4UIcmdWithADoubleAndUnit("/grdm/setTimeWindow", this);
    fCmdT->SetGuidance("Time window to split a radioactive decay chain.");
    fCmdT->SetGuidance("If a daughter nucleus appears after the window,");
    fCmdT->SetGuidance("it is saved in a new entry in the output ntuple.");
    fCmdT->SetGuidance("---Set it to <=0 to disable the splitting---");
    fCmdT->SetParameterName("time window", false, true);
    fCmdT->SetDefaultUnit("s");
    fCmdT->AvailableForStates(G4State_PreInit, G4State_Idle);
  } 
  ~StackingAction() { delete fCmdT; }
  G4ClassificationOfNewTrack ClassifyNewTrack(const G4Track *trk) {
    if (fTimeWindow <= 0)
      return fUrgent; // no need to split
    if (trk->GetGlobalTime() > fT0 + fTimeWindow)
      return fWaiting; // split
    else
      return fUrgent; // too fast to be split
  } 
  void NewStage() { // called after processing urgent trk, before waiting trk
    if (fTimeWindow <= 0)
      return; // do nothing if no time window is specified
    Output *o = ((Output *)G4VSteppingVerbose::GetInstance());
    fT0 = o->t.back();   // update the reference time to the latest decay time
    SaveAndResetEvent(); // end an event manually
  } 
  void SetNewValue(G4UIcommand *cmd, G4String value) {
    if (cmd != fCmdT)
      return;
    fTimeWindow = fCmdT->GetNewDoubleValue(value);
  }
};
//______________________________________________________________________________
//
#include <G4VUserActionInitialization.hh>
class Action : public G4VUserActionInitialization {
  void Build() const {
    SetUserAction(new RunAction);
    SetUserAction(new Generator);
    SetUserAction(new EventAction);
    SetUserAction(new StackingAction);
  }
};
//______________________________________________________________________________
//
#include <G4PhysListFactory.hh>
#include <G4ScoringManager.hh>
#include <G4StepLimiterPhysics.hh>
#include <G4UIExecutive.hh>
#include <G4UImanager.hh> // needed for g4.10 and above
#include <G4VisExecutive.hh>
int main(int argc, char **argv) {
  // inherit G4SteppingVerbose instead of G4UserSteppingAction to record data
  G4VSteppingVerbose::SetInstance(new Output); // must be before run manager
  auto run =
      G4RunManagerFactory::CreateRunManager(G4RunManagerType::SerialOnly);
  G4PhysListFactory factory;
  auto physics = factory.ReferencePhysList();
  physics->RegisterPhysics(new G4StepLimiterPhysics());
  run->SetUserInitialization(physics);      // initialize physics
  run->SetUserInitialization(new Detector); // initialize detector
  run->SetUserInitialization(new Action);   // initialize user actions
  G4ScoringManager::GetScoringManager();    // enable built-in scoring cmds
  G4UIExecutive *ui = nullptr;              // assume batch mode
  if (argc == 1) {
    ui = new G4UIExecutive(argc, argv);
  }                                       // interactive mode
  auto vis = new G4VisExecutive("quiet"); // visualization
  vis->Initialize();                      // do this after ui mode is decided
  if (ui) {                               // interactive mode
    ui->SessionStart();                   // do this after vis
    delete ui;
  } else { // batch mode
    G4String cmd = "/control/execute ";
    G4UImanager::GetUIpointer()->ApplyCommand(cmd + argv[1]);
  }
  delete vis;
  delete run;
  return 0;
}
// -*- C++; indent-tabs-mode:nil; tab-width:2 -*-
// vim: ft=cpp:ts=2:sts=2:sw=2:et