diff --git a/DDG4/include/DDG4/Geant4Converter.h b/DDG4/include/DDG4/Geant4Converter.h index 2a883ba1623bafae782ee13408440f8efe4f9335..c86ed4d01ca42329ef3cc181e15c99879c743de0 100644 --- a/DDG4/include/DDG4/Geant4Converter.h +++ b/DDG4/include/DDG4/Geant4Converter.h @@ -33,7 +33,6 @@ class G4PVPlacement; class G4Region; class G4Field; class G4FieldManager; -class Geant4SensitiveDetector; class G4UserLimits; /* @@ -57,17 +56,24 @@ namespace DD4hep { * @version 1.0 */ struct Geant4Converter : public Geometry::GeoHandler { - typedef std::map<const TNamed*,G4UserLimits*> LimitMap; + typedef std::map<const TGeoElement*,G4Element*> ElementMap; + typedef std::map<const TGeoMedium*, G4Material*> MaterialMap; + typedef std::map<const TNamed*, G4UserLimits*> LimitMap; + typedef std::map<const TGeoNode*, G4PVPlacement*> PlacementMap; + typedef std::map<const TNamed*, G4Region*> RegionMap; + typedef std::map<const TNamed*, Geant4SensitiveDetector*> SensDetMap; + typedef std::map<const TGeoVolume*, G4LogicalVolume*> VolumeMap; + typedef std::map<const TGeoShape*, G4VSolid*> SolidMap; struct G4GeometryInfo : public GeometryInfo { - std::map<const TGeoElement*,G4Element*> g4Elements; - std::map<const TGeoMedium*, G4Material*> g4Materials; - std::map<const TGeoShape*, G4VSolid*> g4Solids; - std::map<const TGeoVolume*, G4LogicalVolume*> g4Volumes; - std::map<const TGeoNode*, G4PVPlacement*> g4Placements; - std::map<const TNamed*, G4Region*> g4Regions; - LimitMap g4Limits; - std::map<const TNamed*, Geant4SensitiveDetector*> g4SensDets; + ElementMap g4Elements; + MaterialMap g4Materials; + SolidMap g4Solids; + VolumeMap g4Volumes; + PlacementMap g4Placements; + RegionMap g4Regions; + LimitMap g4Limits; + SensDetMap g4SensDets; SensitiveVolumes sensitives; RegionVolumes regions; @@ -75,14 +81,21 @@ namespace DD4hep { }; LCDD& m_lcdd; + bool m_checkOverlaps; + G4GeometryInfo* m_dataPtr; G4GeometryInfo& data() const { return *m_dataPtr; } - /// Constructor - Geant4Converter( LCDD& lcdd ) : m_lcdd(lcdd) {} + /// Initializing Constructor + Geant4Converter( LCDD& lcdd ); + /// Standard destructor virtual ~Geant4Converter() {} - /// Create geometry dump + + /// Singleton instance + static Geant4Converter& instance(); + + /// Create geometry conversion void create(DetElement top); /// Convert the geometry type material into the corresponding Geant4 object(s). @@ -91,10 +104,14 @@ namespace DD4hep { virtual void* handleElement(const std::string& name, const TGeoElement* element) const; /// Convert the geometry type solid into the corresponding Geant4 object(s). virtual void* handleSolid(const std::string& name, const TGeoShape* volume) const; + /// Convert the geometry type logical volume into the corresponding Geant4 object(s). virtual void* handleVolume(const std::string& name, const TGeoVolume* volume) const; + virtual void* collectVolume(const std::string& name, const TGeoVolume* volume) const; + /// Convert the geometry type volume placement into the corresponding Geant4 object(s). virtual void* handlePlacement(const std::string& name, const TGeoNode* node) const; + /// Convert the geometry type field into the corresponding Geant4 object(s). ///virtual void* handleField(const std::string& name, Ref_t field) const; @@ -106,6 +123,13 @@ namespace DD4hep { virtual void* handleSensitive(const TNamed* sens_det, const std::set<const TGeoVolume*>& volumes) const; /// Handle the geant 4 specific properties void handleProperties(LCDD::Properties& prp) const; + + /// Print the geometry type SensitiveDetector + virtual void* printSensitive(const TNamed* sens_det, const std::set<const TGeoVolume*>& volumes) const; + /// Print Geant4 placement + virtual void* printPlacement(const std::string& name, const TGeoNode* node) const; + + }; } // End namespace Simulation } // End namespace DD4hep diff --git a/DDG4/include/DDG4/Geant4Hits.h b/DDG4/include/DDG4/Geant4Hits.h index 55ff90837688d5601d1b8c5392f955699dfbc329..71f90ee929a3b0e32dddc7a1b27f31e9ac517142 100644 --- a/DDG4/include/DDG4/Geant4Hits.h +++ b/DDG4/include/DDG4/Geant4Hits.h @@ -11,7 +11,7 @@ // Framework include files #include "DD4hep/Objects.h" -#include "DDG4/Defs.h" +#include "DDG4/Geant4StepHandler.h" // Geant4 include files #include "G4VHit.hh" @@ -32,7 +32,6 @@ namespace DD4hep { // Forward declarations; template<class HIT> struct HitCompare; template<class HIT> struct HitPositionCompare; - class Geant4StepHandler; class Geant4Hit; class Geant4TrackerHit; class Geant4CalorimeterHit; @@ -65,49 +64,6 @@ namespace DD4hep { virtual bool operator()(const HIT* h) const { return pos == h->position; } }; - /** @class Geant4StepHandler Geant4SensitiveDetector.h DDG4/Geant4SensitiveDetector.h - * - * Tiny helper/utility class to easily access Geant4 step information. - * Born by lazyness: Avoid typing millions of statements! - * - * @author M.Frank - * @version 1.0 - */ - class Geant4StepHandler { - public: - G4Step* step; - G4StepPoint* pre; - G4StepPoint* post; - G4Track* track; - Geant4StepHandler(G4Step* s) : step(s) { - pre = s->GetPreStepPoint(); - post = s->GetPostStepPoint(); - track = s->GetTrack(); - } - Position prePos() const { - const G4ThreeVector& p = pre->GetPosition(); - return Position(p.x(),p.y(),p.z()); - } - Position postPos() const { - const G4ThreeVector& p = post->GetPosition(); - return Position(p.x(),p.y(),p.z()); - } - Momentum preMom() const { - const G4ThreeVector& p = pre->GetMomentum(); - return Momentum(p.x(),p.y(),p.z()); - } - Momentum postMom() const { - const G4ThreeVector& p = post->GetMomentum(); - return Momentum(p.x(),p.y(),p.z()); - } - G4VPhysicalVolume* volume(G4StepPoint* p) const { - return p->GetTouchableHandle()->GetVolume(); - } - G4VSensitiveDetector* sd(G4StepPoint* p) const { - return p->GetPhysicalVolume()->GetLogicalVolume()->GetSensitiveDetector(); - } - }; - /** @class Geant4Hit Geant4Hits.h DDG4/Geant4Hits.h * * Geant4 hit base class. Here only the basic diff --git a/DDG4/include/DDG4/Geant4SensitiveDetector.h b/DDG4/include/DDG4/Geant4SensitiveDetector.h index 241f3f3a3b79d655a8c8095466b285d11c32694e..79ad52cc7fb4f2d1f2503ed89c0dbfcbb8dcbbee 100644 --- a/DDG4/include/DDG4/Geant4SensitiveDetector.h +++ b/DDG4/include/DDG4/Geant4SensitiveDetector.h @@ -65,6 +65,9 @@ namespace DD4hep { /// Find hits by position in a collection template <typename T> T* find(const HitCollection* c,const HitCompare<T>& cmp); + /// Dump Step information (careful: very verbose) + void dumpStep(G4Step* step,G4TouchableHistory* history); + public: /// Constructor. The sensitive detector element is identified by the detector name diff --git a/DDG4/include/DDG4/Geant4StepHandler.h b/DDG4/include/DDG4/Geant4StepHandler.h new file mode 100644 index 0000000000000000000000000000000000000000..9cfb9f853b1fc7937f5dcdcdf61909c5ce7ae6d0 --- /dev/null +++ b/DDG4/include/DDG4/Geant4StepHandler.h @@ -0,0 +1,114 @@ +// $Id:$ +//==================================================================== +// AIDA Detector description implementation +//-------------------------------------------------------------------- +// +// Author : M.Frank +// +//==================================================================== +#ifndef DD4HEP_GEANT4STEPHANDLER_H +#define DD4HEP_GEANT4STEPHANDLER_H + +// Framework include files +#include "DDG4/Defs.h" + +// Geant4 include files +#include "G4Step.hh" +#include "G4StepPoint.hh" +#include "G4VTouchable.hh" +#include "G4VSensitiveDetector.hh" + + +/* + * DD4hep namespace declaration + */ +namespace DD4hep { + + /* + * Simulation namespace declaration + */ + namespace Simulation { + + // Forward declarations; + class Geant4StepHandler; + + /** @class Geant4StepHandler Geant4SensitiveDetector.h DDG4/Geant4SensitiveDetector.h + * + * Tiny helper/utility class to easily access Geant4 step information. + * Born by lazyness: Avoid typing millions of statements! + * + * @author M.Frank + * @version 1.0 + */ + class Geant4StepHandler { + public: + const G4Step* step; + G4StepPoint* pre; + G4StepPoint* post; + G4Track* track; + Geant4StepHandler(const G4Step* s) : step(s) { + pre = s->GetPreStepPoint(); + post = s->GetPostStepPoint(); + track = s->GetTrack(); + } + G4ParticleDefinition* trackDef() const { + return track->GetDefinition(); + } + static const char* stepStatus(G4StepStatus status); + const char* preStepStatus() const; + const char* postStepStatus() const; + Position prePos() const { + const G4ThreeVector& p = pre->GetPosition(); + return Position(p.x(),p.y(),p.z()); + } + Position postPos() const { + const G4ThreeVector& p = post->GetPosition(); + return Position(p.x(),p.y(),p.z()); + } + Momentum preMom() const { + const G4ThreeVector& p = pre->GetMomentum(); + return Momentum(p.x(),p.y(),p.z()); + } + Momentum postMom() const { + const G4ThreeVector& p = post->GetMomentum(); + return Momentum(p.x(),p.y(),p.z()); + } + const G4VTouchable* preTouchable() const { + return pre->GetTouchable(); + } + const G4VTouchable* postTouchable() const { + return post->GetTouchable(); + } + const char* volName(G4StepPoint* p, const char* undefined="") const { + G4VPhysicalVolume* v = volume(p); + return v ? v->GetName().c_str() : undefined; + } + G4VPhysicalVolume* volume(G4StepPoint* p) const { + return p->GetTouchableHandle()->GetVolume(); + } + G4VPhysicalVolume* physvol(G4StepPoint* p) const { + return p->GetPhysicalVolume(); + } + G4LogicalVolume* logvol(G4StepPoint* p) const { + G4VPhysicalVolume* pv = physvol(p); + return pv ? pv->GetLogicalVolume() : 0; + } + G4VSensitiveDetector* sd(G4StepPoint* p) const { + G4LogicalVolume* lv = logvol(p); + return lv ? lv->GetSensitiveDetector() : 0; + } + const char* sdName(G4StepPoint* p, const char* undefined="") const { + G4VSensitiveDetector* s = sd(p); + return s ? s->GetName().c_str() : undefined; + } + + G4VPhysicalVolume* preVolume() const { return volume(pre); } + G4VSensitiveDetector* preSD() const { return sd(pre); } + G4VPhysicalVolume* postVolume() const { return volume(post); } + G4VSensitiveDetector* postSD() const { return sd(post); } + }; + + } // End namespace Simulation +} // End namespace DD4hep + +#endif // DD4HEP_GEANT4HITS_H diff --git a/DDG4/src/Geant4Converter.cpp b/DDG4/src/Geant4Converter.cpp index 031fca75f08c0aaa43bd71f61570320d04310dd0..024530453812ac3791984f63a66ef6a16dff06dc 100644 --- a/DDG4/src/Geant4Converter.cpp +++ b/DDG4/src/Geant4Converter.cpp @@ -41,15 +41,16 @@ #include "G4VUserRegionInformation.hh" // Geant4 include files #include "G4Element.hh" +#include "G4SDManager.hh" #include "G4Box.hh" -#include "G4Tubs.hh" #include "G4Trd.hh" -#include "G4Paraboloid.hh" +#include "G4Tubs.hh" +#include "G4Torus.hh" +#include "G4Sphere.hh" #include "G4Polycone.hh" #include "G4Polyhedra.hh" -#include "G4Sphere.hh" -#include "G4Torus.hh" #include "G4UnionSolid.hh" +#include "G4Paraboloid.hh" #include "G4SubtractionSolid.hh" #include "G4IntersectionSolid.hh" @@ -76,10 +77,9 @@ namespace { static string indent = ""; struct MyTransform3D : public G4Transform3D { MyTransform3D(double XX, double XY, double XZ, double DX, - double YX, double YY, double YZ, double DY, - double ZX, double ZY, double ZZ, double DZ) : G4Transform3D() - {} - // { G4Transform3D::setTransform(XX,XY,XZ,DX,YX,YY,YZ,DY,ZX,ZY,ZZ,DZ); } + double YX, double YY, double YZ, double DY, + double ZX, double ZY, double ZZ, double DZ) + : G4Transform3D(XX,XY,XZ,DX,YX,YY,YZ,DY,ZX,ZY,ZZ,DZ) {} }; void handleName(const TGeoNode* n) { @@ -108,6 +108,17 @@ namespace { } +/// Initializing Constructor +Geant4Converter::Geant4Converter( LCDD& lcdd ) : m_lcdd(lcdd) { + m_checkOverlaps = true; +} +#if 0 +/// Dump element in GDML format to output stream +void* Geant4Converter::printElement(const string& name, const TGeoElement* element) const { + G4Element* g4e = G4Element::GetElement(name,false); + +} +#endif /// Dump element in GDML format to output stream void* Geant4Converter::handleElement(const string& name, const TGeoElement* element) const { G4Element* g4e = data().g4Elements[element]; @@ -126,12 +137,9 @@ void* Geant4Converter::handleElement(const string& name, const TGeoElement* elem } } else { - G4double z = element->Z(); - G4double aeff = element->A()*(g/mole); - cout << name << " " << g/mole << " " << element->A() << " " << aeff/(g/mole) - << " " << (aeff/(g/mole) - z) << endl; - g4e = new G4Element(element->GetTitle(),name,z,aeff); + g4e = new G4Element(element->GetTitle(),name,element->Z(),element->A()*(g/mole)); } + cout << "Created G4 " << (*g4e) << endl; } data().g4Elements[element] = g4e; } @@ -144,14 +152,33 @@ void* Geant4Converter::handleMaterial(const string& name, const TGeoMedium* medi if ( !mat ) { mat = G4Material::GetMaterial(name,false); if ( !mat ) { - int nElements = 0; - const char* opt = "----"; TGeoMaterial* m = medium->GetMaterial(); + G4State state = kStateUndefined; + double density = m->GetDensity()*(gram/cm3); + + switch(m->GetState()) { + case TGeoMaterial::kMatStateSolid: + state = kStateSolid; + break; + case TGeoMaterial::kMatStateLiquid: + state = kStateLiquid; + break; + case TGeoMaterial::kMatStateGas: + state = kStateGas; + break; + default: + case TGeoMaterial::kMatStateUndefined: + state = kStateUndefined; + break; + } if ( m->IsMixture() ) { double A_total = 0.0; TGeoMixture* mix = (TGeoMixture*)m; - nElements = mix->GetNelements(); - mat = new G4Material(name,mix->GetDensity(),nElements); + int nElements = mix->GetNelements(); + mat = new G4Material(name,density,nElements, + state, + m->GetTemperature(), + m->GetPressure()); for(int i=0; i<nElements; ++i) A_total += (mix->GetAmixt())[i]; for(int i=0; i<nElements; ++i) { @@ -163,13 +190,12 @@ void* Geant4Converter::handleMaterial(const string& name, const TGeoMedium* medi } mat->AddElement(g4e,(mix->GetAmixt())[i]/A_total); } - opt = "mixture"; } else { - mat = new G4Material(name,m->GetZ(),m->GetA(),m->GetDensity()); - opt = "raw_material"; + mat = new G4Material(name,m->GetZ(),m->GetA(),density,state, + m->GetTemperature(),m->GetPressure()); } - cout << "Created G4 Material " << m->GetName() << " as " << opt << " " << nElements << endl; + cout << "Created G4 " << *mat << endl; } data().g4Materials[medium] = mat; } @@ -182,11 +208,15 @@ void* Geant4Converter::handleSolid(const string& name, const TGeoShape* shape) if ( !solid && shape ) { if ( shape->IsA() == TGeoBBox::Class() ) { const TGeoBBox* s = (const TGeoBBox*)shape; - solid = new G4Box(name,s->GetDX()*CM_2_MM,s->GetDY()*CM_2_MM,s->GetDZ()*CM_2_MM); + G4Box* box = new G4Box(name,s->GetDX()*CM_2_MM,s->GetDY()*CM_2_MM,s->GetDZ()*CM_2_MM); + solid = box; + //::printf("ROOT Box: %s x=%f y=%f z=%f\n",name.c_str(),s->GetDX(),s->GetDY(),s->GetDZ()); + //::printf(" +-->G4 Box: %s x=%f y=%f z=%f\n",name.c_str(),box->GetXHalfLength(),box->GetYHalfLength(),box->GetZHalfLength()); } else if ( shape->IsA() == TGeoTube::Class() ) { const TGeoTube* s = (const TGeoTube*)shape; solid = new G4Tubs(name,s->GetRmin()*CM_2_MM,s->GetRmax()*CM_2_MM,s->GetDz()*CM_2_MM,0,2.*M_PI); + //::printf("Convert Tube: %s r=%f r=%f z=%f\n",name.c_str(),s->GetRmin()*CM_2_MM,s->GetRmax()*CM_2_MM,s->GetDz()*CM_2_MM,0,2.*M_PI); } else if ( shape->IsA() == TGeoTubeSeg::Class() ) { const TGeoTubeSeg* s = (const TGeoTubeSeg*)shape; @@ -257,9 +287,9 @@ void* Geant4Converter::handleSolid(const string& name, const TGeoShape* shape) } if ( m->IsRotation() ) { - MyTransform3D transform(r[0],r[1],r[2],t[0], - r[3],r[4],r[5],t[1], - r[6],r[7],r[8],t[3]); + MyTransform3D transform(r[0],r[1],r[2],t[0]*CM_2_MM, + r[3],r[4],r[5],t[1]*CM_2_MM, + r[6],r[7],r[8],t[2]*CM_2_MM); if ( oper == TGeoBoolNode::kGeoSubtraction ) solid = new G4SubtractionSolid(name,left,right,transform); else if ( oper == TGeoBoolNode::kGeoUnion ) @@ -268,7 +298,7 @@ void* Geant4Converter::handleSolid(const string& name, const TGeoShape* shape) solid = new G4IntersectionSolid(name,left,right,transform); } else { - G4ThreeVector transform(t[0],t[1],t[2]); + G4ThreeVector transform(t[0]*CM_2_MM,t[1]*CM_2_MM,t[2]*CM_2_MM); if ( oper == TGeoBoolNode::kGeoSubtraction ) solid = new G4SubtractionSolid(name,left,right,0,transform); else if ( oper == TGeoBoolNode::kGeoUnion ) @@ -290,12 +320,14 @@ void* Geant4Converter::handleSolid(const string& name, const TGeoShape* shape) /// Dump logical volume in GDML format to output stream void* Geant4Converter::handleVolume(const string& name, const TGeoVolume* volume) const { - G4LogicalVolume* vol = data().g4Volumes[volume]; + G4GeometryInfo& info = data(); + G4LogicalVolume* vol = info.g4Volumes[volume]; if ( !vol ) { const TGeoVolume* v = volume; - string n = v->GetName(); - TGeoMedium* m = v->GetMedium(); - TGeoShape* s = v->GetShape(); + Volume _v = Ref_t(v); + string n = v->GetName(); + TGeoMedium* m = v->GetMedium(); + TGeoShape* s = v->GetShape(); G4VSolid* solid = (G4VSolid*)handleSolid(s->GetName(),s); G4Material* medium = (G4Material*)handleMaterial(m->GetName(),m); @@ -305,70 +337,103 @@ void* Geant4Converter::handleVolume(const string& name, const TGeoVolume* volume if ( !medium ) { throw runtime_error("G4Converter: No Geant4 material present for volume:"+n); } - vol = new G4LogicalVolume(solid, medium, n); - data().g4Volumes[v] = vol; - - // Here we collect all information we need later to handle to volume internals - Volume _v = Ref_t(v); - Region reg = _v.region(); - LimitSet lim = _v.limitSet(); + //Region reg = _v.region(); SensitiveDetector det = _v.sensitiveDetector(); + Geant4SensitiveDetector* sd = 0; + if ( det.isValid() ) { + sd = info.g4SensDets[det.ptr()]; + if ( !sd ) { + throw runtime_error("G4Cnv::volume["+name+"]: + FATAL Failed to " + "access Geant4 sensitive detector."); + } + sd->Activate(true); + } + LimitSet lim = _v.limitSet(); + G4UserLimits* l = 0; + if ( lim.isValid() ) { + l = info.g4Limits[lim.ptr()]; + if ( !l ) { + throw runtime_error("G4Cnv::volume["+name+"]: + FATAL Failed to " + "access Geant4 user limits."); + } + } + vol = new G4LogicalVolume(solid,medium,n,0,sd,l); + info.g4Volumes[v] = vol; + if ( sd ) { + cout << "G4Cnv::volume: + " << name << " <> " << vol->GetName() + << " Solid:" << solid->GetName() << " Mat:" << medium->GetName() + << " SD:" << sd->GetName() + << endl; - if ( lim.isValid() ) data().limits[lim.ptr()].insert(v); - if ( reg.isValid() ) data().regions[reg.ptr()].insert(v); - if ( det.isValid() ) data().sensitives[det.ptr()].insert(v); - + } //cout << "Converted logical volume [" << n << "]:" << v.ptr() << " ---> G4 " << vol << endl; } return vol; } +/// Dump logical volume in GDML format to output stream +void* Geant4Converter::collectVolume(const string& name, const TGeoVolume* volume) const { + G4GeometryInfo& info = data(); + const TGeoVolume* v = volume; + Volume _v = Ref_t(v); + Region reg = _v.region(); + LimitSet lim = _v.limitSet(); + SensitiveDetector det = _v.sensitiveDetector(); + + if ( lim.isValid() ) info.limits[lim.ptr()].insert(v); + if ( reg.isValid() ) info.regions[reg.ptr()].insert(v); + if ( det.isValid() ) info.sensitives[det.ptr()].insert(v); + return (void*)v; +} + /// Dump volume placement in GDML format to output stream void* Geant4Converter::handlePlacement(const string& name, const TGeoNode* node) const { - G4PVPlacement* g4 = data().g4Placements[node]; + G4GeometryInfo& info = data(); + G4PVPlacement* g4 = info.g4Placements[node]; if ( !g4 ) { TGeoMatrix* trafo = node->GetMatrix(); if ( trafo ) { const Double_t* trans = trafo->GetTranslation(); - const Double_t* rot = trafo->GetRotationMatrix(); int copy = node->GetNumber(); - G4LogicalVolume* g4vol = data().g4Volumes[node->GetVolume()]; - G4LogicalVolume* g4mot = data().g4Volumes[node->GetMotherVolume()]; + G4LogicalVolume* vol = info.g4Volumes[node->GetVolume()]; + G4LogicalVolume* mot = info.g4Volumes[node->GetMotherVolume()]; const Value<TGeoNodeMatrix,PlacedVolume::Object>* obj = dynamic_cast<const Value<TGeoNodeMatrix,PlacedVolume::Object>* >(node); - if ( 0 == g4vol ) { + if ( 0 == vol ) { cout << "FATAL: Unknown G4 volume:" << (void*)obj << " " << obj->GetName() << endl; } else if ( trafo->IsRotation() ) { - MyTransform3D transform(rot[0],rot[1],rot[2],trans[0], - rot[3],rot[4],rot[5],trans[1], - rot[6],rot[7],rot[8],trans[3]); + const Double_t* rot = trafo->GetRotationMatrix(); + MyTransform3D transform(rot[0],rot[1],rot[2],trans[0]*CM_2_MM, + rot[3],rot[4],rot[5],trans[1]*CM_2_MM, + rot[6],rot[7],rot[8],trans[2]*CM_2_MM); g4 = new G4PVPlacement(transform, // no rotation - g4vol, // its logical volume + vol, // its logical volume name, // its name - g4mot, // its mother (logical) volume + mot, // its mother (logical) volume false, // no boolean operations - copy); // its copy number + copy, // its copy number + m_checkOverlaps); } else { - G4ThreeVector pos(trans[0],trans[1],trans[2]); + G4ThreeVector pos(trans[0]*CM_2_MM,trans[1]*CM_2_MM,trans[2]*CM_2_MM); g4 = new G4PVPlacement(0, // no rotation pos, // translation position - g4vol, // its logical volume + vol, // its logical volume name, // its name - g4mot, // its mother (logical) volume + mot, // its mother (logical) volume false, // no boolean operations - copy); // its copy number + copy, // its copy number + m_checkOverlaps); } data().g4Placements[node] = g4; - // cout << "Created volume placement:" << name << " No:" << copy << endl; } } else { cout << "Attempt to DOUBLE-place physical volume:" << name << " No:" << node->GetNumber() << endl; } return g4; -} + } /// Convert the geometry type region into the corresponding Geant4 object(s). void* Geant4Converter::handleRegion(const TNamed* region, const set<const TGeoVolume*>& volumes) const { @@ -440,46 +505,28 @@ void* Geant4Converter::handleLimitSet(const TNamed* limitset, const set<const TG /// Convert the geometry type SensitiveDetector into the corresponding Geant4 object(s). void* Geant4Converter::handleSensitive(const TNamed* sens_det, const set<const TGeoVolume*>& volumes) const { - Geant4SensitiveDetector* g4 = data().g4SensDets[sens_det]; + G4GeometryInfo& info = data(); + Geant4SensitiveDetector* g4 = info.g4SensDets[sens_det]; if ( !g4 ) { - G4GeometryInfo& info = data(); SensitiveDetector sd = Ref_t(sens_det); - /// LCDD: Sensitive detector type == item tag name string type = sd.type(), name = sd.name(); + g4 = ROOT::Reflex::PluginService::Create<Geant4SensitiveDetector*>(type,name,&m_lcdd); if ( !g4 ) { throw runtime_error("Geant4Converter<SensitiveDetector>: FATAL Failed to " "create Geant4 sensitive detector "+name+" of type "+type+"."); } + g4->Activate(true); g4->defineCollection(sd.hitsCollection()); - ConstVolumeSet& volset = info.sensitives[sens_det]; - for(ConstVolumeSet::iterator i=volset.begin(); i!=volset.end();++i) { - map<const TGeoVolume*, G4LogicalVolume*>::iterator v = info.g4Volumes.find(*i); - if ( v == info.g4Volumes.end() ) { - throw runtime_error("Geant4Converter<SensitiveDetector>: FATAL Missing converted " - "Geant 4 logical volume"); - } - (*v).second->SetSensitiveDetector(g4); - } - if ( sd.verbose() ) { - cout << "Geant4Converter<SensitiveDetector> +" << setw(18) << left << sd.name() - << setw(20) << left << " ["+sd.type()+"]" - << " Hits:" << setw(16) << left << sd.hitsCollection() << endl; - cout << " | " - << "Cutoff:" << setw(6) << left << sd.energyCutoff() - << setw(5) << right << volset.size() << " volumes "; - if ( sd.region().isValid() ) cout << " Region:" << setw(12) << left << sd.region().name(); - if ( sd.limits().isValid() ) cout << " Limits:" << setw(12) << left << sd.limits().name(); - cout << "." << endl; - } - data().g4SensDets[sens_det] = g4; + G4SDManager::GetSDMpointer()->AddNewDetector(g4); + info.g4SensDets[sens_det] = g4; } return g4; } /// Handle the geant 4 specific properties void Geant4Converter::handleProperties(LCDD::Properties& prp) const { - map<string,string> processors; + map<string,string> processors; static int s_idd = 9999999; string id; for(LCDD::Properties::const_iterator i=prp.begin(); i!=prp.end(); ++i) { @@ -517,6 +564,82 @@ void Geant4Converter::handleProperties(LCDD::Properties& prp) const { } } + +/// Convert the geometry type SensitiveDetector into the corresponding Geant4 object(s). +void* Geant4Converter::printSensitive(const TNamed* sens_det, const set<const TGeoVolume*>& volumes) const { + G4GeometryInfo& info = data(); + Geant4SensitiveDetector* g4 = info.g4SensDets[sens_det]; + ConstVolumeSet& volset = info.sensitives[sens_det]; + SensitiveDetector sd = Ref_t(sens_det); + bool verbose = sd.verbose(); + + if ( verbose ) { + cout << "Geant4Converter<SensitiveDetector> +" << setw(18) << left << sd.name() + << setw(20) << left << " ["+sd.type()+"]" + << " Hits:" << setw(16) << left << sd.hitsCollection() << endl; + cout << " | " + << "Cutoff:" << setw(6) << left << sd.energyCutoff() + << setw(5) << right << volset.size() << " volumes "; + if ( sd.region().isValid() ) cout << " Region:" << setw(12) << left << sd.region().name(); + if ( sd.limits().isValid() ) cout << " Limits:" << setw(12) << left << sd.limits().name(); + cout << "." << endl; + } + for(ConstVolumeSet::iterator i=volset.begin(); i!=volset.end();++i) { + map<const TGeoVolume*, G4LogicalVolume*>::iterator v = info.g4Volumes.find(*i); + G4LogicalVolume* vol = (*v).second; + if ( verbose ) { + cout << " | " + << "Volume:" << setw(24) << left << vol->GetName() + << " " << vol->GetNoDaughters() << " daughters." + << endl; + } + } + return g4; +} + +void printSolid(G4VSolid* sol) { + if ( typeid(*sol) == typeid(G4Box) ) { + const G4Box* b = (G4Box*)sol; + cout << " Box: x=" << b->GetXHalfLength() << " y=" << b->GetYHalfLength() << " z=" << b->GetZHalfLength(); + } + if ( typeid(*sol) == typeid(G4Tubs) ) { + const G4Tubs* t = (const G4Tubs*)sol; + cout << " Tubs: Ri=" << t->GetInnerRadius() << " Ra=" << t->GetOuterRadius() + << " z/2=" << t->GetZHalfLength() << " Phi=" << t->GetStartPhiAngle() + << "..." << t->GetDeltaPhiAngle (); + } +} + +/// Print G4 placement +void* Geant4Converter::printPlacement(const string& name, const TGeoNode* node) const { + G4GeometryInfo& info = data(); + G4PVPlacement* g4 = info.g4Placements[node]; + G4LogicalVolume* vol = info.g4Volumes[node->GetVolume()]; + G4LogicalVolume* mot = info.g4Volumes[node->GetMotherVolume()]; + G4VSolid* sol = vol->GetSolid(); + G4ThreeVector tr = g4->GetObjectTranslation(); + + G4VSensitiveDetector* sd = vol->GetSensitiveDetector(); + if ( !sd ) return g4; + + cout << "G4Cnv::placement: + " << name << " No:" << node->GetNumber() + << " Vol:" << vol->GetName() << " Solid:" << sol->GetName() + << endl; + cout << " |" + << " Loc: x=" << tr.x() << " y=" << tr.y() << " z=" << tr.z(); + printSolid(sol); + cout << endl; + cout << " |" + << " Ndau:" << vol->GetNoDaughters() << " physvols." + << " Mat:" << vol->GetMaterial()->GetName() + << " Mother:" << (char*)(mot ? mot->GetName().c_str() : "---") + << endl; + cout << " |" + << " SD:" << (char*)(sd ? sd->GetName().c_str() : "---") + << endl; + return g4; +} + template <typename O, typename C, typename F> void handle(const O* o, const C& c, F pmf) { for(typename C::const_iterator i=c.begin(); i != c.end(); ++i) { (o->*pmf)((*i)->GetName(),*i); @@ -528,35 +651,60 @@ template <typename O, typename C, typename F> void handleMap(const O* o, const C (o->*pmf)((*i).first, (*i).second); } +template <typename O, typename C, typename F> void handleRMap(const O* o, const C& c, F pmf) { + for(typename C::const_reverse_iterator i=c.rbegin(); i != c.rend(); ++i) + handle(o, (*i).second, pmf); +} + +/// Create geometry conversion void Geant4Converter::create(DetElement top) { LCDD& lcdd = m_lcdd; G4GeometryInfo& geo = *(m_dataPtr=new G4GeometryInfo); m_data->clear(); collect(top,geo); + m_checkOverlaps = false; + // We do not have to handle defines etc. // All positions and the like are not really named. // Hence, start creating the G4 objects for materials, solids and log volumes. handle(this, geo.materials, &Geant4Converter::handleMaterial); cout << "++ Handled " << geo.materials.size() << " materials." << endl; + + handle(this, geo.volumes, &Geant4Converter::collectVolume); + handle(this, geo.solids, &Geant4Converter::handleSolid); cout << "++ Handled " << geo.solids.size() << " solids." << endl; - handle(this, geo.volumes, &Geant4Converter::handleVolume); - cout << "++ Handled " << geo.volumes.size() << " volumes." << endl; - //handleMap(this, lcdd.fields(), &Geant4Converter::handleField); - //cout << "++ Handled " << geo.fields.size() << " field entries." << endl; + + handleMap(this, geo.sensitives, &Geant4Converter::handleSensitive); + cout << "++ Handled " << geo.sensitives.size() << " sensitive detectors." << endl; handleMap(this, geo.limits, &Geant4Converter::handleLimitSet); cout << "++ Handled " << geo.limits.size() << " limit sets." << endl; + handleMap(this, geo.regions, &Geant4Converter::handleRegion); cout << "++ Handled " << geo.regions.size() << " regions." << endl; - handleMap(this, geo.sensitives, &Geant4Converter::handleSensitive); - cout << "++ Handled " << geo.sensitives.size() << " sensitive detectors." << endl; + + handle(this, geo.volumes, &Geant4Converter::handleVolume); + cout << "++ Handled " << geo.volumes.size() << " volumes." << endl; // Now place all this stuff appropriately - for(Data::const_reverse_iterator i=m_data->rbegin(); i != m_data->rend(); ++i) - handle(this, (*i).second, &Geant4Converter::handlePlacement); + handleRMap(this, *m_data, &Geant4Converter::handlePlacement); //==================== Fields handleProperties(m_lcdd.properties()); + + //cout << *(G4Material::GetMaterialTable()) << endl; + //handleMap(this, geo.sensitives, &Geant4Converter::printSensitive); + //handleRMap(this, *m_data, &Geant4Converter::printPlacement); +} + +/// Singleton instance +Geant4Converter& Geant4Converter::instance() { + static Geant4Converter* inst = 0; + if ( 0 == inst ) { + Geometry::LCDD& lcdd = LCDD::getInstance(); + inst = new Geant4Converter(lcdd); + } + return *inst; } diff --git a/DDG4/src/Geant4Hits.cpp b/DDG4/src/Geant4Hits.cpp index 5e652cd1006038a6ee22fd2eee30e5b71ffbd9da..7e43dbd301bd4a3695b0a4672715488c836fd473 100644 --- a/DDG4/src/Geant4Hits.cpp +++ b/DDG4/src/Geant4Hits.cpp @@ -29,7 +29,7 @@ bool Geant4Hit::isGeantino(G4Track* track) { G4ParticleDefinition* def = track->GetDefinition(); if ( def == G4ChargedGeantino::Definition() ) return true; - if ( def != G4Geantino::Definition() ) { + if ( def == G4Geantino::Definition() ) { return true; } } diff --git a/DDG4/src/Geant4SensitiveDetector.cpp b/DDG4/src/Geant4SensitiveDetector.cpp index 4e5abd94a6d0b58582f93b191e85e9259cefd1c5..ca296215a9c4082fae0f44cc786bb8dfd0bfda5f 100644 --- a/DDG4/src/Geant4SensitiveDetector.cpp +++ b/DDG4/src/Geant4SensitiveDetector.cpp @@ -9,9 +9,16 @@ // Framework include files #include "DDG4/Geant4SensitiveDetector.h" +#include "DDG4/Geant4Converter.h" #include "DDG4/Geant4Hits.h" #include "DD4hep/LCDD.h" +#include "G4Step.hh" +#include "G4PVPlacement.hh" + +#include "TGeoNode.h" + + // C/C++ include files #include <iostream> #include <stdexcept> @@ -68,10 +75,10 @@ Geant4SensitiveDetector::find<Geant4CalorimeterHit>(const HitCollection* c,const /// Method invoked at the begining of each event. void Geant4SensitiveDetector::Initialize(G4HCofThisEvent* HCE) { int count = 0; - m_hce = 0; + m_hce = HCE; for(G4CollectionNameVector::const_iterator i=collectionName.begin(); i!=collectionName.end();++i,++count) { G4VHitsCollection* c = createCollection(*i); - HCE->AddHitsCollection(GetCollectionID(count),c); + m_hce->AddHitsCollection(GetCollectionID(count),c); } } @@ -83,8 +90,10 @@ void Geant4SensitiveDetector::EndOfEvent(G4HCofThisEvent* HCE) { /// Method for generating hit(s) using the information of G4Step object. G4bool Geant4SensitiveDetector::ProcessHits(G4Step* step,G4TouchableHistory* hist) { - if ( step->GetTotalEnergyDeposit() > m_sensitive.energyCutoff() ) { + double ene_cut = m_sensitive.energyCutoff(); + if ( step->GetTotalEnergyDeposit() > ene_cut ) { if ( !Geant4Hit::isGeantino(step->GetTrack()) ) { + dumpStep(step, hist); return buildHits(step,hist); } } @@ -100,7 +109,7 @@ Geant4SensitiveDetector::HitCollection* Geant4SensitiveDetector::collectionByID( /// Retrieve the hits collection associated with this detector by its serial number Geant4SensitiveDetector::HitCollection* Geant4SensitiveDetector::collection(int which) { - if ( which < collectionName.size() && which < 0 ) { + if ( which < collectionName.size() && which >= 0 ) { HitCollection* hc = (HitCollection*)m_hce->GetHC(GetCollectionID(which)); if ( hc ) return hc; throw runtime_error("The collection index for subdetector "+name()+" is wrong!"); @@ -111,3 +120,32 @@ Geant4SensitiveDetector::HitCollection* Geant4SensitiveDetector::collection(int /// Method is invoked if the event abortion is occured. void Geant4SensitiveDetector::clear() { } + +/// Dump Step information (careful: very verbose) +void Geant4SensitiveDetector::dumpStep(G4Step* st, G4TouchableHistory* /* history */) { + Geant4StepHandler step(st); + Geant4Converter& cnv = Geant4Converter::instance(); + Geant4Converter::G4GeometryInfo& data = cnv.data(); + + Position pos1 = step.prePos(); + Position pos2 = step.postPos(); + Momentum mom = step.postMom(); + ::printf(" Track:%08ld Pos:(%8f %8f %8f) -> (%f %f %f) Mom:%7.0f %7.0f %7.0f \n", + long(step.track), pos1.x, pos1.y, pos1.z, pos2.x, pos2.y, pos2.z, mom.x, mom.y, mom.z); + ::printf(" pre-Vol: %s Status:%s\n", + step.preVolume()->GetName().c_str(), step.preStepStatus()); + ::printf(" post-Vol:%s Status:%s\n", + step.postVolume()->GetName().c_str(), step.postStepStatus()); + const G4VPhysicalVolume* pv = step.volume(step.post); + + typedef Geant4Converter::PlacementMap Places; + const Places& places = cnv.data().g4Placements; + for(Places::const_iterator i=places.begin(); i!=places.end();++i) { + const G4PVPlacement* pl = (*i).second; + const G4VPhysicalVolume* qv = pl; + if ( qv == pv ) { + const TGeoNode* tpv = (*i).first; + printf(" Found TGeoNode:%s!\n",tpv->GetName()); + } + } +} diff --git a/DDG4/src/Geant4StepHandler.cpp b/DDG4/src/Geant4StepHandler.cpp new file mode 100644 index 0000000000000000000000000000000000000000..ae379aeb3a2ab434954e39690c8933bc5af2a4af --- /dev/null +++ b/DDG4/src/Geant4StepHandler.cpp @@ -0,0 +1,43 @@ +// $Id:$ +//==================================================================== +// AIDA Detector description implementation +//-------------------------------------------------------------------- +// +// Author : M.Frank +// +//==================================================================== + +// Framework include files +#include "DDG4/Geant4StepHandler.h" +using namespace DD4hep::Simulation; + +const char* Geant4StepHandler::stepStatus(G4StepStatus status) { + switch(status) { + // Step reached the world boundary + case fWorldBoundary: return "WorldBoundary"; + // Step defined by a geometry boundary + case fGeomBoundary: return "GeomBoundary"; + // Step defined by a PreStepDoItVector + case fAtRestDoItProc: return "AtRestDoItProc"; + // Step defined by a AlongStepDoItVector + case fAlongStepDoItProc: return "AlongStepDoItProc"; + // Step defined by a PostStepDoItVector + case fPostStepDoItProc: return "PostStepDoItProc"; + // Step defined by the user Step limit in the logical volume + case fUserDefinedLimit: return "UserDefinedLimit"; + // Step defined by an exclusively forced PostStepDoIt process + case fExclusivelyForcedProc: return "ExclusivelyForcedProc"; + // Step not defined yet + case fUndefined: + default: return "Undefined"; + }; +} + +const char* Geant4StepHandler::preStepStatus() const { + return stepStatus(pre ? pre->GetStepStatus() : fUndefined); +} + +const char* Geant4StepHandler::postStepStatus() const { + return stepStatus(post ? post->GetStepStatus() : fUndefined); +} +