athena/src/DIRC_geo.cpp

0001 #include "DD4hep/DetFactoryHelper.h"
0002 #include "DD4hep/OpticalSurfaces.h"
0003 #include "DD4hep/Printout.h"
0004 #include "DDRec/DetectorData.h"
0005 #include "DDRec/Surface.h"
0006 #include <XML/Helper.h>
0007
0008 //////////////////////////////////
0009 // Central Barrel DIRC
0010 //////////////////////////////////
0011
0012 using namespace std;
0013 using namespace dd4hep;
0014
0015 // Fixed Trap constructor. This function is a workaround of this bug:
0016 // https://github.com/AIDASoft/DD4hep/issues/850
0017 // Should be used instead of dd4hep::Trap(pName, pZ, pY, pX, pLTX) constructor
0018 dd4hep::Trap MakeTrap( const std::string& pName, double pZ, double pY, double pX, double pLTX );
0019
0020 static Ref_t createDetector(Detector& desc, xml_h e, SensitiveDetector sens)
0021 {
0022   xml_det_t xml_det   = e;
0023   string    det_name = xml_det.nameStr();
0024   int       det_id   = xml_det.id();
0025
0026   // Main detector xml element
0027   xml_dim_t dirc_dim   = xml_det.dimensions();
0028   xml_dim_t dirc_pos   = xml_det.position();
0029   xml_dim_t dirc_rot   = xml_det.rotation();
0030   double    det_rin   = dirc_dim.rmin();
0031   double    det_rout  = dirc_dim.rmax();
0032   double    SizeZ = dirc_dim.length();
0033
0034   // DEBUG
0035   // double mirror_r1 = x_det.attr<double>(_Unicode(r1));
0036
0037   // DIRC box:
0038   xml_comp_t xml_box_module = xml_det.child(_U(module));
0039
0040   Material Vacuum = desc.material("Vacuum");
0041   Material air = desc.material("AirOptical");
0042   Material quartz = desc.material("Quartz");
0043   Material epotek = desc.material("Epotek");
0044   Material nlak33a = desc.material("Nlak33a");
0045   auto& bar_material = quartz;
0046   auto mirror_material = desc.material("Aluminum");  // mirror material
0047
0048   Tube     det_geo(det_rin, det_rout, SizeZ / 2., 0., 360.0 * deg);
0049   //Volume   det_volume("DIRC", det_geo, Vacuum);
0050   Assembly   det_volume("DIRC");
0051   det_volume.setVisAttributes(desc.visAttributes(xml_det.visStr()));
0052
0053   DetElement   det(det_name, det_id);
0054   Volume       mother_vol = desc.pickMotherVolume(det);
0055
0056   Transform3D  tr(RotationZYX(0, dirc_rot.theta(), 0.0), Position(0.0, 0.0, dirc_pos.z()));
0057   PlacedVolume det_plvol = mother_vol.placeVolume(det_volume, tr);
0058
0059   det_plvol.addPhysVolID("system", det_id);
0060   det.setPlacement(det_plvol);
0061
0062   // Parts Dimentions
0063
0064   int fLensId = 6; // focusing system
0065                    // 0    no lens
0066                    // 1    spherical lens
0067                    // 3    3-layer spherical lens
0068                    // 6    3-layer cylindrical lens
0069                    // 10   ideal lens (thickness = 0, ideal focusing)
0070
0071   int fGeomType = 0;  // Full DIRC - 0, 1 only one plate
0072   int fRunType = 0;   // 0, 10 - simulation, 1, 5 - lookup table, 2,3,4 - reconstruction
0073
0074
0075   double fPrizm[4];
0076   fPrizm[0] = 360 * mm;
0077   fPrizm[1] = 300 * mm;
0078   fPrizm[3] = 50 * mm;
0079   fPrizm[2] = fPrizm[3] + 300 * tan(32 * deg) * mm;
0080   double fBarsGap = 0.15 * mm;
0081
0082   std::cout << "DIRC: fPrizm[2] " << fPrizm[2] << std::endl;
0083
0084   double fdTilt = 80 * deg;
0085   double fPrizmT[6];
0086   fPrizmT[0] = 390 * mm;
0087   fPrizmT[1] = (400 - 290 * cos(fdTilt)) * mm; //
0088   fPrizmT[2] = 290 * sin(fdTilt) * mm;       // hight
0089   fPrizmT[3] = 50 * mm;                      // face
0090   fPrizmT[4] = 290 * mm;
0091   fPrizmT[5] = 290 * cos(fdTilt)* mm;
0092
0093   double fMirror[3];
0094   fMirror[0] = 20 * mm;
0095   fMirror[1] = fPrizm[0];
0096   fMirror[2] = 1 * mm;
0097   //  fPrizm[0] = 170; fPrizm[1] = 300; fPrizm[2] = 50+300*tan(45*deg); fPrizm[3] = 50;
0098
0099 //  double fBar[3];
0100 //  fBar[0] = 17 * mm;
0101 //  fBar[1] = (fPrizm[0] - (fNBar - 1) * fBarsGap) / fNBar;
0102 //  fBar[2] = 1050 * mm; // 4200; //4200
0103
0104   double fMcpTotal[3];
0105   double fMcpActive[3];
0106   fMcpTotal[0] = fMcpTotal[1] = 53 + 4;
0107   fMcpTotal[2] = 1*mm;
0108   fMcpActive[0] = fMcpActive[1] = 53;
0109   fMcpActive[2] = 1*mm;
0110
0111   double fLens[4];
0112   fLens[0] = fLens[1] = 40 * mm;
0113   fLens[2] = 10 * mm;
0114   double fRadius = (det_rin + det_rout)/2;
0115
0116   double fBoxWidth = fPrizm[0];
0117
0118   double fFd[3];
0119   fFd[0] = fBoxWidth;
0120   fFd[1] = fPrizm[2];
0121   fFd[2] = 1 * mm;
0122
0123   fLens[0] = fPrizm[3];
0124   fLens[1] = fPrizm[0];
0125   fLens[2] = 12 * mm;
0126
0127   // Getting box XML
0128   const int fNBoxes = xml_box_module.repeat();
0129   const double box_width = xml_box_module.width();
0130   const double box_height = xml_box_module.height();
0131   const double box_length = xml_box_module.length() + 550*mm;
0132
0133
0134   // The DIRC
0135   Assembly dirc_module("DIRCModule");
0136
0137   //Volume lDirc("lDirc", gDirc, air);
0138   dirc_module.setVisAttributes(desc.visAttributes(xml_box_module.visStr()));
0139
0140
0141   // FD... whatever F and D is
0142   xml_comp_t xml_fd = xml_box_module.child(_Unicode(fd));
0143   Box gFd("gFd", xml_fd.height()/2, xml_fd.width()/2, xml_fd.thickness()/2);
0144   Volume lFd ("lFd", gFd, desc.material(xml_fd.materialStr()));
0145   lFd.setVisAttributes(desc.visAttributes(xml_fd.visStr()));
0146   //lFd.setSensitiveDetector(sens);
0147
0148
0149   // The Bar
0150   xml_comp_t xml_bar = xml_box_module.child(_Unicode(bar));
0151   double bar_height = xml_bar.height();
0152   double bar_width = xml_bar.width();
0153   double bar_length = xml_bar.length();
0154   Box gBar("gBar", bar_height/2, bar_width/2, bar_length/2);
0155   Volume lBar("lBar", gBar, desc.material(xml_bar.materialStr()));
0156   lBar.setVisAttributes(desc.visAttributes(xml_bar.visStr()));
0157
0158   // Glue
0159   xml_comp_t xml_glue = xml_box_module.child(_Unicode(glue));
0160   double glue_thickness = xml_glue.thickness();  // 0.05 * mm;
0161   Box gGlue("gGlue", bar_height/2, bar_width/2, glue_thickness/2);
0162   Volume lGlue("lGlue", gGlue, desc.material(xml_glue.materialStr()));
0163   lGlue.setVisAttributes(desc.visAttributes(xml_glue.visStr()));
0164
0165   sens.setType("tracker");
0166   lBar.setSensitiveDetector(sens);
0167
0168   int bars_repeat_z = 4;    // TODO parametrize!
0169   double bar_assm_length = (bar_length + glue_thickness) * bars_repeat_z;
0170   int fNBar = xml_bar.repeat();
0171   double bar_gap = xml_bar.gap();
0172   for (int y_index = 0; y_index < fNBar; y_index++) {
0173     double shift_y = y_index * (bar_width + bar_gap) - 0.5 * box_width + 0.5 * bar_width;
0174     for (int z_index = 0; z_index < bars_repeat_z; z_index++) {
0175       double z = -0.5 * bar_assm_length + 0.5 * bar_length + (bar_length + glue_thickness) * z_index;
0176       auto placed_bar = dirc_module.placeVolume(lBar, Position(0, shift_y, z));
0177       dirc_module.placeVolume(lGlue, Position(0, shift_y, z + 0.5 * (bar_length + glue_thickness)));
0178       placed_bar.addPhysVolID("section", z_index);
0179       placed_bar.addPhysVolID("bar", y_index);
0180     }
0181   }
0182
0183
0184   // The Mirror
0185   xml_comp_t xml_mirror = xml_box_module.child(_Unicode(mirror));
0186   Box gMirror("gMirror", xml_mirror.height()/2, xml_mirror.width()/2, xml_mirror.thickness()/2);
0187   Volume lMirror("lMirror", gMirror, desc.material(xml_mirror.materialStr()));
0188   dirc_module.placeVolume(lMirror, Position(0, 0, -0.5 * (bar_assm_length - xml_mirror.thickness())));
0189   lMirror.setVisAttributes(desc.visAttributes(xml_mirror.visStr()));
0190
0191   // The mirror optical surface
0192   OpticalSurfaceManager surfMgr = desc.surfaceManager();
0193   auto surf = surfMgr.opticalSurface("MirrorOpticalSurface");
0194   SkinSurface skin(desc, det, Form("dirc_mirror_optical_surface"), surf, lMirror);
0195   skin.isValid();
0196
0197   // LENS
0198   // Lens volumes
0199   Volume lLens1;
0200   Volume lLens2;
0201   Volume lLens3;
0202
0203   double lensMinThikness = 2.0 * mm;
0204   double layer12 = lensMinThikness * 2;
0205
0206   // r1 = (r1==0)? 27.45: r1;
0207   // r2 = (r2==0)? 20.02: r2;
0208
0209   double r1 = 33 * mm;
0210   double r2 = 24 * mm;
0211   double shight = 25 * mm;
0212
0213   Position zTrans1(0, 0, -r1 - fLens[2] / 2. + r1 - sqrt(r1 * r1 - shight / 2. * shight / 2.) + lensMinThikness);
0214   Position zTrans2(0, 0, -r2 - fLens[2] / 2. + r2 - sqrt(r2 * r2 - shight / 2. * shight / 2.) + layer12);
0215
0216   Box gfbox("fbox", 0.5 * fLens[0], 0.5 * fLens[1], 0.5 * fLens[2]);
0217   Box gcbox("cbox", 0.5 * fLens[0], 0.5 * fLens[1] + 1*mm, 0.5 * fLens[2]);
0218
0219 //  Volume gfbox_volume("gfbox_volume", gfbox, bar_material);
0220 //  lDirc.placeVolume(gfbox_volume, Position(0, 0, 0));
0221 //
0222 //  Volume gcbox_volume("gcbox_volume", gcbox, bar_material);
0223 //  lDirc.placeVolume(gcbox_volume, Position(0, 0, 50));
0224
0225
0226   Position tTrans1(0.5 * (fLens[0] + shight), 0, -fLens[2] + layer12);
0227   Position tTrans0(-0.5 * (fLens[0] + shight), 0, -fLens[2] + layer12);
0228   SubtractionSolid tubox("tubox", gfbox, gcbox, tTrans1);
0229   SubtractionSolid gubox("gubox", tubox, gcbox, tTrans0);
0230
0231 //  Volume tubox_volume("tubox_volume", tubox, bar_material);
0232 //  lDirc.placeVolume(tubox_volume, Position(0, 0, 100));
0233 //
0234 //  Volume gubox_volume("gubox_volume", gubox, bar_material);
0235 //  lDirc.placeVolume(gubox_volume, Position(0, 0, 150));
0236
0237   Tube gcylinder1("Cylinder1", 0, r1, 0.5 * fLens[1], 0 * deg, 360 * deg);
0238   Tube gcylinder2("Cylinder2", 0, r2, 0.5 * fLens[1] - 0.5*mm, 0 * deg, 360 * deg);
0239   Tube gcylinder1c("Cylinder1c", 0, r1, 0.5 * fLens[1] + 0.5*mm, 0 * deg, 360 * deg);
0240   Tube gcylinder2c("Cylinder2c", 0, r2, 0.5 * fLens[1] + 0.5*mm, 0 * deg, 360 * deg);
0241   RotationX xRot(-M_PI / 2.);
0242
0243   IntersectionSolid gLens1("Lens1", gubox, gcylinder1, Transform3D(xRot, zTrans1));
0244   SubtractionSolid gLenst("temp", gubox, gcylinder1c, Transform3D(xRot, zTrans1));
0245
0246 //  Volume gLens1_volume("gLens1_volume", gLens1, bar_material);
0247 //  lDirc.placeVolume(gLens1_volume, Position(0, 0, 200));
0248 //
0249 //  Volume gLenst_volume("gLenst_volume", gLenst, bar_material);
0250 //  lDirc.placeVolume(gLenst_volume, Position(0, 0, 250));
0251
0252   IntersectionSolid gLens2("Lens2", gLenst, gcylinder2, Transform3D(xRot, zTrans2));
0253   SubtractionSolid gLens3("Lens3", gLenst, gcylinder2c, Transform3D(xRot, zTrans2));
0254
0255   lLens1 = Volume("lLens1", gLens1, bar_material);
0256   lLens2 = Volume("lLens2", gLens2, nlak33a);
0257   lLens3 = Volume("lLens3", gLens3, bar_material);
0258
0259   lLens1.setVisAttributes(desc.visAttributes("DIRCLens1"));
0260   lLens2.setVisAttributes(desc.visAttributes("DIRCLens2"));
0261   lLens3.setVisAttributes(desc.visAttributes("DIRCLens3"));
0262
0263
0264   double shifth = 0.5 * (bar_assm_length + fLens[2]);
0265   // fmt::print("LENS HERE shifth={}\n", shifth);
0266
0267   lLens1.setVisAttributes(desc.visAttributes("AnlTeal"));
0268   dirc_module.placeVolume(lLens1, Position(0, 0, shifth));
0269   dirc_module.placeVolume(lLens2, Position(0, 0, shifth));
0270   dirc_module.placeVolume(lLens3, Position(0, 0, shifth));
0271
0272   // The Prizm
0273   Trap gPrizm = MakeTrap("gPrizm", fPrizm[0], fPrizm[1], fPrizm[2], fPrizm[3]);
0274   Volume lPrizm("lPrizm", gPrizm, bar_material);
0275   lPrizm.setVisAttributes(desc.visAttributes("DIRCPrism"));
0276
0277   //G4RotationMatrix *fdRot = new G4RotationMatrix();
0278   //G4RotationMatrix *fdrot = new G4RotationMatrix();
0279   double evshiftz = 0.5 * bar_assm_length + fPrizm[1] + fMcpActive[2] / 2. + fLens[2];
0280   double evshiftx = -3*mm;
0281
0282   double prism_shift_x = (fPrizm[2] + fPrizm[3]) / 4. - 0.5 * fPrizm[3] + 1.5*mm;
0283   double prism_shift_z = 0.5 * (bar_assm_length + fPrizm[1]) + fLens[2];
0284
0285   Position  fPrismShift(prism_shift_x, 0, prism_shift_z);
0286   dirc_module.placeVolume(lPrizm, Transform3D(xRot, fPrismShift));
0287   dirc_module.placeVolume(lFd, Position(0.5 * fFd[1] - 0.5 * fPrizm[3] - evshiftx, 0, evshiftz));
0288
0289   double dphi = 2 * M_PI / (double)fNBoxes;
0290   for (int i = 0; i < fNBoxes; i++) {
0291     double phi = dphi * i;
0292     double dx = -fRadius * cos(phi);
0293     double dy = -fRadius * sin(phi);
0294
0295     //G4RotationMatrix *tRot = new G4RotationMatrix();
0296
0297     Transform3D  tr(RotationZ(phi+M_PI), Position(dx, dy, 0));
0298     PlacedVolume box_placement = det_volume.placeVolume(dirc_module, tr);
0299     box_placement.addPhysVolID("module", i);
0300
0301
0302     // fmt::print("placing dircbox # {} -tphi={:.0f} dx={:.0f}, dy={:.0f}\n", i, phi/deg, dx/cm, dy/cm);
0303
0304     //new G4PVPlacement(tRot, G4ThreeVector(dx, dy, 0), lDirc, "wDirc", lExpHall, false, i);
0305   }
0306
0307
0308   //////////////////
0309   // DIRC Bars
0310   //////////////////
0311
0312   // double bar_radius   = 83.65 * cm;
0313   // double bar_length   = SizeZ;
0314   // double bar_width    = 42. * cm;
0315   // double bar_thicknes = 1.7 * cm;
0316   // int    bar_count    = 2 * M_PI * bar_radius / bar_width;
0317   // double bar_dphi     = 2 * 3.1415926 / bar_count;
0318   // Material bar_material = desc.material("Quartz");
0319
0320   // Box    bar_geo(bar_thicknes / 2., bar_width / 2., bar_length / 2.);
0321   // Volume bar_volume("cb_DIRC_bars_Logix", bar_geo, bar_material);
0322   // bar_volume.setVisAttributes(desc.visAttributes(xml_det.visStr()));
0323   // sens.setType("tracker");
0324   // bar_volume.setSensitiveDetector(sens);
0325
0326   // for (int ia = 0; ia < bar_count; ia++) {
0327   //   double phi = (ia * (bar_dphi));
0328   //   double x   = -bar_radius * cos(phi);
0329   //   double y   = -bar_radius * sin(phi);
0330
0331   //   Transform3D  tr(RotationZ(bar_dphi * ia), Position(x, y, 0));
0332   //   PlacedVolume barPV = det_volume.placeVolume(bar_volume, tr);
0333   //   barPV.addPhysVolID("module", ia);
0334   // }
0335   return det;
0336 }
0337
0338 DECLARE_DETELEMENT(cb_DIRC, createDetector)
0339
0340 dd4hep::Trap MakeTrap( const std::string& pName, double pZ, double pY, double pX, double pLTX )
0341 {
0342     // Fixed Trap constructor. This function is a workaround of this bug:
0343     // https://github.com/AIDASoft/DD4hep/issues/850
0344     // Should be used instead of dd4hep::Trap(pName, pZ, pY, pX, pLTX) constructor
0345
0346     double fDz  = 0.5*pZ;
0347     double fTthetaCphi = 0;
0348     double fTthetaSphi = 0;
0349     double fDy1 = 0.5*pY;
0350     double fDx1 = 0.5*pX;
0351     double fDx2 = 0.5*pLTX;
0352     double fTalpha1 = 0.5*(pLTX - pX)/pY;
0353     double fDy2 = fDy1;
0354     double fDx3 = fDx1;
0355     double fDx4 = fDx2;
0356     double fTalpha2 = fTalpha1;
0357
0358     return Trap(pName, fDz,  fTthetaCphi,  fTthetaSphi,
0359                 fDy1,  fDx1,  fDx2,  fTalpha1,
0360                 fDy2,  fDx3,  fDx4,  fTalpha2);
0361 }
0362