Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Make topoclustering work with new theta-module merged readout #63

Merged
merged 9 commits into from
Oct 24, 2023
Merged

Make topoclustering work with new theta-module merged readout #63

Show file tree
Hide file tree
Changes from 8 commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
cbda076
Âcomment unneeded code
giovannimarchiori Sep 27, 2023
e8f9d97
add xml files for calibration of latest FCCee ALLEGRO detector versio…
giovannimarchiori Sep 29, 2023
91fa521
move some printout before the assert to easier debugging
giovannimarchiori Oct 6, 2023
233b960
modify LAr gap thickness to have 1536 modules
giovannimarchiori Oct 6, 2023
d904109
remove extra spaces
giovannimarchiori Oct 6, 2023
dc2f183
remove debug code. Add case of diagonal neighbours
giovannimarchiori Oct 9, 2023
e641849
fix calculation of diagonal cells, simplify and document code
giovannimarchiori Oct 9, 2023
c704544
cleanup
giovannimarchiori Oct 9, 2023
68dc52e
add comment
giovannimarchiori Oct 24, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
26 changes: 23 additions & 3 deletions Detector/DetCommon/include/DetCommon/DetUtils.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@
// FCCSW
#include "DetSegmentation/FCCSWGridPhiEta.h"
#include "DetSegmentation/FCCSWGridPhiTheta.h"

#include "DetSegmentation/FCCSWGridModuleThetaMerged.h"

// DD4hep
#include "DD4hep/DetFactoryHelper.h"
Expand Down Expand Up @@ -84,6 +84,24 @@ std::vector<uint64_t> neighbours(const dd4hep::DDSegmentation::BitFieldCoder& aD
const std::vector<bool>& aFieldCyclic = {false, false, false, false},
bool aDiagonal = true);

/** Special version of the neighbours function for the readout with module and theta merged cells
* Compared to the standard version, it needs a reference to the segmentation class to
* access the number of merged cells per layer. The other parameters and return value are the same
* @param[in] aSeg Reference to the segmentation object.
* @param[in] aDecoder Handle to the bitfield decoder.
* @param[in] aFieldNames Names of the fields for which neighbours are found.
* @param[in] aFieldExtremes Minimal and maximal values for the fields.
* @param[in] aCellId ID of cell.
* @param[in] aDiagonal If diagonal neighbours should be included (all combinations of fields).
* return Vector of neighbours.
*/
std::vector<uint64_t> neighbours_ModuleThetaMerged(const dd4hep::DDSegmentation::FCCSWGridModuleThetaMerged& aSeg,
const dd4hep::DDSegmentation::BitFieldCoder& aDecoder,
const std::vector<std::string>& aFieldNames,
const std::vector<std::pair<int, int>>& aFieldExtremes,
uint64_t aCellId,
bool aDiagonal = false);

/** Get minimal and maximal values that can be decoded in the fields of the bitfield.
* @param[in] aDecoder Handle to the bitfield decoder.
* @param[in] aFieldNames Names of the fields for which extremes are found.
Expand Down Expand Up @@ -146,16 +164,18 @@ std::array<uint, 2> numberOfCells(uint64_t aVolumeId, const dd4hep::DDSegmentati
*/
std::array<uint, 3> numberOfCells(uint64_t aVolumeId, const dd4hep::DDSegmentation::CartesianGridXYZ& aSeg);

/** Get the number of cells for the volume and a given Phi-Eta segmentation.
/** Get the number of cells for the volume and a given Phi-Eta / Phi-Theta / Module-Theta segmentation.
* It is assumed that the volume has a cylindrical shape (and full azimuthal coverage)
* and that it is centred at (0,0,0).
* For an example see: Test/TestReconstruction/tests/options/testcellcountingPhiEta.py.
* @warning No offset in segmentation is currently taken into account.
* @param[in] aVolumeId The volume for which the cells are counted.
* @param[in] aSeg Handle to the segmentation of the volume.
* return Array of the number of cells in (phi, eta) and the minimum eta ID.
* return Array of the number of cells in (phi, eta) / (phi, theta) / (module, theta) and the minimum eta / theta ID.
*/
std::array<uint, 3> numberOfCells(uint64_t aVolumeId, const dd4hep::DDSegmentation::FCCSWGridPhiEta& aSeg);
std::array<uint, 3> numberOfCells(uint64_t aVolumeId, const dd4hep::DDSegmentation::FCCSWGridPhiTheta& aSeg);
std::array<uint, 3> numberOfCells(uint64_t aVolumeId, const dd4hep::DDSegmentation::FCCSWGridModuleThetaMerged& aSeg);

/** Get the number of cells for the volume and a given R-phi segmentation.
* It is assumed that the volume has a cylindrical shape - TGeoTube (and full azimuthal coverage)
Expand Down
252 changes: 240 additions & 12 deletions Detector/DetCommon/src/DetUtils.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,9 @@
#define MM_2_CM 0.1
#endif

#include <iostream>

#include <unordered_set>

namespace det {
namespace utils {
Expand Down Expand Up @@ -78,7 +81,7 @@ std::vector<std::vector<uint>> combinations(int N, int K) {

std::vector<std::vector<int>> permutations(int K) {
std::vector<std::vector<int>> indexes;
int N = pow(2, K); // number of permutations with repetition of 2 numbers (0,1)
int N = pow(2, K); // number of permutations with repetition of 2 numbers (-1,1)
for (int i = 0; i < N; i++) {
// permutation = binary representation of i
std::vector<int> tmp;
Expand All @@ -94,11 +97,13 @@ std::vector<std::vector<int>> permutations(int K) {
return indexes;
}

// use it for field module/phi
int cyclicNeighbour(int aCyclicId, std::pair<int, int> aFieldExtremes) {
int nBins = aFieldExtremes.second - aFieldExtremes.first + 1;
if (aCyclicId < aFieldExtremes.first) {
return aFieldExtremes.second + aCyclicId;
return (aFieldExtremes.second + 1 - ((aFieldExtremes.first - aCyclicId) % nBins) );
} else if (aCyclicId > aFieldExtremes.second) {
return aCyclicId % (aFieldExtremes.second + 1);
return ( ((aCyclicId - aFieldExtremes.first) % nBins) + aFieldExtremes.first) ;
}
return aCyclicId;
}
Expand All @@ -111,23 +116,27 @@ std::vector<uint64_t> neighbours(const dd4hep::DDSegmentation::BitFieldCoder& aD
dd4hep::DDSegmentation::CellID cID = aCellId;
for (uint itField = 0; itField < aFieldNames.size(); itField++) {
const auto& field = aFieldNames[itField];
dd4hep::DDSegmentation::CellID id = aDecoder.get(cID,field);
// note: get(..) returns a FieldID i.e. a int64_t
// while CellID (used previously) is a uint64_t
// similarly, the second argument of set(..)
// is signed (FieldID) rather than unsigned (CellID)
int id = aDecoder.get(cID,field);
if (aFieldCyclic[itField]) {
aDecoder[field].set(cID, cyclicNeighbour(id - 1, aFieldExtremes[itField]));
neighbours.emplace_back(cID);
aDecoder[field].set(cID, cyclicNeighbour(id + 1, aFieldExtremes[itField]));
neighbours.emplace_back(cID);
} else {
if (id > aFieldExtremes[itField].first) {
aDecoder.set(cID, field, id - 1);
aDecoder[field].set(cID, id - 1);
neighbours.emplace_back(cID);
}
if (id < aFieldExtremes[itField].second) {
aDecoder.set(cID, field, id + 1);
aDecoder[field].set(cID, id + 1);
neighbours.emplace_back(cID);
}
}
aDecoder.set(cID, field, id);
aDecoder[field].set(cID, id);
}
if (aDiagonal) {
std::vector<int> fieldIds; // initial IDs
Expand All @@ -150,7 +159,7 @@ std::vector<uint64_t> neighbours(const dd4hep::DDSegmentation::BitFieldCoder& aD
for (uint iField = 0; iField < indexes[iComb].size(); iField++) {
if (aFieldCyclic[indexes[iComb][iField]]) {
aDecoder[aFieldNames[indexes[iComb][iField]]].set(cID, cyclicNeighbour(fieldIds[indexes[iComb][iField]] + calculation[iCalc][iField],
aFieldExtremes[indexes[iComb][iField]]) );
aFieldExtremes[indexes[iComb][iField]]) );
} else if ((calculation[iCalc][iField] > 0 &&
fieldIds[indexes[iComb][iField]] < aFieldExtremes[indexes[iComb][iField]].second) ||
(calculation[iCalc][iField] < 0 &&
Expand All @@ -175,6 +184,176 @@ std::vector<uint64_t> neighbours(const dd4hep::DDSegmentation::BitFieldCoder& aD
return neighbours;
}

// use it for module-theta merged readout (FCCSWGridModuleThetaMerged)
std::vector<uint64_t> neighbours_ModuleThetaMerged(const dd4hep::DDSegmentation::FCCSWGridModuleThetaMerged& aSeg,
const dd4hep::DDSegmentation::BitFieldCoder& aDecoder,
const std::vector<std::string>& aFieldNames,
const std::vector<std::pair<int, int>>& aFieldExtremes,
uint64_t aCellId,
bool aDiagonal) {

std::vector<uint64_t> neighbours;

// check that field names and extremes have the proper length
if (aFieldNames.size() != 3 || aFieldExtremes.size()!=3) {
std::cout << "ERROR: the vectors aFieldNames and aFieldSizes should be of length = 3, corresponding to the theta/module/layer fields" << std::endl;
std::cout << "ERROR: will return empty neighbour map" << std::endl;
return neighbours;
}

// find index of layer, module and theta in field extremes vector
int idModuleField(-1);
int idThetaField(-1);
int idLayerField(-1);
for (uint itField = 0; itField < aFieldNames.size(); itField++) {
if (aFieldNames[itField] == aSeg.fieldNameModule())
idModuleField = (int) itField;
else if (aFieldNames[itField] == aSeg.fieldNameTheta())
idThetaField = (int) itField;
else if (aFieldNames[itField] == aSeg.fieldNameLayer())
idLayerField = (int) itField;

}
if (idModuleField < 0) {
std::cout << "WARNING: module field " << aSeg.fieldNameModule() << " not found in aFieldNames vector" << std::endl;
std::cout << "WARNING: will return empty neighbour map" << std::endl;
return neighbours;
}
if (idThetaField < 0) {
std::cout << "WARNING: theta field " << aSeg.fieldNameTheta() << " not found in aFieldNames vector" << std::endl;
std::cout << "WARNING: will return empty neighbour map" << std::endl;
return neighbours;
}
if (idLayerField < 0) {
std::cout << "WARNING: layer field " << aSeg.fieldNameLayer() << " not found in aFieldNames vector" << std::endl;
std::cout << "WARNING: will return empty neighbour map" << std::endl;
return neighbours;
}

// retrieve layer/module/theta of cell under study
int layer_id = aDecoder.get(aCellId, aFieldNames[idLayerField]);
int module_id = aDecoder.get(aCellId, aFieldNames[idModuleField]);
int theta_id = aDecoder.get(aCellId, aFieldNames[idThetaField]);

// now find the neighbours
dd4hep::DDSegmentation::CellID cID = aCellId;

// for neighbours across different layers, we have to take into
// account that merging along module and/or theta could be different
// so one cell in layer N could be neighbour to several in layer N+-1
// The cells are classified in a different way whether they are
// direct neighbours (common surface), diagonal neighbours (common edge or vertex)
// or neither.
// To decide this, we need to check how the cells are related in both directions:
// neighbours (edge at least partially in common), diagonal neigbours (common vertex),
// none
int neighbourTypeModuleDir; // 0: not neighbour; 1: diagonal neighbour; 2: neighbour in module direction
int neighbourTypeThetaDir; // 0: not neighbour; 1: diagonal neighbour; 2: neighbour in module direction

for (int deltaLayer = -1; deltaLayer<2; deltaLayer+=2) {

// no neighbours in layer N-1 for innermost layer
if (layer_id == aFieldExtremes[idLayerField].first && deltaLayer<0) continue;
// and in layer N+1 for outermost layer
if (layer_id == aFieldExtremes[idLayerField].second && deltaLayer>0) continue;

// set layer field of neighbour cell
aDecoder.set(cID, aSeg.fieldNameLayer(), layer_id + deltaLayer);

// find the neighbour(s) in module and theta
// if the numbers of module (theta) merged cells across 2 layers are the
// same then we just take the same module (theta) ID
// otherwise, we need to do some math to account for the different mergings
// note: even if number of merged cells in layer-1 is larger, a cell
// in layer could neighbour more than one cell in layer-1 if the merged
// cells are not aligned, for example if cells are grouped by 3 in a layer
// and by 4 in the next one, cell 435 in the former (which groups together
// 435-436-437) will be neighbour to cells 432 and 436 of the latter
// this might introduce duplicates, we will remove them later
// another issue is that it could add spurious cells beyond the maximum module number
// to prevent this we would need to know the max module number in layer -1
// which would require modifying this function passing the extrema for all layers
// instead of the extrema only for a certain layer
// this border effect is also present in the original method..
for (int i=-1; i <= aSeg.mergedThetaCells(layer_id); i++) {
int theta_id_neighbour = (theta_id + i) - ((theta_id + i) % aSeg.mergedThetaCells(layer_id+deltaLayer));
if (theta_id_neighbour >= theta_id && theta_id_neighbour < (theta_id + aSeg.mergedThetaCells(layer_id))) neighbourTypeThetaDir = 2;
else if (theta_id_neighbour < theta_id && theta_id_neighbour > (theta_id - aSeg.mergedThetaCells(layer_id+deltaLayer))) neighbourTypeThetaDir = 2;
else if (theta_id_neighbour == (theta_id + aSeg.mergedThetaCells(layer_id))) neighbourTypeThetaDir = 1;
else if (theta_id_neighbour == (theta_id - aSeg.mergedThetaCells(layer_id+deltaLayer))) neighbourTypeThetaDir = 1;
else neighbourTypeThetaDir = 0;

// if there is no point of contact along theta, no need to check also for module direction
if (neighbourTypeThetaDir == 0) continue;
// if we are not considering diagonal neighbours, and cells in theta have only an edge in common, then skip
if (!aDiagonal && neighbourTypeThetaDir == 1) continue;
// otherwise, check status along module direction
for (int j=-1; j <= aSeg.mergedModules(layer_id); j++) {
int module_id_neighbour = (module_id + j) - ((module_id + j) % aSeg.mergedModules(layer_id+deltaLayer));
int module_id_neighbour_cyclic = cyclicNeighbour(module_id_neighbour,
aFieldExtremes[idModuleField]
);

if (module_id_neighbour >= module_id && module_id_neighbour < (module_id + aSeg.mergedModules(layer_id))) neighbourTypeModuleDir = 2;
else if (module_id_neighbour < module_id && module_id_neighbour > (module_id - aSeg.mergedModules(layer_id+deltaLayer))) neighbourTypeModuleDir = 2;
else if (module_id_neighbour == (module_id + aSeg.mergedModules(layer_id))) neighbourTypeModuleDir = 1;
else if (module_id_neighbour == (module_id - aSeg.mergedModules(layer_id+deltaLayer))) neighbourTypeModuleDir = 1;
else neighbourTypeModuleDir = 0;

// if there is no point of contact along module, then skip
if (neighbourTypeModuleDir == 0) continue;
// otherwise: if neighbours along both module and theta, or along one of the two
// and we also consider diagonal neighbours, then add cells to list of neighbours
if ( (neighbourTypeModuleDir == 2 && neighbourTypeThetaDir==2) ||
(aDiagonal && neighbourTypeThetaDir > 0 && neighbourTypeModuleDir>0) ) {
aDecoder.set(cID, aSeg.fieldNameModule(), module_id_neighbour_cyclic);
aDecoder.set(cID, aSeg.fieldNameTheta(), theta_id_neighbour);
neighbours.emplace_back(cID);
}
}
}
}
// reset cellID
// aDecoder.set(cID, aSeg.fieldNameModule(), module_id);
// aDecoder.set(cID, aSeg.fieldNameTheta(), theta_id);

// for neighbours in module/theta direction at same layer_id, do +-nMergedCells instead of +-1
aDecoder.set(cID, aSeg.fieldNameLayer(), layer_id);
// loop over theta cells
for (int i=-1; i<=1; i++) {
// calculate theta_id of neighbour
int theta_id_neighbour = theta_id + i*aSeg.mergedThetaCells(layer_id);
// check that it is within the ranges
if (
(theta_id_neighbour < aFieldExtremes[idThetaField].first) ||
(theta_id_neighbour > aFieldExtremes[idThetaField].second)
) continue;
// set theta_id of cell ID
aDecoder[aSeg.fieldNameTheta()].set(cID, theta_id + i*aSeg.mergedThetaCells(layer_id));
// loop over modules
for (int j=-1; j<=1; j++) {
// skip the cell under study (i==j==0)
if (i==0 && j==0) continue;
// calculate module_id of neighbour
int newid = cyclicNeighbour(module_id + j*aSeg.mergedModules(layer_id), aFieldExtremes[idModuleField]);
newid -= (newid % aSeg.mergedModules(layer_id));
// set module_if of cell ID
aDecoder[aSeg.fieldNameModule()].set(cID, newid);
// add cell to neighbour list
if ( i==0 || j==0 || aDiagonal ) { // first two conditions correspond to non diagonal neighbours
neighbours.emplace_back(cID);
}
}
}

// remove duplicates
std::unordered_set<uint64_t> s;
for (uint64_t i : neighbours)
s.insert(i);
neighbours.assign( s.begin(), s.end() );
return neighbours;
}

std::vector<std::pair<int, int>> bitfieldExtremes(const dd4hep::DDSegmentation::BitFieldCoder& aDecoder,
const std::vector<std::string>& aFieldNames) {
std::vector<std::pair<int, int>> extremes;
Expand Down Expand Up @@ -265,20 +444,32 @@ std::array<double, 2> tubeEtaExtremes(uint64_t aVolumeId) {
// eta segmentation calculate maximum eta from the inner radius (no offset is taken into account)
double maxEta = 0;
double minEta = 0;

double rIn = sizes.x();
double rOut = sizes.y();
double dZ = sizes.z();

// check if it is a cylinder centred at z=0
dd4hep::VolumeManager volMgr = dd4hep::Detector::getInstance().volumeManager();
auto detelement = volMgr.lookupDetElement(aVolumeId);
const auto& transformMatrix = detelement.nominal().worldTransformation();
double outGlobal[3];
double inLocal[] = {0, 0, 0}; // to get middle of the volume
transformMatrix.LocalToMaster(inLocal, outGlobal);
if (outGlobal[2] < 1e-10) {
double zCenter = outGlobal[2];
if (fabs(zCenter) < 1e-10) {
// this assumes cylinder centred at z=0
maxEta = CLHEP::Hep3Vector(sizes.x(), 0, sizes.z()).eta();
maxEta = CLHEP::Hep3Vector(rIn, 0, dZ).eta();
minEta = -maxEta;
} else {
maxEta = CLHEP::Hep3Vector(sizes.x(), 0, std::max(sizes.z() + outGlobal[2], -sizes.z() + outGlobal[2])).eta();
minEta = CLHEP::Hep3Vector(sizes.y(), 0, std::min(sizes.z() + outGlobal[2], -sizes.z() + outGlobal[2])).eta();
maxEta = std::max(
CLHEP::Hep3Vector(rIn, 0, zCenter+dZ).eta(),
CLHEP::Hep3Vector(rOut, 0, zCenter+dZ).eta()
);
minEta = std::min(
CLHEP::Hep3Vector(rIn, 0, zCenter-dZ).eta(),
CLHEP::Hep3Vector(rOut, 0, zCenter-dZ).eta()
);
}
return {minEta, maxEta};
}
Expand Down Expand Up @@ -338,6 +529,43 @@ std::array<uint, 3> numberOfCells(uint64_t aVolumeId, const dd4hep::DDSegmentati
return {phiCellNumber, cellsEta, minEtaID};
}

std::array<uint, 3> numberOfCells(uint64_t aVolumeId, const dd4hep::DDSegmentation::FCCSWGridPhiTheta& aSeg) {
uint phiCellNumber = aSeg.phiBins();
double thetaCellSize = aSeg.gridSizeTheta();
auto etaExtremes = volumeEtaExtremes(aVolumeId);
double thetaMin = 2.*atan(exp(-etaExtremes[1]));
double thetaMax = 2.*atan(exp(-etaExtremes[0]));

uint cellsTheta = ceil(( thetaMax - thetaMin - thetaCellSize ) / 2 / thetaCellSize) * 2 + 1;
uint minThetaID = int(floor((thetaMin + 0.5 * thetaCellSize - aSeg.offsetTheta()) / thetaCellSize));
return {phiCellNumber, cellsTheta, minThetaID};
}

std::array<uint, 3> numberOfCells(uint64_t aVolumeId, const dd4hep::DDSegmentation::FCCSWGridModuleThetaMerged& aSeg) {

const dd4hep::DDSegmentation::BitFieldCoder* aDecoder = aSeg.decoder();
int nLayer = aDecoder->get(aVolumeId, aSeg.fieldNameLayer());
// + 0.5 to avoid integer division
uint nModules = aSeg.nModules() / aSeg.mergedModules(nLayer);
if (aSeg.nModules() % aSeg.mergedModules(nLayer) != 0) nModules++;
// get minimum and maximum theta of volume
auto etaExtremes = volumeEtaExtremes(aVolumeId);
double thetaMin = 2.*atan(exp(-etaExtremes[1]));
double thetaMax = 2.*atan(exp(-etaExtremes[0]));

// convert to minimum and maximum theta bins
double thetaCellSize = aSeg.gridSizeTheta();
double thetaOffset = aSeg.offsetTheta();
uint minThetaID = int(floor((thetaMin + 0.5 * thetaCellSize - thetaOffset) / thetaCellSize));
uint maxThetaID = int(floor((thetaMax + 0.5 * thetaCellSize - thetaOffset) / thetaCellSize));
// correct minThetaID and maxThetaID for merging
uint mergedThetaCells = aSeg.mergedThetaCells(nLayer);
minThetaID -= (minThetaID % mergedThetaCells);
maxThetaID -= (maxThetaID % mergedThetaCells);
uint nThetaCells = 1 + (maxThetaID - minThetaID)/ mergedThetaCells;
return {nModules, nThetaCells, minThetaID};
}

std::array<uint, 2> numberOfCells(uint64_t aVolumeId, const dd4hep::DDSegmentation::PolarGridRPhi& aSeg) {
// get half-widths,
auto tubeSizes = tubeDimensions(aVolumeId);
Expand Down
Loading
Loading