Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2026-09-16 08:27:53 

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

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!