Hadr02/src/HadronicInelasticModelCRMC.cc

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