Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-01-18 09:58:02 

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 // G4ChordFinderDelegate inline methods implementation
0027 //
0028 // Created: D.Sorokin
0029 // --------------------------------------------------------------------
0030 
0031 #include <iomanip>
0032 
0033 template <class Driver>
0034 G4ChordFinderDelegate<Driver>::~G4ChordFinderDelegate()
0035 {
0036 #ifdef G4VERBOSE
0037     if (GetDriver().GetVerboseLevel() > 0)
0038     {
0039         PrintStatistics();
0040     }
0041 #endif
0042 }
0043 
0044 template <class Driver>
0045 void G4ChordFinderDelegate<Driver>::ResetStepEstimate()
0046 {
0047     fLastStepEstimate_Unconstrained = DBL_MAX;  
0048 }
0049 
0050 template <class Driver>
0051 Driver& G4ChordFinderDelegate<Driver>::GetDriver()
0052 {
0053     return static_cast<Driver&>(*this);
0054 }
0055 
0056 template <class Driver>
0057 G4double G4ChordFinderDelegate<Driver>::
0058 AdvanceChordLimitedImpl(G4FieldTrack& yCurrent, G4double stepMax, 
0059                         G4double epsStep, G4double chordDistance)
0060 {
0061     G4double dyErr;
0062     G4FieldTrack yEnd = yCurrent;
0063     G4double nextStep;
0064 
0065     const G4double stepPossible = FindNextChord(yCurrent, stepMax,
0066                                                 epsStep, chordDistance,
0067                                                 yEnd, dyErr, nextStep);
0068     if (dyErr < epsStep * stepPossible)
0069     {
0070          // Accept this accuracy.
0071          //
0072          yCurrent = yEnd;
0073          return stepPossible;
0074     }
0075 
0076     // Advance more accurately to "end of chord"
0077     //
0078     const G4double startCurveLen = yCurrent.GetCurveLength();
0079     const G4bool goodAdvance =
0080           GetDriver().AccurateAdvance(yCurrent,stepPossible,epsStep,nextStep);
0081 
0082     return goodAdvance ? stepPossible
0083                        : yCurrent.GetCurveLength() - startCurveLen;
0084 }
0085 
0086 // Returns Length of Step taken
0087 //
0088 template <class T>
0089 G4double G4ChordFinderDelegate<T>::
0090 FindNextChord(const G4FieldTrack& yStart,
0091               G4double stepMax,
0092               G4double epsStep,
0093               G4double chordDistance,
0094               G4FieldTrack& yEnd, // Endpoint
0095               G4double& dyErrPos, // Error of endpoint
0096               G4double& stepForAccuracy)
0097 {
0098     //  1.)  Try to "leap" to end of interval
0099     //  2.)  Evaluate if resulting chord gives d_chord that is good enough.
0100     // 2a.)  If d_chord is not good enough, find one that is.
0101 
0102     G4double dydx[G4FieldTrack::ncompSVEC];
0103   
0104     G4bool validEndPoint = false;
0105     G4double dChordStep, lastStepLength;
0106 
0107     GetDriver().GetDerivatives(yStart, dydx);
0108 
0109     const G4double safetyFactor = fFirstFraction; //  0.975 or 0.99 ? was 0.999
0110 
0111     G4double stepTrial = std::min(stepMax,
0112                                   safetyFactor*fLastStepEstimate_Unconstrained);
0113 
0114     G4double newStepEst_Uncons = 0.0; 
0115     G4double stepForChord;
0116 
0117     G4int noTrials = 1;
0118     constexpr G4int maxTrials = 75; // Avoid endless loop for bad convergence 
0119     for (; noTrials < maxTrials; ++noTrials)
0120     {
0121         yEnd = yStart; // Always start from initial point  
0122         GetDriver().QuickAdvance(yEnd, dydx, stepTrial, dChordStep, dyErrPos);
0123         lastStepLength = stepTrial; 
0124 
0125         validEndPoint = dChordStep < chordDistance;
0126         stepForChord = NewStep(stepTrial, dChordStep,
0127                                chordDistance, newStepEst_Uncons);
0128         if (validEndPoint)
0129         {
0130             break;
0131         }
0132 
0133         if (stepTrial <= 0.0)
0134         {
0135             stepTrial = stepForChord;
0136         }
0137         else if (stepForChord <= stepTrial)
0138         {
0139             // Reduce by a fraction, possibly up to 20% 
0140             stepTrial = std::min( stepForChord, fFractionLast * stepTrial);
0141         }
0142         else
0143         {
0144             stepTrial *= 0.1;
0145         }
0146     }
0147 
0148     if (noTrials >= maxTrials)
0149     {
0150         std::ostringstream message;
0151         message << "Exceeded maximum number of trials= " << maxTrials << G4endl
0152                 << "Current sagita dist= " << dChordStep << G4endl
0153                 << "Max sagita dist= " << chordDistance << G4endl
0154                 << "Step sizes (actual and proposed): " << G4endl
0155                 << "Last trial =         " << lastStepLength  << G4endl
0156                 << "Next trial =         " << stepTrial  << G4endl
0157                 << "Proposed for chord = " << stepForChord  << G4endl;
0158         G4Exception("G4ChordFinder::FindNextChord()", "GeomField0003",
0159                     JustWarning, message);
0160     }
0161 
0162     if (newStepEst_Uncons > 0.0)
0163     {
0164         fLastStepEstimate_Unconstrained = newStepEst_Uncons;
0165     }
0166 
0167     AccumulateStatistics(noTrials);
0168 
0169 
0170     // Calculate the step size required for accuracy, if it is needed
0171     G4double dyErr_relative = dyErrPos / (epsStep * lastStepLength);
0172     stepForAccuracy = dyErr_relative > 1 ? 
0173         GetDriver().ComputeNewStepSize(dyErr_relative, lastStepLength) : 0;
0174 
0175   return stepTrial; 
0176 }
0177 
0178 // Is called to estimate the next step size, even for successful steps,
0179 // in order to predict an accurate 'chord-sensitive' first step
0180 // which is likely to assist in more performant 'stepping'.
0181 //
0182 template <class T>
0183 G4double G4ChordFinderDelegate<T>::
0184 NewStep(G4double stepTrialOld, 
0185         G4double dChordStep, // Curr. dchord achieved
0186         G4double fDeltaChord,
0187         G4double& stepEstimate_Unconstrained)  
0188 {
0189     G4double stepTrial;
0190 
0191     if (dChordStep > 0.0)
0192     {
0193         stepEstimate_Unconstrained =
0194             stepTrialOld * std::sqrt(fDeltaChord / dChordStep);
0195         stepTrial = fFractionNextEstimate * stepEstimate_Unconstrained;
0196     }
0197     else
0198     {
0199         // Should not update the Unconstrained Step estimate: incorrect!
0200         stepTrial =  stepTrialOld * 2.; 
0201     }
0202 
0203     if (stepTrial <= 0.001 * stepTrialOld)
0204     {
0205         if (dChordStep > 1000.0 * fDeltaChord)
0206         {
0207             stepTrial = stepTrialOld * 0.03;   
0208         }
0209         else
0210         {
0211             if (dChordStep > 100. * fDeltaChord)
0212             {
0213                 stepTrial = stepTrialOld * 0.1;   
0214             }
0215             else   // Try halving the length until dChordStep OK
0216             {
0217                 stepTrial = stepTrialOld * 0.5;   
0218             }
0219         }
0220     }
0221     else if (stepTrial > 1000.0 * stepTrialOld)
0222     {
0223         stepTrial = 1000.0 * stepTrialOld;
0224     }
0225 
0226     if (stepTrial == 0.0)
0227     {
0228         stepTrial= 0.000001;
0229     }
0230 
0231     // A more sophisticated chord-finder could figure out a better
0232     // stepTrial, from dChordStep and the required d_geometry
0233     //   e.g.
0234     //      Calculate R, r_helix (eg at orig point)
0235     //      if( stepTrial < 2 pi  R )
0236     //          stepTrial = R arc_cos( 1 - fDeltaChord / r_helix )
0237     //      else    
0238     //          ??
0239 
0240     return stepTrial;
0241 }
0242 
0243 template <class T>
0244 void G4ChordFinderDelegate<T>::AccumulateStatistics(G4int noTrials) 
0245 {
0246     fTotalNoTrials += noTrials; 
0247     ++fNoCalls; 
0248       
0249     if (noTrials > fmaxTrials) 
0250     { 
0251         fmaxTrials = noTrials; 
0252     }
0253 }
0254 
0255 template <class T>
0256 void G4ChordFinderDelegate<T>::PrintStatistics()
0257 {
0258     // Print Statistics
0259     G4cout << "G4ChordFinder statistics report: \n"
0260            << "  No trials: " << fTotalNoTrials
0261            << "  No Calls: "  << fNoCalls
0262            << "  Max-trial: " <<  fmaxTrials << "\n"
0263            << "  Parameters: " 
0264            << "  fFirstFraction "  << fFirstFraction
0265            << "  fFractionLast "   << fFractionLast
0266            << "  fFractionNextEstimate " << fFractionNextEstimate
0267            << G4endl; 
0268 }
0269 
0270 template <class T>
0271 G4int G4ChordFinderDelegate<T>::GetNoCalls()
0272 {
0273     return fNoCalls;
0274 }
0275 
0276 template <class T>
0277 G4int G4ChordFinderDelegate<T>::GetNoTrials()
0278 {
0279     return fTotalNoTrials;
0280 }
0281 
0282 template <class T>
0283 G4int G4ChordFinderDelegate<T>::GetNoMaxTrials()
0284 {
0285     return fmaxTrials;
0286 }
0287 
0288 template <class T>
0289 void G4ChordFinderDelegate<T>::SetFractions_Last_Next(G4double fractLast, 
0290                                                       G4double fractNext)
0291 {
0292     // Use -1.0 as request for Default.
0293     if (fractLast == -1.0) { fractLast = 1.0; }  // 0.9;
0294     if (fractNext == -1.0) { fractNext = 0.98; } // 0.9; 
0295 
0296     // fFirstFraction  = 0.999; // Safe value, range: ~ 0.95 - 0.999
0297     if (GetDriver().GetVerboseLevel() > 0)
0298     { 
0299         G4cout << " ChordFnd> Trying to set fractions: "
0300                << " first " << fFirstFraction
0301                << " last " <<  fractLast
0302                << " next " <<  fractNext
0303                << G4endl;
0304     } 
0305 
0306     if (fractLast > 0 && fractLast <= 1) 
0307     {
0308         fFractionLast = fractLast;
0309     } else
0310     {
0311         std::ostringstream message;
0312         message << "Invalid fraction Last = " << fractLast
0313                 << "; must be  0 <  fractionLast <= 1 ";
0314         G4Exception("G4ChordFinderDelegate::SetFractions_Last_Next()",
0315                     "GeomField1001", JustWarning, message);
0316     }
0317     if (fractNext > 0. && fractNext < 1)
0318     {
0319         fFractionNextEstimate = fractNext;
0320     } else
0321     {
0322         std::ostringstream message;
0323         message << "Invalid fraction Next = " << fractNext
0324                 << "; must be  0 <  fractionNext < 1 ";
0325         G4Exception("G4ChordFinderDelegate::SetFractions_Last_Next()",
0326                     "GeomField1001", JustWarning, message);
0327     }
0328 }
0329 
0330 template <class T>
0331 void G4ChordFinderDelegate<T>::SetFirstFraction(G4double fractFirst)
0332 {
0333     fFirstFraction = fractFirst;
0334 }
0335 
0336 template <class T>
0337 G4double G4ChordFinderDelegate<T>::GetFirstFraction()
0338 {
0339     return fFirstFraction;
0340 }
0341 
0342 template <class T>
0343 G4double G4ChordFinderDelegate<T>::GetFractionLast()
0344 {
0345     return fFractionLast;
0346 }
0347 
0348 template <class T>
0349 G4double G4ChordFinderDelegate<T>::GetFractionNextEstimate()
0350 {
0351     return fFractionNextEstimate;
0352 }
0353 
0354 template <class T>
0355 G4double G4ChordFinderDelegate<T>::GetLastStepEstimateUnc()
0356 {
0357     return fLastStepEstimate_Unconstrained;   
0358 }
0359 
0360 template <class T>
0361 void G4ChordFinderDelegate<T>::SetLastStepEstimateUnc(G4double stepEst)
0362 {
0363     fLastStepEstimate_Unconstrained = stepEst;
0364 }
0365 
0366 template <class T>
0367 void G4ChordFinderDelegate<T>::TestChordPrint(G4int noTrials, 
0368                                               G4int lastStepTrial, 
0369                                               G4double dChordStep,
0370                                               G4double fDeltaChord,
0371                                               G4double nextStepTrial)
0372 {
0373      G4int oldprec = G4cout.precision(5);
0374      G4cout << " ChF/fnc: notrial " << std::setw( 3) << noTrials 
0375             << " this_step= "       << std::setw(10) << lastStepTrial;
0376      if( std::fabs( (dChordStep / fDeltaChord) - 1.0 ) < 0.001 )
0377      {
0378        G4cout.precision(8);
0379      }
0380      else
0381      {
0382        G4cout.precision(6);
0383      }
0384      G4cout << " dChordStep=  " << std::setw(12) << dChordStep;
0385      if( dChordStep > fDeltaChord ) { G4cout << " d+"; }
0386      else                           { G4cout << " d-"; }
0387      G4cout.precision(5);
0388      G4cout <<  " new_step= "       << std::setw(10)
0389             << fLastStepEstimate_Unconstrained
0390             << " new_step_constr= " << std::setw(10)
0391             << lastStepTrial << G4endl;
0392      G4cout << " nextStepTrial = " << std::setw(10) << nextStepTrial << G4endl;
0393      G4cout.precision(oldprec);
0394 }
0395 
0396 template <class T>
0397 void G4ChordFinderDelegate<T>::StreamDelegateInfo( std::ostream& os ) const
0398 {
0399 // Write out the parameters / state of the driver
0400   os << "State of G4ChordFinderDelegate: " << std::endl;
0401   os << "--Parameters: " << std::endl;
0402   os << "    First Fraction = " << fFirstFraction << std::endl;
0403   os << "    Last  Fraction = " << fFractionLast << std::endl;
0404   os << "    Fract Next est = " << fFractionNextEstimate << std::endl;
0405 
0406   os << "--State (fungible): " << std::endl;
0407   os << "    Maximum No Trials (seen)         = " << fmaxTrials  << std::endl;  
0408   os << "    LastStepEstimate (Unconstrained) = " << fLastStepEstimate_Unconstrained
0409      << std::endl;
0410   // os << "    Statistics NOT printed. " << std::endl;
0411   os << "--Statistics: trials= " << fTotalNoTrials
0412      << "  calls= " << fNoCalls << std::endl;
0413 }
0414
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!