Skip to content

Commit

Permalink
Make topoclustering work with new theta-module merged readout (#63)
Browse files Browse the repository at this point in the history 
* Âcomment unneeded code

* add xml files for calibration of latest FCCee ALLEGRO detector version with 1545 modules and module-theta readout

* move some printout before the assert to easier debugging

* modify LAr gap thickness to have 1536 modules

* remove extra spaces

* remove debug code. Add case of diagonal neighbours

* fix calculation of diagonal cells, simplify and document code

* cleanup

* add comment
  • Loading branch information
@giovannimarchiori
giovannimarchiori committed Oct 24, 2023
1 parent 75c3e84 commit 9cad5e2
Show file tree
Hide file tree
Showing 7 changed files with 551 additions and 29 deletions.
26 changes: 23 additions & 3 deletions Detector/DetCommon/include/DetCommon/DetUtils.h
View file
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
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

0 comments on commit 9cad5e2

Please sign in to comment.