ClientTests/src/LheD_tracker_SiVertexBarrel_geo.cpp

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 //--------------------------------------------------------------------
0011 //
0012 //  Author     : M.Frank
0013 //  Modified   : E.Pilicer ( tube + elliptical vertex layout )
0014 //
0015 //====================================================================
0016 #include "DD4hep/DetFactoryHelper.h"
0017 #include "DD4hep/Printout.h"
0018
0019 using namespace std;
0020 using namespace dd4hep;
0021 using namespace dd4hep::detail;
0022
0023 /*************************************************************
0024  function to calculate path in a given theta
0025 **************************************************************/
0026 static double computeDpt( double ra, double rb, double theta ) {
0027     double dpt_sin = std::pow(ra * std::sin(theta), 2.0);
0028     double dpt_cos = std::pow(rb * std::cos(theta), 2.0);
0029     double dp = std::sqrt(dpt_sin + dpt_cos);
0030     return dp;
0031 }
0032
0033 static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector sens)  {
0034   xml_det_t   x_det     = e;
0035   Material    air       = description.air();
0036   int         det_id    = x_det.id();
0037   string      det_name  = x_det.nameStr();
0038   DetElement  sdet       (det_name,det_id);
0039   Assembly    assembly   (det_name+"_assembly");
0040   map<string, Volume>    volumes;
0041   PlacedVolume pv;
0042
0043   double  mod_width = 0.0;
0044
0045   sens.setType("tracker");
0046   for(xml_coll_t mi(x_det,_U(module)); mi; ++mi)  {
0047     xml_comp_t x_mod  = mi;
0048     xml_comp_t m_env  = x_mod.child(_U(module_envelope));
0049     string     m_nam  = x_mod.nameStr();
0050
0051     mod_width = m_env.width();
0052
0053     Volume     m_vol(det_name+"_"+m_nam,Box(m_env.width()/2,m_env.length()/2,m_env.thickness()/2),air);
0054     int        ncomponents = 0, sensor_number = 1;
0055
0056     if ( volumes.find(m_nam) != volumes.end() )   {
0057       printout(ERROR,"SiTrackerBarrel","Logics error in building modules.");
0058       throw runtime_error("Logics error in building modules.");
0059     }
0060     volumes[m_nam] = m_vol;
0061
0062     for(xml_coll_t ci(x_mod,_U(module_component)); ci; ++ci, ++ncomponents)  {
0063       xml_comp_t x_comp = ci;
0064       xml_comp_t x_pos  = x_comp.position(false);
0065       xml_comp_t x_rot  = x_comp.rotation(false);
0066       string     c_nam  = det_name+"_"+m_nam+_toString(ncomponents,"_component%d");
0067       Box        c_box(x_comp.width()/2,x_comp.length()/2,x_comp.thickness()/2);
0068       Volume     c_vol(c_nam,c_box,description.material(x_comp.materialStr()));
0069       if ( x_pos && x_rot ) {
0070         Position    c_pos(x_pos.x(0),x_pos.y(0),x_pos.z(0));
0071         RotationZYX c_rot(x_rot.z(0),x_rot.y(0),x_rot.x(0));
0072         pv = m_vol.placeVolume(c_vol, Transform3D(c_rot,c_pos));
0073       }
0074       else if ( x_rot ) {
0075         pv = m_vol.placeVolume(c_vol,RotationZYX(x_rot.z(0),x_rot.y(0),x_rot.x(0)));
0076       }
0077       else if ( x_pos ) {
0078         pv = m_vol.placeVolume(c_vol,Position(x_pos.x(0),x_pos.y(0),x_pos.z(0)));
0079       }
0080       else {
0081         pv = m_vol.placeVolume(c_vol);
0082       }
0083       c_vol.setRegion(description, x_comp.regionStr());
0084       c_vol.setLimitSet(description, x_comp.limitsStr());
0085       c_vol.setVisAttributes(description, x_comp.visStr());
0086       if ( x_comp.isSensitive() ) {
0087         pv.addPhysVolID(_U(sensor),sensor_number++);
0088         c_vol.setSensitiveDetector(sens);
0089       }
0090     }
0091     m_vol.setVisAttributes(description.visAttributes(x_mod.visStr()));
0092   }
0093
0094   for(xml_coll_t li(x_det,_U(layer)); li; ++li)  {
0095     xml_comp_t x_layer  = li;
0096     xml_comp_t x_barrel = x_layer.child(_U(barrel_envelope));
0097     xml_comp_t x_layout = x_layer.child(_U(rphi_layout));
0098     xml_comp_t z_layout = x_layer.child(_U(z_layout));      // Get the <z_layout> element.
0099     int        lay_id   = x_layer.id();
0100     string     m_nam    = x_layer.moduleStr();
0101     Volume     m_env    = volumes[m_nam];
0102     string     lay_nam  = det_name+"_"+m_nam+_toString(x_layer.id(),"_layer%d");
0103
0104     double     phi0     = x_layout.phi0();              // Starting phi of first module.
0105     double     phi_tilt = x_layout.phi_tilt();          // Phi tilt of a module.
0106     double     rc       = x_layout.rc();                // Radius of the module center.
0107     double     rphi_dr  = x_layout.dr();                // The delta radius of every other module.
0108     double     phic     = phi0;                         // Phi of the module center.
0109     double     z0       = z_layout.z0();                // Z position of first module in phi.
0110     double     nz       = z_layout.nz();                // Number of modules to place in z.
0111     double     z_dr     = z_layout.dr();                // Radial displacement parameter, of every other module.
0112
0113     // not necessary, calculated here
0114     //int        nphi     = x_layout.nphi();             // Number of modules in phi.
0115     //double     phi_incr = (2*M_PI) / nphi;             // Phi increment for one module.
0116
0117     // Z increment for module placement along Z axis.
0118     // Adjust for z0 at center of module rather than
0119     // the end of cylindrical envelope.
0120     double z_incr   = nz > 1 ? (2.0 * z0) / (nz - 1) : 0.0;
0121     // Starting z for module placement along Z axis.
0122     double module_z = -z0;
0123     int module = 1;
0124
0125     // multiplication factor for ellipse major radius
0126     double c0;
0127     if (x_layer.id() <=1 ) c0 = 3.5;
0128     else if (x_layer.id() == 2 ) c0 = 2.3;
0129     else if (x_layer.id() == 3 ) c0 = 2;
0130     else c0 = 1.8;
0131
0132     // create an envelope
0133     double env_rmin = x_barrel.inner_r();    // inner radius for envelope
0134     double env_rmax = x_barrel.outer_r();    // outer radius for envelope
0135     double env_z    = x_barrel.z_length();   // length of envelope
0136
0137     EllipticalTube envElTubeOut(c0*env_rmax,env_rmax, env_z);
0138     EllipticalTube envElTubeInn(c0*env_rmin,env_rmin, env_z+0.01);
0139     SubtractionSolid envElTube(envElTubeOut,envElTubeInn);
0140
0141     Tube envTube1(env_rmin, env_rmax, env_z+0.01, 3*M_PI/2, M_PI/2);
0142     UnionSolid lay_tub1(envElTube,envTube1);
0143
0144     Tube envTube2(env_rmax, c0*env_rmax, env_z+0.01, 3*M_PI/2, M_PI/2);
0145     SubtractionSolid lay_tub(lay_tub1,envTube2);
0146     Volume     lay_vol   (lay_nam,lay_tub,air);         // Create the layer envelope volume.
0147
0148     /*************************************************************
0149      FIRST-HALF
0150      semi-elliptical, from 90 to 270 degrees
0151     **************************************************************/
0152
0153     double a            = c0 * rc;    // (mm) ellipse major radius
0154     double b            = rc;         // (mm) ellipse minor radius
0155     double theta        = 0.;
0156     double thetaMin     = M_PI / 2.;
0157     double thetaMax     = M_PI * 3. / 2.;
0158     double deltaTheta   = 0.0001;
0159     double numIntegrals = (thetaMax-thetaMin) / deltaTheta;
0160     double cell         = 0.;
0161
0162     // integrate over the elipse to get the circumference
0163     for( int i=0; i < numIntegrals; i++ ) {
0164         if (i==0) theta = thetaMin;
0165         else theta += deltaTheta;
0166         cell += computeDpt( a, b, theta);
0167     }
0168
0169     // number of modules along semi-ellipse path
0170     int    n_ell     = (cell * deltaTheta / mod_width); // - 1 ;
0171     int    nextPoint = 0;
0172     double run       = 0.;
0173     theta            = 0.;
0174
0175     for( int i=0; i < numIntegrals; i++ ) {
0176         if (i==0) theta = thetaMin;
0177         else theta += deltaTheta;
0178         double subIntegral = n_ell*run/cell;
0179         if( (int) subIntegral >= nextPoint ) {
0180             double x = a * std::cos(theta);
0181             double y = b * std::sin(theta);
0182             // Loop over the number of modules in z.
0183             for (int j = 0; j < nz; j++)  {
0184             // Module PhysicalVolume.
0185             Transform3D tr(RotationZYX(0,M_PI/2-theta,-M_PI/2),Position(x,y,module_z));
0186             pv = lay_vol.placeVolume(m_env,tr);
0187             pv.addPhysVolID("module", module++);
0188             // Add z increment to get next z placement pos.
0189             module_z += z_incr;
0190             }
0191             nextPoint++;
0192         }
0193         run += computeDpt(a, b, theta);
0194     }
0195
0196     /*************************************************************
0197      SECOND-HALF
0198      semi-circular, from 270 to 90 degrees
0199     **************************************************************/
0200
0201     // circle circumference
0202     double c_circ       = 2 * M_PI * rc;
0203     int    n_circ       = c_circ / mod_width / 2;  // number of modules for semi-circle
0204     double phi_incr     = (M_PI) / n_circ;         // Phi increment for one module along semi-circle
0205
0206     // Loop over the number of modules in phi.
0207
0208     for (int ii = 0; ii < n_circ; ii++)    {
0209
0210       double dx = z_dr * std::cos(phic + phi_tilt);        // Delta x of module position.
0211       double dy = z_dr * std::sin(phic + phi_tilt);        // Delta y of module position.
0212       double  x = rc * std::cos(phic);                     // Basic x module position.
0213       double  y = rc * std::sin(phic);                     // Basic y module position.
0214
0215       // Loop over the number of modules in z.
0216       for (int j = 0; j < nz; j++)  {
0217          // Module PhysicalVolume.
0218          // Transform3D tr(RotationZYX(0,-((M_PI/2)-phic-phi_tilt),M_PI/2),Position(x,y,module_z));
0219          //NOTE (Nikiforos, 26/08 Rotations needed to be fixed so that component1 (silicon) is on the outside
0220          Transform3D tr(RotationZYX(0,((M_PI/2)-phic-phi_tilt),-M_PI/2),Position(x,y,module_z));
0221          pv = lay_vol.placeVolume(m_env,tr);
0222          pv.addPhysVolID("module", module++);
0223          // Adjust the x and y coordinates of the module.
0224          x += dx;
0225          y += dy;
0226          // Flip sign of x and y adjustments.
0227          dx *= -1;
0228          dy *= -1;
0229          // Add z increment to get next z placement pos.
0230          module_z += z_incr;
0231       }
0232       phic      += phi_incr; // Increment the phi placement of module.
0233       rc        += rphi_dr;  // Increment the center radius according to dr parameter.
0234       rphi_dr   *= -1;       // Flip sign of dr parameter.
0235       module_z   = -z0;      // Reset the Z placement parameter for module.
0236     }
0237
0238     // Create the PhysicalVolume for the layer.
0239     assembly.setVisAttributes(description.invisible());
0240     pv = assembly.placeVolume(lay_vol);     // Place layer in mother
0241     pv.addPhysVolID("layer", lay_id);       // Set the layer ID.
0242     DetElement m_elt(sdet,lay_nam,lay_id);
0243     m_elt.setAttributes(description,lay_vol,x_layer.regionStr(),x_layer.limitsStr(),x_layer.visStr());
0244     m_elt.setPlacement(pv);
0245   }
0246
0247   pv = description.pickMotherVolume(sdet).placeVolume(assembly);
0248   pv.addPhysVolID("system", det_id);      // Set the subdetector system ID.
0249   pv.addPhysVolID("barrel", 0);           // Flag this as a barrel subdetector.
0250   sdet.setPlacement(pv);
0251   return sdet;
0252 }
0253
0254 DECLARE_DETELEMENT(LheD_tracker_SiVertexBarrel,create_detector)