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 // Author: F. Poignant, floriane.poignant@gma 26 // Author: F. Poignant, floriane.poignant@gmail.com 27 // 27 // 28 // file STCyclotronSensitiveTarget.cc 28 // file STCyclotronSensitiveTarget.cc 29 // 29 // 30 #include "STCyclotronRun.hh" 30 #include "STCyclotronRun.hh" 31 #include "STCyclotronSensitiveTarget.hh" 31 #include "STCyclotronSensitiveTarget.hh" 32 32 33 #include "G4AnalysisManager.hh" 33 #include "G4AnalysisManager.hh" 34 #include "G4RunManager.hh" 34 #include "G4RunManager.hh" 35 #include "G4HCofThisEvent.hh" 35 #include "G4HCofThisEvent.hh" 36 #include "G4UnitsTable.hh" 36 #include "G4UnitsTable.hh" 37 #include "G4Step.hh" 37 #include "G4Step.hh" 38 #include "G4SteppingManager.hh" 38 #include "G4SteppingManager.hh" 39 #include "G4ThreeVector.hh" 39 #include "G4ThreeVector.hh" 40 #include "G4SDManager.hh" 40 #include "G4SDManager.hh" 41 #include "G4ios.hh" 41 #include "G4ios.hh" 42 #include "G4SystemOfUnits.hh" 42 #include "G4SystemOfUnits.hh" 43 #include "G4ThreeVector.hh" 43 #include "G4ThreeVector.hh" 44 #include "G4Track.hh" 44 #include "G4Track.hh" 45 #include "G4ParticleDefinition.hh" 45 #include "G4ParticleDefinition.hh" 46 #include "G4DecayTable.hh" 46 #include "G4DecayTable.hh" 47 #include "G4VDecayChannel.hh" 47 #include "G4VDecayChannel.hh" 48 #include "G4TrackVector.hh" 48 #include "G4TrackVector.hh" 49 #include "G4VProcess.hh" 49 #include "G4VProcess.hh" 50 #include "G4Tubs.hh" 50 #include "G4Tubs.hh" 51 #include <map> 51 #include <map> 52 52 53 STCyclotronSensitiveTarget::STCyclotronSensiti 53 STCyclotronSensitiveTarget::STCyclotronSensitiveTarget(G4String name, 54 STCyclotronDetectorConstruc 54 STCyclotronDetectorConstruction* det) 55 : G4VSensitiveDetector(name), 55 : G4VSensitiveDetector(name), 56 fDet(det) 56 fDet(det) 57 { 57 { 58 fTempTrack = 0; 58 fTempTrack = 0; 59 fTempTrack1 = 0; 59 fTempTrack1 = 0; 60 fTempEnergy = 0.; 60 fTempEnergy = 0.; 61 fTempVector = G4ThreeVector(0.,0.,0.); 61 fTempVector = G4ThreeVector(0.,0.,0.); 62 fTrack=0; 62 fTrack=0; 63 } 63 } 64 64 65 STCyclotronSensitiveTarget::~STCyclotronSensit 65 STCyclotronSensitiveTarget::~STCyclotronSensitiveTarget() 66 { 66 { 67 delete fTrack; 67 delete fTrack; 68 } 68 } 69 69 70 G4bool STCyclotronSensitiveTarget::ProcessHits 70 G4bool STCyclotronSensitiveTarget::ProcessHits(G4Step* aStep, G4TouchableHistory*) 71 { 71 { 72 72 73 STCyclotronRun* fRun = static_cast<STCyclotr 73 STCyclotronRun* fRun = static_cast<STCyclotronRun*>(G4RunManager::GetRunManager()->GetNonConstCurrentRun()); 74 fTrack = aStep->GetTrack(); 74 fTrack = aStep->GetTrack(); 75 75 76 auto analysisManager = G4AnalysisManager::In 76 auto analysisManager = G4AnalysisManager::Instance(); 77 77 78 //------------------------------------------ 78 //---------------------------------------------- 79 // Volume info 79 // Volume info 80 //------------------------------------------ 80 //---------------------------------------------- 81 G4double targetHalfDiameter= (fDet->GetTarge 81 G4double targetHalfDiameter= (fDet->GetTargetDiameter())/2.; 82 82 83 //------------------------------------------ 83 //---------------------------------------------- 84 // Step information 84 // Step information 85 //------------------------------------------ 85 //---------------------------------------------- 86 86 87 G4double edep = aStep->GetTotalEnergyDeposit 87 G4double edep = aStep->GetTotalEnergyDeposit(); 88 G4double energy = aStep->GetPreStepPoint()-> 88 G4double energy = aStep->GetPreStepPoint()->GetKineticEnergy(); 89 G4ThreeVector momentumDirection = aStep->Get 89 G4ThreeVector momentumDirection = aStep->GetPreStepPoint()->GetMomentumDirection(); 90 G4ThreeVector vectorPosition = aStep->GetPre 90 G4ThreeVector vectorPosition = aStep->GetPreStepPoint()->GetPosition(); 91 91 92 //------------------------------------------ 92 //---------------------------------------------- 93 // Track 93 // Track 94 //------------------------------------------ 94 //---------------------------------------------- 95 95 96 G4ParticleDefinition* thePartDef = fTrack->G 96 G4ParticleDefinition* thePartDef = fTrack->GetDefinition(); 97 G4String partType= fTrack->GetDefinition()-> 97 G4String partType= fTrack->GetDefinition()->GetParticleType(); 98 G4String name = fTrack->GetDefinition()->Get 98 G4String name = fTrack->GetDefinition()->GetParticleName(); 99 G4double timeLife = fTrack->GetDefinition()- 99 G4double timeLife = fTrack->GetDefinition()->GetPDGLifeTime(); //<--- 100 const G4VProcess* process = fTrack->GetCreat 100 const G4VProcess* process = fTrack->GetCreatorProcess(); 101 101 102 //------------------------------------------ 102 //---------------------------------------------- 103 // Collect general information concerning 103 // Collect general information concerning all of the particles 104 // Collect energy deposition ; separe dec 104 // Collect energy deposition ; separe decay case to beam case 105 //------------------------------------------ 105 //---------------------------------------------- 106 106 107 fRun->EnergyDepositionTarget(edep); 107 fRun->EnergyDepositionTarget(edep); 108 108 109 109 110 //------------------------------------------ 110 //---------------------------------------------- 111 //Collect information about protons and deut 111 //Collect information about protons and deuterons 112 //------------------------------------------ 112 //---------------------------------------------- 113 113 114 if(name == "proton" || name == "deuteron") 114 if(name == "proton" || name == "deuteron") 115 { 115 { 116 116 117 if(fTrack->GetTrackID()!=fTempTrack && ( 117 if(fTrack->GetTrackID()!=fTempTrack && (momentumDirection.getZ()>0.) && 118 vectorPosition.getX()<targetHalfDiameter && 118 vectorPosition.getX()<targetHalfDiameter && 119 vectorPosition.getX()>-targetHalfDiameter & 119 vectorPosition.getX()>-targetHalfDiameter && 120 vectorPosition.getY()<targetHalfDiameter && 120 vectorPosition.getY()<targetHalfDiameter && 121 vectorPosition.getY()>-targetHalfDiameter) 121 vectorPosition.getY()>-targetHalfDiameter) 122 { 122 { 123 analysisManager->FillH2(0,vectorPosition.g 123 analysisManager->FillH2(0,vectorPosition.getX(),vectorPosition.getY()); 124 analysisManager->FillH1(0,energy); 124 analysisManager->FillH1(0,energy); 125 fRun->CountParticlesTarget(); 125 fRun->CountParticlesTarget(); 126 fTempTrack = fTrack->GetTrackID(); 126 fTempTrack = fTrack->GetTrackID(); 127 } 127 } 128 128 129 if(fTempTrack1 == 0) 129 if(fTempTrack1 == 0) 130 { 130 { 131 fTempTrack1 = fTrack->GetTrackID(); 131 fTempTrack1 = fTrack->GetTrackID(); 132 } 132 } 133 133 134 if(fTrack->GetTrackID()!=fTempTrack1 && 134 if(fTrack->GetTrackID()!=fTempTrack1 && (momentumDirection.getZ()>0.) && 135 fTempVector.getX()<targetHalfDiameter && 135 fTempVector.getX()<targetHalfDiameter && 136 fTempVector.getX()>-targetHalfDiameter && 136 fTempVector.getX()>-targetHalfDiameter && 137 fTempVector.getY()<targetHalfDiameter && 137 fTempVector.getY()<targetHalfDiameter && 138 fTempVector.getY()>-targetHalfDiameter ) 138 fTempVector.getY()>-targetHalfDiameter ) 139 { 139 { 140 140 141 analysisManager->FillH2(4,fTempVector.getX 141 analysisManager->FillH2(4,fTempVector.getX(),fTempVector.getY()); 142 analysisManager->FillH1(2,fTempEnergy); 142 analysisManager->FillH1(2,fTempEnergy); 143 fTempTrack1 = fTrack->GetTrackID(); 143 fTempTrack1 = fTrack->GetTrackID(); 144 } 144 } 145 145 146 fTempVector = aStep->GetPostStepPoint() 146 fTempVector = aStep->GetPostStepPoint()->GetPosition(); //vectorPosition; 147 fTempEnergy = aStep->GetPostStepPoint()- 147 fTempEnergy = aStep->GetPostStepPoint()->GetKineticEnergy(); //energy; 148 148 149 analysisManager->FillH2(3,vectorPosition 149 analysisManager->FillH2(3,vectorPosition.getZ(),energy); 150 150 151 } 151 } 152 152 153 //------------------------------------------ 153 //---------------------------------------------- 154 // Store ID for particles that are 154 // Store ID for particles that are 155 // not protons/electrons or deuterons 155 // not protons/electrons or deuterons 156 //------------------------------------------ 156 //---------------------------------------------- 157 157 158 if((name != "proton") && (name != "e-") && ( 158 if((name != "proton") && (name != "e-") && (name != "deuteron")) 159 { 159 { 160 fRun->StoreIsotopeID(fTrack->GetTrackID( 160 fRun->StoreIsotopeID(fTrack->GetTrackID(),name); 161 } 161 } 162 162 163 //------------------------------------------ 163 //---------------------------------------------- 164 // Collect of information for unstable isot 164 // Collect of information for unstable isotopes 165 // generated from an interaction with the t 165 // generated from an interaction with the target 166 //------------------------------------------ 166 //---------------------------------------------- 167 167 168 if (name!="deuteron") 168 if (name!="deuteron") 169 { 169 { 170 if (( partType == "nucleus") && !(thePa 170 if (( partType == "nucleus") && !(thePartDef->GetPDGStable()) && (fTrack->GetCurrentStepNumber()==1) && timeLife!=0.) 171 { 171 { 172 //G4cout << "Saving unstable particles ... 172 //G4cout << "Saving unstable particles ..." << G4endl; 173 G4int Z=thePartDef->GetAtomicNumber(); 173 G4int Z=thePartDef->GetAtomicNumber(); 174 G4int A=thePartDef->GetAtomicMass(); 174 G4int A=thePartDef->GetAtomicMass(); 175 analysisManager->FillH2(2,Z,A); 175 analysisManager->FillH2(2,Z,A); 176 176 177 //---------------------------------------- 177 //---------------------------------------------- 178 // isotopes count 178 // isotopes count 179 //---------------------------------------- 179 //---------------------------------------------- 180 180 181 fRun->PrimaryIsotopeCountTarget(name,timeL 181 fRun->PrimaryIsotopeCountTarget(name,timeLife); 182 analysisManager->FillH1(4,fTrack->GetPosit 182 analysisManager->FillH1(4,fTrack->GetPosition().getZ()); 183 //particle that created the nucleus 183 //particle that created the nucleus 184 std::map<G4int,G4String> parentID = fRun-> 184 std::map<G4int,G4String> parentID = fRun->GetIsotopeID(); 185 G4String nameParent = parentID[fTrack->Get 185 G4String nameParent = parentID[fTrack->GetParentID()]; 186 fRun->ParticleParent(name, process->GetPro 186 fRun->ParticleParent(name, process->GetProcessName()); 187 187 188 //G4cout << name << " : " << process->GetP 188 //G4cout << name << " : " << process->GetProcessName() << " with track ID " << fTrack->GetTrackID() << " and step ID " << fTrack->GetCurrentStepNumber() << G4endl; 189 189 190 } 190 } 191 } 191 } 192 192 193 //------------------------------------------ 193 //---------------------------------------------- 194 // Collect of information for stable isoto 194 // Collect of information for stable isotopes 195 // generated from an interaction with the 195 // generated from an interaction with the target 196 //------------------------------------------ 196 //---------------------------------------------- 197 197 198 if (name!="deuteron") 198 if (name!="deuteron") 199 { 199 { 200 if (( partType == "nucleus") && (thePar 200 if (( partType == "nucleus") && (thePartDef->GetPDGStable()) && (process->GetProcessName() != "RadioactiveDecay") && (fTrack->GetCurrentStepNumber()==1) ) 201 { 201 { 202 //---------------------------------------- 202 //---------------------------------------------- 203 // isotopes count 203 // isotopes count 204 //---------------------------------------- 204 //---------------------------------------------- 205 fRun->CountStableIsotopes(name); 205 fRun->CountStableIsotopes(name); 206 } 206 } 207 } 207 } 208 208 209 209 210 //------------------------------------------ 210 //---------------------------------------------- 211 // Collect unstable isotopes from decay 211 // Collect unstable isotopes from decay 212 //------------------------------------------ 212 //---------------------------------------------- 213 213 214 if (( partType == "nucleus") && !(thePartDe 214 if (( partType == "nucleus") && !(thePartDef->GetPDGStable()) && (process->GetProcessName() == "RadioactiveDecay") && (fTrack->GetCurrentStepNumber()==1) && timeLife!=0) 215 { 215 { 216 std::map<G4int,G4String>::iterator itbis 216 std::map<G4int,G4String>::iterator itbis; 217 std::map<G4int,G4String> parentID = fRun 217 std::map<G4int,G4String> parentID = fRun->GetIsotopeID(); 218 G4String nameParent = parentID[fTrack->G 218 G4String nameParent = parentID[fTrack->GetParentID()]; 219 fRun->DecayIsotopeCountTarget(name,nameP 219 fRun->DecayIsotopeCountTarget(name,nameParent,timeLife); 220 } 220 } 221 221 222 //------------------------------------------ 222 //---------------------------------------------- 223 // Collect any other particles emitted 223 // Collect any other particles emitted 224 //------------------------------------------ 224 //---------------------------------------------- 225 225 226 if((partType!="nucleus")&&(name!="proton")&& 226 if((partType!="nucleus")&&(name!="proton")&&(name!="deuteron")) 227 { 227 { 228 fRun->ParticleCountTarget(name); 228 fRun->ParticleCountTarget(name); 229 //Condition so the particle will be coun 229 //Condition so the particle will be counted for only one step 230 if((fTrack->GetCurrentStepNumber()==1)) 230 if((fTrack->GetCurrentStepNumber()==1)) 231 { 231 { 232 if(process->GetProcessName() != "Radioacti 232 if(process->GetProcessName() != "RadioactiveDecay") 233 { 233 { 234 if(name=="e+"){ 234 if(name=="e+"){ 235 analysisManager->FillH1(5,energy); 235 analysisManager->FillH1(5,energy); 236 } 236 } 237 if(name=="e-"){ 237 if(name=="e-"){ 238 analysisManager->FillH1(6,energy); 238 analysisManager->FillH1(6,energy); 239 } 239 } 240 if(name=="gamma"){ 240 if(name=="gamma"){ 241 analysisManager->FillH1(7,energy); 241 analysisManager->FillH1(7,energy); 242 } 242 } 243 if(name=="neutron"){ 243 if(name=="neutron"){ 244 analysisManager->FillH1(8,energy); 244 analysisManager->FillH1(8,energy); 245 } 245 } 246 } 246 } 247 247 248 if(process->GetProcessName() == "Radioacti 248 if(process->GetProcessName() == "RadioactiveDecay") 249 { 249 { 250 if(name=="e+"){ 250 if(name=="e+"){ 251 analysisManager->FillH1(9,energy); 251 analysisManager->FillH1(9,energy); 252 } 252 } 253 if(name=="e-"){ 253 if(name=="e-"){ 254 analysisManager->FillH1(10,energy); 254 analysisManager->FillH1(10,energy); 255 } 255 } 256 if(name=="gamma"){ 256 if(name=="gamma"){ 257 analysisManager->FillH1(11,energy); 257 analysisManager->FillH1(11,energy); 258 } 258 } 259 if(name=="neutron"){ 259 if(name=="neutron"){ 260 analysisManager->FillH1(12,energy); 260 analysisManager->FillH1(12,energy); 261 } 261 } 262 if(name=="nu_e"){ 262 if(name=="nu_e"){ 263 analysisManager->FillH1(13,energy); 263 analysisManager->FillH1(13,energy); 264 } 264 } 265 if(name=="anti_nu_e"){ 265 if(name=="anti_nu_e"){ 266 analysisManager->FillH1(14,energy); 266 analysisManager->FillH1(14,energy); 267 } 267 } 268 } 268 } 269 } 269 } 270 } 270 } 271 271 272 272 273 fRun->SetTargetVolume(fDet->GetTargetVolume( 273 fRun->SetTargetVolume(fDet->GetTargetVolume()); 274 fRun->SetTargetThickness(fDet->GetTargetThic 274 fRun->SetTargetThickness(fDet->GetTargetThickness()); 275 fRun->SetTargetDiameter(fDet->GetTargetDiame 275 fRun->SetTargetDiameter(fDet->GetTargetDiameter()); 276 276 277 return true; 277 return true; 278 278 279 } 279 } 280 280 281 281 282 282