Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2026-06-01 07:50:55

0001 //==========================================================================
0002 //  AIDA Detector description implementation 
0003 //--------------------------------------------------------------------------
0004 // Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
0005 // All rights reserved.
0006 //
0007 // For the licensing terms see $DD4hepINSTALL/LICENSE.
0008 // For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
0009 //
0010 // Author     : M.Frank
0011 //
0012 //==========================================================================
0013 
0014 // Framework include files
0015 #include <DD4hep/DetFactoryHelper.h>
0016 #include <DD4hep/DD4hepUnits.h>
0017 #include <TGeoArb8.h>
0018 
0019 // C/C++ include files
0020 #include <iomanip>
0021 
0022 using namespace std;
0023 using namespace dd4hep;
0024 using namespace dd4hep::detail;
0025 
0026 static void printCoordinates(char /* sector_type */, const char* /* volume */, 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   // TPC sensitive detector
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   // the TPC mother volume
0065   //Tube    envTub(env.inner,env.outer,env.zhalf);
0066   //Volume  envVol(name+"_envelope",envTub,description.air());
0067   Assembly  envVol(name+"_envelope");
0068   //envVol.setVisAttributes(description.visAttributes(x_envelope.visStr()));
0069   //envVol.setVisAttributes(description.invisible());
0070   //envVol->SetVisibility(kFALSE);
0071   //envVol->SetVisDaughters(kTRUE);
0072 
0073   pv = motherVol.placeVolume(envVol);
0074   pv.addPhysVolID("system",x_det.id());
0075   sdet.setPlacement(pv);
0076 
0077   // TPC Al inner shield 
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   // TPC outer shield 
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   // TPC gas chamber envelope
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   // TPC HV plane
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   // TPC Endcap plane to see sectors and misalignments betters
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     //xml_comp_t x_sector = c;
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   //endCapAVol.setVisAttributes(description.visAttributes(x_outer.visStr()));
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   //endCapBVol.setVisAttributes(description.visAttributes(x_outer.visStr()));
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     //const int  padrows = x_sector.attr<int>(_Unicode(padrows));
0150     //const int  trgrows = x_sector.attr<int>(_Unicode(trgrows));
0151     //const int  nwires  = x_sector.attr<int>(_Unicode(numwires));
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         // For W sectors make the subtraction solid slightly thicker to ensure everything is cut off
0293         // and no left-overs from numerical precision are left.
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)