#include <EDepSimRootGeometryManager.hh>

Classes
struct	DrawAttributes

Public Member Functions
virtual	~RootGeometryManager ()

void	Export (const char *file)
	Export the geometry to a file. More...

void	Update (const G4VPhysicalVolume *aWorld, bool validate)
	Update the root geometry to match the g4 geometry. More...

void	Update (std::string gdmlFile, bool validate)
	Set the root geometry from a gdml file. More...

void	Validate ()
	Make sure that the current geometry passes a bunch of tests. More...

int	GetNodeId (const G4ThreeVector &pos)
	Get a volume ID base on the volume position. More...

void	SetDrawAtt (G4Material *material, int color, double opacity=1.0)
	Set material color. More...

G4Color	GetG4Color (G4Material *material)

virtual void	ShouldPrintMass (std::string name)
	A method to request that a volume mass should be printed;. More...

Static Public Member Functions
static EDepSim::RootGeometryManager *	Get (void)
	If a persistency manager has not been created, create one. More...

Protected Member Functions
	RootGeometryManager ()
	use Get() instead More...

int	GetColor (const G4VPhysicalVolume *vol)

double	GetOpacity (const G4VPhysicalVolume *vol)

Private Types
typedef std::map< G4String, DrawAttributes >	AttributeMap

Private Member Functions
TGeoShape *	CreateShape (const G4VSolid theSolid, TGeoMatrix *mat=NULL)

TGeoVolume *	CreateVolume (const G4VSolid theSolid, std::string theName, TGeoMedium theMedium)
	Create a new ROOT volume object. More...

bool	CreateEnvelope (const G4VPhysicalVolume theVol, TGeoManager theEnvelope, TGeoVolume *theMother)

int	HowManySimilarNodesInVolume (TGeoVolume *theMother, std::string theName)

void	CreateMaterials (const G4VPhysicalVolume *theVol)

virtual bool	IgnoreVolume (const G4VPhysicalVolume *theVol)

virtual bool	PrintMass (const G4VPhysicalVolume *theVol)
	A method to flag that a volume mass should be printed. More...

virtual std::string	MaterialName (const G4VPhysicalVolume *theVol)

virtual TGeoMedium *	AverageMaterial (const G4VPhysicalVolume *theVol)

Private Attributes
std::map< G4String, TGeoMedium * >	fMedium
	A map between G4 material names and Root Material definitions. More...

std::map< G4String, TGeoElement * >	fIsotope
	A map between G4 isotope names and Root Element definitions. More...

std::map< G4String, int >	fNodeCount
	A map between G4 volume names and count of volumes with that name. More...

AttributeMap	fColorMap
	A map between a material name and a color. More...

std::vector< G4String >	fPrintMass
	A vector of volume names to print. More...

std::set< G4String >	fPrintedMass
	A map of which masses have been printed. More...

std::vector< G4String >	fNameStack

std::map< G4LogicalVolume , TGeoVolume >	fKnownVolumes

bool	fCreateAllVolumes

Static Private Attributes
static EDepSim::RootGeometryManager *	fThis = NULL
	The pointer to the instantiation of this object. More...

Detailed Description

Definition at line 32 of file EDepSimRootGeometryManager.hh.

Member Typedef Documentation

typedef std::map<G4String,DrawAttributes> EDepSim::RootGeometryManager::AttributeMap

private

Definition at line 38 of file EDepSimRootGeometryManager.hh.

Constructor & Destructor Documentation

EDepSim::RootGeometryManager::~RootGeometryManager ( )

virtual

Definition at line 69 of file EDepSimRootGeometryManager.cc.

69 {}

EDepSim::RootGeometryManager::RootGeometryManager ( )

protected

use Get() instead

Definition at line 61 of file EDepSimRootGeometryManager.cc.

61 {

62 }

Member Function Documentation

TGeoMedium * EDepSim::RootGeometryManager::AverageMaterial ( const G4VPhysicalVolume * theVol )

privatevirtual

Create a medium that describes the average composition for all the materials in the logical volume. This is used when the geometry is simplified by not including sub-volume (i.e. when the internal structure of the scintillating bars is skipped). The averaged material name is based on the name provided in materialBase;

Definition at line 984 of file EDepSimRootGeometryManager.cc.

                                       {
     G4LogicalVolume* theLog = thePhys->GetLogicalVolume();
     double totalMass = theLog->GetMass(true);
     double totalVolume = theLog->GetSolid()->GetCubicVolume();
     double totalDensity = totalMass/totalVolume;
     std::ostringstream nameStream;
     nameStream << MaterialName(thePhys)
                << "Avg" << totalDensity/(CLHEP::g/CLHEP::cm3);
     std::string averageName = nameStream.str();
     TGeoMedium* avgMedium = fMedium[averageName];
     if (avgMedium) return avgMedium;
     // Create a map of the material pointer and the mass of that material that
     // is contributing to the mixture.  This will be used to calculate the
     // average radiation length and interaction length.
     std::map< const G4Material* , double > materialMass;
     // Create a map of element names and the amount of that element in the
     // mixture.
     std::map<std::string,double> componentMass;
     // Create a stack to walk the geometry tree.
     std::vector<const G4VPhysicalVolume*> stack;
     stack.push_back(thePhys);
     while (!stack.empty()) {
         // Get the new current volume off the stack.
         const G4VPhysicalVolume* currentPhys = stack.back();
         G4LogicalVolume* currentLog = currentPhys->GetLogicalVolume();
         stack.pop_back();
         // Add all of the current children to the stack.
         for (std::size_t child = 0;
              child < currentLog->GetNoDaughters();
              ++child) {
             stack.push_back(currentLog->GetDaughter(child));
         }
         // Add the mass of each element in the current volume.
         double mass = currentLog->GetMass(true,false);
         G4Material *mat = currentLog->GetMaterial();
         materialMass[mat] += mass;
         for (unsigned i = 0; i < mat->GetNumberOfElements(); ++i) {
             const G4Element *elem = mat->GetElement(i);
             for (unsigned j = 0; j < elem->GetNumberOfIsotopes(); ++j) {
                 // Find the isotope for this element and create it if it
                 // doesn't exist
                 std::string isoName = elem->GetIsotope(j)->GetName();
                 componentMass[isoName] += mass
                     *(mat->GetFractionVector()[i])
                     *(elem->GetRelativeAbundanceVector()[j]);
             }
         }
     }
     TGeoMixture *theMix
         = new TGeoMixture(averageName.c_str(),
                           componentMass.size(),
                           totalDensity);
     for (std::map<std::string,double>::iterator c = componentMass.begin();
          c != componentMass.end();
          ++c) {
         TGeoElement* elem = fIsotope[c->first];
         theMix->AddElement(elem,c->second/totalMass);
     }
     double radLen = 0;
     double intLen = 0;
     double massSum = 0;
     for (std::map<const G4Material*,double>::iterator
              mat = materialMass.begin();
          mat != materialMass.end();
          ++mat) {
         const G4Material* material = mat->first;
         double mass = mat->second;
         massSum += mass;
         radLen += mass/material->GetRadlen();
         intLen += mass/material->GetNuclearInterLength();
     }
     if (massSum>0.0) {
         // The minus sign is to make sure this over-rides the approximate ROOT
         // radiation and interaction length calculations.
         theMix->SetRadLen(-massSum/radLen,-massSum/intLen);
     }
     int numed = fMedium.size() + 1;
     TGeoMedium* theMedium = new TGeoMedium(averageName.c_str(),numed,theMix);
     fMedium[averageName] = theMedium;
     return theMedium;
 }

bool EDepSim::RootGeometryManager::CreateEnvelope	(	const G4VPhysicalVolume *	theVol,
		TGeoManager *	theEnvelope,
		TGeoVolume *	theMother
	)

private

Save the detector envelope. This is called recursively to fill in the entire detector. The G4 physical volume, theVol, is used to create a new root TGeoVolume which is added to the existing root TGeoVolume, theMother, as a daughter.

RETURN VALUE: This returns true if one of the immediate daughter volumes was ignored. This tells the mother volume that it will need to adjust it's mass and material to account for the missing detail.

Definition at line 620 of file EDepSimRootGeometryManager.cc.

                            {
 
     if (IgnoreVolume(theG4PhysVol)) return true;
 
     // The new volume that will be added to the mother volume.  This is
     // created in this function.
     TGeoVolume* theVolume = NULL;
     G4LogicalVolume* theLog = theG4PhysVol->GetLogicalVolume();
 
     if (PrintMass(theG4PhysVol)) {
         EDepSimLog("%%% Mass: "
                    << G4BestUnit(theLog->GetMass(true),"Mass")
                    << " Volume: " << theG4PhysVol->GetName());
     }
 
     // Get the name of the expected name of the volume.
     std::string theVolumeName;
     for (std::vector<G4String>::iterator n = fNameStack.begin();
          n != fNameStack.end();
          ++n) {
         theVolumeName += "/";
         theVolumeName += *n;
     }
 
     // Get the volume information from G4.
     const G4VSolid* theSolid = theLog->GetSolid();
     const G4String geometryType = theSolid->GetEntityType();
     std::string theFullName = theG4PhysVol->GetName();
     std::string theShortName = theFullName;
     theShortName.erase(0,theShortName.rfind("/")+1);
 
     // The volume name is built from all the short names in the path to this
     // volume.  This is different than the full name which is built based on
     // the constructor tree.
     theVolumeName += "/";
     theVolumeName += theShortName;
 
     // Check the volume names for validity.  If they don't match then just
     // force it.
     if (theShortName == theFullName) {
         theFullName = theVolumeName;
     }
 
     // Make sure there isn't a bug in the naming.  There should never be a
     // "//" in the name.
     if (theFullName.find("//") != std::string::npos) {
         EDepSimError("Invalid volume name: " << theFullName);
         EDepSimError("   Expected name is: " << theVolumeName);
         theFullName = theVolumeName;
     }
 
     // Find the medium for this volume.
     std::string materialName = theLog->GetMaterial()->GetName();
     TGeoMedium *theMedium = fMedium[materialName];
     if (!theMedium) {
         EDepSimError("MISSING MATERIAL IS " << materialName);
         EDepSimThrow("Material definition is missing");
     }
 
     if (!fCreateAllVolumes) {
         std::map<G4LogicalVolume*, TGeoVolume*>::iterator seenVolume;
         seenVolume = fKnownVolumes.find(theLog);
         if (seenVolume != fKnownVolumes.end()) theVolume = seenVolume->second;
     }
 
     if (!theVolume) {
         // Create the root volume (the solid in G4 lingo).
         theVolume = CreateVolume(theSolid, theShortName, theMedium);
         if (!theVolume) theVolume = theMother;
         theVolume->SetTitle(theVolumeName.c_str());
 
         // Set the visual properties of the new volume
         int color = GetColor(theG4PhysVol);
         double opacity = GetOpacity(theG4PhysVol);
         theVolume->SetVisContainers(true);
         if (color < 0 || opacity < 0.001) {
             theVolume->SetVisibility(false);
         }
         else {
             int i = 100.0*(1.0-opacity);
             if (i>100) i = 100;
             EDepSimInfo("Set color of " << theShortName
                        << " to " << color
                        << " " << opacity
                        << " " << i);
             theVolume->SetLineColor(color);
             theVolume->SetTransparency(i);
             theVolume->SetVisibility(true);
         }
 
         // There is no mother so set this as the top volume.
         if (!theMother) {
             EDepSimLog("Making \"" << theVolume->GetName() << "\" the top\n");
             theEnvelope->SetTopVolume(theVolume);
         }
 
         // Push the name of this volume onto the stack before creating the
         // children.
         fNameStack.push_back(theShortName);
 
         // Add the children to the daughter.
         double missingMass = 0.0;
         for (std::size_t child = 0;
              child < theLog->GetNoDaughters();
              ++child) {
             G4VPhysicalVolume* theChild = theLog->GetDaughter(child);
             if (theLog->GetNoDaughters() > 20000) {
                 G4LogicalVolume *skippedVolume = theChild->GetLogicalVolume();
                 missingMass += skippedVolume->GetMass(true);
             }
             else if (CreateEnvelope(theChild, theEnvelope, theVolume)) {
                 G4LogicalVolume *skippedVolume = theChild->GetLogicalVolume();
                 missingMass += skippedVolume->GetMass(true);
             }
         }
         // Remove this volume from the name stack.
         fNameStack.pop_back();
 
         // A some daughter volumes were ignored and so the density of this
         // volume needs to be adjusted.
         if (missingMass > 0.0) {
             // The correction of the material needs to be implemented.
             double totalMass = theLog->GetMass(true);
             double totalVolume = theLog->GetSolid()->GetCubicVolume();
             double totalDensity = totalMass/totalVolume;
             EDepSimNamedDebug(
                 "ROOTGeom", "Skipping sub-volumes. Correct "
                 << theMedium->GetName() << " density "
                 << theMedium->GetMaterial()->GetDensity()/(CLHEP::g/CLHEP::cm3)
                 << " g/cm3 to "
                 << totalDensity/(CLHEP::g/CLHEP::cm3) << " g/cm3");
             TGeoMedium *avgMedium = AverageMaterial(theG4PhysVol);
             if (avgMedium) theVolume->SetMedium(avgMedium);
         }
 
         // Put this volume into the known volumes.
         fKnownVolumes[theLog] = theVolume;
     }
 
     // Apply the rotation for theMother volume.  This is only done if the
     // theMother is not the top level volume for the envelope so that the
     // detector has the proper orientation.
     if (theMother && (theMother != theVolume)) {
         G4ThreeVector trans = theG4PhysVol->GetObjectTranslation();
         G4RotationMatrix* rot = theG4PhysVol->GetObjectRotation();
         TGeoRotation* rotate =
             new TGeoRotation("rot",
                              TMath::RadToDeg()*rot->thetaX(),
                              TMath::RadToDeg()*rot->phiX(),
                              TMath::RadToDeg()*rot->thetaY(),
                              TMath::RadToDeg()*rot->phiY(),
                              TMath::RadToDeg()*rot->thetaZ(),
                              TMath::RadToDeg()*rot->phiZ());
         if (theG4PhysVol->IsReplicated()) {
             EAxis a; G4int nRep; G4double w; G4double o; G4bool c;
             G4ThreeVector axis;
             theG4PhysVol->GetReplicationData(a,nRep,w,o,c);
             switch (a) {
             case kXAxis: axis.set(1,0,0); break;
             case kYAxis: axis.set(0,1,0); break;
             case kZAxis: axis.set(0,0,1); break;
             default:
                 EDepSimThrow("EDepSim::RootGeometryManager::CreateEnvelope:"
                             "Bad replication data.");
             }
             for (int i=0; i<nRep; ++i) {
                 G4ThreeVector pos = axis*w*(i-0.5*(nRep-1));
                 TGeoCombiTrans* combi
                     = new TGeoCombiTrans(pos.x()/CLHEP::mm,
                                          pos.y()/CLHEP::mm,
                                          pos.z()/CLHEP::mm,
                                          rotate);
                 int index = HowManySimilarNodesInVolume(theMother,
                                                         theVolume->GetName());
                 theMother->AddNode(theVolume,index,combi);
             }
         }
         else {
             TGeoCombiTrans* combi
                 = new TGeoCombiTrans(trans.x()/CLHEP::mm,
                                      trans.y()/CLHEP::mm,
                                      trans.z()/CLHEP::mm,
                                      rotate);
             int index = HowManySimilarNodesInVolume(theMother,
                                                     theVolume->GetName());
             theMother->AddNode(theVolume,index,combi);
         }
     }
 
     return false;
 }

void EDepSim::RootGeometryManager::CreateMaterials ( const G4VPhysicalVolume * theVol )

private

Create the materials needed for the detector. This called recursively to find all of the materials in the detector.

Definition at line 518 of file EDepSimRootGeometryManager.cc.

                                            {
 
     G4LogicalVolume* theLog = theG4PhysVol->GetLogicalVolume();
 
     // Find the medium for this volume and create it if it doesn't exist.
     TGeoMedium *theMedium = fMedium[theLog->GetMaterial()->GetName()];
 
     if (!theMedium) {
         G4Material *mat = theLog->GetMaterial();
         if (mat->GetNumberOfElements()==0) {
             EDepSimError("Material without elements " << mat->GetName());
             EDepSimThrow("Material defined without elements");
         }
         // Make a mixture
         TGeoMixture *theMix
             = new TGeoMixture(mat->GetName(),
                               mat->GetNumberOfElements(),
                               mat->GetDensity());
         unsigned ielement = 0;
         theMix->SetTemperature(mat->GetTemperature());
         // The minus sign is to make sure this over-rides the approximate ROOT
         // radiation and interaction length calculations.
         switch (mat->GetState()) {
         case kStateSolid:
             theMix->SetState(TGeoMaterial::kMatStateSolid);
             break;
         case kStateLiquid:
             theMix->SetState(TGeoMaterial::kMatStateLiquid);
             theMix->SetPressure(mat->GetPressure());
             break;
         case kStateGas:
             theMix->SetState(TGeoMaterial::kMatStateGas);
             theMix->SetPressure(mat->GetPressure());
             break;
         default:
             theMix->SetState(TGeoMaterial::kMatStateUndefined);
         }
         for (unsigned i = 0; i < mat->GetNumberOfElements(); ++i) {
             const G4Element *elem = mat->GetElement(i);
             for (unsigned j = 0; j < elem->GetNumberOfIsotopes(); ++j) {
                 const G4Isotope *iso = elem->GetIsotope(j);
                 theMix->DefineElement(ielement,
                                       iso->GetA()/(CLHEP::g/CLHEP::mole),
                                       iso->GetZ(),
                                       (mat->GetFractionVector()[i])*
                                       (elem->GetRelativeAbundanceVector()[j]));
                 // Find the isotope for this element and create it if it
                 // doesn't exist
                 TGeoElement *theIsotope = fIsotope[iso->GetName()];
                 if (!theIsotope) {
                     // Make an element (isotope)
                     TGeoElement *theEle
                         = new TGeoElement(iso->GetName(),
                                           elem->GetSymbol(),
                                           iso->GetZ(),
                                           iso->GetA()/(CLHEP::g/CLHEP::mole));
                     fIsotope[iso->GetName()] = theEle;
                 }
                 ielement++;
             }
         }
         theMix->SetRadLen(-mat->GetRadlen(), -mat->GetNuclearInterLength());
         int numed = fMedium.size() + 1;
         theMedium = new TGeoMedium(mat->GetName(),numed,theMix);
         fMedium[mat->GetName()] = theMedium;
     }
 
     // Recurse through the children.
     for (std::size_t child = 0;
          child < theLog->GetNoDaughters();
          ++child) {
         G4VPhysicalVolume* theChild = theLog->GetDaughter(child);
         CreateMaterials(theChild);
     }
 }

TGeoShape * EDepSim::RootGeometryManager::CreateShape	(	const G4VSolid *	theSolid,
		TGeoMatrix **	mat = `NULL`
	)

private

Create a new ROOT shape object based on the G4 solid. This might be called recursively if a G4 boolean solid is found.

Definition at line 224 of file EDepSimRootGeometryManager.cc.

                                                                             {
     const G4String geometryType = theSolid->GetEntityType();
     TGeoShape* theShape = NULL;
     if (geometryType == "G4Box") {
         // Create a box
         const G4Box* box = dynamic_cast<const G4Box*>(theSolid);
         theShape = new TGeoBBox(box->GetXHalfLength()/CLHEP::mm,
                                 box->GetYHalfLength()/CLHEP::mm,
                                 box->GetZHalfLength()/CLHEP::mm);
     }
     else if (geometryType == "G4Tubs") {
         const G4Tubs* tube = dynamic_cast<const G4Tubs*>(theSolid);
         // Root takes the angles in degrees so there is no extra
         // conversion.
         double zhalf = tube->GetZHalfLength()/CLHEP::mm;
         double rmin = tube->GetInnerRadius()/CLHEP::mm;
         double rmax = tube->GetOuterRadius()/CLHEP::mm;
         double minPhiDeg = tube->GetStartPhiAngle()/CLHEP::degree;
         double maxPhiDeg = minPhiDeg + tube->GetDeltaPhiAngle()/CLHEP::degree;
         theShape = new TGeoTubeSeg(rmin, rmax,
                                    zhalf,
                                    minPhiDeg, maxPhiDeg);
     }
     else if (geometryType == "G4Sphere") {
         const G4Sphere* sphere = dynamic_cast<const G4Sphere*>(theSolid);
         // Root takes the angles in degrees so there is no extra
         // conversion.
         double minPhiDeg = sphere->GetStartPhiAngle()/CLHEP::degree;
         double maxPhiDeg = minPhiDeg + sphere->GetDeltaPhiAngle()/CLHEP::degree;
         double minThetaDeg = sphere->GetStartThetaAngle()/CLHEP::degree;
         double maxThetaDeg = minThetaDeg + sphere->GetDeltaThetaAngle()/CLHEP::degree;
         theShape = new TGeoSphere(sphere->GetInnerRadius()/CLHEP::mm,
                                   sphere->GetOuterRadius()/CLHEP::mm,
                                   minThetaDeg, maxThetaDeg,
                                   minPhiDeg, maxPhiDeg);
     }
     else if (geometryType == "G4Polyhedra") {
         const G4Polyhedra* polyhedra
             = dynamic_cast<const G4Polyhedra*>(theSolid);
         double phi = polyhedra->GetStartPhi();
         double dPhi = polyhedra->GetEndPhi() - phi;
         if (dPhi>2*M_PI) dPhi -= 2*M_PI;
         if (dPhi<0) dPhi += 2*M_PI;
         int sides = polyhedra->GetNumSide();
         int numZ = polyhedra->GetNumRZCorner()/2;
         // Factor to take into account that ROOT uses the circle that can be
         // inscribed inside the polygon, and G4 uses the corner
         double g4Factor = std::cos(0.5*dPhi/sides);
         TGeoPgon* pgon = new TGeoPgon(phi/CLHEP::degree,
                                       dPhi/CLHEP::degree, sides, numZ);
         for (int i = 0; i< numZ; ++i) {
             double rMin = g4Factor*polyhedra->GetCorner(numZ-i-1).r;
             double rMax = g4Factor*polyhedra->GetCorner(numZ+i).r;
             if (rMax < rMin) std::swap(rMin,rMax);
             pgon->DefineSection(i,
                                 polyhedra->GetCorner(numZ-i-1).z/CLHEP::mm,
                                 rMin/CLHEP::mm,
                                 rMax/CLHEP::mm);
         }
         theShape = pgon;
     }
     else if (geometryType == "G4Polycone") {
         const G4Polycone* polycone
             = dynamic_cast<const G4Polycone*>(theSolid);
         double phi = polycone->GetStartPhi();
         double dPhi = polycone->GetEndPhi() - phi;
         if (dPhi>2*M_PI) dPhi -= 2*M_PI;
         if (dPhi<0) dPhi += 2*M_PI;
 #ifdef G4GEOM_USE_USOLIDS
 #warning GEANT HAS BEEN COMPILED WITH BROKEN USOLIDS.
         int numZ = polycone->GetNumRZCorner()/2;
         TGeoPcon* pcon = new TGeoPcon(phi/CLHEP::degree,
                                       dPhi/CLHEP::degree, numZ);
         // This depends on the (mostly) undocumented order of the corners in
         // the G4Polycone internals.  It's a little unstable...
         for (int i = 0; i< numZ; ++i) {
             pcon->DefineSection(i,
                                 polycone->GetCorner(numZ-i-1).z/CLHEP::mm,
                                 polycone->GetCorner(numZ-i-1).r/CLHEP::mm,
                                 polycone->GetCorner(numZ+i).r/CLHEP::mm);
         }
 #else
         const G4PolyconeHistorical* param = polycone->GetOriginalParameters();
         int numZ = param->Num_z_planes;
         TGeoPcon* pcon = new TGeoPcon(phi/CLHEP::degree, dPhi/CLHEP::degree, numZ);
         // This depends on the older interface.  It's not marked as
         // deprecated, but the documentation discourages it's use.
         for (int i = 0; i< numZ; ++i) {
             pcon->DefineSection(i,
                                 param->Z_values[i]/CLHEP::mm,
                                 param->Rmin[i]/CLHEP::mm,
                                 param->Rmax[i]/CLHEP::mm);
         }
 #endif
         theShape = pcon;
     }
     else if (geometryType == "G4Trap") {
         const G4Trap* trap
             = dynamic_cast<const G4Trap*>(theSolid);
         double dz = trap->GetZHalfLength()/CLHEP::mm;
         double theta = 0;
         double phi = 0;
         double h1 = trap->GetYHalfLength1()/CLHEP::mm;
         double bl1 = trap->GetXHalfLength1()/CLHEP::mm;
         double tl1 = trap->GetXHalfLength2()/CLHEP::mm;
         double alpha1 = std::atan(trap->GetTanAlpha1())/CLHEP::degree;
         double h2 = trap->GetYHalfLength2()/CLHEP::mm;
         double bl2 = trap->GetXHalfLength3()/CLHEP::mm;
         double tl2 = trap->GetXHalfLength4()/CLHEP::mm;
         double alpha2 = std::atan(trap->GetTanAlpha2())/CLHEP::degree;
         theShape = new TGeoTrap(dz, theta, phi,
                                 h1, bl1, tl1, alpha1,
                                 h2, bl2, tl2, alpha2);
     }
     else if (geometryType == "G4Trd") {
         const G4Trd* trd
             = dynamic_cast<const G4Trd*>(theSolid);
         double dz = trd->GetZHalfLength()/CLHEP::mm;
         double dx1 = trd->GetXHalfLength1()/CLHEP::mm;
         double dx2 = trd->GetXHalfLength2()/CLHEP::mm;
         double dy1 = trd->GetYHalfLength1()/CLHEP::mm;
         double dy2 = trd->GetYHalfLength2()/CLHEP::mm;
         theShape = new TGeoTrd2(dx1,dx2,dy1,dy2,dz);
     }
     else if (geometryType == "G4SubtractionSolid") {
         const G4SubtractionSolid* sub
             = dynamic_cast<const G4SubtractionSolid*>(theSolid);
         const G4VSolid* solidA = sub->GetConstituentSolid(0);
         const G4VSolid* solidB = sub->GetConstituentSolid(1);
         // solidA - solidB
         TGeoMatrix* matrixA = NULL;
         TGeoShape* shapeA = CreateShape(solidA, &matrixA);
         TGeoMatrix* matrixB = NULL;
         TGeoShape* shapeB = CreateShape(solidB, &matrixB);
         TGeoSubtraction* subtractNode = new TGeoSubtraction(shapeA,shapeB,
                                                             matrixA, matrixB);
         theShape = new TGeoCompositeShape("name",subtractNode);
     }
     else if (geometryType == "G4DisplacedSolid") {
         const G4DisplacedSolid* disp
             = dynamic_cast<const G4DisplacedSolid*>(theSolid);
         const G4VSolid* movedSolid = disp->GetConstituentMovedSolid();
         G4RotationMatrix rotation = disp->GetObjectRotation();
         G4ThreeVector displacement = disp->GetObjectTranslation();
         theShape = CreateShape(movedSolid);
         if (returnMatrix) {
             TGeoRotation* rotate
                 = new TGeoRotation("rot",
                                    TMath::RadToDeg()*rotation.thetaX(),
                                    TMath::RadToDeg()*rotation.phiX(),
                                    TMath::RadToDeg()*rotation.thetaY(),
                                    TMath::RadToDeg()*rotation.phiY(),
                                    TMath::RadToDeg()*rotation.thetaZ(),
                                    TMath::RadToDeg()*rotation.phiZ());
             *returnMatrix = new TGeoCombiTrans(displacement.x()/CLHEP::mm,
                                                displacement.y()/CLHEP::mm,
                                                displacement.z()/CLHEP::mm,
                                                rotate);
         }
     }
     else if (geometryType == "G4UnionSolid") {
         const G4UnionSolid* sub
             = dynamic_cast<const G4UnionSolid*>(theSolid);
         const G4VSolid* solidA = sub->GetConstituentSolid(0);
         const G4VSolid* solidB = sub->GetConstituentSolid(1);
         // solidA - solidB
         TGeoMatrix* matrixA = NULL;
         TGeoShape* shapeA = CreateShape(solidA, &matrixA);
         TGeoMatrix* matrixB = NULL;
         TGeoShape* shapeB = CreateShape(solidB, &matrixB);
         TGeoUnion* unionNode = new TGeoUnion(shapeA,  shapeB,
                                              matrixA, matrixB);
         theShape = new TGeoCompositeShape("name",unionNode);
     }
     else if (geometryType == "G4IntersectionSolid") {
         const G4IntersectionSolid* sub
           = dynamic_cast<const G4IntersectionSolid*>(theSolid);
         const G4VSolid* solidA = sub->GetConstituentSolid(0);
         const G4VSolid* solidB = sub->GetConstituentSolid(1);
         // solidA - solidB
         TGeoMatrix* matrixA = NULL;
         TGeoShape* shapeA = CreateShape(solidA, &matrixA);
         TGeoMatrix* matrixB = NULL;
         TGeoShape* shapeB = CreateShape(solidB, &matrixB);
         TGeoIntersection* intersectionNode = new TGeoIntersection(shapeA,  shapeB,
                                                                   matrixA, matrixB);
         theShape = new TGeoCompositeShape("name",intersectionNode);
     }
 
     else if (geometryType == "G4ExtrudedSolid"){
         //This following only works when using the 'standard'
         //G4ExtrudedSolid Constructor.
 
         const G4ExtrudedSolid* extr
             = dynamic_cast<const G4ExtrudedSolid*>(theSolid);
 
         //number of z planes
         const G4int nZ = extr->GetNofZSections();
         //number of vertices in the polygon
         const G4int nV = extr->GetNofVertices();
 
         //define and pointers
         const int maxVertices = 1000;
         double vertices_x[maxVertices];
         double vertices_y[maxVertices];
         if (maxVertices < nV) {
             EDepSimThrow("Polygon with more than maxVertices");
         }
 
 
         //define an intermediate extrusion constructor with nZ z planes.
         TGeoXtru *xtru = new TGeoXtru(nZ);
 
         //Get the polygons points.
         std::vector<G4TwoVector> polyPoints = extr->GetPolygon();
 
         //fill the vertices arrays
         for(int i = 0 ; i < nV ; i++){
             vertices_x[i]= polyPoints[i].x();
             vertices_y[i]= polyPoints[i].y();
         }
 
         //Define the polygon
         xtru->DefinePolygon(nV, vertices_x, vertices_y);
 
         double z_pos, x_off, y_off, scale;
 
         //fill the parameters to define the Root extruded solid
         for(int i = 0 ; i < nZ ; i++){
             z_pos = extr->GetZSection(i).fZ;
             x_off = extr->GetZSection(i).fOffset.x() ;
             y_off = extr->GetZSection(i).fOffset.y();
             scale = extr->GetZSection(i).fScale;
             xtru->DefineSection(i, z_pos, x_off, y_off, scale);
         }
         //now assign 'theShape' to this complete extruded object.
         theShape = xtru;
     }
     else {
         EDepSimThrow(geometryType+" --> shape not implemented");
     }
 
     return theShape;
 }

TGeoVolume * EDepSim::RootGeometryManager::CreateVolume	(	const G4VSolid *	theSolid,
		std::string	theName,
		TGeoMedium *	theMedium
	)

private

Create a new ROOT volume object.

Definition at line 471 of file EDepSimRootGeometryManager.cc.

                                                                   {
     TGeoShape* theShape = CreateShape(theSolid);
     TGeoVolume* theVolume = new TGeoVolume(theName.c_str(),
                                            theShape,
                                            theMedium);
     return theVolume;
 }

void EDepSim::RootGeometryManager::Export ( const char * file )

Export the geometry to a file.

Definition at line 71 of file EDepSimRootGeometryManager.cc.

                                                         {
     EDepSimLog("   *** Export to " << file);
     gGeoManager->Export(file);
 }

EDepSim::RootGeometryManager * EDepSim::RootGeometryManager::Get ( void )

static

If a persistency manager has not been created, create one.

Definition at line 64 of file EDepSimRootGeometryManager.cc.

                                                             {
     if (!fThis) fThis = new EDepSim::RootGeometryManager();
     return fThis;
 }

int EDepSim::RootGeometryManager::GetColor ( const G4VPhysicalVolume * vol )

protected

Find the right color for a physical volume. This won't necessarily be the same color as set in G4, but should be a reasonable choice for a detector picture. This will return a negative number if the volume should be invisible.

Definition at line 870 of file EDepSimRootGeometryManager.cc.

                                                                      {
     const G4LogicalVolume* log = vol->GetLogicalVolume();
     std::string theFullName = vol->GetName();
 
     const G4VisAttributes* visAttributes = log->GetVisAttributes();
     if (!visAttributes) {
         return -1;
     }
 
     int index = -1;
     double minDist = 10000;
     G4Color g4Color = visAttributes->GetColor();
 
     // Check the primary color wheel.
     const int rootPrimaryColors[] = {
         kRed, kMagenta, kBlue, kCyan, kGreen, kYellow, -1
     };
     for(const int* rootBaseColor = rootPrimaryColors;
         0 <= *rootBaseColor;
         ++rootBaseColor) {
         for (int i = -10; i<5; ++i) {
             int iColor = *rootBaseColor + i;
             TColor* rootColor = gROOT->GetColor(iColor);
             if (!rootColor) continue;
             double dR = (rootColor->GetRed() - g4Color.GetRed());
             double dG = (rootColor->GetGreen() - g4Color.GetGreen());
             double dB = (rootColor->GetBlue() - g4Color.GetBlue());
             double dist = std::sqrt(dR*dR + dG*dG + dB*dB);
             if (dist > minDist) continue;
             index = iColor;
             minDist = dist;
         }
     }
 
     // Check the secondary colors defined by ROOT.
     const int rootSecondaryColors[] = {
         kPink, kMagenta, kViolet, kTeal, kSpring, kOrange, -1
     };
     for(const int* rootBaseColor = rootSecondaryColors;
         0 <= *rootBaseColor;
         ++rootBaseColor) {
         for (int i = -9; i<=10; ++i) {
             int iColor = *rootBaseColor + i;
             TColor* rootColor = gROOT->GetColor(iColor);
             if (!rootColor) continue;
             double dR = (rootColor->GetRed() - g4Color.GetRed());
             double dG = (rootColor->GetGreen() - g4Color.GetGreen());
             double dB = (rootColor->GetBlue() - g4Color.GetBlue());
             double dist = std::sqrt(dR*dR + dG*dG + dB*dB);
             if (dist > minDist) continue;
             index = iColor;
             minDist = dist;
         }
     }
 
     // Check the different colors of "black"
     const int rootBlackColors[] = {
         kGray, kGray+1, kGray+2, kGray+3, kBlack, -1
     };
     for(const int* rootBaseColor = rootBlackColors;
         0 <= *rootBaseColor;
         ++rootBaseColor) {
         TColor* rootColor = gROOT->GetColor(*rootBaseColor);
         if (!rootColor) continue;
         double dR = (rootColor->GetRed() - g4Color.GetRed());
         double dG = (rootColor->GetGreen() - g4Color.GetGreen());
         double dB = (rootColor->GetBlue() - g4Color.GetBlue());
         double dist = std::sqrt(dR*dR + dG*dG + dB*dB);
         if (dist > minDist) continue;
         index = *rootBaseColor;
         minDist = dist;
     }
 
     // Check the basic colors (the first 50 colors in root).
     for(int rootBaseColor = 1;
         rootBaseColor < 50;
         ++rootBaseColor) {
         TColor* rootColor = gROOT->GetColor(rootBaseColor);
         if (!rootColor) continue;
         double dR = (rootColor->GetRed() - g4Color.GetRed());
         double dG = (rootColor->GetGreen() - g4Color.GetGreen());
         double dB = (rootColor->GetBlue() - g4Color.GetBlue());
         double dist = std::sqrt(dR*dR + dG*dG + dB*dB);
         if (dist > minDist) continue;
         index = rootBaseColor;
         minDist = dist;
     }
 
     return index;
 }

G4Color EDepSim::RootGeometryManager::GetG4Color ( G4Material * material )

Get a color for a material. This can be used to create the G4VisAttributes object.

Definition at line 827 of file EDepSimRootGeometryManager.cc.

                                                                    {
     G4String materialName = material->GetName();
     AttributeMap::const_iterator colorPair = fColorMap.find(materialName);
 
     int index = 0;
     double alpha = 1.0;
 
     if (colorPair == fColorMap.end()) {
         EDepSimError("Missing color for \"" << materialName << "\"");
         fColorMap[materialName].color = kOrange+1;
         fColorMap[materialName].fill = 0;
         index = kOrange+1;
     }
     else {
         index = colorPair->second.color;
         alpha = colorPair->second.alpha;
     }
 
 #include "rootColorToG4ColorMap.hxx"
 
     G4Color color = sRootColorToG4ColorMap[index];
 
     return G4Color(color.GetRed(),color.GetGreen(),color.GetBlue(),
                    color.GetAlpha()*alpha);
 }

int EDepSim::RootGeometryManager::GetNodeId ( const G4ThreeVector & pos )

Get a volume ID base on the volume position.

Definition at line 76 of file EDepSimRootGeometryManager.cc.

                                                                   {
     if (!gGeoManager) return -1;
     gGeoManager->FindNode(pos.x(), pos.y(), pos.z());
     return gGeoManager->GetCurrentNodeId();
 }

double EDepSim::RootGeometryManager::GetOpacity ( const G4VPhysicalVolume * vol )

protected

Get the opacity of the physical volume. This will return an opacity of zero (i.e. transparent) if the volume should be invisible.

Definition at line 853 of file EDepSimRootGeometryManager.cc.

                                                                           {
     const G4LogicalVolume* log = vol->GetLogicalVolume();
     std::string theFullName = vol->GetName();
 
     const G4VisAttributes* visAttributes = log->GetVisAttributes();
     if (!visAttributes) {
         // Invisible if the visual attributes are missing...
         return 0.0;
     }
 
     G4Color g4Color = visAttributes->GetColor();
     double opacity = g4Color.GetAlpha();
     if (opacity > 1.0) opacity = 1.0;
 
     return opacity;
 }

int EDepSim::RootGeometryManager::HowManySimilarNodesInVolume	(	TGeoVolume *	theMother,
		std::string	theName
	)

private

Definition at line 597 of file EDepSimRootGeometryManager.cc.

                                                    {
     daughterName = daughterName + "_";
     int ndaughters = theMother->GetNdaughters();
 
     int ndaughters_samename = 0;
     for(int i = 0; i < ndaughters; i++){
         TGeoNode *node = theMother->GetNode(i);
         if(node){
             std::string name(node->GetName());
             if(name.find(daughterName.c_str()) != std::string::npos)
                 ndaughters_samename++;
         }
     }
 
     return ndaughters_samename;
 
 }

bool EDepSim::RootGeometryManager::IgnoreVolume ( const G4VPhysicalVolume * theVol )

privatevirtual

A method to prune the geometry that gets exported from G4 into root. If this returns true, the volume and all of it's children will not be exported to ROOT. This allows the internal structure of low level objects (i.e. the internal "structure" of scintillating bars) to be hidden.

Definition at line 484 of file EDepSimRootGeometryManager.cc.

                                                                              {
     std::string theFullName = theVol->GetName();
     std::string theShortName = theFullName;
     theShortName.erase(0,theShortName.rfind("/")+1);
 
     // Don't save the internal structure of extruded scintillating bars.  This
     // is required to make sure that hits get assigned to the right geometry
     // volume in the ROOT geometry.
     if (theFullName.find("Bar/Core") != std::string::npos) return true;
     if (theFullName.find("Bar/CrnrOfBar") != std::string::npos) return true;
     if (theFullName.find("Bar/SideOfBar") != std::string::npos) return true;
 
     return false;
 }

std::string EDepSim::RootGeometryManager::MaterialName ( const G4VPhysicalVolume * theVol )

privatevirtual

A method to translate from the local volume material to the right material for the root geometry. This is required for volumes where the inner structure is not saved, and the material of the outer volume does not provide a good base name.

Definition at line 961 of file EDepSimRootGeometryManager.cc.

                                       {
     std::string theFullName = thePhys->GetName();
     G4LogicalVolume* theLog = thePhys->GetLogicalVolume();
     std::string materialName = theLog->GetMaterial()->GetName();
 
     // If this is a scintillator bar, then the base material name should be
     // either FGDScintillator or Scintillator.
     if (theFullName.find("/Bar") != std::string::npos) {
         if (theFullName.find("/FGD") != std::string::npos) {
             materialName = "FGDScintillator";
         }
         else if (theFullName.find("/P0D/") != std::string::npos) {
             materialName = "P0DScintillator";
         }
         else {
             materialName = "Scintillator";
         }
     }
 
     return materialName;
 }

bool EDepSim::RootGeometryManager::PrintMass ( const G4VPhysicalVolume * theVol )

privatevirtual

A method to flag that a volume mass should be printed.

Definition at line 500 of file EDepSimRootGeometryManager.cc.

                                                                           {
     std::string theFullName = theVol->GetName();
 
     for (std::vector<G4String>::iterator n = fPrintMass.begin();
          n != fPrintMass.end();
          ++n) {
         if (theFullName.find(*n) != std::string::npos) {
             if (fPrintedMass.find(*n) != fPrintedMass.end()) continue;
             fPrintedMass.insert(*n);
             return true;
         }
     }
 
     return false;
 }

void EDepSim::RootGeometryManager::SetDrawAtt	(	G4Material *	material,
		int	color,
		double	opacity = `1.0`
	)

Set material color.

Definition at line 815 of file EDepSimRootGeometryManager.cc.

                                                                          {
     G4String materialName = material->GetName();
     fColorMap[materialName].color = color;
     fColorMap[materialName].alpha = opacity;
     if (opacity<0.01) fColorMap[materialName].fill = 0;
     else if (opacity>0.99) fColorMap[materialName].fill = 1;
     else {
         fColorMap[materialName].fill = 4000+100*opacity;
     }
 }

virtual void EDepSim::RootGeometryManager::ShouldPrintMass ( std::string name )

inlinevirtual

A method to request that a volume mass should be printed;.

Definition at line 69 of file EDepSimRootGeometryManager.hh.

                                                  {
         fPrintMass.push_back(name);
     }

void EDepSim::RootGeometryManager::Update	(	const G4VPhysicalVolume *	aWorld,
		bool	validate
	)

Update the root geometry to match the g4 geometry.

Definition at line 93 of file EDepSimRootGeometryManager.cc.

                                                                  {
     EDepSimLog("%%%%%%%%%%%%%%%%%%%%%%%%%% Update ROOT Geometry "
              << "%%%%%%%%%%%%%%%%%%%%%%%%%%" );
     if (gGeoManager) {
         EDepSimLog("%%%%%%%%%%%%%%% Warning: Replacing ROOT Geometry ");
         delete gGeoManager;
         // Clear the node counts definitions.
         fNodeCount.clear();
         // Clear the cached materials.  The material objects are owned by the
         // TGeoManager and are invalid after it has been deleted.
         fMedium.clear();
         // Clear the cached isotopes.  This is a bit complicated since the
         // isotope pointers may be duplicated in the map.
         std::set<TGeoElement*> isotope;
         for (std::map<G4String,TGeoElement*>::iterator i = fIsotope.begin();
              i != fIsotope.end();
              ++i) {
             isotope.insert(i->second);
         }
         fIsotope.clear();
         for (std::set<TGeoElement*>::iterator i = isotope.begin();
              i != isotope.end();
              ++i) {
             delete *i;
         }
     }
     // Create the new geometry.
     gGeoManager = new TGeoManager("EDepSimGeometry",
                                   "Simulated Detector Geometry");
     gGeoManager->SetVisLevel(20);
     // Create all of the materials.
     EDepSimInfo("Create materials");
     CreateMaterials(aWorld);
 
     //Check to see if we can create all of the volumes.  This is done by
     //traversing the GEANT physical volume tree.
     fCreateAllVolumes = false;
     if (CountVolumes(aWorld->GetLogicalVolume()) < 100000) {
         EDepSimInfo("Create all volumes");
         fCreateAllVolumes = true;
     }
 
     // Create the ROOT geometry definitions.
     fPrintedMass.clear();
     fNameStack.clear();
     fKnownVolumes.clear();
     EDepSimInfo("Start defining envelope");
     CreateEnvelope(aWorld,gGeoManager,NULL);
     EDepSimInfo("Geometry is Filled");
 
     gGeoManager->CloseGeometry("di");
 
     EDepSimLog("Geometry initialized and closed");
 
     if (validateGeometry) return Validate();
 }

void EDepSim::RootGeometryManager::Update	(	std::string	gdmlFile,
		bool	validate
	)

Set the root geometry from a gdml file.

Definition at line 151 of file EDepSimRootGeometryManager.cc.

                                                                  {
     EDepSimLog(   "%%%%%%%%%%%%%%%%%% Update ROOT Geometry from GDML"
                << "%%%%%%%%%%%%%%%%%%" );
     if (gGeoManager) {
         EDepSimLog("%%%%%%%%%%%%%%% Warning: Replacing ROOT Geometry ");
         delete gGeoManager;
         // Clear the node counts definitions.
         fNodeCount.clear();
         // Clear the cached materials.  The material objects are owned by the
         // TGeoManager and are invalid after it has been deleted.
         fMedium.clear();
         // Clear the cached isotopes.  This is a bit complicated since the
         // isotope pointers may be duplicated in the map.
         std::set<TGeoElement*> isotope;
         for (std::map<G4String,TGeoElement*>::iterator i = fIsotope.begin();
              i != fIsotope.end();
              ++i) {
             isotope.insert(i->second);
         }
         fIsotope.clear();
         for (std::set<TGeoElement*>::iterator i = isotope.begin();
              i != isotope.end();
              ++i) {
             delete *i;
         }
     }
 
     // Create the new geometry.
     TGeoManager::Import(gdmlFile.c_str());
     gGeoManager->SetName("EDepSimGeometry");
     gGeoManager->SetTitle("Simulated Detector Geometry");
 
     EDepSimLog("####################### Geometry initialized and closed");
 
     if (validateGeometry) return Validate();
 }

void EDepSim::RootGeometryManager::Validate ( )

Make sure that the current geometry passes a bunch of tests.

Definition at line 189 of file EDepSimRootGeometryManager.cc.

                                           {
 
     int count = 0;
     // Check for overlaps at volume corners.  Look for overlaps at 0.1 mm size.
     gGeoManager->CheckOverlaps(0.1*CLHEP::mm);
     {
         TIter next(gGeoManager->GetListOfOverlaps());
         TGeoOverlap* overlap;
         while ((overlap=(TGeoOverlap*)next())) {
             ++count;
             overlap->PrintInfo();
         }
     }
 
     // Check for overlaps with sampling.  Look for overlaps at 0.1 mm size.
     gGeoManager->CheckOverlaps(0.1*CLHEP::mm,"s100000");
     {
         TIter next(gGeoManager->GetListOfOverlaps());
         TGeoOverlap* overlap;
         while ((overlap=(TGeoOverlap*)next())) {
             ++count;
             overlap->PrintInfo();
         }
     }
 
     if (count > 0) {
         EDepSimThrow(
             "The geometry has overlaps and will produce incorrect"
             " results.  To use the geometry, specify the '-C' option"
             " on the command line.");
     }
 
     EDepSimLog("Geometry validated");
 }

Member Data Documentation

AttributeMap EDepSim::RootGeometryManager::fColorMap

private

A map between a material name and a color.

Definition at line 101 of file EDepSimRootGeometryManager.hh.

bool EDepSim::RootGeometryManager::fCreateAllVolumes

private

A flag whether to brute force the geometry or not. If this is true, then a new volume is created for every object in the geometry. This means that replicated volumes, and duplicated volume trees occur multiple times. This is good for small detectors since it means that each volume in the detector will have a unique TGeoNode associated with it. It is bad for large detectors since the ROOT geometry will become impractially large. The decision is made based on traversing the geant geometry and counting the number of unique volumes.

Definition at line 126 of file EDepSimRootGeometryManager.hh.

std::map<G4String,TGeoElement*> EDepSim::RootGeometryManager::fIsotope

private

A map between G4 isotope names and Root Element definitions.

Definition at line 95 of file EDepSimRootGeometryManager.hh.

std::map<G4LogicalVolume*,TGeoVolume*> EDepSim::RootGeometryManager::fKnownVolumes

private

A map between the G4 logical volumes and the root volumes. This keeps track of which volumes have already been created. It needs to be cleared before starting to build a new root geometry.

Definition at line 116 of file EDepSimRootGeometryManager.hh.

std::map<G4String,TGeoMedium*> EDepSim::RootGeometryManager::fMedium

private

A map between G4 material names and Root Material definitions.

Definition at line 92 of file EDepSimRootGeometryManager.hh.

std::vector<G4String> EDepSim::RootGeometryManager::fNameStack

private

A stack of volume names that have been seen. These are the "short" volume names of all of the parents to the current volue.

Definition at line 111 of file EDepSimRootGeometryManager.hh.

std::map<G4String,int> EDepSim::RootGeometryManager::fNodeCount

private

A map between G4 volume names and count of volumes with that name.