From 9da2b6d1c3ae33eca0364757bd6140410c28ec4f Mon Sep 17 00:00:00 2001
From: Markus Frank <Markus.Frank@cern.ch>
Date: Wed, 14 Aug 2019 11:50:37 +0200
Subject: [PATCH] Fix errors in Solid::setDimensions/Solid::dimensions
---
DDCore/src/Shapes.cpp | 388 ++++++++++++++++++++----
DDCore/src/plugins/ShapePlugins.cpp | 10 +-
DDG4/include/DDG4/Geant4SensDetAction.h | 2 +-
DDG4/plugins/Geant4EventReaderHepMC.cpp | 10 +-
DDG4/plugins/Geant4GDMLWriteAction.cpp | 11 +-
5 files changed, 346 insertions(+), 75 deletions(-)
diff --git a/DDCore/src/Shapes.cpp b/DDCore/src/Shapes.cpp
index 8944c7898..465609baa 100644
--- a/DDCore/src/Shapes.cpp
+++ b/DDCore/src/Shapes.cpp
@@ -125,7 +125,7 @@ namespace dd4hep {
}
/// Access Shape dimension parameters (As in TGeo, but angles in radians rather than degrees)
- std::vector<double> get_shape_dimensions(TGeoShape* shape) {
+ vector<double> get_shape_dimensions(TGeoShape* shape) {
if (shape) {
TClass* cl = shape->IsA();
if (cl == TGeoShapeAssembly::Class()) {
@@ -248,25 +248,123 @@ namespace dd4hep {
pars.emplace_back(vtx->xy[i][0]);
pars.emplace_back(vtx->xy[i][1]);
}
-#if 0
- const Double_t* vertices = sh->GetVertices();
- vector<double> pars { sh->GetDz() };
- pars.reserve(17);
- for ( size_t i=0; i<8; ++i )
- pars.emplace_back(vertices[i*2]);
- for ( size_t i=0; i<8; ++i )
- pars.emplace_back(vertices[i*2]+1);
-#endif
return pars;
}
else if (cl == TGeoCompositeShape::Class()) {
+ Solid solid(shape);
const TGeoCompositeShape* sh = (const TGeoCompositeShape*) shape;
const TGeoBoolNode* boolean = sh->GetBoolNode();
- TGeoBoolNode::EGeoBoolType oper = boolean->GetBooleanOperator();
- TGeoMatrix* left_matrix = boolean->GetRightMatrix();
TGeoMatrix* right_matrix = boolean->GetRightMatrix();
TGeoShape* left_solid = boolean->GetLeftShape();
- TGeoShape* right_solid = boolean->GetLeftShape();
+ TGeoShape* right_solid = boolean->GetRightShape();
+ if ( instanceOf<TruncatedTube>(solid) ) {
+ stringstream params(right_matrix->GetTitle());
+ vector<double> pars;
+ pars.reserve(8);
+ for(size_t i=0; i<8; ++i) {
+ double val;
+ params >> val;
+ pars.emplace_back(val);
+ if ( !params.good() ) {
+ except("Solid::dimensions","+++ Invalid parameter stream.");
+ }
+ }
+ return pars;
+#if 0
+ TGeoTubeSeg* tubs = (TGeoTubeSeg*)left_solid;
+ TGeoBBox* box = (TGeoBBox*)right_solid;
+ double rmin = tubs->GetRmin();
+ double rmax = tubs->GetRmax();
+ double zhalf = tubs->GetDz();
+ double startPhi = tubs->GetPhi1()*units::degree;
+ double deltaPhi = tubs->GetPhi2()*units::degree - startPhi;
+ double boxY = box->GetDY();
+ double xBox = right_matrix->GetTranslation()[0];
+ XYZAngles angles = detail::matrix::_xyzAngles(right_matrix);
+ double alpha = -angles.Z();
+ double sin_alpha = std::sin(std::fabs(alpha));
+ double sin2alpha = sin_alpha*sin_alpha;
+ double cos_phi = std::cos(deltaPhi);
+ double cos2phi = cos_phi*cos_phi;
+ double sin2phi = 1.0 - cos2phi;
+ double ss1 = 1.0 - (sin2phi/sin2alpha);
+ double R = 0, r = 0;
+ bool cutInside = true;
+
+ /// Need to select physical solutions: THIS DOES NOT WORK. CANNOT RESOLVE THE AMBIGUITY.
+ double r1 = xBox + boxY/sin_alpha;
+ double r2 = xBox - boxY/sin_alpha;
+ if ( r1 > 0 )
+ r = r1, cutInside = true;
+ else if ( r2 > 0 )
+ r = r2, cutInside = false;
+ else
+ except("","+++ Bad parameters");
+
+ double R1 = -r*cos_phi/ss1*(1.0 + std::sqrt(1.0 + (1.0 - ss1/cos2phi)));
+ double R2 = -r*cos_phi/ss1*(1.0 - std::sqrt(1.0 + (1.0 - ss1/cos2phi)));
+ /// Need to select physical solutions: THIS DOES NOT WORK. CANNOT RESOLVE THE AMBIGUITY.
+ if ( R1 > 0 )
+ R = R1;
+ else if ( R2 > 0 )
+ R = R2;
+ else
+ except("","+++ Bad parameters");
+ double cutAtStart = r;
+ double cutAtDelta = R;
+
+ cout << "dimensions: " << endl;
+ cout << right_matrix->GetTitle() << endl;
+ cout << TRUNCATEDTUBE_TAG << ":" << endl
+ << "\t zhalf: " << zhalf << " " << endl
+ << "\t rmin: " << rmin << " " << endl
+ << "\t rmax: " << rmax << " " << endl
+ << "\t startPhi: " << startPhi/units::degree << " " << endl
+ << "\t deltaPhi: " << deltaPhi/units::degree << " " << endl
+ << "\t cutAtStart:" << cutAtStart << " " << endl
+ << "\t cutAtDelta:" << cutAtDelta << " " << endl
+ << "\t cutInside: " << char(cutInside ? '1' : '0') << endl
+ << "\t\t boxY: " << boxY << endl
+ << "\t\t alpha: " << angles.Z() << " [rad] " << angles.Z()/units::degree << " [degree]" << endl
+ << "\t\t sin_alpha: " << sin_alpha << endl
+ << "\t\t xBox: " << xBox << endl
+ << "\t\t r1: " << r1 << endl
+ << "\t\t r2: " << r2 << endl
+ << "\t\t R1: " << R1 << endl
+ << "\t\t R2: " << R2 << endl
+ ;
+
+ return {
+ zhalf, rmin, rmax,
+ startPhi*units::deg, deltaPhi*units::deg,
+ cutAtStart, cutAtDelta, cutInside ? 1.0 : 0.0 };
+#endif
+ }
+ else if ( instanceOf<PseudoTrap>(solid) ) {
+ stringstream params(right_matrix->GetTitle());
+ vector<double> pars;
+ pars.reserve(7);
+ cout << "dimensions: [" << PSEUDOTRAP_TAG << "]" << endl
+ << right_matrix->GetTitle() << endl;
+ for(size_t i=0; i<7; ++i) {
+ double val;
+ params >> val;
+ pars.emplace_back(val);
+ if ( !params.good() ) {
+ except("Solid::dimensions","+++ Invalid parameter stream.");
+ }
+ }
+ return pars;
+ }
+ else if ( instanceOf<SubtractionSolid>(solid) ) {
+ }
+ else if ( instanceOf<UnionSolid>(solid) ) {
+ }
+ else if ( instanceOf<IntersectionSolid>(solid) ) {
+ }
+
+ TGeoBoolNode::EGeoBoolType oper = boolean->GetBooleanOperator();
+ TGeoMatrix* left_matrix = boolean->GetRightMatrix();
const Double_t* left_tr = left_matrix->GetTranslation();
const Double_t* left_rot = left_matrix->GetRotationMatrix();
const Double_t* right_tr = right_matrix->GetTranslation();
@@ -296,12 +394,12 @@ namespace dd4hep {
}
/// Access Shape dimension parameters (As in TGeo, but angles in radians rather than degrees)
- std::vector<double> get_shape_dimensions(Solid solid) {
+ vector<double> get_shape_dimensions(Solid solid) {
return get_shape_dimensions(solid.ptr());
}
/// Set the shape dimensions. As for the TGeo shape, but angles in rad rather than degrees.
- void set_shape_dimensions(TGeoShape* shape, const std::vector<double>& params) {
+ void set_shape_dimensions(TGeoShape* shape, const vector<double>& params) {
if (shape) {
TClass* cl = shape->IsA();
Solid solid(shape);
@@ -485,10 +583,142 @@ namespace dd4hep {
solid._setDimensions(&pars[0]);
}
else if (cl == TGeoCompositeShape::Class()) {
- //TGeoCompositeShape* sh = (TGeoCompositeShape*) shape;
+ TGeoCompositeShape* sh = (TGeoCompositeShape*) shape;
+ TGeoBoolNode* boolean = sh->GetBoolNode();
+ TGeoMatrix* right_matrix = boolean->GetRightMatrix();
+ TGeoShape* left_solid = boolean->GetLeftShape();
+ TGeoShape* right_solid = boolean->GetRightShape();
if ( instanceOf<TruncatedTube>(solid) ) {
+ TGeoTubeSeg* tubs = (TGeoTubeSeg*)left_solid;
+ TGeoBBox* box = (TGeoBBox*)right_solid;
+ double zhalf = params[0];
+ double rmin = params[1];
+ double rmax = params[2];
+ double startPhi = params[3]/units::deg;
+ double deltaPhi = params[4]/units::deg;
+ double cutAtStart = params[5];
+ double cutAtDelta = params[6];
+ bool cutInside = params[7] > 0.5;
+
+ cout << "setDimensions: [" << TRUNCATEDTUBE_TAG << "]" << endl
+ << right_matrix->GetTitle() << endl;
+ // check the parameters
+ if( rmin <= 0 || rmax <= 0 || cutAtStart <= 0 || cutAtDelta <= 0 )
+ except(TRUNCATEDTUBE_TAG,"++ 0 <= rIn,cutAtStart,rOut,cutAtDelta,rOut violated!");
+ else if( rmin >= rmax )
+ except(TRUNCATEDTUBE_TAG,"++ rIn<rOut violated!");
+ else if( startPhi != 0. )
+ except(TRUNCATEDTUBE_TAG,"++ startPhi != 0 not supported!");
+
+ double r = cutAtStart;
+ double R = cutAtDelta;
+ // angle of the box w.r.t. tubs
+ double cath = r - R * std::cos( deltaPhi*units::deg );
+ double hypo = std::sqrt( r*r + R*R - 2.*r*R * std::cos( deltaPhi*units::deg ));
+ double cos_alpha = cath / hypo;
+ double alpha = std::acos( cos_alpha );
+ double sin_alpha = std::sin( std::fabs(alpha) );
+
+ // exaggerate dimensions - does not matter, it's subtracted!
+ // If we don't, the **edge** of the box would cut into the tube segment
+ // for larger delta-phi values
+ double boxX = 1.1*rmax + rmax/sin_alpha; // Need to adjust for move!
+ double boxY = rmax;
+ // width of the box > width of the tubs
+ double boxZ = 1.1 * zhalf;
+ double xBox; // center point of the box
+ if( cutInside )
+ xBox = r - boxY / sin_alpha;
+ else
+ xBox = r + boxY / sin_alpha;
+
+ // rotationmatrix of box w.r.t. tubs
+ TGeoRotation rot;
+ rot.RotateZ( -alpha/dd4hep::deg );
+ double box_dim[] = {boxX, boxY, boxZ};
+ double tub_dim[] = {rmin, rmax, zhalf, startPhi, deltaPhi};
+ box->SetDimensions(box_dim);
+ tubs->SetDimensions(tub_dim);
+ TGeoCombiTrans* combi = (TGeoCombiTrans*)right_matrix;
+ combi->SetRotation(rot);
+ combi->SetTranslation(xBox, 0, 0);
+ return;
}
else if ( instanceOf<PseudoTrap>(solid) ) {
+ double x1 = params[0];
+ double x2 = params[1];
+ double y1 = params[2];
+ double y2 = params[3];
+ double z = params[4];
+ double r = params[5];
+ bool atMinusZ = params[6] > 0.5;
+ double x = atMinusZ ? x1 : x2;
+ double h = 0;
+ bool intersec = false; // union or intersection solid
+ double displacement = 0;
+ double startPhi = 0;
+ double halfZ = z;
+ double halfOpeningAngle = std::asin( x / std::abs( r ))/units::deg;
+ /// calculate the displacement of the tubs w.r.t. to the trap, determine the opening angle of the tubs
+ double delta = std::sqrt( r * r - x * x );
+
+ cout << "setDimensions: [" << PSEUDOTRAP_TAG << "]" << endl
+ << right_matrix->GetTitle() << endl;
+
+ // Implementation from :
+ // https://cmssdt.cern.ch/lxr/source/Fireworks/Geometry/src/TGeoMgrFromDdd.cc#0538
+ if( r < 0 && std::abs(r) >= x ) {
+ intersec = true; // intersection solid
+ h = y1 < y2 ? y2 : y1; // tubs half height
+ h += h/20.; // enlarge a bit - for subtraction solid
+ if( atMinusZ ) {
+ displacement = - halfZ - delta;
+ startPhi = 90.0 - halfOpeningAngle;
+ }
+ else {
+ displacement = halfZ + delta;
+ startPhi = -90.0 - halfOpeningAngle;
+ }
+ }
+ else if( r > 0 && std::abs(r) >= x ) {
+ if( atMinusZ ) {
+ displacement = - halfZ + delta;
+ startPhi = 270.0 - halfOpeningAngle;
+ h = y1;
+ }
+ else
+ {
+ displacement = halfZ - delta;
+ startPhi = 90.0 - halfOpeningAngle;
+ h = y2;
+ }
+ }
+ else {
+ except(PSEUDOTRAP_TAG,"Check parameters of the PseudoTrap!");
+ }
+ double trd2_dim[] = { x1, x2, y1, y2, halfZ };
+ double tube_dim[] = { 0.0, std::abs(r), h, startPhi, startPhi + halfOpeningAngle*2. };
+ if( intersec && left_solid->IsA() == TGeoTrd2::Class() ) {
+ left_solid->SetDimensions(trd2_dim);
+ right_solid->SetDimensions(tube_dim);
+ }
+ else if ( right_solid->IsA() == TGeoCompositeShape::Class() ) {
+ double box_dim[] = { 1.1*x, 1.1*h, std::sqrt(r*r-x*x) };
+ TGeoCompositeShape* comp = (TGeoCompositeShape*)right_solid;
+ TGeoBoolNode* b_s = comp->GetBoolNode();
+ b_s->GetLeftShape()->SetDimensions(tube_dim);
+ b_s->GetRightShape()->SetDimensions(box_dim);
+ left_solid->SetDimensions(trd2_dim);
+ }
+ else {
+ except("PseudoTrap","+++ Incompatible change of parameters.");
+ }
+ ((TGeoTranslation*)right_matrix)->SetTranslation(0,0,displacement);
+ stringstream params;
+ params << x1 << " " << x2 << " " << y1 << " " << y2 << " " << z << " "
+ << r << " " << char(atMinusZ ? '1' : '0') << " ";
+ right_matrix->SetTitle(params.str().c_str());
+ return;
}
else if ( instanceOf<SubtractionSolid>(solid) ) {
}
@@ -507,7 +737,7 @@ namespace dd4hep {
except("Solid","set_shape_dimensions: Failed to set dimensions [Invalid handle].");
}
/// Set the shape dimensions. As for the TGeo shape, but angles in rad rather than degrees.
- void set_shape_dimensions(Solid solid, const std::vector<double>& params) {
+ void set_shape_dimensions(Solid solid, const vector<double>& params) {
set_shape_dimensions(solid.ptr(), params);
}
}
@@ -645,7 +875,7 @@ string dd4hep::toStringSolid(const TGeoShape* shape, int precision) {
else if (shape->IsA() == TGeoCompositeShape::Class()) {
const TGeoCompositeShape* sh = (const TGeoCompositeShape*) shape;
const TGeoBoolNode* boolean = sh->GetBoolNode();
- const TGeoShape* left = boolean->GetLeftShape();
+ const TGeoShape* left = boolean->GetLeftShape();
const TGeoShape* right = boolean->GetRightShape();
TGeoBoolNode::EGeoBoolType oper = boolean->GetBooleanOperator();
if (oper == TGeoBoolNode::kGeoSubtraction)
@@ -661,7 +891,7 @@ string dd4hep::toStringSolid(const TGeoShape* shape, int precision) {
}
/// Output mesh vertices to string
-std::string dd4hep::toStringMesh(const TGeoShape* shape, int prec) {
+string dd4hep::toStringMesh(const TGeoShape* shape, int prec) {
Solid sol(shape);
stringstream os;
@@ -747,14 +977,14 @@ template <typename T> vector<double> Solid_type<T>::dimensions() {
}
/// Set the shape dimensions. As for the TGeo shape, but angles in rad rather than degrees.
-template <typename T> Solid_type<T>& Solid_type<T>::setDimensions(const std::vector<double>& params) {
+template <typename T> Solid_type<T>& Solid_type<T>::setDimensions(const vector<double>& params) {
set_shape_dimensions(this->access(), params);
return *this;
}
/// Divide volume into subsections (See the ROOT manuloa for details)
template <typename T> TGeoVolume*
-Solid_type<T>::divide(const Volume& voldiv, const std::string& divname,
+Solid_type<T>::divide(const Volume& voldiv, const string& divname,
int iaxis, int ndiv, double start, double step) const {
T* p = this->ptr();
if ( p ) {
@@ -770,7 +1000,7 @@ ShapelessSolid::ShapelessSolid(const string& nam) {
_assign(new TGeoShapeAssembly(), nam, "Assembly", true);
}
-void Box::make(const std::string& nam, double x_val, double y_val, double z_val) {
+void Box::make(const string& nam, double x_val, double y_val, double z_val) {
_assign(new TGeoBBox(nam.c_str(), x_val, y_val, z_val), "", "Box", true);
}
@@ -797,7 +1027,7 @@ double Box::z() const {
}
/// Internal helper method to support object construction
-void HalfSpace::make(const std::string& nam, const double* const point, const double* const normal) {
+void HalfSpace::make(const string& nam, const double* const point, const double* const normal) {
_assign(new TGeoHalfSpace(nam.c_str(),(Double_t*)point, (Double_t*)normal), "", "HalfSpace",true);
}
@@ -848,12 +1078,12 @@ Polycone::Polycone(double startPhi, double deltaPhi, const vector<double>& r, co
}
/// Constructor to be used when creating a new object
-Polycone::Polycone(const std::string& nam, double startPhi, double deltaPhi) {
+Polycone::Polycone(const string& nam, double startPhi, double deltaPhi) {
_assign(new TGeoPcon(nam.c_str(), startPhi/units::deg, deltaPhi/units::deg, 0), "", "Polycone", false);
}
/// Constructor to be used when creating a new polycone object. Add at the same time all Z planes
-Polycone::Polycone(const std::string& nam, double startPhi, double deltaPhi,
+Polycone::Polycone(const string& nam, double startPhi, double deltaPhi,
const vector<double>& rmin, const vector<double>& rmax, const vector<double>& z) {
vector<double> params;
if (rmin.size() < 2) {
@@ -874,7 +1104,7 @@ Polycone::Polycone(const std::string& nam, double startPhi, double deltaPhi,
}
/// Constructor to be used when creating a new polycone object. Add at the same time all Z planes
-Polycone::Polycone(const std::string& nam, double startPhi, double deltaPhi, const vector<double>& r, const vector<double>& z) {
+Polycone::Polycone(const string& nam, double startPhi, double deltaPhi, const vector<double>& r, const vector<double>& z) {
vector<double> params;
if (r.size() < 2) {
throw runtime_error("dd4hep: PolyCone Not enough Z planes. minimum is 2!");
@@ -928,7 +1158,7 @@ ConeSegment::ConeSegment(double dz,
}
/// Constructor to be used when creating a new cone segment object
-ConeSegment::ConeSegment(const std::string& nam,
+ConeSegment::ConeSegment(const string& nam,
double dz,
double rmin1, double rmax1,
double rmin2, double rmax2,
@@ -949,7 +1179,7 @@ ConeSegment& ConeSegment::setDimensions(double dz,
}
/// Constructor to be used when creating a new object with attribute initialization
-void Cone::make(const std::string& nam, double z, double rmin1, double rmax1, double rmin2, double rmax2) {
+void Cone::make(const string& nam, double z, double rmin1, double rmax1, double rmin2, double rmax2) {
_assign(new TGeoCone(nam.c_str(), z, rmin1, rmax1, rmin2, rmax2 ), "", "Cone", true);
}
@@ -979,14 +1209,14 @@ CutTube::CutTube(double rmin, double rmax, double dz, double startPhi, double en
}
/// Constructor to be used when creating a new object with attribute initialization
-CutTube::CutTube(const std::string& nam,
+CutTube::CutTube(const string& nam,
double rmin, double rmax, double dz, double startPhi, double endPhi,
double lx, double ly, double lz, double tx, double ty, double tz) {
make(nam, rmin,rmax,dz,startPhi/units::deg,endPhi/units::deg,lx,ly,lz,tx,ty,tz);
}
/// Constructor to be used when creating a new object with attribute initialization
-void CutTube::make(const std::string& nam, double rmin, double rmax, double dz, double startPhi, double endPhi,
+void CutTube::make(const string& nam, double rmin, double rmax, double dz, double startPhi, double endPhi,
double lx, double ly, double lz, double tx, double ty, double tz) {
_assign(new TGeoCtub(nam.c_str(), rmin,rmax,dz,startPhi,endPhi,lx,ly,lz,tx,ty,tz),"","CutTube",true);
}
@@ -997,13 +1227,13 @@ TruncatedTube::TruncatedTube(double zhalf, double rmin, double rmax, double star
{ make("", zhalf, rmin, rmax, startPhi/units::deg, deltaPhi/units::deg, cutAtStart, cutAtDelta, cutInside); }
/// Constructor to create a truncated tube object with attribute initialization
-TruncatedTube::TruncatedTube(const std::string& nam,
+TruncatedTube::TruncatedTube(const string& nam,
double zhalf, double rmin, double rmax, double startPhi, double deltaPhi,
double cutAtStart, double cutAtDelta, bool cutInside)
{ make(nam, zhalf, rmin, rmax, startPhi/units::deg, deltaPhi/units::deg, cutAtStart, cutAtDelta, cutInside); }
/// Internal helper method to support object construction
-void TruncatedTube::make(const std::string& nam,
+void TruncatedTube::make(const string& nam,
double zhalf, double rmin, double rmax, double startPhi, double deltaPhi,
double cutAtStart, double cutAtDelta, bool cutInside) {
// check the parameters
@@ -1038,15 +1268,39 @@ void TruncatedTube::make(const std::string& nam,
// rotationmatrix of box w.r.t. tubs
TGeoRotation rot;
- //rot.RotateX( 90.0 );
rot.RotateZ( -alpha/dd4hep::deg );
TGeoTranslation trans(xBox, 0., 0.);
TGeoBBox* box = new TGeoBBox((nam+"Box").c_str(), boxX, boxY, boxZ);
TGeoTubeSeg* tubs = new TGeoTubeSeg((nam+"Tubs").c_str(), rmin, rmax, zhalf, startPhi, deltaPhi);
- TGeoSubtraction* sub = new TGeoSubtraction(tubs, box, nullptr, new TGeoCombiTrans(trans, rot));
- // For debugging:
- // TGeoUnion* sub = new TGeoUnion(tubs, box, nullptr, new TGeoCombiTrans(trans, rot));
+ TGeoCombiTrans* combi = new TGeoCombiTrans(trans, rot);
+ TGeoSubtraction* sub = new TGeoSubtraction(tubs, box, nullptr, combi);
_assign(new TGeoCompositeShape(nam.c_str(), sub),"",TRUNCATEDTUBE_TAG,true);
+ stringstream params;
+ params << zhalf << " " << endl
+ << rmin << " " << endl
+ << rmax << " " << endl
+ << startPhi*units::deg << " " << endl
+ << deltaPhi*units::deg << " " << endl
+ << cutAtStart << " " << endl
+ << cutAtDelta << " " << endl
+ << char(cutInside ? '1' : '0') << endl;
+ combi->SetTitle(params.str().c_str());
+ //cout << "Params: " << params.str() << endl;
+#if 0
+ params << TRUNCATEDTUBE_TAG << ":" << endl
+ << "\t zhalf: " << zhalf << " " << endl
+ << "\t rmin: " << rmin << " " << endl
+ << "\t rmax: " << rmax << " " << endl
+ << "\t startPhi: " << startPhi << " " << endl
+ << "\t deltaPhi: " << deltaPhi << " " << endl
+ << "\t r/cutAtStart:" << cutAtStart << " " << endl
+ << "\t R/cutAtDelta:" << cutAtDelta << " " << endl
+ << "\t cutInside: " << char(cutInside ? '1' : '0') << endl
+ << "\t\t alpha: " << alpha << endl
+ << "\t\t sin_alpha: " << sin_alpha << endl
+ << "\t\t boxY: " << boxY << endl
+ << "\t\t xBox: " << xBox << endl;
+#endif
#if 0
cout << "Trans:"; trans.Print(); cout << endl;
cout << "Rot: "; rot.Print(); cout << endl;
@@ -1074,12 +1328,12 @@ void TruncatedTube::make(const std::string& nam,
}
/// Constructor to be used when creating a new object with attribute initialization
-void EllipticalTube::make(const std::string& nam, double a, double b, double dz) {
+void EllipticalTube::make(const string& nam, double a, double b, double dz) {
_assign(new TGeoEltu(nam.c_str(), a, b, dz), "", "EllipticalTube", true);
}
/// Constructor to be used when creating a new object with attribute initialization
-void Trd1::make(const std::string& nam, double x1, double x2, double y, double z) {
+void Trd1::make(const string& nam, double x1, double x2, double y, double z) {
_assign(new TGeoTrd1(nam.c_str(), x1, x2, y, z ), "", "Trd1", true);
}
@@ -1091,7 +1345,7 @@ Trd1& Trd1::setDimensions(double x1, double x2, double y, double z) {
}
/// Constructor to be used when creating a new object with attribute initialization
-void Trd2::make(const std::string& nam, double x1, double x2, double y1, double y2, double z) {
+void Trd2::make(const string& nam, double x1, double x2, double y1, double y2, double z) {
_assign(new TGeoTrd2(nam.c_str(), x1, x2, y1, y2, z ), "", "Trd2", true);
}
@@ -1103,7 +1357,7 @@ Trd2& Trd2::setDimensions(double x1, double x2, double y1, double y2, double z)
}
/// Constructor to be used when creating a new object with attribute initialization
-void Paraboloid::make(const std::string& nam, double r_low, double r_high, double delta_z) {
+void Paraboloid::make(const string& nam, double r_low, double r_high, double delta_z) {
_assign(new TGeoParaboloid(nam.c_str(), r_low, r_high, delta_z ), "", "Paraboloid", true);
}
@@ -1115,7 +1369,7 @@ Paraboloid& Paraboloid::setDimensions(double r_low, double r_high, double delta_
}
/// Constructor to create a new anonymous object with attribute initialization
-void Hyperboloid::make(const std::string& nam, double rin, double stin, double rout, double stout, double dz) {
+void Hyperboloid::make(const string& nam, double rin, double stin, double rout, double stout, double dz) {
_assign(new TGeoHype(nam.c_str(), rin, stin/units::deg, rout, stout/units::deg, dz), "", "Hyperboloid", true);
}
@@ -1127,7 +1381,7 @@ Hyperboloid& Hyperboloid::setDimensions(double rin, double stin, double rout, do
}
/// Constructor function to be used when creating a new object with attribute initialization
-void Sphere::make(const std::string& nam, double rmin, double rmax, double startTheta, double endTheta, double startPhi, double endPhi) {
+void Sphere::make(const string& nam, double rmin, double rmax, double startTheta, double endTheta, double startPhi, double endPhi) {
_assign(new TGeoSphere(nam.c_str(), rmin, rmax,
startTheta/units::deg, endTheta/units::deg,
startPhi/units::deg, endPhi/units::deg), "", "Sphere", true);
@@ -1141,7 +1395,7 @@ Sphere& Sphere::setDimensions(double rmin, double rmax, double startTheta, doubl
}
/// Constructor to be used when creating a new object with attribute initialization
-void Torus::make(const std::string& nam, double r, double rmin, double rmax, double startPhi, double deltaPhi) {
+void Torus::make(const string& nam, double r, double rmin, double rmax, double startPhi, double deltaPhi) {
_assign(new TGeoTorus(nam.c_str(), r, rmin, rmax, startPhi/units::deg, deltaPhi/units::deg), "", "Torus", true);
}
@@ -1162,7 +1416,7 @@ Trap::Trap(double z, double theta, double phi,
}
/// Constructor to be used when creating a new anonymous object with attribute initialization
-Trap::Trap(const std::string& nam,
+Trap::Trap(const string& nam,
double z, double theta, double phi,
double h1, double bl1, double tl1, double alpha1,
double h2, double bl2, double tl2, double alpha2) {
@@ -1172,7 +1426,7 @@ Trap::Trap(const std::string& nam,
}
/// Constructor to be used when creating a new anonymous object with attribute initialization
-void Trap::make(const std::string& nam, double pz, double py, double px, double pLTX) {
+void Trap::make(const string& nam, double pz, double py, double px, double pLTX) {
double z = pz / 2e0;
double theta = 0e0;
double phi = 0e0;
@@ -1197,7 +1451,7 @@ Trap& Trap::setDimensions(double z, double theta, double phi,
}
/// Internal helper method to support object construction
-void PseudoTrap::make(const std::string& nam, double x1, double x2, double y1, double y2, double z, double r, bool atMinusZ) {
+void PseudoTrap::make(const string& nam, double x1, double x2, double y1, double y2, double z, double r, bool atMinusZ) {
double x = atMinusZ ? x1 : x2;
double h = 0;
bool intersec = false; // union or intersection solid
@@ -1285,11 +1539,15 @@ void PseudoTrap::make(const std::string& nam, double x1, double x2, double y1, d
SubtractionSolid sub((nam+"Subs").c_str(), tubs, Box(1.1*x, 1.1*h, std::sqrt(r*r-x*x)), Transform3D(RotationX(M_PI/2.)));
solid = UnionSolid(nam, trap, sub, Transform3D(RotationX(M_PI/2.), Position(0,0,displacement))).ptr();
}
+ stringstream params;
+ params << x1 << " " << x2 << " " << y1 << " " << y2 << " " << z << " "
+ << r << " " << char(atMinusZ ? '1' : '0') << " ";
+ solid->GetBoolNode()->GetRightMatrix()->SetTitle(params.str().c_str());
_assign(solid,"",PSEUDOTRAP_TAG, true);
}
/// Helper function to create poly hedron
-void PolyhedraRegular::make(const std::string& nam, int nsides, double rmin, double rmax,
+void PolyhedraRegular::make(const string& nam, int nsides, double rmin, double rmax,
double zpos, double zneg, double start, double delta) {
if (rmin < 0e0 || rmin > rmax)
throw runtime_error("dd4hep: PolyhedraRegular: Illegal argument rmin:<" + _toString(rmin) + "> is invalid!");
@@ -1301,7 +1559,7 @@ void PolyhedraRegular::make(const std::string& nam, int nsides, double rmin, dou
}
/// Helper function to create poly hedron
-void Polyhedra::make(const std::string& nam, int nsides, double start, double delta,
+void Polyhedra::make(const string& nam, int nsides, double start, double delta,
const vector<double>& z, const vector<double>& rmin, const vector<double>& rmax) {
vector<double> temp;
if ( rmin.size() != z.size() || rmax.size() != z.size() ) {
@@ -1324,7 +1582,7 @@ void Polyhedra::make(const std::string& nam, int nsides, double start, double de
}
/// Helper function to create the polyhedron
-void ExtrudedPolygon::make(const std::string& nam,
+void ExtrudedPolygon::make(const string& nam,
const vector<double>& pt_x,
const vector<double>& pt_y,
const vector<double>& sec_z,
@@ -1341,7 +1599,7 @@ void ExtrudedPolygon::make(const std::string& nam,
}
/// Creator method
-void EightPointSolid::make(const std::string& nam, double dz, const double* vtx) {
+void EightPointSolid::make(const string& nam, double dz, const double* vtx) {
_assign(new TGeoArb8(nam.c_str(), dz, (double*)vtx), "", "EightPointSolid", true);
}
@@ -1376,31 +1634,31 @@ SubtractionSolid::SubtractionSolid(const Solid& shape1, const Solid& shape2, con
}
/// Constructor to be used when creating a new object. Position is identity, Rotation is the identity rotation
-SubtractionSolid::SubtractionSolid(const std::string& nam, const Solid& shape1, const Solid& shape2) {
+SubtractionSolid::SubtractionSolid(const string& nam, const Solid& shape1, const Solid& shape2) {
TGeoSubtraction* sub = new TGeoSubtraction(shape1, shape2, detail::matrix::_identity(), detail::matrix::_identity());
_assign(new TGeoCompositeShape(nam.c_str(), sub), "", SUBTRACTION_TAG, true);
}
/// Constructor to be used when creating a new object. Placement by a generic transformation within the mother
-SubtractionSolid::SubtractionSolid(const std::string& nam, const Solid& shape1, const Solid& shape2, const Transform3D& trans) {
+SubtractionSolid::SubtractionSolid(const string& nam, const Solid& shape1, const Solid& shape2, const Transform3D& trans) {
TGeoSubtraction* sub = new TGeoSubtraction(shape1, shape2, detail::matrix::_identity(), detail::matrix::_transform(trans));
_assign(new TGeoCompositeShape(nam.c_str(), sub), "", SUBTRACTION_TAG, true);
}
/// Constructor to be used when creating a new object. Rotation is the identity rotation
-SubtractionSolid::SubtractionSolid(const std::string& nam, const Solid& shape1, const Solid& shape2, const Position& pos) {
+SubtractionSolid::SubtractionSolid(const string& nam, const Solid& shape1, const Solid& shape2, const Position& pos) {
TGeoSubtraction* sub = new TGeoSubtraction(shape1, shape2, detail::matrix::_identity(), detail::matrix::_translation(pos));
_assign(new TGeoCompositeShape(nam.c_str(), sub), "", SUBTRACTION_TAG, true);
}
/// Constructor to be used when creating a new object
-SubtractionSolid::SubtractionSolid(const std::string& nam, const Solid& shape1, const Solid& shape2, const RotationZYX& rot) {
+SubtractionSolid::SubtractionSolid(const string& nam, const Solid& shape1, const Solid& shape2, const RotationZYX& rot) {
TGeoSubtraction* sub = new TGeoSubtraction(shape1, shape2, detail::matrix::_identity(), detail::matrix::_rotationZYX(rot));
_assign(new TGeoCompositeShape(nam.c_str(), sub), "", SUBTRACTION_TAG, true);
}
/// Constructor to be used when creating a new object
-SubtractionSolid::SubtractionSolid(const std::string& nam, const Solid& shape1, const Solid& shape2, const Rotation3D& rot) {
+SubtractionSolid::SubtractionSolid(const string& nam, const Solid& shape1, const Solid& shape2, const Rotation3D& rot) {
TGeoSubtraction* sub = new TGeoSubtraction(shape1, shape2, detail::matrix::_identity(), detail::matrix::_rotation3D(rot));
_assign(new TGeoCompositeShape(nam.c_str(), sub), "", SUBTRACTION_TAG, true);
}
@@ -1436,31 +1694,31 @@ UnionSolid::UnionSolid(const Solid& shape1, const Solid& shape2, const Rotation3
}
/// Constructor to be used when creating a new object. Position is identity, Rotation is identity rotation
-UnionSolid::UnionSolid(const std::string& nam, const Solid& shape1, const Solid& shape2) {
+UnionSolid::UnionSolid(const string& nam, const Solid& shape1, const Solid& shape2) {
TGeoUnion* uni = new TGeoUnion(shape1, shape2, detail::matrix::_identity(), detail::matrix::_identity());
_assign(new TGeoCompositeShape(nam.c_str(), uni), "", UNION_TAG, true);
}
/// Constructor to be used when creating a new object. Placement by a generic transformation within the mother
-UnionSolid::UnionSolid(const std::string& nam, const Solid& shape1, const Solid& shape2, const Transform3D& trans) {
+UnionSolid::UnionSolid(const string& nam, const Solid& shape1, const Solid& shape2, const Transform3D& trans) {
TGeoUnion* uni = new TGeoUnion(shape1, shape2, detail::matrix::_identity(), detail::matrix::_transform(trans));
_assign(new TGeoCompositeShape(nam.c_str(), uni), "", UNION_TAG, true);
}
/// Constructor to be used when creating a new object. Rotation is identity rotation
-UnionSolid::UnionSolid(const std::string& nam, const Solid& shape1, const Solid& shape2, const Position& pos) {
+UnionSolid::UnionSolid(const string& nam, const Solid& shape1, const Solid& shape2, const Position& pos) {
TGeoUnion* uni = new TGeoUnion(shape1, shape2, detail::matrix::_identity(), detail::matrix::_translation(pos));
_assign(new TGeoCompositeShape(nam.c_str(), uni), "", UNION_TAG, true);
}
/// Constructor to be used when creating a new object
-UnionSolid::UnionSolid(const std::string& nam, const Solid& shape1, const Solid& shape2, const RotationZYX& rot) {
+UnionSolid::UnionSolid(const string& nam, const Solid& shape1, const Solid& shape2, const RotationZYX& rot) {
TGeoUnion *uni = new TGeoUnion(shape1, shape2, detail::matrix::_identity(), detail::matrix::_rotationZYX(rot));
_assign(new TGeoCompositeShape(nam.c_str(), uni), "", UNION_TAG, true);
}
/// Constructor to be used when creating a new object
-UnionSolid::UnionSolid(const std::string& nam, const Solid& shape1, const Solid& shape2, const Rotation3D& rot) {
+UnionSolid::UnionSolid(const string& nam, const Solid& shape1, const Solid& shape2, const Rotation3D& rot) {
TGeoUnion *uni = new TGeoUnion(shape1, shape2, detail::matrix::_identity(), detail::matrix::_rotation3D(rot));
_assign(new TGeoCompositeShape(nam.c_str(), uni), "", UNION_TAG, true);
}
@@ -1496,31 +1754,31 @@ IntersectionSolid::IntersectionSolid(const Solid& shape1, const Solid& shape2, c
}
/// Constructor to be used when creating a new object. Position is identity, Rotation is identity rotation
-IntersectionSolid::IntersectionSolid(const std::string& nam, const Solid& shape1, const Solid& shape2) {
+IntersectionSolid::IntersectionSolid(const string& nam, const Solid& shape1, const Solid& shape2) {
TGeoIntersection* inter = new TGeoIntersection(shape1, shape2, detail::matrix::_identity(), detail::matrix::_identity());
_assign(new TGeoCompositeShape(nam.c_str(), inter), "", INTERSECTION_TAG, true);
}
/// Constructor to be used when creating a new object. Placement by a generic transformation within the mother
-IntersectionSolid::IntersectionSolid(const std::string& nam, const Solid& shape1, const Solid& shape2, const Transform3D& trans) {
+IntersectionSolid::IntersectionSolid(const string& nam, const Solid& shape1, const Solid& shape2, const Transform3D& trans) {
TGeoIntersection* inter = new TGeoIntersection(shape1, shape2, detail::matrix::_identity(), detail::matrix::_transform(trans));
_assign(new TGeoCompositeShape(nam.c_str(), inter), "", INTERSECTION_TAG, true);
}
/// Constructor to be used when creating a new object. Position is identity.
-IntersectionSolid::IntersectionSolid(const std::string& nam, const Solid& shape1, const Solid& shape2, const Position& pos) {
+IntersectionSolid::IntersectionSolid(const string& nam, const Solid& shape1, const Solid& shape2, const Position& pos) {
TGeoIntersection* inter = new TGeoIntersection(shape1, shape2, detail::matrix::_identity(), detail::matrix::_translation(pos));
_assign(new TGeoCompositeShape(nam.c_str(), inter), "", INTERSECTION_TAG, true);
}
/// Constructor to be used when creating a new object
-IntersectionSolid::IntersectionSolid(const std::string& nam, const Solid& shape1, const Solid& shape2, const RotationZYX& rot) {
+IntersectionSolid::IntersectionSolid(const string& nam, const Solid& shape1, const Solid& shape2, const RotationZYX& rot) {
TGeoIntersection* inter = new TGeoIntersection(shape1, shape2, detail::matrix::_identity(), detail::matrix::_rotationZYX(rot));
_assign(new TGeoCompositeShape(nam.c_str(), inter), "", INTERSECTION_TAG, true);
}
/// Constructor to be used when creating a new object
-IntersectionSolid::IntersectionSolid(const std::string& nam, const Solid& shape1, const Solid& shape2, const Rotation3D& rot) {
+IntersectionSolid::IntersectionSolid(const string& nam, const Solid& shape1, const Solid& shape2, const Rotation3D& rot) {
TGeoIntersection* inter = new TGeoIntersection(shape1, shape2, detail::matrix::_identity(), detail::matrix::_rotation3D(rot));
_assign(new TGeoCompositeShape(nam.c_str(), inter), "", INTERSECTION_TAG, true);
}
diff --git a/DDCore/src/plugins/ShapePlugins.cpp b/DDCore/src/plugins/ShapePlugins.cpp
index 5290b8fe7..d464456e3 100644
--- a/DDCore/src/plugins/ShapePlugins.cpp
+++ b/DDCore/src/plugins/ShapePlugins.cpp
@@ -660,7 +660,15 @@ void* shape_mesh_verifier(Detector& description, int argc, char** argv) {
for (Int_t ipv=0, npv=v->GetNdaughters(); ipv < npv; ipv++) {
PlacedVolume place = v->GetNode(ipv);
Solid solid = place.volume().solid();
- if ( solid->IsA() != TGeoCompositeShape::Class() ) {
+ if ( instanceOf<TruncatedTube>(solid) ) {
+ auto params = solid.dimensions();
+ solid.setDimensions(params);
+ }
+ else if ( instanceOf<PseudoTrap>(solid) ) {
+ auto params = solid.dimensions();
+ solid.setDimensions(params);
+ }
+ else if ( solid->IsA() != TGeoCompositeShape::Class() ) {
auto params = solid.dimensions();
solid.setDimensions(params);
}
diff --git a/DDG4/include/DDG4/Geant4SensDetAction.h b/DDG4/include/DDG4/Geant4SensDetAction.h
index 2febc51bd..86c5b5e03 100644
--- a/DDG4/include/DDG4/Geant4SensDetAction.h
+++ b/DDG4/include/DDG4/Geant4SensDetAction.h
@@ -126,7 +126,7 @@ namespace dd4hep {
int m_hitCreationMode = 0;
#if defined(G__ROOT) || defined(__CLING__) || defined(__ROOTCLING__)
/// Reference to the detector description object
- Detector* m_detDesc;
+ Detector* m_detDesc {0};
#else
/// Reference to the detector description object
Detector& m_detDesc;
diff --git a/DDG4/plugins/Geant4EventReaderHepMC.cpp b/DDG4/plugins/Geant4EventReaderHepMC.cpp
index f3842e9a3..ad0afc49d 100644
--- a/DDG4/plugins/Geant4EventReaderHepMC.cpp
+++ b/DDG4/plugins/Geant4EventReaderHepMC.cpp
@@ -309,11 +309,11 @@ void HepMC::fix_particles(EventStream& info) {
for(id=v->out.begin(); id!=v->out.end();++id) {
EventStream::Particles::iterator ipp = parts.find(*id);
Geant4Particle* dau = ipp != parts.end() ? (*ipp).second : 0;
- if ( !dau ) {
+ if ( !dau )
cout << "ERROR: Invalid daughter particle: " << *id << endl;
- }
+ else
+ dau->parents.insert(p->id);
p->daughters.insert(*id);
- dau->parents.insert(p->id);
}
}
}
@@ -486,6 +486,10 @@ int HepMC::read_vertex(EventStream &info, istream& is, istringstream & input)
}
input >> id >> dummy >> v->x >> v->y >> v->z >> v->time
>> num_orphans_in >> num_particles_out >> weights_size;
+ if( !input ) {
+ delete v;
+ return 0;
+ }
#if defined(DD4HEP_DEBUG_HEP_MC_VERTEX)
if ( id == DD4HEP_DEBUG_HEP_MC_VERTEX ) {
printout(ALWAYS,"HepMC","++ Created Vertex ID=%d",id);
diff --git a/DDG4/plugins/Geant4GDMLWriteAction.cpp b/DDG4/plugins/Geant4GDMLWriteAction.cpp
index 2c79155d4..24e66ce20 100644
--- a/DDG4/plugins/Geant4GDMLWriteAction.cpp
+++ b/DDG4/plugins/Geant4GDMLWriteAction.cpp
@@ -134,11 +134,12 @@ void Geant4GDMLWriteAction::writeGDML() {
error("+++ No GDML file name given. Please set the output file (property Output)");
return;
}
- if ( 0 == ::stat(fname.c_str(), &buff) && !m_overWrite ) {
- error("+++ GDML file elready exists. Please set another output file (property Output)");
- return;
- }
- if ( 0 == ::stat(fname.c_str(), &buff) && m_overWrite ) {
+ if ( 0 == ::stat(fname.c_str(), &buff) ) {
+ if ( !m_overWrite ) {
+ error("+++ GDML file %s elready exists. Please set another output file (property Output)",
+ fname.c_str());
+ return;
+ }
warning("+++ GDML file %s already exists. Overwriting existing file.", fname.c_str());
::unlink(fname.c_str());
}
--
GitLab