CaTS/src/RadiatorSD.cc

0001 // ********************************************************************
0002 // * License and Disclaimer                                           *
0003 // *                                                                  *
0004 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0005 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0006 // * conditions of the Geant4 Software License,  included in the file *
0007 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0008 // * include a list of copyright holders.                             *
0009 // *                                                                  *
0010 // * Neither the authors of this software system, nor their employing *
0011 // * institutes,nor the agencies providing financial support for this *
0012 // * work  make  any representation or  warranty, express or implied, *
0013 // * regarding  this  software system or assume any liability for its *
0014 // * use.  Please see the license in the file  LICENSE  and URL above *
0015 // * for the full disclaimer and the limitation of liability.         *
0016 // *                                                                  *
0017 // * This  code  implementation is the result of  the  scientific and *
0018 // * technical work of the GEANT4 collaboration.                      *
0019 // * By using,  copying,  modifying or  distributing the software (or *
0020 // * any work based  on the software)  you  agree  to acknowledge its *
0021 // * use  in  resulting  scientific  publications,  and indicate your *
0022 // * acceptance of all terms of the Geant4 Software license.          *
0023 // ********************************************************************
0024 //
0025 // ********************************************************************
0026 //
0027 //  CaTS (Calorimetry and Tracking Simulation)
0028 //
0029 //  Authors : Hans Wenzel
0030 //            Soon Yung Jun
0031 //            (Fermi National Accelerator Laboratory)
0032 //
0033 // History
0034 //   October 18th, 2021 : first implementation
0035 //
0036 // ********************************************************************
0037 //
0038 /// \file RadiatorSD.cc
0039 /// \brief Implementation of the CaTS::RadiatorSD class
0040
0041 // Geant4 headers
0042 #include "G4Step.hh"
0043 #include "G4Track.hh"
0044 #ifdef WITH_G4OPTICKS
0045 #  include "G4ThreeVector.hh"
0046 #  include "G4ios.hh"
0047 #  include "G4UnitsTable.hh"
0048 #  include "G4SystemOfUnits.hh"
0049 #  include "G4VProcess.hh"
0050 #  include "G4VRestDiscreteProcess.hh"
0051 #  include "G4SDManager.hh"
0052 #  include "G4HCofThisEvent.hh"
0053 #  include "G4Opticks.hh"
0054 #  include "TrackInfo.hh"
0055 #  include "OpticksGenstep.h"
0056 #  include "OpticksFlags.hh"
0057 #  include "G4OpticksHit.hh"
0058 #  include "G4EventManager.hh"
0059 #  include "G4Event.hh"
0060 #  include "G4RunManager.hh"
0061 #  include "G4Version.hh"
0062 #  include "PhotonSD.hh"
0063 #  include "G4Cerenkov.hh"
0064 #  include "G4Scintillation.hh"
0065 #  include <string>
0066 #endif
0067 // project headers
0068 #include "RadiatorSD.hh"
0069 #include "ConfigurationManager.hh"
0070
0071 RadiatorSD::RadiatorSD(G4String name)
0072   : G4VSensitiveDetector(name)
0073 {
0074   verbose = ConfigurationManager::getInstance()->isEnable_verbose();
0075 }
0076
0077 void RadiatorSD::Initialize(G4HCofThisEvent*) {}
0078
0079 G4bool RadiatorSD::ProcessHits(G4Step* aStep, G4TouchableHistory*)
0080 {
0081   G4double edep = aStep->GetTotalEnergyDeposit();
0082   if(edep == 0.)
0083     return false;
0084   // only deal with charged particles
0085   G4Track* aTrack = aStep->GetTrack();
0086   G4double charge = aTrack->GetDynamicParticle()->GetCharge();
0087   if(charge == 0)
0088     return false;
0089 #ifdef WITH_G4OPTICKS
0090   if(ConfigurationManager::getInstance()->isEnable_opticks())
0091   {
0092     G4int materialIndex = 0;
0093     if(first)
0094     {
0095       aMaterial     = aTrack->GetMaterial();
0096       materialIndex = aMaterial->GetIndex();
0097       if(verbose)
0098       {
0099         G4cout << "*******************************" << G4endl;
0100         G4cout << "RadiatorSD::ProcessHits initializing Material:  "
0101                << aMaterial->GetName() << " " << G4endl;
0102         G4cout << "RadiatorSD::ProcessHits: Name "
0103                << aStep->GetPreStepPoint()
0104                     ->GetPhysicalVolume()
0105                     ->GetLogicalVolume()
0106                     ->GetName()
0107                << G4endl;
0108       }
0109       aMaterialPropertiesTable = aMaterial->GetMaterialPropertiesTable();
0110       if(verbose)
0111       {
0112         aMaterialPropertiesTable->DumpTable();
0113       }
0114       //
0115       // properties related to Scintillation
0116       //
0117 #  if(G4VERSION_NUMBER > 1072)
0118       YieldRatio =
0119         aMaterialPropertiesTable->GetConstProperty(kSCINTILLATIONYIELD1) /
0120         aMaterialPropertiesTable->GetConstProperty(
0121           kSCINTILLATIONYIELD2);  // slowerRatio,
0122       FastTimeConstant = aMaterialPropertiesTable->GetConstProperty(
0123         kSCINTILLATIONTIMECONSTANT1);  // TimeConstant,
0124       SlowTimeConstant = aMaterialPropertiesTable->GetConstProperty(
0125         kSCINTILLATIONTIMECONSTANT2);  // slowerTimeConstant,
0126 #  else
0127       Fast_Intensity = aMaterialPropertiesTable->GetProperty(kFASTCOMPONENT);
0128       Slow_Intensity = aMaterialPropertiesTable->GetProperty(kSLOWCOMPONENT);
0129       YieldRatio     = aMaterialPropertiesTable->GetConstProperty(kYIELDRATIO);
0130 #  endif
0131       ScintillationType = Slow;
0132       //
0133       // properties related to Cerenkov
0134       //
0135       Rindex = aMaterialPropertiesTable->GetProperty("RINDEX");
0136 #  if(G4VERSION_NUMBER > 1072)
0137       Pmin = Rindex->GetMinEnergy();
0138       Pmax = Rindex->GetMaxEnergy();
0139 #  else
0140       Pmin           = Rindex->GetMinLowEdgeEnergy();
0141       Pmax           = Rindex->GetMaxLowEdgeEnergy();
0142 #  endif
0143       dp   = Pmax - Pmin;
0144       nMax = Rindex->GetMaxValue();
0145       if(verbose)
0146       {
0147         G4cout << "nMax: " << nMax << "  Pmin: " << Pmin << "  Pmax: " << Pmax
0148                << "  dp: " << dp << G4endl;
0149         Rindex->DumpValues();
0150       }
0151       //
0152       first = false;
0153     }                    // end if first
0154     G4int Sphotons = 0;  // number of scintillation photons this step
0155     G4int Cphotons = 0;  // number of Cerenkov photons this step
0156     //
0157     // info needed for generating Cerenkov photons on the GPU;
0158     //
0159     G4double maxCos                      = 0.0;
0160     G4double maxSin2                     = 0.0;
0161     G4double beta                        = 0.0;
0162     G4double beta1                       = 0.0;
0163     G4double beta2                       = 0.0;
0164     G4double BetaInverse                 = 0.0;
0165     G4double MeanNumberOfPhotons1        = 0.0;
0166     G4double MeanNumberOfPhotons2        = 0.0;
0167     G4SteppingManager* fpSteppingManager = G4EventManager::GetEventManager()
0168                                              ->GetTrackingManager()
0169                                              ->GetSteppingManager();
0170     G4StepStatus stepStatus = fpSteppingManager->GetfStepStatus();
0171     if(stepStatus != fAtRestDoItProc)
0172     {
0173       G4ProcessVector* procPost = fpSteppingManager->GetfPostStepDoItVector();
0174       size_t MAXofPostStepLoops = fpSteppingManager->GetMAXofPostStepLoops();
0175       for(size_t i3 = 0; i3 < MAXofPostStepLoops; i3++)
0176       {
0177         if((*procPost)[i3]->GetProcessName() == "Cerenkov")
0178         {
0179           G4Cerenkov* proc = (G4Cerenkov*) (*procPost)[i3];
0180           thePhysicsTable  = proc->GetPhysicsTable();
0181           CerenkovAngleIntegrals =
0182             (G4PhysicsOrderedFreeVector*) ((*thePhysicsTable)(materialIndex));
0183           Cphotons = proc->GetNumPhotons();
0184           if(Cphotons > 0)
0185           {
0186             beta1       = aStep->GetPreStepPoint()->GetBeta();
0187             beta2       = aStep->GetPostStepPoint()->GetBeta();
0188             beta        = (beta1 + beta2) * 0.5;
0189             BetaInverse = 1. / beta;
0190             maxCos      = BetaInverse / nMax;
0191             maxSin2     = (1.0 - maxCos) * (1.0 + maxCos);
0192             MeanNumberOfPhotons1 =
0193               proc->GetAverageNumberOfPhotons(charge, beta1, aMaterial, Rindex);
0194             MeanNumberOfPhotons2 =
0195               proc->GetAverageNumberOfPhotons(charge, beta2, aMaterial, Rindex);
0196           }
0197         }
0198         if((*procPost)[i3]->GetProcessName() == "Scintillation")
0199         {
0200           G4Scintillation* proc1 = (G4Scintillation*) (*procPost)[i3];
0201           Sphotons               = proc1->GetNumPhotons();
0202         }
0203       }
0204     }
0205     tSphotons += Sphotons;
0206     tCphotons += Cphotons;
0207     G4ThreeVector deltaPosition = aStep->GetDeltaPosition();
0208     G4double ScintillationTime  = 0. * ns;
0209     G4int scntId                = 1;
0210     G4StepPoint* pPreStepPoint  = aStep->GetPreStepPoint();
0211     G4ThreeVector x0            = pPreStepPoint->GetPosition();
0212     G4ThreeVector p0            = aStep->GetDeltaPosition().unit();
0213     //
0214     // harvest the Scintillation photon gensteps:
0215     //
0216     if(Sphotons > 0)
0217     {
0218       G4double ScintillationRiseTime = 0.0;
0219       G4Opticks::Get()->collectGenstep_G4Scintillation_1042(
0220         aTrack, aStep, Sphotons, scntId, ScintillationTime,
0221         ScintillationRiseTime);
0222     }
0223     //
0224     // harvest the Cerenkov photon gensteps:
0225     //
0226     if(Cphotons > 0)
0227     {
0228       G4Opticks::Get()->collectGenstep_G4Cerenkov_1042(
0229         aTrack, aStep, Cphotons, BetaInverse, Pmin, Pmax, maxCos, maxSin2,
0230         MeanNumberOfPhotons1, MeanNumberOfPhotons2);
0231     }
0232     G4Opticks* g4ok      = G4Opticks::Get();
0233     G4RunManager* rm     = G4RunManager::GetRunManager();
0234     const G4Event* event = rm->GetCurrentEvent();
0235     G4int eventid        = event->GetEventID();
0236     G4OpticksHit hit;
0237     unsigned num_photons = g4ok->getNumPhotons();
0238     if(num_photons > ConfigurationManager::getInstance()->getMaxPhotons())
0239     {
0240       g4ok->propagateOpticalPhotons(eventid);
0241       G4HCtable* hctable = G4SDManager::GetSDMpointer()->GetHCtable();
0242       for(G4int i = 0; i < hctable->entries(); ++i)
0243       {
0244         std::string sdn   = hctable->GetSDname(i);
0245         std::size_t found = sdn.find("Photondetector");
0246         if(found != std::string::npos)
0247         {
0248           PhotonSD* aSD =
0249             (PhotonSD*) G4SDManager::GetSDMpointer()->FindSensitiveDetector(
0250               sdn);
0251           aSD->AddOpticksHits();
0252         }
0253       }
0254       g4ok->reset();
0255     }
0256   }
0257 #endif
0258   return true;
0259 }
0260
0261 void RadiatorSD::EndOfEvent(G4HCofThisEvent*)
0262 {
0263   tSphotons = 0;
0264   tCphotons = 0;
0265 }