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