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