Merge pull request #6 from LeonardoPalombini/music10Tev

LP: added hit offset subtraction for CRILIN ECal
MuonColliderSoft · Nov 7, 2024 · 961b89a · 961b89a
2 parents 028769a + 63b7842
commit 961b89a
Show file tree

Hide file tree

Showing 3 changed files with 102 additions and 6 deletions.
diff --git 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
@@ -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
@@ -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<float>({0.,0.,0.,0.,0.,0.}));
+
+    registerProcessorParameter("CrilinEndcapOffsets",
+                               "Energy offsets to subtract in CRILIN ECal endcap in MeV",
+                               m_caloHitCreatorSettings.m_crilinEndcapOffsets,
+                               std::vector<float>({0.,0.,0.,0.,0.,0.}));
 }
 
 //------------------------------------------------------------------------------------------------------------------------------------------