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