Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4TBMagneticFieldSetup.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4TBMagneticFieldSetup.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // $Id: G4TBMagneticFieldSetup.cc,v 1.24 2014/11/14 21:47:38 mccumber Exp $
28 // GEANT4 tag $Name: $
29 //
30 //
31 // User Field class implementation.
32 //
33 
35 #include "G4TBFieldMessenger.hh"
36 #include "PHG4MagneticField.h"
37 
38 #include <Geant4/G4CashKarpRKF45.hh>
39 #include <Geant4/G4ChordFinder.hh>
40 #include <Geant4/G4ClassicalRK4.hh>
41 #include <Geant4/G4ExplicitEuler.hh>
42 #include <Geant4/G4FieldManager.hh>
43 #include <Geant4/G4ImplicitEuler.hh>
44 #include <Geant4/G4MagIntegratorDriver.hh>
45 #include <Geant4/G4MagIntegratorStepper.hh>
46 #include <Geant4/G4Mag_UsualEqRhs.hh>
47 #include <Geant4/G4MagneticField.hh>
48 #include <Geant4/G4SimpleHeum.hh>
49 #include <Geant4/G4SimpleRunge.hh>
50 #include <Geant4/G4SystemOfUnits.hh>
51 #include <Geant4/G4ThreeVector.hh>
52 #include <Geant4/G4TransportationManager.hh>
53 #include <Geant4/G4Types.hh> // for G4double, G4int
54 #include <Geant4/G4UniformMagField.hh>
55 
56 #include <cassert>
57 #include <cstdlib> // for exit, size_t
58 #include <iostream>
59 #include <string>
60 
62 //
63 // Constructors:
64 
66 {
67  assert(phfield);
68 
69  fEMfield = new PHG4MagneticField(phfield);
71  fEquation = new G4Mag_UsualEqRhs(fEMfield);
72  fMinStep = 0.005 * mm; // minimal step of 10 microns
73  fStepperType = 4; // ClassicalRK4 -- the default stepper
74 
76  UpdateField();
77  double point[4] = {0, 0, 0, 0};
78  fEMfield->GetFieldValue(&point[0], &magfield_at_000[0]);
79  for (size_t i = 0; i < sizeof(magfield_at_000) / sizeof(double); i++)
80  {
81  magfield_at_000[i] = magfield_at_000[i] / tesla;
82  }
83  if (verbosity > 0)
84  {
85  std::cout << "field: x" << magfield_at_000[0]
86  << ", y: " << magfield_at_000[1]
87  << ", z: " << magfield_at_000[2]
88  << std::endl;
89  }
90 }
91 
92 //G4TBMagneticFieldSetup::G4TBMagneticFieldSetup(const float magfield)
93 // : verbosity(0), fChordFinder(0), fStepper(0), fIntgrDriver(0)
94 //{
95 // //solenoidal field along the axis of the cyclinders?
96 // fEMfield = new G4UniformMagField(G4ThreeVector(0.0, 0.0, magfield*tesla));
97 // fFieldMessenger = new G4TBFieldMessenger(this) ;
98 // fEquation = new G4Mag_UsualEqRhs(fEMfield);
99 // fMinStep = 0.005 * mm ; // minimal step of 10 microns
100 // fStepperType = 4 ; // ClassicalRK4 -- the default stepper
101 //
102 // fFieldManager = GetGlobalFieldManager();
103 // UpdateField();
104 //
105 // double point[4] = {0,0,0,0};
106 // fEMfield->GetFieldValue(&point[0],&magfield_at_000[0]);
107 // for (size_t i=0; i<sizeof(magfield_at_000)/sizeof(double);i++)
108 // {
109 // magfield_at_000[i] = magfield_at_000[i]/tesla;
110 // }
111 //}
112 //
114 //
115 //G4TBMagneticFieldSetup::G4TBMagneticFieldSetup(const string &fieldmapname, const int dim, const float magfield_rescale)
116 // : verbosity(0),
117 // fChordFinder(0),
118 // fStepper(0),
119 // fIntgrDriver(0)
120 //{
121 //
122 // switch(dim)
123 // {
124 // case 1:
125 // fEMfield = new PHG4FieldsPHENIX(fieldmapname,magfield_rescale);
126 // break;
127 // case 2:
128 // fEMfield = new PHG4Field2D(fieldmapname,0,magfield_rescale);
129 // break;
130 // case 3:
131 // fEMfield = new PHG4Field3D(fieldmapname,0,magfield_rescale);
132 // break;
133 // default:
134 // std::cout << "Invalid dimension, valid is 2 for 2D, 3 for 3D" << std::endl;
135 // exit(1);
136 // }
137 // fFieldMessenger = new G4TBFieldMessenger(this) ;
138 // fEquation = new G4Mag_UsualEqRhs(fEMfield);
139 // fMinStep = 0.005*mm ; // minimal step of 10 microns
140 // fStepperType = 4 ; // ClassicalRK4 -- the default stepper
141 //
142 // fFieldManager = GetGlobalFieldManager();
143 // UpdateField();
144 // double point[4] = {0,0,0,0};
145 // fEMfield->GetFieldValue(&point[0],&magfield_at_000[0]);
146 // for (size_t i=0; i<sizeof(magfield_at_000)/sizeof(double);i++)
147 // {
148 // magfield_at_000[i] = magfield_at_000[i]/tesla;
149 // }
150 // if (verbosity > 0)
151 // {
152 // std::cout << "field: x" << magfield_at_000[0]
153 // << ", y: " << magfield_at_000[1]
154 // << ", z: " << magfield_at_000[2]
155 // << std::endl;
156 // }
157 //}
158 
160 
162 {
163  delete fChordFinder;
164  delete fStepper;
165  delete fFieldMessenger;
166  delete fEquation;
167  delete fEMfield;
168 }
169 
171 //
172 // Register this field to 'global' Field Manager and
173 // Create Stepper and Chord Finder with predefined type, minstep (resp.)
174 //
175 
177 {
178  SetStepper();
179 
180  fFieldManager->SetDetectorField(fEMfield);
181 
182  delete fChordFinder;
183 
184  fIntgrDriver = new G4MagInt_Driver(fMinStep,
185  fStepper,
186  fStepper->GetNumberOfVariables());
187 
188  fChordFinder = new G4ChordFinder(fIntgrDriver);
189 
190  fFieldManager->SetChordFinder(fChordFinder);
191 }
192 
194 //
195 // Set stepper according to the stepper type
196 //
197 
199 {
200  G4int nvar = 8;
201 
202  delete fStepper;
203 
204  std::stringstream message;
205 
206  switch (fStepperType)
207  {
208  case 0:
209  fStepper = new G4ExplicitEuler(fEquation, nvar);
210  message << "Stepper in use: G4ExplicitEuler";
211  break;
212  case 1:
213  fStepper = new G4ImplicitEuler(fEquation, nvar);
214  message << "Stepper in use: G4ImplicitEuler";
215  break;
216  case 2:
217  fStepper = new G4SimpleRunge(fEquation, nvar);
218  message << "Stepper in use: G4SimpleRunge";
219  break;
220  case 3:
221  fStepper = new G4SimpleHeum(fEquation, nvar);
222  message << "Stepper in use: G4SimpleHeum";
223  break;
224  case 4:
225  fStepper = new G4ClassicalRK4(fEquation, nvar);
226  message << "Stepper in use: G4ClassicalRK4 (default)";
227  break;
228  case 5:
229  fStepper = new G4CashKarpRKF45(fEquation, nvar);
230  message << "Stepper in use: G4CashKarpRKF45";
231  break;
232  case 6:
233  fStepper = nullptr; // new G4RKG3_Stepper( fEquation, nvar );
234  message << "G4RKG3_Stepper is not currently working for Magnetic Field";
235  break;
236  case 7:
237  fStepper = nullptr; // new G4HelixExplicitEuler( fEquation );
238  message << "G4HelixExplicitEuler is not valid for Magnetic Field";
239  break;
240  case 8:
241  fStepper = nullptr; // new G4HelixImplicitEuler( fEquation );
242  message << "G4HelixImplicitEuler is not valid for Magnetic Field";
243  break;
244  case 9:
245  fStepper = nullptr; // new G4HelixSimpleRunge( fEquation );
246  message << "G4HelixSimpleRunge is not valid for Magnetic Field";
247  break;
248  default:
249  fStepper = nullptr;
250  }
251 
252  if (verbosity > 0)
253  {
254  std::cout << " ---------- G4TBMagneticFieldSetup::SetStepper() -----------" << std::endl;
255  std::cout << " " << message.str() << std::endl;
256  std::cout << " Minimum step size: " << fMinStep / mm << " mm" << std::endl;
257  std::cout << " -----------------------------------------------------------" << std::endl;
258  }
259 
260  if (!fStepper)
261  {
262  std::cout << "no stepper set, edxiting now" << std::endl;
263  exit(1);
264  }
265 
266  return;
267 }
268 
270 //
271 // Set the value of the Global Field to fieldValue along Z
272 //
273 
274 void G4TBMagneticFieldSetup::SetFieldValue(const G4double fieldValue)
275 {
276  G4ThreeVector fieldVector(0.0, 0.0, fieldValue);
277 
278  SetFieldValue(fieldVector);
279 }
280 
282 //
283 // Set the value of the Global Field value to fieldVector
284 //
285 
286 void G4TBMagneticFieldSetup::SetFieldValue(const G4ThreeVector fieldVector)
287 {
288  // Find the Field Manager for the global field
289  G4FieldManager* fieldMgr = GetGlobalFieldManager();
290 
291  if (fieldVector != G4ThreeVector(0., 0., 0.))
292  {
293  if (fEMfield) delete fEMfield;
294  fEMfield = new G4UniformMagField(fieldVector);
295 
296  fEquation->SetFieldObj(fEMfield); // must now point to the new field
297 
298  // UpdateField();
299 
300  fieldMgr->SetDetectorField(fEMfield);
301  }
302  else
303  {
304  // If the new field's value is Zero, then it is best to
305  // insure that it is not used for propagation.
306  delete fEMfield;
307  fEMfield = nullptr;
308  fEquation->SetFieldObj(fEMfield); // As a double check ...
309  fieldMgr->SetDetectorField(fEMfield);
310  }
311 }
312 
314 //
315 // Utility method
316 
318 {
319  return G4TransportationManager::GetTransportationManager()->GetFieldManager();
320 }