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