Skip to content
Snippets Groups Projects
Commit 8158220c authored by Andy Haas's avatar Andy Haas
Browse files

Tries to remove G4 volume overlaps, so far just in TRT 

Former-commit-id: 7ccbdbf3
parent 9697aabd
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Simulation/G4Utilities/Geo2G4/src/ExtParameterisedVolumeBuilder.cxx
+ 286
− 6
View file @ 8158220c
......@@ -13,6 +13,12 @@
#include "G4ReflectionFactory.hh"
#include "G4VPVParameterisation.hh"
#include "G4PVParameterised.hh"
#include "G4AffineTransform.hh"
#include "G4UnitsTable.hh"
#include "G4LogicalVolume.hh"
#include "G4VSolid.hh"
#include "G4SubtractionSolid.hh"
#include "G4PhysicalVolumeStore.hh"
#include "globals.hh"
#include "SimHelpers/ServiceAccessor.h"
......@@ -35,8 +41,120 @@ ExtParameterisedVolumeBuilder::ExtParameterisedVolumeBuilder(std::string n):
{
}
//stolen from G4bool G4PVPlacement::CheckOverlaps, but revised so we can tell what it overlaps with...
//returns 0 if no overlap, overwise returns overlapping volume pointer...and whether it's with the mother...
G4VPhysicalVolume* check_overlap(G4VPhysicalVolume *me, int &mom){
int res = 1000;//number of points to check per volume!
mom=-1;//will return whether it is the mom that is overlapping or not...
G4VSolid* solid = me->GetLogicalVolume()->GetSolid();
G4LogicalVolume* motherLog = me->GetMotherLogical();
//std::cout<<"check_overlap"<<std::endl;
if (!motherLog) { return 0; }
G4VSolid* motherSolid = motherLog->GetSolid();
// Create the transformation from daughter to mother
G4AffineTransform Tm( me->GetRotation(), me->GetTranslation() );
std::cout<<"check_overlap for "<<me->GetName()<<" with "<<res<<" points"<<std::endl;
for (G4int n=0; n<res; n++)
{
// Generate a random point on the solid's surface
G4ThreeVector point = solid->GetPointOnSurface();
// Transform the generated point to the mother's coordinate system
G4ThreeVector mp = Tm.TransformPoint(point);
// Checking overlaps with the mother volume
if (motherSolid->Inside(mp)==kOutside)
{
G4double distin = motherSolid->DistanceToIn(mp);
if (distin>0)
{
// std::ostringstream message;
// message << "Overlap with mother volume !" << G4endl
// << " Overlap is detected for volume "
// << me->GetName() << G4endl
// << " with its mother volume "
// << motherLog->GetName() << G4endl
// << " at mother local point " << mp << ", "
// << "overlapping by at least: "
// << G4BestUnit(distin, "Length");
//G4Exception("G4PVPlacement::CheckOverlaps()","GeomVol1002", JustWarning, message);
mom=1;
return 0;//motherLog;
}
}
// Checking overlaps with each 'sister' volume
for (G4int i=0; i<motherLog->GetNoDaughters(); i++)
{
G4VPhysicalVolume* daughter = motherLog->GetDaughter(i);
if (daughter == me) { continue; }
// Create the transformation for daughter volume and transform point
G4AffineTransform Td( daughter->GetRotation(),
daughter->GetTranslation() );
G4ThreeVector md = Td.Inverse().TransformPoint(mp);
G4VSolid* daughterSolid = daughter->GetLogicalVolume()->GetSolid();
if (daughterSolid->Inside(md)==kInside)
{
G4double distout = daughterSolid->DistanceToOut(md);
if (distout>0)
{
// std::ostringstream message;
// message << "Overlap with volume already placed !" << G4endl
// << " Overlap is detected for volume "
// << me->GetName() << G4endl
// << " with " << daughter->GetName() << " volume's"
// << G4endl
// << " local point " << md << ", "
// << "overlapping by at least: "
// << G4BestUnit(distout,"Length");
// G4Exception("G4PVPlacement::CheckOverlaps()",
// "GeomVol1002", JustWarning, message);
mom=0;
return daughter;
}
}
// Now checking that 'sister' volume is not totally included and
// overlapping. Do it only once, for the first point generated
if (n==0)
{
// Generate a single point on the surface of the 'sister' volume
// and verify that the point is NOT inside the current volume
G4ThreeVector dPoint = daughterSolid->GetPointOnSurface();
// Transform the generated point to the mother's coordinate system
// and finally to current volume's coordinate system
G4ThreeVector mp2 = Td.TransformPoint(dPoint);
G4ThreeVector msi = Tm.Inverse().TransformPoint(mp2);
if (solid->Inside(msi)==kInside)
{
// std::ostringstream message;
// message << "Overlap with volume already placed !" << G4endl
// << " Overlap is detected for volume "
// << me->GetName() << G4endl
// << " apparently fully encapsulating volume "
// << daughter->GetName() << G4endl
// << " at the same level !";
// G4Exception("G4PVPlacement::CheckOverlaps()",
// "GeomVol1002", JustWarning, message);
mom=2;
return daughter;
}
}
}
}
return 0;//no overlap!
}
G4LogicalVolume* ExtParameterisedVolumeBuilder::Build(const PVConstLink theGeoPhysVolume, OpticalVolumesMap* optical_volumes) const
{
MsgStream log(msgSvc(),GetKey());
PVConstLink theGeoPhysChild;
const GeoSerialTransformer* serialTransformerChild=0;
G4LogicalVolume* theG4LogChild;
......@@ -155,12 +273,173 @@ G4LogicalVolume* ExtParameterisedVolumeBuilder::Build(const PVConstLink theGeoPh
if (nameChild == "ANON") nameChild=theG4LogChild->GetName();
G4PhysicalVolumesPair pvPair = G4ReflectionFactory::Instance()->Place(theG4Position,
nameChild,
theG4LogChild,
theG4LogVolume,
false,
id);
log << MSG::DEBUG << "ACH655 Placing child: "<<nameChild<<"["<<id<<"]"<<" in mother "<<theG4LogVolume->GetName()<<endmsg;
bool surfcheck=true; //whether to check for overlaps
bool doplace=true;
if ("TRT::MainRadiatorA"==nameChild && 16969==id) surfcheck=true;
if ("TRT::ThinRadiatorA"==nameChild && 16969==id) surfcheck=true;
if ("TRT::InnerSupportA"==nameChild && 16969==id) surfcheck=true;
if ("TRT::OuterSupportA"==nameChild && 16969==id) surfcheck=true;
if ("TRT::MiddleRadiatorA"==nameChild && 16969==id) surfcheck=true;
//if ("TRT::Straw"==nameChild) doplace=false;//don't build straws for now
//if ("TRT::StrawPlane"==nameChild && id!=2) doplace=false;//only build one strawplane for now
if (std::string(nameChild).find("LAr:")!=std::string::npos) surfcheck=false;
if (std::string(nameChild).find("Tile:")!=std::string::npos) surfcheck=false;
if (std::string(nameChild).find("DriftTube")!=std::string::npos) surfcheck=false;
if (std::string(nameChild).find("Muon:")!=std::string::npos) surfcheck=false;
if (std::string(theG4LogVolume->GetName()).find("Muon:")!=std::string::npos) surfcheck=false;//or if the mom is muon
if (std::string(nameChild).find("Pixel:")!=std::string::npos) surfcheck=false;
if (std::string(nameChild).find("SCT:")!=std::string::npos) surfcheck=false;
if (std::string(theG4LogVolume->GetName()).find("SCT:")!=std::string::npos) surfcheck=false;//or if the mom is sct
if (std::string(nameChild).find("_cut")!=std::string::npos) {
doplace=false;//don't add again for now, if we already cut it
log<<MSG::WARNING<<"Why is something that was cut already now being placed again?"<<endmsg;
}
bool dosurfcheckfuringplacement=false;//to run the printout in G4 for surfcheck placements
if (doplace) {
G4PhysicalVolumesPair pvPair = G4ReflectionFactory::Instance()->Place(theG4Position,
nameChild,
theG4LogChild,
theG4LogVolume,
false,
id,
surfcheck&&dosurfcheckfuringplacement);
static int trytofixoverlaps=-1;
if (trytofixoverlaps<0) {//just check it once...
char* pTRY; pTRY=getenv("HAAS_TryToFixOverlaps");
if (pTRY!=NULL) trytofixoverlaps=atoi(pTRY); else trytofixoverlaps=0;
log<<MSG::INFO<<"Try to fix overlaps: "<<trytofixoverlaps<<endmsg;
}
static int docheckoverlaps=-1;
if (docheckoverlaps<0) {//just check it once...
char* pTRY; pTRY=getenv("HAAS_CheckOverlaps");
if (pTRY!=NULL) docheckoverlaps=atoi(pTRY); else docheckoverlaps=0;
log<<MSG::INFO<<"Check for overlaps: "<<docheckoverlaps<<endmsg;
}
while (surfcheck){//keep looping until we break out...
bool triedtofix=false;
int mom=-2; G4VPhysicalVolume* other = 0; if (docheckoverlaps||trytofixoverlaps) {other=check_overlap(pvPair.first,mom);}
if (other){
log<<MSG::INFO<<"ACH649 : "<<other->GetName()<<" is overlapping "<<pvPair.first->GetName()<<" and mom is "<<mom<<endmsg;
if (!trytofixoverlaps){
static int ci = 0; std::string cis; std::stringstream cisout;
cisout<<"_cut"; cisout << ++ci; cis = cisout.str();//keep track of increasing cut N
other->SetName(other->GetName()+cis);
pvPair.first->SetName(pvPair.first->GetName()+cis);
std::cout<<"Now have "<<pvPair.first->GetName()<<" and "<<other->GetName()<<std::endl;
break;//
}
}
else {
break;//no overlaps
}
if (other&&trytofixoverlaps){
//do something to try to fix the overlap!
if (0==mom){
static int ci = 0; std::string cis; std::stringstream cisout;
cisout<<"_cut"; cisout << ++ci; cis = cisout.str();//keep track of increasing cut N
//easiest case, overlaping sisters within the same mother
//first make the new solid, via subtraction
G4RotationMatrix *myrot=new G4RotationMatrix();//initially Identity
if (other->GetRotation()) (*myrot)= *(other->GetRotation());//make sure it's not null
(*myrot)*=pvPair.first->GetObjectRotationValue();//straw rot * inverse support rot?
G4SubtractionSolid* cutout = new G4SubtractionSolid(pvPair.first->GetName()+cis, //TRT::InnerSupportB_cut5
pvPair.first->GetLogicalVolume()->GetSolid(), //TRT::InnerSupportB solid
other->GetLogicalVolume()->GetSolid(), //TRT::StrawPlane solid
myrot,//made just above
other->GetTranslation() - pvPair.first->GetTranslation() //straw - support
);
log<<MSG::INFO<<"ACH544 Replacing "<<pvPair.first->GetName()<<" volume of size "
<<pvPair.first->GetLogicalVolume()->GetSolid()->GetCubicVolume()
<<" with cutout volume "<<cutout->GetCubicVolume()<<endmsg;
//then we copy the logical volume, since it will have a different (cutout) solid
G4LogicalVolume *lv=pvPair.first->GetLogicalVolume();
std::cout<<"old name is "<<lv->GetName()<<" and cis is "<<cis<<std::endl;
G4LogicalVolume *lvnew = new G4LogicalVolume(cutout, lv->GetMaterial(),
lv->GetName()+cis, lv->GetFieldManager(), lv->GetSensitiveDetector(), lv->GetUserLimits() );
//copy the daughters of the old volume into the new volume
for (int d=0; d<lv->GetNoDaughters(); ++d){
lvnew->AddDaughter(lv->GetDaughter(d));
}
//really we should copy the daughters too, and then do daughter->SetMotherLogical(pvPair.first->GetLogicalVolume()); !!!
//TODO...
log<<MSG::INFO<<"New LogicalVolume "<<lvnew->GetName()<<" has "<<lvnew->GetNoDaughters()<<" daughters"<<endmsg;
//remove the old and make a new PhysicalVolume
G4String newvolname = pvPair.first->GetName()+cis;
pvPair.first->SetName(pvPair.first->GetName()+"_old");
G4PhysicalVolumeStore::DeRegister(pvPair.first);
delete pvPair.first;//!!!
G4PhysicalVolumesPair pvPair2 = G4ReflectionFactory::Instance()->Place(theG4Position,
newvolname,
lvnew,
theG4LogVolume,
false,
id,
true);//surfcheck&&dosurfcheckfuringplacement);
other->SetName(other->GetName()+cis);
if (pvPair2.first) std::cout<<"Now have first volume "<<pvPair2.first->GetName()<<" and "<<other->GetName()<<std::endl;
if (pvPair2.second) std::cout<<"Now have second volume "<<pvPair2.first->GetName()<<" and "<<other->GetName()<<std::endl;
triedtofix=true;
pvPair.first=pvPair2.first;//for the next time through the while loop checking for more overlaps...
}
else if (2==mom){
//a sister is totally inside me?
}
else {
log<<MSG::ERROR<<" How can mom be "<<mom<<" here ???"<<endmsg;
}
//
}//other
else if (mom==1 && trytofixoverlaps){
log<<MSG::WARNING<<"ACH649 : overlapping pvPair first with mother (not fixing yet), "
<<pvPair.first->GetName()<<" and mom is "<<mom<<endmsg;
}
if (pvPair.second){
mom=-2;
other = check_overlap(pvPair.second,mom);
if (other){
log<<MSG::WARNING<<"ACH649 : "<<other->GetName()<<" is overlapping pvPair second (not fixing yet), "
<<pvPair.second->GetName()<<" and mom is "<<mom<<endmsg;
}
else if (mom==1){
log<<MSG::WARNING<<"ACH649 : overlapping pvPair second (not fixing yet), "
<<pvPair.second->GetName()<<" and mom is "<<mom<<endmsg;
}
}
if (triedtofix){
int mom2=-2; G4VPhysicalVolume* other2 = check_overlap(pvPair.first,mom2);//recheck!
if (other2){
log<<MSG::INFO<<"ACH644 : "<<other2->GetName()<<" is still overlapping "<<pvPair.first->GetName()<<" and mom2 is "<<mom2<<endmsg;
//try again...
}
else if (mom2==1){
log<<MSG::WARNING<<"ACH644 : overlapping pvPair first with mother STILL "<<pvPair.first->GetName()<<" and mom2 is "<<mom2<<endmsg;
break;//not fixing it yet...
}
else{
log<<MSG::INFO<<"ACH644 : no overlaps for "<<pvPair.first->GetName()<<" now and mom2 is "<<mom2<<endmsg;
break;//good jorb!
}
}
}//while surfcheck
// if GeoModel volume is optical store it in the map
if(optical_volumes!=0)
......@@ -170,6 +449,7 @@ G4LogicalVolume* ExtParameterisedVolumeBuilder::Build(const PVConstLink theGeoPh
if(opticalGeoPhysChild)
(*optical_volumes)[opticalGeoPhysChild] = pvPair.first;
}
}//doplace
}
av.next();
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment