From 8a7c23281446e8c9f2d99fa8ed5dbd0219e94ed2 Mon Sep 17 00:00:00 2001 From: varnes <150390087+varnes@users.noreply.github.com> Date: Wed, 18 Sep 2024 00:21:25 -0700 Subject: [PATCH] Fix printout issues in turbine endcap geometry; also add v02 of the geometry (#379) * Fix bug in cell y position for -z endcap * Add v02 of turbine endcal ecal * Fix up some printouts * Fix printouts * Fix printouts * Fix printouts * Resolve conflicts * update for v2 geometry * Fix printouts * Fix printouts * Update to be compatible with v2 segmentation * Fix printout * comment include of DRC in IDEA o1 v3 * IDEA with DRC dedicated test added * fix IDEA with DRC test * apply other SD action to DRC * update IDEA README * still not working... * replace absolute path by cmake variables in ctest * fix test_IDEA_o1_v03 by removing DRC steering file * IDEA+DRC test is now working * Fix a few compiler warnings (#374) * [FCCeeMDI] Use absolute path to import CAD files * ILD_l5_v11: fix a barrel layer ID (#378) * 2nd SIT barrel layer ID corrected; error stemed from outcommenting previous 2nd layer * removed commented layer to avoid any potential issue in the future * copy TrackerBarrel_o1_v06 * copy CLD_o2_v06 to CLD_o2_v07 * update version * Use new TrackerBarrel * TrackerBarrel_o1_v06 remove using namespace std * move NeighbourSurfacesData population into a function * Move LayerLayout population out of the sensor loop * refactor LayerLayout population into function * TrackerBarrel_o1_v06 add assembly stave * add tracking volume * update README * addExtensionsToHCalAllegro * implemented suggestions from discussion * implement comments and improve code readibility * fix failing test and remove unused phi segmentation * fix endcap bitfields and improve readibility * keep ALLEGRO v03 xml file unchanged * adding version nr also for the barrel xml * update hcal xml files to v02 * add back ECalEndcaps_Turbine * remove unused hcal v01 xml files * update READMEs * add HCalTileBarrel_o1_v02_geo.cpp * CLD_o2_v07: start with new lumical to fix overlaps * CLD_o2_v07: fix overlaps in the LumiCal itself Add a lumical_max_z_prime parameter to get the actual end place of the lumical add 7 micron meter to the envelopes. Fix the location and rotation to use the max_z_prime values * CLD_o2_v07: remove old lumical xml * CLD_o2_v07: InnerTracker: fix envelope overlap between lumical and innertracker * Tests: move CLD test to latest model * CLD_o2_v07: account for lumical positioning with difference between primed and unprimed * Fields: correct type to field_type, changed variable name in DD4hep, type still exists as part of handle * Tests: move one test to something with fieldmaps * Do not link against podio and EDM4hep dictionaries (#386) * Add a Key4hepConfig file * Update Key4hepConfig.cmake with fixes from https://github.com/key4hep/key4hep-dev-utils/pull/7 * Use the Key4hepConfig flag to set the standard, compiler flags (#387) * CLD_o2_v07: LumiCalBackShield: change envelope to an assembly instead of algobal union of intersections... No overlaps detected and the finding points on surface does not complain any more * CLD_o2_v07: LumiCal: replace all envelopes with assemblies at least the lumical specific overlap check with /geometry/test/resolution 300000 is much faster now as well??? * ILD models: apply Tracker_limits to the volume inside of beampipe and MDI (#388) * apply Tracker_limits to the volume inside of beampipe and related MDI elements --------- Co-authored-by: JEANS Daniel Thomelin Dietrich * Remove unused variable, and add information about o1_v02 to the README file * Initialize all variables * Initialize all variables --------- Co-authored-by: Erich Varnes Co-authored-by: Brieuc Francois Co-authored-by: Alvaro Tolosa Delgado Co-authored-by: Juan Miguel Carceller <22276694+jmcarcell@users.noreply.github.com> Co-authored-by: BrieucF Co-authored-by: Victor Schwan <162138084+Victor-Schwan@users.noreply.github.com> Co-authored-by: Leonhard Reichenbach Co-authored-by: michaela mlynarikova Co-authored-by: Andre Sailer Co-authored-by: jmcarcell Co-authored-by: Daniel Jeans Co-authored-by: JEANS Daniel Thomelin Dietrich --- .../ALLEGRO_o1_v03/ECalEndcaps_Turbine.xml | 8 +- .../ECalEndcaps_Turbine_calibration.xml | 8 +- .../ECalEndcap_Turbine_o1_v01_geo.cpp | 66 +- .../ECalEndcap_Turbine_o1_v02_geo.cpp | 671 ++++++++++++++++++ detector/calorimeter/README.md | 5 +- .../FCCSWEndcapTurbine_k4geo.h | 4 +- .../src/FCCSWEndcapTurbine_k4geo.cpp | 73 +- 7 files changed, 787 insertions(+), 48 deletions(-) create mode 100644 detector/calorimeter/ECalEndcap_Turbine_o1_v02_geo.cpp diff --git a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/ECalEndcaps_Turbine.xml b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/ECalEndcaps_Turbine.xml index 547053962..ac5b7f2ef 100644 --- a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/ECalEndcaps_Turbine.xml +++ b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/ECalEndcaps_Turbine.xml @@ -43,6 +43,11 @@ + + + + + @@ -51,6 +56,7 @@ + @@ -108,7 +114,7 @@ - + diff --git a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/ECalEndcaps_Turbine_calibration.xml b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/ECalEndcaps_Turbine_calibration.xml index 3d07d3671..1df7dd3af 100644 --- a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/ECalEndcaps_Turbine_calibration.xml +++ b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/ECalEndcaps_Turbine_calibration.xml @@ -35,6 +35,11 @@ + + + + + @@ -43,6 +48,7 @@ + @@ -100,7 +106,7 @@ - + diff --git a/detector/calorimeter/ECalEndcap_Turbine_o1_v01_geo.cpp b/detector/calorimeter/ECalEndcap_Turbine_o1_v01_geo.cpp index 51d79a656..4a4c43822 100644 --- a/detector/calorimeter/ECalEndcap_Turbine_o1_v01_geo.cpp +++ b/detector/calorimeter/ECalEndcap_Turbine_o1_v01_geo.cpp @@ -1,4 +1,5 @@ #include "DD4hep/DetFactoryHelper.h" +#include "DD4hep/Printout.h" #include "TMatrixT.h" // todo: remove gaudi logging and properly capture output @@ -83,18 +84,18 @@ namespace det { float BladeAngle = genericBladeElem.attr(_Unicode(angle)); bool decreaseAnglePerWheel = genericBladeElem.attr(_Unicode(decreaseAnglePerWheel)); - lLog << MSG::DEBUG << "Making wheel with inner, outer radii " << ri << ", " << ro << std:: endl; - lLog << MSG::DEBUG << "Blade angle is " << BladeAngle << "; decrease angle per wheel? " << decreaseAnglePerWheel << endmsg; + dd4hep::printout(dd4hep::DEBUG, "ECalEndcap_Turbine_o1_v01", "Making wheel with inner, outer radii %f, %f", ri, ro); + dd4hep::printout(dd4hep::DEBUG, "ECalEndcap_Turbine_o1_v01", "Blade angle is %f; decrease angle per wheel? ", BladeAngle, decreaseAnglePerWheel); dd4hep::xml::Dimension dim(aXmlElement.child(_Unicode(dimensions))); double grmin = dim.rmin1(); - lLog << MSG::DEBUG << "delZ is " << delZ << endmsg; + dd4hep::printout( dd4hep::DEBUG, "ECalEndcap_Turbine_o1_v01", "delZ is %f", delZ); if (decreaseAnglePerWheel) { float tubeFracCovered = delZ/(2*grmin*TMath::Tan(BladeAngle)); BladeAngle = TMath::ATan(delZ/(2*ri*tubeFracCovered)); } if (TMath::Abs(TMath::Tan(BladeAngle)) < delZ/(2.*ri)) { - lLog << MSG::ERROR << "The requested blade angle is too small for the given delZ and ri values. Please adjust to at least " << TMath::ATan(delZ/(2.*ri))*180./TMath::Pi() << " degrees!" << endmsg; + dd4hep::printout(dd4hep::ERROR, "ECalEndcap_Turbine_o1_v01", "The requested blade angle is too small for the given delZ and ri values. Please adjust to at least %f degrees!", TMath::ATan(delZ/(2.*ri))*180./TMath::Pi() ); return; } @@ -106,7 +107,7 @@ namespace det { float CladdingThick = claddingElem.attr(_Unicode(thickness)); float AbsThickMin = absBladeElem.attr(_Unicode(thickness))-(GlueThick+CladdingThick); if (AbsThickMin < 0.) { - lLog << MSG::ERROR << "Error: requested absorber thickness is negative after accounting for glue and cladding thickness" << endmsg; + dd4hep::printout(dd4hep::ERROR, "ECalEndcap_Turbine_o1_v01", "Error: requested absorber thickness is negative after accounting for glue and cladding thickness"); } float ElectrodeThick = electrodeBladeElem.attr(_Unicode(thickness)); float LArgapi = nobleLiquidElem.attr(_Unicode(gap)); @@ -127,7 +128,7 @@ namespace det { bool scaleBladeThickness = absBladeElem.attr(_Unicode(scaleThickness)); float bladeThicknessScaleFactor = absBladeElem.attr(_Unicode(thicknessScaleFactor)); - lLog << MSG::DEBUG << "nUnitCells: " << nUnitCells << endmsg; + dd4hep::printout(dd4hep::DEBUG, "ECalEndcap_Turbine_o1_v01", "nUnitCells: %d", nUnitCells); float AbsThicki = AbsThickMin; // make volumes for the noble liquid, electrode, and absorber blades @@ -148,7 +149,7 @@ namespace det { double leftoverS = (circ - nUnitCells*delrPhiNoGap); double delrPhiGapOnly = leftoverS/(2*nUnitCells); LArgapi = delrPhiGapOnly*TMath::Sin(BladeAngle); - lLog << MSG::DEBUG << "LArGap at inner radius is " << LArgapi << endmsg; + dd4hep::printout(dd4hep::DEBUG, "ECalEndcap_Turbine_o1_v01", "LArGap at inner radius is %f", LArgapi); // now find gap at outer radius circ = 2*TMath::Pi()*ro; @@ -193,7 +194,7 @@ namespace det { } else if (allNonActiveNotSensitive) { numNonActiveLayers = 1; } else { - lLog << MSG::ERROR << "Some non-active layers are sensitive and others are not -- this is likely a misconfiguration"; + dd4hep::printout(dd4hep::ERROR, "ECalEndcap_Turbine_o1_v01", "Some non-active layers are sensitive and others are not -- this is likely a misconfiguration"); } float delrNonActive = (ro-ri)/numNonActiveLayers; @@ -201,14 +202,14 @@ namespace det { for (unsigned iLayer = 0; iLayer < numNonActiveLayers; iLayer++) { float roLayer = riLayer + delrNonActive; - lLog << MSG::INFO << "Making layer in inner, outer radii " << riLayer << " " << roLayer << endmsg; + dd4hep::printout(dd4hep::INFO, "ECalEndcap_Turbine_o1_v01", "Making layer with inner, outer radii %f, %f", riLayer, roLayer); if (scaleBladeThickness) { AbsThicko = AbsThicki + bladeThicknessScaleFactor*((roLayer/riLayer)-1.)*AbsThicki; } else { AbsThicko = AbsThicki; } - lLog << MSG::DEBUG << "Inner and outer absorber thicknesses " << AbsThicki << " " << AbsThicko << endmsg; + dd4hep::printout(dd4hep::DEBUG, "ECalEndcap_Turbine_o1_v01", "Inner and outer absorber thicknesses %f, %f ", AbsThicki, AbsThicko); dd4hep::Solid claddingLayer = buildOneBlade(AbsThicki+GlueThick+CladdingThick, AbsThicko+GlueThick+CladdingThick, xRange, roLayer, riLayer, BladeAngle, delZ ); dd4hep::Solid glueLayer = buildOneBlade(AbsThicki+GlueThick, AbsThicko+GlueThick, xRange, roLayer, riLayer, BladeAngle, delZ ); @@ -264,8 +265,8 @@ namespace det { leftoverS = (circ - nUnitCells*delrPhiNoGap); delrPhiGapOnly = leftoverS/(2*nUnitCells); LArgapo = delrPhiGapOnly*TMath::Sin(BladeAngle); - lLog << MSG::DEBUG << "Outer LAr gap is " << LArgapo << endmsg ; - lLog << MSG::INFO << "Inner and outer thicknesses of noble liquid volume " << ElectrodeThick+LArgapi*2 << " " << ElectrodeThick+LArgapo*2 << endmsg; + dd4hep::printout(dd4hep::DEBUG, "ECalEndcap_Turbine_o1_v01", "Outer LAr gap is %f", LArgapo) ; + dd4hep::printout(dd4hep::INFO, "ECalEndcap_Turbine_o1_v01", "Inner and outer thicknesses of noble liquid volume %f, %f", ElectrodeThick+LArgapi*2, ElectrodeThick+LArgapo*2); dd4hep::Solid electrodeBladeAndGapLayer = buildOneBlade(ElectrodeThick+LArgapi*2, ElectrodeThick+LArgapo*2, xRange, roLayer, riLayer, BladeAngle, delZ); @@ -289,12 +290,12 @@ namespace det { LArgapi = LArgapo; AbsThicki = AbsThicko; } - lLog << MSG::INFO << "ECal endcap materials: nobleLiquid: " << nobleLiquidElem.materialStr() << " absorber: " << absBladeElem.materialStr() << " electrode: " << electrodeBladeElem.materialStr() << endmsg; + dd4hep::printout(dd4hep::INFO, "ECalEndcap_Turbine_o1_v01", "ECal endcap materials: nobleLiquid: %s absorber %s electrode %s", nobleLiquidElem.materialStr().c_str(), absBladeElem.materialStr().c_str(), electrodeBladeElem.materialStr().c_str() ); int nUnitCellsToDraw = nUnitCells; // nUnitCellsToDraw = 2; - lLog << MSG::INFO << "Number of unit cells "<< nUnitCells << endmsg; + dd4hep::printout(dd4hep::INFO, "ECalEndcap_Turbine_o1_v01", "Number of unit cells %d", nUnitCells); // place all components of the absorber blade inside passive volume @@ -310,7 +311,7 @@ namespace det { dd4hep::PlacedVolume absBladeVol_pv = glueLayerVols[iLayer].placeVolume(absBladeLayerVol, posLayer); absBladeVol_pv.addPhysVolID("subtype", 0); // 0 = absorber, 1 = glue, 2 = cladding - lLog << MSG::DEBUG << "Blade layer, rho is " << iLayer << " " << absBladeVol_pv.position().Rho() << " " << roLayer/2. << endmsg; + dd4hep::printout( dd4hep::DEBUG, "ECalEndcap_Turbine_o1_v01_geo", "Blade layer, rho is %d, %f, %f", iLayer, absBladeVol_pv.position().Rho(), roLayer/2.); absBladeVol_pv.addPhysVolID("layer", iWheel*numNonActiveLayers+iLayer); riLayer = roLayer; @@ -387,10 +388,9 @@ namespace det { double zminLayer = getZmin(riLayer, BladeAngle, delZ); dd4hep::Position posLayer(0,0,(zminLayer-zminri+roLayer-ro)/2.); - std::cout << "for active, riLayer, ri, roLayer, ro = " << riLayer << " " << ri << " " << roLayer << " " << ro << std::endl; dd4hep::PlacedVolume LArVol_pv(activeVol.placeVolume(LArTotalLayerVol, posLayer)); - lLog << MSG::DEBUG << "LAr layer: " << iLayer << endmsg; + dd4hep::printout(dd4hep::DEBUG, "ECalEndcap_Turbine_o1_v01", "LAr layer: %d", iLayer ); LArVol_pv.addPhysVolID("layer", iWheel*ECalEndcapNumCalibLayers+iLayer); riLayer = roLayer; @@ -404,7 +404,7 @@ namespace det { float phi = (iUnitCell-nUnitCellsToDraw/2)*2*TMath::Pi()/nUnitCells; float delPhi = 2*TMath::Pi()/nUnitCells; - lLog << MSG::DEBUG << "Placing blade, ro, ri = " << ro << " " << ri << endmsg; + dd4hep::printout( dd4hep::DEBUG, "ECalEndcap_Turbine_o1_v01", "Placing blade, ro, ri = %f %f", ro, ri); TGeoRotation tgr; tgr.RotateZ(BladeAngle*180/TMath::Pi()); tgr.RotateX(-phi*180/TMath::Pi()); @@ -466,7 +466,7 @@ namespace det { - lLog << MSG::DEBUG << "LArTotalLayerVols.size = " << LArTotalLayerVols.size() << endmsg; + dd4hep::printout(dd4hep::DEBUG, "ECalEndcap_Turbine_o1_v01", "LArTotalLayerVols.size = %d", LArTotalLayerVols.size()); } @@ -510,17 +510,17 @@ namespace det { dd4hep::Tube bathOuterShape(bathRmin, bathRmax, bathDelZ); // make it 4 volumes + 5th for detector envelope dd4hep::Tube bathAndServicesOuterShape(cryoDim.rmin2(), cryoDim.rmax1(), caloDim.dz()); // make it 4 volumes + 5th for detector envelope - lLog << MSG::INFO << "Cryostat front thickness is " << cryoDim.rmin2() << endmsg; + dd4hep::printout( dd4hep::INFO, "ECalEndcap_Turbine_o1_v01", "Cryostat front thickness is %f", cryoDim.rmin2() ); if (cryoThicknessFront > 0) { // 1. Create cryostat dd4hep::Tube cryoFrontShape(cryoDim.rmin1(), cryoDim.rmin2(), cryoDim.dz()); dd4hep::Tube cryoBackShape(cryoDim.rmax1(), cryoDim.rmax2(), cryoDim.dz()); dd4hep::Tube cryoSideOuterShape(cryoDim.rmin2(), cryoDim.rmax1(), cryoDim.dz()); dd4hep::SubtractionSolid cryoSideShape(cryoSideOuterShape, bathAndServicesOuterShape); - lLog << MSG::INFO << "ECAL endcap cryostat: front: rmin (cm) = " << cryoDim.rmin1() << " rmax (cm) = " << cryoDim.rmin2() << " dz (cm) = " << cryoDim.dz() << endmsg; - lLog << MSG::INFO << "ECAL encdap cryostat: back: rmin (cm) = " << cryoDim.rmax1() << " rmax (cm) = " << cryoDim.rmax2() << " dz (cm) = " << cryoDim.dz() << endmsg; - lLog << MSG::INFO << "ECAL endcap cryostat: side: rmin (cm) = " << cryoDim.rmin2() << " rmax (cm) = " << cryoDim.rmax1() << " dz (cm) = " << cryoDim.dz() - caloDim.dz() << endmsg; - lLog << MSG::INFO << "Cryostat is made out of " << cryostat.materialStr() << endmsg; + dd4hep::printout(dd4hep::INFO, "ECalEndcap_Turbine_o1_v01", "ECAL endcap cryostat: front: rmin (cm) = %f rmax (cm) = %f dz (cm) = %f ", cryoDim.rmin1(), cryoDim.rmin2(), cryoDim.dz()); + dd4hep::printout(dd4hep::INFO, "ECalEndcap_Turbine_o1_v01", "ECAL encdap cryostat: back: rmin (cm) = %f rmax (cm) = %f dz (cm) = %f", cryoDim.rmax1(), cryoDim.rmax2(), cryoDim.dz()); + dd4hep::printout( dd4hep::INFO, "ECalEndcap_Turbine_o1_v01", "ECAL endcap cryostat: side: rmin (cm) = %f rmax (cm) = %f dz (cm) = %f", cryoDim.rmin2(), cryoDim.rmax1(), cryoDim.dz() - caloDim.dz()); + dd4hep::printout( dd4hep::INFO, "ECalEndcap_Turbine_o1_v01", "Cryostat is made out of %s", cryostat.materialStr().c_str() ); dd4hep::Volume cryoFrontVol(cryostat.nameStr()+"_front", cryoFrontShape, aLcdd.material(cryostat.materialStr())); dd4hep::Volume cryoBackVol(cryostat.nameStr()+"_back", cryoBackShape, aLcdd.material(cryostat.materialStr())); @@ -533,19 +533,19 @@ namespace det { cryoFrontVol.setSensitiveDetector(aSensDet); cryoFrontPhysVol.addPhysVolID("cryo", 1); cryoFrontPhysVol.addPhysVolID("type", sidetype+1); - lLog << MSG::INFO << "Cryostat front volume set as sensitive" << endmsg; + dd4hep::printout( dd4hep::INFO, "ECalEndcap_Turbine_o1_v01", "Cryostat front volume set as sensitive"); } if (cryoBackSensitive) { cryoBackVol.setSensitiveDetector(aSensDet); cryoBackPhysVol.addPhysVolID("cryo", 1); cryoBackPhysVol.addPhysVolID("type", sidetype+2); - lLog << MSG::INFO << "Cryostat back volume set as sensitive" << endmsg; + dd4hep::printout( dd4hep::INFO, "ECalEndcap_Turbine_o1_v01", "Cryostat back volume set as sensitive" ); } if (cryoSideSensitive) { cryoSideVol.setSensitiveDetector(aSensDet); cryoSidePhysVol.addPhysVolID("cryo", 1); cryoSidePhysVol.addPhysVolID("type", sidetype+3); - lLog << MSG::INFO << "Cryostat front volume set as sensitive" << endmsg; + dd4hep::printout( dd4hep::INFO, "ECalEndcap_Turbine_o1_v01", "Cryostat front volume set as sensitive"); } dd4hep::DetElement cryoFrontDetElem(caloDetElem, "cryo_front", 0); cryoFrontDetElem.setPlacement(cryoFrontPhysVol); @@ -558,9 +558,7 @@ namespace det { // 2. Create noble liquid bath std::string nobleLiquidMaterial = nobleLiquid.materialStr(); dd4hep::Volume bathVol(nobleLiquidMaterial + "_bath", bathOuterShape, aLcdd.material(nobleLiquidMaterial)); - lLog << MSG::INFO << "ECAL endcap bath: material = " << nobleLiquidMaterial << " rmin (cm) = " << bathRmin - << " rmax (cm) = " << bathRmax << " dz (cm) = " << caloDim.dz() << " thickness in front of ECal (cm) = " << caloDim.rmin() - cryoDim.rmin2() - << " thickness behind ECal (cm) = " << cryoDim.rmax1() - caloDim.rmax() << endmsg; + dd4hep::printout( dd4hep::INFO, "ECalEndcap_Turbine_o1_v01", "ECAL endcap bath: material = %s rmin (cm) = %f rmax (cm) = %f, dz (cm) = %f, thickness in front of ECal (cm) = %f, thickness behind ECal (cm) = %f", nobleLiquidMaterial.c_str(), bathRmin, bathRmax, caloDim.dz(), caloDim.rmin() - cryoDim.rmin2(), cryoDim.rmax1() - caloDim.rmax()); dd4hep::DetElement bathDetElem(caloDetElem, "bath", 1); // 3. Create detector structure @@ -569,7 +567,7 @@ namespace det { dd4hep::xml::DetElement supportTubeElem = calo.child(_Unicode(supportTube)); unsigned nWheels = supportTubeElem.attr(_Unicode(nWheels)); - lLog << MSG::INFO << "Will build " << nWheels << " wheels" << endmsg; + dd4hep::printout(dd4hep::INFO, "ECalEndcap_Turbine_o1_v01", "Will build %d wheels", nWheels); double rmin = bathRmin; double rmax = bathRmax; float radiusRatio = pow(rmax/rmin, 1./nWheels); @@ -604,7 +602,7 @@ namespace det { dd4hep::PlacedVolume bathPhysVol = aEnvelope.placeVolume(bathVol); bathDetElem.setPlacement(bathPhysVol); - lLog << MSG::DEBUG << "Total number of modules: " << iModule << endmsg; + dd4hep::printout(dd4hep::DEBUG, "ECalEndcap_Turbine_o1_v01", "Total number of modules: %d", iModule); return; } @@ -634,11 +632,9 @@ createECalEndcapTurbine(dd4hep::Detector& aLcdd, dd4hep::xml::Handle_t aXmlEleme // dd4hep::DetElement caloPositiveDetElem(caloDetElem, "positive", 0); // dd4hep::DetElement caloNegativeDetElem(caloDetElem, "negative", 0); - lLog << MSG::DEBUG << "Placing dector on the positive side: (cm) " << dim.z_offset() << " with min, max radii " << dim.rmin1() << " " << dim.rmax1() << endmsg; + dd4hep::printout(dd4hep::DEBUG, "ECalEndcap_Turbine_o1_v01", "Placing detector on the positive side: (cm) %f with min, max radii %f %f",dim.z_offset(), dim.rmin1(), dim.rmax1() ); unsigned iModule = 0; buildOneSide_Turbine(aLcdd, aSensDet, envelopeVol, aXmlElement, iModule); - // lLog << MSG::DEBUG << "Placing dector on the negative side: (cm) " << -dim.z_offset() << " with min, max radii " << dim.rmin1() << " " << dim.rmax() << endmsg; - // buildOneSide_Turbine(aLcdd, aSensDet, envelopeNegativeVol, aXmlElement, -1, iModule); dd4hep::Assembly endcapsAssembly("ECalEndcaps_turbine"); diff --git a/detector/calorimeter/ECalEndcap_Turbine_o1_v02_geo.cpp b/detector/calorimeter/ECalEndcap_Turbine_o1_v02_geo.cpp new file mode 100644 index 000000000..a8fcafdc5 --- /dev/null +++ b/detector/calorimeter/ECalEndcap_Turbine_o1_v02_geo.cpp @@ -0,0 +1,671 @@ +#include "DD4hep/DetFactoryHelper.h" +#include "DD4hep/Printout.h" +#include "TMatrixT.h" + +// todo: remove gaudi logging and properly capture output +#define endmsg std::endl +#define lLog std::cout +namespace MSG { + const std::string ERROR = " Error: "; + const std::string DEBUG = " Debug: "; + const std::string INFO = " Info: "; +} + +namespace det { + + namespace ECalEndcap_Turbine_o1_v02 { + unsigned ECalEndCapElementCounter = 0; + + unsigned ECalEndcapNumCalibLayers; + + double tForArcLength(double s, double bladeangle, double delZ, double r) { + + // some intermediate constants + double zpos = delZ/2.; + double zp = zpos/TMath::Tan(bladeangle); + double b = zp/(TMath::Sqrt(r*r-zp*zp)); + double c = (TMath::Tan(s/r) +b)/(1.-b*TMath::Tan(s/r)); + double d = c*c*r*r/(1+c*c); + return (TMath::Sqrt(d)-zp)*TMath::Sin(bladeangle); + + // try approximating the arclength as dx. Less accurate, but that + // approximation is used in calculating the LAr gap, so maybe this + // will make it more consistent? + //return s*TMath::Sin(bladeangle); + + } + + // return position of the inner edge of a blade + double getZmin(double r, double bladeangle, double delZ) { + // r: distance from the beamline + // bladeangle: angle of turbine blades wrt xy plane, in radians + // delZ: z extent of the blades + return TMath::Sqrt(r*r - ((delZ/2)/TMath::Tan(bladeangle))*((delZ/2)/TMath::Tan(bladeangle))); + } + + dd4hep::Solid buildOneBlade(double thickness_inner, + double thickness_outer, + double width, + double ro, double ri, + double bladeangle, + double delZ) + { + + dd4hep::Solid shapeBeforeSubtraction; + + // set max and min extent of the blade (along the z axis in the body frame) + double zmax = ro; + double zmin = getZmin(ri, bladeangle, delZ); + + dd4hep::Trd2 tmp1(thickness_inner/2., thickness_outer/2., width/2., width/2., (zmax-zmin)/2. ); + shapeBeforeSubtraction = tmp1; + + dd4hep::Tube allowedTube(ri, ro, delZ); + + return dd4hep::IntersectionSolid (shapeBeforeSubtraction, allowedTube, dd4hep::Transform3D(dd4hep::RotationZYX( 0, TMath::Pi()/2.-bladeangle, TMath::Pi()/2.),dd4hep::Position(0,0, -(zmin+zmax)/2.))); + + } + + void buildWheel(dd4hep::Detector& aLcdd, + dd4hep::SensitiveDetector& aSensDet, + dd4hep::Volume& aEnvelope, + dd4hep::xml::Handle_t& aXmlElement, + dd4hep::DetElement& bathDetElem, + float ri, float ro, float delZ, + unsigned iWheel) { + + + dd4hep::xml::DetElement calorimeterElem = aXmlElement.child(_Unicode(calorimeter)); + dd4hep::xml::DetElement genericBladeElem = calorimeterElem.child(_Unicode(turbineBlade)); + dd4hep::xml::DetElement absBladeElem = genericBladeElem.child(_Unicode(absorberBlade)); + dd4hep::xml::DetElement claddingElem = genericBladeElem.child(_Unicode(cladding)); + dd4hep::xml::DetElement glueElem = genericBladeElem.child(_Unicode(glue)); + dd4hep::xml::DetElement electrodeBladeElem = genericBladeElem.child(_Unicode(electrodeBlade)); + dd4hep::xml::DetElement nobleLiquidElem = genericBladeElem.child(_Unicode(nobleLiquidGap)); + + float BladeAngle = 0.0, AbsThickMin = 0.0, BladeThicknessScaleFactor=0.0; + // hardcode for three wheels + if (iWheel == 0) { + BladeAngle = genericBladeElem.attr(_Unicode(angle1)); + AbsThickMin = absBladeElem.attr(_Unicode(thickness1)); + BladeThicknessScaleFactor = absBladeElem.attr(_Unicode(thicknessScaleFactor1)); + } + if (iWheel == 1) { + BladeAngle = genericBladeElem.attr(_Unicode(angle2)); + AbsThickMin = absBladeElem.attr(_Unicode(thickness2)); + BladeThicknessScaleFactor = absBladeElem.attr(_Unicode(thicknessScaleFactor2)); + } + if (iWheel == 2) { + BladeAngle = genericBladeElem.attr(_Unicode(angle3)); + AbsThickMin = absBladeElem.attr(_Unicode(thickness3)); + BladeThicknessScaleFactor = absBladeElem.attr(_Unicode(thicknessScaleFactor3)); + } + + dd4hep::printout(dd4hep::DEBUG, "ECalEndcap_Turbine_o1_v02", "Making wheel with inner, outer radii %f, %f", ri, ro); + dd4hep::printout(dd4hep::DEBUG, "ECalEndcap_Turbine_o1_v02", "Blade angle is %f ", BladeAngle); + dd4hep::xml::Dimension dim(aXmlElement.child(_Unicode(dimensions))); + dd4hep::printout( dd4hep::DEBUG, "ECalEndcap_Turbine_o1_v02", "delZ is %f", delZ); + + if (TMath::Abs(TMath::Tan(BladeAngle)) < delZ/(2.*ri)) { + dd4hep::printout(dd4hep::ERROR, "ECalEndcap_Turbine_o1_v02", "The requested blade angle is too small for the given delZ and ri values. Please adjust to at least %f degrees!", TMath::ATan(delZ/(2.*ri))*180./TMath::Pi() ); + return; + } + + Float_t xRange = delZ/(TMath::Sin(BladeAngle)); + + double delrPhiNoGap; + + float GlueThick = glueElem.attr(_Unicode(thickness)); + float CladdingThick = claddingElem.attr(_Unicode(thickness)); + + AbsThickMin = AbsThickMin-(GlueThick+CladdingThick); + if (AbsThickMin < 0.) { + dd4hep::printout(dd4hep::ERROR, "ECalEndcap_Turbine_o1_v02", "Error: requested absorber thickness is negative after accounting for glue and cladding thickness"); + } + float ElectrodeThick = electrodeBladeElem.attr(_Unicode(thickness)); + float LArgapi = nobleLiquidElem.attr(_Unicode(gap)); + + bool sameNUnitCells = genericBladeElem.attr(_Unicode(sameNUnitCells)); + char* nUnitCellsStrArr = (char*)genericBladeElem.attr(_Unicode(nUnitCells)).c_str(); + char* nUnitCellsCStr = strtok(nUnitCellsStrArr, " "); + int nUnitCells = -1; + if (!sameNUnitCells) { + for (unsigned i = 0; i < iWheel; i++) { + nUnitCellsCStr = strtok(NULL, " "); + } + std::string nUnitCellsStr = nUnitCellsCStr; + nUnitCells = std::stoi(nUnitCellsStr); + } + int nUnitCellsLeastCommonMultiple = genericBladeElem.attr(_Unicode(nUnitCellsLeastCommonMultiple)); + + dd4hep::printout(dd4hep::DEBUG, "ECalEndcap_Turbine_o1_v02", "nUnitCells: %d", nUnitCells); + + float AbsThicki = AbsThickMin; + // make volumes for the noble liquid, electrode, and absorber blades + float AbsThicko; + AbsThicko = AbsThicki + BladeThicknessScaleFactor*((ro/ri)-1.)*AbsThicki; + + // Calculate gap thickness at inner layer + double circ = 2*TMath::Pi()*ri; + double x2 =(AbsThickMin+(GlueThick+CladdingThick)+ElectrodeThick)/TMath::Sin(BladeAngle); + double y2 = TMath::Sqrt(ri*ri-x2*x2); + double rPhi1 = ri*TMath::Pi()/2.; + double rPhi2 = ri*TMath::ATan(y2/x2); + delrPhiNoGap = TMath::Abs(rPhi1-rPhi2); + double leftoverS = (circ - nUnitCells*delrPhiNoGap); + double delrPhiGapOnly = leftoverS/(2*nUnitCells); + LArgapi = delrPhiGapOnly*TMath::Sin(BladeAngle); + dd4hep::printout(dd4hep::DEBUG, "ECalEndcap_Turbine_o1_v02", "LArGap at inner radius is %f", LArgapi); + + // now find gap at outer radius + circ = 2*TMath::Pi()*ro; + x2 = (AbsThicko+GlueThick+CladdingThick+ElectrodeThick)/TMath::Sin(BladeAngle); + y2 = TMath::Sqrt(ro*ro-x2*x2); + rPhi1 = ro*TMath::Pi()/2.; + rPhi2 = ro*TMath::ATan(y2/x2); + delrPhiNoGap = TMath::Abs(rPhi1-rPhi2); + leftoverS = (circ - nUnitCells*delrPhiNoGap); + delrPhiGapOnly = leftoverS/(2*nUnitCells); + float LArgapo = delrPhiGapOnly*TMath::Sin(BladeAngle); + // LArgapo *= 2.; + + dd4hep::Solid absBlade; + float riLayer = ri; + + std::vector claddingLayerVols; + std::vector glueLayerVols; + std::vector absBladeLayerVols; + std::vector LArTotalLayerVols; + std::vector electrodeBladeLayerVols; + + + dd4hep::Solid passiveShape = buildOneBlade(AbsThicki+GlueThick+CladdingThick, AbsThicko+GlueThick+CladdingThick, xRange, ro, ri, BladeAngle, delZ ); + dd4hep::Volume passiveVol("passive", passiveShape, aLcdd.material("Air")); + + dd4hep::Solid activeShape = buildOneBlade(ElectrodeThick+LArgapi*2, ElectrodeThick+LArgapo*2, xRange, ro, ri, BladeAngle, delZ); + dd4hep::Volume activeVol("active", activeShape, aLcdd.material("Air")); + + unsigned numNonActiveLayers = 1; + // check that either all non-active volumes are set to sensitive (for + // sampling fraction calculations) or none are (for normal running) + bool allNonActiveSensitive = ( claddingElem.isSensitive() && + glueElem.isSensitive() && + absBladeElem.isSensitive() && + electrodeBladeElem.isSensitive() ); + bool allNonActiveNotSensitive = ( !claddingElem.isSensitive() && + !glueElem.isSensitive() && + !absBladeElem.isSensitive() && + !electrodeBladeElem.isSensitive() ); + if (allNonActiveSensitive) { + numNonActiveLayers = ECalEndcapNumCalibLayers; + } else if (allNonActiveNotSensitive) { + numNonActiveLayers = 1; + } else { + dd4hep::printout(dd4hep::ERROR, "ECalEndcap_Turbine_o1_v02", "Some non-active layers are sensitive and others are not -- this is likely a misconfiguration"); + } + + float delrNonActive = (ro-ri)/numNonActiveLayers; + float delrActive = (ro-ri)/ECalEndcapNumCalibLayers; + + for (unsigned iLayer = 0; iLayer < numNonActiveLayers; iLayer++) { + float roLayer = riLayer + delrNonActive; + dd4hep::printout(dd4hep::INFO, "ECalEndcap_Turbine_o1_v02", "Making layer with inner, outer radii %f, %f", riLayer, roLayer); + + AbsThicko = AbsThicki + BladeThicknessScaleFactor*((roLayer/riLayer)-1.)*AbsThicki; + + dd4hep::printout(dd4hep::DEBUG, "ECalEndcap_Turbine_o1_v02", "Inner and outer absorber thicknesses %f, %f ", AbsThicki, AbsThicko); + dd4hep::Solid claddingLayer = buildOneBlade(AbsThicki+GlueThick+CladdingThick, AbsThicko+GlueThick+CladdingThick, xRange, roLayer, riLayer, BladeAngle, delZ ); + + dd4hep::Solid glueLayer = buildOneBlade(AbsThicki+GlueThick, AbsThicko+GlueThick, xRange, roLayer, riLayer, BladeAngle, delZ ); + + // dd4hep::SubtractionSolid claddingLayer(absGlueCladdingLayer, absGlueLayer); + dd4hep::Solid absBladeLayer = buildOneBlade(AbsThicki, AbsThicko, xRange, roLayer, riLayer, BladeAngle, delZ ); + + // dd4hep::SubtractionSolid glueLayer(absGlueLayer, absBladeLayer); + dd4hep::Volume claddingLayerVol("claddingLayer", claddingLayer, aLcdd.material(claddingElem.materialStr())); + if (claddingElem.isSensitive()) { + claddingLayerVol.setSensitiveDetector(aSensDet); + } + claddingLayerVols.push_back(claddingLayerVol); + + dd4hep::Volume glueLayerVol("glueLayer", glueLayer, aLcdd.material(glueElem.materialStr())); + if (glueElem.isSensitive()) { + glueLayerVol.setSensitiveDetector(aSensDet); + } + glueLayerVols.push_back(glueLayerVol); + + dd4hep::Volume absBladeLayerVol("absBladeLayer", absBladeLayer, aLcdd.material(absBladeElem.materialStr())); + if (absBladeElem.isSensitive()) { + absBladeLayerVol.setSensitiveDetector(aSensDet); + } + absBladeLayerVols.push_back(absBladeLayerVol); + + riLayer = roLayer; + AbsThicki = AbsThicko; + } + + riLayer = ri; + + AbsThicki = AbsThickMin; + + for (unsigned iLayer = 0; iLayer < ECalEndcapNumCalibLayers; iLayer++) { + + float roLayer = riLayer + delrActive; + + AbsThicko = AbsThicki + BladeThicknessScaleFactor*((roLayer/riLayer)-1.)*AbsThicki; + + // now find gap at outer layer + circ = 2*TMath::Pi()*roLayer; + x2 = (AbsThicko+GlueThick+CladdingThick+ElectrodeThick)/TMath::Sin(BladeAngle); + y2 = TMath::Sqrt(roLayer*roLayer-x2*x2); + rPhi1 = roLayer*TMath::Pi()/2.; + rPhi2 = roLayer*TMath::ATan(y2/x2); + delrPhiNoGap = TMath::Abs(rPhi1-rPhi2); + leftoverS = (circ - nUnitCells*delrPhiNoGap); + delrPhiGapOnly = leftoverS/(2*nUnitCells); + LArgapo = delrPhiGapOnly*TMath::Sin(BladeAngle); + dd4hep::printout(dd4hep::DEBUG, "ECalEndcap_Turbine_o1_v02", "Outer LAr gap is %f", LArgapo) ; + dd4hep::printout(dd4hep::INFO, "ECalEndcap_Turbine_o1_v02", "Inner and outer thicknesses of noble liquid volume %f, %f", ElectrodeThick+LArgapi*2, ElectrodeThick+LArgapo*2); + + dd4hep::Solid electrodeBladeAndGapLayer = buildOneBlade(ElectrodeThick+LArgapi*2, ElectrodeThick+LArgapo*2, xRange, roLayer, riLayer, BladeAngle, delZ); + + dd4hep::Solid electrodeBladeLayer = buildOneBlade(ElectrodeThick, ElectrodeThick, xRange, roLayer, riLayer, BladeAngle, delZ); + + dd4hep::Volume electrodeBladeLayerVol("electrodeBladeLayer", electrodeBladeLayer, aLcdd.material(electrodeBladeElem.materialStr())); + if (electrodeBladeElem.isSensitive()) { + electrodeBladeLayerVol.setSensitiveDetector(aSensDet); + } + electrodeBladeLayerVols.push_back(electrodeBladeLayerVol); + + // dd4hep::SubtractionSolid LArShapeTotalLayer(electrodeBladeAndGapLayer, electrodeBladeLayer); + dd4hep::Volume LArTotalLayerVol("LArTotalLayerVol", electrodeBladeAndGapLayer, aLcdd.material(nobleLiquidElem.materialStr())); + + if ( nobleLiquidElem.isSensitive() ) { + LArTotalLayerVol.setSensitiveDetector(aSensDet); + } + LArTotalLayerVols.push_back(LArTotalLayerVol); + + riLayer = roLayer; + LArgapi = LArgapo; + AbsThicki = AbsThicko; + } + dd4hep::printout(dd4hep::INFO, "ECalEndcap_Turbine_o1_v02", "ECal endcap materials: nobleLiquid: %s absorber %s electrode %s", nobleLiquidElem.materialStr().c_str(), absBladeElem.materialStr().c_str(), electrodeBladeElem.materialStr().c_str() ); + + int nUnitCellsToDraw = nUnitCells; + // nUnitCellsToDraw = 2; + + dd4hep::printout(dd4hep::INFO, "ECalEndcap_Turbine_o1_v02", "Number of unit cells %d", nUnitCells); + + // place all components of the absorber blade inside passive volume + + unsigned iLayer = 0; + + riLayer = ri; + + for (auto absBladeLayerVol: absBladeLayerVols) { + + float roLayer = riLayer+delrNonActive; + + dd4hep::Position posLayer(0,0,(riLayer-ri+roLayer-ro)/2.); + dd4hep::PlacedVolume absBladeVol_pv = glueLayerVols[iLayer].placeVolume(absBladeLayerVol, posLayer); + + absBladeVol_pv.addPhysVolID("subtype", 0); // 0 = absorber, 1 = glue, 2 = cladding + dd4hep::printout( dd4hep::DEBUG, "ECalEndcap_Turbine_o1_v02_geo", "Blade layer, rho is %d, %f, %f", iLayer, absBladeVol_pv.position().Rho(), roLayer/2.); + absBladeVol_pv.addPhysVolID("layer", iWheel*numNonActiveLayers+iLayer); + + riLayer = roLayer; + iLayer++; + } + + riLayer = ri; + iLayer =0; + + for (auto glueLayerVol: glueLayerVols) { + + float roLayer = riLayer+delrNonActive; + + dd4hep::Position posLayer(0,0,(riLayer-ri+roLayer-ro)/2.); + dd4hep::PlacedVolume glueVol_pv = claddingLayerVols[iLayer].placeVolume(glueLayerVol, posLayer); + + + glueVol_pv.addPhysVolID("subtype", 1); // 0 = absorber, 1 = glue, 2 = cladding + glueVol_pv.addPhysVolID("layer", iWheel*numNonActiveLayers+iLayer); + + // dd4hep::DetElement glueDetElem(passiveDetElem, "glue_",ECalEndCapElementCounter++); + // glueDetElem.setPlacement(glueVol_pv); + + riLayer = roLayer; + iLayer++; + } + + riLayer = ri; + iLayer =0; + + double zminri = getZmin(ri, BladeAngle, delZ); + + for (auto claddingLayerVol: claddingLayerVols) { + + float roLayer = riLayer+delrNonActive; + + double zminLayer = getZmin(riLayer, BladeAngle, delZ); + + dd4hep::Position posLayer(0,0,(zminLayer-zminri+roLayer-ro)/2.); + dd4hep::PlacedVolume claddingVol_pv = passiveVol.placeVolume(claddingLayerVol, posLayer); + + claddingVol_pv.addPhysVolID("subtype", 2); // 0 = absorber, 1 = glue, 2 = cladding + claddingVol_pv.addPhysVolID("layer", iWheel*numNonActiveLayers+iLayer); + + // dd4hep::DetElement claddingDetElem(passiveDetElem, "cladding_", ECalEndCapElementCounter++); + // claddingDetElem.setPlacement(claddingVol_pv); + + riLayer = roLayer; + iLayer++; + } + + + riLayer = ri; + iLayer = 0; + + for (auto electrodeBladeLayerVol: electrodeBladeLayerVols) { + + float roLayer = riLayer+delrActive; + + dd4hep::PlacedVolume electrodeBladeVol_pv = LArTotalLayerVols[iLayer].placeVolume(electrodeBladeLayerVol); + electrodeBladeVol_pv.addPhysVolID("layer", iWheel*numNonActiveLayers+iLayer); + + riLayer = roLayer; + iLayer++; + } + + riLayer = ri; + iLayer = 0; + + for (auto LArTotalLayerVol: LArTotalLayerVols) { + + float roLayer = riLayer+delrActive; + + double zminLayer = getZmin(riLayer, BladeAngle, delZ); + + dd4hep::Position posLayer(0,0,(zminLayer-zminri+roLayer-ro)/2.); + + dd4hep::PlacedVolume LArVol_pv(activeVol.placeVolume(LArTotalLayerVol, posLayer)); + dd4hep::printout(dd4hep::DEBUG, "ECalEndcap_Turbine_o1_v02", "LAr layer: %d", iLayer ); + LArVol_pv.addPhysVolID("layer", iWheel*ECalEndcapNumCalibLayers+iLayer); + + riLayer = roLayer; + iLayer++; + } + + for (int iUnitCell = 0; iUnitCell < nUnitCellsToDraw; iUnitCell++) { + + int modIndex = iUnitCell-nUnitCellsToDraw/2; + if (modIndex < 0) modIndex += nUnitCells; + float phi = (iUnitCell-nUnitCellsToDraw/2)*2*TMath::Pi()/nUnitCells; + float delPhi = 2*TMath::Pi()/nUnitCells; + + dd4hep::printout( dd4hep::DEBUG, "ECalEndcap_Turbine_o1_v02", "Placing blade, ro, ri = %f %f", ro, ri); + + TGeoRotation tgr; + tgr.RotateZ(BladeAngle*180/TMath::Pi()); + tgr.RotateX(-phi*180/TMath::Pi()); + tgr.RotateY(90); + + const Double_t *rotMatPtr; + + rotMatPtr = tgr.GetRotationMatrix(); + TMatrixT rotMat(3,3, rotMatPtr); + dd4hep::Rotation3D r3d; + r3d.SetComponents(rotMat(0,0), rotMat(0,1), rotMat(0,2), + rotMat(1,0), rotMat(1,1), rotMat(1,2), + rotMat(2,0), rotMat(2,1), rotMat(2,2)); + + tgr.Clear(); + tgr.RotateZ(BladeAngle*180/TMath::Pi()); + tgr.RotateX(-(phi+delPhi/2.)*180/TMath::Pi()); + tgr.RotateY(90); + + rotMatPtr = tgr.GetRotationMatrix(); + TMatrixT rotMat2(3,3, rotMatPtr); + dd4hep::Rotation3D r3d2; + r3d2.SetComponents(rotMat2(0,0), rotMat2(0,1), rotMat2(0,2), + rotMat2(1,0), rotMat2(1,1), rotMat2(1,2), + rotMat2(2,0), rotMat2(2,1), rotMat2(2,2)); + + riLayer = ri; + + float xCell = ((ro+zminri)/2.)*TMath::Cos(phi); + float yCell = ((ro+zminri)/2.)*TMath::Sin(phi); //ri*TMath::Sin(phi)/6.; + float zCell = 0.; + + dd4hep::Transform3D comCell(r3d, dd4hep::Translation3D(xCell,yCell,zCell)); + + // place passive volume in LAr bath + dd4hep::PlacedVolume passivePhysVol = aEnvelope.placeVolume(passiveVol, comCell); + passivePhysVol.addPhysVolID("module", modIndex*nUnitCellsLeastCommonMultiple/nUnitCells); + passivePhysVol.addPhysVolID("wheel", iWheel); + passivePhysVol.addPhysVolID("type", 1); // 0 = active, 1 = passive, 2 = readout + dd4hep::DetElement passiveDetElem(bathDetElem, "passive_" + std::to_string(iUnitCell)+"_"+std::to_string(iWheel), ECalEndCapElementCounter++); + passiveDetElem.setPlacement(passivePhysVol); + + // place active volume in LAr bath + + xCell = ((ro+zminri)/2.)*TMath::Cos(phi+delPhi/2.); + yCell = ((ro+zminri)/2.)*TMath::Sin(phi+delPhi/2.); //ri*TMath::Sin(phi)/6.; + zCell = 0.; + dd4hep::Transform3D comCell2(r3d2, dd4hep::Translation3D(xCell,yCell,zCell)); + dd4hep::PlacedVolume activePhysVol = aEnvelope.placeVolume(activeVol, comCell2); + activePhysVol.addPhysVolID("module", modIndex*nUnitCellsLeastCommonMultiple/nUnitCells); + activePhysVol.addPhysVolID("wheel", iWheel); + activePhysVol.addPhysVolID("type", 0); // 0 = active, 1 = passive, 2 = readout + + dd4hep::DetElement activeDetElem(bathDetElem, "active_" + std::to_string(iUnitCell)+"_"+std::to_string(iWheel), ECalEndCapElementCounter++); + activeDetElem.setPlacement(activePhysVol); + + riLayer = ri; + iLayer =0; + + + dd4hep::printout(dd4hep::DEBUG, "ECalEndcap_Turbine_o1_v02", "LArTotalLayerVols.size = %d", LArTotalLayerVols.size()); + + } + + return; + + } + + void buildOneSide_Turbine(dd4hep::Detector& aLcdd, dd4hep::SensitiveDetector& aSensDet, + dd4hep::Volume& aEnvelope, dd4hep::xml::Handle_t& aXmlElement, + unsigned& iModule) { + + dd4hep::xml::DetElement calo = aXmlElement.child(_Unicode(calorimeter)); + dd4hep::xml::Dimension caloDim(calo.dimensions()); + + + dd4hep::xml::DetElement blade = calo.child(_Unicode(turbineBlade)); + dd4hep::xml::DetElement nobleLiquid = blade.child(_Unicode(nobleLiquidGap)); + + dd4hep::xml::DetElement xmlDetElem = aXmlElement; + std::string nameDet = xmlDetElem.nameStr(); + dd4hep::DetElement caloDetElem(nameDet, xmlDetElem.id()); + + dd4hep::xml::Dimension dim(aXmlElement.child(_Unicode(dimensions))); + + //build cryostat + // Retrieve cryostat data + dd4hep::xml::DetElement cryostat = calo.child(_Unicode(cryostat)); + dd4hep::xml::Dimension cryoDim(cryostat.dimensions()); + double cryoThicknessFront = cryoDim.rmin2() - cryoDim.rmin1(); + + dd4hep::xml::DetElement cryoFront = cryostat.child(_Unicode(front)); + dd4hep::xml::DetElement cryoBack = cryostat.child(_Unicode(back)); + dd4hep::xml::DetElement cryoSide = cryostat.child(_Unicode(side)); + bool cryoFrontSensitive = cryoFront.isSensitive(); + bool cryoBackSensitive = cryoBack.isSensitive(); + bool cryoSideSensitive = cryoSide.isSensitive(); + + double bathRmin = caloDim.rmin(); // - margin for inclination + double bathRmax = caloDim.rmax(); // + margin for inclination + double bathDelZ = caloDim.dz(); + dd4hep::Tube bathOuterShape(bathRmin, bathRmax, bathDelZ); // make it 4 volumes + 5th for detector envelope + dd4hep::Tube bathAndServicesOuterShape(cryoDim.rmin2(), cryoDim.rmax1(), caloDim.dz()); // make it 4 volumes + 5th for detector envelope + + dd4hep::printout( dd4hep::INFO, "ECalEndcap_Turbine_o1_v02", "Cryostat front thickness is %f", cryoDim.rmin2() ); + + if (cryoThicknessFront > 0) { + // 1. Create cryostat + dd4hep::Tube cryoFrontShape(cryoDim.rmin1(), cryoDim.rmin2(), cryoDim.dz()); + dd4hep::Tube cryoBackShape(cryoDim.rmax1(), cryoDim.rmax2(), cryoDim.dz()); + dd4hep::Tube cryoSideOuterShape(cryoDim.rmin2(), cryoDim.rmax1(), cryoDim.dz()); + dd4hep::SubtractionSolid cryoSideShape(cryoSideOuterShape, bathAndServicesOuterShape); + dd4hep::printout(dd4hep::INFO, "ECalEndcap_Turbine_o1_v02", "ECAL endcap cryostat: front: rmin (cm) = %f rmax (cm) = %f dz (cm) = %f ", cryoDim.rmin1(), cryoDim.rmin2(), cryoDim.dz()); + dd4hep::printout(dd4hep::INFO, "ECalEndcap_Turbine_o1_v02", "ECAL encdap cryostat: back: rmin (cm) = %f rmax (cm) = %f dz (cm) = %f", cryoDim.rmax1(), cryoDim.rmax2(), cryoDim.dz()); + dd4hep::printout( dd4hep::INFO, "ECalEndcap_Turbine_o1_v02", "ECAL endcap cryostat: side: rmin (cm) = %f rmax (cm) = %f dz (cm) = %f", cryoDim.rmin2(), cryoDim.rmax1(), cryoDim.dz() - caloDim.dz()); + dd4hep::printout( dd4hep::INFO, "ECalEndcap_Turbine_o1_v02", "Cryostat is made out of %s", cryostat.materialStr().c_str() ); + + dd4hep::Volume cryoFrontVol(cryostat.nameStr()+"_front", cryoFrontShape, aLcdd.material(cryostat.materialStr())); + dd4hep::Volume cryoBackVol(cryostat.nameStr()+"_back", cryoBackShape, aLcdd.material(cryostat.materialStr())); + dd4hep::Volume cryoSideVol(cryostat.nameStr()+"_side", cryoSideShape, aLcdd.material(cryostat.materialStr())); + dd4hep::PlacedVolume cryoFrontPhysVol = aEnvelope.placeVolume(cryoFrontVol); + dd4hep::PlacedVolume cryoBackPhysVol = aEnvelope.placeVolume(cryoBackVol); + dd4hep::PlacedVolume cryoSidePhysVol = aEnvelope.placeVolume(cryoSideVol); + unsigned sidetype = 0x4; // probably not needed anymore... + if (cryoFrontSensitive) { + cryoFrontVol.setSensitiveDetector(aSensDet); + cryoFrontPhysVol.addPhysVolID("cryo", 1); + cryoFrontPhysVol.addPhysVolID("type", sidetype+1); + dd4hep::printout( dd4hep::INFO, "ECalEndcap_Turbine_o1_v02", "Cryostat front volume set as sensitive"); + } + if (cryoBackSensitive) { + cryoBackVol.setSensitiveDetector(aSensDet); + cryoBackPhysVol.addPhysVolID("cryo", 1); + cryoBackPhysVol.addPhysVolID("type", sidetype+2); + dd4hep::printout( dd4hep::INFO, "ECalEndcap_Turbine_o1_v02", "Cryostat back volume set as sensitive" ); + } + if (cryoSideSensitive) { + cryoSideVol.setSensitiveDetector(aSensDet); + cryoSidePhysVol.addPhysVolID("cryo", 1); + cryoSidePhysVol.addPhysVolID("type", sidetype+3); + dd4hep::printout( dd4hep::INFO, "ECalEndcap_Turbine_o1_v02", "Cryostat side volume set as sensitive" ); + } + dd4hep::DetElement cryoFrontDetElem(caloDetElem, "cryo_front", 0); + cryoFrontDetElem.setPlacement(cryoFrontPhysVol); + dd4hep::DetElement cryoBackDetElem(caloDetElem, "cryo_back", 0); + cryoBackDetElem.setPlacement(cryoBackPhysVol); + dd4hep::DetElement cryoSideDetElem(caloDetElem, "cryo_side", 0); + cryoSideDetElem.setPlacement(cryoSidePhysVol); + } + + // 2. Create noble liquid bath + std::string nobleLiquidMaterial = nobleLiquid.materialStr(); + dd4hep::Volume bathVol(nobleLiquidMaterial + "_bath", bathOuterShape, aLcdd.material(nobleLiquidMaterial)); + dd4hep::printout( dd4hep::INFO, "ECalEndcap_Turbine_o1_v02", "ECAL endcap bath: material = %s rmin (cm) = %f rmax (cm) = %f, dz (cm) = %f, thickness in front of ECal (cm) = %f, thickness behind ECal (cm) = %f", nobleLiquidMaterial.c_str(), bathRmin, bathRmax, caloDim.dz(), caloDim.rmin() - cryoDim.rmin2(), cryoDim.rmax1() - caloDim.rmax()); + dd4hep::DetElement bathDetElem(caloDetElem, "bath", 1); + + // 3. Create detector structure + double length = dim.dz() * 2.; + double zOffsetEnvelope = -length / 2.; + + dd4hep::xml::DetElement supportTubeElem = calo.child(_Unicode(supportTube)); + unsigned nWheels = supportTubeElem.attr(_Unicode(nWheels)); + dd4hep::printout(dd4hep::INFO, "ECalEndcap_Turbine_o1_v02", "Will build %d wheels", nWheels); + double rmin = bathRmin; + double rmax = bathRmax; + float radiusRatio = pow(rmax/rmin, 1./nWheels); + double ro = rmin*radiusRatio; + double ri = rmin; + + float supportTubeThickness=supportTubeElem.thickness(); + unsigned iSupportTube = 0; + + for (unsigned iWheel = 0; iWheel < nWheels; iWheel++) { + + dd4hep::Tube supportTube(ro, ro+supportTubeThickness, bathDelZ ); + + dd4hep::Volume supportTubeVol("supportTube", supportTube, aLcdd.material(supportTubeElem.materialStr())); + if (supportTubeElem.isSensitive()) { + supportTubeVol.setSensitiveDetector(aSensDet); + } + dd4hep::PlacedVolume supportTube_pv = bathVol.placeVolume(supportTubeVol, dd4hep::Position(0,0,zOffsetEnvelope + dim.dz() )); + supportTube_pv.addPhysVolID("cryo", 1); + // supportTube_pv.addPhysVolID("side",sign); + supportTube_pv.addPhysVolID("wheel", iWheel); + dd4hep::DetElement supportTubeDetElem(bathDetElem, "supportTube_"+std::to_string(iWheel), 0); + supportTubeDetElem.setPlacement(supportTube_pv); + + + buildWheel(aLcdd, aSensDet, bathVol, aXmlElement, bathDetElem, ri+supportTubeThickness, ro, bathDelZ*2, iWheel); + ri = ro; + ro *= radiusRatio; + if (ro > rmax) ro = rmax; + iSupportTube++; + } + + + dd4hep::PlacedVolume bathPhysVol = aEnvelope.placeVolume(bathVol); + bathDetElem.setPlacement(bathPhysVol); + + dd4hep::printout(dd4hep::DEBUG, "ECalEndcap_Turbine_o1_v02", "Total number of modules: %d", iModule); + + + return; + } + + + + static dd4hep::Ref_t + createECalEndcapTurbine(dd4hep::Detector& aLcdd, dd4hep::xml::Handle_t aXmlElement, dd4hep::SensitiveDetector aSensDet) { + + dd4hep::xml::DetElement xmlDetElem = aXmlElement; + std::string nameDet = xmlDetElem.nameStr(); + int idDet = xmlDetElem.id(); + dd4hep::xml::Dimension dim(xmlDetElem.dimensions()); + dd4hep::DetElement caloDetElem(nameDet, idDet); + dd4hep::xml::Dimension sdType = xmlDetElem.child(_U(sensitive)); + aSensDet.setType(sdType.typeStr()); + + ECalEndcapNumCalibLayers = aLcdd.constant("ECalEndcapNumCalibLayers"); + + + // Create air envelope for one endcap (will be copied to make both endcaps) + dd4hep::Tube endcapShape( dim.rmin1(), dim.rmax1(), dim.dz()); + + dd4hep::Volume envelopeVol(nameDet + "_vol", endcapShape, aLcdd.material("Air")); + + + // dd4hep::DetElement caloPositiveDetElem(caloDetElem, "positive", 0); + // dd4hep::DetElement caloNegativeDetElem(caloDetElem, "negative", 0); + + dd4hep::printout(dd4hep::DEBUG, "ECalEndcap_Turbine_o1_v02", "Placing detector on the positive side: (cm) %f with min, max radii %f %f",dim.z_offset(), dim.rmin1(), dim.rmax1() ); + + unsigned iModule = 0; + buildOneSide_Turbine(aLcdd, aSensDet, envelopeVol, aXmlElement, iModule); + + dd4hep::Assembly endcapsAssembly("ECalEndcaps_turbine"); + + // Place the envelope + dd4hep::Transform3D envelopePositiveVolume_tr(dd4hep::RotationZYX( 0 ,0,0), dd4hep::Translation3D(0, 0, dim.z_offset())); + dd4hep::PlacedVolume envelopePositivePhysVol = endcapsAssembly.placeVolume(envelopeVol, envelopePositiveVolume_tr); + envelopePositivePhysVol.addPhysVolID("side", 1); + + dd4hep::DetElement caloPositiveDetElem(caloDetElem, "positive", 0); + caloPositiveDetElem.setPlacement(envelopePositivePhysVol); + + // make another placement for the negative z endcap + dd4hep::Transform3D envelopeNegativeVolume_tr(dd4hep::RotationZYX( 0 ,0,180*dd4hep::deg), dd4hep::Translation3D(0, 0, -dim.z_offset())); + dd4hep::PlacedVolume envelopeNegativePhysVol = + endcapsAssembly.placeVolume(envelopeVol, envelopeNegativeVolume_tr); + envelopeNegativePhysVol.addPhysVolID("side", -1); + + dd4hep::DetElement caloNegativeDetElem(caloDetElem, "negative", 0); + caloNegativeDetElem.setPlacement(envelopeNegativePhysVol); + + dd4hep::Volume motherVol = aLcdd.pickMotherVolume(caloDetElem); + dd4hep::PlacedVolume envelopePhysVol = motherVol.placeVolume(endcapsAssembly); + caloDetElem.setPlacement(envelopePhysVol); + envelopePhysVol.addPhysVolID("system", idDet); + return caloDetElem; + } + } +} // namespace det + +DECLARE_DETELEMENT(ECalEndcap_Turbine_o1_v02, det::ECalEndcap_Turbine_o1_v02::createECalEndcapTurbine) diff --git a/detector/calorimeter/README.md b/detector/calorimeter/README.md index ce9c3b869..1306a0a87 100644 --- a/detector/calorimeter/README.md +++ b/detector/calorimeter/README.md @@ -26,9 +26,12 @@ Original version taken from [FCCDetectors](https://github.com/HEP-FCC/FCCDetecto Sub-detector for the ecal endcaps, with the absorbers and readout boards arranged in a "turbine-like" geometry. -### 01_v01 +### o1_v01 Initial implementation. A custom segmentation that creates readout cells and constant radius and z is also included (FCCSWEndcapTurbine_k4geo). +### o1_v02 +Changes wrt o1_v01: Added flexibility to configure the wheels individual (to allow for example the possibiliity of having different blade angles in each). This is still a work in progress, so o1_v01 should be used for integrated tests. + ## HCalTileBarrel This sub-detector makes calorimeter barrel. It is used in ALLEGRO detector concept. diff --git a/detectorSegmentations/include/detectorSegmentations/FCCSWEndcapTurbine_k4geo.h b/detectorSegmentations/include/detectorSegmentations/FCCSWEndcapTurbine_k4geo.h index b4abef119..3b1d93bb6 100644 --- a/detectorSegmentations/include/detectorSegmentations/FCCSWEndcapTurbine_k4geo.h +++ b/detectorSegmentations/include/detectorSegmentations/FCCSWEndcapTurbine_k4geo.h @@ -129,7 +129,9 @@ class FCCSWEndcapTurbine_k4geo : public Segmentation { protected: /// turbine blade angle - double m_bladeAngle; + std::vector m_bladeAngle; + /// least common multiple of number of unit cells + int m_nUnitCellsLeastCommonMultiple; /// the number of bins in phi int m_phiBins; /// the coordinate offset in phi diff --git a/detectorSegmentations/src/FCCSWEndcapTurbine_k4geo.cpp b/detectorSegmentations/src/FCCSWEndcapTurbine_k4geo.cpp index c25d00625..a3cf00b35 100644 --- a/detectorSegmentations/src/FCCSWEndcapTurbine_k4geo.cpp +++ b/detectorSegmentations/src/FCCSWEndcapTurbine_k4geo.cpp @@ -21,15 +21,41 @@ namespace DDSegmentation { registerIdentifier("identifier_side", "Cell ID identifier for side", m_sideID, "side"); dd4hep::Detector* dd4hepgeo = &(dd4hep::Detector::getInstance()); + m_bladeAngle.clear(); + + try { + m_bladeAngle.push_back(dd4hepgeo->constant("BladeAngle1")); + std::cout << "Blade angle 1 is " << m_bladeAngle[0] << std::endl; + } + catch(...) { + std::cout << "Blade angle 1 not found in detector metadata, exiting..." << std::endl; + exit(1); + } + try { + m_bladeAngle.push_back(dd4hepgeo->constant("BladeAngle2")); + std::cout << "Blade angle 2 is " << m_bladeAngle[1] << std::endl; + } + catch(...) { + std::cout << "Blade angle 2 not found in detector metadata, exiting..." << std::endl; + exit(1); + } try { - m_bladeAngle = dd4hepgeo->constant("BladeAngle"); - std::cout << "Blade angle is " << m_bladeAngle << std::endl; + m_bladeAngle.push_back(dd4hepgeo->constant("BladeAngle3")); + std::cout << "Blade angle 3 is " << m_bladeAngle[2] << std::endl; } catch(...) { - std::cout << "Blade angle not found in detector metadata, exiting..." << std::endl; + std::cout << "Blade angle 3 not found in detector metadata, exiting..." << std::endl; exit(1); } + try { + m_nUnitCellsLeastCommonMultiple = dd4hepgeo->constant("nUnitCellsLeastCommonMultiple"); + } + catch(...) { + std::cout << "nUnitCellsLeastCommonMultiple not found in detector metadata, exiting..." << std::endl; + exit(1); + } + } FCCSWEndcapTurbine_k4geo::FCCSWEndcapTurbine_k4geo(const BitFieldCoder* decoder) : Segmentation(decoder) { @@ -49,15 +75,40 @@ namespace DDSegmentation { registerIdentifier("identifier_wheel", "Cell ID identifier for wheel", m_wheelID, "wheel"); registerIdentifier("identifier_module", "Cell ID identifier for module", m_moduleID, "module"); dd4hep::Detector* dd4hepgeo = &(dd4hep::Detector::getInstance()); + + m_bladeAngle.clear(); + try { - m_bladeAngle = dd4hepgeo->constant("BladeAngle"); - std::cout << "Blade angle is " << m_bladeAngle << std::endl; + m_bladeAngle.push_back(dd4hepgeo->constant("BladeAngle1")); + std::cout << "Blade angle 1 is " << m_bladeAngle[0] << std::endl; } catch(...) { - std::cout << "Blade angle not found in detector metadata, exiting..." << std::endl; + std::cout << "Blade angle 1 not found in detector metadata, exiting..." << std::endl; + exit(1); + } + try { + m_bladeAngle.push_back(dd4hepgeo->constant("BladeAngle2")); + std::cout << "Blade angle 2 is " << m_bladeAngle[1] << std::endl; + } + catch(...) { + std::cout << "Blade angle 2 not found in detector metadata, exiting..." << std::endl; + exit(1); + } + try { + m_bladeAngle.push_back(dd4hepgeo->constant("BladeAngle3")); + std::cout << "Blade angle 3 is " << m_bladeAngle[2] << std::endl; + } + catch(...) { + std::cout << "Blade angle 3 not found in detector metadata, exiting..." << std::endl; + exit(1); + } + try { + m_nUnitCellsLeastCommonMultiple = dd4hepgeo->constant("nUnitCellsLeastCommonMultiple"); + } + catch(...) { + std::cout << "nUnitCellsLeastCommonMultiple not found in detector metadata, exiting..." << std::endl; exit(1); } - } @@ -90,12 +141,16 @@ CellID FCCSWEndcapTurbine_k4geo::cellID(const Vector3D& /* localPosition */, con double FCCSWEndcapTurbine_k4geo::phi(const CellID& cID) const { CellID iModule = _decoder->get(cID, m_moduleID); - double phiCent = twopi*(iModule+0.5)/78336; + CellID iWheel = _decoder->get(cID, m_wheelID); + + double phiCent = twopi*(iModule+0.5)/m_nUnitCellsLeastCommonMultiple; double rhoLoc = rho(cID); double zLoc = TMath::Abs(z(cID))-m_offsetZ - 45; // hard-code midpoint in z for now - double zCotBladeAngle = zLoc/TMath::Tan(m_bladeAngle); + + double zCotBladeAngle = zLoc/TMath::Tan(m_bladeAngle[iWheel]); + double x = zCotBladeAngle; double y = TMath::Sqrt(rhoLoc*rhoLoc - x*x); // rotate about z axis by phiCent