diff --git a/CMakeLists.txt b/CMakeLists.txt index e6f765ad1..cf4c4cb24 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -119,6 +119,8 @@ file(GLOB G4sources ./plugins/DRCaloFastSimModel.h ./plugins/DRTubesSDAction.hh ./plugins/DRTubesSDAction.cpp + ./plugins/SCEPCalSDActionDRCrystalHit.cpp + ./plugins/DRCrystalHit.cpp ) if(DD4HEP_USE_PYROOT) @@ -141,7 +143,7 @@ target_include_directories(${PackageName} PRIVATE ${PROJECT_SOURCE_DIR}/detect target_include_directories(${PackageName}G4 PRIVATE ${PROJECT_SOURCE_DIR}/detector/calorimeter/dual-readout-tubes/include ) target_link_libraries(${PackageName} DD4hep::DDCore DD4hep::DDRec DD4hep::DDParsers ROOT::Core detectorSegmentations) -target_link_libraries(${PackageName}G4 DD4hep::DDCore DD4hep::DDRec DD4hep::DDParsers DD4hep::DDG4 ROOT::Core podio::podioRootIO EDM4HEP::edm4hep ${Geant4_LIBRARIES}) +target_link_libraries(${PackageName}G4 DD4hep::DDCore DD4hep::DDRec DD4hep::DDParsers DD4hep::DDG4 ROOT::Core detectorSegmentations podio::podioRootIO EDM4HEP::edm4hep ${Geant4_LIBRARIES}) if(K4GEO_USE_LCIO) target_link_libraries(${PackageName} LCIO::lcio) diff --git a/FCCee/IDEA/compact/IDEA_o2_v01/DectDimensions_IDEA_o2_v01.xml b/FCCee/IDEA/compact/IDEA_o2_v01/DectDimensions_IDEA_o2_v01.xml index d8e9a9470..ced713ca0 100644 --- a/FCCee/IDEA/compact/IDEA_o2_v01/DectDimensions_IDEA_o2_v01.xml +++ b/FCCee/IDEA/compact/IDEA_o2_v01/DectDimensions_IDEA_o2_v01.xml @@ -269,7 +269,7 @@ - + diff --git a/FCCee/IDEA/compact/IDEA_o2_v01/IDEA_o2_v01.xml b/FCCee/IDEA/compact/IDEA_o2_v01/IDEA_o2_v01.xml index eb4f52957..40edfbd2e 100644 --- a/FCCee/IDEA/compact/IDEA_o2_v01/IDEA_o2_v01.xml +++ b/FCCee/IDEA/compact/IDEA_o2_v01/IDEA_o2_v01.xml @@ -58,6 +58,9 @@ + + + diff --git a/FCCee/IDEA/compact/IDEA_o2_v01/SCEPCal.xml b/FCCee/IDEA/compact/IDEA_o2_v01/SCEPCal.xml new file mode 100644 index 000000000..fdba2e3c1 --- /dev/null +++ b/FCCee/IDEA/compact/IDEA_o2_v01/SCEPCal.xml @@ -0,0 +1,1050 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + system:4,eta:11,phi:11,depth:4 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/detector/calorimeter/SCEPCalConstructor.cpp b/detector/calorimeter/SCEPCalConstructor.cpp new file mode 100644 index 000000000..633f38aee --- /dev/null +++ b/detector/calorimeter/SCEPCalConstructor.cpp @@ -0,0 +1,579 @@ +//=============================== +// Author: Wonyong Chung +// Princeton University +//=============================== + +#include "detectorSegmentations/SCEPCalSegmentation_k4geo.h" +#include "DD4hep/DetFactoryHelper.h" +#include "DD4hep/DetectorTools.h" +#include "DD4hep/Printout.h" +#include "DD4hep/Detector.h" +#include "DDRec/DetectorData.h" +#include "TGeoTrd2.h" +#include + +using dd4hep::Transform3D; +using dd4hep::RotationZYX; +using dd4hep::RotationY; +using ROOT::Math::RotationZ; +using dd4hep::Rotation3D; +using dd4hep::Position; + +static dd4hep::Ref_t +create_detector_SCEPCal(dd4hep::Detector &theDetector,xml_h xmlElement,dd4hep::SensitiveDetector sens) { + xml_det_t detectorXML =xmlElement; + xml_comp_t dimXML =detectorXML.child(_Unicode(dim)); + xml_comp_t timingXML =detectorXML.child(_Unicode(timing)); + xml_comp_t barrelXML =detectorXML.child(_Unicode(barrel)); + xml_comp_t endcapXML =detectorXML.child(_Unicode(endcap)); + xml_comp_t projFXML =detectorXML.child(_Unicode(projF)); + xml_comp_t projRXML =detectorXML.child(_Unicode(projR)); + xml_comp_t crystalFXML =detectorXML.child(_Unicode(crystalF)); + xml_comp_t crystalRXML =detectorXML.child(_Unicode(crystalR)); + xml_comp_t timingTrXML =detectorXML.child(_Unicode(timingLayerTr)); + xml_comp_t timingLgXML =detectorXML.child(_Unicode(timingLayerLg)); + xml_comp_t sipmLgXML =detectorXML.child(_Unicode(sipmLg)); + xml_comp_t sipmTrXML =detectorXML.child(_Unicode(sipmTr)); + // xml_comp_t instXML =detectorXML.child(_Unicode(inst)); + xml_comp_t timingAssemblyGlobalVisXML =detectorXML.child(_Unicode(timingAssemblyGlobalVis)); + xml_comp_t barrelAssemblyGlobalVisXML =detectorXML.child(_Unicode(barrelAssemblyGlobalVis)); + xml_comp_t endcapAssemblyGlobalVisXML =detectorXML.child(_Unicode(endcapAssemblyGlobalVis)); + xml_comp_t scepcalAssemblyXML =detectorXML.child(_Unicode(scepcalAssembly)); + dd4hep::Material crystalFMat =theDetector.material(crystalFXML.materialStr()); + dd4hep::Material crystalRMat =theDetector.material(crystalRXML.materialStr()); + dd4hep::Material timingTrMat =theDetector.material(timingTrXML.materialStr()); + dd4hep::Material timingLgMat =theDetector.material(timingLgXML.materialStr()); + dd4hep::Material sipmLgMat =theDetector.material(sipmLgXML.materialStr()); + dd4hep::Material sipmTrMat =theDetector.material(sipmTrXML.materialStr()); + // dd4hep::Material instMat =theDetector.material(instXML.materialStr()); + const double EBz =dimXML.attr(_Unicode(barrelHalfZ)); + const double Rin =dimXML.attr(_Unicode(barrelInnerR)); + const double nomfw =dimXML.attr(_Unicode(crystalFaceWidthNominal)); + const double Fdz =dimXML.attr(_Unicode(crystalFlength)); + const double Rdz =dimXML.attr(_Unicode(crystalRlength)); + const double nomth =dimXML.attr(_Unicode(crystalTimingThicknessNominal)); + const double sipmth =dimXML.attr(_Unicode(sipmThickness)); + const int PHI_SEGMENTS =dimXML.attr(_Unicode(phiSegments)); + const int N_PROJECTIVE_FILL =dimXML.attr(_Unicode(projectiveFill)); + const bool CONSTRUCT_TIMING =timingXML.attr(_Unicode(construct)); + // const int TIMING_PHI_START =timingXML.attr(_Unicode(phistart)); + // const int TIMING_PHI_END =timingXML.attr(_Unicode(phiend)); + const bool CONSTRUCT_BARREL =barrelXML.attr(_Unicode(construct)); + const int BARREL_PHI_START =barrelXML.attr(_Unicode(phistart)); + const int BARREL_PHI_END =barrelXML.attr(_Unicode(phiend)); + const bool CONSTRUCT_ENDCAP =endcapXML.attr(_Unicode(construct)); + const int ENDCAP_PHI_START =endcapXML.attr(_Unicode(phistart)); + const int ENDCAP_PHI_END =endcapXML.attr(_Unicode(phiend)); + const int ENDCAP_THETA_START =endcapXML.attr(_Unicode(thetastart)); + const double D_PHI_GLOBAL =2*M_PI/PHI_SEGMENTS; + const double PROJECTIVE_GAP =(N_PROJECTIVE_FILL*nomfw)/2; + double THETA_SIZE_BARREL =atan(EBz/Rin); + double THETA_SIZE_ENDCAP =atan(Rin/EBz); + int N_THETA_BARREL =2*floor(EBz/nomfw); + int N_THETA_ENDCAP =floor(Rin/nomfw); + double D_THETA_BARREL =(M_PI-2*THETA_SIZE_ENDCAP)/(N_THETA_BARREL); + double D_THETA_ENDCAP =THETA_SIZE_ENDCAP/N_THETA_ENDCAP; + int N_PHI_BARREL_CRYSTAL =floor(2*M_PI*Rin/(PHI_SEGMENTS*nomfw)); + double D_PHI_BARREL_CRYSTAL =D_PHI_GLOBAL/N_PHI_BARREL_CRYSTAL; + double thC_end =THETA_SIZE_ENDCAP+D_THETA_BARREL/2; + double r0slice_end =Rin/sin(thC_end); + double z0slice_end =r0slice_end*cos(thC_end)+PROJECTIVE_GAP; + double y0slice_end =r0slice_end*tan(D_THETA_BARREL/2.); + double slice_front_jut =y0slice_end*sin(M_PI/2-thC_end); + double slice_side_jut =y0slice_end*cos(M_PI/2-thC_end); + double z1slice =Rin-slice_front_jut; + double z2slice =Rin+Fdz+Rdz+slice_front_jut; + double zheight_slice =(z2slice-z1slice)/2; + double y1slice =z1slice*tan(M_PI/2-THETA_SIZE_ENDCAP)+PROJECTIVE_GAP; + double y2slice =z2slice*tan(M_PI/2-THETA_SIZE_ENDCAP)+PROJECTIVE_GAP; + double rT =z1slice-2*nomth; + double wT =rT *tan(D_PHI_GLOBAL/2); + int nTiles =ceil(y1slice/wT); + double lT =2*y1slice/nTiles; + int nCy =floor(lT/nomth); + double actY =lT/nCy; + double actX =2*wT/nCy; + // double r2slice_end =r0slice_end+Fdz+Rdz; + // double z2slice_end =r2slice_end*cos(thC_end)+PROJECTIVE_GAP; + // double Rin2slice_end =r2slice_end*sin(thC_end); + // double y2slice_end =r2slice_end*tan(D_THETA_BARREL/2.); + // double slice_front_jut2 =y2slice_end*sin(M_PI/2-thC_end); + // double slice_side_jut2 =y2slice_end*cos(M_PI/2-thC_end); + double barrelSlice_z1 =z0slice_end+slice_side_jut; + double barrelSlice_z2 =y2slice; + // double barrelSlice_rmin2 =Rin2slice_end-slice_front_jut2; + double thCEnd =THETA_SIZE_ENDCAP-D_THETA_ENDCAP/2; + double thCBeg =D_THETA_ENDCAP/2+ENDCAP_THETA_START*D_THETA_ENDCAP; + double r0eEnd =EBz/cos(thCEnd); + // double r2eEnd =r0eEnd+Fdz+Rdz; + // double y2eEnd =r2eEnd*tan(D_THETA_ENDCAP/2.); + double r0eBeg =EBz/cos(thCBeg); + double r2eBeg =r0eBeg+Fdz+Rdz; + double y2eBeg =r2eBeg*tan(D_THETA_ENDCAP/2.); + double aEnd =r0eEnd/cos(D_THETA_ENDCAP/2); + // double bEnd =sqrt(r2eEnd*r2eEnd+y2eEnd*y2eEnd); + double z1End =aEnd*cos(thCEnd+D_THETA_ENDCAP/2); + // double z2End =bEnd*cos(thCEnd-D_THETA_ENDCAP/2); + // double aBeg =r0eBeg/cos(D_THETA_ENDCAP/2); + double bBeg =sqrt(r2eBeg*r2eBeg+y2eBeg*y2eBeg); + // double z1Beg =aBeg*cos(thCBeg+D_THETA_ENDCAP/2); + double z2Beg =bBeg*cos(thCBeg-D_THETA_ENDCAP/2); + double z1rmaxE =z1End*tan(thCEnd+D_THETA_ENDCAP/2); + double z2rmaxE =z2Beg*tan(thCEnd+D_THETA_ENDCAP/2); + double z1rminB =z1End*tan(thCBeg-D_THETA_ENDCAP/2); + double z2rminB =z2Beg*tan(thCBeg-D_THETA_ENDCAP/2); + + std::cout << std::endl; + std::cout << "=GEOMETRY INPUTS=" << std::endl; + std::cout << "BARREL_HALF_Z: " << EBz << std::endl; + std::cout << "BARREL_INNER_R: " << Rin << std::endl; + std::cout << "PHI_SEGMENTS: " << PHI_SEGMENTS << std::endl; + std::cout << "N_PROJECTIVE_FILL: " << N_PROJECTIVE_FILL << std::endl; + std::cout << "F_CRYSTAL_LENGTH: " << Fdz << std::endl; + std::cout << "R_CRYSTAL_LENGTH: " << Rdz << std::endl; + std::cout << "XTAL_FACE_NOM: " << nomfw << std::endl; + std::cout << "TIMING_THICK_NOM: " << nomth << std::endl; + std::cout << "SIPM_THICKNESS: " << sipmth << std::endl; + std::cout << std::endl; + std::cout << "=CONTROL=" << std::endl; + std::cout << "CONSTRUCT_TIMING: " << CONSTRUCT_TIMING << std::endl; + std::cout << "CONSTRUCT_BARREL: " << CONSTRUCT_BARREL << std::endl; + std::cout << "CONSTRUCT_ENDCAP: " << CONSTRUCT_ENDCAP << std::endl; + std::cout << "BARREL_PHI_START: " << BARREL_PHI_START << std::endl; + std::cout << "BARREL_PHI_END : " << BARREL_PHI_END << std::endl; + std::cout << "ENDCAP_PHI_START: " << ENDCAP_PHI_START << std::endl; + std::cout << "ENDCAP_PHI_END : " << ENDCAP_PHI_END << std::endl; + std::cout << "ENDCAP_THETA_START: " << ENDCAP_THETA_START << std::endl; + std::cout << std::endl; + std::cout << "=CALCULATED PARAMETERS=" << std::endl; + std::cout << "D_PHI_GLOBAL: " << D_PHI_GLOBAL << std::endl; + std::cout << "PROJECTIVE_GAP: " << PROJECTIVE_GAP << std::endl; + std::cout << std::endl; + std::cout << "=PROJECTIVE LAYER=" << std::endl; + std::cout << "THETA_SIZE_BARREL: " << THETA_SIZE_BARREL << std::endl; + std::cout << "N_THETA_BARREL: " << N_THETA_BARREL << std::endl; + std::cout << "D_THETA_BARREL: " << D_THETA_BARREL << std::endl; + std::cout << "THETA_SIZE_ENDCAP: " << THETA_SIZE_ENDCAP << std::endl; + std::cout << "N_THETA_ENDCAP: " << N_THETA_ENDCAP << std::endl; + std::cout << "D_THETA_ENDCAP: " << D_THETA_ENDCAP << std::endl; + std::cout << "N_PHI_BARREL_CRYSTAL: " << N_PHI_BARREL_CRYSTAL << std::endl; + std::cout << "D_PHI_BARREL_CRYSTAL: " << D_PHI_BARREL_CRYSTAL << std::endl; + std::cout << std::endl; + std::cout << "=TIMING LAYER=" << std::endl; + std::cout << "N_TILES: " << nTiles << std::endl; + std::cout << "N_XTALS_TILE: " << nCy << std::endl; + std::cout << "LENGTH_TILE: " << lT << std::endl; + std::cout << "WIDTH_TILE: " << 2*wT << std::endl; + std::cout << "ACTUAL_X: " << actX << std::endl; + std::cout << "ACTUAL_Y: " << actY << std::endl; + std::cout << std::endl; + + dd4hep::DetElement ScepcalDetElement(detectorXML.nameStr(),detectorXML.id()); + dd4hep::Volume experimentalHall=theDetector.pickMotherVolume(ScepcalDetElement); + dd4hep::xml::Dimension sdType=detectorXML.child(_Unicode(sensitive)); + sens.setType(sdType.typeStr()); + dd4hep::Readout readout=sens.readout(); + dd4hep::Segmentation geomseg=readout.segmentation(); + dd4hep::Segmentation* _geoSeg=&geomseg; + auto segmentation=dynamic_cast(_geoSeg->segmentation()); + segmentation->setGeomParams(Fdz,Rdz,nomfw,nomth,EBz,Rin,sipmth,PHI_SEGMENTS,N_PROJECTIVE_FILL); + + std::vector zTimingPolyhedra ={-barrelSlice_z1,barrelSlice_z1}; + std::vector rminTimingPolyhedra={rT,rT,}; + std::vector rmaxTimingPolyhedra={z1slice,z1slice}; + std::vector zBarrelPolyhedra ={-barrelSlice_z2,-barrelSlice_z1,barrelSlice_z1,barrelSlice_z2}; + std::vector rminBarrelPolyhedra={z2slice,z1slice,z1slice,z2slice}; + std::vector rmaxBarrelPolyhedra={z2slice,z2slice,z2slice,z2slice}; + std::vector zEndcapPolyhedra ={z1End+PROJECTIVE_GAP,z2Beg+PROJECTIVE_GAP}; + std::vector rminEndcapPolyhedra={z1rminB,z2rminB}; + std::vector rmaxEndcapPolyhedra={z1rmaxE,z2rmaxE}; + std::vector zEndcap1Polyhedra ={-(z2Beg+PROJECTIVE_GAP),-(z1End+PROJECTIVE_GAP)}; + std::vector rminEndcap1Polyhedra={z2rminB,z1rminB}; + std::vector rmaxEndcap1Polyhedra={z2rmaxE,z1rmaxE}; + + dd4hep::Polyhedra timingAssemblyShape(PHI_SEGMENTS,D_PHI_GLOBAL/2,2*M_PI,zTimingPolyhedra,rminTimingPolyhedra,rmaxTimingPolyhedra); + dd4hep::Volume timingAssemblyVol("timingAssemblyVol",timingAssemblyShape,theDetector.material("Vacuum")); + timingAssemblyVol.setVisAttributes(theDetector,timingAssemblyGlobalVisXML.visStr()); + dd4hep::Polyhedra barrelAssemblyShape(PHI_SEGMENTS,D_PHI_GLOBAL/2,2*M_PI,zBarrelPolyhedra,rminBarrelPolyhedra,rmaxBarrelPolyhedra); + dd4hep::Volume barrelAssemblyVol("barrelAssemblyVol",barrelAssemblyShape,theDetector.material("Vacuum")); + barrelAssemblyVol.setVisAttributes(theDetector,barrelAssemblyGlobalVisXML.visStr()); + dd4hep::Polyhedra endcapAssemblyShape(PHI_SEGMENTS,D_PHI_GLOBAL/2,2*M_PI,zEndcapPolyhedra,rminEndcapPolyhedra,rmaxEndcapPolyhedra); + dd4hep::Volume endcapAssemblyVol("endcapAssemblyVol",endcapAssemblyShape,theDetector.material("Vacuum")); + endcapAssemblyVol.setVisAttributes(theDetector,endcapAssemblyGlobalVisXML.visStr()); + dd4hep::Polyhedra endcap1AssemblyShape(PHI_SEGMENTS,D_PHI_GLOBAL/2,2*M_PI,zEndcap1Polyhedra,rminEndcap1Polyhedra,rmaxEndcap1Polyhedra); + dd4hep::Volume endcap1AssemblyVol("endcap1AssemblyVol",endcap1AssemblyShape,theDetector.material("Vacuum")); + endcap1AssemblyVol.setVisAttributes(theDetector,endcapAssemblyGlobalVisXML.visStr()); + auto timingAssemblyVolId =segmentation->setVolumeID(6,0,0,0); + int timingAssemblyVolId32=segmentation->getFirst32bits(timingAssemblyVolId); + auto barrelAssemblyVolId =segmentation->setVolumeID(4,0,0,0); + int barrelAssemblyVolId32=segmentation->getFirst32bits(barrelAssemblyVolId); + auto endcapAssemblyVolId =segmentation->setVolumeID(5,0,0,0); + int endcapAssemblyVolId32=segmentation->getFirst32bits(endcapAssemblyVolId); + auto endcap1AssemblyVolId =segmentation->setVolumeID(5,1,0,0); + int endcap1AssemblyVolId32=segmentation->getFirst32bits(endcap1AssemblyVolId); + experimentalHall.placeVolume(timingAssemblyVol,timingAssemblyVolId32); + dd4hep::PlacedVolume barrelPlacedVol =experimentalHall.placeVolume(barrelAssemblyVol,barrelAssemblyVolId32); + experimentalHall.placeVolume(endcapAssemblyVol,endcapAssemblyVolId32); + experimentalHall.placeVolume(endcap1AssemblyVol,endcap1AssemblyVolId32); + + barrelPlacedVol.addPhysVolID("system", detectorXML.id()); + + ScepcalDetElement.setPlacement(barrelPlacedVol); + + int numCrystalsBarrel = 0; + int numCrystalsEndcap = 0; + int numCrystalsTiming = 0; + + for (int iPhi=CONSTRUCT_BARREL? BARREL_PHI_START:BARREL_PHI_END;iPhisetVolumeID(6,nTile*nCy+nC ,iPhi,3); + auto timingTrId64=segmentation->setVolumeID(6,nTile*nCy+nC ,iPhi,6); + int timingLgId32=segmentation->getFirst32bits(timingLgId64); + int timingTrId32=segmentation->getFirst32bits(timingTrId64); + dd4hep::PlacedVolume timingLgp=tileAssemblyVolume.placeVolume(timingCrystalLgVol,timingLgId32,dispLg); + dd4hep::PlacedVolume timingTrp=tileAssemblyVolume.placeVolume(timingCrystalTrVol,timingTrId32,dispTr); + timingLgp.addPhysVolID("system",6); + timingLgp.addPhysVolID("eta",nTile*nCy+nC); + timingLgp.addPhysVolID("phi",iPhi); + timingLgp.addPhysVolID("depth",3); + timingTrp.addPhysVolID("system",6); + timingTrp.addPhysVolID("eta",nTile*nCy+nC); + timingTrp.addPhysVolID("phi",iPhi); + timingTrp.addPhysVolID("depth",6); + auto sipmLgId64_1=segmentation->setVolumeID(6,nTile*nCy+nC ,iPhi,4); + auto sipmLgId64_2=segmentation->setVolumeID(6,nTile*nCy+nC ,iPhi,5); + auto sipmTrId64_1=segmentation->setVolumeID(6,nTile*nCy+nC ,iPhi,7); + auto sipmTrId64_2=segmentation->setVolumeID(6,nTile*nCy+nC ,iPhi,8); + int sipmLgId32_1=segmentation->getFirst32bits(sipmLgId64_1); + int sipmLgId32_2=segmentation->getFirst32bits(sipmLgId64_2); + int sipmTrId32_1=segmentation->getFirst32bits(sipmTrId64_1); + int sipmTrId32_2=segmentation->getFirst32bits(sipmTrId64_2); + dd4hep::PlacedVolume sipmLgp1=tileAssemblyVolume.placeVolume(sipmBoxLgVol,sipmLgId32_1,dispLg+dispSipmLg); + dd4hep::PlacedVolume sipmLgp2=tileAssemblyVolume.placeVolume(sipmBoxLgVol,sipmLgId32_2,dispLg-dispSipmLg); + dd4hep::PlacedVolume sipmTrp1=tileAssemblyVolume.placeVolume(sipmBoxTrVol,sipmTrId32_1,dispTr+dispSipmTr); + dd4hep::PlacedVolume sipmTrp2=tileAssemblyVolume.placeVolume(sipmBoxTrVol,sipmTrId32_2,dispTr-dispSipmTr); + sipmLgp1.addPhysVolID("system",6); + sipmLgp1.addPhysVolID("eta",nTile*nCy+nC); + sipmLgp1.addPhysVolID("phi",iPhi); + sipmLgp1.addPhysVolID("depth",4); + sipmLgp2.addPhysVolID("system",6); + sipmLgp2.addPhysVolID("eta",nTile*nCy+nC); + sipmLgp2.addPhysVolID("phi",iPhi); + sipmLgp2.addPhysVolID("depth",5); + sipmTrp1.addPhysVolID("system",6); + sipmTrp1.addPhysVolID("eta",nTile*nCy+nC); + sipmTrp1.addPhysVolID("phi",iPhi); + sipmTrp1.addPhysVolID("depth",7); + sipmTrp2.addPhysVolID("system",6); + sipmTrp2.addPhysVolID("eta",nTile*nCy+nC); + sipmTrp2.addPhysVolID("phi",iPhi); + sipmTrp2.addPhysVolID("depth",8); + + numCrystalsTiming+=2; + } + } + + for (int nGamma=0; nGamma0?-1:1; + RotationZYX rot(M_PI/2,-M_PI/2+thC,0); + Position dispF(-rF*cos(thC)+projective_sign*PROJECTIVE_GAP,0,-(rSlice-rF*sin(thC))); + Position dispR(-rR*cos(thC)+projective_sign*PROJECTIVE_GAP,0,-(rSlice-rR*sin(thC))); + + dd4hep::EightPointSolid crystalFShape(Fdz/2,verticesF); + dd4hep::EightPointSolid crystalRShape(Rdz/2,verticesR); + dd4hep::Volume crystalFVol("BarrelCrystalF",crystalFShape,crystalFMat); + dd4hep::Volume crystalRVol("BarrelCrystalR",crystalRShape,crystalRMat); + crystalFVol.setVisAttributes(theDetector,crystalFXML.visStr()); + crystalRVol.setVisAttributes(theDetector,crystalRXML.visStr()); + crystalFVol.setSensitiveDetector(sens); + crystalRVol.setSensitiveDetector(sens); + auto crystalFId64=segmentation->setVolumeID(4,N_THETA_ENDCAP+iTheta ,iPhi*N_PHI_BARREL_CRYSTAL+nGamma,1); + auto crystalRId64=segmentation->setVolumeID(4,N_THETA_ENDCAP+iTheta ,iPhi*N_PHI_BARREL_CRYSTAL+nGamma,2); + int crystalFId32=segmentation->getFirst32bits(crystalFId64); + int crystalRId32=segmentation->getFirst32bits(crystalRId64); + dd4hep::PlacedVolume crystalFp=barrelPhiAssemblyVolume.placeVolume(crystalFVol,crystalFId32,Transform3D(rot,dispF)); + dd4hep::PlacedVolume crystalRp=barrelPhiAssemblyVolume.placeVolume(crystalRVol,crystalRId32,Transform3D(rot,dispR)); + crystalFp.addPhysVolID("system",4); + crystalFp.addPhysVolID("eta",N_THETA_ENDCAP+iTheta); + crystalFp.addPhysVolID("phi",iPhi*N_PHI_BARREL_CRYSTAL+nGamma); + crystalFp.addPhysVolID("depth",1); + crystalRp.addPhysVolID("system",4); + crystalRp.addPhysVolID("eta",N_THETA_ENDCAP+iTheta); + crystalRp.addPhysVolID("phi",iPhi*N_PHI_BARREL_CRYSTAL+nGamma); + crystalRp.addPhysVolID("depth",2); + + numCrystalsBarrel+=2; + + } + + for (int iTheta=0; iThetasetVolumeID(7,iTheta,iPhi*N_PHI_BARREL_CRYSTAL+nGamma,1); + auto crystalRId64=segmentation->setVolumeID(7,iTheta,iPhi*N_PHI_BARREL_CRYSTAL+nGamma,2); + int crystalFId32=segmentation->getFirst32bits(crystalFId64); + int crystalRId32=segmentation->getFirst32bits(crystalRId64); + dd4hep::PlacedVolume crystalFp=barrelPhiAssemblyVolume.placeVolume(crystalFVol,crystalFId32,Transform3D(rot,dispF)); + dd4hep::PlacedVolume crystalRp=barrelPhiAssemblyVolume.placeVolume(crystalRVol,crystalRId32,Transform3D(rot,dispR)); + crystalFp.addPhysVolID("system",7); + crystalFp.addPhysVolID("eta",iTheta); + crystalFp.addPhysVolID("phi",iPhi*N_PHI_BARREL_CRYSTAL+nGamma); + crystalFp.addPhysVolID("depth",1); + crystalRp.addPhysVolID("system",7); + crystalRp.addPhysVolID("eta",iTheta); + crystalRp.addPhysVolID("phi",iPhi*N_PHI_BARREL_CRYSTAL+nGamma); + crystalRp.addPhysVolID("depth",2); + + numCrystalsBarrel+=2; + } + } + } + + for (int iTheta=CONSTRUCT_ENDCAP? ENDCAP_THETA_START:N_THETA_ENDCAP;iTheta zPolyhedra={z1+PROJECTIVE_GAP,z2+PROJECTIVE_GAP}; + std::vector rminPolyhedra={r1min,r2min}; + std::vector rmaxPolyhedra={r1max,r2max}; + dd4hep::Polyhedra endcapRingAssemblyShape(PHI_SEGMENTS,D_PHI_GLOBAL/2,2*M_PI,zPolyhedra,rminPolyhedra,rmaxPolyhedra); + dd4hep::Volume endcapRingAssemblyVolume("endcapRingAssembly",endcapRingAssemblyShape,theDetector.material("Vacuum")); + endcapRingAssemblyVolume.setVisAttributes(theDetector,scepcalAssemblyXML.visStr()); + dd4hep::Volume endcap1RingAssemblyVolume("endcapRingAssembly",endcapRingAssemblyShape,theDetector.material("Vacuum")); + endcap1RingAssemblyVolume.setVisAttributes(theDetector,scepcalAssemblyXML.visStr()); + RotationY rotMirror(M_PI); + endcapAssemblyVol.placeVolume(endcapRingAssemblyVolume); + endcap1AssemblyVol.placeVolume(endcap1RingAssemblyVolume,Transform3D(rotMirror)); + + for (int iPhi=ENDCAP_PHI_START;iPhisetVolumeID(5,iTheta,iPhi*nPhiEndcapCrystal+nGamma,1); + auto crystalRId64=segmentation->setVolumeID(5,iTheta,iPhi*nPhiEndcapCrystal+nGamma,2); + int crystalFId32=segmentation->getFirst32bits(crystalFId64); + int crystalRId32=segmentation->getFirst32bits(crystalRId64); + dd4hep::PlacedVolume crystalFp=endcapRingAssemblyVolume.placeVolume(crystalFVol,crystalFId32,Transform3D(rot,dispF)); + dd4hep::PlacedVolume crystalRp=endcapRingAssemblyVolume.placeVolume(crystalRVol,crystalRId32,Transform3D(rot,dispR)); + crystalFp.addPhysVolID("system",5); + crystalFp.addPhysVolID("eta",iTheta); + crystalFp.addPhysVolID("phi",iPhi*nPhiEndcapCrystal+nGamma); + crystalFp.addPhysVolID("depth",1); + crystalRp.addPhysVolID("system",5); + crystalRp.addPhysVolID("eta",iTheta); + crystalRp.addPhysVolID("phi",iPhi*nPhiEndcapCrystal+nGamma); + crystalRp.addPhysVolID("depth",2); + auto crystalFId641=segmentation->setVolumeID(5,N_THETA_ENDCAP+N_THETA_BARREL+N_THETA_ENDCAP-iTheta,iPhi*nPhiEndcapCrystal+nGamma,1); + auto crystalRId641=segmentation->setVolumeID(5,N_THETA_ENDCAP+N_THETA_BARREL+N_THETA_ENDCAP-iTheta,iPhi*nPhiEndcapCrystal+nGamma,2); + int crystalFId321=segmentation->getFirst32bits(crystalFId641); + int crystalRId321=segmentation->getFirst32bits(crystalRId641); + dd4hep::PlacedVolume crystalFp1=endcap1RingAssemblyVolume.placeVolume(crystalFVol,crystalFId321,Transform3D(rot,dispF)); + dd4hep::PlacedVolume crystalRp1=endcap1RingAssemblyVolume.placeVolume(crystalRVol,crystalRId321,Transform3D(rot,dispR)); + crystalFp1.addPhysVolID("system",5); + crystalFp1.addPhysVolID("eta",N_THETA_ENDCAP+N_THETA_BARREL+N_THETA_ENDCAP-iTheta); + crystalFp1.addPhysVolID("phi",iPhi*nPhiEndcapCrystal+nGamma); + crystalFp1.addPhysVolID("depth",1); + crystalRp1.addPhysVolID("system",5); + crystalRp1.addPhysVolID("eta",N_THETA_ENDCAP+N_THETA_BARREL+N_THETA_ENDCAP-iTheta); + crystalRp1.addPhysVolID("phi",iPhi*nPhiEndcapCrystal+nGamma); + crystalRp1.addPhysVolID("depth",2); + + numCrystalsEndcap+=2; + } + } + } + + std::cout << std::endl; + std::cout << std::endl; + std::cout << "NUM_CRYSTALS_BARREL: " << numCrystalsBarrel << std::endl; + std::cout << "NUM_CRYSTALS_ENDCAP: " << numCrystalsEndcap << std::endl; + std::cout << "NUM_CRYSTALS_TIMING: " << numCrystalsTiming << std::endl; + std::cout << std::endl; + std::cout << std::endl; + + return ScepcalDetElement; +} + +DECLARE_DETELEMENT(SegmentedCrystalECAL,create_detector_SCEPCal) diff --git a/detectorSegmentations/include/detectorSegmentations/SCEPCalSegmentationHandle_k4geo.h b/detectorSegmentations/include/detectorSegmentations/SCEPCalSegmentationHandle_k4geo.h new file mode 100644 index 000000000..31104d251 --- /dev/null +++ b/detectorSegmentations/include/detectorSegmentations/SCEPCalSegmentationHandle_k4geo.h @@ -0,0 +1,87 @@ +//=============================== +// Author: Wonyong Chung +// Princeton University +//=============================== +#ifndef SCEPCalSegmentationHandle_k4geo_h +#define SCEPCalSegmentationHandle_k4geo_h 1 +#include "detectorSegmentations/SCEPCalSegmentation_k4geo.h" +#include "DD4hep/Segmentations.h" +#include "DD4hep/detail/SegmentationsInterna.h" + +namespace dd4hep { + +class Segmentation; +template +class SegmentationWrapper; + +typedef Handle> SCEPCalSegmentationHandle_k4geo; + +class SCEPCalSegmentation_k4geo : public SCEPCalSegmentationHandle_k4geo { +public: + typedef SCEPCalSegmentationHandle_k4geo::Object Object; + +public: + SCEPCalSegmentation_k4geo() = default; + SCEPCalSegmentation_k4geo(const SCEPCalSegmentation_k4geo& e) = default; + SCEPCalSegmentation_k4geo(const Segmentation& e) : Handle(e) {} + SCEPCalSegmentation_k4geo(const Handle& e) : Handle(e) {} + template + SCEPCalSegmentation_k4geo(const Handle& e) : Handle(e) {} + SCEPCalSegmentation_k4geo& operator=(const SCEPCalSegmentation_k4geo& seg) = default; + bool operator==(const SCEPCalSegmentation_k4geo& seg) const { return m_element == seg.m_element; } + + inline Position position(const CellID& id) const { + return Position(access()->implementation->position(id)); + } + + inline dd4hep::CellID cellID(const Position& local, const Position& global, const VolumeID& volID) const { + return access()->implementation->cellID(local, global, volID); + } + + inline VolumeID setVolumeID(int System, int Eta, int Phi, int Depth) const { + return access()->implementation->setVolumeID(System, Eta,Phi,Depth); + } + + inline CellID setCellID(int System, int Eta, int Phi, int Depth) const { + return access()->implementation->setCellID(System, Eta, Phi, Depth); + } + + inline int System(const CellID& aCellID) const { return access()->implementation->System(aCellID); } + inline int Eta(const CellID& aCellID) const { return access()->implementation->Eta(aCellID); } + inline int Phi(const CellID& aCellID) const { return access()->implementation->Phi(aCellID); } + inline int Depth(const CellID& aCellID) const { return access()->implementation->Depth(aCellID); } + + inline int getFirst32bits(const CellID& aCellID) const { return access()->implementation->getFirst32bits(aCellID); } + inline int getLast32bits(const CellID& aCellID) const { return access()->implementation->getLast32bits(aCellID); } + + inline CellID convertFirst32to64(const int aId32) const { return access()->implementation->convertFirst32to64(aId32); } + inline CellID convertLast32to64(const int aId32) const { return access()->implementation->convertLast32to64(aId32); } + + inline int System(const int& aId32) const { return access()->implementation->System(aId32); } + inline int Eta(const int& aId32) const { return access()->implementation->Eta(aId32); } + inline int Phi(const int& aId32) const { return access()->implementation->Phi(aId32); } + inline int Depth(const int& aId32) const { return access()->implementation->Depth(aId32); } + + inline void setGeomParams(double Fdz, double Rdz, + double nomfw, double nomth, + double EBz, double Rin, + double sipmth, + int phiSegments, + int NProjectiveFill) const { + return access()->implementation->setGeomParams(Fdz,Rdz,nomfw,nomth,EBz,Rin,sipmth,phiSegments,NProjectiveFill); + } + + inline double getFdz() const{ return access()->implementation->getFdz(); } + inline double getRdz() const{ return access()->implementation->getRdz(); } + inline double getnomfw() const{ return access()->implementation->getnomfw(); } + inline double getnomth() const{ return access()->implementation->getnomth(); } + inline double getEBz() const{ return access()->implementation->getEBz(); } + inline double getRin() const{ return access()->implementation->getRin(); } + inline double getSipmth() const{ return access()->implementation->getSipmth(); } + inline int getphiSegments() const{ return access()->implementation->getphiSegments(); } + inline int getNProjectiveFill() const{ return access()->implementation->getNProjectiveFill(); } + +}; + +} +#endif diff --git a/detectorSegmentations/include/detectorSegmentations/SCEPCalSegmentation_k4geo.h b/detectorSegmentations/include/detectorSegmentations/SCEPCalSegmentation_k4geo.h new file mode 100644 index 000000000..26f184b29 --- /dev/null +++ b/detectorSegmentations/include/detectorSegmentations/SCEPCalSegmentation_k4geo.h @@ -0,0 +1,98 @@ +//=============================== +// Author: Wonyong Chung +// Princeton University +//=============================== +#ifndef SCEPCalSegmentation_k4geo_h +#define SCEPCalSegmentation_k4geo_h 1 +#include "DDSegmentation/Segmentation.h" +#include "TVector3.h" +#include "DD4hep/DetFactoryHelper.h" +#include +#include + +namespace dd4hep { +namespace DDSegmentation { + +class SCEPCalSegmentation_k4geo : public Segmentation { + public: + SCEPCalSegmentation_k4geo(const std::string& aCellEncoding); + SCEPCalSegmentation_k4geo(const BitFieldCoder* decoder); + virtual ~SCEPCalSegmentation_k4geo() override; + + virtual Vector3D position(const CellID& aCellID) const; + + virtual Vector3D myPosition(const CellID& aCellID) const; + + virtual CellID cellID(const Vector3D& aLocalPosition, + const Vector3D& aGlobalPosition, + const VolumeID& aVolumeID) const; + + VolumeID setVolumeID(int System, int Eta, int Phi, int Depth) const; + CellID setCellID(int System, int Eta, int Phi, int Depth) const; + + int System(const CellID& aCellID) const; + int Eta(const CellID& aCellID) const; + int Phi(const CellID& aCellID) const; + int Depth(const CellID& aCellID) const; + + int getFirst32bits(const CellID& aCellID) const { return (int)aCellID; } + int getLast32bits(const CellID& aCellID) const; + + CellID convertFirst32to64(const int aId32) const { return (CellID)aId32; } + CellID convertLast32to64(const int aId32) const; + + int System(const int& aId32) const { return System( convertFirst32to64(aId32) ); } + int Eta(const int& aId32) const { return Eta( convertFirst32to64(aId32) ); } + int Phi(const int& aId32) const { return Phi( convertFirst32to64(aId32) ); } + int Depth(const int& aId32) const { return Depth( convertFirst32to64(aId32) ); } + + inline void setGeomParams(double Fdz, double Rdz, + double nomfw, double nomth, + double EBz, double Rin, + double sipmth, + int phiSegments, int NProjectiveFill) { + f_Fdz = Fdz; + f_Rdz = Rdz; + f_nomfw = nomfw; + f_nomth = nomth; + f_EBz = EBz; + f_Rin = Rin; + f_sipmth = sipmth; + f_phiSegments = phiSegments; + f_NProjectiveFill = NProjectiveFill; + } + + inline double getFdz() const{ return f_Fdz; } + inline double getRdz() const{ return f_Rdz; } + inline double getnomfw() const{ return f_nomfw; } + inline double getnomth() const{ return f_nomth; } + inline double getEBz() const{ return f_EBz; } + inline double getRin() const{ return f_Rin; } + inline double getSipmth() const{ return f_sipmth; } + inline int getphiSegments() const{ return f_phiSegments; } + inline int getNProjectiveFill() const{ return f_NProjectiveFill; } + + protected: + std::string fSystemId; + std::string fEtaId; + std::string fPhiId; + std::string fDepthId; + + double f_Fdz; + double f_Rdz; + double f_nomfw; + double f_nomth; + double f_EBz; + double f_Rin; + double f_sipmth; + int f_phiSegments; + int f_NProjectiveFill; + + private: + std::unordered_map fPositionOf; + +}; +} +} + +#endif diff --git a/detectorSegmentations/src/SCEPCalSegmentationHandle_k4geo.cpp b/detectorSegmentations/src/SCEPCalSegmentationHandle_k4geo.cpp new file mode 100644 index 000000000..cb73ee5d7 --- /dev/null +++ b/detectorSegmentations/src/SCEPCalSegmentationHandle_k4geo.cpp @@ -0,0 +1,4 @@ +#include "detectorSegmentations/SCEPCalSegmentationHandle_k4geo.h" +#include "DD4hep/detail/Handle.inl" + +DD4HEP_INSTANTIATE_HANDLE_UNNAMED(dd4hep::DDSegmentation::SCEPCalSegmentation_k4geo); diff --git a/detectorSegmentations/src/SCEPCalSegmentation_k4geo.cpp b/detectorSegmentations/src/SCEPCalSegmentation_k4geo.cpp new file mode 100644 index 000000000..71d810f35 --- /dev/null +++ b/detectorSegmentations/src/SCEPCalSegmentation_k4geo.cpp @@ -0,0 +1,252 @@ +//=============================== +// Author: Wonyong Chung +// Princeton University +//=============================== +#include "detectorSegmentations/SCEPCalSegmentation_k4geo.h" +#include +#include +#include + +namespace dd4hep { +namespace DDSegmentation { + +SCEPCalSegmentation_k4geo::SCEPCalSegmentation_k4geo(const std::string& cellEncoding) : Segmentation(cellEncoding) { + _type = "SCEPCalSegmentation_k4geo"; + _description = "SCEPCal segmentation based on side/eta/phi/depth/S/C"; + registerIdentifier("identifier_system", "Cell ID identifier for numSystem", fSystemId, "system"); + registerIdentifier("identifier_eta", "Cell ID identifier for numEta", fEtaId, "eta"); + registerIdentifier("identifier_phi", "Cell ID identifier for numPhi", fPhiId, "phi"); + registerIdentifier("identifier_depth", "Cell ID identifier for numDepth", fDepthId, "depth"); +} + +SCEPCalSegmentation_k4geo::SCEPCalSegmentation_k4geo(const BitFieldCoder* decoder) : Segmentation(decoder) { + _type = "SCEPCalSegmentation_k4geo"; + _description = "SCEPCal segmentation based on side/eta/phi/depth/S/C"; + registerIdentifier("identifier_system", "Cell ID identifier for numSystem", fSystemId, "system"); + registerIdentifier("identifier_eta", "Cell ID identifier for Eta", fEtaId, "eta"); + registerIdentifier("identifier_phi", "Cell ID identifier for Phi", fPhiId, "phi"); + registerIdentifier("identifier_depth", "Cell ID identifier for Depth", fDepthId, "depth"); +} + +SCEPCalSegmentation_k4geo::~SCEPCalSegmentation_k4geo() {} + +Vector3D SCEPCalSegmentation_k4geo::position(const CellID& cID) const { + return myPosition(cID); +}; + +Vector3D SCEPCalSegmentation_k4geo::myPosition(const CellID& cID) const { + + int copyNum = (int)cID; + + double Fdz =getFdz(); + double Rdz =getRdz(); + double nomfw =getnomfw(); + double nomth =getnomth(); + double EBz =getEBz(); + double Rin =getRin(); + double sipmth =getSipmth(); + int PHI_SEGMENTS =getphiSegments(); + int N_PROJECTIVE_FILL=getNProjectiveFill(); + int system =System(copyNum); + int nEta_in =Eta(copyNum); + int nPhi_in =Phi(copyNum); + int nDepth_in =Depth(copyNum); + + double D_PHI_GLOBAL =2*M_PI/PHI_SEGMENTS; + double PROJECTIVE_GAP =(N_PROJECTIVE_FILL*nomfw)/2; + // double THETA_SIZE_BARREL =atan(EBz/Rin); + double THETA_SIZE_ENDCAP =atan(Rin/EBz); + int N_THETA_BARREL =2*floor(EBz/nomfw); + int N_THETA_ENDCAP =floor(Rin/nomfw); + double D_THETA_BARREL =(M_PI-2*THETA_SIZE_ENDCAP)/(N_THETA_BARREL); + double D_THETA_ENDCAP =THETA_SIZE_ENDCAP/N_THETA_ENDCAP; + int N_PHI_BARREL_CRYSTAL=floor(2*M_PI*Rin/(PHI_SEGMENTS*nomfw)); + double D_PHI_BARREL_CRYSTAL=D_PHI_GLOBAL/N_PHI_BARREL_CRYSTAL; + + // timing + if (system == 6) { + double thC_end =THETA_SIZE_ENDCAP+D_THETA_BARREL/2; + double r0slice_end =Rin/sin(thC_end); + double y0slice_end =r0slice_end*tan(D_THETA_BARREL/2.); + double slice_front_jut =y0slice_end*sin(M_PI/2-thC_end); + double z1slice =Rin -slice_front_jut; + // double z2slice =Rin +Fdz +Rdz +slice_front_jut; + double y1slice =z1slice*tan(M_PI/2-THETA_SIZE_ENDCAP) +PROJECTIVE_GAP; + double phiTiming =nPhi_in*D_PHI_GLOBAL; + double rT =z1slice -2*nomth; + double wT =rT *tan(D_PHI_GLOBAL/2); + int nTiles =ceil(y1slice/wT); + double lT =2*y1slice/nTiles; + int nCy =floor(lT/nomth); + double actY =lT/nCy; + double actX =2*wT/nCy; + double rTimingAssembly =rT+nomth; + int nTile =int(nEta_in/nCy); + ROOT::Math::XYZVector dispTimingAssembly(rTimingAssembly*cos(0),rTimingAssembly*sin(0),0); + ROOT::Math::XYZVector dispTileAssembly(0,0,-y1slice +nTile*lT +lT/2); + int nC=nEta_in-nTile*nCy; + int phiTimingSign=nPhi_in%2==0? 1:-1; + int sign =nTile%2==0? 1:-1; + ROOT::Math::RotationZ rotZ(phiTiming); + ROOT::Math::XYZVector dispLg(sign*phiTimingSign*(nomth/2),-wT +actX/2 + nC*actX,0); + ROOT::Math::XYZVector dispTr(sign*phiTimingSign*(-nomth/2),0,-lT/2 +actY/2 +nC*actY); + ROOT::Math::XYZVector dispSipmLg(0, 0, lT/2-sipmth/2); + ROOT::Math::XYZVector dispSipmTr(0, wT-sipmth/2, 0); + if (nDepth_in==3) {return rotZ*(dispTimingAssembly +dispTileAssembly +dispLg);} + else if (nDepth_in==4) {return rotZ*(dispTimingAssembly +dispTileAssembly +dispLg +dispSipmLg);} + else if (nDepth_in==5) {return rotZ*(dispTimingAssembly +dispTileAssembly +dispLg -dispSipmLg);} + else if (nDepth_in==6) {return rotZ*(dispTimingAssembly +dispTileAssembly +dispTr);} + else if (nDepth_in==7) {return rotZ*(dispTimingAssembly +dispTileAssembly +dispTr +dispSipmTr);} + else if (nDepth_in==8) {return rotZ*(dispTimingAssembly +dispTileAssembly +dispTr -dispSipmTr);} + } + // endcap + else if (system == 5) { + int nTheta; + if (nEta_inN_THETA_ENDCAP+N_THETA_BARREL) { + nTheta = N_THETA_ENDCAP+N_THETA_BARREL+N_THETA_ENDCAP -nEta_in; + } + double thC =D_THETA_ENDCAP/2+nTheta*D_THETA_ENDCAP; + double RinEndcap =EBz*tan(thC); + int nPhiEndcapCrystal =floor(2*M_PI*RinEndcap/(PHI_SEGMENTS*nomfw)); + double dPhiEndcapCrystal =D_PHI_GLOBAL/nPhiEndcapCrystal; + double r0e =RinEndcap/sin(thC); + double r1e =r0e+Fdz; + int nPhi =int(nPhi_in/nPhiEndcapCrystal); + int nGamma =nPhi_in%nPhiEndcapCrystal; + double phi =nPhi*D_PHI_GLOBAL; + double gamma =-D_PHI_GLOBAL/2+dPhiEndcapCrystal/2+dPhiEndcapCrystal*nGamma; + int mirror =nEta_in0? -1:1; + double r0e =Rin/sin(thC); + double r1e =r0e+Fdz; + double rSlice =Rin+(Fdz+Rdz)/2; + Position dispSlice(rSlice*cos(phi),rSlice*sin(phi),0); + ROOT::Math::RotationZ rotZ(phi); + if (nDepth_in==1) { + double rF=r0e+Fdz/2.; + ROOT::Math::XYZVector dispF(rF*sin(thC)-rSlice,rF*sin(thC)*tan(gamma),rF*cos(thC)-projective_sign*PROJECTIVE_GAP); + return (dispSlice+rotZ*dispF); + } + else if (nDepth_in==2) { + double rR=r1e+Rdz/2.; + ROOT::Math::XYZVector dispR(rR*sin(thC)-rSlice,rR*sin(thC)*tan(gamma),rR*cos(thC)-projective_sign*PROJECTIVE_GAP); + return (dispSlice+rotZ*dispR); + } + } + // projective fill + else if (system == 7) { + int nTheta =nEta_in; + int nPhi =int(nPhi_in/N_PHI_BARREL_CRYSTAL); + int nGamma =nPhi_in%N_PHI_BARREL_CRYSTAL; + double phi =nPhi*D_PHI_GLOBAL; + double gamma =-D_PHI_GLOBAL/2+D_PHI_BARREL_CRYSTAL/2+D_PHI_BARREL_CRYSTAL*nGamma; + double thC =M_PI/2; + double r0e =Rin/sin(thC); + double r1e =r0e+Fdz; + double rSlice=Rin +(Fdz+Rdz)/2; + Position dispSlice(rSlice*cos(phi),rSlice*sin(phi),0); + ROOT::Math::RotationZ rotZ(phi); + if (nDepth_in==1) { + double rF=r0e+Fdz/2.; + ROOT::Math::XYZVector dispF(rF*sin(thC)-rSlice,rF*sin(thC)*tan(gamma),rF*cos(thC)-nomfw*(N_PROJECTIVE_FILL-1)/2+nTheta*nomfw); + return (dispSlice+rotZ*dispF); + } + else if (nDepth_in==2) { + double rR=r1e+Rdz/2.; + ROOT::Math::XYZVector dispR(rR*sin(thC)-rSlice,rR*sin(thC)*tan(gamma),rR*cos(thC)-nomfw*(N_PROJECTIVE_FILL-1)/2+nTheta*nomfw); + return (dispSlice+rotZ*dispR); + } + } + return Vector3D(0,0,0); +} + +CellID SCEPCalSegmentation_k4geo::cellID(const Vector3D& /*localPosition*/, + const Vector3D& /*globalPosition*/, + const VolumeID& vID) const { + return setCellID(System(vID), Eta(vID), Phi(vID), Depth(vID) ); +} + +VolumeID SCEPCalSegmentation_k4geo::setVolumeID(int System, int Eta, int Phi, int Depth) const { + VolumeID SystemId = static_cast(System); + VolumeID EtaId = static_cast(Eta); + VolumeID PhiId = static_cast(Phi); + VolumeID DepthId = static_cast(Depth); + VolumeID vID = 0; + _decoder->set(vID, fSystemId, SystemId); + _decoder->set(vID, fEtaId, EtaId); + _decoder->set(vID, fPhiId, PhiId); + _decoder->set(vID, fDepthId, DepthId); + return vID; +} + +CellID SCEPCalSegmentation_k4geo::setCellID(int System, int Eta, int Phi, int Depth) const { + VolumeID SystemId = static_cast(System); + VolumeID EtaId = static_cast(Eta); + VolumeID PhiId = static_cast(Phi); + VolumeID DepthId = static_cast(Depth); + VolumeID vID = 0; + _decoder->set(vID, fSystemId, SystemId); + _decoder->set(vID, fEtaId, EtaId); + _decoder->set(vID, fPhiId, PhiId); + _decoder->set(vID, fDepthId, DepthId); + return vID; +} + +int SCEPCalSegmentation_k4geo::System(const CellID& aCellID) const { + VolumeID System = static_cast(_decoder->get(aCellID, fSystemId)); + return static_cast(System); +} + +int SCEPCalSegmentation_k4geo::Eta(const CellID& aCellID) const { + VolumeID Eta = static_cast(_decoder->get(aCellID, fEtaId)); + return static_cast(Eta); +} + +int SCEPCalSegmentation_k4geo::Phi(const CellID& aCellID) const { + VolumeID Phi = static_cast(_decoder->get(aCellID, fPhiId)); + return static_cast(Phi); +} + +int SCEPCalSegmentation_k4geo::Depth(const CellID& aCellID) const { + VolumeID Depth = static_cast(_decoder->get(aCellID, fDepthId)); + return static_cast(Depth); +} + +int SCEPCalSegmentation_k4geo::getLast32bits(const CellID& aCellID) const { + CellID aId64 = aCellID >> sizeof(int)*CHAR_BIT; + int aId32 = (int)aId64; + return aId32; +} + +CellID SCEPCalSegmentation_k4geo::convertLast32to64(const int aId32) const { + CellID aId64 = (CellID)aId32; + aId64 <<= sizeof(int)*CHAR_BIT; + return aId64; +} + +} +} diff --git a/detectorSegmentations/src/plugins/SegmentationFactories.cpp b/detectorSegmentations/src/plugins/SegmentationFactories.cpp index a3ac9b0e9..a7081097f 100644 --- a/detectorSegmentations/src/plugins/SegmentationFactories.cpp +++ b/detectorSegmentations/src/plugins/SegmentationFactories.cpp @@ -40,3 +40,6 @@ DECLARE_SEGMENTATION(FCCSWHCalPhiTheta_k4geo, create_segmentation) + +#include "detectorSegmentations/SCEPCalSegmentation_k4geo.h" +DECLARE_SEGMENTATION(SCEPCalSegmentation_k4geo, create_segmentation) diff --git a/example/SteeringFile_IDEA_o2_v01.py b/example/SteeringFile_IDEA_o2_v01.py index 9133c506c..4eef726d7 100644 --- a/example/SteeringFile_IDEA_o2_v01.py +++ b/example/SteeringFile_IDEA_o2_v01.py @@ -103,8 +103,10 @@ ## List of patterns matching sensitive detectors of type Calorimeter. SIM.action.calorimeterSDTypes = ["calorimeter"] -## Replace SDAction for DREndcapTubes subdetector +## Replace SDAction for subdetectors SIM.action.mapActions["DREndcapTubes"] = "DRTubesSDAction" +SIM.action.mapActions["SCEPCal"] = "SCEPCalSDAction_DRHit" + ## Configure the regexSD for DREndcapTubes subdetector SIM.geometry.regexSensitiveDetector["DREndcapTubes"] = { "Match": ["DRETS"], @@ -184,7 +186,7 @@ } ## a map between patterns and filter objects, using patterns to attach filters to sensitive detector -SIM.filter.mapDetFilter = {} +SIM.filter.mapDetFilter = {'SCEPCal' : 'edep1kev'} ## default filter for tracking sensitive detectors; this is applied if no other filter is used for a tracker SIM.filter.tracker = "edep1kev" @@ -440,7 +442,25 @@ ## # arbitrary options can be created and set via the steering file or command line ## SIM.outputConfig.myExtension = '.csv' ## - +def setupEDM4hepOutputDR(dd4hepSimulation): + from DDG4 import EventAction, Kernel + dd = dd4hepSimulation + evt_edm4hep = EventAction(Kernel(), 'Geant4Output2EDM4hep_DRC/' + dd.outputFile, True) + evt_edm4hep.Control = True + output = dd.outputFile + if not dd.outputFile.endswith(dd.outputConfig.myExtension): + output = dd.outputFile + dd.outputConfig.myExtension + evt_edm4hep.Output = output + evt_edm4hep.enableUI() + Kernel().eventAction().add(evt_edm4hep) + eventPars = dd.meta.parseEventParameters() + evt_edm4hep.RunHeader = dd.meta.addParametersToRunHeader(dd) + evt_edm4hep.EventParametersString, evt_edm4hep.EventParametersInt, evt_edm4hep.EventParametersFloat = eventPars + evt_edm4hep.RunNumberOffset = dd.meta.runNumberOffset if dd.meta.runNumberOffset > 0 else 0 + evt_edm4hep.EventNumberOffset = dd.meta.eventNumberOffset if dd.meta.eventNumberOffset > 0 else 0 + return None +SIM.outputConfig.userOutputPlugin = setupEDM4hepOutputDR +SIM.outputConfig.myExtension = '.root' ################################################################################ ## Configuration for the Particle Handler/ MCTruth treatment diff --git a/plugins/DRCrystalHit.cpp b/plugins/DRCrystalHit.cpp new file mode 100644 index 000000000..f2dadbfa1 --- /dev/null +++ b/plugins/DRCrystalHit.cpp @@ -0,0 +1,32 @@ +//=============================== +// Author: Wonyong Chung +// Princeton University +//=============================== +#include +#include +#include +// #include +// #include +// #include +// #include +// #include +#include "DRCrystalHit.h" +// #include "G4Track.hh" + +SCEPCal::DRCrystalHit::DRCrystalHit() +: dd4hep::sim::Geant4HitData(), position(), truth(), energyDeposit(0), nCerenkovProd(0), nScintillationProd(0), tAvgC(0), tAvgS(0) { + + dd4hep::InstanceCount::increment(this); + +} + +SCEPCal::DRCrystalHit::DRCrystalHit(const Position& pos) +: dd4hep::sim::Geant4HitData(), position(pos), truth(), energyDeposit(0), nCerenkovProd(0), nScintillationProd(0), tAvgC(0), tAvgS(0) { + + dd4hep::InstanceCount::increment(this); + +} + +SCEPCal::DRCrystalHit::~DRCrystalHit() { + dd4hep::InstanceCount::decrement(this); +} \ No newline at end of file diff --git a/plugins/DRCrystalHit.h b/plugins/DRCrystalHit.h new file mode 100644 index 000000000..ba22b7f6c --- /dev/null +++ b/plugins/DRCrystalHit.h @@ -0,0 +1,47 @@ +//=============================== +// Author: Wonyong Chung +// Princeton University +//=============================== +#ifndef DRCrystalHit_h +#define DRCrystalHit_h 1 +// #include "G4VHit.hh" +// #include "G4THitsCollection.hh" +// #include "G4Allocator.hh" +// #include "G4ThreeVector.hh" +#include "DDG4/Geant4Data.h" +// #include "G4OpticalPhoton.hh" +// #include "G4VProcess.hh" +// #include "DD4hep/Objects.h" +// #include "DD4hep/Segmentations.h" +// #include "CLHEP/Vector/ThreeVector.h" + +namespace SCEPCal { + + typedef ROOT::Math::XYZVector Position; + + class DRCrystalHit : public dd4hep::sim::Geant4HitData { + + public: + typedef dd4hep::sim::Geant4HitData base_t; + + Position position; + Contributions truth; + double energyDeposit; + + int nCerenkovProd; + int nScintillationProd; + double tAvgC; + double tAvgS; + + public: + DRCrystalHit(); + DRCrystalHit(DRCrystalHit&& c) = delete; + DRCrystalHit(const DRCrystalHit& c) = delete; + DRCrystalHit(const Position& cell_pos); + virtual ~DRCrystalHit(); + DRCrystalHit& operator=(DRCrystalHit&& c) = delete; + DRCrystalHit& operator=(const DRCrystalHit& c) = delete; + }; +}; + +#endif diff --git a/plugins/Geant4Output2EDM4hep_DRC.cpp b/plugins/Geant4Output2EDM4hep_DRC.cpp index 6d408424e..2c74dae07 100644 --- a/plugins/Geant4Output2EDM4hep_DRC.cpp +++ b/plugins/Geant4Output2EDM4hep_DRC.cpp @@ -15,6 +15,8 @@ #include #include #include +#include +#include "DRCrystalHit.h" /// podio include files #include #include @@ -51,12 +53,17 @@ namespace dd4hep { using writer_t = podio::ROOTWriter; using stringmap_t = std::map< std::string, std::string >; using trackermap_t = std::map< std::string, edm4hep::SimTrackerHitCollection >; + using calorimeterpair_t = std::pair< edm4hep::SimCalorimeterHitCollection, edm4hep::CaloHitContributionCollection >; using calorimetermap_t = std::map< std::string, calorimeterpair_t >; + using drcalopair_t = std::pair< edm4hep::RawCalorimeterHitCollection, edm4hep::RawTimeSeriesCollection> ; // Required info for IDEA DRC sim hit using drcalomap_t = std::map< std::string, drcalopair_t >; // Required info for IDEA DRC sim hit using drcaloWavmap_t = std::map< std::string, edm4hep::RawTimeSeriesCollection >; + using scepcalcalopair_t = std::pair< edm4hep::SimDRCalorimeterHitCollection, edm4hep::CaloHitContributionCollection >; + using scepcalcalomap_t = std::map< std::string, scepcalcalopair_t >; + std::unique_ptr m_file { }; podio::Frame m_frame { }; edm4hep::MCParticleCollection m_particles { }; @@ -64,6 +71,7 @@ namespace dd4hep { calorimetermap_t m_calorimeterHits; drcalomap_t m_drcaloHits; drcaloWavmap_t m_drcaloWaves; + scepcalcalomap_t m_scepcalcaloHits; stringmap_t m_runHeader; stringmap_t m_eventParametersInt; stringmap_t m_eventParametersFloat; @@ -290,12 +298,18 @@ void Geant4Output2EDM4hep_DRC::commit( OutputContext& /* ctxt */) { for (auto it = m_drcaloWaves.begin(); it != m_drcaloWaves.end(); ++it) { m_frame.put( std::move(it->second), it->first + "WaveLen"); } + for (auto& [colName, calorimeterHits] : m_scepcalcaloHits) { + m_frame.put( std::move(calorimeterHits.first), colName); + m_frame.put( std::move(calorimeterHits.second), colName + "Contributions"); + } + m_file->writeFrame(m_frame, m_section_name); m_particles.clear(); m_trackerHits.clear(); m_calorimeterHits.clear(); m_drcaloHits.clear(); m_drcaloWaves.clear(); + m_scepcalcaloHits.clear(); m_frame = {}; return; } @@ -335,6 +349,7 @@ void Geant4Output2EDM4hep_DRC::begin(const G4Event* event) { m_calorimeterHits.clear(); m_drcaloHits.clear(); m_drcaloWaves.clear(); + m_scepcalcaloHits.clear(); } /// Data conversion interface for MC particles to EDM4hep format @@ -592,6 +607,52 @@ void Geant4Output2EDM4hep_DRC::saveCollection(OutputContext& /*ctxt*/, } //------------------------------------------------------------------- } + else if( typeid( SCEPCal::DRCrystalHit ) == coll->type().type() ){ + + Geant4Sensitive* sd = coll->sensitive(); + int hit_creation_mode = sd->hitCreationMode(); + + auto& hits = m_scepcalcaloHits[colName]; + + for(unsigned i=0 ; i < nhits ; ++i){ + + auto sch = hits.first->create(); + const SCEPCal::DRCrystalHit* hit = coll->hit(i); + + const auto& pos = hit->position; + edm4hep::Vector3f hitpos( float(pos.x()/CLHEP::mm), float(pos.y()/CLHEP::mm), float(pos.z()/CLHEP::mm) ); + + sch.setCellID( hit->cellID ); + sch.setPosition( hitpos ); + sch.setEnergy( hit->energyDeposit/CLHEP::GeV ); + + sch.setNCerenkovProd( hit->nCerenkovProd ); + sch.setNScintillationProd(hit->nScintillationProd ); + + sch.setTAvgS( hit->tAvgC ); + sch.setTAvgS( hit->tAvgS ); + + for(auto ci=hit->truth.begin(); ci != hit->truth.end(); ++ci){ + + auto sCaloHitCont = hits.second->create(); + sch.addToContributions( sCaloHitCont ); + + const Geant4HitData::MonteCarloContrib& c = *ci; + + int trackID = pm->particleID(c.trackID); + auto mcp = m_particles.at(trackID); + + sCaloHitCont.setEnergy( c.deposit/CLHEP::GeV ); + sCaloHitCont.setTime( c.time/CLHEP::ns ); + sCaloHitCont.setParticle( mcp ); + + edm4hep::Vector3f p(c.x/CLHEP::mm, c.y/CLHEP::mm, c.z/CLHEP::mm); + + sCaloHitCont.setPDG( c.pdgID ); + sCaloHitCont.setStepPosition( p ); + } + } + } else if( typeid( Geant4DRCalorimeter::Hit ) == coll->type().type() ){ Geant4Sensitive* sd = coll->sensitive(); int hit_creation_mode = sd->hitCreationMode(); diff --git a/plugins/SCEPCalSDActionDRCrystalHit.cpp b/plugins/SCEPCalSDActionDRCrystalHit.cpp new file mode 100644 index 000000000..54087c413 --- /dev/null +++ b/plugins/SCEPCalSDActionDRCrystalHit.cpp @@ -0,0 +1,107 @@ +//=============================== +// Author: Wonyong Chung +// Princeton University +//=============================== +#include "DRCrystalHit.h" +#include "detectorSegmentations/SCEPCalSegmentation_k4geo.h" +#include "DDG4/Geant4SensDetAction.inl" +#include "DDG4/Factories.h" +#include "G4OpticalPhoton.hh" +#include "G4VProcess.hh" +#include + +namespace SCEPCal { + G4double convertEvtoNm(G4double energy) + { + return 1239.84187/energy*1000.; //GeV to nm + } + class SegmentedCrystalCalorimeterSD_DRHit { + public: + typedef DRCrystalHit Hit; + }; +} + +namespace dd4hep { + namespace sim { + using namespace SCEPCal; + + template <> void Geant4SensitiveAction::defineCollections() { + m_collectionID = declareReadoutFilteredCollection(); + } + template <> bool + Geant4SensitiveAction::process(const G4Step* step,G4TouchableHistory* /*hist*/ ) { + G4double edep = step->GetTotalEnergyDeposit(); + G4StepPoint *thePrePoint = step->GetPreStepPoint(); + G4StepPoint *thePostPoint = step->GetPostStepPoint(); + G4TouchableHandle thePreStepTouchable = thePrePoint->GetTouchableHandle(); + G4VPhysicalVolume *thePrePV = thePrePoint->GetPhysicalVolume(); + G4double pretime = thePrePoint->GetGlobalTime(); + G4VPhysicalVolume *thePostPV = thePostPoint->GetPhysicalVolume(); + G4double posttime = thePostPoint->GetGlobalTime(); + G4String thePrePVName = ""; + if (thePrePV) thePrePVName = thePrePV->GetName(); + G4String thePostPVName = ""; + if (thePostPV) thePostPVName = thePostPV->GetName(); + + Geant4StepHandler h(step); + Geant4HitData::MonteCarloContrib contrib = Geant4HitData::extractContribution(step); + Geant4HitCollection* coll = collection(m_collectionID); + + dd4hep::Segmentation* _geoSeg = &m_segmentation; + auto segmentation=dynamic_cast(_geoSeg->segmentation()); + auto copyNum64 = segmentation->convertFirst32to64(thePreStepTouchable->GetCopyNumber(0)); + int cellID = (int)copyNum64; + + SegmentedCrystalCalorimeterSD_DRHit::Hit* hit = coll->findByKey(cellID); + + if(!hit) { + DDSegmentation::Vector3D pos = segmentation->myPosition(copyNum64); + Position global(pos.x(),pos.y(),pos.z()); + hit = new SegmentedCrystalCalorimeterSD_DRHit::Hit(global); + hit->cellID = cellID; + coll->add(cellID, hit); + } + + G4Track * track = step->GetTrack(); + + if(track->GetDefinition()==G4OpticalPhoton::OpticalPhotonDefinition()) { + + float avgarrival=(pretime+posttime)/2.; + + // count 1st and kill + // apply scale factor and poisson smearing offline + + int phstep = track->GetCurrentStepNumber(); + + if (track->GetCreatorProcess()->G4VProcess::GetProcessName()=="CerenkovPhys") { + if(phstep==1) { + float tAvgC_new = (((hit->tAvgC)*(hit->nCerenkovProd)) +avgarrival)/(hit->nCerenkovProd+1); + hit->nCerenkovProd+=1; + hit->tAvgC = tAvgC_new; + } + track->SetTrackStatus(fStopAndKill); + } + + else if (track->GetCreatorProcess()->G4VProcess::GetProcessName()=="ScintillationPhys") { + if(phstep==1) { + float tAvgS_new = (((hit->tAvgS)*(hit->nScintillationProd)) +avgarrival)/(hit->nScintillationProd+1); + hit->nScintillationProd+=1; + hit->tAvgS = tAvgS_new; + } + track->SetTrackStatus(fStopAndKill); + } + + } + + hit->truth.emplace_back(contrib); + hit->energyDeposit+=edep; + mark(h.track); + return true; + } + } +} + +namespace dd4hep { namespace sim { + typedef Geant4SensitiveAction SCEPCalSDAction_DRHit; + }} +DECLARE_GEANT4SENSITIVE(SCEPCalSDAction_DRHit) \ No newline at end of file