CHG: Calculate X and Y position using energy-weighted strip positions

include/MGUIOptionsDepthCalibration2024.h

@@ -77,6 +77,8 @@ class MGUIOptionsDepthCalibration2024 : public MGUIOptions
   //! Check button if working with the Card Cage at UCSD
   TGCheckButton* m_UCSDOverride;
 
+  //! Check button if weighting X and Y by energy
+  TGCheckButton* m_WeightedXY;
 
 #ifdef ___CLING___
  public:

include/MModuleDepthCalibration2024.h

@@ -77,6 +77,10 @@ class MModuleDepthCalibration2024 : public MModule
   //! Get whether the data came from the card cage at UCSD
   bool GetUCSDOverride() const {return m_UCSDOverride;}
 
+  //! Set whether X and Y positions should be calculated by weighting strips by energy
+  void SetWeightedXY( bool WeightedXY ) {m_WeightedXY = WeightedXY;}
+  //! Get whether X and Y positions should be calculated by weighting strips by energy
+  bool GetWeightedXY() const {return m_WeightedXY;}
 
   //! Read the XML configuration
   bool ReadXmlConfiguration(MXmlNode* Node);
@@ -101,6 +105,8 @@ class MModuleDepthCalibration2024 : public MModule
   vector<double> norm_pdf(vector<double> x, double mu, double sigma);
 	//! Adds a Depth-to-CTD relation
 	bool AddDepthCTD(vector<double> Depth, vector<vector<double>> CTDArr, int DetID, unordered_map<int, vector<double>>& DepthGrid, unordered_map<int,vector<vector<double>>>& CTDMap, unordered_map<int,vector<TSpline3*>>& SplineMap, unsigned int NPoints);
+  // Calculate the Energy-weighted X and Y strip position
+  bool CalculateEnergyWeightedPosition(vector<MStripHit*> LVStrips, vector<MStripHit*> HVStrips, double& WeightedLVStripID, double& WeightedHVStripID);
   //! Determine the Grade (geometry of charge sharing) of the Hit
   int GetHitGrade(MHit* H);
   //! Load in the specified coefficients file
@@ -153,6 +159,7 @@ class MModuleDepthCalibration2024 : public MModule
 
   // boolean for use with the card cage at UCSD since it tags all events as detector 11
   bool m_UCSDOverride;
+  bool m_WeightedXY;
 
 
 

src/MGUIOptionsDepthCalibration2024.cxx

@@ -80,8 +80,13 @@ void MGUIOptionsDepthCalibration2024::Create()
 
   m_UCSDOverride = new TGCheckButton(m_OptionsFrame, "Check this box if you're using the card cage at UCSD", 1);
   m_UCSDOverride->SetOn(dynamic_cast<MModuleDepthCalibration2024*>(m_Module)->GetUCSDOverride());
+  TGLayoutHints* Label3Layout = new TGLayoutHints(kLHintsTop | kLHintsCenterX | kLHintsExpandX, 10, 10, 10, 10);
+  m_OptionsFrame->AddFrame(m_UCSDOverride, Label3Layout);
+
+  m_WeightedXY = new TGCheckButton(m_OptionsFrame, "Check this box to weight the X and Y positions by energy deposited.", 1);
+  m_WeightedXY->SetOn(dynamic_cast<MModuleDepthCalibration2024*>(m_Module)->GetWeightedXY());
   TGLayoutHints* Label4Layout = new TGLayoutHints(kLHintsTop | kLHintsCenterX | kLHintsExpandX, 10, 10, 10, 10);
-  m_OptionsFrame->AddFrame(m_UCSDOverride, Label4Layout);
+  m_OptionsFrame->AddFrame(m_WeightedXY, Label4Layout);
 
   PostCreate();
 }
@@ -128,6 +133,7 @@ bool MGUIOptionsDepthCalibration2024::OnApply()
   dynamic_cast<MModuleDepthCalibration2024*>(m_Module)->SetCoeffsFileName(m_CoeffsFileSelector->GetFileName());
   dynamic_cast<MModuleDepthCalibration2024*>(m_Module)->SetSplinesFileName(m_SplinesFileSelector->GetFileName());
   dynamic_cast<MModuleDepthCalibration2024*>(m_Module)->SetUCSDOverride(m_UCSDOverride->IsOn());
+  dynamic_cast<MModuleDepthCalibration2024*>(m_Module)->SetWeightedXY(m_WeightedXY->IsOn());
 
   return true;
 }

src/MModuleDepthCalibration2024.cxx

@@ -266,6 +266,19 @@ bool MModuleDepthCalibration2024::AnalyzeEvent(MReadOutAssembly* Event)
         int HVStripID = HVSH->GetStripID();
         int PixelCode = 10000*DetID + 100*LVStripID + HVStripID;
 
+        if (m_WeightedXY==true) {
+
+          // Try to calculate the energy-weighted strip positions. If that doesn't work, place the Hit in the middle of the dominant strips.
+          double WeightedLVStripID = 0;
+          double WeightedHVStripID = 0;
+          bool WeightedPosResult = CalculateEnergyWeightedPosition(LVStrips, HVStrips, WeightedLVStripID, WeightedHVStripID);
+          if (WeightedPosResult == true) {
+            LVStripID = WeightedLVStripID;
+            HVStripID = WeightedHVStripID;
+          }
+
+        }
+
         // TODO: Calculate X and Y positions more rigorously using charge sharing.
         // Somewhat confusing notation: HVStrips run parallel to X-axis, so we calculate X position with LVStrips.
         double Xpos = m_YPitches[DetID]*((m_NYStrips[DetID]/2.0) - ((double)LVStripID));
@@ -544,6 +557,43 @@ bool MModuleDepthCalibration2024::LoadSplinesFile(MString FName)
 
 }
 
+bool MModuleDepthCalibration2024::CalculateEnergyWeightedPosition(vector<MStripHit*> LVStrips, vector<MStripHit*> HVStrips, double& WeightedLVStripID, double& WeightedHVStripID) {
+  // Determine the weighted strip ID positions based on charge sharing.
+  WeightedLVStripID = 0;
+  WeightedHVStripID = 0;
+  double TotalLVEnergy = 0;
+  double TotalHVEnergy = 0;
+
+  // If one side or another doesn't have strips, abort.
+  if ( (LVStrips.size() == 0) || (HVStrips.size()==0) ) {
+    return false;
+  }
+
+  // Calculate the weighted values based on energy per strip.
+  for ( unsigned int i=0; i<LVStrips.size(); ++i ) {
+    double LVSHEnergy = LVStrips[i]->GetEnergy();
+    int LVSHID = LVStrips[i]->GetStripID();
+    TotalLVEnergy += LVSHEnergy;
+    WeightedLVStripID += LVSHEnergy*LVSHID;
+  }
+  for ( unsigned int i=0; i<HVStrips.size(); ++i ) {
+    double HVSHEnergy = HVStrips[i]->GetEnergy();
+    int HVSHID = HVStrips[i]->GetStripID();
+    TotalHVEnergy += HVSHEnergy;
+    WeightedHVStripID += HVSHEnergy*HVSHID;
+  }
+
+  // If one side or another doesn't have energy, abort.
+  if ( (TotalLVEnergy == 0) || (TotalHVEnergy == 0) ) {
+    return false;
+  }
+
+  WeightedLVStripID /= TotalLVEnergy;
+  WeightedHVStripID /= TotalHVEnergy;
+
+  return true;
+}
+
 int MModuleDepthCalibration2024::GetHitGrade(MHit* H){
   // Function for choosing which Depth-to-CTD relation to use for a given event.
   // At time of writing, intention is to choose a CTD based on sub-pixel region determined via charge sharing (Event "grade").
@@ -852,6 +902,11 @@ bool MModuleDepthCalibration2024::ReadXmlConfiguration(MXmlNode* Node)
       m_UCSDOverride = (bool) UCSDOverrideNode->GetValueAsInt();
   }
 
+  MXmlNode* WeightedXYNode = Node->GetNode("WeightedXY");
+  if( WeightedXYNode != NULL ){
+      m_WeightedXY = (bool) WeightedXYNode->GetValueAsInt();
+  }
+
   return true;
 }
 
@@ -866,6 +921,7 @@ MXmlNode* MModuleDepthCalibration2024::CreateXmlConfiguration()
   new MXmlNode(Node, "CoeffsFileName", m_CoeffsFile);
   new MXmlNode(Node, "SplinesFileName", m_SplinesFile);
   new MXmlNode(Node, "UCSDOverride",(unsigned int) m_UCSDOverride);  
+  new MXmlNode(Node, "WeightedXY",(unsigned int) m_WeightedXY);  
 
   return Node;
 }