Skip to content

Commit

Permalink
Optimize St4-5 clustering - preparing comparaison between the 2 algor…
Browse files Browse the repository at this point in the history 
…ithms (AliceO2Group#10457)

[MCHClustering] Optimize St4-5 clustering - Switch to test in similar conditions the current (production) release to the new one.
  • Loading branch information
@grasseau @mwinn2
grasseau authored and mwinn2 committed Jan 17, 2023
1 parent d903634 commit 23ee892
Show file tree
Hide file tree
Showing 26 changed files with 5,029 additions and 904 deletions.
6 changes: 2 additions & 4 deletions Detectors/MUON/MCH/Clustering/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,11 +22,9 @@ o2_target_root_dictionary(MCHClustering
HEADERS include/MCHClustering/ClusterizerParam.h)

o2_add_library(MCHClusteringGEM
SOURCES src/ClusterOriginal.cxx
SOURCES src/ClusterConfig.cxx
src/ClusterDump.cxx
src/ClusterFinderOriginal.cxx
src/ClusterFinderGEM.cxx
src/ClusterizerParam.cxx
src/clusterProcessing.cxx
src/ClusterPEM.cxx
src/mathUtil.cxx
Expand All @@ -36,6 +34,6 @@ o2_add_library(MCHClusteringGEM
src/poissonEM.h
src/poissonEM.cxx
src/InspectModel.cxx
PUBLIC_LINK_LIBRARIES GSL::gsl O2::MCHMappingInterface O2::MCHBase O2::MCHPreClustering
PUBLIC_LINK_LIBRARIES GSL::gsl O2::MCHMappingInterface O2::MCHBase O2::MCHPreClustering O2::MCHClustering
O2::Framework O2::CommonUtils)

49 changes: 33 additions & 16 deletions Detectors/MUON/MCH/Clustering/include/MCHClustering/ClusterConfig.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,14 +29,30 @@ struct ClusterConfig {
//
// Physical-Numerical parameters
//
// Run2
// 4.f * 0.22875f;
double minChargeOfPads; // Lowest Charge of a Pad
double minChargeOfClusterPerCathode; // Lowest Charge of a Pad
// static constexpr double minChargeOfClusterPerCathode = 1.1; // Lowest Charge of a Group
static constexpr double minChargeOfClusterPerCathode = 20.0; // Lowest Charge of a Group
//
// Algorithm limitations
// Run3
// static double minChargeOfPads = 16; // Lowest Charge of a Pad
// static double minChargeOfClusterPerCathode = 1.0 * minChargeOfPads; // Lowest Charge of a Group
//
// ClusterResolution
float SDefaultClusterResolution; ///< default cluster resolution (cm)
float SBadClusterResolution; ///< bad (e.g. mono-cathode) cluster resolution (cm)

// Large Clusters
int nbrPadLimit = 600;
double ratioStepForLargeCluster = 0.05; // increment to find nPads < nbrPadLimit
// Limit of pad number to perform the fitting
static constexpr int nbrOfPadsLimitForTheFitting = 100;
int nbrOfPadsLimitForTheFitting = 100;
// Stop the fitting if small xy shift
double minFittingXYStep = 0.1; // in cm
//
// Algorithm choices
//
int useSpline = 0;
// Logs
//
enum VerboseMode {
Expand All @@ -45,13 +61,13 @@ struct ClusterConfig {
detail = 0x2, ///< Describes in detail
debug = 0x3 ///< Ful details
};
static constexpr VerboseMode fittingLog = no;
static constexpr VerboseMode processingLog = no; // Global
static constexpr VerboseMode padMappingLog = no;
static constexpr VerboseMode groupsLog = no;
static constexpr VerboseMode EMLocalMaxLog = no;
static constexpr VerboseMode inspectModelLog = no;
static constexpr VerboseMode laplacianLocalMaxLog = no;
VerboseMode fittingLog = no;
VerboseMode processingLog = no; // Global
VerboseMode padMappingLog = no;
VerboseMode groupsLog = no;
VerboseMode EMLocalMaxLog = no;
VerboseMode inspectModelLog = no;
VerboseMode laplacianLocalMaxLog = no;
//
// Checks
//
Expand All @@ -60,14 +76,15 @@ struct ClusterConfig {
active = 0x1, ///< Describe default activation
};
// Activate/deactivate InspectModel
static constexpr ActivateMode inspectModel = active;
ActivateMode inspectModel = inactive;
//
static constexpr bool groupsCheck = true;
static constexpr bool padMappingCheck = true;
// TODO ???
// Check, Stat
bool groupsCheck = true;
bool padMappingCheck = true;
bool mathiesonCheck = false;
};

void initClusterConfig();

} // namespace mch
} // end namespace o2

Expand Down
9 changes: 6 additions & 3 deletions Detectors/MUON/MCH/Clustering/include/MCHClustering/ClusterFinderGEM.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -61,7 +61,7 @@ class ClusterFinderGEM
// GG method called by the process workflow ( ClusterFinderGEMSpec )
//

void init(int mode);
void init(int mode, bool run2Config);
void deinit();
void reset();
void fillGEMInputData(gsl::span<const Digit>& digits, uint16_t bunchCrossing, uint32_t orbit, uint32_t iPreCluster);
Expand All @@ -87,8 +87,11 @@ class ClusterFinderGEM
static constexpr int SNFitClustersMax = 3; ///< maximum number of clusters fitted at the same time
static constexpr int SNFitParamMax = 3 * SNFitClustersMax - 1; ///< maximum number of fit parameters
static constexpr double SLowestCoupling = 1.e-2; ///< minimum coupling between clusters of pixels and pads
static constexpr float SDefaultClusterResolution = 0.2f; ///< default cluster resolution (cm)
static constexpr float SBadClusterResolution = 10.f; ///< bad (e.g. mono-cathode) cluster resolution (cm)

// Invalid ???
// static constexpr char statFileName[] = "statistics.csv";
// std::fstream statStream;

// GG Unused
// void resetPreCluster(gsl::span<const Digit>& digits);
// void simplifyPreCluster(std::vector<int>& removedDigits);
Expand Down
34 changes: 30 additions & 4 deletions Detectors/MUON/MCH/Clustering/include/MCHClustering/ClusterPEM.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,13 +18,16 @@
#ifndef O2_MCH_CLUSTER_H_
#define O2_MCH_CLUSTER_H_

#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>

#include "MCHClustering/PadsPEM.h"

namespace o2
{
namespace mch
{
typedef std::pair<int, const double*> DataBlock_t;
typedef std::pair<int, double*> DataBlock_t;

typedef struct {
PadIdx_t i;
Expand All @@ -45,20 +48,32 @@ class ClusterPEM
~ClusterPEM();
inline int getNbrOfPads(int c)
{
return (pads[c] == nullptr ? 0 : pads[c]->getNbrOfPads());
return ((pads[c] == nullptr) ? 0 : pads[c]->getNbrOfPads());
};
inline int getNbrOfObsPads(int c)
{
return ((pads[c] == nullptr) ? 0 : pads[c]->getNbrOfObsPads());
};
inline int getNbrOfPads()
{
return getNbrOfPads(0) + getNbrOfPads(1);
};
inline int getNbrOfObsPads()
{
return getNbrOfObsPads(0) + getNbrOfObsPads(1);
};
inline double getTotalCharge(int c)
{
return (pads[c] == nullptr ? 0 : pads[c]->getTotalCharge());
};

inline const double* getCharges(int c)
{
return (pads[c] == nullptr) ? nullptr : pads[c]->getCharges();
return (pads[c] == nullptr ? nullptr : pads[c]->getCharges());
}

double getMaxCharge();

inline const Pads* getPads(int c) { return pads[c]; };
inline const Groups_t* getCathGroup(int c) { return cathGroup[c]; };
inline Groups_t* getProjPadGroup() { return projPadToGrp; };
Expand All @@ -67,6 +82,9 @@ class ClusterPEM
{
return (projectedPads == nullptr ? -1 : projectedPads->getNbrOfPads());
};
std::pair<double, double> computeChargeBarycenter(int plane);
//
std::pair<int, int> getNxNy(int c);
// Unused - Old version
// double *getProjPadsAsXYdXY( Groups_t group, const Mask_t* maskGrp, int
// nbrProjPadsInTheGroup);
Expand All @@ -85,6 +103,8 @@ class ClusterPEM
void addBoundaryPads();
// Find local maximima with the PET algo
int findLocalMaxWithPEM(double* thetaL, int nbrOfPadsInTheGroupCath);
int findLocalMaxWithPEMFullRefinement(double* thetaL, int nbrOfPadsInTheGroupCath);
int findLocalMaxWithPEM2Lev(double* thetaL, int nbrOfPadsInTheGroupCath);
// Perform the fitting
DataBlock_t fit(double* thetaInit, int K);
// Not used in the Clustering/fitting
Expand Down Expand Up @@ -144,6 +164,9 @@ class ClusterPEM
int renumberGroups(Groups_t* grpToGrp, int nGrp);
// Remove low charged groups
void removeLowChargedGroups(int nGroups);
// Remove smallCharged seeds
int filterFitModelOnSmallChargedSeeds(Pads& pads, double* theta, int K,
Mask_t* maskFilteredTheta);
// Keep the seeds inside the cluster area
// Some fitting cases provide seeds outside of the cluster area
int filterFitModelOnClusterRegion(Pads& pads, double* theta, int K,
Expand All @@ -152,7 +175,8 @@ class ClusterPEM
int filterFitModelOnSpaceVariations(const double* theta0, int K0,
double* theta, int K,
Mask_t* maskFilteredTheta);

Pads* findLocalMaxWithRefinement(double* thetaL, int nbrOfPadsInTheGroupCath);
Pads* findLocalMaxWithoutRefinement(double* thetaL, int nbrOfPadsInTheGroupCath);
// ???
int getIndexByRow(const char* matrix, PadIdx_t N, PadIdx_t M, PadIdx_t* IIdx);
int getIndexByColumns(const char* matrix, PadIdx_t N, PadIdx_t M,
Expand All @@ -171,6 +195,8 @@ class ClusterPEM
PadIdx_t* sortedLocalMax, int kMax, double* smoothQ);
};

gsl_matrix* moore_penrose_pinv(gsl_matrix* A, double rcond);

} // namespace mch
} // namespace o2

Expand Down
2 changes: 2 additions & 0 deletions Detectors/MUON/MCH/Clustering/include/MCHClustering/ClusterizerParam.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,6 +32,8 @@ struct ClusterizerParam : public o2::conf::ConfigurableParamHelper<ClusterizerPa
double defaultClusterResolution = 0.2; ///< default cluster resolution (cm)
double badClusterResolution = 10.; ///< bad (e.g. mono-cathode) cluster resolution (cm)

bool legacy = true; ///< use original (run2) clustering

O2ParamDef(ClusterizerParam, "MCHClustering");
};

Expand Down
54 changes: 45 additions & 9 deletions Detectors/MUON/MCH/Clustering/include/MCHClustering/PadsPEM.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,8 @@
#ifndef O2_MCH_PADSPEM_H_
#define O2_MCH_PADSPEM_H_

#include <vector>

#include "MCHClustering/ClusterConfig.h"

namespace o2
Expand Down Expand Up @@ -45,14 +47,15 @@ inline static PadIdx_t getTheFirstNeighborOf(PadIdx_t* neigh, PadIdx_t i)
class Pads
{
public:
enum padMode {
xydxdyMode = 0x0, ///< x, y, dx, dy pad coordinates
xyInfSupMode = 0x1 ///< xInf=x, xSup=dx, yInf=y, ySup=dy pad coordinates
enum class PadMode {
xydxdyMode, ///< x, y, dx, dy pad coordinates
xyInfSupMode ///< xInf=x, xSup=dx, yInf=y, ySup=dy pad coordinates
};
static constexpr double epsilonGeometry =
1.0e-04; // Uncertainty on pad location (in cm)
// Representation mode (see padMode)
int mode = xydxdyMode;
// PadMode mode;
PadMode mode = PadMode::xydxdyMode;

// Utilities
static void printNeighbors(const PadIdx_t* neigh, int N);
Expand All @@ -63,15 +66,17 @@ class Pads
};

// Allocation constructor
Pads(int N, int chId, int mode = xydxdyMode);
Pads(int N, int chId, PadMode mode = PadMode::xydxdyMode);
// Build a new set of pads with different coordinates
// xydxdy mode or xyInfSup
Pads(const Pads& pads, int mode_);
Pads(const Pads& pads, PadMode mode_);
// Pad object with over allocated pads
Pads(const Pads* pads, int size);
// Build a pads set from those selected by "mask"
// Used to extract sub-clusters
Pads(const Pads& pads, const Groups_t* mask);
// Concatenate the 2 pads sets
Pads(const Pads* pads1, const Pads* pads2, int mode);
Pads(const Pads* pads0, const Pads* pads1, PadMode mode);
// Main constructor
Pads(const double* x_, const double* y_, const double* dx_, const double* dy_,
const double* q_, const short* cathode, const Mask_t* saturate_,
Expand All @@ -82,6 +87,7 @@ class Pads
// Take the ownership of coordinates (x, y, dx, dy)
Pads(double* x_, double* y_, double* dx_, double* dy_, int chId, int nPads_);
inline int getNbrOfPads() const { return nPads; };
inline int getNbrOfObsPads() const { return nObsPads; };
inline const double* getX() const { return x; };
inline const double* getY() const { return y; };
inline const double* getDX() const { return dx; };
Expand All @@ -94,28 +100,51 @@ class Pads
inline const Mask_t* getSaturates() const { return saturate; };
inline const Mask_t* getCathodes() const { return cath; };
inline double getTotalCharge() const { return totalCharge; };
double updateTotalCharge();
// Mean of the 2 cathodes total Charge
double getMeanTotalCharge();
inline int getChamberId() const { return chamberId; };
void setCharges(double c);
void setCharges(double* q_, int n);
void setCathodes(Mask_t cath_);
void setSaturate(Mask_t val);
// Select/keep pads in the list index
Pads* selectPads(int* index, int k);
// Remove pads whos charge is less than qCut
int removePads(double qCut);
// pad coordinates transformations
// xydxyMode <-> xyInfSupMode
void padBoundsToCenter(const Pads& pads);
void padCenterToBounds(const Pads& pads);
void padBoundsToCenter();
void padCenterToBounds();
// Charges normalization
void normalizeCharges();
// Split each pads in 4 smaller pads with the same sizes
Pads* refinePads();
Pads* refineAll();
// refine only at the local max
void refineLocalMax(Pads& localMax, std::vector<PadIdx_t>& localMaxIdx);
// refine only at the local max
void refineLocalMaxAndUpdateCij(const Pads& pads,
std::vector<PadIdx_t>& localMaxIdx, double Cij[]);
// Add zero-charged pads to the neighboring of the pads (cathode cluster)
Pads* addBoundaryPads();
// Building Neighbors
PadIdx_t* buildFirstNeighbors();
// Building K-Neighbors
PadIdx_t* buildKFirstsNeighbors(int kernelSize);
// Extract local maximima
Pads* extractLocalMax();
Pads* extractLocalMax(std::vector<PadIdx_t>& localMaxIdx, double dxMinPadSize, double dyMinPadSize);
Pads* extractLocalMaxOnCoarsePads(std::vector<PadIdx_t>& localMaxIdx);

Pads* extractLocalMaxOnCoarsePads_Remanent(std::vector<PadIdx_t>& localMaxIdx, double dxMinPadSize, double dyMinPadSize);
// Extract local maximima, with of without a neighboring
// Obsolete
Pads* clipOnLocalMax(bool extractLocalMax);
// Groups
int addIsolatedPadInGroups(Mask_t* cathToGrp, Mask_t* grpToGrp, int nGroups);
//
// inv void print(const char* title);
~Pads();

private:
Expand All @@ -132,16 +161,23 @@ class Pads
double* q = nullptr;
double totalCharge = 0;
int nPads = 0;
// n Observable/ measurable pads
// Used to speed-up the fitting
int nObsPads = 0;
int chamberId = -1;
PadIdx_t* neighbors = nullptr;
//
// Memory allocation/deallocation
void allocate();
void allocate(int size);
void release();
void copyPads(const Pads* srcPads, int srcIdx, int destIdx, int N, Mask_t cathValue);
// Utilities
void removePad(int index);
PadIdx_t* buildFirstNeighbors(double* X, double* Y, double* DX, double* DY,
int N);
// Assess or not if xyCheck is a remanent local Max (can be removed)
bool assessRemanent(double xyCheck, double* xy, double precision, int N);
void setToZero();
};

Expand Down
2 changes: 1 addition & 1 deletion Detectors/MUON/MCH/Clustering/include/MCHClustering/clusterProcessing.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@

#include "MCHClustering/ClusterConfig.h"

typedef std::pair<int, const double*> DataBlock_t;
// ??? Inv typedef std::pair<int, double*> DataBlock_t;

namespace o2
{
Expand Down
Loading

0 comments on commit 23ee892

Please sign in to comment.