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 // $Id: F04GlobalField.cc 96981 2016-05-19 09:42:40Z gcosmo $ 26 // 27 // 27 /// \file field/field04/src/F04GlobalField.cc 28 /// \file field/field04/src/F04GlobalField.cc 28 /// \brief Implementation of the F04GlobalFiel 29 /// \brief Implementation of the F04GlobalField class 29 // 30 // 30 31 31 #include "F04GlobalField.hh" << 32 #include <time.h> 32 33 33 #include "F04FocusSolenoid.hh" << 34 #include "Randomize.hh" 34 #include "F04SimpleSolenoid.hh" << 35 #include "G4TransportationManager.hh" 35 36 36 #include "G4CashKarpRKF45.hh" << 37 #include "G4ClassicalRK4.hh" << 38 #include "G4ExplicitEuler.hh" 37 #include "G4ExplicitEuler.hh" 39 #include "G4ImplicitEuler.hh" 38 #include "G4ImplicitEuler.hh" 40 #include "G4SimpleHeum.hh" << 41 #include "G4SimpleRunge.hh" 39 #include "G4SimpleRunge.hh" >> 40 #include "G4SimpleHeum.hh" >> 41 #include "G4ClassicalRK4.hh" >> 42 #include "G4CashKarpRKF45.hh" 42 #include "G4SystemOfUnits.hh" 43 #include "G4SystemOfUnits.hh" 43 #include "G4TransportationManager.hh" << 44 #include "Randomize.hh" << 45 44 46 #include <ctime> << 45 #include "F04GlobalField.hh" >> 46 #include "F04SimpleSolenoid.hh" >> 47 #include "F04FocusSolenoid.hh" 47 48 48 //....oooOO0OOooo........oooOO0OOooo........oo 49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 49 50 50 G4ThreadLocal F04GlobalField* F04GlobalField:: << 51 G4ThreadLocal F04GlobalField* F04GlobalField::fObject = 0; 51 52 52 //....oooOO0OOooo........oooOO0OOooo........oo 53 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 53 54 54 F04GlobalField::F04GlobalField(F04DetectorCons << 55 F04GlobalField::F04GlobalField(F04DetectorConstruction* det) >> 56 : G4ElectroMagneticField(), >> 57 fMinStep(0.01*mm), fDeltaChord(3.0*mm), >> 58 fDeltaOneStep(0.01*mm), fDeltaIntersection(0.1*mm), >> 59 fEpsMin(2.5e-7*mm), fEpsMax(0.05*mm), >> 60 fEquation(0), fFieldManager(0), >> 61 fFieldPropagator(0), fStepper(0), fChordFinder(0), >> 62 fDetectorConstruction(det) >> 63 //F04GlobalField::F04GlobalField(F04DetectorConstruction* det) >> 64 // : G4MagneticField(), >> 65 // fMinStep(0.01*mm), fDeltaChord(3.0*mm), >> 66 // fDeltaOneStep(0.01*mm), fDeltaIntersection(0.1*mm), >> 67 // fEpsMin(2.5e-7*mm), fEpsMax(0.05*mm), >> 68 // fEquation(0), fFieldManager(0), >> 69 // fFieldPropagator(0), fStepper(0), fChordFinder(0), >> 70 // fDetectorConstruction(det) 55 { 71 { 56 fFieldMessenger = new F04FieldMessenger(this << 72 fFieldMessenger = new F04FieldMessenger(this,det); 57 73 58 fFields = new FieldList(); 74 fFields = new FieldList(); 59 75 >> 76 fStepperType = 4 ; // ClassicalRK4 is default stepper >> 77 60 // set object 78 // set object >> 79 61 fObject = this; 80 fObject = this; >> 81 fFirst = true; >> 82 >> 83 fNfp = 0; >> 84 fFp = NULL; 62 85 63 ConstructField(); 86 ConstructField(); 64 } 87 } 65 88 66 //....oooOO0OOooo........oooOO0OOooo........oo 89 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 67 90 68 F04GlobalField::~F04GlobalField() 91 F04GlobalField::~F04GlobalField() 69 { 92 { 70 Clear(); 93 Clear(); 71 94 72 delete fFields; 95 delete fFields; 73 96 74 delete fFieldMessenger; 97 delete fFieldMessenger; 75 98 76 delete fEquation; << 99 if (fEquation) delete fEquation; 77 delete fStepper; << 100 if (fStepper) delete fStepper; 78 delete fChordFinder; << 101 if (fChordFinder) delete fChordFinder; 79 } 102 } 80 103 81 //....oooOO0OOooo........oooOO0OOooo........oo 104 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 82 105 83 void F04GlobalField::ConstructField() 106 void F04GlobalField::ConstructField() 84 { 107 { 85 Clear(); 108 Clear(); 86 109 87 // Construct equ. of motion of particles th 110 // Construct equ. of motion of particles through B fields 88 // fEquation = new G4Mag_EqRhs(this); << 111 // fEquation = new G4Mag_EqRhs(this); 89 // Construct equ. of motion of particles th 112 // Construct equ. of motion of particles through e.m. fields 90 // fEquation = new G4EqMagElectricField(thi << 113 // fEquation = new G4EqMagElectricField(this); 91 // Construct equ. of motion of particles in 114 // Construct equ. of motion of particles including spin through B fields 92 // fEquation = new G4Mag_SpinEqRhs(this); << 115 // fEquation = new G4Mag_SpinEqRhs(this); 93 // Construct equ. of motion of particles in 116 // Construct equ. of motion of particles including spin through e.m. fields 94 fEquation = new G4EqEMFieldWithSpin(this); 117 fEquation = new G4EqEMFieldWithSpin(this); 95 118 96 // Get transportation, field, and propagato 119 // Get transportation, field, and propagator managers 97 G4TransportationManager* transportManager = << 120 G4TransportationManager* transportManager = >> 121 G4TransportationManager::GetTransportationManager(); 98 122 99 fFieldManager = GetGlobalFieldManager(); 123 fFieldManager = GetGlobalFieldManager(); 100 124 101 fFieldPropagator = transportManager->GetProp 125 fFieldPropagator = transportManager->GetPropagatorInField(); 102 126 103 // Need to SetFieldChangesEnergy to account 127 // Need to SetFieldChangesEnergy to account for a time varying electric 104 // field (r.f. fields) 128 // field (r.f. fields) 105 fFieldManager->SetFieldChangesEnergy(true); 129 fFieldManager->SetFieldChangesEnergy(true); 106 130 107 // Set the field 131 // Set the field 108 fFieldManager->SetDetectorField(this); 132 fFieldManager->SetDetectorField(this); 109 133 110 // Choose a stepper for integration of the 134 // Choose a stepper for integration of the equation of motion 111 SetStepper(); 135 SetStepper(); 112 136 113 // Create a cord finder providing the (glob 137 // Create a cord finder providing the (global field, min step length, 114 // a pointer to the stepper) 138 // a pointer to the stepper) 115 fChordFinder = new G4ChordFinder((G4Magnetic << 139 fChordFinder = new G4ChordFinder((G4MagneticField*)this,fMinStep,fStepper); 116 140 117 // Set accuracy parameters 141 // Set accuracy parameters 118 fChordFinder->SetDeltaChord(fDeltaChord); << 142 fChordFinder->SetDeltaChord( fDeltaChord ); 119 143 120 fFieldManager->SetAccuraciesWithDeltaOneStep 144 fFieldManager->SetAccuraciesWithDeltaOneStep(fDeltaOneStep); 121 145 122 fFieldManager->SetDeltaIntersection(fDeltaIn 146 fFieldManager->SetDeltaIntersection(fDeltaIntersection); 123 147 124 fFieldPropagator->SetMinimumEpsilonStep(fEps 148 fFieldPropagator->SetMinimumEpsilonStep(fEpsMin); 125 fFieldPropagator->SetMaximumEpsilonStep(fEps 149 fFieldPropagator->SetMaximumEpsilonStep(fEpsMax); 126 150 127 G4cout << "Accuracy Parameters:" << 151 G4cout << "Accuracy Parameters:" << 128 << " MinStep=" << fMinStep << " Delta << 152 " MinStep=" << fMinStep << 129 << " DeltaOneStep=" << fDeltaOneStep << 153 " DeltaChord=" << fDeltaChord << 130 G4cout << " " << 154 " DeltaOneStep=" << fDeltaOneStep << G4endl; 131 << " DeltaIntersection=" << fDeltaInt << 155 G4cout << " " << 132 << " EpsMax=" << fEpsMax << G4endl; << 156 " DeltaIntersection=" << fDeltaIntersection << >> 157 " EpsMin=" << fEpsMin << >> 158 " EpsMax=" << fEpsMax << G4endl; 133 159 134 fFieldManager->SetChordFinder(fChordFinder); 160 fFieldManager->SetChordFinder(fChordFinder); 135 161 136 G4double l = 0.0; 162 G4double l = 0.0; 137 G4double B1 = fDetectorConstruction->GetCapt 163 G4double B1 = fDetectorConstruction->GetCaptureMgntB1(); 138 G4double B2 = fDetectorConstruction->GetCapt 164 G4double B2 = fDetectorConstruction->GetCaptureMgntB2(); 139 165 140 G4LogicalVolume* logicCaptureMgnt = fDetecto 166 G4LogicalVolume* logicCaptureMgnt = fDetectorConstruction->GetCaptureMgnt(); 141 G4ThreeVector captureMgntCenter = fDetectorC << 167 G4ThreeVector captureMgntCenter = >> 168 fDetectorConstruction->GetCaptureMgntCenter(); 142 169 143 auto focusSolenoid = new F04FocusSolenoid(B1 << 170 F04FocusSolenoid* focusSolenoid = 144 focusSolenoid->SetHalf(true); << 171 new F04FocusSolenoid(B1, B2, l, logicCaptureMgnt,captureMgntCenter); >> 172 focusSolenoid -> SetHalf(true); 145 173 146 G4double B = fDetectorConstruction->GetTrans 174 G4double B = fDetectorConstruction->GetTransferMgntB(); 147 175 148 G4LogicalVolume* logicTransferMgnt = fDetect << 176 G4LogicalVolume* logicTransferMgnt = 149 G4ThreeVector transferMgntCenter = fDetector << 177 fDetectorConstruction->GetTransferMgnt(); >> 178 G4ThreeVector transferMgntCenter = >> 179 fDetectorConstruction->GetTransferMgntCenter(); 150 180 151 auto simpleSolenoid = new F04SimpleSolenoid( << 181 F04SimpleSolenoid* simpleSolenoid = >> 182 new F04SimpleSolenoid(B, l, logicTransferMgnt,transferMgntCenter); 152 183 153 simpleSolenoid->SetColor("1,0,1"); 184 simpleSolenoid->SetColor("1,0,1"); 154 simpleSolenoid->SetColor("0,1,1"); 185 simpleSolenoid->SetColor("0,1,1"); 155 simpleSolenoid->SetMaxStep(1.5 * mm); << 186 simpleSolenoid->SetMaxStep(1.5*mm); 156 simpleSolenoid->SetMaxStep(2.5 * mm); << 187 simpleSolenoid->SetMaxStep(2.5*mm); 157 188 158 if (fFields) { 189 if (fFields) { 159 if (fFields->size() > 0) { << 190 if (fFields->size()>0) { 160 FieldList::iterator i; << 191 FieldList::iterator i; 161 for (i = fFields->begin(); i != fFields- << 192 for (i=fFields->begin(); i!=fFields->end(); ++i){ 162 (*i)->Construct(); << 193 (*i)->Construct(); 163 } << 194 } 164 } << 195 } 165 } 196 } 166 } 197 } 167 198 168 //....oooOO0OOooo........oooOO0OOooo........oo 199 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 169 200 170 F04GlobalField* F04GlobalField::GetObject(F04D 201 F04GlobalField* F04GlobalField::GetObject(F04DetectorConstruction* det) 171 { 202 { 172 if (!fObject) new F04GlobalField(det); 203 if (!fObject) new F04GlobalField(det); 173 return fObject; 204 return fObject; 174 } 205 } 175 206 176 //....oooOO0OOooo........oooOO0OOooo........oo 207 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 177 208 178 F04GlobalField* F04GlobalField::GetObject() 209 F04GlobalField* F04GlobalField::GetObject() 179 { 210 { 180 if (fObject) return fObject; 211 if (fObject) return fObject; 181 return nullptr; << 212 return NULL; 182 } 213 } 183 214 184 //....oooOO0OOooo........oooOO0OOooo........oo 215 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 185 216 186 void F04GlobalField::SetStepper() 217 void F04GlobalField::SetStepper() 187 { 218 { 188 delete fStepper; << 219 if(fStepper) delete fStepper; 189 220 190 switch (fStepperType) { << 221 switch ( fStepperType ) >> 222 { 191 case 0: 223 case 0: 192 // fStepper = new G4ExplicitEuler( << 224 // fStepper = new G4ExplicitEuler( fEquation, 8 ); // no spin tracking 193 fStepper = new G4ExplicitEuler(fEquation << 225 fStepper = new G4ExplicitEuler( fEquation, 12 ); // with spin tracking 194 G4cout << "G4ExplicitEuler is called" << 226 G4cout << "G4ExplicitEuler is called" << G4endl; 195 break; 227 break; 196 case 1: 228 case 1: 197 // fStepper = new G4ImplicitEuler( << 229 // fStepper = new G4ImplicitEuler( fEquation, 8 ); // no spin tracking 198 fStepper = new G4ImplicitEuler(fEquation << 230 fStepper = new G4ImplicitEuler( fEquation, 12 ); // with spin tracking 199 G4cout << "G4ImplicitEuler is called" << 231 G4cout << "G4ImplicitEuler is called" << G4endl; 200 break; 232 break; 201 case 2: 233 case 2: 202 // fStepper = new G4SimpleRunge( fE << 234 // fStepper = new G4SimpleRunge( fEquation, 8 ); // no spin tracking 203 fStepper = new G4SimpleRunge(fEquation, << 235 fStepper = new G4SimpleRunge( fEquation, 12 ); // with spin tracking 204 G4cout << "G4SimpleRunge is called" << G 236 G4cout << "G4SimpleRunge is called" << G4endl; 205 break; 237 break; 206 case 3: 238 case 3: 207 // fStepper = new G4SimpleHeum( fEq << 239 // fStepper = new G4SimpleHeum( fEquation, 8 ); // no spin tracking 208 fStepper = new G4SimpleHeum(fEquation, 1 << 240 fStepper = new G4SimpleHeum( fEquation, 12 ); // with spin tracking 209 G4cout << "G4SimpleHeum is called" << G4 241 G4cout << "G4SimpleHeum is called" << G4endl; 210 break; 242 break; 211 case 4: 243 case 4: 212 // fStepper = new G4ClassicalRK4( f << 244 // fStepper = new G4ClassicalRK4( fEquation, 8 ); // no spin tracking 213 fStepper = new G4ClassicalRK4(fEquation, << 245 fStepper = new G4ClassicalRK4( fEquation, 12 ); // with spin tracking 214 G4cout << "G4ClassicalRK4 (default) is c 246 G4cout << "G4ClassicalRK4 (default) is called" << G4endl; 215 break; 247 break; 216 case 5: 248 case 5: 217 // fStepper = new G4CashKarpRKF45( << 249 // fStepper = new G4CashKarpRKF45( fEquation, 8 ); // no spin tracking 218 fStepper = new G4CashKarpRKF45(fEquation << 250 fStepper = new G4CashKarpRKF45( fEquation, 12 ); // with spin tracking 219 G4cout << "G4CashKarpRKF45 is called" << 251 G4cout << "G4CashKarpRKF45 is called" << G4endl; 220 break; 252 break; 221 default: << 253 default: fStepper = 0; 222 fStepper = nullptr; << 223 } 254 } 224 } 255 } 225 256 226 //....oooOO0OOooo........oooOO0OOooo........oo 257 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 227 258 228 G4FieldManager* F04GlobalField::GetGlobalField 259 G4FieldManager* F04GlobalField::GetGlobalFieldManager() 229 { 260 { 230 return G4TransportationManager::GetTransport << 261 return G4TransportationManager::GetTransportationManager() >> 262 ->GetFieldManager(); 231 } 263 } 232 264 233 //....oooOO0OOooo........oooOO0OOooo........oo 265 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 234 266 235 void F04GlobalField::GetFieldValue(const G4dou 267 void F04GlobalField::GetFieldValue(const G4double* point, G4double* field) const 236 { 268 { 237 // NOTE: this routine dominates the CPU time 269 // NOTE: this routine dominates the CPU time for tracking. 238 // Using the simple array fFp[] instead of f << 270 // Using the simple array fFp[] instead of fields[] 239 // directly sped it up 271 // directly sped it up 240 272 241 field[0] = field[1] = field[2] = field[3] = 273 field[0] = field[1] = field[2] = field[3] = field[4] = field[5] = 0.0; 242 274 243 // protect against Geant4 bug that calls us 275 // protect against Geant4 bug that calls us with point[] NaN. 244 if (point[0] != point[0]) return; << 276 if(point[0] != point[0]) return; 245 277 246 // (can't use fNfp or fFp, as they may chang 278 // (can't use fNfp or fFp, as they may change) 247 if (fFirst) ((F04GlobalField*)this)->SetupAr << 279 if (fFirst) ((F04GlobalField*)this)->SetupArray(); // (cast away const) 248 280 249 for (int i = 0; i < fNfp; ++i) { << 281 for (int i=0; i<fNfp; ++i) { 250 const F04ElementField* p = fFp[i]; << 282 const F04ElementField* p = fFp[i]; 251 if (p->IsInBoundingBox(point)) { << 283 if (p->IsInBoundingBox(point)) { 252 p->AddFieldValue(point, field); << 284 p->AddFieldValue(point,field); 253 } << 285 } 254 } 286 } >> 287 255 } 288 } 256 289 257 //....oooOO0OOooo........oooOO0OOooo........oo 290 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 258 291 259 void F04GlobalField::Clear() 292 void F04GlobalField::Clear() 260 { 293 { 261 if (fFields) { 294 if (fFields) { 262 if (fFields->size() > 0) { << 295 if (fFields->size()>0) { 263 FieldList::iterator i; << 296 FieldList::iterator i; 264 for (i = fFields->begin(); i != fFields- << 297 for (i=fFields->begin(); i!=fFields->end(); ++i) delete *i; 265 delete *i; << 298 fFields->clear(); 266 fFields->clear(); << 299 } 267 } << 268 } 300 } 269 301 270 delete[] fFp; << 302 if (fFp) { delete [] fFp; } 271 fFirst = true; 303 fFirst = true; 272 fNfp = 0; 304 fNfp = 0; 273 fFp = nullptr; << 305 fFp = NULL; 274 } 306 } 275 307 276 //....oooOO0OOooo........oooOO0OOooo........oo 308 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 277 309 278 void F04GlobalField::SetupArray() 310 void F04GlobalField::SetupArray() 279 { 311 { 280 fFirst = false; 312 fFirst = false; 281 fNfp = fFields->size(); 313 fNfp = fFields->size(); 282 fFp = new const F04ElementField*[fNfp + 1]; << 314 fFp = new const F04ElementField* [fNfp+1]; // add 1 so it's never 0 283 for (int i = 0; i < fNfp; ++i) << 315 for (int i=0; i<fNfp; ++i) fFp[i] = (*fFields)[i]; 284 fFp[i] = (*fFields)[i]; << 285 } 316 } 286 317