GGSUniqueTouchableIDComputer.cpp

/*
 * GGSUniqueTouchableIDComputer.cpp
 *
 *  Created on: 3 Jun 2020
 *      Author: Nicola Mori
 */

#include "montecarlo/scoring/GGSUniqueTouchableIDComputer.h"

#include "G4Step.hh"

RegisterTIDC(GGSUniqueTouchableIDComputer);

int GGSUniqueTouchableIDComputer::ComputeTouchableID(G4Step *aStep) {

  G4TouchableHandle touchable = aStep->GetPreStepPoint()->GetTouchableHandle();
  G4ThreeVector volPos = touchable->GetTranslation(); // Absolute volume translation

  // Build the absolute touchable ID
  // Definition: the multiplicity of a physical volume is defined as the number of that volumes inside a given mother
  // logical volume. Definition: given a history, the multiplicity at depth i is the multiplicity of the volume in
  // history at depth i Definition: given a history, the copy number at depth i is the copy number of the volume in
  // history at depth i Being N_i and M_i the copy number and multiplicity at depth i, the absolute ID of the touchable
  // at history depth 0 is defined as:
  //   ID = N_0 + M_0*N_1 + M_0*M_1*N_2 + ...
  // At each step, the product M_0*M_1*N_2*... is stored in the variable totMult.

  int historyDepth = touchable->GetHistoryDepth();
  int totMult = 1;
  int touchableID = 0;
  // Map to store the multiplicity of each physical volume

  for (int iDepth = 0; iDepth < historyDepth; iDepth++) { // This loop will go up in history up to the daughter of the
                                                          // world volume (i.e. iDepth == historyDepth - 1)
    touchableID += totMult * touchable->GetCopyNumber(iDepth);
    // Compute total multiplicity for next step
    if (iDepth != historyDepth - 1) { // No need to compute the total multiplicity on last step
      int currVolMult = 0;
      if (!(touchable->GetVolume(iDepth)->IsReplicated())) {
        // See if the multiplicity for the current non-replicated volume has been computed before
        auto storedMult = _multMap.find(touchable->GetVolume(iDepth));
        if (storedMult == _multMap.end()) {
          // Compute the multiplicity for non-replicated volumes placed inside the current mother logical
          G4LogicalVolume *currLogVol = touchable->GetVolume(iDepth)->GetLogicalVolume();
          G4LogicalVolume *motherLogVol = touchable->GetVolume(iDepth)->GetMotherLogical();
          int nDaughters = motherLogVol->GetNoDaughters();
          for (int iDaughter = 0; iDaughter < nDaughters; iDaughter++) {
            if (motherLogVol->GetDaughter(iDaughter)->GetLogicalVolume() == currLogVol)
              currVolMult++;
          }
          // Store the computed multiplicity
          _multMap[touchable->GetVolume(iDepth)] = currVolMult;
        } else {
          // Use the stored multiplicity
          currVolMult = storedMult->second;
        }
      } else {
        // Use the multiplicity of the replicated volume
        EAxis axis;
        G4int nReplicas;
        G4double width, offset;
        G4bool consuming;
        touchable->GetVolume(iDepth)->GetReplicationData(axis, nReplicas, width, offset, consuming);
        currVolMult = nReplicas;
      }
      totMult *= currVolMult;
    }
  }

  return touchableID;
}