* Add new method to fill DSSSD Scorers (Rectangle only for the momen)

 Please enter the commit message for your changes. Lines starting

NPSimulation/Scorers/DSSDScorers.cc

@@ -155,7 +155,7 @@ PS_Rectangle::PS_Rectangle(G4String name, G4int Level, G4double StripPlaneLength
   m_StripPlaneWidth = StripPlaneWidth;
   m_NumberOfStripLength = NumberOfStripLength;
   m_NumberOfStripWidth = NumberOfStripWidth;
-  m_StripPitchLength = m_StripPlaneLength / m_NumberOfStripLength ;
+  m_StripPitchLength = m_StripPlaneLength / m_NumberOfStripLength;
   m_StripPitchWidth = m_StripPlaneWidth / m_NumberOfStripWidth;
   m_Level = Level;
   if (axis == "xy")
@@ -170,60 +170,179 @@ PS_Rectangle::PS_Rectangle(G4String name, G4int Level, G4double StripPlaneLength
 PS_Rectangle::~PS_Rectangle() {}
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-G4bool PS_Rectangle::ProcessHits(G4Step* aStep, G4TouchableHistory*) {
+G4bool PS_Rectangle::ProcessHitsNew(G4Step* aStep) {
   // contain Energy Time, DetNbr, StripFront and StripBack
   t_Energy = aStep->GetTotalEnergyDeposit();
   t_Time = aStep->GetPreStepPoint()->GetGlobalTime();
-
   t_DetectorNumber = aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(m_Level);
-  t_Position = aStep->GetPreStepPoint()->GetPosition();
 
+  // Transform Pre/Post-StepPoint position from world coordinate system to local coordinate system of the volume
+  t_Position = aStep->GetPreStepPoint()->GetPosition();
   t_Position =
       aStep->GetPreStepPoint()->GetTouchableHandle()->GetHistory()->GetTopTransform().TransformPoint(t_Position);
+  post_Position = aStep->GetPostStepPoint()->GetPosition();
+  post_Position =
+      aStep->GetPostStepPoint()->GetTouchableHandle()->GetHistory()->GetTopTransform().TransformPoint(post_Position);
+
+  double Win = 0, Wout = 0;
+  double Wstart = 0, Wend = 0;
+  double Lin = 0, Lout = 0;
+  double Lstart = 0, Lend = 0;
   if (m_Axis == ps_xy) {
-
-    t_StripLengthNumber = (int)((t_Position.x() + m_StripPlaneLength / 2.) / m_StripPitchLength) + 1;
-    t_StripWidthNumber = (int)((t_Position.y() + m_StripPlaneWidth / 2.) / m_StripPitchWidth) + 1;
-    }
+    Lin = t_Position.x() + m_StripPlaneLength / 2.;
+    Lout = post_Position.x() + m_StripPlaneLength / 2.;
+    Win = t_Position.y() + m_StripPlaneWidth / 2.;
+    Wout = post_Position.y() + m_StripPlaneWidth / 2.;
+  }
   else if (m_Axis == ps_yz) {
-    t_StripLengthNumber = (int)((t_Position.y() + m_StripPlaneLength / 2.) / m_StripPitchLength) + 1;
-    t_StripWidthNumber = (int)((t_Position.z() + m_StripPlaneWidth / 2.) / m_StripPitchWidth) + 1;
+    Lin = t_Position.y() + m_StripPlaneLength / 2.;
+    Lout = post_Position.y() + m_StripPlaneLength / 2.;
+    Win = t_Position.z() + m_StripPlaneWidth / 2.;
+    Wout = post_Position.z() + m_StripPlaneWidth / 2.;
   }
   else if (m_Axis == ps_xz) {
-    t_StripLengthNumber = (int)((t_Position.x() + m_StripPlaneLength / 2.) / m_StripPitchLength) + 1;
-    t_StripWidthNumber = (int)((t_Position.z() + m_StripPlaneWidth / 2.) / m_StripPitchWidth) + 1;
+    Lin = t_Position.x() + m_StripPlaneLength / 2.;
+    Lout = post_Position.x() + m_StripPlaneLength / 2.;
+    Win = t_Position.z() + m_StripPlaneWidth / 2.;
+    Wout = post_Position.z() + m_StripPlaneWidth / 2.;
   }
-  // Rare case where particle is close to edge of silicon plan
-  if (t_StripLengthNumber > m_NumberOfStripLength)
-    t_StripLengthNumber = m_NumberOfStripLength;
-  if (t_StripWidthNumber > m_NumberOfStripWidth)
-    t_StripWidthNumber = m_NumberOfStripWidth;
 
+  if (Win < Wout) {
+    Wstart = Win;
+    Wend = Wout;
+  }
+  else {
+    Wstart = Wout;
+    Wend = Win;
+  }
+  if (Lin < Lout) {
+    Lstart = Lin;
+    Lend = Lout;
+  }
+  else {
+    Lstart = Lout;
+    Lend = Lin;
+  }
 
-  // Check if the particle has interact before, if yes, add up the energies.
+  // Width
+  // cout << "--------Width--------" << endl;
+  // cout << "totalE:\t" << t_Energy << endl;
+  int Nstart = (int)(Wstart / m_StripPitchWidth);
+  int Nend = (int)(Wend / m_StripPitchWidth);
   vector<DSSDData>::iterator it;
-  // Length
-  it = m_HitLength.find(DSSDData::CalculateIndex(t_StripLengthNumber, t_DetectorNumber));
-
-  if (it != m_HitLength.end()) {
-    it->Add(t_Energy);
+  double dW = 0;
+  for (int i = Nstart; i <= Nend; i++) {
+    if (i == Nstart)
+      dW = (Nstart + 1) * m_StripPitchWidth - Wstart;
+    else if (i == Nend)
+      dW = Wend - Nend * m_StripPitchWidth;
+    else
+      dW = m_StripPitchWidth;
+    t_StripWidthNumber = i + 1;
+
+    // Rare case where particle is close to edge of silicon plan
+    if (t_StripWidthNumber > m_NumberOfStripWidth)
+      t_StripWidthNumber = m_NumberOfStripWidth;
+
+    it = m_HitWidth.find(DSSDData::CalculateIndex(t_StripWidthNumber, t_DetectorNumber));
+    if (it != m_HitWidth.end()) {
+      it->Add(t_Energy * dW / (Wend - Wstart));
+    }
+    else {
+      m_HitWidth.Set(t_Energy * dW / (Wend - Wstart), t_Time, t_StripWidthNumber, t_DetectorNumber);
+    }
   }
-  else
-  {
-    m_HitLength.Set(t_Energy, t_Time, t_StripLengthNumber, t_DetectorNumber);
+
+  // Length
+  // cout << "--------Length--------" << endl;
+  // cout << "totalE:\t" << t_Energy << endl;
+  Nstart = (int)(Lstart / m_StripPitchLength);
+  Nend = (int)(Lend / m_StripPitchLength);
+  double dL = 0;
+  for (int i = Nstart; i <= Nend; i++) {
+    if (i == Nstart)
+      dL = (Nstart + 1) * m_StripPitchLength - Lstart;
+    else if (i == Nend)
+      dL = Lend - Nend * m_StripPitchLength;
+    else
+      dL = m_StripPitchLength;
+    t_StripLengthNumber = i + 1;
+
+    // Rare case where particle is close to edge of silicon plan
+    if (t_StripLengthNumber > m_NumberOfStripLength)
+      t_StripLengthNumber = m_NumberOfStripLength;
+
+    it = m_HitLength.find(DSSDData::CalculateIndex(t_StripLengthNumber, t_DetectorNumber));
+    if (it != m_HitLength.end()) {
+      it->Add(t_Energy * dL / (Lend - Lstart));
+    }
+    else {
+      m_HitLength.Set(t_Energy * dL / (Lend - Lstart), t_Time, t_StripLengthNumber, t_DetectorNumber);
+    }
   }
 
-  it = m_HitWidth.find(DSSDData::CalculateIndex(t_StripWidthNumber, t_DetectorNumber));
+  return TRUE;
+}
 
-  if (it != m_HitWidth.end()) {
-    it->Add(t_Energy);
+G4bool PS_Rectangle::ProcessHits(G4Step* aStep, G4TouchableHistory*) {
+
+  bool UseNewMethod = false;
+
+  if (UseNewMethod) {
+    ProcessHitsNew(aStep);
   }
-  else
-    {
-    m_HitWidth.Set(t_Energy, t_Time, t_StripWidthNumber, t_DetectorNumber);
+  else {
+    t_Energy = aStep->GetTotalEnergyDeposit();
+    t_Time = aStep->GetPreStepPoint()->GetGlobalTime();
+
+    t_DetectorNumber = aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(m_Level);
+    t_Position = aStep->GetPreStepPoint()->GetPosition();
+
+    // Transform PreStepPoint position from world coordinate system to local coordinate system of the volume
+    t_Position =
+        aStep->GetPreStepPoint()->GetTouchableHandle()->GetHistory()->GetTopTransform().TransformPoint(t_Position);
+
+    if (m_Axis == ps_xy) {
+      t_StripLengthNumber = (int)((t_Position.x() + m_StripPlaneLength / 2.) / m_StripPitchLength) + 1;
+      t_StripWidthNumber = (int)((t_Position.y() + m_StripPlaneWidth / 2.) / m_StripPitchWidth) + 1;
+    }
+    else if (m_Axis == ps_yz) {
+      t_StripLengthNumber = (int)((t_Position.y() + m_StripPlaneLength / 2.) / m_StripPitchLength) + 1;
+      t_StripWidthNumber = (int)((t_Position.z() + m_StripPlaneWidth / 2.) / m_StripPitchWidth) + 1;
+    }
+    else if (m_Axis == ps_xz) {
+      t_StripLengthNumber = (int)((t_Position.x() + m_StripPlaneLength / 2.) / m_StripPitchLength) + 1;
+      t_StripWidthNumber = (int)((t_Position.z() + m_StripPlaneWidth / 2.) / m_StripPitchWidth) + 1;
+    }
+    // Rare case where particle is close to edge of silicon plan
+    if (t_StripLengthNumber > m_NumberOfStripLength)
+      t_StripLengthNumber = m_NumberOfStripLength;
+    if (t_StripWidthNumber > m_NumberOfStripWidth)
+      t_StripWidthNumber = m_NumberOfStripWidth;
+
+    // Check if the particle has interact before, if yes, add up the energies.
+    vector<DSSDData>::iterator it;
+    // Length
+    it = m_HitLength.find(DSSDData::CalculateIndex(t_StripLengthNumber, t_DetectorNumber));
+
+    if (it != m_HitLength.end()) {
+      it->Add(t_Energy);
+    }
+    else {
+      m_HitLength.Set(t_Energy, t_Time, t_StripLengthNumber, t_DetectorNumber);
+    }
+
+    it = m_HitWidth.find(DSSDData::CalculateIndex(t_StripWidthNumber, t_DetectorNumber));
+
+    if (it != m_HitWidth.end()) {
+      it->Add(t_Energy);
+    }
+    else {
+      m_HitWidth.Set(t_Energy, t_Time, t_StripWidthNumber, t_DetectorNumber);
+    }
+    //}
+    return TRUE;
   }
-  //}
-  return TRUE;
 }
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

NPSimulation/Scorers/DSSDScorers.hh

@@ -176,6 +176,7 @@ namespace DSSDScorers {
 
    protected: // with description
     G4bool ProcessHits(G4Step*, G4TouchableHistory*);
+    G4bool ProcessHitsNew(G4Step*);
 
    public:
     void Initialize(G4HCofThisEvent*);
@@ -200,6 +201,7 @@ namespace DSSDScorers {
 
    private: // Needed for intermediate calculation (avoid multiple instantiation in Processing Hit)
     G4ThreeVector t_Position;
+    G4ThreeVector post_Position;
     double t_Energy;
     double t_Time;
     unsigned int t_DetectorNumber;