From 96a31054078f37870349067199d48adf6a3d552c Mon Sep 17 00:00:00 2001
From: Chengdong Fu <fucd@ihep.ac.cn>
Date: Tue, 21 Jul 2020 13:51:12 +0800
Subject: [PATCH] add interface/convertor to get detector parameters/materials
---
Detector/DetInterface/DetInterface/IGeoSvc.h | 29 +-
Detector/GeoSvc/src/GeoSvc.cpp | 533 ++++++++++++++++++-
Detector/GeoSvc/src/GeoSvc.h | 66 ++-
3 files changed, 600 insertions(+), 28 deletions(-)
diff --git a/Detector/DetInterface/DetInterface/IGeoSvc.h b/Detector/DetInterface/DetInterface/IGeoSvc.h
index 58df046d..7e26a243 100644
--- a/Detector/DetInterface/DetInterface/IGeoSvc.h
+++ b/Detector/DetInterface/DetInterface/IGeoSvc.h
@@ -11,12 +11,28 @@
#define IGEOSVC_H
#include "GaudiKernel/IService.h"
+#include "DDRec/DetectorData.h"
+#include <map>
namespace dd4hep {
-class Detector;
-class DetElement;
+ class Detector;
+ class DetElement;
+ // if not include DDRec/DetectorData.h, can not compile, why?
+ //namespace rec{
+ //struct StructExtension;
+ //struct ZPlanarStruct;
+ //typedef StructExtension<ZPlanarStruct> ZPlanarData;
+ //struct ConicalSupportStruct;
+ //typedef StructExtension<ConicalSupportStruct> ConicalSupportData;
}
+class StructExtension;
+
+namespace gear{
+ class ZPlanarParametersImpl;
+ class GearParametersImpl;
+}
+class TMaterial;
// class G4VUserDetectorConstruction;
class GAUDI_API IGeoSvc : virtual public IService {
@@ -30,6 +46,15 @@ public:
// receive Geant4 Geometry
// virtual G4VUserDetectorConstruction* getGeant4Geo() = 0;
+ // obsolete parameter format, will remove once StructExtension<> validated
+ virtual const gear::ZPlanarParametersImpl* getVXDParameters() = 0;
+
+ virtual const dd4hep::rec::ZPlanarData* getVXDData() = 0;
+ virtual const dd4hep::rec::ConicalSupportData* getBeamPipeData() =0;
+
+ virtual const std::map<std::string,double>& getDetParameters(std::string s) = 0;
+ virtual const double getDetParameter(std::string set_name, std::string par_name) = 0;
+ virtual TMaterial* getMaterial(std::string s) = 0;
virtual ~IGeoSvc() {}
};
diff --git a/Detector/GeoSvc/src/GeoSvc.cpp b/Detector/GeoSvc/src/GeoSvc.cpp
index 305896c4..42990b30 100644
--- a/Detector/GeoSvc/src/GeoSvc.cpp
+++ b/Detector/GeoSvc/src/GeoSvc.cpp
@@ -1,14 +1,28 @@
#include "GeoSvc.h"
+#include "gearimpl/GearParametersImpl.h"
+#include "TMath.h"
+#include "TMaterial.h"
+#include "CLHEP/Units/SystemOfUnits.h"
+
+#include "G4Box.hh"
+#include "G4Tubs.hh"
+#include "G4LogicalVolume.hh"
+#include "G4PVPlacement.hh"
+#include "G4NistManager.hh"
#include "DD4hep/Detector.h"
#include "DD4hep/Plugins.h"
#include "DDG4/Geant4Converter.h"
#include "DDG4/Geant4Mapping.h"
+#include "DDRec/DetectorData.h"
+
+#include <iomanip>
+#include <iostream>
DECLARE_COMPONENT(GeoSvc)
GeoSvc::GeoSvc(const std::string& name, ISvcLocator* svc)
-: base_class(name, svc), m_dd4hep_geo(nullptr) {
+: base_class(name, svc), m_dd4hep_geo(nullptr), m_vxdData(nullptr), m_beamData(nullptr){
}
@@ -18,20 +32,40 @@ GeoSvc::~GeoSvc() {
StatusCode
GeoSvc::initialize() {
- StatusCode sc = Service::initialize();
+ StatusCode sc = Service::initialize();
- m_dd4hep_geo = &(dd4hep::Detector::getInstance());
- // if failed to load the compact, a runtime error will be thrown.
- m_dd4hep_geo->fromCompact(m_dd4hep_xmls.value());
-
- return sc;
+ m_dd4hep_geo = &(dd4hep::Detector::getInstance());
+ // if failed to load the compact, a runtime error will be thrown.
+ m_dd4hep_geo->fromCompact(m_dd4hep_xmls.value());
+
+ dd4hep::DetElement world = m_dd4hep_geo->world();
+ //info() << world.type() << " " << world.path() << " " << world.placementPath() << endmsg;
+ const std::map<std::string, dd4hep::DetElement>& subs = world.children();
+ for(std::map<std::string, dd4hep::DetElement>::const_iterator it=subs.begin();it!=subs.end();it++){
+ dd4hep::DetElement sub = it->second;
+ info() << it->first << " " << sub.path() << " " << sub.placementPath() << endmsg;
+ if(it->first=="Tube"){
+ dd4hep::Volume vol = sub.volume();
+ dd4hep::Solid solid = vol.solid();
+ //info() << " " << solid.type() << " " << solid.name() << endmsg;
+ const std::map<std::string, dd4hep::DetElement>& pipes = sub.children();
+ for(std::map<std::string, dd4hep::DetElement>::const_iterator it=pipes.begin();it!=pipes.end();it++){
+ dd4hep::DetElement pipe = it->second;
+ //info() << " " << it->first << " " << pipe.id() << " " << pipe.path() << " " << pipe.placementPath() << endmsg;
+ }
+ }
+ if(it->first=="VXD"){
+ sc = convertVXD(sub);
+ }
+ }
+ return sc;
}
StatusCode
GeoSvc::finalize() {
- StatusCode sc;
-
- return sc;
+ StatusCode sc;
+ if(m_vxdParameters) delete m_vxdParameters;
+ return sc;
}
dd4hep::DetElement
@@ -46,3 +80,482 @@ dd4hep::Detector*
GeoSvc::lcdd() {
return m_dd4hep_geo;
}
+
+const std::map<std::string,double>& GeoSvc::getDetParameters(std::string name){
+ if(m_detParameters.find(name)!=m_detParameters.end()) return m_detParameters[name];
+ else{
+ char message[200];
+ sprintf(message,"GeoSvc has not the parameter set named %s", name);
+ throw std::runtime_error(message);
+ }
+}
+
+const double GeoSvc::getDetParameter(std::string set_name, std::string par_name){
+ std::map<std::string, std::map<std::string,double> >::iterator it=m_detParameters.find(set_name);
+ if(it!=m_detParameters.end()){
+ if(it->second.find(par_name)!=it->second.end()) return it->second[par_name];
+ }
+ char message[200];
+ sprintf(message,"GeoSvc has not the parameter named %s in set %s", par_name, set_name);
+ throw std::runtime_error(message);
+}
+
+StatusCode GeoSvc::convertVXD(dd4hep::DetElement& vxd){
+ StatusCode sc;
+ //fucd: another method to obtain parameters, but not fully for KalDet
+ bool extensionDataValid = true;
+ try{
+ m_vxdData = vxd.extension<dd4hep::rec::ZPlanarData>();
+ }
+ catch(std::runtime_error& e){
+ extensionDataValid = false;
+ info() << e.what() << " " << m_vxdData << endmsg;
+ }
+
+ std::vector<helpLayer> helpSensitives;
+ std::vector<helpLayer> helpLadders;
+ std::vector<int> helpNumberLadders;
+ std::vector<double> helpPhi0;
+ int helpCount=0;
+ int type=0;
+ double shellInnerRadius, shellOuterRadius, shellHalfLength, gap, shellRadLength;
+ int nLadders=0;
+ double phi0=0;
+ helpLayer thisLadder;
+ double beryllium_ladder_block_length=0,end_electronics_half_z=0,side_band_electronics_width=0;
+ double rAlu, drAlu, rSty, drSty, dzSty, rInner, aluEndcapZ, aluHalfZ, alu_RadLen, Cryostat_dEdx;
+ double VXDSupportDensity, VXDSupportZeff, VXDSupportAeff, VXDSupportRadLen;
+ dd4hep::Volume vxd_vol = vxd.volume();
+ for(int i=0;i<vxd_vol->GetNdaughters();i++){
+ TGeoNode* daughter = vxd_vol->GetNode(i);
+ std::string nodeName = daughter->GetName();
+ //info() << daughter->GetName() << endmsg;
+ if(nodeName=="VXD_support_assembly_0"){
+ TGeoNode* shell = FindNode(daughter, "SupportShell");
+ if(shell){
+ const TGeoShape* shape_shell = shell->GetVolume()->GetShape();
+ //fucd: IsA() method does not always work for TGeoTube, sometimes, strange?
+ //if(shape_shell->IsA()==TGeoTube::Class()){
+ if(shape_shell->TestShapeBit(TGeoTube::kGeoTube)){
+ const TGeoTube* tube = (const TGeoTube*) shape_shell;
+ shellInnerRadius = tube->GetRmin()*CLHEP::cm;
+ shellOuterRadius = tube->GetRmax()*CLHEP::cm;
+ shellHalfLength = tube->GetDz()*CLHEP::cm;
+ }
+ else{
+ error() << shell->GetName() << " is not a TGeoTube! Shape bits = " << shape_shell->TestShapeBits(0xFFFFFFFF) << endmsg;
+ }
+ TGeoMaterial* mat = shell->GetMedium()->GetMaterial();
+ shellRadLength = mat->GetRadLen()*CLHEP::cm;
+ }
+ TGeoNode* block = FindNode(daughter, "BerylliumAnnulusBlock");
+ if(block){
+ const TGeoShape* shape_block = block->GetVolume()->GetShape();
+ if(shape_block->IsA()==TGeoBBox::Class()){
+ const TGeoBBox* box = (const TGeoBBox*) shape_block;
+ beryllium_ladder_block_length = box->GetDY()*CLHEP::cm;
+ }
+ else{
+ error() << block->GetName() << " is not a TGeoTube! Shape bits = " << shape_block->TestShapeBits(0xFFFFFFFF) << endmsg;
+ }
+ }
+ TGeoNode* skin = FindNode(daughter, "CryostatAluSkinBarrel");
+ if(skin){
+ const TGeoShape* shape_skin = skin->GetVolume()->GetShape();
+ if(shape_skin->TestShapeBit(TGeoTube::kGeoTube)){
+ const TGeoTube* tube = (const TGeoTube*) shape_skin;
+ rAlu = tube->GetRmin()*CLHEP::cm;
+ drAlu = tube->GetRmax()*CLHEP::cm - rAlu;
+ aluHalfZ = tube->GetDz()*CLHEP::cm;
+ //info() << rmin << "," << rmax << "," << zhalf << endmsg;
+ }
+ else{
+ error() << skin->GetName() << " is not a TGeoTube! Shape bits = " << shape_skin->TestShapeBits(0xFFFFFFFF) << endmsg;
+ }
+ }
+ TGeoNode* foam = FindNode(daughter, "CryostatFoamBarrel");
+ if(foam){
+ const TGeoShape* shape_foam = foam->GetVolume()->GetShape();
+ if(shape_foam->TestShapeBit(TGeoTube::kGeoTube)){
+ const TGeoTube* tube = (const TGeoTube*) shape_foam;
+ rSty = tube->GetRmin()*CLHEP::cm;
+ drSty = tube->GetRmax()*CLHEP::cm - rSty;
+ dzSty = tube->GetDz()*CLHEP::cm;
+ //info() << rmin << "," << rmax << "," << zhalf << endmsg;
+ }
+ else{
+ error() << foam->GetName() << " is not a TGeoTube! Shape bits = " << shape_foam->TestShapeBits(0xFFFFFFFF) << endmsg;
+ }
+ }
+ TGeoNode* skinEnd = FindNode(daughter, "CryostatAluSkinEndPlateInner");
+ if(skinEnd){
+ const TGeoShape* shape_skinEnd = skinEnd->GetVolume()->GetShape();
+ if(shape_skinEnd->TestShapeBit(TGeoTube::kGeoTube)){
+ const TGeoTube* tube = (const TGeoTube*) shape_skinEnd;
+ rInner = tube->GetRmin()*CLHEP::cm;
+ double rmax = tube->GetRmax()*CLHEP::cm;
+ drAlu = tube->GetDz()*CLHEP::cm*2;
+ //info() << rmin << "," << rmax << "," << zhalf << endmsg;
+ }
+ else{
+ error() << skinEnd->GetName() << " is not a TGeoTube! Shape bits = " << shape_skinEnd->TestShapeBits(0xFFFFFFFF) << endmsg;
+ }
+ }
+ TGeoNode* shellEnd = FindNode(daughter, "EndPlateShell_outer");
+ if(shellEnd){
+ const TGeoShape* shape_shellEnd = shellEnd->GetVolume()->GetShape();
+ if(shape_shellEnd->TestShapeBit(TGeoTube::kGeoTube)){
+ const TGeoTube* tube = (const TGeoTube*) shape_shellEnd;
+ double rmin = tube->GetRmin()*CLHEP::cm;
+ double rmax = tube->GetRmax()*CLHEP::cm;
+ double zhalf = tube->GetDz()*CLHEP::cm;
+ //info() << rmin << "," << rmax << "," << zhalf << endmsg;
+ }
+ else{
+ error() << shellEnd->GetName() << " is not a TGeoTube! Shape bits = " << shape_shellEnd->TestShapeBits(0xFFFFFFFF) << endmsg;
+ }
+ }
+
+ }
+ else if(nodeName=="layer_assembly_0_1"){
+ if(TGeoNode* side_band = FindNode(daughter, "ElectronicsBand")){
+ const TGeoShape* shape_band = side_band->GetVolume()->GetShape();
+ if(shape_band->IsA()==TGeoBBox::Class()){
+ const TGeoBBox* box = (const TGeoBBox*) shape_band;
+ side_band_electronics_width = box->GetDX()*CLHEP::cm*2;
+ //info() << "fucd: "<< box->GetDX() << " " << box->GetDY() << " " << box->GetDZ() <<endmsg;
+ }
+ else{
+ error() << "ElectronicsBand is not a TGeoBBox!!!"<< endmsg;
+ }
+ }
+ if(TGeoNode* end = FindNode(daughter, "ElectronicsEnd")){
+ const TGeoShape* shape_end = end->GetVolume()->GetShape();
+ if(shape_end->IsA()==TGeoBBox::Class()){
+ const TGeoBBox* box = (const TGeoBBox*) shape_end;
+ end_electronics_half_z = box->GetDY()*CLHEP::cm;
+ //info() << "fucd: " << box->GetDX() << " " << box->GetDY() << " " << box->GetDZ() << endmsg;
+ }
+ else{
+ error() << "ElectronicsEnd is not a TGeoBBox!!!"<< endmsg;
+ }
+ }
+ }
+
+ /*
+ for(int j=0;j<daughter->GetNdaughters();j++){
+ TGeoNode* next = daughter->GetDaughter(j);
+ info() << "fucd: " << next->GetName() << endmsg;
+ }
+ */
+ }
+
+ const std::map<std::string, dd4hep::DetElement>& components = vxd.children();
+ for(std::map<std::string, dd4hep::DetElement>::const_iterator it=components.begin();it!=components.end();it++){
+ dd4hep::DetElement component = it->second;
+ dd4hep::Volume vol = component.volume();
+ dd4hep::PlacedVolume phys = component.placement();
+ TGeoMaterial* mat = vol->GetMaterial();
+ const TGeoShape* shape = vol->GetShape();
+ const dd4hep::PlacedVolumeExtension::VolIDs& ids = phys.volIDs();
+ //info() << " " << it->first << " " << vol->GetName() << " " << component.id() << " " << component.path() << " " << component.placementPath() << endmsg;
+ //info() << " " << shape->GetName() << " " << vol.solid().name() << endmsg;
+ //info() << " " << mat->GetName() << " " << mat->GetRadLen() << endmsg;
+ //info() << " " << ids.str() << endmsg;
+ //info() << " " << vol->GetNdaughters() << endmsg;
+ //double radL = mat->GetRadLen();
+ //dd4hep::Solid solid = vol.solid();
+ //info() << " " << sh->TestShapeBit(TGeoShape::kGeoBox) << " " << sh->GetName() << " " << phys.material().radLength() << endmsg;
+ if(vol.isSensitive()&&shape->IsA()==TGeoBBox::Class()){
+ int iLayer = ids.find("layer")->second;
+ int iModule = ids.find("module")->second;
+ int iSide = ids.find("side")->second;
+ //info() << " layer=" << iLayer << " module=" << iModule << mat->GetName() << endmsg;
+ if(iModule==0&&iLayer==helpCount+1){
+ helpCount++;
+ helpSensitives.push_back(thisLadder);
+ helpLadders.push_back(thisLadder);
+ helpNumberLadders.push_back(nLadders);
+ helpPhi0.push_back(phi0);
+ nLadders = 0;
+ thisLadder.length = 0;
+ }
+ if(iLayer == helpCount){
+ if(iModule == 0){
+ const TGeoBBox* box = (const TGeoBBox*) shape;
+ double width = box->GetDX()*CLHEP::cm;
+ double length = box->GetDY()*CLHEP::cm;
+ double thickness = box->GetDZ()*CLHEP::cm;
+ TGeoMatrix* matrix = phys->GetMatrix();
+ const double* pos = matrix->GetTranslation();
+ const double* rot_data = matrix->GetRotationMatrix();
+ TGeoRotation rot;
+ rot.SetMatrix(rot_data);
+ double theta,phi,psi;
+ rot.GetAngles(phi,theta,psi);
+ phi *= TMath::DegToRad();
+ theta *= TMath::DegToRad();
+ psi *= TMath::DegToRad();
+ phi0 = -TMath::PiOver2()+phi;
+ double distance = fabs(cos(phi0)*sin(theta)*pos[0]+sin(phi0)*sin(theta)*pos[1]+cos(theta)*pos[2]);
+ double offset = sqrt(pos[0]*pos[0]+pos[1]*pos[1]-distance*distance)*pos[0]/fabs(pos[0])*CLHEP::cm;
+ distance *= CLHEP::cm;
+ distance -= thickness;
+ double radL = mat->GetRadLen()*CLHEP::cm;
+ //info() << " -> " << helpCount << ": " << distance << " " << cos(atan2(pos[1],pos[0])-phi)*sqrt(pos[0]*pos[0]+pos[1]*pos[1]) << endmsg;
+ thisLadder.distance = distance;
+ thisLadder.offset = offset;
+ thisLadder.thickness = thickness;
+ thisLadder.length += length;
+ thisLadder.width = width;
+ thisLadder.radLength = radL;
+ thisLadder.z = pos[2]*CLHEP::cm;
+ }
+ if(iModule==nLadders) nLadders++;
+ }
+ }
+ //info() << " " << vol.solid().type() << " " << vol.solid().name() << " " << vol->GetNdaughters() << endmsg;
+ else if(it->first=="VXD_support"){
+ helpCount++;
+ helpSensitives.push_back(thisLadder);
+ helpLadders.push_back(thisLadder);
+ helpNumberLadders.push_back(nLadders);
+ helpPhi0.push_back(phi0);
+ nLadders = 0;
+ if(vol->GetNdaughters()==0) error() << "!!!!!!!!!" << endmsg;
+
+ int nFlexCable = 0, nFoamSpacer=0, nMetalTraces=0;
+ int currentLayer = -1;
+ double tFlexCable, tFoamSpacer, tMetalTraces;
+ double radLFlexCable, radLFoamSpacer, radLMetalTraces;
+ double intLFlexCable, intLFoamSpacer, intLMetalTraces;
+ double dFlexCable, dFoamSpacer, dMetalTraces;
+ double metalZeff, metalZAeff, foamZeff, foamZAeff, flexZeff, flexZAeff;
+ for(int i=0;i<vol->GetNdaughters();i++){
+ TGeoNode* daughter = vol->GetNode(i);
+ TGeoMaterial* matDaughter = daughter->GetMedium()->GetMaterial();
+ const TGeoShape* shape_sup = daughter->GetVolume()->GetShape();
+ TGeoMatrix* matrix = daughter->GetMatrix();
+ const double* pos = matrix->GetTranslation();
+ const double* rot_data = matrix->GetRotationMatrix();
+ TGeoRotation rot;
+ rot.SetMatrix(rot_data);
+ double theta,phi,psi;
+ rot.GetAngles(phi,theta,psi);
+ phi *= TMath::DegToRad();
+ theta *= TMath::DegToRad();
+ psi *= TMath::DegToRad();
+ phi0 = -TMath::PiOver2()+phi;
+ std::string phy_name = daughter->GetName();
+ if(phy_name.find("FoamSpacer")==-1&&phy_name.find("FlexCable")==-1&&phy_name.find("MetalTraces")==-1){
+ //info() << phy_name <<endmsg;
+ continue;
+ }
+ int iLayer = atoi(phy_name.substr(phy_name.find("_")+1,2).c_str());
+ if(iLayer!=currentLayer){
+ //info() << tFoamSpacer << "," << tFlexCable << "," << tMetalTraces << endmsg;
+ helpLadders[currentLayer].thickness = tFoamSpacer+tFlexCable+tMetalTraces;
+ helpLadders[currentLayer].radLength = helpLadders[currentLayer].thickness / (tFoamSpacer/radLFoamSpacer+tFlexCable/radLFlexCable+tMetalTraces/radLMetalTraces);
+ nFlexCable = 0;
+ nFoamSpacer=0;
+ nMetalTraces=0;
+ currentLayer=iLayer;
+ }
+ //info() << "ss pos=" << pos[0] << "," << pos[1] << "," << pos[2] << " distance=" << sqrt(pos[0]*pos[0]+pos[1]*pos[1]) << endmsg;
+ //info() << "ss rot=" << phi << "," << theta << "," << psi << endmsg;
+ //info() << "ss " << daughter->GetName() << " " << daughter->GetVolume()->GetName() << " " << endmsg;
+ if(shape_sup->IsA()==TGeoBBox::Class()&&(nFoamSpacer==0||nFlexCable==0||nMetalTraces==0)){
+ const TGeoBBox* box = (const TGeoBBox*) shape_sup;
+ //info() << phy_name.substr(phy_name.find("_")+1,2) << " " << iLayer << " " << box->GetDX() << "," << box->GetDY() << "," << box->GetDZ() << endmsg;
+ //info() << "fucd: pos " << pos[0] << " " << pos[1] << " " << pos[2] << endmsg;
+ double distance = fabs(cos(phi0)*sin(theta)*pos[0]+sin(phi0)*sin(theta)*pos[1]+cos(theta)*pos[2]);
+ double offset = sqrt(pos[0]*pos[0]+pos[1]*pos[1]-distance*distance)*pos[0]/fabs(pos[0])*CLHEP::cm;
+ distance -= box->GetDZ();
+ distance *= CLHEP::cm;
+ if(helpLadders[iLayer].distance == helpSensitives[iLayer].distance) helpLadders[iLayer].distance = distance;
+ else helpLadders[iLayer].distance = TMath::Min(helpLadders[iLayer].distance, distance);
+ //info() << phy_name << " " << distance << " " << offset << endmsg;
+ if(phy_name.find("FoamSpacer")!=-1&&nFoamSpacer==0){
+ helpLadders[iLayer].offset = offset;
+ tFoamSpacer = box->GetDZ()*CLHEP::cm;
+ radLFoamSpacer = matDaughter->GetRadLen()*CLHEP::cm;
+ intLFoamSpacer = matDaughter->GetIntLen()*CLHEP::cm;
+ dFoamSpacer = matDaughter->GetDensity()*CLHEP::g/CLHEP::cm3;
+ //fucd: A calculated by TGeoMaterial class is not equal to Zeff/ZAeff, Zeff = sum(Zi*Ai/sumA), ZAeff = sum(Zi/Ai*A/totalA)
+ // use Zeff and ZAeff to keep same with Mokka case
+ //foamZ = matDaughter->GetZ();
+ //foamA = matDaughter->GetA();
+ double totalA = 0, Zeff = 0, ZAeff = 0;
+ for(int iEle = 0; iEle<matDaughter->GetNelements(); iEle++){
+ totalA += matDaughter->GetElement(iEle)->A();
+ }
+ for(int iEle = 0; iEle<matDaughter->GetNelements(); iEle++){
+ double A, Z, w;
+ // by fucd: w!=A/totalA, strange! to fix
+ matDaughter->GetElementProp(A,Z,w,iEle);
+ Zeff += Z*w;
+ ZAeff += Z/A*w;
+ //info() << std::setprecision(16) << Z << " " << A << " " << A/totalA << " " << w << endmsg;
+ }
+ foamZeff = Zeff;
+ foamZAeff = ZAeff;
+ nFoamSpacer++;
+ }
+ if(phy_name.find("FlexCable")!=-1&&nFlexCable==0){
+ helpLadders[iLayer].width = box->GetDX()*CLHEP::cm;
+ helpLadders[iLayer].length = box->GetDY()*CLHEP::cm-beryllium_ladder_block_length*2-end_electronics_half_z*2;
+ tFlexCable = box->GetDZ()*CLHEP::cm;
+ radLFlexCable = matDaughter->GetRadLen()*CLHEP::cm;
+ intLFlexCable = matDaughter->GetIntLen()*CLHEP::cm;
+ dFlexCable = matDaughter->GetDensity()*CLHEP::g/CLHEP::cm3;
+ double Zeff = 0, ZAeff = 0;
+ for(int iEle = 0; iEle<matDaughter->GetNelements(); iEle++){
+ double A, Z, w;
+ matDaughter->GetElementProp(A,Z,w,iEle);
+ Zeff += Z*w;
+ ZAeff += Z/A*w;
+ //std::cout << std::setprecision(16) << Z << " " << A << " " << w << std::endl;
+ }
+ flexZeff = Zeff;
+ flexZAeff = ZAeff;
+ nFlexCable++;
+ }
+ if(phy_name.find("MetalTraces")!=-1&&nMetalTraces==0){
+ tMetalTraces = box->GetDZ()*CLHEP::cm;
+ radLMetalTraces = matDaughter->GetRadLen()*CLHEP::cm;
+ intLMetalTraces = matDaughter->GetIntLen()*CLHEP::cm;
+ dMetalTraces = matDaughter->GetDensity()*CLHEP::g/CLHEP::cm3;
+ double totalA = 0, Zeff = 0, ZAeff = 0;
+ for(int iEle = 0; iEle<matDaughter->GetNelements(); iEle++){
+ totalA += matDaughter->GetElement(iEle)->A();
+ }
+ for(int iEle = 0; iEle<matDaughter->GetNelements(); iEle++){
+ double A, Z, w;
+ matDaughter->GetElementProp(A,Z,w,iEle);
+ Zeff += Z*w;
+ ZAeff += Z/A*w;
+ //info() << Z << " " << A << " " << w << endmsg;
+ }
+ metalZeff = Zeff;
+ metalZAeff = ZAeff;
+ nMetalTraces++;
+ }
+ }
+ }
+ {
+ //info() << tFoamSpacer << "," << tFlexCable << "," << tMetalTraces << endmsg;
+ double tSupport = tMetalTraces + tFoamSpacer + tFlexCable;
+ helpLadders[currentLayer].thickness = tSupport;
+ helpLadders[currentLayer].radLength = helpLadders[currentLayer].thickness / (tFoamSpacer/radLFoamSpacer+tFlexCable/radLFlexCable+tMetalTraces/radLMetalTraces);
+ nFlexCable = 0;
+ nFoamSpacer=0;
+ nMetalTraces=0;
+
+ //calculations of thickness fractions of each layer of the support
+ double metalTF = tMetalTraces / tSupport;
+ double foamTF = tFoamSpacer / tSupport;
+ double flexTF = tFlexCable / tSupport;
+ //info() << foamTF << "," << flexTF << "," << metalTF << endmsg;
+ //info() << dFoamSpacer/(CLHEP::kg/CLHEP::cm3) << "," << dFlexCable/(CLHEP::kg/CLHEP::cm3) << "," << dMetalTraces/(CLHEP::kg/CLHEP::cm3) << endmsg;
+ //info() << foamZeff << " " << flexZeff << " " << metalZeff << endmsg;
+ //info() << foamZAeff << " " << flexZAeff << " " << metalZAeff << endmsg;
+ G4double elemVol = 1*CLHEP::cm3;
+ G4double VXDSupportMass = (foamTF*dFoamSpacer + flexTF*dFlexCable + metalTF*dMetalTraces)*elemVol;
+ VXDSupportDensity = VXDSupportMass/elemVol;
+
+ G4double foamFM = 100. * ((foamTF*(elemVol)*dFoamSpacer) / VXDSupportMass) ;
+ G4double kaptonFM = 100. * ((flexTF*(elemVol)*dFlexCable) / VXDSupportMass) ;
+ G4double metalFM = 100. * ((metalTF*(elemVol)*dMetalTraces) / VXDSupportMass) ;
+
+ VXDSupportRadLen = helpLadders[currentLayer].radLength;
+
+ VXDSupportZeff = (metalFM/100.)*metalZeff + (kaptonFM/100.)*flexZeff + (foamFM/100.)*foamZeff;
+ G4double VXDSupportZAeff = (metalFM/100.)*metalZAeff + (kaptonFM/100.)*flexZAeff + (foamFM/100.)*foamZAeff;
+ VXDSupportAeff = VXDSupportZeff / VXDSupportZAeff;
+ G4double VXDSupportIntLength = 1. / ((metalTF/intLMetalTraces) + (flexTF/intLFlexCable) + (foamTF/intLFoamSpacer));
+ VXDSupportDensity = 1000*VXDSupportDensity/(CLHEP::g/CLHEP::cm3);
+ //info() << "fucd: " << VXDSupportZeff << " " << VXDSupportAeff << " " << VXDSupportRadLen << " " << VXDSupportIntLength << " " << VXDSupportDensity << endmsg;
+ //info() << intLMetalTraces << " " << intLFlexCable << " " << intLFoamSpacer <<endmsg;
+ }
+ }
+ //info() << it->first << endmsg;
+ }
+ if(end_electronics_half_z>0 && side_band_electronics_width==0) type = gear::ZPlanarParametersImpl::CCD ;
+ if(side_band_electronics_width>0 && end_electronics_half_z==0 ) type = gear::ZPlanarParametersImpl::CMOS ;
+ if(side_band_electronics_width>0 && end_electronics_half_z>0) type = gear::ZPlanarParametersImpl::HYBRID ;
+
+ m_vxdParameters = new gear::ZPlanarParametersImpl(type, shellInnerRadius, shellOuterRadius, shellHalfLength, gap, shellRadLength );
+ // by fucd: debug info, if validated enough, merge them in future
+ info() << "=====================from convertor==============================" << endmsg;
+ info() << type << " " << shellInnerRadius << " " << shellOuterRadius << " " << shellHalfLength << " " << gap << " " << shellRadLength << endmsg;
+ for(int i=0;i<helpCount;i++){
+ m_vxdParameters->addLayer(helpNumberLadders[i] , helpPhi0[i] ,
+ helpLadders[i].distance , helpLadders[i].offset, helpLadders[i].thickness*2 ,
+ helpLadders[i].length , helpLadders[i].width*2 , helpLadders[i].radLength ,
+ helpSensitives[i].distance, helpSensitives[i].offset , helpSensitives[i].thickness*2 ,
+ helpSensitives[i].length , helpSensitives[i].width*2 , helpSensitives[i].radLength ) ;
+ info() << "fucd " << i << ": " << helpNumberLadders[i] << ", " << helpPhi0[i] << ", "
+ << helpLadders[i].distance << ", " << helpLadders[i].offset << ", " << helpLadders[i].thickness*2 << ", " << helpLadders[i].length << ", "
+ << helpLadders[i].width*2 << ", " << helpLadders[i].radLength << ", " << helpSensitives[i].distance << ", " << helpSensitives[i].offset << ", "
+ << helpSensitives[i].thickness*2 << ", " << helpSensitives[i].length << ", " << helpSensitives[i].width*2 << ", " << helpSensitives[i].radLength << endmsg;
+ }
+ //m_vxdInfra = new gear::GearParametersImpl;
+ //CryostatAlRadius, CryostatAlThickness, CryostatAlInnerR, CryostatAlZEndCap, CryostatAlHalfZ
+ //m_vxdInfra->setDoubleVal("CryostatAlRadius",rAlu);
+ //m_vxdInfra->setDoubleVal("CryostatAlThickness",drAlu);
+ //m_vxdInfra->setDoubleVal("CryostatAlInnerR",rInner);
+ //m_vxdInfra->setDoubleVal("CryostatAlZEndCap",aluEndcapZ=dzSty + drSty + drAlu / 2);
+ //m_vxdInfra->setDoubleVal("CryostatAlHalfZ",aluHalfZ= dzSty + drSty);
+ // change GearParametersImpl to map
+ std::map<std::string,double> vxdInfra;
+ vxdInfra["CryostatAlRadius"] = rAlu;
+ vxdInfra["CryostatAlThickness"] = drAlu;
+ vxdInfra["CryostatAlInnerR"] = rInner;
+ vxdInfra["CryostatAlZEndCap"] = dzSty+drSty+drAlu/2;
+ vxdInfra["CryostatAlHalfZ"] = dzSty+drSty;
+ m_detParameters["VXDInfra"] = vxdInfra;
+ //effective A different with what in Mokka, fix them as Mokka's
+ //m_materials["VXDSupportMaterial"] = new TMaterial("VXDSupportMaterial", "", VXDSupportAeff, VXDSupportZeff, VXDSupportDensity, VXDSupportRadLen, 0.);
+ m_materials["VXDSupportMaterial"] = new TMaterial("VXDSupportMaterial", "", 2.075865162e+01, 1.039383117e+01, 2.765900000e+02/1000, 1.014262421e+02, 0.);
+
+ info() << "=====================from ZPlanarData==============================" << endmsg;
+ if(m_vxdData){
+ info() << m_vxdData->rInnerShell << " " << m_vxdData->rOuterShell << " " << m_vxdData->zHalfShell << " " << m_vxdData->gapShell << endmsg;
+ const std::vector<dd4hep::rec::ZPlanarData::LayerLayout>& layers = m_vxdData->layers;
+ for(int i=0;i<layers.size();i++){
+ const dd4hep::rec::ZPlanarData::LayerLayout& thisLayer = layers[i];
+ info() << i << ": " << thisLayer.ladderNumber << "," << thisLayer.phi0 << "," << thisLayer.distanceSupport << "," << thisLayer.offsetSupport << ","
+ << thisLayer.thicknessSupport << "," << thisLayer.zHalfSupport << "," << thisLayer.widthSupport << "," << "NULL,"
+ << thisLayer.distanceSensitive << "," << thisLayer.offsetSensitive << "," << thisLayer.thicknessSensitive << "," << thisLayer.zHalfSensitive << ","
+ << thisLayer.widthSensitive << ",NULL" << endmsg;
+ }
+ }
+ info() << rAlu << " " << drAlu << " " << rInner << " " << aluEndcapZ << " " << aluHalfZ << endmsg;
+ //info() << m_materials["VXDSupportMaterial"] << endmsg;
+ return sc;
+}
+
+TGeoNode* GeoSvc::FindNode(TGeoNode* mother, char* name){
+ TGeoNode* next = 0;
+ if(mother->GetNdaughters()!=0){
+ for(int i=0;i<mother->GetNdaughters();i++){
+ TGeoNode* daughter = mother->GetDaughter(i);
+ std::string s = daughter->GetName();
+ //info() << "current: " << s << " search for" << name << endmsg;
+ if(s.find(name)!=-1){
+ next = daughter;
+ break;
+ }
+ else{
+ next = FindNode(daughter, name);
+ }
+ }
+ }
+ return next;
+}
+
+TMaterial* GeoSvc::getMaterial(std::string name){
+ std::map<std::string, TMaterial*>::const_iterator it = m_materials.find(name);
+ if(it!=m_materials.end()) return it->second;
+ else return 0;
+}
diff --git a/Detector/GeoSvc/src/GeoSvc.h b/Detector/GeoSvc/src/GeoSvc.h
index 9c98daa9..cc54142f 100644
--- a/Detector/GeoSvc/src/GeoSvc.h
+++ b/Detector/GeoSvc/src/GeoSvc.h
@@ -15,27 +15,61 @@
// DD4Hep
#include "DD4hep/Detector.h"
-class GeoSvc: public extends<Service, IGeoSvc> {
-public:
- GeoSvc(const std::string& name, ISvcLocator* svc);
- ~GeoSvc();
+#include <gear/GEAR.h>
+#include <gearimpl/ZPlanarParametersImpl.h>
+#include <gearimpl/GearParametersImpl.h>
+
+class dd4hep::DetElement;
+class TGeoNode;
- // Service
- StatusCode initialize() override;
- StatusCode finalize() override;
+class GeoSvc: public extends<Service, IGeoSvc> {
+ public:
+ GeoSvc(const std::string& name, ISvcLocator* svc);
+ ~GeoSvc();
+
+ // Service
+ StatusCode initialize() override;
+ StatusCode finalize() override;
+
+ // IGeoSvc
+ dd4hep::DetElement getDD4HepGeo() override;
+ dd4hep::Detector* lcdd() override;
+
+ const gear::ZPlanarParametersImpl* getVXDParameters() override {return m_vxdParameters;};
+ const dd4hep::rec::ZPlanarData* getVXDData() override {return m_vxdData;};
+ const dd4hep::rec::ConicalSupportData* getBeamPipeData() override {return m_beamData;};
- // IGeoSvc
- dd4hep::DetElement getDD4HepGeo() override;
- dd4hep::Detector* lcdd() override;
+ const std::map<std::string,double>& getDetParameters(std::string name) override;
+ const double getDetParameter(std::string set_name, std::string par_name) override;
+ TMaterial* getMaterial(std::string name);
+
+ private:
+ StatusCode convertVXD(dd4hep::DetElement& sub);
-private:
+ TGeoNode* FindNode(TGeoNode* mother, char* name);
+ // DD4hep XML compact file path
+ Gaudi::Property<std::string> m_dd4hep_xmls{this, "compact"};
+
+ //
+ dd4hep::Detector* m_dd4hep_geo;
- // DD4hep XML compact file path
- Gaudi::Property<std::string> m_dd4hep_xmls{this, "compact"};
+ gear::ZPlanarParametersImpl* m_vxdParameters;
+ dd4hep::rec::ZPlanarData* m_vxdData;
+ dd4hep::rec::ConicalSupportData* m_beamData;
- //
- dd4hep::Detector* m_dd4hep_geo;
+ //gear::GearParametersImpl* m_vxdInfra;
+ std::map<std::string, std::map<std::string,double> > m_detParameters;
+ std::map<std::string, TMaterial*> m_materials;
+ struct helpLayer {
+ double distance =0;
+ double offset =0;
+ double thickness =0;
+ double length =0;
+ double width =0;
+ double radLength =0;
+ double z =0;
+ double foam_spacer_radLength =0;
+ };
};
-
#endif GeoSvc_h
--
GitLab