Hello, I am experiencing a weird issue where some of my events completely stop. I am running Geant4 Version v11.0.2.
I am simulating scintillation photons in Liquid Argon (LAr) being WLS by TPB. I have two concentric cylinders of LAr, and the inner one is coated by TPB to WLS the scintillation photons created in it (TPB is a daughter of innerCell).
around the innerCell I have a cylindrical layer of Aluminium which is covered with a skin surface to reflect the incident photons back to the innerCell LAr, and also top/bottom bases made of Aluminium covered with the same skin surface for reflection. I then measure the hit positions of the photons on the bases and sides of the innerCell by looking at the boundary steps of photons from TPB→innerCell and save them to a .csv file.
When running a large amount of events the simulation can sometimes come to a complete stop, even though most events process as intended,and Im not sure why this happens.
I did notice alot of ping-pong steps between the TPB and innerCell / TPB layers themselves which are of 0 length, but not sure if this is related, and I have not been able to solve it by say moving the TPB layers so they are well inside the innerCell instead of sitting exactly on the edges. I know the event stops since Im printing every step taken by every track and the prints just stop (not related to lag).
This is the relevant part of my DetectorConstruction.cc:
// inner cell
innerDiameter = 10.*cm;
innerHeight = 10.*cm;innerCellSolid = new G4Tubs(“innerCell”,0., innerDiameter/2, innerHeight/2, 0., 360. * deg);
innerCellLogic = new G4LogicalVolume(innerCellSolid, //its solid
fAbsorberMaterial, //its material
“innerCell”); //its nameinnerCellPhysi = new G4PVPlacement(0, //no rotation
G4ThreeVector(0,0.,0.),
innerCellLogic, //its logical volume
“innerCell”, //its name
outerCellLogic, //its mother
false, //no boulean operat
0); //copy numberstd::cout << “after inner cell” << std::endl;
innerLayerSolid = new G4Tubs(“innerLayer”,innerDiameter/2, innerDiameter/2 + 1.5*mm, innerHeight/2, 0., 360. * deg);
innerLayerLogic = new G4LogicalVolume(innerLayerSolid, //its solid
Al, //its material
“innerLayer”); //its nameinnerLayerPhysi = new G4PVPlacement(0,
G4ThreeVector(),
innerLayerLogic, //its logical volume
“innerLayer”, //its name
outerCellLogic, //its mother
false, //no boulean operat
0);std::cout << “after inner layer” << std::endl;
G4OpticalSurface* OpSurface = new G4OpticalSurface(“name”);
//G4LogicalBorderSurface* Surface = new G4LogicalBorderSurface(“name”,innerCellPhysi,innerLayerPhysi,OpSurface);OpSurface->SetType(dielectric_metal);
OpSurface->SetModel(unified);
OpSurface->SetFinish(polished);
std::vector pp = {1.5eV, 4.0eV,7.0eV,10.0eV};
//rindex = {1.35, 1.40};
std::vector reflectivity = {0.98, 0.98,0.2,0.1};G4MaterialPropertiesTable* SMPT = new G4MaterialPropertiesTable();
//SMPT->AddProperty(“RINDEX”, pp, rindex);
SMPT->AddProperty(“REFLECTIVITY”, pp, reflectivity);OpSurface->SetMaterialPropertiesTable(SMPT);
new G4LogicalSkinSurface(“skinInnerLayer”,innerLayerLogic,OpSurface);// Top/Bottom optical planes (same radius as inner cell), placed flush
auto planeMat = G4NistManager::Instance()->FindOrBuildMaterial(“G4_Al”); // any bulk; detection is on surface
auto r = innerDiameter/2;
auto t = 1.5*mm;auto topSolid = new G4Tubs(“TopPlane”, 0, r, t/2, 0, CLHEP::twopi);
auto botSolid = new G4Tubs(“BotPlane”, 0, r, t/2, 0, CLHEP::twopi);
auto topLV = new G4LogicalVolume(topSolid, planeMat, “TopPlaneLV”);
auto botLV = new G4LogicalVolume(botSolid, planeMat, “BotPlaneLV”);auto topPV = new G4PVPlacement(nullptr, {0, 0, +innerHeight/2 + t/2}, topLV, “TopPlanePV”, outerCellLogic, false, 0);
auto botPV = new G4PVPlacement(nullptr, {0, 0, -innerHeight/2 - t/2}, botLV, “BotPlanePV”, outerCellLogic, false, 0);new G4LogicalSkinSurface(“skinTop”,topLV,OpSurface);
new G4LogicalSkinSurface(“skinBot”,botLV,OpSurface);auto mkReflector = [&](const char* name){
auto s = new G4OpticalSurface(name);
s->SetType(dielectric_metal);
s->SetModel(unified);
s->SetFinish(polished);
auto mpt = new G4MaterialPropertiesTable();
// Example: flat 95% reflectivity in visible; add your spectrum
G4double e[2] = {1.5eV, 4.0eV};
G4double R[2] = {0.98, 0.98};
mpt->AddProperty(“REFLECTIVITY”, e, R, 2);
s->SetMaterialPropertiesTable(mpt);
return s;
};
auto reflTop = mkReflector(“ReflTop”);
auto reflBot = mkReflector(“ReflBot”);// TPB Volumes
const G4double tTPB = 1.5um; // ~1–3 µm film is typical
const G4double inset = 0.1mm; // instead of 5 um// mother: LAr logical volume is LArLV (you already have it)
// Side coating: a thin cylindrical shell hugging the wall
auto tpbSideSolid = new G4Tubs(“TPB_Side”, r - tTPB - inset, r - inset, 0.5*innerHeight - tTPB - inset, 0.*deg, 360.*deg);
//auto tpbSideLV = new G4LogicalVolume(tpbSideSolid, tpbMat, “TPB_Side_LV”);
//auto tpbSidePV = new G4PVPlacement(nullptr, G4ThreeVector(), tpbSideLV, “TPB_Side_PV”, innerCellLogic, true, 0);// Bottom disk
auto tpbBotSolid = new G4Tubs(“TPB_Bot”, 0., r, 0.5*tTPB, 0.deg, 360.deg);
//auto tpbBotLV = new G4LogicalVolume(tpbBotSolid, tpbMat, “TPB_Bot_LV”);
//auto tpbBotPV = new G4PVPlacement(nullptr, G4ThreeVector(0,0,-0.5innerHeight + 0.5tTPB + inset ), tpbBotLV,“TPB_Bot_PV”, innerCellLogic, true, 0);// Top disk
auto tpbTopSolid = new G4Tubs(“TPB_Top”, 0., r, 0.5*tTPB, 0.deg, 360.deg);
//auto tpbTopLV = new G4LogicalVolume(tpbTopSolid, tpbMat, “TPB_Top_LV”);
//auto tpbTopPV = new G4PVPlacement(nullptr, G4ThreeVector(0,0, +0.5innerHeight - 0.5tTPB - inset), tpbTopLV,“TPB_Top_PV”, innerCellLogic, true, 0);auto tpbUnion1 = new G4UnionSolid(“TPB_union”, tpbSideSolid, tpbTopSolid, nullptr, G4ThreeVector(0,0, +0.5innerHeight - 0.5tTPB - inset));
auto tpbUnion = new G4UnionSolid(“TPB_union2”, tpbUnion1, tpbBotSolid, nullptr, G4ThreeVector(0,0,-0.5innerHeight + 0.5tTPB + inset ));auto tpbLV = new G4LogicalVolume(tpbUnion, tpbMat, “TPB_LV”);
auto tpbPV = new G4PVPlacement(nullptr, {}, tpbLV, “TPB_PV”, innerCellLogic, true, 0);
Note that I previously did not use unionSolid, I just tried it as a method of solving the ping-pong but this didn’t work.
This is the relevant part of my SteppingAction.cc:
if(particleDef == opticalphoton){
if(creator && creator->GetProcessName() == "OpWLS" ){
double t = post->GetGlobalTime()/CLHEP::ns; // Time of flight at boundary
if (post->GetStepStatus() == fGeomBoundary) {
const bool fromTPB = (preName.rfind("TPB", 0) == 0); // any TPB_* PV
const bool toInner = (postName == "innerCell");
// accept only crossings from tpb to the inner cell
if (fromTPB && toInner){
//G4cout << "from tpb to inner" << G4endl;
//fEventAction->numPhotons += 1;
// classify surface based on which TPB piece we just left (prePV)
char surfC = 'O';
const G4String& preName = prePV->GetName();
if (preName.find("Side") != std::string::npos) surfC = 'S';
else if (preName.find("Top") != std::string::npos) surfC = 'T';
else if (preName.find("Bot") != std::string::npos) surfC = 'B';
// else 'O' (other) by default
const auto p = post->GetPosition();
if (surfC != 'O') {
if (t < fEventAction->tZeroEx){
fEventAction->tZeroEx = t;
}
}
// Helpers
auto wrap = [](double a, double w) {
double r = std::fmod(a, w);
if (r < 0) r += w;
return r;
};
int u_int=0, v_int=0; // quantized positions
// Bases: (u,v)=(x_cm,z_cm) → mm ticks
if (surfC=='T' || surfC=='B') {
u_int = static_cast<int>(std::llround( (p.x()/CLHEP::cm) * kPosScale_cm_to_mm ));
v_int = static_cast<int>(std::llround( (p.z()/CLHEP::cm) * kPosScale_cm_to_mm ));
}
// Side: (u,v)=(phi_rad,y_cm) → tile index & mm ticks
else if (surfC=='S') {
//tile count for R=5 cm, tile=0.6 cm
constexpr double R_cm = 5.0;
constexpr double tile_cm = 0.6;
const int nphi = static_cast<int>(std::floor((kTwoPi*R_cm)/tile_cm)); // 52
const double dphi = kTwoPi / nphi;
double phi = std::atan2(p.z(), p.x()); // [-pi, pi]
phi = wrap(phi, kTwoPi); // [0, 2pi)
int k = static_cast<int>(std::floor(phi / dphi));
if (k >= nphi) k = nphi - 1; // guard
u_int = k; // store tile index instead of φ (massive savings)
v_int = static_cast<int>(std::llround( (p.y()/CLHEP::cm) * kPosScale_cm_to_mm ));
//u_int = static_cast<int>(std::llround( (p.x()/CLHEP::cm) * kPosScale_cm_to_mm ));
//v_int = static_cast<int>(std::llround( (p.z()/CLHEP::cm) * kPosScale_cm_to_mm ));
}
// Time in 10 ps ticks
int t_10ps = static_cast<int>(std::llround(t * kTimeScale_10ps ));
if (fEventAction->ShouldLogWLS(trackID, surfC, u_int, v_int, t_10ps)) {
//G4cout << "should log WLS" << G4endl;
//if (track->GetCurrentStepNumber() > 2500) G4cout << "when logging wls, step number is : " << track->GetCurrentStepNumber() << G4endl;
fEventAction->BufferHitRow({std::to_string(eventID),std::to_string(trackID),std::string(1,surfC),std::to_string(u_int),std::to_string(v_int),std::to_string(t_10ps)});
}
}
}
}
}
}
Any help would be much appreciated, hinting as to what could be causing this, because I’ve been spinning my wheels trying to solve this.
Thanks!