From 7dd79da83d61ce4100ac4538c13bf392006c9934 Mon Sep 17 00:00:00 2001
From: lintao <lintao51@gmail.com>
Date: Mon, 10 Aug 2020 11:26:26 +0800
Subject: [PATCH] Import SEcal05 detector constructor from Chengdong Fu.
---
Detector/DetCEPCv4/CMakeLists.txt | 4 +
.../src/calorimeter/SEcal05_Barrel.cpp | 444 +++++++++++
.../src/calorimeter/SEcal05_ECRing.cpp | 689 +++++++++++++++++
.../src/calorimeter/SEcal05_Endcaps.cpp | 390 ++++++++++
.../src/calorimeter/SEcal05_Helpers.cpp | 729 ++++++++++++++++++
.../src/calorimeter/SEcal05_Helpers.h | 309 ++++++++
6 files changed, 2565 insertions(+)
create mode 100644 Detector/DetCEPCv4/src/calorimeter/SEcal05_Barrel.cpp
create mode 100644 Detector/DetCEPCv4/src/calorimeter/SEcal05_ECRing.cpp
create mode 100644 Detector/DetCEPCv4/src/calorimeter/SEcal05_Endcaps.cpp
create mode 100644 Detector/DetCEPCv4/src/calorimeter/SEcal05_Helpers.cpp
create mode 100644 Detector/DetCEPCv4/src/calorimeter/SEcal05_Helpers.h
diff --git a/Detector/DetCEPCv4/CMakeLists.txt b/Detector/DetCEPCv4/CMakeLists.txt
index 7211736c..ff3a8a70 100644
--- a/Detector/DetCEPCv4/CMakeLists.txt
+++ b/Detector/DetCEPCv4/CMakeLists.txt
@@ -26,6 +26,10 @@ set(DetCEPCv4_src
src/tracker/SIT_Simple_Pixel_geo.cpp
src/tracker/TPC10_geo.cpp
src/tracker/SET_Simple_Planar_geo.cpp
+ src/calorimeter/SEcal05_Helpers.cpp
+ src/calorimeter/SEcal05_Barrel.cpp
+ src/calorimeter/SEcal05_Endcaps.cpp
+ src/calorimeter/SEcal05_ECRing.cpp
src/other/BoxSupport_o1_v01_geo.cpp
src/other/TubeSupport_o1_v01_geo.cpp
)
diff --git a/Detector/DetCEPCv4/src/calorimeter/SEcal05_Barrel.cpp b/Detector/DetCEPCv4/src/calorimeter/SEcal05_Barrel.cpp
new file mode 100644
index 00000000..a0329b4e
--- /dev/null
+++ b/Detector/DetCEPCv4/src/calorimeter/SEcal05_Barrel.cpp
@@ -0,0 +1,444 @@
+//====================================================================
+// lcgeo - LC detector models in DD4hep
+//--------------------------------------------------------------------
+// DD4hep Geometry driver for SiWEcalBarrel
+// Ported from Mokka
+//--------------------------------------------------------------------
+// S.Lu, DESY
+// D.Jeans, Tokyo <-- to 05: add possibility to remove preshower layer
+//====================================================================
+
+#include "DD4hep/DetFactoryHelper.h"
+#include "DD4hep/DetType.h"
+
+#include "XML/Layering.h"
+#include "TGeoTrd2.h"
+
+#include "XML/Utilities.h"
+#include "DDRec/DetectorData.h"
+
+#include "DDSegmentation/MegatileLayerGridXY.h"
+#include "DDSegmentation/WaferGridXY.h"
+
+#include "SEcal05_Helpers.h"
+
+#include <sstream>
+
+#undef NDEBUG
+#include <assert.h>
+
+using namespace std;
+
+using dd4hep::BUILD_ENVELOPE;
+using dd4hep::DetElement;
+using dd4hep::Detector;
+using dd4hep::Layering;
+using dd4hep::Material;
+using dd4hep::PlacedVolume;
+using dd4hep::Position;
+using dd4hep::Readout;
+using dd4hep::Ref_t;
+using dd4hep::RotationZYX;
+using dd4hep::Segmentation;
+using dd4hep::SensitiveDetector;
+using dd4hep::Transform3D;
+using dd4hep::Translation3D;
+using dd4hep::Trapezoid;
+using dd4hep::Volume;
+using dd4hep::_toString;
+
+using dd4hep::rec::LayeredCalorimeterData;
+
+#define VERBOSE 1
+
+// workaround for DD4hep v00-14 (and older)
+#ifndef DD4HEP_VERSION_GE
+#define DD4HEP_VERSION_GE(a,b) 0
+#endif
+
+/** SEcal05.cc
+ *
+ * new SEcal05 barrel driver: allows removal of preshower layer. DJeans UTokyo, sep/2016
+ * based on SEcal04
+ *
+ * @author: Shaojun Lu, DESY
+ * @version $Id: SEcal04_Barrel.cpp 1060 2016-09-05 07:48:41Z /C=JP/O=KEK/OU=CRC/CN=JEANS Daniel Thomelin Dietrich $
+ * Ported from Mokka SEcal04 Barrel part. Read the constants from XML
+ * instead of the DB. Then build the Barrel in the same way with DD4hep
+ * construct.
+ *
+ * @history: F.Gaede, CERN/DESY, Nov. 10, 2014
+ * added information for reconstruction: LayeringExtension and surfaces (experimental)
+ * removed DetElement for slices (not needed) increased multiplicity for layer DetElement
+ * along tower index
+ * F.Gaede: 03/2015:
+ * create the envelope volume with create_placed_envelope() using the xml
+ *
+ * D. Jeans: 03/2015:
+ * generalise to non-octagonal barrel shape
+ * allow dead region between end of each slab and the module edge (e.g. for the slab fastening system)
+ */
+
+/*
+
+
+
+
+ y
+ ^
+ | z is perpendicular to page (along beamline)
+ |
+ -----> x
+
+ original coord system was very mixed up.
+ Daniel tried to rationalise it to make it more understandable (hopefully)
+
+ absorber layers are wrapped in a number of layers of carbon fibre composite (CF)
+
+ structure of slab (which is inserted into the structure) is defined in the compact description
+
+ general structure (from +ve y)
+ -------------------------------
+
+ if preshower, front face (CF) \\\\\\\\\\\\\\\\\\
+
+ if no preshower, wrapped absorber ==================
+
+ __________________
+ alveolus, containing the slab |__________________|
+
+ wrapped absorber ===================
+ __________________
+ alveolus, containing the slab |__________________|
+
+ .
+ repeat as required
+ .
+
+ wrapped absorber ===================
+ __________________
+ alveolus, containing the slab |__________________|
+
+ back support plate (CF) \\\\\\\\\\\\\\\\\\\\
+ ////////////////////
+
+ */
+
+
+
+static Ref_t create_detector(Detector& theDetector, xml_h element, SensitiveDetector sens) {
+
+ cout << "-----------------------" << endl;
+ cout << "creating SEcal05_Barrel" << endl;
+ cout << "-----------------------" << endl;
+
+ xml_det_t x_det = element;
+ string det_name = x_det.nameStr();
+ Layering layering (element);
+
+ Material carbon_fibre = theDetector.material("CarbonFiber");
+
+ int det_id = x_det.id();
+ DetElement sdet (det_name,det_id);
+
+ xml_comp_t x_dim = x_det.dimensions();
+ int nsides = x_dim.numsides();
+
+ double dphi = (2*M_PI/nsides);
+ double hphi = dphi/2;
+
+ // --- create an envelope volume and position it into the world ---------------------
+
+ Volume envelope = dd4hep::xml::createPlacedEnvelope( theDetector, element , sdet ) ;
+
+ dd4hep::xml::setDetectorTypeFlag( element, sdet ) ;
+
+ if( theDetector.buildType() == BUILD_ENVELOPE ) return sdet ;
+
+ //-----------------------------------------------------------------------------------
+
+ sens.setType("calorimeter");
+
+ DetElement stave_det("module1stave1",det_id);
+
+ //====================================================================
+ //
+ // Read all the constant from ILD_o1_v05.xml
+ // Use them to build HcalBarrel
+ //
+ //====================================================================
+
+ // some hardcoded values!
+ // const int N_FIBERS_W_STRUCTURE = 2; // number of CF layers around absorber layers in the structure
+ // const int N_FIBERS_ALVEOLUS = 3; // number of CF layers used to make alveolus
+
+ // read other parameters from compact.xml file
+
+ // overall size
+ double Ecal_inner_radius = theDetector.constant<double>("Ecal_inner_radius");
+
+ double Ecal_Barrel_halfZ = theDetector.constant<double>("Ecal_Barrel_halfZ");
+ int Ecal_barrel_z_modules = theDetector.constant<int>("Ecal_barrel_z_modules");
+ int Ecal_barrel_number_of_towers = theDetector.constant<int>("Ecal_barrel_number_of_towers");
+
+ double Ecal_barrel_thickness = theDetector.constant<double>("Ecal_barrel_thickness"); // what's assumed in the compact description
+
+ // absorber layers
+ int Ecal_nlayers1 = theDetector.constant<int>("Ecal_nlayers1");
+ int Ecal_nlayers2 = theDetector.constant<int>("Ecal_nlayers2");
+ int Ecal_nlayers3 = theDetector.constant<int>("Ecal_nlayers3");
+
+ double Ecal_radiator_thickness1 = theDetector.constant<double>("Ecal_radiator_layers_set1_thickness");
+ double Ecal_radiator_thickness2 = theDetector.constant<double>("Ecal_radiator_layers_set2_thickness");
+ double Ecal_radiator_thickness3 = theDetector.constant<double>("Ecal_radiator_layers_set3_thickness");
+
+ // CF support dimensions
+ double Ecal_support_thickness = theDetector.constant<double>("Ecal_support_thickness");
+ double Ecal_front_face_thickness = theDetector.constant<double>("Ecal_front_face_thickness");
+ double Ecal_lateral_face_thickness = theDetector.constant<double>("Ecal_lateral_face_thickness");
+ double Ecal_Slab_H_fiber_thickness = theDetector.constant<double>("Ecal_Slab_H_fiber_thickness");
+
+
+ // double Ecal_fiber_thickness = theDetector.constant<double>("Ecal_fiber_thickness");
+ // some hardcoded values!
+ // const int N_FIBERS_W_STRUCTURE = 2; // number of CF layers around absorber layers in the structure
+ // const int N_FIBERS_ALVEOLUS = 3; // number of CF layers used to make alveolus
+ double Ecal_fiber_thickness_structure = theDetector.constant<double>("Ecal_fiber_thickness_structure"); // absorber wrapping thickness
+ double Ecal_fiber_thickness_alveolus = theDetector.constant<double>("Ecal_fiber_thickness_alveolus"); // alveolar wall thickness
+
+ double Ecal_Slab_shielding = theDetector.constant<double>("Ecal_Slab_shielding");
+
+ // first layer is preshower?
+ bool Ecal_Barrel_PreshowerLayer = theDetector.constant<int>("Ecal_Barrel_Preshower") > 0;
+
+ // internal dimensions of slab
+ double Ecal_guard_ring_size = theDetector.constant<double>("Ecal_guard_ring_size");
+ int Ecal_n_wafers_per_tower = theDetector.constant<int>("Ecal_n_wafers_per_tower");
+
+ std::string Ecal_layerConfig = theDetector.constant<string>("Ecal_layer_pattern");
+
+ int Ecal_end_of_slab_strategy = theDetector.constant<int>("Ecal_end_of_slab_strategy");
+
+ int Ecal_cells_across_megatile = theDetector.constant <int> ("Ecal_cells_across_megatile" );
+ int Ecal_strips_across_megatile = theDetector.constant <int> ("Ecal_strips_across_megatile");
+ int Ecal_strips_along_megatile = theDetector.constant <int> ("Ecal_strips_along_megatile" );
+
+ float Ecal_plugLength = 0.;
+ try {
+ Ecal_plugLength = theDetector.constant<double>("Ecal_Slab_Plug_length");
+ } catch (std::runtime_error&e) {
+ cout << "SEcal05_Barrel: Ecal_plugLength not found, using " << Ecal_plugLength << endl;
+ }
+
+ //---------------------------------
+
+ // set up the helper, which will make a module for us
+
+ SEcal05_Helpers helper;
+ helper.setDet( & x_det );
+
+ // absorber layer structure
+ helper.setPreshower( Ecal_Barrel_PreshowerLayer );
+
+ cout << "SEcal05_Barrel: Preshower ? " << Ecal_Barrel_PreshowerLayer << endl;
+
+ helper.setAbsLayers( Ecal_nlayers1, Ecal_radiator_thickness1,
+ Ecal_nlayers2, Ecal_radiator_thickness2,
+ Ecal_nlayers3, Ecal_radiator_thickness3 );
+
+ cout << "SEcal05_Barrel: absorber layers " <<
+ Ecal_nlayers1 << "*" << Ecal_radiator_thickness1/dd4hep::mm << " mm + " <<
+ Ecal_nlayers2 << "*" << Ecal_radiator_thickness2/dd4hep::mm << " mm + " <<
+ Ecal_nlayers3 << "*" << Ecal_radiator_thickness3/dd4hep::mm << " mm " << endl;
+
+ helper.checkLayerConsistency();
+
+ // structural CF thicknesses
+ helper.setCFthickness( Ecal_fiber_thickness_structure, // was N_FIBERS_W_STRUCTURE*Ecal_fiber_thickness, : updated DJ
+ Ecal_fiber_thickness_alveolus, // N_FIBERS_ALVEOLUS*Ecal_fiber_thickness,
+ Ecal_front_face_thickness,
+ Ecal_support_thickness);
+
+ helper.setPlugLength( Ecal_plugLength );
+
+ // check resulting thickness is consistent with what's in compact description
+ float module_thickness = helper.getTotalThickness();
+ if ( fabs( Ecal_barrel_thickness - module_thickness ) > 0.01 ) {
+ cout << "ERROR SEcal05_Barrel : ECAL barrel thickness in comapct decription not consistent with what I calculate!" << endl;
+ cout << " calculated = " << module_thickness << " compact description: " << Ecal_barrel_thickness << endl;
+ assert( 0 ); // exit
+ }
+
+ cout << "SEcal05_Barrel : module thickness = " << module_thickness << endl;
+
+ // set up the sensitive layer segmentation
+
+ Readout readout = sens.readout();
+ Segmentation seg = readout.segmentation();
+ helper.setSegmentation( &seg );
+
+ helper.setNCells( Ecal_cells_across_megatile, Ecal_strips_across_megatile, Ecal_strips_along_megatile);
+ helper.setMagicMegatileStrategy ( Ecal_end_of_slab_strategy );
+
+ // layer configuration
+ int ntemp;
+ stringstream stream(Ecal_layerConfig);
+ std::vector < int > layerConfig;
+ while ( stream >> ntemp ) {
+ assert (ntemp>=0 );
+ layerConfig.push_back( ntemp );
+ }
+ cout << "SEcal05_Barrel : layer config: ";
+ for (size_t i=0; i<layerConfig.size(); i++) cout << layerConfig[i] << " ";
+ cout << endl;
+
+ helper.setLayerConfig( layerConfig );
+
+
+ // set the alveolar structure
+
+ // get width of alveolus
+ double Ecal_Barrel_module_dim_z = 2 * Ecal_Barrel_halfZ / Ecal_barrel_z_modules ;
+ double alv_width = ( Ecal_Barrel_module_dim_z - 2.*Ecal_lateral_face_thickness ) / Ecal_barrel_number_of_towers;
+
+ cout << "SEcal05_Barrel : width of module, alveolus = " << Ecal_Barrel_module_dim_z << " " << alv_width << endl;
+
+ std::vector < int > ntowers; ntowers.push_back(Ecal_barrel_number_of_towers);
+
+ helper.setTowersUnits( ntowers,
+ alv_width,
+ Ecal_n_wafers_per_tower,
+ Ecal_lateral_face_thickness,
+ Ecal_fiber_thickness_alveolus + Ecal_Slab_H_fiber_thickness + Ecal_Slab_shielding,
+ Ecal_guard_ring_size );
+
+ // shape of barrel module
+ // to get alignment as in ECAL interface design doc fig10
+ double module_thickness_noSupport = module_thickness - Ecal_support_thickness;
+
+ // double max_dim_x = 2. * tan(M_PI/8.) * Ecal_inner_radius + module_thickness_noSupport/sin(M_PI/4.); // longest side assumes oct
+ // double min_dim_x = max_dim_x - 2*module_thickness; // shorter one (assumes octagon)
+
+ // double max_dim_x = 2. * tan(M_PI/nsides) * Ecal_inner_radius + module_thickness_noSupport / tan( 2*M_PI/nsides ); // longest side
+ double max_dim_x = 2. * tan(M_PI/nsides) * Ecal_inner_radius + module_thickness_noSupport / sin( 2*M_PI/nsides ); // longest side : fix DJeans 03/2017
+ double min_dim_x = max_dim_x - 2*module_thickness/tan( 2*M_PI/nsides ) ; // shorter one
+
+ if ( min_dim_x<0 || max_dim_x<0 ) {
+ std::cout << "SEcal05_Barrel ERROR : requesting too many sides! barrel modules have max/min extent " << max_dim_x << " " << min_dim_x << std::endl;
+ std::cout << "exiting" << std::endl;
+ assert(0); // exit gracefully
+ }
+
+ Trapezoid trd(max_dim_x / 2,
+ min_dim_x / 2,
+ Ecal_Barrel_module_dim_z / 2,
+ Ecal_Barrel_module_dim_z / 2,
+ module_thickness/2);
+
+ Volume mod_vol(det_name+"_module", trd, carbon_fibre); // DJeans 5-sep-2016
+
+ //========== fill data for reconstruction ============================
+ LayeredCalorimeterData* caloData = new LayeredCalorimeterData ;
+
+ helper.setModuleDimensions( 0, // int XYtype, // module shape in XY
+ 1, // int XZtype, // module shape in XZ
+ max_dim_x, // double dX_max, // maximum extent in X
+ -999, // dummy
+ 2.*M_PI/nsides
+ );
+
+ // traslation to get layers within the envelope
+ helper.setTranslation ( Position ( -max_dim_x/2. , -Ecal_Barrel_module_dim_z/2., -module_thickness/2. ) );
+
+ // make the module
+
+ helper.makeModule( mod_vol, stave_det,
+ *caloData,
+ theDetector, sens );
+
+ for (size_t i=0; i<caloData->layers.size(); i++) {
+ caloData->layers[i].distance += Ecal_inner_radius; // add IP->front face distance
+ }
+
+ // cout << "SEcal05_Barrel : cell sizes: " << endl;
+ // for (size_t i=0; i<caloData->layers.size(); i++) {
+ // cout << "sensitive layer " << i << " : x,y = " << caloData->layers[i].cellSize0 << " " << caloData->layers[i].cellSize1 << endl;
+ // }
+
+ // ------------- create extension objects for reconstruction -----------------
+
+ caloData->layoutType = LayeredCalorimeterData::BarrelLayout ;
+ caloData->inner_symmetry = nsides ;
+ caloData->outer_symmetry = nsides ;
+ caloData->phi0 = 0 ; // hardcoded
+
+ // extent of the calorimeter in the r-z-plane [ rmin, rmax, zmin, zmax ] in mm.
+ caloData->extent[0] = Ecal_inner_radius ;
+ caloData->extent[1] = ( Ecal_inner_radius + module_thickness );
+ caloData->extent[2] = 0. ;
+ caloData->extent[3] = Ecal_Barrel_halfZ ;
+ //-------------------------------------------------------
+
+ //====================================================================
+ // Place ECAL Barrel stave module into the envelope volume
+ //====================================================================
+
+ double X = Ecal_inner_radius + module_thickness / 2.;
+
+ // altered to get alignment as in ECAL interface design doc fig10
+ // end of slabs aligned to inner face of support plate in next stave (not the outer surface)
+ // double Y = (module_thickness/2.) / sin(M_PI/4.);
+ double Y = (module_thickness_noSupport/2.) / sin(2.*M_PI/nsides);
+
+ // stave numbering from 1->8
+ // stave = 1 is in +ve x direction
+ // stave = 3 is in +ve y direction (ie at top of barrel)
+ // module index from 1->5
+ // module = 1 is at -ve z
+
+ for (int istave = 0; istave < nsides ; istave++) {
+ //for (int istave = 0; istave < 1 ; istave++) {
+ int stave_id = istave+1;
+
+ // dstave is the change in stave index from the top module to the one with smallest positive phi (which has istave=0)
+ int dstave = int( nsides/4. );
+
+
+ // rotations around the center to get the modules in the correct orientation
+ // top module (dstave) shits by hphi + pi/4
+ double phirot = hphi + M_PI/2.; //
+ phirot += (istave - dstave)*dphi;
+
+
+ // this angle is used to get the stave in the right position wrt the origin
+ double phirot2 = (istave - dstave ) * dphi + hphi;
+
+ for (int imodule = 0; imodule < Ecal_barrel_z_modules; imodule++) {
+ int module_id = imodule+1;
+ Transform3D tr( RotationZYX( 0 , phirot, M_PI/2.), // magic rotation!
+ Translation3D(
+ ( X*cos(phirot2)-Y*sin(phirot2) ) ,
+ ( X*sin(phirot2)+Y*cos(phirot2) ) ,
+ ( imodule+0.5 )*Ecal_Barrel_module_dim_z - Ecal_Barrel_halfZ )
+ );
+ PlacedVolume pv = envelope.placeVolume(mod_vol,tr);
+ pv.addPhysVolID("module",module_id);
+ pv.addPhysVolID("stave",stave_id);
+
+ DetElement sd = (imodule==0 && istave==0) ? stave_det : stave_det.clone(_toString(module_id,"module%d")+_toString(stave_id,"stave%d"));
+ sd.setPlacement(pv);
+ sdet.add(sd);
+ }
+ }
+
+ // Set envelope volume attributes.
+ envelope.setAttributes(theDetector,x_det.regionStr(),x_det.limitsStr(),x_det.visStr());
+
+ sdet.addExtension< LayeredCalorimeterData >( caloData ) ;
+
+ // cout << "finished SEcal05_Barrel" << endl;
+
+ return sdet;
+}
+
+DECLARE_DETELEMENT(SEcal05_Barrel,create_detector)
diff --git a/Detector/DetCEPCv4/src/calorimeter/SEcal05_ECRing.cpp b/Detector/DetCEPCv4/src/calorimeter/SEcal05_ECRing.cpp
new file mode 100644
index 00000000..9499c1fa
--- /dev/null
+++ b/Detector/DetCEPCv4/src/calorimeter/SEcal05_ECRing.cpp
@@ -0,0 +1,689 @@
+//====================================================================
+// lcgeo - LC detector models in DD4hep
+//--------------------------------------------------------------------
+// DD4hep Geometry driver for SiWEcalEndcaps
+// Ported from Mokka
+//--------------------------------------------------------------------
+// S.Lu, DESY
+// $Id$
+//====================================================================
+
+/* History:
+
+// *******************************************************
+// * *
+// * Mokka *
+// * - the detailed geant4 simulation for ILC - *
+// * *
+// * For more information about Mokka, visit the *
+// * *
+// * Mokka.in2p3.fr Mokka home page. *
+// * *
+// *******************************************************
+//
+// $Id$
+// $Name: mokka-07-00 $
+//
+//
+//
+// SEcal04.cc
+//
+Shaojun Lu: Ported from Mokka SEcal04 Endcaps part. Read the constants from XML
+instead of the DB. Then build the Endcap in the same way with DD4hep
+construct.
+Inside SEcal04, some parameters, which used by Ecal Endcaps, come from
+Ecal Barrel. They can be seen here again.
+Start ECRing ...
+*/
+
+
+/*
+
+ ==============
+ SEcal05_ECRing
+ ===============
+ SEcal04_ECRing driver originally assumend 2 identical modules in +/- z
+ this is not true, because lumical (and therefore hole in EC plug) is not centred, but offset in +z position
+
+ SEcal05_ECRing driver corrects this (applies offset according to the crossing algle, to center it on the outgoing beampipe).
+ it also deals with the presence or not of a preshower layer
+
+ this is a hack of the 04 version, not a clean rewrite as was done for barrel and endcap, so it's not very clean & tidy....but I think it works
+
+ DJeans jan 2017
+
+ ========
+
+
+ - fixed dd4hep::rec::LayeredCalorimeterData for case with no preshower
+ - small fixes to layer thicknesses (now take from compact description) to get exactly same structure as main endcaps
+
+ DJeans July 2017
+
+
+ ======
+
+ - fixes to LayeredCalorimeterData: "distance" defined to start of layer, some other issues
+ - fix to use correct absorber thickness at transition between stacks
+ - LayeredCalorimeterData material describtion: overall structure made in CF (previously air)
+
+ DJeans Sep 2017
+
+*/
+
+
+#include "DD4hep/DetFactoryHelper.h"
+#include "DD4hep/DetType.h"
+#include "XML/Layering.h"
+#include "XML/Utilities.h"
+#include "DDRec/DetectorData.h"
+
+using namespace std;
+
+using dd4hep::BUILD_ENVELOPE;
+using dd4hep::Box;
+using dd4hep::DetElement;
+using dd4hep::Detector;
+using dd4hep::Layering;
+using dd4hep::Material;
+using dd4hep::PlacedVolume;
+using dd4hep::Position;
+using dd4hep::Readout;
+using dd4hep::Ref_t;
+using dd4hep::Rotation3D;
+using dd4hep::RotationZYX;
+using dd4hep::Segmentation;
+using dd4hep::SensitiveDetector;
+using dd4hep::SubtractionSolid;
+using dd4hep::Transform3D;
+using dd4hep::Tube;
+using dd4hep::Volume;
+using dd4hep::_toString;
+
+using dd4hep::rec::LayeredCalorimeterData;
+
+//#define VERBOSE 1
+
+// workaround for DD4hep v00-14 (and older)
+#ifndef DD4HEP_VERSION_GE
+#define DD4HEP_VERSION_GE(a,b) 0
+#endif
+
+static Ref_t create_detector(Detector& theDetector, xml_h element, SensitiveDetector sens) {
+ static double tolerance = 0e0;
+
+ cout << "---------------------------------" << endl;
+ cout << " creating Ecal ECRing ( SEcal05_ECRing ) " << endl;
+ cout << "---------------------------------" << endl;
+
+ xml_det_t x_det = element;
+ string det_name = x_det.nameStr();
+ Layering layering (element);
+
+ Material air = theDetector.air();
+ Material CF = theDetector.material("CarbonFiber");
+
+ int det_id = x_det.id();
+ xml_comp_t x_staves = x_det.staves();
+ DetElement sdet (det_name,det_id);
+
+ // --- create an envelope volume and position it into the world ---------------------
+
+ Volume envelope = dd4hep::xml::createPlacedEnvelope( theDetector, element , sdet ) ;
+
+ dd4hep::xml::setDetectorTypeFlag( element, sdet ) ;
+
+ if( theDetector.buildType() == BUILD_ENVELOPE ) return sdet ;
+
+ //-----------------------------------------------------------------------------------
+
+ sens.setType("calorimeter");
+
+ Material stave_material = theDetector.material(x_staves.materialStr());
+
+ Readout readout = sens.readout();
+ Segmentation seg = readout.segmentation();
+
+ std::vector<double> cellSizeVector = seg.segmentation()->cellDimensions(0); //Assume uniform cell sizes, provide dummy cellID
+ double cell_sizeX = cellSizeVector[0];
+ double cell_sizeY = cellSizeVector[1];
+
+ //====================================================================
+ //
+ // Read all the constant from ILD_o1_v05.xml
+ // Use them to build Ecal ECRing
+ //
+ //====================================================================
+
+ // read parametere from compact.xml file
+ double Ecal_fiber_thickness_alveolus = theDetector.constant<double>("Ecal_fiber_thickness_alveolus");
+ double Ecal_fiber_thickness_structure = theDetector.constant<double>("Ecal_fiber_thickness_structure");
+
+ double Ecal_radiator_thickness1 = theDetector.constant<double>("Ecal_radiator_layers_set1_thickness");
+ double Ecal_radiator_thickness2 = theDetector.constant<double>("Ecal_radiator_layers_set2_thickness");
+ double Ecal_radiator_thickness3 = theDetector.constant<double>("Ecal_radiator_layers_set3_thickness");
+ double Ecal_Barrel_halfZ = theDetector.constant<double>("Ecal_Barrel_halfZ");
+
+ double Ecal_support_thickness = theDetector.constant<double>("Ecal_support_thickness");
+ double Ecal_front_face_thickness = theDetector.constant<double>("Ecal_front_face_thickness");
+ double Ecal_lateral_face_thickness = theDetector.constant<double>("Ecal_lateral_face_thickness");
+
+ double Ecal_EC_Ring_gap = theDetector.constant<double>("Ecal_EC_Ring_gap");
+
+ double Ecal_cables_gap = theDetector.constant<double>("Ecal_cables_gap");
+ double Lcal_outer_radius = theDetector.constant<double>("Lcal_outer_radius");
+ double Ecal_Lcal_ring_gap = theDetector.constant<double>("Ecal_Lcal_ring_gap");
+ double Ecal_endcap_center_box_size = theDetector.constant<double>("Ecal_endcap_center_box_size");
+
+ int Ecal_nlayers1 = theDetector.constant<int>("Ecal_nlayers1");
+ int Ecal_nlayers2 = theDetector.constant<int>("Ecal_nlayers2");
+ int Ecal_nlayers3 = theDetector.constant<int>("Ecal_nlayers3");
+
+ double EcalEndcapRing_inner_radius = theDetector.constant<double>("EcalEndcapRing_inner_radius");
+ double EcalEndcapRing_outer_radius = theDetector.constant<double>("EcalEndcapRing_outer_radius");
+ double EcalEndcapRing_min_z = theDetector.constant<double>("EcalEndcapRing_min_z");
+ double EcalEndcapRing_max_z = theDetector.constant<double>("EcalEndcapRing_max_z");
+
+
+ double crossing_angle = theDetector.constant<double>("CepC_Main_Crossing_Angle");
+ bool Ecal_Barrel_PreshowerLayer = theDetector.constant<int>("Ecal_Barrel_Preshower") > 0;
+
+ double Ecal_barrel_thickness = theDetector.constant<double>("Ecal_barrel_thickness");
+
+ //========== fill data for reconstruction ============================
+ dd4hep::rec::LayeredCalorimeterData* caloData = new dd4hep::rec::LayeredCalorimeterData ;
+ caloData->layoutType = dd4hep::rec::LayeredCalorimeterData::EndcapLayout ;
+ caloData->inner_symmetry = 0 ; // hard code cernter pipe hole
+ caloData->outer_symmetry = 4 ; // outer box
+ caloData->phi0 = 0 ;
+
+ /// extent of the calorimeter in the r-z-plane [ rmin, rmax, zmin, zmax ] in mm.
+ caloData->extent[0] = EcalEndcapRing_inner_radius ;
+ caloData->extent[1] = EcalEndcapRing_outer_radius ;
+ caloData->extent[2] = EcalEndcapRing_min_z ;
+ caloData->extent[3] = EcalEndcapRing_max_z ;
+
+
+ //====================================================================
+ //
+ // general calculated parameters
+ //
+ //====================================================================
+
+ int Number_of_Si_Layers_in_Barrel = 0;
+
+#ifdef VERBOSE
+ std::cout << " Ecal total number of Silicon layers = " << Number_of_Si_Layers_in_Barrel << std::endl;
+#endif
+
+ int n_total_layers = Ecal_nlayers1 + Ecal_nlayers2 + Ecal_nlayers3;
+ Number_of_Si_Layers_in_Barrel = n_total_layers; // take account of preshower on/off (djeans)
+ if ( Ecal_Barrel_PreshowerLayer ) Number_of_Si_Layers_in_Barrel++;
+
+ // calculate total thickness of ECAL
+ // updated to take info from xml description - djeans: 12 July 2017
+ double module_thickness = Ecal_support_thickness + Ecal_front_face_thickness;// front and back supports
+
+ // the absorber in the structure
+ for (int i=0; i<Ecal_nlayers1+Ecal_nlayers2+Ecal_nlayers3; i++) {
+ bool inStructure = Ecal_Barrel_PreshowerLayer==1 ? i%2==1 : i%2==0 ;
+ if ( inStructure ) {
+ double thickness (Ecal_radiator_thickness1);
+ if ( i>=Ecal_nlayers1 ) thickness = Ecal_radiator_thickness2;
+ if ( i>=Ecal_nlayers1+Ecal_nlayers2 ) thickness = Ecal_radiator_thickness3;
+ module_thickness += thickness + 2*Ecal_fiber_thickness_structure; // the absorber and its wrapping
+ }
+ }
+
+ // the slabs
+ int l_num2(0);
+ for(xml_coll_t li(x_det,_U(layer)); li; ++li) { // types of layers (i.e. thin/thick absorber) or "stack"
+ xml_comp_t x_layer = li;
+ // Loop over number of repeats for this layer.
+ for (int j=0; j< x_layer.repeat(); j++) { // layers within this type (or "stack")
+ float thisthick = layering.layer(l_num2)->thickness();
+ module_thickness+=thisthick + 2*Ecal_fiber_thickness_alveolus; // slab thickness, and the alveolar wall around it
+ l_num2++;
+ }
+ }
+
+
+ if ( fabs( Ecal_barrel_thickness - module_thickness) > 0.1 ) {
+ cout << "ERROR EC ecal thickness inconsistency: calculated, nominal = " << module_thickness << " " << Ecal_barrel_thickness << endl;
+ assert(0);
+ }
+
+#ifdef VERBOSE
+ std::cout << " module_thickness = " << module_thickness << std::endl;
+#endif
+
+
+ double ECRingSiplateSize = Ecal_endcap_center_box_size
+ - 2 * Ecal_EC_Ring_gap
+ - 2 * Ecal_lateral_face_thickness;
+
+
+ // central hole is not centred on detector axis, but on outgoing beampipe (as is the lumical)
+ double hole_x_offset = tan(crossing_angle/2.) * (EcalEndcapRing_min_z + EcalEndcapRing_max_z)/2.;
+
+
+ // ========= Create Ecal end cap ring ====================================
+ // It will be the volume for palcing the Ecal endcaps alveolus(i.e. Layers).
+ // And the structure W plate.
+ // Itself will be placed into the world volume.
+ // ==========================================================================
+
+ //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ //~ EndcapRing ~
+ //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+ for(int module_num=0;module_num<2;module_num++) { // move module loop to here (modules not identical)
+
+ DetElement module_det(_toString(module_num,"module%d"),det_id);
+
+ double thismodule_hole_dx = hole_x_offset;
+ if ( module_num==0 ) thismodule_hole_dx*=-1;
+ Position hole_position(thismodule_hole_dx, 0, 0);
+
+ double Ecal_endcap_Tube_rmax = Lcal_outer_radius + Ecal_Lcal_ring_gap;
+
+ Tube CenterECTub(0., Ecal_endcap_Tube_rmax, module_thickness);
+ Box CenterECBox((Ecal_endcap_center_box_size/2. - Ecal_EC_Ring_gap),
+ (Ecal_endcap_center_box_size/2. - Ecal_EC_Ring_gap),
+ module_thickness/2.);
+
+ SubtractionSolid ECRingSolid( CenterECBox, CenterECTub, hole_position);
+
+ // Volume EnvLogECRing("ECRing",ECRingSolid,theDetector.material("g10"));
+ Volume EnvLogECRing("ECRing",ECRingSolid,theDetector.material("CarbonFiber")); // DJeans 5-sep-2016
+
+ //==============================================================
+ // build the layer and place into the ECRing
+ //==============================================================
+
+
+ //-------------------------------------------------------
+ // Radiator and towers placements inside the ECRing module
+ //-------------------------------------------------------
+
+ // We count the layers starting from IP and from 1,
+ // so odd layers should be inside slabs and
+ // even ones on the structure.
+ //
+
+
+ double distance_ip_to_front_face = Ecal_Barrel_halfZ + Ecal_cables_gap;
+ double dist_from_front_face_tally(0);
+
+ double z_floor = - module_thickness/2;
+
+ double l_pos_z = z_floor;
+
+ dd4hep::rec::LayeredCalorimeterData::Layer caloLayer ;
+ caloLayer.cellSize0 = cell_sizeX;
+ caloLayer.cellSize1 = cell_sizeY;
+
+ //-------------------- start loop over ECAL layers ----------------------
+ // Loop over the sets of layer elements in the detector.
+ double radiator_dim_y = Ecal_radiator_thickness1; //to be updated with slice radiator thickness
+
+ double nRadiationLengths=0.;
+ double nInteractionLengths=0.;
+ double thickness_sum=0;
+
+ int l_num = 1;
+ bool isFirstSens = true;
+ int myLayerNum = 0 ;
+ bool isFrontFace=true;
+ int absorber_index(0);
+
+ for(xml_coll_t li(x_det,_U(layer)); li; ++li) { // types of layer
+ xml_comp_t x_layer = li;
+ int repeat = x_layer.repeat();
+
+ // Loop over number of repeats for this layer.
+ for (int j=0; j<repeat; j++) {
+
+ // move structural layer to before the slabs - djeans
+
+ // deal with preshower or lack thereof ------------------
+ double radiator_dim_Z(0);
+ //double rad_pos_Z(0);
+ double this_struct_CFthick_beforeAbs(0);
+ double this_struct_CFthick_afterAbs(0);
+ double _CF_absWrap = Ecal_fiber_thickness_structure;
+ double _CF_alvWall = Ecal_fiber_thickness_alveolus;
+ if ( isFrontFace ) { // the first part of the module depends on whether we have preshwer or not
+ if ( Ecal_Barrel_PreshowerLayer==1 ) { // don't include W+CF wrapping; only one side of alveolus
+ this_struct_CFthick_beforeAbs = 0;
+ this_struct_CFthick_afterAbs = Ecal_front_face_thickness + _CF_alvWall;
+ radiator_dim_Z = 0;
+ } else { // include W+CF wrapping; only one side of alveolus
+ this_struct_CFthick_beforeAbs = Ecal_front_face_thickness + _CF_absWrap;
+ this_struct_CFthick_afterAbs = _CF_absWrap + _CF_alvWall;
+ radiator_dim_Z = Ecal_radiator_thickness1; // this line implicitly assumes Ecal_nlayers1>0...
+ //rad_pos_Z = radiator_dim_Z/2. + this_struct_CFthick_beforeAbs; // distance from top surface of structure to centre of radiator
+ absorber_index++;
+ }
+
+ thickness_sum += radiator_dim_Z;
+ nRadiationLengths += radiator_dim_Z/stave_material.radLength();
+ nInteractionLengths += radiator_dim_Z/stave_material.intLength();
+
+ thickness_sum += ( this_struct_CFthick_beforeAbs + this_struct_CFthick_afterAbs );
+ nRadiationLengths += ( this_struct_CFthick_beforeAbs + this_struct_CFthick_afterAbs ) / CF.radLength();
+ nInteractionLengths += ( this_struct_CFthick_beforeAbs + this_struct_CFthick_afterAbs ) / CF.intLength();
+
+ isFrontFace=false;
+ } else { // internal layer: include W+CF wrapping; both sides of alveolus
+ this_struct_CFthick_beforeAbs = _CF_alvWall+_CF_absWrap;
+ this_struct_CFthick_afterAbs = _CF_alvWall+_CF_absWrap;
+
+ if ( absorber_index<Ecal_nlayers1 ) radiator_dim_Z=Ecal_radiator_thickness1;
+ else if ( absorber_index<Ecal_nlayers1+Ecal_nlayers2 ) radiator_dim_Z=Ecal_radiator_thickness2;
+ else if ( absorber_index<Ecal_nlayers1+Ecal_nlayers2+Ecal_nlayers3 ) radiator_dim_Z=Ecal_radiator_thickness3;
+ else {
+ assert(0 && "ERROR cannot determine absorber thickness for layer ");
+ }
+
+ absorber_index++;
+
+ assert( radiator_dim_Z>0 && "no radiator!" );
+ }
+ //-------------------------------
+
+ radiator_dim_y = radiator_dim_Z; // ahh different axis conventions....
+
+ l_pos_z += this_struct_CFthick_beforeAbs + radiator_dim_y/2;
+
+ // #########################
+ // BuildECRingStructureLayer
+ // #########################
+
+ Tube CenterECTubForSi(0.,
+ Lcal_outer_radius + Ecal_Lcal_ring_gap + 0.001, // + tolerance,
+ module_thickness,
+ 0.,
+ 2 * M_PI);
+
+
+ string l_name = _toString(l_num,"layer%d");
+ RotationZYX rot(0,0,0); // a null rotation
+
+ if ( radiator_dim_y>0 ) {
+
+ string bs_name="bs";
+ Box EndcapStructureLayer_box(ECRingSiplateSize/ 2. - tolerance, ECRingSiplateSize/ 2. - tolerance, radiator_dim_y/2.);
+ SubtractionSolid EndcapStructureLayerSolid( EndcapStructureLayer_box, CenterECTubForSi, hole_position);
+ Volume EndcapStructureLayer_vol(det_name+"_"+l_name+"_"+bs_name,EndcapStructureLayerSolid,theDetector.material(x_staves.materialStr()));
+ EndcapStructureLayer_vol.setVisAttributes(theDetector.visAttributes( x_staves.visStr()));
+
+ // this is where we place the tungsten into the envelope
+ double bsl_pos_z = l_pos_z;
+ Position bsl_pos(0,0,bsl_pos_z);
+ Transform3D bsl_tran3D(rot,bsl_pos);
+ EnvLogECRing.placeVolume(EndcapStructureLayer_vol,bsl_tran3D);
+ }
+
+ // Increment to next layer Z position.
+ l_pos_z += radiator_dim_y/2 + this_struct_CFthick_afterAbs;
+
+ // now move to the alveolus and what's inside
+ double l_thickness = layering.layer(l_num-1)->thickness(); // Layer's thickness. this is the internal thickness of the alveolus
+
+ l_pos_z += l_thickness/2.; // add in half the alveolus thickness
+
+ int EC_Number_of_towers = 0;
+
+ // We use the same method able to create the Barrel
+ // Slabs, so we have to rotate it later when placing
+ // into the EC modules.
+ //
+ // We use the same method able to create the Barrel
+ // radiator plates between slabs, but with the good
+ // dimensions to avoid to rotate it later when placing
+ // into the EC modules.
+
+ // While the towers have the same shape use the same
+ // logical volumes and parameters.
+
+ string tower_name = _toString(EC_Number_of_towers,"tower%d");
+
+ // build the ECRing layer, a box with hole in the middle
+ Box ECRingSiBox( ECRingSiplateSize/ 2. - tolerance, ECRingSiplateSize/ 2. - tolerance, l_thickness/2.0-tolerance);
+
+ SubtractionSolid ECRingSiSolid( ECRingSiBox, CenterECTubForSi, hole_position);
+
+ Volume l_vol(det_name+"_"+l_name+"_"+tower_name,ECRingSiSolid,air);
+ DetElement layer(module_det, l_name+tower_name, det_id);
+
+ l_vol.setVisAttributes(theDetector.visAttributes( x_layer.visStr() ));
+
+ // Loop over the sublayers or slices for this layer.
+ int s_num = 1;
+ double s_pos_z = -(l_thickness / 2);
+
+ //--------------------------------------------------------------------------------
+ // BuildECRing keep the Alveolus structure, the same as the Barrel and Endcap
+ //--------------------------------------------------------------------------------
+
+ for(xml_coll_t si(x_layer,_U(slice)); si; ++si) {
+ xml_comp_t x_slice = si;
+ string s_name = _toString(s_num,"slice%d");
+ double s_thick = x_slice.thickness();
+ Material slice_material = theDetector.material(x_slice.materialStr());
+
+ double slab_dim_x = ECRingSiplateSize/2.-tolerance;
+ double slab_dim_y = s_thick/2.-tolerance;
+ double slab_dim_z = ECRingSiplateSize/2.-tolerance;
+
+ Box s_box(slab_dim_x,slab_dim_z,slab_dim_y);
+
+ SubtractionSolid ECRingSiSliceSolid( s_box, CenterECTubForSi, hole_position);
+
+ Volume s_vol(det_name+"_"+l_name+"_"+s_name,ECRingSiSliceSolid,slice_material);
+ DetElement slice(layer,s_name,det_id);
+
+ s_vol.setVisAttributes(theDetector.visAttributes(x_slice.visStr()));
+
+#ifdef VERBOSE
+ std::cout<<"x_slice.materialStr(): "<< x_slice.materialStr() <<std::endl;
+#endif
+ if (x_slice.materialStr().compare(x_staves.materialStr()) == 0){
+ radiator_dim_y = s_thick;
+
+ absorber_index++;
+
+#if DD4HEP_VERSION_GE( 0, 15 )
+ caloLayer.outer_nRadiationLengths = nRadiationLengths;
+ caloLayer.outer_nInteractionLengths = nInteractionLengths;
+ caloLayer.outer_thickness = thickness_sum;
+ caloLayer.distance = distance_ip_to_front_face + dist_from_front_face_tally;
+#endif
+ if (Ecal_Barrel_PreshowerLayer==0 || !isFirstSens ){ // add this layer to the caloData (unless its the first absorber layer of a no-preshower calo)
+ if ( module_num==0 ) {
+ caloData->layers.push_back( caloLayer ) ;
+#ifdef VERBOSE
+#if DD4HEP_VERSION_GE( 0, 15 )
+ std::cout<<" caloLayer.distance: "<< caloLayer.distance <<std::endl;
+ std::cout<<" caloLayer.inner_nRadiationLengths: "<< caloLayer.inner_nRadiationLengths <<std::endl;
+ std::cout<<" caloLayer.inner_nInteractionLengths: "<< caloLayer.inner_nInteractionLengths <<std::endl;
+ std::cout<<" caloLayer.inner_thickness: "<< caloLayer.inner_thickness <<std::endl;
+ std::cout<<" caloLayer.sensitive_thickness: "<< caloLayer.sensitive_thickness <<std::endl;
+ std::cout<<" caloLayer.cellSize0, 1: " << caloLayer.cellSize0 << " " << caloLayer.cellSize1 << std::endl;
+ std::cout<<" caloLayer.outer_nRadiationLengths: "<< caloLayer.outer_nRadiationLengths <<std::endl;
+ std::cout<<" caloLayer.outer_nInteractionLengths: "<< caloLayer.outer_nInteractionLengths <<std::endl;
+ std::cout<<" caloLayer.outer_thickness: "<< caloLayer.outer_thickness <<std::endl;
+ std::cout<<" EcalECRing[1]==>caloLayer.inner_thickness + caloLayer.outer_thickness: "
+ << caloLayer.inner_thickness + caloLayer.outer_thickness <<std::endl;
+#endif
+#endif
+ dist_from_front_face_tally += caloLayer.inner_thickness+caloLayer.outer_thickness;
+ }
+ }
+
+ nRadiationLengths = 0. ;
+ nInteractionLengths = 0. ;
+ thickness_sum = 0. ;
+ isFirstSens = false;
+ } // if radiator
+
+ nRadiationLengths += s_thick/(2.*slice_material.radLength());
+ nInteractionLengths += s_thick/(2.*slice_material.intLength());
+ thickness_sum += s_thick/2;
+
+ if ( x_slice.isSensitive() ) {
+ s_vol.setSensitiveDetector(sens);
+
+#if DD4HEP_VERSION_GE( 0, 15 )
+ //Store "inner" quantities
+ caloLayer.inner_nRadiationLengths = nRadiationLengths;
+ caloLayer.inner_nInteractionLengths = nInteractionLengths;
+ caloLayer.inner_thickness = thickness_sum;
+ //Store sensitive slice thickness
+ caloLayer.sensitive_thickness = s_thick;
+#ifdef VERBOSE
+ std::cout<<" l_num: "<<l_num <<std::endl;
+ std::cout<<" s_num: "<<s_num <<std::endl;
+ std::cout<<" module_thickness: "<< module_thickness <<std::endl;
+ std::cout<<" l_pos_z: "<< l_pos_z <<std::endl;
+ std::cout<<" l_thickness: "<< l_thickness <<std::endl;
+ std::cout<<" s_pos_z: "<< s_pos_z <<std::endl;
+ std::cout<<" s_thick: "<< s_thick <<std::endl;
+ std::cout<<" radiator_dim_y: "<< radiator_dim_y <<std::endl;
+#endif
+ caloLayer.absorberThickness = radiator_dim_y ;
+#endif
+ nRadiationLengths = 0. ;
+ nInteractionLengths = 0. ;
+ thickness_sum = 0. ;
+ } // if sensitive
+
+ nRadiationLengths += s_thick/(2.*slice_material.radLength());
+ nInteractionLengths += s_thick/(2.*slice_material.intLength());
+ thickness_sum += s_thick/2;
+
+ slice.setAttributes(theDetector,s_vol,x_slice.regionStr(),x_slice.limitsStr(),x_slice.visStr());
+
+ // Slice placement.
+ PlacedVolume slice_phv = l_vol.placeVolume(s_vol,Position(0,0,s_pos_z+s_thick/2));
+
+ if ( x_slice.isSensitive() ) {
+ slice_phv.addPhysVolID("layer", myLayerNum++);
+ }
+
+ slice.setPlacement(slice_phv);
+ // Increment Z position of slice.
+ s_pos_z += s_thick;
+
+ // Increment slice number.
+ ++s_num;
+ }
+
+ // add material between the "slabs"
+
+ // we need to use correct thickness here! especially at the interface between stacks.
+ if ( absorber_index < Ecal_nlayers1 ) radiator_dim_y=Ecal_radiator_thickness1;
+ else if ( absorber_index < Ecal_nlayers1+Ecal_nlayers2 ) radiator_dim_y=Ecal_radiator_thickness2;
+ else if ( absorber_index < Ecal_nlayers1+Ecal_nlayers2+Ecal_nlayers3 ) radiator_dim_y=Ecal_radiator_thickness3;
+ else {
+ radiator_dim_y=0;
+ }
+
+ // deal with final support plate
+ float cf_thick = Ecal_fiber_thickness_alveolus;
+ if ( absorber_index<Ecal_nlayers1+Ecal_nlayers2+Ecal_nlayers3 ) {
+ cf_thick += Ecal_fiber_thickness_structure;
+ } else {
+ cf_thick += Ecal_support_thickness;
+ }
+
+ if ( module_num==0 ) {
+#if DD4HEP_VERSION_GE( 0, 15 )
+ caloLayer.distance = distance_ip_to_front_face + dist_from_front_face_tally;
+ caloLayer.outer_nRadiationLengths = nRadiationLengths + cf_thick/CF.radLength();
+ caloLayer.outer_nInteractionLengths = nInteractionLengths + cf_thick/CF.intLength();
+ caloLayer.outer_thickness = thickness_sum + cf_thick;
+#endif
+ caloData->layers.push_back( caloLayer ) ;
+
+#ifdef VERBOSE
+#if DD4HEP_VERSION_GE( 0, 15 )
+ std::cout<<" caloLayer.distance: "<< caloLayer.distance <<std::endl;
+ std::cout<<" caloLayer.inner_nRadiationLengths: "<< caloLayer.inner_nRadiationLengths <<std::endl;
+ std::cout<<" caloLayer.inner_nInteractionLengths: "<< caloLayer.inner_nInteractionLengths <<std::endl;
+ std::cout<<" caloLayer.inner_thickness: "<< caloLayer.inner_thickness <<std::endl;
+ std::cout<<" caloLayer.sensitive_thickness: "<< caloLayer.sensitive_thickness <<std::endl;
+ std::cout<<" caloLayer.cellSize0, 1: " << caloLayer.cellSize0 << " " << caloLayer.cellSize1 << std::endl;
+ std::cout<<" caloLayer.outer_nRadiationLengths: "<< caloLayer.outer_nRadiationLengths <<std::endl;
+ std::cout<<" caloLayer.outer_nInteractionLengths: "<< caloLayer.outer_nInteractionLengths <<std::endl;
+ std::cout<<" caloLayer.outer_thickness: "<< caloLayer.outer_thickness <<std::endl;
+
+ std::cout<<" EcalECRing[2]==>caloLayer.inner_thickness + caloLayer.outer_thickness: "
+ << caloLayer.inner_thickness + caloLayer.outer_thickness <<std::endl;
+#endif
+#endif
+ dist_from_front_face_tally += caloLayer.inner_thickness+caloLayer.outer_thickness;
+ }
+
+ nRadiationLengths = radiator_dim_y/stave_material.radLength() + cf_thick/CF.radLength();
+ nInteractionLengths = radiator_dim_y/stave_material.intLength() + cf_thick/CF.intLength();
+ thickness_sum = radiator_dim_y + cf_thick;
+
+ // this is where we place the slab (l_vol) into the envelope
+ int i_stave = 1;
+ int i_tower = 1;
+ Position l_pos(0,0,l_pos_z);
+ Transform3D tran3D(rot,l_pos);
+ PlacedVolume layer_phv = EnvLogECRing.placeVolume(l_vol,tran3D);
+ layer_phv.addPhysVolID("tower", i_tower);
+ layer_phv.addPhysVolID("stave", i_stave);
+ layer.setPlacement(layer_phv);
+
+ l_pos_z += l_thickness/2.; // add in the other half of the alveolus
+
+ ++l_num;
+
+ }
+ }
+
+
+ // Set stave visualization.
+ if (x_staves) {
+ EnvLogECRing.setVisAttributes(theDetector.visAttributes(x_staves.visStr()));
+ }
+
+
+ //====================================================================
+ // Place Ecal Endcap module into the assembly envelope volume
+ //====================================================================
+
+ double EC_module_z_offset = (EcalEndcapRing_min_z + EcalEndcapRing_max_z)/2.;
+
+ int module_id = ( module_num == 0 ) ? 0:6;
+ double this_module_z_offset = ( module_id == 0 ) ? - EC_module_z_offset : EC_module_z_offset;
+ double this_module_rotY = ( module_id == 0 ) ? M_PI:0;
+
+ Position xyzVec(0,0,this_module_z_offset);
+ RotationZYX rot(0,this_module_rotY,0);
+ Rotation3D rot3D(rot);
+ Transform3D tran3D(rot3D,xyzVec);
+
+ PlacedVolume pv = envelope.placeVolume(EnvLogECRing,tran3D);
+ pv.addPhysVolID("module",module_id); // z: -/+ 0/6
+
+ DetElement sd = module_det;
+ sd.setPlacement(pv);
+
+ }
+
+ sdet.addExtension< dd4hep::rec::LayeredCalorimeterData >( caloData ) ;
+
+ return sdet;
+
+}
+
+
+
+DECLARE_DETELEMENT(SEcal05_ECRing, create_detector)
+
diff --git a/Detector/DetCEPCv4/src/calorimeter/SEcal05_Endcaps.cpp b/Detector/DetCEPCv4/src/calorimeter/SEcal05_Endcaps.cpp
new file mode 100644
index 00000000..a0f4c69a
--- /dev/null
+++ b/Detector/DetCEPCv4/src/calorimeter/SEcal05_Endcaps.cpp
@@ -0,0 +1,390 @@
+//====================================================================
+// lcgeo - LC detector models in DD4hep
+//--------------------------------------------------------------------
+// DD4hep Geometry driver for SiWEcalEndcaps
+// Ported from Mokka
+//--------------------------------------------------------------------
+// S.Lu, DESY
+// $Id: SEcal04_Endcaps.cpp 1060 2016-09-05 07:48:41Z /C=JP/O=KEK/OU=CRC/CN=JEANS Daniel Thomelin Dietrich $
+//====================================================================
+
+ /* History:
+
+// *******************************************************
+// * *
+// * Mokka *
+// * - the detailed geant4 simulation for ILC - *
+// * *
+// * For more information about Mokka, visit the *
+// * *
+// * Mokka.in2p3.fr Mokka home page. *
+// * *
+// *******************************************************
+//
+// $Id: SEcal04_Endcaps.cpp 1060 2016-09-05 07:48:41Z /C=JP/O=KEK/OU=CRC/CN=JEANS Daniel Thomelin Dietrich $
+// $Name: mokka-07-00 $
+//
+//
+//
+// SEcal04.cc
+//
+ Shaojun Lu: Ported from Mokka SEcal04 Endcaps part. Read the constants from XML
+ instead of the DB. Then build the Endcap in the same way with DD4hep
+ construct.
+ Inside SEcal04, some parameters, which used by Ecal Endcaps, come from
+ Ecal Barrel. They can be ssen here again.
+ */
+
+#include "DD4hep/DetFactoryHelper.h"
+#include "XML/Layering.h"
+#include "DD4hep/Shapes.h"
+#include "XML/Utilities.h"
+#include "DDRec/DetectorData.h"
+#include "DDSegmentation/WaferGridXY.h"
+#include <sstream>
+
+using namespace std;
+
+using dd4hep::BUILD_ENVELOPE;
+using dd4hep::Box;
+using dd4hep::DetElement;
+using dd4hep::Detector;
+using dd4hep::IntersectionSolid;
+using dd4hep::Layering;
+using dd4hep::Material;
+using dd4hep::PlacedVolume;
+using dd4hep::PolyhedraRegular;
+using dd4hep::Position;
+using dd4hep::Readout;
+using dd4hep::Ref_t;
+using dd4hep::Rotation3D;
+using dd4hep::RotationZYX;
+using dd4hep::Segmentation;
+using dd4hep::SensitiveDetector;
+using dd4hep::Transform3D;
+using dd4hep::Volume;
+using dd4hep::_toString;
+
+using dd4hep::rec::LayeredCalorimeterData;
+
+#include "SEcal05_Helpers.h"
+
+#undef NDEBUG
+#include <assert.h>
+
+
+//#define VERBOSE 1
+
+// workaround for DD4hep v00-14 (and older)
+#ifndef DD4HEP_VERSION_GE
+#define DD4HEP_VERSION_GE(a,b) 0
+#endif
+
+static Ref_t create_detector(Detector& theDetector, xml_h element, SensitiveDetector sens) {
+
+ cout << "------------------------" << endl;
+ cout << "creating SEcal05_Endcaps" << endl;
+ cout << "------------------------" << endl;
+
+ xml_det_t x_det = element;
+ string det_name = x_det.nameStr();
+ Layering layering (element);
+
+ int det_id = x_det.id();
+ // xml_comp_t x_staves = x_det.staves();
+ DetElement sdet (det_name,det_id);
+ // Volume motherVol = theDetector.pickMotherVolume(sdet);
+
+ // --- create an envelope volume and position it into the world ---------------------
+
+ Volume envelope = dd4hep::xml::createPlacedEnvelope( theDetector, element , sdet ) ;
+
+ dd4hep::xml::setDetectorTypeFlag( element, sdet ) ;
+
+ if( theDetector.buildType() == BUILD_ENVELOPE ) return sdet ;
+ //-----------------------------------------------------------------------------------
+
+ sens.setType("calorimeter");
+
+ Readout readout = sens.readout();
+ Segmentation seg = readout.segmentation();
+
+ //====================================================================
+ //
+ // Read all the constant from ILD_o1_v05.xml
+ // Use them to build Ecal endcaps
+ //
+ //====================================================================
+
+ // read parameters from compact.xml file
+ double Ecal_Slab_shielding = theDetector.constant<double>("Ecal_Slab_shielding");
+ double Ecal_fiber_thickness_structure = theDetector.constant<double>("Ecal_fiber_thickness_structure"); // absorber wrapping thickness
+ double Ecal_fiber_thickness_alveolus = theDetector.constant<double>("Ecal_fiber_thickness_alveolus"); // alveolar wall thickness
+
+ double EcalBarrel_inner_radius = theDetector.constant<double>("TPC_outer_radius") +theDetector.constant<double>("Ecal_Tpc_gap");
+ int Ecal_barrel_z_modules = theDetector.constant<int>("Ecal_barrel_z_modules");
+
+ double Ecal_radiator_thickness1 = theDetector.constant<double>("Ecal_radiator_layers_set1_thickness");
+ double Ecal_radiator_thickness2 = theDetector.constant<double>("Ecal_radiator_layers_set2_thickness");
+ double Ecal_radiator_thickness3 = theDetector.constant<double>("Ecal_radiator_layers_set3_thickness");
+ double Ecal_Barrel_halfZ = theDetector.constant<double>("Ecal_Barrel_halfZ");
+
+ int Ecal_barrel_number_of_towers = theDetector.constant<int>("Ecal_barrel_number_of_towers");
+
+ double Ecal_support_thickness = theDetector.constant<double>("Ecal_support_thickness");
+ double Ecal_front_face_thickness = theDetector.constant<double>("Ecal_front_face_thickness");
+ double Ecal_lateral_face_thickness = theDetector.constant<double>("Ecal_lateral_face_thickness");
+ double Ecal_Slab_H_fiber_thickness = theDetector.constant<double>("Ecal_Slab_H_fiber_thickness");
+
+ double Ecal_endcap_extra_size = theDetector.constant<double>("Ecal_endcap_extra_size");
+ double Ecal_cables_gap = theDetector.constant<double>("Ecal_cables_gap");
+
+ int Ecal_nlayers1 = theDetector.constant<int>("Ecal_nlayers1");
+ int Ecal_nlayers2 = theDetector.constant<int>("Ecal_nlayers2");
+ int Ecal_nlayers3 = theDetector.constant<int>("Ecal_nlayers3");
+
+ // first layer is preshower?
+ bool Ecal_Endcap_PreshowerLayer = theDetector.constant<int>("Ecal_Endcap_Preshower") > 0;
+
+ std::string Ecal_endcap_number_of_towers = theDetector.constant<string>("Ecal_endcap_number_of_towers");
+
+ int Ecal_end_of_slab_strategy = theDetector.constant<int>("Ecal_end_of_slab_strategy");
+
+
+ int Ecal_n_wafers_per_tower = theDetector.constant<int>("Ecal_n_wafers_per_tower");
+
+ double Ecal_guard_ring_size = theDetector.constant<double>("Ecal_guard_ring_size");
+
+ double EcalEndcap_inner_radius = theDetector.constant<double>("EcalEndcap_inner_radius");
+ double EcalEndcap_outer_radius = theDetector.constant<double>("EcalEndcap_outer_radius");
+ double EcalEndcap_min_z = theDetector.constant<double>("EcalEndcap_min_z");
+ double EcalEndcap_max_z = theDetector.constant<double>("EcalEndcap_max_z");
+
+ std::string Ecal_layerConfig = theDetector.constant<string>("Ecal_layer_pattern");
+
+ int Ecal_cells_across_megatile = theDetector.constant <int> ("Ecal_cells_across_megatile" );
+ int Ecal_strips_across_megatile= theDetector.constant <int> ("Ecal_strips_across_megatile");
+ int Ecal_strips_along_megatile = theDetector.constant <int> ("Ecal_strips_along_megatile" );
+
+ double Ecal_barrel_thickness = theDetector.constant<double>("Ecal_barrel_thickness"); // what's assumed in the compact description
+
+
+ //====================================================================
+ //
+ // set up the helper
+ //
+ //====================================================================
+
+ SEcal05_Helpers helper;
+
+ helper.setDet( & x_det );
+
+ // layer configuration
+ helper.setPreshower( Ecal_Endcap_PreshowerLayer );
+
+ cout << "Preshower ? " << Ecal_Endcap_PreshowerLayer << endl;
+
+ helper.setAbsLayers( Ecal_nlayers1, Ecal_radiator_thickness1,
+ Ecal_nlayers2, Ecal_radiator_thickness2,
+ Ecal_nlayers3, Ecal_radiator_thickness3 );
+
+ cout << "absorber layers " <<
+ Ecal_nlayers1 << "*" << Ecal_radiator_thickness1 << "mm + " <<
+ Ecal_nlayers2 << "*" << Ecal_radiator_thickness2 << "mm + " <<
+ Ecal_nlayers3 << "*" << Ecal_radiator_thickness3 << "mm " << endl;
+
+ // check this setup is self-consistent
+ helper.checkLayerConsistency();
+
+ // set structural CF thicknesses
+ helper.setCFthickness( Ecal_fiber_thickness_structure,
+ Ecal_fiber_thickness_alveolus,
+ Ecal_front_face_thickness,
+ Ecal_support_thickness);
+
+ // check that resulting total thickness is consistent with what's in the compact file
+ // n.b. this assumes same thickness in barrel and endcap!
+ float module_thickness = helper.getTotalThickness();
+ if ( fabs( Ecal_barrel_thickness - module_thickness ) > 0.1*dd4hep::mm ) {
+ cout << "ERROR : thickness in comapct decription not consistent with what I calculate!" << endl;
+ cout << " calculated = " << module_thickness << " compact description: " << Ecal_barrel_thickness << endl;
+ assert( 0 ); // exit
+ }
+
+ cout << "module thickness = " << module_thickness << endl;
+
+ // set up the sensitive layer segmentation
+ helper.setSegmentation( &seg );
+
+ helper.setNCells( Ecal_cells_across_megatile, Ecal_strips_across_megatile, Ecal_strips_along_megatile);
+ helper.setMagicMegatileStrategy ( Ecal_end_of_slab_strategy );
+
+
+ int ntemp;
+ std::stringstream stream(Ecal_layerConfig);
+ std::vector < int > layerConfig;
+ while ( stream >> ntemp ) {
+ assert( ntemp>=0 );
+ layerConfig.push_back( ntemp );
+ }
+ cout << "layer config: ";
+ for (size_t i=0; i<layerConfig.size(); i++) cout << layerConfig[i] << " ";
+ cout << endl;
+
+ helper.setLayerConfig( layerConfig );
+
+ // set the number of towers/modules
+ std::vector < int > ntowers;
+ std::stringstream stream2( Ecal_endcap_number_of_towers );
+ while ( stream2 >> ntemp ) {
+ ntowers.push_back( ntemp );
+ }
+
+ cout << "ECAL endcap tower configuration " << Ecal_endcap_number_of_towers << " : ";
+ for (size_t i=0; i<ntowers.size(); i++)
+ cout << ntowers[i] << " ";
+ cout << endl;
+
+ // check that resulting quadrant size is consistent with compact description
+ // here we assume that the width of an alveolus in the endcaps is the same as that in the barrel.
+ double barrel_alv_width = ( Ecal_Barrel_halfZ*2./Ecal_barrel_z_modules - 2.*Ecal_lateral_face_thickness ) / Ecal_barrel_number_of_towers;
+ double calc_endcap_rout(EcalEndcap_inner_radius);
+ for (size_t i=0; i<ntowers.size(); i++) {
+ calc_endcap_rout+=ntowers[i]*barrel_alv_width + 2.*Ecal_lateral_face_thickness;
+ }
+
+ cout << "alveolus width = " << barrel_alv_width << endl;
+
+ // compare calculated size vs that in compact description
+ if ( fabs( calc_endcap_rout - EcalEndcap_outer_radius ) > 0.1*dd4hep::mm ) {
+ cout << "WARNING, inconsistent ECAL endcap radial extent! calculated: " << calc_endcap_rout <<
+ " cm , nominal (from compact descrition): " << EcalEndcap_outer_radius << endl;
+ cout << "consider changing Ecal_endcap_extra_size from " << Ecal_endcap_extra_size <<
+ " to " << calc_endcap_rout - EcalBarrel_inner_radius - Ecal_barrel_thickness << endl;
+ assert(0);
+ }
+
+
+ helper.setTowersUnits( ntowers,
+ barrel_alv_width,
+ Ecal_n_wafers_per_tower,
+ Ecal_lateral_face_thickness,
+ Ecal_fiber_thickness_alveolus + Ecal_Slab_H_fiber_thickness + Ecal_Slab_shielding,
+ Ecal_guard_ring_size );
+
+
+
+ // ========= Create Ecal endcaps ====================================
+ // It will be the volume for palcing the Ecal endcaps alveolus(i.e. Layers).
+ // And the structure W plate.
+ // Itself will be placed into the world volume.
+ // ==========================================================================
+
+ //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ //~ EndcapStandardModule ~
+ //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+ // octagon
+ PolyhedraRegular ECPolyHedra(8, M_PI/8., 0., calc_endcap_rout, module_thickness);
+
+ // take just one quadrant of the octagon, and remove the centre box
+ double quadr = calc_endcap_rout; // anything which is big enough
+ Box quadrant( quadr, quadr , module_thickness);
+ IntersectionSolid EndCapSolid( ECPolyHedra,
+ quadrant,
+ Position( quadr - EcalEndcap_inner_radius,
+ quadr + EcalEndcap_inner_radius, 0 ) );
+
+ Volume EnvLogEndCap("EcalEndcapQuadrant",EndCapSolid,theDetector.material("CarbonFiber"));
+
+ // kink position wrt bottom of quadrant (inside the lateral support)
+ // Y==0 is defined as outer edge of lateral face of inner module of quadrant
+ double kink_y = calc_endcap_rout*tan(M_PI/8.) - EcalEndcap_inner_radius;
+
+ helper.setModuleDimensions(1, // xytype
+ 0, // xztype
+ calc_endcap_rout + EcalEndcap_inner_radius, // dxMax
+ kink_y
+ );
+
+ helper.setTranslation ( Position ( -EcalEndcap_inner_radius , EcalEndcap_inner_radius, -module_thickness/2. ) );
+
+ // make the module
+
+ LayeredCalorimeterData* caloData = new LayeredCalorimeterData ;
+ caloData->layoutType = LayeredCalorimeterData::EndcapLayout ;
+
+ DetElement mod_det ("quad0",det_id);
+
+ helper.makeModule( EnvLogEndCap,
+ mod_det,
+ *caloData,
+ theDetector,
+ sens );
+
+ for (size_t i=0; i<caloData->layers.size(); i++) {
+ caloData->layers[i].distance += Ecal_Barrel_halfZ + Ecal_cables_gap; // add IP->front face distance
+ }
+
+ // cout << "cell sizes: " << endl;
+ // for (size_t i=0; i<caloData->layers.size(); i++) {
+ // cout << "sensitive layer " << i << " : x,y = " << caloData->layers[i].cellSize0 << " " << caloData->layers[i].cellSize1 << endl;
+ // }
+
+ caloData->layoutType = LayeredCalorimeterData::EndcapLayout ;
+ caloData->inner_symmetry = 4 ; // hard code cernter box hole
+ caloData->outer_symmetry = 8 ; // outer Octagun
+ caloData->phi0 = 0 ;
+
+ /// extent of the calorimeter in the r-z-plane [ rmin, rmax, zmin, zmax ] in mm.
+ caloData->extent[0] = EcalEndcap_inner_radius ;
+ caloData->extent[1] = EcalEndcap_outer_radius ;
+ caloData->extent[2] = EcalEndcap_min_z ;
+ caloData->extent[3] = EcalEndcap_max_z ;
+
+
+
+ //====================================================================
+ // Place Ecal Endcap modules into the assembly envelope volume
+ //====================================================================
+
+ for (int iend=0; iend<2; iend++) { // the 2 endcaps
+ int module_id = ( iend == 0 ) ? 0:6;
+
+ double this_module_z_offset = (EcalEndcap_min_z + EcalEndcap_max_z)/2.;
+ if ( iend == 0 ) this_module_z_offset*=-1;
+
+ double this_module_rotY = iend==0 ? M_PI : 0;
+ double rotZ_offset = iend==0 ? M_PI/8. - M_PI/2. : M_PI/8. + 3.*M_PI/4.; // magic rotation to get modules in right place
+ for (int iquad=0; iquad<4; iquad++) { // the 4 quadrants per endcap
+ int stave_id=iquad+1;
+ double this_module_rotZ(0);
+ if ( iend==0 ) {
+ this_module_rotZ = rotZ_offset - (iquad-2) * M_PI/2.;
+ } else {
+ this_module_rotZ = rotZ_offset + (iquad+1) * M_PI/2.;
+ }
+ Position xyzVec(0,0,this_module_z_offset);
+ RotationZYX rot( this_module_rotZ ,this_module_rotY, 0);
+ Rotation3D rot3D(rot);
+ Transform3D tran3D(rot3D,xyzVec);
+ PlacedVolume pv = envelope.placeVolume(EnvLogEndCap,tran3D);
+ pv.addPhysVolID("module",module_id); // z: -/+ 0/6
+ pv.addPhysVolID("stave",stave_id);
+ DetElement sd = (iend==0 && iquad==0) ? mod_det : mod_det.clone(_toString(module_id,"module%d")+_toString(stave_id,"stave%d"));
+ sd.setPlacement(pv);
+ }
+ }
+
+ sdet.addExtension< LayeredCalorimeterData >( caloData ) ;
+
+ // cout << "finished SEcal05_Endcaps" << endl;
+
+ return sdet;
+
+}
+
+
+
+DECLARE_DETELEMENT(SEcal05_Endcaps, create_detector)
+
diff --git a/Detector/DetCEPCv4/src/calorimeter/SEcal05_Helpers.cpp b/Detector/DetCEPCv4/src/calorimeter/SEcal05_Helpers.cpp
new file mode 100644
index 00000000..f5b57dcf
--- /dev/null
+++ b/Detector/DetCEPCv4/src/calorimeter/SEcal05_Helpers.cpp
@@ -0,0 +1,729 @@
+#include "SEcal05_Helpers.h"
+
+SEcal05_Helpers::SEcal05_Helpers() {
+ // constructor, initialise internal variables
+ _x_det=NULL;
+ _layering=NULL;
+ _preshower=-1;
+
+ _CF_absWrap=-1;
+ _CF_alvWall=-1;
+ _CF_front=-1;
+ _CF_back=-1;
+
+ _ntowers.clear();
+ _towerGap=-999;
+ _unitsPerTower=0;
+ _unitDeadEdge=-999;
+
+ _cells_across_megatile=0;
+ _strips_across_megatile=0;
+ _strips_along_megatile=0;
+
+ _constantSlabXYDimensions.clear();
+
+ _caloLayer.absorberThickness = -999;
+ _caloLayer.sensitive_thickness = -999;
+ _caloLayer.distance = 0;
+
+ _totThick=0;
+
+ _plugLength=0;
+
+ _magicMegatileStrategy=-1;
+}
+
+
+void SEcal05_Helpers::setAbsLayers( int nl1, double th1, int nl2, double th2, int nl3, double th3 ) {
+ // set the number and thicknesses of absorber layers
+ _nlayers1=nl1;
+ _nlayers2=nl2;
+ _nlayers3=nl3;
+ _radiator_thickness1=th1;
+ _radiator_thickness2=th2;
+ _radiator_thickness3=th3;
+ return;
+}
+
+void SEcal05_Helpers::checkLayerConsistency() {
+ // check requested number of absorber layers is cinsistent with preshower status
+ // we are constrained to have an even number of sensitive layers [ 2 sens. layers / slab ]
+ assert( (_preshower==0 || _preshower==1) && "_preshower not set" );
+ int n_total_abs_layers = _nlayers1 + _nlayers2 + _nlayers3; // total number of Absorber layers
+ // check that number of requested absober layers is consistent
+ // if we want a preshower layer, total number of Absorber layers should be odd; otherwise even
+ if ( _preshower==1 && n_total_abs_layers%2==0 ) {
+ std::cout << "SEcal05_Helpers ERROR: inconsistent ECAL model !! if you request a preshower layer, the number of absorber layers = _nlayers1 + _nlayers2 + _nlayers3 must be odd" << std::endl;
+ std::cout << " Ecal_PreshowerLayer = " << _preshower << " ; _nlayers1/2/3 = " << _nlayers1 << " " << _nlayers2 << " " << _nlayers3 << std::endl;
+ assert(0);
+ } else if ( _preshower==0 && n_total_abs_layers%2==1 ) {
+ std::cout << "SEcal05_Helpers ERROR: inconsistent ECAL model !! if you request no preshower layer, the number of absorber layers = _nlayers1 + _nlayers2 + _nlayers3 must be even" << std::endl;
+ std::cout << " Ecal_PreshowerLayer = " << _preshower << " ; _nlayers1/2/3 = " << _nlayers1 << " " << _nlayers2 << " " << _nlayers3 << std::endl;
+ assert(0);
+ }
+ return;
+}
+
+float SEcal05_Helpers::getTotalThickness() {
+ // calculate total thickness of ECAL
+ checkLayerConsistency();
+ assert( _x_det && _layering && "_x_det or _layering not set" );
+ assert( (_preshower==0 || _preshower==1) && "_preshower not set" );
+ assert( _CF_absWrap>=0. && _CF_alvWall>=0. && _CF_front>=0. && _CF_back>=0. && "CF thicknesses not set" );
+ float totalThickness = _CF_front+_CF_back;// front and back supports
+
+ // the absorber in the structure
+ for (unsigned int i=0; i<_nlayers1+_nlayers2+_nlayers3; i++) {
+ bool inStructure = _preshower ? i%2==1 : i%2==0 ;
+ if ( inStructure ) {
+ double thickness (_radiator_thickness1);
+ if ( i>=_nlayers1 ) thickness = _radiator_thickness2;
+ if ( i>=_nlayers1+_nlayers2 ) thickness = _radiator_thickness3;
+ totalThickness += thickness + 2*_CF_absWrap; // the absorber and its wrapping
+ }
+ }
+ // the slabs
+ int l_num(0);
+ for(xml_coll_t li(*_x_det,_U(layer)); li; ++li) { // types of layers (i.e. thin/thick absorber) or "stack"
+ xml_comp_t x_layer = li;
+ // Loop over number of repeats for this layer.
+ for (int j=0; j< x_layer.repeat(); j++) { // layers within this type (or "stack")
+ float thisthick = _layering->layer(l_num)->thickness();
+ totalThickness+=thisthick + 2*_CF_alvWall; // slab thickness, and the alveolar wall around it
+ l_num++;
+ }
+ }
+ return totalThickness;
+}
+
+void SEcal05_Helpers::printSEcal05LayerInfo( dd4hep::rec::LayeredCalorimeterData::Layer & caloLayer) {
+ std::cout<<"SEcal05_Helpers === CALOLAYER printout: " << std::endl;
+ std::cout<<" caloLayer.distance: " << caloLayer.distance <<std::endl;
+ std::cout<<" caloLayer.inner_nRadiationLengths: " << caloLayer.inner_nRadiationLengths <<std::endl;
+ std::cout<<" caloLayer.inner_nInteractionLengths: "<< caloLayer.inner_nInteractionLengths <<std::endl;
+ std::cout<<" caloLayer.inner_thickness: " << caloLayer.inner_thickness <<std::endl;
+ std::cout<<" caloLayer.sensitive_thickness: " << caloLayer.sensitive_thickness <<std::endl;
+ std::cout<<" caloLayer.cellSize0, 1: " << caloLayer.cellSize0 << " " << caloLayer.cellSize1 << std::endl;
+ std::cout<<" caloLayer.outer_nRadiationLengths: " << caloLayer.outer_nRadiationLengths <<std::endl;
+ std::cout<<" caloLayer.outer_nInteractionLengths: "<< caloLayer.outer_nInteractionLengths <<std::endl;
+ std::cout<<" caloLayer.outer_thickness: " << caloLayer.outer_thickness <<std::endl;
+ return;
+}
+
+double SEcal05_Helpers::getAbsThickness( unsigned int iAbsLay ) { //, int n1, int n2, int n3, double t1, double t2, double t3) {
+ // get the thickness of a given absorber layer
+ if ( iAbsLay < _nlayers1 ) return _radiator_thickness1;
+ else if ( iAbsLay - _nlayers1 < _nlayers2 ) return _radiator_thickness2;
+ else if ( iAbsLay - _nlayers1 - _nlayers2 < _nlayers3 ) return _radiator_thickness3;
+ assert(0 && "impossible layer number"); // should never get here
+ return -999.;
+}
+
+
+
+std::vector <SEcal05_Helpers::dimposXYStruct> SEcal05_Helpers::getAbsPlateXYDimensions( double ztop ) {
+ // get the dimensions of absorber plates
+ // for trapeziodal case, this depends on the Z position
+ // CF wrap not taken into account, so this width constains both absorber and its wrapping
+
+ if ( _module_XZtype==1 ) assert( ztop>=0 && "getAbsPlateXYDimensions: ztop not specified" ); // must specify Z position for case with non-uniform layers
+
+ std::vector <dimposXYStruct> absorbersheets;
+
+ if ( _module_XYtype == 1 ) { // each slab is different: calculate slab-by-slab (eg in endcap)
+ absorbersheets = getSlabXYDimensions( ztop );
+ } else { // single sheet over all slabs (eg in barrel)
+ dimposXYStruct ss;
+ if ( _module_XZtype==0 ) { // no taper
+ ss.sizeX = _module_dX_max;
+ ss.posX = ss.sizeX/2.;
+ } else { // size varies by laer
+
+ // ss.sizeX = _module_dX_max - 2.*ztop; // this assumes octagon
+
+ ss.sizeX = _module_dX_max - 2.*ztop/tan(_module_angle); // for general shape
+
+ ss.posX = _module_dX_max/2.; // keep it centered
+ }
+ ss.sizeY = _module_dY_total;
+ ss.posY = _module_dY_total/2.;
+ absorbersheets.push_back(ss);
+ }
+ return absorbersheets;
+}
+
+std::vector <SEcal05_Helpers::dimposXYStruct> SEcal05_Helpers::getSlabXYDimensions( double ztop ) {
+ // calculate size and position of slab volumes
+ // size includes slab and dead area around it
+
+ // ztop is Z of upper face of this slab (the face shortest in X)
+ if ( _module_XZtype==1 ) assert( ztop>=0 && "getSlabXYDimensions: ztop not specified" ); // must specify Z position for case with non-uniform layers
+
+ std::vector <dimposXYStruct> layerslabdims;
+ if ( _module_XZtype == 0 && _constantSlabXYDimensions.size()>0 ) { // every layer is the same, and has already been calculated
+ layerslabdims = _constantSlabXYDimensions;
+ } else {
+ dimposXYStruct ss;
+ int totTow(0);
+ for (size_t imod=0; imod<_ntowers.size(); imod++) { // loop over "modules" [nb, for barrel, we only make a single module; for endcaps, typically have 3 per quadrant]
+ for (int itow = 0; itow<_ntowers[imod]; itow++) { // towers within the module
+ // each slab has same width
+ ss.sizeY = _alveolus_total_dim_Y; // this size include edge-of-tower dead space
+ ss.posY =
+ _moduleGap * ( 1 + 2*imod ) + // edge-of-module gaps from Y=0 to this slab
+ _alveolus_total_dim_Y * (totTow+0.5); // position of slab centre w.r.t (X,Y) = (0,0)
+
+ if ( _module_XYtype==0 ) { // all slabs in a layer have same length
+
+ if ( _module_XZtype == 0 ) { // all layers have same length
+ ss.sizeX = _module_dX_max;
+ ss.posX = ss.sizeX/2.;
+ } else { // layers vary in length : barrel module
+
+ // ss.sizeX = _module_dX_max - 2.*ztop; // assumes octagon
+
+ ss.sizeX = _module_dX_max - 2.*ztop/tan(_module_angle); // for general shape
+
+ ss.posX = _module_dX_max/2.; // slabs all centred
+ }
+
+ ss.sizeX -= _plugLength; // DANIELHACK
+ ss.posX += _plugLength/2.; // DANIELHACK
+
+
+
+ } else if ( _module_XYtype==1 ) { // slabs within layer have different lengths : this is for endcap
+ // placing in Y
+ double upperY = ss.posY + 0.5*_alveolus_total_dim_Y; // Y position of upper edge
+
+ if ( upperY<=_module_dY_kink ) { // straight edge part under kink
+ ss.sizeX = _module_dX_max;
+ } else { // sloping edge: take length at upperY, which is the shortest one
+ ss.sizeX = _module_dX_max - ( upperY - _module_dY_kink );
+ }
+
+ ss.posX = ss.sizeX/2.; // this aligns the -X end of slab
+
+ ss.sizeX -= _plugLength; // DANIELHACK
+ ss.posX += _plugLength/2.; // DANIELHACK
+
+ } else {
+ cout << " SEcal05_Helpers ERROR _module_XYtype = " << _module_XYtype << "!!!" << endl;
+ assert(0);
+ }
+ layerslabdims.push_back(ss);
+ totTow++;
+ } // towers
+ } // modules
+
+ // if slab dimensions do not change layer-by-layer, memorise for next time
+ if ( _module_XZtype == 0 ) _constantSlabXYDimensions = layerslabdims;
+ }
+ return layerslabdims;
+}
+
+
+void SEcal05_Helpers::updateCaloLayers(double thickness,
+ dd4hep::Material mat,
+ bool isAbsorber,
+ bool isSensitive,
+ double cell_size_x, double cell_size_y,
+ bool isFinal
+ ) {
+
+ if ( isFinal ) { // add material before saving layer
+ _layer_thickness += (thickness);
+ _layer_nRadiationLengths += (thickness)/mat.radLength() ;
+ _layer_nInteractionLengths += (thickness)/mat.intLength() ;
+ }
+
+ // should we finish off the current layer?
+ if ( isFinal || // last slice of calo
+ (isAbsorber && _caloLayer.sensitive_thickness > 0 ) // end of a caloLayer
+ ) {
+
+ // finalise caloLayer entry, add to caloData
+ _caloLayer.outer_thickness = _layer_thickness;
+ _caloLayer.outer_nRadiationLengths = _layer_nRadiationLengths;
+ _caloLayer.outer_nInteractionLengths = _layer_nInteractionLengths;
+ //_caloLayer.thickness = _caloLayer.inner_thickness + _caloLayer.outer_thickness ;
+
+ // push it back
+ _caloData->layers.push_back( _caloLayer ) ;
+ // printSEcal05LayerInfo( _caloLayer );
+
+ // reset layer thicknesses
+ _layer_thickness=0;
+ _layer_nRadiationLengths=0;
+ _layer_nInteractionLengths=0;
+
+ // update the calolayer.distance ( DJeans 12 sep 2017)
+ // here this is distance from ECAL start to the layer start; the Barrel and Endcap drivers add the distance from IP
+ _caloLayer.distance += _caloLayer.inner_thickness+_caloLayer.outer_thickness;
+
+
+ }
+
+ if (!isFinal) {
+
+ // first add half the material
+ _layer_thickness += (thickness/2.);
+ _layer_nRadiationLengths += (thickness/2.)/mat.radLength() ;
+ _layer_nInteractionLengths += (thickness/2.)/mat.intLength() ;
+
+ // then we store material before and after centre of this layer
+ if ( isSensitive ) {
+ _caloLayer.cellSize0 = cell_size_x;
+ _caloLayer.cellSize1 = cell_size_y;
+ _caloLayer.sensitive_thickness = thickness ;
+ _caloLayer.inner_nRadiationLengths = _layer_nRadiationLengths ;
+ _caloLayer.inner_nInteractionLengths = _layer_nInteractionLengths ;
+ _caloLayer.inner_thickness = _layer_thickness ;
+
+ // reset layer thicknesses
+ _layer_thickness=0;
+ _layer_nRadiationLengths=0;
+ _layer_nInteractionLengths=0;
+ }
+
+ // add in remaining half of layer
+ _layer_thickness += (thickness/2.);
+ _layer_nRadiationLengths += (thickness/2.)/mat.radLength() ;
+ _layer_nInteractionLengths += (thickness/2.)/mat.intLength() ;
+ }
+
+ _totThick+=thickness;
+
+ return;
+}
+
+
+SEcal05_Helpers::dxinfo SEcal05_Helpers::getNormalMagicUnitsInX( double dx_total, double dx_unit, double dx_cell, double dx_dead ,
+ int magicStrategy ) {
+
+ dxinfo dxInf;
+ dxInf.normal_nX=-1;
+ dxInf.magic1_unitDX=-1;
+ dxInf.magic1_ncellsX=-1;
+ dxInf.magic2_unitDX=-1;
+
+ int nNormalUnit = int( floor( dx_total / dx_unit ) );
+
+ if ( magicStrategy==0 ) { // just integer number of standard units (wafers/megatiles)
+ dxInf.normal_nX = nNormalUnit;
+ } else {
+
+ double extraSpace = dx_total - nNormalUnit*dx_unit;
+ double extraSensitiveSpace = extraSpace - 2.*dx_dead;
+ double extraNCells = extraSensitiveSpace / dx_cell;
+
+ if ( magicStrategy==1 ) { // last magic megatile has integer number of standard-sized cells
+ dxInf.normal_nX = nNormalUnit;
+ if ( extraNCells > 1.0 ) {
+ int iext = int( floor( extraNCells ) );
+ dxInf.magic1_unitDX = iext*dx_cell + 2.*dx_dead;
+ dxInf.magic1_ncellsX = iext;
+ }
+ } else if ( magicStrategy==2 ) {
+
+ // one magic megatile with integer number of standard cells,
+ // second with non-standard cell size to fill space exactly
+ // last magic cell size in range shortestMacigCell<magicCellsSize/standardCellSize<shortestMagicCell
+
+ // the shortest magic cell (in terms of the usual cell size)
+ const double shortestMagicCell = 1.0; // this means that the last magic cell is at least as long as the standard cell
+
+ extraSensitiveSpace = extraSpace - 4.*dx_dead; // because both the magic tiles may have a dead area
+ extraNCells = extraSensitiveSpace / dx_cell;
+
+ if ( extraNCells < shortestMagicCell ) { // too small to make magic unit. remove one normal unit
+ nNormalUnit-=1;
+ extraSpace = dx_total - nNormalUnit*dx_unit;
+ extraSensitiveSpace = extraSpace - 4.*dx_dead;
+ extraNCells = extraSensitiveSpace / dx_cell;
+ }
+
+ int imagic1 = int( floor( extraNCells ) ); // rounded down number of standard-sized cells
+
+ double magic1size = imagic1>0 ? imagic1*dx_cell + 2*dx_dead : 0.;
+
+ double magic2cellsize = extraSpace - magic1size - 2*dx_dead;
+
+ if ( magic2cellsize/dx_cell < shortestMagicCell ) { // too small
+ imagic1 -= 1; // remove last standard cell
+ magic1size = imagic1>0 ? imagic1*dx_cell + 2*dx_dead : 0;
+ magic2cellsize = extraSpace - magic1size - 2*dx_dead;
+ }
+
+ assert( magic2cellsize/dx_cell >= shortestMagicCell && magic2cellsize/dx_cell <= shortestMagicCell+1.0 && "problem in deciding magic unit size" );
+
+ dxInf.normal_nX = nNormalUnit;
+ dxInf.magic1_ncellsX = imagic1;
+ dxInf.magic1_unitDX = imagic1>0 ? imagic1*dx_cell + 2*dx_dead : 0;
+
+ dxInf.magic2_unitDX = magic2cellsize + 2*dx_dead;
+
+ }
+ }
+
+ if ( magicStrategy==2 ) {
+ // check it fills exactly
+ if ( fabs( dxInf.normal_nX*dxInf.magic1_unitDX + dxInf.magic1_unitDX + dxInf.magic2_unitDX - dx_total ) < 0.01*dd4hep::mm ) {
+ cout << " SEcal05_Helpers ERROR : " << endl;
+ cout << dxInf.normal_nX*dxInf.magic1_unitDX << " " << dxInf.magic1_unitDX << " " << dxInf.magic2_unitDX << endl;
+ cout << dxInf.normal_nX*dxInf.magic1_unitDX + dxInf.magic1_unitDX + dxInf.magic2_unitDX << " " << dx_total << endl;
+ cout << dxInf.normal_nX*dxInf.magic1_unitDX + dxInf.magic1_unitDX + dxInf.magic2_unitDX - dx_total << endl;
+ assert(0 && "magic unit does not fill space exactly");
+ }
+
+ }
+
+ return dxInf;
+}
+
+
+void SEcal05_Helpers::makeModule( dd4hep::Volume & mod_vol, // the volume we'll fill
+ dd4hep::DetElement & stave_det, // the detector element
+ dd4hep::rec::LayeredCalorimeterData & caloData, // the reco data we'll fill
+ dd4hep::Detector & theDetector,
+ dd4hep::SensitiveDetector & sens
+ ) {
+ // make the module
+ _caloData = &caloData;
+
+ // calculate widths of module/tower/unit (in z-direction for barrel: across the slab/module
+ assert ( _ntowers.size()>0 && _unitsPerTower>0 && "_ntowers or _unitsPerTower not set" );
+
+ _module_dY_total=0;
+ for (size_t i=0; i<_ntowers.size(); i++) {
+ _module_dY_total = _ntowers[i]*_alveolus_total_dim_Y + 2.*_moduleGap;
+ }
+
+ double alveolus_active_dim_Y = _alveolus_total_dim_Y - 2.*_towerGap;
+ double unit_dim_Y = alveolus_active_dim_Y/_unitsPerTower;
+ double unit_sensitive_dim_Y = unit_dim_Y - 2*_unitDeadEdge; // width of the sensitive area in each sensor
+ assert( unit_sensitive_dim_Y>0 && "negative-sized sensitive area..." ); // otherwise really weird!
+
+ // get detector stuff
+ assert ( _x_det && "_x_det not set");
+ int det_id = _x_det->id();
+ xml_comp_t x_staves = _x_det->staves();
+
+ _carbon_fibre_material = theDetector.material("CarbonFiber");
+ _radiator_material = theDetector.material(x_staves.materialStr());
+ _air_material = theDetector.air();
+
+ dd4hep::VisAttr _radiator_visatt = theDetector.visAttributes( x_staves.visStr() );
+
+ // get segmentation stuff
+ assert (_geomseg && "segmentation not set");
+
+ dd4hep::DDSegmentation::WaferGridXY* waferSeg = dynamic_cast< dd4hep::DDSegmentation::WaferGridXY* > ( _geomseg->segmentation() ) ;
+
+
+ dd4hep::DDSegmentation::MegatileLayerGridXY* megatileSeg = dynamic_cast< dd4hep::DDSegmentation::MegatileLayerGridXY* > ( _geomseg->segmentation() ) ;
+ assert( (waferSeg || megatileSeg) && "no segmentation found" );
+
+ // set up the standard megatile size and offset
+ if ( megatileSeg ) {
+ megatileSeg->setMegaTileSizeXY( unit_sensitive_dim_Y, unit_sensitive_dim_Y );
+ megatileSeg->setMegaTileOffsetXY( -unit_sensitive_dim_Y/2., -unit_sensitive_dim_Y/2. );
+ }
+
+ _module_thickness = getTotalThickness();
+
+ double currentLayerBase_pos_Z = 0; // this gets updated: it's the position of the bottom face of the next detector slice
+
+ int layer_index = 0; // 1;// layer number - Daniel changes to c-type counting from 0...lets get rid of fortran habits!
+
+ // to deal with initial layers
+ bool isFrontFace = true;
+ int myLayerNum = 0 ; // this is the sensitive layer number (one per sensitive)
+ unsigned int absorber_index(0); // keep track of which absorber layer we're in
+
+ // keep track of thickness within a layer
+ // initialise
+ _layer_thickness = 0. ;
+ _layer_nRadiationLengths = 0. ;
+ _layer_nInteractionLengths = 0. ;
+
+ //--------------------------------
+ // loop over the layers
+ //---------------------------------
+
+ for(xml_coll_t li(*_x_det,_U(layer)); li; ++li) { // types of layers (i.e. thin/thick absorber) or "stack"
+ xml_comp_t x_layer = li;
+ // cout << " ---- NEW LAYER TYPE " << layer_index << " repeat " << x_layer.repeat() << endl;
+
+ // Loop over number of repeats for this layer type
+ for (int j=0; j< x_layer.repeat(); j++) { // layers within this type (or "stack")
+ std::string l_name = dd4hep::_toString(layer_index,"layer%d");
+
+ // #########################
+ // Build Structure Layer
+ // #########################
+
+ //----------------------------------------------------
+ // position and thickness of absorber plates in the structure
+ //----------------------------------------------------
+ double radiator_dim_Z(0);
+ double rad_pos_Z(0);
+ double this_struct_CFthick_beforeAbs(0);
+ double this_struct_CFthick_afterAbs(0);
+
+ if ( isFrontFace ) { // the first part of the module depends on whether we have preshwer or not
+ if ( _preshower==1 ) { // don't include W+CF wrapping; only one side of alveolus
+ this_struct_CFthick_beforeAbs = 0;
+ this_struct_CFthick_afterAbs = _CF_front + _CF_alvWall;
+ radiator_dim_Z = 0;
+ } else { // include W+CF wrapping; only one side of alveolus
+ // cout << " -- no preshower, including absorber" << endl;
+ this_struct_CFthick_beforeAbs = _CF_front + _CF_absWrap; // allow for initial CF front plate also in no preshower case - djeans 6 july 2017
+ this_struct_CFthick_afterAbs = _CF_absWrap + _CF_alvWall;
+ radiator_dim_Z = getAbsThickness( absorber_index++ );
+ rad_pos_Z = radiator_dim_Z/2. + this_struct_CFthick_beforeAbs; // distance from top surface of structure to centre of radiator
+ }
+ isFrontFace=false;
+ } else { // internal layer: include W+CF wrapping; both sides of alveolus
+ this_struct_CFthick_beforeAbs = _CF_alvWall+_CF_absWrap;
+ this_struct_CFthick_afterAbs = _CF_alvWall+_CF_absWrap;
+ radiator_dim_Z = getAbsThickness( absorber_index++ );
+ assert( radiator_dim_Z>0 && "no radiator!" );
+ rad_pos_Z = this_struct_CFthick_beforeAbs + radiator_dim_Z/2.; // distance from top surface of structure to centre of radiator
+ }
+
+
+ // create the radiator volume
+ if ( radiator_dim_Z>0 ) { // only create the volume if we have radiator
+ std::vector < dimposXYStruct > absorbersheets = getAbsPlateXYDimensions( currentLayerBase_pos_Z + this_struct_CFthick_beforeAbs + radiator_dim_Z ); // add CF before abs. added djeans 21 nov 2016
+ // create and place the absorber sheets
+ for ( size_t ipl=0; ipl<absorbersheets.size(); ipl++) {
+ dimposXYStruct plSize = absorbersheets[ipl];
+ dd4hep::Box barrelStructureLayer_box( plSize.sizeX/2.,
+ plSize.sizeY/2. - _CF_absWrap, // remove CF wrapping
+ radiator_dim_Z/2.);
+
+ dd4hep::Volume barrelStructureLayer_vol( _det_name+"_"+l_name+"_"+dd4hep::_toString(int(ipl),"bs%02d"),
+ barrelStructureLayer_box, _radiator_material);
+
+ barrelStructureLayer_vol.setVisAttributes( _radiator_visatt );
+
+ // Position the layer.
+ dd4hep::Position bsl_pos = getTranslatedPosition(plSize.posX, plSize.posY, currentLayerBase_pos_Z + rad_pos_Z );
+
+ // dd4hep::PlacedVolume barrelStructureLayer_phv =
+ mod_vol.placeVolume(barrelStructureLayer_vol, bsl_pos);
+
+ } // loop over sheets
+ }
+
+ // update layer thickness until front face of absorber
+ updateCaloLayers( this_struct_CFthick_beforeAbs, _carbon_fibre_material, false, false);
+ updateCaloLayers( radiator_dim_Z, _radiator_material, true, false );
+ updateCaloLayers( this_struct_CFthick_afterAbs, _carbon_fibre_material, false, false);
+
+ // update position within module
+ currentLayerBase_pos_Z += this_struct_CFthick_beforeAbs + radiator_dim_Z + this_struct_CFthick_afterAbs;
+
+ // #########################
+ // Build Slab
+ // #########################
+
+ //-------------------------------
+ // first sum the various materials for the caloData/caloLayer
+ //-------------------------------
+ int myLayerNumTemp = myLayerNum;
+ for(xml_coll_t si(x_layer,_U(slice)); si; ++si) {
+ xml_comp_t x_slice = si;
+ double s_thick = x_slice.thickness();
+ dd4hep::Material slice_material = theDetector.material( x_slice.materialStr() );
+ if (x_slice.materialStr().compare(x_staves.materialStr()) == 0){
+ // this is absorber material
+ // check it's consistent with what we expect from the detector parameters
+ assert ( fabs( s_thick - getAbsThickness( absorber_index++ ) ) < 1e-5 && "inconsistent radiator thickness" );
+ updateCaloLayers( s_thick, slice_material, true, false ); // absorber
+ } else if ( x_slice.isSensitive() ) {
+ double cell_size_x(0), cell_size_y(0);
+ if ( waferSeg ) {
+ cell_size_x = waferSeg->cellDimensions(0)[0];
+ cell_size_y = waferSeg->cellDimensions(0)[1];
+ } else if ( megatileSeg ) {
+ // setup megatile
+ int laytype = _layerConfig [ myLayerNumTemp%_layerConfig.size() ];
+ if ( laytype==0 ) {
+ megatileSeg->setMegaTileCellsXY( myLayerNumTemp, _cells_across_megatile , _cells_across_megatile );
+ } else if ( laytype==1 ) {
+ megatileSeg->setMegaTileCellsXY( myLayerNumTemp, _strips_across_megatile , _strips_along_megatile ); // strips in one orientation
+ } else if ( laytype==2 ) {
+ megatileSeg->setMegaTileCellsXY( myLayerNumTemp, _strips_along_megatile , _strips_across_megatile ); // and in the other
+ } else {
+ assert(0 && "unknown layer type");
+ }
+ cell_size_x = megatileSeg->cellDimensions(myLayerNumTemp, 0)[0]; // dummy wafer
+ cell_size_y = megatileSeg->cellDimensions(myLayerNumTemp, 0)[1];
+ }
+ updateCaloLayers( s_thick, slice_material, false, true, cell_size_x, cell_size_y ); // sensitive
+ myLayerNumTemp++;
+ } else {
+ updateCaloLayers( s_thick, slice_material, false, false );
+ }
+ }
+
+ //------------------------------------
+ // then actually construct the slabs
+ //------------------------------------
+ double slab_dim_Z = _layering->layer(layer_index)->thickness();
+ double slab_pos_Z = currentLayerBase_pos_Z + slab_dim_Z/2.; // centre position of slab
+
+ std::vector <dimposXYStruct> slabDims = getSlabXYDimensions( currentLayerBase_pos_Z + slab_dim_Z );
+
+ for (size_t islab = 0; islab<slabDims.size(); islab++) {
+ myLayerNumTemp = myLayerNum;
+ double slab_dim_X = slabDims[islab].sizeX;
+
+ // make an air volume for the alveolus
+ dd4hep::Box l_box( slab_dim_X/2. ,
+ slabDims[islab].sizeY/2. - _CF_alvWall,
+ slab_dim_Z/2. );
+
+ dd4hep::Volume l_vol( _det_name+"_alveolus_"+l_name, l_box, _air_material);
+ l_vol.setVisAttributes(theDetector.visAttributes( "GrayVis" ) );
+
+ dd4hep::DetElement l_det( stave_det, l_name+dd4hep::_toString(int(islab),"tower%02d") , det_id );
+ dd4hep::Position l_pos = getTranslatedPosition(slabDims[islab].posX, slabDims[islab].posY, slab_pos_Z );
+ dd4hep::PlacedVolume l_phv = mod_vol.placeVolume(l_vol,l_pos);
+ l_phv.addPhysVolID("tower", int(islab) );
+ l_det.setPlacement(l_phv);
+
+ // then fill it with the slab sublayers
+ int s_num(0);
+ double s_pos_Z = -slab_dim_Z / 2.; // position of sub-layer with respect to centre of this slab
+
+ for(xml_coll_t si(x_layer,_U(slice)); si; ++si) { // the sub-layers
+ xml_comp_t x_slice = si;
+ std::string s_name = dd4hep::_toString(s_num,"slice%d");
+ double s_thick = x_slice.thickness();
+
+ dd4hep::Material slice_material = theDetector.material( x_slice.materialStr() );
+
+ std::string vis_str = x_slice.visStr();
+
+ if ( !x_slice.isSensitive() ) { // not the sensitive slice: just a layer of stuff
+
+ dd4hep::Box s_box( slab_dim_X/2. , slabDims[islab].sizeY/2. - _CF_alvWall, s_thick/2. );
+ dd4hep::Volume s_vol(_det_name+"_"+l_name+"_"+s_name, s_box, slice_material);
+ s_vol.setVisAttributes(theDetector.visAttributes( vis_str ));
+
+ dd4hep::Position s_pos( 0, 0, s_pos_Z + s_thick/2. );
+ // dd4hep::PlacedVolume slice_phv =
+ l_vol.placeVolume(s_vol, s_pos );
+
+ } else { // sensitive slice
+
+ // Normal squared wafers - this is just the sensitive part
+ // square piece of silicon, not including guard ring. guard ring material is not included
+
+ dd4hep::Box WaferSiSolid( unit_sensitive_dim_Y/2., unit_sensitive_dim_Y/2., s_thick/2.);
+
+ // get the standard cell size in X for this layer
+ double cell_size_x = waferSeg ? waferSeg->cellDimensions(0)[0] : megatileSeg->cellDimensions(myLayerNumTemp, 0)[0];
+ double cell_size_y = waferSeg ? waferSeg->cellDimensions(0)[1] : megatileSeg->cellDimensions(myLayerNumTemp, 0)[1];
+
+ // work out how to make the magic units, if requested
+ dxinfo xseg = getNormalMagicUnitsInX( slab_dim_X,
+ unit_dim_Y,
+ cell_size_x,
+ _unitDeadEdge,
+ _magicMegatileStrategy );
+
+ int n_wafers_x = xseg.normal_nX;
+ int ncellsy = int( unit_sensitive_dim_Y/cell_size_y + 0.5 ); // should be exactly integer. take nearest integer to account for possible rounding errors.
+
+ int wafer_num = 0;
+
+ double wafer_pos_X = -slab_dim_X/2.; // keep track of wafer position in X
+
+ for (int n_wafer_x = 0; n_wafer_x < n_wafers_x+2 ; n_wafer_x++) { // loop along slab, including magic megatile/wafer
+
+ double megatile_size_x = unit_dim_Y;
+ int ncellsx(-1);
+
+ bool isMagic = n_wafer_x >= n_wafers_x;
+
+ if ( n_wafer_x==n_wafers_x ) { // first magic unit
+ megatile_size_x = xseg.magic1_unitDX;
+ ncellsx = xseg.magic1_ncellsX;
+ } else if ( n_wafer_x==n_wafers_x+1 ) { // second magic unit
+ megatile_size_x = xseg.magic2_unitDX;
+ ncellsx = 1;
+ }
+
+ // cout << "ismagic " << isMagic << " wafer " << n_wafer_x << " / " << n_wafers_x << " : sizex " << megatile_size_x << " nx " << ncellsx << endl;
+
+ if ( megatile_size_x<=0 ) continue;
+
+ // cout << " PASSED" << endl;
+
+ double megatile_sensitive_size_x = megatile_size_x - 2*_unitDeadEdge;
+
+ for (int n_wafer_Y = 0; n_wafer_Y < _unitsPerTower; n_wafer_Y++) { // loop across slab (usually 2 wafers, or 1 EBU) [ie along beam dir for barrel module]
+
+ double wafer_pos_Y = -alveolus_active_dim_Y/2.0 + (n_wafer_Y+0.5)*unit_dim_Y;
+ wafer_num++;
+ std::string Wafer_name;
+ if ( isMagic ) Wafer_name="magic";
+ Wafer_name += dd4hep::_toString(wafer_num,"wafer%d");
+
+ dd4hep::Box* box = isMagic ? new dd4hep::Box( megatile_sensitive_size_x/2,unit_sensitive_dim_Y/2,s_thick/2.) : &WaferSiSolid;
+ dd4hep::Volume WaferSiLog(_det_name+"_"+l_name+"_"+s_name+"_"+Wafer_name,*box,slice_material);
+
+ std::string wafer_vis_str = isMagic ? "YellowVis" : vis_str;
+
+ WaferSiLog.setVisAttributes(theDetector.visAttributes( wafer_vis_str ));
+ WaferSiLog.setSensitiveDetector(sens);
+
+ dd4hep::Position w_pos(wafer_pos_X + megatile_size_x/2., wafer_pos_Y, s_pos_Z + s_thick/2. );
+ dd4hep::PlacedVolume wafer_phv = l_vol.placeVolume(WaferSiLog, w_pos );
+ wafer_phv.addPhysVolID("wafer", wafer_num);
+ wafer_phv.addPhysVolID("layer", myLayerNumTemp );
+
+ if ( isMagic ) {
+ if ( megatileSeg ) { // define the special megatile
+ megatileSeg->setSpecialMegaTile( myLayerNumTemp, wafer_num,
+ megatile_sensitive_size_x, unit_sensitive_dim_Y,
+ -megatile_size_x/2., -unit_sensitive_dim_Y/2., // the offset
+ ncellsx, ncellsy ); // the segmentation
+ }
+ } else { // not magic
+ if ( waferSeg ) { // Normal squared wafers, this waferOffsetX is 0.0 // if its an odd number of cells, need to do something?
+ waferSeg->setWaferOffsetX(myLayerNumTemp, wafer_num, 0.0);
+ }
+ } // isMagic
+
+ } // y-wafers
+
+ wafer_pos_X += megatile_size_x;
+
+ } // x-wafers
+ myLayerNumTemp++;
+ } // end sensitive slice
+
+ // Increment Z position of slice.
+ s_pos_Z += s_thick;
+
+ // Increment slice number.
+ ++s_num;
+
+ } // end slice within one slab
+ } // slabs
+ myLayerNum = myLayerNumTemp;
+ currentLayerBase_pos_Z += slab_dim_Z;
+ layer_index++; // this is the structural layer counter (one per slab, not per sensitive layer)
+ } // layers in compact file
+ } // layer types in compact file
+
+ // add material after last slab. Just CF
+ updateCaloLayers( _CF_alvWall + _CF_back, _carbon_fibre_material, false, false, -1, -1, true ); // the last layer
+
+ return;
+}
diff --git a/Detector/DetCEPCv4/src/calorimeter/SEcal05_Helpers.h b/Detector/DetCEPCv4/src/calorimeter/SEcal05_Helpers.h
new file mode 100644
index 00000000..80e8e829
--- /dev/null
+++ b/Detector/DetCEPCv4/src/calorimeter/SEcal05_Helpers.h
@@ -0,0 +1,309 @@
+#ifndef __SECAL05_HELPERS_H__
+#define __SECAL05_HELPERS_H__ 1
+
+#include "DD4hep/DetFactoryHelper.h"
+
+#include "DDRec/DetectorData.h"
+
+#include "XML/Layering.h"
+#include "XML/Utilities.h"
+
+#include "DD4hep/Segmentations.h"
+
+#include "DDSegmentation/MegatileLayerGridXY.h"
+
+#include "DDSegmentation/WaferGridXY.h"
+
+#include <iostream>
+
+#undef NDEBUG
+#include <assert.h>
+
+using std::cout;
+using std::endl;
+
+/*
+
+we build general ECAL modules on a horizontal table.
+
+X-Y are on the table
+Z is vertically upwards
+slabs are aligned along X
+the lower Z side will face the IP
+
+2 shapes when looking down on table:
+
+XYtype = 0: (barrel and endcap)
+ ----------------------------- ^ Y
+ | | |
+ | | |
+ ----------------------------- ----> X
+
+XYtype = 1: (endcap)
+ ----------------------
+ | \
+ | \
+ | \
+ | | ^
+ | | | dY_kink
+ ----------------------------- |
+
+2 shapes when looking from side:
+
+XZtype = 0: (endcap)
+ ---------------------------------
+ | |
+ | |
+ ---------------------------------
+
+XZtype = 1: (barrel)
+ -------------------------- ^ Z
+ / \ |
+ / \ |
+ / \ |
+ -------------------------------- ----> X
+
+the slabs are prefectly aligned at -ve X side
+"magic" wafers are at +ve X
+
+local position: 0,0 defined as -veX,Y,Z corner
+
+we start building from Z=0, the face nearest IP, moving in the +ve Z direction
+
+
+D.Jeans update: 03/2015
+dead space between slab end and module edge
+other than 8-fold symmetry for barrel
+
+ */
+
+class SEcal05_Helpers {
+
+ public:
+
+ SEcal05_Helpers();
+
+ ~SEcal05_Helpers() {
+ delete _layering;
+ }
+
+ // SET THE PARAMETERS
+
+ void setDet( xml_det_t* x_det ) {
+ _x_det = x_det;
+ _layering = new dd4hep::Layering(*x_det);
+ _det_name = x_det->nameStr();
+ }
+
+ void setAbsLayers( int nl1, double th1, int nl2, double th2, int nl3=0, double th3=0 );
+
+ void setLayerConfig( std::vector < int > layerConfig ) {
+ _layerConfig=layerConfig;
+ }
+
+ void setSegmentation( dd4hep::Segmentation & seg ) {
+ _geomseg = &seg;
+ }
+
+ void setSegmentation( dd4hep::Segmentation * seg ) {
+ _geomseg = seg;
+ }
+
+ void setNCells( int Ecal_cells_across_megatile, int Ecal_strips_across_megatile, int Ecal_strips_along_megatile) {
+ _cells_across_megatile = Ecal_cells_across_megatile ;
+ _strips_across_megatile = Ecal_strips_across_megatile ;
+ _strips_along_megatile = Ecal_strips_along_megatile ;
+ }
+
+ void setMagicMegatileStrategy ( int i ) {
+ // magic megatile strategy
+ // 0: no magic megatile at lend of slab
+ // 1: magic megatile includes integer number of standard-sized cells
+ // 2: last cell of magic megatile adjusted to fill up to end of slab (size between 1 and 2 times normal cell size)
+ assert ( i>=0 && i<=2 );
+ _magicMegatileStrategy = i;
+ }
+
+ void setPreshower( bool p ) {
+ _preshower = p? 1 : 0; // do we have preshower layer? (if 1, first layer is sensitive; if 0, first layer is absorber)
+ }
+
+ void setCFthickness( double absWrap, double alvWall, double front, double back ) {
+ _CF_absWrap = absWrap; // wrapping around absorber
+ _CF_alvWall = alvWall; // alveolar wall
+ _CF_front = front; // front (IP side) support plate
+ _CF_back = back; // back support plate
+ }
+
+ void setModuleDimensions( int XYtype, // module shape in XY
+ int XZtype, // module shape in XZ
+ double dX_max, // maximum extent in X
+ double dY_kink = -999, // distance from lowerY edge to kink for XYtype=2
+ double angle = M_PI/4. // angle, if not rectangular. pi/4 for octagon
+ ) {
+ if ( XYtype>0 ) assert ( XZtype==0 ); // endcap module has vertical edges
+ if ( XZtype>0 ) assert ( XYtype==0 ); // barrel module should have rectangular shadow
+
+ _module_XYtype = XYtype;
+ _module_XZtype = XZtype;
+ _module_dX_max = dX_max;
+ _module_dY_kink = dY_kink;
+ _module_angle = angle;
+
+ }
+
+ void setTowersUnits( std::vector <int> ntowers, // a vector of modules, containing the specified # of towers/alveolii
+ double towerWidth,
+ int unitsPerTower, // # megatiles/units per tower
+ double moduleDeadEdge, // dead region at edge of module
+ double towerDeadEdge, // dead region at edge of tower
+ double unitDeadEdge // width of dead region at edge of each unit (e.g. wafer's guard ring)
+ ) {
+ _ntowers = ntowers; // number of towers (or alveoli) across Y
+ _alveolus_total_dim_Y = towerWidth; // width of tower (including dead area)
+ _unitsPerTower = unitsPerTower; // how many units (megatiles/wafers) across one tower
+ _moduleGap = moduleDeadEdge; // distance from tower boundary to edge of unit
+ _towerGap = towerDeadEdge; // distance from tower boundary to edge of unit
+ _unitDeadEdge = unitDeadEdge; // width of dead area around edge of unit
+ }
+
+ void checkLayerConsistency();
+
+ float getTotalThickness();
+
+ void setTranslation( dd4hep::Position trans ) {_trans = trans;}
+
+ // ---- this is the main workhorse
+ void makeModule( dd4hep::Volume & mod_vol, // the volume we'll fill
+ dd4hep::DetElement & stave_det, // the detector element
+ dd4hep::rec::LayeredCalorimeterData & caloData, // the reco data we'll fill
+ dd4hep::Detector & theDetector,
+ dd4hep::SensitiveDetector & sens
+ );
+
+ void setPlugLength( float ll ) { _plugLength = ll; }
+
+
+ private:
+
+ void printSEcal05LayerInfo( dd4hep::rec::LayeredCalorimeterData::Layer & caloLayer);
+
+ double getAbsThickness( unsigned int iAbsLay );
+
+ int getNumberOfAbsorberLayers() {
+ // total number of absorber layers
+ return _nlayers1+_nlayers2+_nlayers2;
+ }
+
+ int getNumberOfStructuralAbsorberLayers() {
+ // number of absorber layers in the module (not in the slabs)
+ assert( _preshower==0 || _preshower==1 );
+ return _preshower ? getNumberOfAbsorberLayers()/2 : (getNumberOfAbsorberLayers()-1)/2 ;
+ }
+
+ struct dimposXYStruct {
+ double sizeX, sizeY;
+ double posX, posY;
+ };
+
+ struct dxinfo {
+ int normal_nX;
+
+ double magic1_unitDX;
+ int magic1_ncellsX;
+
+ double magic2_unitDX;
+ };
+
+ xml_det_t* _x_det;
+ dd4hep::Layering* _layering=NULL;
+ std::string _det_name;
+
+ std::vector <dimposXYStruct> getAbsPlateXYDimensions( double ztop=-999 );
+
+ std::vector <dimposXYStruct> getSlabXYDimensions( double ztop=-999 );
+
+
+ dd4hep::Position getTranslatedPosition(double x, double y, double z) {
+ return dd4hep::Position ( x, y, z ) + _trans;
+ }
+
+ dxinfo getNormalMagicUnitsInX( double dx_total, double dx_unit, double dx_cell, double dx_dead ,
+ int magicStrategy );
+
+ void updateCaloLayers(double thickness,
+ dd4hep::Material mat,
+ bool isAbsorber,
+ bool isSensitive,
+ double cell_size_x=0, double cell_size_y=0,
+ bool isFinal=false
+ );
+
+
+ dd4hep::Segmentation* _geomseg;
+
+ dd4hep::Material _air_material;
+ dd4hep::Material _carbon_fibre_material;
+ dd4hep::Material _radiator_material;
+
+ unsigned int _cells_across_megatile;
+ unsigned int _strips_across_megatile;
+ unsigned int _strips_along_megatile;
+
+ int _preshower;
+
+ unsigned int _nlayers1;
+ unsigned int _nlayers2;
+ unsigned int _nlayers3;
+
+ double _radiator_thickness1;
+ double _radiator_thickness2;
+ double _radiator_thickness3;
+
+ std::vector < int > _layerConfig; // square or X or Y strips
+
+ double _CF_absWrap;
+ double _CF_alvWall;
+ double _CF_front;
+ double _CF_back;
+
+ // overall module size
+ int _module_XYtype;
+ int _module_XZtype;
+ double _module_dX_max;
+ double _module_dY_total;
+ double _module_dY_kink;
+ double _module_angle;
+
+ // internal details
+ std::vector <int> _ntowers;
+ double _towerGap;
+ double _moduleGap;
+ int _unitsPerTower;
+ double _unitDeadEdge;
+
+ dd4hep::rec::LayeredCalorimeterData* _caloData;
+ dd4hep::rec::LayeredCalorimeterData::Layer _caloLayer ; // this is the output info which is passed to reconstruction
+
+ double _layer_thickness;
+ double _layer_nRadiationLengths;
+ double _layer_nInteractionLengths;
+
+ double _totThick;
+
+ double _alveolus_total_dim_Y;
+ double _module_thickness;
+
+ int _magicMegatileStrategy;
+
+ dd4hep::Position _trans;
+
+ std::vector <dimposXYStruct> _constantSlabXYDimensions;
+
+
+ float _plugLength;
+
+
+};
+
+#endif
--
GitLab