Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2026-10-05 09:14:01 

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