Skip to content

Commit

Permalink
CalorimeterClusterRecoCoG: set intrinsic{Phi,Theta} (#1511)
Browse files Browse the repository at this point in the history 
This partially reverts effect of #1391. The clusterShape has become a
way to extend EDM without going through extending it. We should sort
those variable out into proper members, this is a first step.

---------

Co-authored-by: Wouter Deconinck <[email protected]>
  • Loading branch information
@veprbl @wdconinc
veprbl and wdconinc authored Jul 26, 2024
1 parent 260f0d9 commit 80866b3
Show file tree
Hide file tree
Showing 7 changed files with 53 additions and 47 deletions.
49 changes: 25 additions & 24 deletions src/algorithms/calorimetry/CalorimeterClusterRecoCoG.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -228,12 +228,6 @@ std::optional<edm4eic::MutableCluster> CalorimeterClusterRecoCoG::reconstruct(co

// Additional convenience variables

// best estimate on the cluster direction is the cluster position
// for simple 2D CoG clustering
cl.setIntrinsicTheta(edm4hep::utils::anglePolar(cl.getPosition()));
cl.setIntrinsicPhi(edm4hep::utils::angleAzimuthal(cl.getPosition()));
// TODO errors

//_______________________________________
// Calculate cluster profile:
// radius,
Expand All @@ -248,12 +242,11 @@ std::optional<edm4eic::MutableCluster> CalorimeterClusterRecoCoG::reconstruct(co
Eigen::Vector3f sum1_3D = Eigen::Vector3f::Zero();
Eigen::Vector2cf eigenValues_2D = Eigen::Vector2cf::Zero();
Eigen::Vector3cf eigenValues_3D = Eigen::Vector3cf::Zero();
//the axis is the direction of the eigenvalue corresponding to the largest eigenvalue.
double axis_x=0, axis_y=0, axis_z=0;
if (cl.getNhits() > 1) {
// the axis is the direction of the eigenvalue corresponding to the largest eigenvalue.
edm4hep::Vector3f axis;

if (cl.getNhits() > 1) {
for (const auto& hit : pcl.getHits()) {

float w = weightFunc(hit.getEnergy(), totalE, logWeightBase, 0);

// theta, phi
Expand Down Expand Up @@ -306,19 +299,15 @@ std::optional<edm4eic::MutableCluster> CalorimeterClusterRecoCoG::reconstruct(co
return std::real(a) < std::real(b);
}
);
auto axis = eigenvectors.col(std::distance(
auto axis_eigen = eigenvectors.col(std::distance(
eigenValues_3D.begin(),
max_eigenvalue_it
));
axis_x=axis(0,0).real();
axis_y=axis(1,0).real();
axis_z=axis(2,0).real();
double norm=sqrt(axis_x*axis_x+axis_y*axis_y+axis_z*axis_z);
if (norm!=0){
axis_x/=norm;
axis_y/=norm;
axis_z/=norm;
}
axis = {
axis_eigen(0,0).real(),
axis_eigen(1,0).real(),
axis_eigen(2,0).real(),
};
}
}

Expand All @@ -329,10 +318,22 @@ std::optional<edm4eic::MutableCluster> CalorimeterClusterRecoCoG::reconstruct(co
cl.addToShapeParameters( eigenValues_3D[0].real() ); // 3D x-y-z cluster width 1
cl.addToShapeParameters( eigenValues_3D[1].real() ); // 3D x-y-z cluster width 2
cl.addToShapeParameters( eigenValues_3D[2].real() ); // 3D x-y-z cluster width 3
//last 3 shape parameters are the components of the axis direction
cl.addToShapeParameters( axis_x );
cl.addToShapeParameters( axis_y );
cl.addToShapeParameters( axis_z );


double dot_product = cl.getPosition() * axis;
if (dot_product < 0) {
axis = -1 * axis;
}

if (m_cfg.longitudinalShowerInfoAvailable) {
cl.setIntrinsicTheta(edm4hep::utils::anglePolar(axis));
cl.setIntrinsicPhi(edm4hep::utils::angleAzimuthal(axis));
// TODO intrinsicDirectionError
} else {
cl.setIntrinsicTheta(NAN);
cl.setIntrinsicPhi(NAN);
}

return std::move(cl);
}

Expand Down
2 changes: 2 additions & 0 deletions src/algorithms/calorimetry/CalorimeterClusterRecoCoGConfig.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,8 @@ namespace eicrecon {
std::vector<double> logWeightBaseCoeffs{};
double logWeightBase_Eref = 50 * dd4hep::MeV;

bool longitudinalShowerInfoAvailable = false;

// Constrain the cluster position eta to be within
// the eta of the contributing hits. This is useful to avoid edge effects
// for endcaps.
Expand Down
1 change: 1 addition & 0 deletions src/detectors/BEMC/BEMC.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -97,6 +97,7 @@ extern "C" {
.energyWeight = "log",
.sampFrac = 1.0,
.logWeightBase = 6.2,
.longitudinalShowerInfoAvailable = true,
.enableEtaBounds = false
},
app // TODO: Remove me once fixed
Expand Down
4 changes: 4 additions & 0 deletions src/detectors/FHCAL/FHCAL.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -117,6 +117,7 @@ extern "C" {
.energyWeight = "log",
.sampFrac = 0.0257,
.logWeightBase = 3.6,
.longitudinalShowerInfoAvailable = true,
.enableEtaBounds = true
},
app // TODO: Remove me once fixed
Expand All @@ -134,6 +135,7 @@ extern "C" {
.energyWeight = "log",
.sampFrac = 0.0257,
.logWeightBase = 6.2,
.longitudinalShowerInfoAvailable = true,
.enableEtaBounds = false,
},
app // TODO: Remove me once fixed
Expand Down Expand Up @@ -230,6 +232,7 @@ extern "C" {
.energyWeight = "log",
.sampFrac = 1.0,
.logWeightBase = 4.5,
.longitudinalShowerInfoAvailable = true,
.enableEtaBounds = false
},
app // TODO: Remove me once fixed
Expand All @@ -247,6 +250,7 @@ extern "C" {
.energyWeight = "log",
.sampFrac = 1.0,
.logWeightBase = 4.5,
.longitudinalShowerInfoAvailable = true,
.enableEtaBounds = false,
},
app // TODO: Remove me once fixed
Expand Down
5 changes: 5 additions & 0 deletions src/detectors/ZDC/ZDC.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -85,6 +85,7 @@ extern "C" {
.energyWeight = "log",
.sampFrac = 1.0,
.logWeightBase = 3.6,
.longitudinalShowerInfoAvailable = true,
.enableEtaBounds = false
},
app // TODO: Remove me once fixed
Expand All @@ -102,6 +103,7 @@ extern "C" {
.energyWeight = "log",
.sampFrac = 1.0,
.logWeightBase = 6.2,
.longitudinalShowerInfoAvailable = true,
.enableEtaBounds = false,
},
app // TODO: Remove me once fixed
Expand Down Expand Up @@ -193,6 +195,7 @@ extern "C" {
.sampFrac = 0.0203,
.logWeightBaseCoeffs={5.0,0.65,0.31},
.logWeightBase_Eref=50*dd4hep::GeV,
.longitudinalShowerInfoAvailable = true,
},
app // TODO: Remove me once fixed
));
Expand Down Expand Up @@ -228,6 +231,7 @@ extern "C" {
.energyWeight = "log",
.sampFrac = 1.0,
.logWeightBase = 3.6,
.longitudinalShowerInfoAvailable = true,
.enableEtaBounds = false
},
app // TODO: Remove me once fixed
Expand All @@ -244,6 +248,7 @@ extern "C" {
.energyWeight = "log",
.sampFrac = 0.0203,
.logWeightBase = 6.2,
.longitudinalShowerInfoAvailable = true,
.enableEtaBounds = false,
},
app // TODO: Remove me once fixed
Expand Down
1 change: 1 addition & 0 deletions src/factories/calorimetry/CalorimeterClusterRecoCoG_factory.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,7 @@ class CalorimeterClusterRecoCoG_factory : public JOmniFactory<CalorimeterCluster
ParameterRef<double> m_logWeightBase {this, "logWeightBase", config().logWeightBase};
ParameterRef<std::vector<double>> m_logWeightBaseCoeffs {this, "logWeightBaseCoeffs", config().logWeightBaseCoeffs};
ParameterRef<double> m_logWeightBase_Eref {this, "logWeightBase_Eref", config().logWeightBase_Eref};
ParameterRef<bool> m_longitudinalShowerInfoAvailable {this, "longitudinalShowerInfoAvailable", config().longitudinalShowerInfoAvailable};
ParameterRef<bool> m_enableEtaBounds {this, "enableEtaBounds", config().enableEtaBounds};

Service<AlgorithmsInit_service> m_algorithmsInit {this};
Expand Down
38 changes: 15 additions & 23 deletions src/tests/algorithms_test/calorimetry_CalorimeterClusterRecoCoG.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,14 +4,15 @@

#include <Evaluator/DD4hepUnits.h> // for MeV, mm, keV, ns
#include <catch2/catch_test_macros.hpp> // for AssertionHandler, operator""_catch_sr, StringRef, REQUIRE, operator<, operator==, operator>, TEST_CASE
#include <catch2/matchers/catch_matchers.hpp>
#include <catch2/matchers/catch_matchers_floating_point.hpp>
#include <edm4eic/CalorimeterHitCollection.h> // for CalorimeterHitCollection, MutableCalorimeterHit, CalorimeterHitMutableCollectionIterator
#include <edm4eic/ClusterCollection.h>
#include <edm4eic/MCRecoClusterParticleAssociationCollection.h>
#include <edm4eic/ProtoClusterCollection.h>
#include <edm4hep/SimCalorimeterHitCollection.h>
#include <edm4hep/Vector3f.h> // for Vector3f
#include <math.h>
#include <podio/RelationRange.h>
#include <spdlog/common.h> // for level_enum
#include <spdlog/logger.h> // for logger
#include <spdlog/spdlog.h> // for default_logger
Expand All @@ -28,16 +29,19 @@ using eicrecon::CalorimeterClusterRecoCoGConfig;
using edm4eic::CalorimeterHit;

TEST_CASE( "the calorimeter CoG algorithm runs", "[CalorimeterClusterRecoCoG]" ) {
const float EPSILON = 1e-5;

CalorimeterClusterRecoCoG algo("CalorimeterClusterRecoCoG");

std::shared_ptr<spdlog::logger> logger = spdlog::default_logger()->clone("CalorimeterClusterRecoCoG");
logger->set_level(spdlog::level::trace);

CalorimeterClusterRecoCoGConfig cfg;
cfg.energyWeight = "log";
cfg.sampFrac = 0.0203,
cfg.logWeightBaseCoeffs={5.0,0.65,0.31},
cfg.logWeightBase_Eref=50*dd4hep::GeV,
cfg.sampFrac = 0.0203;
cfg.logWeightBaseCoeffs = {5.0,0.65,0.31};
cfg.logWeightBase_Eref = 50*dd4hep::GeV;
cfg.longitudinalShowerInfoAvailable = true;


algo.applyConfig(cfg);
Expand All @@ -51,7 +55,7 @@ TEST_CASE( "the calorimeter CoG algorithm runs", "[CalorimeterClusterRecoCoG]" )

//create a protocluster with 3 hits
auto pclust = pclust_coll.create();
edm4hep::Vector3f position({0,0,0*dd4hep::mm});
edm4hep::Vector3f position({0, 0, 1 *dd4hep::mm});

auto hit1 = hits_coll.create();
hit1.setCellID(0);
Expand All @@ -63,8 +67,9 @@ TEST_CASE( "the calorimeter CoG algorithm runs", "[CalorimeterClusterRecoCoG]" )
hit1.setDimension({0,0,0});
hit1.setLocal(position);
pclust.addToHits(hit1);
pclust.addToWeights(1);

position={0,0, 1*dd4hep::mm};
position={-1 * dd4hep::mm, 0, 2 * dd4hep::mm};
auto hit2 = hits_coll.create();
hit2.setCellID(0);
hit2.setEnergy(0.1*dd4hep::GeV);
Expand All @@ -75,19 +80,7 @@ TEST_CASE( "the calorimeter CoG algorithm runs", "[CalorimeterClusterRecoCoG]" )
hit2.setDimension({0,0,0});
hit2.setLocal(position);
pclust.addToHits(hit2);

position={0,0, 2*dd4hep::mm};
auto hit3 = hits_coll.create();//0, 0.1*dd4hep::GeV, 0,0,0,position, {0,0,0}, 0,0, position);
hit3.setCellID(0);
hit3.setEnergy(0.1*dd4hep::GeV);
hit3.setEnergyError(0);
hit3.setTime(0);
hit3.setTimeError(0);
hit3.setPosition(position);
hit3.setDimension({0,0,0});
hit3.setLocal(position);
pclust.addToHits(hit3);
pclust.addToWeights(1);pclust.addToWeights(1);pclust.addToWeights(1);
pclust.addToWeights(1);

// Constructing input and output as per the algorithm's expected signature
auto input = std::make_tuple(&pclust_coll, &simhits);
Expand All @@ -97,10 +90,9 @@ TEST_CASE( "the calorimeter CoG algorithm runs", "[CalorimeterClusterRecoCoG]" )


for (auto clust : *clust_coll){
// require that this cluster's axis is 0,0,1
REQUIRE(clust.getShapeParameters()[7] == 0);
REQUIRE(clust.getShapeParameters()[8] == 0);
REQUIRE(abs(clust.getShapeParameters()[9]) == 1);
REQUIRE_THAT(clust.getIntrinsicTheta(), Catch::Matchers::WithinAbs(M_PI / 4, EPSILON));
// std::abs() checks if we land on -M_PI
REQUIRE_THAT(std::abs(clust.getIntrinsicPhi()), Catch::Matchers::WithinAbs(M_PI, EPSILON));
}


Expand Down

0 comments on commit 80866b3

Please sign in to comment.