Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // 26 /// \file Par01/src/Par01EMShowerModel.cc << 27 /// \brief Implementation of the Par01EMShower << 28 // << 29 // 26 // >> 27 // $Id: Par01EMShowerModel.cc 77940 2013-11-29 15:15:17Z gcosmo $ 30 // 28 // 31 #include "Par01EMShowerModel.hh" 29 #include "Par01EMShowerModel.hh" 32 << 33 #include "Par01EnergySpot.hh" 30 #include "Par01EnergySpot.hh" 34 31 >> 32 #include "Randomize.hh" >> 33 35 #include "G4Electron.hh" 34 #include "G4Electron.hh" 36 #include "G4Gamma.hh" << 37 #include "G4NistManager.hh" << 38 #include "G4PhysicalConstants.hh" << 39 #include "G4Positron.hh" 35 #include "G4Positron.hh" 40 #include "G4SystemOfUnits.hh" << 36 #include "G4Gamma.hh" 41 #include "G4TouchableHandle.hh" << 42 #include "G4TransportationManager.hh" 37 #include "G4TransportationManager.hh" 43 #include "G4VSensitiveDetector.hh" 38 #include "G4VSensitiveDetector.hh" 44 #include "Randomize.hh" << 39 #include "G4TouchableHandle.hh" 45 << 40 #include "G4PhysicalConstants.hh" 46 //....oooOO0OOooo........oooOO0OOooo........oo << 41 #include "G4SystemOfUnits.hh" >> 42 #include "G4NistManager.hh" 47 43 48 Par01EMShowerModel::Par01EMShowerModel(G4Strin 44 Par01EMShowerModel::Par01EMShowerModel(G4String modelName, G4Region* envelope) 49 : G4VFastSimulationModel(modelName, envelope << 45 : G4VFastSimulationModel(modelName, envelope) 50 { 46 { 51 fFakeStep = new G4Step(); << 47 fFakeStep = new G4Step(); 52 fFakePreStepPoint = fFakeStep->GetPreStepPoi << 48 fFakePreStepPoint = fFakeStep->GetPreStepPoint(); 53 fFakePostStepPoint = fFakeStep->GetPostStepP 49 fFakePostStepPoint = fFakeStep->GetPostStepPoint(); 54 fTouchableHandle = new G4TouchableHistory(); << 50 fTouchableHandle = new G4TouchableHistory(); 55 fpNavigator = new G4Navigator(); << 51 fpNavigator = new G4Navigator(); 56 fNaviSetup = false; << 52 fNaviSetup = false; 57 fCsI = nullptr; << 53 fCsI = 0; 58 } 54 } 59 55 60 //....oooOO0OOooo........oooOO0OOooo........oo << 56 Par01EMShowerModel::Par01EMShowerModel(G4String modelName) 61 << 57 : G4VFastSimulationModel(modelName) 62 Par01EMShowerModel::Par01EMShowerModel(G4Strin << 63 { 58 { 64 fFakeStep = new G4Step(); << 59 fFakeStep = new G4Step(); 65 fFakePreStepPoint = fFakeStep->GetPreStepPoi << 60 fFakePreStepPoint = fFakeStep->GetPreStepPoint(); 66 fFakePostStepPoint = fFakeStep->GetPostStepP 61 fFakePostStepPoint = fFakeStep->GetPostStepPoint(); 67 fTouchableHandle = new G4TouchableHistory(); << 62 fTouchableHandle = new G4TouchableHistory(); 68 fpNavigator = new G4Navigator(); << 63 fpNavigator = new G4Navigator(); 69 fNaviSetup = false; << 64 fNaviSetup = false; 70 fCsI = nullptr; << 65 fCsI = 0; 71 } 66 } 72 67 73 //....oooOO0OOooo........oooOO0OOooo........oo << 74 << 75 Par01EMShowerModel::~Par01EMShowerModel() 68 Par01EMShowerModel::~Par01EMShowerModel() 76 { 69 { 77 delete fFakeStep; 70 delete fFakeStep; 78 delete fpNavigator; 71 delete fpNavigator; 79 } 72 } 80 73 81 //....oooOO0OOooo........oooOO0OOooo........oo << 82 << 83 G4bool Par01EMShowerModel::IsApplicable(const 74 G4bool Par01EMShowerModel::IsApplicable(const G4ParticleDefinition& particleType) 84 { 75 { 85 return &particleType == G4Electron::Electron << 76 return 86 || &particleType == G4Positron::Posit << 77 &particleType == G4Electron::ElectronDefinition() || 87 || &particleType == G4Gamma::GammaDef << 78 &particleType == G4Positron::PositronDefinition() || >> 79 &particleType == G4Gamma::GammaDefinition(); 88 } 80 } 89 81 90 //....oooOO0OOooo........oooOO0OOooo........oo << 91 << 92 G4bool Par01EMShowerModel::ModelTrigger(const 82 G4bool Par01EMShowerModel::ModelTrigger(const G4FastTrack& fastTrack) 93 { 83 { 94 // Applies the parameterisation above 100 Me 84 // Applies the parameterisation above 100 MeV: 95 return fastTrack.GetPrimaryTrack()->GetKinet << 85 return fastTrack.GetPrimaryTrack()->GetKineticEnergy() > 100*MeV; 96 } 86 } 97 87 98 //....oooOO0OOooo........oooOO0OOooo........oo << 88 void Par01EMShowerModel::DoIt(const G4FastTrack& fastTrack, 99 << 89 G4FastStep& fastStep) 100 void Par01EMShowerModel::DoIt(const G4FastTrac << 101 { 90 { 102 // Kill the parameterised particle: 91 // Kill the parameterised particle: 103 fastStep.KillPrimaryTrack(); 92 fastStep.KillPrimaryTrack(); 104 fastStep.ProposePrimaryTrackPathLength(0.0); 93 fastStep.ProposePrimaryTrackPathLength(0.0); 105 fastStep.ProposeTotalEnergyDeposited(fastTra 94 fastStep.ProposeTotalEnergyDeposited(fastTrack.GetPrimaryTrack()->GetKineticEnergy()); 106 95 107 // split into "energy spots" energy accordin 96 // split into "energy spots" energy according to the shower shape: 108 Explode(fastTrack); 97 Explode(fastTrack); 109 << 98 110 // and put those energy spots into the cryst 99 // and put those energy spots into the crystals: 111 BuildDetectorResponse(); 100 BuildDetectorResponse(); >> 101 112 } 102 } 113 103 114 //....oooOO0OOooo........oooOO0OOooo........oo << 115 << 116 void Par01EMShowerModel::Explode(const G4FastT 104 void Par01EMShowerModel::Explode(const G4FastTrack& fastTrack) 117 { 105 { 118 //------------------------------------------ 106 //----------------------------------------------------- 119 // 107 // 120 //------------------------------------------ 108 //----------------------------------------------------- 121 109 122 // Reduced quantities: 110 // Reduced quantities: 123 // -- critical energy in CsI: 111 // -- critical energy in CsI: 124 G4double Ec = 800 * MeV / (54. + 1.2); // 5 << 112 G4double Ec = 800*MeV/(54. + 1.2); // 54 = mean Z of CsI 125 G4double Energy = fastTrack.GetPrimaryTrack( 113 G4double Energy = fastTrack.GetPrimaryTrack()->GetKineticEnergy(); 126 G4double y = Energy / Ec; << 114 G4double y = Energy/Ec; 127 115 128 // compute value of parameter "a" of longitu 116 // compute value of parameter "a" of longitudinal profile, b assumed = 0.5 129 G4double a, tmax, b(0.5), C; 117 G4double a, tmax, b(0.5), C; 130 if (fastTrack.GetPrimaryTrack()->GetDefiniti << 118 if (fastTrack.GetPrimaryTrack()->GetDefinition() == G4Gamma::GammaDefinition()) C = 0.5; 131 C = 0.5; << 119 else C = -0.5; 132 else << 133 C = -0.5; << 134 tmax = 1.0 * (std::log(y) + C); 120 tmax = 1.0 * (std::log(y) + C); 135 a = 1.0 + b * tmax; << 121 a = 1.0 + b*tmax; 136 122 137 // t : reduced quantity = z/X0: 123 // t : reduced quantity = z/X0: 138 G4double t, bt; << 124 G4double t, bt; 139 if (fCsI == nullptr) fCsI = G4NistManager::I << 125 if ( fCsI == 0 ) fCsI = G4NistManager::Instance()->FindOrBuildMaterial("G4_CESIUM_IODIDE"); 140 G4double X0 = fCsI->GetRadlen(); << 126 G4double X0 = fCsI->GetRadlen(); 141 // Moliere radius: 127 // Moliere radius: 142 G4double Es = 21 * MeV; << 128 G4double Es = 21*MeV; 143 G4double Rm = X0 * Es / Ec; << 129 G4double Rm = X0*Es/Ec; 144 130 145 // axis of the shower, in global reference f 131 // axis of the shower, in global reference frame: 146 G4ThreeVector xShower, yShower, zShower; 132 G4ThreeVector xShower, yShower, zShower; 147 zShower = fastTrack.GetPrimaryTrack()->GetMo 133 zShower = fastTrack.GetPrimaryTrack()->GetMomentumDirection(); 148 xShower = zShower.orthogonal(); 134 xShower = zShower.orthogonal(); 149 yShower = zShower.cross(xShower); 135 yShower = zShower.cross(xShower); 150 // starting point of the shower: 136 // starting point of the shower: 151 G4ThreeVector sShower = fastTrack.GetPrimary 137 G4ThreeVector sShower = fastTrack.GetPrimaryTrack()->GetPosition(); 152 138 153 // We shoot 100 spots of energy: 139 // We shoot 100 spots of energy: 154 G4int nSpots = 100; << 140 G4int nSpots = 100; 155 G4double deposit = Energy / double(nSpots); << 141 G4double deposit = Energy/double(nSpots); 156 Par01EnergySpot eSpot; 142 Par01EnergySpot eSpot; 157 eSpot.SetEnergy(deposit); 143 eSpot.SetEnergy(deposit); 158 G4ThreeVector ePoint; 144 G4ThreeVector ePoint; 159 G4double z, r, phi; 145 G4double z, r, phi; 160 146 161 feSpotList.clear(); 147 feSpotList.clear(); 162 for (int i = 0; i < nSpots; i++) { << 148 for (int i = 0; i < nSpots; i++) 163 // Longitudinal profile: << 149 { 164 // -- shoot z according to Gamma distribut << 150 // Longitudinal profile: 165 bt = G4RandGamma::shoot(a, 1.0); << 151 // -- shoot z according to Gamma distribution: 166 t = bt / b; << 152 bt = G4RandGamma::shoot(a,1.0); 167 z = t * X0; << 153 t = bt/b; 168 << 154 z = t*X0; 169 // transverse profile: << 155 170 // we set 90% of energy in one Rm, << 156 // transverse profile: 171 // the rest between 1 and 3.5 Rm: << 157 // we set 90% of energy in one Rm, 172 G4double xr = G4UniformRand(); << 158 // the rest between 1 and 3.5 Rm: 173 if (xr < 0.9) << 159 G4double xr = G4UniformRand(); 174 r = xr / 0.9 * Rm; << 160 if (xr < 0.9) r = xr/0.9*Rm; 175 else << 161 else r = ((xr - 0.9)/0.1*2.5 + 1.0)*Rm; 176 r = ((xr - 0.9) / 0.1 * 2.5 + 1.0) * Rm; << 162 phi = G4UniformRand()*twopi; 177 phi = G4UniformRand() * twopi; << 163 178 << 164 // build the position: 179 // build the position: << 165 ePoint = sShower + 180 ePoint = sShower + z * zShower + r * std:: << 166 z*zShower + 181 << 167 r*std::cos(phi)*xShower + r*std::sin(phi)*yShower; 182 // and the energy spot: << 168 183 eSpot.SetPosition(ePoint); << 169 // and the energy spot: 184 << 170 eSpot.SetPosition(ePoint); 185 // Records the eSpot: << 171 186 feSpotList.push_back(eSpot); << 172 // Records the eSpot: 187 } << 173 feSpotList.push_back(eSpot); >> 174 } 188 } 175 } 189 176 190 //....oooOO0OOooo........oooOO0OOooo........oo << 191 177 192 void Par01EMShowerModel::BuildDetectorResponse 178 void Par01EMShowerModel::BuildDetectorResponse() 193 { 179 { 194 // Does the assignation of the energy spots 180 // Does the assignation of the energy spots to the sensitive volumes: 195 for (size_t i = 0; i < feSpotList.size(); i+ << 181 for (size_t i = 0; i < feSpotList.size(); i++) 196 // Draw the energy spot: << 182 { 197 // feSpotList[i].Draw(); << 183 // Draw the energy spot: 198 // feSpotList[i].Print(); << 184 feSpotList[i].Draw(); 199 << 185 // feSpotList[i].Print(); 200 // "converts" the energy spot into the fak << 186 201 // G4Step to pass to sensitive detector: << 187 // "converts" the energy spot into the fake 202 AssignSpotAndCallHit(feSpotList[i]); << 188 // G4Step to pass to sensitive detector: 203 } << 189 AssignSpotAndCallHit(feSpotList[i]); >> 190 } 204 } 191 } 205 192 206 //....oooOO0OOooo........oooOO0OOooo........oo << 207 193 208 void Par01EMShowerModel::AssignSpotAndCallHit( << 194 void Par01EMShowerModel::AssignSpotAndCallHit(const Par01EnergySpot &eSpot) 209 { 195 { 210 // 196 // 211 // "converts" the energy spot into the fake 197 // "converts" the energy spot into the fake 212 // G4Step to pass to sensitive detector: 198 // G4Step to pass to sensitive detector: 213 // 199 // 214 FillFakeStep(eSpot); 200 FillFakeStep(eSpot); 215 201 216 // 202 // 217 // call sensitive part: taken/adapted from t 203 // call sensitive part: taken/adapted from the stepping: 218 // Send G4Step information to Hit/Dig if the 204 // Send G4Step information to Hit/Dig if the volume is sensitive 219 // 205 // 220 G4VPhysicalVolume* pCurrentVolume = fFakeSte << 206 G4VPhysicalVolume* pCurrentVolume = >> 207 fFakeStep->GetPreStepPoint()->GetPhysicalVolume(); 221 G4VSensitiveDetector* pSensitive; 208 G4VSensitiveDetector* pSensitive; 222 << 209 223 if (pCurrentVolume != nullptr) { << 210 if( pCurrentVolume != 0 ) 224 pSensitive = pCurrentVolume->GetLogicalVol << 211 { 225 if (pSensitive != nullptr) { << 212 pSensitive = pCurrentVolume->GetLogicalVolume()-> 226 pSensitive->Hit(fFakeStep); << 213 GetSensitiveDetector(); >> 214 if( pSensitive != 0 ) >> 215 { >> 216 pSensitive->Hit(fFakeStep); >> 217 } 227 } 218 } 228 } << 229 } 219 } 230 220 231 //....oooOO0OOooo........oooOO0OOooo........oo << 232 221 233 void Par01EMShowerModel::FillFakeStep(const Pa << 222 void Par01EMShowerModel::FillFakeStep(const Par01EnergySpot &eSpot) 234 { 223 { 235 //------------------------------------------ 224 //----------------------------------------------------------- 236 // find in which volume the spot is. 225 // find in which volume the spot is. 237 //------------------------------------------ 226 //----------------------------------------------------------- 238 if (!fNaviSetup) { << 227 if (!fNaviSetup) 239 fpNavigator->SetWorldVolume(G4Transportati << 228 { 240 ->GetNavigat << 229 fpNavigator-> 241 ->GetWorldVo << 230 SetWorldVolume(G4TransportationManager::GetTransportationManager()-> 242 fpNavigator->LocateGlobalPointAndUpdateTou << 231 GetNavigatorForTracking()->GetWorldVolume()); 243 eSpot.GetPosition(), G4ThreeVector(0., 0 << 232 fpNavigator-> 244 fNaviSetup = true; << 233 LocateGlobalPointAndUpdateTouchableHandle(eSpot.GetPosition(), 245 } << 234 G4ThreeVector(0.,0.,0.), 246 else { << 235 fTouchableHandle, 247 fpNavigator->LocateGlobalPointAndUpdateTou << 236 false); 248 eSpot.GetPosition(), G4ThreeVector(0., 0 << 237 fNaviSetup = true; 249 } << 238 } >> 239 else >> 240 { >> 241 fpNavigator-> >> 242 LocateGlobalPointAndUpdateTouchableHandle(eSpot.GetPosition(), >> 243 G4ThreeVector(0.,0.,0.), >> 244 fTouchableHandle); >> 245 } 250 //-------------------------------------- 246 //-------------------------------------- 251 // Fills attribute of the G4Step needed 247 // Fills attribute of the G4Step needed 252 // by our sensitive detector: 248 // by our sensitive detector: 253 //------------------------------------- 249 //------------------------------------- 254 // set touchable volume at PreStepPoint: 250 // set touchable volume at PreStepPoint: 255 fFakePreStepPoint->SetTouchableHandle(fTouch 251 fFakePreStepPoint->SetTouchableHandle(fTouchableHandle); 256 // set total energy deposit: 252 // set total energy deposit: 257 fFakeStep->SetTotalEnergyDeposit(eSpot.GetEn 253 fFakeStep->SetTotalEnergyDeposit(eSpot.GetEnergy()); 258 } 254 } 259 255