Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/magneticfield/src/G4BFieldIntegrationDriver.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /geometry/magneticfield/src/G4BFieldIntegrationDriver.cc (Version 11.3.0) and /geometry/magneticfield/src/G4BFieldIntegrationDriver.cc (Version 10.6.p1)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
 26 // G4BFieldIntegrationDriver implementation        26 // G4BFieldIntegrationDriver implementation
 27 //                                                 27 //
 28 // Specialized integration driver for pure mag     28 // Specialized integration driver for pure magnetic field
 29 //                                                 29 //
 30 // Author: D.Sorokin                               30 // Author: D.Sorokin
 31 // -------------------------------------------     31 // --------------------------------------------------------------------
 32                                                    32 
 33 #include "G4BFieldIntegrationDriver.hh"            33 #include "G4BFieldIntegrationDriver.hh"
 34                                                    34 
 35 #include "G4FieldTrack.hh"                         35 #include "G4FieldTrack.hh"
 36 #include "G4FieldUtils.hh"                         36 #include "G4FieldUtils.hh"
 37 #include "G4Exception.hh"                          37 #include "G4Exception.hh"
 38 #include "G4PhysicalConstants.hh"                  38 #include "G4PhysicalConstants.hh"
 39 #include "G4SystemOfUnits.hh"                      39 #include "G4SystemOfUnits.hh"
 40 #include "templates.hh"                            40 #include "templates.hh"
 41                                                    41 
 42                                                    42 
 43 namespace {                                        43 namespace {
 44                                                    44 
 45 G4Mag_EqRhs* toMagneticEquation(G4EquationOfMo     45 G4Mag_EqRhs* toMagneticEquation(G4EquationOfMotion* equation)
 46 {                                                  46 {
 47     auto e = dynamic_cast<G4Mag_EqRhs*>(equati     47     auto e = dynamic_cast<G4Mag_EqRhs*>(equation);
 48                                                    48 
 49     if (e == nullptr)                          <<  49     if (!e) 
 50     {                                              50     {
 51         G4Exception("G4BFieldIntegrationDriver     51         G4Exception("G4BFieldIntegrationDriver::G4BFieldIntegrationDriver",
 52                     "GeomField0003", FatalErro     52                     "GeomField0003", FatalErrorInArgument,
 53                     "Works only with G4Mag_EqR     53                     "Works only with G4Mag_EqRhs");
 54     }                                              54     }
 55                                                    55 
 56     return e;                                      56     return e;
 57 }                                                  57 }
 58                                                    58 
 59 } // namespace                                     59 } // namespace
 60                                                    60 
 61                                                    61 
 62 G4BFieldIntegrationDriver::G4BFieldIntegration     62 G4BFieldIntegrationDriver::G4BFieldIntegrationDriver(
 63     std::unique_ptr<G4VIntegrationDriver> smal     63     std::unique_ptr<G4VIntegrationDriver> smallStepDriver, 
 64     std::unique_ptr<G4VIntegrationDriver> larg     64     std::unique_ptr<G4VIntegrationDriver> largeStepDriver)
 65     : fSmallStepDriver(std::move(smallStepDriv     65     : fSmallStepDriver(std::move(smallStepDriver)),
 66       fLargeStepDriver(std::move(largeStepDriv     66       fLargeStepDriver(std::move(largeStepDriver)),
 67       fCurrDriver(fSmallStepDriver.get()),         67       fCurrDriver(fSmallStepDriver.get()),
 68       fEquation(toMagneticEquation(fCurrDriver     68       fEquation(toMagneticEquation(fCurrDriver->GetEquationOfMotion()))
 69 {                                                  69 {
 70     if (fSmallStepDriver->GetEquationOfMotion(     70     if (fSmallStepDriver->GetEquationOfMotion()
 71      != fLargeStepDriver->GetEquationOfMotion(     71      != fLargeStepDriver->GetEquationOfMotion())
 72     {                                              72     {
 73         G4Exception("G4BFieldIntegrationDriver     73         G4Exception("G4BFieldIntegrationDriver Constructor:",
 74                     "GeomField1001", FatalExce     74                     "GeomField1001", FatalException, "different EoM");  
 75     }                                              75     }
 76 }                                                  76 }
 77                                                    77 
 78 G4double G4BFieldIntegrationDriver::AdvanceCho     78 G4double G4BFieldIntegrationDriver::AdvanceChordLimited(G4FieldTrack& yCurrent, 
 79                                                    79                                                         G4double stepMax, 
 80                                                    80                                                         G4double epsStep, 
 81                                                    81                                                         G4double chordDistance)
 82 {                                                  82 {
 83     const G4double radius = CurvatureRadius(yC     83     const G4double radius = CurvatureRadius(yCurrent);
 84                                                    84 
 85     G4VIntegrationDriver* driver = nullptr;        85     G4VIntegrationDriver* driver = nullptr;
 86     if (chordDistance < 2 * radius)                86     if (chordDistance < 2 * radius)
 87     {                                              87     {
 88         stepMax = std::min(stepMax, twopi * ra     88         stepMax = std::min(stepMax, twopi * radius);
 89         driver = fSmallStepDriver.get();           89         driver = fSmallStepDriver.get();
 90         ++fSmallDriverSteps;                       90         ++fSmallDriverSteps;
 91     } else                                         91     } else
 92     {                                              92     {
 93         driver = fLargeStepDriver.get();           93         driver = fLargeStepDriver.get();
 94         ++fLargeDriverSteps;                       94         ++fLargeDriverSteps;
 95     }                                              95     }
 96                                                    96 
 97     if (driver != fCurrDriver)                     97     if (driver != fCurrDriver)
 98     {                                              98     {
 99         driver->OnComputeStep();                   99         driver->OnComputeStep();
100     }                                             100     }
101                                                   101 
102     fCurrDriver = driver;                         102     fCurrDriver = driver;
103                                                   103 
104     return fCurrDriver->AdvanceChordLimited(yC    104     return fCurrDriver->AdvanceChordLimited(yCurrent, stepMax,
105                                             ep    105                                             epsStep, chordDistance);
106 }                                                 106 }
107                                                   107 
108 void                                              108 void
109 G4BFieldIntegrationDriver::SetEquationOfMotion    109 G4BFieldIntegrationDriver::SetEquationOfMotion(G4EquationOfMotion* equation)
110 {                                                 110 {
111     fEquation = toMagneticEquation(equation);     111     fEquation = toMagneticEquation(equation);
112     fSmallStepDriver->SetEquationOfMotion(equa    112     fSmallStepDriver->SetEquationOfMotion(equation);
113     fLargeStepDriver->SetEquationOfMotion(equa    113     fLargeStepDriver->SetEquationOfMotion(equation);
114 }                                                 114 }
115                                                   115 
116 G4double                                          116 G4double
117 G4BFieldIntegrationDriver::CurvatureRadius(con    117 G4BFieldIntegrationDriver::CurvatureRadius(const G4FieldTrack& track) const
118 {                                                 118 {
119     G4double field[G4Field::MAX_NUMBER_OF_COMP    119     G4double field[G4Field::MAX_NUMBER_OF_COMPONENTS];
120                                                   120     
121     GetFieldValue(track, field);                  121     GetFieldValue(track, field);
122                                                   122     
123     const G4double Bmag2 =   field[0] * field[    123     const G4double Bmag2 =   field[0] * field[0]
124                            + field[1] * field[    124                            + field[1] * field[1]
125                            + field[2] * field[    125                            + field[2] * field[2] ;
126     if (Bmag2 == 0.0 )                            126     if (Bmag2 == 0.0 )
127     {                                             127     {
128        return DBL_MAX;                            128        return DBL_MAX;
129     }                                             129     }
130                                                   130 
131     const G4double momentum2 = track.GetMoment    131     const G4double momentum2 = track.GetMomentum().mag2();
132     const G4double fCof_inv = eplus / std::abs    132     const G4double fCof_inv = eplus / std::abs(fEquation->FCof());
133                                                   133 
134     return std::sqrt(momentum2 / Bmag2) * fCof    134     return std::sqrt(momentum2 / Bmag2) * fCof_inv;
135 }                                                 135 }
136                                                   136 
137 void                                              137 void
138 G4BFieldIntegrationDriver::GetFieldValue(const    138 G4BFieldIntegrationDriver::GetFieldValue(const G4FieldTrack& track,
139                                          G4dou    139                                          G4double Field[] ) const
140 {                                                 140 {
141     G4ThreeVector pos= track.GetPosition();       141     G4ThreeVector pos= track.GetPosition();
142     G4double positionTime[4]= { pos.x(), pos.y    142     G4double positionTime[4]= { pos.x(), pos.y(), pos.z(),
143                                 track.GetLabTi    143                                 track.GetLabTimeOfFlight() } ;
144                                                   144     
145     fEquation->GetFieldValue(positionTime, Fie    145     fEquation->GetFieldValue(positionTime, Field);
146 }                                                 146 }
147                                                   147 
148 void G4BFieldIntegrationDriver::PrintStatistic    148 void G4BFieldIntegrationDriver::PrintStatistics() const
149 {                                                 149 {
150     const auto totSteps = fSmallDriverSteps +     150     const auto totSteps = fSmallDriverSteps + fLargeDriverSteps;
151     const auto toFraction = [&](double value)     151     const auto toFraction = [&](double value) { return value / totSteps * 100; };
152                                                   152 
153     G4cout << "============= G4BFieldIntegrati    153     G4cout << "============= G4BFieldIntegrationDriver statistics ===========\n"
154            << "total steps " << totSteps << "     154            << "total steps " << totSteps << " "
155            << "smallDriverSteps " << toFractio    155            << "smallDriverSteps " << toFraction(fSmallDriverSteps) << " "
156            << "largeDriverSteps " << toFractio    156            << "largeDriverSteps " << toFraction(fLargeDriverSteps) << "\n"
157            << "===============================    157            << "======================================\n";
158 }                                                 158 }
159                                                   159