| [807] | 1 |
|
|---|
| 2 | =========================================================
|
|---|
| 3 | Geant4 - an Object-Oriented Toolkit for Simulation in HEP
|
|---|
| 4 | =========================================================
|
|---|
| 5 |
|
|---|
| 6 |
|
|---|
| 7 |
|
|---|
| 8 | field04 Example
|
|---|
| 9 | ---------------
|
|---|
| 10 |
|
|---|
| 11 | This example shows how to define/use OVERLAPPING field elements
|
|---|
| 12 | in Geant4. Fields might be either magnetic, electric or both.
|
|---|
| 13 |
|
|---|
| 14 | Credit goes to Tom Roberts and Muons Inc. since much of the code
|
|---|
| 15 | and ideas were taken at liberty from the (GNU GPL) source of
|
|---|
| 16 | G4BEAMLINE release 1.12.
|
|---|
| 17 |
|
|---|
| 18 | http://g4beamline.muonsinc.com
|
|---|
| 19 |
|
|---|
| 20 | **************
|
|---|
| 21 | *Classes Used*
|
|---|
| 22 | **************
|
|---|
| 23 |
|
|---|
| 24 | 1 - main()
|
|---|
| 25 |
|
|---|
| 26 | Note:
|
|---|
| 27 |
|
|---|
| 28 | // runManager->Initialize();
|
|---|
| 29 |
|
|---|
| 30 | is commented out and must be executed with: /run/initialize
|
|---|
| 31 | either in a macro or on the command line. Note: One can change
|
|---|
| 32 | the geometry (see below) only BEFORE the initialization!
|
|---|
| 33 |
|
|---|
| 34 | Is is possible to assign the PhysicsList via an argument to
|
|---|
| 35 | program submission:
|
|---|
| 36 |
|
|---|
| 37 | field04 -p QGSP_BERT
|
|---|
| 38 |
|
|---|
| 39 | and an initial random number seed with:
|
|---|
| 40 |
|
|---|
| 41 | field04 field04.in 12345
|
|---|
| 42 |
|
|---|
| 43 | It is also possible to specify more than one macro file which are
|
|---|
| 44 | then executed in succession:
|
|---|
| 45 |
|
|---|
| 46 | field04 field04.in field04.in 12345
|
|---|
| 47 |
|
|---|
| 48 | (Note, in this case the initial random # must also be provided)
|
|---|
| 49 |
|
|---|
| 50 | (1) argc holds the number of arguments on the command line. Since
|
|---|
| 51 | this MUST include at least the program name, argc value is >= 1!
|
|---|
| 52 | (2) argv[0] is always the name of the program
|
|---|
| 53 | argv[1] points to the first argument, and so on ...
|
|---|
| 54 | (3) argv[argc-1] is assumed to be the intial random # seed
|
|---|
| 55 | (4) the loop over macro files goes from optind = 1 to (argc-1)
|
|---|
| 56 |
|
|---|
| 57 |
|
|---|
| 58 | 2- GEOMETRY DEFINITION
|
|---|
| 59 |
|
|---|
| 60 | The geometry consists of two solenoidal magnets: a "CaptureMgnt"
|
|---|
| 61 | followed by a (bloe-colored "TransferMgnt". By definition, the
|
|---|
| 62 | axis and center of the "CaptureMgnt" coincide with the "World". The
|
|---|
| 63 | position of the "TransferMgnt" relative to the downstream end of the
|
|---|
| 64 | "CaptureMgnt", as well as its axis angle, both may vary. A cylindrical
|
|---|
| 65 | "Target" is positioned inside the "CaptureMgnt". Its axis can vary
|
|---|
| 66 | from 0 to 180 deg, and hence also the direction of the incoming
|
|---|
| 67 | proton beam wrt the "CaptureMgnt"'s axis. A "Degrader" is located
|
|---|
| 68 | inside the "TransferMgnt", its default position being at the
|
|---|
| 69 | upstream end of the "TransferMgnt". Finally, also a "TestPlane" is
|
|---|
| 70 | located inside the "TransferMgnt", by default at its downstream end.
|
|---|
| 71 |
|
|---|
| 72 |
|
|---|
| 73 | The "World" consists of a solid cylinder made of a given material.
|
|---|
| 74 | (It is the responsibility of the user to make the world
|
|---|
| 75 | large enough to contain the rest of the geometry!)
|
|---|
| 76 |
|
|---|
| 77 | Three parameters define the world :
|
|---|
| 78 | - the material of the world,
|
|---|
| 79 | - the world radius,
|
|---|
| 80 | - the world length.
|
|---|
| 81 |
|
|---|
| 82 | Example (default values):
|
|---|
| 83 | /field04/SetWorldMat G4_AIR
|
|---|
| 84 | /field04/SetWorldR 5.0 m
|
|---|
| 85 | /field04/SetWorldZ 50.0 m
|
|---|
| 86 |
|
|---|
| 87 |
|
|---|
| 88 | The "Target" is a solid cylinder made of a given material.
|
|---|
| 89 |
|
|---|
| 90 | Five parameters define the target:
|
|---|
| 91 | - the material of the target,
|
|---|
| 92 | - the target radius,
|
|---|
| 93 | - the target thickness,
|
|---|
| 94 | - the target position inside the "CaptureMgnt",
|
|---|
| 95 | - the target axis angle relative to that of the "CaptureMgnt".
|
|---|
| 96 |
|
|---|
| 97 | Example (default values):
|
|---|
| 98 | /field04/SetTgtMat G4_W
|
|---|
| 99 | /field04/SetTgtRad 0.4 cm
|
|---|
| 100 | /field04/SetTgtThick 16.0 cm
|
|---|
| 101 | /field04/SetTgtPos 0.0 cm
|
|---|
| 102 | /field04/SetTgtAng 170
|
|---|
| 103 |
|
|---|
| 104 |
|
|---|
| 105 | The "Degrader" is a solid cylinder made of a given material.
|
|---|
| 106 |
|
|---|
| 107 | Four parameters define the degrader:
|
|---|
| 108 | - the material of the degrader,
|
|---|
| 109 | - the degrader radius,
|
|---|
| 110 | - the degrader thickness,
|
|---|
| 111 | - the degrader position relative to the "TransferMgnt" center.
|
|---|
| 112 |
|
|---|
| 113 | Example (default values):
|
|---|
| 114 | /field04/SetDgrMat G4_Pb
|
|---|
| 115 | /field04/SetDgrRad 30.0 cm
|
|---|
| 116 | /field04/SetDgrThick 0.1 cm
|
|---|
| 117 | #/field04/SetDgrPos -7.4 m
|
|---|
| 118 |
|
|---|
| 119 |
|
|---|
| 120 | The "CaptureMgnt" is a solenoid (vacuum cylinder). It is either
|
|---|
| 121 | a two-sided or a one-sided magnetic bottle with the B field
|
|---|
| 122 | varying linearly from the center value B1 to the edge value B2.
|
|---|
| 123 | The one-sided 'FocusSolenoid' has the open end at +z and focuses
|
|---|
| 124 | on the z < 0 side.
|
|---|
| 125 |
|
|---|
| 126 | Four parameters define the "CaptureMgnt":
|
|---|
| 127 | - the magnet radius,
|
|---|
| 128 | - the magnet length,
|
|---|
| 129 | - the weaker magnetic field at the center B1
|
|---|
| 130 | - the stronger magnetic field at the edge B2
|
|---|
| 131 |
|
|---|
| 132 | Example (default values):
|
|---|
| 133 | /field04/SetCaptureR 0.6 m
|
|---|
| 134 | /field04/SetCaptureZ 4.0 m
|
|---|
| 135 | /field04/SetCaptureB1 2.5 tesla
|
|---|
| 136 | /field04/SetCaptureB2 5.0 tesla
|
|---|
| 137 |
|
|---|
| 138 |
|
|---|
| 139 | The "TransferMgnt" is a solenoid (vacuum cylinder) with a
|
|---|
| 140 | constant B-field. When the "TransferMgnt" follows immediately
|
|---|
| 141 | the "CaptureMgnt", its relative position is at 0 cm.
|
|---|
| 142 |
|
|---|
| 143 | Four parameters define the "TransferMgnt":
|
|---|
| 144 | - the magnet radius,
|
|---|
| 145 | - the magnet length,
|
|---|
| 146 | - the magnet field,
|
|---|
| 147 | - the magnet relative position
|
|---|
| 148 | (its upstream face wrt the downstream face of the "CaptureMgnt".)
|
|---|
| 149 |
|
|---|
| 150 | Example (default values):
|
|---|
| 151 | /field04/SetTransferR 0.3 m
|
|---|
| 152 | /field04/SetTransferZ 15.0 m
|
|---|
| 153 | /field04/SetTransferB 5.0 tesla
|
|---|
| 154 | /field04/SetTransferP 0.0 m
|
|---|
| 155 |
|
|---|
| 156 | The default geometry is constructed in DetectorConstruction class,
|
|---|
| 157 | but all the parameters can be changed via the commands defined in
|
|---|
| 158 | the DetectorMessenger class.
|
|---|
| 159 |
|
|---|
| 160 |
|
|---|
| 161 | 3- MATERIAL DEFINITION
|
|---|
| 162 |
|
|---|
| 163 | Material definitions are done through the singleton class Materials
|
|---|
| 164 | which keeps a pointer to the G4NistManager. It has a method
|
|---|
| 165 | GetMaterial by name (G4String) which in turn invokes the
|
|---|
| 166 | G4NistManager::FindOrBuildMaterial, and/or G4Material::GetMaterial
|
|---|
| 167 | methods. It has also a method CreateMaterials which, for materials
|
|---|
| 168 | absent from the NIST data base, shows how to create them using the
|
|---|
| 169 | G4NistManager::ConstructNewMaterial method.
|
|---|
| 170 |
|
|---|
| 171 |
|
|---|
| 172 | 4- AN EVENT: THE PRIMARY GENERATOR
|
|---|
| 173 |
|
|---|
| 174 | The primary kinematic consists of a single particle which hits the
|
|---|
| 175 | target perpendicular to its upstream face. The type of the particle
|
|---|
| 176 | and its energy are set in the PrimaryGeneratorAction class, and can
|
|---|
| 177 | be changed via the G4 build-in commands of the ParticleGun class.
|
|---|
| 178 | In addition, there is a rndmFlag, which once set allows the beam to
|
|---|
| 179 | explore randomly the whole cross section of the target. The default
|
|---|
| 180 | beam consists of 500 MeV protons, starting at the upstream face of
|
|---|
| 181 | the target, directed along dx = dy = 0, dz = 1 wrt the target frame.
|
|---|
| 182 | The default direction should NOT be changed! The arguments of the
|
|---|
| 183 | x/y/zvertex commands are relative to the target center.
|
|---|
| 184 |
|
|---|
| 185 | Example:
|
|---|
| 186 | /gun/random on
|
|---|
| 187 | #/gun/xvertex 0 mm
|
|---|
| 188 | #/gun/yvertex 0 mm
|
|---|
| 189 | #/gun/zvertex -100 mm
|
|---|
| 190 |
|
|---|
| 191 |
|
|---|
| 192 | 5- DETECTOR RESPONSE
|
|---|
| 193 |
|
|---|
| 194 | Information is extracted from the program via SteppingAction
|
|---|
| 195 | at the TestPlane.
|
|---|
| 196 |
|
|---|
| 197 |
|
|---|
| 198 | 6- PHYSICS
|
|---|
| 199 |
|
|---|
| 200 | The PhysicsList is adopted from examples/extended/hadronic/Hadr01.
|
|---|
| 201 | It uses components which are distributed with Geant4 in
|
|---|
| 202 | /geant4/physics_lists subdirectory. So, before compiling field04
|
|---|
| 203 | it is necessary to compile physics_lists.
|
|---|
| 204 |
|
|---|
| 205 | There are two std::vector<G4VPhysicsConstructor*> PhysicsListVectors;
|
|---|
| 206 | one for EM- and one for HadronPHysics and a user defined process
|
|---|
| 207 | "StepMax".
|
|---|
| 208 |
|
|---|
| 209 | The choice of the physics is provided by the UI command:
|
|---|
| 210 |
|
|---|
| 211 | /exp/phys/addPhysics emstandard_opt1
|
|---|
| 212 | /exp/phys/addPhysics QGSP_BERT
|
|---|
| 213 |
|
|---|
| 214 | To see the list of available configurations one can use:
|
|---|
| 215 |
|
|---|
| 216 | /exp/phys/list
|
|---|
| 217 |
|
|---|
| 218 | The command:
|
|---|
| 219 |
|
|---|
| 220 | /exp/phys/addPhysics PHYSLIST
|
|---|
| 221 |
|
|---|
| 222 | allows to download a physics configuration defined by an
|
|---|
| 223 | environment variable PHYSLIST.
|
|---|
| 224 |
|
|---|
| 225 | The cuts for electromagnetic phsyics can be established via:
|
|---|
| 226 |
|
|---|
| 227 | /exp/phys/allCuts 1 mm
|
|---|
| 228 | /exp/phys/gammaCut 0.1 mm
|
|---|
| 229 | /exp/phys/electronCut 0.2 mm
|
|---|
| 230 | /exp/phys/positronCut 0.3 mm
|
|---|
| 231 |
|
|---|
| 232 | The cut for the StepMax process via:
|
|---|
| 233 |
|
|---|
| 234 | /exp/phys/stepMax
|
|---|
| 235 |
|
|---|
| 236 | The two PhysicsListVectors can be cleared via:
|
|---|
| 237 |
|
|---|
| 238 | /exp/phys/clearEMPhysics
|
|---|
| 239 | /exp/phys/clearHadronPhysics
|
|---|
| 240 |
|
|---|
| 241 | and individual processes can be removed from the vectors via,
|
|---|
| 242 | for example:
|
|---|
| 243 |
|
|---|
| 244 | /exp/phys/removeEMPhysics msc
|
|---|
| 245 | /exp/phys/removeHadronPhysics gamma_nuc
|
|---|
| 246 |
|
|---|
| 247 | Furthermore, the following commands are also available, but
|
|---|
| 248 | may only be used AFTER /run/initialize
|
|---|
| 249 |
|
|---|
| 250 | /process/inactivate msc
|
|---|
| 251 | /process/activate msc
|
|---|
| 252 |
|
|---|
| 253 | The decay of pions can be assigned via (pi -> e nu, pi -> mu nu):
|
|---|
| 254 |
|
|---|
| 255 | /decay/pienu
|
|---|
| 256 | /decay/pimunu
|
|---|
| 257 |
|
|---|
| 258 | The pienu assignment includes a small fraction of radiative decay:
|
|---|
| 259 | e nu gamma (G4PionRadiativeDecayChannel).
|
|---|
| 260 |
|
|---|
| 261 | The standard/default muon decay chain is modified to be 98.6%
|
|---|
| 262 | G4MuonDecayChannelWithSpin and 1.4% G4MuonRadiativeDecayChannelWithSpin
|
|---|
| 263 | in ConstructParticle().
|
|---|
| 264 |
|
|---|
| 265 | The pion decay process G4PolDecay inherits from G4Decay and implements
|
|---|
| 266 | the virtual method - empty in the base class - DaughterPolarization
|
|---|
| 267 |
|
|---|
| 268 | The muon decay process is G4DecayWithSpin
|
|---|
| 269 |
|
|---|
| 270 |
|
|---|
| 271 | 7- Overlapping Fields
|
|---|
| 272 |
|
|---|
| 273 | The GlobalField (a singleton) is instantiated in DetectorConstruction()
|
|---|
| 274 | and assigned to the global field manager in updateField():
|
|---|
| 275 |
|
|---|
| 276 | fFieldManager = GetGlobalFieldManager();
|
|---|
| 277 | fFieldManager->SetDetectorField(this);
|
|---|
| 278 |
|
|---|
| 279 | The GlobalField has a std::vector<ElementField*> FieldList
|
|---|
| 280 |
|
|---|
| 281 | The field from each individual beamline element is given by a
|
|---|
| 282 | ElementField object. Any number of overlapping ElementField
|
|---|
| 283 | objects can be added to the global field. Any element that
|
|---|
| 284 | represents an element with an EM field must add the appropriate
|
|---|
| 285 | ElementField to the global GlobalField object.
|
|---|
| 286 |
|
|---|
| 287 | Of course, the GlobalField has the method GetFieldValue implemented.
|
|---|
| 288 |
|
|---|
| 289 | Before /run/initialize in the macro file or command, the update
|
|---|
| 290 | field command must have been issued if any of the other following
|
|---|
| 291 | field commands was employed:
|
|---|
| 292 |
|
|---|
| 293 | /field/update
|
|---|
| 294 |
|
|---|
| 295 | Other options are:
|
|---|
| 296 |
|
|---|
| 297 | /field/setStepperType 4
|
|---|
| 298 | /field/setMinStep 10 mm
|
|---|
| 299 | /field/setDeltaChord 3.0 mm
|
|---|
| 300 | /field/setDeltaOneStep 0.01 mm
|
|---|
| 301 | /field/setDeltaIntersection 0.1 mm
|
|---|
| 302 | /field/setEpsMin 2.5e-7 mm
|
|---|
| 303 | /field/setEpsMax 0.05 mm
|
|---|
| 304 |
|
|---|
| 305 | Each field element has a rectilinear bounding box in global
|
|---|
| 306 | coordinate space which is checked before a point is verified to
|
|---|
| 307 | actually be inside the ElementField (isWithin and isOutside).
|
|---|
| 308 | setGlobalPoint is called 8 times for the corners of the local
|
|---|
| 309 | bounding box, after a local->global coordinate transform.
|
|---|
| 310 |
|
|---|
| 311 | The ElementField is the interface class used by GlobalField to
|
|---|
| 312 | compute the field value at a given point[].
|
|---|
| 313 |
|
|---|
| 314 | A beamline element, for example the SimpleSolenoid, will derive
|
|---|
| 315 | from ElementField and implement the computation for the element.
|
|---|
| 316 |
|
|---|
| 317 | simpleSolenoid = new SimpleSolenoid(B, l,
|
|---|
| 318 | logicTransferMgnt,TransferMgntCenter);
|
|---|
| 319 |
|
|---|
| 320 | Besides the magnetic field and the length of the simple solenoid,
|
|---|
| 321 | the constructor needs the knowledge of the G4LogicalVolume for
|
|---|
| 322 | the beamline element and where its center is located in the
|
|---|
| 323 | 'World'.
|
|---|
| 324 |
|
|---|
| 325 | The ElementField has a G4AffineTransform "global2local" which
|
|---|
| 326 | allows the quick computation of coordinate transformations. It can
|
|---|
| 327 | only be determined by knowing the element's coordinate origin in
|
|---|
| 328 | the global frame and after all of the geometry has been defined.
|
|---|
| 329 | For this reason, the object is prepared in two stages, through the
|
|---|
| 330 | constructor providing it with the coordinate center and a pointer
|
|---|
| 331 | to the G4LogicalVolume. Later the construct() method is called to
|
|---|
| 332 | calculate the global2local and the bounding box. This can be done
|
|---|
| 333 | from the RunAction::BeginOfRunAction method, for only then are we
|
|---|
| 334 | certain that the geometry has been completely built:
|
|---|
| 335 |
|
|---|
| 336 | FieldList* fields = GlobalField::getObject()->getFields();
|
|---|
| 337 |
|
|---|
| 338 | if (fields) {
|
|---|
| 339 | if (fields->size()>0) {
|
|---|
| 340 | FieldList::iterator i;
|
|---|
| 341 | for (i=fields->begin(); i!=fields->end(); ++i)(*i)->construct();
|
|---|
| 342 | }
|
|---|
| 343 | }
|
|---|
| 344 |
|
|---|
| 345 | The ElementField constructor will also add the derived object into
|
|---|
| 346 | GlobalField. Finally, its addFieldValue() will add the field value
|
|---|
| 347 | for this element to field[].
|
|---|
| 348 |
|
|---|
| 349 |
|
|---|
| 350 | 8- User Action Classes
|
|---|
| 351 |
|
|---|
| 352 | RunActionMessenger:
|
|---|
| 353 |
|
|---|
| 354 | /rndm/save freq - to save rndm status in external files
|
|---|
| 355 | 0 not saved
|
|---|
| 356 | >0 saved on: beginOfRun.rndm
|
|---|
| 357 | 1 saved on: endOfRun.rndm
|
|---|
| 358 | 2 saved on: endOfEvent.rndm
|
|---|
| 359 | /rndm/read random/run0evt8268.rndm
|
|---|
| 360 |
|
|---|
| 361 | RunAction:
|
|---|
| 362 | BeginOfRunAction: Deal with random number storage,
|
|---|
| 363 | initialization etc. Call the construct() method of
|
|---|
| 364 | ElementFields in the FieldList of GlobalField object.
|
|---|
| 365 | EndOfRunAction: random number storage/status printing.
|
|---|
| 366 |
|
|---|
| 367 | EventActionMessenger:
|
|---|
| 368 | /event/setverbose
|
|---|
| 369 | /event/drawTracks none/charged/all
|
|---|
| 370 | /event/printModulo
|
|---|
| 371 |
|
|---|
| 372 | EventAction(RunAction* RA):
|
|---|
| 373 | Customized BeginOfEvent printing (frequency: printModulo)
|
|---|
| 374 | EndofEvent: condition depending drawing of trajectories,
|
|---|
| 375 | saveEngingStatus and showEngineStatus according to flag
|
|---|
| 376 | in RunAction
|
|---|
| 377 |
|
|---|
| 378 | TrackingAction:
|
|---|
| 379 | PreUserTrackingAction: Instantiate UserTrackInformation
|
|---|
| 380 | and set the application TrackStatus.
|
|---|
| 381 | PostUserTrackingAction: Retreive UserTrackInformation
|
|---|
| 382 | and decide to save random number status accordingly.
|
|---|
| 383 |
|
|---|
| 384 | SteppingActionMessenger:
|
|---|
| 385 |
|
|---|
| 386 | SteppingAction:
|
|---|
| 387 | UserSteppingAction: Kill primary if/when outside Target
|
|---|
| 388 | volume. Diagnostic/histogram filling for particles at a
|
|---|
| 389 | TestPlane. Find decay position and when particle
|
|---|
| 390 | FIRST reverses z-momentum component via using a
|
|---|
| 391 | UserTrackInformation object.
|
|---|
| 392 |
|
|---|
| 393 | StackingAction:
|
|---|
| 394 | Track only primaries, pi+ or mu+
|
|---|
| 395 |
|
|---|
| 396 | UserTrackInformation:
|
|---|
| 397 | Keep a application TrackStatus for the track:
|
|---|
| 398 | undefined, left, right, reverse
|
|---|
| 399 |
|
|---|
| 400 | SteppingVerbose:
|
|---|
| 401 | Only print track header and step information for
|
|---|
| 402 | pi+ and mu+
|
|---|
| 403 |
|
|---|
| 404 | Trajectory, TrajectoryPoint:
|
|---|
| 405 | Example of application specific implementations
|
|---|
| 406 |
|
|---|
| 407 | 9- HOW TO START ?
|
|---|
| 408 |
|
|---|
| 409 | - compile and link to generate an executable
|
|---|
| 410 | % cd $G4INSTALL/example/extended/field/field04
|
|---|
| 411 | % gmake
|
|---|
| 412 |
|
|---|
| 413 | - execute field04 in 'batch' mode from macro files e.g.
|
|---|
| 414 | % $(G4INSTALL)/bin/$(G4SYSTEM)/field04 field04.in
|
|---|
| 415 |
|
|---|
| 416 | - execute field04 in 'interactive' mode with visualization e.g.
|
|---|
| 417 | % $(G4INSTALL)/bin/$(G4SYSTEM)/field04
|
|---|
| 418 | ....
|
|---|
| 419 | Idle> type your commands, for example:
|
|---|
| 420 | Idle> control/execute field04.in
|
|---|
| 421 | ....
|
|---|