Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-02-25 09:20:18

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   
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   // TPC sensitive detector
0059   sens_det.setType("tracker");
0060 
0061   // the TPC mother volume
0062   //Tube    envTub(env.inner,env.outer,env.zhalf);
0063   //Volume  envVol(name+"_envelope",envTub,description.air());
0064   Assembly  envVol(name+"_envelope");
0065   //envVol.setVisAttributes(description.visAttributes(x_envelope.visStr()));
0066   //envVol.setVisAttributes(description.invisible());
0067   //envVol->SetVisibility(kFALSE);
0068   //envVol->SetVisDaughters(kTRUE);
0069 
0070   pv = motherVol.placeVolume(envVol);
0071   pv.addPhysVolID("system",x_det.id());
0072   sdet.setPlacement(pv);
0073 
0074   // TPC Al inner shield 
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   // TPC outer shield 
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   // TPC gas chamber envelope
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   // TPC HV plane
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   // TPC Endcap plane to see sectors and misalignments betters
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     //xml_comp_t x_sector = c;
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   //endCapAVol.setVisAttributes(description.visAttributes(x_outer.visStr()));
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   //endCapBVol.setVisAttributes(description.visAttributes(x_outer.visStr()));
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   //envVol.setVisAttributes(description.invisible());
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     //const int  padrows = x_sector.attr<int>(_Unicode(padrows));
0146     //const int  trgrows = x_sector.attr<int>(_Unicode(trgrows));
0147     //const int  nwires  = x_sector.attr<int>(_Unicode(numwires));
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         // For W sectors make the subtraction solid slightly thicker to ensure everything is cut off
0289         // and no left-overs from numerical precision are left.
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)