Skip to content
GitLab
About GitLab
GitLab: the DevOps platform
Explore GitLab
Install GitLab
How GitLab compares
Get started
GitLab docs
GitLab Learn
Pricing
Talk to an expert
/
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Projects
Groups
Snippets
Sign up now
Login
Sign in / Register
Toggle navigation
Menu
Open sidebar
EIC
Project Juggler
Commits
37c4cff4
Commit
37c4cff4
authored
Aug 03, 2021
by
Sylvester Joosten
Browse files
Fix issues introduced by switching off exposePODMembers in data model
parent
5f8cd5d6
Changes
17
Hide whitespace changes
Inline
Side-by-side
JugReco/src/components/CalorimeterHitsEtaPhiProjector.cpp
View file @
37c4cff4
...
...
@@ -97,7 +97,7 @@ namespace Jug::Reco {
std
::
unordered_map
<
std
::
pair
<
int64_t
,
int64_t
>
,
std
::
vector
<
eic
::
ConstCalorimeterHit
>
,
pair_hash
>
merged_hits
;
for
(
const
auto
h
:
*
m_inputHitCollection
.
get
())
{
Point3D
p
(
h
.
x
(),
h
.
y
(),
h
.
z
()
);
Point3D
p
(
h
.
position
().
x
,
h
.
position
().
y
,
h
.
position
().
z
);
auto
bins
=
std
::
make_pair
(
static_cast
<
int64_t
>
(
pos2bin
(
p
.
eta
(),
gridSizes
[
0
],
0.
)),
static_cast
<
int64_t
>
(
pos2bin
(
p
.
phi
(),
gridSizes
[
1
],
0.
)));
merged_hits
[
bins
].
push_back
(
h
);
...
...
@@ -111,7 +111,7 @@ namespace Jug::Reco {
hit
.
layerID
(
ref
.
layerID
());
// TODO, we can do timing cut to reject noises
hit
.
time
(
ref
.
time
());
double
r
=
std
::
hypot
(
ref
.
x
(),
ref
.
y
(),
ref
.
z
()
);
double
r
=
std
::
hypot
(
ref
.
position
().
x
,
ref
.
position
().
y
,
ref
.
position
().
z
);
double
eta
=
bin2pos
(
bins
.
first
,
gridSizes
[
0
],
0.
);
double
phi
=
bin2pos
(
bins
.
second
,
gridSizes
[
1
],
1.
);
double
theta
=
std
::
atan
(
std
::
exp
(
-
eta
))
*
2.
;
...
...
JugReco/src/components/CalorimeterIslandCluster.cpp
View file @
37c4cff4
...
...
@@ -44,24 +44,24 @@ typedef ROOT::Math::XYZPoint Point3D;
namespace
Jug
::
Reco
{
// helper functions to get distance between hits
static
Point
localDistXY
(
eic
::
ConstCalorimeterHit
h1
,
eic
::
ConstCalorimeterHit
h2
)
{
return
Point
(
h1
.
local
_x
()
-
h2
.
local
_x
(),
h1
.
local_y
()
-
h2
.
local
_y
());
return
Point
(
h1
.
local
()
.
local_x
-
h2
.
local
().
local_x
,
h1
.
local
().
local_y
-
h2
.
local
()
.
local_y
);
}
static
Point
localDistXZ
(
eic
::
ConstCalorimeterHit
h1
,
eic
::
ConstCalorimeterHit
h2
)
{
return
Point
(
h1
.
local
_x
()
-
h2
.
local
_x
(),
h1
.
local_z
()
-
h2
.
local
_z
());
return
Point
(
h1
.
local
()
.
local_x
-
h2
.
local
().
local_x
,
h1
.
local
().
local_z
-
h2
.
local
()
.
local_z
);
}
static
Point
localDistYZ
(
eic
::
ConstCalorimeterHit
h1
,
eic
::
ConstCalorimeterHit
h2
)
{
return
Point
(
h1
.
local
_y
()
-
h2
.
local
_y
(),
h1
.
local_z
()
-
h2
.
local
_z
());
return
Point
(
h1
.
local
()
.
local_y
-
h2
.
local
().
local_y
,
h1
.
local
().
local_z
-
h2
.
local
()
.
local_z
);
}
static
Point
dimScaledLocalDistXY
(
eic
::
ConstCalorimeterHit
h1
,
eic
::
ConstCalorimeterHit
h2
)
{
return
Point
(
2.
*
(
h1
.
local
_x
()
-
h2
.
local
_x
())
/
(
h1
.
dim_x
()
+
h2
.
dim_x
()
),
2.
*
(
h1
.
local
_y
()
-
h2
.
local
_y
())
/
(
h1
.
dim_y
()
+
h2
.
dim_y
()
));
return
Point
(
2.
*
(
h1
.
local
()
.
local_x
-
h2
.
local
().
local_x
)
/
(
h1
.
dimension
().
dim_x
+
h2
.
dimension
().
dim_x
),
2.
*
(
h1
.
local
()
.
local_y
-
h2
.
local
().
local_y
)
/
(
h1
.
dimension
().
dim_y
+
h2
.
dimension
().
dim_y
));
}
static
Point
globalDistRPhi
(
eic
::
ConstCalorimeterHit
h1
,
eic
::
ConstCalorimeterHit
h2
)
{
Point3D
p1
(
h1
.
x
(),
h1
.
y
(),
h1
.
z
()),
p2
(
h2
.
x
(),
h2
.
y
(),
h2
.
z
()
);
Point3D
p1
(
h1
.
position
().
x
,
h1
.
position
().
y
,
h1
.
position
().
z
),
p2
(
h2
.
position
().
x
,
h2
.
position
().
y
,
h2
.
position
().
z
);
return
Point
(
p1
.
r
()
-
p2
.
r
(),
p1
.
phi
()
-
p2
.
phi
());
}
static
Point
globalDistEtaPhi
(
eic
::
ConstCalorimeterHit
h1
,
eic
::
ConstCalorimeterHit
h2
)
{
Point3D
p1
(
h1
.
x
(),
h1
.
y
(),
h1
.
z
()),
p2
(
h2
.
x
(),
h2
.
y
(),
h2
.
z
()
);
Point3D
p1
(
h1
.
position
().
x
,
h1
.
position
().
y
,
h1
.
position
().
z
),
p2
(
h2
.
position
().
x
,
h2
.
position
().
y
,
h2
.
position
().
z
);
return
Point
(
p1
.
eta
()
-
p2
.
eta
(),
p1
.
phi
()
-
p2
.
phi
());
}
// name: {method, units}
...
...
@@ -187,8 +187,8 @@ namespace Jug::Reco {
for
(
size_t
i
=
0
;
i
<
hits
.
size
();
++
i
)
{
debug
()
<<
fmt
::
format
(
"hit {:d}: energy = {:.4f} MeV, local = ({:.4f}, {:.4f}) mm, "
"global=({:.4f}, {:.4f}, {:.4f}) mm, layer = {:d}, sector = {:d}."
,
i
,
hits
[
i
].
energy
()
*
1000.
,
hits
[
i
].
local
_x
(),
hits
[
i
].
local
_y
(),
hits
[
i
].
x
(),
hits
[
i
].
y
(),
hits
[
i
].
z
()
,
i
,
hits
[
i
].
energy
()
*
1000.
,
hits
[
i
].
local
()
.
local_x
,
hits
[
i
].
local
()
.
local_y
,
hits
[
i
].
position
().
x
,
hits
[
i
].
position
().
y
,
hits
[
i
].
position
().
z
,
hits
[
i
].
layerID
(),
hits
[
i
].
sectorID
())
<<
endmsg
;
// already in a group
if
(
visits
[
i
])
{
...
...
@@ -230,7 +230,7 @@ namespace Jug::Reco {
// different sector, local coordinates do not work, using global coordinates
}
else
{
// sector may have rotation (barrel), so z is included
return
dist3d
(
h1
.
x
()
-
h2
.
x
(),
h1
.
y
()
-
h2
.
y
(),
h1
.
z
()
-
h2
.
z
()
)
<=
sectorDist
;
return
dist3d
(
h1
.
position
().
x
-
h2
.
position
().
x
,
h1
.
position
().
y
-
h2
.
position
().
y
,
h1
.
position
().
z
-
h2
.
position
().
z
)
<=
sectorDist
;
}
}
...
...
@@ -347,7 +347,7 @@ namespace Jug::Reco {
size_t
j
=
0
;
// calculate weights for local maxima
for
(
auto
cit
=
maxima
.
begin
();
cit
!=
maxima
.
end
();
++
cit
,
++
j
)
{
double
dist_ref
=
cit
->
dim_x
()
;
double
dist_ref
=
cit
->
dim
ension
().
dim
_x
;
double
energy
=
cit
->
energy
();
double
dist
=
hitsDist
(
*
cit
,
*
it
).
r
();
weights
[
j
]
=
std
::
exp
(
-
dist
/
dist_ref
)
*
energy
;
...
...
JugReco/src/components/ClusterRecoCoG.cpp
View file @
37c4cff4
...
...
@@ -194,13 +194,13 @@ namespace Jug::Reco {
// info() << std::log(hit.energy()/totalE) << endmsg;
float
w
=
weightFunc
(
hit
.
energy
(),
totalE
,
m_logWeightBase
.
value
(),
0
);
tw
+=
w
;
x
+=
hit
.
x
()
*
w
;
y
+=
hit
.
y
()
*
w
;
z
+=
hit
.
z
()
*
w
;
x
+=
hit
.
position
().
x
*
w
;
y
+=
hit
.
position
().
y
*
w
;
z
+=
hit
.
position
().
z
*
w
;
/*
debug() << hit.cellID() << ": (" << hit.local_x() << ", " << hit.local_y() << ", "
<< hit.local_z() << "), "
<< "(" << hit.
x()
<< ", " << hit.
y()
<< ", " << hit.
z()
<< "), " << endmsg;
<< "(" << hit.
position().x
<< ", " << hit.
position().y
<< ", " << hit.
position().z
<< "), " << endmsg;
*/
}
if
(
tw
==
0.
)
{
...
...
@@ -210,7 +210,7 @@ namespace Jug::Reco {
// convert global position to local position, use the cell with max edep as a reference
const
auto
volman
=
m_geoSvc
->
detector
()
->
volumeManager
();
const
auto
alignment
=
volman
.
lookupDetElement
(
centerID
).
nominal
();
const
auto
lpos
=
alignment
.
worldToLocal
(
dd4hep
::
Position
(
res
.
x
(),
res
.
y
(),
res
.
z
()
));
const
auto
lpos
=
alignment
.
worldToLocal
(
dd4hep
::
Position
(
res
.
position
().
x
,
res
.
position
().
y
,
res
.
position
().
z
));
// TODO: may need convert back to have depthCorrection in global positions
res
.
local
({
lpos
.
x
(),
lpos
.
y
(),
lpos
.
z
()
+
m_depthCorrection
});
...
...
JugReco/src/components/ImagingClusterReco.cpp
View file @
37c4cff4
...
...
@@ -158,21 +158,19 @@ namespace Jug::Reco {
double
mz
=
0.
;
double
edep
=
0.
;
for
(
const
auto
&
hit
:
hits
)
{
mx
+=
hit
.
x
()
;
my
+=
hit
.
y
()
;
mz
+=
hit
.
z
()
;
mx
+=
hit
.
position
().
x
;
my
+=
hit
.
position
().
y
;
mz
+=
hit
.
position
().
z
;
edep
+=
hit
.
edep
();
}
layer
.
x
(
mx
/
layer
.
nhits
());
layer
.
y
(
my
/
layer
.
nhits
());
layer
.
z
(
mz
/
layer
.
nhits
());
layer
.
position
({
mx
/
layer
.
nhits
(),
my
/
layer
.
nhits
(),
mz
/
layer
.
nhits
()});
layer
.
edep
(
edep
);
double
radius
=
0.
;
for
(
const
auto
&
hit
:
hits
)
{
radius
+=
std
::
sqrt
(
pow2
(
hit
.
x
()
-
layer
.
x
())
+
pow2
(
hit
.
y
()
-
layer
.
y
()
)
+
pow2
(
hit
.
z
()
-
layer
.
z
()
));
radius
+=
std
::
sqrt
(
pow2
(
hit
.
position
().
x
-
layer
.
position
().
x
)
+
pow2
(
hit
.
position
().
y
-
layer
.
position
().
y
)
+
pow2
(
hit
.
position
().
z
-
layer
.
position
().
z
));
}
layer
.
radius
(
radius
/
layer
.
nhits
());
return
layer
;
...
...
@@ -190,9 +188,9 @@ namespace Jug::Reco {
double
r
=
9999
*
cm
;
for
(
const
auto
&
hit
:
hits
)
{
meta
+=
hit
.
eta
()
*
hit
.
edep
();
mphi
+=
hit
.
phi
()
*
hit
.
edep
();
mphi
+=
hit
.
p
olar
().
p
hi
*
hit
.
edep
();
edep
+=
hit
.
edep
();
r
=
std
::
min
(
hit
.
r
(),
r
);
r
=
std
::
min
(
hit
.
pola
r
()
.
r
,
r
);
}
const
double
eta
=
meta
/
edep
;
const
double
phi
=
mphi
/
edep
;
...
...
@@ -201,20 +199,15 @@ namespace Jug::Reco {
cluster
.
edep
(
edep
);
cluster
.
energy
(
edep
/
m_sampFrac
);
// simple energy reconstruction
cluster
.
eta
(
eta
);
cluster
.
phi
(
phi
);
cluster
.
theta
(
theta
);
// project it to the inner surface (the shallowest hit)
cluster
.
r
(
r
);
cluster
.
polar
({
r
,
phi
,
theta
});
// cartesian coordinates
ROOT
::
Math
::
Polar3DVectorD
polar
{
r
,
theta
,
(
phi
>
M_PI
?
phi
-
M_PI
:
phi
)};
cluster
.
x
(
polar
.
x
());
cluster
.
y
(
polar
.
y
());
cluster
.
z
(
polar
.
z
());
cluster
.
position
({
polar
.
x
(),
polar
.
y
(),
polar
.
z
()});
// shower radius estimate (eta-phi plane)
double
radius
=
0.
;
for
(
auto
hit
:
cluster
.
hits
())
{
radius
+=
std
::
sqrt
(
pow2
(
hit
.
eta
()
-
cluster
.
eta
())
+
pow2
(
hit
.
phi
()
-
cluster
.
phi
()
));
radius
+=
std
::
sqrt
(
pow2
(
hit
.
eta
()
-
cluster
.
eta
())
+
pow2
(
hit
.
p
olar
().
p
hi
-
cluster
.
p
olar
().
p
hi
));
}
cluster
.
radius
(
radius
/
cluster
.
nhits
());
...
...
@@ -230,9 +223,9 @@ namespace Jug::Reco {
double
mz
=
0.
;
for
(
const
auto
&
layer
:
layers
)
{
if
((
layer
.
layerID
()
<=
stop_layer
)
&&
(
layer
.
nhits
()
>
0
))
{
mx
+=
layer
.
x
()
;
my
+=
layer
.
y
()
;
mz
+=
layer
.
z
()
;
mx
+=
layer
.
position
().
x
;
my
+=
layer
.
position
().
y
;
mz
+=
layer
.
position
().
z
;
nrows
+=
1
;
}
}
...
...
@@ -249,9 +242,9 @@ namespace Jug::Reco {
int
ir
=
0
;
for
(
const
auto
&
layer
:
layers
)
{
if
((
layer
.
layerID
()
<=
stop_layer
)
&&
(
layer
.
nhits
()
>
0
))
{
pos
(
ir
,
0
)
=
layer
.
x
()
-
mx
;
pos
(
ir
,
1
)
=
layer
.
y
()
-
my
;
pos
(
ir
,
2
)
=
layer
.
z
()
-
mz
;
pos
(
ir
,
0
)
=
layer
.
position
().
x
-
mx
;
pos
(
ir
,
1
)
=
layer
.
position
().
y
-
my
;
pos
(
ir
,
2
)
=
layer
.
position
().
z
-
mz
;
ir
+=
1
;
}
}
...
...
JugReco/src/components/ImagingPixelMerger.cpp
View file @
37c4cff4
...
...
@@ -94,10 +94,10 @@ namespace Jug::Reco {
if
((
int
)
k
>=
m_nLayers
)
{
continue
;
}
double
r
=
std
::
sqrt
(
h
.
x
()
*
h
.
x
()
+
h
.
y
()
*
h
.
y
()
+
h
.
z
()
*
h
.
z
()
);
double
th
=
std
::
acos
(
h
.
z
()
/
r
);
double
r
=
std
::
sqrt
(
h
.
position
().
x
*
h
.
position
().
x
+
h
.
position
().
y
*
h
.
position
().
y
+
h
.
position
().
z
*
h
.
position
().
z
);
double
th
=
std
::
acos
(
h
.
position
().
z
/
r
);
double
eta
=
-
std
::
log
(
std
::
tan
(
th
/
2.
));
double
phi
=
std
::
atan2
(
h
.
y
(),
h
.
x
()
);
double
phi
=
std
::
atan2
(
h
.
position
().
y
,
h
.
position
().
x
);
// debug() << th << ", " << eta << ", " << phi << endmsg;
auto
&
layer
=
group_hits
[
k
];
...
...
@@ -135,7 +135,7 @@ namespace Jug::Reco {
auto
h
=
mhits
.
create
();
h
.
edep
(
grid
.
energy
);
h
.
eta
(
grid
.
eta
);
h
.
p
hi
(
grid
.
phi
);
h
.
p
olar
({
0.
,
0.
,
grid
.
phi
}
);
h
.
layerID
(
i
);
h
.
hitID
(
k
);
}
...
...
JugReco/src/components/ImagingTopoCluster.cpp
View file @
37c4cff4
...
...
@@ -135,8 +135,8 @@ namespace Jug::Reco {
/*
debug() << fmt::format("hit {:d}: local position = ({}, {}, {}), global position = ({}, {}, {})",
i + 1,
hits[i].local
_x
(), hits[i].local
_y
(), hits[i].local_z(),
hits[i].
x(), hits[i].y(), hits[i].z()
)
hits[i].local()
.local_x
, hits[i].local()
.local_y
, hits[i].local_z(),
hits[i].
position().x, hits[i].position().y, hits[i].position().z
)
<< endmsg;
*/
// already in a group, or not energetic enough to form a cluster
...
...
@@ -189,7 +189,7 @@ namespace Jug::Reco {
{
// different sectors, simple distance check
if
(
h1
.
sectorID
()
!=
h2
.
sectorID
())
{
return
std
::
sqrt
(
pow2
(
h1
.
x
()
-
h2
.
x
())
+
pow2
(
h1
.
y
()
-
h2
.
y
())
+
pow2
(
h1
.
z
()
-
h2
.
z
()
))
<=
return
std
::
sqrt
(
pow2
(
h1
.
position
().
x
-
h2
.
position
().
x
)
+
pow2
(
h1
.
position
().
y
-
h2
.
position
().
y
)
+
pow2
(
h1
.
position
().
z
-
h2
.
position
().
z
))
<=
sectorDist
;
}
...
...
@@ -197,11 +197,11 @@ namespace Jug::Reco {
int
ldiff
=
std
::
abs
(
h1
.
layerID
()
-
h2
.
layerID
());
// same layer, check local positions
if
(
!
ldiff
)
{
return
(
std
::
abs
(
h1
.
local
_x
()
-
h2
.
local
_x
())
<=
localDistXY
[
0
])
&&
(
std
::
abs
(
h1
.
local
_y
()
-
h2
.
local
_y
())
<=
localDistXY
[
1
]);
return
(
std
::
abs
(
h1
.
local
()
.
local_x
-
h2
.
local
()
.
local_x
)
<=
localDistXY
[
0
])
&&
(
std
::
abs
(
h1
.
local
()
.
local_y
-
h2
.
local
()
.
local_y
)
<=
localDistXY
[
1
]);
}
else
if
(
ldiff
<=
m_neighbourLayersRange
)
{
return
(
std
::
abs
(
h1
.
eta
()
-
h2
.
eta
())
<=
layerDistEtaPhi
[
0
])
&&
(
std
::
abs
(
h1
.
p
hi
()
-
h2
.
phi
()
)
<=
layerDistEtaPhi
[
1
]);
(
std
::
abs
(
h1
.
p
olar
().
phi
-
h2
.
polar
().
phi
)
<=
layerDistEtaPhi
[
1
]);
}
// not in adjacent layers
...
...
JugReco/src/components/PhotoRingClusters.cpp
View file @
37c4cff4
...
...
@@ -9,47 +9,43 @@
#include
"Gaudi/Property.h"
#include
"GaudiAlg/GaudiAlgorithm.h"
#include
"GaudiKernel/ToolHandle.h"
#include
"GaudiAlg/Transformer.h"
#include
"GaudiAlg/GaudiTool.h"
#include
"Gaudi
Kernel/RndmGenerators
.h"
#include
"Gaudi
Alg/Transformer
.h"
#include
"GaudiKernel/PhysicalConstants.h"
#include
"GaudiKernel/RndmGenerators.h"
#include
"GaudiKernel/ToolHandle.h"
#include
"DDRec/CellIDPositionConverter.h"
#include
"DDRec/SurfaceManager.h"
#include
"DDRec/Surface.h"
#include
"DDRec/SurfaceManager.h"
// FCCSW
#include
"JugBase/DataHandle.h"
#include
"JugBase/IGeoSvc.h"
// Event Model related classes
#include
"FuzzyKClusters.h"
#include
"eicd/PMTHitCollection.h"
#include
"eicd/RIChClusterCollection.h"
#include
"FuzzyKClusters.h"
using
namespace
Gaudi
::
Units
;
using
namespace
Eigen
;
namespace
Jug
::
Reco
{
/** Clustering Algorithm for Ring Imaging Cherenkov (RICH) events.
*
* \ingroup reco
*/
class
PhotoRingClusters
:
public
GaudiAlgorithm
{
public:
DataHandle
<
eic
::
PMTHitCollection
>
m_inputHitCollection
{
"inputHitCollection"
,
Gaudi
::
DataHandle
::
Reader
,
this
};
DataHandle
<
eic
::
RIChClusterCollection
>
m_outputClusterCollection
{
"outputClusterCollection"
,
Gaudi
::
DataHandle
::
Writer
,
this
};
/** Clustering Algorithm for Ring Imaging Cherenkov (RICH) events.
*
* \ingroup reco
*/
class
PhotoRingClusters
:
public
GaudiAlgorithm
{
public:
DataHandle
<
eic
::
PMTHitCollection
>
m_inputHitCollection
{
"inputHitCollection"
,
Gaudi
::
DataHandle
::
Reader
,
this
};
DataHandle
<
eic
::
RIChClusterCollection
>
m_outputClusterCollection
{
"outputClusterCollection"
,
Gaudi
::
DataHandle
::
Writer
,
this
};
// @TODO
// A more realistic way is to have tracker info as the input to determine how much clusters should be found
Gaudi
::
Property
<
int
>
m_nRings
{
this
,
"nRings"
,
1
};
Gaudi
::
Property
<
int
>
m_nIters
{
this
,
"nIters"
,
1000
};
Gaudi
::
Property
<
int
>
m_nRings
{
this
,
"nRings"
,
1
};
Gaudi
::
Property
<
int
>
m_nIters
{
this
,
"nIters"
,
1000
};
Gaudi
::
Property
<
double
>
m_q
{
this
,
"q"
,
2.0
};
Gaudi
::
Property
<
double
>
m_eps
{
this
,
"epsilon"
,
1e-4
};
Gaudi
::
Property
<
double
>
m_minNpe
{
this
,
"minNpe"
,
0.5
};
...
...
@@ -57,62 +53,59 @@ public:
SmartIF
<
IGeoSvc
>
m_geoSvc
;
// ill-formed: using GaudiAlgorithm::GaudiAlgorithm;
PhotoRingClusters
(
const
std
::
string
&
name
,
ISvcLocator
*
svcLoc
)
:
GaudiAlgorithm
(
name
,
svcLoc
)
PhotoRingClusters
(
const
std
::
string
&
name
,
ISvcLocator
*
svcLoc
)
:
GaudiAlgorithm
(
name
,
svcLoc
)
{
declareProperty
(
"inputHitCollection"
,
m_inputHitCollection
,
""
);
declareProperty
(
"outputClusterCollection"
,
m_outputClusterCollection
,
""
);
declareProperty
(
"inputHitCollection"
,
m_inputHitCollection
,
""
);
declareProperty
(
"outputClusterCollection"
,
m_outputClusterCollection
,
""
);
}
StatusCode
initialize
()
override
{
if
(
GaudiAlgorithm
::
initialize
().
isFailure
())
{
return
StatusCode
::
FAILURE
;
}
m_geoSvc
=
service
(
"GeoSvc"
);
if
(
!
m_geoSvc
)
{
error
()
<<
"Unable to locate Geometry Service. "
<<
"Make sure you have GeoSvc and SimSvc in the right order in the configuration."
<<
endmsg
;
return
StatusCode
::
FAILURE
;
}
return
StatusCode
::
SUCCESS
;
if
(
GaudiAlgorithm
::
initialize
().
isFailure
())
{
return
StatusCode
::
FAILURE
;
}
m_geoSvc
=
service
(
"GeoSvc"
);
if
(
!
m_geoSvc
)
{
error
()
<<
"Unable to locate Geometry Service. "
<<
"Make sure you have GeoSvc and SimSvc in the right order in the configuration."
<<
endmsg
;
return
StatusCode
::
FAILURE
;
}
return
StatusCode
::
SUCCESS
;
}
StatusCode
execute
()
override
{
// input collections
const
auto
&
rawhits
=
*
m_inputHitCollection
.
get
();
// Create output collections
auto
&
clusters
=
*
m_outputClusterCollection
.
createAndPut
();
// algorithm
auto
alg
=
fkc
::
KRings
();
// fill data
MatrixXd
data
(
rawhits
.
size
(),
2
);
for
(
int
i
=
0
;
i
<
data
.
rows
();
++
i
)
{
if
(
rawhits
[
i
].
npe
()
>
m_minNpe
)
{
data
.
row
(
i
)
<<
rawhits
[
i
].
local_x
(),
rawhits
[
i
].
local_y
();
}
// input collections
const
auto
&
rawhits
=
*
m_inputHitCollection
.
get
();
// Create output collections
auto
&
clusters
=
*
m_outputClusterCollection
.
createAndPut
();
// algorithm
auto
alg
=
fkc
::
KRings
();
// fill data
MatrixXd
data
(
rawhits
.
size
(),
2
);
for
(
int
i
=
0
;
i
<
data
.
rows
();
++
i
)
{
if
(
rawhits
[
i
].
npe
()
>
m_minNpe
)
{
data
.
row
(
i
)
<<
rawhits
[
i
].
local
().
local_x
,
rawhits
[
i
].
local
().
local_y
;
}
}
// clustering
auto
res
=
alg
.
Fit
(
data
,
m_nRings
,
m_q
,
m_eps
,
m_nIters
);
// clustering
auto
res
=
alg
.
Fit
(
data
,
m_nRings
,
m_q
,
m_eps
,
m_nIters
);
// local position
for
(
int
i
=
0
;
i
<
res
.
rows
();
++
i
)
{
auto
cl
=
clusters
.
create
();
cl
.
x
(
res
(
i
,
0
));
cl
.
y
(
res
(
i
,
1
));
cl
.
radius
(
res
(
i
,
2
));
}
// local position
for
(
int
i
=
0
;
i
<
res
.
rows
();
++
i
)
{
auto
cl
=
clusters
.
create
();
cl
.
position
({
res
(
i
,
0
),
res
(
i
,
1
)});
cl
.
radius
(
res
(
i
,
2
));
}
return
StatusCode
::
SUCCESS
;
return
StatusCode
::
SUCCESS
;
}
};
};
DECLARE_COMPONENT
(
PhotoRingClusters
)
DECLARE_COMPONENT
(
PhotoRingClusters
)
}
// namespace Jug::Reco
JugReco/src/components/SimpleClustering.cpp
View file @
37c4cff4
...
...
@@ -34,7 +34,7 @@ namespace Jug::Reco {
using
RecHits
=
eic
::
CalorimeterHitCollection
;
using
Clusters
=
eic
::
ClusterCollection
;
DataHandle
<
RecHits
>
m_inputHitCollection
{
"inputHitCollection"
,
Gaudi
::
DataHandle
::
Reader
,
this
};
DataHandle
<
RecHits
>
m_inputHitCollection
{
"inputHitCollection"
,
Gaudi
::
DataHandle
::
Reader
,
this
};
DataHandle
<
Clusters
>
m_outputClusters
{
"outputClusters"
,
Gaudi
::
DataHandle
::
Writer
,
this
};
Gaudi
::
Property
<
double
>
m_minModuleEdep
{
this
,
"minModuleEdep"
,
5.0
*
MeV
};
Gaudi
::
Property
<
double
>
m_maxDistance
{
this
,
"maxDistance"
,
20.0
*
cm
};
...
...
@@ -56,8 +56,7 @@ namespace Jug::Reco {
m_geoSvc
=
service
(
"GeoSvc"
);
if
(
!
m_geoSvc
)
{
error
()
<<
"Unable to locate Geometry Service. "
<<
"Make sure you have GeoSvc and SimSvc in the right order in the configuration."
<<
endmsg
;
<<
"Make sure you have GeoSvc and SimSvc in the right order in the configuration."
<<
endmsg
;
return
StatusCode
::
FAILURE
;
}
return
StatusCode
::
SUCCESS
;
...
...
@@ -79,9 +78,8 @@ namespace Jug::Reco {
eic
::
CalorimeterHit
ref_hit
;
bool
have_ref
=
false
;
for
(
const
auto
&
ch
:
hits
)
{
const
eic
::
CalorimeterHit
h
{
ch
.
cellID
(),
ch
.
clusterID
(),
ch
.
layerID
(),
ch
.
sectorID
(),
ch
.
energy
(),
ch
.
time
(),
ch
.
position
(),
ch
.
local
(),
ch
.
dimension
(),
1
};
const
eic
::
CalorimeterHit
h
{
ch
.
cellID
(),
ch
.
clusterID
(),
ch
.
layerID
(),
ch
.
sectorID
(),
ch
.
energy
(),
ch
.
time
(),
ch
.
position
(),
ch
.
local
(),
ch
.
dimension
(),
1
};
if
(
!
have_ref
||
h
.
energy
()
>
ref_hit
.
energy
())
{
ref_hit
=
h
;
have_ref
=
true
;
...
...
@@ -96,8 +94,8 @@ namespace Jug::Reco {
std
::
vector
<
eic
::
CalorimeterHit
>
cluster_hits
;
for
(
const
auto
&
h
:
the_hits
)
{
if
(
std
::
hypot
(
h
.
x
()
-
ref_hit
.
x
(),
h
.
y
()
-
ref_hit
.
y
(),
h
.
z
()
-
ref_hit
.
z
())
<
max_dist
)
{
if
(
std
::
hypot
(
h
.
position
().
x
-
ref_hit
.
position
().
x
,
h
.
position
().
y
-
ref_hit
.
position
().
y
,
h
.
position
().
z
-
ref_hit
.
position
().
z
)
<
max_dist
)
{
cluster_hits
.
push_back
(
h
);
}
else
{
remaining_hits
.
push_back
(
h
);
...
...
@@ -105,17 +103,17 @@ namespace Jug::Reco {
}
debug
()
<<
" cluster size "
<<
cluster_hits
.
size
()
<<
endmsg
;
double
total_energy
=
std
::
accumulate
(
std
::
begin
(
cluster_hits
),
std
::
end
(
cluster_hits
),
0.0
,
[](
double
t
,
const
eic
::
CalorimeterHit
&
h1
)
{
return
(
t
+
h1
.
energy
());
});
double
total_energy
=
std
::
accumulate
(
std
::
begin
(
cluster_hits
),
std
::
end
(
cluster_hits
),
0.0
,
[](
double
t
,
const
eic
::
CalorimeterHit
&
h1
)
{
return
(
t
+
h1
.
energy
());
});
debug
()
<<
" total_energy = "
<<
total_energy
<<
endmsg
;
eic
::
Cluster
cluster0
;
for
(
const
auto
&
h
:
cluster_hits
)
{
cluster0
.
energy
(
cluster0
.
energy
()
+
h
.
energy
());
cluster0
.
x
(
cluster0
.
x
()
+
h
.
energy
()
*
h
.
x
()
/
total_energy
);
cluster0
.
y
(
cluster0
.
y
()
+
h
.
energy
()
*
h
.
y
()
/
total_energy
);
cluster0
.
z
(
cluster0
.
z
()
+
h
.
energy
()
*
h
.
z
()
/
total_energy
);
cluster0
.
position
({
cluster0
.
position
().
x
+
h
.
energy
()
*
h
.
position
().
x
/
total_energy
,
cluster0
.
position
().
y
+
h
.
energy
()
*
h
.
position
().
y
/
total_energy
,
cluster0
.
position
().
z
+
h
.
energy
()
*
h
.
position
().
z
/
total_energy
}
);
}
clusters
.
push_back
(
cluster0
);
...
...
JugReco/src/components/TopologicalCellCluster.cpp
View file @
37c4cff4
...
...
@@ -190,7 +190,7 @@ namespace Jug::Reco {
int
s1
=
id_dec
->
get
(
h1
.
cellID
(),
sector_idx
);
int
s2
=
id_dec
->
get
(
h2
.
cellID
(),
sector_idx
);
if
(
s1
!=
s2
)
{
return
std
::
sqrt
(
pow2
(
h1
.
x
()
-
h2
.
x
())
+
pow2
(
h1
.
y
()
-
h2
.
y
())
+
pow2
(
h1
.
z
()
-
h2
.
z
()
))
<=
return
std
::
sqrt
(
pow2
(
h1
.
position
().
x
-
h2
.
position
().
x
)
+
pow2
(
h1
.
position
().
y
-
h2
.
position
().
y
)
+
pow2
(
h1
.
position
().
z
-
h2
.
position
().
z
))
<=
m_adjSectorDist
;
}
...
...
@@ -201,12 +201,12 @@ namespace Jug::Reco {
int
ldiff
=
std
::
abs
(
l1
-
l2
);
// same layer, check local positions
if
(
!
ldiff
)
{
return
(
std
::
abs
(
h1
.
local
_x
()
-
h2
.
local
_x
())
<=
u_localRanges
[
0
])
&&
(
std
::
abs
(
h1
.
local
_y
()
-
h2
.
local
_y
())
<=
u_localRanges
[
1
]);
return
(
std
::
abs
(
h1
.
local
()
.
local_x
-
h2
.
local
()
.
local_x
)
<=
u_localRanges
[
0
])
&&
(
std
::
abs
(
h1
.
local
()
.
local_y
-
h2
.
local
()
.
local_y
)
<=
u_localRanges
[
1
]);
}
else
if
(
ldiff
<=
m_adjLayerDiff
)
{
double
eta1
,
phi1
,
r1
,
eta2
,
phi2
,
r2
;
calc_eta_phi_r
(
h1
.
x
(),
h1
.
y
(),
h1
.
z
()
,
eta1
,
phi1
,
r1
);
calc_eta_phi_r
(
h2
.
x
(),
h2
.
y
(),
h2
.
z
()
,
eta2
,
phi2
,
r2
);
calc_eta_phi_r
(
h1
.
position
().
x
,
h1
.
position
().
y
,
h1
.
position
().
z
,
eta1
,
phi1
,
r1
);
calc_eta_phi_r
(
h2
.
position
().
x
,
h2
.
position
().
y
,
h2
.
position
().
z
,
eta2
,
phi2
,
r2
);
return
(
std
::
abs
(
eta1
-
eta2
)
<=
u_adjLayerRanges
[
0
])
&&
(
std
::
abs
(
phi1
-
phi2
)
<=
u_adjLayerRanges
[
1
]);
...
...
JugTrack/src/components/GenFitTrackFitter.cpp
View file @
37c4cff4
#include
"GenFitTrackFitter.h"
// Gaudi
#include
"Gaudi/Property.h"
#include
"GaudiAlg/GaudiAlgorithm.h"
#include
"GaudiKernel/ToolHandle.h"
#include
"GaudiAlg/Transformer.h"
#include
"GaudiAlg/GaudiTool.h"
#include
"GaudiAlg/Transformer.h"
#include
"GaudiKernel/RndmGenerators.h"
#include
"Gaudi
/Property
.h"
#include
"Gaudi
Kernel/ToolHandle
.h"
#include
"DDRec/CellIDPositionConverter.h"
#include
"DDRec/SurfaceManager.h"
#include
"DDRec/Surface.h"
#include
"DDRec/SurfaceManager.h"
#include
"JugBase/DataHandle.h"
#include
"JugBase/IGeoSvc.h"
#include
"JugTrack/BField.h"
#include
"JugTrack/GeometryContainers.hpp"
#include
"JugTrack/IndexSourceLink.hpp"
#include
"JugTrack/Track.hpp"
#include
"JugTrack/BField.h"
#include
"JugTrack/Measurement.hpp"
#include
"JugTrack/Track.hpp"
#include
"eicd/TrackerHitCollection.h"