File indexing completed on 2026-06-01 07:50:55
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015 #include <DD4hep/DetFactoryHelper.h>
0016 #include <DD4hep/DD4hepUnits.h>
0017 #include <TGeoArb8.h>
0018
0019
0020 #include <iomanip>
0021
0022 using namespace std;
0023 using namespace dd4hep;
0024 using namespace dd4hep::detail;
0025
0026 static void printCoordinates(char , const char* , double [8][2]) {
0027 }
0028 #if 0
0029 static void printCoordinates(char sector_type, const char* volume, double v[8][2]) {
0030 cout << "Sector type " << sector_type << " Volume " << volume << endl;
0031 cout << "Coordinate 1: x:" << setw(10) << v[0][0] << " y:" << setw(10) << v[0][1] << endl;
0032 cout << "Coordinate 2: x:" << setw(10) << v[1][0] << " y:" << setw(10) << v[1][1] << endl;
0033 cout << "Coordinate 3: x:" << setw(10) << v[2][0] << " y:" << setw(10) << v[2][1] << endl;
0034 cout << "Coordinate 4: x:" << setw(10) << v[3][0] << " y:" << setw(10) << v[3][1] << endl;
0035 }
0036 #endif
0037
0038 static Ref_t create_element(Detector& description, xml_h e, SensitiveDetector sens_det) {
0039 xml_det_t x_det = e;
0040 string name = x_det.nameStr();
0041 DetElement sdet(name,x_det.id());
0042 Volume motherVol = description.pickMotherVolume(sdet);
0043 xml_comp_t x_envelope = x_det.child(_Unicode(envelope));
0044 xml_comp_t x_sectors = x_det.child(_Unicode(sectors));
0045 xml_comp_t x_inner = x_det.child(_Unicode(inner_wall));
0046 xml_comp_t x_outer = x_det.child(_Unicode(outer_wall));
0047 xml_comp_t x_gas = x_det.child(_Unicode(gas));
0048 xml_comp_t x_cathode = x_det.child(_Unicode(cathode));
0049 xml_comp_t x_sens = x_det.child(_Unicode(sensitive), false);
0050
0051 PlacedVolume pv;
0052
0053 struct cylinder_t { double inner, outer, zhalf; };
0054 cylinder_t env = { x_envelope.inner_r(), x_envelope.outer_r(), x_envelope.zhalf() };
0055 cylinder_t inner_wall = { x_inner.inner_r(), x_inner.inner_r()+x_inner.thickness(), env.zhalf };
0056 cylinder_t outer_wall = { x_outer.outer_r()-x_outer.thickness(), x_outer.outer_r(), env.zhalf };
0057 cylinder_t gas = { inner_wall.outer+5*dd4hep::cm, outer_wall.inner-5*dd4hep::cm, x_gas.zhalf()-5*dd4hep::cm };
0058
0059
0060 string sens_typ = x_sens.ptr() ? x_sens.typeStr("tracker") : string("tracker");
0061 sens_det.setType(sens_typ);
0062 cout << "Detector: " << name << " id: " << x_det.id() << " Sensitive type: " << sens_typ << endl;
0063
0064
0065
0066
0067 Assembly envVol(name+"_envelope");
0068
0069
0070
0071
0072
0073 pv = motherVol.placeVolume(envVol);
0074 pv.addPhysVolID("system",x_det.id());
0075 sdet.setPlacement(pv);
0076
0077
0078 Tube innerTub(inner_wall.inner, inner_wall.outer, inner_wall.zhalf);
0079 Volume innerVol(name+"_inner",innerTub,description.material(x_inner.materialStr()));
0080 innerVol.setVisAttributes(description.visAttributes(x_inner.visStr()));
0081 envVol.placeVolume(innerVol);
0082
0083
0084 Tube outerTub(outer_wall.inner, outer_wall.outer, outer_wall.zhalf);
0085 Volume outerVol(name+"_outer",outerTub,description.material(x_outer.materialStr()));
0086 outerVol.setVisAttributes(description.visAttributes(x_outer.visStr()));
0087 envVol.placeVolume(outerVol);
0088
0089
0090 Material gasMat = description.material(x_gas.materialStr());
0091 Tube gasTub(gas.inner,gas.outer,gas.zhalf);
0092 Volume gasVol(name+"_chamber",gasTub,gasMat);
0093 gasVol.setVisAttributes(description.visAttributes(x_gas.visStr()));
0094 gasVol.setSensitiveDetector(sens_det);
0095 pv = envVol.placeVolume(gasVol);
0096 pv.addPhysVolID("layer", 10);
0097
0098
0099 Tube hvTub(gas.inner,gas.outer,x_cathode.thickness()/2);
0100 Volume hvVol(name+"_cathode",hvTub,description.material(x_cathode.materialStr()));
0101 hvVol.setVisAttributes(description.visAttributes(x_cathode.visStr()));
0102 envVol.placeVolume(hvVol);
0103
0104
0105 Tube endCapPlane(env.inner,env.outer,0.0001);
0106 Volume endCapPlaneVol(name+"_plane",endCapPlane,description.material("Copper"));
0107 endCapPlaneVol.setVisAttributes(description.visAttributes("TPCEndcapVis"));
0108
0109 double endcap_thickness = 0;
0110 for(xml_coll_t c(x_sectors,_Unicode(sector)); c; ++c) {
0111
0112 double thickness = endCapPlane->GetDz();
0113 for(xml_coll_t l(x_sectors.child(_Unicode(layers)),_Unicode(layer)); l; ++l) {
0114 xml_comp_t x_layer = l;
0115 thickness += 2*x_layer.thickness();
0116 }
0117 if ( thickness > endcap_thickness ) endcap_thickness = thickness;
0118 }
0119
0120 Tube endCapTub(env.inner,env.outer,endcap_thickness/2);
0121 Volume endCapAVol(name+"_end_A",endCapTub,description.material("Copper"));
0122 endCapAVol.setVisAttributes(description.invisible());
0123
0124 endCapAVol.placeVolume(endCapPlaneVol,Position(0,0,-endcap_thickness/2));
0125 pv = envVol.placeVolume(endCapAVol,Position(0,0,env.zhalf+endcap_thickness/2));
0126 pv.addPhysVolID("side",1);
0127 DetElement detA(sdet,name+"_SideA",1);
0128 detA.setPlacement(pv);
0129
0130 Volume endCapBVol(name+"_end_B",endCapTub,description.material("Copper"));
0131 endCapBVol.setVisAttributes(description.invisible());
0132
0133 endCapBVol.placeVolume(endCapPlaneVol,Position(0,0,-endcap_thickness/2));
0134
0135 Transform3D trEndB(RotationY(M_PI),Position(0,0,-env.zhalf-endcap_thickness/2));
0136 pv = envVol.placeVolume(endCapBVol, trEndB);
0137 pv.addPhysVolID("side",2);
0138 DetElement detB(sdet,name+"_SideB",2);
0139 detB.setPlacement(pv);
0140
0141 envVol.setVisAttributes(description.visAttributes(x_det.visStr()));
0142
0143 int sector_count = 0;
0144 for(xml_coll_t c(x_sectors,_Unicode(sector)); c; ++c) {
0145 xml_comp_t x_sector = c;
0146 const char sector_type = x_sector.typeStr()[0];
0147 const double rmin0 = x_sector.rmin();
0148 const double rmax0 = x_sector.rmax();
0149
0150
0151
0152 const int num_sectors = sector_type == 'K' ? 6 : 12;
0153 const double shift = sector_type == 'K' ? 0 : M_PI/num_sectors;
0154 const double dphi = 2*M_PI/double(num_sectors);
0155 string sector_vis = x_sector.visStr();
0156 Solid tm;
0157 double z_start = 0.0;
0158 Assembly sector(name+"_sector_"+sector_type);
0159 int i_layer = 0;
0160 for( xml_coll_t l(x_sectors.child(_Unicode(layers)),_Unicode(layer)); l; ++l, ++i_layer) {
0161 xml_comp_t x_layer = l;
0162 double layer_thickness = x_layer.thickness();
0163 string layer_vis = x_layer.visStr();
0164 string layer_mat = x_layer.materialStr();
0165 double gap_half = 1;
0166 double rmin = rmin0;
0167 double rmax = rmax0;
0168
0169 if ( layer_vis.empty() ) layer_vis = sector_vis;
0170
0171 if ( sector_type == 'K' ) {
0172 double angle = M_PI/12.0;
0173 double angle1 = std::tan(angle);
0174 double v[8][2];
0175 v[0][0] = rmin;
0176 v[0][1] = 0;
0177 v[1][0] = rmin;
0178 v[1][1] = rmin*angle1;
0179 v[2][0] = rmax;
0180 v[2][1] = rmax*angle1;
0181 v[3][0] = rmax;
0182 v[3][1] = 0;
0183 ::memcpy(&v[4][0],&v[0][0],8*sizeof(double));
0184 EightPointSolid upper(layer_thickness/2,&v[0][0]);
0185
0186 v[0][0] = rmin;
0187 v[0][1] = 0;
0188 v[1][0] = rmax;
0189 v[1][1] = 0;
0190 v[2][0] = rmax;
0191 v[2][1] = -rmax*angle1;
0192 v[3][0] = rmin;
0193 v[3][1] = -rmin*angle1;
0194 ::memcpy(&v[4][0],&v[0][0],8*sizeof(double));
0195 EightPointSolid lower(layer_thickness/2,&v[0][0]);
0196
0197 v[0][0] = rmin;
0198 v[0][1] = gap_half;
0199 v[1][0] = rmin;
0200 v[1][1] = rmin*angle1;
0201 v[2][0] = rmax;
0202 v[2][1] = rmax*angle1;
0203 v[3][0] = rmax;
0204 v[3][1] = gap_half;
0205 ::memcpy(&v[4][0],&v[0][0],8*sizeof(double));
0206 EightPointSolid top(layer_thickness/2,&v[0][0]);
0207
0208 v[0][0] = rmin;
0209 v[0][1] = -gap_half;
0210 v[1][0] = rmax;
0211 v[1][1] = -gap_half;
0212 v[2][0] = rmax;
0213 v[2][1] = -rmax*angle1;
0214 v[3][0] = rmin;
0215 v[3][1] = -rmin*angle1;
0216 ::memcpy(&v[4][0],&v[0][0],8*sizeof(double));
0217 EightPointSolid bottom(layer_thickness/2,&v[0][0]);
0218
0219 UnionSolid u(UnionSolid(upper,lower),top,Rotation3D(RotationZ(2*(angle))));
0220 tm = UnionSolid(u,bottom,Rotation3D(RotationZ(-2*(angle))));
0221 }
0222 else {
0223 double overlap = sector_type=='W' ? 20 : -20;
0224 double angle = M_PI/12.0;
0225 double angle1 = std::tan(angle);
0226 double v[8][2], dr = 0;
0227
0228 if ( sector_type == 'W' ) {
0229 rmax += overlap*std::tan(angle/2);
0230 dr = overlap*std::tan(angle/2);
0231 }
0232
0233 v[0][0] = rmin;
0234 v[0][1] = -rmin*angle1+gap_half;
0235 v[1][0] = rmin;
0236 v[1][1] = rmin*angle1-gap_half;
0237 v[2][0] = rmax;
0238 v[2][1] = rmax*angle1-gap_half;
0239 v[3][0] = rmax;
0240 v[3][1] = -rmax*angle1+gap_half;
0241 ::memcpy(&v[4][0],&v[0][0],8*sizeof(double));
0242 printCoordinates(sector_type,"t2:",v);
0243 EightPointSolid sectorSolid(layer_thickness/2,&v[0][0]);
0244
0245 if ( sector_type=='W' ) {
0246 v[0][0] = (rmax+rmin)/2-dr;
0247 v[0][1] = (rmax+rmin)/2*angle1-gap_half;
0248 v[1][0] = rmax+0.0001;
0249 v[1][1] = rmax*angle1-gap_half;
0250 v[2][0] = rmax+0.0001;
0251 v[2][1] = (rmax-overlap)*angle1-gap_half;
0252 v[3][0] = (rmax+rmin)/2-dr;
0253 v[3][1] = (rmax+rmin-overlap)/2*angle1-gap_half;
0254 }
0255 else {
0256 v[0][0] = (rmax+rmin)/2-dr;
0257 v[0][1] = (rmax+rmin)/2*angle1-gap_half;
0258 v[3][0] = rmax+0.0001;
0259 v[3][1] = rmax*angle1-gap_half;
0260 v[2][0] = rmax+0.0001;
0261 v[2][1] = (rmax-overlap)*angle1-gap_half;
0262 v[1][0] = (rmax+rmin)/2-dr;
0263 v[1][1] = (rmax+rmin-overlap)/2*angle1-gap_half;
0264 }
0265 printCoordinates(sector_type,"upper_boolean:",v);
0266 ::memcpy(&v[4][0],&v[0][0],8*sizeof(double));
0267 EightPointSolid upper_boolean((sector_type == 'W' ? 1.1 : 1.0) * layer_thickness/2,&v[0][0]);
0268
0269 if ( sector_type=='W' ) {
0270 v[0][0] = (rmax+rmin)/2-dr;
0271 v[0][1] = -((rmax+rmin)/2*angle1-gap_half);
0272 v[1][0] = (rmax+rmin)/2-dr;
0273 v[1][1] = -((rmax+rmin-overlap)/2*angle1-gap_half);
0274 v[2][0] = rmax+0.0001;
0275 v[2][1] = -((rmax-overlap)*angle1-gap_half);
0276 v[3][0] = rmax+0.0001;
0277 v[3][1] = -(rmax*angle1-gap_half);
0278 }
0279 else {
0280 v[0][0] = (rmax+rmin)/2-dr;
0281 v[0][1] = -((rmax+rmin)/2*angle1-gap_half);
0282 v[3][0] = (rmax+rmin)/2-dr;
0283 v[3][1] = -((rmax+rmin-overlap)/2*angle1-gap_half);
0284 v[2][0] = rmax+0.0001;
0285 v[2][1] = -((rmax-overlap)*angle1-gap_half);
0286 v[1][0] = rmax+0.0001;
0287 v[1][1] = -(rmax*angle1-gap_half);
0288 }
0289 printCoordinates(sector_type,"lower_boolean:",v);
0290 ::memcpy(&v[4][0],&v[0][0],8*sizeof(double));
0291
0292
0293
0294 EightPointSolid lower_boolean((sector_type == 'W' ? 1.1 : 1.0) * layer_thickness/2,&v[0][0]);
0295 if ( sector_type == 'W' )
0296 tm = SubtractionSolid(SubtractionSolid(sectorSolid,upper_boolean),lower_boolean);
0297 else
0298 tm = UnionSolid(UnionSolid(sectorSolid,upper_boolean),lower_boolean);
0299 }
0300 Material lmat = description.material(layer_mat);
0301 Volume secVol( name+"_sector_"+sector_type+_toString(i_layer, "_layer%d"), tm, lmat );
0302 secVol.setVisAttributes( description.visAttributes(layer_vis) );
0303 if ( x_layer.isSensitive() ) {
0304 secVol.setSensitiveDetector(sens_det);
0305 pv.addPhysVolID("layer", i_layer+1);
0306 }
0307
0308 sector.placeVolume(secVol,Position(0,0,z_start+layer_thickness/2));
0309 z_start += layer_thickness;
0310 }
0311 for(int i=0; i<num_sectors;++i) {
0312 int j = i + (sector_type=='W' ? 1 : 0);
0313 double phi = dphi*j+shift + (sector_type=='K' ? 0 : M_PI/12.0);
0314 if ( sector_type == 'K' || (i%2)==0 ) {
0315 Transform3D trA( RotationZYX(phi,0,0), Position(0,0,0.00001) );
0316 Transform3D trB( RotationZYX(phi,0,0), Position(0,0,0.00001) );
0317
0318 pv = endCapAVol.placeVolume( sector, trA );
0319 pv.addPhysVolID( "type", sector_type=='K' ? 1 : sector_type=='M' ? 2 : 3 );
0320 pv.addPhysVolID( "sector", j );
0321 DetElement sectorA( detA, detA.name()+_toString(sector_count,"_sector%02d"),1 );
0322 sectorA.setPlacement(pv);
0323
0324 pv = endCapBVol.placeVolume( sector, trB );
0325 pv.addPhysVolID( "type", sector_type=='K' ? 1 : sector_type=='M' ? 2 : 3 );
0326 pv.addPhysVolID( "sector",j );
0327 DetElement sectorB(detB, detB.name()+_toString(sector_count,"_sector%02d"), 1 );
0328 sectorB.setPlacement(pv);
0329 cout << "Placed " << sector_type << " sector at phi=" << phi << endl;
0330 ++sector_count;
0331 }
0332 }
0333 }
0334 return sdet;
0335 }
0336
0337 DECLARE_SUBDETECTOR(AlephTPC,create_element)