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 SteppingAction.cc 26 /// \file SteppingAction.cc 27 /// \brief Implementation of the SteppingActio 27 /// \brief Implementation of the SteppingAction class 28 // 28 // 29 // << 29 // $Id: SteppingAction.cc 98748 2016-08-09 13:42:11Z gcosmo $ >> 30 // 30 //....oooOO0OOooo........oooOO0OOooo........oo 31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 31 //....oooOO0OOooo........oooOO0OOooo........oo 32 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 32 33 33 #include "SteppingAction.hh" 34 #include "SteppingAction.hh" 34 << 35 #include "HistoManager.hh" << 36 #include "Run.hh" 35 #include "Run.hh" >> 36 #include "HistoManager.hh" 37 37 38 #include "G4HadronicProcess.hh" << 39 #include "G4ParticleTypes.hh" 38 #include "G4ParticleTypes.hh" 40 #include "G4RunManager.hh" 39 #include "G4RunManager.hh" >> 40 #include "G4HadronicProcess.hh" >> 41 >> 42 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... >> 43 >> 44 SteppingAction::SteppingAction() >> 45 : G4UserSteppingAction() >> 46 { } 41 47 42 //....oooOO0OOooo........oooOO0OOooo........oo 48 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 43 49 44 // SteppingAction::SteppingAction() << 50 SteppingAction::~SteppingAction() 45 //: G4UserSteppingAction() << 51 { } 46 //{ } << 47 52 48 //....oooOO0OOooo........oooOO0OOooo........oo 53 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 49 54 50 void SteppingAction::UserSteppingAction(const 55 void SteppingAction::UserSteppingAction(const G4Step* aStep) 51 { 56 { 52 // check trackID and stepNumber << 57 Run* run 53 G4int trackID = aStep->GetTrack()->GetTrackI << 58 = static_cast<Run*>(G4RunManager::GetRunManager()->GetNonConstCurrentRun()); 54 G4int stepNb = aStep->GetTrack()->GetCurrent << 59 55 if (trackID * stepNb != 1) return; << 56 // ok, we are at first interaction of the pr << 57 << 58 Run* run = static_cast<Run*>(G4RunManager::G << 59 << 60 // count processes 60 // count processes 61 // << 61 // 62 const G4StepPoint* endPoint = aStep->GetPost 62 const G4StepPoint* endPoint = aStep->GetPostStepPoint(); 63 G4VProcess* process = const_cast<G4VProcess* << 63 G4VProcess* process = >> 64 const_cast<G4VProcess*>(endPoint->GetProcessDefinedStep()); 64 run->CountProcesses(process); 65 run->CountProcesses(process); 65 << 66 66 // check that an real interaction occured (e 67 // check that an real interaction occured (eg. not a transportation) 67 G4StepStatus stepStatus = endPoint->GetStepS 68 G4StepStatus stepStatus = endPoint->GetStepStatus(); 68 G4bool transmit = (stepStatus == fGeomBounda << 69 G4bool transmit = (stepStatus==fGeomBoundary || stepStatus==fWorldBoundary); 69 if (transmit) return; 70 if (transmit) return; 70 << 71 71 // real processes : sum track length << 72 //real processes : sum track length 72 // 73 // 73 G4double stepLength = aStep->GetStepLength() 74 G4double stepLength = aStep->GetStepLength(); 74 run->SumTrack(stepLength); 75 run->SumTrack(stepLength); 75 << 76 76 // energy-momentum balance initialisation << 77 //energy-momentum balance initialisation 77 // 78 // 78 const G4StepPoint* prePoint = aStep->GetPreS 79 const G4StepPoint* prePoint = aStep->GetPreStepPoint(); 79 G4double Q = -prePoint->GetKineticEnergy(); << 80 G4double Q = - prePoint->GetKineticEnergy(); 80 G4ThreeVector Pbalance = -prePoint->GetMomen << 81 G4ThreeVector Pbalance = - prePoint->GetMomentum(); 81 << 82 82 // initialisation of the nuclear channel ide << 83 //initialisation of the nuclear channel identification 83 // 84 // 84 G4ParticleDefinition* particle = aStep->GetT 85 G4ParticleDefinition* particle = aStep->GetTrack()->GetDefinition(); 85 G4String partName = particle->GetParticleNam 86 G4String partName = particle->GetParticleName(); 86 G4String nuclearChannel = partName; 87 G4String nuclearChannel = partName; 87 G4HadronicProcess* hproc = dynamic_cast<G4Ha 88 G4HadronicProcess* hproc = dynamic_cast<G4HadronicProcess*>(process); 88 const G4Isotope* target = NULL; 89 const G4Isotope* target = NULL; 89 if (hproc) target = hproc->GetTargetIsotope( 90 if (hproc) target = hproc->GetTargetIsotope(); 90 G4String targetName = "XXXX"; << 91 G4String targetName = "XXXX"; 91 if (target) targetName = target->GetName(); 92 if (target) targetName = target->GetName(); 92 nuclearChannel += " + " + targetName + " --> 93 nuclearChannel += " + " + targetName + " --> "; 93 if (targetName == "XXXX") run->SetTargetXXX( 94 if (targetName == "XXXX") run->SetTargetXXX(true); 94 << 95 95 // scattered primary particle (if any) << 96 //scattered primary particle (if any) 96 // 97 // 97 G4AnalysisManager* analysis = G4AnalysisMana 98 G4AnalysisManager* analysis = G4AnalysisManager::Instance(); 98 G4int ih = 1; 99 G4int ih = 1; 99 if (aStep->GetTrack()->GetTrackStatus() == f 100 if (aStep->GetTrack()->GetTrackStatus() == fAlive) { 100 G4double energy = endPoint->GetKineticEner << 101 G4double energy = endPoint->GetKineticEnergy(); 101 analysis->FillH1(ih, energy); << 102 analysis->FillH1(ih,energy); 102 // 103 // 103 G4ThreeVector momentum = endPoint->GetMome 104 G4ThreeVector momentum = endPoint->GetMomentum(); 104 Q += energy; << 105 Q += energy; 105 Pbalance += momentum; 106 Pbalance += momentum; 106 // 107 // 107 nuclearChannel += partName + " + "; 108 nuclearChannel += partName + " + "; 108 } << 109 } 109 << 110 110 // secondaries << 111 //secondaries 111 // << 112 // 112 const std::vector<const G4Track*>* secondary << 113 const std::vector<const G4Track*>* secondary 113 for (size_t lp = 0; lp < (*secondary).size() << 114 = aStep->GetSecondaryInCurrentStep(); 114 particle = (*secondary)[lp]->GetDefinition << 115 for (size_t lp=0; lp<(*secondary).size(); lp++) { 115 G4String name = particle->GetParticleName( << 116 particle = (*secondary)[lp]->GetDefinition(); 116 G4String type = particle->GetParticleType( << 117 G4String name = particle->GetParticleName(); >> 118 G4String type = particle->GetParticleType(); 117 G4double energy = (*secondary)[lp]->GetKin 119 G4double energy = (*secondary)[lp]->GetKineticEnergy(); 118 run->ParticleCount(name, energy); << 120 run->ParticleCount(name,energy); 119 // energy spectrum << 121 //energy spectrum 120 ih = 0; << 122 ih = 0; 121 if (particle == G4Gamma::Gamma()) << 123 if (particle == G4Gamma::Gamma()) ih = 2; 122 ih = 2; << 124 else if (particle == G4Neutron::Neutron()) ih = 3; 123 else if (particle == G4Electron::Electron( << 125 else if (particle == G4Proton::Proton()) ih = 4; 124 ih = 3; << 126 else if (particle == G4Deuteron::Deuteron()) ih = 5; 125 else if (particle == G4Neutron::Neutron()) << 127 else if (particle == G4Alpha::Alpha()) ih = 6; 126 ih = 4; << 128 else if (type == "nucleus") ih = 7; 127 else if (particle == G4Proton::Proton()) << 129 else if (type == "meson") ih = 8; 128 ih = 5; << 130 else if (type == "baryon") ih = 9; 129 else if (particle == G4Deuteron::Deuteron( << 131 if (ih > 0) analysis->FillH1(ih,energy); 130 ih = 6; << 132 //atomic mass 131 else if (particle == G4Alpha::Alpha()) << 132 ih = 7; << 133 else if (type == "nucleus") << 134 ih = 8; << 135 else if (type == "meson") << 136 ih = 9; << 137 else if (type == "baryon") << 138 ih = 10; << 139 if (ih > 0) analysis->FillH1(ih, energy); << 140 // atomic mass << 141 if (type == "nucleus") { 133 if (type == "nucleus") { 142 G4int A = particle->GetAtomicMass(); 134 G4int A = particle->GetAtomicMass(); 143 analysis->FillH1(13, A); << 135 analysis->FillH1(12, A); 144 } 136 } 145 // energy-momentum balance << 137 //energy-momentum balance 146 G4ThreeVector momentum = (*secondary)[lp]- 138 G4ThreeVector momentum = (*secondary)[lp]->GetMomentum(); 147 Q += energy; << 139 Q += energy; 148 Pbalance += momentum; 140 Pbalance += momentum; 149 // count e- from internal conversion toget << 141 //particle flag 150 if (particle == G4Electron::Electron()) pa << 151 // particle flag << 152 fParticleFlag[particle]++; 142 fParticleFlag[particle]++; 153 } 143 } 154 << 144 155 // energy-momentum balance << 145 //energy-momentum balance 156 G4double Pbal = Pbalance.mag(); 146 G4double Pbal = Pbalance.mag(); 157 run->Balance(Pbal); 147 run->Balance(Pbal); >> 148 ih = 10; >> 149 analysis->FillH1(ih,Q); 158 ih = 11; 150 ih = 11; 159 analysis->FillH1(ih, Q); << 151 analysis->FillH1(ih,Pbal); 160 ih = 12; << 152 161 analysis->FillH1(ih, Pbal); << 162 << 163 // nuclear channel 153 // nuclear channel 164 const G4int kMax = 16; << 154 const G4int kMax = 16; 165 const G4String conver[] = {"0", "", "2 " << 155 const G4String conver[] = {"0","","2 ","3 ","4 ","5 ","6 ","7 ","8 ","9 ", 166 "9 ", "10 ", "11 << 156 "10 ","11 ","12 ","13 ","14 ","15 ","16 "}; 167 std::map<G4ParticleDefinition*, G4int>::iter << 157 std::map<G4ParticleDefinition*,G4int>::iterator ip; 168 for (ip = fParticleFlag.begin(); ip != fPart 158 for (ip = fParticleFlag.begin(); ip != fParticleFlag.end(); ip++) { 169 particle = ip->first; 159 particle = ip->first; 170 G4String name = particle->GetParticleName( << 160 G4String name = particle->GetParticleName(); 171 G4int nb = ip->second; 161 G4int nb = ip->second; 172 if (nb > kMax) nb = kMax; << 162 if (nb > kMax) nb = kMax; 173 G4String Nb = conver[nb]; << 163 G4String Nb = conver[nb]; 174 if (particle == G4Gamma::Gamma()) { 164 if (particle == G4Gamma::Gamma()) { 175 run->CountGamma(nb); << 165 run->CountGamma(nb); 176 Nb = "N "; << 166 Nb = "N "; 177 name = "gamma or e-"; << 167 } 178 } << 179 if (ip != fParticleFlag.begin()) nuclearCh 168 if (ip != fParticleFlag.begin()) nuclearChannel += " + "; 180 nuclearChannel += Nb + name; 169 nuclearChannel += Nb + name; 181 } 170 } 182 << 171 183 /// G4cout << "\n nuclear channel: " << nucl << 172 ///G4cout << "\n nuclear channel: " << nuclearChannel << G4endl; 184 run->CountNuclearChannel(nuclearChannel, Q); 173 run->CountNuclearChannel(nuclearChannel, Q); 185 << 174 186 fParticleFlag.clear(); 175 fParticleFlag.clear(); 187 << 176 188 // kill event after first interaction 177 // kill event after first interaction 189 // 178 // 190 G4RunManager::GetRunManager()->AbortEvent(); << 179 G4RunManager::GetRunManager()->AbortEvent(); 191 } 180 } 192 181 193 //....oooOO0OOooo........oooOO0OOooo........oo 182 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... >> 183 194 184