g4simulation/g4main/G4TBMagneticFieldSetup.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 // $Id: G4TBMagneticFieldSetup.cc,v 1.24 2014/11/14 21:47:38 mccumber Exp $
0028 // GEANT4 tag $Name:  $
0029 //
0030 //
0031 //   User Field class implementation.
0032 //
0033
0034 #include "G4TBMagneticFieldSetup.hh"
0035
0036 #include "G4TBFieldMessenger.hh"
0037 #include "PHG4MagneticField.h"
0038
0039 #include <Geant4/G4CashKarpRKF45.hh>
0040 #include <Geant4/G4ChordFinder.hh>
0041 #include <Geant4/G4ClassicalRK4.hh>
0042 #include <Geant4/G4ExplicitEuler.hh>
0043 #include <Geant4/G4FieldManager.hh>
0044 #include <Geant4/G4ImplicitEuler.hh>
0045 #include <Geant4/G4MagIntegratorDriver.hh>
0046 #include <Geant4/G4MagIntegratorStepper.hh>
0047 #include <Geant4/G4Mag_UsualEqRhs.hh>
0048 #include <Geant4/G4MagneticField.hh>
0049 #include <Geant4/G4SimpleHeum.hh>
0050 #include <Geant4/G4SimpleRunge.hh>
0051 #include <Geant4/G4SystemOfUnits.hh>
0052 #include <Geant4/G4ThreeVector.hh>
0053 #include <Geant4/G4TransportationManager.hh>
0054 #include <Geant4/G4Types.hh>  // for G4double, G4int
0055 #include <Geant4/G4UniformMagField.hh>
0056
0057 #include <cassert>
0058 #include <cstdlib>  // for exit, size_t
0059 #include <iostream>
0060 #include <string>
0061
0062 //////////////////////////////////////////////////////////////////////////
0063 //
0064 //  Constructors:
0065
0066 G4TBMagneticFieldSetup::G4TBMagneticFieldSetup(PHField* phfield)
0067 {
0068   assert(phfield);
0069
0070   fEMfield = new PHG4MagneticField(phfield);
0071   fFieldMessenger = new G4TBFieldMessenger(this);
0072   fEquation = new G4Mag_UsualEqRhs(fEMfield);
0073   fMinStep = 0.005 * mm;  // minimal step of 5 microns
0074   fStepperType = 4;       // ClassicalRK4 -- the default stepper
0075
0076   fFieldManager = GetGlobalFieldManager();
0077   UpdateField();
0078   double point[4] = {0, 0, 0, 0};
0079   fEMfield->GetFieldValue(&point[0], &magfield_at_000[0]);
0080   for (double& i : magfield_at_000)
0081   {
0082     i = i / tesla;
0083   }
0084   if (verbosity > 0)
0085   {
0086     std::cout << "field: x" << magfield_at_000[0]
0087               << ", y: " << magfield_at_000[1]
0088               << ", z: " << magfield_at_000[2]
0089               << std::endl;
0090   }
0091 }
0092
0093 // G4TBMagneticFieldSetup::G4TBMagneticFieldSetup(const float magfield)
0094 //   : verbosity(0), fChordFinder(0), fStepper(0), fIntgrDriver(0)
0095 //{
0096 //   //solenoidal field along the axis of the cyclinders?
0097 //   fEMfield = new G4UniformMagField(G4ThreeVector(0.0, 0.0, magfield*tesla));
0098 //   fFieldMessenger = new G4TBFieldMessenger(this) ;
0099 //   fEquation = new  G4Mag_UsualEqRhs(fEMfield);
0100 //   fMinStep     = 0.005 * mm ; // minimal step of 10 microns
0101 //   fStepperType = 4 ;        // ClassicalRK4 -- the default stepper
0102 //
0103 //   fFieldManager = GetGlobalFieldManager();
0104 //   UpdateField();
0105 //
0106 //   double point[4] = {0,0,0,0};
0107 //   fEMfield->GetFieldValue(&point[0],&magfield_at_000[0]);
0108 //   for (size_t i=0; i<sizeof(magfield_at_000)/sizeof(double);i++)
0109 //     {
0110 //       magfield_at_000[i] = magfield_at_000[i]/tesla;
0111 //     }
0112 // }
0113 //
0114 ///////////////////////////////////////////////////////////////////////////////////
0115 //
0116 // G4TBMagneticFieldSetup::G4TBMagneticFieldSetup(const string &fieldmapname, const int dim, const float magfield_rescale)
0117 //  : verbosity(0),
0118 //    fChordFinder(0),
0119 //    fStepper(0),
0120 //    fIntgrDriver(0)
0121 //{
0122 //
0123 //  switch(dim)
0124 //    {
0125 //    case 1:
0126 //      fEMfield = new PHG4FieldsPHENIX(fieldmapname,magfield_rescale);
0127 //      break;
0128 //    case 2:
0129 //      fEMfield = new PHG4Field2D(fieldmapname,0,magfield_rescale);
0130 //      break;
0131 //    case 3:
0132 //      fEMfield = new PHG4Field3D(fieldmapname,0,magfield_rescale);
0133 //      break;
0134 //    default:
0135 //      std::cout << "Invalid dimension, valid is 2 for 2D, 3 for 3D" << std::endl;
0136 //      exit(1);
0137 //    }
0138 //  fFieldMessenger = new G4TBFieldMessenger(this) ;
0139 //  fEquation = new  G4Mag_UsualEqRhs(fEMfield);
0140 //  fMinStep     = 0.005*mm ; // minimal step of 10 microns
0141 //  fStepperType = 4 ;        // ClassicalRK4 -- the default stepper
0142 //
0143 //  fFieldManager = GetGlobalFieldManager();
0144 //  UpdateField();
0145 //  double point[4] = {0,0,0,0};
0146 //  fEMfield->GetFieldValue(&point[0],&magfield_at_000[0]);
0147 //  for (size_t i=0; i<sizeof(magfield_at_000)/sizeof(double);i++)
0148 //    {
0149 //      magfield_at_000[i] = magfield_at_000[i]/tesla;
0150 //    }
0151 //  if (verbosity > 0)
0152 //    {
0153 //      std::cout << "field: x" << magfield_at_000[0]
0154 //     << ", y: " << magfield_at_000[1]
0155 //     << ", z: " << magfield_at_000[2]
0156 //     << std::endl;
0157 //    }
0158 //}
0159
0160 ////////////////////////////////////////////////////////////////////////////////
0161
0162 G4TBMagneticFieldSetup::~G4TBMagneticFieldSetup()
0163 {
0164   delete fChordFinder;
0165   delete fStepper;
0166   delete fFieldMessenger;
0167   delete fEquation;
0168   delete fEMfield;
0169 }
0170
0171 /////////////////////////////////////////////////////////////////////////////
0172 //
0173 // Register this field to 'global' Field Manager and
0174 // Create Stepper and Chord Finder with predefined type, minstep (resp.)
0175 //
0176
0177 void G4TBMagneticFieldSetup::UpdateField()
0178 {
0179   SetStepper();
0180
0181   fFieldManager->SetDetectorField(fEMfield);
0182
0183   delete fChordFinder;
0184
0185   fIntgrDriver = new G4MagInt_Driver(fMinStep,
0186                                      fStepper,
0187                                      fStepper->GetNumberOfVariables());
0188
0189   fChordFinder = new G4ChordFinder(fIntgrDriver);
0190
0191   fFieldManager->SetChordFinder(fChordFinder);
0192 }
0193
0194 /////////////////////////////////////////////////////////////////////////////
0195 //
0196 // Set stepper according to the stepper type
0197 //
0198
0199 void G4TBMagneticFieldSetup::SetStepper()
0200 {
0201   G4int nvar = 8;
0202
0203   delete fStepper;
0204
0205   std::stringstream message;
0206
0207   switch (fStepperType)
0208   {
0209   case 0:
0210     fStepper = new G4ExplicitEuler(fEquation, nvar);
0211     message << "Stepper in use: G4ExplicitEuler";
0212     break;
0213   case 1:
0214     fStepper = new G4ImplicitEuler(fEquation, nvar);
0215     message << "Stepper in use: G4ImplicitEuler";
0216     break;
0217   case 2:
0218     fStepper = new G4SimpleRunge(fEquation, nvar);
0219     message << "Stepper in use: G4SimpleRunge";
0220     break;
0221   case 3:
0222     fStepper = new G4SimpleHeum(fEquation, nvar);
0223     message << "Stepper in use: G4SimpleHeum";
0224     break;
0225   case 4:
0226     fStepper = new G4ClassicalRK4(fEquation, nvar);
0227     message << "Stepper in use: G4ClassicalRK4 (default)";
0228     break;
0229   case 5:
0230     fStepper = new G4CashKarpRKF45(fEquation, nvar);
0231     message << "Stepper in use: G4CashKarpRKF45";
0232     break;
0233   case 6:
0234     fStepper = nullptr;  // new G4RKG3_Stepper( fEquation, nvar );
0235     message << "G4RKG3_Stepper is not currently working for Magnetic Field";
0236     break;
0237   case 7:
0238     fStepper = nullptr;  // new G4HelixExplicitEuler( fEquation );
0239     message << "G4HelixExplicitEuler is not valid for Magnetic Field";
0240     break;
0241   case 8:
0242     fStepper = nullptr;  // new G4HelixImplicitEuler( fEquation );
0243     message << "G4HelixImplicitEuler is not valid for Magnetic Field";
0244     break;
0245   case 9:
0246     fStepper = nullptr;  // new G4HelixSimpleRunge( fEquation );
0247     message << "G4HelixSimpleRunge is not valid for Magnetic Field";
0248     break;
0249   default:
0250     fStepper = nullptr;
0251   }
0252
0253   if (verbosity > 0)
0254   {
0255     std::cout << " ---------- G4TBMagneticFieldSetup::SetStepper() -----------" << std::endl;
0256     std::cout << "  " << message.str() << std::endl;
0257     std::cout << "  Minimum step size: " << fMinStep / mm << " mm" << std::endl;
0258     std::cout << " -----------------------------------------------------------" << std::endl;
0259   }
0260
0261   if (!fStepper)
0262   {
0263     std::cout << "no stepper set, edxiting now" << std::endl;
0264     exit(1);
0265   }
0266
0267   return;
0268 }
0269
0270 /////////////////////////////////////////////////////////////////////////////
0271 //
0272 // Set the value of the Global Field to fieldValue along Z
0273 //
0274
0275 void G4TBMagneticFieldSetup::SetFieldValue(const G4double fieldValue)
0276 {
0277   G4ThreeVector fieldVector(0.0, 0.0, fieldValue);
0278
0279   SetFieldValue(fieldVector);
0280 }
0281
0282 ///////////////////////////////////////////////////////////////////////////////
0283 //
0284 // Set the value of the Global Field value to fieldVector
0285 //
0286
0287 void G4TBMagneticFieldSetup::SetFieldValue(const G4ThreeVector& fieldVector)
0288 {
0289   // Find the Field Manager for the global field
0290   G4FieldManager* fieldMgr = GetGlobalFieldManager();
0291
0292   if (fieldVector != G4ThreeVector(0., 0., 0.))
0293   {
0294     delete fEMfield;
0295
0296     fEMfield = new G4UniformMagField(fieldVector);
0297
0298     fEquation->SetFieldObj(fEMfield);  // must now point to the new field
0299
0300     // UpdateField();
0301
0302     fieldMgr->SetDetectorField(fEMfield);
0303   }
0304   else
0305   {
0306     // If the new field's value is Zero, then it is best to
0307     //  insure that it is not used for propagation.
0308     delete fEMfield;
0309     fEMfield = nullptr;
0310     fEquation->SetFieldObj(fEMfield);  // As a double check ...
0311     fieldMgr->SetDetectorField(fEMfield);
0312   }
0313 }
0314
0315 ////////////////////////////////////////////////////////////////////////////////
0316 //
0317 //  Utility method
0318
0319 G4FieldManager* G4TBMagneticFieldSetup::GetGlobalFieldManager()
0320 {
0321   return G4TransportationManager::GetTransportationManager()->GetFieldManager();
0322 }