diff --git a/include/DDCaloHitCreator.h b/include/DDCaloHitCreator.h index 7acd489..b4c6e8b 100644 --- a/include/DDCaloHitCreator.h +++ b/include/DDCaloHitCreator.h @@ -87,6 +87,12 @@ class DDCaloHitCreator float m_eCalScToHadGeVBarrel; ///< The calibration from deposited Sc-layer energy on the endcaps to hadronic energy float m_eCalSiToHadGeVEndCap; ///< The calibration from deposited Si-layer energy on the enecaps to hadronic energy float m_eCalScToHadGeVEndCap; ///< The calibration from deposited Sc-layer energy on the endcaps to hadronic energy + + + /// LEOP == HIT ENERGY OFFSET SUBTRACTION BY ECAL REGION (SAME ORG AS CALODIGI_BIB) + /// It works for a 6-layers CRILIN ECal, the geometry is hardcoded + FloatVector m_crilinBarrelOffsets; + FloatVector m_crilinEndcapOffsets; ///ADDED BY NIKIFOROS @@ -252,6 +258,10 @@ class DDCaloHitCreator dd4hep::VolumeManager m_volumeManager; ///< DD4hep volume manager + ///LEOP == function to subtract energy offset inCRIILIN ECal + float GetOffsetCrilinBarrel(float &x, float &y, float &z); + float GetOffsetCrilinEndcap(float &x, float &y, float &z); + float SubtractedEnergy(float x, float y, float z, float en, bool useEndcap); }; //------------------------------------------------------------------------------------------------------------------------------------------ diff --git a/src/DDCaloHitCreator.cc b/src/DDCaloHitCreator.cc index 61b0500..cfa4569 100644 --- a/src/DDCaloHitCreator.cc +++ b/src/DDCaloHitCreator.cc @@ -88,6 +88,8 @@ pandora::StatusCode DDCaloHitCreator::CreateCaloHits(const EVENT::LCEvent *const pandora::StatusCode DDCaloHitCreator::CreateECalCaloHits(const EVENT::LCEvent *const pLCEvent) { + float subtrEne = 0.; + for (StringVector::const_iterator iter = m_settings.m_eCalCaloHitCollections.begin(), iterEnd = m_settings.m_eCalCaloHitCollections.end(); iter != iterEnd; ++iter) { @@ -176,21 +178,28 @@ pandora::StatusCode DDCaloHitCreator::CreateECalCaloHits(const EVENT::LCEvent *c << std::endl; this->GetBarrelCaloHitProperties(pCaloHit, barrelLayers, m_settings.m_eCalBarrelInnerSymmetry, caloHitParameters, m_settings.m_eCalBarrelNormalVector, absorberCorrection); - - caloHitParameters.m_hadronicEnergy = eCalToHadGeVBarrel * pCaloHit->getEnergy(); + + subtrEne = SubtractedEnergy(pCaloHit->getPosition()[0], pCaloHit->getPosition()[1], pCaloHit->getPosition()[2], pCaloHit->getEnergy(), false); + if(subtrEne < 0.) continue; + + caloHitParameters.m_hadronicEnergy = eCalToHadGeVBarrel * subtrEne; } else { this->GetEndCapCaloHitProperties(pCaloHit, endcapLayers, caloHitParameters, absorberCorrection); - caloHitParameters.m_hadronicEnergy = eCalToHadGeVEndCap * pCaloHit->getEnergy(); + + subtrEne = SubtractedEnergy(pCaloHit->getPosition()[0], pCaloHit->getPosition()[1], pCaloHit->getPosition()[2], pCaloHit->getEnergy(), true); + if(subtrEne < 0.) continue; + + caloHitParameters.m_hadronicEnergy = eCalToHadGeVEndCap * subtrEne; } - caloHitParameters.m_mipEquivalentEnergy = pCaloHit->getEnergy() * eCalToMip * absorberCorrection; + caloHitParameters.m_mipEquivalentEnergy = subtrEne * eCalToMip * absorberCorrection; if (caloHitParameters.m_mipEquivalentEnergy.Get() < eCalMipThreshold) continue; - caloHitParameters.m_electromagneticEnergy = eCalToEMGeV * pCaloHit->getEnergy(); + caloHitParameters.m_electromagneticEnergy = eCalToEMGeV * subtrEne; // ATTN If using strip splitting, must correct cell sizes for use in PFA to minimum of strip width and strip length if (m_settings.m_stripSplittingOn) @@ -809,6 +818,68 @@ float DDCaloHitCreator::GetMaximumRadius(const EVENT::CalorimeterHit *const pCal return maximumRadius; } +//------------------------------------------------------------------------------------------------------------------------------------------ +///LEOP == functions to subtract offset from ecal hits, only for CRILIN ecal +float DDCaloHitCreator::GetOffsetCrilinBarrel(float &x, float &y, float &z){ + + float r = sqrt(x*x+y*y); + + float phi_rel= acos(x/r); + if (y<0) phi_rel = -phi_rel+2*3.14159; + phi_rel = phi_rel - round( phi_rel / ( 2 * 3.14159 / 12 ) ) * 2 * 3.14159 / 12; + float x_rel = r*cos(phi_rel); + + float nZones = ((float)m_settings.m_crilinBarrelOffsets.size()) / 6.; + + int nx = 0; + + if (x_rel<1720.) nx=0; + else if (x_rel>1720. && x_rel<1770.) nx=1; + else if (x_rel>1770. && x_rel<1810.) nx=2; + else if (x_rel>1810. && x_rel<1860.) nx=3; + else if (x_rel>1860. && x_rel<1900.) nx=4; + else if (x_rel>1900.) nx=5; + + int nz = (int)abs( z / ( 2210. / nZones ) ); + if (nz > nZones - 1) nz = nZones - 1; + + return m_settings.m_crilinBarrelOffsets[ nx*(int)nZones+nz ] * 0.001; +} + +float DDCaloHitCreator::GetOffsetCrilinEndcap(float &x, float &y, float &z){ + + float r_rel = sqrt(x*x+y*y) - 305.; //relative to endcap inner radius (-5mm tolerance) + + float nZones = ((float)m_settings.m_crilinEndcapOffsets.size()) / 6.; + + int nx = 0; + if (abs(z)<2350.) nx=0; + else if (abs(z)>2350. && abs(z)<2400.) nx=1; + else if (abs(z)>2400. && abs(z)<2450.) nx=2; + else if (abs(z)>2450. && abs(z)<2490.) nx=3; + else if (abs(z)>2490. && abs(z)<2530.) nx=4; + else if (abs(z)>2530.) nx=5; + + int nz = (int) r_rel / ( 1395. / nZones ); + if( nz < 0 ) nz = 0; + else if( nz > nZones - 1 ) nz = nZones - 1; + + return m_settings.m_crilinEndcapOffsets[ nx*(int)nZones+nz ] * 0.001; +} + +float DDCaloHitCreator::SubtractedEnergy(float x, float y, float z, float en, bool useEndcap){ + + float offset = 0.; + + if(useEndcap) offset = GetOffsetCrilinEndcap(x, y, z); + else offset = GetOffsetCrilinBarrel(x, y, z); + + float energy = en - offset; + energy = energy>0. ? energy : -1.; + + return energy; +} + //------------------------------------------------------------------------------------------------------------------------------------------ //------------------------------------------------------------------------------------------------------------------------------------------ @@ -865,6 +936,10 @@ DDCaloHitCreator::Settings::Settings() m_hCalBarrelOuterSymmetry(0.f), m_eCalBarrelNormalVector({0.0, 0.0, 1.0}), m_hCalBarrelNormalVector({0.0, 0.0, 1.0}), - m_muonBarrelNormalVector({0.0, 0.0, 1.0}) + m_muonBarrelNormalVector({0.0, 0.0, 1.0}), + + ///LEOP == offset subtraction parameters for ECal + m_crilinBarrelOffsets({0.,0.,0.,0.,0.,0.}), + m_crilinEndcapOffsets({0.,0.,0.,0.,0.,0.}) { } diff --git a/src/DDPandoraPFANewProcessor.cc b/src/DDPandoraPFANewProcessor.cc index 4b261a0..4aa0a82 100644 --- a/src/DDPandoraPFANewProcessor.cc +++ b/src/DDPandoraPFANewProcessor.cc @@ -854,6 +854,17 @@ void DDPandoraPFANewProcessor::ProcessSteeringFile() "The minimum correction to on ecal hit in Pandora energy correction", m_settings.m_minCleanCorrectedHitEnergy, softwareCompensationParameters.m_minCleanCorrectedHitEnergy); + + ///LEOP == offset subtraction in ECal + registerProcessorParameter("CrilinBarrelOffsets", + "Energy offsets to subtract in CRILIN ECal barrel in MeV", + m_caloHitCreatorSettings.m_crilinBarrelOffsets, + std::vector({0.,0.,0.,0.,0.,0.})); + + registerProcessorParameter("CrilinEndcapOffsets", + "Energy offsets to subtract in CRILIN ECal endcap in MeV", + m_caloHitCreatorSettings.m_crilinEndcapOffsets, + std::vector({0.,0.,0.,0.,0.,0.})); } //------------------------------------------------------------------------------------------------------------------------------------------