Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2026-09-15 08:28:57 

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 /// \file SteppingAction.cc
0027 /// \brief Implementation of the SteppingAction class
0028 
0029 #include "SteppingAction.hh"
0030 
0031 #include "Run.hh"
0032 
0033 #include "G4DecayProducts.hh"
0034 #include "G4DecayTable.hh"
0035 #include "G4LossTableManager.hh"
0036 #include "G4ParticleDefinition.hh"
0037 #include "G4ParticleTypes.hh"
0038 #include "G4Step.hh"
0039 #include "G4StepPoint.hh"
0040 #include "G4SystemOfUnits.hh"
0041 #include "G4TouchableHistory.hh"
0042 #include "G4Track.hh"
0043 #include "G4VDecayChannel.hh"
0044 #include "G4VPhysicalVolume.hh"
0045 #include "G4VTouchable.hh"
0046 
0047 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0048 
0049 SteppingAction::SteppingAction() : G4UserSteppingAction()
0050 {
0051   Initialize();
0052 }
0053 
0054 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0055 
0056 SteppingAction::~SteppingAction() {}
0057 
0058 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0059 
0060 void SteppingAction::Initialize()
0061 {
0062   // Initialization needed at the beginning of each Run
0063   fRunPtr = nullptr;
0064   fToleranceEPviolations = 1.0 * CLHEP::eV;  //***LOOKHERE***
0065   fPrimaryParticleId = 0;
0066   fPrimaryParticleInitialKineticEnergy = 0.0;
0067   fPrimaryParticleInitialTotalEnergy = 0.0;
0068   fPrimaryParticleInitialMomentum = 0.0;
0069   fPrimaryParticleInitialBeta = 1.0;
0070   fPrimaryParticleInitialGamma = 1.0;
0071   fPrimaryParticleInitial3Momentum = G4ThreeVector(0.0, 0.0, 0.0);
0072   fPrimaryParticleInitialPosition = G4ThreeVector(0.0, 0.0, 0.0);
0073   fMaxEkin_deltaMax = 0.0;
0074   fMaxEtot_deltaMax = 0.0;
0075   fMaxP_deltaMax = 0.0;
0076   fMaxPdir_deltaMax = 0.0;
0077   fMaxMass_deltaMax1 = 0.0;
0078   fMaxMass_deltaMax2 = 0.0;
0079   fMaxMass_deltaMax3 = 0.0;
0080   fMeanMass_deltaMax3 = 0.0;
0081   fMaxBeta_deltaMax1 = 0.0;
0082   fMaxBeta_deltaMax2 = 0.0;
0083   fMaxGamma_deltaMax1 = 0.0;
0084   fMaxGamma_deltaMax2 = 0.0;
0085   fMaxGamma_deltaMax3 = 0.0;
0086   fMaxT_proper_deltaMax = 0.0;
0087   fMaxT_lab_deltaMax = 0.0;
0088   fMaxMc_truth_rPos_deltaMax = 0.0;
0089   fMeanMc_truth_rPos_deltaMax = 0.0;
0090   fMeanDeltaR_primaryDecay = 0.0;
0091   fMinDeltaR_primaryDecay = 9999999.9;
0092   fMaxDeltaR_primaryDecay = -9999999.9;
0093   fMeanR_primaryDecay = 0.0;
0094   fMinR_primaryDecay = 9999999.9;
0095   fMaxR_primaryDecay = -9999999.9;
0096   fMeanX_primaryDecay = 0.0;
0097   fMinX_primaryDecay = 9999999.9;
0098   fMaxX_primaryDecay = -9999999.9;
0099   fMeanY_primaryDecay = 0.0;
0100   fMinY_primaryDecay = 9999999.9;
0101   fMaxY_primaryDecay = -9999999.9;
0102   fMeanZ_primaryDecay = 0.0;
0103   fMinZ_primaryDecay = 9999999.9;
0104   fMaxZ_primaryDecay = -9999999.9;
0105   fMeanDeltaAngle_primaryDecay = 0.0;
0106   fMinDeltaAngle_primaryDecay = 9999999.9;
0107   fMaxDeltaAngle_primaryDecay = -9999999.9;
0108   fMeanDeltaEkin_primaryDecay = 0.0;
0109   fMinDeltaEkin_primaryDecay = 9999999.9;
0110   fMaxDeltaEkin_primaryDecay = -9999999.9;
0111   fMeanEkin_primaryDecay = 0.0;
0112   fMinEkin_primaryDecay = 9999999.9;
0113   fMaxEkin_primaryDecay = -9999999.9;
0114   fMeanPx_primaryDecay = 0.0;
0115   fMinPx_primaryDecay = 9999999.9;
0116   fMaxPx_primaryDecay = -9999999.9;
0117   fMeanPy_primaryDecay = 0.0;
0118   fMinPy_primaryDecay = 9999999.9;
0119   fMaxPy_primaryDecay = -9999999.9;
0120   fMeanPz_primaryDecay = 0.0;
0121   fMinPz_primaryDecay = 9999999.9;
0122   fMaxPz_primaryDecay = -9999999.9;
0123   fMinUnderestimated_mc_truth_rPos_delta = 9999999.9;
0124   fMaxOverestimated_mc_truth_rPos_delta = -9999999.9;
0125   fMeanUnderestimated_mc_truth_rPos_delta = 0.0;
0126   fMeanOverestimated_mc_truth_rPos_delta = 0.0;
0127   fMinUnderestimated_rDeltaPos = 9999999.9;
0128   fMaxOverestimated_rDeltaPos = -9999999.9;
0129   fMeanUnderestimated_rDeltaPos = 0.0;
0130   fMeanOverestimated_rDeltaPos = 0.0;
0131   fMaxFloat_rDeltaPos_deltaMax = -9999999.9;
0132   fMeanViolationE_primaryDecay = 0.0;
0133   fMinViolationE_primaryDecay = 9999999.9;
0134   fMaxViolationE_primaryDecay = -9999999.9;
0135   fMeanViolationPx_primaryDecay = 0.0;
0136   fMinViolationPx_primaryDecay = 9999999.9;
0137   fMaxViolationPx_primaryDecay = -9999999.9;
0138   fMeanViolationPy_primaryDecay = 0.0;
0139   fMinViolationPy_primaryDecay = 9999999.9;
0140   fMaxViolationPy_primaryDecay = -9999999.9;
0141   fMeanViolationPz_primaryDecay = 0.0;
0142   fMinViolationPz_primaryDecay = 9999999.9;
0143   fMaxViolationPz_primaryDecay = -9999999.9;
0144 }
0145 
0146 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0147 
0148 void SteppingAction::UserSteppingAction(const G4Step* theStep)
0149 {
0150   // Store the information about the ID and the kinetic energy of the primary particle,
0151   // at the first step of the first event.
0152   // Note that for the kinetic energy, we are considering the "pre-step" point of such first step.
0153   if (theStep->GetTrack()->GetParentID() == 0 && theStep->GetTrack()->GetCurrentStepNumber() == 1) {
0154     fPrimaryParticleId = theStep->GetTrack()->GetDefinition()->GetPDGEncoding();
0155     fPrimaryParticleInitialKineticEnergy = theStep->GetPreStepPoint()->GetKineticEnergy();
0156     fPrimaryParticleInitialTotalEnergy = theStep->GetPreStepPoint()->GetTotalEnergy();
0157     fPrimaryParticleInitial3Momentum = theStep->GetPreStepPoint()->GetMomentum();
0158     fPrimaryParticleInitialMomentum = fPrimaryParticleInitial3Momentum.mag();
0159     fPrimaryParticleInitialPosition = theStep->GetPreStepPoint()->GetPosition();
0160     fPrimaryParticleInitialBeta = theStep->GetPreStepPoint()->GetBeta();
0161     fPrimaryParticleInitialGamma = theStep->GetPreStepPoint()->GetGamma();
0162     // As tolerance for EP violations, consider the max value between the default value
0163     // and 1 billionth of the initial, primary particle kinetic energy.
0164     if (fToleranceEPviolations < fPrimaryParticleInitialKineticEnergy * 1.0e-9) {
0165       fToleranceEPviolations = fPrimaryParticleInitialKineticEnergy * 1.0e-9;
0166     }
0167     // Set the values of this run to the Run object
0168     if (fRunPtr) {
0169       fRunPtr->SetPrimaryParticleId(fPrimaryParticleId);
0170       fRunPtr->SetPrimaryParticleInitialKineticEnergy(fPrimaryParticleInitialKineticEnergy);
0171       fRunPtr->SetPrimaryParticleInitialTotalEnergy(fPrimaryParticleInitialTotalEnergy);
0172       fRunPtr->SetPrimaryParticleInitialMomentum(fPrimaryParticleInitialMomentum);
0173       fRunPtr->SetPrimaryParticleInitialBeta(fPrimaryParticleInitialBeta);
0174       fRunPtr->SetPrimaryParticleInitialGamma(fPrimaryParticleInitialGamma);
0175       fRunPtr->SetPrimaryParticleInitial3Momentum(fPrimaryParticleInitial3Momentum);
0176       fRunPtr->SetPrimaryParticleInitialPosition(fPrimaryParticleInitialPosition);
0177       fRunPtr->SetToleranceEPviolations(ToleranceEPviolations());
0178       fRunPtr->SetToleranceDeltaDecayRadius(ToleranceDeltaDecayRadius());
0179       fRunPtr->SetIsPreassignedDecayEnabled(IsPreassignedDecayEnabled());
0180       fRunPtr->SetIsBoostToLabEnabled(IsBoostToLabEnabled());
0181     }
0182     // Use the preassigned decay is enabled
0183     if (IsPreassignedDecayEnabled() && (!theStep->GetTrack()->GetDefinition()->GetPDGStable())) {
0184       G4DynamicParticle* dynamicParent =
0185         const_cast<G4DynamicParticle*>(theStep->GetTrack()->GetDynamicParticle());
0186       if (dynamicParent != nullptr) {
0187         G4DecayProducts* decayProducts =
0188           (G4DecayProducts*)(dynamicParent->GetPreAssignedDecayProducts());
0189         if (decayProducts == nullptr) {
0190           G4ParticleDefinition* parentDef = theStep->GetTrack()->GetDefinition();
0191           G4DecayTable* decayTable = (parentDef == nullptr ? nullptr : parentDef->GetDecayTable());
0192           if (decayTable != nullptr) {
0193             G4double parentMass = dynamicParent->GetMass();
0194             G4VDecayChannel* decayChannel = decayTable->SelectADecayChannel(parentMass);
0195             if (decayChannel != nullptr) {
0196               decayProducts = decayChannel->DecayIt(parentMass);
0197               if (!decayProducts->IsChecked()) decayProducts->DumpInfo();
0198               if (IsBoostToLabEnabled()) {
0199                 // boost all decay products to laboratory frame
0200                 decayProducts->Boost(dynamicParent->GetTotalEnergy(),
0201                                      dynamicParent->GetMomentumDirection());
0202               }
0203             }
0204             else {
0205               decayProducts = new G4DecayProducts(*dynamicParent);
0206             }
0207             dynamicParent->SetPreAssignedDecayProducts(decayProducts);
0208           }
0209         }
0210         else {
0211           G4cout << "WARNING : already present preassign decay !" << G4endl;
0212         }
0213       }
0214     }
0215   }
0216 
0217   // G4cout << theStep->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName() << G4endl;
0218 
0219   // If the primary decays somewhere inside the World volume, get the information about the decay
0220   if (theStep->GetTrack()->GetParentID() == 0
0221       && theStep->GetPostStepPoint()->GetProcessDefinedStep() != nullptr
0222       && theStep->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName().find("Decay")
0223            != std::string::npos)
0224   {
0225     // Get properties of the primary particle when it decays
0226 
0227     //--- Get values in different ways and check their consistency ---
0228     // Kinetic energy of the primary particle at the decay
0229     const G4double ekin_dynamicParticle =
0230       theStep->GetTrack()->GetDynamicParticle()->GetKineticEnergy();
0231     const G4double ekin_track = theStep->GetTrack()->GetKineticEnergy();
0232     const G4double ekin_postStepPoint = theStep->GetPostStepPoint()->GetKineticEnergy();
0233     const G4double ekin_deltaMax = std::max(std::abs(ekin_dynamicParticle - ekin_track),
0234                                             std::abs(ekin_dynamicParticle - ekin_postStepPoint));
0235     // G4cout << "\t ekin_deltaMax [eV] = " << ekin_deltaMax / CLHEP::eV << G4endl;
0236     const G4double ekin_val = ekin_dynamicParticle;  // To be used later
0237     // Total energy of the primary particle at the decay
0238     const G4double etot_dynamicParticle =
0239       theStep->GetTrack()->GetDynamicParticle()->GetTotalEnergy();
0240     const G4double etot_track = theStep->GetTrack()->GetTotalEnergy();
0241     const G4double etot_postStepPoint = theStep->GetPostStepPoint()->GetTotalEnergy();
0242     const G4double etot_deltaMax = std::max(std::abs(etot_dynamicParticle - etot_track),
0243                                             std::abs(etot_dynamicParticle - etot_postStepPoint));
0244     // G4cout << "\t etot_deltaMax [eV] = " << etot_deltaMax / CLHEP::eV << G4endl;
0245     const G4double etot_val = etot_dynamicParticle;  // To be used later
0246     // Module of the 3-momentum of the primary particle at the decay
0247     const G4double p_dynamicParticle =
0248       theStep->GetTrack()->GetDynamicParticle()->GetMomentum().mag();
0249     const G4double p_track = theStep->GetTrack()->GetMomentum().mag();
0250     const G4double p_postStepPoint = theStep->GetPostStepPoint()->GetMomentum().mag();
0251     const G4double p_deltaMax = std::max(std::abs(p_dynamicParticle - p_track),
0252                                          std::abs(p_dynamicParticle - p_postStepPoint));
0253     // G4cout << "\t p_deltaMax [eV] = " << p_deltaMax / CLHEP::eV << G4endl;
0254     const G4double p_val = p_dynamicParticle;  // To be used later
0255     // 3-momentum direction (adimensional) of the primary particle at the decay
0256     const G4ThreeVector pdir_dynamicParticle =
0257       theStep->GetTrack()->GetDynamicParticle()->GetMomentumDirection();
0258     const G4ThreeVector pdir_track = theStep->GetTrack()->GetMomentumDirection();
0259     const G4ThreeVector pdir_postStepPoint = theStep->GetPostStepPoint()->GetMomentumDirection();
0260     const G4double pdir_x_deltaMax =
0261       std::max(std::abs(pdir_dynamicParticle.x() - pdir_track.x()),
0262                std::abs(pdir_dynamicParticle.x() - pdir_postStepPoint.x()));
0263     const G4double pdir_y_deltaMax =
0264       std::max(std::abs(pdir_dynamicParticle.y() - pdir_track.y()),
0265                std::abs(pdir_dynamicParticle.y() - pdir_postStepPoint.y()));
0266     const G4double pdir_z_deltaMax =
0267       std::max(std::abs(pdir_dynamicParticle.z() - pdir_track.z()),
0268                std::abs(pdir_dynamicParticle.z() - pdir_postStepPoint.z()));
0269     const G4double pdir_deltaMax =
0270       std::max(std::max(pdir_x_deltaMax, pdir_y_deltaMax), pdir_z_deltaMax);
0271     // G4cout << "\t pdir_deltaMax = " << pdir_deltaMax << G4endl;
0272     //  Mass of the primary particle at the decay
0273     const G4double mass_dynamicParticle = theStep->GetTrack()->GetDynamicParticle()->GetMass();
0274     const G4double mass_preStepPoint = theStep->GetPreStepPoint()->GetMass();
0275     const G4double mass_postStepPoint = theStep->GetPostStepPoint()->GetMass();
0276     const G4double mass_from_etot_ekin = etot_val - ekin_val;
0277     const G4double mass_from4mom = std::sqrt(etot_val * etot_val - p_val * p_val);
0278     G4double mass_deltaMax1 = std::max(std::abs(mass_dynamicParticle - mass_preStepPoint),
0279                                        std::abs(mass_dynamicParticle - mass_postStepPoint));
0280     G4double mass_deltaMax2 = std::abs(mass_dynamicParticle - mass_from_etot_ekin);
0281     G4double mass_deltaMax3 = std::abs(mass_dynamicParticle - mass_from4mom);
0282     fMeanMass_deltaMax3 += mass_deltaMax3;
0283     // G4cout << "\t mass_deltaMax{1,2,3} [eV] = " << mass_deltaMax1 / CLHEP::eV << "\t"
0284     //        << mass_deltaMax2 / CLHEP::eV << "\t" << mass_deltaMax3 / CLHEP::eV << G4endl;
0285     const G4double mass_val = mass_dynamicParticle;  // To be used later
0286     // Lorentz beta of the primary particle at the decay
0287     // The following line works only for G4 versions >= 10.7
0288     const G4double beta_dynamicParticle = theStep->GetTrack()->GetDynamicParticle()->GetBeta();
0289     const G4double beta_postStepPoint = theStep->GetPostStepPoint()->GetBeta();
0290     // Before-10.7  const G4double beta_dynamicParticle = beta_postStepPoint;
0291     const G4double beta_velocity_track = theStep->GetTrack()->GetVelocity() / CLHEP::c_light;
0292     const G4double beta_velocity_postStepPoint =
0293       theStep->GetPostStepPoint()->GetVelocity() / CLHEP::c_light;
0294     const G4double beta_p_over_etot = p_val / etot_val;
0295     G4double beta_deltaMax1 = std::max(std::abs(beta_dynamicParticle - beta_postStepPoint),
0296                                        std::abs(beta_dynamicParticle - beta_velocity_track));
0297     beta_deltaMax1 =
0298       std::max(beta_deltaMax1, std::abs(beta_dynamicParticle - beta_velocity_postStepPoint));
0299     const G4double beta_deltaMax2 = std::abs(beta_dynamicParticle - beta_p_over_etot);
0300     // G4cout << "\t beta_deltaMax{1,2} = " << beta_deltaMax1 << " , " << beta_deltaMax2 << G4endl;
0301     const G4double beta_val = beta_dynamicParticle;  // To be used later
0302     // Lorentz gamma of the primary particle at the decay
0303     const G4double gamma_postStepPoint = theStep->GetPostStepPoint()->GetGamma();
0304     const G4double gamma_from_e_over_m = etot_val / mass_val;
0305     const G4double gamma_deltaMax1 = std::abs(gamma_postStepPoint - gamma_from_e_over_m);
0306     G4double gamma_from_beta = 0.0;
0307     G4double gamma_deltaMax2 = 0.0;
0308     G4double gamma_deltaMax3 = 0.0;
0309     if (beta_val < 1.0) {
0310       gamma_from_beta = 1.0 / std::sqrt(1.0 - beta_val * beta_val);
0311       gamma_deltaMax2 = std::abs(gamma_postStepPoint - gamma_from_beta);
0312       gamma_deltaMax3 = std::abs(gamma_from_e_over_m - gamma_from_beta);
0313     }
0314     const G4double gamma_val = gamma_postStepPoint;  // To be used later;
0315     // G4cout << "\t gamma_deltaMax{1,2,3} = " << gamma_deltaMax1 << " , " << gamma_deltaMax2
0316     //        << " , " << gamma_deltaMax3  << " ; gamma_postStepPoint = " << gamma_postStepPoint
0317     //        << " ; gamma_from_e_over_m = " << gamma_from_e_over_m << " ; gamma_from_beta = "
0318     //        << gamma_from_beta << G4endl;
0319     //  Proper time of the primary particle at the decay
0320     const G4double t_proper_track = theStep->GetTrack()->GetProperTime();
0321     const G4double t_proper_postStepPoint = theStep->GetPostStepPoint()->GetProperTime();
0322     const G4double t_proper_deltaMax = std::abs(t_proper_track - t_proper_postStepPoint);
0323     // G4cout << "\t t_proper_deltaMax [fs] = " << t_proper_deltaMax / femtosecond << G4endl;
0324     const G4double t_proper_val = t_proper_track;  // To be used later
0325     // Lab time of the primary particle at the decay
0326     // (Note: it would be wrong to trying to compute this lab time from the
0327     //        above proper time via the simple formula:
0328     //          const G4double t_lab_from_gamma = t_proper_val * gamma_val;
0329     //        because the gamma value of the primary particle has changed
0330     //        during its lifetime.)
0331     const G4double t_local_track = theStep->GetTrack()->GetLocalTime();
0332     const G4double t_local_postStepPoint = theStep->GetPostStepPoint()->GetLocalTime();
0333     const G4double t_global_track = theStep->GetTrack()->GetGlobalTime();
0334     const G4double t_global_postStepPoint = theStep->GetPostStepPoint()->GetGlobalTime();
0335     G4double t_lab_deltaMax = std::max(std::abs(t_local_track - t_local_postStepPoint),
0336                                        std::abs(t_local_track - t_global_track));
0337     t_lab_deltaMax = std::max(t_lab_deltaMax, std::abs(t_local_track - t_global_postStepPoint));
0338     // G4cout << "\t t_lab_deltaMax [fs] = " << t_lab_deltaMax / femtosecond << G4endl;
0339     const G4double t_lab_val = t_local_track;  // To be used later
0340     // "MC-truth" decay radius of the primary particle at the decay
0341     // (defined as the one that would happen if there are neither magnetic field effects
0342     // nor interactions with matter).
0343     const G4double primaryBeta =
0344       fPrimaryParticleInitialMomentum / fPrimaryParticleInitialTotalEnergy;
0345     const G4double mc_truth_rPos1 = t_lab_val * fPrimaryParticleInitialBeta * CLHEP::c_light;
0346     const G4double mc_truth_rPos2 = t_lab_val * primaryBeta * CLHEP::c_light;
0347     const G4double mc_truth_rPos_deltaMax = std::abs(mc_truth_rPos1 - mc_truth_rPos2);
0348     fMeanMc_truth_rPos_deltaMax += mc_truth_rPos_deltaMax;
0349     // G4cout << "\t mc_truth_rPos_deltaMax [mum] = "
0350     //        << mc_truth_rPos_deltaMax / CLHEP::micrometer << G4endl;
0351     if (mc_truth_rPos_deltaMax > ToleranceDeltaDecayRadius()) {
0352       // G4cout << std::setprecision(6)
0353       //        << " Large : mc_truth_rPos_deltaMax [mum]="
0354       //        << mc_truth_rPos_deltaMax / CLHEP::micrometer
0355       //        << " ; " << mc_truth_rPos1 << " , " << mc_truth_rPos2 << " mm" << G4endl;
0356       if (fRunPtr) fRunPtr->IncrementNumber_mc_truth_rPos_deltaMax_above();
0357     }
0358     const G4double mc_truth_rPos_val = mc_truth_rPos1;  // To be used later
0359     // Keep note of the biggest discrepancies
0360     fMaxEkin_deltaMax = std::max(fMaxEkin_deltaMax, ekin_deltaMax);
0361     fMaxEtot_deltaMax = std::max(fMaxEtot_deltaMax, etot_deltaMax);
0362     fMaxP_deltaMax = std::max(fMaxP_deltaMax, p_deltaMax);
0363     fMaxPdir_deltaMax = std::max(fMaxPdir_deltaMax, pdir_deltaMax);
0364     fMaxMass_deltaMax1 = std::max(fMaxMass_deltaMax1, mass_deltaMax1);
0365     fMaxMass_deltaMax2 = std::max(fMaxMass_deltaMax2, mass_deltaMax2);
0366     fMaxMass_deltaMax3 = std::max(fMaxMass_deltaMax3, mass_deltaMax3);
0367     fMaxBeta_deltaMax1 = std::max(fMaxBeta_deltaMax1, beta_deltaMax1);
0368     fMaxBeta_deltaMax2 = std::max(fMaxBeta_deltaMax2, beta_deltaMax2);
0369     fMaxGamma_deltaMax1 = std::max(fMaxGamma_deltaMax1, gamma_deltaMax1);
0370     fMaxGamma_deltaMax2 = std::max(fMaxGamma_deltaMax2, gamma_deltaMax2);
0371     fMaxGamma_deltaMax3 = std::max(fMaxGamma_deltaMax3, gamma_deltaMax3);
0372     fMaxT_lab_deltaMax = std::max(fMaxT_lab_deltaMax, t_lab_deltaMax);
0373     fMaxT_proper_deltaMax = std::max(fMaxT_proper_deltaMax, t_proper_deltaMax);
0374     fMaxMc_truth_rPos_deltaMax = std::max(fMaxMc_truth_rPos_deltaMax, mc_truth_rPos_deltaMax);
0375     //--- End consistency checks ---
0376 
0377     // Global position
0378     const G4double xPos = theStep->GetPostStepPoint()->GetPosition().x();
0379     const G4double yPos = theStep->GetPostStepPoint()->GetPosition().y();
0380     const G4double zPos = theStep->GetPostStepPoint()->GetPosition().z();
0381     const G4double rPos = std::sqrt(xPos * xPos + yPos * yPos + zPos * zPos);
0382     // I have verified that for this case in which only primaries are considered, the
0383     // "GetGlobalTime()" is the same as "GetLocalTime()" (the one we use).
0384     // Moreover, this value is also the same as "GetProperTime()"*gamma .
0385     G4double tPos = theStep->GetPostStepPoint()->GetLocalTime();
0386     // The "MC-truth" decay radius is defined as the one that would happen if there are
0387     // neither magnetic field effects nor interactions with matter.
0388     const G4double mc_truth_rPos = tPos * fPrimaryParticleInitialBeta * CLHEP::c_light;
0389     const G4double rDeltaPos = mc_truth_rPos - rPos;
0390     const G4double eKin = theStep->GetPostStepPoint()->GetKineticEnergy();
0391     const G4double xMom = theStep->GetPostStepPoint()->GetMomentum().x();
0392     const G4double yMom = theStep->GetPostStepPoint()->GetMomentum().y();
0393     const G4double zMom = theStep->GetPostStepPoint()->GetMomentum().z();
0394     // The compute here the angular deflection, in degrees, between the initial direction of
0395     // the primary particle - which is along the x-axis, and its direction when it decays.
0396     G4double xDirection = std::min(theStep->GetPostStepPoint()->GetMomentumDirection().x(), 1.0);
0397     if (xDirection < -1.0) xDirection = -1.0;
0398     const G4double deflection_angle_in_degrees = 57.29 * std::acos(xDirection);
0399     const G4double delta_ekin = fPrimaryParticleInitialKineticEnergy - eKin;
0400     // G4cout << std::setprecision(6)
0401     //        << " Decay: tPos[ns]=" << tPos << " ; rPos[mm]=" << rPos << " ; deltaR[mum]="
0402     //        << rDeltaPos /CLHEP::micrometer << " ; deltaEkin[MeV]=" << delta_ekin
0403     //        << " ; deltaAngle(deg)=" << deflection_angle_in_degrees << G4endl;
0404     //  If the absolute difference between the "MC-truth" decay radius and the real one is above a
0405     //  given threshold, then we notify this special situation in the output, with "LARGE_DELTA_R"
0406     //  for post-processing evaluation. Moreover, in this case, if the "MC-truth" decay radius is
0407     //  smaller than the real one, then we count this unexpected occurrence and we further notify
0408     //  this special situation in the output with "***UNEXPECTED***" for post-processing
0409     //  evaluation.
0410     if (std::abs(rDeltaPos) > ToleranceDeltaDecayRadius()) {
0411       // G4cout << "\t LARGE_DELTA_R : mc_truth_rPos[mm]=" << mc_truth_rPos
0412       //        << " ; rPos[mm]=" << rPos;
0413       if (rDeltaPos < 0.0) {
0414         // G4cout << "\t ***UNEXPECTED***";
0415         if (fRunPtr) fRunPtr->IncrementNumberUnexpectedDecays();
0416       }
0417       // G4cout << G4endl;
0418     }
0419     fMeanDeltaR_primaryDecay += rDeltaPos;
0420     fMinDeltaR_primaryDecay = std::min(fMinDeltaR_primaryDecay, rDeltaPos);
0421     fMaxDeltaR_primaryDecay = std::max(fMaxDeltaR_primaryDecay, rDeltaPos);
0422     fMeanR_primaryDecay += rPos;
0423     fMinR_primaryDecay = std::min(fMinR_primaryDecay, rPos);
0424     fMaxR_primaryDecay = std::max(fMaxR_primaryDecay, rPos);
0425     fMeanX_primaryDecay += xPos;
0426     fMinX_primaryDecay = std::min(fMinX_primaryDecay, xPos);
0427     fMaxX_primaryDecay = std::max(fMaxX_primaryDecay, xPos);
0428     fMeanY_primaryDecay += yPos;
0429     fMinY_primaryDecay = std::min(fMinY_primaryDecay, yPos);
0430     fMaxY_primaryDecay = std::max(fMaxY_primaryDecay, yPos);
0431     fMeanZ_primaryDecay += zPos;
0432     fMinZ_primaryDecay = std::min(fMinZ_primaryDecay, zPos);
0433     fMaxZ_primaryDecay = std::max(fMaxZ_primaryDecay, zPos);
0434     fMeanDeltaAngle_primaryDecay += deflection_angle_in_degrees;
0435     fMinDeltaAngle_primaryDecay =
0436       std::min(fMinDeltaAngle_primaryDecay, deflection_angle_in_degrees);
0437     fMaxDeltaAngle_primaryDecay =
0438       std::max(fMaxDeltaAngle_primaryDecay, deflection_angle_in_degrees);
0439     fMeanDeltaEkin_primaryDecay += delta_ekin;
0440     fMinDeltaEkin_primaryDecay = std::min(fMinDeltaEkin_primaryDecay, delta_ekin);
0441     fMaxDeltaEkin_primaryDecay = std::max(fMaxDeltaEkin_primaryDecay, delta_ekin);
0442     fMeanEkin_primaryDecay += eKin;
0443     fMinEkin_primaryDecay = std::min(fMinEkin_primaryDecay, eKin);
0444     fMaxEkin_primaryDecay = std::max(fMaxEkin_primaryDecay, eKin);
0445     fMeanPx_primaryDecay += xMom;
0446     fMinPx_primaryDecay = std::min(fMinPx_primaryDecay, xMom);
0447     fMaxPx_primaryDecay = std::max(fMaxPx_primaryDecay, xMom);
0448     fMeanPy_primaryDecay += yMom;
0449     fMinPy_primaryDecay = std::min(fMinPy_primaryDecay, yMom);
0450     fMaxPy_primaryDecay = std::max(fMaxPy_primaryDecay, yMom);
0451     fMeanPz_primaryDecay += zMom;
0452     fMinPz_primaryDecay = std::min(fMinPz_primaryDecay, zMom);
0453     fMaxPz_primaryDecay = std::max(fMaxPz_primaryDecay, zMom);
0454 
0455     //--- Extra checks ---
0456     // Compute the "MC-truth" decay radius using the proper time of the primary particle when
0457     // it decays.
0458     // To do this, we would need an "effective" or "average" Lorentz beta and gamma of the primary
0459     // particle during its lifetime, whereas in practice we have only the Lorentz beta and gamma
0460     // values at the beginning and at the end when it decays. So, we can get only either an
0461     // overestimate of the "MC-truth" decay radius - by using the initial Lorentz beta and gamma -
0462     // or an underestimate of it - by using the Lorentz beta and gamma at the decay.
0463     // We want to check the average values and the largest values of these wrong estimates.
0464     const G4double underestimated_mc_truth_rPos =
0465       t_proper_val * gamma_val * beta_val * CLHEP::c_light;
0466     const G4double overestimated_mc_truth_rPos =
0467       t_proper_val * fPrimaryParticleInitialGamma * fPrimaryParticleInitialBeta * CLHEP::c_light;
0468     const G4double underestimated_mc_truth_rPos_delta =
0469       underestimated_mc_truth_rPos - mc_truth_rPos_val;
0470     const G4double overestimated_mc_truth_rPos_delta =
0471       overestimated_mc_truth_rPos - mc_truth_rPos_val;
0472     fMeanUnderestimated_mc_truth_rPos_delta += underestimated_mc_truth_rPos_delta;
0473     fMeanOverestimated_mc_truth_rPos_delta += overestimated_mc_truth_rPos_delta;
0474     // G4cout << "\t underestimated_mc_truth_rPos_delta [mum] = "
0475     //        << underestimated_mc_truth_rPos_delta / CLHEP::micrometer
0476     //        << " ; overestimated_mc_truth_rPos_delta [mum] = "
0477     //        << overestimated_mc_truth_rPos_delta / CLHEP::micrometer << G4endl;
0478     if (-underestimated_mc_truth_rPos_delta > ToleranceDeltaDecayRadius()) {
0479       // G4cout << std::setprecision(6)
0480       //        << " Large : underestimated_mc_truth_rPos_delta [mum]="
0481       //        << underestimated_mc_truth_rPos_delta / CLHEP::micrometer
0482       //        << " ; " << underestimated_mc_truth_rPos << " , "
0483       //        << mc_truth_rPos_val << " mm" << G4endl;
0484       if (fRunPtr) fRunPtr->IncrementNumber_underestimated_mc_truth_rPos_delta_above();
0485     }
0486     if (overestimated_mc_truth_rPos_delta > ToleranceDeltaDecayRadius()) {
0487       // G4cout << std::setprecision(6)
0488       //        << " Large : overestimated_mc_truth_rPos_delta [mum]="
0489       //        << overestimated_mc_truth_rPos_delta / CLHEP::micrometer
0490       //        << " ; " << overestimated_mc_truth_rPos << " , "
0491       //        << mc_truth_rPos_val << " mm" << G4endl;
0492       if (fRunPtr) fRunPtr->IncrementNumber_overestimated_mc_truth_rPos_delta_above();
0493     }
0494     const G4double underestimated_rDeltaPos = underestimated_mc_truth_rPos - rPos;
0495     const G4double overestimated_rDeltaPos = overestimated_mc_truth_rPos - rPos;
0496     fMeanUnderestimated_rDeltaPos += underestimated_rDeltaPos;
0497     fMeanOverestimated_rDeltaPos += overestimated_rDeltaPos;
0498     // G4cout << std::setprecision(6)
0499     //        << "\t underestimated_rDeltaPos=" << underestimated_rDeltaPos/CLHEP::micrometer
0500     //        << " ; overestimated_rDeltaPos=" << overestimated_rDeltaPos/CLHEP::micrometer
0501     //        << " mum" << G4endl;
0502     if (-underestimated_rDeltaPos > ToleranceDeltaDecayRadius()) {
0503       if (fRunPtr) fRunPtr->IncrementNumberLargeUnderestimates();
0504     }
0505     if (overestimated_rDeltaPos > ToleranceDeltaDecayRadius()) {
0506       if (fRunPtr) fRunPtr->IncrementNumberLargeOverestimates();
0507     }
0508     // Keep note of the biggest discrepancies
0509     fMinUnderestimated_mc_truth_rPos_delta =
0510       std::min(fMinUnderestimated_mc_truth_rPos_delta, underestimated_mc_truth_rPos_delta);
0511     fMaxOverestimated_mc_truth_rPos_delta =
0512       std::max(fMaxOverestimated_mc_truth_rPos_delta, overestimated_mc_truth_rPos_delta);
0513     fMinUnderestimated_rDeltaPos = std::min(fMinUnderestimated_rDeltaPos, underestimated_rDeltaPos);
0514     fMaxOverestimated_rDeltaPos = std::max(fMaxOverestimated_rDeltaPos, overestimated_rDeltaPos);
0515     //--- End extra checks ---
0516 
0517     // Check numerical errors due to the use of  float  instead of  double : try out
0518     // several, equivalent computations, taking the one with the largest numerical error.
0519     const G4float float_xPos = static_cast<G4float>(theStep->GetPostStepPoint()->GetPosition().x());
0520     const G4float float_yPos = static_cast<G4float>(theStep->GetPostStepPoint()->GetPosition().y());
0521     const G4float float_zPos = static_cast<G4float>(theStep->GetPostStepPoint()->GetPosition().z());
0522     const G4float float_rPos =
0523       std::sqrt(float_xPos * float_xPos + float_yPos * float_yPos + float_zPos * float_zPos);
0524     const G4float float_tPos = static_cast<G4float>(theStep->GetPostStepPoint()->GetLocalTime());
0525     const G4float float_initialBeta1 = static_cast<G4float>(fPrimaryParticleInitialBeta);
0526     const G4float float_initialBeta2 = static_cast<G4float>(fPrimaryParticleInitialMomentum)
0527                                        / static_cast<G4float>(fPrimaryParticleInitialTotalEnergy);
0528     const G4float float_initialGamma = static_cast<G4float>(fPrimaryParticleInitialGamma);
0529     const G4float float_initialBeta3 =
0530       std::sqrt(float_initialGamma * float_initialGamma - 1.0) / float_initialGamma;
0531     const G4float float_c_light = static_cast<G4float>(CLHEP::c_light);
0532     const G4float float_mc_truth_rPos1 = float_tPos * float_initialBeta1 * float_c_light;
0533     const G4float float_mc_truth_rPos2 = float_tPos * float_initialBeta2 * float_c_light;
0534     const G4float float_mc_truth_rPos3 = float_tPos * float_initialBeta3 * float_c_light;
0535     const G4float float_rDeltaPos_0 = static_cast<G4float>(rDeltaPos);
0536     const G4float float_rDeltaPos_1 = float_mc_truth_rPos1 - float_rPos;
0537     const G4float float_rDeltaPos_2 = float_mc_truth_rPos2 - float_rPos;
0538     const G4float float_rDeltaPos_3 = float_mc_truth_rPos3 - float_rPos;
0539     const G4float float_rDeltaPos_4 = static_cast<G4float>(mc_truth_rPos) - float_rPos;
0540     const G4float float_rDeltaPos_5 = float_mc_truth_rPos1 - static_cast<G4float>(rPos);
0541     const G4float float_rDeltaPos_6 = float_mc_truth_rPos2 - static_cast<G4float>(rPos);
0542     const G4float float_rDeltaPos_7 = float_mc_truth_rPos3 - static_cast<G4float>(rPos);
0543     G4double rDeltaPos_deltaMax =
0544       std::max(std::abs(float_rDeltaPos_0 - rDeltaPos), std::abs(float_rDeltaPos_1 - rDeltaPos));
0545     rDeltaPos_deltaMax = std::max(rDeltaPos_deltaMax, std::abs(float_rDeltaPos_2 - rDeltaPos));
0546     rDeltaPos_deltaMax = std::max(rDeltaPos_deltaMax, std::abs(float_rDeltaPos_3 - rDeltaPos));
0547     rDeltaPos_deltaMax = std::max(rDeltaPos_deltaMax, std::abs(float_rDeltaPos_4 - rDeltaPos));
0548     rDeltaPos_deltaMax = std::max(rDeltaPos_deltaMax, std::abs(float_rDeltaPos_5 - rDeltaPos));
0549     rDeltaPos_deltaMax = std::max(rDeltaPos_deltaMax, std::abs(float_rDeltaPos_6 - rDeltaPos));
0550     rDeltaPos_deltaMax = std::max(rDeltaPos_deltaMax, std::abs(float_rDeltaPos_7 - rDeltaPos));
0551     // G4cout << std::setprecision(6) << " rDeltaPos_deltaMax[mum]="
0552     //        << rDeltaPos_deltaMax / CLHEP::micrometer << G4endl;
0553     fMaxFloat_rDeltaPos_deltaMax = std::max(fMaxFloat_rDeltaPos_deltaMax, rDeltaPos_deltaMax);
0554 
0555     // Get properties of the decay products and check the energy-momentum conservation of the
0556     // decay
0557     std::size_t nSec = theStep->GetNumberOfSecondariesInCurrentStep();
0558     const std::vector<const G4Track*>* ptrVecSecondaries = theStep->GetSecondaryInCurrentStep();
0559     G4double deltaE = 0.0, deltaPx = 0.0, deltaPy = 0.0, deltaPz = 0.0;
0560     if (nSec > 0 && ptrVecSecondaries != nullptr) {
0561       G4double sumEsecondaries = 0.0;
0562       G4ThreeVector sumPsecondaries(0.0, 0.0, 0.0);
0563       for (std::size_t i = 0; i < nSec; ++i) {
0564         if ((*ptrVecSecondaries)[i]) {
0565           sumEsecondaries += (*ptrVecSecondaries)[i]->GetTotalEnergy();
0566           sumPsecondaries += (*ptrVecSecondaries)[i]->GetMomentum();
0567         }
0568       }
0569       deltaE = sumEsecondaries - theStep->GetPostStepPoint()->GetTotalEnergy();
0570       fMeanViolationE_primaryDecay += deltaE;
0571       fMinViolationE_primaryDecay = std::min(fMinViolationE_primaryDecay, deltaE);
0572       fMaxViolationE_primaryDecay = std::max(fMaxViolationE_primaryDecay, deltaE);
0573       if (std::abs(deltaE) > ToleranceEPviolations()) {
0574         if (fRunPtr) fRunPtr->IncrementNumberEviolations();
0575       }
0576       deltaPx = sumPsecondaries.x() - xMom;
0577       fMeanViolationPx_primaryDecay += deltaPx;
0578       fMinViolationPx_primaryDecay = std::min(fMinViolationPx_primaryDecay, deltaPx);
0579       fMaxViolationPx_primaryDecay = std::max(fMaxViolationPx_primaryDecay, deltaPx);
0580       deltaPy = sumPsecondaries.y() - yMom;
0581       fMeanViolationPy_primaryDecay += deltaPy;
0582       fMinViolationPy_primaryDecay = std::min(fMinViolationPy_primaryDecay, deltaPy);
0583       fMaxViolationPy_primaryDecay = std::max(fMaxViolationPy_primaryDecay, deltaPy);
0584       deltaPz = sumPsecondaries.z() - zMom;
0585       fMeanViolationPz_primaryDecay += deltaPz;
0586       fMinViolationPz_primaryDecay = std::min(fMinViolationPz_primaryDecay, deltaPz);
0587       fMaxViolationPz_primaryDecay = std::max(fMaxViolationPz_primaryDecay, deltaPz);
0588       if (std::abs(deltaPx) > ToleranceEPviolations() || std::abs(deltaPy) > ToleranceEPviolations()
0589           || std::abs(deltaPz) > ToleranceEPviolations())
0590       {
0591         if (fRunPtr) fRunPtr->IncrementNumberPviolations();
0592       }
0593     }
0594     else {
0595       if (fRunPtr) fRunPtr->IncrementNumberBadPrimaryDecays();
0596     }
0597 
0598     if (fRunPtr) {
0599       fRunPtr->IncrementNumberDecays();
0600       fRunPtr->SetDecayT(tPos);
0601       fRunPtr->SetDecayR_mc_truth(mc_truth_rPos);
0602       fRunPtr->SetDecayR(rPos);
0603       fRunPtr->SetDecayX(xPos);
0604       fRunPtr->SetDecayY(yPos);
0605       fRunPtr->SetDecayZ(zPos);
0606       fRunPtr->SetDeltaDecayR(rDeltaPos);
0607       fRunPtr->SetDeflectionAngle(deflection_angle_in_degrees);
0608       fRunPtr->SetDeltaEkin(delta_ekin);
0609       fRunPtr->SetDecayEkin(eKin);
0610       fRunPtr->SetDecayPx(xMom);
0611       fRunPtr->SetDecayPy(yMom);
0612       fRunPtr->SetDecayPz(zMom);
0613       fRunPtr->SetDecayEtotViolation(deltaE);
0614       fRunPtr->SetDecayPxViolation(deltaPx);
0615       fRunPtr->SetDecayPyViolation(deltaPy);
0616       fRunPtr->SetDecayPzViolation(deltaPz);
0617       fRunPtr->SetMaxEkin_deltaMax(ekin_deltaMax);
0618       fRunPtr->SetMaxEtot_deltaMax(etot_deltaMax);
0619       fRunPtr->SetMaxP_deltaMax(p_deltaMax);
0620       fRunPtr->SetMaxPdir_deltaMax(pdir_deltaMax);
0621       fRunPtr->SetMaxMass_deltaMax1(mass_deltaMax1);
0622       fRunPtr->SetMaxMass_deltaMax2(mass_deltaMax2);
0623       fRunPtr->SetMaxMass_deltaMax3(mass_deltaMax3);
0624       fRunPtr->SetMaxBeta_deltaMax1(beta_deltaMax1);
0625       fRunPtr->SetMaxBeta_deltaMax2(beta_deltaMax2);
0626       fRunPtr->SetMaxGamma_deltaMax1(gamma_deltaMax1);
0627       fRunPtr->SetMaxGamma_deltaMax2(gamma_deltaMax2);
0628       fRunPtr->SetMaxGamma_deltaMax3(gamma_deltaMax3);
0629       fRunPtr->SetMaxT_proper_deltaMax(t_proper_deltaMax);
0630       fRunPtr->SetMaxT_lab_deltaMax(t_lab_deltaMax);
0631       fRunPtr->SetMaxMc_truth_rPos_deltaMax(mc_truth_rPos_deltaMax);
0632       fRunPtr->SetMinUnderestimated_mc_truth_rPos_delta(underestimated_mc_truth_rPos_delta);
0633       fRunPtr->SetMaxOverestimated_mc_truth_rPos_delta(overestimated_mc_truth_rPos_delta);
0634       fRunPtr->SetMinUnderestimated_rDeltaPos(underestimated_rDeltaPos);
0635       fRunPtr->SetMaxOverestimated_rDeltaPos(overestimated_rDeltaPos);
0636       fRunPtr->SetMaxFloat_rDeltaPos_deltaMax(fMaxFloat_rDeltaPos_deltaMax);
0637     }
0638   }
0639 }
0640 
0641 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!