Add spatial resolution 1D and 2D

source/digits_hits/include/GateDistributionFile.hh

@@ -32,7 +32,8 @@ class GateDistributionFile : public GateVDistributionArray
     inline G4int GetColumnX() const {return m_column_for_X;}
     inline G4int GetColumnY() const {return m_column_for_Y;}
 
-    void Read();
+     void Read();
+     void ReadMatrix2d();
 
     virtual void DescribeMyself(size_t indent);
   private:
@@ -41,6 +42,11 @@ class GateDistributionFile : public GateVDistributionArray
     G4int m_column_for_X;
     G4int m_column_for_Y;
     GateDistributionFileMessenger* m_messenger;
+    std::map<std::pair<double, double>, double> stddevMap;
+
+    std::vector<double> xValues;
+        std::vector<double> yValues;
+
 };
 
 

source/digits_hits/include/GateDistributionFileMessenger.hh

@@ -33,6 +33,7 @@ class GateDistributionFileMessenger: public GateDistributionArrayMessenger
     G4UIcmdWithAnInteger      *setColXCmd;
     G4UIcmdWithAnInteger      *setColYCmd;
     G4UIcmdWithoutParameter   *readCmd;
+    G4UIcmdWithoutParameter   *read2DCmd;
     G4UIcmdWithoutParameter   *autoXCmd;
 };
 

source/digits_hits/include/GateDistributionMessenger.hh

@@ -30,6 +30,7 @@ class GateDistributionMessenger: public GateNamedObjectMessenger
     void SetNewValue(G4UIcommand* aCommand, G4String aString);
     void SetUnitX(const G4String& unitX);
     void SetUnitY(const G4String& unitY);
+
     inline G4String UnitCategoryX() const {return m_unitX.empty()?"":G4UIcommand::CategoryOf(m_unitX);}
     inline G4String UnitCategoryY() const {return m_unitY.empty()?"":G4UIcommand::CategoryOf(m_unitY);}
 
@@ -41,8 +42,10 @@ class GateDistributionMessenger: public GateNamedObjectMessenger
     G4UIcmdWithoutParameter        *getMaxY_Cmd ;
     G4UIcmdWithoutParameter        *getRandom_Cmd ;
     G4UIcmdWithADoubleAndUnit      *getValueCmd ;
+
     G4String	    	    	   m_unitX;
     G4String	    	    	   m_unitY;
+
 };
 
 #endif

source/digits_hits/include/GateSpatialResolution.hh

@@ -31,7 +31,10 @@ See LICENSE.md for further details
 #include "G4TouchableHistoryHandle.hh"
 #include "GateSinglesDigitizer.hh"
 
+class GateVDistribution;
+
 class GateSpatialResolution : public GateVDigitizerModule
+
 {
 public:
 
@@ -40,9 +43,15 @@ public:
 
   void Digitize() override;
 
+
+
   //! These functions return the resolution in use.
     G4double GetFWHM()   	       { return m_fwhm; }
-    G4double GetFWHMx()   	       { return m_fwhmX; }
+    GateVDistribution* GetFWHMxdistrib()    { return m_fwhmXdistrib; }
+    GateVDistribution* GetFWHMydistrib()    { return m_fwhmYdistrib; }
+    GateVDistribution* GetFWHMxydistrib2D()    { return m_fwhmXYdistrib2D; }
+
+    G4double GetFWHMx()            { return m_fwhmX; }
     G4double GetFWHMy()   	       { return m_fwhmY; }
     G4double GetFWHMz()   	       { return m_fwhmZ; }
 
@@ -51,11 +60,14 @@ public:
       If you want a resolution of 10%, SetSpresolution(0.1)
     */
     void SetFWHM(G4double val)   { m_fwhm = val;  }
+    void SetFWHMxdistrib(GateVDistribution* dist)  { m_fwhmXdistrib= dist; }
+    void SetFWHMydistrib(GateVDistribution* dist)  { m_fwhmYdistrib = dist; }
+    void SetFWHMxydistrib2D(GateVDistribution* dist)  { m_fwhmXYdistrib2D= dist; }
+
     void SetFWHMx(G4double val)   { m_fwhmX = val;  }
     void SetFWHMy(G4double val)   { m_fwhmY = val;  }
     void SetFWHMz(G4double val)   { m_fwhmZ = val;  }
-
-
+    void SetSpatialResolutionParameters();
     inline void ConfineInsideOfSmallestElement(const G4bool& value) { m_IsConfined = value; };
     inline G4bool IsConfinedInsideOfSmallestElement() const  	      	{ return m_IsConfined; }
 
@@ -64,17 +76,23 @@ public:
 
     void UpdateVolumeID();
 
-
     //! Implementation of the pure virtual method declared by the base class GateClockDependent
     //! print-out the attributes specific of the blurring
     void DescribeMyself(size_t );
 
 protected:
     G4double m_fwhm;
 
+
     G4double m_fwhmX;
+
     G4double m_fwhmY;
     G4double m_fwhmZ;
+
+    GateVDistribution*  m_fwhmXdistrib;
+    GateVDistribution* m_fwhmYdistrib;
+    GateVDistribution*  m_fwhmXYdistrib2D;
+
     G4bool m_IsConfined;
     G4Navigator* m_Navigator;
     G4TouchableHistoryHandle m_Touchable;
@@ -83,16 +101,15 @@ protected:
 
 private:
 
-    G4int m_systemDepth;
-
-  GateDigi* m_outputDigi;
+   G4int m_systemDepth;
 
+  GateDigi* m_outputDigi;;
   GateSpatialResolutionMessenger *m_Messenger;
 
   GateDigiCollection*  m_OutputDigiCollection;
 
   GateSinglesDigitizer *m_digitizer;
-
+  G4bool   m_IsFirstEntrance;
   G4VoxelLimits limits;
   G4double Xmin, Xmax, Ymin, Ymax, Zmin, Zmax;
   G4AffineTransform at;

source/digits_hits/include/GateSpatialResolutionMessenger.hh

@@ -27,7 +27,6 @@ See LICENSE.md for further details
 #include "GateClockDependentMessenger.hh"
 class GateSpatialResolution;
 class G4UIcmdWithAString;
-
 class GateSpatialResolutionMessenger : public GateClockDependentMessenger
 {
 public:
@@ -41,10 +40,13 @@ public:
 private:
   GateSpatialResolution* m_SpatialResolution;
 
-  G4UIcmdWithADouble*  	spresolutionCmd;
-  G4UIcmdWithADouble*   spresolutionXCmd;
-  G4UIcmdWithADouble*   spresolutionYCmd;
-  G4UIcmdWithADouble*   spresolutionZCmd;
+  G4UIcmdWithADoubleAndUnit*  	spresolutionCmd;
+  G4UIcmdWithADoubleAndUnit*   spresolutionXCmd;
+  G4UIcmdWithADoubleAndUnit*   spresolutionYCmd;
+  G4UIcmdWithADoubleAndUnit*   spresolutionZCmd;
+  G4UIcmdWithAString   *spresolutionXdistribCmd;// Command declaration for 1D X-resolution distribution
+  G4UIcmdWithAString   *spresolutionYdistribCmd;// Command declaration for 1D Y-resolution distribution
+  G4UIcmdWithAString   *spresolutionXYdistrib2DCmd; // Command declaration for 2D XY-resolution distribution
   G4UIcmdWithABool* 	confineCmd;
 
 

source/digits_hits/include/GateVDistribution.hh

@@ -41,6 +41,8 @@ class GateVDistribution : public GateNamedObject
     virtual G4double MaxX() const=0;
     virtual G4double MaxY() const=0;
     virtual G4double Value(G4double x) const=0;
+    virtual G4double Value2D(G4double  x, G4double y) const;
+
     // Returns a random number following the current distribution
     // should be optimised according to each distrbution type
     virtual G4double ShootRandom() const=0;

source/digits_hits/include/GateVDistributionArray.hh

@@ -27,43 +27,57 @@ class GateVDistributionArray : public GateVDistribution
     inline void SetFactorX(G4double factor) {m_factorX=factor;}
     inline void SetFactorY(G4double factor) {m_factorY=factor;}
 
-    G4double Integral() const {return m_arrayRepartition.back();}
 
+    G4double Integral() const {return m_arrayRepartition.back();}
     virtual G4double MinX() const;
     virtual G4double MinY() const;
     virtual G4double MaxX() const;
     virtual G4double MaxY() const;
     virtual G4double Value(G4double x) const;
+
     virtual G4double RepartitionValue(G4double x) const;
+   G4double BilinearInterpolation(G4double x, G4double y) const ;
     // Returns a random number following the current distribution
     // should be optimised according to each distrbution type
     virtual G4double ShootRandom() const;
     size_t GetSize() const {return m_arrayX.size();}
-
     void Clear();
     void SetAutoStart(G4int start) {m_autoStart=start;}
+    virtual G4double Value2D(G4double x, G4double y) const;
   protected:
-    void InsertPoint(G4double x,G4double y);
-    void InsertPoint(G4double y);
+     void InsertPoint(G4double x,G4double y, G4double sigma );
+     void InsertPoint(G4double x, G4double y );
+    void InsertPoint(G4double y );
     //! private function
     G4int FindIdxBefore(G4double x
             ,const std::vector<G4double>& array) const;
-    void  FillRepartition();
+
+
+    void FillRepartition() ;
+
     std::vector<G4double>& GetArrayX() {return m_arrayX;}
     std::vector<G4double>& GetArrayY() {return m_arrayY;}
+
     std::vector<G4double>& GetArrayRepartition() {return m_arrayRepartition;}
   private:
     //! private members
     std::vector<G4double> m_arrayX;
     std::vector<G4double> m_arrayY;
+
     G4double m_minX;
     G4double m_minY;
     G4double m_maxX;
     G4double m_maxY;
     std::vector<G4double> m_arrayRepartition; //! repartition function calculus
     G4double m_factorX;
     G4double m_factorY;
+
     G4int m_autoStart;
+
+    std::map<std::pair<double, double>, double> stddevMap;
+    std::vector<double> xValues; // Stocker les valeurs x
+    std::vector<double> yValues; // Stocker les valeurs y
+
 };
 
 

source/digits_hits/src/GateDistributionFile.cc

@@ -20,6 +20,7 @@ GateDistributionFile::GateDistributionFile(const G4String& itsName)
   , m_FileName()
   , m_column_for_X(0)
   , m_column_for_Y(1)
+
 {
     m_messenger = new GateDistributionFileMessenger(this,itsName);
 }
@@ -30,16 +31,17 @@ GateDistributionFile::~GateDistributionFile()
 //___________________________________________________________________
 void GateDistributionFile::DescribeMyself(size_t indent)
 {
-  G4cout << GateTools::Indent(indent)
-    	 <<"File : "         << m_FileName
-         <<'{'    << m_column_for_X<<':'<<m_column_for_Y<<'}'
-	 << Gateendl;
+  //G4cout << GateTools::Indent(indent)
+    //	 <<"File : "         << m_FileName
+      //   <<'{'    << m_column_for_X<<':'<<m_column_for_Y<<'}'
+	 //<< Gateendl;
 }
 //___________________________________________________________________
 void GateDistributionFile::Read()
 {
     Clear();
     G4cout<<"OPENING FILE "<<m_FileName<< Gateendl;
+    G4cout << "This must be a 1D distribution" << Gateendl;
     std::ifstream f(m_FileName,std::ios::in);
     if (!f){
        G4cerr<<"[GateDistributionFile::Read] WARNING : File "<<m_FileName<<" can't be opened\n";
@@ -84,14 +86,92 @@ void GateDistributionFile::Read()
 	  else {
 	    ok = (sscanf(line,pattern.c_str(),addrFirst,addrSecond)==2);
 	    if (ok)
-	      InsertPoint(x,y);	    
+	      InsertPoint(x,y);
 	  }
-	  
+
       	  if (!ok){
 	    G4cerr<<"[GateDistributionFile::Read] WARNING : Line format unrecognised (expected == '" << pattern << "' )"
 	          << Gateendl<<line<< Gateendl;
 	  }
-	  
+
     }
     FillRepartition();
 }
+void GateDistributionFile::ReadMatrix2d() {
+    Clear();
+    G4cout << "OPENING FILE " << m_FileName << Gateendl;
+    G4cout << "This must be a 2D distribution" << Gateendl;
+    std::ifstream f(m_FileName);
+    if (!f.is_open()) {
+        G4cerr << "[GateDistributionFile::ReadMatrix2d] WARNING: File " << m_FileName << " can't be opened\n";
+        return;
+    }
+    std::string line;
+    std::vector<G4double> xValues;
+    bool isFirstLine = true;
+    bool is2D = true;
+
+    while (std::getline(f, line)) {
+        std::stringstream iss(line);
+
+        if (isFirstLine) {
+            // Process x values from the first line
+            G4double x;
+            size_t xCount = 0;
+
+            while (iss >> x) {
+                xValues.push_back(x);
+                ++xCount;
+                if (iss.peek() == ',') iss.ignore();
+            }
+
+            // Check if we have at least two x values to be considered 2D
+            if (xCount < 3) {
+                is2D = false;
+                break;
+            }
+
+            isFirstLine = false;
+        } else {
+            // Process y and stddev values for subsequent lines
+            G4double y;
+            iss >> y;
+            if (iss.peek() == ',') iss.ignore();
+            size_t stddevCount = 0;
+
+            // Read stddev values for each x value
+            for (size_t i = 0; i < xValues.size(); ++i) {
+                G4double stddev;
+                if (iss >> stddev) {
+                    ++stddevCount;
+                    // Insert the (x, y) -> stddev pair into the map using InsertPoint
+                    InsertPoint(xValues[i], y, stddev);
+                    if (iss.peek() == ',') iss.ignore();
+                }
+            }
+
+            // Check if the number of stddev values matches the number of x values
+            if (stddevCount != xValues.size()) {
+                is2D = false;
+                break;
+            }
+        }
+    }
+
+    f.close();
+    //G4cout << "Content of stddevMap:\n";
+  	    //for (const auto& entry : stddevMap) {
+  	       // std::pair<double, double> coordinates = entry.first;
+  	       // G4double stddev = entry.second;
+  	        //std::cout << "Coordinates: (" << coordinates.first << ", " << coordinates.second << "), stddev: " << stddev << std::endl;
+  	    //}
+
+    if (!is2D) {
+        G4cerr << "[GateDistributionFile::ReadMatrix2d] ERROR: File " << m_FileName << " is not a valid 2D distribution file\n";
+        Clear();
+        return;
+    }
+
+	    FillRepartition();
+	}
+

source/digits_hits/src/GateDistributionFileMessenger.cc

@@ -45,7 +45,10 @@ GateDistributionFileMessenger::GateDistributionFileMessenger(GateDistributionFil
   guidance = "Do read the file";
   readCmd = new G4UIcmdWithoutParameter(cmdName,this);
   readCmd->SetGuidance(guidance);
-
+  cmdName = GetDirectoryName()+"readMatrix2d";
+  guidance = "Do read the file";
+  read2DCmd = new G4UIcmdWithoutParameter(cmdName,this);
+  read2DCmd->SetGuidance(guidance);
   cmdName = GetDirectoryName()+"autoX";
   guidance = "Set automatic mode for X";
   autoXCmd = new G4UIcmdWithoutParameter(cmdName,this);
@@ -58,6 +61,7 @@ GateDistributionFileMessenger::GateDistributionFileMessenger(GateDistributionFil
 GateDistributionFileMessenger::~GateDistributionFileMessenger()
 {
     delete readCmd;
+    delete read2DCmd;
     delete autoXCmd;
     delete setColYCmd;
     delete setColXCmd;
@@ -83,6 +87,9 @@ void GateDistributionFileMessenger::SetNewValue(G4UIcommand* command,G4String ne
   }   else if( command==readCmd ) {
     GetDistributionFile()->Read();
   }
+ else if( command==read2DCmd ){
+	  GetDistributionFile()->ReadMatrix2d();
+  }
   else
     GateDistributionArrayMessenger::SetNewValue(command,newValue);
 }