ECCE @ EIC 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 
61 using namespace std;
62 
64 //
65 // Constructors:
66 
68  : verbosity(0)
69  , fChordFinder(0)
70  , fStepper(0)
71  , fIntgrDriver(0)
72 {
73  assert(phfield);
74 
75  fEMfield = new PHG4MagneticField(phfield);
78  fMinStep = 0.005 * mm; // minimal step of 10 microns
79  fStepperType = 4; // ClassicalRK4 -- the default stepper
80 
82  UpdateField();
83  double point[4] = {0, 0, 0, 0};
84  fEMfield->GetFieldValue(&point[0], &magfield_at_000[0]);
85  for (size_t i = 0; i < sizeof(magfield_at_000) / sizeof(double); i++)
86  {
88  }
89  if (verbosity > 0)
90  {
91  cout << "field: x" << magfield_at_000[0]
92  << ", y: " << magfield_at_000[1]
93  << ", z: " << magfield_at_000[2]
94  << endl;
95  }
96 }
97 
98 //G4TBMagneticFieldSetup::G4TBMagneticFieldSetup(const float magfield)
99 // : verbosity(0), fChordFinder(0), fStepper(0), fIntgrDriver(0)
100 //{
101 // //solenoidal field along the axis of the cyclinders?
102 // fEMfield = new G4UniformMagField(G4ThreeVector(0.0, 0.0, magfield*tesla));
103 // fFieldMessenger = new G4TBFieldMessenger(this) ;
104 // fEquation = new G4Mag_UsualEqRhs(fEMfield);
105 // fMinStep = 0.005 * mm ; // minimal step of 10 microns
106 // fStepperType = 4 ; // ClassicalRK4 -- the default stepper
107 //
108 // fFieldManager = GetGlobalFieldManager();
109 // UpdateField();
110 //
111 // double point[4] = {0,0,0,0};
112 // fEMfield->GetFieldValue(&point[0],&magfield_at_000[0]);
113 // for (size_t i=0; i<sizeof(magfield_at_000)/sizeof(double);i++)
114 // {
115 // magfield_at_000[i] = magfield_at_000[i]/tesla;
116 // }
117 //}
118 //
120 //
121 //G4TBMagneticFieldSetup::G4TBMagneticFieldSetup(const string &fieldmapname, const int dim, const float magfield_rescale)
122 // : verbosity(0),
123 // fChordFinder(0),
124 // fStepper(0),
125 // fIntgrDriver(0)
126 //{
127 //
128 // switch(dim)
129 // {
130 // case 1:
131 // fEMfield = new PHG4FieldsPHENIX(fieldmapname,magfield_rescale);
132 // break;
133 // case 2:
134 // fEMfield = new PHG4Field2D(fieldmapname,0,magfield_rescale);
135 // break;
136 // case 3:
137 // fEMfield = new PHG4Field3D(fieldmapname,0,magfield_rescale);
138 // break;
139 // default:
140 // cout << "Invalid dimension, valid is 2 for 2D, 3 for 3D" << endl;
141 // exit(1);
142 // }
143 // fFieldMessenger = new G4TBFieldMessenger(this) ;
144 // fEquation = new G4Mag_UsualEqRhs(fEMfield);
145 // fMinStep = 0.005*mm ; // minimal step of 10 microns
146 // fStepperType = 4 ; // ClassicalRK4 -- the default stepper
147 //
148 // fFieldManager = GetGlobalFieldManager();
149 // UpdateField();
150 // double point[4] = {0,0,0,0};
151 // fEMfield->GetFieldValue(&point[0],&magfield_at_000[0]);
152 // for (size_t i=0; i<sizeof(magfield_at_000)/sizeof(double);i++)
153 // {
154 // magfield_at_000[i] = magfield_at_000[i]/tesla;
155 // }
156 // if (verbosity > 0)
157 // {
158 // cout << "field: x" << magfield_at_000[0]
159 // << ", y: " << magfield_at_000[1]
160 // << ", z: " << magfield_at_000[2]
161 // << endl;
162 // }
163 //}
164 
166 
168 {
169  delete fChordFinder;
170  delete fStepper;
171  delete fFieldMessenger;
172  delete fEquation;
173  delete fEMfield;
174 }
175 
177 //
178 // Register this field to 'global' Field Manager and
179 // Create Stepper and Chord Finder with predefined type, minstep (resp.)
180 //
181 
183 {
184  SetStepper();
185 
187 
188  delete fChordFinder;
189 
191  fStepper,
193 
195 
196  fFieldManager->SetChordFinder(fChordFinder);
197 }
198 
200 //
201 // Set stepper according to the stepper type
202 //
203 
205 {
206  G4int nvar = 8;
207 
208  delete fStepper;
209 
210  std::stringstream message;
211 
212  switch (fStepperType)
213  {
214  case 0:
215  fStepper = new G4ExplicitEuler(fEquation, nvar);
216  message << "Stepper in use: G4ExplicitEuler";
217  break;
218  case 1:
219  fStepper = new G4ImplicitEuler(fEquation, nvar);
220  message << "Stepper in use: G4ImplicitEuler";
221  break;
222  case 2:
223  fStepper = new G4SimpleRunge(fEquation, nvar);
224  message << "Stepper in use: G4SimpleRunge";
225  break;
226  case 3:
227  fStepper = new G4SimpleHeum(fEquation, nvar);
228  message << "Stepper in use: G4SimpleHeum";
229  break;
230  case 4:
231  fStepper = new G4ClassicalRK4(fEquation, nvar);
232  message << "Stepper in use: G4ClassicalRK4 (default)";
233  break;
234  case 5:
235  fStepper = new G4CashKarpRKF45(fEquation, nvar);
236  message << "Stepper in use: G4CashKarpRKF45";
237  break;
238  case 6:
239  fStepper = nullptr; // new G4RKG3_Stepper( fEquation, nvar );
240  message << "G4RKG3_Stepper is not currently working for Magnetic Field";
241  break;
242  case 7:
243  fStepper = nullptr; // new G4HelixExplicitEuler( fEquation );
244  message << "G4HelixExplicitEuler is not valid for Magnetic Field";
245  break;
246  case 8:
247  fStepper = nullptr; // new G4HelixImplicitEuler( fEquation );
248  message << "G4HelixImplicitEuler is not valid for Magnetic Field";
249  break;
250  case 9:
251  fStepper = nullptr; // new G4HelixSimpleRunge( fEquation );
252  message << "G4HelixSimpleRunge is not valid for Magnetic Field";
253  break;
254  default:
255  fStepper = nullptr;
256  }
257 
258  if (verbosity > 0)
259  {
260  std::cout << " ---------- G4TBMagneticFieldSetup::SetStepper() -----------" << std::endl;
261  std::cout << " " << message.str() << endl;
262  std::cout << " Minimum step size: " << fMinStep / mm << " mm" << std::endl;
263  std::cout << " -----------------------------------------------------------" << std::endl;
264  }
265 
266  if (!fStepper)
267  {
268  cout << "no stepper set, edxiting now" << endl;
269  exit(1);
270  }
271 
272  return;
273 }
274 
276 //
277 // Set the value of the Global Field to fieldValue along Z
278 //
279 
281 {
282  G4ThreeVector fieldVector(0.0, 0.0, fieldValue);
283 
284  SetFieldValue(fieldVector);
285 }
286 
288 //
289 // Set the value of the Global Field value to fieldVector
290 //
291 
293 {
294  // Find the Field Manager for the global field
296 
297  if (fieldVector != G4ThreeVector(0., 0., 0.))
298  {
299  if (fEMfield) delete fEMfield;
300  fEMfield = new G4UniformMagField(fieldVector);
301 
302  fEquation->SetFieldObj(fEMfield); // must now point to the new field
303 
304  // UpdateField();
305 
306  fieldMgr->SetDetectorField(fEMfield);
307  }
308  else
309  {
310  // If the new field's value is Zero, then it is best to
311  // insure that it is not used for propagation.
312  delete fEMfield;
313  fEMfield = nullptr;
314  fEquation->SetFieldObj(fEMfield); // As a double check ...
315  fieldMgr->SetDetectorField(fEMfield);
316  }
317 }
318 
320 //
321 // Utility method
322 
324 {
326 }