Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-11-25 09:25:09 

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2025 Akio Ogawa
0003 
0004 //==========================================================================
0005 //  Implementation of forward calorimeter with 2025 design and
0006 //  the insert shape cut out, with ScFi
0007 //--------------------------------------------------------------------------
0008 //  Author: Akio Ogawa (BNL)
0009 //==========================================================================
0010 
0011 #include "DD4hep/DetFactoryHelper.h"
0012 #include "DD4hep/Printout.h"
0013 #include <XML/Helper.h>
0014 #include <XML/Layering.h>
0015 #include "forwardEcalMap.h"
0016 
0017 using namespace dd4hep;
0018 
0019 double blocksize, blockgap, offsetX[2], offsetY[2], insert_dx[2];
0020 
0021 double xBlock(int ns, int row, int col) {
0022   if (row >= 18 && row <= 20) {
0023     return (1 - 2 * ns) * (offsetX[ns] + insert_dx[ns] + (blocksize + blockgap) * (col + 0.5));
0024   } else {
0025     return (1 - 2 * ns) * (offsetX[ns] + (blocksize + blockgap) * (col + 0.5));
0026   }
0027 }
0028 
0029 double yBlock(int ns, int row) {
0030   return offsetY[ns] + (blocksize + blockgap) * (mMaxRowBlock / 2.0 - row - 0.5);
0031 }
0032 
0033 static Ref_t createDetector(Detector& desc, xml_h handle, SensitiveDetector sens) {
0034   blocksize    = desc.constant<double>("EcalEndcapP_blockSize");
0035   blockgap     = desc.constant<double>("EcalEndcapP_spaceBetweenBlock");
0036   double nsgap = desc.constant<double>("EcalEndcapP_xOffsetNorth") +
0037                  desc.constant<double>("EcalEndcapP_xOffsetSouth");
0038   double rmin             = 0.0; // Dummy variable. Set to 0 since cutting out insert
0039   double rmax             = desc.constant<double>("EcalEndcapP_rmax");
0040   double rmaxWithGap      = desc.constant<double>("EcalEndcapP_rmaxWithGap");
0041   double zmax             = desc.constant<double>("EcalEndcapP_zmax");
0042   double length           = desc.constant<double>("EcalEndcapP_length");
0043   double zmin             = desc.constant<double>("EcalEndcapP_zmin");
0044   offsetX[0]              = desc.constant<double>("EcalEndcapP_xOffsetNorth");
0045   offsetX[1]              = desc.constant<double>("EcalEndcapP_xOffsetSouth");
0046   offsetY[0]              = desc.constant<double>("EcalEndcapP_yOffsetNorth");
0047   offsetY[1]              = desc.constant<double>("EcalEndcapP_yOffsetSouth");
0048   insert_dx[0]            = desc.constant<double>("EcalEndcapP_xOffsetBeamPipeNorth");
0049   insert_dx[1]            = desc.constant<double>("EcalEndcapP_xOffsetBeamPipeSouth");
0050   double insert_dy        = desc.constant<double>("EcalEndcapP_insert_dy");
0051   double insert_dz        = desc.constant<double>("EcalEndcapP_insert_dz");
0052   double insert_thickness = desc.constant<double>("EcalEndcapP_insert_thickness");
0053   double insert_x         = desc.constant<double>("EcalEndcapP_insert_center_x");
0054   int nx                  = desc.constant<int>("EcalEndcapP_nfiber_x");
0055   int ny                  = desc.constant<int>("EcalEndcapP_nfiber_y");
0056   double rFiber           = desc.constant<double>("EcalEndcapP_rfiber");
0057   double rScfi            = desc.constant<double>("EcalEndcapP_rscfi");
0058 
0059   const double phi1[2]  = {-M_PI / 2.0, M_PI / 2.0};
0060   const double phi2[2]  = {M_PI / 2.0, 3.0 * M_PI / 2.0};
0061   const char* nsName[2] = {"N", "S"};
0062   const double pm[2]    = {1.0, -1.0}; //positive x for north, and negative for south
0063 
0064   //from compact files
0065   xml_det_t detElem   = handle;
0066   std::string detName = detElem.nameStr();
0067   int detID           = detElem.id();
0068 
0069   int Homogeneous_Scfi = 0;
0070   Homogeneous_Scfi     = desc.constant<int>("EcalEndcapP_Homogeneous_ScFi");
0071   if (Homogeneous_Scfi <= 1)
0072     printout(INFO, "FEMC", "Making Homogeneous geometry model\n");
0073   else
0074     printout(INFO, "FEMC", "Making ScFi geometry model\n");
0075 
0076   //xml_dim_t dim = detElem.dimensions();
0077   //xml_dim_t pos = detElem.position();
0078   if (zmax != mBackPlateZ)
0079     printf("WARNING!!! forwardEcal_geo.cpp detect inconsistent Z pos %f(compact) %f(map)\n", zmax,
0080            mBackPlateZ);
0081   //printf("forwardEcal_geo : zmax= %f vs %f\n",zmax,mBackPlateZ);
0082 
0083   PlacedVolume pv;
0084   Material air     = desc.material("Air");
0085   Material alumi   = desc.material("Aluminum5083"); //actually using 6061... does it matter?
0086   Material Wpowder = desc.material("WPowderplusEpoxy");
0087   Material PMMA    = desc.material("Plexiglass");
0088   Material ScFi    = desc.material("Polystyrene");
0089 
0090   // Defining envelope with full phi,with slightly increased radius for NS gap
0091   Tube envelope(rmin, rmaxWithGap, length / 2.0);
0092 
0093   // Removing insert shape from envelope
0094   Box insert((insert_dx[0] + insert_dx[1] - insert_thickness * 2 + nsgap) / 2.0,
0095              (insert_dy - insert_thickness * 2) / 2.0, length / 2.);
0096   SubtractionSolid envelope_with_inserthole(envelope, insert, Position(insert_x, 0.0, 0.0));
0097   Volume envelopeVol(detName, envelope_with_inserthole, air);
0098   envelopeVol.setAttributes(desc, detElem.regionStr(), detElem.limitsStr(), detElem.visStr());
0099 
0100   //double thickness=0.0;
0101   //int slice_num  = 1;
0102   double slice_z = -length / 2.0; // Keeps track of slices' z locations in each layer
0103   // Looping over each layer's slices
0104   for (xml_coll_t sl(detElem, _U(slice)); sl; ++sl) {
0105     xml_comp_t x_slice     = sl;
0106     double slice_thickness = x_slice.thickness();
0107     //thickness+=slice_thickness;
0108     //printf("forwardEcal slice=%1d %8.4f %s \n",slice_num,slice_thickness,x_slice.materialStr().c_str());
0109     std::string slice_name = "fEcal" + x_slice.nameStr();
0110     Material slice_mat     = desc.material(x_slice.materialStr());
0111     slice_z += slice_thickness / 2.; // Going to slice halfway point
0112     Tube slice(rmin, rmaxWithGap, slice_thickness / 2.);
0113 
0114     // Removing insert shape from each slice
0115     Box slice_insert((insert_dx[0] + insert_dx[1] + nsgap) / 2.0, insert_dy / 2.0,
0116                      slice_thickness / 2.0);
0117     SubtractionSolid slice_with_inserthole(slice, slice_insert, Position(insert_x, 0.0, 0.0));
0118     Volume slice_vol(slice_name, slice_with_inserthole, air); //Still air
0119     slice_vol.setAttributes(desc, x_slice.regionStr(), x_slice.limitsStr(), x_slice.visStr());
0120     pv = envelopeVol.placeVolume(slice_vol,
0121                                  Transform3D(RotationZYX(0, 0, 0), Position(0., 0., slice_z)));
0122 
0123     //Loop over north and south halves
0124     for (int ns = 0; ns < 2; ns++) {
0125       Tube half(rmin, rmax, slice_thickness / 2.0, phi1[ns], phi2[ns]);
0126       std::string half_name = slice_name + "_" + nsName[ns];
0127 
0128       // Removing insert shape from each slice & halves
0129       Box half_insert(insert_dx[ns] / 2.0, insert_dy / 2.0, slice_thickness / 2.0);
0130       SubtractionSolid half_with_inserthole(half, half_insert,
0131                                             Position(pm[ns] * insert_dx[ns] / 2.0, 0.0, 0.0));
0132 
0133       Material mat = slice_mat;
0134       if (x_slice.isSensitive())
0135         mat = air; //for calorimeter itself, still air to place blocks inside
0136       Volume half_vol(half_name, half_with_inserthole, mat);
0137       half_vol.setAttributes(desc, x_slice.regionStr(), x_slice.limitsStr(), x_slice.visStr());
0138       pv = slice_vol.placeVolume(
0139           half_vol, Transform3D(RotationZYX(0, 0, 0), Position(pm[ns] * nsgap / 2.0, 0.0, 0.0)));
0140       pv.addPhysVolID("northsouth", ns);
0141 
0142       // For detector (sensitive) slice, placing detector blocks in col and row
0143       double bsize = blocksize + blockgap;
0144       if (x_slice.isSensitive()) {
0145 
0146         //Define WSiFi block (4x4 towers)
0147         if (Homogeneous_Scfi <= 1)
0148           mat = slice_mat;
0149         Box block(blocksize / 2.0, blocksize / 2.0, slice_thickness / 2.0);
0150         Volume block_vol("FEMCBlock", block, mat);
0151         block_vol.setAttributes(desc, x_slice.regionStr(), x_slice.limitsStr(), x_slice.visStr());
0152         if (Homogeneous_Scfi <= 1) {
0153           sens.setType("calorimeter");
0154           block_vol.setSensitiveDetector(sens);
0155         } // end if Homogeneous_Scfi<=1
0156 
0157         if (Homogeneous_Scfi == 2) {
0158           //4 rows of towers
0159           Box trow(blocksize / 2.0, blocksize / 8.0, slice_thickness / 2.0);
0160           Volume trow_vol("FEMCTowerRow", trow, air);
0161           trow_vol.setAttributes(desc, x_slice.regionStr(), x_slice.limitsStr(), x_slice.visStr());
0162           for (int tr = 0; tr < 4; tr++) {
0163             pv = block_vol.placeVolume(
0164                 trow_vol,
0165                 Transform3D(RotationZYX(0, 0, 0), Position(0.0, (tr - 1.5) * blocksize / 4, 0)));
0166             pv.addPhysVolID("y", tr);
0167           }
0168 
0169           //4 towers in a row - finally a W powder volume, not air
0170           Box tower(blocksize / 8.0, blocksize / 8.0, slice_thickness / 2.0);
0171           Volume tower_vol("FEMCTower", tower, Wpowder);
0172           tower_vol.setAttributes(desc, x_slice.regionStr(), x_slice.limitsStr(), x_slice.visStr());
0173           for (int tc = 0; tc < 4; tc++) {
0174             pv = trow_vol.placeVolume(
0175                 tower_vol,
0176                 Transform3D(RotationZYX(0, 0, 0), Position((tc - 1.5) * blocksize / 4, 0, 0)));
0177             pv.addPhysVolID("x", tc);
0178           }
0179 
0180           //rows of fibers
0181           double fiberDistanceX =
0182               blocksize / 4.0 / (nx + 0.5); //exrea 0.5 for even/odd rows shifted by 1/2
0183           Box frow(blocksize / 8.0 - fiberDistanceX / 2.0, blocksize / 8.0 / ny,
0184                    slice_thickness / 2.0);
0185           Volume frow_vol("FEMCFiberRow", frow, Wpowder);
0186           frow_vol.setAttributes(desc, x_slice.regionStr(), x_slice.limitsStr(), x_slice.visStr());
0187           for (int iy = 0; iy < ny; iy++) {
0188             double xx = 0;
0189             if (iy % 2 == 1)
0190               xx += fiberDistanceX / 2.0;
0191             pv = tower_vol.placeVolume(
0192                 frow_vol,
0193                 Transform3D(RotationZYX(0, 0, 0),
0194                             Position(xx, (iy - ny / 2.0 + 0.5) * blocksize / 4.0 / ny, 0)));
0195             //printf("iy=%2d dy=%8.4f fiberRx2=%8.4f xx=%8.4f\n",iy,blocksize/4.0/ny,rFiber*2,xx);
0196             pv.addPhysVolID("fiber_y", iy);
0197           }
0198 
0199           //columns of fibers, with 1/2 fiber distance shifted each row
0200           Box fcol(fiberDistanceX / 2.0, blocksize / 8.0 / ny, slice_thickness / 2.0);
0201           Volume fcol_vol("FEMCFiberCol", fcol, Wpowder);
0202           fcol_vol.setAttributes(desc, x_slice.regionStr(), x_slice.limitsStr(), x_slice.visStr());
0203           for (int ix = 0; ix < nx; ix++) {
0204             double xx = (ix - nx / 2.0 + 0.5) * fiberDistanceX;
0205             pv        = frow_vol.placeVolume(fcol_vol,
0206                                              Transform3D(RotationZYX(0, 0, 0), Position(xx, 0, 0)));
0207             //printf("ix=%2d dx=%8.4f xx=%8.4f x0=%8.4f x1=%8.4f\n",ix,fiberDistanceX,xx,xx-fiberDistanceX/2,xx+fiberDistanceX/2);
0208             pv.addPhysVolID("fiber_x", ix);
0209           }
0210 
0211           //a fiber (with coating material, not sensitive yet)
0212           Tube fiber(0, rFiber, slice_thickness / 2.0);
0213           Volume fiber_vol("FEMCFiber", fiber, PMMA);
0214           fiber_vol.setAttributes(desc, x_slice.regionStr(), x_slice.limitsStr(), x_slice.visStr());
0215           pv =
0216               fcol_vol.placeVolume(fiber_vol, Transform3D(RotationZYX(0, 0, 0), Position(0, 0, 0)));
0217 
0218           //scintillating fiber core - and finally a sensitive volume
0219           Tube scfi(0, rScfi, slice_thickness / 2.0);
0220           Volume scfi_vol("FEMCScFi", scfi, ScFi);
0221           scfi_vol.setAttributes(desc, x_slice.regionStr(), x_slice.limitsStr(), x_slice.visStr());
0222           pv =
0223               fiber_vol.placeVolume(scfi_vol, Transform3D(RotationZYX(0, 0, 0), Position(0, 0, 0)));
0224           sens.setType("calorimeter");
0225           scfi_vol.setSensitiveDetector(sens);
0226         } //end if Homogeneous_Scfi==2
0227 
0228         //rows of blocks
0229         //int nRowBlock = map->maxRowBlock(); //# of rows
0230         int nRowBlock = mMaxRowBlock;
0231         for (int r = 0; r < nRowBlock; r++) {
0232           //int nColBlock = map->nColBlock(ns, r);
0233           int nColBlock = mNColBlock[ns][r];
0234           double dxrow  = bsize * nColBlock;
0235           Box row(dxrow / 2.0, bsize / 2.0, slice_thickness / 2.0);
0236           std::string row_name = half_name + _toString(r, "_R%02d");
0237           Volume row_vol(row_name, row, air);
0238           row_vol.setAttributes(desc, x_slice.regionStr(), x_slice.limitsStr(), x_slice.visStr());
0239           //double xrow = (map->xBlock(ns, r, 0) + map->xBlock(ns, r, nColBlock - 1)) / 2.0 -
0240           //              pm[ns] * nsgap / 2.0;
0241           //double yrow = map->yBlock(ns, r);
0242           double xrow =
0243               (xBlock(ns, r, 0) + xBlock(ns, r, nColBlock - 1)) / 2.0 - pm[ns] * nsgap / 2.0;
0244           double yrow = yBlock(ns, r);
0245           pv          = half_vol.placeVolume(row_vol,
0246                                              Transform3D(RotationZYX(0, 0, 0), Position(xrow, yrow, 0)));
0247           pv.addPhysVolID("blockrow", r);
0248 
0249           //column of blocks
0250           double xcol = -pm[ns] * (dxrow / 2.0 - bsize / 2.0);
0251           for (int c = 0; c < nColBlock; c++) {
0252             Box col(bsize / 2.0, bsize / 2.0, slice_thickness / 2.0);
0253             std::string col_name = row_name + _toString(c, "C%02d");
0254             Volume col_vol(col_name, col, air);
0255             col_vol.setAttributes(desc, x_slice.regionStr(), x_slice.limitsStr(), x_slice.visStr());
0256             pv = row_vol.placeVolume(col_vol,
0257                                      Transform3D(RotationZYX(0, 0, 0), Position(xcol, 0, 0)));
0258             pv.addPhysVolID("blockcol", c);
0259             //printf("r=%2d dx=%8.3f x=%8.3f y=%8.3f   c=%2d %8.3f  mapx=%8.3f\n",r,dxrow,xrow,yrow,c,xcol,map->xBlock(ns,r,c));
0260             xcol += pm[ns] * bsize;
0261 
0262             //a block inside with air gap
0263             pv = col_vol.placeVolume(block_vol,
0264                                      Transform3D(RotationZYX(0, 0, 0), Position(0, 0, 0)));
0265           } //end loop block col
0266         } //end loop block raw
0267       } //end if isSensitive
0268     } //end loop over ns
0269     slice_z += slice_thickness / 2.; // Going to end of slice
0270     //++slice_num;
0271   } //end loop over slice
0272   //printf("forwardEcal Total thickness=%f Slice end at %f\n",thickness,slice_z);
0273 
0274   //Al Beampipe Protector placed in envelope volume outside slices
0275   for (int ns = 0; ns < 2; ns++) {
0276     Box bpp(insert_dx[ns] / 2.0, insert_dy / 2.0, insert_dz / 2.0);
0277     std::string bpp_name = detName + "_BeamPipeProtector_" + nsName[ns];
0278     Box bpp_hole((insert_dx[ns] - insert_thickness) / 2.0, insert_dy / 2.0 - insert_thickness,
0279                  insert_dz / 2.0);
0280     SubtractionSolid bpp_with_hole(bpp, bpp_hole,
0281                                    Position(-pm[ns] * insert_thickness / 2.0, 0.0, 0.0));
0282     Volume bpp_vol(bpp_name, bpp_with_hole, alumi);
0283     bpp_vol.setAttributes(desc, detElem.regionStr(), detElem.limitsStr(), detElem.visStr());
0284     pv = envelopeVol.placeVolume(
0285         bpp_vol, Transform3D(RotationZYX(0, 0, 0), Position(pm[ns] * (insert_dx[ns] + nsgap) / 2.0,
0286                                                             0.0, (length - insert_dz) / 2.0)));
0287   }
0288 
0289   // Placing in the world volume
0290   DetElement det(detName, detID);
0291   Volume motherVol = desc.pickMotherVolume(det);
0292   auto tr          = Transform3D(Position(0.0, 0.0, zmin + length / 2.0));
0293   pv               = motherVol.placeVolume(envelopeVol, tr);
0294   pv.addPhysVolID("system", detID);
0295   det.setPlacement(pv);
0296 
0297   //printf("forwardEcal_geo Done\n");
0298   return det;
0299 }
0300 DECLARE_DETELEMENT(epic_ForwardEcal, createDetector)
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!