diff --git a/DDCore/src/DetectorImp.cpp b/DDCore/src/DetectorImp.cpp index 7cc325f5006d67ce36478ee0a5d36cb04548e13d..18ca12268eb98bf642eb99486ff15d54a6818052 100644 --- a/DDCore/src/DetectorImp.cpp +++ b/DDCore/src/DetectorImp.cpp @@ -670,6 +670,9 @@ void DetectorImp::endDocument(bool close_geometry) { add(trackingVis); #endif m_worldVol.solid()->ComputeBBox(); + // Propagating reflections + ReflectionBuilder rb(*this); + rb.execute(); // Since we allow now for anonymous shapes, // we will rename them to use the name of the volume they are assigned to mgr->CloseGeometry(); @@ -677,9 +680,6 @@ void DetectorImp::endDocument(bool close_geometry) { // Patching shape names of anaonymous shapes ShapePatcher patcher(m_volManager, m_world); patcher.patchShapes(); - // Propagating reflections - ReflectionBuilder rb(*this); - rb.execute(); mapDetectorTypes(); m_state = READY; diff --git a/DDCore/src/Volumes.cpp b/DDCore/src/Volumes.cpp index dd23f31eab01c0e639e00da245ef1f77fd9c68e9..0bfeaf13d686b887fbb5926ffa46a8e9769afc4b 100644 --- a/DDCore/src/Volumes.cpp +++ b/DDCore/src/Volumes.cpp @@ -185,6 +185,7 @@ namespace { static TMap map(100); TGeoVolume* vol = (TGeoVolume*)map.GetValue(v); if ( vol ) { + if (newname && newname[0]) v->SetName(newname); return vol; } vol = v->CloneVolume(); @@ -203,19 +204,17 @@ namespace { vol->SetName((nam+"_refl").c_str()); } delete vol->GetNodes(); - vol->SetNodes(NULL); + vol->SetNodes(nullptr); vol->SetBit(TGeoVolume::kVolumeImportNodes, kFALSE); v->CloneNodesAndConnect(vol); // The volume is now properly cloned, but with the same shape. // Reflect the shape (if any) and connect it. -#if 0 if (v->GetShape()) { TGeoScale* scale = new TGeoScale( 1., 1.,-1.); TGeoShape* reflected_shape = TGeoScaledShape::MakeScaledShape((nam+"_shape_refl").c_str(), v->GetShape(), scale); vol->SetShape(reflected_shape); } -#endif // Reflect the daughters. Int_t nd = vol->GetNdaughters(); if ( !nd ) return vol; @@ -233,8 +232,8 @@ namespace { if (!reflected) { // We need to reflect only the translation and propagate to daughters. // H' = Sz * H * Sz - local_cloned->ReflectZ(kTRUE,kFALSE); - local_cloned->ReflectZ(kFALSE,kFALSE); + local_cloned->ReflectZ(kTRUE); + local_cloned->ReflectZ(kFALSE); // printf("%s after\n", node->GetName()); // node->GetMatrix()->Print(); new_vol = MakeReflection(node->GetVolume()); @@ -281,18 +280,29 @@ namespace { /// Perform scan void ReflectionBuilder::execute() const { TGeoIterator next(detector.manager().GetTopVolume()); + bool print_active = isActivePrintLevel(DEBUG); TGeoNode *node; while ( (node=next()) ) { TGeoMatrix* m = node->GetMatrix(); if (m->IsReflection()) { + if ( print_active ) { + printout(INFO,"ReflectionBuilder","Reflection matrix:"); + m->Print(); + } Volume vol(node->GetVolume()); TGeoMatrix* mclone = new TGeoCombiTrans(*m); mclone->RegisterYourself(); // Reflect just the rotation component mclone->ReflectZ(kFALSE, kTRUE); + if ( print_active ) { + printout(INFO,"ReflectionBuilder","CLONE matrix:"); + mclone->Print(); + } TGeoNodeMatrix* nodematrix = (TGeoNodeMatrix*)node; nodematrix->SetMatrix(mclone); - printout(INFO,"Detector","Reflecting volume: %s ",vol.name()); + if ( print_active ) { + printout(INFO,"ReflectionBuilder","Reflected volume: %s ",vol.name()); + } Volume refl = vol.reflect(vol.sensitiveDetector()); node->SetVolume(refl.ptr()); } diff --git a/DDG4/python/DDG4.py b/DDG4/python/DDG4.py index 07eee9f27250b0217ac8e3bce58d4cde10380d8d..aa6276a2710fb06b7210b9fd8e80f7e012706f1f 100644 --- a/DDG4/python/DDG4.py +++ b/DDG4/python/DDG4.py @@ -561,6 +561,11 @@ class Geant4: return self def printDetectors(self): + """ + Scan the list of detectors and print detector name and sensitive type + + \author M.Frank + """ logger.info('+++ List of sensitive detectors:') for i in self.description.detectors(): o = DetElement(i.second.ptr()) # noqa: F405 @@ -572,7 +577,36 @@ class Geant4: sdtyp = self.sensitive_types[typ] logger.info('+++ %-32s type:%-12s --> Sensitive type: %s', o.name(), typ, sdtyp) + def setupDetectors(self): + """ + Scan the list of detectors and assign the proper sensitive actions + + \author M.Frank + """ + seq = None + actions = [] + logger.info('+++ Setting up sensitive detectors:') + for i in self.description.detectors(): + o = DetElement(i.second.ptr()) # noqa: F405 + sd = self.description.sensitiveDetector(str(o.name())) + if sd.isValid(): + typ = sd.type() + sdtyp = 'Unknown' + if typ in self.sensitive_types: + sdtyp = self.sensitive_types[typ] + seq, act = self.setupDetector(o.name(), sdtyp, collections=None) + logger.info('+++ %-32s type:%-12s --> Sensitive type: %s', o.name(), typ, sdtyp) + actions.append(act) + continue + logger.info('+++ %-32s --> UNKNOWN Sensitive type: %s', o.name(), typ) + return (seq, actions) + def setupDetector(self, name, action, collections=None): + """ + Setup single subdetector and assign the proper sensitive action + + \author M.Frank + """ # fg: allow the action to be a tuple with parameter dictionary sensitive_type = "" parameterDict = {} @@ -625,6 +659,11 @@ class Geant4: return (seq, acts[0]) def setupCalorimeter(self, name, type=None, collections=None): + """ + Setup subdetector of type 'calorimeter' and assign the proper sensitive action + + \author M.Frank + """ self.description.sensitiveDetector(str(name)) # sd.setType('calorimeter') if type is None: @@ -632,6 +671,11 @@ class Geant4: return self.setupDetector(name, type, collections) def setupTracker(self, name, type=None, collections=None): + """ + Setup subdetector of type 'tracker' and assign the proper sensitive action + + \author M.Frank + """ self.description.sensitiveDetector(str(name)) # sd.setType('tracker') if type is None: diff --git a/examples/ClientTests/compact/Check_Shape_Eightpoint_Reflect_DetElement.xml b/examples/ClientTests/compact/Check_Shape_Eightpoint_Reflect_DetElement.xml index c0d94c8435293378d13061f94871301cb9d025bc..f57e381206d4daa9e4a86c38bb2a050446aa400e 100644 --- a/examples/ClientTests/compact/Check_Shape_Eightpoint_Reflect_DetElement.xml +++ b/examples/ClientTests/compact/Check_Shape_Eightpoint_Reflect_DetElement.xml @@ -20,7 +20,8 @@ <rotation x="0" y="0" z="0"/> </check> <reflect vis="Shape2_vis"/> - <test type="DD4hep_Mesh_Verifier" ref="${DD4hepExamplesINSTALL}/examples/ClientTests/ref/Ref_Eightpoint_Reflect_DetElement.txt" create="CheckShape_create"/> + <test type="DD4hep_Mesh_Verifier" ref="${DD4hepExamplesINSTALL}/examples/ClientTests/ref/Ref_Eightpoint_Reflect_DetElement.txt" create="CheckShape_create"/> + <testn type="DD4hep_Mesh_Verifier" ref="${DD4hepExamplesINSTALL}/examples/ClientTests/ref/Ref_Eightpoint_Reflect_DetElement.txt" create="1"/> </detector> </detectors> </lccdd> diff --git a/examples/ClientTests/compact/NestedBoxReflection.xml b/examples/ClientTests/compact/NestedBoxReflection.xml index af560b628af42db338ae5468da69df7deb00aa49..f722a62bccca182ea760d34d6e89811bc5c2ba13 100644 --- a/examples/ClientTests/compact/NestedBoxReflection.xml +++ b/examples/ClientTests/compact/NestedBoxReflection.xml @@ -31,10 +31,9 @@ </limitset> </limits> - <geometry reflect="true"/> <detectors> - <detector id="1" name="NestedBox" type="NestedBoxReflection" readout="NestedBoxHits" vis="VisibleGreen" limits="cal_limits"> - <comment>A box with 3 boxes inside spanning a coordinate system</comment> + + <detector id="1" name="NestedBoxes" type="DD4hep_Test_NestedBoxReflection" readout="NestedBoxHits" vis="VisibleGreen" limits="cal_limits"> <assembly/> <dimensions x="50*cm" y="70*cm" z="90*cm" level="2"> <position x="400*cm" y="0*cm" z="0*cm"/> @@ -53,6 +52,19 @@ <reflect_z x="0*cm" y="0*cm" z="200*cm"/> <reflect_z x="0*cm" y="0*cm" z="-200*cm"/> + </detector> + + <!-- + --> + + <detector id="2" name="NestedTube" type="DD4hep_Test_TubeReflection" readout="NestedTubeHits" vis="VisibleGreen" limits="cal_limits"> + <assembly/> + <dimensions r="100*cm" dz="50*cm" level="2"> + <position x="0*cm" y="200*cm" z="200*cm"/> + </dimensions> + <reflect_z x="0*cm" y="200*cm" z="-200*cm"/> + </detector> + <!-- <reflect_x x="-200*cm" y="200*cm" z="0*cm"/> <reflect_x x="-200*cm" y="200*cm" z="200*cm"/> @@ -77,12 +89,15 @@ <reflect_z x="200*cm" y="-200*cm" z="-200*cm"/> <reflect_z x="200*cm" y="-200*cm" z="0*cm"/> --> - </detector> + </detectors> <readouts> <readout name="NestedBoxHits"> - <id>system:8,lvl4:5,lvl3:5,lvl2:5,lvl1:5,lvl0:5</id> + <id>system:8,lvl4:4,lvl3:4,lvl2:4,lvl1:4,lvl0:4</id> + </readout> + <readout name="NestedTubeHits"> + <id>system:8,lvl4:4,lvl3:4,lvl2:4,lvl1:4,lvl0:4</id> </readout> </readouts> diff --git a/examples/ClientTests/ref/Ref_Eightpoint_Reflect_DetElement.txt b/examples/ClientTests/ref/Ref_Eightpoint_Reflect_DetElement.txt index 6b6d2a102725249669f1e59d5619ae95c24939ac..54d1dc479047a7f343b8485fb22de5a3be555564 100644 --- a/examples/ClientTests/ref/Ref_Eightpoint_Reflect_DetElement.txt +++ b/examples/ClientTests/ref/Ref_Eightpoint_Reflect_DetElement.txt @@ -10,15 +10,15 @@ TGeoArb8 6 Local ( -23.00, 27.00, 30.00) Global ( -23.00, 27.00 TGeoArb8 7 Local ( 13.00, -27.00, 30.00) Global ( 13.00, -27.00, 130.00) TGeoArb8 Bounding box: dx= 27.50 dy= 27.00 dz= 30.00 Origin: x= -2.50 y= 0.00 z= 0.00 -ShapeCheck[1] TGeoArb8 8 Mesh-points: -TGeoArb8 EightPointSolid N(mesh)=8 N(vert)=8 N(seg)=12 N(pols)=6 -TGeoArb8 0 Local ( -30.00, -25.00, -30.00) Global ( -30.00, -25.00, -130.00) -TGeoArb8 1 Local ( -25.00, 25.00, -30.00) Global ( -25.00, 25.00, -130.00) -TGeoArb8 2 Local ( 5.00, 25.00, -30.00) Global ( 5.00, 25.00, -130.00) -TGeoArb8 3 Local ( 25.00, -25.00, -30.00) Global ( 25.00, -25.00, -130.00) -TGeoArb8 4 Local ( -28.00, -23.00, 30.00) Global ( -28.00, -23.00, -70.00) -TGeoArb8 5 Local ( -23.00, 27.00, 30.00) Global ( -23.00, 27.00, -70.00) -TGeoArb8 6 Local ( -23.00, 27.00, 30.00) Global ( -23.00, 27.00, -70.00) -TGeoArb8 7 Local ( 13.00, -27.00, 30.00) Global ( 13.00, -27.00, -70.00) -TGeoArb8 Bounding box: dx= 27.50 dy= 27.00 dz= 30.00 Origin: x= -2.50 y= 0.00 z= 0.00 +ShapeCheck[1] TGeoScaledShape 8 Mesh-points: +TGeoScaledShape Shape_Arbitrary_Eightpoint_vol_0_shape_refl N(mesh)=8 N(vert)=8 N(seg)=12 N(pols)=6 +TGeoScaledShape 0 Local ( -30.00, -25.00, 30.00) Global ( -30.00, -25.00, -70.00) +TGeoScaledShape 1 Local ( -25.00, 25.00, 30.00) Global ( -25.00, 25.00, -70.00) +TGeoScaledShape 2 Local ( 5.00, 25.00, 30.00) Global ( 5.00, 25.00, -70.00) +TGeoScaledShape 3 Local ( 25.00, -25.00, 30.00) Global ( 25.00, -25.00, -70.00) +TGeoScaledShape 4 Local ( -28.00, -23.00, -30.00) Global ( -28.00, -23.00, -130.00) +TGeoScaledShape 5 Local ( -23.00, 27.00, -30.00) Global ( -23.00, 27.00, -130.00) +TGeoScaledShape 6 Local ( -23.00, 27.00, -30.00) Global ( -23.00, 27.00, -130.00) +TGeoScaledShape 7 Local ( 13.00, -27.00, -30.00) Global ( 13.00, -27.00, -130.00) +TGeoScaledShape Bounding box: dx= 27.50 dy= 27.00 dz= 30.00 Origin: x= -2.50 y= 0.00 z= -0.00 diff --git a/examples/ClientTests/scripts/NestedBoxReflection.py b/examples/ClientTests/scripts/NestedBoxReflection.py new file mode 100644 index 0000000000000000000000000000000000000000..25281731826c90d2db7c3dd4dd0f4233b0514441 --- /dev/null +++ b/examples/ClientTests/scripts/NestedBoxReflection.py @@ -0,0 +1,182 @@ +# +# +from __future__ import absolute_import, unicode_literals +import os +import sys +import math +import time +import logging +import DDG4 +from DDG4 import OutputLevel as Output +from g4units import rad, keV, MeV, GeV, TeV, m, cm, mm, ns +# +# +logging.basicConfig(format='%(levelname)s: %(message)s', level=logging.INFO) +logger = logging.getLogger(__name__) +""" + + dd4hep simulation example setup using the python configuration + + @author M.Frank + @version 1.0 + +""" + + +def help(): + logging.info("Check_shape.py -option [-option] ") + logging.info(" -vis Enable visualization ") + logging.info(" -batch Batch execution ") + + +def run(): + hlp = False + vis = False + dump = False + batch = False + install_dir = os.environ['DD4hepINSTALL'] + # + for i in list(range(len(sys.argv))): + c = sys.argv[i].upper() + if c.find('BATCH') < 2 and c.find('BATCH') >= 0: + batch = True + elif c[:4] == '-VIS': + vis = True + elif c[:4] == '-DUM': + dump = True + elif c[:2] == '-H': + hlp = True + + if hlp: + help() + sys.exit(1) + + kernel = DDG4.Kernel() + description = kernel.detectorDescription() + install_dir = os.environ['DD4hepExamplesINSTALL'] + geant4 = DDG4.Geant4(kernel, tracker='Geant4TrackerCombineAction') + # + logger.info("# Configure UI") + ui = None + if batch: + geant4.setupCshUI(ui=None, vis=None) + kernel.UI = 'UI' + else: + ui = geant4.setupCshUI(vis=vis) + + kernel.loadGeometry(str("file:" + install_dir + "/examples/ClientTests/compact/NestedBoxReflection.xml")) + DDG4.importConstants(description) + + geant4.printDetectors() + if dump: + seq, act = geant4.addDetectorConstruction("Geant4DetectorGeometryConstruction/ConstructGeo") + act.DebugReflections = True + act.DebugMaterials = False + act.DebugElements = False + act.DebugVolumes = False + act.DebugShapes = False + act.DumpHierarchy = ~0x0 + + logger.info("# Configure G4 magnetic field tracking") + geant4.setupTrackingField() + + logger.info("# Setup random generator") + rndm = DDG4.Action(kernel, 'Geant4Random/Random') + rndm.Seed = 987654321 + rndm.initialize() + # + logger.info("# Configure Event actions") + prt = DDG4.EventAction(kernel, 'Geant4ParticlePrint/ParticlePrint') + prt.OutputType = 3 # Print both: table and tree + prt.OutputLevel = Output.INFO + kernel.eventAction().adopt(prt) + # + logger.info("# Configure I/O") + geant4.setupROOTOutput('RootOutput', 'BoxReflect_' + time.strftime('%Y-%m-%d_%H-%M')) + # + gen = DDG4.GeneratorAction(kernel, "Geant4GeneratorActionInit/GenerationInit") + gen.enableUI() + kernel.generatorAction().adopt(gen) + # + logger.info("# Generation of isotrope tracks of a given multiplicity with overlay:") + gen = DDG4.GeneratorAction(kernel, "Geant4ParticleGun/IsotropE+") + gen.mask = 4 + gen.isotrop = True + gen.particle = 'e+' + gen.energy = 100 * GeV + gen.multiplicity = 200 + gen.position = (0*m, 0*m, 0*m) + gen.direction = (0, 0, 1.) + gen.distribution = 'uniform' + gen.standalone = False + #gen.PhiMin = 0.0*rad + #gen.PhiMax = 2.0*math.pi*rad + #gen.ThetaMin = 0.0*math.pi*rad + #gen.ThetaMax = 1.0*math.pi*rad + gen.enableUI() + kernel.generatorAction().adopt(gen) + # + logger.info("# Merge all existing interaction records") + gen = DDG4.GeneratorAction(kernel, "Geant4InteractionMerger/InteractionMerger") + gen.OutputLevel = 4 # generator_output_level + gen.enableUI() + kernel.generatorAction().adopt(gen) + # + logger.info("# Finally generate Geant4 primaries") + gen = DDG4.GeneratorAction(kernel, "Geant4PrimaryHandler/PrimaryHandler") + gen.OutputLevel = 4 # generator_output_level + gen.enableUI() + kernel.generatorAction().adopt(gen) + # + logger.info("# ....and handle the simulation particles.") + part = DDG4.GeneratorAction(kernel, "Geant4ParticleHandler/ParticleHandler") + kernel.generatorAction().adopt(part) + # part.SaveProcesses = ['conv','Decay'] + part.SaveProcesses = ['Decay'] + part.MinimalKineticEnergy = 100 * MeV + part.OutputLevel = 5 # generator_output_level + part.enableUI() + user = DDG4.Action(kernel, "Geant4TCUserParticleHandler/UserParticleHandler") + user.TrackingVolume_Zmax = 3.0 * m + user.TrackingVolume_Rmax = 3.0 * m + user.enableUI() + part.adopt(user) + # + logger.info("# Now setup the calorimeters") + seq, actions = geant4.setupDetectors() + # + logger.info("# Now build the physics list:") + phys = geant4.setupPhysics('QGSP_BERT') + ph = geant4.addPhysics(str('Geant4PhysicsList/Myphysics')) + ph.addPhysicsConstructor(str('G4StepLimiterPhysics')) + # + # Add special particle types from specialized physics constructor + part = geant4.addPhysics('Geant4ExtraParticles/ExtraParticles') + part.pdgfile = os.path.join(install_dir, 'examples/DDG4/examples/particle.tbl') + # + # Add global range cut + rg = geant4.addPhysics('Geant4DefaultRangeCut/GlobalRangeCut') + rg.RangeCut = 0.7 * mm + # + # + if ui and vis: + cmds = [] + cmds.append('/control/verbose 2') + cmds.append('/run/initialize') + cmds.append('/vis/open OGL') + cmds.append('/vis/verbose errors') + cmds.append('/vis/drawVolume') + cmds.append('/vis/viewer/set/viewpointThetaPhi 55. 45.') + cmds.append('/vis/scene/add/axes 0 0 0 3 m') + ui.Commands = cmds + + kernel.configure() + kernel.initialize() + + # DDG4.setPrintLevel(Output.DEBUG) + kernel.run() + kernel.terminate() + + +if __name__ == "__main__": + run() diff --git a/examples/ClientTests/src/NestedBoxReflection_geo.cpp b/examples/ClientTests/src/NestedBoxReflection_geo.cpp index c8cf2264dcf35803e8e6bf874f3c240845f2db1e..6964950a9fe885feadc922c5282c600be12afafc 100644 --- a/examples/ClientTests/src/NestedBoxReflection_geo.cpp +++ b/examples/ClientTests/src/NestedBoxReflection_geo.cpp @@ -23,148 +23,140 @@ using namespace std; using namespace dd4hep; using namespace dd4hep::detail; +namespace { + const char* col(int idx) { + const char* cols[5] = {"VisibleRed","VisibleBlue","VisibleGreen","VisibleYellow","VisibleCyan"}; + return cols[idx%(sizeof(cols)/sizeof(cols[0]))]; + } + Rotation3D makeRotReflect(double thetaX, double phiX, double thetaY, double phiY, double thetaZ, double phiZ) { + // define 3 unit std::vectors forming the new left-handed axes + Position x(cos(phiX) * sin(thetaX), sin(phiX) * sin(thetaX), cos(thetaX)); + Position y(cos(phiY) * sin(thetaY), sin(phiY) * sin(thetaY), cos(thetaY)); + Position z(cos(phiZ) * sin(thetaZ), sin(phiZ) * sin(thetaZ), cos(thetaZ)); + + constexpr double tol = 1.0e-3; // Geant4 compatible + double check = (x.Cross(y)).Dot(z); // in case of a LEFT-handed orthogonal system this must be -1 + if (abs(1. + check) > tol) { + except("NestedBoxReflection", "+++ FAILED to construct Rotation is not LEFT-handed!"); + } + printout(INFO, "NestedBoxReflection", "+++ Constructed LEFT-handed reflection rotation."); + Rotation3D rotation(x.x(), y.x(), z.x(), x.y(), y.y(), z.y(), x.z(), y.z(), z.z()); + return rotation; + } + + TGeoCombiTrans* createPlacement(const Rotation3D& iRot, const Position& iTrans) { + TGeoRotation r; + double elements[9]; + iRot.GetComponents(elements); + r.SetMatrix(elements); + return new TGeoCombiTrans(TGeoTranslation(iTrans.x(), iTrans.y(), iTrans.z()), r); + } + + //static Transform3D transform_reflect(xml_h element) { + TGeoCombiTrans* transform_reflect(xml_h element) { + xml_dim_t xrot(element.child(_U(rotation))); + xml_dim_t xpos(element.child(_U(position))); + double thetaX = xrot.attr<double>(Unicode("thetaX")); + double phiX = xrot.attr<double>(Unicode("phiX")); + double thetaY = xrot.attr<double>(Unicode("thetaY")); + double phiY = xrot.attr<double>(Unicode("phiY")); + double thetaZ = xrot.attr<double>(Unicode("thetaZ")); + double phiZ = xrot.attr<double>(Unicode("phiZ")); + printout(INFO, "NestedBoxReflection", + "+++ Adding reflection rotation \"%s\": " + "(theta/phi)[rad] X: %6.3f %6.3f Y: %6.3f %6.3f Z: %6.3f %6.3f", + element.attr<string>(_U(name)).c_str(), thetaX, phiX, thetaY, phiY, thetaZ, phiZ); + Rotation3D rot = makeRotReflect(thetaX, phiX, thetaY, phiY, thetaZ, phiZ); + Position pos = Position(xpos.x(),xpos.y(),xpos.z()); + // return Transform3D(rot, pos); + return createPlacement(rot, pos); + } +} -// makes sure that the RotationMatrix built is -// LEFT-handed coordinate system (i.e. reflected) namespace { - struct NestedBoxReflection : public xml::tools::VolumeBuilder { + struct TubeReflection : public xml::tools::VolumeBuilder { + Material silicon; + Material iron; + /// using VolumeBuilder::VolumeBuilder; - - Rotation3D makeRotReflect(double thetaX, double phiX, double thetaY, double phiY, double thetaZ, double phiZ) { - // define 3 unit std::vectors forming the new left-handed axes - Position x(cos(phiX) * sin(thetaX), sin(phiX) * sin(thetaX), cos(thetaX)); - Position y(cos(phiY) * sin(thetaY), sin(phiY) * sin(thetaY), cos(thetaY)); - Position z(cos(phiZ) * sin(thetaZ), sin(phiZ) * sin(thetaZ), cos(thetaZ)); - - constexpr double tol = 1.0e-3; // Geant4 compatible - double check = (x.Cross(y)).Dot(z); // in case of a LEFT-handed orthogonal system this must be -1 - if (abs(1. + check) > tol) { - except("NestedBoxReflection", "+++ FAILED to construct Rotation is not LEFT-handed!"); - } - printout(INFO, "NestedBoxReflection", "+++ Constructed LEFT-handed reflection rotation."); - Rotation3D rotation(x.x(), y.x(), z.x(), x.y(), y.y(), z.y(), x.z(), y.z(), z.z()); - return rotation; - } + /// + void place_quadrants(int /* level */, Volume vol) { + Tube tube = vol.solid(); + double rmin = tube.rMin(), rmax = tube.rMax(), dz = tube.dZ(); + Tube tsol1(rmin, rmax, dz*0.05, 0, M_PI*2.); + Tube tsol2(rmin, rmax, dz*0.15, 0, M_PI/2.); + Volume tvol0("Tube0", tsol1, silicon); + Volume tvol1("Tube1", tsol2, silicon); + Volume tvol2("Tube2", tsol2, silicon); + Volume tvol3("Tube3", tsol2, silicon); + Volume tvol4("Tube4", tsol2, silicon); + TGeoHMatrix* mat; + Transform3D tr; + + double dz1 = dz-2.0*tsol2.dZ()-tsol1.dZ(); + tvol0.setVisAttributes(description.visAttributes("VisibleCyan")); + tvol1.setVisAttributes(description.visAttributes("VisibleRed")); + tvol2.setVisAttributes(description.visAttributes("VisibleBlue")); + tvol3.setVisAttributes(description.visAttributes("VisibleYellow")); + tvol4.setVisAttributes(description.visAttributes("VisibleGreen")); - TGeoCombiTrans* createPlacement(const Rotation3D& iRot, const Position& iTrans) { - TGeoRotation r; - double elements[9]; - iRot.GetComponents(elements); - r.SetMatrix(elements); - return new TGeoCombiTrans(TGeoTranslation(iTrans.x(), iTrans.y(), iTrans.z()), r); - } + mat = detail::matrix::_transform(Transform3D(Position(0, 0, dz1))); + vol.placeVolume(tvol0, mat); - //static Transform3D transform_reflect(xml_h element) { - TGeoCombiTrans* transform_reflect(xml_h element) { - xml_dim_t xrot(element.child(_U(rotation))); - xml_dim_t xpos(element.child(_U(position))); - double thetaX = xrot.attr<double>(Unicode("thetaX")); - double phiX = xrot.attr<double>(Unicode("phiX")); - double thetaY = xrot.attr<double>(Unicode("thetaY")); - double phiY = xrot.attr<double>(Unicode("phiY")); - double thetaZ = xrot.attr<double>(Unicode("thetaZ")); - double phiZ = xrot.attr<double>(Unicode("phiZ")); - printout(INFO, "NestedBoxReflection", - "+++ Adding reflection rotation \"%s\": " - "(theta/phi)[rad] X: %6.3f %6.3f Y: %6.3f %6.3f Z: %6.3f %6.3f", - element.attr<string>(_U(name)).c_str(), thetaX, phiX, thetaY, phiY, thetaZ, phiZ); - Rotation3D rot = makeRotReflect(thetaX, phiX, thetaY, phiY, thetaZ, phiZ); - Position pos = Position(xpos.x(),xpos.y(),xpos.z()); - // return Transform3D(rot, pos); - return createPlacement(rot, pos); - } + mat = detail::matrix::_transform(Transform3D(Position(0, 0,-dz1))); + mat->ReflectZ(kTRUE, kTRUE); + vol.placeVolume(tvol0, mat); - void place_boxes(int level, Volume vol) { - if ( level >= 0 ) { - Box box = vol.solid(); - double line = 0.015; - double bx = box.x(); - double by = box.y(); - double bz = box.z(); - Material mat = description.material("Si"); - bool sens = level == 2 || level == 0; - Box long_box(bx*0.2, by*0.2, bz*0.4); - Box flat_box(bx*1.0, by*1.0, bz*0.1); - Box mini_box(bx*0.2, by*0.2, bz*0.2); - const char* cols[5] = {"VisibleRed","VisibleBlue","VisibleGreen","VisibleYellow","VisibleCyan"}; - const char* c; - PlacedVolume pv; - Volume v; - - c = cols[(0+level)%5]; - v = Volume(_toString(1,"box%d"), mini_box, mat); - v.setRegion(vol.region()); - v.setLimitSet(vol.limitSet()); - if ( !sens ) v.setSensitiveDetector(sensitive); - v.setVisAttributes(description, c); - pv = vol.placeVolume(v, Position(0,0,0)); - pv.addPhysVolID(_toString(level,"lvl%d"), 1); - printout(INFO,"NestedBoxReflection","++ Volume: %s Color: %s", v.name(), c); - place_boxes(level-1, v); - - c = cols[(1+level)%5]; - v = Volume(_toString(2,"box%d"), long_box, mat); - v.setRegion(vol.region()); - v.setLimitSet(vol.limitSet()); - if ( !sens ) v.setSensitiveDetector(sensitive); - v.setVisAttributes(description, c); - pv = vol.placeVolume(v, Position(0.8*bx,0,0)); - pv.addPhysVolID(_toString(level,"lvl%d"), 2); - printout(INFO,"NestedBoxReflection","++ Volume: %s Color: %s", v.name(), c); - place_boxes(level-1, v); - - v = Volume(_toString(1,"axis_x"), Box(bx*0.2, by*line, bz*line), mat); - v.setVisAttributes(description, c); - pv = vol.placeVolume(v, Position(0.4*bx,0,0)); - printout(INFO,"NestedBoxReflection","++ Volume: %s Color: %s", v.name(), c); - - c = cols[(2+level)%5]; - v = Volume(_toString(3,"box%d"), mini_box, mat); - v.setRegion(vol.region()); - v.setLimitSet(vol.limitSet()); - if ( !sens ) v.setSensitiveDetector(sensitive); - v.setVisAttributes(description, c); - pv = vol.placeVolume(v, Position(0,0.8*by,0)); - pv.addPhysVolID(_toString(level,"lvl%d"), 3); - printout(INFO,"NestedBoxReflection","++ Volume: %s Color: %s", v.name(), c); - place_boxes(level-1, v); - - v = Volume(_toString(1,"axis_y"), Box(bx*line, by*0.2, bz*line), mat); - v.setVisAttributes(description, c); - pv = vol.placeVolume(v, Position(0,0.4*by,0)); - printout(INFO,"NestedBoxReflection","++ Volume: %s Color: %s", v.name(), c); - - c = cols[(3+level)%5]; - v = Volume(_toString(4,"box%d"), mini_box, mat); - v.setRegion(vol.region()); - v.setLimitSet(vol.limitSet()); - if ( !sens ) v.setSensitiveDetector(sensitive); - v.setVisAttributes(description, c); - pv = vol.placeVolume(v, Position(0,0,0.8*bz)); - pv.addPhysVolID(_toString(level,"lvl%d"), 4); - printout(INFO,"NestedBoxReflection","++ Volume: %s Color: %s", v.name(), c); - place_boxes(level-1, v); - - v = Volume(_toString(1,"axis_z"), Box(bx*line, by*line, bz*0.2), mat); - v.setVisAttributes(description, c); - pv = vol.placeVolume(v, Position(0,0,0.4*bz)); - printout(INFO,"NestedBoxReflection","++ Volume: %s Color: %s", v.name(), c); - - c = cols[(4+level)%5]; - v = Volume(_toString(5,"box%d"), flat_box, mat); - v.setRegion(vol.region()); - v.setLimitSet(vol.limitSet()); - v.setVisAttributes(description, c); - if ( !sens ) v.setSensitiveDetector(sensitive); - TGeoHMatrix* mtx = detail::matrix::_transform(Position(0,0,0.9*bz)); - mtx->ReflectZ(kTRUE, kFALSE); - pv = vol.placeVolume(v, mtx); - pv.addPhysVolID(_toString(level,"lvl%d"), 5); - printout(INFO,"NestedBoxReflection","++ Volume: %s Color: %s", v.name(), c); - place_boxes(level-1, v); - } +#if 0 +#endif + double dz2 = dz-tsol2.dZ(); + // OK + mat = detail::matrix::_transform(Transform3D(Position(0, 0, dz2))); + vol.placeVolume(tvol1, mat); + // OK + mat = detail::matrix::_transform(Transform3D(Position(0, 0, dz2))); + mat->RotateZ(1.0*M_PI/2.0/dd4hep::degree); + vol.placeVolume(tvol2, mat); + // OK + mat = detail::matrix::_transform(Transform3D(Position(0, 0, dz2))); + mat->RotateZ(2.0*M_PI/2.0/dd4hep::degree); + vol.placeVolume(tvol3, mat); + // OK + mat = detail::matrix::_transform(Transform3D(Position(0, 0, dz2))); + mat->RotateZ(3.0*M_PI/2.0/dd4hep::degree); + vol.placeVolume(tvol4, mat); + + + /** Now eflect the quadrants to the other endcap: */ + + // OK + mat = detail::matrix::_transform(Transform3D(Position(0, 0,-dz2))); + mat->ReflectY(kTRUE,kTRUE); + //mat->RotateZ(1.0*M_PI/2.0/dd4hep::degree); + vol.placeVolume(tvol1, mat); + + // OK + mat = detail::matrix::_transform(Transform3D(Position(0, 0,-dz2))); + mat->RotateZ(1.0*M_PI/2.0/dd4hep::degree); + mat->ReflectY(kTRUE,kTRUE); + //mat->RotateZ(-1.0*M_PI/2.0/dd4hep::degree); + vol.placeVolume(tvol2, mat); + + // OK + mat = detail::matrix::_transform(Transform3D(Position(0, 0,-dz2))); + mat->RotateZ(2.0*M_PI/2.0/dd4hep::degree); + mat->ReflectY(kTRUE,kTRUE); + //mat->RotateZ(1.0*M_PI/2.0/dd4hep::degree); + vol.placeVolume(tvol3, mat); + + // OK + mat = detail::matrix::_transform(Transform3D(Position(0, 0,-dz2))); + mat->RotateZ(3.0*M_PI/2.0/dd4hep::degree); + mat->ReflectY(kTRUE,kTRUE); + //mat->RotateZ(-1.0*M_PI/2.0/dd4hep::degree); + vol.placeVolume(tvol4, mat); } - + /// PlacedVolume place(Volume mother, Volume vol, xml_elt_t e, int level, int copyNo, char reflection) { Position pos; RotationZYX rot; @@ -191,53 +183,45 @@ namespace { pv.addPhysVolID(_toString(level,"lvl%d"), copyNo); return pv; } - - Ref_t create() { - xml_dim_t x_box(x_det.dimensions()); - int levels = x_box.level(2); - double bx = x_box.x(); - double by = x_box.y(); - double bz = x_box.z(); + /// + DetElement create() { + xml_dim_t x_tube(x_det.dimensions()); + double r = x_tube.r(); + double dz = x_tube.dz(); + int levels = x_tube.level(2); + Volume tube_vol; Volume v_det; - Box box_solid(bx,by,bz); - Volume box_vol("nested_box",box_solid,description.air()); PlacedVolume pv; sensitive.setType("tracker"); + silicon = description.material("Si"); + iron = description.material("Iron"); + tube_vol = Volume("tube", Tube(0,r,dz,0,M_PI*2.0), silicon); if ( levels != 0 && x_det.hasChild(_U(assembly)) ) - v_det = Assembly("envelope"); + v_det = Assembly("envelope"); else - v_det = Volume("envelope",Box(4.5*bx,4.5*by,4.5*bz),description.air()); + v_det = Volume("envelope",Tube(4.5*r,4.5*r,4.5*r),description.air()); - if ( levels == 0 ) { - v_det.setSensitiveDetector(sensitive); - } - else if ( levels == 1 ) { - box_vol.setSensitiveDetector(sensitive); - box_vol.setAttributes(description,x_det.regionStr(),x_det.limitsStr(),"VisibleGrey"); - } - else { - int cnt = 1; - Transform3D tr = xml::createTransformation(x_box); - box_vol.setAttributes(description,x_det.regionStr(),x_det.limitsStr(),"VisibleGrey"); - place_boxes(levels-1, box_vol); - pv = v_det.placeVolume(box_vol, tr); - pv.addPhysVolID(_toString(levels,"lvl%d"), ++cnt); - - for(xml_coll_t c(x_det,_U(reflect_x)); c; ++c) - place(v_det, box_vol, c, levels, ++cnt, 'X'); - for(xml_coll_t c(x_det,_U(reflect_y)); c; ++c) - place(v_det, box_vol, c, levels, ++cnt, 'Y'); - for(xml_coll_t c(x_det,_U(reflect_z)); c; ++c) - place(v_det, box_vol, c, levels, ++cnt, 'Z'); - - if ( x_det.hasChild(_U(reflect)) ) { - Volume reflect_vol = box_vol; - for(xml_coll_t c(x_det,_U(reflect)); c; ++c) { - TGeoCombiTrans* reflect_tr = transform_reflect(c); - pv = v_det.placeVolume(reflect_vol.ptr(), reflect_tr); - pv.addPhysVolID(_toString(levels,"lvl%d"), ++cnt); - } + int cnt = 0; + Transform3D tr = xml::createTransformation(x_tube); + tube_vol.setAttributes(description, x_det.regionStr(), x_det.limitsStr(), "VisibleGrey"); + place_quadrants(levels-1, tube_vol); + pv = v_det.placeVolume(tube_vol, tr); + pv.addPhysVolID(_toString(levels,"lvl%d"), ++cnt); + + for(xml_coll_t c(x_det,_U(reflect_x)); c; ++c) + place(v_det, tube_vol, c, levels, ++cnt, 'X'); + for(xml_coll_t c(x_det,_U(reflect_y)); c; ++c) + place(v_det, tube_vol, c, levels, ++cnt, 'Y'); + for(xml_coll_t c(x_det,_U(reflect_z)); c; ++c) + place(v_det, tube_vol, c, levels, ++cnt, 'Z'); + + if ( x_det.hasChild(_U(reflect)) ) { + Volume reflect_vol = tube_vol; + for(xml_coll_t c(x_det,_U(reflect)); c; ++c) { + TGeoCombiTrans* reflect_tr = transform_reflect(c); + pv = v_det.placeVolume(reflect_vol.ptr(), reflect_tr); + pv.addPhysVolID(_toString(levels,"lvl%d"), ++cnt); } } // Place the calo inside the world @@ -245,11 +229,167 @@ namespace { return detector; } }; + + static Ref_t create_refl_tube(Detector& description, xml_dim_t x_det, SensitiveDetector sens) { + TubeReflection builder(description, x_det, sens); + return builder.create(); + } +} +DECLARE_DETELEMENT(DD4hep_Test_TubeReflection,create_refl_tube) + + +// makes sure that the RotationMatrix built is +// LEFT-handed coordinate system (i.e. reflected) +namespace { + struct NestedBoxReflection : public xml::tools::VolumeBuilder { + using VolumeBuilder::VolumeBuilder; + Material silicon; + Material iron; + + void place_box(int level, int ident, const char* vis, Volume par, Solid sol, TGeoHMatrix* mtx); + void place_box(int level, int ident, const char* vis, Volume par, Solid sol, const Transform3D& tr) { + TGeoHMatrix* mtx = detail::matrix::_transform(tr); + place_box(level, ident, vis, par, sol, mtx); + } + void place_box(int level, int ident, const char* vis, Volume par, Solid sol, Position pos) { + place_box(level, ident, vis, par, sol, Transform3D(pos)); + } + void place_boxes(int level, Volume vol); + PlacedVolume place(Volume mother, Volume vol, xml_elt_t e, int level, int copyNo, char reflection); + DetElement create(); + }; + /// + void NestedBoxReflection::place_boxes(int level, Volume vol) { + if ( level >= 0 ) { + Box box = vol.solid(); + double line = 0.015, bx = box.x(), by = box.y(), bz = box.z(); + Box mbox(bx*0.2, by*0.2, bz*0.2); + + printout(INFO,"NestedBoxReflection","+++ Placing boxes at level %d",level); + place_box(level, 1, col(0+level), vol, mbox, Position(0,0,0)); + place_box(level, 2, col(1+level), vol, Box(bx*0.2, by*0.2, bz*0.4), Position(0.8*bx,0,0)); + place_box(level, 3, col(2+level), vol, mbox, Position(0,0.8*by,0)); + place_box(level, 4, col(3+level), vol, mbox, Position(0,0,0.8*bz)); + auto mtx = make_unique<TGeoHMatrix>(TGeoTranslation(0,0,-0.9*bz)); + mtx->ReflectZ(kTRUE,kTRUE); + place_box(level, 5, col(4+level), vol, Box(bx*1.0, by*1.0, bz*0.1), mtx.release()); + place_box(level, 0, col(1+level), vol, Box(bx*0.2, by*line, bz*line), Position(0.4*bx,0,0)); + place_box(level, 0, col(2+level), vol, Box(bx*line, by*0.2, bz*line), Position(0,0.4*by,0)); + place_box(level, 0, col(3+level), vol, Box(bx*line, by*line, bz*0.2), Position(0,0,0.4*bz)); + } + } + + /// + void NestedBoxReflection::place_box(int level, int ident, const char* vis, Volume par, Solid sol, TGeoHMatrix* tr) { + PlacedVolume pv; + bool sens = level == 2 || level == 0; + string idnam = _toString(ident,"box%d"); + Volume vol = Volume(idnam, sol, sens ? silicon : iron); + vol.setRegion(par.region()); + vol.setLimitSet(par.limitSet()); + vol.setVisAttributes(description, vis); + pv = par.placeVolume(vol, tr); + if ( ident > 0 ) { + string vidn = _toString(level,"lvl%d"); + if ( sens ) { + vol.setSensitiveDetector(sensitive); + } + pv.addPhysVolID(vidn, ident); + string n = par.name()+string("/")+vol.name(); + printout(INFO,"NestedBoxReflection","++ Level: %3d Volume:%-24s Sensitive:%s Color:%s vid:%s=%d", + level, n.c_str(), yes_no(sens), vis, vidn.c_str(), ident); + place_boxes(level-1, vol); + } + } + + /// + PlacedVolume NestedBoxReflection::place(Volume mother, Volume vol, xml_elt_t e, int level, int copyNo, char reflection) { + Position pos; + RotationZYX rot; + xml_dim_t xpos = xml_dim_t(e).child(_U(position), false); + xml_dim_t xrot = xml_dim_t(e).child(_U(rotation), false); + if ( !xpos.ptr() ) xpos = e; + if ( xpos.ptr() ) pos = Position(xpos.x(0),xpos.y(0),xpos.z(0)); + if ( xrot.ptr() ) rot = RotationZYX(xrot.x(0),xrot.y(0),xrot.z(0)); + + TGeoHMatrix* mat = detail::matrix::_transform(pos, rot); + switch(reflection) { + case 'X': + mat->ReflectX(kTRUE, kTRUE); + break; + case 'Y': + mat->ReflectY(kTRUE, kTRUE); + break; + case 'Z': + default: + mat->ReflectZ(kTRUE, kTRUE); + break; + } + PlacedVolume pv = mother.placeVolume(vol, mat); + pv.addPhysVolID(_toString(level,"lvl%d"), copyNo); + return pv; + } + + /// + DetElement NestedBoxReflection::create() { + xml_dim_t x_box(x_det.dimensions()); + int levels = x_box.level(2); + double bx = x_box.x(); + double by = x_box.y(); + double bz = x_box.z(); + Volume v_det; + Box box_solid(bx,by,bz); + Volume box_vol("nested_box",box_solid,description.air()); + PlacedVolume pv; + + sensitive.setType("tracker"); + iron = description.material("Iron"); + silicon = description.material("Si"); + if ( levels != 0 && x_det.hasChild(_U(assembly)) ) + v_det = Assembly("envelope"); + else + v_det = Volume("envelope",Box(4.5*bx,4.5*by,4.5*bz),description.air()); + + if ( levels == 0 ) { + v_det.setSensitiveDetector(sensitive); + } + else if ( levels == 1 ) { + box_vol.setSensitiveDetector(sensitive); + box_vol.setAttributes(description,x_det.regionStr(),x_det.limitsStr(),"VisibleGrey"); + } + else { + int cnt = 0; + Transform3D tr = xml::createTransformation(x_box); + box_vol.setAttributes(description,x_det.regionStr(),x_det.limitsStr(),"VisibleGrey"); + place_boxes(levels-1, box_vol); + pv = v_det.placeVolume(box_vol, tr); + pv.addPhysVolID(_toString(levels,"lvl%d"), ++cnt); + + for(xml_coll_t c(x_det,_U(reflect_x)); c; ++c) + place(v_det, box_vol, c, levels, ++cnt, 'X'); + for(xml_coll_t c(x_det,_U(reflect_y)); c; ++c) + place(v_det, box_vol, c, levels, ++cnt, 'Y'); + for(xml_coll_t c(x_det,_U(reflect_z)); c; ++c) + place(v_det, box_vol, c, levels, ++cnt, 'Z'); + + if ( x_det.hasChild(_U(reflect)) ) { + Volume reflect_vol = box_vol; + for(xml_coll_t c(x_det,_U(reflect)); c; ++c) { + TGeoCombiTrans* reflect_tr = transform_reflect(c); + pv = v_det.placeVolume(reflect_vol.ptr(), reflect_tr); + pv.addPhysVolID(_toString(levels,"lvl%d"), ++cnt); + } + } + } + // Place the calo inside the world + placeDetector(v_det, x_det).addPhysVolID("system",x_det.id()); + return detector; + } } -static Ref_t create_detector(Detector& description, xml_dim_t x_det, SensitiveDetector sens) { +static Ref_t create_nested_box(Detector& description, xml_dim_t x_det, SensitiveDetector sens) { NestedBoxReflection builder(description, x_det, sens); return builder.create(); } -DECLARE_DETELEMENT(NestedBoxReflection,create_detector) +DECLARE_DETELEMENT(DD4hep_Test_NestedBoxReflection,create_nested_box) diff --git a/examples/DDCMS/CMakeLists.txt b/examples/DDCMS/CMakeLists.txt index eda773331ec0ea7325ec01187ec409aac8c84a2d..a1411d462b2aee812a91a6b114981e30a02c8991 100644 --- a/examples/DDCMS/CMakeLists.txt +++ b/examples/DDCMS/CMakeLists.txt @@ -35,6 +35,12 @@ else() dd4hep_print("|++> CLHEP is not present. NOT building DDCMS examples.") return() endif() +#========================================================================== +if(NOT TARGET Geant4::Interface) + dd4hep_print("Not creating DDCMS tests [No Geant4 found]") + return() +endif() +#========================================================================== # dd4hep_configure_output () # @@ -49,9 +55,10 @@ dd4hep_configure_output () dd4hep_add_plugin ( DDCMS SOURCES src/*.cpp src/plugins/*.cpp USES DD4hep::DDCore DD4hep::DDAlign DD4hep::DDCond - ROOT::Core ROOT::Geom ROOT::GenVector CLHEP::CLHEP + ROOT::Core ROOT::Geom ROOT::GenVector CLHEP::CLHEP ${CLHEP} Geant4::Interface ) target_include_directories(DDCMS PUBLIC include) +target_link_options(DDCMS PRIVATE -L${Geant4_DIR}/.. -lG4clhep) # # install(TARGETS DDCMS LIBRARY DESTINATION lib) diff --git a/examples/DDCMS/src/plugins/DDEcalEndcapAlgo.cpp b/examples/DDCMS/src/plugins/DDEcalEndcapAlgo.cpp index 8ca14e83b6c5518f06c4f701d2d0e79ec14729c3..1437b27c50434f3866b82cf9e1d3fae2d4577037 100644 --- a/examples/DDCMS/src/plugins/DDEcalEndcapAlgo.cpp +++ b/examples/DDCMS/src/plugins/DDEcalEndcapAlgo.cpp @@ -5,8 +5,8 @@ #include <Math/Vector3D.h> #include <CLHEP/Geometry/Transform3D.h> -#include <CLHEP/Units/GlobalPhysicalConstants.h> -#include <CLHEP/Units/GlobalSystemOfUnits.h> +#include <CLHEP/Units/PhysicalConstants.h> +#include <CLHEP/Units/SystemOfUnits.h> #include <CLHEP/Geometry/Point3D.h> #include <CLHEP/Geometry/Plane3D.h> #include <CLHEP/Geometry/Vector3D.h>