Skip to content

Commit

Permalink
fix bug in calculation of detector geometry information (#433)
Browse files Browse the repository at this point in the history 
* fix bug in calculation of detector geometry information saved for reconstruction

* fix comment
  • Loading branch information
@giovannimarchiori
giovannimarchiori authored Feb 20, 2025
1 parent 91e5123 commit 6a8703e
Show file tree
Hide file tree
Showing 2 changed files with 82 additions and 77 deletions.
68 changes: 35 additions & 33 deletions detector/calorimeter/ECalBarrel_NobleLiquid_InclinedTrapezoids_o1_v03_geo.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -135,17 +135,17 @@ static dd4hep::detail::Ref_t createECalBarrelInclined(dd4hep::Detector& aLcdd,
dd4hep::Tube cryoSideOuterShape(cryoDim.rmin2(), cryoDim.rmax1(), cryoDim.dz());
dd4hep::SubtractionSolid cryoSideShape(cryoSideOuterShape, bathAndServicesOuterShape);
lLog << MSG::INFO
<< "ECAL cryostat: front: rmin (cm) = " << cryoDim.rmin1()/dd4hep::cm
<< " rmax (cm) = " << cryoDim.rmin2()/dd4hep::cm
<< " dz (cm) = " << cryoDim.dz()/dd4hep::cm << endmsg;
<< "ECAL cryostat: front: rmin (cm) = " << cryoDim.rmin1()/dd4hep::cm
<< " rmax (cm) = " << cryoDim.rmin2()/dd4hep::cm
<< " dz (cm) = " << cryoDim.dz()/dd4hep::cm << endmsg;
lLog << MSG::INFO
<< "ECAL cryostat: back: rmin (cm) = " << cryoDim.rmax1()/dd4hep::cm
<< " rmax (cm) = " << cryoDim.rmax2()/dd4hep::cm
<< " dz (cm) = " << cryoDim.dz()/dd4hep::cm << endmsg;
<< "ECAL cryostat: back: rmin (cm) = " << cryoDim.rmax1()/dd4hep::cm
<< " rmax (cm) = " << cryoDim.rmax2()/dd4hep::cm
<< " dz (cm) = " << cryoDim.dz()/dd4hep::cm << endmsg;
lLog << MSG::INFO
<< "ECAL cryostat: side: rmin (cm) = " << cryoDim.rmin2()/dd4hep::cm
<< " rmax (cm) = " << cryoDim.rmax1()/dd4hep::cm
<< " dz (cm) = " << (cryoDim.dz() - caloDim.dz())/dd4hep::cm << endmsg;
<< "ECAL cryostat: side: rmin (cm) = " << cryoDim.rmin2()/dd4hep::cm
<< " rmax (cm) = " << cryoDim.rmax1()/dd4hep::cm
<< " dz (cm) = " << (cryoDim.dz() - caloDim.dz())/dd4hep::cm << endmsg;
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()));
Expand Down Expand Up @@ -180,13 +180,13 @@ static dd4hep::detail::Ref_t createECalBarrelInclined(dd4hep::Detector& aLcdd,
dd4hep::Tube servicesFrontShape(cryoDim.rmin2(), bathRmin, caloDim.dz());
dd4hep::Tube servicesBackShape(bathRmax, cryoDim.rmax1(), caloDim.dz());
lLog << MSG::INFO
<< "ECAL services: front: rmin (cm) = " << cryoDim.rmin2()/dd4hep::cm
<< " rmax (cm) = " << bathRmin/dd4hep::cm/dd4hep::cm
<< " dz (cm) = " << caloDim.dz()/dd4hep::cm << endmsg;
<< "ECAL services: front: rmin (cm) = " << cryoDim.rmin2()/dd4hep::cm
<< " rmax (cm) = " << bathRmin/dd4hep::cm/dd4hep::cm
<< " dz (cm) = " << caloDim.dz()/dd4hep::cm << endmsg;
lLog << MSG::INFO
<< "ECAL services: back: rmin (cm) = " << bathRmax/dd4hep::cm
<< " rmax (cm) = " << cryoDim.rmax1()/dd4hep::cm
<< " dz (cm) = " << caloDim.dz()/dd4hep::cm << endmsg;
<< "ECAL services: back: rmin (cm) = " << bathRmax/dd4hep::cm
<< " rmax (cm) = " << cryoDim.rmax1()/dd4hep::cm
<< " dz (cm) = " << caloDim.dz()/dd4hep::cm << endmsg;
dd4hep::Volume servicesFrontVol("services_front", servicesFrontShape, aLcdd.material(activeMaterial));
dd4hep::Volume servicesBackVol("services_back", servicesBackShape, aLcdd.material(activeMaterial));
dd4hep::PlacedVolume servicesFrontPhysVol = envelopeVol.placeVolume(servicesFrontVol);
Expand Down Expand Up @@ -376,7 +376,7 @@ static dd4hep::detail::Ref_t createECalBarrelInclined(dd4hep::Detector& aLcdd,
//dd4hep::Box layerPassiveInnerShape(passiveInnerThickness / 2., caloDim.dz(), layerHeight[iLayer] / 2.);
dd4hep::Trd1 layerPassiveInnerShape(passiveInnerThicknessLayer[iLayer] / 2., passiveInnerThicknessLayer[iLayer+1] / 2., caloDim.dz(), layerHeight[iLayer] / 2.);
dd4hep::Volume layerPassiveInnerVol(passiveInnerMaterial, layerPassiveInnerShape,
aLcdd.material(passiveInnerMaterial));
aLcdd.material(passiveInnerMaterial));
layerPassiveInnerVol.setSensitiveDetector(aSensDet);
dd4hep::PlacedVolume layerPassiveInnerPhysVol =
passiveInnerVol.placeVolume(layerPassiveInnerVol, dd4hep::Position(0, 0, layerOffset));
Expand Down Expand Up @@ -699,26 +699,29 @@ static dd4hep::detail::Ref_t createECalBarrelInclined(dd4hep::Detector& aLcdd,
dd4hep::xml::setDetectorTypeFlag(xmlDetElem, caloDetElem);
dd4hep::rec::MaterialManager matMgr(envelopeVol);
dd4hep::rec::LayeredCalorimeterData::Layer caloLayer;

double rad_first = Rmin;
double rad_last = 0;
double scale_fact = dR / (-Rmin * cos(angle) + sqrt(pow(Rmax, 2) - pow(Rmin * sin(angle), 2)));
// since the layer height is given along the electrode and not along the radius it needs to be scaled to get the values of layer height radially
lLog << MSG::DEBUG << "Scaling factor " << scale_fact << endmsg;
for (size_t il = 0; il < layerHeight.size(); il++) {

runningHeight = 0.0;
for (uint iLay = 0; iLay < numLayers; iLay++) {
double Lmin = runningHeight;
double Lmax = runningHeight + layerHeight[iLay];
double Lmid = runningHeight + layerHeight[iLay]/2.0;
runningHeight += layerHeight[iLay];
double rad_first = sqrt(Rmin*Rmin + Lmin*Lmin + 2*Rmin*Lmin*cos(angle));
double rad_last = sqrt(Rmin*Rmin + Lmax*Lmax + 2*Rmin*Lmax*cos(angle));
double rad_mid = sqrt(Rmin*Rmin + Lmid*Lmid + 2*Rmin*Lmid*cos(angle));

double thickness_sen = 0.;
double absorberThickness = 0.;

rad_last = rad_first + (layerHeight[il] * scale_fact);
dd4hep::rec::Vector3D ivr1 = dd4hep::rec::Vector3D(0., rad_first, 0); // defining starting vector points of the given layer
dd4hep::rec::Vector3D ivr2 = dd4hep::rec::Vector3D(0., rad_last, 0); // defining end vector points of the given layer

lLog << MSG::DEBUG << "radius first " << rad_first << " radius last " << rad_last << endmsg;
lLog << MSG::DEBUG << "radius first " << rad_first << " radius last " << rad_last << " radius middle " << rad_mid << endmsg;
const dd4hep::rec::MaterialVec &materials = matMgr.materialsBetween(ivr1, ivr2); // calling material manager to get material info between two points
auto mat = matMgr.createAveragedMaterial(materials); // creating average of all the material between two points to calculate X0 and lambda of averaged material
const double nRadiationLengths = mat.radiationLength();
const double nInteractionLengths = mat.interactionLength();
const double difference_bet_r1r2 = (ivr1 - ivr2).r(); // equal to layerHeight[il]*scale_fact
const double difference_bet_r1r2 = (ivr1 - ivr2).r();
const double value_of_x0 = difference_bet_r1r2 / nRadiationLengths;
const double value_of_lambda = difference_bet_r1r2 / nInteractionLengths;
std::string str1("LAr");
Expand All @@ -733,20 +736,19 @@ static dd4hep::detail::Ref_t createECalBarrelInclined(dd4hep::Detector& aLcdd,
absorberThickness += materials.at(imat).second;
}
}
rad_first = rad_last;
lLog << MSG::DEBUG << "The sensitive thickness is " << thickness_sen << endmsg;
lLog << MSG::DEBUG << "The absorber thickness is " << absorberThickness << endmsg;
lLog << MSG::DEBUG << "The radiation length is " << value_of_x0 << " and the interaction length is " << value_of_lambda << endmsg;

caloLayer.distance = rad_first;
caloLayer.sensitive_thickness = thickness_sen;
caloLayer.inner_nRadiationLengths = value_of_x0 / 2.0;
caloLayer.inner_nInteractionLengths = value_of_lambda / 2.0;
caloLayer.inner_thickness = difference_bet_r1r2 / 2.0;
caloLayer.inner_nRadiationLengths = value_of_x0 * (rad_mid - rad_first)/(rad_last - rad_first);
caloLayer.inner_nInteractionLengths = value_of_lambda * (rad_mid - rad_first)/(rad_last - rad_first);
caloLayer.inner_thickness = rad_mid - rad_first;

caloLayer.outer_nRadiationLengths = value_of_x0 / 2.0;
caloLayer.outer_nInteractionLengths = value_of_lambda / 2.0;
caloLayer.outer_thickness = difference_bet_r1r2 / 2;
caloLayer.outer_nRadiationLengths = value_of_x0 * (rad_last - rad_mid)/(rad_last - rad_first);
caloLayer.outer_nInteractionLengths = value_of_lambda * (rad_last - rad_mid)/(rad_last - rad_first);
caloLayer.outer_thickness = rad_last - rad_mid;

caloLayer.absorberThickness = absorberThickness;
caloLayer.cellSize0 = 20 * dd4hep::mm;
Expand Down
Loading

0 comments on commit 6a8703e

Please sign in to comment.