source: trunk/source/geometry/magneticfield/test/field06/exampleS5.cc@ 1199

Last change on this file since 1199 was 1199, checked in by garnier, 16 years ago

nvx fichiers dans CVS
File size: 11.2 KB
Line 
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: exampleS5.cc,v 1.2 2007/05/02 14:59:26 gunter Exp $
28// GEANT4 tag $Name: HEAD $
29//
30//
31// --------------------------------------------------------------
32// GEANT 4 - exampleN01
33// --------------------------------------------------------------
34
35#include "G4RunManager.hh"
36#include "G4UImanager.hh"
37
38#include "ExN01DetectorConstruction.hh"
39#include "ExN01PhysicsList.hh"
40#include "ExN01PrimaryGeneratorAction.hh"
41#include "ExN01RunAction.hh"
42#include "ExN01EventAction.hh"
43#include "ExN01SteppingAction.hh"
44
45#include "G4UIterminal.hh"
46#include "G4UItcsh.hh"
47#include "G4UniformMagField.hh"
48#include "G4VPhysicalVolume.hh"
49#include "G4ChordFinder.hh"
50#include "G4PropagatorInField.hh"
51#include "G4MagneticField.hh"
52#include "G4FieldManager.hh"
53#include "G4TransportationManager.hh"
54#include "G4Navigator.hh"
55#include "G4HelixExplicitEuler.hh"
56#include "G4HelixSimpleRunge.hh"
57#include "G4HelixImplicitEuler.hh"
58#include "G4ExactHelixStepper.hh"
59#include "G4ExplicitEuler.hh"
60#include "G4ImplicitEuler.hh"
61#include "G4SimpleRunge.hh"
62#include "G4SimpleHeum.hh"
63#include "G4ClassicalRK4.hh"
64#include "G4Mag_UsualEqRhs.hh"
65#include "G4MagIntegratorStepper.hh"
66#include "G4CashKarpRKF45.hh"
67#include "G4RKG3_Stepper.hh"
68#include "G4Timer.hh"
69#include "G4ios.hh"
70#include "G4QuadrupoleMagField.hh"
71#include "G4MagIntegratorDriver.hh"
72#include <fstream>
73///
74///MagField
75///
76
77///
78///EndMagField
79///
80
81///
82///SetupStepper
83///
84G4MagIntegratorStepper* SetupStepper(G4int type,G4Mag_UsualEqRhs* pE)
85{
86 G4MagIntegratorStepper* pStepper;
87 switch ( type )
88 {
89 case 0: pStepper = new G4ExplicitEuler( pE ); break;
90 case 1: pStepper = new G4ImplicitEuler( pE ); break;
91 case 2: pStepper = new G4SimpleRunge( pE ); break;
92 case 3: pStepper = new G4SimpleHeum( pE ); break;
93 case 4: pStepper = new G4ClassicalRK4( pE ); break;
94 case 5: pStepper = new G4HelixExplicitEuler( pE ); break;
95 case 6: pStepper = new G4HelixImplicitEuler( pE ); break;
96 case 7: pStepper = new G4HelixSimpleRunge( pE ); break;
97 case 8: pStepper = new G4CashKarpRKF45( pE ); break;
98 case 9: pStepper = new G4ExactHelixStepper( pE ); break;
99 case 10: pStepper = new G4RKG3_Stepper( pE ); break;
100 default: pStepper = 0;
101 }
102 return pStepper;
103}
104///
105///EndSetupStepper
106///
107
108///
109///CompareStepStepper
110///
111void CompareStep(G4double y[],G4double dydx[],G4double hhh,
112 G4double yout1[],G4double yout2[],G4double ydif[],
113 G4double yerr1[],G4double yerr2[],
114 G4MagIntegratorStepper* pSt1,G4MagIntegratorStepper* pSt2)
115{
116 //G4double yerr1[G4FieldTrack::ncompSVEC];
117 // G4double yerr2[G4FieldTrack::ncompSVEC];
118 G4double ydiferr[G4FieldTrack::ncompSVEC];
119 G4double r[10];
120 //
121 //r0=r;r1=Pr;r2=difR;r3=difPr;r4=diferrR;r5=diferrPr;r6=err1r;r7=err1Pr;
122 //r8=err2r;r9=err2Pr;
123 //
124 r[0]=std::sqrt(y[0]*y[0]+y[1]*y[1]+y[2]*y[2]);
125 r[1]=std::sqrt(y[3]*y[3]+y[4]*y[4]+y[5]*y[5]);
126 pSt1->Stepper(y,dydx,hhh,yout1,yerr1);
127 //G4cout << " One step with stepper type 1"<<G4endl;
128 //G4cout << " Step Length= "<<hhh/mm<<"mm"<<G4endl;
129 // G4cout << "Y DyDx Yout "<<G4endl;
130 //G4cout.precision(10);
131 //for (G4int i=0;i<6;i++)
132 //G4cout << y[i]<<" "<<dydx[i]<<" "<<yout1[i]<<G4endl;
133 pSt2->Stepper(y,dydx,hhh,yout2,yerr2);
134 //G4cout << " One step with stepper type 2 "<<G4endl;
135 //G4cout << "Y DyDx Yout "<<G4endl;
136 //G4cout.precision(10);
137 //for (G4int i=0;i<6;i++)
138 // G4cout << y[i]<<" "<<dydx[i]<<" "<<yout2[i]<<G4endl;
139 for (G4int i=0;i<6;i++) ydif[i]= std::abs(yout1[i]-yout2[i]);
140 for (G4int i=0;i<6;i++) ydiferr[i]=std::sqrt(yerr1[i]*yerr1[i]+yerr2[i]*yerr2[i]);
141 r[2]=std::sqrt(ydif[0]*ydif[0]+ydif[1]*ydif[1]+ydif[2]*ydif[2]);
142 r[3]=std::sqrt(ydif[3]*ydif[3]+ydif[4]*ydif[4]+ydif[5]*ydif[5]);
143 r[4]=std::sqrt(ydiferr[0]*ydiferr[0]+ydiferr[1]*ydiferr[1]+ydiferr[2]*ydiferr[2]);
144 r[5]=std::sqrt(ydiferr[3]*ydiferr[3]+ydiferr[4]*ydiferr[4]+ydiferr[5]*ydiferr[5]);
145 r[6]=std::sqrt(yerr1[0]*yerr1[0]+yerr1[1]*yerr1[1]+yerr1[2]*yerr1[2]);
146 r[7]=std::sqrt(yerr1[3]*yerr1[3]+yerr1[4]*yerr1[4]+yerr1[5]*yerr1[5]);
147 r[8]=std::sqrt(yerr2[0]*yerr2[0]+yerr2[1]*yerr2[1]+yerr2[2]*yerr2[2]);
148 r[9]=std::sqrt(yerr2[3]*yerr2[3]+yerr2[4]*yerr2[4]+yerr2[5]*yerr2[5]);
149 G4cout.precision(10);
150 for (G4int i=0;i<6;i++)
151 // G4cout << "deltaY["<<i<<"]="<<ydif[i]/mm<<"mm"<<" "
152 // <<ydiferr[i]/mm<<" "
153 // <<yerr1[i]/mm<<" "<<yerr2[i]/mm<<G4endl;
154 G4cout <<hhh<<" "<<y[i]<<" "<<ydif[i]/mm <<" "
155 <<ydiferr[i]/mm<<" "
156 <<yerr1[i]/mm<<" "<<yerr2[i]/mm<<G4endl;
157 // G4cout <<hhh<<" "<<y[i]<<" "<<ydif[i]/mm <<" "
158 // <<yout1[i]/mm<<" "
159 // <<yout2[i]/mm<<" "<<yerr2[i]/mm<<G4endl;
160}
161///
162///OneStepStepper
163///
164///
165void OneStep(G4double y[],G4double dydx[],G4double hhh,
166 G4double yout1[],G4double yerr1[],
167 G4MagIntegratorStepper* pSt1)
168{
169
170 pSt1->Stepper(y,dydx,hhh,yout1,yerr1);
171 //G4cout << " One step with stepper type 1"<<G4endl;
172 //G4cout << " Step Length= "<<hhh/mm<<"mm"<<G4endl;
173 // G4cout << "Y DyDx Yout "<<G4endl;
174 //G4cout.precision(10);
175 //for (G4int i=0;i<6;i++)
176 //G4cout << y[i]<<" "<<dydx[i]<<" "<<yout1[i]<<G4endl;
177
178}
179///
180///OneStepStepper
181///
182 #ifdef G4VIS_USE
183 #include "G4VisExecutive.hh"
184 #endif
185
186int main(int argc,char** argv)
187{
188 // Construct the default run manager
189 //
190 G4RunManager* runManager = new G4RunManager;
191
192 // set mandatory initialization classes
193 //
194
195 ExN01DetectorConstruction* detector = new ExN01DetectorConstruction;
196 runManager->SetUserInitialization(detector);
197 //
198 G4VUserPhysicsList* physics = new ExN01PhysicsList;
199 runManager->SetUserInitialization(physics);
200
201 #ifdef G4VIS_USE
202 // Visualization, if you choose to have it!
203 //
204 G4VisManager* visManager = new G4VisExecutive;
205 visManager->Initialize();
206 #endif
207
208 // set mandatory user action class
209 //
210 G4VUserPrimaryGeneratorAction* gen_action = new ExN01PrimaryGeneratorAction;
211 runManager->SetUserAction(gen_action);
212
213
214 ExN01RunAction* run_action = new ExN01RunAction;
215 runManager->SetUserAction(run_action);
216 //
217 ExN01EventAction* event_action = new ExN01EventAction(run_action);
218 runManager->SetUserAction(event_action);
219 //
220 G4UserSteppingAction* stepping_action =
221 new ExN01SteppingAction(detector,event_action);
222 runManager->SetUserAction(stepping_action);
223 //
224 //Mag Field
225 //
226 G4int type=4;
227 G4ThreeVector Bconst(0.,0.,-2.*tesla);
228
229 G4UniformMagField* pmagField;
230 G4FieldManager* pfieldMgr;
231 G4MagIntegratorStepper* pStepper1;
232 G4MagIntegratorStepper* pStepper2;
233 G4ChordFinder* pChordFinder;
234 G4Mag_UsualEqRhs* pEquation;
235
236 pmagField = new G4UniformMagField(Bconst);
237
238 //// Get the global field manager
239 pfieldMgr= G4TransportationManager::GetTransportationManager()->
240 GetFieldManager();
241 // // Set this field to the global field manager
242 pfieldMgr->SetDetectorField(pmagField );
243 // //stepper choise
244 pEquation=new G4Mag_UsualEqRhs(pmagField);
245 pStepper1=SetupStepper(type,pEquation);
246 pChordFinder = new G4ChordFinder(pmagField,
247 1.0e-2 * mm,
248 pStepper1);
249 pfieldMgr->SetChordFinder( pChordFinder );
250
251 //
252 //Propagator
253 //
254 G4PropagatorInField* pPropagator;
255 pPropagator = G4TransportationManager::GetTransportationManager()->
256 GetPropagatorInField ();
257 pPropagator -> SetMinimumEpsilonStep( 1.0e-5*mm ) ;
258 pPropagator -> SetMaximumEpsilonStep( 1.0e-5*mm ) ;
259 G4cout << " Using values for "
260 << " Min Eps = " << pPropagator->GetMinimumEpsilonStep()
261 << " and "
262 << " MaxEps = " << pPropagator->GetMaximumEpsilonStep()
263 << G4endl;
264 pPropagator ->SetLargestAcceptableStep(5.0*cm);
265
266 //
267 //End MagField
268
269 //
270
271 // Initialize G4 kernel
272 //
273 runManager->Initialize();
274// Get the pointer to the User Interface manager
275 //
276 G4UImanager * UI = G4UImanager::GetUIpointer();
277
278 if(argc==1) // Define (G)UI terminal for interactive mode
279 {
280 // G4UIterminal is a (dumb) terminal
281 //
282 G4UIsession * session = 0;
283 #ifdef G4UI_USE_TCSH
284 session = new G4UIterminal(new G4UItcsh);
285 #else
286 session = new G4UIterminal();
287 #endif
288 UI->ApplyCommand("/control/execute vis.mac");
289
290 // UI->ApplyCommand("/gun/energy 1 GeV");
291 G4int numberOfEvent = 1;
292 runManager->BeamOn(numberOfEvent);
293 UI->ApplyCommand("/run/verbose 1");
294
295 session->SessionStart();
296
297 delete session;
298 }
299 else // Batch mode
300 {
301 G4String command = "/control/execute ";
302 G4String fileName = argv[1];
303 UI->ApplyCommand(command+fileName);
304 }
305
306 // Free the store: user actions, physics_list and detector_description are
307 // owned and deleted by the run manager, so they should not
308 // be deleted in the main() program !
309
310 // Get the pointer to the UI manager and set verbosities
311 //
312 //G4UImanager* UI = G4UImanager::GetUIpointer();
313 //UI->ApplyCommand("/run/verbose 1");
314 //UI->ApplyCommand("/event/verbose 1");
315 // UI->ApplyCommand("/tracking/verbose 1");
316
317 // Start a run
318 //
319 //G4int numberOfEvent = 3;
320 //runManager->BeamOn(numberOfEvent);
321
322 // Job termination
323 //
324 // Free the store: user actions, physics_list and detector_description are
325 // owned and deleted by the run manager, so they should not
326 // be deleted in the main() program !
327 //
328#ifdef G4VIS_USE
329 delete visManager;
330 #endif
331 delete runManager;
332
333 return 0;
334}
335
336
Note: See TracBrowser for help on using the repository browser.