| Geant4 Cross Reference |
>> 1 // This code implementation is the intellectual property of
>> 2 // the GEANT4 collaboration.
1 // 3 //
2 // ******************************************* << 4 // By copying, distributing or modifying the Program (or any work
3 // * License and Disclaimer << 5 // based on the Program) you indicate your acceptance of this statement,
4 // * << 6 // and all its terms.
5 // * The Geant4 software is copyright of th <<
6 // * the Geant4 Collaboration. It is provided <<
7 // * conditions of the Geant4 Software License <<
8 // * LICENSE and available at http://cern.ch/ <<
9 // * include a list of copyright holders. <<
10 // * <<
11 // * Neither the authors of this software syst <<
12 // * institutes,nor the agencies providing fin <<
13 // * work make any representation or warran <<
14 // * regarding this software system or assum <<
15 // * use. Please see the license in the file <<
16 // * for the full disclaimer and the limitatio <<
17 // * <<
18 // * This code implementation is the result <<
19 // * technical work of the GEANT4 collaboratio <<
20 // * By using, copying, modifying or distri <<
21 // * any work based on the software) you ag <<
22 // * use in resulting scientific publicati <<
23 // * acceptance of all terms of the Geant4 Sof <<
24 // ******************************************* <<
25 // 7 //
26 // G4ChordFinder implementation << 8 // $Id: G4ChordFinder.cc,v 1.18 2001/03/23 18:50:33 japost Exp $
>> 9 // GEANT4 tag $Name: geant4-03-01 $
27 // 10 //
28 // Author: J.Apostolakis - Design and implemen << 11 //
29 // ------------------------------------------- << 12 // 25.02.97 John Apostolakis, design and implimentation
30 << 13 // 05.03.97 V. Grichine , style modification
31 #include <iomanip> <<
32 14
33 #include "G4ChordFinder.hh" 15 #include "G4ChordFinder.hh"
34 #include "G4SystemOfUnits.hh" << 16 #include "G4MagIntegratorDriver.hh"
35 #include "G4MagneticField.hh" <<
36 #include "G4Mag_UsualEqRhs.hh" 17 #include "G4Mag_UsualEqRhs.hh"
>> 18 #include "G4ClassicalRK4.hh"
>> 19 // #include "G4Field.hh"
>> 20 // #include "G4MagIntegratorStepper.hh"
37 #include "G4MagIntegratorDriver.hh" 21 #include "G4MagIntegratorDriver.hh"
38 // #include "G4ClassicalRK4.hh" << 22 #include "g4std/iomanip"
39 // #include "G4CashKarpRKF45.hh" <<
40 // #include "G4NystromRK4.hh" <<
41 // #include "G4BogackiShampine23.hh" <<
42 // #include "G4BogackiShampine45.hh" <<
43 <<
44 #include "G4DormandPrince745.hh" <<
45 23
46 // New templated stepper(s) -- avoid virtual c << 24 // For the moment fDeltaChord is a constant!
47 #include "G4TDormandPrince45.hh" <<
48 25
49 // FSAL type driver / steppers ----- << 26 const G4double G4ChordFinder::fDefaultDeltaChord = 3. * mm;
50 #include "G4FSALIntegrationDriver.hh" <<
51 #include "G4VFSALIntegrationStepper.hh" <<
52 #include "G4RK547FEq1.hh" <<
53 // #include "G4RK547FEq2.hh" <<
54 // #include "G4RK547FEq3.hh" <<
55 // #include "G4FSALBogackiShampine45.hh" <<
56 // #include "G4FSALDormandPrince745.hh" <<
57 27
58 // Templated type drivers ----- << 28 // ..........................................................................
59 #include "G4IntegrationDriver.hh" <<
60 #include "G4InterpolationDriver.hh" <<
61 29
62 #include "G4HelixHeum.hh" << 30 G4ChordFinder::G4ChordFinder( G4MagneticField* theMagField,
63 #include "G4BFieldIntegrationDriver.hh" << 31 G4double stepMinimum,
>> 32 G4MagIntegratorStepper* pItsStepper ) // A default one
>> 33 : fDeltaChord( fDefaultDeltaChord )
>> 34 {
>> 35 // Construct the Chord Finder
>> 36 // by creating in inverse order the Driver, the Stepper and EqRhs ...
>> 37 G4Mag_EqRhs *pEquation = new G4Mag_UsualEqRhs(theMagField);
>> 38 fEquation = pEquation;
>> 39 fLastStepEstimate_Unconstrained = DBL_MAX; // Should move q, p to
>> 40 // G4FieldTrack ??
>> 41 // --->> Charge Q = 0
>> 42 // --->> Momentum P = 1 NOMINAL VALUES !!!!!!!!!!!!!!!!!!
64 43
65 #include "G4QSSDriverCreator.hh" << 44 if( pItsStepper == 0 )
>> 45 {
>> 46 pItsStepper = fDriversStepper = new G4ClassicalRK4(pEquation);
>> 47 fAllocatedStepper= true;
>> 48 }
>> 49 else
>> 50 {
>> 51 fAllocatedStepper= false;
>> 52 }
>> 53 fIntgrDriver = new G4MagInt_Driver(stepMinimum,
>> 54 pItsStepper,
>> 55 pItsStepper->GetNumberOfVariables() );
>> 56 }
66 57
67 #include "G4CachedMagneticField.hh" << 58 // ......................................................................
68 59
69 #include <cassert> << 60 G4ChordFinder::~G4ChordFinder()
70 #include <memory> << 61 {
>> 62 delete fEquation; // fIntgrDriver->pIntStepper->theEquation_Rhs;
>> 63 if( fAllocatedStepper)
>> 64 {
>> 65 delete fDriversStepper;
>> 66 } // fIntgrDriver->pIntStepper;}
>> 67 delete fIntgrDriver;
>> 68 }
71 69
72 G4bool G4ChordFinder::gVerboseCtor = false; << 70 // ......................................................................
73 // ........................................... <<
74 71
75 G4ChordFinder::G4ChordFinder(G4VIntegrationDri << 72 G4double
76 : fDefaultDeltaChord(0.25 * mm), fIntgrDrive << 73 G4ChordFinder::AdvanceChordLimited( G4FieldTrack& yCurrent,
>> 74 G4double stepMax,
>> 75 G4double epsStep )
77 { 76 {
78 // Simple constructor -- it does not create << 77 G4double stepPossible;
79 if( gVerboseCtor ) << 78 G4double dyErr;
>> 79 G4FieldTrack yEnd( yCurrent);
>> 80 G4double startCurveLen= yCurrent.GetCurveLength();
>> 81 G4bool dbg= false;
>> 82
>> 83 #ifdef G4VERBOSE
>> 84 if( dbg )
>> 85 G4cerr << "Entered AdvanceChordLimited with:\n yCurrent: " << yCurrent
>> 86 << " and initial Step=stepMax=" << stepMax << " mm. " << G4endl;
>> 87 #endif
>> 88
>> 89 stepPossible= FindNextChord(yCurrent, stepMax, yEnd, dyErr, epsStep);
>> 90 G4bool good_advance;
>> 91 if ( dyErr < epsStep * stepPossible )
>> 92 {
>> 93 // Accept this accuracy.
>> 94 yCurrent = yEnd;
>> 95 good_advance = true;
>> 96 }
>> 97 else
80 { 98 {
81 G4cout << "G4ChordFinder: Simple construct << 99 // Advance more accurately to "end of chord"
>> 100 good_advance = fIntgrDriver->AccurateAdvance(yCurrent, stepPossible, epsStep);
>> 101 #ifdef G4VERBOSE
>> 102 if (dbg) G4cerr << "Accurate advance to end of chord attemped"
>> 103 << "with result " << good_advance << G4endl ;
>> 104 #endif
>> 105 if ( ! good_advance ){
>> 106 // In this case the driver could not do the full distance
>> 107 stepPossible= yCurrent.GetCurveLength()-startCurveLen;
>> 108 }
82 } 109 }
83 << 110
84 fDeltaChord = fDefaultDeltaChord; // P << 111 #ifdef G4VERBOSE
>> 112 if( dbg ) G4cerr << "Exiting FindNextChord Limited with:\n yCurrent: "
>> 113 << yCurrent<< G4endl;
>> 114 #endif
>> 115
>> 116 return stepPossible;
85 } 117 }
86 118
87 // ........................................... << 119 // #define TEST_CHORD_PRINT 1
88 120
89 G4ChordFinder::G4ChordFinder( G4MagneticField* << 121 // ..............................................................................
90 G4double << 122
91 G4MagIntegratorS << 123 G4double
92 G4int << 124 G4ChordFinder::FindNextChord( const G4FieldTrack yStart,
93 : fDefaultDeltaChord(0.25 * mm) << 125 G4double stepMax,
>> 126 G4FieldTrack& yEnd, // Endpoint
>> 127 G4double& dyErr, // Error of endpoint
>> 128 G4double epsStep )
>> 129
>> 130 // Returns Length of Step taken
94 { 131 {
95 // Construct the Chord Finder << 132 // G4int stepRKnumber=0;
96 // by creating in inverse order the Driver, << 133 G4FieldTrack yCurrent= yStart;
97 constexpr G4int nVar6 = 6; // Components i << 134 G4double stepTrial;
>> 135 G4double dydx[G4FieldTrack::ncompSVEC];
>> 136
>> 137 // 1.) Try to "leap" to end of interval
>> 138 // 2.) Evaluate if resulting chord gives d_chord that is good enough.
>> 139 // 2a.) If d_chord is not good enough, find one that is.
98 140
99 fDeltaChord = fDefaultDeltaChord; // P << 141 G4bool validEndPoint= false;
>> 142 G4double dChordStep, oldStepTrial, stepOfLastGoodChord;
100 143
101 G4cout << " G4ChordFinder: stepperDriverId: << 144 fIntgrDriver-> GetDerivatives( yCurrent, dydx ) ;
102 145
103 G4bool useFSALstepper = (stepperDriverId << 146 G4int noTrials=0;
104 G4bool useTemplatedStepper= (stepperDriverId <<
105 G4bool useRegularStepper = (stepperDriverId <<
106 G4bool useBfieldDriver = (stepperDriverId <<
107 G4bool useG4QSSDriver = (stepperDriverId <<
108 <<
109 if( stepperDriverId == kQss3DriverType) <<
110 { <<
111 stepperDriverId = kQss2DriverType; <<
112 G4cout << " G4ChordFinder: QSS 3 is curren <<
113 } <<
114 147
115 using EquationType = G4Mag_UsualEqRhs; << 148 stepTrial = G4std::min( stepMax,
116 << 149 (1-perThousand)*fLastStepEstimate_Unconstrained );
117 using TemplatedStepperType = <<
118 G4TDormandPrince45<EquationType,nVar6 <<
119 const char* TemplatedStepperName = <<
120 "G4TDormandPrince745 (templated Dormand- <<
121 <<
122 using RegularStepperType = <<
123 G4DormandPrince745; // 5th order embe <<
124 // G4ClassicalRK4; // The old <<
125 // G4CashKarpRKF45; // First em <<
126 // G4BogackiShampine45; // High eff <<
127 // G4NystromRK4; // Nystrom <<
128 // G4RK547FEq1; // or 2 or 3 <<
129 const char* RegularStepperName = <<
130 "G4DormandPrince745 (aka DOPRI5): 5th/4t <<
131 // "BogackiShampine 45 (Embedded 5th/4th <<
132 // "Nystrom stepper 4th order"; <<
133 <<
134 using NewFsalStepperType = G4DormandPrince74 <<
135 <<
136 const char* NewFSALStepperName = <<
137 "G4RK574FEq1> FSAL 4th/5th order 7-stage <<
138 <<
139 #ifdef G4DEBUG_FIELD <<
140 static G4bool verboseDebug = true; <<
141 if( verboseDebug ) <<
142 { <<
143 G4cout << "G4ChordFinder 2nd Constructor <<
144 G4cout << " Arguments: " << G4endl <<
145 << " - min step = " << stepMinimum <<
146 << " - stepper ptr provided : " <<
147 << ( pItsStepper==nullptr ? " no <<
148 if( pItsStepper==nullptr ) <<
149 G4cout << " - stepper/driver Id = " << <<
150 << " useFSAL = " << useFSALste <<
151 << " , useTemplated = " << use <<
152 << " , useRegular = " << useRe <<
153 << " , useFSAL = " << useFSALs <<
154 << G4endl; <<
155 } <<
156 #endif <<
157 150
158 // useHigherStepper = forceHigherEffiencySte << 151 do
>> 152 {
>> 153 G4double stepForChord; // , stepForAccuracy;
>> 154
>> 155 yCurrent = yStart; // Always start from initial point
159 156
160 auto pEquation = new G4Mag_UsualEqRhs(theMa << 157 fIntgrDriver->QuickAdvance( yCurrent, dydx, stepTrial, dChordStep, dyErr);
161 fEquation = pEquation; <<
162 158
163 // G4MagIntegratorStepper* regularStepper = << 159 // First debug print
164 // G4VFSALIntegrationStepper* fsalStepper = <<
165 // G4MagIntegratorStepper* oldFSALStepper = <<
166 160
167 G4bool errorInStepperCreation = false; << 161 // We check whether the criterion is met here.
>> 162 validEndPoint = AcceptableMissDist(dChordStep);
>> 163 // && (dyErr < eps) ;
168 164
169 std::ostringstream message; // In case of f << 165 oldStepTrial = stepTrial;
170 166
171 if( pItsStepper != nullptr ) << 167 // This method estimates to step size for a good chord.
172 { << 168 stepForChord = NewStep(stepTrial, dChordStep, fLastStepEstimate_Unconstrained );
173 if( gVerboseCtor ) <<
174 { <<
175 G4cout << " G4ChordFinder: Creating G4I <<
176 << " stepMinimum = " << stepMini <<
177 << " numVar= " << pItsStepper->G <<
178 } <<
179 169
180 // Stepper type is not known - so must us << 170 if( ! validEndPoint ) {
181 if(pItsStepper->isQSS()) << 171 stepTrial = stepForChord;
182 { << 172 #if 0
183 // fIntgrDriver = pItsStepper->build_ << 173 // Possible complementary approach:
184 G4Exception("G4ChordFinder::G4ChordFi << 174 // Get the driver to calculate the new step size, if it is needed
185 "GeomField1001", FatalEx << 175 stepForAccuracy = fIntgrDriver->ComputeNewStepSize( dyErr/(epsStep*oldStepTrial),
186 "Cannot provide QSS ste << 176 stepTrial);
187 } << 177 stepTrial = G4std::min(stepForChord, stepForAccuracy);
188 else << 178 #endif
189 { << 179
190 fIntgrDriver = new G4IntegrationDrive << 180 // if(dbg) G4cerr<<"Dchord too big. Try new hstep="<<stepTrial<<G4endl;
191 pItsStepper, <<
192 // Stepper type is not known - so mus <<
193 // Non-interpolating driver used by d <<
194 // WAS: fIntgrDriver = pItsStepper-> <<
195 } <<
196 // -- Older: <<
197 // G4cout << " G4ChordFinder: Creating G4 <<
198 // Type is not known - so must use old cl <<
199 // fIntgrDriver = new G4MagInt_Driver( st <<
200 // pItsSt <<
201 } <<
202 else if ( useTemplatedStepper ) <<
203 { <<
204 if( gVerboseCtor ) <<
205 { <<
206 G4cout << " G4ChordFinder: Creating Te <<
207 << TemplatedStepperName << G4en <<
208 } <<
209 // RegularStepperType* regularStepper = n <<
210 auto templatedStepper = new TemplatedStep <<
211 // *** *************** <<
212 // <<
213 // Alternative - for G4NystromRK4: <<
214 // = new G4NystromRK4(pEquation, 0.1*mm ) <<
215 fRegularStepperOwned = templatedStepper; <<
216 if( templatedStepper == nullptr ) <<
217 { <<
218 message << "Templated Stepper instanti <<
219 message << "G4ChordFinder: Attempted t <<
220 << TemplatedStepperName << " t <<
221 errorInStepperCreation = true; <<
222 } <<
223 else <<
224 { <<
225 fIntgrDriver = new G4IntegrationDriver <<
226 stepMinimum, templatedStepper, nVar <<
227 if( gVerboseCtor ) <<
228 { <<
229 G4cout << " G4ChordFinder: Using G4 <<
230 } <<
231 } <<
232 <<
233 } <<
234 else if ( useRegularStepper ) // Plain st <<
235 { <<
236 auto regularStepper = new RegularStepperT <<
237 // *** *************** <<
238 fRegularStepperOwned = regularStepper; <<
239 <<
240 if( gVerboseCtor ) <<
241 { <<
242 G4cout << " G4ChordFinder: Creating Dr <<
243 } <<
244 <<
245 if( regularStepper == nullptr ) <<
246 { <<
247 message << "Regular Stepper instantiat <<
248 message << "G4ChordFinder: Attempted t <<
249 << RegularStepperName << " typ <<
250 errorInStepperCreation = true; <<
251 } <<
252 else <<
253 { <<
254 auto dp5= dynamic_cast<G4DormandPrince <<
255 if( dp5 != nullptr ) <<
256 { <<
257 fIntgrDriver = new G4InterpolationD <<
258 stepMinimum, <<
259 if( gVerboseCtor ) <<
260 { <<
261 G4cout << " Using InterpolationD <<
262 } <<
263 } <<
264 else <<
265 { <<
266 fIntgrDriver = new G4IntegrationDri <<
267 stepMinimum, <<
268 if( gVerboseCtor ) <<
269 { <<
270 G4cout << " Using IntegrationDri <<
271 } <<
272 } <<
273 } <<
274 } <<
275 else if ( useBfieldDriver ) <<
276 { <<
277 auto regularStepper = new G4DormandPrince <<
278 // *** *************** <<
279 // <<
280 fRegularStepperOwned = regularStepper; <<
281 <<
282 { <<
283 using SmallStepDriver = G4Interpolatio <<
284 using LargeStepDriver = G4IntegrationD <<
285 <<
286 fLongStepper = std::make_unique<G4Heli <<
287 <<
288 fIntgrDriver = new G4BFieldIntegration <<
289 std::make_unique<SmallStepDriver>(st <<
290 regularStepper, regularStepper-> <<
291 std::make_unique<LargeStepDriver>(st <<
292 fLongStepper.get(), regularStepp <<
293 <<
294 if( fIntgrDriver == nullptr) <<
295 { <<
296 message << "Using G4BFieldIntegrati <<
297 << RegularStepperName << " <<
298 message << "Driver instantiation FA <<
299 G4Exception("G4ChordFinder::G4Chord <<
300 "GeomField1001", JustWa <<
301 } <<
302 } <<
303 } <<
304 else if( useG4QSSDriver ) <<
305 { <<
306 if( stepperDriverId == kQss2DriverType ) <<
307 { <<
308 auto qssStepper2 = G4QSSDriverCreator:: <<
309 if( gVerboseCtor ) <<
310 { <<
311 G4cout << "-- Created QSS-2 stepper" <<
312 } <<
313 fIntgrDriver = G4QSSDriverCreator::Crea <<
314 } <<
315 else <<
316 { <<
317 auto qssStepper3 = G4QSSDriverCreator:: <<
318 if( gVerboseCtor ) <<
319 { <<
320 G4cout << "-- Created QSS-3 stepper" <<
321 } <<
322 fIntgrDriver = G4QSSDriverCreator::Crea <<
323 } <<
324 if( gVerboseCtor ) <<
325 { <<
326 G4cout << "-- G4ChordFinder: Using QSS <<
327 } 181 }
328 } << 182 #ifdef TEST_CHORD_PRINT
329 else << 183 G4cout.precision(5);
330 { << 184 G4cout << " ChF/fnc: notrial " << G4std::setw( 3) << noTrials
331 auto fsalStepper= new NewFsalStepperType << 185 << " this_step= " << G4std::setw(10) << oldStepTrial;
332 // *** ****************** << 186 if( fabs( (dChordStep / fDeltaChord) - 1.0 ) < 0.001 ){
333 fNewFSALStepperOwned = fsalStepper; << 187 G4cout.precision(8);
334 << 188 G4cout << " dChordStep= " << G4std::setw(12) << dChordStep;
335 if( fsalStepper == nullptr ) << 189 }else{
336 { << 190 G4cout.precision(6);
337 message << "Stepper instantiation FAIL << 191 G4cout << " dChordStep= " << G4std::setw(12) << dChordStep;
338 message << "Attempted to instantiate " <<
339 << NewFSALStepperName << " typ <<
340 G4Exception("G4ChordFinder::G4ChordFin <<
341 "GeomField1001", JustWarni <<
342 errorInStepperCreation = true; <<
343 } 192 }
>> 193 if( dChordStep > fDeltaChord )
>> 194 G4cout << " d+";
344 else 195 else
345 { << 196 G4cout << " d-";
346 fIntgrDriver = new << 197 G4cout.precision(5);
347 G4FSALIntegrationDriver<NewFsalStep << 198 G4cout << " new_step= " << G4std::setw(10) << fLastStepEstimate_Unconstrained
348 fsal << 199 << " new_step_constr= " << G4std::setw(10) << stepTrial << endl;
349 // ==== Create the driver which k << 200 #endif
350 << 201 noTrials++;
351 if( fIntgrDriver == nullptr ) <<
352 { <<
353 message << "Using G4FSALIntegration <<
354 << NewFSALStepperName << G4 <<
355 message << "Integration Driver inst <<
356 G4Exception("G4ChordFinder::G4Chord <<
357 "GeomField1001", JustWa <<
358 } <<
359 } <<
360 } 202 }
>> 203 while( ! validEndPoint ); // End of do-while RKD
361 204
362 // -- Main work is now done << 205 stepOfLastGoodChord = stepTrial;
363 << 206 #ifdef TEST_CHORD_PRINT
364 // Now check that no error occured, and r << 207 G4cout << "ChordF/FindNextChord: NoTrials= " << noTrials
365 << 208 << " StepForGoodChord=" << G4std::setw(10) << stepTrial << endl;
366 // To test failure to create driver << 209 #endif
367 // delete fIntgrDriver; <<
368 // fIntgrDriver = nullptr; <<
369 <<
370 // Detect and report Error conditions <<
371 // <<
372 if( errorInStepperCreation || (fIntgrDriver <<
373 { <<
374 std::ostringstream errmsg; <<
375 <<
376 if( errorInStepperCreation ) <<
377 { <<
378 errmsg << "ERROR> Failure to create S <<
379 << " ------------------- <<
380 } <<
381 if (fIntgrDriver == nullptr ) <<
382 { <<
383 errmsg << "ERROR> Failure to create I <<
384 << G4endl <<
385 << " ------------------- <<
386 << G4endl; <<
387 } <<
388 const std::string BoolName[2]= { "False", <<
389 errmsg << " Configuration: (constructor <<
390 << " provided Stepper = " << pI <<
391 << " stepper/driver Id = " << step <<
392 << " useTemplated = " << BoolNam <<
393 << " useRegular = " << BoolName[ <<
394 << " useFSAL = " << BoolName[use <<
395 << " using combo BField Driver = <<
396 BoolName[ ! (useFSALstepper <<
397 || useRegularSt <<
398 << G4endl; <<
399 errmsg << message.str(); <<
400 errmsg << "Aborting."; <<
401 G4Exception("G4ChordFinder::G4ChordFinder <<
402 "GeomField0003", FatalExcepti <<
403 } <<
404 <<
405 assert( ( pItsStepper != nullptr ) <<
406 || ( fRegularStepperOwned != nullptr <<
407 || ( fNewFSALStepperOwned != nullptr <<
408 || useG4QSSDriver <<
409 ); <<
410 assert( fIntgrDriver != nullptr ); <<
411 } <<
412 <<
413 // ........................................... <<
414 210
415 G4ChordFinder::~G4ChordFinder() << 211 yEnd= yCurrent;
416 { << 212 return stepTrial;
417 delete fEquation; <<
418 delete fRegularStepperOwned; <<
419 delete fNewFSALStepperOwned; <<
420 delete fCachedField; <<
421 delete fIntgrDriver; <<
422 } 213 }
423 214
>> 215 // ----------------------------------------------------------------------------
>> 216 #if 0
>> 217 // First debug print // older OPTIONAL code
>> 218 // #ifdef G4VERBOSE
>> 219 if( dbg ) {
>> 220 G4cerr << "Returned from QuickAdvance with: yCur=" << yCurrent <<G4endl;
>> 221 G4cerr << " dChordStep= "<< dChordStep <<" dyErr=" << dyErr << G4endl;
>> 222 }
>> 223 #endif
>> 224 // ----------------------------------------------------------------------------
>> 225
424 // ........................................... 226 // ...........................................................................
425 227
426 G4FieldTrack << 228 G4double G4ChordFinder::NewStep(G4double stepTrialOld,
427 G4ChordFinder::ApproxCurvePointS( const G4Fiel << 229 G4double dChordStep, // Current dchord achieved.
428 const G4Fiel << 230 G4double& stepEstimate_Unconstrained )
429 const G4Fiel << 231
430 const G4Thre <<
431 const G4Thre <<
432 const G4Thre <<
433 G4bool <<
434 { 232 {
435 // ApproxCurvePointS is 2nd implementation o << 233 G4double stepTrial;
436 // Use Brent Algorithm (or InvParabolic) whe << 234 static G4double lastStepTrial = 1., lastDchordStep= 1.;
437 // Given a starting curve point A (CurveA_Po <<
438 // (CurveB_PointVelocity), a point E which i <<
439 // and a point F which is on the curve (fir <<
440 // point S on the curve closer to point E. <<
441 // While advancing towards S utilise 'eps_st <<
442 // relative accuracy of each Step. <<
443 <<
444 G4FieldTrack EndPoint(CurveA_PointVelocity); <<
445 if(!first) { EndPoint = ApproxCurveV; } <<
446 <<
447 G4ThreeVector Point_A,Point_B; <<
448 Point_A=CurveA_PointVelocity.GetPosition(); <<
449 Point_B=CurveB_PointVelocity.GetPosition(); <<
450 235
451 G4double xa,xb,xc,ya,yb,yc; << 236 #if 1
452 << 237 const G4double threshold = 1.21, multiplier = 0.9; // 0.9 < 1 / sqrt(1.21)
453 // InverseParabolic. AF Intersects (First Pa <<
454 238
455 if(first) <<
456 { <<
457 xa=0.; <<
458 ya=(PointG-Point_A).mag(); <<
459 xb=(Point_A-CurrentF_Point).mag(); <<
460 yb=-(PointG-CurrentF_Point).mag(); <<
461 xc=(Point_A-Point_B).mag(); <<
462 yc=-(CurrentE_Point-Point_B).mag(); <<
463 } <<
464 else <<
465 { <<
466 xa=0.; <<
467 ya=(Point_A-CurrentE_Point).mag(); <<
468 xb=(Point_A-CurrentF_Point).mag(); <<
469 yb=(PointG-CurrentF_Point).mag(); <<
470 xc=(Point_A-Point_B).mag(); <<
471 yc=-(Point_B-PointG).mag(); <<
472 if(xb==0.) <<
473 { <<
474 EndPoint = ApproxCurvePointV(CurveA_Poin <<
475 CurrentE_Po <<
476 return EndPoint; <<
477 } <<
478 } <<
479 239
480 const G4double tolerance = 1.e-12; << 240 stepEstimate_Unconstrained = stepTrialOld * sqrt( fDeltaChord / dChordStep );
481 if(std::abs(ya)<=tolerance||std::abs(yc)<=to << 241 stepTrial = 0.98 * stepEstimate_Unconstrained;
482 { << 242
483 ; // What to do for the moment: return the << 243 if ( dChordStep < threshold * fDeltaChord ){
484 // then PropagatorInField will take care << 244 stepTrial= stepTrialOld * multiplier;
485 } 245 }
486 else <<
487 { <<
488 G4double test_step = InvParabolic(xa,ya,xb <<
489 G4double curve; <<
490 if(first) <<
491 { <<
492 curve=std::abs(EndPoint.GetCurveLength() <<
493 -ApproxCurveV.GetCurveLeng <<
494 } <<
495 else <<
496 { <<
497 test_step = test_step - xb; <<
498 curve=std::abs(EndPoint.GetCurveLength() <<
499 -CurveB_PointVelocity.GetC <<
500 xb = (CurrentF_Point-Point_B).mag(); <<
501 } <<
502 <<
503 if(test_step<=0) { test_step=0.1*xb; } <<
504 if(test_step>=xb) { test_step=0.5*xb; } <<
505 if(test_step>=curve){ test_step=0.5*curve; <<
506 <<
507 if(curve*(1.+eps_step)<xb) // Similar to R <<
508 { // G4VIntersect <<
509 test_step=0.5*curve; <<
510 } <<
511 246
512 fIntgrDriver->AccurateAdvance(EndPoint,tes << 247 lastStepTrial = stepTrialOld;
513 << 248 lastDchordStep= dChordStep;
514 #ifdef G4DEBUG_FIELD << 249 #else
515 // Printing Brent and Linear Approximation << 250 if ( dChordStep > 1000. * fDeltaChord ){
516 // << 251 stepTrial= stepTrialOld * 0.03;
517 G4cout << "G4ChordFinder::ApproxCurvePoint << 252 }else{
518 << test_step << " EndPoint = " << << 253 if ( dChordStep > 100. * fDeltaChord ){
519 << 254 stepTrial= stepTrialOld * 0.1;
520 // Test Track << 255 }else{
521 // << 256 // Keep halving the length until dChordStep OK
522 G4FieldTrack TestTrack( CurveA_PointVeloci << 257 stepTrial= stepTrialOld * 0.5;
523 TestTrack = ApproxCurvePointV( CurveA_Poin << 258 }
524 CurveB_Poin <<
525 CurrentE_Po <<
526 G4cout.precision(14); <<
527 G4cout << "G4ChordFinder::BrentApprox = " <<
528 G4cout << "G4ChordFinder::LinearApprox= " <<
529 #endif <<
530 } 259 }
531 return EndPoint; << 260 #endif
532 } <<
533 261
>> 262 // A more sophisticated chord-finder could figure out a better
>> 263 // stepTrial, from dChordStep and the required d_geometry
>> 264 // eg
>> 265 // Calculate R, r_helix (eg at orig point)
>> 266 // if( stepTrial < 2 pi R )
>> 267 // stepTrial = R arc_cos( 1 - fDeltaChord / r_helix )
>> 268 // else
>> 269 // ??
534 270
535 // ........................................... << 271 return stepTrial;
>> 272 }
536 273
537 G4FieldTrack G4ChordFinder:: << 274 //
538 ApproxCurvePointV( const G4FieldTrack& CurveA_ << 275 // Given a starting curve point A (CurveA_PointVelocity), a later
539 const G4FieldTrack& CurveB_ << 276 // curve point B (CurveB_PointVelocity) and a point E which is (generally)
540 const G4ThreeVector& Curren << 277 // not on the curve, find and return a point F which is on the curve and
541 G4double eps_step) << 278 // which is close to E. While advancing towards F utilise eps_step
>> 279 // as a measure of the relative accuracy of each Step.
>> 280
>> 281 G4FieldTrack G4ChordFinder::ApproxCurvePointV(
>> 282 const G4FieldTrack& CurveA_PointVelocity,
>> 283 const G4FieldTrack& CurveB_PointVelocity,
>> 284 const G4ThreeVector& CurrentE_Point,
>> 285 G4double eps_step)
542 { 286 {
543 // If r=|AE|/|AB|, and s=true path lenght (A << 287 // 1st implementation:
544 // return the point that is r*s along the cu << 288 // if r=|AE|/|AB|, and s=true path lenght (AB)
545 << 289 // return the point that is r*s along the curve!
546 G4FieldTrack Current_PointVelocity = Curve << 290
>> 291 G4FieldTrack Current_PointVelocity= CurveA_PointVelocity;
547 292
548 G4ThreeVector CurveA_Point= CurveA_PointVel 293 G4ThreeVector CurveA_Point= CurveA_PointVelocity.GetPosition();
549 G4ThreeVector CurveB_Point= CurveB_PointVel 294 G4ThreeVector CurveB_Point= CurveB_PointVelocity.GetPosition();
550 295
551 G4ThreeVector ChordAB_Vector= CurveB_Point 296 G4ThreeVector ChordAB_Vector= CurveB_Point - CurveA_Point;
552 G4ThreeVector ChordAE_Vector= CurrentE_Poin 297 G4ThreeVector ChordAE_Vector= CurrentE_Point - CurveA_Point;
553 298
554 G4double ABdist= ChordAB_Vector.mag(); 299 G4double ABdist= ChordAB_Vector.mag();
555 G4double curve_length; // A curve length 300 G4double curve_length; // A curve length of AB
556 G4double AE_fraction; 301 G4double AE_fraction;
557 302
558 curve_length= CurveB_PointVelocity.GetCurveL << 303 curve_length=
559 - CurveA_PointVelocity.GetCurveL << 304 CurveB_PointVelocity.GetCurveLength() - CurveA_PointVelocity.GetCurveLength();
560 << 305
561 G4double integrationInaccuracyLimit= std::ma << 306 // const
562 if( curve_length < ABdist * (1. - integratio << 307 G4double integrationInaccuracyLimit= G4std::max( perMillion, 0.5*eps_step );
563 { << 308 if( curve_length < ABdist * (1. - integrationInaccuracyLimit) ){
564 #ifdef G4DEBUG_FIELD << 309 // #ifdef G4DEBUG
565 G4cerr << " Warning in G4ChordFinder::Appr << 310 G4cerr << " Warning in G4ChordFinder::ApproxCurvePoint: " << G4endl <<
566 << G4endl << 311 " The two points are further apart than the curve length " << G4endl <<
567 << " The two points are further apa << 312 " Dist = " << ABdist <<
568 << G4endl << 313 " curve length = " << curve_length
569 << " Dist = " << ABdist << 314 << " relativeDiff = " << (curve_length-ABdist)/ABdist
570 << " curve length = " << curve_leng << 315 << G4endl;
571 << " relativeDiff = " << (curve_len << 316 // #endif
572 << G4endl; << 317 if( curve_length < ABdist * (1. - 10*eps_step) ) {
573 if( curve_length < ABdist * (1. - 10*eps_s << 318 // #ifdef G4DEBUG
574 { << 319 G4cerr << " ERROR: the size of the above difference exceeds allowed limits. Aborting."
575 std::ostringstream message; << 320 << G4endl;
576 message << "Unphysical curve length." << << 321 // #endif
577 << "The size of the above differ << 322 G4Exception("G4ChordFinder::ApproxCurvePoint> Unphysical curve length.");
578 << G4endl <<
579 << "Aborting."; <<
580 G4Exception("G4ChordFinder::ApproxCurveP <<
581 FatalException, message); <<
582 } 323 }
583 #endif << 324 // Take default corrective action:
584 // Take default corrective action: adjust << 325 // --> adjust the maximum curve length.
585 // NOTE: this case only happens for relati << 326 // NOTE: this case only happens for relatively straight paths.
586 // curve_length = ABdist; << 327 curve_length = ABdist;
587 } 328 }
588 329
589 G4double new_st_length; << 330 G4double new_st_length;
590 331
591 if ( ABdist > 0.0 ) << 332 if ( ABdist > 0.0 ){
592 { <<
593 AE_fraction = ChordAE_Vector.mag() / ABdi 333 AE_fraction = ChordAE_Vector.mag() / ABdist;
594 } << 334 }else{
595 else << 335 G4cerr << " Error in G4ChordFinder::ApproxCurvePoint: A and B are the same point\n" <<
596 { << 336 " Chord AB length = " << ChordAE_Vector.mag() << G4endl << G4endl;
597 AE_fraction = 0.5; 337 AE_fraction = 0.5; // Guess .. ?;
598 #ifdef G4DEBUG_FIELD <<
599 G4cout << "Warning in G4ChordFinder::Appr <<
600 << " A and B are the same point!" <<
601 << " Chord AB length = " << ChordA <<
602 << G4endl; <<
603 #endif <<
604 } 338 }
605 339
606 if( (AE_fraction> 1.0 + perMillion) || (AE_f << 340 if( (AE_fraction> 1.0 + perMillion) || (AE_fraction< 0.) ){
607 { << 341 G4cerr << " G4ChordFinder::ApproxCurvePointV: Warning: Anomalous condition:AE > AB or AE/AB <= 0 " << G4endl <<
608 #ifdef G4DEBUG_FIELD << 342 " AE_fraction = " << AE_fraction << G4endl <<
609 G4cerr << " G4ChordFinder::ApproxCurvePoin << 343 " Chord AE length = " << ChordAE_Vector.mag() << G4endl <<
610 << " Anomalous condition:AE > AB or << 344 " Chord AB length = " << ABdist << G4endl << G4endl;
611 << " AE_fraction = " << AE_fract << 345 G4cerr << " OK if this condition occurs after a recalculation of 'B'" << G4endl
612 << " Chord AE length = " << Chord << 346 << " Otherwise it is an error. " << G4endl ;
613 << " Chord AB length = " << ABdis <<
614 G4cerr << " OK if this condition occurs af <<
615 << G4endl << " Otherwise it is an e <<
616 #endif <<
617 // This course can now result if B has be 347 // This course can now result if B has been re-evaluated,
618 // without E being recomputed (1 July 99) << 348 // without E being recomputed (1 July 99)
619 // In this case this is not a "real error << 349 // In this case this is not a "real error" - but it undesired
620 // and we cope with it by a default corre << 350 // and we cope with it by a default corrective action ...
621 // <<
622 AE_fraction = 0.5; 351 AE_fraction = 0.5; // Default value
623 } 352 }
624 353
625 new_st_length = AE_fraction * curve_length; << 354 new_st_length= AE_fraction * curve_length;
626 355
627 if ( AE_fraction > 0.0 ) << 356 G4bool good_advance;
628 { << 357 if ( AE_fraction > 0.0 ) {
629 fIntgrDriver->AccurateAdvance(Current_Poi << 358 good_advance =
630 new_st_leng << 359 fIntgrDriver->AccurateAdvance(Current_PointVelocity,
631 // << 360 new_st_length,
>> 361 eps_step ); // Relative accuracy
632 // In this case it does not matter if it 362 // In this case it does not matter if it cannot advance the full distance
633 } 363 }
634 364
635 // If there was a memory of the step_length << 365 // If there was a memory of the step_length actually require at the start
636 // of the integration Step, this could be re 366 // of the integration Step, this could be re-used ...
637 367
638 G4cout.precision(14); <<
639 <<
640 return Current_PointVelocity; 368 return Current_PointVelocity;
641 } 369 }
642 370
643 // ........................................... <<
644 371
645 std::ostream& operator<<( std::ostream& os, co <<
646 { <<
647 // Dumping the state of G4ChordFinder <<
648 os << "State of G4ChordFinder : " << std::e <<
649 os << " delta_chord = " << cf.fDeltaCh <<
650 os << " Default d_c = " << cf.fDefault <<
651 <<
652 os << " stats-verbose = " << cf.fStatsVe <<
653 <<
654 return os; <<
655 } <<
656 372