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
8aebe6a3
Commit
8aebe6a3
authored
Jul 08, 2021
by
Chao Peng
Browse files
Ensure unit consistency among clustering algorithms
parent
0a8cdc20
Changes
7
Hide whitespace changes
Inline
Side-by-side
JugDigi/src/components/CalorimeterHitDigi.cpp
View file @
8aebe6a3
...
...
@@ -33,10 +33,6 @@ public:
Gaudi
::
Property
<
std
::
vector
<
double
>>
u_eRes
{
this
,
"energyResolutions"
,
{}};
// a/sqrt(E/GeV) + b + c/(E/GeV)
Gaudi
::
Property
<
double
>
m_tRes
{
this
,
"timineResolution"
,
0.0
*
ns
};
// input units, should be fixed
Gaudi
::
Property
<
double
>
m_eUnit
{
this
,
"inputEnergyUnit"
,
GeV
};
Gaudi
::
Property
<
double
>
m_tUnit
{
this
,
"inputTimeUnit"
,
ns
};
// digitization settings
Gaudi
::
Property
<
int
>
m_capADC
{
this
,
"capacityADC"
,
8096
};
Gaudi
::
Property
<
double
>
m_dyRangeADC
{
this
,
"dynamicRangeADC"
,
100
*
MeV
};
...
...
@@ -48,7 +44,8 @@ public:
m_inputHitCollection
{
"inputHitCollection"
,
Gaudi
::
DataHandle
::
Reader
,
this
};
DataHandle
<
eic
::
RawCalorimeterHitCollection
>
m_outputHitCollection
{
"outputHitCollection"
,
Gaudi
::
DataHandle
::
Writer
,
this
};
double
res
[
3
]
=
{
0.
,
0.
,
0.
};
// unitless counterparts of inputs
double
dyRangeADC
,
tRes
,
eRes
[
3
]
=
{
0.
,
0.
,
0.
};
// ill-formed: using GaudiAlgorithm::GaudiAlgorithm;
CalorimeterHitDigi
(
const
std
::
string
&
name
,
ISvcLocator
*
svcLoc
)
:
GaudiAlgorithm
(
name
,
svcLoc
)
...
...
@@ -70,8 +67,13 @@ public:
}
// set energy resolution numbers
for
(
size_t
i
=
0
;
i
<
u_eRes
.
size
()
&&
i
<
3
;
++
i
)
{
r
es
[
i
]
=
u_eRes
[
i
];
eR
es
[
i
]
=
u_eRes
[
i
];
}
// using juggler internal units (GeV, mm, radian, ns)
dyRangeADC
=
m_dyRangeADC
.
value
()
/
GeV
;
tRes
=
m_tRes
.
value
()
/
ns
;
return
StatusCode
::
SUCCESS
;
}
...
...
@@ -82,15 +84,16 @@ public:
// Create output collections
auto
rawhits
=
m_outputHitCollection
.
createAndPut
();
for
(
const
auto
&
ahit
:
*
simhits
)
{
double
eres
=
std
::
sqrt
(
std
::
pow
(
m_normDist
()
*
res
[
0
]
/
sqrt
(
ahit
.
energyDeposit
()
*
m_eUnit
/
GeV
),
2
)
+
std
::
pow
(
m_normDist
()
*
res
[
1
],
2
)
+
std
::
pow
(
m_normDist
()
*
res
[
2
]
/
(
ahit
.
energyDeposit
()
*
m_eUnit
/
GeV
),
2
));
// Note: juggler internal unit of energy is GeV
double
eResRel
=
std
::
sqrt
(
std
::
pow
(
m_normDist
()
*
eRes
[
0
]
/
sqrt
(
ahit
.
energyDeposit
()),
2
)
+
std
::
pow
(
m_normDist
()
*
eRes
[
1
],
2
)
+
std
::
pow
(
m_normDist
()
*
eRes
[
2
]
/
(
ahit
.
energyDeposit
()),
2
));
double
ped
=
m_pedMeanADC
+
m_normDist
()
*
m_pedSigmaADC
;
long
long
adc
=
std
::
llround
(
ped
+
ahit
.
energyDeposit
()
*
(
1.
+
e
r
es
)
*
m_eUnit
/
m_
dyRangeADC
*
m_capADC
);
long
long
adc
=
std
::
llround
(
ped
+
ahit
.
energyDeposit
()
*
(
1.
+
e
R
es
Rel
)
/
dyRangeADC
*
m_capADC
);
eic
::
RawCalorimeterHit
rawhit
(
(
long
long
)
ahit
.
cellID
(),
(
adc
>
m_capADC
?
m_capADC
.
value
()
:
adc
),
(
double
)
ahit
.
truth
().
time
*
m_tUnit
/
ns
+
m_normDist
()
*
m_
tRes
/
ns
(
adc
>
m_capADC
.
value
()
?
m_capADC
.
value
()
:
adc
),
(
double
)
ahit
.
truth
().
time
+
m_normDist
()
*
tRes
);
rawhits
->
push_back
(
rawhit
);
}
...
...
JugReco/src/components/CalorimeterHitReco.cpp
View file @
8aebe6a3
...
...
@@ -44,7 +44,9 @@ namespace Jug::Reco {
// zero suppression values
Gaudi
::
Property
<
double
>
m_thresholdADC
{
this
,
"thresholdFactor"
,
3.0
};
Gaudi
::
Property
<
double
>
m_minEdep
{
this
,
"minimumHitEdep"
,
0.
};
// unitless counterparts of the input parameters
double
dyRangeADC
;
DataHandle
<
eic
::
RawCalorimeterHitCollection
>
m_inputHitCollection
{
"inputHitCollection"
,
Gaudi
::
DataHandle
::
Reader
,
this
};
...
...
@@ -87,6 +89,9 @@ namespace Jug::Reco {
return
StatusCode
::
FAILURE
;
}
// unitless conversion
dyRangeADC
=
m_dyRangeADC
.
value
()
/
GeV
;
// do not get the layer/sector ID if no readout class provided
if
(
m_readout
.
value
().
empty
())
{
return
StatusCode
::
SUCCESS
;
...
...
@@ -97,9 +102,11 @@ namespace Jug::Reco {
id_dec
=
id_spec
.
decoder
();
if
(
m_sectorField
.
value
().
size
())
{
sector_idx
=
id_dec
->
index
(
m_sectorField
);
info
()
<<
"Find sector field "
<<
m_sectorField
.
value
()
<<
", index = "
<<
sector_idx
<<
endmsg
;
}
if
(
m_layerField
.
value
().
size
())
{
layer_idx
=
id_dec
->
index
(
m_layerField
);
info
()
<<
"Find layer field "
<<
m_layerField
.
value
()
<<
", index = "
<<
sector_idx
<<
endmsg
;
}
}
catch
(...)
{
error
()
<<
"Failed to load ID decoder for "
<<
m_readout
<<
endmsg
;
...
...
@@ -117,7 +124,7 @@ namespace Jug::Reco {
<<
m_localDetElement
.
value
()
<<
endmsg
;
return
StatusCode
::
FAILURE
;
}
// or get from fields
// or get from fields
}
else
{
std
::
vector
<
std
::
pair
<
std
::
string
,
int
>>
fields
;
for
(
auto
&
f
:
u_localDetFields
.
value
())
{
...
...
@@ -151,20 +158,16 @@ namespace Jug::Reco {
}
// convert ADC -> energy
float
energy
=
(
rh
.
amplitude
()
-
m_pedMeanADC
)
/
(
float
)
m_capADC
*
m_dyRangeADC
;
if
(
energy
<
m_minEdep
)
{
continue
;
}
float
energy
=
(
rh
.
amplitude
()
-
m_pedMeanADC
)
/
static_cast
<
float
>
(
m_capADC
.
value
())
*
dyRangeADC
;
float
time
=
rh
.
timeStamp
();
// ns
auto
id
=
rh
.
cellID
();
int
lid
=
((
id_dec
!=
nullptr
)
&
m_layerField
.
value
().
size
())
int
lid
=
((
id_dec
!=
nullptr
)
&
&
m_layerField
.
value
().
size
())
?
static_cast
<
int
>
(
id_dec
->
get
(
id
,
layer_idx
))
:
-
1
;
int
sid
=
((
id_dec
!=
nullptr
)
&
m_sectorField
.
value
().
size
())
int
sid
=
((
id_dec
!=
nullptr
)
&
&
m_sectorField
.
value
().
size
())
?
static_cast
<
int
>
(
id_dec
->
get
(
id
,
sector_idx
))
:
-
1
;
// global positions
auto
gpos
=
m_geoSvc
->
cellIDPositionConverter
()
->
position
(
id
);
// local positions
...
...
@@ -173,9 +176,8 @@ namespace Jug::Reco {
local
=
volman
.
lookupDetElement
(
id
&
local_mask
).
nominal
();
}
auto
pos
=
local
.
worldToLocal
(
dd4hep
::
Position
(
gpos
.
x
(),
gpos
.
y
(),
gpos
.
z
()));
// auto pos =
// m_geoSvc->cellIDPositionConverter()->findContext(id)->volumePlacement().position(); cell
// dimension
// auto pos = m_geoSvc->cellIDPositionConverter()->findContext(id)->volumePlacement().position();
// cell dimension
auto
cdim
=
m_geoSvc
->
cellIDPositionConverter
()
->
cellDimensions
(
id
);
double
dim
[
3
]
=
{
0.
,
0.
,
0.
};
for
(
size_t
i
=
0
;
i
<
cdim
.
size
()
&&
i
<
3
;
++
i
)
{
...
...
JugReco/src/components/CalorimeterHitsMerger.cpp
View file @
8aebe6a3
...
...
@@ -4,8 +4,9 @@
*
* Author: Chao Peng (ANL), 03/31/2021
*/
#include
<
algorithm
>
#include
<
tuple
>
#include
<bitset>
#include
<algorithm>
#include
<unordered_map>
#include
"Gaudi/Property.h"
...
...
@@ -102,50 +103,60 @@ namespace Jug::Reco {
StatusCode
execute
()
override
{
// input collections
const
auto
&
h
its
=
*
m_inputHitCollection
.
get
();
const
auto
&
i
npu
ts
=
*
m_inputHitCollection
.
get
();
// Create output collections
auto
&
mhits
=
*
m_outputHitCollection
.
createAndPut
();
// dd4hep decoders
auto
poscon
=
m_geoSvc
->
cellIDPositionConverter
();
auto
volman
=
m_geoSvc
->
detector
()
->
volumeManager
();
auto
&
outputs
=
*
m_outputHitCollection
.
createAndPut
();
//
sum energies that has the same id
std
::
unordered_map
<
long
long
,
s
ize_t
>
merge_map
;
for
(
const
auto
&
h
:
h
its
)
{
//
find the hits that belong to the same group (for merging)
std
::
unordered_map
<
long
long
,
s
td
::
vector
<
eic
::
ConstCalorimeterHit
>
>
merge_map
;
for
(
const
auto
&
h
:
i
npu
ts
)
{
int64_t
id
=
h
.
cellID
()
&
id_mask
;
// use the reference field position
int64_t
ref_id
=
id
|
ref_mask
;
// debug() << h.cellID() << " - " << std::bitset<64>(h.cellID()) << endmsg;
auto
it
=
merge_map
.
find
(
id
);
debug
()
<<
fmt
::
format
(
"{:#064b} - {:#064b}, ref: {:#064b}"
,
h
.
cellID
(),
id
,
ref_id
)
<<
endmsg
;
if
(
it
==
merge_map
.
end
())
{
merge_map
[
id
]
=
mhits
.
size
();
// global positions
auto
gpos
=
poscon
->
position
(
ref_id
);
// local positions
auto
alignment
=
volman
.
lookupDetElement
(
ref_id
).
nominal
();
// debug() << volman.lookupDetElement(ref_id).path() << ", " <<
// volman.lookupDetector(ref_id).path() << endmsg;
auto
pos
=
alignment
.
worldToLocal
(
dd4hep
::
Position
(
gpos
.
x
(),
gpos
.
y
(),
gpos
.
z
()));
mhits
.
push_back
({
ref_id
,
h
.
clusterID
(),
h
.
layerID
(),
h
.
sectorID
(),
h
.
energy
(),
h
.
time
(),
{
gpos
.
x
()
/
dd4hep
::
mm
,
gpos
.
y
()
/
dd4hep
::
mm
,
gpos
.
z
()
/
dd4hep
::
mm
},
{
pos
.
x
()
/
dd4hep
::
mm
,
pos
.
y
()
/
dd4hep
::
mm
,
pos
.
z
()
/
dd4hep
::
mm
},
h
.
dimension
(),
h
.
type
()});
merge_map
[
id
]
=
{
h
};
}
else
{
mhits
[
it
->
second
].
energy
(
mhits
[
it
->
second
].
energy
()
+
h
.
energy
());
it
->
second
.
push_back
(
h
);
}
}
// reconstruct info for merged hits
// dd4hep decoders
auto
poscon
=
m_geoSvc
->
cellIDPositionConverter
();
auto
volman
=
m_geoSvc
->
detector
()
->
volumeManager
();
for
(
auto
&
[
id
,
hits
]
:
merge_map
)
{
// reference fields id
int64_t
ref_id
=
id
|
ref_mask
;
// global positions
auto
gpos
=
poscon
->
position
(
ref_id
);
// local positions
auto
alignment
=
volman
.
lookupDetElement
(
ref_id
).
nominal
();
auto
pos
=
alignment
.
worldToLocal
(
dd4hep
::
Position
(
gpos
.
x
(),
gpos
.
y
(),
gpos
.
z
()));
// debug() << volman.lookupDetElement(ref_id).path() << ", "
// << volman.lookupDetector(ref_id).path() << endmsg;
// sum energy
float
energy
=
0.
;
for
(
auto
&
hit
:
hits
)
{
energy
+=
hit
.
energy
();
// debug() << fmt::format("{:#064b} - {:#064b}, ref: {:#064b}", hit.cellID(), id, ref_id)
// << endmsg;
}
const
auto
&
href
=
hits
.
front
();
outputs
.
push_back
(
eic
::
CalorimeterHit
{
ref_id
,
href
.
clusterID
(),
href
.
layerID
(),
href
.
sectorID
(),
energy
,
href
.
time
(),
{
gpos
.
x
()
/
dd4hep
::
mm
,
gpos
.
y
()
/
dd4hep
::
mm
,
gpos
.
z
()
/
dd4hep
::
mm
},
{
pos
.
x
()
/
dd4hep
::
mm
,
pos
.
y
()
/
dd4hep
::
mm
,
pos
.
z
()
/
dd4hep
::
mm
},
href
.
dimension
(),
href
.
type
()});
}
debug
()
<<
"Size before = "
<<
h
its
.
size
()
<<
", after = "
<<
mhi
ts
.
size
()
<<
endmsg
;
debug
()
<<
"Size before = "
<<
i
npu
ts
.
size
()
<<
", after = "
<<
outpu
ts
.
size
()
<<
endmsg
;
return
StatusCode
::
SUCCESS
;
}
...
...
JugReco/src/components/CalorimeterIslandCluster.cpp
View file @
8aebe6a3
/*
* Island Clustering Algorithm for Calorimeter Blocks
* 1. group all the adjacent modules
* 2. split the groups between their local maxima with the energy deposit above
*
<minClusterCenterEdep>
* 2. split the groups between their local maxima with the energy deposit above
<minClusterCenterEdep>
*
Output hits collection with clusterID
*
* Author: Chao Peng (ANL), 09/27/2020
* References:
...
...
@@ -11,6 +11,7 @@
*/
#include
<algorithm>
#include
"fmt/format.h"
#include
"Gaudi/Property.h"
#include
"GaudiAlg/GaudiAlgorithm.h"
#include
"GaudiAlg/GaudiTool.h"
...
...
@@ -37,14 +38,18 @@ namespace Jug::Reco {
class
CalorimeterIslandCluster
:
public
GaudiAlgorithm
{
public:
Gaudi
::
Property
<
bool
>
m_splitCluster
{
this
,
"splitCluster"
,
true
};
Gaudi
::
Property
<
double
>
m_groupRange
{
this
,
"groupRange"
,
1.8
};
Gaudi
::
Property
<
std
::
vector
<
double
>>
u_groupRanges
{
this
,
"groupRanges"
,
{}};
Gaudi
::
Property
<
double
>
m_minClusterCenterEdep
{
this
,
"minClusterCenterEdep"
,
50.0
*
MeV
};
DataHandle
<
eic
::
CalorimeterHitCollection
>
m_inputHitCollection
{
"inputHitCollection"
,
Gaudi
::
DataHandle
::
Reader
,
this
};
DataHandle
<
eic
::
CalorimeterHitCollection
>
m_splitHitCollection
{
"outputHitCollection"
,
Gaudi
::
DataHandle
::
Writer
,
this
};
Gaudi
::
Property
<
bool
>
m_splitCluster
{
this
,
"splitCluster"
,
true
};
Gaudi
::
Property
<
double
>
m_dimScaledDist
{
this
,
"dimScaledDist"
,
1.8
};
Gaudi
::
Property
<
std
::
vector
<
double
>>
u_localDistXY
{
this
,
"localDistXY"
,
{}};
Gaudi
::
Property
<
double
>
m_minClusterHitEdep
{
this
,
"minClusterHitEdep"
,
0.
};
Gaudi
::
Property
<
double
>
m_minClusterCenterEdep
{
this
,
"minClusterCenterEdep"
,
50.0
*
MeV
};
DataHandle
<
eic
::
CalorimeterHitCollection
>
m_inputHitCollection
{
"inputHitCollection"
,
Gaudi
::
DataHandle
::
Reader
,
this
};
DataHandle
<
eic
::
CalorimeterHitCollection
>
m_splitHitCollection
{
"outputHitCollection"
,
Gaudi
::
DataHandle
::
Writer
,
this
};
// unitless counterparts of the input parameters
double
minClusterHitEdep
,
minClusterCenterEdep
,
localDistXY
[
2
]
=
{
0.
,
0.
};
CalorimeterIslandCluster
(
const
std
::
string
&
name
,
ISvcLocator
*
svcLoc
)
:
GaudiAlgorithm
(
name
,
svcLoc
)
...
...
@@ -58,6 +63,21 @@ namespace Jug::Reco {
if
(
GaudiAlgorithm
::
initialize
().
isFailure
())
{
return
StatusCode
::
FAILURE
;
}
// unitless conversion, keep consistency with juggler internal units (GeV, mm, ns, rad)
minClusterHitEdep
=
m_minClusterHitEdep
.
value
()
/
GeV
;
minClusterCenterEdep
=
m_minClusterCenterEdep
.
value
()
/
GeV
;
// sanity check
if
(
u_localDistXY
.
size
()
==
2
)
{
localDistXY
[
0
]
=
u_localDistXY
.
value
()[
0
]
/
mm
;
localDistXY
[
1
]
=
u_localDistXY
.
value
()[
1
]
/
mm
;
info
()
<<
fmt
::
format
(
"Clustering using distance [x, y] <= [{:.4f} mm, {:.4f} mm]"
,
localDistXY
[
0
],
localDistXY
[
1
])
<<
endmsg
;
}
else
if
(
u_localDistXY
.
size
()
>
0
)
{
warning
()
<<
fmt
::
format
(
"Expect two values (x, y) from localDistXY, got {}. "
"Using dimScaledDist = {} instead"
,
u_localDistXY
.
size
(),
m_dimScaledDist
.
value
())
<<
endmsg
;
}
return
StatusCode
::
SUCCESS
;
}
...
...
@@ -73,6 +93,11 @@ namespace Jug::Reco {
std
::
vector
<
bool
>
visits
(
hits
.
size
(),
false
);
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
(),
hits
[
i
].
layerID
(),
hits
[
i
].
sectorID
())
<<
endmsg
;
// already in a group
if
(
visits
[
i
])
{
continue
;
...
...
@@ -81,14 +106,16 @@ namespace Jug::Reco {
// create a new group, and group all the neighboring hits
dfs_group
(
groups
.
back
(),
i
,
hits
,
visits
);
}
debug
()
<<
"we have "
<<
groups
.
size
()
<<
" groups of hits"
<<
endmsg
;
size_t
clusterID
=
0
;
for
(
auto
&
group
:
groups
)
{
auto
maxima
=
find_maxima
(
group
,
!
m_splitCluster
);
if
(
group
.
empty
())
{
continue
;
}
auto
maxima
=
find_maxima
(
group
,
!
m_splitCluster
.
value
());
split_group
(
group
,
maxima
,
clusterID
,
clustered_hits
);
debug
()
<<
"hits in a group: "
<<
group
.
size
()
<<
", "
<<
"local maxima: "
<<
maxima
.
hits_
size
()
<<
endmsg
;
<<
"local maxima: "
<<
maxima
.
size
()
<<
endmsg
;
debug
()
<<
"total number of clusters so far: "
<<
clusterID
<<
", "
<<
endmsg
;
}
...
...
@@ -112,23 +139,23 @@ namespace Jug::Reco {
{
// in the same sector
if
(
h1
.
sectorID
()
==
h2
.
sectorID
())
{
if
(
u_groupRanges
.
size
()
>=
2
)
{
return
(
std
::
abs
(
h1
.
local_x
()
-
h2
.
local_x
())
<=
u_groupRanges
[
0
])
&&
(
std
::
abs
(
h1
.
local_y
()
-
h2
.
local_y
())
<=
u_groupRanges
[
1
]);
// use absolute value
if
(
u_localDistXY
.
size
()
==
2
)
{
return
(
std
::
abs
(
h1
.
local_x
()
-
h2
.
local_x
())
<=
localDistXY
[
0
])
&&
(
std
::
abs
(
h1
.
local_y
()
-
h2
.
local_y
())
<=
localDistXY
[
1
]);
// use scaled value (module size as the scales)
}
else
{
return
(
std
::
abs
(
h1
.
local_x
()
-
h2
.
local_x
())
<=
m_
groupRange
*
(
h1
.
dim_x
()
+
h2
.
dim_x
())
/
2.
)
&&
m_
dimScaledDist
*
(
h1
.
dim_x
()
+
h2
.
dim_x
())
/
2.
)
&&
(
std
::
abs
(
h1
.
local_y
()
-
h2
.
local_y
())
<=
m_
groupRange
*
(
h1
.
dim_y
()
+
h2
.
dim_y
())
/
2.
);
m_
dimScaledDist
*
(
h1
.
dim_y
()
+
h2
.
dim_y
())
/
2.
);
}
// different sector
// different sector
, local coordinates do not work, using global coordinates
}
else
{
double
range
=
(
u_groupRanges
.
size
()
>=
2
)
?
std
::
sqrt
(
std
::
accumulate
(
u_groupRanges
.
begin
(),
u_groupRanges
.
end
(),
0.
,
[](
double
s
,
double
x
)
{
return
s
+
x
*
x
;
}))
:
m_groupRange
*
dist2d
((
h1
.
dim_x
()
+
h2
.
dim_x
())
/
2.
,
(
h1
.
dim_y
()
+
h2
.
dim_y
())
/
2.
);
(
u_localDistXY
.
size
()
==
2
)
?
dist2d
(
localDistXY
[
0
],
localDistXY
[
1
])
:
m_dimScaledDist
*
dist2d
((
h1
.
dim_x
()
+
h2
.
dim_x
())
/
2.
,
(
h1
.
dim_y
()
+
h2
.
dim_y
())
/
2.
);
// sector may have rotation (barrel), so z is included
return
dist3d
(
h1
.
x
()
-
h2
.
x
(),
h1
.
y
()
-
h2
.
y
(),
h1
.
z
()
-
h2
.
z
())
<=
range
;
}
...
...
@@ -138,6 +165,12 @@ namespace Jug::Reco {
void
dfs_group
(
std
::
vector
<
eic
::
CalorimeterHit
>&
group
,
int
idx
,
const
eic
::
CalorimeterHitCollection
&
hits
,
std
::
vector
<
bool
>&
visits
)
const
{
// not a qualified hit to particpate clustering, stop here
if
(
hits
[
idx
].
energy
()
<
minClusterHitEdep
)
{
visits
[
idx
]
=
true
;
return
;
}
eic
::
CalorimeterHit
hit
{
hits
[
idx
].
cellID
(),
hits
[
idx
].
clusterID
(),
hits
[
idx
].
layerID
(),
hits
[
idx
].
sectorID
(),
hits
[
idx
].
energy
(),
hits
[
idx
].
time
(),
...
...
@@ -154,10 +187,10 @@ namespace Jug::Reco {
}
// find local maxima that above a certain threshold
eic
::
Cluster
find_maxima
(
const
std
::
vector
<
eic
::
CalorimeterHit
>&
group
,
bool
global
=
false
)
const
std
::
vector
<
eic
::
ConstCalorimeterHit
>
find_maxima
(
const
std
::
vector
<
eic
::
CalorimeterHit
>&
group
,
bool
global
=
false
)
const
{
eic
::
Cluster
maxima
;
std
::
vector
<
eic
::
ConstCalorimeterHit
>
maxima
;
if
(
group
.
empty
())
{
return
maxima
;
}
...
...
@@ -169,13 +202,15 @@ namespace Jug::Reco {
mpos
=
i
;
}
}
maxima
.
addhits
(
group
[
mpos
]);
if
(
group
[
mpos
].
energy
()
>=
minClusterCenterEdep
)
{
maxima
.
push_back
(
group
[
mpos
]);
}
return
maxima
;
}
for
(
auto
&
hit
:
group
)
{
// not a qualified center
if
(
hit
.
energy
()
<
m_
minClusterCenterEdep
/
GeV
)
{
if
(
hit
.
energy
()
<
minClusterCenterEdep
)
{
continue
;
}
...
...
@@ -191,7 +226,7 @@ namespace Jug::Reco {
}
if
(
maximum
)
{
maxima
.
addhits
(
hit
);
maxima
.
push_back
(
hit
);
}
}
...
...
@@ -212,13 +247,14 @@ namespace Jug::Reco {
// split a group of hits according to the local maxima
// split_hits is used to persistify the data
void
split_group
(
std
::
vector
<
eic
::
CalorimeterHit
>&
group
,
const
eic
::
Cluster
&
maxima
,
void
split_group
(
std
::
vector
<
eic
::
CalorimeterHit
>&
group
,
const
std
::
vector
<
eic
::
ConstCalorimeterHit
>&
maxima
,
size_t
&
clusterID
,
eic
::
CalorimeterHitCollection
&
clustered_hits
)
const
{
// special cases
if
(
maxima
.
hits_
size
()
==
0
)
{
if
(
maxima
.
size
()
==
0
)
{
return
;
}
else
if
(
maxima
.
hits_
size
()
==
1
)
{
}
else
if
(
maxima
.
size
()
==
1
)
{
for
(
auto
&
hit
:
group
)
{
hit
.
clusterID
(
clusterID
);
clustered_hits
.
push_back
(
hit
);
...
...
@@ -228,15 +264,15 @@ namespace Jug::Reco {
}
// split between maxima
std
::
vector
<
double
>
weights
(
maxima
.
hits_
size
());
std
::
vector
<
eic
::
Cluster
>
splits
(
maxima
.
hits_
size
());
std
::
vector
<
double
>
weights
(
maxima
.
size
());
std
::
vector
<
eic
::
Cluster
>
splits
(
maxima
.
size
());
size_t
n_clus
=
clusterID
+
1
;
size_t
i
=
0
;
for
(
auto
it
=
group
.
begin
();
it
!=
group
.
end
();
++
it
,
++
i
)
{
auto
hedep
=
it
->
energy
();
size_t
j
=
0
;
// calculate weights for local maxima
for
(
auto
cit
=
maxima
.
hits_
begin
();
cit
!=
maxima
.
hits_
end
();
++
cit
,
++
j
)
{
for
(
auto
cit
=
maxima
.
begin
();
cit
!=
maxima
.
end
();
++
cit
,
++
j
)
{
double
dist_ref
=
cit
->
dim_x
();
double
energy
=
cit
->
energy
();
double
dist
=
std
::
sqrt
(
std
::
pow
(
it
->
local_x
()
-
cit
->
local_x
(),
2
)
+
...
...
JugReco/src/components/ImagingClusterReco.cpp
View file @
8aebe6a3
...
...
@@ -102,7 +102,7 @@ namespace Jug::Reco {
// debug output
for
(
const
auto
&
cl
:
clusters
)
{
debug
()
<<
fmt
::
format
(
"Cluster {:d}: Edep = {:.3f} MeV, Dir = ({:.3f}, {:.3f}) deg"
,
cl
.
clusterID
(),
cl
.
edep
()
/
MeV
,
cl
.
cl_theta
()
/
M_PI
*
180.
,
cl
.
clusterID
(),
cl
.
edep
()
*
1000.
,
cl
.
cl_theta
()
/
M_PI
*
180.
,
cl
.
cl_phi
()
/
M_PI
*
180.
)
<<
endmsg
;
}
...
...
JugReco/src/components/ImagingPixelReco.cpp
View file @
8aebe6a3
...
...
@@ -44,6 +44,10 @@ namespace Jug::Reco {
Gaudi
::
Property
<
double
>
m_dyRangeADC
{
this
,
"dynamicRangeADC"
,
100
*
MeV
};
Gaudi
::
Property
<
double
>
m_pedSigmaADC
{
this
,
"pedestalSigma"
,
3.2
};
Gaudi
::
Property
<
double
>
m_thresholdADC
{
this
,
"thresholdFactor"
,
3.0
};
// unitless counterparts for the input parameters
double
dyRangeADC
;
// hits containers
DataHandle
<
eic
::
RawCalorimeterHitCollection
>
m_inputHitCollection
{
"inputHitCollection"
,
Gaudi
::
DataHandle
::
Reader
,
this
};
...
...
@@ -90,6 +94,9 @@ namespace Jug::Reco {
return
StatusCode
::
FAILURE
;
}
// unitless conversion
dyRangeADC
=
m_dyRangeADC
.
value
()
/
GeV
;
return
StatusCode
::
SUCCESS
;
}
...
...
@@ -106,8 +113,7 @@ namespace Jug::Reco {
if
((
rh
.
amplitude
()
-
m_pedMeanADC
)
<
m_thresholdADC
*
m_pedSigmaADC
)
{
continue
;
}
double
edep
=
(
rh
.
amplitude
()
-
m_pedMeanADC
)
/
(
double
)
m_capADC
*
m_dyRangeADC
;
// convert ADC -> energy
double
edep
=
(
rh
.
amplitude
()
-
m_pedMeanADC
)
/
(
double
)
m_capADC
*
dyRangeADC
;
// convert ADC -> energy
double
time
=
rh
.
timeStamp
();
// ns
auto
id
=
rh
.
cellID
();
int
lid
=
(
int
)
id_dec
->
get
(
id
,
layer_idx
);
...
...
JugReco/src/components/ImagingTopoCluster.cpp
View file @
8aebe6a3
...
...
@@ -37,20 +37,22 @@ namespace Jug::Reco {
class
ImagingTopoCluster
:
public
GaudiAlgorithm
{
public:
// maximum difference in layer numbers that can be considered as neighbours
Gaudi
::
Property
<
int
>
m_
adjLayerDiff
{
this
,
"adjLayerDiff
"
,
1
};
Gaudi
::
Property
<
int
>
m_
neighbourLayersRange
{
this
,
"neighbourLayersRange
"
,
1
};
// maximum distance of local (x, y) to be considered as neighbors at the same layer
Gaudi
::
Property
<
std
::
vector
<
double
>>
u_local
Ranges
{
this
,
"local
Ranges
"
,
{
1.0
*
mm
,
1.0
*
mm
}};
Gaudi
::
Property
<
std
::
vector
<
double
>>
u_local
DistXY
{
this
,
"local
DistXY
"
,
{
1.0
*
mm
,
1.0
*
mm
}};
// maximum distance of global (eta, phi) to be considered as neighbors at different layers
Gaudi
::
Property
<
std
::
vector
<
double
>>
u_adjLayerRanges
{
this
,
"adjLayerRanges"
,
{
0.01
*
M_PI
,
0.01
*
M_PI
}};
Gaudi
::
Property
<
std
::
vector
<
double
>>
u_layerDistEtaPhi
{
this
,
"layerDistEtaPhi"
,
{
0.01
,
0.01
}};
// maximum global distance to be considered as neighbors in different sectors
Gaudi
::
Property
<
double
>
m_adjSectorDist
{
this
,
"adjSectorDist"
,
1.0
*
cm
};
Gaudi
::
Property
<
double
>
m_sectorDist
{
this
,
"sectorDist"
,
1.0
*
cm
};
// minimum hit energy to participate clustering
Gaudi
::
Property
<
double
>
m_minClusterHitEdep
{
this
,
"minClusterHitEdep"
,
0.
};
// minimum cluster center energy (to be considered as a seed for cluster)
Gaudi
::
Property
<
double
>
m_minClusterCenterEdep
{
this
,
"minClusterCenterEdep"
,
0.
};
// minimum cluster energy (to save this cluster)
Gaudi
::
Property
<
double
>
m_minEdep
{
this
,
"minClusterEdep"
,
0.5
*
MeV
};
Gaudi
::
Property
<
double
>
m_min
Cluster
Edep
{
this
,
"minClusterEdep"
,
0.5
*
MeV
};
// minimum number of hits (to save this cluster)
Gaudi
::
Property
<
int
>
m_minNhits
{
this
,
"minClusterNhits"
,
10
};
Gaudi
::
Property
<
int
>
m_min
Cluster
Nhits
{
this
,
"minClusterNhits"
,
10
};
// input hits collection
DataHandle
<
eic
::
ImagingPixelCollection
>
m_inputHitCollection
{
"inputHitCollection"
,
Gaudi
::
DataHandle
::
Reader
,
this
};
...
...
@@ -58,6 +60,10 @@ namespace Jug::Reco {
DataHandle
<
eic
::
ImagingPixelCollection
>
m_outputHitCollection
{
"outputHitCollection"
,
Gaudi
::
DataHandle
::
Writer
,
this
};
// unitless counterparts of the input parameters
double
localDistXY
[
2
],
layerDistEtaPhi
[
2
],
sectorDist
;
double
minClusterHitEdep
,
minClusterCenterEdep
,
minClusterEdep
,
minClusterNhits
;
ImagingTopoCluster
(
const
std
::
string
&
name
,
ISvcLocator
*
svcLoc
)
:
GaudiAlgorithm
(
name
,
svcLoc
)
{
declareProperty
(
"inputHitCollection"
,
m_inputHitCollection
,
""
);
...
...
@@ -70,25 +76,37 @@ namespace Jug::Reco {
return
StatusCode
::
FAILURE
;
}
//
check local clustering range
if
(
u_localRanges
.
size
()
<
2
)
{
error
()
<<
"Need 2-dimensional ranges for same-layer clustering, only have "
<<
u_localRanges
.
size
()
<<
endmsg
;
//
unitless conversion
// sanity checks
if
(
u_localDistXY
.
size
()
!=
2
)
{
error
()
<<
"Expected 2 values (x_dist, y_dist) for localDistXY"
<<
endmsg
;
return
StatusCode
::
FAILURE
;
}
info
()
<<
"Same layer hits group ranges"
<<
" ("
<<
u_localRanges
[
0
]
/
mm
<<
" mm, "
<<
u_localRanges
[
1
]
/
mm
<<
" mm)"
<<
endmsg
;
// check adjacent layer clustering range
if
(
u_adjLayerRanges
.
size
()
<
2
)
{
error
()
<<
"Need 2-dimensional ranges for adjacent-layer clustering, only have "
<<
u_adjLayerRanges
.
size
()
<<
endmsg
;
if
(
u_layerDistEtaPhi
.
size
()
!=
2
)
{
error
()
<<
"Expected 2 values (eta_dist, phi_dist) for layerDistEtaPhi"
<<
endmsg
;
return
StatusCode
::
FAILURE
;
}
info
()
<<
"Same layer hits group ranges"
<<
" ("
<<
u_adjLayerRanges
[
0
]
/
M_PI
*
1000.
<<
" mrad, "
<<
u_adjLayerRanges
[
1
]
/
M_PI
*
1000.
<<
" mrad)"
<<
endmsg
;
// using juggler internal units (GeV, mm, ns, rad)