[MCHClustering] Optimize St4-5 clustering - preparing PR - option Tim…

…ingStat

Detectors/MUON/MCH/Clustering/CMakeLists.txt

@@ -23,6 +23,7 @@ o2_target_root_dictionary(MCHClustering
 
 o2_add_library(MCHClusteringGEM
                SOURCES src/ClusterOriginal.cxx
+                       src/ClusterConfig.cxx
                        src/ClusterDump.cxx
                        src/ClusterFinderOriginal.cxx
                        src/ClusterFinderGEM.cxx

Detectors/MUON/MCH/Clustering/include/MCHClustering/ClusterConfig.h

@@ -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 {
@@ -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
   //
@@ -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
 

Detectors/MUON/MCH/Clustering/include/MCHClustering/ClusterFinderGEM.h

@@ -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);
@@ -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);

Detectors/MUON/MCH/Clustering/include/MCHClustering/ClusterPEM.h

@@ -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;
@@ -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; };
@@ -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);
@@ -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
@@ -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,
@@ -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,
@@ -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
 

Detectors/MUON/MCH/Clustering/include/MCHClustering/PadsPEM.h

@@ -17,6 +17,8 @@
 #ifndef O2_MCH_PADSPEM_H_
 #define O2_MCH_PADSPEM_H_
 
+#include <vector>
+
 #include "MCHClustering/ClusterConfig.h"
 
 namespace o2
@@ -45,14 +47,15 @@ inline static PadIdx_t getTheFirstNeighborOf(PadIdx_t* neigh, PadIdx_t i)
 class Pads
 {
  public:
-  enum padMode {
+  enum PadMode {
     xydxdyMode = 0x0,  ///< x, y, dx, dy pad coordinates
     xyInfSupMode = 0x1 ///< 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 = xydxdyMode;
 
   // Utilities
   static void printNeighbors(const PadIdx_t* neigh, int N);
@@ -63,15 +66,17 @@ class Pads
   };
 
   // Allocation constructor
-  Pads(int N, int chId, int mode = xydxdyMode);
+  Pads(int N, int chId, PadMode mode = 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_,
@@ -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; };
@@ -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:
@@ -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();
 };
 

Detectors/MUON/MCH/Clustering/include/MCHClustering/clusterProcessing.h

@@ -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
 {

Detectors/MUON/MCH/Clustering/src/ClusterConfig.cxx

@@ -0,0 +1,42 @@
+// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
+// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
+// All rights not expressly granted are reserved.
+//
+// This software is distributed under the terms of the GNU General Public
+// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
+//
+// In applying this license CERN does not waive the privileges and immunities
+// granted to it by virtue of its status as an Intergovernmental Organization
+// or submit itself to any jurisdiction.
+
+// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
+// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
+// All rights not expressly granted are reserved.
+//
+// This software is distributed under the terms of the GNU General Public
+// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
+//
+// In applying this license CERN does not waive the privileges and immunities
+// granted to it by virtue of its status as an Intergovernmental Organization
+// or submit itself to any jurisdiction.
+
+/// \file ClusterConfig.cxx
+/// \brief Clustering and fitting parameters
+/// \author Gilles Grasseau, Subatech
+
+#include "MCHClustering/ClusterConfig.h"
+
+namespace o2
+{
+namespace mch
+{
+
+ClusterConfig clusterConfig;
+
+void initClusterConfig()
+{
+  clusterConfig.minChargeOfClusterPerCathode = -1;
+}
+
+} // namespace mch
+} // end namespace o2