Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-30 09:17:48

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 includes
0015 #include "DD4hep/Printout.h"
0016 #include "DD4hep/MatrixHelpers.h"
0017 #include "DD4hep/DetFactoryHelper.h"
0018 #include "XML/VolumeBuilder.h"
0019 #include "XML/Utilities.h"
0020 #include <cmath>
0021 
0022 using namespace std;
0023 using namespace dd4hep;
0024 using namespace dd4hep::detail;
0025 
0026 namespace   {
0027   const char* col(int idx)  {
0028     const char*  cols[5] = {"VisibleRed","VisibleBlue","VisibleGreen","VisibleYellow","VisibleCyan"};
0029     return cols[idx%(sizeof(cols)/sizeof(cols[0]))];
0030   }
0031   Rotation3D makeRotReflect(double thetaX, double phiX, double thetaY, double phiY, double thetaZ, double phiZ) {
0032     // define 3 unit std::vectors forming the new left-handed axes
0033     Position x(cos(phiX) * sin(thetaX), sin(phiX) * sin(thetaX), cos(thetaX));
0034     Position y(cos(phiY) * sin(thetaY), sin(phiY) * sin(thetaY), cos(thetaY));
0035     Position z(cos(phiZ) * sin(thetaZ), sin(phiZ) * sin(thetaZ), cos(thetaZ));
0036 
0037     constexpr double tol = 1.0e-3;       // Geant4 compatible
0038     double check = (x.Cross(y)).Dot(z);  // in case of a LEFT-handed orthogonal system this must be -1
0039     if (std::abs(1. + check) > tol) {
0040       except("NestedBoxReflection", "+++ FAILED to construct Rotation is not LEFT-handed!");
0041     }
0042     printout(INFO, "NestedBoxReflection", "+++ Constructed LEFT-handed reflection rotation.");
0043     Rotation3D rotation(x.x(), y.x(), z.x(), x.y(), y.y(), z.y(), x.z(), y.z(), z.z());
0044     return rotation;
0045   }
0046 
0047   TGeoCombiTrans* createPlacement(const Rotation3D& iRot, const Position& iTrans) {
0048     TGeoRotation r;
0049     double elements[9];
0050     iRot.GetComponents(elements);
0051     r.SetMatrix(elements);
0052     return new TGeoCombiTrans(TGeoTranslation(iTrans.x(), iTrans.y(), iTrans.z()), r);
0053   }
0054 
0055   //static Transform3D transform_reflect(xml_h element)   {
0056   TGeoCombiTrans* transform_reflect(xml_h element)   {
0057     xml_dim_t xrot(element.child(_U(rotation)));
0058     xml_dim_t xpos(element.child(_U(position)));
0059     double thetaX = xrot.attr<double>(Unicode("thetaX"));
0060     double phiX   = xrot.attr<double>(Unicode("phiX"));
0061     double thetaY = xrot.attr<double>(Unicode("thetaY"));
0062     double phiY   = xrot.attr<double>(Unicode("phiY"));
0063     double thetaZ = xrot.attr<double>(Unicode("thetaZ"));
0064     double phiZ   = xrot.attr<double>(Unicode("phiZ"));
0065     printout(INFO, "NestedBoxReflection",
0066          "+++ Adding reflection rotation \"%s\": "
0067          "(theta/phi)[rad] X: %6.3f %6.3f Y: %6.3f %6.3f Z: %6.3f %6.3f",
0068          element.attr<string>(_U(name)).c_str(), thetaX, phiX, thetaY, phiY, thetaZ, phiZ);
0069     Rotation3D rot = makeRotReflect(thetaX, phiX, thetaY, phiY, thetaZ, phiZ);
0070     Position   pos = Position(xpos.x(),xpos.y(),xpos.z());
0071     //  return Transform3D(rot, pos);
0072     return createPlacement(rot, pos);
0073   }
0074 }
0075 
0076 namespace   {
0077   struct TubeReflection : public xml::tools::VolumeBuilder  {
0078     Material   silicon;
0079     Material   iron;
0080     ///
0081     using VolumeBuilder::VolumeBuilder;
0082     ///
0083     void place_quadrants(int /* level */, Volume vol)    {
0084       Tube         tube = vol.solid();
0085       double       rmin = tube.rMin(), rmax = tube.rMax(), dz = tube.dZ();
0086       Tube         tsol1(rmin, rmax, dz*0.05, 0, M_PI*2.);
0087       Tube         tsol2(rmin, rmax, dz*0.15, 0, M_PI/2.);
0088       Volume       tvol0("Tube0", tsol1, silicon);
0089       Volume       tvol1("Tube1", tsol2, silicon);
0090       Volume       tvol2("Tube2", tsol2, silicon);
0091       Volume       tvol3("Tube3", tsol2, silicon);
0092       Volume       tvol4("Tube4", tsol2, silicon);
0093       TGeoHMatrix* mat;
0094       Transform3D  tr;
0095       
0096       double dz1 = dz-2.0*tsol2.dZ()-tsol1.dZ();
0097       tvol0.setVisAttributes(description.visAttributes("VisibleCyan"));
0098       tvol1.setVisAttributes(description.visAttributes("VisibleRed"));
0099       tvol2.setVisAttributes(description.visAttributes("VisibleBlue"));
0100       tvol3.setVisAttributes(description.visAttributes("VisibleYellow"));
0101       tvol4.setVisAttributes(description.visAttributes("VisibleGreen"));
0102 
0103       mat = detail::matrix::_transform(Transform3D(Position(0, 0, dz1)));
0104       vol.placeVolume(tvol0, mat);
0105 
0106       mat = detail::matrix::_transform(Transform3D(Position(0, 0,-dz1)));
0107       mat->ReflectZ(kTRUE, kTRUE);
0108       vol.placeVolume(tvol0, mat);
0109 
0110       double dz2 = dz-tsol2.dZ();
0111       // OK
0112       mat = detail::matrix::_transform(Transform3D(Position(0, 0, dz2)));
0113       vol.placeVolume(tvol1, mat);
0114       // OK
0115       mat = detail::matrix::_transform(Transform3D(Position(0, 0, dz2)));
0116       mat->RotateZ(1.0*M_PI/2.0/dd4hep::degree);
0117       vol.placeVolume(tvol2, mat);
0118       // OK
0119       mat = detail::matrix::_transform(Transform3D(Position(0, 0, dz2)));
0120       mat->RotateZ(2.0*M_PI/2.0/dd4hep::degree);
0121       vol.placeVolume(tvol3, mat);
0122       // OK
0123       mat = detail::matrix::_transform(Transform3D(Position(0, 0, dz2)));
0124       mat->RotateZ(3.0*M_PI/2.0/dd4hep::degree);
0125       vol.placeVolume(tvol4, mat);
0126 
0127 
0128       /** Now eflect the quadrants to the other endcap:  */      
0129       // OK
0130       mat = detail::matrix::_transform(Transform3D(Position(0, 0,-dz2)));
0131       mat->ReflectY(kTRUE,kTRUE);
0132       //mat->RotateZ(1.0*M_PI/2.0/dd4hep::degree);
0133       vol.placeVolume(tvol1, mat);
0134 
0135       // OK
0136       mat = detail::matrix::_transform(Transform3D(Position(0, 0,-dz2)));
0137       mat->RotateZ(1.0*M_PI/2.0/dd4hep::degree);
0138       mat->ReflectY(kTRUE,kTRUE);
0139       //mat->RotateZ(-1.0*M_PI/2.0/dd4hep::degree);
0140       vol.placeVolume(tvol2, mat);
0141 
0142       // OK
0143       mat = detail::matrix::_transform(Transform3D(Position(0, 0,-dz2)));
0144       mat->RotateZ(2.0*M_PI/2.0/dd4hep::degree);
0145       mat->ReflectY(kTRUE,kTRUE);
0146       //mat->RotateZ(1.0*M_PI/2.0/dd4hep::degree);
0147       vol.placeVolume(tvol3, mat);
0148 
0149       // OK
0150       mat = detail::matrix::_transform(Transform3D(Position(0, 0,-dz2)));
0151       mat->RotateZ(3.0*M_PI/2.0/dd4hep::degree);
0152       mat->ReflectY(kTRUE,kTRUE);
0153       //mat->RotateZ(-1.0*M_PI/2.0/dd4hep::degree);
0154       vol.placeVolume(tvol4, mat);
0155     }
0156     /// 
0157     PlacedVolume place(Volume mother, Volume vol, xml_elt_t e, int level, int copyNo, char reflection)   {
0158       Position    pos;
0159       RotationZYX rot;
0160       xml_dim_t   xpos = xml_dim_t(e).child(_U(position), false);
0161       xml_dim_t   xrot = xml_dim_t(e).child(_U(rotation), false);
0162       if ( !xpos.ptr() ) xpos = e;
0163       if ( xpos.ptr()  ) pos = Position(xpos.x(0),xpos.y(0),xpos.z(0));
0164       if ( xrot.ptr()  ) rot = RotationZYX(xrot.x(0),xrot.y(0),xrot.z(0));
0165 
0166       TGeoHMatrix* mat = detail::matrix::_transform(pos, rot);
0167       switch(reflection)  {
0168       case 'X':
0169     mat->ReflectX(kTRUE, kTRUE);
0170     break;
0171       case 'Y':
0172     mat->ReflectY(kTRUE, kTRUE);
0173     break;
0174       case 'Z':
0175       default:
0176     mat->ReflectZ(kTRUE, kTRUE);
0177     break;
0178       }
0179       PlacedVolume pv = mother.placeVolume(vol, mat);
0180       pv.addPhysVolID(_toString(level,"lvl%d"), copyNo);
0181       return pv;
0182     }
0183     ///
0184     DetElement create()    {
0185       xml_dim_t    x_tube(x_det.dimensions()); 
0186       double       r      = x_tube.r();
0187       double       dz     = x_tube.dz();
0188       int          levels = x_tube.level(2);
0189       Volume       tube_vol;
0190       Volume       v_det;
0191       PlacedVolume pv;
0192 
0193       sensitive.setType("tracker");
0194       silicon = description.material("Si");
0195       iron    = description.material("Iron");
0196       tube_vol = Volume("tube", Tube(0,r,dz,0,M_PI*2.0), silicon);  
0197       if ( levels != 0 && x_det.hasChild(_U(assembly)) )
0198     v_det = Assembly("envelope");
0199       else
0200     v_det = Volume("envelope",Tube(4.5*r,4.5*r,4.5*r),description.air());
0201 
0202       int cnt = 0;
0203       Transform3D tr = xml::createTransformation(x_tube);
0204       tube_vol.setAttributes(description, x_det.regionStr(), x_det.limitsStr(), "VisibleGrey");
0205       place_quadrants(levels-1, tube_vol);
0206       pv = v_det.placeVolume(tube_vol, tr);
0207       pv.addPhysVolID(_toString(levels,"lvl%d"), ++cnt);
0208 
0209       for(xml_coll_t c(x_det,_U(reflect_x)); c; ++c)
0210     place(v_det, tube_vol, c, levels, ++cnt, 'X');
0211       for(xml_coll_t c(x_det,_U(reflect_y)); c; ++c)
0212     place(v_det, tube_vol, c, levels, ++cnt, 'Y');
0213       for(xml_coll_t c(x_det,_U(reflect_z)); c; ++c)
0214     place(v_det, tube_vol, c, levels, ++cnt, 'Z');
0215 
0216       if ( x_det.hasChild(_U(reflect)) )   {
0217     Volume reflect_vol = tube_vol;
0218     for(xml_coll_t c(x_det,_U(reflect)); c; ++c)   {
0219       TGeoCombiTrans* reflect_tr = transform_reflect(c);
0220       pv = v_det.placeVolume(reflect_vol.ptr(), reflect_tr);
0221       pv.addPhysVolID(_toString(levels,"lvl%d"), ++cnt);
0222     }
0223       }
0224       // Place the calo inside the world
0225       placeDetector(v_det, x_det).addPhysVolID("system",x_det.id());
0226       return detector;
0227     }
0228   };
0229 
0230   static Ref_t create_refl_tube(Detector& description, xml_dim_t x_det, SensitiveDetector sens)  {
0231     TubeReflection builder(description, x_det, sens);
0232     return builder.create();
0233   }
0234 }
0235 DECLARE_DETELEMENT(DD4hep_Test_TubeReflection,create_refl_tube)
0236 
0237 
0238 // makes sure that the RotationMatrix built is
0239 // LEFT-handed coordinate system (i.e. reflected)
0240 namespace   {
0241   struct NestedBoxReflection : public xml::tools::VolumeBuilder  {
0242     using VolumeBuilder::VolumeBuilder;
0243     Material   silicon;
0244     Material   iron;
0245     
0246     void place_box(int level, int ident, const char* vis, Volume par, Solid sol, TGeoHMatrix* mtx);
0247     void place_box(int level, int ident, const char* vis, Volume par, Solid sol, const Transform3D& tr)   {
0248       TGeoHMatrix* mtx = detail::matrix::_transform(tr);
0249       place_box(level, ident, vis, par, sol, mtx);
0250     }
0251     void place_box(int level, int ident, const char* vis, Volume par, Solid sol, Position pos)    {
0252       place_box(level, ident, vis, par, sol, Transform3D(pos));
0253     }
0254     void place_boxes(int level, Volume vol);
0255     PlacedVolume place(Volume mother, Volume vol, xml_elt_t e, int level, int copyNo, char reflection);
0256     DetElement   create();
0257   };
0258   /// 
0259   void NestedBoxReflection::place_boxes(int level, Volume vol)    {
0260     if ( level >= 0 )   {
0261       Box          box  = vol.solid();
0262       double       line = 0.015, bx = box.x(), by = box.y(), bz = box.z();
0263       Box          mbox(bx*0.2, by*0.2, bz*0.2);
0264 
0265       printout(INFO,"NestedBoxReflection","+++ Placing boxes at level %d",level);
0266       place_box(level, 1, col(0+level), vol, mbox,                           Position(0,0,0));
0267       place_box(level, 2, col(1+level), vol, Box(bx*0.2,  by*0.2,  bz*0.4),  Position(0.8*bx,0,0));
0268       place_box(level, 3, col(2+level), vol, mbox,                           Position(0,0.8*by,0));
0269       place_box(level, 4, col(3+level), vol, mbox,                           Position(0,0,0.8*bz));
0270       auto mtx = make_unique<TGeoHMatrix>(TGeoTranslation(0,0,-0.9*bz));
0271       mtx->ReflectZ(kTRUE,kTRUE);
0272       place_box(level, 5, col(4+level), vol, Box(bx*1.0,  by*1.0,  bz*0.1),  mtx.release());
0273       place_box(level, 0, col(1+level), vol, Box(bx*0.2,  by*line, bz*line), Position(0.4*bx,0,0));
0274       place_box(level, 0, col(2+level), vol, Box(bx*line, by*0.2,  bz*line), Position(0,0.4*by,0));
0275       place_box(level, 0, col(3+level), vol, Box(bx*line, by*line, bz*0.2),  Position(0,0,0.4*bz));
0276     }
0277   }
0278   
0279   /// 
0280   void NestedBoxReflection::place_box(int level, int ident, const char* vis, Volume par, Solid sol, TGeoHMatrix* tr)   {
0281     PlacedVolume pv;
0282     bool         sens  = level == 2 || level == 0;
0283     string       idnam = _toString(ident,"box%d");
0284     Volume       vol   = Volume(idnam, sol, sens ? silicon : iron);
0285     vol.setRegion(par.region());
0286     vol.setLimitSet(par.limitSet());
0287     vol.setVisAttributes(description, vis);
0288     pv = par.placeVolume(vol, tr);
0289     if ( ident > 0 )   {
0290       string vidn  = _toString(level,"lvl%d");
0291       if ( sens )   {
0292     vol.setSensitiveDetector(sensitive);
0293       }
0294       pv.addPhysVolID(vidn, ident);
0295       string n = par.name()+string("/")+vol.name();
0296       printout(INFO,"NestedBoxReflection","++ Level: %3d Volume:%-24s Sensitive:%s Color:%s vid:%s=%d",
0297            level, n.c_str(), yes_no(sens), vis, vidn.c_str(), ident);
0298       place_boxes(level-1, vol);
0299     }
0300   }
0301 
0302   /// 
0303   PlacedVolume NestedBoxReflection::place(Volume mother, Volume vol, xml_elt_t e, int level, int copyNo, char reflection)   {
0304     Position    pos;
0305     RotationZYX rot;
0306     xml_dim_t   xpos = xml_dim_t(e).child(_U(position), false);
0307     xml_dim_t   xrot = xml_dim_t(e).child(_U(rotation), false);
0308     if ( !xpos.ptr() ) xpos = e;
0309     if ( xpos.ptr()  ) pos = Position(xpos.x(0),xpos.y(0),xpos.z(0));
0310     if ( xrot.ptr()  ) rot = RotationZYX(xrot.x(0),xrot.y(0),xrot.z(0));
0311 
0312     TGeoHMatrix* mat = detail::matrix::_transform(pos, rot);
0313     switch(reflection)  {
0314     case 'X':
0315       mat->ReflectX(kTRUE, kTRUE);
0316       break;
0317     case 'Y':
0318       mat->ReflectY(kTRUE, kTRUE);
0319       break;
0320     case 'Z':
0321     default:
0322       mat->ReflectZ(kTRUE, kTRUE);
0323       break;
0324     }
0325     PlacedVolume pv = mother.placeVolume(vol, mat);
0326     pv.addPhysVolID(_toString(level,"lvl%d"), copyNo);
0327     return pv;
0328   }
0329   
0330   /// 
0331   DetElement NestedBoxReflection::create()    {
0332     xml_dim_t    x_box(x_det.dimensions()); 
0333     int          levels = x_box.level(2);
0334     double       bx = x_box.x();
0335     double       by = x_box.y();
0336     double       bz = x_box.z();
0337     Volume       v_det;
0338     Box          box_solid(bx,by,bz);
0339     Volume       box_vol("nested_box",box_solid,description.air());  
0340     PlacedVolume pv;
0341 
0342     sensitive.setType("tracker");
0343     iron    = description.material("Iron");
0344     silicon = description.material("Si");
0345     if ( levels != 0 && x_det.hasChild(_U(assembly)) )
0346       v_det = Assembly("envelope");
0347     else
0348       v_det = Volume("envelope",Box(4.5*bx,4.5*by,4.5*bz),description.air());
0349 
0350     if ( levels == 0 )   {
0351       v_det.setSensitiveDetector(sensitive);
0352     }
0353     else if ( levels == 1 )   {
0354       box_vol.setSensitiveDetector(sensitive);
0355       box_vol.setAttributes(description,x_det.regionStr(),x_det.limitsStr(),"VisibleGrey");
0356     }
0357     else  {
0358       int cnt = 0;
0359       Transform3D tr = xml::createTransformation(x_box);
0360       box_vol.setAttributes(description,x_det.regionStr(),x_det.limitsStr(),"VisibleGrey");
0361       place_boxes(levels-1, box_vol);
0362       pv = v_det.placeVolume(box_vol, tr);
0363       pv.addPhysVolID(_toString(levels,"lvl%d"), ++cnt);
0364 
0365       for(xml_coll_t c(x_det,_U(reflect_x)); c; ++c)
0366     place(v_det, box_vol, c, levels, ++cnt, 'X');
0367       for(xml_coll_t c(x_det,_U(reflect_y)); c; ++c)
0368     place(v_det, box_vol, c, levels, ++cnt, 'Y');
0369       for(xml_coll_t c(x_det,_U(reflect_z)); c; ++c)
0370     place(v_det, box_vol, c, levels, ++cnt, 'Z');
0371 
0372       if ( x_det.hasChild(_U(reflect)) )   {
0373     Volume reflect_vol = box_vol;
0374     for(xml_coll_t c(x_det,_U(reflect)); c; ++c)   {
0375       TGeoCombiTrans* reflect_tr = transform_reflect(c);
0376       pv = v_det.placeVolume(reflect_vol.ptr(), reflect_tr);
0377       pv.addPhysVolID(_toString(levels,"lvl%d"), ++cnt);
0378     }
0379       }
0380     }
0381     // Place the calo inside the world
0382     placeDetector(v_det, x_det).addPhysVolID("system",x_det.id());
0383     return detector;
0384   }
0385 }
0386 
0387 static Ref_t create_nested_box(Detector& description, xml_dim_t x_det, SensitiveDetector sens)  {
0388   NestedBoxReflection builder(description, x_det, sens);
0389   return builder.create();
0390 }
0391 
0392 DECLARE_DETELEMENT(DD4hep_Test_NestedBoxReflection,create_nested_box)