diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index fbfec5b1..6a4fa49d 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -17,6 +17,7 @@ add_subdirectory(Example11) add_subdirectory(Example12) add_subdirectory(Example13) add_subdirectory(Example14) +add_subdirectory(Example15) add_subdirectory(Example16) add_subdirectory(Example17) add_subdirectory(Example18) diff --git a/examples/Example11/README.md b/examples/Example11/README.md index 91c3e8ff..ebd0e79d 100644 --- a/examples/Example11/README.md +++ b/examples/Example11/README.md @@ -19,7 +19,7 @@ Results are reproducible using one RANLUX++ state per track. ### Kernels This example uses one stream per particle type to launch kernels asynchronously. -They are synchronized via a forth stream using CUDA events. +They are synchronized via a fourth stream using CUDA events. #### `TransportElectrons` diff --git a/examples/Example12/README.md b/examples/Example12/README.md index d282f472..6a55a190 100644 --- a/examples/Example12/README.md +++ b/examples/Example12/README.md @@ -22,7 +22,7 @@ Additionally, the kernels score energy deposit and the charged track length per ### Kernels This example uses one stream per particle type to launch kernels asynchronously. -They are synchronized via a forth stream using CUDA events. +They are synchronized via a fourth stream using CUDA events. #### `TransportElectrons` diff --git a/examples/Example13/README.md b/examples/Example13/README.md index 68a7d59c..b522e03c 100644 --- a/examples/Example13/README.md +++ b/examples/Example13/README.md @@ -24,7 +24,7 @@ Additionally, the kernels score energy deposit and the charged track length per ### Kernels This example uses one stream per particle type to launch kernels asynchronously. -They are synchronized via a forth stream using CUDA events. +They are synchronized via a fourth stream using CUDA events. #### `TransportElectrons` diff --git a/examples/Example13/electrons.cu b/examples/Example13/electrons.cu index a6b9bce1..dcd0b5a4 100644 --- a/examples/Example13/electrons.cu +++ b/examples/Example13/electrons.cu @@ -123,9 +123,12 @@ static __device__ __forceinline__ void TransportElectrons(Track *electrons, cons bool propagated = true; double geometryStepLength; vecgeom::NavStateIndex nextState; + constexpr int max_iters = 10; + // printf("max_iters= %3d\n", max_iters); if (BzFieldValue != 0) { geometryStepLength = fieldPropagatorBz.ComputeStepAndNextVolume( - energy, Mass, Charge, geometricalStepLengthFromPhysics, pos, dir, navState, nextState, propagated, safety); + energy, Mass, Charge, geometricalStepLengthFromPhysics, pos, dir, navState, nextState, + propagated, safety, max_iters ); } else { geometryStepLength = BVHNavigator::ComputeStepAndNextVolume(pos, dir, geometricalStepLengthFromPhysics, navState, nextState, kPush); diff --git a/examples/Example13/example13.cuh b/examples/Example13/example13.cuh index 7345c96f..bfcb3682 100644 --- a/examples/Example13/example13.cuh +++ b/examples/Example13/example13.cuh @@ -135,7 +135,7 @@ extern __constant__ __device__ struct G4HepEmData g4HepEmData; extern __constant__ __device__ int *MCIndex; -// constexpr vecgeom::Precision BzFieldValue = 0.1 * copcore::units::tesla; -constexpr vecgeom::Precision BzFieldValue = 0; +constexpr vecgeom::Precision BzFieldValue = 3.8 * copcore::units::tesla; +// constexpr vecgeom::Precision BzFieldValue = 0; #endif diff --git a/examples/Example15/CMakeLists.txt b/examples/Example15/CMakeLists.txt new file mode 100644 index 00000000..4ffd29dc --- /dev/null +++ b/examples/Example15/CMakeLists.txt @@ -0,0 +1,42 @@ +# SPDX-FileCopyrightText: 2021 CERN +# SPDX-License-Identifier: Apache-2.0 + +if(NOT TARGET G4HepEm::g4HepEm) + message(STATUS "Disabling example15 (needs G4HepEm)") + return() +endif() + +if(Geant4_FOUND) + if(NOT Geant4_gdml_FOUND) + message(STATUS "Disabling example15 (needs Geant4 with GDML support)") + return() + endif() +else() + message(STATUS "Disabling example15 (needs Geant4)") + return() +endif() + +# example15 is the AdePT demo example using material cuts as defined in the input gdml file +add_executable(example15 + example15.cpp + example15.cu + electrons.cu + gammas.cu) +target_link_libraries(example15 + PRIVATE + AdePT + CopCore::CopCore + VecGeom::vecgeom + VecGeom::vecgeomcuda_static + VecGeom::vgdml + ${Geant4_LIBRARIES} + G4HepEm::g4HepEmData + G4HepEm::g4HepEmInit + G4HepEm::g4HepEmRun + CUDA::cudart) + +set_target_properties(example15 PROPERTIES CUDA_SEPARABLE_COMPILATION ON CUDA_RESOLVE_DEVICE_SYMBOLS ON) + +# Tests +add_test(NAME example15 + COMMAND $ -gdml_file ${CMAKE_BINARY_DIR}/cms2018.gdml) diff --git a/examples/Example15/README.md b/examples/Example15/README.md new file mode 100644 index 00000000..a68335b7 --- /dev/null +++ b/examples/Example15/README.md @@ -0,0 +1,64 @@ + + +## Example 14 + +Example demonstrating particle transportation on GPUs with a non-uniform magnetic field - in arbitrary geometry read from a GDML file. + +New feature is + - definition of a non-uniform magnetic field + - integration of the equation of motion of charged particles using embedded Runge-Kutta method, Dormand Prince 4(5). + Caveats (at present) : + - fixed accuracy of integration (constant in source) + +The features carried over from Example13 are: + * arbitrary geometry via gdml file (tested with cms2018.gdml from VecGeom persistency/gdml/gdmls folder) and optionally a magnetic field with constant Bz, + * geometry read as Geant4 geometry, reading in regions and cuts, to initialize G4HepEm data + * geometry read then into VecGeom, and synchronized to GPU + * G4HepEm material-cuts couple indices mapped to VecGeom logical volume id's + * physics processes for e-/e+ (including MSC) and gammas using G4HepEm. + * scoring per placed volume, no sensitive detector feature + +Electrons, positrons, and gammas are stored in separate containers in device memory. +Free positions in the storage are handed out with monotonic slot numbers, slots are not reused. +Active tracks are passed via three queues per particle type (see `struct ParticleQueues`). +Results are reproducible using one RANLUX++ state per track. + +Additionally, the kernels score energy deposit and the charged track length per volume. + +### Kernels + +This example uses one stream per particle type to launch kernels asynchronously. +They are synchronized via a fourth stream using CUDA events. + +#### `TransportElectrons` + +1. Obtain safety unless the track is currently on a boundary. +2. Determine physics step limit, including conversion to geometric step length according to MSC. +3. Query geometry (or optionally magnetic field) to get geometry step length. +4. Convert geometry to true step length according to MSC, apply net direction change and discplacement. +5. Apply continuous effects; kill track if stopped. +6. If the particle reaches a boundary, perform relocation. +7. If not, and if there is a discrete process: + 1. Sample the final state. + 2. Update the primary and produce secondaries. + +#### `TransportGammas` + +1. Determine the physics step limit. +2. Query VecGeom to get geometry step length (no magnetic field for neutral particles!). +3. If the particle reaches a boundary, perform relocation. +4. If not, and if there is a discrete process: + 1. Sample the final state. + 2. Update the primary and produce secondaries. + +#### `FinishIteration` + +Clear the queues and return the tracks in flight. +This kernel runs after all secondary particles were produced. + +#### `InitPrimaries` and `InitParticleQueues` + +Used to initialize multiple primary particles with separate seeds. diff --git a/examples/Example15/electrons.cu b/examples/Example15/electrons.cu new file mode 100644 index 00000000..3a90c780 --- /dev/null +++ b/examples/Example15/electrons.cu @@ -0,0 +1,491 @@ +// SPDX-FileCopyrightText: 2022 CERN +// SPDX-License-Identifier: Apache-2.0 + +#include "example15.cuh" + +#include + +#include + +#define NOFLUCTUATION + +// Classes for Runge-Kutta integration +#include "MagneticFieldEquation.h" +#include "DormandPrinceRK45.h" +#include "fieldPropagatorRungeKutta.h" +#include "fieldPropagatorConstBz.h" + +#include +#include +#include +#include +#include +#include +// Pull in implementation. +#include +#include +#include +#include +#include +#include +#include + +#ifdef CHECK_RESULTS +#include "CompareResponses.h" +#endif + +__device__ float *gPtrBzFieldValue_dev = nullptr; + +// Transfer pointer to memory address of BzFieldValue_dev to device ... +// +__global__ void SetBzFieldPtr( float* pBzFieldValue_dev ) +{ + gPtrBzFieldValue_dev = pBzFieldValue_dev; +} + +// Compute the physics and geometry step limit, transport the electrons while +// applying the continuous effects and maybe a discrete process that could +// generate secondaries. +template +static __device__ __forceinline__ void TransportElectrons(Track *electrons, const adept::MParray *active, + Secondaries &secondaries, adept::MParray *activeQueue, + GlobalScoring *globalScoring, + ScoringPerVolume *scoringPerVolume) +{ +#ifdef VECGEOM_FLOAT_PRECISION + const Precision kPush = 10 * vecgeom::kTolerance; +#else + const Precision kPush = 0.; +#endif + constexpr int Charge = IsElectron ? -1 : 1; + constexpr double Mass = copcore::units::kElectronMassC2; + + constexpr int Nvar = 6; + using Field_t = UniformMagneticField; // ToDO: Change to non-uniform type !! + using Equation_t = MagneticFieldEquation; + using Stepper_t = DormandPrinceRK45; + using DoPri5Driver_t = RkIntegrationDriver; + + Field_t magField( vecgeom::Vector3D(0.0, 0.0, *gPtrBzFieldValue_dev) ); + // 2.0*copcore::units::tesla) ); // -> Obtain it from object ? + +#ifdef REPORT_OPTION + static bool ReportOption = true; + static const char* RunType= "Runge-Kutta field propagation"; + if( ReportOption && blockIdx.x == 0 && threadIdx.x == 0 ) { + printf( "-- Run type: %s .\n\n", RunType ); + ReportOption= false; + } +#endif + + // DoPri5Driver_t + // Static method fieldPropagatorRungeKutta + // no object fieldPropagatorRK() + + int activeSize = active->size(); + for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < activeSize; i += blockDim.x * gridDim.x) { + const int slot = (*active)[i]; + Track ¤tTrack = electrons[slot]; + auto energy = currentTrack.energy; + auto pos = currentTrack.pos; + auto dir = currentTrack.dir; + auto navState = currentTrack.navState; + const auto volume = navState.Top(); + const int volumeID = volume->id(); + // the MCC vector is indexed by the logical volume id + const int theMCIndex = MCIndex[volume->GetLogicalVolume()->id()]; + + auto survive = [&] { + currentTrack.energy = energy; + currentTrack.pos = pos; + currentTrack.dir = dir; + currentTrack.navState = navState; + activeQueue->push_back(slot); + }; + + // Init a track with the needed data to call into G4HepEm. + G4HepEmElectronTrack elTrack; + G4HepEmTrack *theTrack = elTrack.GetTrack(); + theTrack->SetEKin(energy); + theTrack->SetMCIndex(theMCIndex); + theTrack->SetOnBoundary(navState.IsOnBoundary()); + theTrack->SetCharge(Charge); + G4HepEmMSCTrackData *mscData = elTrack.GetMSCTrackData(); + mscData->fIsFirstStep = currentTrack.initialRange < 0; + mscData->fInitialRange = currentTrack.initialRange; + mscData->fDynamicRangeFactor = currentTrack.dynamicRangeFactor; + mscData->fTlimitMin = currentTrack.tlimitMin; + + // Initial value of magnetic field + // Equation_t::EvaluateDerivativesReturnB(magField, pos, + // , momentum* dir, + // , charge, dy_ds, magFieldStart); + + // Prepare a branched RNG state while threads are synchronized. Even if not + // used, this provides a fresh round of random numbers and reduces thread + // divergence because the RNG state doesn't need to be advanced later. + RanluxppDouble newRNG(currentTrack.rngState.BranchNoAdvance()); + + // Compute safety, needed for MSC step limit. + double safety = 0; + if (!navState.IsOnBoundary()) { + safety = BVHNavigator::ComputeSafety(pos, navState); + } + theTrack->SetSafety(safety); + + G4HepEmRandomEngine rnge(¤tTrack.rngState); + + // Sample the `number-of-interaction-left` and put it into the track. + for (int ip = 0; ip < 3; ++ip) { + double numIALeft = currentTrack.numIALeft[ip]; + if (numIALeft <= 0) { + numIALeft = -std::log(currentTrack.Uniform()); + } + theTrack->SetNumIALeft(numIALeft, ip); + } + + // Call G4HepEm to compute the physics step limit. + G4HepEmElectronManager::HowFar(&g4HepEmData, &g4HepEmPars, &elTrack, &rnge); + + // Remember MSC values for the next step(s). + currentTrack.initialRange = mscData->fInitialRange; + currentTrack.dynamicRangeFactor = mscData->fDynamicRangeFactor; + currentTrack.tlimitMin = mscData->fTlimitMin; + + // Get result into variables. + double geometricalStepLengthFromPhysics = theTrack->GetGStepLength(); + // The phyiscal step length is the amount that the particle experiences + // which might be longer than the geometrical step length due to MSC. As + // long as we call PerformContinuous in the same kernel we don't need to + // care, but we need to make this available when splitting the operations. + // double physicalStepLength = elTrack.GetPStepLength(); + int winnerProcessIndex = theTrack->GetWinnerProcessIndex(); + // Leave the range and MFP inside the G4HepEmTrack. If we split kernels, we + // also need to carry them over! + + // Check if there's a volume boundary in between. + bool propagated = true; + double geometryStepLength; + vecgeom::NavStateIndex nextState; + + float BzFieldValue= *gPtrBzFieldValue_dev; // Use vecgeom::Precision ? + if (BzFieldValue != 0.0) { + UniformMagneticField magneticFieldB( vecgeom::Vector3D(0.0, 0.0, BzFieldValue ) ); + // using fieldPropagatorRK = fieldPropagatorRungetKutta; + + // Set up the Integration using Runge-Kutta DoPri5 method + // + constexpr unsigned int Nvar = 6; // Number of integration variables + using Field_t = UniformMagneticField; + using Equation_t = MagneticFieldEquation; + using Stepper_t = DormandPrinceRK45; + using RkDriver_t = RkIntegrationDriver; + + constexpr int max_iterations= 10; + +#ifdef CHECK_RESULTS + // Store starting values + const vecgeom::Vector3D startPosition= pos; + const vecgeom::Vector3D startDirection= dir; + + // Start Baseline reply + vecgeom::Vector3D positionHx= startPosition; + vecgeom::Vector3D directionHx= startDirection; + vecgeom::NavStateIndex nextStateHx; + bool propagatedHx; + + fieldPropagatorConstBz fieldPropagatorBz(BzFieldValue); + Precision helixStepLength = fieldPropagatorBz.ComputeStepAndNextVolume( + energy, Mass, Charge, geometricalStepLengthFromPhysics, + positionHx, directionHx, navState, nextStateHx, propagatedHx, safety, + max_iterations ); + // activeSize < 100 ? max_iterations : max_iters_tail ); + // End Baseline reply +#endif + int iterDone= -1; + geometryStepLength = + fieldPropagatorRungeKutta::ComputeStepAndNextVolume( + magneticFieldB, + energy, Mass, Charge, geometricalStepLengthFromPhysics, + pos, dir, navState, nextState, + propagated, /*lengthDone,*/ safety, + // activeSize < 100 ? max_iterations : max_iters_tail ), // Was + max_iterations, + iterDone, slot + ); +#ifdef CHECK_RESULTS +#define formatBool(b) ((b) ? "yes " : "no") + constexpr Precision thresholdDiff=3.0e-3; + bool diffLength = false, badPosition = false, badDirection = false; + vecgeom::NavStateIndex& currNavState= navState; + bool sameLevel = nextState.GetLevel() == nextStateHx.GetLevel(); + bool sameIndex = nextState.GetNavIndex() == nextStateHx.GetNavIndex(); + + if( std::fabs( helixStepLength - geometryStepLength ) > 1.0e-4 * helixStepLength ) { + bool sameNextVol= (nextState.GetLevel() == nextStateHx.GetLevel()) && (nextStateHx.GetNavIndex() == nextStateHx.GetNavIndex() ) + && ( nextState.IsOnBoundary() == nextStateHx.IsOnBoundary() ); + printf ("\ns-len diff: id= %3d kinE= %12.7g phys-request= %11.6g helix-did= %11.6g rk-did= %11.6g (l-diff= %7.4g)" + " -- NavStates (curr/next RK/next Hx) : Levels %1d %1d %1d NavIdx: %5u %5u %5u OnBoundary: %3s %3s %3s Agree? %9s \n", + slot, energy, geometricalStepLengthFromPhysics, helixStepLength, geometryStepLength, geometryStepLength-helixStepLength, + navState.GetLevel(), nextState.GetLevel(), nextStateHx.GetLevel(), + navState.GetNavIndex(), nextState.GetNavIndex(), nextStateHx.GetNavIndex(), + // navState.IsOnBoundary()), nextState.IsOnBoundary(), nextStateHx.IsOnBoundary(), + formatBool(navState.IsOnBoundary()), formatBool(nextState.IsOnBoundary()), formatBool(nextStateHx.IsOnBoundary()), + ( sameNextVol ? "-Same-" : "-NotSame-" ) + ); + diffLength = true; + } else { + badPosition = + CompareResponseVector3D( slot, startPosition, positionHx, pos, "Position", thresholdDiff ); + badDirection = + CompareResponseVector3D( slot, startDirection, directionHx, dir, "Direction", thresholdDiff ); + } + const char* Outcome[2]={ "Good", " Bad" }; + bool problem = diffLength || badPosition || badDirection; + if( problem ) { + printf("%4s track (id= %3d) e_kin= %8.4g stepReq= %9.5g (did: RK= %9.5g vs hlx= %9.5g , diff= %9.5g) iters= %5d\n ", + Outcome[problem], // [diffLength||badPosition||badDirection], + slot, energy, geometricalStepLengthFromPhysics, geometryStepLength, helixStepLength, + geometryStepLength - helixStepLength, iterDone); + currentTrack.print(slot, /* verbose= */ true ); + } +#endif + } else { + geometryStepLength = BVHNavigator::ComputeStepAndNextVolume(pos, dir, geometricalStepLengthFromPhysics, navState, + nextState, kPush); + pos += geometryStepLength * dir; + } + + // Set boundary state in navState so the next step and secondaries get the + // correct information (navState = nextState only if relocated + // in case of a boundary; see below) + navState.SetBoundaryState(nextState.IsOnBoundary()); + + // Propagate information from geometrical step to MSC. + theTrack->SetDirection(dir.x(), dir.y(), dir.z()); + theTrack->SetGStepLength(geometryStepLength); + theTrack->SetOnBoundary(nextState.IsOnBoundary()); + + // Apply continuous effects. + bool stopped = G4HepEmElectronManager::PerformContinuous(&g4HepEmData, &g4HepEmPars, &elTrack, &rnge); + + // Collect the direction change and displacement by MSC. + const double *direction = theTrack->GetDirection(); + dir.Set(direction[0], direction[1], direction[2]); + if (!nextState.IsOnBoundary()) { + const double *mscDisplacement = mscData->GetDisplacement(); + vecgeom::Vector3D displacement(mscDisplacement[0], mscDisplacement[1], mscDisplacement[2]); + const double dLength2 = displacement.Length2(); + constexpr double kGeomMinLength = 5 * copcore::units::nm; // 0.05 [nm] + constexpr double kGeomMinLength2 = kGeomMinLength * kGeomMinLength; // (0.05 [nm])^2 + if (dLength2 > kGeomMinLength2) { + const double dispR = std::sqrt(dLength2); + // Estimate safety by subtracting the geometrical step length. + safety -= geometryStepLength; + constexpr double sFact = 0.99; + double reducedSafety = sFact * safety; + + // Apply displacement, depending on how close we are to a boundary. + // 1a. Far away from geometry boundary: + if (reducedSafety > 0.0 && dispR <= reducedSafety) { + pos += displacement; + } else { + // Recompute safety. + safety = BVHNavigator::ComputeSafety(pos, navState); + reducedSafety = sFact * safety; + + // 1b. Far away from geometry boundary: + if (reducedSafety > 0.0 && dispR <= reducedSafety) { + pos += displacement; + // 2. Push to boundary: + } else if (reducedSafety > kGeomMinLength) { + pos += displacement * (reducedSafety / dispR); + } + // 3. Very small safety: do nothing. + } + } + } + + // Collect the charged step length (might be changed by MSC). + atomicAdd(&globalScoring->chargedSteps, 1); + atomicAdd(&scoringPerVolume->chargedTrackLength[volumeID], elTrack.GetPStepLength()); + + // Collect the changes in energy and deposit. + energy = theTrack->GetEKin(); + double energyDeposit = theTrack->GetEnergyDeposit(); + atomicAdd(&globalScoring->energyDeposit, energyDeposit); + atomicAdd(&scoringPerVolume->energyDeposit[volumeID], energyDeposit); + + // Save the `number-of-interaction-left` in our track. + for (int ip = 0; ip < 3; ++ip) { + double numIALeft = theTrack->GetNumIALeft(ip); + currentTrack.numIALeft[ip] = numIALeft; + } + + if (stopped) { + if (!IsElectron) { + // Annihilate the stopped positron into two gammas heading to opposite + // directions (isotropic). + Track &gamma1 = secondaries.gammas.NextTrack(); + Track &gamma2 = secondaries.gammas.NextTrack(); + atomicAdd(&globalScoring->numGammas, 2); + + const double cost = 2 * currentTrack.Uniform() - 1; + const double sint = sqrt(1 - cost * cost); + const double phi = k2Pi * currentTrack.Uniform(); + double sinPhi, cosPhi; + sincos(phi, &sinPhi, &cosPhi); + + gamma1.InitAsSecondary(pos, navState); + newRNG.Advance(); + gamma1.rngState = newRNG; + gamma1.energy = copcore::units::kElectronMassC2; + gamma1.dir.Set(sint * cosPhi, sint * sinPhi, cost); + + gamma2.InitAsSecondary(pos, navState); + // Reuse the RNG state of the dying track. + gamma2.rngState = currentTrack.rngState; + gamma2.energy = copcore::units::kElectronMassC2; + gamma2.dir = -gamma1.dir; + } + // Particles are killed by not enqueuing them into the new activeQueue. + continue; + } + + if (nextState.IsOnBoundary()) { + // For now, just count that we hit something. + atomicAdd(&globalScoring->hits, 1); + + // Kill the particle if it left the world. + if (nextState.Top() != nullptr) { + BVHNavigator::RelocateToNextVolume(pos, dir, nextState); + + // Move to the next boundary. + navState = nextState; + survive(); + } + continue; + } else if (!propagated) { + // Did not yet reach the interaction point due to error in the magnetic + // field propagation. Try again next time. + survive(); + continue; + } else if (winnerProcessIndex < 0) { + // No discrete process, move on. + survive(); + continue; + } + + // Reset number of interaction left for the winner discrete process. + // (Will be resampled in the next iteration.) + currentTrack.numIALeft[winnerProcessIndex] = -1.0; + + // Check if a delta interaction happens instead of the real discrete process. + if (G4HepEmElectronManager::CheckDelta(&g4HepEmData, theTrack, currentTrack.Uniform())) { + // A delta interaction happened, move on. + survive(); + continue; + } + + // Perform the discrete interaction, make sure the branched RNG state is + // ready to be used. + newRNG.Advance(); + // Also advance the current RNG state to provide a fresh round of random + // numbers after MSC used up a fair share for sampling the displacement. + currentTrack.rngState.Advance(); + + const double theElCut = g4HepEmData.fTheMatCutData->fMatCutData[theMCIndex].fSecElProdCutE; + + switch (winnerProcessIndex) { + case 0: { + // Invoke ionization (for e-/e+): + double deltaEkin = (IsElectron) ? G4HepEmElectronInteractionIoni::SampleETransferMoller(theElCut, energy, &rnge) + : G4HepEmElectronInteractionIoni::SampleETransferBhabha(theElCut, energy, &rnge); + + double dirPrimary[] = {dir.x(), dir.y(), dir.z()}; + double dirSecondary[3]; + G4HepEmElectronInteractionIoni::SampleDirections(energy, deltaEkin, dirSecondary, dirPrimary, &rnge); + + Track &secondary = secondaries.electrons.NextTrack(); + atomicAdd(&globalScoring->numElectrons, 1); + + secondary.InitAsSecondary(pos, navState); + secondary.rngState = newRNG; + secondary.energy = deltaEkin; + secondary.dir.Set(dirSecondary[0], dirSecondary[1], dirSecondary[2]); + + energy -= deltaEkin; + dir.Set(dirPrimary[0], dirPrimary[1], dirPrimary[2]); + survive(); + break; + } + case 1: { + // Invoke model for Bremsstrahlung: either SB- or Rel-Brem. + double logEnergy = std::log(energy); + double deltaEkin = energy < g4HepEmPars.fElectronBremModelLim + ? G4HepEmElectronInteractionBrem::SampleETransferSB(&g4HepEmData, energy, logEnergy, + theMCIndex, &rnge, IsElectron) + : G4HepEmElectronInteractionBrem::SampleETransferRB(&g4HepEmData, energy, logEnergy, + theMCIndex, &rnge, IsElectron); + + double dirPrimary[] = {dir.x(), dir.y(), dir.z()}; + double dirSecondary[3]; + G4HepEmElectronInteractionBrem::SampleDirections(energy, deltaEkin, dirSecondary, dirPrimary, &rnge); + + Track &gamma = secondaries.gammas.NextTrack(); + atomicAdd(&globalScoring->numGammas, 1); + + gamma.InitAsSecondary(pos, navState); + gamma.rngState = newRNG; + gamma.energy = deltaEkin; + gamma.dir.Set(dirSecondary[0], dirSecondary[1], dirSecondary[2]); + + energy -= deltaEkin; + dir.Set(dirPrimary[0], dirPrimary[1], dirPrimary[2]); + survive(); + break; + } + case 2: { + // Invoke annihilation (in-flight) for e+ + double dirPrimary[] = {dir.x(), dir.y(), dir.z()}; + double theGamma1Ekin, theGamma2Ekin; + double theGamma1Dir[3], theGamma2Dir[3]; + G4HepEmPositronInteractionAnnihilation::SampleEnergyAndDirectionsInFlight( + energy, dirPrimary, &theGamma1Ekin, theGamma1Dir, &theGamma2Ekin, theGamma2Dir, &rnge); + + Track &gamma1 = secondaries.gammas.NextTrack(); + Track &gamma2 = secondaries.gammas.NextTrack(); + atomicAdd(&globalScoring->numGammas, 2); + + gamma1.InitAsSecondary(pos, navState); + gamma1.rngState = newRNG; + gamma1.energy = theGamma1Ekin; + gamma1.dir.Set(theGamma1Dir[0], theGamma1Dir[1], theGamma1Dir[2]); + + gamma2.InitAsSecondary(pos, navState); + // Reuse the RNG state of the dying track. + gamma2.rngState = currentTrack.rngState; + gamma2.energy = theGamma2Ekin; + gamma2.dir.Set(theGamma2Dir[0], theGamma2Dir[1], theGamma2Dir[2]); + + // The current track is killed by not enqueuing into the next activeQueue. + break; + } + } + } +} + +// Instantiate kernels for electrons and positrons. +__global__ void TransportElectrons(Track *electrons, const adept::MParray *active, Secondaries secondaries, + adept::MParray *activeQueue, GlobalScoring *globalScoring, + ScoringPerVolume *scoringPerVolume) +{ + TransportElectrons(electrons, active, secondaries, activeQueue, globalScoring, scoringPerVolume); +} +__global__ void TransportPositrons(Track *positrons, const adept::MParray *active, Secondaries secondaries, + adept::MParray *activeQueue, GlobalScoring *globalScoring, + ScoringPerVolume *scoringPerVolume) +{ + TransportElectrons(positrons, active, secondaries, activeQueue, globalScoring, + scoringPerVolume); +} diff --git a/examples/Example15/example15.cpp b/examples/Example15/example15.cpp new file mode 100644 index 00000000..22ab6066 --- /dev/null +++ b/examples/Example15/example15.cpp @@ -0,0 +1,301 @@ +// SPDX-FileCopyrightText: 2021 CERN +// SPDX-License-Identifier: Apache-2.0 + +#include "example15.h" + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include + +static constexpr double DefaultCut = 0.7 * mm; + +float BzFieldValue_host = 1.0 * copcore::units::tesla; +float *BzFieldValue_dev = nullptr; + +const G4VPhysicalVolume *InitGeant4(const std::string &gdml_file) +{ + auto defaultRegion = new G4Region("DefaultRegionForTheWorld"); // deleted by store + auto pcuts = G4ProductionCutsTable::GetProductionCutsTable()->GetDefaultProductionCuts(); + pcuts->SetProductionCut(DefaultCut, "gamma"); + pcuts->SetProductionCut(DefaultCut, "e-"); + pcuts->SetProductionCut(DefaultCut, "e+"); + pcuts->SetProductionCut(DefaultCut, "proton"); + defaultRegion->SetProductionCuts(pcuts); + + // Read the geometry, regions and cuts from the GDML file + std::cout << "reading " << gdml_file << " transiently on CPU for Geant4 ...\n"; + G4GDMLParser parser; + parser.Read(gdml_file, false); // turn off schema checker + G4VPhysicalVolume *world = parser.GetWorldVolume(); + + if (world == nullptr) { + std::cerr << "example15: World volume not set properly check your setup selection criteria or GDML input!\n"; + return world; + } + + // - PHYSICS + // --- Create particles that have secondary production threshold. + G4Gamma::Gamma(); + G4Electron::Electron(); + G4Positron::Positron(); + G4Proton::Proton(); + G4ParticleTable *partTable = G4ParticleTable::GetParticleTable(); + partTable->SetReadiness(); + + // - REGIONS + if (world->GetLogicalVolume()->GetRegion() == nullptr) { + // Add default region if none available + defaultRegion->AddRootLogicalVolume(world->GetLogicalVolume()); + } + for (auto region : *G4RegionStore::GetInstance()) { + region->UsedInMassGeometry(true); // make sure all regions are marked as used + region->UpdateMaterialList(); + } + + // - UPDATE COUPLES + std::cout << "updating material-cut couples based on " << G4RegionStore::GetInstance()->size() << " regions ...\n"; + G4ProductionCutsTable *theCoupleTable = G4ProductionCutsTable::GetProductionCutsTable(); + theCoupleTable->UpdateCoupleTable(world); + + // --- Set the world volume to fix initialization of G4SafetyHelper (used by G4UrbanMscModel) + G4TransportationManager::GetTransportationManager()->SetWorldForTracking(world); + + // --- Set MSC range factor to match G4HepEm physics lists. + G4EmParameters *param = G4EmParameters::Instance(); + param->SetDefaults(); + param->SetMscRangeFactor(0.04); + + return world; +} + +const vecgeom::cxx::VPlacedVolume *InitVecGeom(const std::string &gdml_file, int cache_depth) +{ + // Import the gdml file into VecGeom + vecgeom::GeoManager::Instance().SetTransformationCacheDepth(cache_depth); + vgdml::Parser vgdmlParser; + auto middleWare = vgdmlParser.Load(gdml_file.c_str(), false, copcore::units::mm); + if (middleWare == nullptr) { + std::cerr << "Failed to read geometry from GDML file '" << gdml_file << "'" << std::endl; + return nullptr; + } + + const vecgeom::VPlacedVolume *world = vecgeom::GeoManager::Instance().GetWorld(); + if (world == nullptr) { + std::cerr << "GeoManager world volume is nullptr" << std::endl; + return nullptr; + } + return world; +} + +int *CreateMCCindex(const G4VPhysicalVolume *g4world, const vecgeom::VPlacedVolume *world, + const G4HepEmState &hepEmState) +{ + const int *g4tohepmcindex = hepEmState.fData->fTheMatCutData->fG4MCIndexToHepEmMCIndex; + + // - FIND vecgeom::LogicalVolume corresponding to each and every G4LogicalVolume + int nphysical = 0; + + int nvolumes = vecgeom::GeoManager::Instance().GetRegisteredVolumesCount(); + int *mcc_index = new int[nvolumes]; + memset(mcc_index, 0, nvolumes * sizeof(int)); + + // recursive geometry visitor lambda matching one by one Geant4 and VecGeom logical volumes + // (we need to make sure we set the right MCC index to the right volume) + typedef std::function func_t; + func_t visitAndSetMCindex = [&](G4LogicalVolume const *g4vol, vecgeom::LogicalVolume const *vol) { + int nd = g4vol->GetNoDaughters(); + auto daughters = vol->GetDaughters(); + if (nd != daughters.size()) throw std::runtime_error("Mismatch in number of daughters"); + // Check the couples + if (g4vol->GetMaterialCutsCouple() == nullptr) + throw std::runtime_error("G4LogicalVolume " + std::string(g4vol->GetName()) + + std::string(" has no material-cuts couple")); + int g4mcindex = g4vol->GetMaterialCutsCouple()->GetIndex(); + int hepemmcindex = g4tohepmcindex[g4mcindex]; + // Check consistency with G4HepEm data + if (hepEmState.fData->fTheMatCutData->fMatCutData[hepemmcindex].fG4MatCutIndex != g4mcindex) + throw std::runtime_error("Mismatch between Geant4 mcindex and corresponding G4HepEm index"); + if (vol->id() >= nvolumes) throw std::runtime_error("Volume id larger than number of volumes"); + + // All OK, now fill the index in the array + mcc_index[vol->id()] = hepemmcindex; + nphysical++; + + // Now do the daughters + for (int id = 0; id < nd; ++id) { + auto g4pvol = g4vol->GetDaughter(id); + auto pvol = daughters[id]; + // VecGeom does not strip pointers from logical volume names + if (std::string(pvol->GetLogicalVolume()->GetName()).rfind(g4pvol->GetLogicalVolume()->GetName(), 0) != 0) + throw std::runtime_error("Volume names " + std::string(pvol->GetLogicalVolume()->GetName()) + " and " + + std::string(g4pvol->GetLogicalVolume()->GetName()) + " mismatch"); + visitAndSetMCindex(g4pvol->GetLogicalVolume(), pvol->GetLogicalVolume()); + } + }; + + visitAndSetMCindex(g4world->GetLogicalVolume(), world->GetLogicalVolume()); + std::cout << "Visited " << nphysical << " matching physical volumes\n"; + return mcc_index; +} + +void FreeG4HepEm(G4HepEmState *state) +{ + FreeG4HepEmData(state->fData); +} + +void InitBVH() +{ + vecgeom::cxx::BVHManager::Init(); + vecgeom::cxx::BVHManager::DeviceInit(); +} + +void PrintScoringPerVolume(const vecgeom::VPlacedVolume *placed, const ScoringPerVolume *scoring, int level = 0) +{ + for (auto *daughter : placed->GetDaughters()) { + std::cout << std::setw(level * 2) << ""; + auto id = daughter->id(); + std::cout << " ID " << id << " Charged-TrakL " << scoring->chargedTrackLength[id] / copcore::units::mm + << " mm; Energy-Dep " << scoring->energyDeposit[id] / copcore::units::MeV << " MeV" << std::endl; + PrintScoringPerVolume(daughter, scoring, level + 1); + } +} + +int main(int argc, char *argv[]) +{ + // Only outputs are the data file(s) + // Separate for now, but will want to unify + OPTION_STRING(gdml_file, "cms2018.gdml"); // "trackML.gdml"); // "cms2018.gdml"); + OPTION_INT(cache_depth, 0); // 0 = full depth + OPTION_INT(particles, 100); + OPTION_DOUBLE(energy, 10); // entered in GeV + energy *= copcore::units::GeV; + OPTION_INT(batch, -1); + OPTION_BOOL(rotatingParticleGun, false); + + OPTION_DOUBLE(bzValue, -1.00); // entered in tesla + BzFieldValue_host = bzValue * copcore::units::tesla; + + std::cout << "INFO: Particles to simulate = " << particles << " \n"; + std::cout << "INFO: Bz field value = " << bzValue << " T (tesla) \n"; + + constexpr int stackLimit= 8 * 1024; // 8192 + printf("example15.cpp/main(): Setting Cuda Device Stack Limit to %6d \n", stackLimit); + cudaError_t error = cudaDeviceSetLimit(cudaLimitStackSize, stackLimit); + if (error != cudaSuccess) { + printf("cudaDeviceSetLimit failed with %d, line(%d)\n", error, __LINE__); + exit(EXIT_FAILURE); + } + + vecgeom::Stopwatch timer; + timer.Start(); + + // Initialize Geant4 + auto g4world = InitGeant4(gdml_file); + if (!g4world) return 3; + + // Initialize VecGeom + std::cout << "reading " << gdml_file << " transiently on CPU for VecGeom ...\n"; + auto world = InitVecGeom(gdml_file, cache_depth); + if (!world) return 3; + + // Construct and initialize the G4HepEmState data/tables + std::cout << "initializing G4HepEm state ...\n"; + G4HepEmState hepEmState; + InitG4HepEmState(&hepEmState); + + // Initialize G4HepEm material-cut couple array indexed by VecGeom volume id. + // (In future we should connect the index directly to the VecGeom logical volume) + std::cout << "initializing material-cut couple indices ...\n"; + int *MCCindex = nullptr; + + try { + MCCindex = CreateMCCindex(g4world, world, hepEmState); + } catch (const std::runtime_error &ex) { + std::cerr << "*** CreateMCCindex: " << ex.what() << "\n"; + return 1; + } + + // Load and synchronize the geometry on the GPU + std::cout << "synchronizing VecGeom geometry to GPU ...\n"; + auto &cudaManager = vecgeom::cxx::CudaManager::Instance(); + cudaManager.LoadGeometry(world); + cudaManager.Synchronize(); + InitBVH(); + + auto time_cpu = timer.Stop(); + std::cout << "Initialization took: " << time_cpu << " sec\n"; + + int NumVolumes = vecgeom::GeoManager::Instance().GetRegisteredVolumesCount(); + int NumPlaced = vecgeom::GeoManager::Instance().GetPlacedVolumesCount(); + + // Scoring is done per placed volume (for now...) + double *chargedTrackLength = new double[NumPlaced]; + double *energyDeposit = new double[NumPlaced]; + ScoringPerVolume scoringPerVolume; + scoringPerVolume.chargedTrackLength = chargedTrackLength; + scoringPerVolume.energyDeposit = energyDeposit; + GlobalScoring globalScoring; + + example15(particles, energy, batch, MCCindex, &scoringPerVolume, &globalScoring, NumVolumes, NumPlaced, &hepEmState, + rotatingParticleGun); + + std::cout << std::endl; + std::cout << std::endl; + double meanEnergyDeposit = globalScoring.energyDeposit / particles; + std::cout << "Mean energy deposit " << (meanEnergyDeposit / copcore::units::GeV) << " GeV\n" + << "Mean number of gamma " << ((double)globalScoring.numGammas / particles) << "\n" + << "Mean number of e- " << ((double)globalScoring.numElectrons / particles) << "\n" + << "Mean number of e+ " << ((double)globalScoring.numPositrons / particles) << "\n" + << "Mean number of charged steps " << ((double)globalScoring.chargedSteps / particles) << "\n" + << "Mean number of neutral steps " << ((double)globalScoring.neutralSteps / particles) << "\n" + << "Mean number of hits " << ((double)globalScoring.hits / particles) << "\n"; + if (globalScoring.numKilled > 0) { + std::cout << "Total killed particles " << globalScoring.numKilled << "\n"; + } + std::cout << std::endl; + + // Average charged track length and energy deposit per particle. + for (int i = 0; i < NumPlaced; i++) { + chargedTrackLength[i] /= particles; + energyDeposit[i] /= particles; + } + + std::cout << std::scientific; +#ifdef VERBOSE + std::cout << std::endl; + const int id = world->id(); + std::cout << "World (ID " << world->id() << ") Charged-TrakL " << chargedTrackLength[id] / copcore::units::mm + << " Energy-Dep " << energyDeposit[id] / copcore::units::MeV << " MeV" << std::endl; + PrintScoringPerVolume(world, &scoringPerVolume); +#endif + FreeG4HepEm(&hepEmState); + delete[] MCCindex; + delete[] chargedTrackLength; + delete[] energyDeposit; + return 0; +} diff --git a/examples/Example15/example15.cu b/examples/Example15/example15.cu new file mode 100644 index 00000000..1fd91dc2 --- /dev/null +++ b/examples/Example15/example15.cu @@ -0,0 +1,552 @@ +// SPDX-FileCopyrightText: 2021 CERN +// SPDX-License-Identifier: Apache-2.0 + +#include "example15.h" +#include "example15.cuh" + +#include +#include +#include + +#include +#include +#include + +#include +#include +#include +#ifdef VECGEOM_ENABLE_CUDA +#include +#endif + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +__constant__ __device__ struct G4HepEmParameters g4HepEmPars; +__constant__ __device__ struct G4HepEmData g4HepEmData; + +__constant__ __device__ int *MCIndex = nullptr; + +void InitG4HepEmGPU(G4HepEmState *state) +{ + // Copy to GPU. + CopyG4HepEmDataToGPU(state->fData); + COPCORE_CUDA_CHECK(cudaMemcpyToSymbol(g4HepEmPars, state->fParameters, sizeof(G4HepEmParameters))); + + // Create G4HepEmData with the device pointers. + G4HepEmData dataOnDevice; + dataOnDevice.fTheMatCutData = state->fData->fTheMatCutData_gpu; + dataOnDevice.fTheMaterialData = state->fData->fTheMaterialData_gpu; + dataOnDevice.fTheElementData = state->fData->fTheElementData_gpu; + dataOnDevice.fTheElectronData = state->fData->fTheElectronData_gpu; + dataOnDevice.fThePositronData = state->fData->fThePositronData_gpu; + dataOnDevice.fTheSBTableData = state->fData->fTheSBTableData_gpu; + dataOnDevice.fTheGammaData = state->fData->fTheGammaData_gpu; + // The other pointers should never be used. + dataOnDevice.fTheMatCutData_gpu = nullptr; + dataOnDevice.fTheMaterialData_gpu = nullptr; + dataOnDevice.fTheElementData_gpu = nullptr; + dataOnDevice.fTheElectronData_gpu = nullptr; + dataOnDevice.fThePositronData_gpu = nullptr; + dataOnDevice.fTheSBTableData_gpu = nullptr; + dataOnDevice.fTheGammaData_gpu = nullptr; + + COPCORE_CUDA_CHECK(cudaMemcpyToSymbol(g4HepEmData, &dataOnDevice, sizeof(G4HepEmData))); +} + +// A bundle of queues per particle type: +// * Two for active particles, one for the current iteration and the second for the next. +struct ParticleQueues { + adept::MParray *currentlyActive; + adept::MParray *nextActive; + + void SwapActive() { std::swap(currentlyActive, nextActive); } +}; + +struct ParticleType { + Track *tracks; + SlotManager *slotManager; + ParticleQueues queues; + cudaStream_t stream; + cudaEvent_t event; + + enum { + Electron = 0, + Positron = 1, + Gamma = 2, + + NumParticleTypes, + }; +}; + +// A bundle of queues for the three particle types. +struct AllParticleQueues { + ParticleQueues queues[ParticleType::NumParticleTypes]; +}; + +// Kernel to initialize the set of queues per particle type. +__global__ void InitParticleQueues(ParticleQueues queues, size_t Capacity) +{ + adept::MParray::MakeInstanceAt(Capacity, queues.currentlyActive); + adept::MParray::MakeInstanceAt(Capacity, queues.nextActive); +} + +// Kernel function to initialize a set of primary particles. +__global__ void InitPrimaries(ParticleGenerator generator, int startEvent, int numEvents, double energy, + const vecgeom::VPlacedVolume *world, GlobalScoring *globalScoring, bool rotatingParticleGun) +{ + for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < numEvents; i += blockDim.x * gridDim.x) { + Track &track = generator.NextTrack(); + + track.rngState.SetSeed(startEvent + i); + track.energy = energy; + track.numIALeft[0] = -1.0; + track.numIALeft[1] = -1.0; + track.numIALeft[2] = -1.0; + + track.initialRange = -1.0; + track.dynamicRangeFactor = -1.0; + track.tlimitMin = -1.0; + + track.pos = {0, 0, 0}; + if (rotatingParticleGun) { + // Generate particles flat in phi and in eta between -5 and 5. We'll lose the far forwards ones, so no need to simulate. + const double phi = 2. * M_PI * track.rngState.Rndm(); + const double eta = -5. + 10. * track.rngState.Rndm(); + track.dir.x() = static_cast(cos(phi) / cosh(eta)); + track.dir.y() = static_cast(sin(phi) / cosh(eta)); + track.dir.z() = static_cast(tanh(eta)); + } else { + track.dir = {1.0, 0, 0}; + } + track.navState.Clear(); + BVHNavigator::LocatePointIn(world, track.pos, track.navState, true); + + atomicAdd(&globalScoring->numElectrons, 1); + } +} + +// A data structure to transfer statistics after each iteration. +struct Stats { + int inFlight[ParticleType::NumParticleTypes]; +}; + +// Finish iteration: clear queues and fill statistics. +__global__ void FinishIteration(AllParticleQueues all, Stats *stats) +{ + for (int i = 0; i < ParticleType::NumParticleTypes; i++) { + all.queues[i].currentlyActive->clear(); + stats->inFlight[i] = all.queues[i].nextActive->size(); + } +} + +__global__ void ClearQueue(adept::MParray *queue) +{ + queue->clear(); +} + +__global__ void printTrackV2(const Track* trkArray_dev, + const adept::MParray *active_dev, + int i, + bool verbose) +{ + const int slot = (*active_dev)[i]; + trkArray_dev[slot].print(slot, verbose); +} + +__global__ +void get_size( const adept::MParray *arr_dev, int* capacity_dev ) +{ + *capacity_dev = arr_dev->size(); + // printf( " Size of arr_dev %p = %d \n", arr_dev, *capacity_dev ); +} + +__host__ +int obtain_size( const adept::MParray *arr_dev ) +{ + static int* pCapacity_dev = nullptr; + int capacity_host; + + if( ! pCapacity_dev ) { + COPCORE_CUDA_CHECK(cudaMalloc(&pCapacity_dev, sizeof(int))); + std::cout << " Created memory space for capacity_dev. Ptr = " << pCapacity_dev << std::endl; + } + + // Fill it on device + get_size<<<1,1>>>(arr_dev, pCapacity_dev); + + // Fetch it from the device + COPCORE_CUDA_CHECK(cudaMemcpy( &capacity_host, pCapacity_dev, + sizeof(int), cudaMemcpyDeviceToHost) ); + return capacity_host; +} + +__host__ +void printTracks( const Track *trackArr_dev, + const adept::MParray *active_dev, + short particleTypeId, + bool verbose, + int numTracks ) +{ + static const char* particleTypeName[ParticleType::NumParticleTypes] = { "Electrons", "Positrons", "Gammas" }; + static const char* particleShortName[] = { "e-", "e+", "g" }; + // int numPrinted = 0; + int capacity_host= obtain_size( active_dev ); + + if( numTracks > 0 ) { + printf("%s : %d\n", particleTypeName[particleTypeId], numTracks ); + + for( int i= 0; i < capacity_host; i++ ){ + // Track &trk = *(trackArr+i); + printf(" %2s (pid= %2d) ", particleShortName[particleTypeId], particleTypeId ); + //if (numPrinted++ < numTracks ) + //{ + printTrackV2<<<1, 1>>>(trackArr_dev, active_dev, i, verbose); + // std::cout << "\n"; + // } + cudaDeviceSynchronize(); + } + // cudaDeviceSynchronize(); + } +} + +__host__ +void printActiveTracks( ParticleType& electrons, ParticleType& positrons, ParticleType& gammas, bool verbTrk ) +{ + int numElectrons = obtain_size(electrons.queues.currentlyActive); + int numPositrons = obtain_size(positrons.queues.currentlyActive); + int numGammas = obtain_size(gammas.queues.currentlyActive); + // printf("Electrons: \n"); + printTracks( electrons.tracks, electrons.queues.currentlyActive, ParticleType::Electron, /* "Electrons", */ verbTrk, numElectrons ); + // printf("Positrons: \n"); + printTracks( positrons.tracks, positrons.queues.currentlyActive, ParticleType::Positron, verbTrk, numPositrons ); + // printf("Gammas: \n"); + printTracks( gammas.tracks, gammas.queues.currentlyActive, ParticleType::Gamma, verbTrk, numGammas ); + std::cout << std::endl; +} + +// #define VERBOSE_REPORT 1 +// To obtain information both intermediately in an iteration and at the end + +void example15(int numParticles, double energy, int batch, const int *MCIndex_host, + ScoringPerVolume *scoringPerVolume_host, GlobalScoring *globalScoring_host, int numVolumes, + int numPlaced, G4HepEmState *state, bool rotatingParticleGun) +{ + // NVTXTracer tracer("InitG4HepEM"); + InitG4HepEmGPU(state); + + // Transfer MC indices. + int *MCIndex_dev = nullptr; + COPCORE_CUDA_CHECK(cudaMalloc(&MCIndex_dev, sizeof(int) * numVolumes)); + COPCORE_CUDA_CHECK(cudaMemcpy(MCIndex_dev, MCIndex_host, sizeof(int) * numVolumes, cudaMemcpyHostToDevice)); + COPCORE_CUDA_CHECK(cudaMemcpyToSymbol(MCIndex, &MCIndex_dev, sizeof(int *))); + + // Capacity of the different containers aka the maximum number of particles. + constexpr int Capacity = 1024 * 1024; + + std::cout << "INFO: capacity of containers set to " << Capacity << std::endl; + std::cout << "INFO: number of particles = " << numParticles << std::endl; + if (batch == -1) { + // Rule of thumb: at most 2000 particles of one type per GeV primary. + batch = Capacity / ((int)energy / copcore::units::GeV) / 2000; + } else if (batch < 1) { + batch = 1; + } + std::cout << "INFO: batching " << batch << " particles for transport on the GPU" << std::endl; + if (BzFieldValue_host != 0) { + std::cout << "INFO: running with field Bz = " << BzFieldValue_host / copcore::units::tesla << " T" << std::endl; + } else { + std::cout << "INFO: running with magnetic field OFF" << std::endl; + } + std::cout<>>(particles[i].queues, Capacity); + + COPCORE_CUDA_CHECK(cudaStreamCreate(&particles[i].stream)); + COPCORE_CUDA_CHECK(cudaEventCreate(&particles[i].event)); + } + COPCORE_CUDA_CHECK(cudaDeviceSynchronize()); + + ParticleType &electrons = particles[ParticleType::Electron]; + ParticleType &positrons = particles[ParticleType::Positron]; + ParticleType &gammas = particles[ParticleType::Gamma]; + + // Create a stream to synchronize kernels of all particle types. + cudaStream_t stream; + COPCORE_CUDA_CHECK(cudaStreamCreate(&stream)); + + // Allocate memory to score charged track length and energy deposit per volume. + double *chargedTrackLength = nullptr; + COPCORE_CUDA_CHECK(cudaMalloc(&chargedTrackLength, sizeof(double) * numPlaced)); + COPCORE_CUDA_CHECK(cudaMemset(chargedTrackLength, 0, sizeof(double) * numPlaced)); + double *energyDeposit = nullptr; + COPCORE_CUDA_CHECK(cudaMalloc(&energyDeposit, sizeof(double) * numPlaced)); + COPCORE_CUDA_CHECK(cudaMemset(energyDeposit, 0, sizeof(double) * numPlaced)); + + // Allocate and initialize scoring and statistics. + GlobalScoring *globalScoring = nullptr; + COPCORE_CUDA_CHECK(cudaMalloc(&globalScoring, sizeof(GlobalScoring))); + COPCORE_CUDA_CHECK(cudaMemset(globalScoring, 0, sizeof(GlobalScoring))); + + ScoringPerVolume *scoringPerVolume = nullptr; + ScoringPerVolume scoringPerVolume_devPtrs; + scoringPerVolume_devPtrs.chargedTrackLength = chargedTrackLength; + scoringPerVolume_devPtrs.energyDeposit = energyDeposit; + COPCORE_CUDA_CHECK(cudaMalloc(&scoringPerVolume, sizeof(ScoringPerVolume))); + COPCORE_CUDA_CHECK( + cudaMemcpy(scoringPerVolume, &scoringPerVolume_devPtrs, sizeof(ScoringPerVolume), cudaMemcpyHostToDevice)); + + Stats *stats_dev = nullptr; + COPCORE_CUDA_CHECK(cudaMalloc(&stats_dev, sizeof(Stats))); + Stats *stats = nullptr; + COPCORE_CUDA_CHECK(cudaMallocHost(&stats, sizeof(Stats))); + + // Allocate memory to hold a "vanilla" SlotManager to initialize for each batch. + SlotManager slotManagerInit(Capacity); + SlotManager *slotManagerInit_dev = nullptr; + COPCORE_CUDA_CHECK(cudaMalloc(&slotManagerInit_dev, sizeof(SlotManager))); + COPCORE_CUDA_CHECK(cudaMemcpy(slotManagerInit_dev, &slotManagerInit, sizeof(SlotManager), cudaMemcpyHostToDevice)); + + COPCORE_CUDA_CHECK(cudaMalloc(&BzFieldValue_dev, sizeof(BzFieldValue_host))); + COPCORE_CUDA_CHECK(cudaMemcpy(BzFieldValue_dev, &BzFieldValue_host, sizeof(BzFieldValue_host), cudaMemcpyHostToDevice)); + std::cout << " Host: passed value of BzField to device at " << BzFieldValue_dev << " value = " << BzFieldValue_host << "\n"; + + // Pass the location of the BzFieldValue_dev ! + SetBzFieldPtr<<<1,1>>>(BzFieldValue_dev); + COPCORE_CUDA_CHECK(cudaStreamSynchronize(stream)); + + vecgeom::Stopwatch timer; + timer.Start(); + + std::cout << "Entering loop to simulate " << numParticles << " particles. " << std::flush; + + std::cout << std::endl << "Simulating particles "; + const bool detailed = true; // (numParticles / batch) < 50; + if (!detailed) { + std::cout << "... " << std::flush; + } + + constexpr bool printStats= false; + + unsigned long long killed = 0; + + for (int startEvent = 1; startEvent <= numParticles; startEvent += batch) { + if (detailed) { + std::cout << startEvent << " ... " << std::flush; + } + int left = numParticles - startEvent + 1; + int chunk = std::min(left, batch); + + for (int i = 0; i < ParticleType::NumParticleTypes; i++) { + COPCORE_CUDA_CHECK(cudaMemcpyAsync(particles[i].slotManager, slotManagerInit_dev, ManagerSize, + cudaMemcpyDeviceToDevice, stream)); + } + + // Initialize primary particles. + constexpr int InitThreads = 32; + int initBlocks = (chunk + InitThreads - 1) / InitThreads; + ParticleGenerator electronGenerator(electrons.tracks, electrons.slotManager, electrons.queues.currentlyActive); + auto world_dev = vecgeom::cxx::CudaManager::Instance().world_gpu(); + InitPrimaries<<>>(electronGenerator, startEvent, chunk, energy, world_dev, + globalScoring, rotatingParticleGun); + COPCORE_CUDA_CHECK(cudaStreamSynchronize(stream)); + + stats->inFlight[ParticleType::Electron] = chunk; + stats->inFlight[ParticleType::Positron] = 0; + stats->inFlight[ParticleType::Gamma] = 0; + + constexpr int MaxBlocks = 1024; + constexpr int TransportThreads = 32; + int transportBlocks; + + int inFlight; + int loopingNo = 0; + int previousElectrons = -1, previousPositrons = -1; + + // constexpr int mod_it=25; + int iteration= 0; + + do { + Secondaries secondaries = { + .electrons = {electrons.tracks, electrons.slotManager, electrons.queues.nextActive}, + .positrons = {positrons.tracks, positrons.slotManager, positrons.queues.nextActive}, + .gammas = {gammas.tracks, gammas.slotManager, gammas.queues.nextActive}, + }; + // printf("\n\n-- Top of loop ----------------------\n"); + + // *** ELECTRONS *** + int numElectrons = stats->inFlight[ParticleType::Electron]; + if (numElectrons > 0) { + transportBlocks = (numElectrons + TransportThreads - 1) / TransportThreads; + transportBlocks = std::min(transportBlocks, MaxBlocks); + + TransportElectrons<<>>( + electrons.tracks, electrons.queues.currentlyActive, secondaries, electrons.queues.nextActive, globalScoring, + scoringPerVolume); + + COPCORE_CUDA_CHECK(cudaEventRecord(electrons.event, electrons.stream)); + COPCORE_CUDA_CHECK(cudaStreamWaitEvent(stream, electrons.event, 0)); + } + + // *** POSITRONS *** + int numPositrons = stats->inFlight[ParticleType::Positron]; + if (numPositrons > 0) { + transportBlocks = (numPositrons + TransportThreads - 1) / TransportThreads; + transportBlocks = std::min(transportBlocks, MaxBlocks); + + TransportPositrons<<>>( + positrons.tracks, positrons.queues.currentlyActive, secondaries, positrons.queues.nextActive, globalScoring, + scoringPerVolume); + + COPCORE_CUDA_CHECK(cudaEventRecord(positrons.event, positrons.stream)); + COPCORE_CUDA_CHECK(cudaStreamWaitEvent(stream, positrons.event, 0)); + } + // *** GAMMAS *** + int numGammas = stats->inFlight[ParticleType::Gamma]; + if (numGammas > 0) { + transportBlocks = (numGammas + TransportThreads - 1) / TransportThreads; + transportBlocks = std::min(transportBlocks, MaxBlocks); + + TransportGammas<<>>( + gammas.tracks, gammas.queues.currentlyActive, secondaries, gammas.queues.nextActive, globalScoring, + scoringPerVolume); + + COPCORE_CUDA_CHECK(cudaEventRecord(gammas.event, gammas.stream)); + COPCORE_CUDA_CHECK(cudaStreamWaitEvent(stream, gammas.event, 0)); + } + + // *** END OF TRANSPORT *** + + // Extra sync for debugging !!? JA 2022.09.19 + // COPCORE_CUDA_CHECK(cudaStreamSynchronize(stream)); + + // The events ensure synchronization before finishing this iteration and + // copying the Stats back to the host. + AllParticleQueues queues = {{electrons.queues, positrons.queues, gammas.queues}}; + FinishIteration<<<1, 1, 0, stream>>>(queues, stats_dev); + COPCORE_CUDA_CHECK(cudaMemcpyAsync(stats, stats_dev, sizeof(Stats), cudaMemcpyDeviceToHost, stream)); + + // Finally synchronize all kernels. + COPCORE_CUDA_CHECK(cudaStreamSynchronize(stream)); + + // Count the number of particles in flight. + inFlight = 0; + for (int i = 0; i < ParticleType::NumParticleTypes; i++) { + inFlight += stats->inFlight[i]; + // int npi= stats->inFlight[i]; inFlight += npi; std::cout << " [" << i << "]=" << npi << " "; + } + + // Swap the queues for the next iteration. + electrons.queues.SwapActive(); + positrons.queues.SwapActive(); + gammas.queues.SwapActive(); + + // Check if only charged particles are left that are looping. + numElectrons = stats->inFlight[ParticleType::Electron]; + numPositrons = stats->inFlight[ParticleType::Positron]; + numGammas = stats->inFlight[ParticleType::Gamma]; + + if (numElectrons == previousElectrons && numPositrons == previousPositrons && numGammas == 0) { + loopingNo++; + } else { + previousElectrons = numElectrons; + previousPositrons = numPositrons; + loopingNo = 0; + } + + iteration++; + const int iterModPrint = 100; // Every how many to print info + if( printStats && ( iteration % iterModPrint == 0 ) ) { + printf("Iteration = %6d event = %4d inflight= %6d ( e- : %6d g : %5d e+ : %4d) iter looping=%3d \n" , + iteration, startEvent, inFlight, + numElectrons, numGammas, numPositrons, loopingNo ); + } + + } while (inFlight > 0 && loopingNo < 200); + + // if( verbose ) std::cout << std::endl; + + if( printStats ) { + printf( "-- Tracks: %5d InFlight at the end (killed). Iterations needed = %4d \n", inFlight, iteration ); + } + + if (inFlight > 0) { + killed += inFlight; + for (int i = 0; i < ParticleType::NumParticleTypes; i++) { + ParticleType &pType = particles[i]; + int inFlightParticles = stats->inFlight[i]; + if (inFlightParticles == 0) { + continue; + } + + ClearQueue<<<1, 1, 0, stream>>>(pType.queues.currentlyActive); + } + COPCORE_CUDA_CHECK(cudaStreamSynchronize(stream)); + } + } + std::cout << "\ndone!" << std::endl; + + auto time = timer.Stop(); + std::cout << "Run time: " << time << "\n"; + + // Transfer back scoring. + COPCORE_CUDA_CHECK(cudaMemcpy(globalScoring_host, globalScoring, sizeof(GlobalScoring), cudaMemcpyDeviceToHost)); + globalScoring_host->numKilled = killed; + + // Transfer back the scoring per volume (charged track length and energy deposit). + COPCORE_CUDA_CHECK(cudaMemcpy(scoringPerVolume_host->chargedTrackLength, scoringPerVolume_devPtrs.chargedTrackLength, + sizeof(double) * numPlaced, cudaMemcpyDeviceToHost)); + COPCORE_CUDA_CHECK(cudaMemcpy(scoringPerVolume_host->energyDeposit, scoringPerVolume_devPtrs.energyDeposit, + sizeof(double) * numPlaced, cudaMemcpyDeviceToHost)); + + // Free resources. + COPCORE_CUDA_CHECK(cudaFree(MCIndex_dev)); + COPCORE_CUDA_CHECK(cudaFree(chargedTrackLength)); + COPCORE_CUDA_CHECK(cudaFree(energyDeposit)); + + COPCORE_CUDA_CHECK(cudaFree(globalScoring)); + COPCORE_CUDA_CHECK(cudaFree(scoringPerVolume)); + COPCORE_CUDA_CHECK(cudaFree(stats_dev)); + COPCORE_CUDA_CHECK(cudaFreeHost(stats)); + COPCORE_CUDA_CHECK(cudaFree(slotManagerInit_dev)); + COPCORE_CUDA_CHECK(cudaFree(BzFieldValue_dev)); + COPCORE_CUDA_CHECK(cudaStreamDestroy(stream)); + + for (int i = 0; i < ParticleType::NumParticleTypes; i++) { + COPCORE_CUDA_CHECK(cudaFree(particles[i].tracks)); + COPCORE_CUDA_CHECK(cudaFree(particles[i].slotManager)); + + COPCORE_CUDA_CHECK(cudaFree(particles[i].queues.currentlyActive)); + COPCORE_CUDA_CHECK(cudaFree(particles[i].queues.nextActive)); + + COPCORE_CUDA_CHECK(cudaStreamDestroy(particles[i].stream)); + COPCORE_CUDA_CHECK(cudaEventDestroy(particles[i].event)); + } +} diff --git a/examples/Example15/example15.cuh b/examples/Example15/example15.cuh new file mode 100644 index 00000000..2b01ee5b --- /dev/null +++ b/examples/Example15/example15.cuh @@ -0,0 +1,171 @@ +// SPDX-FileCopyrightText: 2021 CERN +// SPDX-License-Identifier: Apache-2.0 + +#ifndef EXAMPLE15_CUH +#define EXAMPLE15_CUH + +#include "example15.h" + +#include +#include +#include + +#include +#include +#include + +#include +#include + +// A data structure to represent a particle track. The particle type is implicit +// by the queue and not stored in memory. +struct Track { + using Precision = vecgeom::Precision; + RanluxppDouble rngState; + double energy; + double numIALeft[3]; + double initialRange; + double dynamicRangeFactor; + double tlimitMin; + + vecgeom::Vector3D pos; + vecgeom::Vector3D dir; + vecgeom::NavStateIndex navState; + + __host__ __device__ double Uniform() { return rngState.Rndm(); } + + __host__ __device__ void InitAsSecondary(const vecgeom::Vector3D &parentPos, + const vecgeom::NavStateIndex &parentNavState) + { + // The caller is responsible to branch a new RNG state and to set the energy. + this->numIALeft[0] = -1.0; + this->numIALeft[1] = -1.0; + this->numIALeft[2] = -1.0; + + this->initialRange = -1.0; + this->dynamicRangeFactor = -1.0; + this->tlimitMin = -1.0; + + // A secondary inherits the position of its parent; the caller is responsible + // to update the directions. + this->pos = parentPos; + this->navState = parentNavState; + } + + __host__ __device__ void print(int id = -1, bool verbose = false) const + { + using copcore::units::MeV; + using copcore::units::mm; + // static const char *particleName[3] = {"e-", "g", "e+"}; + printf( " id= %3d Pos[mm]: " + "%10.7g, %10.7g, %10.7g " + " kE[MeV]= %10.5g Dir: %9.6f %9.6f %9.6f " + " - intl= %6.3f %6.3f %6.3f r0= %10.5g\n", + id, pos[0] / mm, pos[1] / mm, pos[2] / mm , + energy / MeV, dir[0], dir[1], dir[2], + numIALeft[0] , numIALeft[1], numIALeft[2], + initialRange / mm ); + +#ifdef COPCORE_DEVICE_COMPILATION + if (verbose) { + auto currentLevel = navState.GetLevel(); + auto currentIndex = navState.GetNavIndex(currentLevel); + printf(" id= %3d current: (lv= %3d ind= %8u bnd= %1d) ", id, currentLevel, currentIndex, + (int) navState.IsOnBoundary() ); + printf("\n"); + } +#endif + } +}; + +// Defined in example15.cu +// extern __constant__ __device__ int Zero; + +#ifdef __CUDA_ARCH__ +// Define inline implementations of the RNG methods for the device. +// (nvcc ignores the __device__ attribute in definitions, so this is only to +// communicate the intent.) +inline __device__ double G4HepEmRandomEngine::flat() +{ + return ((RanluxppDouble *)fObject)->Rndm(); +} + +inline __device__ void G4HepEmRandomEngine::flatArray(const int size, double *vect) +{ + for (int i = 0; i < size; i++) { + vect[i] = ((RanluxppDouble *)fObject)->Rndm(); + } +} +#endif + +// A data structure to manage slots in the track storage. +class SlotManager { + adept::Atomic_t fNextSlot; + const int fMaxSlot; + +public: + __host__ __device__ SlotManager(int maxSlot) : fMaxSlot(maxSlot) { fNextSlot = 0; } + + __host__ __device__ int NextSlot() + { + int next = fNextSlot.fetch_add(1); + if (next >= fMaxSlot) return -1; + return next; + } +}; + +// A bundle of pointers to generate particles of an implicit type. +class ParticleGenerator { + Track *fTracks; + SlotManager *fSlotManager; + adept::MParray *fActiveQueue; + +public: + __host__ __device__ ParticleGenerator(Track *tracks, SlotManager *slotManager, adept::MParray *activeQueue) + : fTracks(tracks), fSlotManager(slotManager), fActiveQueue(activeQueue) + { + } + + __host__ __device__ Track &NextTrack() + { + int slot = fSlotManager->NextSlot(); + if (slot == -1) { + COPCORE_EXCEPTION("No slot available in ParticleGenerator::NextTrack"); + } + fActiveQueue->push_back(slot); + return fTracks[slot]; + } +}; + +// A bundle of generators for the three particle types. +struct Secondaries { + ParticleGenerator electrons; + ParticleGenerator positrons; + ParticleGenerator gammas; +}; + +// Kernels in different TUs. +__global__ void TransportElectrons(Track *electrons, const adept::MParray *active, Secondaries secondaries, + adept::MParray *activeQueue, GlobalScoring *globalScoring, + ScoringPerVolume *scoringPerVolume); +__global__ void TransportPositrons(Track *positrons, const adept::MParray *active, Secondaries secondaries, + adept::MParray *activeQueue, GlobalScoring *globalScoring, + ScoringPerVolume *scoringPerVolume); + +__global__ void TransportGammas(Track *gammas, const adept::MParray *active, Secondaries secondaries, + adept::MParray *activeQueue, GlobalScoring *globalScoring, + ScoringPerVolume *scoringPerVolume); + +// Constant data structures from G4HepEm accessed by the kernels. +// (defined in TestEm3.cu) +extern __constant__ __device__ struct G4HepEmParameters g4HepEmPars; +extern __constant__ __device__ struct G4HepEmData g4HepEmData; + +extern __constant__ __device__ int *MCIndex; + +extern float BzFieldValue_host; // = 1.0 * copcore::units::tesla; +extern /*__device__*/ float *BzFieldValue_dev; + +__global__ void SetBzFieldPtr( float* pBzFieldValue_dev ); + +#endif diff --git a/examples/Example15/example15.h b/examples/Example15/example15.h new file mode 100644 index 00000000..229b651e --- /dev/null +++ b/examples/Example15/example15.h @@ -0,0 +1,33 @@ +// SPDX-FileCopyrightText: 2021 CERN +// SPDX-License-Identifier: Apache-2.0 + +#ifndef EXAMPLE15_H +#define EXAMPLE15_H + +struct G4HepEmState; + +// Data structures for scoring. The accessors must make sure to use atomic operations if needed. +struct GlobalScoring { + double energyDeposit; + // Not int to avoid overflows for more than 100,000 events; unsigned long long + // is the only other data type available for atomicAdd(). + unsigned long long chargedSteps; + unsigned long long neutralSteps; + unsigned long long hits; + unsigned long long numGammas; + unsigned long long numElectrons; + unsigned long long numPositrons; + // Not used on the device, filled in by the host. + unsigned long long numKilled; +}; + +struct ScoringPerVolume { + double *energyDeposit; + double *chargedTrackLength; +}; + +// Interface between C++ and CUDA. +void example15(int numParticles, double energy, int batch, const int *MCCindex, ScoringPerVolume *scoringPerVolume, + GlobalScoring *globalScoring, int numVolumes, int numPlaced, G4HepEmState *state, bool rotatingParticleGun); + +#endif diff --git a/examples/Example15/gammas.cu b/examples/Example15/gammas.cu new file mode 100644 index 00000000..5561ce99 --- /dev/null +++ b/examples/Example15/gammas.cu @@ -0,0 +1,237 @@ +// SPDX-FileCopyrightText: 2022 CERN +// SPDX-License-Identifier: Apache-2.0 + +#include "example15.cuh" + +#include + +#include + +#include +#include +#include +#include +#include +#include +// Pull in implementation. +#include +#include +#include +#include + +__global__ void TransportGammas(Track *gammas, const adept::MParray *active, Secondaries secondaries, + adept::MParray *activeQueue, GlobalScoring *globalScoring, + ScoringPerVolume *scoringPerVolume) +{ +#ifdef VECGEOM_FLOAT_PRECISION + const Precision kPush = 10 * vecgeom::kTolerance; +#else + const Precision kPush = 0.; +#endif + int activeSize = active->size(); + for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < activeSize; i += blockDim.x * gridDim.x) { + const int slot = (*active)[i]; + Track ¤tTrack = gammas[slot]; + auto energy = currentTrack.energy; + auto pos = currentTrack.pos; + auto dir = currentTrack.dir; + auto navState = currentTrack.navState; + const auto volume = navState.Top(); + const int volumeID = volume->id(); + // the MCC vector is indexed by the logical volume id + const int theMCIndex = MCIndex[volume->GetLogicalVolume()->id()]; + + auto survive = [&] { + currentTrack.energy = energy; + currentTrack.pos = pos; + currentTrack.dir = dir; + currentTrack.navState = navState; + activeQueue->push_back(slot); + }; + + // Init a track with the needed data to call into G4HepEm. + G4HepEmGammaTrack gammaTrack; + G4HepEmTrack *theTrack = gammaTrack.GetTrack(); + theTrack->SetEKin(energy); + theTrack->SetMCIndex(theMCIndex); + + // Sample the `number-of-interaction-left` and put it into the track. + for (int ip = 0; ip < 3; ++ip) { + double numIALeft = currentTrack.numIALeft[ip]; + if (numIALeft <= 0) { + numIALeft = -std::log(currentTrack.Uniform()); + } + theTrack->SetNumIALeft(numIALeft, ip); + } + + // Call G4HepEm to compute the physics step limit. + G4HepEmGammaManager::HowFar(&g4HepEmData, &g4HepEmPars, &gammaTrack); + + // Get result into variables. + double geometricalStepLengthFromPhysics = theTrack->GetGStepLength(); + int winnerProcessIndex = theTrack->GetWinnerProcessIndex(); + // Leave the range and MFP inside the G4HepEmTrack. If we split kernels, we + // also need to carry them over! + + // Check if there's a volume boundary in between. + vecgeom::NavStateIndex nextState; + double geometryStepLength = + BVHNavigator::ComputeStepAndNextVolume(pos, dir, geometricalStepLengthFromPhysics, navState, nextState, kPush); + pos += geometryStepLength * dir; + atomicAdd(&globalScoring->neutralSteps, 1); + + // Set boundary state in navState so the next step and secondaries get the + // correct information (navState = nextState only if relocated + // in case of a boundary; see below) + navState.SetBoundaryState(nextState.IsOnBoundary()); + + // Propagate information from geometrical step to G4HepEm. + theTrack->SetGStepLength(geometryStepLength); + theTrack->SetOnBoundary(nextState.IsOnBoundary()); + + G4HepEmGammaManager::UpdateNumIALeft(theTrack); + + // Save the `number-of-interaction-left` in our track. + for (int ip = 0; ip < 3; ++ip) { + double numIALeft = theTrack->GetNumIALeft(ip); + currentTrack.numIALeft[ip] = numIALeft; + } + + if (nextState.IsOnBoundary()) { + // For now, just count that we hit something. + atomicAdd(&globalScoring->hits, 1); + + // Kill the particle if it left the world. + if (nextState.Top() != nullptr) { + BVHNavigator::RelocateToNextVolume(pos, dir, nextState); + + // Move to the next boundary. + navState = nextState; + survive(); + } + continue; + } else if (winnerProcessIndex < 0) { + // No discrete process, move on. + survive(); + continue; + } + + // Reset number of interaction left for the winner discrete process. + // (Will be resampled in the next iteration.) + currentTrack.numIALeft[winnerProcessIndex] = -1.0; + + // Perform the discrete interaction. + G4HepEmRandomEngine rnge(¤tTrack.rngState); + // We might need one branched RNG state, prepare while threads are synchronized. + RanluxppDouble newRNG(currentTrack.rngState.Branch()); + + switch (winnerProcessIndex) { + case 0: { + // Invoke gamma conversion to e-/e+ pairs, if the energy is above the threshold. + if (energy < 2 * copcore::units::kElectronMassC2) { + survive(); + continue; + } + + double logEnergy = std::log(energy); + double elKinEnergy, posKinEnergy; + G4HepEmGammaInteractionConversion::SampleKinEnergies(&g4HepEmData, energy, logEnergy, theMCIndex, elKinEnergy, + posKinEnergy, &rnge); + + double dirPrimary[] = {dir.x(), dir.y(), dir.z()}; + double dirSecondaryEl[3], dirSecondaryPos[3]; + G4HepEmGammaInteractionConversion::SampleDirections(dirPrimary, dirSecondaryEl, dirSecondaryPos, elKinEnergy, + posKinEnergy, &rnge); + + Track &electron = secondaries.electrons.NextTrack(); + Track &positron = secondaries.positrons.NextTrack(); + atomicAdd(&globalScoring->numElectrons, 1); + atomicAdd(&globalScoring->numPositrons, 1); + + electron.InitAsSecondary(pos, navState); + electron.rngState = newRNG; + electron.energy = elKinEnergy; + electron.dir.Set(dirSecondaryEl[0], dirSecondaryEl[1], dirSecondaryEl[2]); + + positron.InitAsSecondary(pos, navState); + // Reuse the RNG state of the dying track. + positron.rngState = currentTrack.rngState; + positron.energy = posKinEnergy; + positron.dir.Set(dirSecondaryPos[0], dirSecondaryPos[1], dirSecondaryPos[2]); + + // The current track is killed by not enqueuing into the next activeQueue. + break; + } + case 1: { + // Invoke Compton scattering of gamma. + constexpr double LowEnergyThreshold = 100 * copcore::units::eV; + if (energy < LowEnergyThreshold) { + survive(); + continue; + } + const double origDirPrimary[] = {dir.x(), dir.y(), dir.z()}; + double dirPrimary[3]; + const double newEnergyGamma = + G4HepEmGammaInteractionCompton::SamplePhotonEnergyAndDirection(energy, dirPrimary, origDirPrimary, &rnge); + vecgeom::Vector3D newDirGamma(dirPrimary[0], dirPrimary[1], dirPrimary[2]); + + const double energyEl = energy - newEnergyGamma; + if (energyEl > LowEnergyThreshold) { + // Create a secondary electron and sample/compute directions. + Track &electron = secondaries.electrons.NextTrack(); + atomicAdd(&globalScoring->numElectrons, 1); + + electron.InitAsSecondary(pos, navState); + electron.rngState = newRNG; + electron.energy = energyEl; + electron.dir = energy * dir - newEnergyGamma * newDirGamma; + electron.dir.Normalize(); + } else { + atomicAdd(&globalScoring->energyDeposit, energyEl); + atomicAdd(&scoringPerVolume->energyDeposit[volumeID], energyEl); + } + + // Check the new gamma energy and deposit if below threshold. + if (newEnergyGamma > LowEnergyThreshold) { + energy = newEnergyGamma; + dir = newDirGamma; + survive(); + } else { + atomicAdd(&globalScoring->energyDeposit, newEnergyGamma); + atomicAdd(&scoringPerVolume->energyDeposit[volumeID], newEnergyGamma); + // The current track is killed by not enqueuing into the next activeQueue. + } + break; + } + case 2: { + // Invoke photoelectric process. + const double theLowEnergyThreshold = 1 * copcore::units::eV; + + const double bindingEnergy = G4HepEmGammaInteractionPhotoelectric::SelectElementBindingEnergy( + &g4HepEmData, theMCIndex, gammaTrack.GetPEmxSec(), energy, &rnge); + + double edep = bindingEnergy; + const double photoElecE = energy - edep; + if (photoElecE > theLowEnergyThreshold) { + // Create a secondary electron and sample directions. + Track &electron = secondaries.electrons.NextTrack(); + atomicAdd(&globalScoring->numElectrons, 1); + + double dirGamma[] = {dir.x(), dir.y(), dir.z()}; + double dirPhotoElec[3]; + G4HepEmGammaInteractionPhotoelectric::SamplePhotoElectronDirection(photoElecE, dirGamma, dirPhotoElec, &rnge); + + electron.InitAsSecondary(pos, navState); + electron.rngState = newRNG; + electron.energy = photoElecE; + electron.dir.Set(dirPhotoElec[0], dirPhotoElec[1], dirPhotoElec[2]); + } else { + edep = energy; + } + atomicAdd(&globalScoring->energyDeposit, edep); + // The current track is killed by not enqueuing into the next activeQueue. + break; + } + } + } +} diff --git a/examples/TestEm3/README.md b/examples/TestEm3/README.md index abdd6cb3..31fd25cd 100644 --- a/examples/TestEm3/README.md +++ b/examples/TestEm3/README.md @@ -20,7 +20,7 @@ Additionally, the kernels score energy deposit and the charged track length per ### Kernels This example uses one stream per particle type to launch kernels asynchronously. -They are synchronized via a forth stream using CUDA events. +They are synchronized via a fourth stream using CUDA events. #### `TransportElectrons` diff --git a/examples/TestEm3MT/README.md b/examples/TestEm3MT/README.md index 36a17caf..9cf53b29 100644 --- a/examples/TestEm3MT/README.md +++ b/examples/TestEm3MT/README.md @@ -22,7 +22,7 @@ This version uses multiple threads on the host to submit work to the GPU. ### Kernels This example uses one stream per particle type to launch kernels asynchronously. -They are synchronized via a forth stream using CUDA events. +They are synchronized via a fourth stream using CUDA events. #### `TransportElectrons` diff --git a/magneticfield/inc/CompareResponses.h b/magneticfield/inc/CompareResponses.h new file mode 100644 index 00000000..3a03e4c5 --- /dev/null +++ b/magneticfield/inc/CompareResponses.h @@ -0,0 +1,71 @@ +// SPDX-FileCopyrightText: 2020-3 CERN +// SPDX-License-Identifier: Apache-2.0 + +// Author: J. Apostolakis, 22 June 2022 + +#ifndef CompareResponses_hh +#define CompareResponses_hh 1 + +#include + +static __device__ __host__ +bool CompareResponseVector3D( + int id, + vecgeom::Vector3D const & originalVec, + vecgeom::Vector3D const & baselineVec, + vecgeom::Vector3D const & resultVec, // Output of new method + const char * vecName, + Precision thresholdRel // fraction difference allowed + ) +// Returns 'true' if values are 'bad'... +{ + bool bad = false; // Good .. + Precision magOrig= originalVec.Mag(); + vecgeom::Vector3D moveBase = baselineVec-originalVec; + vecgeom::Vector3D moveRes = resultVec-originalVec; + Precision magMoveBase = moveBase.Mag(); + Precision magDiffRes = moveRes.Mag(); + + if ( std::fabs( magDiffRes / magMoveBase) - 1.0 > thresholdRel + || + ( resultVec - baselineVec ).Mag() > thresholdRel * magMoveBase + ){ + // printf("Difference seen in vector %s : ", vecName ); + printf("\n id %3d - Diff in %s: " + " new-base= %14.9g %14.9g %14.9g (mag= %14.9g) " + " mv_Res/mv_Base-1 = %7.3g | mv/base: 3v= %14.9g %14.9g %14.9g (mag= %9.4g)" + " | mv/new: 3v= %14.9g %14.9g %14.9g (mag = %14.9g)" + " || origVec= %14.9f %14.9f %14.9f (mag=%14.9f) | base= %14.9f %14.9f %14.9f (mag=%9.4g) \n", + id, vecName, + resultVec[0]-baselineVec[0], resultVec[1]-baselineVec[1], resultVec[2]-baselineVec[2], (resultVec-baselineVec).Mag(), + (moveRes.Mag() / moveBase.Mag() - 1.0), + moveBase[0], moveBase[1], moveBase[2], moveBase.Mag(), +// printf(" new-original: mag= %20.16g , new_vec= %14.9f , %14.9f , %14.9f \n", + moveRes[0], moveRes[1], moveRes[2], moveRes.Mag(), + originalVec[0], originalVec[1], originalVec[2], originalVec.Mag(), + baselineVec[0], baselineVec[1], baselineVec[2], baselineVec.Mag() + ); + bad= true; + } + return bad; +}; + +static __device__ __host__ +void ReportSameMoveVector3D( int id, + vecgeom::Vector3D const & originalVec, + vecgeom::Vector3D const & resultVec, + const char * vecName ) +{ + printf( " id %3d - Same %s: " + " mv/base: 3v= %14.9f %14.9f %14.9f (mag= %9.4g)" + " || origVec= %14.9f %14.9f %14.9f (mag=%14.9f) | result= %14.9f %14.9f %14.9f (mag=%9.4g) \n", + id, vecName, + resultVec[0]-originalVec[0], resultVec[1]-originalVec[1], resultVec[2]-originalVec[2], + (resultVec-originalVec).Mag(), + originalVec[0], originalVec[1], originalVec[2], originalVec.Mag(), + resultVec[0], resultVec[1], resultVec[2], resultVec.Mag() + ); + +} + +#endif /* CompareResponses_hh */ diff --git a/magneticfield/inc/ConstBzFieldStepper.h b/magneticfield/inc/ConstBzFieldStepper.h index 75450e0b..e3e6c128 100644 --- a/magneticfield/inc/ConstBzFieldStepper.h +++ b/magneticfield/inc/ConstBzFieldStepper.h @@ -15,6 +15,8 @@ #include "VecGeom/base/Global.h" +#include "fieldConstants.h" + // namespace adept { // inline namespace ADEPT_IMPL_NAMESPACE { @@ -30,14 +32,14 @@ class ConstBzFieldStepper { Precision fBz; public: - VECCORE_ATT_HOST_DEVICE - ConstBzFieldStepper(Precision Bz = 0) : fBz(Bz) {} + __host__ __device__ + ConstBzFieldStepper(float Bz = 0.) : fBz(Bz) {} void SetBz(Precision Bz) { fBz = Bz; } Precision GetBz() const { return fBz; } - static constexpr Precision kB2C = - -0.299792458 * (copcore::units::GeV / (copcore::units::tesla * copcore::units::meter)); + // static constexpr Precision kB2C = + // -0.299792458 * (copcore::units::GeV / (copcore::units::tesla * copcore::units::meter)); /* template @@ -45,7 +47,7 @@ class ConstBzFieldStepper { double const charge, double const momentum ) const { if (charge == 0) return RT(0.); - return abs( kB2C * fBz * dir.FastInverseScaledXYLength( momentum ) ); + return abs( fieldConstants::kB2C * fBz * dir.FastInverseScaledXYLength( momentum ) ); } */ @@ -55,7 +57,7 @@ class ConstBzFieldStepper { * output: new position, new direction of particle */ template - inline __attribute__((always_inline)) VECCORE_ATT_HOST_DEVICE void DoStep( + inline __attribute__((always_inline)) __host__ __device__ void DoStep( BaseType const & /*posx*/, BaseType const & /*posy*/, BaseType const & /*posz*/, BaseType const & /*dirx*/, BaseType const & /*diry*/, BaseType const & /*dirz*/, BaseIType const & /*charge*/, BaseType const & /*momentum*/, BaseType const & /*step*/, BaseType & /*newsposx*/, BaseType & /*newposy*/, BaseType & /*newposz*/, @@ -68,7 +70,7 @@ class ConstBzFieldStepper { * output: new position, new direction of particle */ template - VECCORE_ATT_HOST_DEVICE void DoStep(Vector3D const &pos, Vector3D const &dir, BaseIType const &charge, + __host__ __device__ void DoStep(Vector3D const &pos, Vector3D const &dir, BaseIType const &charge, BaseType const &momentum, BaseType const &step, Vector3D &newpos, Vector3D &newdir) const { @@ -94,11 +96,11 @@ inline __attribute__((always_inline)) void ConstBzFieldStepper::DoStep( BaseDType dt = sqrt((dx0 * dx0) + (dy0 * dy0)) + kSmall; BaseDType invnorm = 1. / dt; // radius has sign and determines the sense of rotation - BaseDType R = momentum * dt / ((BaseDType(kB2C) * BaseDType(charge)) * (fBz)); + BaseDType R = momentum * dt / ((BaseDType(fieldConstants::kB2C) * BaseDType(charge)) * (fBz)); BaseDType cosa = dx0 * invnorm; BaseDType sina = dy0 * invnorm; - BaseDType phi = step * BaseDType(charge) * fBz * BaseDType(kB2C) / momentum; + BaseDType phi = step * BaseDType(charge) * fBz * BaseDType(fieldConstants::kB2C) / momentum; BaseDType cosphi; BaseDType sinphi; diff --git a/magneticfield/inc/ConstFieldHelixStepper.h b/magneticfield/inc/ConstFieldHelixStepper.h index 3dd04b67..b5f3bd1e 100644 --- a/magneticfield/inc/ConstFieldHelixStepper.h +++ b/magneticfield/inc/ConstFieldHelixStepper.h @@ -12,9 +12,11 @@ #ifndef CONSTFIELDHELIXSTEPPER_H_ #define CONSTFIELDHELIXSTEPPER_H_ -#include "VecGeom/base/Global.h" +#include // Needed for Vector3D +#include "fieldConstants.h" + // namespace adept { // inline namespace ADEPT_IMPL_NAMESPACE { /** @@ -28,31 +30,29 @@ class ConstFieldHelixStepper { using Vector3D = vecgeom::Vector3D; public: - // VECCORE_ATT_HOST_DEVICE + // __host__ __device__ // ConstFieldHelixStepper(); // For default initialisation only - VECCORE_ATT_HOST_DEVICE + __host__ __device__ ConstFieldHelixStepper(Precision Bx, Precision By, Precision Bz); - VECCORE_ATT_HOST_DEVICE + __host__ __device__ ConstFieldHelixStepper(Precision Bfield[3]); - VECCORE_ATT_HOST_DEVICE + __host__ __device__ ConstFieldHelixStepper(Vector3D const &Bfield); - void VECCORE_ATT_HOST_DEVICE SetB(Precision Bx, Precision By, Precision Bz) + void __host__ __device__ + SetB(Precision Bx, Precision By, Precision Bz) { // fB.Set(Bx, By, Bz); Vector3D Bfield(Bx, By, Bz); CalculateDerived(Bfield); } - VECCORE_ATT_HOST_DEVICE + __host__ __device__ Vector3D const GetFieldVec() const { return fBmag * fUnit; } - static constexpr Precision kB2C = - -0.299792458 * (copcore::units::GeV / (copcore::units::tesla * copcore::units::meter)); - /* template RT GetCurvature(Vector3D const & dir, @@ -69,7 +69,7 @@ class ConstFieldHelixStepper { * output: new position, new direction of particle */ template - VECCORE_ATT_HOST_DEVICE void DoStep(Real_t const &posx, Real_t const &posy, Real_t const &posz, Real_t const &dirx, + __host__ __device__ void DoStep(Real_t const &posx, Real_t const &posy, Real_t const &posz, Real_t const &dirx, Real_t const &diry, Real_t const &dirz, Int_t const &charge, Real_t const &momentum, Real_t const &step, Real_t &newposx, Real_t &newposy, Real_t &newposz, Real_t &newdirx, Real_t &newdiry, Real_t &newdirz) const; @@ -80,7 +80,7 @@ class ConstFieldHelixStepper { * output: new position, new direction of particle */ template - inline VECCORE_ATT_HOST_DEVICE void DoStep(Vector3D const &position, Vector3D const &direction, + inline __host__ __device__ void DoStep(Vector3D const &position, Vector3D const &direction, Int_t const &charge, Real_t const &momentum, Real_t const &step, Vector3D &endPosition, Vector3D &endDirection) const; @@ -91,10 +91,10 @@ class ConstFieldHelixStepper { vecgeom::Vector3D &endPosition, vecgeom::Vector3D &endDirection) const; protected: - inline VECCORE_ATT_HOST_DEVICE void CalculateDerived(Vector3D Bvec); + inline __host__ __device__ void CalculateDerived(Vector3D Bvec); template - inline VECCORE_ATT_HOST_DEVICE bool CheckModulus(Real_t &newdirX_v, Real_t &newdirY_v, Real_t &newdirZ_v) const; + inline __host__ __device__ bool CheckModulus(Real_t &newdirX_v, Real_t &newdirY_v, Real_t &newdirZ_v) const; private: // Values below used for speed, code simplicity @@ -103,9 +103,8 @@ class ConstFieldHelixStepper { }; // end class declaration -inline // __host__ __device__ - void - ConstFieldHelixStepper::CalculateDerived(Vector3D Bvec) +inline __host__ __device__ +void ConstFieldHelixStepper::CalculateDerived(Vector3D Bvec) { fBmag = Bvec.Mag(); Precision bMagInv = (1 / fBmag); @@ -113,17 +112,20 @@ inline // __host__ __device__ if (fBmag > 0) fUnit = bMagInv * Bvec; } -inline ConstFieldHelixStepper::ConstFieldHelixStepper(Precision Bx, Precision By, Precision Bz) // : fB(Bx, By, gBz) +inline __host__ __device__ +ConstFieldHelixStepper::ConstFieldHelixStepper(Precision Bx, Precision By, Precision Bz) // : fB(Bx, By, gBz) { CalculateDerived({Bx, By, Bz}); } -inline ConstFieldHelixStepper::ConstFieldHelixStepper(Precision B[3]) // : fB(B[0], B[1], B[2]) +inline __host__ __device__ +ConstFieldHelixStepper::ConstFieldHelixStepper(Precision B[3]) // : fB(B[0], B[1], B[2]) { CalculateDerived({B[0], B[1], B[2]}); } -inline VECCORE_ATT_HOST_DEVICE ConstFieldHelixStepper::ConstFieldHelixStepper( + +inline __host__ __device__ ConstFieldHelixStepper::ConstFieldHelixStepper( Vector3D const &Bfield) // : fB(Bfield) { CalculateDerived(Bfield); @@ -160,7 +162,7 @@ inline void ConstFieldHelixStepper::DoStep(Real_t const &x0, Real_t const &y0, R } template -inline VECCORE_ATT_HOST_DEVICE void ConstFieldHelixStepper::DoStep(vecgeom::Vector3D const &startPosition, +inline __host__ __device__ void ConstFieldHelixStepper::DoStep(vecgeom::Vector3D const &startPosition, vecgeom::Vector3D const &startDirection, Int_t const &charge, Real_t const &momentum, Real_t const &step, @@ -178,11 +180,11 @@ inline VECCORE_ATT_HOST_DEVICE void ConstFieldHelixStepper::DoStep(vecgeom::Vect Vector3D dir1Field(fUnit); Real_t UVdotUB = startDirection.Dot(dir1Field); // Limit cases 0.0 and 1.0 - Real_t dt2 = max(startDirection.Mag2() - UVdotUB * UVdotUB, Real_t(0.0)); + Real_t dt2 = vecCore::Max(startDirection.Mag2() - UVdotUB * UVdotUB, Real_t(0.0)); Real_t sinVB = sqrt(dt2) + kSmall; // radius has sign and determines the sense of rotation - Real_t R = momentum * sinVB / (kB2C * Real_t(charge) * fBmag); + Real_t R = momentum * sinVB / (fieldConstants::kB2C * Real_t(charge) * fBmag); Vector3D restVelX = startDirection - UVdotUB * dir1Field; @@ -205,14 +207,17 @@ inline VECCORE_ATT_HOST_DEVICE void ConstFieldHelixStepper::DoStep(vecgeom::Vect assert(fabs(dirVelX.Dot(dirCrossVB)) < 1.e-6); assert(fabs(dirCrossVB.Dot(dir1Field)) < 1.e-6); - Real_t phi = -step * Real_t(charge) * fBmag * kB2C / momentum; + Real_t phi = -step * Real_t(charge) * fBmag * fieldConstants::kB2C / momentum; // printf("CVFHS> phi= %g \n", vecCore::Get(phi,0) ); // phi (scalar) or phi[0] (vector) - Real_t cosphi; // = cos(phi); - Real_t sinphi; // = sin(phi); - sincos(phi, &sinphi, &cosphi); + double cosphiD, sinphiD; + sincos(phi, &sinphiD, &cosphiD); + + Real_t cosphi = cosphiD; // = cos(phi); + Real_t sinphi = sinphiD; // = sin(phi); + endPosition = startPosition + R * (cosphi - 1) * dirCrossVB - R * sinphi * dirVelX + step * UVdotUB * dir1Field; // 'Drift' along field direction diff --git a/magneticfield/inc/DormandPrinceRK45.h b/magneticfield/inc/DormandPrinceRK45.h new file mode 100644 index 00000000..14690171 --- /dev/null +++ b/magneticfield/inc/DormandPrinceRK45.h @@ -0,0 +1,238 @@ +// SPDX-FileCopyrightText: 2021 CERN +// SPDX-License-Identifier: Apache-2.0 +// +// Author: J. Apostolakis, 15 Nov 2021 +// +// Implementation of the Dormand Price 5(4) 7-stage Runge-Kutta integrator for AdePT. +// +// Notes: +// - Current version is restricted to Magnetic fields (see EvaluateDerivatives.) +// - It provides the next value of dy/ds in 'next_dydx' +// - It uses a large number of registers and/or stack locations - 7 derivatives + In + Out + Err + +template +class DormandPrinceRK45 // : public VScalarIntegrationStepper +{ + // using ThreeVector = vecgeom::Vector3D; + +public: + static constexpr unsigned int kMethodOrder = 4; + // inline DormandPrinceRK45(Equation_t *EqRhs, bool verbose = false); + // DormandPrinceRK45(const DormandPrinceRK45 &) = delete; + // ~DormandPrinceRK45() {} + + static __host__ __device__ + void EvaluateDerivatives(const T_Field & field, const Real_t y[], int charge, Real_t dydx[]) ; + + static __host__ __device__ + void StepWithErrorEstimate( const T_Field & field, + const Real_t * yIn, + const Real_t * dydx, + int charge, + Real_t Step, + Real_t * yOut, // Output: y values at end, + Real_t * yerr, // estimated errors, + Real_t * next_dydx); // next value of dydx +}; + +template +__host__ __device__ void +DormandPrinceRK45:: + EvaluateDerivatives( const T_Field& field, const Real_t yIn[], int charge, Real_t dy_ds[] ) +{ +/* #ifdef VERBOSE_RHS + using geant::units::tesla; + std::cout << "DormandPrinceRK45::EvaluateDerivatives called with q= " << charge + << " at Position = " << yIn[0] << " y= " << yIn[1] << " z= " << yIn[2] + << " with Momentum = " << yIn[3] << " y= " << yIn[4] << " z= " << yIn[5] << " "; +#endif */ + + // Vector3D Bfield; + // Equation_t::EvaluateDerivativesReturnB( field, yIn, charge, dy_ds, Bfield ); + + Equation_t::EvaluateDerivatives( /* const T_Field& */ field, yIn, charge, dy_ds ); + + /********* + using copcore::units::tesla; + using std::setw; + constexpr int prec= 5; + constexpr int nf= prec+5; + int old_prec = std::cout.precision(prec); + std::cout << " DoPri5: Evaluate Derivatives - using B-field, Bx= " << Bfield.x() / tesla << " By= " << Bfield.y() / tesla << " Bz= " << Bfield.z() / tesla << " "; + std::cout << " gives Derivs dy_ds= : " + << " x = " << setw(nf) << dy_ds[0] << " y = " << setw(nf) << dy_ds[1] << " z = " << setw(nf) << dy_ds[2] + << " px= " << setw(nf) << dy_ds[3] << " py= " << setw(nf) << dy_ds[4] << " pz= " << setw(nf) << dy_ds[5] + << std::endl; + std::cout.precision(old_prec); + ********/ +} + +template +inline +__host__ __device__ void +DormandPrinceRK45:: + StepWithErrorEstimate( const T_Field & field, + const Real_t * yIn, + const Real_t * dydx, + int charge, + Real_t Step, + Real_t * yOut, + Real_t * yErr, + Real_t * next_dydx ) +// yIn and yOut MUST NOT be aliases for same array +{ + assert( yIn != yOut ); + + static constexpr Real_t + b21 = 0.2 , + + b31 = 3.0/40.0, b32 = 9.0/40.0 , + + b41 = 44.0/45.0, b42 = -56.0/15.0, b43 = 32.0/9.0, + + b51 = 19372.0/6561.0, b52 = -25360.0/2187.0, b53 = 64448.0/6561.0, + b54 = -212.0/729.0 , + + b61 = 9017.0/3168.0 , b62 = -355.0/33.0, b63 = 46732.0/5247.0, + b64 = 49.0/176.0 , b65 = -5103.0/18656.0 , + + b71 = 35.0/384.0, b72 = 0., b73 = 500.0/1113.0, b74 = 125.0/192.0, + b75 = -2187.0/6784.0, b76 = 11.0/84.0; + + static constexpr Real_t + dc1 = -( b71 - 5179.0/57600.0), + dc2 = -( b72 - .0), + dc3 = -( b73 - 7571.0/16695.0), + dc4 = -( b74 - 393.0/640.0), + dc5 = -( b75 + 92097.0/339200.0), + dc6 = -( b76 - 187.0/2100.0), + dc7 = -( - 1.0/40.0 ); + + // Initialise time to t0, needed when it is not updated by the integration. + // [ Note: Only for time dependent fields (usually electric) + // is it neccessary to integrate the time.] + // yOut[7] = yTemp[7] = yIn[7]; + + // EvaluateDerivatives( field, yIn, charge, dydx) ; // 1st Step + + Real_t ak2[Nvar]; + { + Real_t yTemp2[Nvar]; + for (unsigned int i = 0; i < Nvar; i++) { + yTemp2[i] = yIn[i] + b21 * Step * dydx[i]; + } + EvaluateDerivatives( field, yTemp2, charge, ak2); // 2nd Step + } + + Real_t ak3[Nvar]; + { + Real_t yTemp3[Nvar]; + for (unsigned int i = 0; i < Nvar; i++) { + yTemp3[i] = yIn[i] + Step * (b31 * dydx[i] + b32 * ak2[i]); + } + EvaluateDerivatives( field, yTemp3, charge, ak3); // 3rd Step + } + + Real_t ak4[Nvar]; + { + Real_t yTemp4[Nvar]; + for (unsigned int i = 0; i < Nvar; i++) { + yTemp4[i] = yIn[i] + Step * (b41 * dydx[i] + b42 * ak2[i] + b43 * ak3[i]); + } + EvaluateDerivatives( field, yTemp4, charge, ak4); // 4th Step + } + + Real_t ak5[Nvar]; + { + Real_t yTemp5[Nvar]; + for (unsigned int i = 0; i < Nvar; i++) { + yTemp5[i] = yIn[i] + Step * (b51 * dydx[i] + b52 * ak2[i] + b53 * ak3[i] + b54 * ak4[i]); + } + EvaluateDerivatives( field, yTemp5, charge, ak5); // 5th Step + } + + Real_t ak6[Nvar]; + { + Real_t yTemp6[Nvar]; + for (unsigned int i = 0; i < Nvar; i++) { + yTemp6[i] = + yIn[i] + Step * (b61 * dydx[i] + b62 * ak2[i] + b63 * ak3[i] + b64 * ak4[i] + b65 * ak5[i]); + } + EvaluateDerivatives( field, yTemp6, charge, ak6); // 6th Step + } + + // Real_t ak7[Nvar]; // -> Replaced by next_dydx + for(unsigned int i=0;i< Nvar;i++) + { + yOut[i] = yIn[i] + Step*(b71*dydx[i] + b72*ak2[i] + b73*ak3[i] + + b74*ak4[i] + b75*ak5[i] + b76*ak6[i] ); + } + EvaluateDerivatives( field, yOut, charge, next_dydx); //7th and Final stage + + for (unsigned int i = 0; i < Nvar; i++) { + // Estimate error as difference between 4th and 5th order methods + // + yErr[i] = Step * ( dc1 * dydx[i] + dc2 * ak2[i] + dc3 * ak3[i] + dc4 * ak4[i] + + dc5 * ak5[i] + dc6 * ak6[i] + dc7 * next_dydx[i] ); + // std::cout<< "----In Stepper, yerrr is: "< +inline Real_t DormandPrinceRK45::DistChord() const +{ + // Coefficients were taken from Some Practical Runge-Kutta Formulas by Lawrence F. Shampine, page 149, c* + static constexpr Real_t hf1 = 6025192743.0 / 30085553152.0, + hf2 = 0.0, + hf3 = 51252292925.0 / 65400821598.0, + hf4 = - 2691868925.0 / 45128329728.0, + hf5 = 187940372067.0 / 1594534317056.0, + hf6 = - 1776094331.0 / 19743644256.0, + hf7 = 11237099.0 / 235043384.0; + + Real_t midVector[3]; + + for(int i = 0; i < 3; ++i) { + midVector[i] = fLastInitialVector[i] + 0.5 * fLastStepLength * + (hf1 * fInitialDyDx[i] + hf2 * ak2[i] + hf3 * ak3[i] + + hf4 * ak4[i] + hf5 * ak5[i] + hf6 * ak6[i] + hf7 * next_dydx[i]); + } + Real_t distChord; + ThreeVector initialPoint, finalPoint, midPoint; + + initialPoint = ThreeVector(fLastInitialVector[0], fLastInitialVector[1], fLastInitialVector[2]); + finalPoint = ThreeVector(fLastFinalVector[0], fLastFinalVector[1], fLastFinalVector[2]); + midPoint = ThreeVector(midVector[0], midVector[1], midVector[2]); + + // Use stored values of Initial and Endpoint + new Midpoint to evaluate + // distance of Chord + distChord = GULineSection::Distline(midPoint, initialPoint, finalPoint); + + return distChord; +} +#endif +**********************************************************************/ + +// template +// inline void DormandPrinceRK45::PrintField(const char *label, const Real_t y[Nvar], +// const vecgeom::Vector3D &Bfield) const + + +// template +// inline void DormandPrinceRK45::PrintDyDx(const char *label, const Real_t dydx[Nvar], +// const Real_t y[Nvar]) const diff --git a/magneticfield/inc/ErrorEstimatorRK.h b/magneticfield/inc/ErrorEstimatorRK.h new file mode 100644 index 00000000..0c7d26f4 --- /dev/null +++ b/magneticfield/inc/ErrorEstimatorRK.h @@ -0,0 +1,82 @@ +// SPDX-FileCopyrightText: 2021 CERN +// SPDX-License-Identifier: Apache-2.0 +// +// ErrorEstimatorRK: Simple class to Estimate the overall relative error +// from the RK method's error per (integration) variable +// +// Created on: November 15, 2021 +// Author: J. Apostolakis + +#ifndef ErrorEstimatorRK_h +#define ErrorEstimatorRK_h + +class ErrorEstimatorRK +{ + public: + __host__ __device__ + ErrorEstimatorRK( float eps_rel_max, + float minimumStep, + int noComponents = 6 + ) : + fEpsRelMax (eps_rel_max ), + fInvEpsilonRelSq( 1.0 / (eps_rel_max * eps_rel_max) ), + fMinimumStep( minimumStep ) // , fNoComponents( 6 ) + + {} + + template + __host__ __device__ + Real_t EstimateSquareError( const Real_t yError[], + const Real_t & hStep, + // const Real_t yValue[fNoComponents], + const Real_t & magMomentumSq // Initial momentum square (used for rel. error) + ) const; + // Returns the Maximum Square Error: the maximum between + // - position relative error square (magnitude^2): (momentum error vec)^2 / (initial momentum)^2 + // - momentum relative error square (magnitude^2): (position error vec)^2 / (step length)^2 + // + // Last argument enables the use of initial momentum square in calculating the relative error + + __host__ __device__ + float GetMaxRelativeError() const { return fEpsRelMax; } + + public: + static constexpr const double tinyValue = 1.0e-30; // Just to ensure there is no division by zero + + private: + const float fEpsRelMax; + const float fInvEpsilonRelSq; // = 1.0 / (eps_rel_max * eps_rel_max); + const float fMinimumStep; + // constexpr int fNoComponents = 6; + +}; + +template +__host__ __device__ +Real_t ErrorEstimatorRK::EstimateSquareError( + const Real_t yEstError[], // [fNoComponents] + const Real_t & hStep, + const Real_t & magInitMomentumSq // (Initial) momentum square (used for rel. error) + ) const +{ + Real_t invMagMomentumSq = 1.0 / (magInitMomentumSq + tinyValue); + Real_t epsPosition; + Real_t errpos_sq, errmom_sq; + + epsPosition = fEpsRelMax * vecCore::math::Max(hStep, Real_t(fMinimumStep)); + // Note: it uses the remaining step 'hStep' + // Could change it to use full step size ==> move it outside loop !! 2017.11.10 JA + + Real_t invEpsPositionSq = 1.0 / (epsPosition * epsPosition); + + // Evaluate accuracy + errpos_sq = yEstError[0] * yEstError[0] + yEstError[1] * yEstError[1] + yEstError[2] * yEstError[2]; + errpos_sq *= invEpsPositionSq; // Scale relative to required tolerance + // Accuracy for momentum + + Real_t sumerr_sq = yEstError[3] * yEstError[3] + yEstError[4] * yEstError[4] + yEstError[5] * yEstError[5]; + errmom_sq = fInvEpsilonRelSq * invMagMomentumSq * sumerr_sq ; + + return vecCore::math::Max(errpos_sq, errmom_sq); // Maximum Square Error +} +#endif diff --git a/magneticfield/inc/MagneticFieldEquation.h b/magneticfield/inc/MagneticFieldEquation.h new file mode 100644 index 00000000..0e8b8e4b --- /dev/null +++ b/magneticfield/inc/MagneticFieldEquation.h @@ -0,0 +1,312 @@ +// SPDX-FileCopyrightText: 2021 CERN +// SPDX-License-Identifier: Apache-2.0 +// +// Author: J. Apostolakis, 12-17 Nov 2021 +// +// Equation of motion for pure magnetic field for +// use in solving ODEs of motion using Runge-Kutta methods. +// + +#ifndef MagneticFieldEquation_H_ +#define MagneticFieldEquation_H_ + +#include // For cout only + +// #include +#include + +#include + +template +class MagneticFieldEquation +{ +public: + + // Key methods + // ----------- + template + static __host__ __device__ void + EvaluateDerivativesGivenB(const Real_t y[], const Real_t Bfield[3], int charge, Real_t dy_ds[]) ; + // Evaluate + // dy_ds = d/ds ( position(x, y, z), momentum (px, py, pz) ) + // given charge and + // y = ( x, y, z, p_x, p_y, p_z ) + // Bfield = ( B_x, B_y, B_z ) + + template + static __host__ __device__ void + EvaluateDerivatives(const Real_t y[], + vecgeom::Vector3D const & Bvec, + int charge, + Real_t dy_ds[] ) + { EvaluateRhsGivenB(y, charge, Bvec, dy_ds); } + // Same with Vector3D arguments + + template + static __host__ __device__ void + EvaluateDerivatives( MagneticField_t const & magField, + const Real_t y[], + int charge, + Real_t dy_ds[]); + // Same, with MagneticField_t object, used to obtain B-field. + + template + static __host__ __device__ void + EvaluateDerivativesReturnB( MagneticField_t const & magField, + const Real_t y[], + int charge, + Real_t dy_ds[], + vecgeom::Vector3D & BfieldVal + ); + // Same, with MagneticField_t object, used to obtain B-field. + + template + inline __host__ __device__ void + EvaluateDerivativesReturnB( MagneticField_t const & magField, + vecgeom::Vector3D const & position, + vecgeom::Vector3D const & momentum, + int charge, + Real_t dy_ds[], + vecgeom::Vector3D& BfieldVal + ); + + // Implementation method + // --------------------- + template + static + __host__ __device__ void + EvaluateRhsGivenB(const Real_t y[/*Nvar*/], + const Int_t & charge, + const vecgeom::Vector3D & Bvalue, + Real_t dy_ds[/*Nvar*/]) + { + vecgeom::Vector3D momentum = { y[3], y[4], y[5] }; + Force( momentum, charge, Bvalue, + dy_ds[0], dy_ds[1], dy_ds[2], dy_ds[3], dy_ds[4], dy_ds[5] ); + // Force( y[4], y[5], B[0], B[1], B[2], charge, + // dy_ds[0], dy_ds[1], dy_ds[2], dy_ds[3], dy_ds[4], dy_ds[5] ); + } + + static constexpr double gCof = copcore::units::kCLight; // / fieldUnits::meter ; + static constexpr int Nvar = 6; + + private: + + template + static __host__ __device__ void + Force( // Vector3D const & position, + vecgeom::Vector3D const & momentum, + Int_t const & charge, + vecgeom::Vector3D const & Bfield, + // Real_t const & momentumMag, // Magnitude of initial momentum + // Vector3D dy_ds + Real_t & dx_ds, Real_t & dy_ds, Real_t & dz_ds, + Real_t & dpx_ds, Real_t & dpy_ds, Real_t & dpz_ds ); + /****** + template + static + __host__ __device__ void + Force(Real_t const & momx, Real_t const & momy, Real_t const & momz, + Real_t const & Bx, Real_t const & By, Real_t const & Bz, + Int_t const & charge, + // Real_t const & momentumMag, // Magnitude of initial momentum + Real_t & dx_ds, Real_t & dy_ds, Real_t & dz_ds, + Real_t & dpx_ds, Real_t & dpy_ds, Real_t & dpz_ds + ) ; + *********/ + // Implements the Lorentz force for RK integration of d/ds ( x, y, z, px, py, pz ) + +}; +// ------------------------------------------------------------------------------------ + +template +template +inline +__host__ __device__ void +MagneticFieldEquation:: +Force( vecgeom::Vector3D const & momentum, + Int_t const & charge, + vecgeom::Vector3D const & Bfield, + // Real_t const & momentumMag, // Magnitude of initial momentum + // Vector3D dy_ds + // Real_t dydsArr[/*Nvar*/] + Real_t & dx_ds, Real_t & dy_ds, Real_t & dz_ds, + Real_t & dpx_ds, Real_t & dpy_ds, Real_t & dpz_ds + ) +{ + Real_t inv_momentum_mag = Real_t(1.) / + // vecCore::math::Sqrt + std::sqrt + (momentum[0] * momentum[0] + momentum[1] * momentum[1] + momentum[2] * momentum[2]); + // = vdt::fast_isqrt_general( momentum_sqr, 2); // Alternative + + dpx_ds = (momentum[1] * Bfield[2] - momentum[2] * Bfield[1]); // Ax = a*(Vy*Bz - Vz*By) + dpy_ds = (momentum[2] * Bfield[0] - momentum[0] * Bfield[2]); // Ay = a*(Vz*Bx - Vx*Bz) + dpz_ds = (momentum[0] * Bfield[1] - momentum[1] * Bfield[0]); // Az = a*(Vx*By - Vy*Bx) + + Real_t cof = charge * Real_t(gCof) * inv_momentum_mag; + + dx_ds = momentum[0] * inv_momentum_mag; // (d/ds)x = Vx/V + dy_ds = momentum[1] * inv_momentum_mag; // (d/ds)y = Vy/V + dz_ds = momentum[2] * inv_momentum_mag; // (d/ds)z = Vz/V + + dpx_ds *= cof; + dpy_ds *= cof; + dpz_ds *= cof; + // std::cout << "Force(v2) - mom = " << momentum[0] << " " << momentum[1] << " " << momentum[2] << " B= " << Bfield[0] << " " << Bfield[1] << " " << Bfield[2] << "\n"; + // std::cout << "Force(v2) - d/ds [px] = " << dpx_ds << " [py] = " << dpy_ds << " [pz] = " << dpz_ds << "\n"; +} + +// ------------------------------------------------------------------------------------ + +/*********** +template +template +inline +__host__ __device__ void +MagneticFieldEquation:: + Force(Real_t const & momx, Real_t const & momy, Real_t const & momz, + Real_t const & Bx, Real_t const & By, Real_t const & Bz, + Int_t const & charge, + // Real_t const & momentumMag, // Magnitude of initial momentum + Real_t & dx_ds, Real_t & dy_ds, Real_t & dz_ds, + Real_t & dpx_ds, Real_t & dpy_ds, Real_t & dpz_ds + ) +{ + // using vecgeom::Vector3D; + vecgeom::Vector3D Momentum = { momx, momy, momz } ; + vecgeom::Vector3D Bfield = { Bx, By, Bz }; + // Real_t dydsArr[Nvar]; + Force( Momentum, charge, Bfield, dx_ds, dy_ds, dz_ds, dpx_ds, dpy_ds, dpz_ds ); + + // dx_ds= dydsArr[0]; dy_ds= dydsArr[1]; dz_ds= dydsArr[2]; + // dpx_ds= dydsArr[3]; dpy_ds= dydsArr[4]; dpz_ds= dydsArr[5]; +} +*************/ + +/************************************************************** +{ + Real_t inv_momentum_mag = Real_t(1.) / + // vecCore::math::Sqrt + std::sqrt + (momx * momx + momy * momy + momz * momz); + // = vdt::fast_isqrt_general( momentum_sqr, 2); // Alternative + + dpx_ds = (momy * Bz - momz * By); // Ax = a*(Vy*Bz - Vz*By) + dpy_ds = (momz * Bx - momx * Bz); // Ay = a*(Vz*Bx - Vx*Bz) + dpz_ds = (momx * By - momy * Bx); // Az = a*(Vx*By - Vy*Bx) + + Real_t cof = charge * Real_t(gCof) * inv_momentum_mag; + + dx_ds = momx * inv_momentum_mag; // (d/ds)x = Vx/V + dy_ds = momy * inv_momentum_mag; // (d/ds)y = Vy/V + dz_ds = momz * inv_momentum_mag; // (d/ds)z = Vz/V + + dpx_ds *= cof; + dpy_ds *= cof; + dpz_ds *= cof; + std::cout << "Force - mom = " << momx << " " << momy << " " << momz << " B= " << Bx << " " << By << " " << Bz << "\n"; + std::cout << "Force - d/ds [px] = " << dpx_ds << " [py] = " << dpy_ds << " [pz] = " << dpz_ds << "\n"; +} +************************************************************************/ + +// ------------------------------------------------------------------------------------ + +template +template +inline __host__ __device__ void +MagneticFieldEquation:: +EvaluateDerivatives( MagField_t const & magField, + const Real_t y[], + int charge, + Real_t dy_ds[]) +{ + float Bx, By, Bz; + magField.Evaluate( y[0], y[1], y[2], Bx, By, Bz ); + Real_t Bfield[3] = { Bx, By, Bz } ; + EvaluateDerivativesGivenB(y, Bfield, charge, dy_ds ); +} + +// ------------------------------------------------------------------------------------ + +template +template +inline __host__ __device__ void +MagneticFieldEquation:: +EvaluateDerivativesReturnB( MagField_t const & magField, + const Real_t y[], + int charge, + Real_t dy_ds[], + vecgeom::Vector3D& BfieldVec + // float BfieldValue[3] + ) +{ + // vecgeom::Vector3D position( y[0], y[1], y[2] ); + // magField.Evaluate( position, BfieldVec ); + // float Bfield[3] = { BfieldVec[0], BfieldVec[1], BfieldVec[2] } ; + float Bx, By, Bz; + magField.Evaluate( y[0], y[1], y[2], Bx, By, Bz ); + // std::cout << "EvalDerivRetB: Bx= " << Bx << " By= " << By << " Bz=" << Bz << std::endl; + Real_t Bfield[3] = { Bx, By, Bz } ; + EvaluateDerivativesGivenB(y, Bfield, charge, dy_ds ); + BfieldVec = vecgeom::Vector3D( Bfield[0], Bfield[1], Bfield[2] ); // Bx, By, Bz); +} + +// ------------------------------------------------------------------------------------ + +template +template +inline __host__ __device__ void +MagneticFieldEquation:: +EvaluateDerivativesReturnB( MagField_t const & magField, + vecgeom::Vector3D const & position, + vecgeom::Vector3D const & momentum, + int charge, + Real_t dy_ds[], + vecgeom::Vector3D & BfieldVec + // float BfieldValue[3] + ) +{ + // vecgeom::Vector3D position( y[0], y[1], y[2] ); + // magField.Evaluate( position, BfieldVec ); + // float Bfield[3] = { BfieldVec[0], BfieldVec[1], BfieldVec[2] } ; + const Real_t y[Nvar] = { position[0], position[1], position[2], momentum[0], momentum[1], momentum[2] }; + float Bx, By, Bz; + magField.Evaluate( y[0], y[1], y[2], Bx, By, Bz ); + // std::cout << "EvalDerivRetB: Bx= " << Bx << " By= " << By << " Bz=" << Bz << std::endl; + Real_t Bfield[3] = { Bx, By, Bz } ; + EvaluateDerivativesGivenB(y, Bfield, charge, dy_ds ); + BfieldVec = vecgeom::Vector3D( Bfield[0], Bfield[1], Bfield[2] ); // Bx, By, Bz); +} + + +// ------------------------------------------------------------------------------------ + +template +template +inline __host__ __device__ void +MagneticFieldEquation:: + EvaluateDerivativesGivenB( const Real_t y[], + const Real_t Bfield[3], + int charge, + Real_t dy_ds[] ) +{ + const vecgeom::Vector3D Bvec= { Bfield[0], Bfield[1], Bfield[2] }; + EvaluateRhsGivenB(y, charge, Bvec, dy_ds); +} + +// template +// MagneticFieldEquation:: + +/****** + template< class Real_t> + static __host__ __device__ void + EvaluateField(Real_t const & posx, Real_t const & posy, Real_t const & posz, + MagField_t const & magField, + Real_t & Bx, Real_t & By, Real_t & Bz) + { + magField.Evaluated (posx, posy, posz, Bx, By, Bz ); + } + ******/ + +#endif diff --git a/magneticfield/inc/PrintFieldVectors.h b/magneticfield/inc/PrintFieldVectors.h new file mode 100644 index 00000000..6d5ce4e6 --- /dev/null +++ b/magneticfield/inc/PrintFieldVectors.h @@ -0,0 +1,153 @@ +// SPDX-FileCopyrightText: 2021 CERN +// SPDX-License-Identifier: Apache-2.0 +// +// Author: J. Apostolakis, 19 Nov 2021 +// +#ifndef PRINT_FIELD_VECTORS_H__ +#define PRINT_FIELD_VECTORS_H__ + +static constexpr int NvarPrint= 6; + +#include +#include + +#include "VecGeom/base/Vector3D.h" + +namespace PrintFieldVectors { + + constexpr int prec = 6; + constexpr int nd = prec + 4; + // constexpr char* numFormat= "10.6f"; + +template +static inline __host__ __device__ +void Print3vec(const Real_t y[3], const char* name ) +{ + char numFormat[12]; + char format[64]; + + sprintf( numFormat, "%%%d.%df", nd, prec ); + sprintf( format, " %%-12s = %s %s %s ", numFormat, numFormat, numFormat ); + // std::cout << " Position = " << setw(nd) << y[0] << " " << setw(nd) << y[1] << " " << setw(nd) << y[2]; + printf( format, name, y[0], y[1], y[2] ); +} + +template +static inline __host__ __device__ +void PrintScalar(const Real_t yVal, const char* name[10] ) +{ + char numFormat[12]; + // char format[64]; + + sprintf( numFormat, "%%%d.%df", nd, prec ); + // sprintf( format, " %%s = %s ", numFormat ); + // std::cout << " Position = " << setw(nd) << y[0] << " " << setw(nd) << y[1] << " " << setw(nd) << y[2]; + // printf( format, name, yVal ); + printf( "%12s", name ); + printf( numFormat, yVal ); +} + +template +static inline __host__ __device__ +void PrintSixvec(const Real_t y[NvarPrint], Real_t originalMomentum = -1.0, bool newline=true ) +{ + char numFormat[12]; + sprintf( numFormat, "%%%d.%df", nd, prec ); + + // printf( "\n# Six-Vector & B field \n"); + Print3vec( y, "Position" ); if( newline ) { printf("\n"); } + Print3vec( &y[3], "Momentum" ); if( newline ) { printf("\n"); } + + if( (originalMomentum != -1.0) && (originalMomentum != 0.0) ) { + Real_t mag = sqrt( y[3] * y[3] + y[4] * y[4] + y[5] * y[5] ); + // Vector3D( y[3], y[4], y[5]).Mag(); + // std::cout << " (|p|-o)/o = " << setw(nd) << (mag - originalMomentum) / originalMomentum; + printf( " (|p|-o)/o = " ); + printf( numFormat, (mag - originalMomentum) / originalMomentum ); + } + if( newline ) printf("\n"); +} + +template +static inline +void PrintSixvec( vecgeom::Vector3D const & position, vecgeom::Vector3D const & momentum, Real_t origMomentumMag = -1.0 ) +{ + Real_t y[NvarPrint]= { position[0], position[1], position[2], + momentum[0], momentum[1], momentum[2] }; + PrintSixVec( y, origMomentumMag ); +} + +template +static inline __host__ __device__ +void PrintLineSixvec(const Real_t y[NvarPrint] ) +{ + Print3vec( y, "x:" ); // std::cout << " x: " << setw(nd) << y[0] << " " << setw(nd) << y[1] << " " << setw(nd) << y[2]; + Print3vec( &y[3], "p:" ); // std::cout << " p: " << setw(nd) << y[3] << " " << setw(nd) << y[4] << " " << setw(nd) << y[5]; + printf("\n"); +} + +template +static inline __host__ __device__ +void PrintSixvecAndDyDx(const Real_t y[NvarPrint], int charge, const Real_t Bfield[3], Real_t const dydx[NvarPrint]) +{ + // RightHandSide(y, charge, dydx); + + // Obtain the field value + // Vector3D Bfield; + // FieldFromY(y, Bfield); + // EvaluateRhsGivenB(y, charge, Bfield, dydx); + + PrintSixvec( y ); + // char bFieldName[12] = " B field [0-2]"; + Print3vec( Bfield, /* bFieldName ); */ " B field [0-2]"); printf("\n"); + + // std::cout << "# 'Force' from B field \n"; + Print3vec( dydx , " dy/dx [0-2] (=dX/ds) = " ); printf("\n"); + Print3vec( dydx+3, " dy/dx [3-5] (=dP/ds) = " ); printf("\n"); +} + +// std::cout.unsetf(std::ios_base::scientific); + + +// template static inline // __host__ __device__ +// void PrintSixvecAndDyDx(const Real_t y[NvarPrint], int charge, const Real_t Bfield[3], const Real_t dydx[NvarPrint]) { } + + +template +static inline // __host__ __device__ +void PrintLineSixvecDyDx(const Real_t y[NvarPrint], int charge, const Real_t Bfield[3], Real_t const dydx[NvarPrint]) +{ + using std::cout; + using std::endl; + using std::setw; + constexpr int linPrec = 3; // Output precision - number of significant digits + + int oldprec= std::cout.precision(linPrec); + + PrintLineSixvec( y ); // Was PrintSixvec( y, false ); + + cout.setf(std::ios_base::fixed); + // cout << " B field [0-2] = "; + // cout << " B: "; + // cout << setw(nd) << Bfield[0] << " " << setw(nd) << Bfield[1] << " " << setw(nd) << Bfield[2]; + + // cout << "# 'Force' from B field \n"; + // cout << " dy/dx [0-2] (=dX/ds) = "; + cout << " dy_ds "; + cout << setw(nd) << dydx[0] << " " << setw(nd) << dydx[1] << " " << setw(nd) << dydx[2] << " | " ; + // cout << " dy/dx [3-5] (=dP/ds) = "; + cout << setw(nd+2) << dydx[3] << " " << setw(nd+2) << dydx[4] << " " << setw(nd+2) << dydx[5]; + cout.unsetf(std::ios_base::fixed); + cout << std::endl; + cout.precision(oldprec); +} + +// cout.unsetf(std::ios_base::scientific); + + +// template static inline // __host__ __device__ +// void PrintSixvecAndDyDx(const Real_t y[NvarPrint], int charge, const Real_t Bfield[3], const Real_t dydx[NvarPrint]) { } + +}; + +#endif diff --git a/magneticfield/inc/RkIntegrationDriver.h b/magneticfield/inc/RkIntegrationDriver.h new file mode 100644 index 00000000..3eb357f9 --- /dev/null +++ b/magneticfield/inc/RkIntegrationDriver.h @@ -0,0 +1,456 @@ +// SPDX-FileCopyrightText: 2021 CERN +// SPDX-License-Identifier: Apache-2.0 + +/* + * Created on: November 11-15, 2021 + * + * Author: J. Apostolakis + */ + +#ifndef RKINTEGRATION_DRIVER_H_ +#define RKINTEGRATION_DRIVER_H_ + +// #include "VecGeom/base/Global.h" + +#include +#include // Needed for Vector3D +// #include + +#include "ErrorEstimatorRK.h" + +// namespace adept { +// inline namespace ADEPT_IMPL_NAMESPACE { + +/** + * A simple stepper treating the propagation of particles in a constant magnetic field + * ( not just along the z-axis-- for that there is ConstBzFieldHelixStepper ) + */ + +template + +class RkIntegrationDriver { + +public: + // No constructors! + // ---------------- + // RkIntegrationDriver() = delete; + // RkIntegrationDriver( const RkIntegrationDriver& ) = delete; + // ~RkIntegrationDriver() = delete; + + // template + // Real_t GetCurvature(Vector3D const & dir, float const charge, float const momentum) const; + + + // Propagate the track along the numerical solution of the ODE by the requested length + // Inputs: current position, current direction, some particle properties, max # of steps + // Output: new position, new direction of particle, number of integration steps (current integration) + // + // + // Note: care needed to 'tune' max number of steps - compromise between getting work done & divergence + // + // Versions: + // 1. Vector3D version + // template + template + static inline __host__ __device__ bool Advance(vecgeom::Vector3D & position, + vecgeom::Vector3D & momentumVec, + Int_t const &charge, + // Real_t const &momentum, + Real_t const &step, + MagField_t const &magField, + Real_t &htry, // Suggested integration step -- from previous stages + Real_t dydx_next[], // dy_ds[Nvar] at final point (return only !! ) + Real_t &lengthDone, + unsigned int &totalTrials, + unsigned int maxTrials = 5 + ) ; + // 2. per variable version + // template + static inline __host__ __device__ bool + Advance(Real_t const &posx, Real_t const &posy, Real_t const &posz, + Real_t const &dirx, Real_t const &diry, Real_t const &dirz, + Int_t const &charge, Real_t const &momentum, Real_t const &step, + MagField_t const &magField, + Real_t &htry, // Suggested integration step -- from previous stages + Real_t &newposx, Real_t &newposy, Real_t &newposz, + Real_t &newdirx, Real_t &newdiry, Real_t &newdirz, + Real_t dydx_next[], // dy_ds[] at final point (return only !! ) + Real_t &lengthDone, + unsigned int &totalTrials, + unsigned int maxTrials = 5 + ) ; + + // Versions: + // 1. Original Vector3D version + // template + static inline __host__ __device__ bool + AdvanceV1(vecgeom::Vector3D const &position, + vecgeom::Vector3D const &direction, + Int_t const &charge, + Real_t const &momentum, + Real_t const &step, + MagField_t const &magField, + Real_t &htry, // Suggested integration step, from previous stages + vecgeom::Vector3D &endPosition, + vecgeom::Vector3D &endDirection, + Real_t dydx_next[], // dy_ds[Nvar] at final point (return only !! ) + Real_t &lengthDone, + unsigned int &totalTrials, + unsigned int maxTrials = 5 + ) ; + + // Invariants + // ---------- + static constexpr int Nvar = 6; // For now, adequate to integrate over x, y, z, px, py, pz + static constexpr Real_t fEpsilonRelativeMax = 1.0e-6; // For now .. to be a parameter + static constexpr Real_t fMinimumStep= 2.0e-4 * copcore::units::millimeter; + static constexpr Real_t kSmall = 1.0e-30; // amount to add to vector magnitude to avoid div by zero + + // Auxiliary methods + // ----------------- + // static void SetMaxRelativeError( Real_t epsMax ) { fEpsilonRelativeMax = epsMax; } + + // template + + static inline __host__ __device__ + void PrintStep( vecgeom::Vector3D const &startPosition, + vecgeom::Vector3D const &startDirection, + Int_t const &charge, Real_t const &momentum, Real_t const &step, + vecgeom::Vector3D &endPosition, vecgeom::Vector3D &endDirection); + + static inline __host__ __device__ bool + IntegrateStep( const Real_t yStart[], + const Real_t dydx[], + int charge, + Real_t & xCurrent, // InOut + Real_t htry, + MagField_t magField, + Real_t yEnd[], // Out - values + Real_t next_dydx[], // - next derivative + Real_t & hnext ); + + + // static inline __host__ __device__ EquationType() { return Equation_t; } + +protected: + // template + static inline __host__ __device__ bool CheckModulus(Real_t &newdirX_v, Real_t &newdirY_v, Real_t &newdirZ_v); + +}; // end class declaration + +// ------------------------------------------------------------------------------------------------ + +#include +#include + +template +template +inline __host__ __device__ +bool RkIntegrationDriver +::Advance( vecgeom::Vector3D & position, // In/Out + vecgeom::Vector3D & momentumVec, // In/Out + Int_t const &chargeInt, + // Real_t const &momentum, + Real_t const &length, + MagField_t const &magField, + Real_t &htry, // Suggested integration step -- from previous stages + Real_t dydx_next[Nvar], // dy_ds[] at final point (return only !! ) + Real_t &lenAdvanced, // how far it advanced + unsigned int &totalTrials, // total overall steps for this integration + unsigned int maxTrials // max to do now + ) +{ + using vecgeom::Vector3D; + + Real_t yStart[Nvar] = { position[0], position[1], position[2], + momentumVec[0], momentumVec[1], momentumVec[2] }; + + // vecgeom::Vector3D &B0fieldVal; + // EvaluateField( position, B0fieldVal ); + + Real_t dydx[Nvar]; + Real_t yEnd[Nvar]; + // Real_t charge= Real_t( chargeInt ); + + // Stepper_t::EquationType ??? ToDO + Equation_t::EvaluateDerivatives( magField, yStart, chargeInt, dydx ); + + // std::cout << " RK: i 00 s= " << setw(10) << lenAdvanced; + // PrintFieldVectors::PrintLineSixvecDyDx( yStart, charge, magFieldDummy, dydx ); + + Real_t x = 0.0 ; + Real_t hgood = 0.0; // Length possible with adequate accuracy -- could return this for next step! + + if( htry <= 0.0 || htry >= length ){ + htry = length; + } + + bool done= false; + int numSteps = 0; + + bool allFailed = true; // Have all steps until now failed + bool goodStep= false; // last step was good + do + { + Real_t hnext; // , hdid; + + goodStep= IntegrateStep(yStart, dydx, chargeInt, x, htry, magField, + yEnd, dydx_next, hnext ); + + Real_t hdid = goodStep ? htry : 0.0 ; + allFailed = allFailed && !goodStep; + + lenAdvanced += hdid; + done = (x >= length ); + hgood= vecCore::Max( hnext, Real_t(fMinimumStep) ); + +#ifdef RK_VERBOSE + Real_t htryOld = htry; +#endif + htry = hgood; + if( goodStep && !done ) { + for( int i=0; i< Nvar; i++ ) { + yStart[i] = yEnd[i]; + dydx[i] = dydx_next[i]; // Using FSAL property ! + } + htry = vecCore::Min( hgood , length - x ); + } +#ifdef RK_VERBOSE + if( Verbose > 1 ) { + printf( " n = %4d x = %10.5f ltot = %10.5f ", numSteps, x, lenAdvanced ); + printf( " h: suggested (input try) = %10.7f good? %1s next= %10.7f \n", // did= %10.7f good= %10.7f + htryOld, /* did, */ ( goodStep ? "Y" : "N" ), hnext /* , good */ ); + } +#endif + + ++numSteps; + + } while ( !done && numSteps < maxTrials ); + + totalTrials += numSteps; + + // In case of failed last step, we must not use it's 'yEnd' values !! + if( goodStep ) { + // Real_t invM = 1.0 / (momentum+kSmall); + position.Set( yEnd[0], yEnd[1], yEnd[2] ); + momentumVec.Set( yEnd[3], yEnd[4], yEnd[5] ); + // endDirection.Set( invM*yEnd[3], invM*yEnd[4], invM*yEnd[5] ); + } else { + if( !allFailed ) { // numSteps == maxTrials ){ + position.Set( yStart[0], yStart[1], yStart[2] ); + momentumVec.Set( yStart[3], yStart[4], yStart[5] ); + } + } + + return done; +} + +// ---------------------------------------------------------------------------------------- + + // Versions: + // 1. Vector3D version + // template + +template +inline __host__ __device__ bool +RkIntegrationDriver:: + AdvanceV1( vecgeom::Vector3D const &startPosition, + vecgeom::Vector3D const &startDirection, + Int_t const &charge, + Real_t const &momentum, + Real_t const &step, + MagField_t const &magField, + Real_t &htry, // Suggested integration step -- from previous stages + vecgeom::Vector3D &endPosition, + vecgeom::Vector3D &endDirection, + Real_t dydx_next[], // dy_ds[Nvar] at final point (return only !! ) + Real_t &lengthDone, + unsigned int &totalTrials, + unsigned int maxTrials + ) +{ + vecgeom::Vector3D positionVec = startPosition; + vecgeom::Vector3D momentumVec = momentum * startDirection; + bool done = + Advance( positionVec, momentumVec, charge, step, magField, htry, dydx_next, lengthDone, totalTrials, maxTrials ); + + Real_t invM = 1.0 / (momentum+kSmall); + endPosition = positionVec; + endDirection = invM * momentumVec; + + return done; +} + +// ---------------------------------------------------------------------------------------- + +template +inline __host__ __device__ bool +RkIntegrationDriver:: + IntegrateStep( const Real_t yStart[], + const Real_t dydx[], + int charge, + Real_t & xCurrent, // InOut + Real_t htry, + MagField_t magField, + // Real_t eps_rel_max, + Real_t yEnd[], // Out - values + Real_t next_dydx[], // - next derivative + Real_t & hnext ) +{ + constexpr Real_t safetyFactor= 0.9; + constexpr Real_t shrinkPower = -1.0 / Real_t(Stepper_t::kMethodOrder); + constexpr Real_t growPower = -1.0 / (Real_t(Stepper_t::kMethodOrder+1)); + constexpr Real_t max_step_increase = 10.0; // Step size must not grow more than 10x + constexpr Real_t max_step_decrease = 0.1; // Step size must not shrink more than 10x + // const Real_t errcon = std::pow( max_step_increase / safetyFactor , 1.0/growPower ); + + Real_t yErr[Nvar]; + + // static constexpr iMom0 = 3; // Momentum starts at [3] + // = Mag2 ( Vector3D( yStart[iMom0], yStart[imom0+1], yStart[iMom0+2]) ); + Real_t magMomentumSq = yStart[3] * yStart[3] + yStart[4] * yStart[4] + yStart[5] * yStart[5]; + + Stepper_t::StepWithErrorEstimate(magField, + yStart, + dydx, + charge, + htry, + yEnd, // Output: y values at end, + yErr, // estimated errors, + next_dydx); // next value of dydx + + ErrorEstimatorRK errorEstimator(fEpsilonRelativeMax, fMinimumStep ); + Real_t errmax_sq = errorEstimator.EstimateSquareError( yErr, htry, magMomentumSq ); + + bool goodStep = errmax_sq <= 1.0; + if ( goodStep ) + { + xCurrent += htry; + // Store next_dydx ... +#if 1 + // New code - all threads run the same code, even power ... + hnext = htry * vecCore::Min( safetyFactor*vecCore::Pow(errmax_sq, Real_t(0.5)*growPower ), max_step_increase ); +#else + // Compute size of next Step + hnext = ( errmax_sq > errcon*errcon ) ? h * safetyFactor*std::pow(errmax_sq, Real_t(0.5)*growPower ) + : hnext = max_step_increase*h ; +#endif + } + else + { + // Step failed; compute the size of retrial Step. + Real_t htemp = safetyFactor * htry * std::pow( errmax_sq, 0.5* shrinkPower ); + hnext = vecCore::Max( htemp, max_step_decrease * htry ); + // no more than than a factor of 10 smaller + + if( xCurrent + hnext == xCurrent ) { + // Serious Problem --- under FLOW !!! Report it ?????????????????????????? + hnext = 2.0 * vecCore::math::Max( htemp, htry ); + } + } + return goodStep; +} + +// ---------------------------------------------------------------------------------------- + +/** + * this function propagates the track along the "Runge-Kutta" solution by a step step + * input: current position (x0, y0, z0), current direction ( dirX0, dirY0, dirZ0 ), some particle properties + * output: new position, new direction of particle + */ +template +inline __host__ __device__ bool +RkIntegrationDriver + ::Advance(Real_t const &x0, Real_t const &y0, Real_t const &z0, + Real_t const &dirX0, Real_t const &dirY0, Real_t const &dirZ0, + Int_t const &charge, Real_t const &momentum, Real_t const &step, + MagField_t const &magField, + Real_t &htry, // Suggested integration step -- from previous stages + Real_t &x, Real_t &y, Real_t &z, + Real_t &dx, Real_t &dy, Real_t &dz, + Real_t dydx_next[Nvar], + Real_t &lenAdvanced, // how far it advanced + unsigned int &totalTrials, + unsigned int maxTrials + ) +{ + vecgeom::Vector3D position(x0, y0, z0); + vecgeom::Vector3D startDirection(dirX0, dirY0, dirZ0); + vecgeom::Vector3D endPosition, endDirection; + + // position.Set( x0, y0, z0); + // startDirection.Set( dirX0, dirY0, dirZ0); + + bool done= + Advance(position, startDirection, charge, momentum, step, magField, htry, + endPosition, endDirection, dydx_next, lenAdvanced, maxTrials, totalTrials); + x = endPosition.x(); + y = endPosition.y(); + z = endPosition.z(); + dx = endDirection.x(); + dy = endDirection.y(); + dz = endDirection.z(); + + // PrintStep(position, startDirection, charge, momentum, step, endPosition, endDirection); + return done; +} + +// ---------------------------------------------------------------------------------------- + +template +bool RkIntegrationDriver::CheckModulus(Real_t &newdirX, Real_t &newdirY, Real_t &newdirZ) +{ + constexpr float perMillion = 1.0e-6; + + Real_t modulusDir = newdirX * newdirX + newdirY * newdirY + newdirZ * newdirZ; + typename vecCore::Mask goodDir; + goodDir = vecCore::Abs(modulusDir - Real_t(1.0)) < perMillion; + + bool allGood = vecCore::MaskFull(goodDir); + assert(allGood && "Not all Directions are nearly 1"); + + return allGood; +} + +// ---------------------------------------------------------------------------------------- + +template +inline __host__ __device__ +void RkIntegrationDriver:: + PrintStep(vecgeom::Vector3D const &startPosition, + vecgeom::Vector3D const &startDirection, Int_t const &charge, + Real_t const &momentum, Real_t const &step, + vecgeom::Vector3D &endPosition, + vecgeom::Vector3D &endDirection) +{ + // Debug printing of input & output + printf(" HelixSteper::PrintStep \n"); + const int vectorSize = vecCore::VectorSize(); + Real_t x0, y0, z0, dirX0, dirY0, dirZ0; + Real_t x, y, z, dx, dy, dz; + x0 = startPosition.x(); + y0 = startPosition.y(); + z0 = startPosition.z(); + dirX0 = startDirection.x(); + dirY0 = startDirection.y(); + dirZ0 = startDirection.z(); + x = endPosition.x(); + y = endPosition.y(); + z = endPosition.z(); + dx = endDirection.x(); + dy = endDirection.y(); + dz = endDirection.z(); + for (int i = 0; i < vectorSize; i++) { + printf("Start> Lane= %1d Pos= %8.5f %8.5f %8.5f Dir= %8.5f %8.5f %8.5f ", i, x0, y0, z0, dirX0, dirY0, dirZ0); + printf(" s= %10.6f ", step); // / units::mm ); + printf(" q= %d ", charge); // in e+ units ? + printf(" p= %10.6f ", momentum); // / units::GeV ); + // printf(" ang= %7.5f ", angle ); + printf(" End> Pos= %9.6f %9.6f %9.6f Mom= %9.6f %9.6f %9.6f\n", x, y, z, dx, dy, dz); + } +} + +// } // namespace ADEPT_IMPL_NAMESPACE +// } // namespace adept + +#endif /* RKINTEGRATION_DRIVER_H_ */ diff --git a/magneticfield/inc/UniformMagneticField.h b/magneticfield/inc/UniformMagneticField.h new file mode 100644 index 00000000..a0cc776c --- /dev/null +++ b/magneticfield/inc/UniformMagneticField.h @@ -0,0 +1,70 @@ +// SPDX-FileCopyrightText: 2021 CERN +// SPDX-License-Identifier: Apache-2.0 +// +// Author: J. Apostolakis, 16 Nov 2021 + +#ifndef UniformMagneticField_H__ +#define UniformMagneticField_H__ + +// #include +#include + +class UniformMagneticField // : public VVectorField +{ +public: + // static constexpr int gNumFieldComponents = 3; + // static constexpr bool gFieldChangesEnergy = false; + + /** @brief Constructor providing the constant field value (cartesian) */ + __host__ __device__ + UniformMagneticField(const vecgeom::Vector3D &fieldVector) : fFieldComponents(fieldVector) {} + //: VVectorField(gNumFieldComponents, gFieldChangesEnergy), + + /** @brief Constructor providing the constant field value (spherical) */ + __host__ __device__ + UniformMagneticField(char, double vField, double vTheta, double vPhi); + + /** @brief Destructor */ + __host__ __device__ + ~UniformMagneticField() {} + + /** @brief Copy constructor */ + __host__ __device__ + UniformMagneticField(const UniformMagneticField &p) : fFieldComponents(p.fFieldComponents) {} + // : VVectorField(gNumFieldComponents, gFieldChangesEnergy), + + /** Assignment operator */ + __host__ __device__ + UniformMagneticField &operator=(const UniformMagneticField &p); + + /** @brief Templated field interface */ + template + __host__ __device__ + void Evaluate(const vecgeom::Vector3D & /*position*/, vecgeom::Vector3D &fieldValue) const + { + fieldValue.Set(Real_t(fFieldComponents.x()), Real_t(fFieldComponents.y()), Real_t(fFieldComponents.z())); + } + + /** @brief Templated field interface */ + template + __host__ __device__ + void Evaluate(Real_tp1 /*x*/, Real_tp1 /*y*/, Real_tp1 /*z*/, Real_tp2 &Bx, Real_tp2 &By, Real_tp2 &Bz ) const + { + Bx = Real_tp2(fFieldComponents.x()); + By = Real_tp2(fFieldComponents.y()); + Bz = Real_tp2(fFieldComponents.z()); + } + + /** @brief Interface to evaluate field at location */ + __host__ __device__ + void GetFieldValue(const vecgeom::Vector3D &position, vecgeom::Vector3D &fieldValue) + { + Evaluate(position, fieldValue); + } + +private: + vecgeom::Vector3D fFieldComponents; + +}; + +#endif diff --git a/magneticfield/inc/fieldConstants.h b/magneticfield/inc/fieldConstants.h new file mode 100644 index 00000000..0db5103a --- /dev/null +++ b/magneticfield/inc/fieldConstants.h @@ -0,0 +1,21 @@ +// SPDX-FileCopyrightText: 2020 CERN +// SPDX-License-Identifier: Apache-2.0 + +// Author: J. Apostolakis, 29 Nov 2021 + +#ifndef FIELD_CONSTANTS_H__ +#define FIELD_CONSTANTS_H__ + +namespace fieldConstants { + + static constexpr double kB2C = -0.299792458 * (copcore::units::GeV / (copcore::units::tesla * copcore::units::meter)); + // + static constexpr float deltaIntersection = 1.0e-4 * copcore::units::millimeter; + // Accuracy required for intersection of curved trajectory with surface(s) + + static constexpr float gEpsilonDeflect = 1.E-2 * copcore::units::cm; // Allowable deflection during an integration step + // The difference between the endpoint and the projection of the straight-line path after such a step + // Used to ensure the accuracy of (intersecting with volumes) the curved path is well approximated by chords. +}; + +#endif diff --git a/magneticfield/inc/fieldPropagatorConstBany.h b/magneticfield/inc/fieldPropagatorConstBany.h index 3d74e1a0..f663de7f 100644 --- a/magneticfield/inc/fieldPropagatorConstBany.h +++ b/magneticfield/inc/fieldPropagatorConstBany.h @@ -15,17 +15,18 @@ class fieldPropagatorConstBany { using Precision = vecgeom::Precision; public: - __host__ __device__ void stepInField(ConstFieldHelixStepper &helixAnyB, double kinE, double mass, int charge, - Precision step, vecgeom::Vector3D &position, - vecgeom::Vector3D &direction); + inline __host__ __device__ + void stepInField(ConstFieldHelixStepper &helixAnyB, double kinE, double mass, int charge, + Precision step, vecgeom::Vector3D &position, + vecgeom::Vector3D &direction); }; // ---------------------------------------------------------------------------- -__host__ __device__ void fieldPropagatorConstBany::stepInField(ConstFieldHelixStepper &helixAnyB, double kinE, - double mass, int charge, Precision step, - vecgeom::Vector3D &position, - vecgeom::Vector3D &direction) +inline __host__ __device__ void fieldPropagatorConstBany::stepInField(ConstFieldHelixStepper &helixAnyB, double kinE, + double mass, int charge, Precision step, + vecgeom::Vector3D &position, + vecgeom::Vector3D &direction) { using Precision = vecgeom::Precision; if (charge != 0) { diff --git a/magneticfield/inc/fieldPropagatorConstBz.h b/magneticfield/inc/fieldPropagatorConstBz.h index c8c4f91f..5b98990d 100644 --- a/magneticfield/inc/fieldPropagatorConstBz.h +++ b/magneticfield/inc/fieldPropagatorConstBz.h @@ -33,8 +33,8 @@ class fieldPropagatorConstBz { Vector3D &position, Vector3D &direction, vecgeom::NavStateIndex const ¤t_state, vecgeom::NavStateIndex &new_state, bool &propagated, - Precision safety = 0, const int max_iteration = 100); - + const Precision safety = 0.0, + const int max_iteration = 100); private: Precision BzValue; }; @@ -72,7 +72,7 @@ __host__ __device__ Precision fieldPropagatorConstBz::ComputeSafeLength(Precisio // Direction projection in plane perpendicular to field vector Precision dirxy = sqrt((1 - direction[2]) * (1 + direction[2])); - Precision bend = std::fabs(ConstBzFieldStepper::kB2C * charge * BzValue) / momentumMag; + Precision bend = std::fabs(fieldConstants::kB2C * charge * BzValue) / momentumMag; // R = helix radius, curv = 1./R = curvature in plane perpendicular to the field //Precision curv = bend / (dirxy + 1.e-30); @@ -87,7 +87,7 @@ template __host__ __device__ Precision fieldPropagatorConstBz::ComputeStepAndNextVolume( double kinE, double mass, int charge, Precision physicsStep, vecgeom::Vector3D &position, vecgeom::Vector3D &direction, vecgeom::NavStateIndex const ¤t_state, - vecgeom::NavStateIndex &next_state, bool &propagated, Precision safety, const int max_iterations) + vecgeom::NavStateIndex &next_state, bool &propagated, const vecgeom::Precision safetyIn, const int max_iterations) { using Precision = vecgeom::Precision; #ifdef VECGEOM_FLOAT_PRECISION @@ -99,7 +99,7 @@ __host__ __device__ Precision fieldPropagatorConstBz::ComputeStepAndNextVolume( Precision momentumMag = sqrt(kinE * (kinE + 2 * mass)); // Distance along the track direction to reach the maximum allowed error - Precision safeLength = ComputeSafeLength(momentumMag, charge, direction); + const Precision safeLength = ComputeSafeLength(momentumMag, charge, direction); ConstBzFieldStepper helixBz(BzValue); @@ -114,6 +114,7 @@ __host__ __device__ Precision fieldPropagatorConstBz::ComputeStepAndNextVolume( } else { bool continueIteration = false; + Precision safety = safetyIn; Vector3D safetyOrigin = position; // Prepare next_state in case we skip navigation inside the safety sphere. current_state.CopyTo(&next_state); @@ -121,6 +122,7 @@ __host__ __device__ Precision fieldPropagatorConstBz::ComputeStepAndNextVolume( Precision maxNextSafeMove = safeLength; + bool lastWasZero = false; // Locate the intersection of the curved trajectory and the boundaries of the current // volume (including daughters). do { @@ -153,6 +155,13 @@ __host__ __device__ Precision fieldPropagatorConstBz::ComputeStepAndNextVolume( } } + static constexpr Precision ReduceFactor = 0.1; + static constexpr int ReduceIters = 6; + + if( lastWasZero && chordIters >= ReduceIters ) { + lastWasZero = false; + } + if (move == chordLen) { position = endPosition; direction = endDirection; @@ -164,13 +173,18 @@ __host__ __device__ Precision fieldPropagatorConstBz::ComputeStepAndNextVolume( // FIXME: Even for zero steps, the Navigator will return kPush + possibly // Navigator::kBoundaryPush instead of a real 0. move = 0; + lastWasZero = true; - static constexpr Precision ReduceFactor = 0.5; - static constexpr Precision ReduceIters = 5; // Reduce the step attempted in the next iteration to navigate around // boundaries where the chord step may end in a volume we just left. maxNextSafeMove = ReduceFactor * safeMove; continueIteration = chordIters < ReduceIters; + + if( ! continueIteration ) { + // Let's move to the other side of this boundary -- this side we cannot progress !! + move = Navigator::kBoundaryPush; + // printf("fieldProp-ConstBz: pushing by %10.4g \n ", move ); + } } else { // Accept the intersection point on the surface. This means that // the point at the boundary will be on the 'straight'-line chord, diff --git a/magneticfield/inc/fieldPropagatorRungeKutta.h b/magneticfield/inc/fieldPropagatorRungeKutta.h new file mode 100644 index 00000000..691665b2 --- /dev/null +++ b/magneticfield/inc/fieldPropagatorRungeKutta.h @@ -0,0 +1,385 @@ +// SPDX-FileCopyrightText: 2020 CERN +// SPDX-License-Identifier: Apache-2.0 + +// Author: J. Apostolakis 15 Nov 2021 + +#ifndef FIELD_PROPAGATOR_RUNGEKUTTA_H +#define FIELD_PROPAGATOR_RUNGEKUTTA_H + +#include + +#include "fieldConstants.h" // For kB2C factor with units + +#include "UniformMagneticField.h" +#include "RkIntegrationDriver.h" + +template +class fieldPropagatorRungeKutta { +public: + static + inline __host__ __device__ + void stepInField( Field_t const & magneticField, + Real_t kinE, + Real_t mass, + int charge, + Real_t step, + vecgeom::Vector3D &position, + vecgeom::Vector3D &direction + , int id // For debugging + ); + + static + inline __host__ __device__ + __host__ __device__ Real_t ComputeStepAndNextVolume( + Field_t const & magneticField, + double kinE, double mass, int charge, double physicsStep, + vecgeom::Vector3D &position, + vecgeom::Vector3D &direction, + vecgeom::NavStateIndex const ¤t_state, + vecgeom::NavStateIndex &next_state, + bool & propagated, + const Real_t & /*safety*/, + const int max_iterations + , int & iterDone + , int threadId + ); + // Move the track, + // updating 'position', 'direction', the next state and returning the length moved. +protected: + static inline __host__ __device__ + void IntegrateTrackToEnd( Field_t const & magField, // RkDriver_t &integrationDriverRK, + vecgeom::Vector3D & position, vecgeom::Vector3D & momentumVec, + int charge, Real_t stepLength + , int id // Temporary - for debugging + , bool verbose = true // >>2 + ); + + static inline __host__ __device__ + bool IntegrateTrackForProgress( Field_t const & magField, // RkDriver_t &integrationDriverRK, + vecgeom::Vector3D & position, + vecgeom::Vector3D & momentumVec, + int charge, + Real_t & stepLength // In/Out - In = requested; Out = last trial / next value ?? + // unsigned int & totalTrials + ); + + static constexpr unsigned int fMaxTrials= 100; + static constexpr unsigned int Nvar = 6; // For position (3) and momentum (3) -- invariant + +#ifdef VECGEOM_FLOAT_PRECISION + static constexpr Real_t kPush = 10 * vecgeom::kTolerance; +#else + static constexpr Real_t kPush = 0.; +#endif + + // Cannot change the energy (or momentum magnitude) -- currently usable only for pure magnetic fields +}; + +#define VERBOSE_STEP_IN_THREAD 1 // 2022.09.05 14:00 -- look for failed RK integration hAdvanced = 0.0 + +// ---------------------------------------------------------------------------- +template +inline __host__ __device__ // __host__ __device_ +void fieldPropagatorRungeKutta + ::stepInField( Field_t const & magField, // RkDriver_t &integrationDriverRK, + Real_t kinE, + Real_t mass, + int charge, + Real_t step, + vecgeom::Vector3D &position, + vecgeom::Vector3D &direction + , int id // For debugging + ) +{ + Real_t momentumMag = sqrt(kinE * (kinE + 2.0 * mass)); + Real_t invMomentumMag = 1.0 / momentumMag; + + // Only charged particles ( e-, e+, any .. ) can be propagated + vecgeom::Vector3D positionVec = position; + vecgeom::Vector3D momentumVec = momentumMag * direction; + IntegrateTrackToEnd( magField, + positionVec, + momentumVec, + charge, + step + , id, true // For debugging + ); + position = positionVec; + direction = invMomentumMag * momentumVec; + // Deviation of magnitude of direction from unit indicates integration error +} + +// ---------------------------------------------------------------------------- + +template +inline __host__ __device__ +void fieldPropagatorRungeKutta::IntegrateTrackToEnd( + // RkDriver_t & integrationDriverRK, + Field_t const & magField, + vecgeom::Vector3D & position, + vecgeom::Vector3D & momentumVec, + int charge, + Real_t stepLength + , int id + , bool verbose + ) +{ + // Version 1. Finish the integration of lanes ... + // Future alternative (ToDo): return unfinished intergration, in order to interleave loading of other 'lanes' + + const unsigned int trialsPerCall = vecCore::Min( 30U, fMaxTrials / 2) ; // Parameter that can be tuned + unsigned int totalTrials=0; + static int callNum = 0; + callNum ++; + + Real_t lenRemains = stepLength; + + Real_t hTry = stepLength; // suggested 'good' length per integration step + bool unfinished = true; + + Real_t totLen = 0.0; + unsigned int loopCt=0; + do { + Real_t hAdvanced = 0; // length integrated this iteration (of do-while) + Real_t dydx_end[Nvar]; + + bool done= + RkDriver_t::Advance( position, momentumVec, charge, lenRemains, magField, hTry, dydx_end, + hAdvanced, totalTrials, + // id, // Temporary ? + trialsPerCall); + // Runge-Kutta single call ( number of steps <= trialsPerCall ) + + lenRemains -= hAdvanced; + unfinished = lenRemains > 0.0; /// Was = !done ... for debugging ???? + + totLen+= hAdvanced; + loopCt++; + + // sumAdvanced += hAdvanced; // Gravy .. + + } while ( unfinished && (totalTrials < fMaxTrials) ); + + + if( loopCt > 1 && verbose ) { printf( " fieldPropagatorRK: id %3d call %4d --- LoopCt reached %d ", id, callNum, loopCt ); } + +} + +// ---------------------------------------------------------------------------- + +template +inline __host__ __device__ +bool fieldPropagatorRungeKutta::IntegrateTrackForProgress( + Field_t const & magField, + vecgeom::Vector3D & position, + vecgeom::Vector3D & momentumVec, + int charge, + Real_t & stepLength // In/Out - In = requested; Out = last trial / next value ?? + // , unsigned int & totalTrials + ) +{ + // Version 2. Try to get some progress in the integration of this threads - but not more than trialsPerCall ... + // Future alternative (ToDo): return unfinished intergration, in order to interleave loading of other 'lanes' + const unsigned int trialsPerCall = vecCore::Min( 6U, fMaxTrials / 2) ; // Parameter that can be tuned + + Real_t lenRemains = stepLength; + + Real_t hTry = stepLength; // suggested 'good' length per integration step + bool done= false; + + int totalTrials= 0; + Real_t hAdvanced = 0; // length integrated this iteration (of do-while) + do { + + Real_t dydx_end[Nvar]; + + done= RkDriver_t::Advance( position, momentumVec, charge, lenRemains, magField, + hTry, dydx_end, hAdvanced, totalTrials, trialsPerCall); + // Runge-Kutta one call for 1+ iterations ( number of steps <= trialsPerCall ) + + stepLength -= hAdvanced; + // unfinished = !done; // (lenRemains > 0.0); + + } while ( hAdvanced == 0.0 && totalTrials < fMaxTrials ); + + return done; +} + +template +inline __host__ __device__ Real_t +inverseCurvature( + vecgeom::Vector3D & momentumVec, + vecgeom::Vector3D & BfieldVec, + int charge + ) +{ + Real_t bmag2 = BfieldVec.Mag2(); + Real_t ratioOverFld = (bmag2>0) ? momentumVec.Dot(BfieldVec) / bmag2 : 0.0; + vecgeom::Vector3D PtransB = momentumVec - ratioOverFld * BfieldVec; + + Real_t bmag = sqrt(bmag2); + + // Real_t curv = fabs(Track::kB2C * charge * bmag / ( PtransB.Mag() + tiny)); + + // Calculate inverse curvature instead - save a division + Real_t inv_curv = fabs( PtransB.Mag() + / ( fieldConstants::kB2C * Real_t(charge) * bmag + 1.0e-30) + ); + return inv_curv; +} + +// Determine the step along curved trajectory for charged particles in a field. +// ( Same name as as navigator method. ) + +template +inline __host__ __device__ Real_t +fieldPropagatorRungeKutta ::ComputeStepAndNextVolume( + Field_t const & magField, + double kinE, double mass, int charge, double physicsStep, + vecgeom::Vector3D & position, + vecgeom::Vector3D & direction, + vecgeom::NavStateIndex const & current_state, + vecgeom::NavStateIndex & next_state, + bool & propagated, + const Real_t & /*safety*/, // eventually In/Out ? + const int max_iterations + , int & itersDone // useful for now - to monitor and report -- unclear if needed later + , int indx + ) +{ + // using copcore::units::MeV; + + const Real_t momentumMag = sqrt(kinE * (kinE + 2.0 * mass)); + vecgeom::Vector3D momentumVec = momentumMag * direction; + + vecgeom::Vector3D B0fieldVec = { 0.0, 0.0, 0.0 }; // Field value at starting point + magField.Evaluate( position, B0fieldVec ); + + Real_t inv_curv = inverseCurvature /**/ ( momentumVec, B0fieldVec, charge ); + + // acceptable lateral error from field ~ related to delta_chord sagital distance + const Real_t safeLength = + sqrt( Real_t(2.0) * fieldConstants::gEpsilonDeflect * inv_curv); // max length along curve for deflectionn + // = sqrt( 2.0 / ( invEpsD * curv) ); // Candidate for fast inv-sqrt + + Precision maxNextSafeMove = safeLength; // It can be reduced if, at the start, a boundary is encountered + + Real_t stepDone = 0.0; + Real_t remains = physicsStep; + const Real_t tiniest_step = 1.0e-7 * physicsStep; // Ignore remainder if < e_s * PhysicsStep + int chordIters = 0; + + constexpr bool inZeroFieldRegion= false; // This could be a per-region flag ... - better depend on template parameter? + bool found_end = false; + + if ( inZeroFieldRegion ) { + stepDone = Navigator_t::ComputeStepAndNextVolume(position, direction, remains, current_state, next_state, kPush); + position += stepDone * direction; + } else { + bool continueIteration = false; + bool fullChord = false; + // vecgeom::Vector3D momentumVec = momentumMag * direction; + const Real_t inv_momentumMag = 1.0 / momentumMag; + + bool lastWasZero = false; // Debug only ? JA 2022.09.05 + + // Locate the intersection of the curved trajectory and the boundaries of the current + // volume (including daughters). + do { + static constexpr Precision ReduceFactor = 0.1; + static constexpr int ReduceIters = 6; + + vecgeom::Vector3D endPosition = position; + vecgeom::Vector3D endMomentumVec = momentumVec; // momentumMag * direction; + const Real_t safeArc = min(remains, maxNextSafeMove); // safeLength); + + IntegrateTrackToEnd( magField, endPosition, endMomentumVec, charge, safeArc, indx); + //----------------- + vecgeom::Vector3D chordVec = endPosition - position; + Real_t chordDist = chordVec.Length(); + vecgeom::Vector3D endDirection = inv_momentumMag * endMomentumVec; + chordVec *= (1.0 / chordDist); // Now the direction of the chord! + + // Check Intersection + //-- vecgeom::Vector3D ChordDir= (1.0/chordDist) * ChordVec; + Real_t linearStep = Navigator_t::ComputeStepAndNextVolume(position, chordVec, chordDist, current_state, next_state, kPush); + Real_t curvedStep; + + if( lastWasZero && chordIters >= ReduceIters ) { + lastWasZero = false; + } + + fullChord = (linearStep == chordDist); + if (fullChord) { + position = endPosition; + momentumVec = endMomentumVec; + + direction = endDirection; + curvedStep = safeArc; + + maxNextSafeMove = safeArc; // Reset it, once a step succeeds!! + continueIteration = true; + } else if (linearStep <= kPush + Navigator_t::kBoundaryPush && stepDone == 0) { + // Cope with a track at a boundary that wants to bend back into the previous + // volume in the first step (by reducing the attempted distance.) + // FIXME: Even for zero steps, the Navigator will return kPush + possibly + // Navigator::kBoundaryPush instead of a real 0. + curvedStep = 0; + lastWasZero = true; + + // Reduce the step attempted in the next iteration to navigate around + // boundaries where the chord step may end in a volume we just left. + maxNextSafeMove = ReduceFactor * safeArc; + continueIteration = chordIters < ReduceIters; + + if( ! continueIteration ) { + // Let's move to the other side of this boundary -- this side we cannot progress !! + curvedStep = Navigator_t::kBoundaryPush; + } + } + else + { + assert( next_state.IsOnBoundary() ); + // assert( linearStep == chordDist ); + + // USE the intersection point on the chord & surface as the 'solution', ie. instead + // of the (potential) true point on the intersection of the curve and the boundary. + // ( This involves a bias -- typically important only for muons in trackers. + // Currently it's controlled/limited by the acceptable step size ie. 'safeLength' ) + Real_t fraction = chordDist > 0 ? linearStep / chordDist : 0.0; + curvedStep = fraction * safeArc; +#ifndef ENDPOINT_ON_CURVE + // Primitive approximation of end direction and linearStep to the crossing point ... + position = position + linearStep * chordVec; + direction = direction * (1.0 - fraction) + endDirection * fraction; + direction = direction.Unit(); + momentumVec = momentumMag * direction; + // safeArc is how much the track would have been moved if not hitting the boundary + // We approximate the actual reduction along the curved trajectory to be the same + // as the reduction of the full chord due to the boundary crossing. +#else + // Alternative approximation of end position & direction -- calling RK again + // Better accuracy (e.g. for comparing with Helix) -- but the point will not be on the surface !! + IntegrateTrackToEnd( magField, position, momentumVec, charge, curvedStep, indx); + direction = inv_momentumMag * momentumVec; // momentumVec.Unit(); +#endif + continueIteration = false; + } + + stepDone += curvedStep; + remains -= curvedStep; + chordIters++; + + found_end = ( (curvedStep > 0) && next_state.IsOnBoundary() ) // Fix 2022.09.05 JA + || (remains <= tiniest_step); + + } while ( !found_end && continueIteration && (chordIters < max_iterations) ); + } + + propagated = found_end; + itersDone += chordIters; + // = (chordIters < max_iterations); // ---> Misses success on the last step! + return stepDone; +} + +#endif diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index e28fd996..93db24f7 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -30,6 +30,7 @@ set(ADEPT_UNIT_TESTS_BASE test_track_block.cu # Unit test for BlockData test_sparsevector.cu # Unit test for SparseVector on GPU test_sparsevector_cpu.cpp # Unit test for SparseVector on CPU + test_magfieldRK.cpp # Unit test for Mag-Field integration classes ) add_compile_options("$<$:--extended-lambda;>") diff --git a/test/test_magfieldRK.cpp b/test/test_magfieldRK.cpp new file mode 100644 index 00000000..2f262675 --- /dev/null +++ b/test/test_magfieldRK.cpp @@ -0,0 +1,522 @@ +// SPDX-FileCopyrightText: 2021 CERN +// SPDX-License-Identifier: Apache-2.0 +// +// Author: J. Apostolakis, 17 Nov 2021 +// + +// #include + +#define __device__ +#define __host__ + +#include + +#include "MagneticFieldEquation.h" + +#include "UniformMagneticField.h" + +#include "PrintFieldVectors.h" + +using Real_t = float; +// using Real_t = double; +using Field_t = UniformMagneticField; +using MagFieldEquation_t= MagneticFieldEquation; + +template +using Vector3D = vecgeom::Vector3D; + +#include + +using std::cout; +using std::cerr; +using std::endl; +using std::setw; + +using namespace vecCore::math; + +// MagFieldEquation_t *CreateUniformFieldAndEquation(Vector3D const &fieldValue); +// MagFieldEquation_t *CreateFieldAndEquation(const char *filename); + +// CreateUniformField(); + +template // , typename Equation_t> +bool TestEquation( Field_t const & magField); + +constexpr unsigned int gNposmom = 6; // Position 3-vec + Momentum 3-vec + + +int gVerbose = 1; + +template // , typename Equation_t> +bool TestEquation( Field_t const & magField) +{ + using Equation_t = MagneticFieldEquation; + + const Real_t perMillion = Real_t(1e-6); + bool hasError = false; // Return value + using Bool_v = vecCore::Mask; + + // 1. Initialize meaningful state + // + Vector3D PositionVec(1., 2., 3.); // initial + Vector3D MomentumVec(0., 0.1, 1.); + Vector3D FieldVec(0., 0., 1.); // Magnetic field value (constant) + + // double PositionTime[4] = { PositionVec.x(), PositionVec.y(), PositionVec.z(), 0.0}; + + // 2. Initialise object for Equation (of motion) + // + Real_t PositionMomentum[gNposmom] = { PositionVec[0], PositionVec[1], PositionVec[2], + MomentumVec[0], MomentumVec[1], MomentumVec[2] } ; + Real_t dydx[gNposmom]; + int charge = -1; + + // 3. Evaluate RHS / ie Derivatives of Equation of motion + // + + // Option 3a) use a field value + Equation_t::EvaluateDerivatives(PositionMomentum, FieldVec, charge, dydx); + // **************************** + + // Option 3b) use a magnetic field class + // Equation_t::EvaluateDerivatives(PositionMomentum, charge, magField, dydx); + //-- - Better form of method: + //-- Equation_t::EvaluateDerivatives(magField, PositionMomentum, charge, dydx); + + // **************************** + + Real_t FieldArr[3] = { FieldVec[0], FieldVec[1], FieldVec[2] }; + std::cout << " Values obtained: " << std::endl; + PrintFieldVectors::PrintSixvecAndDyDx( PositionMomentum, charge, FieldArr, dydx ); + + // 4. Test output values + // + Vector3D dPos_ds(dydx[0], dydx[1], dydx[2]); + Real_t mag = dPos_ds.Mag(); + Real_t magDiff = fabs ( mag - 1.0 ); + assert ( magDiff < 1.0e-6 ); + if( magDiff > 2.0e-7 ) { + std::cerr << "WARNING> dPos_ps (direction) is not a unit vector. | Mag - 1 | = " << magDiff << std::endl; + } + + Vector3D ForceVec(dydx[3], dydx[4], dydx[5]); + + // Check result + Real_t MdotF = MomentumVec.Dot(ForceVec); + Real_t BdotF = FieldVec.Dot(ForceVec); + + Real_t momentumMag = MomentumVec.Mag(); + Real_t fieldMag = FieldVec.Mag(); + Real_t sineAngle = FieldVec.Cross(MomentumVec).Mag() / (momentumMag * fieldMag); + + Real_t ForceMag = ForceVec.Mag(); + const Real_t c = Real_t(copcore::units::kCLight); + + // Tolerance of difference in values (used below) + Real_t tolerance = perMillion; + + if (gVerbose) { + std::cout << " Output: " << std::endl; + } + Bool_v error = (Abs(ForceMag - c * Abs(charge) * fieldMag * sineAngle)) > (tolerance * ForceMag); + if (!vecCore::MaskEmpty(error)) { + cerr << "ERROR: Force magnitude is not equal to c * |charge| * |field| * sin( p, B )." << endl; + cerr << " Force magnitude = " << ForceMag << endl; + cerr << " other side = " << c * Abs(charge) * fieldMag * sineAngle; + cerr << " charge = " << charge << " field-Mag= " << fieldMag << std::endl; + cerr << " Force = " << ForceVec[0] << " " << ForceVec[1] << " " << ForceVec[2] << " " << endl; + hasError = true; + } + + error = (ForceMag == momentumMag * fieldMag); + assert(vecCore::MaskEmpty(error)); // Must add coefficient !! + + error = Abs(MdotF) > tolerance * MomentumVec.Mag() * ForceVec.Mag(); + if (!vecCore::MaskEmpty(error)) { + cerr << "ERROR: Force due to magnetic field is not perpendicular to momentum!!" << endl; + hasError = true; + } else if (gVerbose) { + cout << " Success: Good (near zero) dot product momentum . force " << endl; + } + + error = Abs(BdotF) > tolerance * FieldVec.Mag() * ForceVec.Mag(); + if (!vecCore::MaskEmpty(error)) { + cerr << "ERROR: Force due to magnetic field is not perpendicular to B field!" << std::endl; + cerr << " Vectors: BField Force " << std::endl; + for (int i = 0; i < 3; i++) + cerr << " [" << i << "] " << FieldVec[i] << " " << ForceVec[i] << std::endl; + + hasError = true; + } else if (gVerbose) { + cout << " Success: Good (near zero) dot product magnetic-field . force " << std::endl; + } + + return !hasError; +} + +// NEW Development Thurs 18 December 2021 +// --------------- +#include "DormandPrinceRK45.h" +#include + +template // , typename Equation_t> +bool TestStepper( Field_t const & magField) +// ----------- +{ + constexpr unsigned int Nvar = 6; // Number of integration variables + using Equation_t = MagneticFieldEquation; + using Stepper_t = DormandPrinceRK45; + + Real_t yPosMom[Nvar] = { 0, 1, 2, // PositionVec[0], PositionVec[1], PositionVec[2], + 0, 0, 3.0 * copcore::units::GeV } ; // MomentumVec[0], MomentumVec[1], MomentumVec[2] } ; + + Real_t dy_ds[Nvar]; + + Real_t yOut[Nvar]; // Output: y values at end, + Real_t yerr[Nvar]; // estimated errors, + Real_t next_dyds[Nvar]; // next value of dydx + + int charge = -1; + Real_t step = 0.1; + + // Evaluate dy_ds + // Equation_t::EvaluateDerivatives(magField, yPosMom, charge, dy_ds); + Vector3D magFieldVec( 1, 2, 3); + Equation_t::EvaluateDerivativesReturnB(magField, yPosMom, charge, dy_ds, magFieldVec); + std::cout << " magField = " << magFieldVec[0] << " , " << magFieldVec[1] << " , " << magFieldVec[2] << std::endl; + + Stepper_t::StepWithErrorEstimate(magField, yPosMom, dy_ds, charge, step, yOut, yerr, next_dyds); + + float magFieldOut[3] = { magFieldVec[0], magFieldVec[1], magFieldVec[2] }; + + PrintFieldVectors::PrintSixvecAndDyDx( yPosMom, charge, magFieldOut, dy_ds); + + return true; +} + + +#include "RkIntegrationDriver.h" + +// NEW Development Friday 19 November 2021 +// --------------- + +const unsigned int trialsPerCall = 6; // Try 100 later + +template // , typename Equation_t> +bool TestDriverAdvance( Field_t const & magField, Real_t hLength = 300 ) +// ----------- +{ + constexpr unsigned int Nvar = 6; // Number of integration variables + using Equation_t = MagneticFieldEquation; + using Stepper_t = DormandPrinceRK45; + using Driver_t = RkIntegrationDriver; + + Vector3D Position = { 0, 0.001, 0.002 }; // PositionVec[0], PositionVec[1], PositionVec[2], + Vector3D Direction = { 0, 1.0, 0.0 }; + Real_t momentumMag = 3.0 * copcore::units::GeV ; + + Vector3D Momentum = momentumMag * Direction; + + Real_t dy_ds[Nvar]; + Real_t yOut[Nvar]; // Output: y values at end, + Real_t yerr[Nvar]; // estimated errors, + Real_t next_dyds[Nvar]; // next value of dydx + + int charge = -1; + + Real_t yStart[Nvar] = { Position[0], Position[1], Position[2], + Momentum[0], Momentum[1], Momentum[2] } ; + // std::cout << "yStart: " << std::endl; + // PrintFieldVectors::PrintSixvec( yStart ); + + Vector3D magFieldStart( 0, 0, 0); + Equation_t::EvaluateDerivativesReturnB(magField, yStart, charge, dy_ds, magFieldStart); + // std::cout << " magField = " << magFieldStart[0] << " , " << magFieldStart[1] << " , " << magFieldStart[2] << std::endl; + + // Do a couple of steps + unsigned int totalTrials = 0; + + // Test the simple Advance method + bool verbose = false; + double sumAdvanced = 0.0; + + bool done=false; + Real_t hTry = hLength; // suggested 'good' length per integration step + Real_t dydx_end[Nvar]; + + std::cout << " t: " << setw(3) << totalTrials << " s= " << setw(9) << sumAdvanced << " "; + Real_t magFieldStartArr[3] = { magFieldStart[0], magFieldStart[1], magFieldStart[2] }; + PrintFieldVectors::PrintLineSixvecDyDx( yStart, charge, magFieldStartArr, dy_ds ); + + Vector3D momentumVec = momentumMag * Direction; + + do { + Real_t hAdvanced= 0.0; + + done = + Driver_t::Advance( Position, momentumVec, charge, hLength, + magField, hTry, dydx_end, + hAdvanced, totalTrials, trialsPerCall); + // Runge-Kutta single call + + std::cout << "Advanced returned: done= " << ( done ? "Yes" : " No" ) + << " hAdvanced = " << hAdvanced + << " hNext = " << hTry + << std::endl; + + Real_t yPosMom[Nvar] = { Position[0], Position[1], Position[2], + momentumVec[0], momentumVec[1], momentumVec[2] } ; + sumAdvanced += hAdvanced; + + if( verbose ) + std::cout << "- Trials: max/call = " << trialsPerCall << " total done = " << totalTrials + << ( (totalTrials == trialsPerCall) ? " MAXimum reached !! " : " " ) + << std::endl; + else + std::cout << " t: " << setw(4) << totalTrials << " s= " << setw(9) << sumAdvanced << " "; + + Vector3D magFieldEnd; + Equation_t::EvaluateDerivativesReturnB(magField, yPosMom, charge, dy_ds, magFieldEnd); + Real_t magFieldEndArr[3] = { magFieldEnd[0], magFieldEnd[1], magFieldEnd[2] }; + + //--- PrintFieldVectors::PrintSixvecAndDyDx( + PrintFieldVectors::PrintLineSixvecDyDx( yPosMom, charge, magFieldEndArr, // dydx_end ); + dy_ds); + + hLength -= hAdvanced; + done = (hLength <= 0.0); + + } while ( ! done ); + + Vector3D magFieldEnd; + magField.Evaluate( Position, magFieldEnd ); + + std::cout << " Mag field at end = " << magFieldEnd[0] << " , " << magFieldEnd[1] << " , " << magFieldEnd[2] << std::endl; + + // Test the 'full' DoStep methodx + + + // Finish the integration ... + + return true; +} + +// ------------------------------------------------------------------------------------------------ + +#include "ConstFieldHelixStepper.h" + +// NEW test Code Friday 26 November 2021 +// ------------- + +template // , typename Equation_t> +bool CheckDriverVsHelix( Field_t const & magField, const Real_t stepLength = 300.0 ) +// ----------- +{ + constexpr unsigned int Nvar = 6; // Number of integration variables + using Equation_t = MagneticFieldEquation; + using Stepper_t = DormandPrinceRK45; + using Driver_t = RkIntegrationDriver; + + Vector3D const startPosition = { 0, 0.001, 0.002 }; // PositionVec[0], PositionVec[1], PositionVec[2], + Vector3D const startDirection = { 0, 1.0, 0.0 }; + const Real_t momentumMag = 30.0 * copcore::units::MeV; + + Vector3D const startMomentum = momentumMag * startDirection; + + Real_t dy_ds[Nvar]; + Real_t yOut[Nvar]; // Output: y values at end, + Real_t yerr[Nvar]; // estimated errors, + Real_t next_dyds[Nvar]; // next value of dydx + + int charge = -1; + + Vector3D magFieldStart( 1, 2, 3); + Real_t yStart[Nvar] = { startPosition[0], startPosition[1], startPosition[2], + startMomentum[0], startMomentum[1], startMomentum[2] } ; + Equation_t::EvaluateDerivativesReturnB(magField, yStart, charge, dy_ds, magFieldStart); + // std::cout << " magField = " << magFieldStart[0] << " , " << magFieldStart[1] << " , " << magFieldStart[2] << std::endl; + + // For use by Runge Kutta + Real_t hLength = stepLength; + Vector3D Position = startPosition; + Vector3D Direction= startDirection; + + // -- Test the simple Advance method + unsigned int totalTrials = 0; + bool verbose = false; + Real_t dydx_end[Nvar]; + + if( verbose ) { + // std::cout << "yStart: " << std::endl; + // PrintFieldVectors::PrintSixvec( yStart ); + std::cout << " t: " << setw(3) << totalTrials << " s= " << setw(9) << 0.0 << " "; + Real_t magFieldStartArr[3] = { magFieldStart[0], magFieldStart[1], magFieldStart[2] }; + PrintFieldVectors::PrintLineSixvecDyDx( yStart, charge, magFieldStartArr, dy_ds ); + } + + bool unfinished = true; + Real_t hTry = hLength; // suggested 'good' length per integration step + Real_t sumAdvanced = 0; // length integrated + + constexpr int maxTrials = 500; + + Vector3D momentumVec= momentumMag * Direction; + + do { + Real_t hAdvanced = 0; // length integrated + + bool done = + Driver_t::Advance( Position, momentumVec, charge, hLength, magField, hTry, dydx_end, + hAdvanced, totalTrials, trialsPerCall); + // Runge-Kutta single call ( number of steps <= trialsPerCall ) + + hLength -= hAdvanced; + unfinished = !done; // (hLength > 0.0); + + sumAdvanced += hAdvanced; // Gravy .. + + } while ( unfinished && (totalTrials < maxTrials) ); + // Ensure that the 'lane' integration is done .... or at least we tried hard + + // Cast the solution into simple variables - for comparison + Vector3D endPositionRK = Position; + Vector3D endDirectionRK = (1.0 / momentumMag) * momentumVec; + Real_t magDirectionRK = endDirectionRK.Mag(); + + // For use by Helix method ... + Vector3D endPositionHx, endDirectionHx; + + ConstFieldHelixStepper helixBvec( magFieldStart ); + //******************** + + // Helix solution --- for comparison +#if 1 + // fieldPropagatorConstBany fieldPropagatorBany; + helixBvec.DoStep(startPosition, startDirection, charge, momentumMag, stepLength, + /**************/ endPositionHx, endDirectionHx); +#else + helixBvec.DoStep(startPosition, startDirection, charge, momentumMag, + stepLength, endPositionHx, endDirectionHx ); +#endif + + // Compare the results HERE + Vector3D shiftVec = ( endPositionRK - endPositionHx ); + Vector3D deltaDir = ( endDirectionRK - endDirectionHx ); + + if( verbose ) { + std::cout << "- Trials: max/call= " << trialsPerCall + << " total done = " << totalTrials + << ( (totalTrials == trialsPerCall) ? " MAXimum reached !!":" ") + << std::endl; + } else { + std::cout << " t: " << setw(4) << totalTrials + << " s= " << setw(7) << sumAdvanced << " " + << " l= " << setw(7) << stepLength << " " + << " d||*1e6 = " << setw(10) << 1.e6 * ( magDirectionRK - 1.0 ) << " | " ; + } + Real_t yPosMom[Nvar] = { endPositionHx[0], endPositionHx[1], endPositionHx[2], + momentumVec[0], momentumVec[1], momentumVec[2] } ; + Vector3D magFieldEnd; + Equation_t::EvaluateDerivativesReturnB(magField, yPosMom, charge, dy_ds, magFieldEnd ); + Real_t magFieldEndArr[3] = { magFieldEnd[0], magFieldEnd[1], magFieldEnd[2] }; + + constexpr Real_t magFactorPos= 1.0 / Driver_t:: // RkIntegrationDriver<> + fEpsilonRelativeMax; + constexpr Real_t magFactorDir= magFactorPos ; // 1.0e+6; + + Real_t mulFactorPos = magFactorPos / stepLength; + Real_t deltaPosMom[Nvar]= + { mulFactorPos * shiftVec[0], mulFactorPos * shiftVec[1], mulFactorPos * shiftVec[2], + magFactorDir * deltaDir[0], magFactorDir * deltaDir[1], magFactorDir * deltaDir[2] } ; + + //--- PrintFieldVectors::PrintSixvecAndDyDx( + PrintFieldVectors::PrintLineSixvecDyDx( yPosMom, charge, magFieldEndArr, deltaPosMom ); + + // magField.Evaluate( endPositionHx, magFieldEnd ); + // std::cout << " Mag field at end = " << magFieldEnd[0] << " , " << magFieldEnd[1] + // << " , " << magFieldEnd[2] << std::endl; + + // Test the 'full' DoStep methodx + + // Finish the integration ... + + return true; +} + +//-------------------------------------------------------------------------------------------------- + +int main(int argc, char **argv) +{ + using copcore::units::tesla; + Vector3D fieldValueZ(0.0, 0.0, 1.0*tesla); + bool allStepGood; + + UniformMagneticField *pConstantBfield = new UniformMagneticField(fieldValueZ); + UniformMagneticField Bx( Vector3D(10.0*tesla, 0.0, 0.0) ); + UniformMagneticField By( Vector3D(0.0, 10.0*tesla, 0.0) ); + UniformMagneticField Bz( Vector3D(0.0, 0.0, 10.0*tesla) ); + + cout << " -- Testing mag-field equation with float. " << endl; + cout << " ---------------------------------------------------------------------" << endl; + bool okUniformFloat = TestEquation(*pConstantBfield); + bool okBxFloat = TestEquation(Bx); + bool okByFloat = TestEquation(By); + bool okBzFloat = TestEquation(Bz); + bool allEqGood = okUniformFloat && okBxFloat && okByFloat && okBzFloat; + cout << endl; + + cout << " ---------------------------------------------------------------------" << endl; + cout << " -- Testing Dormand-Prince stepper with magfield and equation (float). " << endl; + cout << " ---------------------------------------------------------------------" << endl; + bool okStep1Float = TestStepper(*pConstantBfield); + allStepGood = okStep1Float; +#ifdef MORE + bool okStep2Float = TestStepper(Bx); + bool okStep3Float = TestStepper(By); + bool okStep4Float = TestStepper(Bz); + allStepGood &&= okStep2Float && okStep3Float && okStep4Float ; +#endif + cout << endl; + + // float lenInc= 250; + cout << " -------------------------------------------------------------------------------" << endl; + cout << " -- Testing Integration driver 'Advance' with Dormand-Prince stepper (float). --" << endl; + cout << " -------------------------------------------------------------------------------" << endl; + + using copcore::units::millimeter; + bool okDoPri1Float = TestDriverAdvance(*pConstantBfield, 1000.0 * millimeter); + + bool allDoPriGood = okDoPri1Float; + for ( float len = 250; len < 5000 ; len *= 2.0 ) + { + cout << " -- Bx field -- \n"; + bool okDoPri2Float = TestDriverAdvance(Bx, len); + + cout << " -- By field -- \n"; + bool okDoPri3Float = TestDriverAdvance(By, 2.0* len); + + cout << " -- Bz field -- \n"; + bool okDoPri4Float = TestDriverAdvance(Bz, 3.0* len); + + allDoPriGood |= okDoPri1Float && okDoPri2Float && okDoPri3Float && okDoPri4Float; + } + + cout << " -------------------------------------------------------------------------------" << endl; + cout << " -- Testing Integration driver 'Advance' vs Helix stepper (float). --" << endl; + cout << " -------------------------------------------------------------------------------" << endl; + bool allChecksGood = true; + + for ( float len = 10 ; len < 1000000 ; len *= 2.0 ) { + bool okCheck1= CheckDriverVsHelix( *pConstantBfield, len ); + allChecksGood = allChecksGood && okCheck1; + } + cout << " Checks vs Helix result = " << ( allChecksGood ? " OK ! " : " Errors Found " ) << "\n"; + + for ( double len = 10 ; len < 1000000 ; len *= 2.0 ) { + bool okCheck1= CheckDriverVsHelix ( *pConstantBfield, len ); + allChecksGood = allChecksGood && okCheck1; + } + + + return ! ( allStepGood && allDoPriGood && allChecksGood ); +}