Changeset 1340 for trunk/source/processes/hadronic/models/incl
- Timestamp:
- Nov 5, 2010, 3:45:55 PM (14 years ago)
- Location:
- trunk/source/processes/hadronic/models/incl
- Files:
-
- 1 added
- 21 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/source/processes/hadronic/models/incl/GNUmakefile
r1315 r1340 1 # $Id: GNUmakefile,v 1. 9 2010/04/27 16:02:37kaitanie Exp $1 # $Id: GNUmakefile,v 1.11 2010/10/26 02:47:58 kaitanie Exp $ 2 2 # ----------------------------------------------------------- 3 3 # GNUmakefile for hadronic library. Gabriele Cosmo, 18/9/96. … … 8 8 ifndef G4INSTALL 9 9 G4INSTALL = ../../../../.. 10 endif 11 12 ifdef G4INCLDEBUG 13 CPPFLAGS += -DG4INCLDEBUG 10 14 endif 11 15 … … 33 37 -I$(G4BASE)/processes/hadronic/models/util/include \ 34 38 -I$(G4BASE)/processes/hadronic/models/de_excitation/evaporation/include \ 39 -I$(G4BASE)/processes/hadronic/models/de_excitation/fermi_breakup/include \ 35 40 -I$(G4BASE)/particles/management/include \ 36 41 -I$(G4BASE)/particles/leptons/include \ -
trunk/source/processes/hadronic/models/incl/History
r1337 r1340 4 4 Geant4 - an Object-Oriented Toolkit for Physics Simulation 5 5 ========================================================== 6 $Id: History,v 1. 28 2010/06/14 16:10:01 gcosmoExp $6 $Id: History,v 1.32 2010/10/29 06:54:23 gunter Exp $ 7 7 --------------------------------------------------------------------- 8 8 … … 16 16 * Please list in reverse chronological order (last date on top) 17 17 --------------------------------------------------------------- 18 19 29 October 2010 - Gunter Folger (hadr-incl-V09-03-05) 20 - Fixed several more compilation warnings for conversion of double to int 21 in G4Incl.cc, G4InclAblaCascadeInterface.cc, and G4InclAblaLightIonInterface.cc. 22 23 28 October 2010 - Gabriele Cosmo (hadr-incl-V09-03-04) 24 ----------------------------------------------------- 25 - Fixed compilation warning in G4InclDataDefs.hh for conversion of double 26 to integer in initialization. 27 28 25 October 2010 - Pekka Kaitaniemi (hadr-incl-V09-03-03) 29 -------------------------------------------------------- 30 - Refactored INCL input data structure to its own class G4InclInput 31 o Made data members private 32 o Access to input data only through accessor methods 33 o INCL now uses integer A and Z internally as well as in the interface 34 35 14 September 2010 - Pekka Kaitaniemi (hadr-incl-V09-03-02) 36 ---------------------------------------------------------- 37 - Migrated to integer A and Z when using G4Nucleus 38 - Bugfixes: 39 o Make sure we use Geant4 random number generators everywhere 40 o Silenced some unnecessary warnings 41 o Make sure we use INCL for all light ion projectiles up to Carbon-12 42 43 16 August 2010 - Pekka Kaitaniemi 44 --------------------------------- 45 - Moved the INCL input initialization to the INCL input class 46 - Implemented inverse-kinematics treatment for ion-hydrogen collisions 47 48 01 June 2010 - Pekka Kaitaniemi 49 ------------------------------- 50 - Improved the Fermi break-up configuration 51 52 26 February 2010 - Pekka Kaitaniemi 53 ----------------------------------- 54 - Merged INCL light ion projectile support to the integration branch 55 o Support infrastructure for light ions up to Carbon 56 o Use Geant4 Fermi break-up for projectile spectators and light cascade 57 remnants 18 58 19 59 14 June 2010 - Gabriele Cosmo (hadr-incl-V09-03-01) … … 98 138 ----------------------------------------------------- 99 139 - use GetNuclearMass() instead of GetAtomicMass() in G4AblaEvaporation.cc 140 141 10 October 2008 - Pekka Kaitaniemi 142 ------------------------------------ 143 - Added ability to use Geant4 Fermi break-up 144 for light cascade remnants 100 145 101 146 12 September 2008 - Pekka Kaitaniemi (hadr-incl-V09-01-02) -
trunk/source/processes/hadronic/models/incl/include/G4Abla.hh
r962 r1340 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4Abla.hh,v 1.1 1 2008/06/25 17:20:03kaitanie Exp $26 // $Id: G4Abla.hh,v 1.13 2010/10/26 02:47:59 kaitanie Exp $ 27 27 // Translation of INCL4.2/ABLA V3 28 28 // Pekka Kaitaniemi, HIP (translation) … … 33 33 #include "globals.hh" 34 34 35 #include "G4VInclLogger.hh" 35 36 #include "G4InclRandomNumbers.hh" 36 37 #include "G4AblaDataDefs.hh" … … 75 76 76 77 /** 78 * Register the INCL/ABLA internal variable logger. 79 */ 80 void registerLogger(G4VInclLogger *theLogger); 81 82 /** 77 83 * Set verbosity level. 78 84 */ 79 void setVerboseLevel(G4int level) { 80 verboseLevel = level; 81 } 85 void setVerboseLevel(G4int level); 82 86 83 87 /** … … 102 106 * @param eventnumber number of the event 103 107 */ 104 void breakItUp(G4 double nucleusA, G4doublenucleusZ, G4double nucleusMass, G4double excitationEnergy,108 void breakItUp(G4int nucleusA, G4int nucleusZ, G4double nucleusMass, G4double excitationEnergy, 105 109 G4double angularMomentum, G4double recoilEnergy, G4double momX, G4double momY, G4double momZ, 106 110 G4int eventnumber); … … 322 326 G4Volant *volant; 323 327 G4VarNtp *varntp; 328 329 G4VInclLogger *theLogger; 324 330 }; 325 331 -
trunk/source/processes/hadronic/models/incl/include/G4AblaDataDefs.hh
r962 r1340 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4AblaDataDefs.hh,v 1. 9 2008/06/25 17:20:04kaitanie Exp $26 // $Id: G4AblaDataDefs.hh,v 1.11 2010/10/26 02:47:59 kaitanie Exp $ 27 27 // Translation of INCL4.2/ABLA V3 28 28 // Pekka Kaitaniemi, HIP (translation) … … 189 189 190 190 //#define VOLANTSIZE 200 191 #define VOLANTSIZE 2000191 #define VOLANTSIZE 301 192 192 /** 193 193 * Evaporation and fission output data. -
trunk/source/processes/hadronic/models/incl/include/G4AblaEvaporation.hh
r819 r1340 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4AblaEvaporation.hh,v 1. 1 2007/09/11 13:18:42 miheikkiExp $26 // $Id: G4AblaEvaporation.hh,v 1.3 2010/10/26 02:47:59 kaitanie Exp $ 27 27 // Defines an interface to evaporation models of Bertini cascase (BERT) 28 28 // based on INUCL code. … … 35 35 #include "G4Fragment.hh" 36 36 #include "G4DynamicParticle.hh" 37 38 37 #include "G4Abla.hh" 39 38 … … 96 95 G4double CoulombBarrier; 97 96 97 /** 98 * ABLA evaporation 99 */ 100 G4Abla *abla; 101 G4VarNtp *varntp; 102 98 103 #ifdef DEBUG 99 104 -
trunk/source/processes/hadronic/models/incl/include/G4AblaFissionBase.hh
r1337 r1340 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4AblaFissionBase.hh,v 1. 3 2010/06/14 16:10:01 gcosmoExp $26 // $Id: G4AblaFissionBase.hh,v 1.5 2010/10/26 02:47:59 kaitanie Exp $ 27 27 // Translation of INCL4.2/ABLA V3 28 28 // Pekka Kaitaniemi, HIP (translation) … … 35 35 36 36 #include "globals.hh" 37 #include "G4InclUtils.hh" 37 38 38 39 /* … … 50 51 G4double &A2, G4double &Z2, G4double &E2, G4double &K2) = 0; 51 52 53 void setVerboseLevel(G4int level) { 54 verboseLevel = level; 55 if(verboseLevel > G4InclUtils::silent) { 56 about(); 57 G4cout <<";; Fission model verbosity level set to " << verboseLevel << G4endl; 58 } 59 } 60 52 61 void about() { 53 G4cout << aboutModel << G4endl;62 G4cout << ";; " << aboutModel << G4endl; 54 63 } 55 64 … … 59 68 60 69 private: 70 G4int verboseLevel; 61 71 G4String aboutModel; 62 72 }; -
trunk/source/processes/hadronic/models/incl/include/G4Incl.hh
r1315 r1340 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4Incl.hh,v 1.1 6 2010/04/27 16:02:37kaitanie Exp $26 // $Id: G4Incl.hh,v 1.18 2010/10/26 02:47:59 kaitanie Exp $ 27 27 // Translation of INCL4.2/ABLA V3 28 28 // Pekka Kaitaniemi, HIP (translation) … … 39 39 #include "G4Abla.hh" 40 40 #include <fstream> 41 #include "G4VInclLogger.hh" 42 #include "G4InclInput.hh" 41 43 42 44 using namespace std; … … 63 65 * @param varntp a pointer to G4VarNtp structure. 64 66 */ 65 G4Incl(G4Hazard *hazard, G4 Calincl*calincl, G4Ws *ws, G4Mat *mat, G4VarNtp *varntp);67 G4Incl(G4Hazard *hazard, G4InclInput *calincl, G4Ws *ws, G4Mat *mat, G4VarNtp *varntp); 66 68 67 69 /** … … 85 87 void dumpBl2(std::ofstream& dumpOut); // Dump the contents of G4Bl2. 86 88 void dumpBl3(std::ofstream& dumpOut); // Dump the contents of G4Bl3. 89 90 /** 91 * Set Fermi break-up use flag. 92 */ 93 void setUseFermiBreakUp(G4bool useIt); 94 95 /** 96 * Set projectile spectator use flag. 97 * 98 * Select whether or not to produce so-called spectator nucleus from 99 * the projectile nucleons that do not hit the target. 100 */ 101 void setUseProjectileSpectators(G4bool useIt); 87 102 88 103 /** … … 101 116 void setSpl2Data(G4Spl2 *newSpl2); 102 117 void setMatData(G4Mat *newMat); 103 void set CalinclData(G4Calincl*newCalincl);118 void setInput(G4InclInput *newCalincl); 104 119 void setLightNucData(G4LightNuc *newLightNuc); 105 120 void setLightGausNucData(G4LightGausNuc *newLightGausNuc); … … 115 130 void setKindData(G4Kind *newKind); 116 131 132 /** 133 * Register the logger. 134 */ 135 void registerLogger(G4VInclLogger *aLogger) { 136 if(aLogger != 0) { 137 theLogger = aLogger; 138 abla->registerLogger(aLogger); 139 } 140 } 141 117 142 public: 118 143 /** 119 144 * Process one event with INCL4 only. 120 145 */ 121 void processEventIncl( );146 void processEventIncl(G4InclInput *input); 122 147 123 148 /** 124 149 * Process one event with INCL4 and built-in ABLA evaporation and fission. 125 150 */ 126 void processEventInclAbla(G4 int eventnumber);151 void processEventInclAbla(G4InclInput *input, G4int eventnumber); 127 152 128 153 public: // Methods used to initialize INCL … … 280 305 void pnu(G4int *ibert_p, G4int *nopart_p, G4int *izrem_p, G4int *iarem_p, G4double *esrem_p, 281 306 G4double *erecrem_p, G4double *alrem_p, G4double *berem_p, G4double *garem_p, 282 G4double *bimpact_p, G4int *l_p); 307 G4double *bimpact_p, G4int *l_p, G4double *xjrem, G4double *yjrem, G4double *zjrem); 308 309 //C projection de JREM sur Z: 310 // MREM=ANINT(ZJREM/197.328) 311 // IF (MREM.GT.JREM) MREM=JREM 312 // IF (MREM.LT.-JREM) MREM=-JREM 283 313 284 314 /** … … 460 490 * @return a double value 461 491 */ 462 G4double transmissionProb(G4double E, G4double iz, G4double izn, G4double r, G4double v0); 492 // G4double transmissionProb(G4double E, G4double iz, G4double izn, G4double r, G4double v0); 493 // G4double transmissionProb(G4double E, G4int iz, G4int izn, G4double r, G4double v0) 494 G4double transmissionProb(G4double E, G4int iz, G4int ia, G4int izn, G4double r, G4double v0); 495 // G4double transmissionProb(G4double E, G4int iz, G4int ia, G4int izn,G4double R, G4double v0); 496 // G4double transmissionProb(G4double E, G4double iz, G4double izn, G4double r, G4double v0); 497 498 void projo_spec(G4int ia1, G4int ips, 499 G4double fmpinc, G4double pinc, G4double tlab); 500 501 void ordered(G4double t, G4int nb); 463 502 464 503 /** … … 536 575 * @param A mass number (double parameter) 537 576 */ 538 G4double radius(G4 doubleA);577 G4double radius(G4int A); 539 578 540 579 /** Parametrisation de la section efficace de réaction calculée par incl4.1 … … 686 725 687 726 /** 688 * G4Calincl689 */ 690 G4 Calincl*calincl;727 * INCL input data structure 728 */ 729 G4InclInput *calincl; 691 730 692 731 /** … … 744 783 */ 745 784 G4Kind *kindstruct; 785 786 /** 787 * Projectile properties. 788 */ 789 G4Bev *bev; 746 790 747 791 /** … … 820 864 G4int densFunction; 821 865 866 /** 867 * Use fermi break-up? 868 */ 869 G4bool useFermiBreakup; 870 871 /** 872 * Projectile spectator nucleus support? 873 */ 874 G4bool useProjSpect; 875 876 /** 877 * Type of the particle at index i. 878 * 879 * The extension to large composite projectiles the value is 880 * negative for projectile spectators. Otherwise the particle types 881 * are given in exactly the same way as in G4Calincl::f[6]. 882 * @see G4Calincl::f 883 */ 822 884 G4int kind[300]; //= (*kind_p); 823 885 G4double ep[300]; // = (*ep_p); … … 826 888 G4double gam[300]; // = (*gam_p); 827 889 890 G4VBe *be; 891 G4InclProjSpect *ps; 892 G4InclFermi *fermi; 893 G4QuadvectProjo *qvp; 828 894 G4Volant *volant; 829 895 G4Abla *abla; 830 896 G4InclRandomInterface *randomGenerator; 831 897 898 G4VInclLogger *theLogger; 832 899 G4int inside_step; // Flag to determine whether we are inside or outside a simulation step 833 900 }; -
trunk/source/processes/hadronic/models/incl/include/G4InclAblaCascadeInterface.hh
r1196 r1340 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4InclAblaCascadeInterface.hh,v 1. 8 2009/11/18 10:43:14kaitanie Exp $26 // $Id: G4InclAblaCascadeInterface.hh,v 1.9 2010/10/26 02:47:59 kaitanie Exp $ 27 27 // Translation of INCL4.2/ABLA V3 28 28 // Pekka Kaitaniemi, HIP (translation) … … 129 129 G4Hazard *hazard; // The random seeds used by INCL. 130 130 G4VarNtp *varntp; 131 G4 Calincl*calincl;131 G4InclInput *calincl; 132 132 G4Ws *ws; 133 133 G4Mat *mat; -
trunk/source/processes/hadronic/models/incl/include/G4InclAblaLightIonInterface.hh
r819 r1340 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4InclAblaLightIonInterface.hh,v 1. 5 2007/10/31 10:44:22 miheikkiExp $26 // $Id: G4InclAblaLightIonInterface.hh,v 1.7 2010/10/26 02:47:59 kaitanie Exp $ 27 27 // Translation of INCL4.2/ABLA V3 28 28 // Pekka Kaitaniemi, HIP (translation) … … 46 46 #define G4INCLABLALIGHTIONINTERFACE_H 1 47 47 48 #include "G4Nucleon.hh"49 #include "G4Nucleus.hh"50 #include "G4HadronicInteraction.hh"51 48 #include "G4VIntraNuclearTransportModel.hh" 52 49 #include "G4KineticTrackVector.hh" … … 108 105 G4Hazard *hazard; // The random seeds used by INCL. 109 106 G4VarNtp *varntp; 110 G4 Calincl*calincl;107 G4InclInput *calincl; 111 108 G4Ws *ws; 112 109 G4Mat *mat; … … 119 116 G4double previousTargetA; 120 117 G4double previousTargetZ; 118 G4bool useProjectileSpectator; 119 G4bool useFermiBreakup; 121 120 }; 122 121 -
trunk/source/processes/hadronic/models/incl/include/G4InclCascadeInterface.hh
r819 r1340 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4InclCascadeInterface.hh,v 1. 5 2007/10/31 10:44:22 miheikkiExp $26 // $Id: G4InclCascadeInterface.hh,v 1.6 2010/10/26 02:47:59 kaitanie Exp $ 27 27 // Translation of INCL4.2/ABLA V3 28 28 // Pekka Kaitaniemi, HIP (translation) … … 59 59 #include "G4AblaDataDefs.hh" 60 60 #include "G4Incl.hh" 61 #include "G4InclInput.hh" 61 62 62 63 #include <fstream> … … 107 108 G4Hazard *hazard; // The random seeds used by INCL. 108 109 G4VarNtp *varntp; 109 G4 Calincl *calincl;110 G4InclInput *inclInput; 110 111 G4Ws *ws; 111 112 G4Mat *mat; -
trunk/source/processes/hadronic/models/incl/include/G4InclDataDefs.hh
r1196 r1340 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4InclDataDefs.hh,v 1. 7 2009/11/18 10:43:14 kaitanieExp $26 // $Id: G4InclDataDefs.hh,v 1.10 2010/10/28 15:35:50 gcosmo Exp $ 27 27 // Translation of INCL4.2/ABLA V3 28 28 // Pekka Kaitaniemi, HIP (translation) … … 36 36 #define InclDataDefs_hh 1 37 37 38 #include "G4Nucleus.hh" 39 #include "G4HadProjectile.hh" 40 #include "G4ParticleTable.hh" 41 #include "G4Track.hh" 42 43 class G4InclFermi { 44 public: 45 G4InclFermi() { 46 G4double hc = 197.328; 47 G4double fmp = 938.2796; 48 pf=1.37*hc; 49 pf2=pf*pf; 50 tf=std::sqrt(pf*pf+fmp*fmp)-fmp; 51 }; 52 ~G4InclFermi() {}; 53 54 G4double tf,pf,pf2; 55 }; 56 57 #define max_a_proj 61 58 59 /** 60 * (eps_c,p1_s,p2_s,p3_s,eps_c used to store the kinematics of 61 * nucleons for composit projectiles before entering the potential) 62 */ 63 class G4QuadvectProjo { 64 public: 65 G4QuadvectProjo() {}; 66 ~G4QuadvectProjo() {}; 67 68 G4double eps_c[max_a_proj],p3_c[max_a_proj], 69 p1_s[max_a_proj],p2_s[max_a_proj],p3_s[max_a_proj], 70 t_c[max_a_proj]; 71 }; 72 73 class G4VBe { 74 public: 75 G4VBe() {}; 76 ~G4VBe() {}; 77 78 G4int ia_be, iz_be; 79 G4double rms_be, pms_be, bind_be; 80 }; 81 82 /** 83 * Projectile spectator 84 */ 85 class G4InclProjSpect { 86 public: 87 G4InclProjSpect() { 88 // G4cout <<"Projectile spectator data structure created!" << G4endl; 89 }; 90 ~G4InclProjSpect() {}; 91 92 void clear() { 93 for(G4int i = 0; i < 21; i++) tab[i] = 0.0; 94 for(G4int i = 0; i < 61; i++) n_projspec[i] = 0; 95 a_projspec = 0; 96 z_projspec = 0; 97 m_projspec = 0.0; 98 ex_projspec = 0.0; 99 p1_projspec = 0.0; 100 p2_projspec = 0.0; 101 p3_projspec = 0.0; 102 }; 103 104 G4double tab[21]; 105 G4int n_projspec[61]; 106 G4int a_projspec,z_projspec; 107 G4double ex_projspec,t_projspec, p1_projspec, p2_projspec, p3_projspec, m_projspec; 108 }; 109 38 110 #define FSIZE 15 39 111 /** … … 42 114 class G4Calincl { 43 115 public: 44 G4Calincl() {}; 116 G4Calincl() { 117 isExtendedProjectile = false; 118 }; 119 120 G4Calincl(const G4HadProjectile &aTrack, const G4Nucleus &theNucleus, G4bool inverseKinematics = false) { 121 for(int i = 0; i < 15; i++) { 122 f[i] = 0.0; // Initialize INCL input data 123 } 124 125 usingInverseKinematics = inverseKinematics; 126 127 f[0] = theNucleus.GetA_asInt(); // Target mass number 128 f[1] = theNucleus.GetZ_asInt(); // Target charge number 129 f[6] = getBulletType(aTrack.GetDefinition()); // Projectile type (INCL particle code) 130 f[2] = aTrack.GetKineticEnergy() / MeV; // Projectile energy (total, in MeV) 131 f[5] = 1.0; // Time scaling 132 f[4] = 45.0; // Nuclear potential 133 setExtendedProjectileInfo(aTrack.GetDefinition()); 134 // G4cout <<"Projectile type = " << f[6] << G4endl; 135 // G4cout <<"Energy = " << f[2] << G4endl; 136 // G4cout <<"Target A = " << f[0] << " Z = " << f[1] << G4endl; 137 }; 138 45 139 ~G4Calincl() {}; 46 140 141 static void printProjectileTargetInfo(const G4HadProjectile &aTrack, const G4Nucleus &theNucleus) { 142 G4cout <<"Projectile = " << aTrack.GetDefinition()->GetParticleName() << G4endl; 143 G4cout <<" four-momentum: " << aTrack.Get4Momentum() << G4endl; 144 G4cout <<"Energy = " << aTrack.GetKineticEnergy() / MeV << G4endl; 145 G4cout <<"Target A = " << theNucleus.GetA_asInt() << " Z = " << theNucleus.GetZ_asInt() << G4endl; 146 } 147 148 static G4bool canUseInverseKinematics(const G4HadProjectile &aTrack, const G4Nucleus &theNucleus) { 149 G4int targetA = theNucleus.GetA_asInt(); 150 const G4ParticleDefinition *projectileDef = aTrack.GetDefinition(); 151 G4int projectileA = projectileDef->GetAtomicMass(); 152 // G4int projectileZ = projectileDef->GetAtomicNumber(); 153 if(targetA > 0 && targetA < 18 && (projectileDef != G4Proton::Proton() && 154 projectileDef != G4Neutron::Neutron() && 155 projectileDef != G4PionPlus::PionPlus() && 156 projectileDef != G4PionZero::PionZero() && 157 projectileDef != G4PionMinus::PionMinus()) && 158 projectileA > 1) { 159 return true; 160 } else { 161 return false; 162 } 163 } 164 165 G4double bulletE() { 166 return f[2]; 167 } 168 169 G4int getBulletType() { 170 return G4int(f[6]); 171 }; 172 173 void setExtendedProjectileInfo(const G4ParticleDefinition *pd) { 174 if(getBulletType(pd) == -666) { 175 extendedProjectileA = pd->GetAtomicMass(); 176 extendedProjectileZ = pd->GetAtomicNumber(); 177 isExtendedProjectile = true; 178 } else { 179 isExtendedProjectile = false; 180 } 181 }; 182 183 G4int getBulletType(const G4ParticleDefinition *pd) { 184 G4ParticleTable *pt = G4ParticleTable::GetParticleTable(); 185 186 if(pd == G4Proton::Proton()) { 187 return 1; 188 } else if(pd == G4Neutron::Neutron()) { 189 return 2; 190 } else if(pd == G4PionPlus::PionPlus()) { 191 return 3; 192 } else if(pd == G4PionMinus::PionMinus()) { 193 return 5; 194 } else if(pd == G4PionZero::PionZero()) { 195 return 4; 196 } else if(pd == G4Deuteron::Deuteron()) { 197 return 6; 198 } else if(pd == G4Triton::Triton()) { 199 return 7; 200 } else if(pd == G4He3::He3()) { 201 return 8; 202 } else if(pd == G4Alpha::Alpha()) { 203 return 9; 204 } else if(pd == pt->GetIon(6, 12, 0.0)) { // C12 special case. This should be phased-out in favor of "extended projectile" 205 return -12; 206 } else { // Is this extended projectile? 207 G4int A = pd->GetAtomicMass(); 208 G4int Z = pd->GetAtomicNumber(); 209 if(A > 4 && A <= 16 && Z > 2 && Z <= 8) { // Ions from Lithium to Oxygen 210 return -666; // Code of an extended projectile 211 } 212 } 213 G4cout <<"Error! Projectile " << pd->GetParticleName() << " not defined!" << G4endl; 214 return 0; 215 }; 216 217 G4bool isInverseKinematics() { 218 return usingInverseKinematics; 219 }; 220 221 G4double targetA() { return f[0]; }; 222 G4double targetZ() { return f[1]; }; 223 224 G4int extendedProjectileA; 225 G4int extendedProjectileZ; 226 G4bool isExtendedProjectile; 227 47 228 /** 48 229 * Here f is an array containing the following initial values: … … 69 250 */ 70 251 G4int icoup; 252 253 G4bool usingInverseKinematics; 71 254 }; 72 255 … … 206 389 class G4Ws { 207 390 public: 208 G4Ws() {}; 391 G4Ws() { 392 fneck = 0.0; 393 }; 209 394 ~G4Ws() {}; 210 395 … … 256 441 */ 257 442 G4double bmax; 443 444 G4double fneck; 258 445 }; 259 446 … … 311 498 G4int ind1[BL1SIZE],ind2[BL1SIZE]; 312 499 G4double ta; 500 501 void dump(G4int numberOfParticles) { 502 static G4int dumpNumber = 0; 503 G4cout <<"Dump number" << dumpNumber << " of particle 4-momenta (G4Bl1):" << G4endl; 504 G4cout <<"ta = " << ta << G4endl; 505 for(G4int i = 0; i < numberOfParticles; i++) { 506 G4cout <<"i = " << i << " p1 = " << p1[i] << " p2 = " << p2[i] << " p3 = " << p3[i] << " eps = " << eps[i] << G4endl; 507 } 508 dumpNumber++; 509 } 313 510 }; 314 511 … … 383 580 */ 384 581 G4double rab2; 582 583 void dump() { 584 static G4int dumpNumber = 0; 585 G4cout <<"Dump number" << dumpNumber << " of particle positions (G4Bl3):" << G4endl; 586 G4cout <<" ia1 = " << ia1 << G4endl; 587 G4cout <<" ia2 = " << ia2 << G4endl; 588 G4cout <<" rab2 = " << rab2 << G4endl; 589 for(G4int i = 0; i <= (ia1 + ia2); i++) { 590 G4cout <<"i = " << i << " x1 = " << x1[i] << " x2 = " << x2[i] << " x3 = " << x3[i] << G4endl; 591 } 592 dumpNumber++; 593 } 385 594 }; 386 595 … … 510 719 }; 511 720 721 /** 722 * Projectile parameters. 723 */ 724 class G4Bev { 725 public: 726 /** 727 * Initialize all variables to zero. 728 */ 729 G4Bev() { 730 ia_be = 0; 731 iz_be = 0; 732 rms_be = 0.0; 733 pms_be = 0.0; 734 bind_be = 0.0; 735 }; 736 ~G4Bev() {}; 737 738 /** 739 * Mass number. 740 */ 741 G4int ia_be; 742 743 /** 744 * Charge number. 745 */ 746 G4int iz_be; 747 748 /** 749 * rms 750 */ 751 G4double rms_be; 752 753 /** 754 * pms 755 */ 756 G4double pms_be; 757 758 /** 759 * bind 760 */ 761 G4double bind_be; 762 }; 763 512 764 #define VARSIZE 3 513 765 #define VAEPSSIZE 250 … … 567 819 }; 568 820 569 #define VARNTPSIZE 255821 #define VARNTPSIZE 301 570 822 class G4VarNtp { 571 823 public: … … 582 834 targetA = 0; 583 835 targetZ = 0; 836 masp = 0.0; mzsp = 0.0; exsp = 0.0; 837 // To be deleted? 838 spectatorA = 0; 839 spectatorZ = 0; 840 spectatorEx = 0.0; 841 spectatorM = 0.0; 842 spectatorT = 0.0; 843 spectatorP1 = 0.0; 844 spectatorP2 = 0.0; 845 spectatorP3 = 0.0; 584 846 massini = 0; 585 847 mzini = 0; … … 590 852 pyrem = 0; 591 853 pzrem = 0; 854 erecrem = 0; 592 855 mulncasc = 0; 593 856 mulnevap = 0; … … 600 863 iafis = 0; 601 864 ntrack = 0; 865 needsFermiBreakup = false; 602 866 for(G4int i = 0; i < VARNTPSIZE; i++) { 603 867 itypcasc[i] = 0; … … 612 876 } 613 877 878 /** 879 * Add a particle to the INCL/ABLA final output. 880 */ 614 881 void addParticle(G4double A, G4double Z, G4double E, G4double P, G4double theta, G4double phi) { 615 882 if(full[particleIndex]) { … … 752 1019 753 1020 /** 1021 * Projectile spectator A, Z, Eex; 1022 */ 1023 G4double masp, mzsp, exsp, mrem; 1024 1025 /** 1026 * Spectator nucleus mass number for light ion projectile support. 1027 */ 1028 G4int spectatorA; 1029 1030 /** 1031 * Spectator nucleus charge number for light ion projectile support. 1032 */ 1033 G4int spectatorZ; 1034 1035 /** 1036 * Spectator nucleus excitation energy for light ion projectile support. 1037 */ 1038 G4double spectatorEx; 1039 1040 /** 1041 * Spectator nucleus mass. 1042 */ 1043 G4double spectatorM; 1044 1045 /** 1046 * Spectator nucleus kinetic energy. 1047 */ 1048 G4double spectatorT; 1049 1050 /** 1051 * Spectator nucleus momentum x-component. 1052 */ 1053 G4double spectatorP1; 1054 1055 /** 1056 * Spectator nucleus momentum y-component. 1057 */ 1058 G4double spectatorP2; 1059 1060 /** 1061 * Spectator nucleus momentum z-component. 1062 */ 1063 G4double spectatorP3; 1064 1065 /** 754 1066 * A of the remnant. 755 1067 */ … … 766 1078 G4double exini; 767 1079 768 G4double pcorem, mcorem, pxrem, pyrem, pzrem ;1080 G4double pcorem, mcorem, pxrem, pyrem, pzrem, erecrem; 769 1081 770 1082 /** … … 826 1138 827 1139 /** 1140 * Does this nucleus require Fermi break-up treatment? Only 1141 * applicable when used together with Geant4. 1142 * true = do fermi break-up (and skip ABLA part) 1143 * false = use ABLA 1144 */ 1145 G4bool needsFermiBreakup; 1146 1147 /** 828 1148 * emitted in cascade (0) or evaporation (1). 829 1149 */ -
trunk/source/processes/hadronic/models/incl/include/G4InclLightIonInterface.hh
r819 r1340 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4InclLightIonInterface.hh,v 1. 5 2007/10/31 10:44:22 miheikkiExp $26 // $Id: G4InclLightIonInterface.hh,v 1.6 2010/10/26 02:47:59 kaitanie Exp $ 27 27 // Translation of INCL4.2/ABLA V3 28 28 // Pekka Kaitaniemi, HIP (translation) … … 108 108 G4Hazard *hazard; // The random seeds used by INCL. 109 109 G4VarNtp *varntp; 110 G4 Calincl*calincl;110 G4InclInput *calincl; 111 111 G4Ws *ws; 112 112 G4Mat *mat; -
trunk/source/processes/hadronic/models/incl/include/G4InclRandomNumbers.hh
r1337 r1340 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4InclRandomNumbers.hh,v 1. 4 2010/06/14 16:10:01 gcosmoExp $26 // $Id: G4InclRandomNumbers.hh,v 1.6 2010/10/26 02:47:59 kaitanie Exp $ 27 27 // Translation of INCL4.2/ABLA V3 28 28 // Pekka Kaitaniemi, HIP (translation) … … 54 54 */ 55 55 virtual G4double getRandom() = 0; 56 56 virtual void printSeeds() = 0; 57 57 private: 58 58 G4long seed; … … 77 77 return 0.5; 78 78 } 79 80 void printSeeds() {}; 79 81 }; 80 82 … … 96 98 return G4UniformRand(); 97 99 } 100 101 void printSeeds() { 102 G4cout <<"Using Geant4 random number generator." << G4endl; 103 }; 98 104 }; 99 105 -
trunk/source/processes/hadronic/models/incl/include/G4Ranecu.hh
r967 r1340 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4Ranecu.hh,v 1. 3 2008/06/25 17:20:04kaitanie Exp $26 // $Id: G4Ranecu.hh,v 1.5 2010/10/26 02:47:59 kaitanie Exp $ 27 27 // Translation of INCL4.2/ABLA V3 28 28 // Pekka Kaitaniemi, HIP (translation) … … 39 39 40 40 G4double getRandom(); 41 void printSeeds() { 42 G4cout <<"Seed1 = " << iseed1 << G4endl; 43 G4cout <<"Seed2 = " << iseed2 << G4endl; 44 }; 41 45 42 46 private: -
trunk/source/processes/hadronic/models/incl/src/G4Abla.cc
r1196 r1340 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4Abla.cc,v 1.2 0 2009/11/18 10:43:14kaitanie Exp $26 // $Id: G4Abla.cc,v 1.22 2010/10/26 02:47:59 kaitanie Exp $ 27 27 // Translation of INCL4.2/ABLA V3 28 28 // Pekka Kaitaniemi, HIP (translation) … … 98 98 } 99 99 100 void G4Abla::setVerboseLevel(G4int level) 101 { 102 verboseLevel = level; 103 if(verboseLevel > G4InclUtils::silent) { 104 G4cout <<";; G4Abla: Setting verbose level to " << verboseLevel << G4endl; 105 } 106 fissionModel->setVerboseLevel(verboseLevel); 107 } 108 100 109 G4Abla::~G4Abla() 101 110 { … … 111 120 } 112 121 122 void G4Abla::registerLogger(G4VInclLogger *theLogger) { 123 if(theLogger != NULL) { 124 this->theLogger = theLogger; 125 } 126 } 127 113 128 // Main interface to the evaporation 114 129 … … 118 133 // G4Fragment? 119 134 120 void G4Abla::breakItUp(G4 double nucleusA, G4doublenucleusZ, G4double nucleusMass, G4double excitationEnergy,135 void G4Abla::breakItUp(G4int nucleusA, G4int nucleusZ, G4double nucleusMass, G4double excitationEnergy, 121 136 G4double angularMomentum, G4double recoilEnergy, G4double momX, G4double momY, G4double momZ, 122 137 G4int eventnumber) … … 187 202 G4double esrem = excitationEnergy; 188 203 189 G4double aprf = nucleusA;190 G4double zprf = nucleusZ;204 G4double aprf = (double) nucleusA; 205 G4double zprf = (double) nucleusZ; 191 206 G4double mcorem = nucleusMass; 192 207 G4double ee = excitationEnergy; … … 208 223 209 224 G4double pcorem = std::sqrt(erecrem*(erecrem +2.*938.2796*nucleusA)); 225 #ifdef G4INCLDEBUG 226 theLogger->fillHistogram1D("pcorem", pcorem); 227 #endif 210 228 // G4double pcorem = std::sqrt(std::pow(momX,2) + std::pow(momY,2) + std::pow(momZ,2)); 211 229 if(pcorem != 0) { // Guard against division by zero. -
trunk/source/processes/hadronic/models/incl/src/G4AblaEvaporation.cc
r962 r1340 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4AblaEvaporation.cc,v 1. 4 2008/10/24 21:07:40 dennisExp $26 // $Id: G4AblaEvaporation.cc,v 1.6 2010/10/26 02:47:59 kaitanie Exp $ 27 27 // 28 28 #include <numeric> … … 82 82 hazard->igraine[17] = 33759; 83 83 hazard->igraine[18] = 13227; 84 85 G4VarNtp *evaporationResult = new G4VarNtp(); 86 G4Volant *volant = new G4Volant(); 87 88 // Initialize evaporation. 89 abla = new G4Abla(hazard, volant, evaporationResult); 90 abla->initEvapora(); 84 91 } 85 92 … … 108 115 } 109 116 110 G4FragmentVector * G4AblaEvaporation::BreakItUp(const G4Fragment &theNucleus) { 111 112 113 G4VarNtp *varntp = new G4VarNtp(); 114 G4Volant *volant = new G4Volant(); 115 116 G4Abla *abla = new G4Abla(hazard, volant, varntp); 117 G4cout <<"Initializing evaporation..." << G4endl; 118 abla->initEvapora(); 119 G4cout <<"Initialization complete!" << G4endl; 120 121 G4double nucleusA = theNucleus.GetA(); 122 G4double nucleusZ = theNucleus.GetZ(); 117 G4FragmentVector * G4AblaEvaporation::BreakItUp(const G4Fragment &theNucleus) { 118 G4int nucleusA = theNucleus.GetA_asInt(); 119 G4int nucleusZ = theNucleus.GetZ_asInt(); 123 120 G4double nucleusMass = G4NucleiProperties::GetNuclearMass(nucleusA, nucleusZ); 124 121 G4double excitationEnergy = theNucleus.GetExcitationEnergy(); … … 137 134 138 135 varntp->ntrack = -1; 139 varntp->massini = theNucleus.GetA ();140 varntp->mzini = theNucleus.GetZ ();136 varntp->massini = theNucleus.GetA_asInt(); 137 varntp->mzini = theNucleus.GetZ_asInt(); 141 138 142 139 std::vector<G4DynamicParticle*> cascadeParticles; -
trunk/source/processes/hadronic/models/incl/src/G4Incl.cc
r1315 r1340 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4Incl.cc,v 1.3 0 2010/04/27 16:02:37 kaitanieExp $26 // $Id: G4Incl.cc,v 1.33 2010/10/29 06:48:43 gunter Exp $ 27 27 // Translation of INCL4.2/ABLA V3 28 28 // Pekka Kaitaniemi, HIP (translation) … … 49 49 densFunction = 5; 50 50 51 // Default: no Fermi break-up 52 useFermiBreakup = false; 53 54 // Default: no projectile spectators 55 useProjSpect = false; 56 51 57 randomGenerator = new G4InclGeant4Random(); 52 // randomGenerator = new G4Ranecu();58 // randomGenerator = new G4Ranecu(); 53 59 } 54 60 … … 56 62 { 57 63 verboseLevel = 0; 64 65 // Default: no Fermi break-up 66 useFermiBreakup = false; 67 68 // Default: no projectile spectators 69 useProjSpect = false; 58 70 59 71 // Set functions to be used for integration routine. … … 71 83 ws = aWs; 72 84 73 // randomGenerator = new G4Ranecu();85 // randomGenerator = new G4Ranecu(); 74 86 randomGenerator = new G4InclGeant4Random(); 75 87 } 76 88 77 G4Incl::G4Incl(G4Hazard *aHazard, G4 Calincl*aCalincl, G4Ws *aWs, G4Mat *aMat, G4VarNtp *aVarntp)89 G4Incl::G4Incl(G4Hazard *aHazard, G4InclInput *aCalincl, G4Ws *aWs, G4Mat *aMat, G4VarNtp *aVarntp) 78 90 { 79 91 verboseLevel = 0; 92 93 // Default: no Fermi break-up 94 useFermiBreakup = false; 95 96 // Default: no projectile spectators 97 useProjSpect = false; 80 98 81 99 // Set functions to be used for integration routine. … … 95 113 96 114 randomGenerator = new G4InclGeant4Random(); 97 // randomGenerator = new G4Ranecu();115 // randomGenerator = new G4Ranecu(); 98 116 light_gaus_nuc = new G4LightGausNuc(); 99 117 light_nuc = new G4LightNuc(); … … 111 129 bl10 = new G4Bl10(); 112 130 kindstruct = new G4Kind(); 131 bev = new G4Bev(); 113 132 paul = new G4Paul(); 114 133 varavat = new G4VarAvat(); 115 134 varavat->kveux = 0; 135 136 be = new G4VBe(); 137 ps = new G4InclProjSpect(); 138 fermi = new G4InclFermi(); 139 qvp = new G4QuadvectProjo(); 116 140 117 141 volant = new G4Volant(); … … 123 147 abla = new G4Abla(hazard, volant, evaporationResult); 124 148 abla->initEvapora(); 149 150 theLogger = 0; 125 151 } 126 152 … … 143 169 delete bl10; 144 170 delete kindstruct; 171 delete bev; 145 172 delete paul; 146 173 delete varavat; … … 149 176 delete evaporationResult; 150 177 delete volant; 178 } 179 180 void G4Incl::setUseFermiBreakUp(G4bool useIt) 181 { 182 useFermiBreakup = useIt; 183 } 184 185 void G4Incl::setUseProjectileSpectators(G4bool useIt) 186 { 187 useProjSpect = useIt; 151 188 } 152 189 … … 236 273 { 237 274 verboseLevel = level; 275 if(verboseLevel > G4InclUtils::silent) { 276 G4cout <<";; G4Incl: Setting verbose level to " << verboseLevel << G4endl; 277 } 278 abla->setVerboseLevel(level); 238 279 } 239 280 … … 268 309 } 269 310 270 void G4Incl::set CalinclData(G4Calincl*newCalincl)311 void G4Incl::setInput(G4InclInput *newCalincl) 271 312 { 272 313 calincl = newCalincl; … … 342 383 */ 343 384 344 void G4Incl::processEventIncl() 345 { 385 void G4Incl::processEventIncl(G4InclInput *input) 386 { 387 if(input == 0) { 388 G4cerr <<"G4Incl fatal error: NULL pointer passed as input!" << G4endl; 389 return; 390 } 391 calincl = input; 392 346 393 const G4double uma = 931.4942; 347 394 const G4double melec = 0.511; … … 355 402 356 403 varntp->clear(); 357 if(calincl-> f[6] == 3.0) { // pi+404 if(calincl->bulletType() == 3) { // pi+ 358 405 mprojo = 139.56995; 359 406 ap = 0.0; … … 361 408 } 362 409 363 if(calincl-> f[6] == 4.0) { // pi0410 if(calincl->bulletType() == 4) { // pi0 364 411 mprojo = 134.9764; 365 412 ap = 0.0; … … 367 414 } 368 415 369 if(calincl-> f[6] == 5.0) { // pi-416 if(calincl->bulletType() == 5) { // pi- 370 417 mprojo = 139.56995; 371 418 ap = 0.0; … … 374 421 375 422 // Coulomb en entree seulement pour les particules ci-dessous. 376 if(calincl-> f[6] == 1.0) { // proton423 if(calincl->bulletType() == 1) { // proton 377 424 mprojo = 938.27231; 378 425 ap = 1.0; … … 380 427 } 381 428 382 if(calincl-> f[6] == 2.0) { // neutron429 if(calincl->bulletType() == 2) { // neutron 383 430 mprojo = 939.56563; 384 431 ap = 1.0; … … 386 433 } 387 434 388 if(calincl-> f[6] == 6.0) { // deuteron435 if(calincl->bulletType() == 6) { // deuteron 389 436 mprojo = 1875.61276; 390 437 ap = 2.0; … … 392 439 } 393 440 394 if(calincl-> f[6] == 7.0) { // triton441 if(calincl->bulletType() == 7) { // triton 395 442 mprojo = 2808.95; 396 443 ap = 3.0; … … 398 445 } 399 446 400 if(calincl-> f[6] == 8.0) { // He3447 if(calincl->bulletType() == 8) { // He3 401 448 mprojo = 2808.42; 402 449 ap = 3.0; … … 404 451 } 405 452 406 if(calincl-> f[6] == 9.0) { // Alpha453 if(calincl->bulletType() == 9) { // Alpha 407 454 mprojo = 3727.42; 408 455 ap = 4.0; … … 410 457 } 411 458 412 pbeam = std::sqrt(calincl->f[2]*(calincl->f[2] + 2.0*mprojo)); 413 414 G4double at = calincl->f[0]; 459 if(calincl->bulletType() == -12) { // Carbon 460 mprojo = 12.0 * 938.27231; 461 ap = 12.0; 462 zp = 6.0; 463 } else if(calincl->bulletType() == -666) { // Generic extended ion projectile 464 mprojo = calincl->extendedProjectileA() * 938.27231; 465 ap = calincl->extendedProjectileA(); 466 zp = calincl->extendedProjectileZ(); 467 } 468 469 pbeam = std::sqrt(calincl->bulletE()*(calincl->bulletE() + 2.0*mprojo)); 470 471 G4double at = calincl->targetA(); 415 472 416 calincl->f[3] = 0.0; // seuil sortie proton417 calincl->f[7] = 0.0; // seuil sortie neutron418 419 473 G4int ibert = 1; 420 474 … … 435 489 G4double probaTrans = 0.0; 436 490 G4double rndm = 0.0; 437 if((calincl->f[6] == 1.0) || (calincl->f[6] >= 6.0)) { 438 probaTrans = coulombTransm(calincl->f[2],ap,zp,calincl->f[0],calincl->f[1]); 491 492 if((calincl->bulletType() == 1) || (calincl->bulletType() >= 6)) { 493 probaTrans = coulombTransm(calincl->bulletE(),ap,zp,calincl->targetA(),calincl->targetZ()); 439 494 standardRandom(&rndm, &(hazard->ial)); 440 495 if(rndm <= (1.0 - probaTrans)) { … … 447 502 * Call the actual INCL routine. 448 503 */ 449 pnu(&ibert, &nopart,&izrem,&iarem,&esrem,&erecrem,&alrem,&berem,&garem,&bimpac,&jrem); 504 505 // randomGenerator->printSeeds(); 506 G4double jremx, jremy, jremz; 507 pnu(&ibert, &nopart,&izrem,&iarem,&esrem,&erecrem,&alrem,&berem,&garem,&bimpac,&jrem, 508 &jremx, &jremy, &jremz); 450 509 forceAbsor(&nopart, &iarem, &izrem, &esrem, &erecrem, &alrem, &berem, &garem, &jrem); 451 510 G4double aprf = double(iarem); // mass number of the prefragment … … 480 539 * -following the approximations of the cugnon code (esrem...) 481 540 */ 482 G4double f0 = calincl-> f[0];483 G4double f1 = calincl-> f[1];484 G4double f2 = calincl-> f[2];541 G4double f0 = calincl->targetA(); 542 G4double f1 = calincl->targetZ(); 543 G4double f2 = calincl->bulletE(); 485 544 G4double mcorem = mprojo + f2 + abla->pace2(f0, f1) + f0 * uma - f1 * melec; 486 545 … … 680 739 681 740 682 void G4Incl::processEventInclAbla(G4int eventnumber) 683 { 741 void G4Incl::processEventInclAbla(G4InclInput *input, G4int eventnumber) 742 { 743 //G4cout <<"Starting event " << eventnumber << G4endl; 744 if(input == 0) { 745 G4cerr <<"G4Incl fatal error: NULL pointer passed as input!" << G4endl; 746 return; 747 } 748 calincl = input; 749 684 750 const G4double uma = 931.4942; 685 751 const G4double melec = 0.511; 752 const G4double fmp = 938.2796; 686 753 687 754 G4double pcorem = 0.0; … … 694 761 varntp->clear(); 695 762 763 if(calincl->bulletType() == -12) { 764 be->ia_be = std::abs(calincl->bulletType()); 765 be->iz_be = 6; 766 } else if(calincl->bulletType() == -666) { 767 be->iz_be = calincl->extendedProjectileZ(); 768 be->ia_be = calincl->extendedProjectileA(); 769 } 770 771 if(calincl->isExtendedProjectile() == false && calincl->bulletType() < -max_a_proj) { 772 // if(calincl->bulletType() < -max_a_proj) { 773 G4cout <<"max a of composite projectile is: " << max_a_proj << G4endl; 774 exit(0); 775 } 776 if(calincl->bulletType() < 0) { 777 // calincl->bulletType() = std::floor(calincl->bulletType() + 0.1); WTF??? 778 be->pms_be=100.; 779 G4int i_tabled=0; 780 if(be->iz_be == 3 && be->ia_be == 6) { 781 be->rms_be=2.56; 782 be->bind_be=32.0; 783 i_tabled=1; 784 } else if(be->iz_be == 3 && be->ia_be == 7) { // TODO: Check the values! 785 be->rms_be=2.56; 786 be->bind_be=32.0; 787 i_tabled=1; 788 } else if(be->iz_be == 3 && be->ia_be == 8) { 789 be->rms_be=2.40; 790 be->bind_be=39.25; 791 i_tabled=1; 792 } else if(be->iz_be == 4 && be->ia_be == 7) { 793 be->rms_be=2.51; 794 be->bind_be=58.17; 795 i_tabled=1; 796 } else if(be->iz_be == 4 && be->ia_be == 9) { 797 be->rms_be=2.51; 798 be->bind_be=58.17; 799 i_tabled=1; 800 } else if(be->iz_be == 4 && be->ia_be == 10) { 801 be->rms_be=2.45; 802 be->bind_be=64.75; 803 i_tabled=1; 804 } else if(be->iz_be == 5 && be->ia_be == 10) { 805 be->rms_be=2.45; 806 be->bind_be=64.75; 807 i_tabled=1; 808 } else if(be->iz_be == 5 && be->ia_be == 11) { 809 be->rms_be=2.40; 810 be->bind_be=76.21; 811 i_tabled=1; 812 } else if(be->iz_be == 6 && be->ia_be == 9) { // TODO: Check the values! 813 be->rms_be=2.44; 814 be->bind_be=92.17; 815 i_tabled=1; 816 } else if(be->iz_be == 6 && be->ia_be == 10) { // TODO: Check the values! 817 be->rms_be=2.44; 818 be->bind_be=92.17; 819 i_tabled=1; 820 } else if(be->iz_be == 6 && be->ia_be == 11) { // TODO: Check the values! 821 be->rms_be=2.44; 822 be->bind_be=92.17; 823 i_tabled=1; 824 } else if(be->iz_be == 6 && calincl->bulletType() == -12) { // Special Carbon case 825 G4cout <<"Carbon 12 (special) selected." << G4endl; 826 be->rms_be=2.44; 827 be->bind_be=92.17; 828 i_tabled=1; 829 } else if(be->iz_be == 6 && be->ia_be == 12) { 830 be->rms_be=2.44; 831 be->bind_be=92.17; 832 i_tabled=1; 833 } else if(be->iz_be == 7 && be->ia_be == 16) { 834 be->rms_be=2.73; 835 be->bind_be=127.62; 836 i_tabled=1; 837 } else { 838 G4cout <<"Warning: No rms and binding for projectile ion A = " << be->ia_be << " Z = " << be->iz_be << G4endl; 839 be->rms_be=2.44; 840 be->bind_be=92.17; 841 G4cout <<"Warning: Using probably bad values rms = " << be->rms_be << " binding = " << be->bind_be << G4endl; 842 i_tabled=1; 843 } 844 845 if(i_tabled == 0) { 846 G4cout <<"This heavy ion (a,z)= " << be->ia_be << " " << be->iz_be << " is not defined as beam in INCL" << G4endl; 847 exit(0); 848 } 849 850 // G4cout <<"z projectile, rms_r, rms_p (gaussian model)" << be->iz_be << " " << be->rms_be << " " << be->pms_be << G4endl; 851 // G4cout <<"binding energy (mev):" << be->bind_be << G4endl; 852 // G4cout <<"fermi-breakup dresner below a=" << calincl->f[11] << G4endl; 853 } 854 // G4cout <<"Target Mass and Charge: " << calincl->targetA() << " " << calincl->targetZ() << G4endl; 855 // calincl->f[10] = 0; // No clusters 856 857 if(calincl->bulletType() == -12) { // C12 special case 858 mprojo=fmp*std::abs(calincl->bulletType()) - be->bind_be; 859 pbeam=std::sqrt(calincl->bulletE()*(calincl->bulletE()+2.*mprojo)); 860 ap=std::abs(calincl->bulletType()); 861 zp=be->iz_be; 862 } else if(calincl->bulletType() == -666) { // Generic extended projectile 863 mprojo=fmp*be->ia_be - be->bind_be; 864 pbeam=std::sqrt(calincl->bulletE()*(calincl->bulletE()+2.*mprojo)); 865 ap=be->ia_be; 866 zp=be->iz_be; 867 } 696 868 // pi+ 697 if(calincl-> f[6] == 3.0) {869 if(calincl->bulletType() == 3) { 698 870 mprojo = 139.56995; 699 871 ap = 0.0; … … 702 874 703 875 // pi0 704 if(calincl-> f[6] == 4.0) {876 if(calincl->bulletType() == 4) { 705 877 mprojo = 134.9764; 706 878 ap = 0.0; … … 709 881 710 882 // pi- 711 if(calincl-> f[6] == 5.0) {883 if(calincl->bulletType() == 5) { 712 884 mprojo = 139.56995; 713 885 ap = 0.0; … … 718 890 719 891 // proton 720 if(calincl-> f[6] == 1.0) {892 if(calincl->bulletType() == 1) { 721 893 mprojo = 938.27231; 722 894 ap = 1.0; … … 725 897 726 898 // neutron 727 if(calincl-> f[6] == 2.0) {899 if(calincl->bulletType() == 2) { 728 900 mprojo = 939.56563; 729 901 ap = 1.0; … … 732 904 733 905 // deuteron 734 if(calincl-> f[6] == 6.0) {906 if(calincl->bulletType() == 6) { 735 907 mprojo = 1875.61276; 736 908 ap = 2.0; … … 739 911 740 912 // triton 741 if(calincl-> f[6] == 7.0) {913 if(calincl->bulletType() == 7) { 742 914 mprojo = 2808.95; 743 915 ap = 3.0; … … 746 918 747 919 // He3 748 if(calincl-> f[6] == 8.0) {920 if(calincl->bulletType() == 8) { 749 921 mprojo = 2808.42; 750 922 ap = 3.0; … … 753 925 754 926 // Alpha 755 if(calincl-> f[6] == 9.0) {927 if(calincl->bulletType() == 9) { 756 928 mprojo = 3727.42; 757 929 ap = 4.0; … … 759 931 } 760 932 761 pbeam = std::sqrt(calincl->f[2]*(calincl->f[2] + 2.0*mprojo)); 762 763 G4double at = calincl->f[0]; 933 // Carbon 934 if(calincl->bulletType() == -12) { 935 mprojo = 6.0*938.27231 + 6.0*939.56563; 936 ap = 12.0; 937 zp = 6.0; 938 } 939 940 pbeam = std::sqrt(calincl->bulletE()*(calincl->bulletE() + 2.0*mprojo)); 941 942 G4double at = calincl->targetA(); 764 943 765 calincl->f[3] = 0.0; // !seuil sortie proton766 calincl->f[7] = 0.0; // !seuil sortie neutron767 768 944 G4int ibert = 1; 769 945 … … 777 953 G4double bimpac = 0.0; 778 954 G4int jrem = 0; 955 G4double xjrem = 0.0, yjrem = 0.0, zjrem = 0.0; 779 956 G4double alrem = 0.0; 780 957 … … 784 961 G4double rndm = 0.0; 785 962 786 if((calincl-> f[6] == 1.0) || (calincl->f[6] >= 6.0)) {787 // probaTrans = coulombTransm(calincl-> f[2],apro,zpro,calincl->f[0],calincl->f[1]);788 probaTrans = coulombTransm(calincl-> f[2],ap,zp,calincl->f[0],calincl->f[1]);963 if((calincl->bulletType() == 1) || (calincl->bulletType() >= 6)) { 964 // probaTrans = coulombTransm(calincl->bulletE(),apro,zpro,calincl->targetA(),calincl->targetZ()); 965 probaTrans = coulombTransm(calincl->bulletE(),ap,zp,calincl->targetA(),calincl->targetZ()); 789 966 standardRandom(&rndm, &(hazard->ial)); 790 967 if(rndm <= (1.0 - probaTrans)) { … … 794 971 } 795 972 973 // G4cout <<"Before PNU:" << G4endl; 974 // randomGenerator->printSeeds(); 796 975 // Call the actual INCL routine: 797 pnu(&ibert, &nopart,&izrem,&iarem,&esrem,&erecrem,&alrem,&berem,&garem,&bimpac,&jrem); 976 pnu(&ibert, &nopart,&izrem,&iarem,&esrem,&erecrem,&alrem,&berem,&garem,&bimpac, 977 &jrem, &xjrem, &yjrem, &zjrem); 978 // G4cout <<"After PNU:" << G4endl; 979 // randomGenerator->printSeeds(); 980 G4double mrem = int(zjrem/197.328); // CHECK 981 if (mrem > jrem) mrem=jrem; 982 if (mrem < -jrem) mrem=-jrem; 983 798 984 // nopart=1; 799 985 // kind[0]=1; … … 846 1032 // -from energy balance (input - all emitted energies) 847 1033 // -following the approximations of the cugnon code (esrem...) 848 G4double mcorem = mprojo + calincl-> f[2] + abla->pace2(double(calincl->f[0]),double(calincl->f[1]))849 + calincl-> f[0]*uma - calincl->f[1]*melec;1034 G4double mcorem = mprojo + calincl->bulletE() + abla->pace2(double(calincl->targetA()),double(calincl->targetZ())) 1035 + calincl->targetA()*uma - calincl->targetZ()*melec; 850 1036 851 1037 G4double pxbil = 0.0; … … 854 1040 855 1041 if(nopart > -1) { 856 for(G4int j = 0; j < nopart; j++) { 857 if(ep[j] < 0.0) continue; // Workaround to avoid negative energies (and taking std::sqrt of a negative number). 858 varntp->itypcasc[j] = 1; 1042 // Fill the projectile spectator variables 1043 varntp->masp = ps->a_projspec; 1044 varntp->mzsp = ps->z_projspec; 1045 varntp->exsp = ps->ex_projspec; 1046 varntp->spectatorP1 = ps->p1_projspec; 1047 varntp->spectatorP2 = ps->p2_projspec; 1048 varntp->spectatorP3 = ps->p3_projspec; 1049 varntp->spectatorT = ps->t_projspec; 1050 for(G4int j = 0; j <= nopart; j++) { 1051 if(ep[j] < 0.0) { // Workaround to avoid negative energies (and taking std::sqrt of a negative number). 1052 G4cout <<"G4Incl: Not registering particle with energy: " << ep[j] << G4endl; 1053 continue; 1054 } 1055 if(kind[j] == 0) continue; // Empty particle rows are sometimes produced by lurking indexing problems. We can simply skip those "bad" entries... 1056 if(gam[j] > CLHEP::pi) { 1057 if(verboseLevel > 2) { 1058 G4cout <<"G4Incl: Just avoided floating point exception by using an ugly hack..." << G4endl; 1059 } 1060 continue; // Avoid floating point exception 1061 } 1062 1063 varntp->itypcasc[j] = 1; // Particle was produced by the cascade 1064 // Spectators of composite projectiles (7/2006, AB) 1065 // (kind is negative in that case) 1066 if(kind[j] <= 0) { // Particle is a projectile spectator that comes directly from the cascade 1067 kind[j] *= -1; 1068 varntp->itypcasc[j]=-1; 1069 // G4cout <<"Spectator registered!" << G4endl; 1070 // continue; 1071 } 1072 859 1073 // kind(): 1=proton, 2=neutron, 3=pi+, 4=pi0, 5=pi - 860 1074 if(kind[j] == 1) { … … 998 1212 pzbil = pzbil + pzrem; 999 1213 1214 // If on purpose, add here the spectator nuclei: 1215 if(calincl->bulletType() < 0 && ps->a_projspec != 0) { 1216 pxbil=pxbil+ps->p1_projspec; 1217 pybil=pybil+ps->p2_projspec; 1218 pzbil=pzbil+ps->p3_projspec; 1219 } 1220 1000 1221 if((std::fabs(pzbil - pbeam) > 5.0) || (std::sqrt(std::pow(pxbil,2) + std::pow(pybil,2)) >= 3.0)) { 1001 1222 if(verboseLevel > 3) { … … 1013 1234 1014 1235 // on recopie le remnant dans le ntuple 1015 //varntp->ntrack = varntp->ntrack + 1;1236 varntp->ntrack = varntp->ntrack + 1; 1016 1237 varntp->massini = iarem; 1017 1238 varntp->mzini = izrem; 1018 1239 varntp->exini = esrem; 1240 varntp->pxrem = pxrem; 1241 varntp->pyrem = pyrem; 1242 varntp->pzrem = pzrem; 1243 varntp->mcorem = mcorem; 1244 varntp->erecrem = pcorem; 1245 varntp->erecrem = erecrem; 1019 1246 1020 1247 if(verboseLevel > 2) { … … 1022 1249 varntp->dump(); 1023 1250 } 1024 // Evaporation/fission: 1025 evaporationResult->ntrack = 0; 1026 abla->breakItUp(varntp->massini, varntp->mzini, mcorem, varntp->exini, varntp->jremn, 1027 erecrem, pxrem, pyrem, pzrem, eventnumber); 1028 1029 if(verboseLevel > 2) { 1030 G4cout << __FILE__ << ":" << __LINE__ << "Dump evaporationResult after evaporation: " << G4endl; 1031 evaporationResult->dump(); 1032 } 1033 for(G4int evaporatedParticle = 0; evaporatedParticle < evaporationResult->ntrack; evaporatedParticle++) { 1034 if(evaporationResult->avv[evaporatedParticle] == 0 && evaporationResult->zvv[evaporatedParticle] == 0) { //Fix: Skip "empty" particles with A = 0 and Z = 0 1035 // continue; 1036 } 1037 varntp->kfis = evaporationResult->kfis; 1038 varntp->itypcasc[varntp->ntrack] = evaporationResult->itypcasc[evaporatedParticle]; 1039 varntp->avv[varntp->ntrack] = evaporationResult->avv[evaporatedParticle]; 1040 varntp->zvv[varntp->ntrack]= evaporationResult->zvv[evaporatedParticle]; 1041 varntp->plab[varntp->ntrack] = evaporationResult->plab[evaporatedParticle]; 1042 varntp->enerj[varntp->ntrack] = evaporationResult->enerj[evaporatedParticle]; 1043 varntp->tetlab[varntp->ntrack] = evaporationResult->tetlab[evaporatedParticle]; 1044 varntp->philab[varntp->ntrack] = evaporationResult->philab[evaporatedParticle]; 1045 varntp->ntrack++; 1046 if(verboseLevel > 3) { 1047 G4cout <<"G4Incl: Returning evaporation result" << G4endl; 1048 G4cout <<"G4Incl: A = " << varntp->avv[varntp->ntrack] << " Z = " << varntp->zvv[varntp->ntrack] << G4endl; 1049 G4cout <<"Energy: " << varntp->enerj[varntp->ntrack] << G4endl; 1050 G4cout <<"Momentum: " << varntp->plab[varntp->ntrack] << G4endl; 1051 G4cout <<"Theta: " << varntp->tetlab[varntp->ntrack] << G4endl; 1052 G4cout <<"Phi: " << varntp->philab[varntp->ntrack] << G4endl; 1053 } 1054 } 1055 if(verboseLevel > 2) { 1056 G4cout <<"G4Incl: ntrack = " << varntp->ntrack << G4endl; 1057 G4cout <<"G4Incl: Done extracting..." << G4endl; 1058 } 1251 // Maximum remnant for Geant4 Fermi break-up: A = 17, Z = 8 1252 if(varntp->massini < 17 && varntp->exini < 1000 && varntp->massini > 0 && varntp->mzini > 0 && useFermiBreakup) { // Choose between Fermi break-up and ABLA 1253 varntp->needsFermiBreakup = true; 1254 } else { 1255 varntp->needsFermiBreakup = false; 1256 // Evaporation/fission: 1257 evaporationResult->ntrack = 0; 1258 abla->breakItUp(G4int(varntp->massini), G4int(varntp->mzini), mcorem, varntp->exini, varntp->jremn, 1259 erecrem, pxrem, pyrem, pzrem, eventnumber); 1260 1261 if(verboseLevel > 2) { 1262 G4cout << __FILE__ << ":" << __LINE__ << "Dump evaporationResult after evaporation: " << G4endl; 1263 evaporationResult->dump(); 1264 } 1265 for(G4int evaporatedParticle = 0; evaporatedParticle < evaporationResult->ntrack; evaporatedParticle++) { 1266 if(evaporationResult->avv[evaporatedParticle] == 0 && evaporationResult->zvv[evaporatedParticle] == 0) { //Fix: Skip "empty" particles with A = 0 and Z = 0 1267 // continue; 1268 } 1269 varntp->kfis = evaporationResult->kfis; 1270 varntp->itypcasc[varntp->ntrack] = evaporationResult->itypcasc[evaporatedParticle]; 1271 varntp->avv[varntp->ntrack] = evaporationResult->avv[evaporatedParticle]; 1272 varntp->zvv[varntp->ntrack]= evaporationResult->zvv[evaporatedParticle]; 1273 varntp->plab[varntp->ntrack] = evaporationResult->plab[evaporatedParticle]; 1274 varntp->enerj[varntp->ntrack] = evaporationResult->enerj[evaporatedParticle]; 1275 varntp->tetlab[varntp->ntrack] = evaporationResult->tetlab[evaporatedParticle]; 1276 varntp->philab[varntp->ntrack] = evaporationResult->philab[evaporatedParticle]; 1277 varntp->ntrack++; 1278 if(verboseLevel > 3) { 1279 G4cout <<"G4Incl: Returning evaporation result" << G4endl; 1280 G4cout <<"G4Incl: A = " << varntp->avv[varntp->ntrack] << " Z = " << varntp->zvv[varntp->ntrack] << G4endl; 1281 G4cout <<"Energy: " << varntp->enerj[varntp->ntrack] << G4endl; 1282 G4cout <<"Momentum: " << varntp->plab[varntp->ntrack] << G4endl; 1283 G4cout <<"Theta: " << varntp->tetlab[varntp->ntrack] << G4endl; 1284 G4cout <<"Phi: " << varntp->philab[varntp->ntrack] << G4endl; 1285 } 1286 } 1287 if(verboseLevel > 2) { 1288 G4cout <<"G4Incl: ntrack = " << varntp->ntrack << G4endl; 1289 G4cout <<"G4Incl: Done extracting..." << G4endl; 1290 } 1291 } 1292 if(calincl->bulletType() < 0 && useProjSpect && !useFermiBreakup) { // If we use projectile spectators for carbon but no fermi breakup 1293 // Evaporation/fission: 1294 evaporationResult->ntrack = 0; 1295 // G4cout <<"Warning: Using ABLA to de-excite projectile spectator..." << G4endl; 1296 abla->breakItUp(G4int(varntp->masp), G4int(varntp->mzsp), ps->m_projspec, varntp->exsp, 0.0, 1297 ps->t_projspec, ps->p1_projspec, ps->p1_projspec, ps->p1_projspec, eventnumber); 1298 1299 if(verboseLevel > 2) { 1300 G4cout << __FILE__ << ":" << __LINE__ << "Dump evaporationResult after evaporation: " << G4endl; 1301 evaporationResult->dump(); 1302 } 1303 for(G4int evaporatedParticle = 0; evaporatedParticle < evaporationResult->ntrack; evaporatedParticle++) { 1304 if(evaporationResult->avv[evaporatedParticle] == 0 && evaporationResult->zvv[evaporatedParticle] == 0) { //Fix: Skip "empty" particles with A = 0 and Z = 0 1305 } 1306 varntp->kfis = evaporationResult->kfis; 1307 varntp->itypcasc[varntp->ntrack] = evaporationResult->itypcasc[evaporatedParticle]; 1308 varntp->avv[varntp->ntrack] = evaporationResult->avv[evaporatedParticle]; 1309 varntp->zvv[varntp->ntrack]= evaporationResult->zvv[evaporatedParticle]; 1310 varntp->plab[varntp->ntrack] = evaporationResult->plab[evaporatedParticle]; 1311 varntp->enerj[varntp->ntrack] = evaporationResult->enerj[evaporatedParticle]; 1312 varntp->tetlab[varntp->ntrack] = evaporationResult->tetlab[evaporatedParticle]; 1313 varntp->philab[varntp->ntrack] = evaporationResult->philab[evaporatedParticle]; 1314 varntp->ntrack++; 1315 if(verboseLevel > 3) { 1316 G4cout <<"G4Incl: Returning evaporation result" << G4endl; 1317 G4cout <<"G4Incl: A = " << varntp->avv[varntp->ntrack] << " Z = " << varntp->zvv[varntp->ntrack] << G4endl; 1318 G4cout <<"Energy: " << varntp->enerj[varntp->ntrack] << G4endl; 1319 G4cout <<"Momentum: " << varntp->plab[varntp->ntrack] << G4endl; 1320 G4cout <<"Theta: " << varntp->tetlab[varntp->ntrack] << G4endl; 1321 G4cout <<"Phi: " << varntp->philab[varntp->ntrack] << G4endl; 1322 } 1323 } 1324 } 1325 #ifdef G4INCLDEBUG 1326 theLogger->fillHistogram1D("bimpact", varntp->bimpact); 1327 #endif 1328 1329 #ifdef G4INCLDEBUG 1330 theLogger->fillHistogram1D("mzini", varntp->mzini); 1331 #endif 1059 1332 } 1060 1333 if(nopart == -2) { … … 1403 1676 1404 1677 if(tz < 0) { 1678 #ifdef G4INCLDEBUG 1679 theLogger->fillHistogram1D("interpolationResult", (saxw->y[0][saxw->imat] + saxw->s[0][saxw->imat]*tz)); 1680 theLogger->fillHistogram2D("interpolationPoints", (saxw->y[0][saxw->imat] + saxw->s[0][saxw->imat]*tz), xv); 1681 #endif 1405 1682 return (saxw->y[0][saxw->imat] + saxw->s[0][saxw->imat]*tz); 1406 1683 } 1407 1684 else if(tz == 0) { 1408 return (saxw->y[0][saxw->imat]); 1685 #ifdef G4INCLDEBUG 1686 theLogger->fillHistogram1D("interpolationResult", (saxw->y[0][saxw->imat])); 1687 theLogger->fillHistogram2D("interpolationPoints", (saxw->y[0][saxw->imat]), xv); 1688 #endif 1689 return (saxw->y[0][saxw->imat]); 1409 1690 } 1410 1691 else { // tz > 0 … … 1417 1698 } 1418 1699 if(tz >= 0) { 1700 #ifdef G4INCLDEBUG 1701 theLogger->fillHistogram1D("interpolationResult", saxw->y[j][saxw->imat]); 1702 theLogger->fillHistogram2D("interpolationPoints", saxw->y[j][saxw->imat], xv); 1703 #endif 1419 1704 return saxw->y[j][saxw->imat]; 1420 1705 } else if(tz < 0.0) { 1421 1706 j = j - 1; 1422 1707 G4double dgx = xv - saxw->x[j][saxw->imat]; 1708 #ifdef G4INCLDEBUG 1709 theLogger->fillHistogram1D("interpolationResult", (saxw->y[j][saxw->imat] + saxw->s[j][saxw->imat]*dgx)); 1710 theLogger->fillHistogram2D("interpolationPoints", (saxw->y[j][saxw->imat] + saxw->s[j][saxw->imat]*dgx), xv); 1711 #endif 1423 1712 return(saxw->y[j][saxw->imat] + saxw->s[j][saxw->imat]*dgx); 1424 1713 } … … 1673 1962 void G4Incl::pnu(G4int *ibert_p, G4int *nopart_p, G4int *izrem_p, G4int *iarem_p, G4double *esrem_p, 1674 1963 G4double *erecrem_p, G4double *alrem_p, G4double *berem_p, G4double *garem_p, 1675 G4double *bimpact_p, G4int *l_p) 1676 { 1964 G4double *bimpact_p, G4int *l_p, G4double *xjrem, G4double *yjrem, G4double *zjrem) 1965 { 1966 // G4cout <<"Now in pnu..." << G4endl; 1967 // G4cout <<"(clean-up)" << G4endl; 1968 G4int npenter = 0; 1969 G4int nnenter = 0; 1970 G4int n_enter_pot = 0; 1971 G4int avatarCounter = 0; 1972 ps->clear(); // For projectile spectators 1973 // G4int i_c = 0; 1974 G4double p1spec=0.; 1975 G4double p2spec=0.; 1976 G4double p3spec=0.; 1977 // G4double p_spec2 = 0.0; 1978 // G4double s_spec = 0.0; 1979 // G4double e_spec = 0.0; 1980 G4double xl1_i, xl2_i, xl3_i; 1981 G4int idq = 0; 1982 G4double destar = 0.0; 1983 // G4cout <<"(checkpoint 'setipszero0)" << G4endl; 1984 G4int ips = 0; 1985 // G4double p1_spec = 0.0; 1986 // G4double p2_spec = 0.0; 1987 // G4double p3_spec = 0.0; 1988 G4int n_activnuc = 0; 1677 1989 G4int ibert = (*ibert_p); 1678 1990 // float f[15]; // = (*f_p); … … 1702 2014 G4double b3 = 0.0; 1703 2015 G4double bb2 = 0.0; 1704 G4double be = 0.0;2016 G4double bevar = 0.0; 1705 2017 G4double bmass[2000]; 1706 2018 for(G4int init_i = 0; init_i < 2000; init_i++) { … … 1736 2048 G4double ener_max = 0.0; 1737 2049 G4double eout = 0.0; 1738 G4double eps_c[BL1SIZE]; 1739 for(G4int init_i = 0; init_i < BL1SIZE; init_i++) { 1740 eps_c[init_i] = 0.0; 1741 } 2050 // G4double eps_c[BL1SIZE]; 2051 // G4double p3_c[BL1SIZE]; 2052 // for(G4int init_i = 0; init_i < BL1SIZE; init_i++) { 2053 // eps_c[init_i] = 0.0; 2054 // p3_c[init_i] = 0.0; 2055 // } 1742 2056 G4double epsv = 0.0; 1743 2057 G4double erecg = 0.0; … … 1778 2092 G4int imin = 0; 1779 2093 G4int indic[3000]; 2094 G4int nb_transprojo=0; 2095 G4double v_proj = 0.0; 1780 2096 for(G4int init_i = 0; init_i < 3000; init_i++) { 1781 2097 indic[init_i] = 0; … … 1825 2141 G4int nbtest = 0; 1826 2142 G4int nc[300]; 2143 G4bool isPartOfSpectatorNucleus[300]; 1827 2144 G4int npproj[300]; 1828 2145 for(G4int init_i = 0; init_i < 300; init_i++) { … … 1942 2259 G4double xleng = 0.0; 1943 2260 G4double xlengm = 0.0; 1944 G4double xmodp = 0.0;1945 2261 G4double xpb = 0.0; 1946 2262 G4double xq = 0.0; … … 1960 2276 G4double xye = 0.0; 1961 2277 G4double y = 0.0; 1962 G4double p3_c[BL1SIZE];1963 2278 G4double q1[BL1SIZE]; 1964 2279 G4double q2[BL1SIZE]; … … 2220 2535 for(G4int i = 0; i < 300; i++) { 2221 2536 npproj[i] = 0; 2537 isPartOfSpectatorNucleus[i] = false; 2222 2538 nc[i] = 0; 2223 2539 } … … 2270 2586 // 427 c k6=0 no angular momentum conservation p-n02070 2271 2587 // 428 c k6=1 angular momentum conservation p-n02080 2588 // Kclst=F(10) Light clusters computed if not =0 2272 2589 // 429 c p-n02090 2273 2590 // 430 c b=impact parameter p-n02100 … … 2312 2629 // end for logging 2313 2630 2631 G4int ia_breakup = 0; 2314 2632 G4int jparticip[BL1SIZE]; 2315 2633 for(G4int i = 0; i < BL1SIZE; i++) { … … 2318 2636 2319 2637 G4double beproj = 0.; 2320 bl3->ia2 = G4int(std::floor(calincl->f[0] + 0.1)); // f(1)->f[0] and so on..., calincl added2321 G4int iz2 = G4int(std::floor(calincl->f[1] + 0.1));2638 bl3->ia2 = calincl->targetA(); // f(1)->f[0] and so on..., calincl added 2639 G4int iz2 = calincl->targetZ(); 2322 2640 G4double r02 = 1.12; 2323 kindstruct->kindf7 = int(std::floor(calincl->f[6] + 0.1)); 2324 2325 bl3->ia1 = ia1t[G4int(kindstruct->kindf7)-1]; 2326 G4int iz1 = iz1t[G4int(kindstruct->kindf7)-1]; 2327 G4double fmpinc = fmpinct[G4int(kindstruct->kindf7)-1]; 2328 G4double rms1 = light_gaus_nuc->rms1t[G4int(kindstruct->kindf7)-1]; 2329 G4double pf1 = light_gaus_nuc->pf1t[G4int(kindstruct->kindf7)-1]; 2330 G4double tlab = calincl->f[2]; 2641 kindstruct->kindf7 = calincl->bulletType(); 2642 if(kindstruct->kindf7 < 0) kindstruct->kindf7 = kindstruct->kindf7-1; 2643 G4int kclst = calincl->getClusterOption(); 2644 2645 // G4cout <<"Projectile kind = " << kindstruct->kindf7 << G4endl; 2646 2647 G4int iz1 = 0; 2648 G4double fmpinc = 0.0; 2649 G4double rms1 = 0.0; 2650 G4double pf1 = 0.0; 2651 G4double tlab = calincl->bulletE();; 2652 2653 if(kindstruct->kindf7 > 0) { 2654 bl3->ia1 = ia1t[G4int(kindstruct->kindf7)-1]; // In the C++ version indices start from 0 2655 iz1 = iz1t[G4int(kindstruct->kindf7)-1]; 2656 fmpinc = fmpinct[G4int(kindstruct->kindf7)-1]; 2657 rms1 = light_gaus_nuc->rms1t[G4int(kindstruct->kindf7)-1]; 2658 pf1 = light_gaus_nuc->pf1t[G4int(kindstruct->kindf7)-1]; 2659 2660 } else { 2661 // Starting values for carbon beams: 2662 if(kindstruct->kindf7 == -12) { // Handle Carbon-12 as a special case 2663 bl3->ia1 = std::abs(kindstruct->kindf7 + 1); 2664 // check here that kindf7 is -a and ia1 is a of the projectile. 2665 // write(6,*) 'kindf7,ia1:',kindf7,ia1 2666 be->ia_be=bl3->ia1; 2667 iz1=be->iz_be; 2668 fmpinc=bl3->ia1*fmp - be->bind_be; 2669 rms1=be->rms_be; 2670 pf1=be->pms_be; 2671 ia_breakup = 10; // Stopping threshold for the cascade in case of light nuclei 2672 // ia_breakup=calincl->getCascadeStoppingAThreshold(); 2673 } else { // For extended projectiles 2674 bl3->ia1 = calincl->extendedProjectileA(); 2675 be->ia_be = calincl->extendedProjectileA(); 2676 iz1=be->iz_be; 2677 fmpinc=bl3->ia1*fmp - be->bind_be; 2678 rms1=be->rms_be; 2679 pf1=be->pms_be; 2680 ia_breakup = 10; // Stopping threshold for the cascade in case of light nuclei 2681 // ia_breakup=calincl->getCascadeStoppingAThreshold(); 2682 } 2683 } 2684 2685 if(verboseLevel > 1) { 2686 G4cout <<"rms1 = " << rms1 << G4endl; 2687 G4cout <<"pf1 = " << pf1 << G4endl; 2688 } 2331 2689 2332 2690 G4int k2 = 0; … … 2338 2696 2339 2697 // material number: 2340 saxw->imat = G4int(std::floor(calincl->f[8] + 0.5)); // f(9) -> f[8]2698 saxw->imat = 0; 2341 2699 // espace de phases test (r et p) pour pauli: 2342 2700 // valeur recommandee par j.c. v-test=0.589 h**3: … … 2356 2714 // temfin (time at which the inc is stopped), tmax5 defined after chosing b 2357 2715 2358 G4double v0 = calincl-> f[4]; // f(5)->f[4]2716 G4double v0 = calincl->getPotential(); 2359 2717 G4double v1 = v0; 2360 2718 bl8->rathr = 0.; … … 2428 2786 G4cout <<"Radius bl3->r2 = " << bl3->r2 << G4endl; 2429 2787 } 2430 2788 2431 2789 G4double tnr = tlab; 2432 2790 G4double binc = std::sqrt(tnr*tnr + 2.0*tlab*fmpinc)/(tnr+fmpinc); 2433 G4double ginc =1.0/std::sqrt(1.0 - binc*binc);2791 G4double ginc = 1.0/std::sqrt(1.0 - binc*binc); 2434 2792 G4double pinc = fmpinc*binc*ginc; 2435 2793 … … 2463 2821 // 710 c if (kindf7.gt.2) bmax=bmax ! a.b. (avec w.s., idem les nucleons) 2464 2822 // deutons cv 22/01/2001 2465 if (kindstruct->kindf7 < = 2) { //then2823 if (kindstruct->kindf7 < 6 && kindstruct->kindf7 > 0) { //then 2466 2824 ws->bmax = ws->bmax; // comme alain 2467 2825 } 2468 2826 else { 2469 if (kindstruct->kindf7 < 6) { // then 2470 ws->bmax = ws->bmax; // comme alain 2471 } 2472 else { 2473 beproj = fmpinc - bl3->ia1*fmp; 2474 ws->bmax = ws->rmaxws + rms1; // maximum extension of the nucleus ( w.s.) 2475 } 2827 beproj = fmpinc - bl3->ia1*fmp; 2828 ws->bmax = ws->rmaxws + rms1; // maximum extension of the nucleus ( w.s.) 2476 2829 } 2477 2830 … … 2496 2849 bred = 0.; 2497 2850 } 2498 if (kindstruct->kindf7 <= 2 ) {2851 if (kindstruct->kindf7 <= 2 && kindstruct->kindf7 > 0) { 2499 2852 if (tlab < 400.0) { 2500 2853 cb0 = 6.86 - 0.0035 * tlab; … … 2515 2868 } 2516 2869 else { 2517 if (kindstruct->kindf7 < 6 ) {2870 if (kindstruct->kindf7 < 6 && kindstruct->kindf7 > 0) { 2518 2871 // here for pions: 2519 2872 temfin = 30.0*std::pow((float(bl3->ia2)/208.0),0.25)*(1.0 - 0.2*bred)*(1.0 - tlab/1250.0); … … 2557 2910 // temfin=60.*(float(ia2)/208.)**0.25*(1.-0.2*bred)*(1.-tlab/1250.) 2558 2911 // modified in april 2003 (more reasonable but not yet checked!) 2559 temfin = 25.5*std::pow(G4double(bl3->ia2),0.16); // pb208->60fm/c 2912 // temfin = 25.5*std::pow(G4double(bl3->ia2),0.16); // pb208->60fm/c 2913 2914 // c modified in June 2005 (function of A and TLAB) 2915 temfin = 30.18*std::pow(bl3->ia2,0.17)*(1.0-5.7*std::pow(10.0,(-5))*tlab); 2560 2916 } 2561 2917 else { … … 2567 2923 // deutons 2568 2924 // a way to change stopping time f[5] not used here 2569 factemp = calincl-> f[5]; // f(6)->f[5]2925 factemp = calincl->getTimeScale(); 2570 2926 // attention !!! 30/04/2001 scaling time is now a multiplication factor 2571 2927 temfin = temfin*factemp; … … 2604 2960 else { 2605 2961 npion = 1; 2606 ipi[1] = 8 - 2*(kindstruct->kindf7); 2962 if(kindstruct->kindf7 == -666) { // Extended projectile... 2963 ipi[1] = 8 - 2*be->ia_be; 2964 } else { 2965 ipi[1] = 8 - 2*(kindstruct->kindf7); 2966 } 2607 2967 y1[1] = 0.0; 2608 2968 y2[1] = 0.0; … … 2617 2977 goto pnu9; 2618 2978 // 850 2619 pnu7: 2620 s1t1 = 0.0; 2621 s2t1 = 0.0; 2622 s3t1 = 0.0; 2623 sp1t1 = 0.0; 2624 sp2t1 = 0.0; 2625 sp3t1 = 0.0; 2626 for(G4int i = 1; i <= bl3->ia1; i++) { 2627 //for(G4int i = 1; i <= (bl3->ia1-1); i++) { 2628 // for(G4int i = 1; i < (bl3->ia1-1); i++) { 2629 bl9->hel[i] = 0; 2630 bl5->nesc[i] = 0; 2631 bl1->ind2[i] = 1; 2632 if (i > iz1) { 2633 bl1->ind2[i] = -1; 2634 } 2635 bl1->ind1[i] = 0; 2636 // deutons 2637 jparticip[i] = 1; 2638 2639 // deutons 2979 pnu7: 2980 // G4cout <<"(checkpoint 'pnu7)" << G4endl; 2981 if(kindstruct->kindf7 == 6) { // Deuteron 2640 2982 gaussianRandom(&xga); 2641 bl3->x1[i] = xga*rms1*0.57735; 2642 s1t1 = s1t1 + bl3->x1[i]; 2983 bl3->x1[1] = xga*rms1*0.5775; 2643 2984 gaussianRandom(&xga); 2644 bl3->x2[i] = xga*rms1*0.57735; 2645 s2t1 = s2t1 + bl3->x2[i]; 2985 bl3->x2[1] = xga*rms1*0.5775; 2646 2986 gaussianRandom(&xga); 2647 bl3->x3[i] = xga*rms1*0.57735; 2648 s3t1 = s3t1 + bl3->x3[i]; 2649 2650 if(kindstruct->kindf7 == 6) { //then 2651 // deuteron density from paris potential in q space: 2652 standardRandom(&xq, &(hazard->igraine[9])); 2653 qdeut = splineab(xq) * 197.3289; 2654 standardRandom(&u, &(hazard->igraine[10])); 2655 cstet = u*2 - 1; 2656 sitet = std::sqrt(1.0 - std::pow(cstet,2)); 2657 standardRandom(&v, &(hazard->igraine[11])); 2658 phi = 2.0*3.141592654*v; 2659 2660 bl1->p1[i] = qdeut*sitet*std::cos(phi); 2661 bl1->p2[i] = qdeut*sitet*std::sin(phi); 2662 bl1->p3[i] = qdeut*cstet; 2663 } 2664 else { 2665 // density of composite as a gaussien in q space: 2987 bl3->x3[1] = xga*rms1*0.5775; 2988 2989 bl3->x1[2] = -bl3->x1[1]; 2990 bl3->x2[2] = -bl3->x2[1]; 2991 bl3->x3[2] = -bl3->x3[1]; 2992 2993 // Deuteron density from Paris potential in q space: 2994 standardRandom(&xq, &(hazard->igraine[9])); 2995 qdeut = splineab(xq) * 197.3289; 2996 standardRandom(&u, &(hazard->igraine[11])); 2997 cstet = u * 2 - 1; 2998 sitet = std::sqrt(1.0 - std::pow(cstet, 2)); 2999 standardRandom(&v, &(hazard->igraine[11])); 3000 phi = 2.0 * 3.141592654 * v; 3001 bl1->p1[1] = qdeut * sitet * std::cos(phi); 3002 bl1->p2[1] = qdeut * sitet * std::sin(phi); 3003 bl1->p3[1] = qdeut * cstet; 3004 bl1->eps[1] = w(bl1->p1[1], bl1->p2[1], bl1->p3[1], fmp); 3005 bl1->p1[2] = - bl1->p1[1]; 3006 bl1->p2[2] = - bl1->p2[1]; 3007 bl1->p3[2] = - bl1->p3[1]; 3008 bl1->eps[2] = bl1->eps[1]; 3009 // bl1->eps[2] = w(bl1->p1[2], bl1->p2[2], bl1->p3[2], fmp); 3010 3011 jparticip[1] = 1; 3012 bl9->hel[1] = 0; 3013 bl5->nesc[1] = 0; 3014 bl1->ind2[1] = 1; 3015 bl1->ind1[1] = 0; 3016 jparticip[2] = 1; 3017 bl9->hel[2] = 0; 3018 bl5->nesc[2] = 0; 3019 bl1->ind2[2] = -1; 3020 bl1->ind1[2] = 0; 3021 } else { 3022 // Composite heavier than a deuteron 3023 s1t1 = 0.0; 3024 s2t1 = 0.0; 3025 s3t1 = 0.0; 3026 sp1t1 = 0.0; 3027 sp2t1 = 0.0; 3028 sp3t1 = 0.0; 3029 G4double sei = 0.0; 3030 3031 bl9->hel[0] = 0; 3032 bl5->nesc[0] = 0; 3033 for(G4int i = 1; i <= bl3->ia1; i++) { 3034 bl9->hel[i] = 0; 3035 bl5->nesc[i] = 0; 3036 bl1->ind2[i] = 1; 3037 if(i > iz1) { 3038 bl1->ind2[i] = -1; 3039 } 3040 bl1->ind1[i] = 0; 3041 jparticip[i] = 1; 2666 3042 gaussianRandom(&xga); 2667 bl1->p1[i] = xga*pf1*0.57735; 3043 bl3->x1[i] = xga * rms1 * 0.57735; 3044 s1t1 = s1t1 + bl3->x1[i]; 2668 3045 gaussianRandom(&xga); 2669 bl1->p2[i] = xga*pf1*0.57735; 3046 bl3->x2[i] = xga * rms1 * 0.57735; 3047 s2t1 = s2t1 + bl3->x2[i]; 2670 3048 gaussianRandom(&xga); 2671 bl1->p3[i] = xga*pf1*0.57735; 2672 } 2673 2674 bl1->eps[i] = w(bl1->p1[i],bl1->p2[i],bl1->p3[i],fmp); 2675 2676 sp1t1 = sp1t1 + bl1->p1[i]; 2677 sp2t1 = sp2t1 + bl1->p2[i]; 2678 sp3t1 = sp3t1 + bl1->p3[i]; 2679 } 2680 2681 bl9->hel[bl3->ia1] = 0; 2682 bl5->nesc[bl3->ia1] = 0; 2683 bl1->ind2[bl3->ia1] = -1; 2684 bl1->ind1[bl3->ia1] = 0; 2685 bl3->x1[bl3->ia1] = -s1t1; 2686 bl3->x2[bl3->ia1] = -s2t1; 2687 bl3->x3[bl3->ia1] = -s3t1; 2688 bl1->p1[bl3->ia1] = -sp1t1; 2689 bl1->p2[bl3->ia1] = -sp2t1; 2690 bl1->p3[bl3->ia1] = -sp3t1; 2691 bl1->eps[bl3->ia1] = w(bl1->p1[bl3->ia1],bl1->p2[bl3->ia1],bl1->p3[bl3->ia1],fmp); 2692 2693 // deutons 2694 jparticip[bl3->ia1] = 1; 2695 if(verboseLevel > 3) { 2696 G4cout <<"Particle " << bl3->ia1-1 << " is now participant." << G4endl; 2697 } 3049 bl3->x3[i] = xga * rms1 * 0.57735; 3050 s3t1 = s3t1 + bl3->x3[i]; 3051 3052 G4int igaussfermi = 2; 3053 if(igaussfermi == 1) { // Gaussian 3054 // Density of composite as a gaussian in q 3055 // space: 3056 gaussianRandom(&xga); 3057 bl1->p1[i] = xga * pf1 * 0.57735; 3058 gaussianRandom(&xga); 3059 bl1->p2[i] = xga * pf1 * 0.57735; 3060 gaussianRandom(&xga); 3061 bl1->p3[i] = xga * pf1 * 0.57735; 3062 } else { // Fermi 3063 // Density of composite as a Fermi sphere in q space: 3064 for(int iloc=1; iloc <= 3; iloc++) { 3065 standardRandom(&t[iloc],&(hazard->igraine[11])); 3066 } 3067 t[2]=-1+2.*t[2]; 3068 t[3]=6.283185*t[3]; 3069 t2=t[2]; 3070 t3=std::sqrt(1.-t2*t2); 3071 t4=std::cos(t[3]); 3072 t5=std::sin(t[3]); 3073 y=fermi->pf*std::pow(t[1],0.33333333); 3074 bl1->p1[i]=y*t3*t4; 3075 bl1->p2[i]=y*t3*t5; 3076 bl1->p3[i]=y*t2; 3077 } 3078 bl1->eps[i] = w(bl1->p1[i], bl1->p2[i], bl1->p3[i], fmp); 3079 sei = sei + bl1->eps[i]; 3080 3081 if(verboseLevel > 1) { 3082 G4cout <<"projectile nucleon = " << i << G4endl; 3083 G4cout <<"x1 = " << bl3->x1[i] << G4endl; 3084 G4cout <<"x2 = " << bl3->x2[i] << G4endl; 3085 G4cout <<"x3 = " << bl3->x3[i] << G4endl; 3086 G4cout <<"p1 = " << bl1->p1[i] << G4endl; 3087 G4cout <<"p2 = " << bl1->p2[i] << G4endl; 3088 G4cout <<"p3 = " << bl1->p3[i] << G4endl; 3089 G4cout <<"eps = " << bl1->eps[i] << G4endl; 3090 } 3091 sp1t1 = sp1t1 + bl1->p1[i]; 3092 sp2t1 = sp2t1 + bl1->p2[i]; 3093 sp3t1 = sp3t1 + bl1->p3[i]; 3094 } 3095 3096 s1t1 = s1t1/double(bl3->ia1); 3097 s2t1 = s2t1/double(bl3->ia1); 3098 s3t1 = s3t1/double(bl3->ia1); 3099 sp1t1 = sp1t1/double(bl3->ia1); 3100 sp2t1 = sp2t1/double(bl3->ia1); 3101 sp3t1 = sp3t1/double(bl3->ia1); 3102 sei = 0.0; 3103 3104 for(G4int i = 1; i <= bl3->ia1; i++) { 3105 bl3->x1[i] = bl3->x1[i] - s1t1; 3106 bl3->x2[i] = bl3->x2[i] - s2t1; 3107 bl3->x3[i] = bl3->x3[i] - s3t1; 3108 bl1->p1[i] = bl1->p1[i] - sp1t1; 3109 bl1->p2[i] = bl1->p2[i] - sp2t1; 3110 bl1->p3[i] = bl1->p3[i] - sp3t1; 3111 bl1->eps[i] = w(bl1->p1[i], bl1->p2[i], bl1->p3[i], fmp); 3112 sei = sei + bl1->eps[i]; 3113 } 3114 v_proj = 0.0; 3115 if(bl3->ia1 > 4) v_proj = (sei - fmpinc)/bl3->ia1; 3116 3117 // bl1->dump(26); 3118 // bl3->dump(); 3119 // G4cout <<"(checkpoint 'endheavyprojinit)" << G4endl; 3120 } // End of: Composites heavier than deuteron 3121 3122 2698 3123 pnu9: // continue 3124 // G4cout <<"(checkpoint 'pnu9)" << G4endl; 2699 3125 // deutons 2700 3126 // target preparation for 1 < a < 5 (with sum of momentum =0) … … 2860 3286 // random azimuthal direction of the impact parameter (sept 99) 2861 3287 2862 if (kindstruct->kindf7 <= 2) { 2863 G4long testseed = hazard->igraine[13]; 2864 // standardRandom(&tbid, &(hazard->igraine[13])); 2865 standardRandom(&tbid, &testseed); 2866 hazard->igraine[13] = testseed; 3288 if (kindstruct->kindf7 <= 2 && kindstruct->kindf7 > 0) { 3289 standardRandom(&tbid, &(hazard->igraine[13])); 2867 3290 tbid = tbid*6.283185; 2868 3291 bl3->x1[1] = bl3->x1[1] + b*std::cos(tbid); //x1(1)->x1[1] 2869 3292 bl3->x2[1] = bl3->x2[1] + b*std::sin(tbid); //x2(1)->x2[1] 2870 3293 bl3->x3[1] = bl3->x3[1] - z; 3294 // Counter of p(n) entering the potential 4/2008 AB 3295 if(bl1->ind2[1] == 1) { // then 3296 npenter=1; 3297 nnenter=0; 3298 } else { 3299 npenter=0; 3300 nnenter=1; 3301 } // endif 3302 2871 3303 // pour le ntuple des avatars: 2872 3304 if(varavat->kveux == 1) { … … 2877 3309 } 2878 3310 else { 2879 if (kindstruct->kindf7 < 6 ) { //then ! pour les pions on laisse3311 if (kindstruct->kindf7 < 6 && kindstruct->kindf7 > 0) { //then ! pour les pions on laisse 2880 3312 //call standardRandom(tbid,iy14) 2881 3313 standardRandom(&tbid, &(hazard->igraine[13])); … … 2916 3348 // pnu21: 2917 3349 if (nmiss == bl3->ia1) { //then 3350 // G4cout <<"(checkpoint 'projectilemissedtarget)" << G4endl; 2918 3351 nopart = -1; 2919 3352 // return; … … 2925 3358 } 2926 3359 else { 3360 // Counter of p(n) entering the potential 4/2008 AB 3361 if(bl1->ind2[ilm]==1) { //then 3362 npenter=1; 3363 nnenter=0; 3364 } else { 3365 npenter=0; 3366 nnenter=1; 3367 } // endif 3368 2927 3369 zshif = bl3->x3[ilm] - ztouch; 2928 3370 standardRandom(&tbid, &(hazard->igraine[13])); … … 2943 3385 } 2944 3386 3387 // for rho(r),rho(q) checking 3388 #ifdef G4INCLDEBUG 3389 for(G4int rho_i = bl3->ia1 + 1; rho_i <= ia; rho_i++) { // DO I=IA1+1,IA 3390 G4double r_dist = std::sqrt(bl3->x1[rho_i] * bl3->x1[rho_i] + bl3->x2[rho_i] * bl3->x2[rho_i] + bl3->x3[rho_i] * bl3->x3[rho_i]); 3391 G4double p_mod = std::sqrt(std::pow(bl1->p1[rho_i],2)+std::pow(bl1->p2[rho_i],2) + std::pow(bl1->p3[rho_i],2)); 3392 theLogger->fillHistogram1D("r_distrib", r_dist); 3393 theLogger->fillHistogram1D("p_distrib", p_mod); 3394 theLogger->fillHistogram2D("r-p_correl", r_dist, p_mod); 3395 } 3396 #endif 3397 2945 3398 // initial momentum for all type of incident particles: 2946 3399 xl1 = b*pinc*std::sin(tbid); … … 2951 3404 // (here,=lab frame) 2952 3405 2953 be = 0.0;3406 bevar = 0.0; 2954 3407 ge = 1.0; 2955 b1 = (binc - be )/(1.0 - be*binc);2956 b2 = -be ;3408 b1 = (binc - bevar)/(1.0 - bevar*binc); 3409 b2 = -bevar; 2957 3410 g1 = 1.0/std::sqrt(1.0 - b1*b1); 2958 3411 g2 = 1.0; 2959 3412 // deutons 2960 3413 // here for nucleons 2961 if (kindstruct->kindf7 <= 2 ) {3414 if (kindstruct->kindf7 <= 2 && kindstruct->kindf7 > 0) { 2962 3415 bl1->eps[1] = g1*fmp + v0; 2963 3416 bl1->p3[1] = std::sqrt(std::pow(bl1->eps[1],2) - std::pow(fmp,2)); … … 2965 3418 else { 2966 3419 // here for pions 2967 if (kindstruct->kindf7 < 6 ) { //then3420 if (kindstruct->kindf7 < 6 && kindstruct->kindf7 > 0) { //then 2968 3421 q4[1] = g1*fmpi; // q4(1)->q4[0] 2969 3422 q3[1] = b1*q4[1]; 2970 3423 } 2971 else { 3424 else { // Composite projectiles 2972 3425 // here for composite projectiles: 2973 3426 // the kinetic energy is below the threshold. put all … … 2985 3438 // here the composit is above threshold 2986 3439 for(G4int i = 1; i <= bl3->ia1; i++) { //do i=1,bl3->ia1 !save e,p in the composit rest frame 2987 eps_c[i] = bl1->eps[i]; 2988 p3_c[i] = bl1->p3[i]; 3440 qvp->t_c[i] = bl1->eps[i] - fmp; 3441 bl1->eps[i] = bl1->eps[i] - v_proj; 3442 qvp->eps_c[i] = bl1->eps[i]; 3443 qvp->p3_c[i] = bl1->p3[i]; 2989 3444 } //enddo 2990 3445 … … 2995 3450 iflag = 0; 2996 3451 pnu1870: 3452 // G4cout <<"(checkpoint 'pnu1870)" << G4endl; 2997 3453 sueps = 0.0; 2998 3454 … … 3000 3456 3001 3457 for(G4int i = 1; i <= bl3->ia1; i++) { //do i=1,bl3->ia1 3002 tte = eps_c[i];3003 bl1->eps[i] = g1*( eps_c[i] + b1*p3_c[i]);3004 bl1->p3[i] = g1*(b1*tte + p3_c[i]);3458 tte = qvp->eps_c[i]; 3459 bl1->eps[i] = g1*(qvp->eps_c[i] + b1*qvp->p3_c[i]); 3460 bl1->p3[i] = g1*(b1*tte + qvp->p3_c[i]); 3005 3461 sueps = sueps + bl1->eps[i]; 3006 3462 } //enddo 3007 3463 3464 goto pnu987; 3465 // With a potential for the projectile it is nonsens to search for 3466 // nucleons ON shell 31/05/2010. 3467 3008 3468 cobe = (tlab + fmpinc)/sueps; 3009 3469 … … 3012 3472 if(iflag == nbtest) { // too much..all momentum to 0 3013 3473 for(G4int klm = 1; klm <= bl3->ia1; klm++) { //do klm=1,bl3->ia1 3014 eps_c[klm] = fmp;3474 qvp->eps_c[klm] = fmp; 3015 3475 bl1->p1[klm] = 0.0; 3016 3476 bl1->p2[klm] = 0.0; 3017 p3_c[klm] = 0;3477 qvp->p3_c[klm] = 0; 3018 3478 } 3019 3479 goto pnu1870; … … 3030 3490 } 3031 3491 } 3032 eps_c[i_emax] = fmp;3492 qvp->eps_c[i_emax] = fmp; 3033 3493 bl1->p1[i_emax] = 0.0; 3034 3494 bl1->p2[i_emax] = 0.0; 3035 p3_c[i_emax] = 0.0;3495 qvp->p3_c[i_emax] = 0.0; 3036 3496 3037 3497 if(i_emax == bl3->ia1) { // circular permut if the last one 3038 epsv = eps_c[bl3->ia1]; // permutation circulaire3498 epsv = qvp->eps_c[bl3->ia1]; // permutation circulaire 3039 3499 p1v = bl1->p1[bl3->ia1]; 3040 3500 p2v = bl1->p2[bl3->ia1]; 3041 p3v = p3_c[bl3->ia1];3501 p3v = qvp->p3_c[bl3->ia1]; 3042 3502 for(bl2->k = bl3->ia1-1; bl2->k >= 1; bl2->k = bl2->k - 1) { //do k=bl3->ia1-1,1,-1 3043 eps_c[bl2->k+1] =eps_c[bl2->k];3503 qvp->eps_c[bl2->k+1] = qvp->eps_c[bl2->k]; 3044 3504 bl1->p1[bl2->k+1] = bl1->p1[bl2->k]; 3045 3505 bl1->p2[bl2->k+1] = bl1->p2[bl2->k]; 3046 p3_c[bl2->k+1] =p3_c[bl2->k];3506 qvp->p3_c[bl2->k+1] = qvp->p3_c[bl2->k]; 3047 3507 } 3048 eps_c[1] = epsv;3508 qvp->eps_c[1] = epsv; 3049 3509 bl1->p1[1] = p1v; 3050 3510 bl1->p2[1] = p2v; 3051 p3_c[1] = p3v; // fin permut.3511 qvp->p3_c[1] = p3v; // fin permut. 3052 3512 } 3053 3513 sp1t1 = 0.0; // re-compute the last one … … 3057 3517 sp1t1 = sp1t1 + bl1->p1[j]; 3058 3518 sp2t1 = sp2t1 + bl1->p2[j]; 3059 sp3t1 = sp3t1 + p3_c[j];3519 sp3t1 = sp3t1 + qvp->p3_c[j]; 3060 3520 } 3061 3521 bl1->p1[bl3->ia1] = -sp1t1; 3062 3522 bl1->p2[bl3->ia1] = -sp2t1; 3063 p3_c[bl3->ia1] = -sp3t1;3064 eps_c[bl3->ia1] = w(bl1->p1[bl3->ia1],bl1->p2[bl3->ia1],p3_c[bl3->ia1],fmp);3523 qvp->p3_c[bl3->ia1] = -sp3t1; 3524 qvp->eps_c[bl3->ia1] = w(bl1->p1[bl3->ia1],bl1->p2[bl3->ia1],qvp->p3_c[bl3->ia1],fmp); 3065 3525 3066 3526 goto pnu1870; // ..and boost all of them. … … 3083 3543 } 3084 3544 bl1->eps[ilm] = bl1->eps[ilm] + v0; 3085 } 3086 } 3545 3546 pnu987: 3547 // G4cout <<"(checkpoint 'pnu987)" << G4endl; 3548 // ici sauvegarde des energies et impulsions des nucleons du projectile i.l.: 3549 // (hors potentiel) 3550 for(G4int i = 1; i <= bl3->ia1; i++) { 3551 qvp->eps_c[i]=bl1->eps[i]; 3552 qvp->p1_s[i]=bl1->p1[i]; 3553 qvp->p2_s[i]=bl1->p2[i]; 3554 qvp->p3_s[i]=bl1->p3[i]; 3555 } 3556 3557 // first particle entering in the target potential: 3558 bl1->eps[ilm]=bl1->eps[ilm]+v0; 3559 3560 // correction en test 16/04/2009 (pourquoi pas la correction aussi de l'impulsion 3561 // comme pour un nucleon in ???? ....on essaye de traiter les nucleons hors couche?) 3562 // 29/4/2009: je pense que c'est une erreur et je transforme aussi maintenant l'impulsion: 3563 var_ab = std::pow(bl1->p1[ilm],2) 3564 + std::pow(bl1->p2[ilm],2) 3565 + std::pow(bl1->p3[ilm],2); 3566 if(var_ab > 0.) { 3567 gpsg=std::sqrt((std::pow(bl1->eps[ilm],2) 3568 -fmp*fmp)/var_ab); 3569 } 3570 bl1->p1[ilm]=bl1->p1[ilm]*gpsg; 3571 bl1->p2[ilm]=bl1->p2[ilm]*gpsg; 3572 bl1->p3[ilm]=bl1->p3[ilm]*gpsg; 3573 } // For pion-in 3574 } // For nucleon-in 3087 3575 3088 3576 pnu1871: 3577 // G4cout <<"(checkpoint 'pnu1871)" << G4endl; 3578 // bl1->dump(27); 3579 // bl3->dump(); 3580 // randomGenerator->printSeeds(); 3089 3581 // evaluation of the times t(a,b) 3090 3582 bl2->k = 0; 3091 3583 kcol = 0; 3092 if (kindstruct->kindf7 <= 2 ) {3584 if (kindstruct->kindf7 <= 2 && kindstruct->kindf7 > 0) { 3093 3585 // modif s.vuillier tient compte propagation projectile,1e collision 3094 3586 // imposee pour lui (c'est une maniere de faire!!) … … 3155 3647 else { 3156 3648 // deutons 3157 if (kindstruct->kindf7 < 6 ) { //then3649 if (kindstruct->kindf7 < 6 && kindstruct->kindf7 > 0) { //then 3158 3650 // here for incoming pions: 3159 3651 for(G4int i = bl3->ia1+1; i <= ia; i++) { //do i=ia1+1,ia … … 3176 3668 } 3177 3669 else { 3670 // Counting the spectators and transparents of composit projectiles: 3671 // G4cout <<"(checkpoint 'setipszero1)" << G4endl; 3672 ips=0; 3673 nb_transprojo=0; 3178 3674 for(G4int i = 1; i <= bl3->ia1; i++) { //do 38 i=1,ia1 3179 3675 bl5->nesc[i] = 1; 3676 // Spectators of composite projectiles (7/2006, AB) 3677 // Pour garder trace des spectateurs du projectile: 3678 npproj[i]=1; 3679 3180 3680 if (i != ilm) { 3181 3681 goto pnu36; … … 3193 3693 goto pnu37; 3194 3694 pnu36: 3695 // G4cout <<"(checkpoint 'pnu36)" << G4endl; 3195 3696 t1 = bl3->x1[i]*(bl1->p1[i])+bl3->x2[i]*(bl1->p2[i])+bl3->x3[i]*(bl1->p3[i]); 3196 3697 t2 = bl1->p1[i]*(bl1->p1[i])+bl1->p2[i]*(bl1->p2[i])+bl1->p3[i]*(bl1->p3[i]); … … 3199 3700 // 1379 c incoming nucleons enter potential at maximum radius (modif. 13/06/01) 3200 3701 t5 = t3*t3 + ((ws->rmaxws)*(ws->rmaxws) - t4)/t2; 3702 if (t5 < 0.) { 3703 // this is a projectile spectator: 3704 // G4cout <<"(checkpoint 'incripspnu36)" << G4endl; 3705 ips=ips+1; 3706 isPartOfSpectatorNucleus[i] = true; 3707 ps->n_projspec[ips] = i; 3708 continue; 3709 //goto pnu38; 3710 } // endif 3711 3201 3712 if(verboseLevel > 3) { 3202 3713 G4cout <<"x1 = " << bl3->x1[i] <<" x2 = " << bl3->x2[i] <<" x3 = " << bl3->x3[i] << G4endl; … … 3208 3719 G4cout <<"t5 = " << t5 << G4endl; 3209 3720 } 3210 if (t5 < 0.) { 3211 continue; 3212 } 3721 // if (t5 < 0.) { 3722 // ips++; 3723 // ps->n_projspec[ips] = i; 3724 // continue; 3725 // } 3213 3726 tref = (-1.0*t3 - std::sqrt(t5))*(bl1->eps[i]); 3214 if (tref > bl4->tmax5) { 3215 continue; 3216 } 3727 3728 // Drop reflection avatars with too large tref and the ones 3729 // with negative tref (spectator, outside the nucleus) 3730 // if (tref > bl4->tmax5 || tref <= 0.0) { 3731 // continue; 3732 // } 3733 if (tref > bl4->tmax5) { //then 3734 // This is also a projectile spectator: 3735 ips=ips+1; 3736 ps->n_projspec[ips] = i; 3737 continue; // go to 38 3738 } // endif 3739 3740 // Test to avoid a nucleon entering in the past (3/06/2010 as in incl4.6 a.b.) 3741 if(tref <= 0.0) { // then 3742 // This is also a projectile spectator: 3743 ips=ips+1; 3744 ps->n_projspec[ips] = i; 3745 continue; // go to 38 3746 } // endif 3747 3217 3748 npproj[i] = 1; 3218 3749 pnu37: 3750 // G4cout <<"(checkpoint 'pnu37)" << G4endl; 3219 3751 bl2->k = bl2->k + 1; 3220 3752 bl2->crois[bl2->k] = tref; 3221 3753 bl2->ind[bl2->k] = i; 3222 3754 bl2->jnd[bl2->k] = -1; 3223 } 3755 // pnu38: 3756 } 3757 3758 // For composit beams, cannot decide that there is no interaction, 3759 // (another nucleon can enter later at RMAX). 3224 3760 kcol = 1; 3225 3761 n_activnuc = ia-ips; //number of active nucleons (not missing the target) 3762 n_enter_pot = ia-ips; //number of projectile nucleons that can enter 3763 3764 /* 3765 if(useProjSpect) { // Treat the projectile spectators: 3766 n_activnuc = ia - ips; // number of active nucleons (not missing the target) 3767 if(ips != 0) { //then 3768 e_spec = 0.0; 3769 p1_spec = 0.0; 3770 p2_spec = 0.0; 3771 p3_spec = 0.0; 3772 z_projspec = 0; 3773 for(G4int i_spec = 1; i_spec <= ips; i_spec++) { // do i_spec=1,ips 3774 G4int i_c = n_projspec[i_spec]; 3775 e_spec = e_spec + bl1->eps[i_c]; 3776 p1_spec = p1_spec + bl1->p1[i_c]; 3777 p2_spec = p2_spec + bl1->p2[i_c]; 3778 p3_spec = p3_spec + bl1->p3[i_c]; 3779 if(bl1->ind2[i_c] == 1) z_projspec = z_projspec + 1; 3780 } 3781 G4double p_spec2 = std::pow(p1_spec,2) + std::pow(p2_spec,2) + std::pow(p3_spec,2); 3782 G4double s_spec = std::sqrt(std::pow(e_spec,2) - p_spec2); 3783 // write(6,*) 'e,p,sqs:',e_spec,sqrt(p_spec2),s_spec 3784 3785 // c masses from dresner code: 3786 G4int fja = ips; 3787 G4int fjz = z_projspec; 3788 G4double ex2 = abla->eflmac(fja, fjz, 0, 0); // dnrgy(fja,fjz); 3789 G4double rnmass = ex2+um*fja; 3790 3791 // c write(6,*) fja,fjz,ex2,um 3792 ex_projspec = s_spec - rnmass; 3793 if(ex_projspec <= 0.0) { //then 3794 // c no spectator nucleus compatible with available energy: 3795 ips = 0; 3796 // c n_activnuc = ia-ips!number of active nucleons (not missing the target) 3797 goto pnu239; 3798 } 3799 t_projspec = e_spec - s_spec; 3800 G4double coef = std::sqrt(t_projspec*(t_projspec + 2.0*s_spec)/p_spec2); 3801 p1_projspec = p1_spec*coef; 3802 p2_projspec = p2_spec*coef; 3803 p3_projspec = p3_spec*coef; 3804 m_projspec = rnmass; 3805 3806 // c check module... 3807 // c write(6,*) 'projectile spectator a,z:',a_projspec,z_projspec 3808 // c write(6,*) 'm,t,e*:',m_projspec,t_projspec,ex_projspec 3809 // c write(6,*) 'p1,p2,p3:',p1_projspec,p2_projspec,p3_projspec 3810 // c write(6,*) ' constituants:',ips 3811 // c do i_spec=1,ips 3812 // c i_c = n_projspec(i_spec) 3813 // c write(6,*) i_c,eps(i_c),p1(i_c),p2(i_c),p3(i_c) 3814 // c enddo 3815 // c ...end check module 3816 } // endif 3817 pnu239: 3818 a_projspec = ips; // number of nucleons in the spectator nucleus 3819 } // End of projectile spectator treatment 3820 */ 3821 3226 3822 // for(G4int i = bl3->ia1+1; i < ia; i++) { //do 39 i=ia1+1,ia 3227 3823 for(G4int i = bl3->ia1+1; i <= ia; i++) { //do 39 i=ia1+1,ia … … 3298 3894 goto pnu48; 3299 3895 } 3896 // Here transparent event (no interaction avatar found) 3897 // G4cout <<"(checkpoint 'heretransparent)" << G4endl; 3300 3898 nopart = -1; 3301 3899 // Pour eviter renvoi des resultats du run precedent cv 7/7/98 … … 3316 3914 // Initialization at the beginning of the run 3317 3915 pnu48: 3916 // G4cout <<"(checkpoint 'pnu48)" << G4endl; 3318 3917 if(verboseLevel > 3) { 3319 3918 G4cout <<"Beginning a run..." << G4endl; … … 3364 3963 } 3365 3964 pnu449: 3965 // G4cout <<"(checkpoint 'pnu449)" << G4endl; 3366 3966 if(verboseLevel > 3) { 3367 3967 G4cout <<"Now at 449" << G4endl; … … 3372 3972 3373 3973 pnu44: 3974 // G4cout <<"(checkpoint 'pnu44)" << G4endl; 3374 3975 if(verboseLevel > 3) { 3375 3976 G4cout <<"Now at 44" << G4endl; … … 3406 4007 } 3407 4008 pnu448: 4009 // G4cout <<"(checkpoint 'pnu448)" << G4endl; 4010 // bl2->dump(); 3408 4011 imin = indic[next]; 3409 4012 bl9->l1 = bl2->ind[imin]; //NOTE: l1 changed to bl9->l1. … … 3419 4022 // test le 20/3/2003: tue sinon le dernier avatar? 3420 4023 if (bl2->k == 0) { 3421 if(verboseLevel > 2) {4024 if(verboseLevel > -2) { 3422 4025 G4cout <<"k == 0. Going to the end of the avatar." << G4endl; 3423 4026 } … … 3438 4041 } 3439 4042 pnu46: 4043 // G4cout <<"(checkpoint 'pnu46)" << G4endl; 3440 4044 if(verboseLevel > 3) { 3441 4045 G4cout <<"G4Incl: Now at pnu46." << G4endl; 3442 4046 } 3443 4047 tim = timi + tau; 4048 avatarCounter++; 4049 // G4cout <<"FOUND AVATAR " << avatarCounter << " : tau = " << tau; 4050 // G4cout <<" l1 = " << bl9->l1 << " l2 = " << bl9->l2 << G4endl; 3444 4051 3445 4052 // tableau des energies a t=20,40,60 fm/c … … 3496 4103 } 3497 4104 } 4105 // G4cout <<"(checkpoint 'pnu645)" << G4endl; 4106 // Here, stoping condition for heavy (light) ions: 4107 if(kindstruct->kindf7 < 0) { 4108 if(n_activnuc-nbquit <= ia_breakup) { 4109 // G4cout <<"stopped cascade for fermi breakup " << G4endl; 4110 goto pnu255; 4111 } 4112 } 3498 4113 3499 4114 // modif: pas de reflexions avant au moins un avatar du (des) nucleon incident 3500 4115 // celui-ci ne peut etre qu'une collision nn (ou pin) 3501 4116 4117 /* 3502 4118 if((irst_avatar == 0) && (bl9->l2 == -1)) { 3503 4119 if(verboseLevel > 3) { … … 3508 4124 3509 4125 irst_avatar = irst_avatar+1; 4126 */ 4127 // if((bl3->ia1 + bl3->ia2 - nbquit) <= 10) goto pnu255; // Stop the cascade if only very few nucleons are left. 3510 4128 3511 4129 if (tim < temfin) { … … 3518 4136 goto pnu255; 3519 4137 pnu49: 4138 // G4cout <<"(checkpoint 'pnu49)" << G4endl; 3520 4139 if(verboseLevel > 3) { 3521 4140 G4cout <<"G4Incl: Now at pnu49. " << G4endl; … … 3551 4170 // l1 est un delta: 3552 4171 pnu830: 4172 // G4cout <<"(checkpoint 'pnu830)" << G4endl; 3553 4173 if(verboseLevel > 3) { 3554 4174 G4cout <<"G4Incl: Now at pnu830." << G4endl; … … 3566 4186 } 3567 4187 pnu803: 4188 // G4cout <<"(checkpoint 'pnu803)" << G4endl; 3568 4189 // pas de collision entre 2 non participants: 3569 4190 if(jparticip[bl9->l1] == 0 && jparticip[bl9->l2] == 0) { … … 3632 4253 goto pnu261; 3633 4254 pnu260: 4255 // G4cout <<"(checkpoint 'pnu260)" << G4endl; 3634 4256 bmax2 = totalCrossSection(sq,mg,isos)/31.41592; 3635 4257 pnu261: 4258 // G4cout <<"(checkpoint 'pnu261)" << G4endl; 3636 4259 if (bb2 < bmax2) { 3637 4260 goto pnu220; … … 3648 4271 // evaluation of the positions at time = tim 3649 4272 pnu220: 4273 // G4cout <<"(checkpoint 'pnu220)" << G4endl; 3650 4274 timi = tim; 3651 4275 if(verboseLevel > 3) { … … 3682 4306 3683 4307 // gel des nucleons non participants sur le premier avatar (nn)=(l1,1) 4308 /* 3684 4309 if (irst_avatar == 1) { 3685 4310 for(G4int i = 1; i <= bl9->l1; i = i + bl9->l1 - 1) { // bugfix! … … 3697 4322 } 3698 4323 else { 4324 */ 3699 4325 for(G4int i = 1; i <= ia; i++) { 3700 4326 bl1->ta = tau/bl1->eps[i]; … … 3703 4329 bl3->x3[i] = bl3->x3[i] + bl1->p3[i]*(bl1->ta); 3704 4330 } 3705 }4331 // } 3706 4332 3707 4333 // if(npion != 0) { … … 3714 4340 } 3715 4341 } 4342 // G4cout <<"(checkpoint 'pnu840)" << G4endl; 3716 4343 3717 4344 if(bl9->l2 == 0) { … … 3828 4455 } 3829 4456 pnu243: 4457 // G4cout <<"(checkpoint 'pnu243)" << G4endl; 3830 4458 if (bl1->ind1[bl9->l2] == 1) { // bugfix 1 -> 0 3831 4459 goto pnu241; … … 3838 4466 goto pnu241; 3839 4467 pnu248: 4468 // G4cout <<"(checkpoint 'pnu248)" << G4endl; 3840 4469 if(verboseLevel > 3) { 3841 4470 G4cout <<"Pauli blocked transition!" << G4endl; … … 3874 4503 3875 4504 pnu241: 3876 4505 // G4cout <<"(checkpoint 'pnu241)" << G4endl; 3877 4506 // la premiere collision a deux corps ne peut pas baisser l'energie 3878 4507 // du nucleon de recul (bloque par pauli dans un noyau cible froid). … … 3974 4603 q4[npion] = t[14]; //t(15)->t[14] 3975 4604 pnu240: 4605 // G4cout <<"(checkpoint 'pnu240)" << G4endl; 3976 4606 ncol = ncol + 1; 3977 4607 if (bl9->l2 != 1) { … … 4045 4675 4046 4676 pnu870: 4677 // G4cout <<"(checkpoint 'pnu870)" << G4endl; 4047 4678 bl5->tlg[bl9->l1] = th*(bl1->eps[bl9->l1])/std::sqrt(std::pow(bl1->eps[bl9->l1],2)-std::pow(bl1->p1[bl9->l1],2)-std::pow(bl1->p2[bl9->l1],2)-std::pow(bl1->p3[bl9->l1],2)); 4048 4679 bl5->tlg[bl9->l2] = th*(bl1->eps[bl9->l2])/std::sqrt(std::pow(bl1->eps[bl9->l2],2)-std::pow(bl1->p1[bl9->l2],2)-std::pow(bl1->p2[bl9->l2],2)-std::pow(bl1->p3[bl9->l2],2)); … … 4250 4881 4251 4882 pnu848: 4883 // G4cout <<"(checkpoint 'pnu848)" << G4endl; 4252 4884 if (npion == 0) { 4253 4885 goto pnu844; … … 4269 4901 4270 4902 pnu843: 4903 // G4cout <<"(checkpoint 'pnu843)" << G4endl; 4271 4904 if(bl1->ind1[bl9->l2] != 1) { 4272 4905 for(G4int k20 = 1; k20 <= npion; k20++) { … … 4276 4909 } 4277 4910 pnu844: 4278 4911 // G4cout <<"(checkpoint 'pnu844)" << G4endl; 4279 4912 if(bl1->ind1[bl9->l1]+bl1->ind1[bl9->l2] <= ich1+ich2) { 4280 4913 goto pnu849; … … 4286 4919 goto pnu821; 4287 4920 pnu820: 4921 // G4cout <<"(checkpoint 'pnu820)" << G4endl; 4288 4922 if(bl1->ind1[bl9->l2]-ich2 != 1) { 4289 4923 goto pnu849; … … 4292 4926 4293 4927 pnu821: 4928 // G4cout <<"(checkpoint 'pnu821)" << G4endl; 4294 4929 standardRandom(&rndm,&(hazard->igraine[16])); 4295 4930 // largeur variable du delta (phase space factor G4introduced 4/2001) … … 4320 4955 // decay of the delta particle p-n09780 4321 4956 pnu805: 4957 // G4cout <<"(checkpoint 'pnu805)" << G4endl; 4322 4958 npion = npion + 1; 4323 4959 ichd = bl1->ind2[bl9->l1]; … … 4372 5008 4373 5009 pnu837: 5010 // G4cout <<"(checkpoint 'pnu837)" << G4endl; 4374 5011 ipi[npion]=bl1->ind2[bl9->l1]*2; 4375 5012 bl1->ind2[bl9->l1]=-1*(bl1->ind2[bl9->l1]); 4376 5013 goto pnu809; 4377 5014 pnu806: 5015 // G4cout <<"(checkpoint 'pnu806)" << G4endl; 4378 5016 bl1->ind2[bl9->l1]=bl1->ind2[bl9->l1]/3; 4379 5017 ipi[npion]=2*(bl1->ind2[bl9->l1]); 4380 5018 pnu809: // continue 5019 // G4cout <<"(checkpoint 'pnu809)" << G4endl; 4381 5020 bl1->ind1[bl9->l1]=0; 4382 5021 bl5->tlg[bl9->l1]=0.; … … 4422 5061 // it is blocked! 4423 5062 pnu1848: 5063 // G4cout <<"(checkpoint 'pnu1848)" << G4endl; 4424 5064 mpaul2 = mpaul2 + 1; 4425 5065 if(varavat->kveux == 1) { … … 4467 5107 // valid decay of the delta 4468 5108 pnu850: 5109 // G4cout <<"(checkpoint 'pnu850)" << G4endl; 4469 5110 if(varavat->kveux == 1) { 4470 5111 varavat->bloc_paul[iavat] = 0; … … 4514 5155 4515 5156 pnu4047: 5157 // G4cout <<"(checkpoint 'pnu4047)" << G4endl; 4516 5158 if (bl5->nesc[bl9->l1] != 0) { 4517 5159 goto pnu845; … … 4538 5180 4539 5181 pnu845: 5182 // G4cout <<"(checkpoint 'pnu845)" << G4endl; 4540 5183 new2(y1[npion], y2[npion], y3[npion], q1[npion], q2[npion], q3[npion], q4[npion], npion, bl9->l1); 4541 5184 if(bl2->k == 0) { … … 4550 5193 // pion-nucleon collision 4551 5194 pnu801: 5195 // G4cout <<"(checkpoint 'pnu801)" << G4endl; 4552 5196 if(verboseLevel > 3) { 4553 5197 G4cout <<"Pion-nucleon collision!" << G4endl; … … 4588 5232 4589 5233 pnu832: 5234 // G4cout <<"(checkpoint 'pnu832)" << G4endl; 4590 5235 if (bl2->k == 0) { 4591 5236 goto pnu230; … … 4597 5242 goto pnu44; 4598 5243 pnu831: 5244 // G4cout <<"(checkpoint 'pnu831)" << G4endl; 4599 5245 standardRandom(&rndm, &(hazard->igraine[18])); 4600 5246 geff = t[9]/sq; //t(10)->t[9] … … 4714 5360 // reflection on or transmission through the potential wall 4715 5361 pnu600: 5362 // G4cout <<"(checkpoint 'pnu600)" << G4endl; 4716 5363 // deutons pas bien compris ici cv ? 4717 5364 if (npproj[bl9->l1] == 0) { … … 4736 5383 } 4737 5384 5385 if (bl1->ind2[bl9->l1] == 1) { //then 5386 npenter=npenter+1; 5387 } else { 5388 nnenter=nnenter+1; 5389 } // endif 5390 4738 5391 var_ab = std::pow(bl1->p1[bl9->l1],2) + std::pow(bl1->p2[bl9->l1],2) + std::pow(bl1->p3[bl9->l1],2); 4739 5392 gpsg = 0.0; … … 4756 5409 // pour un non participant la transmission est impossible: 4757 5410 pnu608: 5411 // G4cout <<"(checkpoint 'pnu608)" << G4endl; 4758 5412 if(varavat->kveux == 1) { 4759 5413 varavat->del1avat[iavat] = bl1->ind1[bl9->l1]; … … 4777 5431 4778 5432 pnu605: 5433 // G4cout <<"(checkpoint 'pnu605)" << G4endl; 4779 5434 fm = fmp; 4780 5435 pot = v0; 4781 5436 4782 5437 pnu606: 5438 // G4cout <<"(checkpoint 'pnu606)" << G4endl; 4783 5439 if(verboseLevel > 3) { 4784 5440 G4cout <<"G4Incl: Now at pnu606. Calculating transmission probability." << G4endl; 4785 5441 } 4786 tp = transmissionProb(bl1->eps[bl9->l1]-fm,bl1->ind2[bl9->l1], itch,bl3->r2,v0);5442 tp = transmissionProb(bl1->eps[bl9->l1]-fm,bl1->ind2[bl9->l1],1,itch,bl3->r2,v0); 4787 5443 if(varavat->kveux == 1) { 4788 5444 varavat->energyavat[iavat] = bl1->eps[bl9->l1] - fm; … … 4796 5452 goto pnu601; 4797 5453 } 4798 5454 // Commande clusters ou non (1/0)! 5455 if(kclst == 0) goto pnu1610; 5456 // Cluster goes here! 5457 pnu1610: 5458 // G4cout <<"(checkpoint 'pnu1610)" << G4endl; 4799 5459 // ici la particule l1 s'échappe du noyau: 4800 5460 bl5->nesc[bl9->l1] = 1; … … 4811 5471 bl1->eps[bl9->l1] = bl1->eps[bl9->l1] - pot; 4812 5472 5473 // Test for a transparent nucleon for composit beams and stored for FermiBreakup 5474 if(bl3->ia1 > 1 && nc[bl9->l1] == 0) { 5475 // G4cout <<"(checkpoint 'incripspnu1610)" << G4endl; 5476 ips=ips+1; 5477 isPartOfSpectatorNucleus[bl9->l1] = true; 5478 ps->n_projspec[ips]=bl9->l1; 5479 nb_transprojo = nb_transprojo + 1; 5480 } // endif 5481 5482 // pnu610: goes here 5483 4813 5484 // comptage des particules hors du noyau (7/6/2002): 4814 5485 // (remnant minimum=1 nucleon) … … 4827 5498 // here no transmission possible 4828 5499 pnu601: 5500 // G4cout <<"(checkpoint 'pnu601)" << G4endl; 4829 5501 pspr=bl3->x1[bl9->l1]*(bl1->p1[bl9->l1])+bl3->x2[bl9->l1]*(bl1->p2[bl9->l1])+bl3->x3[bl9->l1]*(bl1->p3[bl9->l1]); 4830 5502 if(varavat->kveux == 1) { … … 4841 5513 4842 5514 pnu602: 5515 // G4cout <<"(checkpoint 'pnu602)" << G4endl; 4843 5516 if(verboseLevel > 3) { 4844 5517 G4cout <<"G4Incl: Now at pnu602." << G4endl; 4845 5518 } 4846 5519 4847 if(bl2->k != 0) { 5520 if(bl2->k != 0) { // GOTO elimination? 4848 5521 kd = 0; 4849 5522 ccr = tau; … … 4862 5535 bl2->k = bl2->k - kd; 4863 5536 4864 if (bl5->nesc[bl9->l1] == 1) {5537 if (bl5->nesc[bl9->l1] != 0) { // modif 10/02 pour logique des clusters (nesc()=2!) 4865 5538 goto pnu613; 4866 5539 } … … 4898 5571 4899 5572 pnu615: 5573 // G4cout <<"(checkpoint 'pnu615)" << G4endl; 4900 5574 if(verboseLevel > 3) { 4901 5575 G4cout <<"G4Incl: Now at pnu615." << G4endl; … … 4912 5586 } 4913 5587 pnu613: 5588 // G4cout <<"(checkpoint 'pnu613)" << G4endl; 4914 5589 if(verboseLevel > 3) { 4915 5590 G4cout <<"G4Incl: Now at pnu613." << G4endl; … … 4948 5623 // decay of the surviving deltas 4949 5624 pnu230: 5625 // G4cout <<"(checkpoint 'pnu230)" << G4endl; 4950 5626 if(verboseLevel > 3) { 4951 5627 G4cout <<"G4Incl: Now at pnu230." << G4endl; … … 4953 5629 // decay of the surviving deltas 4954 5630 pnu255: 5631 // G4cout <<"(checkpoint 'pnu255)" << G4endl; 4955 5632 if(verboseLevel > 3) { 4956 5633 G4cout <<"G4Incl: Now at pnu6255." << G4endl; … … 4965 5642 4966 5643 npidir = npion; 5644 5645 idq = 0; 5646 destar = 0.; 5647 4967 5648 for(G4int i = 1; i <= ia; i++) { 5649 G4int iblcdpp=0; 4968 5650 if (bl1->ind1[i] == 0) { 4969 5651 continue; … … 4980 5662 xy3 = bl1->p3[i]; 4981 5663 xye = bl1->eps[i]; 5664 ichd=bl1->ind2[i]; 4982 5665 if(varavat->kveux == 1) { 4983 5666 iavat = iavat + 1; … … 5000 5683 G4cout <<"i = " << i << " bl1->p1[i] = " << bl1->p1[i] << " bl1->p2[i] = " << bl1->p2[i] << " bl1->p3[i] = " << bl1->p3[i] << " bl1->eps[i] = " << bl1->eps[i] << G4endl; 5001 5684 } 5002 5685 5686 // C 6/08/09; Don't have to check CDPP for the forced decay of a DELTA OUTSIDE! 5687 // IF(NESC(I).EQ.1) GO TO 1850 5688 5689 // C------ 5690 // C On teste CDPP apres la decroissance forcee June/2005 (AB-CV) 5691 // C Si rejet de la decroissance -> delta devient nucleon sans pion 5692 // C et toute l'energie en exces est mise en energie d'excitation du noyau 5693 // C On respecte aussi la conservation de charge (Z du remnant!) 5694 5695 // C Le decay ne peut conduire a un noyau de A nucleons 5696 // C sous l'energie de Fermi et dans une config. d'energie inferieure a 5697 // C EFER-(IA2-NBALTTF)*TF). 5698 // EGS=0. 5699 // NBALTTF=0 5700 // DO II=1,IA 5701 // IF(NESC(II).EQ.0) THEN 5702 // IF(SQRT(P1(II)**2+P2(II)**2+P3(II)**2).LT.PF) THEN 5703 // NBALTTF=NBALTTF+1 5704 // EGS=EGS+EPS(II)-FMP 5705 // ENDIF 5706 // ENDIF 5707 // END DO 5708 // IF(EGS.GE.(EFER-(IA2-NBALTTF)*TF)) GO TO 1850 5709 // C ATTENTION, logique negative!!! Liberer le goto si on veut supprimer la 5710 // C sequence precedente (CDPP sur Delta-> pi N) 5711 // C GO TO 1850 5712 5713 // C Ici ce decay est interdit par CDPP. 5714 // C Restitution du Delta: un neutron au repos et correction conservee pour 5715 // C la charge IDQ et l'energie d'excitation DESTAR du remnant. On viole donc 5716 // C un peu la conservation d'impulsion, mais semble moins grave que de convertir 5717 // C la masse d'un pion en impulsion d'un nucleon unique.... 5718 // IBLCDPP=1 !flag de blocage par CDPP 5719 // IND1(i) = 0 !identif delta (nucleon) 5720 // IND2(i) = -1 !isospin du neutron 5721 // IDQ = IDQ + (ICHD+1)/2 5722 // p1(i)=0. !impulsion energie 5723 // p2(i)=0. 5724 // p3(i)=0. 5725 // eps(i)=FMP 5726 // DESTAR = DESTAR + xye - FMP 5727 // NPION=NPION-1 !pas le dernier pion 5728 5729 // 1850 CONTINUE 5730 // C------ 5003 5731 if(bl5->nesc[i] == 0) { 5004 5732 idecf = 1; … … 5022 5750 // fin surface 5023 5751 } 5752 if(iblcdpp == 1) continue; // goto pnu257; 5024 5753 5025 5754 if (bl1->ind2[i]*(bl1->ind2[i]) == 9) { … … 5036 5765 5037 5766 pnu283: 5767 // G4cout <<"(checkpoint 'pnu283)" << G4endl; 5038 5768 if(verboseLevel > 3) { 5039 5769 G4cout <<"G4Incl: Now at pnu283." << G4endl; … … 5045 5775 5046 5776 pnu280: 5777 // G4cout <<"(checkpoint 'pnu285)" << G4endl; 5047 5778 if(verboseLevel > 3) { 5048 5779 G4cout <<"G4Incl: Now at pnu280." << G4endl; … … 5053 5784 5054 5785 pnu285: 5786 // G4cout <<"(checkpoint 'pnu285)" << G4endl; 5055 5787 if(verboseLevel > 3) { 5056 5788 G4cout <<"G4Incl: Now at pnu285." << G4endl; … … 5094 5826 5095 5827 pnu256: 5828 // G4cout <<"(checkpoint 'pnu256)" << G4endl; 5096 5829 if(verboseLevel > 3) { 5097 5830 G4cout <<"G4Incl: Now at pnu256." << G4endl; … … 5102 5835 npx = 0; 5103 5836 erem = 0.; 5837 // Excitation energy of the final delta of blocked decay: (AB CV 6/2005) 5838 erem = erem + destar; 5104 5839 izrem = 0; 5840 // Correct the charge of the remnant if the final delta is blocked: (A.B., C.V. 6/2005) 5841 izrem = izrem + idq; 5105 5842 inrem = 0; 5106 5843 iarem = 0; … … 5115 5852 pout3 = 0.0; 5116 5853 eout = 0.0; 5854 p1spec=0.; 5855 p2spec=0.; 5856 p3spec=0.; 5117 5857 cmultn = 0.0; 5118 5858 5119 if (kindstruct->kindf7 <= 2 ) {5859 if (kindstruct->kindf7 <= 2 && kindstruct->kindf7 > 0) { 5120 5860 if (ncol == 0 || nc[1] == 0) { // then nc(1)->nc[0] 5121 5861 if(verboseLevel > 3) { … … 5126 5866 } 5127 5867 else { 5128 if (kindstruct->kindf7 <= 5 ) {5868 if (kindstruct->kindf7 <= 5 && kindstruct->kindf7 > 0) { 5129 5869 if (ncol == 0) { 5130 5870 if(verboseLevel > 3) { … … 5149 5889 // pour eviter renvoi des resultats du run precedent cv 20/11/98 5150 5890 pnu9100: 5891 // G4cout <<"(checkpoint 'pnu9100)" << G4endl; 5151 5892 iarem = bl3->ia2; 5152 5893 izrem = iz2; … … 5161 5902 5162 5903 pnu9101: 5904 // G4cout <<"(checkpoint 'pnu9101)" << G4endl; 5163 5905 if(verboseLevel > 3) { 5164 5906 G4cout <<"G4Incl: Now at pnu9101." << G4endl; … … 5168 5910 ekout = 0.0; 5169 5911 5170 for(G4int i = 1; i <= ia; i++) { 5171 if(bl5->nesc[i] != 0) { 5172 xl1 = xl1-bl3->x2[i]*(bl1->p3[i]) + (bl3->x3[i])*(bl1->p2[i]); 5173 xl2 = xl2-bl3->x3[i]*(bl1->p1[i]) + bl3->x1[i]*(bl1->p3[i]); 5174 xl3 = xl3-bl3->x1[i]*(bl1->p2[i]) + bl3->x2[i]*(bl1->p1[i]); 5175 5176 // ici ajout de pout cv le 5/7/95 5177 pout1 = pout1 + bl1->p1[i]; 5178 pout2 = pout2 + bl1->p2[i]; 5179 pout3 = pout3 + bl1->p3[i]; 5180 eout = eout+bl1->eps[i] - fmp; 5181 // ic33 = int((std::floor(double(bl1->ind2[i]+3))/double(2))); 5182 ic33 = (bl1->ind2[i]+3)/2; 5183 //nopart = nopart + 1; 5184 kind[nopart] = 3 - ic33; 5185 if(verboseLevel > 3) { 5186 G4cout <<"kind[" << nopart << "] = " << kind[nopart] << G4endl; 5187 } 5188 ep[nopart] = bl1->eps[i] - fmp; 5189 if(verboseLevel > 2) { 5190 if(ep[nopart] > calincl->f[2]) { 5191 G4cout <<"ep[" << nopart << "] = " << ep[nopart] << " > " << calincl->f[2] << G4endl; 5192 G4cout <<"E = " << bl1->eps[nopart] << G4endl; 5193 G4cout <<"px = " << bl1->p1[nopart] << " py = " << bl1->p2[nopart] << " pz = " << bl1->p3[nopart] << G4endl; 5194 G4cout <<"Particle mass = " << fmp << G4endl; 5195 } 5196 } 5197 bmass[nopart] = fmp; 5198 ekout = ekout + ep[nopart]; 5199 ptotl = std::sqrt(std::pow(bl1->p1[i],2) + std::pow(bl1->p2[i],2) + std::pow(bl1->p3[i],2)); 5200 alpha[nopart] = bl1->p1[i]/ptotl; 5201 beta[nopart] = bl1->p2[i]/ptotl; 5202 gam[nopart] = bl1->p3[i]/ptotl; 5203 nopart = nopart + 1; 5204 continue; 5205 } 5206 5207 t[3] = bl3->x1[i]*(bl3->x1[i]) + bl3->x2[i]*(bl3->x2[i]) + bl3->x3[i]*(bl3->x3[i]); //t(4)->t[3] 5208 erem = erem + bl1->eps[i] - fmp; 5209 rcm1 = rcm1 + bl3->x1[i]; 5210 rcm2 = rcm2 + bl3->x2[i]; 5211 rcm3 = rcm3 + bl3->x3[i]; 5212 prem1 = prem1 + bl1->p1[i]; 5213 prem2 = prem2 + bl1->p2[i]; 5214 prem3 = prem3 + bl1->p3[i]; 5215 izrem = izrem + (1 + bl1->ind2[i])/2; 5216 iarem = iarem + 1; 5217 } 5218 5912 if(kclst != 0) { // Clusters go here! 5913 } 5914 5915 // Treat here the projectile spectators (including transparents): 5916 projo_spec(bl3->ia1,ips,fmpinc,pinc,tlab); 5917 5918 for(G4int i = 1; i <= ia; i++) { // do 258 i=1,ia 5919 // pnu259: 5920 if (bl5->nesc[i] == 0) goto pnu254; 5921 if (bl5->nesc[i] > 1) continue; // on evite les nucleons des clusters 5922 5923 // here nesc=1 or negative value (emitted nucleons, spectators included) 5924 xl1_i=bl3->x2[i]*bl1->p3[i]-bl3->x3[i]*bl1->p2[i]; // moment angulaire emporte 5925 xl2_i=bl3->x3[i]*bl1->p1[i]-bl3->x1[i]*bl1->p3[i]; 5926 xl3_i=bl3->x1[i]*bl1->p2[i]-bl3->x2[i]*bl1->p1[i]; 5927 5928 xl1=xl1-xl1_i; 5929 xl2=xl2-xl2_i; 5930 xl3=xl3-xl3_i; 5931 5932 // ici ajout de pout cv le 5/7/95 5933 pout1=pout1+bl1->p1[i]; 5934 pout2=pout2+bl1->p2[i]; 5935 pout3=pout3+bl1->p3[i]; 5936 // write(6,*)'eout nucleon',eout 5937 eout = eout + bl1->eps[i]-fmp; 5938 5939 ic33=(bl1->ind2[i]+3)/2; 5940 5941 // here we don't copy the spectators of the projectile bound in a 5942 // spectator nucleus (treated by the deexcitation module): 5943 if(kindstruct->kindf7 < 0 && i <= bl3->ia1 && ps->a_projspec != 0) { 5944 for(G4int iloc = 1; iloc <= ps->a_projspec; iloc++) { //do iloc=1,a_projspec 5945 if(i == ps->n_projspec[iloc]) { 5946 p1spec=p1spec+bl1->p1[i]; 5947 p2spec=p2spec+bl1->p2[i]; 5948 p3spec=p3spec+bl1->p3[i]; 5949 // for balance, a transparent from projectile is out: 5950 if(i > ips-nb_transprojo) ekout=ekout+bl1->eps[i]-fmp; 5951 goto pnu258_end_of_outer_loop; 5952 // continue; 5953 } // endif 5954 } // enddo 5955 } // endif 5956 // if (isPartOfSpectatorNucleus[i]) goto pnu258_end_of_outer_loop; 5957 5958 nopart=nopart+1; 5959 5960 ekout=ekout+bl1->eps[i]-fmp; 5961 5962 kind[nopart]=3-ic33; 5963 // spectators of composite projectiles (7/2006, ab) 5964 if(npproj[i] == 1) kind[nopart] = -kind[nopart]; 5965 5966 ep[nopart]=bl1->eps[i]-fmp; 5967 bmass[nopart]=fmp; 5968 ptotl=std::sqrt(std::pow(bl1->p1[i],2) + std::pow(bl1->p2[i],2)+std::pow(bl1->p3[i],2)); 5969 alpha[nopart]=bl1->p1[i]/ptotl; 5970 beta[nopart]=bl1->p2[i]/ptotl; 5971 gam[nopart]=bl1->p3[i]/ptotl; 5972 5973 continue; 5974 pnu254: 5975 t[4]=bl3->x1[i]*bl3->x1[i]+bl3->x2[i]*bl3->x2[i]+bl3->x3[i]*bl3->x3[i]; 5976 erem=erem+bl1->eps[i]-fmp; 5977 rcm1=rcm1+bl3->x1[i]; 5978 rcm2=rcm2+bl3->x2[i]; 5979 rcm3=rcm3+bl3->x3[i]; 5980 prem1=prem1+bl1->p1[i]; 5981 prem2=prem2+bl1->p2[i]; 5982 prem3=prem3+bl1->p3[i]; 5983 izrem=izrem+(1+bl1->ind2[i])/2; 5984 iarem=iarem+1; 5985 pnu258_end_of_outer_loop: 5986 continue; 5987 } 5988 // 258 continue 5989 5990 //pnu258: 5991 // G4cout <<"(checkpoint 'pnu258)" << G4endl; 5219 5992 // correction pions 21/3/95 jc 5220 5993 ichpion = 0; … … 5229 6002 xl3 = xl3 - y1[ipion]*q2[ipion] + y2[ipion]*q1[ipion]; 5230 6003 ichpion = ichpion + ipi[ipion]/2; 5231 //nopart = nopart + 1;6004 nopart = nopart + 1; 5232 6005 kind[nopart] = 4 - ipi[ipion]/2; 5233 6006 ptotl = std::sqrt(std::pow(q1[ipion],2) + std::pow(q2[ipion],2) + std::pow(q3[ipion],2)); … … 5238 6011 beta[nopart] = q2[ipion]/ptotl; 5239 6012 gam[nopart] = q3[ipion]/ptotl; 5240 nopart++;6013 // nopart++; 5241 6014 } 5242 6015 } … … 5354 6127 } 5355 6128 5356 for(G4int ipart = 0; ipart < nopart; ipart++) {6129 for(G4int ipart = 0; ipart <= nopart; ipart++) { 5357 6130 ep[ipart] = ep[ipart]*fffc; 5358 6131 } … … 5361 6134 // modif boudard juillet 99 (il faut tenir compte de la renormalisation 5362 6135 // des energies pour les impulsions.) 5363 pfrem1 = 0.0; 5364 pfrem2 = 0.0; 5365 pfrem3 = pinc; 5366 for(G4int ipart = 0; ipart < nopart; ipart++) { 5367 xmodp = std::sqrt(ep[ipart]*(2.0*bmass[ipart] + ep[ipart])); 6136 pfrem1= -p1spec; 6137 pfrem2= -p2spec; 6138 pfrem3=pinc-p3spec; 6139 // if(useProjSpect) { 6140 // pfrem1 = -p1_spec; 6141 // pfrem2 = -p2_spec; 6142 // pfrem3 = pinc - p3_spec; 6143 // } else { 6144 // pfrem1 = 0.0; 6145 // pfrem2 = 0.0; 6146 // pfrem3 = pinc; 6147 // } 6148 for(G4int ipart = 0; ipart <= nopart; ipart++) { 6149 G4double xmodp = 0.0; 6150 if(ep[ipart] < 0.0) { // Safeguard against particles with negative energy 6151 xmodp = 0.0; 6152 G4cout <<"The energy of particle " << ipart << " is negative (" << ep[ipart] << "). Forcing xmodp = " << xmodp << "." << G4endl; 6153 } else { 6154 xmodp = std::sqrt(ep[ipart]*(2.0*bmass[ipart] + ep[ipart])); 6155 } 5368 6156 pfrem1 = pfrem1 - alpha[ipart]*xmodp; 5369 6157 pfrem2 = pfrem2 - beta[ipart]*xmodp; … … 5394 6182 esrem = 0.0; 5395 6183 } 6184 goto pnureturn; 5396 6185 5397 6186 if(verboseLevel > 3) { … … 5404 6193 } 5405 6194 pnureturn: 6195 // G4cout <<"(checkpoint 'pnureturn)" << G4endl; 6196 // if(useProjSpect) { 6197 // varntp->spectatorA = a_projspec; 6198 // varntp->spectatorZ = z_projspec; 6199 // varntp->spectatorEx = ex_projspec; 6200 // varntp->spectatorT = t_projspec; 6201 // varntp->spectatorM = m_projspec; 6202 // varntp->spectatorP1 = p1_projspec; 6203 // varntp->spectatorP2 = p1_projspec; 6204 // varntp->spectatorP3 = p1_projspec; 6205 // } 5406 6206 (*ibert_p) = ibert; 5407 6207 (*nopart_p) = nopart; … … 5415 6215 (*bimpact_p) = bimpact; 5416 6216 (*l_p) = l; 5417 } 5418 6217 (*xjrem) = 0; 6218 (*yjrem) = 0; 6219 (*zjrem) = 0; 6220 } 5419 6221 5420 6222 void G4Incl::collis(G4double *p1_p, G4double *p2_p, G4double *p3_p, G4double *e1_p, G4double *pout11_p, G4double *pout12_p, … … 6567 7369 // pauli strict 6568 7370 pmod = std::sqrt(std::pow(bl1->p1[l],2) + std::pow(bl1->p2[l],2) + std::pow(bl1->p3[l],2)); 6569 if (pmod < 270.0) {7371 if (pmod < fermi->pf) { 6570 7372 return 1.0; 6571 7373 } else { … … 6674 7476 G4double sine = 0.0; 6675 7477 6676 if((m-1) < 0) { 7478 if((m-1) < 0) { // Nucleon - Nucleon 6677 7479 sine = deltaProductionCrossSection(E,int(i)); 6678 7480 } 6679 7481 6680 if((m-1) == 0) { 7482 if((m-1) == 0) { // Nucleon - Delta 6681 7483 sine = srec(E,(bl6->xx10),i,int(bl6->isa)); 6682 7484 } 6683 7485 6684 if((m-1) > 0) { 7486 if((m-1) > 0) { // Delta - Delta 6685 7487 sine = 0.0; 6686 7488 } 6687 7489 6688 7490 stotResult = sine + lowEnergy(E,m,i); 7491 #ifdef G4INCLDEBUG 7492 theLogger->fillHistogram1D("totalCXResult", stotResult); 7493 theLogger->fillHistogram1D("totalLowECX", lowEnergy(E, m, i)); 7494 if(m == 0) { // Nucleon - Nucleon 7495 if(i == 2) { // pp 7496 theLogger->fillHistogram2D("ppCX", E, stotResult); 7497 } else if(i == 0) { // pn 7498 theLogger->fillHistogram2D("pnCX", E, stotResult); 7499 } else if(i == -2) { // nn 7500 theLogger->fillHistogram2D("nnCX", E, stotResult); 7501 } 7502 } else if(m == 1) { // Nucleon - Delta 7503 theLogger->fillHistogram2D("NDCX", E, stotResult); 7504 } else if(m == 2) { // Delta - Delta 7505 theLogger->fillHistogram2D("DDCX", E, stotResult); 7506 } 7507 #endif 6689 7508 return stotResult; 6690 7509 } … … 6858 7677 } 6859 7678 6860 G4double G4Incl::transmissionProb(G4double E, G4double iz, G4double izn, G4double r, G4double v0) 6861 { 6862 // transmission probability for a nucleon of kinetic energy 6863 // E on the edge of the well of depth v0 (nr approximation) 6864 // iz is the isospin of the nucleon,izn the instanteneous charge 6865 // of the nucleus and r is the target radius 6866 6867 G4double x = 0.0; 7679 // FUNCTION BARR(E,IZ,IA,IZN,R,V0) 7680 // tp = transmissionProb(bl1->eps[bl9->l1]-fm,bl1->ind2[bl9->l1],1,itch,bl3->r2,v0); 7681 G4double G4Incl::transmissionProb(G4double E, G4int iz, G4int ia, G4int izn,G4double R, G4double v0) 7682 { 7683 // C 7684 // C 7685 // C NOUVEAU BARR RELATIVISTE ATTENTION AUX ENTREES SORTIES! 7686 // C (06/2005) 7687 // C P-N24450 7688 // CCC BARR=TRANSMISSION PROBABILITY FOR A PARTICLE (Nucleon, cluster or 7689 // CCC pion) of kinetic energy E on the edge of the (attractive) well of 7690 // CCC depth V0 felt by the particle. 7691 // CCC IZ is the isospin of the particle, 7692 // CCC IZN the instantaneous charge of the nucleus and R the radius of 7693 // CCC the well. 7694 // CCC IA is the mass number of the particle, and MRE its mass energy. 7695 // C P-N24500 7696 // C Modified 9/10/2002 for clusters (d,t,3He,4He) (IZ=isospin,IA=A) 7697 // C IZ must be correct so that charge Q=(IA+IZ)/2. 7698 // C Modified 4/2004 for relativistic expressions and pions. 7699 7700 G4double mre; 6868 7701 G4double barr = 0.0; 6869 6870 // We need enough energy to escape from the potential well. 6871 if (E > v0) { 6872 x = std::sqrt(E*(E - v0)); 6873 barr = 4.0*x/(E + E - v0 + x + x); 6874 if (iz > 0 && izn != 0) { // izn = 0 causes division by zero 6875 G4double b = izn*1.44/r; 6876 G4double px = std::sqrt((E - v0)/b); 6877 6878 if (px < 1.0) { 6879 G4double g = izn/137.03*std::sqrt(2.0*938.3/(E - v0))*(std::acos(px) - px*std::sqrt(1.0 - px*px)); 6880 if (g > 35.){ 6881 barr=0.0; 6882 } else { 6883 barr = barr*std::exp(-2.0*g); 6884 } 6885 return barr; 6886 } else { 6887 return barr; 6888 } 6889 } else { 6890 return barr; 6891 } 7702 G4double v0ia = v0*ia; 7703 mre = ia*938.3; 7704 7705 if (E > v0ia && izn > 0 && R != 0) goto tp1; 7706 barr = 0.0; 7707 return barr; 7708 7709 tp1: 7710 G4double x=std::sqrt((2.*mre*E+E*E)*(2.*mre*(E-v0ia)+std::pow((E-v0ia),2))); 7711 barr=4.*x/(2.*mre*(2.*E-v0ia)+std::pow(E,2)+std::pow(E-v0ia,2)+2.*x); 7712 G4double iq=(ia+iz)/2.; 7713 if (iq > 0 && E > v0ia) goto tp2; 7714 return barr; 7715 7716 tp2: 7717 G4double b=iq*izn*1.44/R; 7718 G4double px=std::sqrt((E-v0ia)/b); 7719 if (px < 1.) goto tp3; 7720 return barr; 7721 7722 tp3: 7723 G4double g=iq*izn/137.03*std::sqrt(2.*mre/(E-v0ia)/(1.+(E-v0ia)/2./mre)) 7724 *(std::acos(px)-px*std::sqrt(1.-px*px)); 7725 if (g > 35.) { 7726 barr=0.0; 6892 7727 } else { 6893 return barr; 6894 } 7728 barr=barr*std::exp(-2.*g); 7729 } // endif 7730 return barr; 6895 7731 } 6896 7732 … … 6908 7744 G4double s_l = 0.0; 6909 7745 G4double t1 = 0.0, t3 = 0.0, t4 = 0.0, t5 = 0.0; 6910 if (ws->nosurf <= 0) { // modif pour w.s.: 6911 xv = p/pf; 6912 if(t2 <= pf2) { 6913 r = interpolateFunction(xv); 6914 } else { 6915 r = ws->rmaxws; 6916 } 6917 r = r*r; 6918 } 7746 if (ws->nosurf <= 0) { // modif pour w.s.: 7747 xv = p/pf; 7748 if(t2 <= pf2) { 7749 r = interpolateFunction(xv); 7750 #ifdef G4INCLDEBUG 7751 theLogger->fillHistogram1D("refInterpolXV", xv); 7752 theLogger->fillHistogram1D("refInterpolRes", r); 7753 if(xv > 1.0) theLogger->fillHistogram1D("refInterpolResWhenXVgt1", r); 7754 #endif 7755 } else { 7756 r = ws->rmaxws; 7757 } 7758 r = r*r; 7759 } 7760 /* 7761 if (ws->nosurf <= 0) { 7762 // modif pour w.s.: 7763 xv=p/pf; 7764 r=interpolateFunction(xv); 7765 r=r*r; 7766 if (t2 > pf2) r=std::pow(ws->rmaxws,2); 7767 } //endif 7768 */ 6919 7769 ref21: 6920 7770 t4 = x1*x1 + x2*x2 + x3*x3; … … 6969 7819 } 6970 7820 6971 bl3->ia2 = int(std::floor(calincl->f[0] + 0.1)); // f(1) -> f[0] 7821 // bl3->ia2 = int(std::floor(calincl->targetA() + 0.1)); // f(1) -> f[0] 7822 bl3->ia2 = calincl->targetA(); 6972 7823 sep = 6.8309; 6973 7824 … … 6976 7827 itg = 6 - 1; 6977 7828 } 6978 if(bl3->ia2 == 3 && calincl-> f[1]== 1) {7829 if(bl3->ia2 == 3 && calincl->targetZ() == 1) { 6979 7830 itg = 7 - 1; 6980 7831 } 6981 if(bl3->ia2 == 3 && calincl-> f[1]== 2) {7832 if(bl3->ia2 == 3 && calincl->targetZ() == 2) { 6982 7833 itg = 8 - 1; 6983 7834 } … … 6988 7839 } 6989 7840 6990 if((calincl-> f[2] >= 10.0) && (calincl->f[2]<= 100.0)) {6991 if(calincl-> f[6] == 1.0) {7841 if((calincl->bulletE() >= 10.0) && (calincl->bulletE() <= 100.0)) { 7842 if(calincl->bulletType() == 1) { 6992 7843 bl3->ia1 = int(1.0); 6993 7844 iz1 = 1.0; 6994 7845 G4double fmpinc = 938.2796; 6995 G4double pbeam2 = calincl-> f[2]*(calincl->f[2]+ 2.0*fmpinc);7846 G4double pbeam2 = calincl->bulletE()*(calincl->bulletE() + 2.0*fmpinc); 6996 7847 bmaxt = ws->bmax; 6997 proba_trans = coulombTransm(calincl-> f[2],bl3->ia1,iz1,calincl->f[0],calincl->f[1]);6998 proba = forceAbs(1,calincl-> f[0],calincl->f[1],calincl->f[2],bmaxt,proba_trans);7848 proba_trans = coulombTransm(calincl->bulletE(),bl3->ia1,iz1,calincl->targetA(),calincl->targetZ()); 7849 proba = forceAbs(1,calincl->targetA(),calincl->targetZ(),calincl->bulletE(),bmaxt,proba_trans); 6999 7850 7000 7851 standardRandom(&alea,&(hazard->igraine[4])); … … 7003 7854 } 7004 7855 7005 (*iarem) = int(calincl-> f[0]) + bl3->ia1;7006 (*izrem) = int(calincl-> f[1]) + int(iz1);7856 (*iarem) = int(calincl->targetA()) + bl3->ia1; 7857 (*izrem) = int(calincl->targetZ()) + int(iz1); 7007 7858 7008 del = std::sqrt(std::pow(((calincl-> f[0] + 1.0)*fmpinc + calincl->f[2]),2) - pbeam2);7859 del = std::sqrt(std::pow(((calincl->targetA() + 1.0)*fmpinc + calincl->bulletE()),2) - pbeam2); 7009 7860 7010 (*erecrem) = pbeam2/((calincl-> f[0] + 1.0)*fmpinc+calincl->f[2]+ del);7011 7012 (*esrem) = calincl-> f[2]+ sep - (*erecrem);7861 (*erecrem) = pbeam2/((calincl->targetA() + 1.0)*fmpinc+calincl->bulletE() + del); 7862 7863 (*esrem) = calincl->bulletE() + sep - (*erecrem); 7013 7864 7014 7865 (*alrem) = 0.00001; … … 7019 7870 return; 7020 7871 } 7021 else if((calincl-> f[6] == 2) && (calincl->f[2]>= 20.0)) {7872 else if((calincl->bulletType() == 2) && (calincl->bulletE() >= 20.0)) { 7022 7873 bl3->ia1 = int(1.0); 7023 7874 iz1 = 0.0; 7024 7875 G4double fmpinc = 938.2796; 7025 G4double pbeam2 = calincl-> f[2]*(calincl->f[2]+ 2.0*fmpinc);7876 G4double pbeam2 = calincl->bulletE()*(calincl->bulletE() + 2.0*fmpinc); 7026 7877 bmaxt = ws->bmax; 7027 proba_trans = coulombTransm(calincl-> f[2],bl3->ia1,iz1,calincl->f[0],calincl->f[1]);7028 proba = forceAbs(1,calincl-> f[0],calincl->f[1],calincl->f[2],bmaxt,proba_trans);7878 proba_trans = coulombTransm(calincl->bulletE(),bl3->ia1,iz1,calincl->targetA(),calincl->targetZ()); 7879 proba = forceAbs(1,calincl->targetA(),calincl->targetZ(),calincl->bulletE(),bmaxt,proba_trans); 7029 7880 7030 7881 standardRandom(&alea,&(hazard->igraine[4])); … … 7033 7884 } 7034 7885 7035 (*iarem) = int(calincl-> f[0]) + bl3->ia1;7036 (*izrem) = int(calincl-> f[1]) + int(iz1);7886 (*iarem) = int(calincl->targetA()) + bl3->ia1; 7887 (*izrem) = int(calincl->targetZ()) + int(iz1); 7037 7888 7038 del = std::sqrt(std::pow(((calincl-> f[0]+1.)*fmpinc+calincl->f[2]),2)-pbeam2);7889 del = std::sqrt(std::pow(((calincl->targetA()+1.)*fmpinc+calincl->bulletE()),2)-pbeam2); 7039 7890 7040 (*erecrem) = pbeam2/((calincl-> f[0] + 1.0)*fmpinc + calincl->f[2]+ del);7041 7042 (*esrem) = calincl-> f[2]+ sep - (*erecrem);7891 (*erecrem) = pbeam2/((calincl->targetA() + 1.0)*fmpinc + calincl->bulletE() + del); 7892 7893 (*esrem) = calincl->bulletE() + sep - (*erecrem); 7043 7894 7044 7895 (*alrem) = 0.00001; … … 7108 7959 return dp0; 7109 7960 } 7110 G4double rp = radius( ap);7111 G4double rt = radius( at);7961 G4double rp = radius(G4int(ap)); 7962 G4double rt = radius(G4int(at)); 7112 7963 G4double vp = (dp1 + dpth)*dppi*std::pow(rp,3); 7113 7964 G4double vt = (dp1 + dpth)*dppi*std::pow(rt,3); … … 7310 8161 { 7311 8162 // Gaussian random number generator 7312 8163 G4double xg, xxg, xshuf; 8164 gaussianRandom1: 8165 xg=0.; 8166 for(G4int i = 1; i <= 12; i++) { //do 6 j=1,12 8167 standardRandom(&xxg,&(hazard->ial)); 8168 // gaussianRandom6: 8169 xg=xg+xxg; 8170 } 8171 xg=xg-6.; 8172 if (xg*xg > 9.0) goto gaussianRandom1; 8173 standardRandom(&xshuf,&(hazard->igraine[10])); 8174 if (xshuf > 0.5) standardRandom(&xxg,&(hazard->ial)); 8175 // gaussianRandom2: 8176 (*rndm) = xg; 8177 return; 8178 /* 7313 8179 G4double tempRandom = 0.0, random = 0.0, randomShuffle = 0.0; 7314 8180 … … 7332 8198 7333 8199 (*rndm) = random; 8200 */ 7334 8201 } 7335 8202 … … 7346 8213 } 7347 8214 7348 G4double G4Incl::radius(G4 doubleA)8215 G4double G4Incl::radius(G4int A) 7349 8216 { 7350 8217 const G4double dp1 = 1.0, dp3 = 3.0; … … 7360 8227 7361 8228 G4double fact = std::sqrt(dp5/dp3); 7362 G4int ia = int(std::floor(A+0.4)); 8229 // G4int ia = int(std::floor(A+0.4)); 8230 G4int ia = A; 7363 8231 G4double result = fact * (0.84 * std::pow(A,dpth) + 0.55); 7364 8232 for(G4int i = 0; i < naSize; i++) { … … 7609 8477 } 7610 8478 8479 // C**************************************************************************** 8480 8481 void G4Incl::projo_spec(G4int, G4int ips, 8482 G4double, G4double, G4double) // max_a_proj is a global #define'd as a global constant 8483 { 8484 // subroutine projo_spec(ia1,ips,max_a_proj, 8485 // sn_projspec,fmpinc,pinc,tlab) 8486 // common/bl1/p1(300),p2(300),p3(300),eps(300),ind1(300),ind2(300),ta 8487 // spectators of the projectile treated here (nucleus for fermi-breakup): 8488 // ips is the number of spect, their label is in n_projspec(ips). 8489 // e,p are outside the target potential for these nucleons. 8490 // la transmission de max_a_proj=60 ne marche pas....explicite ici! 8491 // integer a_projspec,z_projspec,n_projspec(60) 8492 // real*4 m_projspec 8493 // common/proj_spect/a_projspec,z_projspec,ex_projspec,t_projspec, 8494 // s p1_projspec,p2_projspec,p3_projspec,m_projspec 8495 /// common/quatvectprojo/eps_c(60),p3_c(60), 8496 // s p1_s(60),p2_s(60),p3_s(60), 8497 // s t_c(60) 8498 // (eps_c,p1_s,p2_s,p3_s,eps_c used to store the kinematics of nucleons 8499 // for composit projectiles before entering the potential) 8500 // common/fermi/tf,pf,pf2 8501 8502 // G4double ttrou[21]; 8503 G4double fja = 0.0; 8504 G4double fjz = 0.0; 8505 // G4double ex2; 8506 // G4double rmmass; 8507 // const G4double um = 931.20793; 8508 8509 // c verifs: ground state 8510 // c es=0. 8511 // c p1s=0. 8512 // c p2s=0. 8513 // c p3s=0. 8514 // c gsmass=0. 8515 // c do iloc=1,ia1 8516 // c es=es+eps_c(iloc) 8517 // c p1s=p1s+p1_s(iloc) 8518 // c p2s=p2s+p2_s(iloc) 8519 // c p3s=p3s+p3_s(iloc) 8520 // c gsmass=gsmass+sqrt(eps_c(iloc)**2 8521 // c s -p1_s(iloc)**2-p2_s(iloc)**2-p3_s(iloc)**2) 8522 // c enddo 8523 // c sgs=sqrt(es**2-p1s**2-p2s**2-p3s**2) 8524 // c write(6,*) 'sqs, gsmass :',sgs,gsmass 8525 // c write(6,*) 'es,p1s,p2s,p3s:',es,p1s,p2s,p3s 8526 // c end verifs 8527 8528 // G4cout <<"proj_spect called with: " << G4endl; 8529 // G4cout <<"ips = " << ips << G4endl; 8530 ps->a_projspec=0; 8531 ps->z_projspec=0; 8532 ps->t_projspec=0.0; 8533 G4double e_spec=0.; 8534 G4double p1_spec=0.; 8535 G4double p2_spec=0.; 8536 G4double p3_spec=0.; 8537 G4double gs_mass=0.; 8538 8539 if(ips > 1) { 8540 for(G4int i_spec=1; i_spec <= ips; i_spec++) { 8541 G4int i_c = ps->n_projspec[i_spec]; 8542 e_spec = e_spec+bl1->eps[i_c]; 8543 p1_spec = p1_spec+bl1->p1[i_c]; 8544 p2_spec = p2_spec+bl1->p2[i_c]; 8545 p3_spec = p3_spec+bl1->p3[i_c]; 8546 G4double pi_spec2= std::pow(bl1->p1[i_c],2) + std::pow(bl1->p2[i_c],2) + std::pow(bl1->p3[i_c],2); 8547 // current ground state mass: 8548 gs_mass=gs_mass+std::sqrt(std::pow(bl1->eps[i_c],2)-pi_spec2); 8549 if(bl1->ind2[i_c] == 1) ps->z_projspec=ps->z_projspec+1; 8550 } // enddo 8551 G4double p_spec2=std::pow(p1_spec,2)+std::pow(p2_spec,2)+std::pow(p3_spec,2); 8552 G4double s_spec = sqrt(std::pow(e_spec,2)-p_spec2); 8553 8554 // no projectile spectator if a>=4 and a=z or a=n (no dresner breakup) 8555 if(ips >= 4 && (ps->z_projspec == ips || ps->z_projspec == 0)) { 8556 for(G4int i_spec = 1; i_spec <= ips; i_spec++) { //do i_spec=1,ips 8557 G4int i_c = ps->n_projspec[i_spec]; 8558 // single nucleon forced on shell! momentum kept. 8559 // single nucleon forced on shell! energy kept. 8560 G4double arg=(std::pow(bl1->eps[i_c],2)-std::pow(938.2796,2))/ 8561 (std::pow(bl1->p1[i_c],2)+std::pow(bl1->p2[i_c], 2)+std::pow(bl1->p3[i_c],2)); 8562 G4double coef=0.0; 8563 if(arg < 0.) { // Single nucleon is forced on shell. We have 8564 // to choose between energy and momentum 8565 // conservation... so we choose energy: 8566 coef = 1.0; 8567 bl1->eps[i_c]=std::sqrt(std::pow(bl1->p1[i_c],2) + std::pow(bl1->p2[i_c], 2) + 8568 std::pow(bl1->p3[i_c], 2) + std::pow(938.2796,2)); 8569 } else { 8570 coef=std::sqrt(arg); 8571 bl1->p1[i_c]=bl1->p1[i_c]*coef; 8572 bl1->p2[i_c]=bl1->p2[i_c]*coef; 8573 bl1->p3[i_c]=bl1->p3[i_c]*coef; 8574 } 8575 } // enddo 8576 return; 8577 } // endif 8578 8579 // excitation energy from particle-hole excitation: 8580 // seconde methode de calcul ex4: 8581 for(G4int i = 1; i <= bl3->ia1; i++) { // do i=1,ia1 8582 G4int nb_p=i; 8583 // call ordered(t_c(i),nb_p,tpart) 8584 ordered(qvp->t_c[i], nb_p); 8585 } // enddo 8586 G4double ex4=0.; 8587 for(G4int i = 1; i <= ips; i++) { //do i=1,ips 8588 G4int j=ps->n_projspec[i]; 8589 ex4=ex4+qvp->t_c[j]-ps->tab[i]; //ex4=ex4+qvp->t_c[j]-ps->tpart[i]; 8590 } // enddo 8591 8592 ps->ex_projspec =ex4; 8593 8594 // coherent rest mass: 8595 ps->m_projspec = s_spec - ps->ex_projspec; 8596 8597 ps->a_projspec=ips; 8598 8599 if((ps->ex_projspec < 0.) || (ps->m_projspec < 0.)) { 8600 G4cout <<"ex_projspec is negative,a,z: " << ps->ex_projspec << " " << fja << " " << fjz << G4endl; 8601 G4cout <<"...or m is negative,a,z:" << ps->m_projspec << " " << fja << " " << fjz << G4endl; 8602 } // endif 8603 8604 ps->t_projspec = e_spec-s_spec; 8605 ps->p1_projspec = p1_spec; 8606 ps->p2_projspec = p2_spec; 8607 ps->p3_projspec = p3_spec; 8608 8609 // c check module... 8610 // c write(6,*) 'projectile spectator a,z:',a_projspec,z_projspec 8611 // c write(6,*) 'm,t,e*:',m_projspec,t_projspec,ex_projspec 8612 // c write(6,*) 'p1,p2,p3:',p1_projspec,p2_projspec,p3_projspec 8613 // c write(6,*) ' constituants:',ips 8614 // c do i_spec=1,ips 8615 // c i_c = n_projspec(i_spec) 8616 // c write(6,*) i_c,eps(i_c),p1(i_c),p2(i_c),p3(i_c) 8617 // c enddo 8618 // c ...end check module 8619 8620 return; 8621 } else if(ips == 1) { // then 8622 ps->a_projspec = 0; // number of nucleons in the spectator nucleus 8623 G4int i_c = ps->n_projspec[1]; 8624 //c single nucleon forced on shell! momentum kept. 8625 //c eps(i_c)=sqrt(p1(i_c)**2+p2(i_c)**2+p3(i_c)**2+938.2796**2) 8626 //c single nucleon forced on shell! energy kept. 8627 G4double arg=(std::pow(bl1->eps[i_c],2)-std::pow(938.2796,2))/ 8628 (std::pow(bl1->p1[i_c],2)+std::pow(bl1->p2[i_c],2)+std::pow(bl1->p3[i_c],2)); 8629 G4double coef=0.0; // Initialize the coefficieant 8630 if(arg < 0.) { 8631 G4cout <<"off shell problems in projo_spec, forced coef = 1.0 to recover from the error" << G4endl; 8632 coef = 1.0; 8633 } else { 8634 coef = std::sqrt(arg); 8635 } 8636 8637 bl1->p1[i_c]=bl1->p1[i_c]*coef; 8638 bl1->p2[i_c]=bl1->p2[i_c]*coef; 8639 bl1->p3[i_c]=bl1->p3[i_c]*coef; 8640 } //endif 8641 8642 ps->a_projspec = 0; // number of nucleons in the spectator nucleus 8643 return; 8644 } 8645 8646 void G4Incl::ordered(G4double t, G4int nb) 8647 { 8648 // SUBROUTINE ORDERED(T,nb,TAB) 8649 //C This subroutine will ordered in TAB() from min TAB(1) to max 8650 //C TAB(nb), nb values of T sent successively. 8651 //C CALL ORDERED (T1,1), CALL ORDERED (T2,2),..CALL ORDERED (Tnb,nb) 8652 // DIMENSION TAB(1) 8653 8654 // G4cout <<"Ordered was called" << G4endl; 8655 if(nb == 1) { 8656 ps->tab[1]=t; 8657 return; 8658 } // endif 8659 8660 G4int i=0; 8661 8662 ordered100: 8663 i=i+1; 8664 if(t < ps->tab[i]) { 8665 // for(G4int l = nb-1; l != i; l--) { //do l=nb-1,i,-1 8666 for(G4int l = nb-1; l >= i; l--) { //do l=nb-1,i,-1 8667 // G4cout <<"Ordered: L = " << l << " nb = " << nb << " i = " << i << G4endl; 8668 ps->tab[l+1]=ps->tab[l]; 8669 } // enddo 8670 ps->tab[i]=t; 8671 return; 8672 } else { 8673 if(i <= nb-1) goto ordered100; 8674 ps->tab[i]=t; 8675 return; 8676 } // endif 8677 return; 8678 } 8679 7611 8680 // Utilities 7612 8681 … … 7903 8972 G4cout <<"(particle " << index << G4endl; 7904 8973 if(bl1->ind1[index] == 0) { // Nucleon 7905 if(bl1->ind1[index] == 1) G4cout << "(quote proton)" << G4endl; 7906 else if(bl1->ind1[index] == -1) G4cout << "(quote neutron)" << G4endl; 7907 else G4cout << "(quote unidentified-nucleon)" << G4endl; 8974 if(bl1->ind2[index] == 1) { 8975 G4cout << "(quote proton)" << G4endl; 8976 } else if(bl1->ind2[index] == -1) { 8977 G4cout << "(quote neutron)" << G4endl; 8978 } else { 8979 G4cout << "(quote unidentified-nucleon-ind1-" << bl1->ind1[index] <<"-ind2-" << bl1->ind2[index] << G4endl; 8980 } 7908 8981 } else if(bl1->ind1[index] == 1) { // Delta 7909 if(bl1->ind2[index] >= -2 && bl1->ind2[index] <= 2) G4cout << "(quote delta-" 7910 << bl1->ind2[index] 7911 << ")" << G4endl; 7912 else G4cout << "(quote unidentified-delta)" << G4endl; 8982 if(bl1->ind2[index] >= -2 && bl1->ind2[index] <= 2) { 8983 G4cout << "(quote delta-" 8984 << bl1->ind2[index] 8985 << ")" << G4endl; 8986 } else { 8987 G4cout << "(quote unidentified-delta)" << G4endl; 8988 } 7913 8989 } else { 7914 8990 G4cout <<"(quote unidentified-particle)" << G4endl; … … 7938 9014 G4cout <<"(list ;; Particles at the beginning of the avatar" << G4endl; 7939 9015 G4int ia = bl3->ia1 + bl3->ia2; 7940 for(G4int i = 1; i != ia; ++i) { // do i=1,IA9016 for(G4int i = 1; i <= ia; ++i) { // do i=1,IA 7941 9017 print_one_particle(i); 7942 9018 } -
trunk/source/processes/hadronic/models/incl/src/G4InclAblaCascadeInterface.cc
r1228 r1340 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4InclAblaCascadeInterface.cc,v 1.1 3 2009/12/04 13:16:57 kaitanieExp $26 // $Id: G4InclAblaCascadeInterface.cc,v 1.16 2010/10/29 06:48:43 gunter Exp $ 27 27 // Translation of INCL4.2/ABLA V3 28 28 // Pekka Kaitaniemi, HIP (translation) … … 34 34 35 35 #include "G4InclAblaCascadeInterface.hh" 36 #include "G4FermiBreakUp.hh" 36 37 #include "math.h" 37 38 #include "G4GenericIon.hh" … … 46 47 47 48 varntp = new G4VarNtp(); 48 calincl = new G4Calincl();49 calincl = 0; 49 50 ws = new G4Ws(); 50 51 mat = new G4Mat(); 51 52 incl = new G4Incl(hazard, calincl, ws, mat, varntp); 53 if(!getenv("G4INCLABLANOFERMIBREAKUP")) { // Use Fermi Break-up by default if it is NOT explicitly disabled 54 incl->setUseFermiBreakUp(true); 55 } 52 56 53 57 verboseLevel = 0; … … 58 62 delete hazard; 59 63 delete varntp; 60 delete calincl;61 64 delete ws; 62 65 delete mat; … … 68 71 G4int maxTries = 200; 69 72 70 G4int particleI, n = 0; 71 73 G4int particleI; 72 74 G4int bulletType = 0; 73 75 … … 86 88 } 87 89 88 // INCL4 needs the energy in units MeV89 G4double bulletE = aTrack.GetKineticEnergy() * MeV;90 91 90 #ifdef DEBUGINCL 92 91 G4cout <<"Bullet energy = " << bulletE / MeV << G4endl; 93 92 #endif 94 95 G4double targetA = theNucleus.GetN();96 G4double targetZ = theNucleus.GetZ();97 93 98 94 G4double eKin; … … 100 96 G4DynamicParticle *cascadeParticle = 0; 101 97 G4ParticleDefinition *aParticleDefinition = 0; 98 99 G4FermiBreakUp *fermiBreakUp = new G4FermiBreakUp(); 100 G4FragmentVector *theFermiBreakupResult = 0; 101 G4ParticleTable *theTableOfParticles = G4ParticleTable::GetParticleTable(); 102 102 103 103 // INCL assumes the projectile particle is going in the direction of … … 113 113 theResult.Clear(); // Make sure the output data structure is clean. 114 114 115 // Map Geant4 particle types to corresponding INCL4 types. 116 enum bulletParticleType {nucleus = 0, proton = 1, neutron = 2, pionPlus = 3, pionZero = 4, 117 pionMinus = 5, deuteron = 6, triton = 7, he3 = 8, he4 = 9}; 118 119 // Coding particles for use with INCL4 and ABLA 120 if (aTrack.GetDefinition() == G4Proton::Proton() ) bulletType = proton; 121 if (aTrack.GetDefinition() == G4Neutron::Neutron() ) bulletType = neutron; 122 if (aTrack.GetDefinition() == G4PionPlus::PionPlus() ) bulletType = pionPlus; 123 if (aTrack.GetDefinition() == G4PionMinus::PionMinus() ) bulletType = pionMinus; 124 if (aTrack.GetDefinition() == G4PionZero::PionZero() ) bulletType = pionZero; 115 calincl = new G4InclInput(aTrack, theNucleus, false); 116 incl->setInput(calincl); 117 118 G4InclInput::printProjectileTargetInfo(aTrack, theNucleus); 119 calincl->printInfo(); 125 120 126 121 #ifdef DEBUGINCL … … 139 134 #endif 140 135 141 for(int i = 0; i < 15; i++) {142 calincl->f[i] = 0.0; // Initialize INCL input data143 }144 145 136 // Check wheter the input is acceptable. 146 if((bulletType != 0) && ((targetA != 1) && (targetZ != 1))) { 147 calincl->f[0] = targetA; // Target mass number 148 calincl->f[1] = targetZ; // Charge number 149 calincl->f[6] = bulletType; // Type 150 calincl->f[2] = bulletE; // Energy [MeV] 151 calincl->f[5] = 1.0; // Time scaling 152 calincl->f[4] = 45.0; // Nuclear potential 153 137 if((calincl->bulletType() != 0) && ((calincl->targetA() != 1) && (calincl->targetZ() != 1))) { 154 138 ws->nosurf = -2; // Nucleus surface, -2 = Woods-Saxon 155 139 ws->xfoisa = 8; … … 160 144 161 145 mat->nbmat = 1; 162 mat->amat[0] = int(calincl-> f[0]);163 mat->zmat[0] = int(calincl-> f[1]);146 mat->amat[0] = int(calincl->targetA()); 147 mat->zmat[0] = int(calincl->targetZ()); 164 148 165 149 incl->initIncl(true); … … 170 154 G4cout <<"G4InclAblaCascadeInterface: Try number = " << nTries << G4endl; 171 155 } 172 incl->processEventInclAbla( eventNumber);156 incl->processEventInclAbla(calincl, eventNumber); 173 157 174 158 if(verboseLevel > 1) { … … 181 165 * Diagnostic output 182 166 */ 183 G4cout <<"G4InclAblaCascadeInterface: Bullet type: " << bulletType<< G4endl;184 G4cout <<"G4Incl4AblaCascadeInterface: Bullet energy: " << bulletE<< " MeV" << G4endl;185 186 G4cout <<"G4InclAblaCascadeInterface: Target A: " << targetA<< G4endl;187 G4cout <<"G4InclAblaCascadeInterface: Target Z: " << targetZ<< G4endl;167 G4cout <<"G4InclAblaCascadeInterface: Bullet type: " << calincl->bulletType() << G4endl; 168 G4cout <<"G4Incl4AblaCascadeInterface: Bullet energy: " << calincl->bulletE() << " MeV" << G4endl; 169 170 G4cout <<"G4InclAblaCascadeInterface: Target A: " << calincl->targetA() << G4endl; 171 G4cout <<"G4InclAblaCascadeInterface: Target Z: " << calincl->targetZ() << G4endl; 188 172 189 173 if(verboseLevel > 3) { 190 diagdata <<"G4InclAblaCascadeInterface: Bullet type: " << bulletType<< G4endl;191 diagdata <<"G4InclAblaCascadeInterface: Bullet energy: " << bulletE<< " MeV" << G4endl;174 diagdata <<"G4InclAblaCascadeInterface: Bullet type: " << calincl->bulletType() << G4endl; 175 diagdata <<"G4InclAblaCascadeInterface: Bullet energy: " << calincl->bulletE() << " MeV" << G4endl; 192 176 193 diagdata <<"G4InclAblaCascadeInterface: Target A: " << targetA << G4endl; 194 diagdata <<"G4InclAblaCascadeInterface: Target Z: " << targetZ << G4endl; 195 } 196 197 for(particleI = 0; particleI < varntp->ntrack; particleI++) { 198 G4cout << n << " " << calincl->f[6] << " " << calincl->f[2] << " "; 199 G4cout << varntp->massini << " " << varntp->mzini << " "; 200 G4cout << varntp->exini << " " << varntp->mulncasc << " " << varntp->mulnevap << " " << varntp->mulntot << " "; 201 G4cout << varntp->bimpact << " " << varntp->jremn << " " << varntp->kfis << " " << varntp->estfis << " "; 202 G4cout << varntp->izfis << " " << varntp->iafis << " " << varntp->ntrack << " " << varntp->itypcasc[particleI] << " "; 203 G4cout << varntp->avv[particleI] << " " << varntp->zvv[particleI] << " " << varntp->enerj[particleI] << " "; 204 G4cout << varntp->plab[particleI] << " " << varntp->tetlab[particleI] << " " << varntp->philab[particleI] << G4endl; 205 // For diagnostic output 206 if(verboseLevel > 3) { 207 diagdata << n << " " << calincl->f[6] << " " << calincl->f[2] << " "; 208 diagdata << varntp->massini << " " << varntp->mzini << " "; 209 diagdata << varntp->exini << " " << varntp->mulncasc << " " << varntp->mulnevap << " " << varntp->mulntot << " "; 210 diagdata << varntp->bimpact << " " << varntp->jremn << " " << varntp->kfis << " " << varntp->estfis << " "; 211 diagdata << varntp->izfis << " " << varntp->iafis << " " << varntp->ntrack << " "; 212 diagdata << varntp->itypcasc[particleI] << " "; 213 diagdata << varntp->avv[particleI] << " " << varntp->zvv[particleI] << " " << varntp->enerj[particleI] << " "; 214 diagdata << varntp->plab[particleI] << " " << varntp->tetlab[particleI] << " " << varntp->philab[particleI] << G4endl; 215 } 216 } 177 diagdata <<"G4InclAblaCascadeInterface: Target A: " << calincl->targetA() << G4endl; 178 diagdata <<"G4InclAblaCascadeInterface: Target Z: " << calincl->targetZ() << G4endl; 179 } 180 217 181 } 218 182 … … 227 191 theResult.SetStatusChange(stopAndKill); 228 192 229 if(bulletType == proton) { 230 aParticleDefinition = G4Proton::ProtonDefinition(); 231 } else if(bulletType == neutron) { 232 aParticleDefinition = G4Neutron::NeutronDefinition(); 233 } else if(bulletType == pionPlus) { 234 aParticleDefinition = G4PionPlus::PionPlusDefinition(); 235 } else if(bulletType == pionZero) { 236 aParticleDefinition = G4PionZero::PionZeroDefinition(); 237 } else if(bulletType == pionMinus) { 238 aParticleDefinition = G4PionMinus::PionMinusDefinition(); 239 } else { // Projectile was not regognized 240 aParticleDefinition = 0; 241 } 193 G4int bulletType = calincl->bulletType(); 194 aParticleDefinition = G4InclInput::getParticleDefinition(bulletType); 242 195 243 196 if(aParticleDefinition != 0) { … … 328 281 if((varntp->avv[particleI] > 1) && (varntp->zvv[particleI] >= 1)) { // Nucleus fragment 329 282 G4ParticleDefinition * aIonDef = 0; 330 G4ParticleTable *theTableOfParticles = G4ParticleTable::GetParticleTable();331 283 332 284 G4int A = G4int(varntp->avv[particleI]); … … 404 356 } 405 357 } 358 } 359 360 // Finally do Fermi break-up if needed 361 if(varntp->needsFermiBreakup) { 362 // baryonNumberBalanceInINCL -= varntp->massini; 363 // chargeNumberBalanceInINCL -= varntp->mzini; 364 // Call Fermi Break-up 365 G4double nuclearMass = G4NucleiProperties::GetNuclearMass(G4int(varntp->massini), G4int(varntp->mzini)) + varntp->exini * MeV; 366 G4LorentzVector fragmentMomentum(varntp->pxrem * MeV, varntp->pyrem * MeV, varntp->pzrem * MeV, 367 varntp->erecrem * MeV + nuclearMass); 368 G4double momentumScaling = G4InclUtils::calculate4MomentumScaling(G4int(varntp->massini), G4int(varntp->mzini), 369 varntp->exini, 370 varntp->erecrem, 371 varntp->pxrem, 372 varntp->pyrem, 373 varntp->pzrem); 374 G4LorentzVector p4(momentumScaling * varntp->pxrem * MeV, momentumScaling * varntp->pyrem * MeV, 375 momentumScaling * varntp->pzrem * MeV, 376 varntp->erecrem + nuclearMass); 377 378 // For four-momentum, baryon number and charge conservation check: 379 G4LorentzVector fourMomentumBalance = p4; 380 G4int baryonNumberBalance = G4int(varntp->massini); 381 G4int chargeBalance = G4int(varntp->mzini); 382 383 G4LorentzRotation toFragmentZ; 384 toFragmentZ.rotateZ(-p4.theta()); 385 toFragmentZ.rotateY(-p4.phi()); 386 G4LorentzRotation toFragmentLab = toFragmentZ.inverse(); 387 // p4 *= toFragmentZ; 388 389 G4LorentzVector p4rest = p4; 390 // p4rest.boost(-p4.boostVector()); 391 if(verboseLevel > 0) { 392 G4cout <<"Cascade remnant nucleus:" << G4endl; 393 G4cout <<"p4: " << G4endl; 394 G4cout <<" px: " << p4.px() <<" py: " << p4.py() <<" pz: " << p4.pz() << G4endl; 395 G4cout <<" E = " << p4.e() << G4endl; 396 397 G4cout <<"p4rest: " << G4endl; 398 G4cout <<" px: " << p4rest.px() <<" py: " << p4rest.py() <<" pz: " << p4rest.pz() << G4endl; 399 G4cout <<" E = " << p4rest.e() << G4endl; 400 } 401 402 G4Fragment theCascadeRemnant(G4int(varntp->massini), G4int(varntp->mzini), p4rest); 403 theFermiBreakupResult = fermiBreakUp->BreakItUp(theCascadeRemnant); 404 if(theFermiBreakupResult != 0) { 405 G4FragmentVector::iterator fragment; 406 for(fragment = theFermiBreakupResult->begin(); fragment != theFermiBreakupResult->end(); fragment++) { 407 G4ParticleDefinition *theFragmentDefinition = 0; 408 if((*fragment)->GetA_asInt() == 1 && (*fragment)->GetZ_asInt() == 0) { // Neutron 409 theFragmentDefinition = G4Neutron::NeutronDefinition(); 410 } else if ((*fragment)->GetA_asInt() == 1 && (*fragment)->GetZ_asInt() == 1) { 411 theFragmentDefinition = G4Proton::ProtonDefinition(); 412 } else { 413 theFragmentDefinition = theTableOfParticles->GetIon((*fragment)->GetZ_asInt(), (*fragment)->GetA_asInt(), (*fragment)->GetExcitationEnergy()); 414 } 415 416 if(theFragmentDefinition != 0) { 417 G4DynamicParticle *theFragment = new G4DynamicParticle(theFragmentDefinition, (*fragment)->GetMomentum()); 418 G4LorentzVector labMomentum = theFragment->Get4Momentum(); 419 // labMomentum.boost(p4.boostVector()); 420 // labMomentum *= toFragmentLab; 421 // labMomentum *= toLabFrame; 422 theFragment->Set4Momentum(labMomentum); 423 fourMomentumBalance -= theFragment->Get4Momentum(); 424 baryonNumberBalance -= theFragmentDefinition->GetAtomicMass(); 425 chargeBalance -= theFragmentDefinition->GetAtomicNumber(); 426 if(verboseLevel > 0) { 427 G4cout <<"Resulting fragment: " << G4endl; 428 G4cout <<" kinetic energy = " << theFragment->GetKineticEnergy() / MeV << " MeV" << G4endl; 429 G4cout <<" momentum = " << theFragment->GetMomentum().mag() / MeV << " MeV" << G4endl; 430 } 431 theResult.AddSecondary(theFragment); 432 } else { 433 G4cout <<"G4InclAblaCascadeInterface: Error. Fragment produced by Fermi break-up does not exist." << G4endl; 434 G4cout <<"Resulting fragment: " << G4endl; 435 G4cout <<" Z = " << (*fragment)->GetZ_asInt() << G4endl; 436 G4cout <<" A = " << (*fragment)->GetA_asInt() << G4endl; 437 G4cout <<" Excitation : " << (*fragment)->GetExcitationEnergy() / MeV << " MeV" << G4endl; 438 G4cout <<" momentum = " << (*fragment)->GetMomentum().mag() / MeV << " MeV" << G4endl; 439 } 440 } 441 if(std::abs(fourMomentumBalance.mag() / MeV) > 0.1 * MeV) { 442 G4cout <<"Four-momentum balance after remnant nucleus Fermi break-up:" << G4endl; 443 G4cout <<"Magnitude: " << fourMomentumBalance.mag() / MeV << " MeV" << G4endl; 444 G4cout <<"Vector components (px, py, pz, E) = (" 445 << fourMomentumBalance.px() << ", " 446 << fourMomentumBalance.py() << ", " 447 << fourMomentumBalance.pz() << ", " 448 << fourMomentumBalance.e() << ")" << G4endl; 449 } 450 if(baryonNumberBalance != 0) { 451 G4cout <<"Baryon number balance after remnant nucleus Fermi break-up: " << baryonNumberBalance << G4endl; 452 } 453 if(chargeBalance != 0) { 454 G4cout <<"Charge balance after remnant nucleus Fermi break-up: " << chargeBalance << G4endl; 455 } 456 } 406 457 } 407 458 … … 449 500 theResult.SetStatusChange(stopAndKill); 450 501 451 G4ParticleTable *theTableOfParticles = G4ParticleTable::GetParticleTable();452 502 cascadeParticle = new G4DynamicParticle(theTableOfParticles->FindParticle(aTrack.GetDefinition()), aTrack.Get4Momentum()); 453 503 … … 472 522 } 473 523 474 if(( targetA == 1) && (targetZ== 1)) { // Unsupported target524 if((calincl->targetA() == 1) && (calincl->targetZ() == 1)) { // Unsupported target 475 525 if(verboseLevel > 1) { 476 526 G4cout <<"Unsupported target: " << G4endl; 477 G4cout <<"Target A: " << targetA<< G4endl;478 G4cout <<"TargetZ: " << targetZ<< G4endl;527 G4cout <<"Target A: " << calincl->targetA() << G4endl; 528 G4cout <<"TargetZ: " << calincl->targetZ() << G4endl; 479 529 } 480 530 if(verboseLevel > 3) { 481 531 diagdata <<"Unsupported target: " << G4endl; 482 diagdata <<"Target A: " << targetA<< G4endl;483 diagdata <<"TargetZ: " << targetZ<< G4endl;484 } 485 } 486 487 if( bulletE< 100) { // INCL does not support E < 100 MeV.532 diagdata <<"Target A: " << calincl->targetA() << G4endl; 533 diagdata <<"TargetZ: " << calincl->targetZ() << G4endl; 534 } 535 } 536 537 if(calincl->bulletE() < 100) { // INCL does not support E < 100 MeV. 488 538 if(verboseLevel > 1) { 489 G4cout <<"Unsupported bullet energy: " << bulletE<< " MeV. (Lower limit is 100 MeV)." << G4endl;539 G4cout <<"Unsupported bullet energy: " << calincl->bulletE() << " MeV. (Lower limit is 100 MeV)." << G4endl; 490 540 G4cout <<"WARNING: Returning the original bullet with original energy back to Geant4." << G4endl; 491 541 } 492 542 if(verboseLevel > 3) { 493 diagdata <<"Unsupported bullet energy: " << bulletE<< " MeV. (Lower limit is 100 MeV)." << G4endl;543 diagdata <<"Unsupported bullet energy: " << calincl->bulletE() << " MeV. (Lower limit is 100 MeV)." << G4endl; 494 544 } 495 545 } … … 500 550 } 501 551 552 delete fermiBreakUp; 553 delete calincl; 554 calincl = 0; 502 555 return &theResult; 503 556 } -
trunk/source/processes/hadronic/models/incl/src/G4InclAblaLightIonInterface.cc
r1228 r1340 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4InclAblaLightIonInterface.cc,v 1.1 1 2009/12/04 13:16:57 kaitanieExp $26 // $Id: G4InclAblaLightIonInterface.cc,v 1.14 2010/10/29 06:48:43 gunter Exp $ 27 27 // Translation of INCL4.2/ABLA V3 28 28 // Pekka Kaitaniemi, HIP (translation) … … 31 31 // Aatos Heikkinen, HIP (project coordination) 32 32 33 #include <vector> 34 33 35 #include "G4InclAblaLightIonInterface.hh" 36 #include "G4FermiBreakUp.hh" 34 37 #include "math.h" 35 38 #include "G4GenericIon.hh" … … 44 47 45 48 varntp = new G4VarNtp(); 46 calincl = new G4Calincl();49 calincl = 0; 47 50 ws = new G4Ws(); 48 51 mat = new G4Mat(); 49 52 incl = new G4Incl(hazard, calincl, ws, mat, varntp); 50 53 useProjectileSpectator = true; 54 useFermiBreakup = true; 55 incl->setUseProjectileSpectators(useProjectileSpectator); 56 if(!getenv("G4INCLABLANOFERMIBREAKUP")) { // Use Fermi Break-up by default if it is NOT explicitly disabled 57 incl->setUseFermiBreakUp(true); 58 useFermiBreakup = true; 59 } 51 60 verboseLevel = 0; 61 if(getenv("G4INCLVERBOSE")) { 62 verboseLevel = 1; 63 } 52 64 } 53 65 … … 64 76 G4HadFinalState* G4InclAblaLightIonInterface::ApplyYourself(const G4HadProjectile& aTrack, G4Nucleus& theNucleus) 65 77 { 78 // const G4bool useFermiBreakup = false; 66 79 G4int maxTries = 200; 67 80 68 G4int particleI , n = 0;69 70 G4int b ulletType= 0;71 72 // Print diagnostic messages: 0 = silent, 1 and 2 = verbose 73 verboseLevel = 0;81 G4int particleI; 82 83 G4int baryonNumberBalanceInINCL = 0; 84 G4int chargeNumberBalanceInINCL = 0; 85 86 G4ParticleTable *theTableOfParticles = G4ParticleTable::GetParticleTable(); 74 87 75 88 // Increase the event number: 76 89 eventNumber++; 77 90 91 // Clean up the INCL input 92 if(calincl != 0) { 93 delete calincl; 94 calincl = 0; 95 } 96 78 97 if (verboseLevel > 1) { 79 98 G4cout << " >>> G4InclAblaLightIonInterface::ApplyYourself called" << G4endl; … … 84 103 } 85 104 86 // INCL4 needs the energy in units MeV 87 G4double bulletE = aTrack.GetKineticEnergy() / MeV; 88 89 G4double targetA = theNucleus.GetN(); 90 G4double targetZ = theNucleus.GetZ(); 105 if(verboseLevel > 1) { 106 G4Calincl::printProjectileTargetInfo(aTrack, theNucleus); 107 } 108 109 // Inverse kinematics for targets with Z = 1 and A = 1 110 // if(false) { 111 G4LorentzRotation toBreit = aTrack.Get4Momentum().boostVector(); 112 113 if(theNucleus.GetZ_asInt() == 1 && theNucleus.GetA_asInt() == 1 && G4Calincl::canUseInverseKinematics(aTrack, theNucleus)) { 114 G4ParticleDefinition *oldTargetDef = theTableOfParticles->GetIon(theNucleus.GetA_asInt(), theNucleus.GetZ_asInt(), 0.0); 115 const G4ParticleDefinition *oldProjectileDef = aTrack.GetDefinition(); 116 117 G4int oldTargetA = oldTargetDef->GetAtomicMass(); 118 G4int newTargetA = oldProjectileDef->GetAtomicMass(); 119 G4int newTargetZ = oldProjectileDef->GetAtomicNumber(); 120 121 if(newTargetA > 0 && newTargetZ > 0 && oldTargetDef != 0 && oldProjectileDef != 0) { 122 G4Nucleus swappedTarget(oldProjectileDef->GetAtomicMass(), oldProjectileDef->GetAtomicNumber()); 123 124 // G4cout <<"Original projectile kinE = " << aTrack.GetKineticEnergy() / MeV << G4endl; 125 126 // We need the same energy/nucleon. 127 G4double projectileE = ((aTrack.GetKineticEnergy() / MeV) / newTargetA) * oldTargetA * MeV; 128 129 // G4cout <<"projectileE = " << projectileE << G4endl; 130 G4DynamicParticle swappedProjectileParticle(oldTargetDef, G4ThreeVector(0.0, 0.0, 1.0), projectileE); 131 const G4LorentzVector swapped4Momentum = (swappedProjectileParticle.Get4Momentum()*=toBreit); 132 swappedProjectileParticle.Set4Momentum(swapped4Momentum); 133 const G4HadProjectile swappedProjectile(swappedProjectileParticle); 134 // G4cout <<"New projectile kinE = " << swappedProjectile.GetKineticEnergy() / MeV << G4endl; 135 calincl = new G4InclInput(swappedProjectile, swappedTarget, true); 136 } else { 137 G4cout <<"Badly defined target after swapping. Falling back to normal (non-swapped) mode." << G4endl; 138 calincl = new G4InclInput(aTrack, theNucleus, false); 139 } 140 } else { 141 calincl = new G4InclInput(aTrack, theNucleus, false); 142 } 91 143 92 144 G4double eKin; … … 105 157 G4LorentzRotation toLabFrame = toZ.inverse(); 106 158 159 /* 160 G4cout <<"Projectile theta = " << projectileMomentum.theta() << " phi = " << projectileMomentum.phi() << G4endl; 161 G4cout <<"Projectile momentum " 162 << "(px = " << projectileMomentum.px() 163 << ", py = " << projectileMomentum.py() 164 << ", pz = " << projectileMomentum.pz() << ")" << G4endl; 165 G4cout << "Projectile energy = " << bulletE << " MeV" << G4endl; 166 */ 167 168 G4FermiBreakUp *fermiBreakUp = new G4FermiBreakUp(); 169 G4FragmentVector *theSpectatorFermiBreakupResult = 0; 170 G4FragmentVector *theFermiBreakupResult = 0; 171 107 172 theResult.Clear(); // Make sure the output data structure is clean. 173 174 std::vector<G4DynamicParticle*> result; // Temporary list for the results 108 175 109 176 // Map Geant4 particle types to corresponding INCL4 types. 110 177 enum bulletParticleType {nucleus = 0, proton = 1, neutron = 2, pionPlus = 3, pionZero = 4, 111 pionMinus = 5, deuteron = 6, triton = 7, he3 = 8, he4 = 9}; 112 113 // Coding particles for use with INCL4 and ABLA 114 if (aTrack.GetDefinition() == G4Deuteron::Deuteron() ) bulletType = deuteron; 115 if (aTrack.GetDefinition() == G4Triton::Triton() ) bulletType = triton; 116 if (aTrack.GetDefinition() == G4He3::He3() ) bulletType = he3; 117 if (aTrack.GetDefinition() == G4Alpha::Alpha() ) bulletType = he4; 118 119 for(int i = 0; i < 15; i++) { 120 calincl->f[i] = 0.0; // Initialize INCL input data 178 pionMinus = 5, deuteron = 6, triton = 7, he3 = 8, he4 = 9, 179 c12 = -12}; // Carbon beam support. 180 181 G4int bulletType = calincl->bulletType(); 182 chargeNumberBalanceInINCL = calincl->targetZ(); 183 baryonNumberBalanceInINCL = calincl->targetA(); 184 185 // G4cout <<"Type of the projectile (INCL projectile code): " << bulletType << G4endl; 186 187 if(bulletType == proton) { 188 chargeNumberBalanceInINCL += 1; 189 baryonNumberBalanceInINCL += 1; 190 } else if(bulletType == neutron) { 191 baryonNumberBalanceInINCL += 1; 192 } else if(bulletType == pionPlus) { //Note: positive pion doesn't contribute to the baryon and charge number counters 193 chargeNumberBalanceInINCL += 1; 194 } else if(bulletType == pionMinus) { 195 chargeNumberBalanceInINCL -= 1; 196 } else if(bulletType == deuteron) { 197 chargeNumberBalanceInINCL += 1; 198 baryonNumberBalanceInINCL += 2; 199 } else if(bulletType == triton) { 200 chargeNumberBalanceInINCL += 1; 201 baryonNumberBalanceInINCL += 3; 202 } else if(bulletType == he3) { 203 chargeNumberBalanceInINCL += 2; 204 baryonNumberBalanceInINCL += 3; 205 } else if(bulletType == he4) { 206 chargeNumberBalanceInINCL += 2; 207 baryonNumberBalanceInINCL += 4; 208 } if(bulletType == c12) { 209 chargeNumberBalanceInINCL += 6; 210 baryonNumberBalanceInINCL += 12; 211 } if(bulletType == -666) { 212 chargeNumberBalanceInINCL += calincl->extendedProjectileZ(); 213 baryonNumberBalanceInINCL += calincl->extendedProjectileA(); 121 214 } 122 215 123 216 // Check wheter the input is acceptable. 124 if((bulletType != 0) && ((targetA != 1) && (targetZ != 1))) { 125 calincl->f[0] = targetA; // Target mass number 126 calincl->f[1] = targetZ; // Charge number 127 calincl->f[6] = bulletType; // Type 128 calincl->f[2] = bulletE; // Energy [MeV] 129 calincl->f[5] = 1.0; // Time scaling 130 calincl->f[4] = 45.0; // Nuclear potential 131 217 if((bulletType != 0) && ((calincl->targetA() != 1) && (calincl->targetZ() != 1))) { 132 218 ws->nosurf = -2; // Nucleus surface, -2 = Woods-Saxon 133 219 ws->xfoisa = 8; … … 138 224 139 225 mat->nbmat = 1; 140 mat->amat[0] = int(calincl->f[0]); 141 mat->zmat[0] = int(calincl->f[1]); 142 226 mat->amat[0] = int(calincl->targetA()); 227 mat->zmat[0] = int(calincl->targetA()); 228 229 incl->setInput(calincl); 143 230 incl->initIncl(true); 144 231 … … 148 235 G4cout <<"G4InclAblaLightIonInterface: Try number = " << nTries << G4endl; 149 236 } 150 incl->processEventInclAbla( eventNumber);237 incl->processEventInclAbla(calincl, eventNumber); 151 238 152 239 if(verboseLevel > 1) { … … 159 246 * Diagnostic output 160 247 */ 161 G4cout <<"G4InclAblaLightIonInterface: Bullet type: " << bulletType << G4endl; 162 G4cout <<"G4Incl4AblaCascadeInterface: Bullet energy: " << bulletE << " MeV" << G4endl; 163 164 G4cout <<"G4InclAblaLightIonInterface: Target A: " << targetA << G4endl; 165 G4cout <<"G4InclAblaLightIonInterface: Target Z: " << targetZ << G4endl; 248 G4cout <<"G4InclAblaLightIonInterface: Bullet type: " << calincl->bulletType() << G4endl; 249 G4cout <<"G4Incl4AblaCascadeInterface: Bullet energy: " << calincl->bulletE() << " MeV" << G4endl; 250 if(bulletType == -666) { 251 G4cout <<" Extended projectile: A = " << calincl->extendedProjectileA() 252 <<" Z = " << calincl->extendedProjectileZ() << G4endl; 253 } 254 255 G4cout <<"G4InclAblaLightIonInterface: Target A: " << calincl->targetA() << G4endl; 256 G4cout <<"G4InclAblaLightIonInterface: Target Z: " << calincl->targetZ() << G4endl; 166 257 167 258 if(verboseLevel > 3) { 168 diagdata <<"G4InclAblaLightIonInterface: Bullet type: " << bulletType<< G4endl;169 diagdata <<"G4InclAblaLightIonInterface: Bullet energy: " << bulletE<< " MeV" << G4endl;259 diagdata <<"G4InclAblaLightIonInterface: Bullet type: " << calincl->bulletType() << G4endl; 260 diagdata <<"G4InclAblaLightIonInterface: Bullet energy: " << calincl->bulletE() << " MeV" << G4endl; 170 261 171 diagdata <<"G4InclAblaLightIonInterface: Target A: " << targetA << G4endl; 172 diagdata <<"G4InclAblaLightIonInterface: Target Z: " << targetZ << G4endl; 173 } 174 175 for(particleI = 0; particleI < varntp->ntrack; particleI++) { 176 G4cout << n << " " << calincl->f[6] << " " << calincl->f[2] << " "; 177 G4cout << varntp->massini << " " << varntp->mzini << " "; 178 G4cout << varntp->exini << " " << varntp->mulncasc << " " << varntp->mulnevap << " " << varntp->mulntot << " "; 179 G4cout << varntp->bimpact << " " << varntp->jremn << " " << varntp->kfis << " " << varntp->estfis << " "; 180 G4cout << varntp->izfis << " " << varntp->iafis << " " << varntp->ntrack << " " << varntp->itypcasc[particleI] << " "; 181 G4cout << varntp->avv[particleI] << " " << varntp->zvv[particleI] << " " << varntp->enerj[particleI] << " "; 182 G4cout << varntp->plab[particleI] << " " << varntp->tetlab[particleI] << " " << varntp->philab[particleI] << G4endl; 183 // For diagnostic output 184 if(verboseLevel > 3) { 185 diagdata << n << " " << calincl->f[6] << " " << calincl->f[2] << " "; 186 diagdata << varntp->massini << " " << varntp->mzini << " "; 187 diagdata << varntp->exini << " " << varntp->mulncasc << " " << varntp->mulnevap << " " << varntp->mulntot << " "; 188 diagdata << varntp->bimpact << " " << varntp->jremn << " " << varntp->kfis << " " << varntp->estfis << " "; 189 diagdata << varntp->izfis << " " << varntp->iafis << " " << varntp->ntrack << " "; 190 diagdata << varntp->itypcasc[particleI] << " "; 191 diagdata << varntp->avv[particleI] << " " << varntp->zvv[particleI] << " " << varntp->enerj[particleI] << " "; 192 diagdata << varntp->plab[particleI] << " " << varntp->tetlab[particleI] << " " << varntp->philab[particleI] << G4endl; 193 } 262 diagdata <<"G4InclAblaLightIonInterface: Target A: " << calincl->targetA() << G4endl; 263 diagdata <<"G4InclAblaLightIonInterface: Target Z: " << calincl->targetZ() << G4endl; 194 264 } 195 265 } … … 231 301 cascadeParticle->SetDefinition(aParticleDefinition); 232 302 cascadeParticle->Set4Momentum(aTrack.Get4Momentum()); 233 theResult.AddSecondary(cascadeParticle);303 result.push_back(cascadeParticle); 234 304 } 235 305 } … … 239 309 theResult.SetStatusChange(stopAndKill); 240 310 241 for(particleI = 0; particleI < varntp->ntrack; particleI++) { // Loop through the INCL4+ABLA output.311 for(particleI = 0; particleI <= varntp->ntrack; particleI++) { // Loop through the INCL4+ABLA output. 242 312 // Get energy/momentum and construct momentum vector in INCL4 coordinates. 313 // if(varntp->itypcasc[particleI] == -1) continue; // Avoid nucleons that are part of the spectator 314 if(varntp->avv[particleI] == 0 && varntp->zvv[particleI] == 0) continue; 243 315 momx = varntp->plab[particleI]*std::sin(varntp->tetlab[particleI]*CLHEP::pi/180.0)*std::cos(varntp->philab[particleI]*CLHEP::pi/180.0)*MeV; 244 316 momy = varntp->plab[particleI]*std::sin(varntp->tetlab[particleI]*CLHEP::pi/180.0)*std::sin(varntp->philab[particleI]*CLHEP::pi/180.0)*MeV; … … 262 334 new G4DynamicParticle(G4Proton::ProtonDefinition(), momDirection, eKin); 263 335 particleIdentified++; 336 baryonNumberBalanceInINCL -= 1; 337 chargeNumberBalanceInINCL -= 1; 264 338 } 265 339 … … 268 342 new G4DynamicParticle(G4Neutron::NeutronDefinition(), momDirection, eKin); 269 343 particleIdentified++; 344 baryonNumberBalanceInINCL -= 1; 270 345 } 271 346 … … 274 349 new G4DynamicParticle(G4PionPlus::PionPlusDefinition(), momDirection, eKin); 275 350 particleIdentified++; 351 chargeNumberBalanceInINCL -= 1; 276 352 } 277 353 … … 280 356 new G4DynamicParticle(G4PionZero::PionZeroDefinition(), momDirection, eKin); 281 357 particleIdentified++; 358 chargeNumberBalanceInINCL -= 0; 282 359 } 283 360 … … 286 363 new G4DynamicParticle(G4PionMinus::PionMinusDefinition(), momDirection, eKin); 287 364 particleIdentified++; 365 chargeNumberBalanceInINCL -= -1; 288 366 } 289 367 290 368 if((varntp->avv[particleI] > 1) && (varntp->zvv[particleI] >= 1)) { // Nucleus fragment 291 G4ParticleDefinition * aIonDef = 0; 292 G4ParticleTable *theTableOfParticles = G4ParticleTable::GetParticleTable(); 369 G4ParticleDefinition * aIonDef = 0; 293 370 294 371 G4int A = G4int(varntp->avv[particleI]); … … 313 390 new G4DynamicParticle(aIonDef, momDirection, eKin); 314 391 particleIdentified++; 392 baryonNumberBalanceInINCL -= A; 393 chargeNumberBalanceInINCL -= Z; 315 394 } 316 395 } 317 396 318 397 if(particleIdentified == 1) { // Particle identified properly. 319 320 theResult.AddSecondary(cascadeParticle); // Put data into G4HadFinalState.398 cascadeParticle->Set4Momentum(cascadeParticle->Get4Momentum()*=toLabFrame); 399 result.push_back(cascadeParticle); 321 400 } 322 401 else { // Particle identification failed. … … 332 411 } 333 412 413 // Spectator nucleus Fermi break-up 414 if(useFermiBreakup && useProjectileSpectator && varntp->masp > 1) { 415 baryonNumberBalanceInINCL -= G4int(varntp->masp); 416 G4double nuclearMass = G4NucleiProperties::GetNuclearMass(G4int(varntp->masp), G4int(varntp->mzsp)) + varntp->exsp * MeV; 417 // Use momentum scaling to compensate for different masses in G4 and INCL: 418 G4double momentumScaling = G4InclUtils::calculate4MomentumScaling(G4int(varntp->masp), 419 G4int(varntp->mzsp), 420 varntp->exsp, 421 varntp->spectatorT, 422 varntp->spectatorP1, 423 varntp->spectatorP2, 424 varntp->spectatorP3); 425 G4LorentzVector p4(momentumScaling * varntp->spectatorP1 * MeV, momentumScaling * varntp->spectatorP2 * MeV, 426 momentumScaling * varntp->spectatorP3 * MeV, 427 varntp->spectatorT * MeV + nuclearMass); 428 // Four-momentum, baryon number and charge balance: 429 G4LorentzVector fourMomentumBalance = p4; 430 G4int baryonNumberBalance = G4int(varntp->masp); 431 chargeNumberBalanceInINCL -= G4int(varntp->mzsp); 432 G4int chargeBalance = G4int(varntp->mzsp); 433 434 G4LorentzRotation toFragmentZ; 435 // Assume that Fermi breakup uses Z as the direction of the projectile 436 toFragmentZ.rotateZ(-p4.theta()); 437 toFragmentZ.rotateY(-p4.phi()); 438 G4LorentzRotation toFragmentLab = toFragmentZ.inverse(); 439 // p4 *= toFragmentZ; 440 441 G4LorentzVector p4rest = p4; 442 // p4rest.boost(-p4.boostVector()); 443 if(verboseLevel > 0) { 444 G4cout <<"Spectator nucleus:" << G4endl; 445 G4cout <<"p4: " << G4endl; 446 G4cout <<" px: " << p4.px() <<" py: " << p4.py() <<" pz: " << p4.pz() << G4endl; 447 G4cout <<" E = " << p4.e() << G4endl; 448 G4cout <<"p4rest: " << G4endl; 449 G4cout <<" px: " << p4rest.px() <<" py: " << p4rest.py() <<" pz: " << p4rest.pz() << G4endl; 450 G4cout <<" E = " << p4rest.e() << G4endl; 451 } 452 G4Fragment theSpectatorNucleus(G4int(varntp->masp), G4int(varntp->mzsp), p4rest); 453 theSpectatorFermiBreakupResult = fermiBreakUp->BreakItUp(theSpectatorNucleus); 454 if(theSpectatorFermiBreakupResult != 0) { 455 G4FragmentVector::iterator fragment; 456 for(fragment = theSpectatorFermiBreakupResult->begin(); fragment != theSpectatorFermiBreakupResult->end(); fragment++) { 457 G4ParticleDefinition *theFragmentDefinition = 0; 458 if((*fragment)->GetA_asInt() == 1 && (*fragment)->GetZ_asInt() == 0) { // Neutron 459 theFragmentDefinition = G4Neutron::NeutronDefinition(); 460 } else if ((*fragment)->GetA_asInt() == 1 && (*fragment)->GetZ_asInt() == 1) { 461 theFragmentDefinition = G4Proton::ProtonDefinition(); 462 } else { 463 theFragmentDefinition = theTableOfParticles->GetIon((*fragment)->GetZ_asInt(), (*fragment)->GetA_asInt(), (*fragment)->GetExcitationEnergy()); 464 } 465 if(theFragmentDefinition != 0) { 466 G4DynamicParticle *theFragment = new G4DynamicParticle(theFragmentDefinition, (*fragment)->GetMomentum()); 467 G4LorentzVector labMomentum = theFragment->Get4Momentum(); 468 // labMomentum.boost(p4.boostVector()); 469 // labMomentum *= toFragmentLab; 470 // labMomentum *= toLabFrame; 471 theFragment->Set4Momentum(labMomentum); 472 fourMomentumBalance -= theFragment->Get4Momentum(); 473 baryonNumberBalance -= theFragmentDefinition->GetAtomicMass(); 474 chargeBalance -= theFragmentDefinition->GetAtomicNumber(); 475 if(verboseLevel > 0) { 476 G4cout <<"Resulting fragment: " << G4endl; 477 G4cout <<" kinetic energy = " << theFragment->GetKineticEnergy() / MeV << " MeV" << G4endl; 478 G4cout <<" momentum = " << theFragment->GetMomentum().mag() / MeV << " MeV" << G4endl; 479 } 480 result.push_back(theFragment); 481 } else { 482 G4cout <<"G4InclAblaCascadeInterface: Error. Fragment produced by Fermi break-up does not exist." 483 << G4endl; 484 G4cout <<"Resulting fragment: " << G4endl; 485 G4cout <<" Z = " << (*fragment)->GetZ_asInt() << G4endl; 486 G4cout <<" A = " << (*fragment)->GetA_asInt() << G4endl; 487 G4cout <<" Excitation : " << (*fragment)->GetExcitationEnergy() / MeV << " MeV" << G4endl; 488 G4cout <<" momentum = " << (*fragment)->GetMomentum().mag() / MeV << " MeV" << G4endl; 489 } 490 } 491 if(std::abs(fourMomentumBalance.mag() / MeV) > 0.1 * MeV) { 492 G4cout <<"Four-momentum balance after spectator nucleus Fermi break-up:" << G4endl; 493 G4cout <<"Magnitude: " << fourMomentumBalance.mag() / MeV << " MeV" << G4endl; 494 G4cout <<"Vector components (px, py, pz, E) = (" 495 << fourMomentumBalance.px() << ", " 496 << fourMomentumBalance.py() << ", " 497 << fourMomentumBalance.pz() << ", " 498 << fourMomentumBalance.e() << ")" << G4endl; 499 } 500 if(baryonNumberBalance != 0) { 501 G4cout <<"Event " << eventNumber << ": Baryon number balance after spectator nucleus Fermi break-up: " << baryonNumberBalance << G4endl; 502 } 503 if(chargeBalance != 0) { 504 G4cout <<"Event " << eventNumber <<": Charge balance after spectator nucleus Fermi break-up: " << chargeBalance << G4endl; 505 } 506 } 507 } 508 509 // Finally do Fermi break-up if needed 510 if(varntp->needsFermiBreakup && varntp->massini > 0) { 511 baryonNumberBalanceInINCL -= G4int(varntp->massini); 512 chargeNumberBalanceInINCL -= G4int(varntp->mzini); 513 // Call Fermi Break-up 514 G4double nuclearMass = G4NucleiProperties::GetNuclearMass(G4int(varntp->massini), G4int(varntp->mzini)) + varntp->exini * MeV; 515 G4LorentzVector fragmentMomentum(varntp->pxrem * MeV, varntp->pyrem * MeV, varntp->pzrem * MeV, 516 varntp->erecrem * MeV + nuclearMass); 517 G4double momentumScaling = G4InclUtils::calculate4MomentumScaling(G4int(varntp->massini), G4int(varntp->mzini), 518 varntp->exini, 519 varntp->erecrem, 520 varntp->pxrem, 521 varntp->pyrem, 522 varntp->pzrem); 523 G4LorentzVector p4(momentumScaling * varntp->pxrem * MeV, momentumScaling * varntp->pyrem * MeV, 524 momentumScaling * varntp->pzrem * MeV, 525 varntp->erecrem + nuclearMass); 526 527 // For four-momentum, baryon number and charge conservation check: 528 G4LorentzVector fourMomentumBalance = p4; 529 G4int baryonNumberBalance = G4int(varntp->massini); 530 G4int chargeBalance = G4int(varntp->mzini); 531 532 G4LorentzRotation toFragmentZ; 533 toFragmentZ.rotateZ(-p4.theta()); 534 toFragmentZ.rotateY(-p4.phi()); 535 G4LorentzRotation toFragmentLab = toFragmentZ.inverse(); 536 // p4 *= toFragmentZ; 537 538 G4LorentzVector p4rest = p4; 539 // p4rest.boost(-p4.boostVector()); 540 if(verboseLevel > 0) { 541 G4cout <<"Cascade remnant nucleus:" << G4endl; 542 G4cout <<"p4: " << G4endl; 543 G4cout <<" px: " << p4.px() <<" py: " << p4.py() <<" pz: " << p4.pz() << G4endl; 544 G4cout <<" E = " << p4.e() << G4endl; 545 546 G4cout <<"p4rest: " << G4endl; 547 G4cout <<" px: " << p4rest.px() <<" py: " << p4rest.py() <<" pz: " << p4rest.pz() << G4endl; 548 G4cout <<" E = " << p4rest.e() << G4endl; 549 } 550 551 G4Fragment theCascadeRemnant(G4int(varntp->massini), G4int(varntp->mzini), p4rest); 552 theFermiBreakupResult = fermiBreakUp->BreakItUp(theCascadeRemnant); 553 if(theFermiBreakupResult != 0) { 554 G4FragmentVector::iterator fragment; 555 for(fragment = theFermiBreakupResult->begin(); fragment != theFermiBreakupResult->end(); fragment++) { 556 G4ParticleDefinition *theFragmentDefinition = 0; 557 if((*fragment)->GetA_asInt() == 1 && (*fragment)->GetZ_asInt() == 0) { // Neutron 558 theFragmentDefinition = G4Neutron::NeutronDefinition(); 559 } else if ((*fragment)->GetA_asInt() == 1 && (*fragment)->GetZ_asInt() == 1) { 560 theFragmentDefinition = G4Proton::ProtonDefinition(); 561 } else { 562 theFragmentDefinition = theTableOfParticles->GetIon((*fragment)->GetZ_asInt(), (*fragment)->GetA_asInt(), (*fragment)->GetExcitationEnergy()); 563 } 564 565 if(theFragmentDefinition != 0) { 566 G4DynamicParticle *theFragment = new G4DynamicParticle(theFragmentDefinition, (*fragment)->GetMomentum()); 567 G4LorentzVector labMomentum = theFragment->Get4Momentum(); 568 // labMomentum.boost(p4.boostVector()); 569 // labMomentum *= toFragmentLab; 570 // labMomentum *= toLabFrame; 571 theFragment->Set4Momentum(labMomentum); 572 fourMomentumBalance -= theFragment->Get4Momentum(); 573 baryonNumberBalance -= theFragmentDefinition->GetAtomicMass(); 574 chargeBalance -= theFragmentDefinition->GetAtomicNumber(); 575 if(verboseLevel > 0) { 576 G4cout <<"Resulting fragment: " << G4endl; 577 G4cout <<" kinetic energy = " << theFragment->GetKineticEnergy() / MeV << " MeV" << G4endl; 578 G4cout <<" momentum = " << theFragment->GetMomentum().mag() / MeV << " MeV" << G4endl; 579 } 580 result.push_back(theFragment); 581 } else { 582 G4cout <<"G4InclAblaCascadeInterface: Error. Fragment produced by Fermi break-up does not exist." << G4endl; 583 G4cout <<"Resulting fragment: " << G4endl; 584 G4cout <<" Z = " << (*fragment)->GetZ_asInt() << G4endl; 585 G4cout <<" A = " << (*fragment)->GetA_asInt() << G4endl; 586 G4cout <<" Excitation : " << (*fragment)->GetExcitationEnergy() / MeV << " MeV" << G4endl; 587 G4cout <<" momentum = " << (*fragment)->GetMomentum().mag() / MeV << " MeV" << G4endl; 588 } 589 } 590 if(std::abs(fourMomentumBalance.mag() / MeV) > 0.1 * MeV) { 591 G4cout <<"Four-momentum balance after remnant nucleus Fermi break-up:" << G4endl; 592 G4cout <<"Magnitude: " << fourMomentumBalance.mag() / MeV << " MeV" << G4endl; 593 G4cout <<"Vector components (px, py, pz, E) = (" 594 << fourMomentumBalance.px() << ", " 595 << fourMomentumBalance.py() << ", " 596 << fourMomentumBalance.pz() << ", " 597 << fourMomentumBalance.e() << ")" << G4endl; 598 } 599 if(baryonNumberBalance != 0) { 600 G4cout <<"Baryon number balance after remnant nucleus Fermi break-up: " << baryonNumberBalance << G4endl; 601 } 602 if(chargeBalance != 0) { 603 G4cout <<"Charge balance after remnant nucleus Fermi break-up: " << chargeBalance << G4endl; 604 } 605 } 606 } 607 334 608 varntp->ntrack = 0; // Clean up the number of generated particles in the event. 609 610 if(baryonNumberBalanceInINCL != 0 && verboseLevel > 1) { 611 G4cout <<"Event " << eventNumber <<": G4InclAblaLightIonInterface: Baryon number conservation problem in INCL detected!" << G4endl; 612 G4cout <<"Baryon number balance: " << baryonNumberBalanceInINCL << G4endl; 613 if(baryonNumberBalanceInINCL < 0) { 614 G4cout <<"Event " << eventNumber <<": Too many outcoming baryons!" << G4endl; 615 } else if(baryonNumberBalanceInINCL > 0) { 616 G4cout <<"Event " << eventNumber <<": Too few outcoming baryons!" << G4endl; 617 } 618 } 619 620 if(chargeNumberBalanceInINCL != 0 && verboseLevel > 1) { 621 G4cout <<"Event " << eventNumber <<": G4InclAblaLightIonInterface: Charge number conservation problem in INCL detected!" << G4endl; 622 G4cout <<"Event " << eventNumber <<": Charge number balance: " << chargeNumberBalanceInINCL << G4endl; 623 } 335 624 } 336 625 /** … … 344 633 cascadeParticle = new G4DynamicParticle(theTableOfParticles->FindParticle(aTrack.GetDefinition()), aTrack.Get4Momentum()); 345 634 346 theResult.AddSecondary(cascadeParticle);635 result.push_back(cascadeParticle); 347 636 348 637 if(verboseLevel > 1) { … … 364 653 } 365 654 366 if(( targetA == 1) && (targetZ== 1)) { // Unsupported target655 if((calincl->targetA() == 1) && (calincl->targetZ() == 1)) { // Unsupported target 367 656 if(verboseLevel > 1) { 368 657 G4cout <<"Unsupported target: " << G4endl; 369 G4cout <<"Target A: " << targetA<< G4endl;370 G4cout <<"TargetZ: " << targetZ<< G4endl;658 G4cout <<"Target A: " << calincl->targetA() << G4endl; 659 G4cout <<"TargetZ: " << calincl->targetZ() << G4endl; 371 660 } 372 661 if(verboseLevel > 3) { 373 662 diagdata <<"Unsupported target: " << G4endl; 374 diagdata <<"Target A: " << targetA<< G4endl;375 diagdata <<"TargetZ: " << targetZ<< G4endl;376 } 377 } 378 379 if( bulletE< 100) { // INCL does not support E < 100 MeV.663 diagdata <<"Target A: " << calincl->targetA() << G4endl; 664 diagdata <<"TargetZ: " << calincl->targetZ() << G4endl; 665 } 666 } 667 668 if(calincl->bulletE() < 100) { // INCL does not support E < 100 MeV. 380 669 if(verboseLevel > 1) { 381 G4cout <<"Unsupported bullet energy: " << bulletE<< " MeV. (Lower limit is 100 MeV)." << G4endl;670 G4cout <<"Unsupported bullet energy: " << calincl->bulletE() << " MeV. (Lower limit is 100 MeV)." << G4endl; 382 671 G4cout <<"WARNING: Returning the original bullet with original energy back to Geant4." << G4endl; 383 672 } 384 673 if(verboseLevel > 3) { 385 diagdata <<"Unsupported bullet energy: " << bulletE<< " MeV. (Lower limit is 100 MeV)." << G4endl;674 diagdata <<"Unsupported bullet energy: " << calincl->bulletE() << " MeV. (Lower limit is 100 MeV)." << G4endl; 386 675 } 387 676 } … … 392 681 } 393 682 683 // Finally copy the accumulated secondaries into the result collection: 684 G4ThreeVector boostVector = aTrack.Get4Momentum().boostVector(); 685 G4LorentzRotation boostBack = toBreit.inverse(); 686 687 for(std::vector<G4DynamicParticle*>::iterator i = result.begin(); i != result.end(); ++i) { 688 // If the calculation was performed in inverse kinematics we have to 689 // convert the result back... 690 if(calincl->isInverseKinematics()) { 691 G4LorentzVector mom = (*i)->Get4Momentum(); 692 mom.setPz(-1.0 * mom.pz()); // Reverse the z-component of the momentum vector 693 mom *= boostBack; 694 (*i)->Set4Momentum(mom); 695 } 696 theResult.AddSecondary((*i)); 697 } 698 699 delete calincl; 700 calincl = 0; 394 701 return &theResult; 395 702 } -
trunk/source/processes/hadronic/models/incl/src/G4InclCascadeInterface.cc
r819 r1340 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4InclCascadeInterface.cc,v 1.1 0 2007/12/10 16:32:02 gunterExp $26 // $Id: G4InclCascadeInterface.cc,v 1.12 2010/10/26 02:47:59 kaitanie Exp $ 27 27 // Translation of INCL4.2/ABLA V3 28 28 // Pekka Kaitaniemi, HIP (translation) … … 45 45 46 46 varntp = new G4VarNtp(); 47 calincl = new G4Calincl();47 inclInput = 0; 48 48 ws = new G4Ws(); 49 49 mat = new G4Mat(); 50 incl = new G4Incl(hazard, calincl, ws, mat, varntp);50 incl = new G4Incl(hazard, inclInput, ws, mat, varntp); 51 51 52 52 verboseLevel = 0; … … 57 57 delete hazard; 58 58 delete varntp; 59 delete calincl;59 delete inclInput; 60 60 delete ws; 61 61 delete mat; … … 67 67 G4int maxTries = 200; 68 68 69 G4int particleI, n = 0; 70 71 G4int bulletType = 0; 69 G4int particleI; 72 70 73 71 // Print diagnostic messages: 0 = silent, 1 and 2 = verbose … … 85 83 } 86 84 87 // INCL4 needs the energy in units MeV88 G4double bulletE = aTrack.GetKineticEnergy() * MeV;89 90 85 #ifdef DEBUGINCL 91 86 G4cout <<"Bullet energy = " << bulletE / MeV << G4endl; 92 87 #endif 93 94 G4double targetA = theNucleus.GetN();95 G4double targetZ = theNucleus.GetZ();96 88 97 89 G4double eKin; … … 110 102 G4LorentzRotation toLabFrame = toZ.inverse(); 111 103 104 inclInput = new G4InclInput(aTrack, theNucleus, false); 105 112 106 theResult.Clear(); // Make sure the output data structure is clean. 113 114 // Map Geant4 particle types to corresponding INCL4 types.115 enum bulletParticleType {nucleus = 0, proton = 1, neutron = 2, pionPlus = 3, pionZero = 4,116 pionMinus = 5, deuteron = 6, triton = 7, he3 = 8, he4 = 9};117 118 // Coding particles for use with INCL4 and ABLA119 if (aTrack.GetDefinition() == G4Proton::Proton() ) bulletType = proton;120 if (aTrack.GetDefinition() == G4Neutron::Neutron() ) bulletType = neutron;121 if (aTrack.GetDefinition() == G4PionPlus::PionPlus() ) bulletType = pionPlus;122 if (aTrack.GetDefinition() == G4PionMinus::PionMinus() ) bulletType = pionMinus;123 if (aTrack.GetDefinition() == G4PionZero::PionZero() ) bulletType = pionZero;124 107 125 108 #ifdef DEBUGINCL … … 137 120 #endif 138 121 139 for(int i = 0; i < 15; i++) {140 calincl->f[i] = 0.0; // Initialize INCL input data141 }142 143 122 // Check wheter the input is acceptable. 144 if((bulletType != 0) && ((targetA != 1) && (targetZ != 1))) { 145 calincl->f[0] = targetA; // Target mass number 146 calincl->f[1] = targetZ; // Charge number 147 calincl->f[6] = bulletType; // Type 148 calincl->f[2] = bulletE; // Energy [MeV] 149 calincl->f[5] = 1.0; // Time scaling 150 calincl->f[4] = 45.0; // Nuclear potential 151 123 if((inclInput->bulletType() != 0) && ((inclInput->targetA() != 1) && (inclInput->targetZ() != 1))) { 152 124 ws->nosurf = -2; // Nucleus surface, -2 = Woods-Saxon 153 125 ws->xfoisa = 8; … … 158 130 159 131 mat->nbmat = 1; 160 mat->amat[0] = int( calincl->f[0]);161 mat->zmat[0] = int( calincl->f[1]);132 mat->amat[0] = int(inclInput->targetA()); 133 mat->zmat[0] = int(inclInput->targetZ()); 162 134 163 135 incl->initIncl(true); … … 168 140 G4cout <<"G4InclCascadeInterface: Try number = " << nTries << G4endl; 169 141 } 170 incl->processEventIncl( );142 incl->processEventIncl(inclInput); 171 143 172 144 if(verboseLevel > 1) { … … 179 151 * Diagnostic output 180 152 */ 181 G4cout <<"G4InclCascadeInterface: Bullet type: " << bulletType<< G4endl;182 G4cout <<"G4Incl4AblaCascadeInterface: Bullet energy: " << bulletE<< " MeV" << G4endl;183 184 G4cout <<"G4InclCascadeInterface: Target A: " << targetA<< G4endl;185 G4cout <<"G4InclCascadeInterface: Target Z: " << targetZ<< G4endl;153 G4cout <<"G4InclCascadeInterface: Bullet type: " << inclInput->bulletType() << G4endl; 154 G4cout <<"G4Incl4AblaCascadeInterface: Bullet energy: " << inclInput->bulletE() << " MeV" << G4endl; 155 156 G4cout <<"G4InclCascadeInterface: Target A: " << inclInput->targetA() << G4endl; 157 G4cout <<"G4InclCascadeInterface: Target Z: " << inclInput->targetZ() << G4endl; 186 158 187 159 if(verboseLevel > 3) { 188 diagdata <<"G4InclCascadeInterface: Bullet type: " << bulletType<< G4endl;189 diagdata <<"G4InclCascadeInterface: Bullet energy: " << bulletE<< " MeV" << G4endl;160 diagdata <<"G4InclCascadeInterface: Bullet type: " << inclInput->bulletType() << G4endl; 161 diagdata <<"G4InclCascadeInterface: Bullet energy: " << inclInput->bulletE() << " MeV" << G4endl; 190 162 191 diagdata <<"G4InclCascadeInterface: Target A: " << targetA<< G4endl;192 diagdata <<"G4InclCascadeInterface: Target Z: " << targetZ<< G4endl;193 } 194 195 for(particleI = 0; particleI < varntp->ntrack; particleI++) {196 G4cout << n << " " << calincl->f[6] << " " << calincl->f[2] << " ";197 G4cout << varntp->massini << " " << varntp->mzini << " ";198 G4cout << varntp->exini << " " << varntp->mulncasc << " " << varntp->mulnevap << " " << varntp->mulntot << " ";199 G4cout << varntp->bimpact << " " << varntp->jremn << " " << varntp->kfis << " " << varntp->estfis << " ";200 G4cout << varntp->izfis << " " << varntp->iafis << " " << varntp->ntrack << " " << varntp->itypcasc[particleI] << " ";201 G4cout << varntp->avv[particleI] << " " << varntp->zvv[particleI] << " " << varntp->enerj[particleI] << " ";202 G4cout << varntp->plab[particleI] << " " << varntp->tetlab[particleI] << " " << varntp->philab[particleI] << G4endl;203 // For diagnostic output204 if(verboseLevel > 3) {205 diagdata << n << " " << calincl->f[6] << " " << calincl->f[2] << " ";206 diagdata << varntp->massini << " " << varntp->mzini << " ";207 diagdata << varntp->exini << " " << varntp->mulncasc << " " << varntp->mulnevap << " " << varntp->mulntot << " ";208 diagdata << varntp->bimpact << " " << varntp->jremn << " " << varntp->kfis << " " << varntp->estfis << " ";209 diagdata << varntp->izfis << " " << varntp->iafis << " " << varntp->ntrack << " ";210 diagdata << varntp->itypcasc[particleI] << " ";211 diagdata << varntp->avv[particleI] << " " << varntp->zvv[particleI] << " " << varntp->enerj[particleI] << " ";212 diagdata << varntp->plab[particleI] << " " << varntp->tetlab[particleI] << " " << varntp->philab[particleI] << G4endl;213 }214 }163 diagdata <<"G4InclCascadeInterface: Target A: " << inclInput->targetA() << G4endl; 164 diagdata <<"G4InclCascadeInterface: Target Z: " << inclInput->targetZ() << G4endl; 165 } 166 167 // for(particleI = 0; particleI < varntp->ntrack; particleI++) { 168 // G4cout << n << " " << inclInput->f[6] << " " << inclInput->f[2] << " "; 169 // G4cout << varntp->massini << " " << varntp->mzini << " "; 170 // G4cout << varntp->exini << " " << varntp->mulncasc << " " << varntp->mulnevap << " " << varntp->mulntot << " "; 171 // G4cout << varntp->bimpact << " " << varntp->jremn << " " << varntp->kfis << " " << varntp->estfis << " "; 172 // G4cout << varntp->izfis << " " << varntp->iafis << " " << varntp->ntrack << " " << varntp->itypcasc[particleI] << " "; 173 // G4cout << varntp->avv[particleI] << " " << varntp->zvv[particleI] << " " << varntp->enerj[particleI] << " "; 174 // G4cout << varntp->plab[particleI] << " " << varntp->tetlab[particleI] << " " << varntp->philab[particleI] << G4endl; 175 // // For diagnostic output 176 // if(verboseLevel > 3) { 177 // diagdata << n << " " << inclInput->f[6] << " " << inclInput->f[2] << " "; 178 // diagdata << varntp->massini << " " << varntp->mzini << " "; 179 // diagdata << varntp->exini << " " << varntp->mulncasc << " " << varntp->mulnevap << " " << varntp->mulntot << " "; 180 // diagdata << varntp->bimpact << " " << varntp->jremn << " " << varntp->kfis << " " << varntp->estfis << " "; 181 // diagdata << varntp->izfis << " " << varntp->iafis << " " << varntp->ntrack << " "; 182 // diagdata << varntp->itypcasc[particleI] << " "; 183 // diagdata << varntp->avv[particleI] << " " << varntp->zvv[particleI] << " " << varntp->enerj[particleI] << " "; 184 // diagdata << varntp->plab[particleI] << " " << varntp->tetlab[particleI] << " " << varntp->philab[particleI] << G4endl; 185 // } 186 // } 215 187 } 216 188 … … 224 196 225 197 theResult.SetStatusChange(stopAndKill); 226 227 if(bulletType == proton) { 228 aParticleDefinition = G4Proton::ProtonDefinition(); 229 } 230 if(bulletType == neutron) { 231 aParticleDefinition = G4Neutron::NeutronDefinition(); 232 } 233 if(bulletType == pionPlus) { 234 aParticleDefinition = G4PionPlus::PionPlusDefinition(); 235 } 236 if(bulletType == pionZero) { 237 aParticleDefinition = G4PionZero::PionZeroDefinition(); 238 } 239 if(bulletType == pionMinus) { 240 aParticleDefinition = G4PionMinus::PionMinusDefinition(); 241 } 198 aParticleDefinition = G4InclInput::getParticleDefinition(inclInput->bulletType()); 242 199 243 200 cascadeParticle = new G4DynamicParticle(); … … 434 391 } 435 392 436 if( bulletType== 0) {393 if(inclInput->bulletType() == 0) { 437 394 if(verboseLevel > 1) { 438 395 G4cout <<"G4InclCascadeInterface: Unknown bullet type" << G4endl; … … 445 402 } 446 403 447 if(( targetA == 1) && (targetZ== 1)) { // Unsupported target404 if((inclInput->targetA() == 1) && (inclInput->targetZ() == 1)) { // Unsupported target 448 405 if(verboseLevel > 1) { 449 406 G4cout <<"Unsupported target: " << G4endl; 450 G4cout <<"Target A: " << targetA<< G4endl;451 G4cout <<"TargetZ: " << targetZ<< G4endl;407 G4cout <<"Target A: " << inclInput->targetA() << G4endl; 408 G4cout <<"TargetZ: " << inclInput->targetZ() << G4endl; 452 409 } 453 410 if(verboseLevel > 3) { 454 411 diagdata <<"Unsupported target: " << G4endl; 455 diagdata <<"Target A: " << targetA<< G4endl;456 diagdata <<"TargetZ: " << targetZ<< G4endl;457 } 458 } 459 460 if( bulletE< 100) { // INCL does not support E < 100 MeV.461 if(verboseLevel > 1) { 462 G4cout <<"Unsupported bullet energy: " << bulletE<< " MeV. (Lower limit is 100 MeV)." << G4endl;412 diagdata <<"Target A: " << inclInput->targetA() << G4endl; 413 diagdata <<"TargetZ: " << inclInput->targetZ() << G4endl; 414 } 415 } 416 417 if(inclInput->bulletE() < 100) { // INCL does not support E < 100 MeV. 418 if(verboseLevel > 1) { 419 G4cout <<"Unsupported bullet energy: " << inclInput->bulletE() << " MeV. (Lower limit is 100 MeV)." << G4endl; 463 420 G4cout <<"WARNING: Returning the original bullet with original energy back to Geant4." << G4endl; 464 421 } 465 422 if(verboseLevel > 3) { 466 diagdata <<"Unsupported bullet energy: " << bulletE<< " MeV. (Lower limit is 100 MeV)." << G4endl;423 diagdata <<"Unsupported bullet energy: " << inclInput->bulletE() << " MeV. (Lower limit is 100 MeV)." << G4endl; 467 424 } 468 425 } -
trunk/source/processes/hadronic/models/incl/src/G4InclLightIonInterface.cc
r819 r1340 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4InclLightIonInterface.cc,v 1.1 0 2007/12/10 16:32:07 gunterExp $26 // $Id: G4InclLightIonInterface.cc,v 1.12 2010/10/26 02:47:59 kaitanie Exp $ 27 27 // Translation of INCL4.2/ABLA V3 28 28 // Pekka Kaitaniemi, HIP (translation) … … 43 43 44 44 varntp = new G4VarNtp(); 45 calincl = new G4Calincl();45 calincl = 0; 46 46 ws = new G4Ws(); 47 47 mat = new G4Mat(); … … 65 65 G4int maxTries = 200; 66 66 67 G4int particleI, n = 0;68 69 67 G4int bulletType = 0; 68 G4int particleI; 70 69 71 70 // Print diagnostic messages: 0 = silent, 1 and 2 = verbose … … 86 85 G4double bulletE = aTrack.GetKineticEnergy() / MeV; 87 86 88 G4 double targetA = theNucleus.GetN();89 G4 double targetZ = theNucleus.GetZ();87 G4int targetA = theNucleus.GetA_asInt(); 88 G4int targetZ = theNucleus.GetZ_asInt(); 90 89 91 90 G4double eKin; … … 116 115 if (aTrack.GetDefinition() == G4Alpha::Alpha() ) bulletType = he4; 117 116 118 for(int i = 0; i < 15; i++) { 119 calincl->f[i] = 0.0; // Initialize INCL input data 120 } 117 calincl = new G4InclInput(); 121 118 122 119 // Check wheter the input is acceptable. 123 120 if((bulletType != 0) && ((targetA != 1) && (targetZ != 1))) { 124 calincl->f[0] = targetA; // Target mass number125 calincl->f[1] = targetZ; // Charge number126 calincl->f[6] = bulletType; // Type127 calincl->f[2] = bulletE; // Energy [MeV]128 calincl->f[5] = 1.0; // Time scaling129 calincl->f[4] = 45.0; // Nuclear potential130 131 121 ws->nosurf = -2; // Nucleus surface, -2 = Woods-Saxon 132 122 ws->xfoisa = 8; … … 137 127 138 128 mat->nbmat = 1; 139 mat->amat[0] = int(calincl-> f[0]);140 mat->zmat[0] = int(calincl-> f[1]);129 mat->amat[0] = int(calincl->targetA()); 130 mat->zmat[0] = int(calincl->targetZ()); 141 131 142 132 incl->initIncl(true); … … 147 137 G4cout <<"G4InclLightIonInterface: Try number = " << nTries << G4endl; 148 138 } 149 incl->processEventIncl( );139 incl->processEventIncl(calincl); 150 140 151 141 if(verboseLevel > 1) { … … 170 160 diagdata <<"G4InclLightIonInterface: Target A: " << targetA << G4endl; 171 161 diagdata <<"G4InclLightIonInterface: Target Z: " << targetZ << G4endl; 172 }173 174 for(particleI = 0; particleI < varntp->ntrack; particleI++) {175 G4cout << n << " " << calincl->f[6] << " " << calincl->f[2] << " ";176 G4cout << varntp->massini << " " << varntp->mzini << " ";177 G4cout << varntp->exini << " " << varntp->mulncasc << " " << varntp->mulnevap << " " << varntp->mulntot << " ";178 G4cout << varntp->bimpact << " " << varntp->jremn << " " << varntp->kfis << " " << varntp->estfis << " ";179 G4cout << varntp->izfis << " " << varntp->iafis << " " << varntp->ntrack << " " << varntp->itypcasc[particleI] << " ";180 G4cout << varntp->avv[particleI] << " " << varntp->zvv[particleI] << " " << varntp->enerj[particleI] << " ";181 G4cout << varntp->plab[particleI] << " " << varntp->tetlab[particleI] << " " << varntp->philab[particleI] << G4endl;182 // For diagnostic output183 if(verboseLevel > 3) {184 diagdata << n << " " << calincl->f[6] << " " << calincl->f[2] << " ";185 diagdata << varntp->massini << " " << varntp->mzini << " ";186 diagdata << varntp->exini << " " << varntp->mulncasc << " " << varntp->mulnevap << " " << varntp->mulntot << " ";187 diagdata << varntp->bimpact << " " << varntp->jremn << " " << varntp->kfis << " " << varntp->estfis << " ";188 diagdata << varntp->izfis << " " << varntp->iafis << " " << varntp->ntrack << " ";189 diagdata << varntp->itypcasc[particleI] << " ";190 diagdata << varntp->avv[particleI] << " " << varntp->zvv[particleI] << " " << varntp->enerj[particleI] << " ";191 diagdata << varntp->plab[particleI] << " " << varntp->tetlab[particleI] << " " << varntp->philab[particleI] << G4endl;192 }193 162 } 194 163 }
Note: See TracChangeset
for help on using the changeset viewer.