field04/src/F04GlobalField.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 //
0027 /// \file field/field04/src/F04GlobalField.cc
0028 /// \brief Implementation of the F04GlobalField class
0029 //
0030
0031 #include "F04GlobalField.hh"
0032
0033 #include "F04FocusSolenoid.hh"
0034 #include "F04SimpleSolenoid.hh"
0035
0036 #include "G4CashKarpRKF45.hh"
0037 #include "G4ClassicalRK4.hh"
0038 #include "G4ExplicitEuler.hh"
0039 #include "G4ImplicitEuler.hh"
0040 #include "G4SimpleHeum.hh"
0041 #include "G4SimpleRunge.hh"
0042 #include "G4SystemOfUnits.hh"
0043 #include "G4TransportationManager.hh"
0044 #include "Randomize.hh"
0045
0046 #include <ctime>
0047
0048 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0049
0050 G4ThreadLocal F04GlobalField* F04GlobalField::fObject = nullptr;
0051
0052 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0053
0054 F04GlobalField::F04GlobalField(F04DetectorConstruction* det) : fDetectorConstruction(det)
0055 {
0056   fFieldMessenger = new F04FieldMessenger(this, det);
0057
0058   fFields = new FieldList();
0059
0060   //  set object
0061   fObject = this;
0062
0063   ConstructField();
0064 }
0065
0066 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0067
0068 F04GlobalField::~F04GlobalField()
0069 {
0070   Clear();
0071
0072   delete fFields;
0073
0074   delete fFieldMessenger;
0075
0076   delete fEquation;
0077   delete fStepper;
0078   delete fChordFinder;
0079 }
0080
0081 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0082
0083 void F04GlobalField::ConstructField()
0084 {
0085   Clear();
0086
0087   //  Construct equ. of motion of particles through B fields
0088   //  fEquation = new G4Mag_EqRhs(this);
0089   //  Construct equ. of motion of particles through e.m. fields
0090   //  fEquation = new G4EqMagElectricField(this);
0091   //  Construct equ. of motion of particles including spin through B fields
0092   //  fEquation = new G4Mag_SpinEqRhs(this);
0093   //  Construct equ. of motion of particles including spin through e.m. fields
0094   fEquation = new G4EqEMFieldWithSpin(this);
0095
0096   //  Get transportation, field, and propagator managers
0097   G4TransportationManager* transportManager = G4TransportationManager::GetTransportationManager();
0098
0099   fFieldManager = GetGlobalFieldManager();
0100
0101   fFieldPropagator = transportManager->GetPropagatorInField();
0102
0103   //  Need to SetFieldChangesEnergy to account for a time varying electric
0104   //  field (r.f. fields)
0105   fFieldManager->SetFieldChangesEnergy(true);
0106
0107   //  Set the field
0108   fFieldManager->SetDetectorField(this);
0109
0110   //  Choose a stepper for integration of the equation of motion
0111   SetStepper();
0112
0113   //  Create a cord finder providing the (global field, min step length,
0114   //  a pointer to the stepper)
0115   fChordFinder = new G4ChordFinder((G4MagneticField*)this, fMinStep, fStepper);
0116
0117   // Set accuracy parameters
0118   fChordFinder->SetDeltaChord(fDeltaChord);
0119
0120   fFieldManager->SetAccuraciesWithDeltaOneStep(fDeltaOneStep);
0121
0122   fFieldManager->SetDeltaIntersection(fDeltaIntersection);
0123
0124   fFieldPropagator->SetMinimumEpsilonStep(fEpsMin);
0125   fFieldPropagator->SetMaximumEpsilonStep(fEpsMax);
0126
0127   G4cout << "Accuracy Parameters:"
0128          << " MinStep=" << fMinStep << " DeltaChord=" << fDeltaChord
0129          << " DeltaOneStep=" << fDeltaOneStep << G4endl;
0130   G4cout << "                    "
0131          << " DeltaIntersection=" << fDeltaIntersection << " EpsMin=" << fEpsMin
0132          << " EpsMax=" << fEpsMax << G4endl;
0133
0134   fFieldManager->SetChordFinder(fChordFinder);
0135
0136   G4double l = 0.0;
0137   G4double B1 = fDetectorConstruction->GetCaptureMgntB1();
0138   G4double B2 = fDetectorConstruction->GetCaptureMgntB2();
0139
0140   G4LogicalVolume* logicCaptureMgnt = fDetectorConstruction->GetCaptureMgnt();
0141   G4ThreeVector captureMgntCenter = fDetectorConstruction->GetCaptureMgntCenter();
0142
0143   auto focusSolenoid = new F04FocusSolenoid(B1, B2, l, logicCaptureMgnt, captureMgntCenter);
0144   focusSolenoid->SetHalf(true);
0145
0146   G4double B = fDetectorConstruction->GetTransferMgntB();
0147
0148   G4LogicalVolume* logicTransferMgnt = fDetectorConstruction->GetTransferMgnt();
0149   G4ThreeVector transferMgntCenter = fDetectorConstruction->GetTransferMgntCenter();
0150
0151   auto simpleSolenoid = new F04SimpleSolenoid(B, l, logicTransferMgnt, transferMgntCenter);
0152
0153   simpleSolenoid->SetColor("1,0,1");
0154   simpleSolenoid->SetColor("0,1,1");
0155   simpleSolenoid->SetMaxStep(1.5 * mm);
0156   simpleSolenoid->SetMaxStep(2.5 * mm);
0157
0158   if (fFields) {
0159     if (fFields->size() > 0) {
0160       FieldList::iterator i;
0161       for (i = fFields->begin(); i != fFields->end(); ++i) {
0162         (*i)->Construct();
0163       }
0164     }
0165   }
0166 }
0167
0168 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0169
0170 F04GlobalField* F04GlobalField::GetObject(F04DetectorConstruction* det)
0171 {
0172   if (!fObject) new F04GlobalField(det);
0173   return fObject;
0174 }
0175
0176 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0177
0178 F04GlobalField* F04GlobalField::GetObject()
0179 {
0180   if (fObject) return fObject;
0181   return nullptr;
0182 }
0183
0184 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0185
0186 void F04GlobalField::SetStepper()
0187 {
0188   delete fStepper;
0189
0190   switch (fStepperType) {
0191     case 0:
0192       //      fStepper = new G4ExplicitEuler( fEquation, 8 ); // no spin tracking
0193       fStepper = new G4ExplicitEuler(fEquation, 12);  // with spin tracking
0194       G4cout << "G4ExplicitEuler is called" << G4endl;
0195       break;
0196     case 1:
0197       //      fStepper = new G4ImplicitEuler( fEquation, 8 ); // no spin tracking
0198       fStepper = new G4ImplicitEuler(fEquation, 12);  // with spin tracking
0199       G4cout << "G4ImplicitEuler is called" << G4endl;
0200       break;
0201     case 2:
0202       //      fStepper = new G4SimpleRunge( fEquation, 8 ); // no spin tracking
0203       fStepper = new G4SimpleRunge(fEquation, 12);  // with spin tracking
0204       G4cout << "G4SimpleRunge is called" << G4endl;
0205       break;
0206     case 3:
0207       //      fStepper = new G4SimpleHeum( fEquation, 8 ); // no spin tracking
0208       fStepper = new G4SimpleHeum(fEquation, 12);  // with spin tracking
0209       G4cout << "G4SimpleHeum is called" << G4endl;
0210       break;
0211     case 4:
0212       //      fStepper = new G4ClassicalRK4( fEquation, 8 ); // no spin tracking
0213       fStepper = new G4ClassicalRK4(fEquation, 12);  // with spin tracking
0214       G4cout << "G4ClassicalRK4 (default) is called" << G4endl;
0215       break;
0216     case 5:
0217       //      fStepper = new G4CashKarpRKF45( fEquation, 8 ); // no spin tracking
0218       fStepper = new G4CashKarpRKF45(fEquation, 12);  // with spin tracking
0219       G4cout << "G4CashKarpRKF45 is called" << G4endl;
0220       break;
0221     default:
0222       fStepper = nullptr;
0223   }
0224 }
0225
0226 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0227
0228 G4FieldManager* F04GlobalField::GetGlobalFieldManager()
0229 {
0230   return G4TransportationManager::GetTransportationManager()->GetFieldManager();
0231 }
0232
0233 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0234
0235 void F04GlobalField::GetFieldValue(const G4double* point, G4double* field) const
0236 {
0237   // NOTE: this routine dominates the CPU time for tracking.
0238   // Using the simple array fFp[] instead of fields[]
0239   // directly sped it up
0240
0241   field[0] = field[1] = field[2] = field[3] = field[4] = field[5] = 0.0;
0242
0243   // protect against Geant4 bug that calls us with point[] NaN.
0244   if (point[0] != point[0]) return;
0245
0246   // (can't use fNfp or fFp, as they may change)
0247   if (fFirst) ((F04GlobalField*)this)->SetupArray();  // (cast away const)
0248
0249   for (int i = 0; i < fNfp; ++i) {
0250     const F04ElementField* p = fFp[i];
0251     if (p->IsInBoundingBox(point)) {
0252       p->AddFieldValue(point, field);
0253     }
0254   }
0255 }
0256
0257 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0258
0259 void F04GlobalField::Clear()
0260 {
0261   if (fFields) {
0262     if (fFields->size() > 0) {
0263       FieldList::iterator i;
0264       for (i = fFields->begin(); i != fFields->end(); ++i)
0265         delete *i;
0266       fFields->clear();
0267     }
0268   }
0269
0270   delete[] fFp;
0271   fFirst = true;
0272   fNfp = 0;
0273   fFp = nullptr;
0274 }
0275
0276 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0277
0278 void F04GlobalField::SetupArray()
0279 {
0280   fFirst = false;
0281   fNfp = fFields->size();
0282   fFp = new const F04ElementField*[fNfp + 1];  // add 1 so it's never 0
0283   for (int i = 0; i < fNfp; ++i)
0284     fFp[i] = (*fFields)[i];
0285 }