Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-04-04 08:05:12 

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 hadronic/Hadr02/src/HadronicInelasticModelCRMC.cc
0027 /// \brief Implementation of the HadronicInelasticModelCRMC class methods
0028 //
0029 //
0030 //---------------------------------------------------------------------------
0031 //
0032 #ifdef G4_USE_CRMC
0033 
0034 #  include "HadronicInelasticModelCRMC.hh"
0035 
0036 #  include "G4IonTable.hh"
0037 #  include "G4NucleiProperties.hh"
0038 #  include "G4ParticleDefinition.hh"
0039 #  include "G4ParticleTable.hh"
0040 #  include "G4SystemOfUnits.hh"
0041 #  include "G4ThreeVector.hh"
0042 #  include "Randomize.hh"
0043 
0044 #  include <cmath>
0045 #  include <cstdlib>
0046 #  include <iostream>
0047 #  include <string>
0048 
0049 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0050 
0051 #  define MAX_ENERGY_LAB_GEV 10000000.
0052 #  define MAX_ENERGY_CMS_GEV \
0053     30000.  // assuming that the target is <=100 times heavier than the projectile
0054 
0055 #  define IGNORE_PARTICLE_UNKNOWN_PDGID false
0056 #  define USE_ENERGY_CORR false
0057 #  define ENERGY_NON_CONSERVATION_RESAMPLE false
0058 #  define ENERGY_NON_CONSERVATION_EMAX_GEV 0.999
0059 #  define ENERGY_NON_CONSERVATION_FRACTION_MAX 0.00001
0060 #  define ENERGY_NON_CONSERVATION_FRACTION_MAX_ATTEMPT 10
0061 #  define ENERGY_NON_CONSERVATION_FRACTION_MAX_ENERGYTRY ENERGY_NON_CONSERVATION_EMAX_GEV
0062 
0063 #  define SPLIT_MULTI_NEUTRONS_MAXN 10
0064 #  define PARTICLE_MULTI_NEUTRONS_ERRORCODE -1
0065 
0066 #  define ERROR_REPORT_EMAIL "andrii.tykhonov@SPAMNOTcern.ch"
0067 #  define CRMC_CONFIG_FILE_ENV_VARIABLE "CRMC_CONFIG_FILE"
0068 
0069 //***********************************
0070 // CRMC ION DEFINITION
0071 // ID = CRMC_ION_COEF_0 +
0072 //      CRMC_ION_COEF_Z * Z +
0073 //      CRMC_ION_COEF_A * A
0074 //
0075 #  define CRMC_ION_COEF_0 1000000000
0076 #  define CRMC_ION_COEF_Z 10000
0077 #  define CRMC_ION_COEF_A 10
0078 //***********************************
0079 
0080 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0081 
0082 HadronicInelasticModelCRMC::HadronicInelasticModelCRMC(int model, const G4String& modelName)
0083   : G4HadronicInteraction(modelName), fPrintDebug(false)
0084 {
0085   SetMaxEnergy(MAX_ENERGY_LAB_GEV * GeV);
0086 
0087   // int model = 1; // Epos (use temporary), it is faster
0088   // int model = 12; // Dpmjet
0089   int seed = 123456789;
0090   // int seed = CLHEP::HepRandom::getTheSeed();  // Returns 0 which is invalid
0091   int produce_tables = 0;  // CRMC default, see CRMCoptions.cc in the CRMC package
0092   fTypeOutput = 0;  // CRMC default, see CRMCoptions.cc in the CRMC package
0093   static std::string crmc_param =
0094     GetCrmcParamPath();  //"crmc.param"; // CRMC default, see CRMCoptions.cc in the CRMC package
0095 
0096   fInterface = new CRMCinterface();
0097   fInterface->init(model);
0098 
0099   // open FORTRAN IO at first call
0100   fInterface->crmc_init(MAX_ENERGY_CMS_GEV, seed, model, produce_tables, fTypeOutput,
0101                         crmc_param.c_str(), "", 0);
0102 
0103   // final state
0104   finalState = new G4HadFinalState();
0105 
0106   // geant4 particle helpers (tables)
0107   fParticleTable = G4ParticleTable::GetParticleTable();
0108   fIonTable = fParticleTable->GetIonTable();
0109 }
0110 
0111 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0112 
0113 std::string HadronicInelasticModelCRMC::GetCrmcParamPath()
0114 {
0115   std::string crmcParamPath = std::getenv(CRMC_CONFIG_FILE_ENV_VARIABLE);
0116   if (crmcParamPath == "") {
0117     std::ostringstream errorstr;
0118     errorstr << "CRMC ERROR: could not find crmc param file, please check "
0119              << CRMC_CONFIG_FILE_ENV_VARIABLE << " envornoment variable!";
0120     std::string error(errorstr.str());
0121     std::cout << error << std::endl;
0122     throw error;
0123   }
0124   std::cout << "Using CRMC parameter file: " << crmcParamPath << std::endl;
0125   return crmcParamPath;
0126 }
0127 
0128 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0129 
0130 HadronicInelasticModelCRMC::~HadronicInelasticModelCRMC() {}
0131 
0132 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0133 
0134 G4HadFinalState* HadronicInelasticModelCRMC::ApplyYourself(const G4HadProjectile& aTrack,
0135                                                            G4Nucleus& targetNucleus)
0136 {
0137   //* leanup data vectors
0138   gCRMC_data.Clean();
0139 
0140   //* cleanup geant4 final state vector
0141   finalState->Clear();
0142   finalState->SetStatusChange(
0143     G4HadFinalStateStatus::stopAndKill);  // TODO: check: inelastic collisions kills previos
0144                                           // particles?
0145 
0146   //* git input particles parameters
0147   int id_proj = aTrack.GetDefinition()->GetPDGEncoding();
0148   int id_targ = targetNucleus.GetZ_asInt() * 10000 + targetNucleus.GetA_asInt() * 10;
0149   double p_proj = aTrack.Get4Momentum().pz() / GeV;
0150   double e_proj = aTrack.Get4Momentum().e() / GeV;
0151   double p_targ = 0.;
0152   double e_targ =
0153     G4NucleiProperties::GetNuclearMass(targetNucleus.GetA_asInt(), targetNucleus.GetZ_asInt())
0154     / GeV;
0155   double e_initial = e_proj + e_targ;
0156   // ... bug fix (March 2, 2020 - momentum per nucleon!)
0157   double a_proj = (double)(aTrack.GetDefinition()->GetAtomicMass());  // GetAtomicNumber());
0158   if (a_proj < 1.0)
0159     a_proj = 1.0;  // explanation: if particle is not an ion/proton, the GetAtomicMass returns 0
0160   double a_targ = (double)(targetNucleus.GetA_asInt());
0161 
0162   //* DEBUG messages
0163   if (fPrintDebug) {
0164     std::cout << "\n\n\n\n\n\n\n==============================================" << std::endl;
0165     std::cout << "Start interaction" << std::endl;
0166     std::cout << "id_proj=" << id_proj << std::endl;
0167     std::cout << "id_targ=" << id_targ << std::endl;
0168     std::cout << "p_proj=" << p_proj << std::endl;
0169     std::cout << "p_targ=" << p_targ << std::endl;
0170   }
0171 
0172   // set up input particle type and energy
0173   fInterface->crmc_set(1,  // fNCollision,
0174                        p_proj / a_proj,  // fCfg.fProjectileMomentum (per nucleon!!!),
0175                        p_targ / a_targ,  // fCfg.fTargetMomentum (per nucleon!!!),
0176                        id_proj,  // fCfg.fProjectileId,
0177                        id_targ);  // fCfg.fTargetId);
0178 
0179   //=================================================
0180   // sample 1 interaction until the energy
0181   // conservation is fulfilled
0182   int resample_attampts = 1;
0183   double max_energy_diff = ENERGY_NON_CONSERVATION_EMAX_GEV;
0184   double energy_diff_coef = 1.;
0185   double forbid_energy_corr = false;
0186   while (true) {
0187     // run one interaction
0188     fInterface->crmc_generate(fTypeOutput,  // fCfg.fTypoaut,
0189                               1,  // iColl+1,
0190                               gCRMC_data.fNParticles, gCRMC_data.fImpactParameter,
0191                               gCRMC_data.fPartId[0], gCRMC_data.fPartPx[0], gCRMC_data.fPartPy[0],
0192                               gCRMC_data.fPartPz[0], gCRMC_data.fPartEnergy[0],
0193                               gCRMC_data.fPartMass[0], gCRMC_data.fPartStatus[0]);
0194 
0195     // split Z=0 A>1 "particles" into multiple neutrons
0196     SplitMultiNeutrons(gCRMC_data);
0197 
0198     // energy check
0199     double e_final = 0;
0200     for (int i = 0; i < gCRMC_data.fNParticles; i++) {
0201       if (gCRMC_data.fPartStatus[i] != 1) continue;  // only final state particles
0202       G4ParticleDefinition* pdef;
0203       int Z_test = (gCRMC_data.fPartId[i] - CRMC_ION_COEF_0) / CRMC_ION_COEF_Z;
0204       int A_test =
0205         (gCRMC_data.fPartId[i] - CRMC_ION_COEF_0 - CRMC_ION_COEF_Z * Z_test) / CRMC_ION_COEF_A;
0206       if (fPrintDebug) {
0207         std::cout << std::endl;
0208         std::cout << "**********************************************************************"
0209                   << std::endl;
0210         std::cout << "PDG test: " << gCRMC_data.fPartId[i] << std::endl;
0211         std::cout << "fIonTable->GetIon(Z_test, A_test)                  = "
0212                   << fIonTable->GetIon(Z_test, A_test) << std::endl;
0213         std::cout << "ParticleTable->FindParticle(gCRMC_data.fPartId[i]) = "
0214                   << fParticleTable->FindParticle(gCRMC_data.fPartId[i]) << std::endl;
0215         std::cout << "**********************************************************************"
0216                   << std::endl;
0217       }
0218 
0219       // pdef = fParticleTable->FindParticle(gCRMC_data.fPartId[i]);
0220       int pdef_errorcode;
0221       pdef = GetParticleDefinition(gCRMC_data.fPartId[i], pdef_errorcode);
0222       if (!pdef && IGNORE_PARTICLE_UNKNOWN_PDGID) {
0223         continue;
0224       }
0225 
0226       double p2 = std::pow(gCRMC_data.fPartPx[i], 2) + std::pow(gCRMC_data.fPartPy[i], 2)
0227                   + std::pow(gCRMC_data.fPartPz[i], 2);
0228       double mass = pdef->GetPDGMass() / GeV;
0229       e_final += std::sqrt(mass * mass + p2);
0230     }
0231 
0232     // Check if we need to resample again...
0233     double diff = std::fabs(e_final - e_initial);
0234     if (e_final != 0. && e_initial != 0. && USE_ENERGY_CORR) energy_diff_coef = e_final / e_initial;
0235     if (fPrintDebug) {
0236       std::cout << "# e_initial = " << e_initial << " GeV" << std::endl;
0237       std::cout << "# e_final   = " << e_final << " GeV" << std::endl;
0238       std::cout << "# energy_diff_coef = " << energy_diff_coef << std::endl;
0239     }
0240 
0241     // energy conservation check, if yes
0242     if (!ENERGY_NON_CONSERVATION_RESAMPLE) {
0243       // ===== NOCHECK ========== NOCHECK ============== NOCHECK ========
0244       break;
0245       // ===== NOCHECK ========== NOCHECK ============== NOCHECK ========
0246     }
0247     else if (diff < max_energy_diff || diff / e_initial < ENERGY_NON_CONSERVATION_FRACTION_MAX) {
0248       // ===== OK ========== OK ============== OK ========
0249       forbid_energy_corr = true;
0250       break;  // everything is fine, no need to resample, break the re-sampling loop
0251       // ===== OK ========== OK ============== OK ========
0252     }
0253     else if (resample_attampts < ENERGY_NON_CONSERVATION_FRACTION_MAX_ATTEMPT) {
0254       resample_attampts++;
0255       std::cout << std::endl;
0256       std::cout
0257         << "#==== WARNING ==== WARNING ==== WARNING ==== WARNING ==== WARNING ==== WARNING ====#"
0258         << std::endl;
0259       std::cout
0260         << "#                                                                                  #"
0261         << std::endl;
0262       std::cout
0263         << "# [HadronicInelasticModelCRMC::ApplyYourself]: Energy non conservation detected: #"
0264         << std::endl;
0265       std::cout << "# e_initial = " << e_initial << " GeV" << std::endl;
0266       std::cout << "# e_final   = " << e_final << " GeV" << std::endl;
0267       std::cout << "# diff      = " << diff << " GeV" << std::endl;
0268       std::cout << "# Running attempt #" << resample_attampts << std::endl;
0269       std::cout
0270         << "#                                                                                  #"
0271         << std::endl;
0272       std::cout
0273         << "#==== WARNING ==== WARNING ==== WARNING ==== WARNING ==== WARNING ==== WARNING ====#"
0274         << std::endl;
0275       std::cout << std::endl;
0276     }
0277     else if (max_energy_diff < ENERGY_NON_CONSERVATION_FRACTION_MAX_ENERGYTRY) {
0278       std::cout << std::endl;
0279       std::cout
0280         << "#==== WARNING ==== WARNING ==== WARNING ==== WARNING ==== WARNING ==== WARNING ====#"
0281         << std::endl;
0282       std::cout << "# reached maximum number of attempts = "
0283                 << ENERGY_NON_CONSERVATION_FRACTION_MAX_ATTEMPT
0284                 << " ==> increasing twice the energy threshold!" << std::endl;
0285       max_energy_diff *= 2.;
0286       std::cout << "# max_energy_diff = " << max_energy_diff << std::endl;
0287       std::cout
0288         << "#==== WARNING ==== WARNING ==== WARNING ==== WARNING ==== WARNING ==== WARNING ====#"
0289         << std::endl;
0290       std::cout << std::endl;
0291       resample_attampts = 1;
0292     }
0293     else {
0294       std::cout << std::endl;
0295       std::cout
0296         << "#==== WARNING ==== WARNING ==== WARNING ==== WARNING ==== WARNING ==== WARNING ====#"
0297         << std::endl;
0298       std::cout << "# reached maximum number of attempts = "
0299                 << ENERGY_NON_CONSERVATION_FRACTION_MAX_ATTEMPT << "not resampling any more!"
0300                 << std::endl;
0301       std::cout
0302         << "#==== WARNING ==== WARNING ==== WARNING ==== WARNING ==== WARNING ==== WARNING ====#"
0303         << std::endl;
0304       std::cout << std::endl;
0305       // ===== FAIL ========== FAIL ============== FAIl ========
0306       break;
0307       // ===== FAIL ========== FAIL ============== FAIl ========
0308     }
0309   }
0310   // ... finished sampling one interaction
0311   //=================================================
0312 
0313   // ... for DEBUG messages
0314   double totalenergy = 0;
0315   double totalz = 0;
0316   double eaftertest0 = 0.;
0317   double eaftertest1 = 0.;
0318   double eaftertest2 = 0.;
0319 
0320   // save secondary particles for outputa
0321   for (int i = 0; i < gCRMC_data.fNParticles; i++) {
0322     //* Keep only final state particles
0323     // .. (-9 is the beam, 2 is a particle which decayed and 1 is final)
0324     if (gCRMC_data.fPartStatus[i] != 1) continue;
0325 
0326     if (fPrintDebug) {
0327       std::cout << "\n\nSecondary:" << std::endl
0328                 << gCRMC_data.fPartId[i] << std::endl
0329                 << gCRMC_data.fPartPx[i] << std::endl
0330                 << gCRMC_data.fPartPy[i] << std::endl
0331                 << gCRMC_data.fPartPz[i] << std::endl
0332                 << gCRMC_data.fPartEnergy[i] << "  ENERGY  " << std::endl;
0333     }
0334 
0335     // G4ParticleDefinition* pdef = fParticleTable->FindParticle(gCRMC_data.fPartId[i]);
0336     int pdef_errorcode;
0337     G4ParticleDefinition* pdef = GetParticleDefinition(gCRMC_data.fPartId[i], pdef_errorcode);
0338     if (!pdef) {
0339       if (IGNORE_PARTICLE_UNKNOWN_PDGID) {
0340         std::cout << std::endl;
0341         std::cout << "*****************************************************************************"
0342                      "***************************"
0343                   << std::endl;
0344         std::cout << " -- WARNING "
0345                      "-----------------------------------------------------------------------------"
0346                      "------ WARNING --  "
0347                   << std::endl;
0348         std::cout
0349           << " [HadronicInelasticModelCRMC] Geant4 could not find particle definition for PDG ID = "
0350           << gCRMC_data.fPartId[i] << std::endl;
0351         std::cout << " [HadronicInelasticModelCRMC] Ignoring this particle. This might cause "
0352                      "energy non-conservation!"
0353                   << std::endl;
0354         std::cout << " -- WARNING "
0355                      "-----------------------------------------------------------------------------"
0356                      "------ WARNING --  "
0357                   << std::endl;
0358         std::cout << "*****************************************************************************"
0359                      "***************************"
0360                   << std::endl;
0361         continue;
0362       }
0363       else {
0364         std::cout << std::endl;
0365         std::cout << "*****************************************************************************"
0366                      "***************************"
0367                   << std::endl;
0368         std::cout << " -- ERROR "
0369                      "-----------------------------------------------------------------------------"
0370                      "------ ERROR --  "
0371                   << std::endl;
0372         std::cout
0373           << " [HadronicInelasticModelCRMC] Geant4 could not find particle definition for PDG ID = "
0374           << gCRMC_data.fPartId[i] << std::endl;
0375         std::cout << " [HadronicInelasticModelCRMC] Throwing exception! Please report to: "
0376                   << ERROR_REPORT_EMAIL << std::endl;
0377         std::cout << " -- ERROR "
0378                      "-----------------------------------------------------------------------------"
0379                      "------ ERROR --  "
0380                   << std::endl;
0381         std::cout << "*****************************************************************************"
0382                      "***************************"
0383                   << std::endl;
0384         throw;
0385       }
0386     }
0387 
0388     double part_e_corr = 1.;
0389     double part_p_corr = 1.;
0390     if (USE_ENERGY_CORR && !forbid_energy_corr && energy_diff_coef != 0) {
0391       part_e_corr = 1. / energy_diff_coef;
0392       double pbefore2 = std::pow(gCRMC_data.fPartPx[i], 2) + std::pow(gCRMC_data.fPartPy[i], 2)
0393                         + std::pow(gCRMC_data.fPartPz[i], 2);
0394       double mass2 =
0395         std::pow(pdef->GetPDGMass() / GeV, 2);  // std::pow(gCRMC_data.fPartEnergy[i],2) - pbefore2;
0396       double ebefore2 = pbefore2 + mass2;
0397       double pafter2 = ebefore2 * part_e_corr * part_e_corr - mass2;
0398       if (pbefore2) part_p_corr = std::sqrt(std::fabs(pafter2 / pbefore2));
0399       if (fPrintDebug) std::cout << "part_p_corr=" << part_p_corr << std::endl;
0400       eaftertest0 += std::sqrt(mass2 + pbefore2);
0401       eaftertest1 += std::sqrt(mass2 + pafter2);
0402     }
0403 
0404     G4DynamicParticle* part =
0405       new G4DynamicParticle(pdef, G4ThreeVector(gCRMC_data.fPartPx[i] * GeV * part_p_corr,
0406                                                 gCRMC_data.fPartPy[i] * GeV * part_p_corr,
0407                                                 gCRMC_data.fPartPz[i] * GeV * part_p_corr));
0408     eaftertest2 += part->GetTotalEnergy();
0409     finalState->AddSecondary(part);
0410     totalenergy += gCRMC_data.fPartEnergy[i];
0411     totalz += gCRMC_data.fPartPz[i];
0412   }
0413 
0414   if (fPrintDebug) {
0415     std::cout << "totalenergy (GeV) = " << totalenergy << std::endl;
0416     std::cout << "totalz (GeV)      = " << totalz << std::endl;
0417     std::cout << "initialz (GeV)    = " << p_proj + p_targ << std::endl;
0418     std::cout << "eaftertest0       = " << eaftertest0 << std::endl;
0419     std::cout << "eaftertest1       = " << eaftertest1 << std::endl;
0420     std::cout << "eaftertest2       = " << eaftertest2 << std::endl;
0421     std::cout << "Finishing interaction: " << std::endl;
0422     const G4LorentzVector& p1 = aTrack.Get4Momentum();
0423     std::cout << "e=" << p1.e() << " px=" << p1.px() << " py=" << p1.py() << " pz=" << p1.pz()
0424               << std::endl;
0425     std::cout << aTrack.GetDefinition()->GetAtomicNumber() << std::endl;
0426     std::cout << aTrack.GetDefinition()->GetPDGCharge() << std::endl;
0427     std::cout << targetNucleus.GetA_asInt() << std::endl;
0428     std::cout << targetNucleus.GetZ_asInt() << std::endl;
0429     std::cout << "Stop interaction" << std::endl;
0430     std::cout << "==============================================\n\n\n\n\n\n" << std::endl;
0431   }
0432   // std::cout<<"finalState->GetNumberOfSecondaries()="<<finalState->GetNumberOfSecondaries()<<
0433   // std::endl; // Debugging info
0434   return finalState;
0435 }
0436 
0437 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0438 
0439 G4bool HadronicInelasticModelCRMC::IsApplicable(const G4HadProjectile&, G4Nucleus&)
0440 {
0441   return true;
0442 }
0443 
0444 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0445 
0446 G4ParticleDefinition* HadronicInelasticModelCRMC::GetParticleDefinition(long particle_id,
0447                                                                         int& error_code)
0448 {
0449   G4ParticleDefinition* pdef = fParticleTable->FindParticle(particle_id);
0450   if (!pdef && particle_id > CRMC_ION_COEF_0) {
0451     int Z = (particle_id - CRMC_ION_COEF_0) / CRMC_ION_COEF_Z;
0452     int A = (particle_id - CRMC_ION_COEF_0 - CRMC_ION_COEF_Z * Z) / CRMC_ION_COEF_A;
0453     if (IsMultiNeutron(Z, A)) {
0454       error_code = PARTICLE_MULTI_NEUTRONS_ERRORCODE;
0455       pdef = NULL;
0456     }
0457     else {
0458       pdef = fIonTable->GetIon(Z, A);
0459     }
0460   }
0461   return pdef;
0462 }
0463 
0464 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0465 
0466 bool HadronicInelasticModelCRMC::IsMultiNeutron(int Z, int A)
0467 {
0468   bool result = false;
0469   if (!Z && A > 1) {
0470     if (A <= SPLIT_MULTI_NEUTRONS_MAXN) {
0471       result = true;
0472     }
0473     else {
0474       std::cout << " [HadronicInelasticModelCRMC::IsMultiNeutron] ERROR A=" << A
0475                 << " is higher than " << SPLIT_MULTI_NEUTRONS_MAXN << " throwing exception!"
0476                 << std::endl;
0477       throw;
0478     }
0479   }
0480   return result;
0481 }
0482 
0483 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0484 
0485 void HadronicInelasticModelCRMC::SplitMultiNeutrons(CRMCdata& CRMC_data)
0486 {
0487   for (int i = 0; i < CRMC_data.fNParticles; i++) {
0488     // check if it is a final-state secondary particle
0489     if (CRMC_data.fPartStatus[i] != 1) continue;
0490 
0491     int pdef_errorcode;
0492     GetParticleDefinition(CRMC_data.fPartId[i], pdef_errorcode);
0493     if (pdef_errorcode != PARTICLE_MULTI_NEUTRONS_ERRORCODE) continue;
0494 
0495     //
0496     int particle_id = gCRMC_data.fPartId[i];
0497     int Z = (particle_id - CRMC_ION_COEF_0) / CRMC_ION_COEF_Z;
0498     int A = (particle_id - CRMC_ION_COEF_0 - CRMC_ION_COEF_Z * Z) / CRMC_ION_COEF_A;
0499     if (Z != 0 || A < 2) {
0500       std::cout << " [HadronicInelasticModelCRMC::SplitMultiNeutrons] ERROR consistency check "
0501                    "failed! Throwing exception! "
0502                 << std::endl;
0503       throw;
0504     }
0505 
0506     //
0507     std::cout << std::endl;
0508     std::cout << " [HadronicInelasticModelCRMC::SplitMultiNeutrons] INFO splitting the floowing "
0509                  "particle into neutrons: "
0510               << std::endl;
0511     std::cout << " [HadronicInelasticModelCRMC::SplitMultiNeutrons] INFO    Z = " << Z << std::endl;
0512     std::cout << " [HadronicInelasticModelCRMC::SplitMultiNeutrons] INFO    A = " << A << std::endl;
0513     std::cout << " [HadronicInelasticModelCRMC::SplitMultiNeutrons] INFO    fPartId = "
0514               << CRMC_data.fPartId[i] << std::endl;
0515     std::cout << " [HadronicInelasticModelCRMC::SplitMultiNeutrons] INFO    fPartPx = "
0516               << CRMC_data.fPartPx[i] << std::endl;
0517     std::cout << " [HadronicInelasticModelCRMC::SplitMultiNeutrons] INFO    fPartPy = "
0518               << CRMC_data.fPartPy[i] << std::endl;
0519     std::cout << " [HadronicInelasticModelCRMC::SplitMultiNeutrons] INFO    fPartPz = "
0520               << CRMC_data.fPartPz[i] << std::endl;
0521     std::cout << " [HadronicInelasticModelCRMC::SplitMultiNeutrons] INFO    fPartEnergy = "
0522               << CRMC_data.fPartEnergy[i] << std::endl;
0523     std::cout << " [HadronicInelasticModelCRMC::SplitMultiNeutrons] INFO    fPartMass   = "
0524               << CRMC_data.fPartMass[i] << std::endl;
0525     std::cout << " [HadronicInelasticModelCRMC::SplitMultiNeutrons] INFO    fPartStatus = "
0526               << CRMC_data.fPartStatus[i] << std::endl;
0527 
0528     //
0529     int NEUTRON_PDG_ID = 2112;
0530     G4ParticleDefinition* p_n_def = fParticleTable->FindParticle(NEUTRON_PDG_ID);
0531     double m_n = p_n_def->GetPDGMass() / GeV;
0532     double e_n = CRMC_data.fPartEnergy[i] / A;
0533     int status_n = CRMC_data.fPartStatus[i];
0534     if (e_n < m_n) {
0535       std::cout << " [HadronicInelasticModelCRMC::SplitMultiNeutrons] WARNING neutron energy "
0536                 << e_n << " lower than neutron mass " << m_n << " assigning e_n = m_n     "
0537                 << std::endl;
0538       e_n = m_n;
0539     }
0540     double p_tot_before = std::sqrt(CRMC_data.fPartPx[i] * CRMC_data.fPartPx[i]
0541                                     + CRMC_data.fPartPy[i] * CRMC_data.fPartPy[i]
0542                                     + CRMC_data.fPartPz[i] * CRMC_data.fPartPz[i]);
0543     double p_tot_after = std::sqrt(e_n * e_n - m_n * m_n);
0544     double px_n = 0;
0545     double py_n = 0;
0546     double pz_n = 0;
0547     if (p_tot_before > 0. && p_tot_after > 0.) {
0548       px_n = CRMC_data.fPartPx[i] * p_tot_after / p_tot_before;
0549       py_n = CRMC_data.fPartPy[i] * p_tot_after / p_tot_before;
0550       pz_n = CRMC_data.fPartPz[i] * p_tot_after / p_tot_before;
0551     }
0552     for (int j = 0; j < A; j++) {
0553       int i_neutron = j ? CRMC_data.fNParticles + j : i;
0554       CRMC_data.fPartId[i_neutron] = NEUTRON_PDG_ID;
0555       CRMC_data.fPartPx[i_neutron] = px_n;
0556       CRMC_data.fPartPy[i_neutron] = py_n;
0557       CRMC_data.fPartPz[i_neutron] = pz_n;
0558       CRMC_data.fPartEnergy[i_neutron] = e_n;
0559       CRMC_data.fPartMass[i_neutron] = m_n;
0560       CRMC_data.fPartStatus[i_neutron] = status_n;
0561     }
0562     CRMC_data.fNParticles += A - 1;
0563 
0564     //
0565     std::cout << " [HadronicInelasticModelCRMC::SplitMultiNeutrons] done for a particle. "
0566               << std::endl;
0567     std::cout << std::endl;
0568   }
0569 }
0570 
0571 #endif  // G4_USE_CRMC
[ 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!