LArFCS_StepInfoSD.cxx

Go to the documentation of this file.

/*
  Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
*/
 
// class header
#include "LArFCS_StepInfoSD.h"
#include "G4Step.hh"
#include "G4ThreeVector.hh"
#include "G4TouchableHistory.hh"

#include <string>
#include <utility>
 
#include "CaloDetDescr/CaloDetDescrElement.h"
#include "CaloDetDescr/CaloDetDescrManager.h"
#include "GaudiKernel/MsgStream.h"
#include "LArG4Code/ILArCalculatorSvc.h"
 
LArFCS_StepInfoSD::LArFCS_StepInfoSD(G4String a_name,
                                     const FCS_Param::Config& config)
    : FCS_StepInfoSD(std::move(a_name), config),
      m_calculator(m_config.m_LArCalculator) {}
 
LArFCS_StepInfoSD::~LArFCS_StepInfoSD() {}
 
G4bool LArFCS_StepInfoSD::ProcessHits(G4Step* a_step, G4TouchableHistory*) {
  G4bool result(false);
  // If there's no energy, there's no hit.  (Aside: Isn't this energy
  // the same as the energy from the calculator?  Not necessarily.
  // The calculator may include detector effects such as
  // charge-collection which are not modeled by Geant4.)
  if (a_step->GetTotalEnergyDeposit() <= 0.) {
    return result;
  }
 
  if (m_calculator) {
    if (!m_calo_dd_man.get()) {
      getCaloDDManager();
    }
 
    const double StepLength = a_step->GetStepLength() / CLHEP::mm;
    const G4ThreeVector preStepPosition =
        a_step->GetPreStepPoint()
            ->GetPosition();  // pre step is the position we're interested in
    const G4ThreeVector postStepPosition =
        a_step->GetPostStepPoint()->GetPosition();
    std::vector<const G4Step*> steps;
    bool madeSubSteps(false);
    if (m_config.shorten_lar_step && StepLength > m_config.substpsize) {
      // create smaller substeps instead
      G4int nsub_step = (int)(StepLength / m_config.substpsize) + 1;
      G4double delta = 1. / ((double)nsub_step);
      // G4cout <<"Orig prestep: "<<preStepPosition<<std::endl;
      for (G4int i = 0; i < nsub_step; i++) {
        G4double fraction1 = ((G4double)i) * delta;
        G4double fraction2 = (((G4double)i) + 1.) * delta;
        G4ThreeVector subpoint1 =
            preStepPosition * (1 - fraction1) + postStepPosition * fraction1;
        G4ThreeVector subpoint2 =
            preStepPosition * (1 - fraction2) + postStepPosition * fraction2;
        G4StepPoint* startpoint = new G4StepPoint(*(a_step->GetPreStepPoint()));
        G4StepPoint* endpoint = new G4StepPoint(*(a_step->GetPostStepPoint()));
        startpoint->SetPosition(subpoint1);
        endpoint->SetPosition(subpoint2);
 
        G4Step* newstep = new G4Step(*a_step);
        newstep->SetPreStepPoint(startpoint);
        newstep->SetPostStepPoint(endpoint);
        newstep->SetStepLength((subpoint1 - subpoint2).mag());
        newstep->SetTotalEnergyDeposit(a_step->GetTotalEnergyDeposit() /
                                       nsub_step);
        steps.push_back(newstep);
        madeSubSteps = true;
      }
    } else {
      steps.push_back(a_step);
    }
 
    if (m_config.verboseLevel > 4) {
      G4cout << "LArFCS_StepInfoSD::ProcessHits(): original step in Volume: "
             << a_step->GetPreStepPoint()->GetPhysicalVolume()->GetName()
             << " position: " << a_step->GetPreStepPoint()->GetPosition()
             << " Length=" << a_step->GetStepLength() / CLHEP::mm
             << " E=" << a_step->GetTotalEnergyDeposit() << G4endl;
      std::vector<LArHitData> processedHits;
      if (m_calculator->Process(a_step, processedHits)) {
        if (m_config.verboseLevel > -1) {
          G4cout << "  #LArHitData=" << processedHits.size();
          for (const auto& lhd : processedHits) {
            G4cout << " ; id=" << (std::string)lhd.id << " E=" << lhd.energy
                   << G4endl;
          }
          G4cout << G4endl;
        }
      }
      G4cout << "LArFCS_StepInfoSD::ProcessHits(): #substep=" << steps.size()
             << G4endl;
    }
 
    double et_all = 0;  // Total collected charge in all substeps
    for (const G4Step* substep : steps) {
      double et(0.);  // Total collected charge in this substep
      G4ThreeVector stepPosition =
          0.5 * (substep->GetPreStepPoint()->GetPosition() +
                 substep->GetPostStepPoint()->GetPosition());
      std::vector<LArHitData> processedHits;
      if (m_config.verboseLevel > 4) {
        G4cout << "  LArFCS_StepInfoSD::ProcessHits(): substep in Volume: "
               << substep->GetPreStepPoint()->GetPhysicalVolume()->GetName()
               << " position: " << substep->GetPreStepPoint()->GetPosition()
               << " Length=" << substep->GetStepLength() / CLHEP::mm
               << " E=" << substep->GetTotalEnergyDeposit() << G4endl;
      }
      if (m_calculator->Process(substep, processedHits)) {
        for (const auto& larhit : processedHits) {
          et += larhit.energy;
          et_all += larhit.energy;
        }
        if (m_config.verboseLevel > 4) {
          G4cout << "    substep #LArHitData=" << processedHits.size();
          for (const auto& lhd : processedHits) {
            G4cout << " ; id=" << (std::string)lhd.id << " E=" << lhd.energy
                   << G4endl;
          }
          G4cout << G4endl;
        }
      } else {
        if (m_config.verboseLevel > 4) {
          // Maybe 0 hits or something like that...
          G4cout << this->GetName()
                 << " WARNING ProcessHits: Call to ILArCalculatorSvc::Process "
                    "failed! Details:"
                 << G4endl << "          " << "Volume: "
                 << a_step->GetPreStepPoint()->GetPhysicalVolume()->GetName()
                 << " " << m_calculator << " position: " << stepPosition
                 << " SL: " << StepLength << G4endl << "          "
                 << "Orig position: "
                 << substep->GetPreStepPoint()->GetPosition() << "  /  "
                 << substep->GetPostStepPoint()->GetPosition() << G4endl
                 << "          " << "StepLength: " << StepLength
                 << " step: " << preStepPosition << " / " << postStepPosition
                 << G4endl;
        }
        continue;
      }
 
      // drop hits with zero deposited energy (could still happen with negative
      // corrections from calculator)
      // Or if total energy is <0
      if (et <= 0.) {
        if (m_config.verboseLevel > 4) {
          G4cout << this->GetName()
                 << " WARNING ProcessHits: Total negative energy: " << et
                 << " not processing..." << G4endl;
        }
        continue;
      }
 
      const size_t numberOfProcessedHits = processedHits.size();
      const G4ThreeVector originalStepPosition = stepPosition;
      double maxSubHitEnergy = -999.;
      int maxSubHitEnergyindex = -1;
      if (numberOfProcessedHits > 0) {
        maxSubHitEnergy = processedHits[0].energy;
        maxSubHitEnergyindex = 0;
      }
      // Figure out the subhit with most energy
      for (size_t i = 1; i < numberOfProcessedHits; ++i) {
        if (maxSubHitEnergy < processedHits[i].energy) {
          maxSubHitEnergy = processedHits[i].energy;
          maxSubHitEnergyindex = i;
        }
      }
      if (maxSubHitEnergyindex ==
          -1) {  // because there were no hits; numberOfProcessedHits ==0
        G4cout << this->GetName()
               << " WARNING ProcessHits: numberOfProcessedHits is zero"
               << G4endl;
        continue;
      }
      // Identifier for the subhit with max energy
      Identifier maxEnergyIdentifier =
          this->ConvertID(processedHits[maxSubHitEnergyindex].id);
      const CaloDetDescrElement* maxEnergyCell =
          m_calo_dd_man.get()->get_element(maxEnergyIdentifier);
 
      Identifier invalidIdentifier;
      // for (size_t i=0; i<numberOfProcessedHits; ++i) {
      for (const auto& larhit : processedHits) {
        Identifier id = this->ConvertID(larhit.id);
        if (larhit.id[0] == 10) {
          if (m_config.verboseLevel > 9) {
            G4cout << this->GetName()
                   << " VERBOSE ProcessHits: Dead Material LArG4Identifier: "
                   << (std::string)larhit.id << G4endl << "          " << id
                   << G4endl << "          " << id.getString() << G4endl
                   << "          "
                   << a_step->GetPreStepPoint()
                          ->GetPhysicalVolume()
                          ->GetLogicalVolume()
                          ->GetName()
                   << G4endl << "          numberOfProcessedHits: "
                   << numberOfProcessedHits << G4endl;
            G4cout << "          " << invalidIdentifier << G4endl;
          }
        } else if (id == invalidIdentifier) {
          G4cout
              << this->GetName()
              << " WARNING ProcessHits: Something wrong with LArG4Identifier: "
              << (std::string)larhit.id << G4endl << "          " << id
              << G4endl << "          " << id.getString() << G4endl
              << "          "
              << a_step->GetPreStepPoint()
                     ->GetPhysicalVolume()
                     ->GetLogicalVolume()
                     ->GetName()
              << G4endl
              << "          numberOfProcessedHits: " << numberOfProcessedHits
              << G4endl;
          G4cout << "          " << invalidIdentifier << G4endl;
        }
        // need to get the cell information
        if (numberOfProcessedHits > 1) {
          if (!m_larEmID->is_em_barrel(id)) {
            // It didn't seem to happen outside em_barrel, so flag up if it
            // does:
            G4cout << this->GetName()
                   << " WARNING ProcessHits: Outside LAr barrel, but "
                      "numberOfProcessedHits="
                   << numberOfProcessedHits
                   << ", LArG4Identifier: " << (std::string)larhit.id << G4endl;
          } else {
            if (m_config.shift_lar_subhit) {
              // find subhit with largest energy
              if (maxSubHitEnergyindex == -1) {
                G4cout
                    << this->GetName()
                    << " WARNING ProcessHits: no subhit index with e>-999??? "
                    << G4endl;
                continue;
              }
              if (m_config.verboseLevel > 9) {
                G4cout << this->GetName()
                       << " VERBOSE ProcessHits: shifting subhits: largest "
                          "energy subhit index is "
                       << maxSubHitEnergyindex << " E: " << maxSubHitEnergy
                       << " identifier: " << maxEnergyIdentifier.getString()
                       << G4endl;
              }
              // from identifier
              const CaloDetDescrElement* thiscell =
                  m_calo_dd_man.get()->get_element(id);
              if (!maxEnergyCell) {
                // How often does this happen? Do not shift.
                G4cout << this->GetName()
                       << " WARNING ProcessHits: maxEnergyCell failed: "
                       << maxEnergyIdentifier.getString() << G4endl
                       << "          "
                       << m_calo_dd_man.get()->get_element(id)->getSampling()
                       << G4endl << "          " << originalStepPosition.eta()
                       << G4endl << "          " << originalStepPosition.phi()
                       << G4endl;
                stepPosition = originalStepPosition;
              } else if (maxEnergyCell == thiscell) {
                // The cells match, so do not shift this hit.
                if (m_config.verboseLevel > 9) {
                  G4cout << this->GetName()
                         << " VERBOSE ProcessHits: Original step position: "
                         << originalStepPosition.x() << " "
                         << originalStepPosition.y() << " "
                         << originalStepPosition.z() << G4endl << "          "
                         << "This cell: " << thiscell->x() << " "
                         << thiscell->y() << " " << thiscell->z() << G4endl
                         << "          " << "No shift" << G4endl;
                }
                stepPosition = originalStepPosition;
              } else {
                // the two cells do not match => shift
                G4ThreeVector diff(thiscell->x() - maxEnergyCell->x(),
                                   thiscell->y() - maxEnergyCell->y(),
                                   thiscell->z() - maxEnergyCell->z());
                stepPosition = originalStepPosition + diff;
                if (m_config.verboseLevel > 9) {
                  const CaloDetDescrElement* bestcell =
                      m_calo_dd_man.get()->get_element(
                          m_calo_dd_man.get()->get_element(id)->getSampling(),
                          originalStepPosition.eta(),
                          originalStepPosition.phi());
                  G4cout << this->GetName()
                         << " VERBOSE ProcessHits: Original step position: "
                         << originalStepPosition.x() << " "
                         << originalStepPosition.y() << " "
                         << originalStepPosition.z() << G4endl << "          "
                         << "This cell: " << thiscell->x() << " "
                         << thiscell->y() << " " << thiscell->z() << G4endl
                         << "          "
                         << "Highest E cell: " << maxEnergyCell->x() << " "
                         << maxEnergyCell->y() << " " << maxEnergyCell->z()
                         << G4endl << "          "
                         << "(Best cell: " << bestcell->x() << " "
                         << bestcell->y() << " " << bestcell->z() << ")"
                         << G4endl << "          "
                         << "Shifted step position: " << stepPosition.x() << " "
                         << stepPosition.y() << " " << stepPosition.z()
                         << G4endl;
                }
              }
            }
          }
        }
        // Finalize time for LAr hits?: NO
        // double time = larhit.energy)==0 ? 0. : (double)
        // larhit.time/larhit.energy/CLHEP::ns;
        double time = larhit.time;
        double energy = larhit.energy / CLHEP::MeV;
        // Drop any hits that don't have a good identifier attached
        if (!m_calo_dd_man.get()->get_element(id)) {
          if (m_config.verboseLevel > 4) {
            G4cout << this->GetName()
                   << " DEBUG update_map: bad identifier: " << id.getString()
                   << " skipping this hit." << G4endl;
          }
          continue;
        }
        this->update_map(
            stepPosition, id, energy, time, true,
            numberOfProcessedHits);  // store numberOfProcessedHits as info
      }  // numberOfProcessedHits
    }  // istep
    if (madeSubSteps) {
      // only delete steps when doing substeps. Do not delete the original
      // G4Step!
      while (!steps.empty()) {
        delete steps.back();
        steps.pop_back();
      }
    }
    if (m_config.verboseLevel > 4) {
      G4cout << "LArFCS_StepInfoSD::ProcessHits(): Etotal substeps=" << et_all
             << G4endl << G4endl << G4endl;
    }
  }
 
  return result;
}
 
Identifier LArFCS_StepInfoSD::ConvertID(const LArG4Identifier& a_ident) const {
  Identifier id;
  if (a_ident[0] == 4) {
    // is LAr
    if (a_ident[1] == 1) {
      // is LAr EM
      try {
        id = m_larEmID->channel_id(a_ident[2],   // barrel_ec
                                   a_ident[3],   // sampling
                                   a_ident[4],   // region
                                   a_ident[5],   // eta
                                   a_ident[6]);  // phi
      } catch (LArID_Exception& e) {
        G4ExceptionDescription description;
        description << "ConvertID: LArEM_ID error code " << e.code() << " "
                    << (std::string)e;
        G4Exception("LArFCS_StepInfoSD", "ConvertIDEM", FatalException,
                    description);
        abort();
      }
    } else if (a_ident[1] == 2) {
      // is EM HEC
      try {
        id = m_larHecID->channel_id(a_ident[2],   // zSide
                                    a_ident[3],   // sampling
                                    a_ident[4],   // region
                                    a_ident[5],   // eta
                                    a_ident[6]);  // phi
      } catch (LArID_Exception& e) {
        G4ExceptionDescription description;
        description << "ConvertID: LArHEC_ID error code " << e.code() << " "
                    << (std::string)e;
        G4Exception("LArFCS_StepInfoSD", "ConvertIDHEC", FatalException,
                    description);
        abort();
      }
    } else if (a_ident[1] == 3) {
      // FCAL
      if (a_ident[3] > 0) {
        // is EM FCAL
        try {
          id = m_larFcalID->channel_id(a_ident[2],   // zSide
                                       a_ident[3],   // sampling
                                       a_ident[4],   // eta
                                       a_ident[5]);  // phi
        } catch (LArID_Exception& e) {
          G4ExceptionDescription description;
          description << "ConvertID: LArFCAL_ID error code " << e.code() << " "
                      << (std::string)e;
          G4Exception("LArFCS_StepInfoSD", "ConvertIDFCAL", FatalException,
                      description);
          abort();
        }
      } else {
        G4ExceptionDescription description;
        description << "ConvertID: Unsupported ID.  ";
        G4Exception("LArFCS_StepInfoSD", "ConvertIDMiniFCAL", FatalException,
                    description);
        abort();
      }
    }
  }
  return id;
}