Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 // 27 /// \file field/field04/src/F04GlobalField.cc 28 /// \brief Implementation of the F04GlobalFiel 29 // 30 31 #include "F04GlobalField.hh" 32 33 #include "F04FocusSolenoid.hh" 34 #include "F04SimpleSolenoid.hh" 35 36 #include "G4CashKarpRKF45.hh" 37 #include "G4ClassicalRK4.hh" 38 #include "G4ExplicitEuler.hh" 39 #include "G4ImplicitEuler.hh" 40 #include "G4SimpleHeum.hh" 41 #include "G4SimpleRunge.hh" 42 #include "G4SystemOfUnits.hh" 43 #include "G4TransportationManager.hh" 44 #include "Randomize.hh" 45 46 #include <ctime> 47 48 //....oooOO0OOooo........oooOO0OOooo........oo 49 50 G4ThreadLocal F04GlobalField* F04GlobalField:: 51 52 //....oooOO0OOooo........oooOO0OOooo........oo 53 54 F04GlobalField::F04GlobalField(F04DetectorCons 55 { 56 fFieldMessenger = new F04FieldMessenger(this 57 58 fFields = new FieldList(); 59 60 // set object 61 fObject = this; 62 63 ConstructField(); 64 } 65 66 //....oooOO0OOooo........oooOO0OOooo........oo 67 68 F04GlobalField::~F04GlobalField() 69 { 70 Clear(); 71 72 delete fFields; 73 74 delete fFieldMessenger; 75 76 delete fEquation; 77 delete fStepper; 78 delete fChordFinder; 79 } 80 81 //....oooOO0OOooo........oooOO0OOooo........oo 82 83 void F04GlobalField::ConstructField() 84 { 85 Clear(); 86 87 // Construct equ. of motion of particles th 88 // fEquation = new G4Mag_EqRhs(this); 89 // Construct equ. of motion of particles th 90 // fEquation = new G4EqMagElectricField(thi 91 // Construct equ. of motion of particles in 92 // fEquation = new G4Mag_SpinEqRhs(this); 93 // Construct equ. of motion of particles in 94 fEquation = new G4EqEMFieldWithSpin(this); 95 96 // Get transportation, field, and propagato 97 G4TransportationManager* transportManager = 98 99 fFieldManager = GetGlobalFieldManager(); 100 101 fFieldPropagator = transportManager->GetProp 102 103 // Need to SetFieldChangesEnergy to account 104 // field (r.f. fields) 105 fFieldManager->SetFieldChangesEnergy(true); 106 107 // Set the field 108 fFieldManager->SetDetectorField(this); 109 110 // Choose a stepper for integration of the 111 SetStepper(); 112 113 // Create a cord finder providing the (glob 114 // a pointer to the stepper) 115 fChordFinder = new G4ChordFinder((G4Magnetic 116 117 // Set accuracy parameters 118 fChordFinder->SetDeltaChord(fDeltaChord); 119 120 fFieldManager->SetAccuraciesWithDeltaOneStep 121 122 fFieldManager->SetDeltaIntersection(fDeltaIn 123 124 fFieldPropagator->SetMinimumEpsilonStep(fEps 125 fFieldPropagator->SetMaximumEpsilonStep(fEps 126 127 G4cout << "Accuracy Parameters:" 128 << " MinStep=" << fMinStep << " Delta 129 << " DeltaOneStep=" << fDeltaOneStep 130 G4cout << " " 131 << " DeltaIntersection=" << fDeltaInt 132 << " EpsMax=" << fEpsMax << G4endl; 133 134 fFieldManager->SetChordFinder(fChordFinder); 135 136 G4double l = 0.0; 137 G4double B1 = fDetectorConstruction->GetCapt 138 G4double B2 = fDetectorConstruction->GetCapt 139 140 G4LogicalVolume* logicCaptureMgnt = fDetecto 141 G4ThreeVector captureMgntCenter = fDetectorC 142 143 auto focusSolenoid = new F04FocusSolenoid(B1 144 focusSolenoid->SetHalf(true); 145 146 G4double B = fDetectorConstruction->GetTrans 147 148 G4LogicalVolume* logicTransferMgnt = fDetect 149 G4ThreeVector transferMgntCenter = fDetector 150 151 auto simpleSolenoid = new F04SimpleSolenoid( 152 153 simpleSolenoid->SetColor("1,0,1"); 154 simpleSolenoid->SetColor("0,1,1"); 155 simpleSolenoid->SetMaxStep(1.5 * mm); 156 simpleSolenoid->SetMaxStep(2.5 * mm); 157 158 if (fFields) { 159 if (fFields->size() > 0) { 160 FieldList::iterator i; 161 for (i = fFields->begin(); i != fFields- 162 (*i)->Construct(); 163 } 164 } 165 } 166 } 167 168 //....oooOO0OOooo........oooOO0OOooo........oo 169 170 F04GlobalField* F04GlobalField::GetObject(F04D 171 { 172 if (!fObject) new F04GlobalField(det); 173 return fObject; 174 } 175 176 //....oooOO0OOooo........oooOO0OOooo........oo 177 178 F04GlobalField* F04GlobalField::GetObject() 179 { 180 if (fObject) return fObject; 181 return nullptr; 182 } 183 184 //....oooOO0OOooo........oooOO0OOooo........oo 185 186 void F04GlobalField::SetStepper() 187 { 188 delete fStepper; 189 190 switch (fStepperType) { 191 case 0: 192 // fStepper = new G4ExplicitEuler( 193 fStepper = new G4ExplicitEuler(fEquation 194 G4cout << "G4ExplicitEuler is called" << 195 break; 196 case 1: 197 // fStepper = new G4ImplicitEuler( 198 fStepper = new G4ImplicitEuler(fEquation 199 G4cout << "G4ImplicitEuler is called" << 200 break; 201 case 2: 202 // fStepper = new G4SimpleRunge( fE 203 fStepper = new G4SimpleRunge(fEquation, 204 G4cout << "G4SimpleRunge is called" << G 205 break; 206 case 3: 207 // fStepper = new G4SimpleHeum( fEq 208 fStepper = new G4SimpleHeum(fEquation, 1 209 G4cout << "G4SimpleHeum is called" << G4 210 break; 211 case 4: 212 // fStepper = new G4ClassicalRK4( f 213 fStepper = new G4ClassicalRK4(fEquation, 214 G4cout << "G4ClassicalRK4 (default) is c 215 break; 216 case 5: 217 // fStepper = new G4CashKarpRKF45( 218 fStepper = new G4CashKarpRKF45(fEquation 219 G4cout << "G4CashKarpRKF45 is called" << 220 break; 221 default: 222 fStepper = nullptr; 223 } 224 } 225 226 //....oooOO0OOooo........oooOO0OOooo........oo 227 228 G4FieldManager* F04GlobalField::GetGlobalField 229 { 230 return G4TransportationManager::GetTransport 231 } 232 233 //....oooOO0OOooo........oooOO0OOooo........oo 234 235 void F04GlobalField::GetFieldValue(const G4dou 236 { 237 // NOTE: this routine dominates the CPU time 238 // Using the simple array fFp[] instead of f 239 // directly sped it up 240 241 field[0] = field[1] = field[2] = field[3] = 242 243 // protect against Geant4 bug that calls us 244 if (point[0] != point[0]) return; 245 246 // (can't use fNfp or fFp, as they may chang 247 if (fFirst) ((F04GlobalField*)this)->SetupAr 248 249 for (int i = 0; i < fNfp; ++i) { 250 const F04ElementField* p = fFp[i]; 251 if (p->IsInBoundingBox(point)) { 252 p->AddFieldValue(point, field); 253 } 254 } 255 } 256 257 //....oooOO0OOooo........oooOO0OOooo........oo 258 259 void F04GlobalField::Clear() 260 { 261 if (fFields) { 262 if (fFields->size() > 0) { 263 FieldList::iterator i; 264 for (i = fFields->begin(); i != fFields- 265 delete *i; 266 fFields->clear(); 267 } 268 } 269 270 delete[] fFp; 271 fFirst = true; 272 fNfp = 0; 273 fFp = nullptr; 274 } 275 276 //....oooOO0OOooo........oooOO0OOooo........oo 277 278 void F04GlobalField::SetupArray() 279 { 280 fFirst = false; 281 fNfp = fFields->size(); 282 fFp = new const F04ElementField*[fNfp + 1]; 283 for (int i = 0; i < fNfp; ++i) 284 fFp[i] = (*fFields)[i]; 285 } 286