Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/em_dissociation/src/G4EMDissociation.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /processes/hadronic/models/em_dissociation/src/G4EMDissociation.cc (Version 11.3.0) and /processes/hadronic/models/em_dissociation/src/G4EMDissociation.cc (Version 10.0.p2)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // *                                               20 // *                                                                  *
 21 // * Parts of this code which have been  devel     21 // * Parts of this code which have been  developed by QinetiQ Ltd     *
 22 // * under contract to the European Space Agen     22 // * under contract to the European Space Agency (ESA) are the        *
 23 // * intellectual property of ESA. Rights to u     23 // * intellectual property of ESA. Rights to use, copy, modify and    *
 24 // * redistribute this software for general pu     24 // * redistribute this software for general public use are granted    *
 25 // * in compliance with any licensing, distrib     25 // * in compliance with any licensing, distribution and development   *
 26 // * policy adopted by the Geant4 Collaboratio     26 // * policy adopted by the Geant4 Collaboration. This code has been   *
 27 // * written by QinetiQ Ltd for the European S     27 // * written by QinetiQ Ltd for the European Space Agency, under ESA  *
 28 // * contract 17191/03/NL/LvH (Aurora Programm     28 // * contract 17191/03/NL/LvH (Aurora Programme).                     *
 29 // *                                               29 // *                                                                  *
 30 // * By using,  copying,  modifying or  distri     30 // * By using,  copying,  modifying or  distributing the software (or *
 31 // * any work based  on the software)  you  ag     31 // * any work based  on the software)  you  agree  to acknowledge its *
 32 // * use  in  resulting  scientific  publicati     32 // * use  in  resulting  scientific  publications,  and indicate your *
 33 // * acceptance of all terms of the Geant4 Sof     33 // * acceptance of all terms of the Geant4 Software license.          *
 34 // *******************************************     34 // ********************************************************************
 35 //                                                 35 //
 36 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%     36 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 37 //                                                 37 //
 38 // MODULE:    G4EMDissociation.cc                  38 // MODULE:    G4EMDissociation.cc
 39 //                                                 39 //
 40 // Version:   B.1                                  40 // Version:   B.1
 41 // Date:    15/04/04                               41 // Date:    15/04/04
 42 // Author:    P R Truscott                         42 // Author:    P R Truscott
 43 // Organisation:  QinetiQ Ltd, UK                  43 // Organisation:  QinetiQ Ltd, UK
 44 // Customer:    ESA/ESTEC, NOORDWIJK               44 // Customer:    ESA/ESTEC, NOORDWIJK
 45 // Contract:    17191/03/NL/LvH                    45 // Contract:    17191/03/NL/LvH
 46 //                                                 46 //
 47 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%     47 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 48 //                                                 48 //
 49 // CHANGE HISTORY                                  49 // CHANGE HISTORY
 50 // --------------                                  50 // --------------
 51 //                                                 51 //
 52 // 17 October 2003, P R Truscott, QinetiQ Ltd,     52 // 17 October 2003, P R Truscott, QinetiQ Ltd, UK
 53 // Created.                                        53 // Created.
 54 //                                                 54 //
 55 // 15 March 2004, P R Truscott, QinetiQ Ltd, U     55 // 15 March 2004, P R Truscott, QinetiQ Ltd, UK
 56 // Beta release                                    56 // Beta release
 57 //                                                 57 //
 58 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%     58 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 59 //////////////////////////////////////////////     59 ////////////////////////////////////////////////////////////////////////////////
 60 //                                                 60 //
 61 #include "G4EMDissociation.hh"                     61 #include "G4EMDissociation.hh"
 62 #include "G4PhysicalConstants.hh"                  62 #include "G4PhysicalConstants.hh"
 63 #include "G4SystemOfUnits.hh"                      63 #include "G4SystemOfUnits.hh"
                                                   >>  64 #include "G4Evaporation.hh"
                                                   >>  65 #include "G4FermiBreakUp.hh"
                                                   >>  66 #include "G4StatMF.hh"
 64 #include "G4ParticleDefinition.hh"                 67 #include "G4ParticleDefinition.hh"
 65 #include "G4LorentzVector.hh"                      68 #include "G4LorentzVector.hh"
 66 #include "G4PhysicsFreeVector.hh"                  69 #include "G4PhysicsFreeVector.hh"
 67 #include "G4EMDissociationCrossSection.hh"         70 #include "G4EMDissociationCrossSection.hh"
 68 #include "G4Proton.hh"                             71 #include "G4Proton.hh"
 69 #include "G4Neutron.hh"                            72 #include "G4Neutron.hh"
                                                   >>  73 #include "G4ParticleTable.hh"
 70 #include "G4IonTable.hh"                           74 #include "G4IonTable.hh"
                                                   >>  75 #include "G4GeneralPhaseSpaceDecay.hh"
 71 #include "G4DecayProducts.hh"                      76 #include "G4DecayProducts.hh"
 72 #include "G4DynamicParticle.hh"                    77 #include "G4DynamicParticle.hh"
 73 #include "G4Fragment.hh"                           78 #include "G4Fragment.hh"
 74 #include "G4ReactionProductVector.hh"              79 #include "G4ReactionProductVector.hh"
 75 #include "Randomize.hh"                            80 #include "Randomize.hh"
 76 #include "globals.hh"                              81 #include "globals.hh"
 77 #include "G4PhysicsModelCatalog.hh"            << 
 78                                                    82 
 79 G4EMDissociation::G4EMDissociation() :         <<  83 G4EMDissociation::G4EMDissociation():G4HadronicInteraction("EMDissociation") {
 80   G4HadronicInteraction("EMDissociation"),     <<  84 
 81   secID_projectileDissociation(-1), secID_targ << 
 82 {                                              << 
 83   // Send message to stdout to advise that the     85   // Send message to stdout to advise that the G4EMDissociation model is being
 84   // used.                                         86   // used.
 85   PrintWelcomeMessage();                           87   PrintWelcomeMessage();
 86                                                    88 
 87   // No de-excitation handler has been supplie     89   // No de-excitation handler has been supplied - define the default handler.
 88   theExcitationHandler            = new G4Exci     90   theExcitationHandler            = new G4ExcitationHandler;
                                                   >>  91   G4Evaporation* theEvaporation   = new G4Evaporation;
                                                   >>  92   G4FermiBreakUp* theFermiBreakUp = new G4FermiBreakUp;
                                                   >>  93   G4StatMF* theMF                 = new G4StatMF;
                                                   >>  94   theExcitationHandler->SetEvaporation(theEvaporation);
                                                   >>  95   theExcitationHandler->SetFermiModel(theFermiBreakUp);
                                                   >>  96   theExcitationHandler->SetMultiFragmentation(theMF);
                                                   >>  97   theExcitationHandler->SetMaxAandZForFermiBreakUp(12, 6);
 89   theExcitationHandler->SetMinEForMultiFrag(5.     98   theExcitationHandler->SetMinEForMultiFrag(5.0*MeV);
 90   handlerDefinedInternally = true;                 99   handlerDefinedInternally = true;
 91                                                   100 
 92   // This EM dissociation model needs access t    101   // This EM dissociation model needs access to the cross-sections held in
 93   // G4EMDissociationCrossSection.                102   // G4EMDissociationCrossSection.
 94   dissociationCrossSection = new G4EMDissociat    103   dissociationCrossSection = new G4EMDissociationCrossSection;
 95   thePhotonSpectrum = new G4EMDissociationSpec    104   thePhotonSpectrum = new G4EMDissociationSpectrum;
 96                                                   105 
 97   // Set the minimum and maximum range for the    106   // Set the minimum and maximum range for the model (despite nomanclature, this
 98   // is in energy per nucleon number).            107   // is in energy per nucleon number).    
 99   SetMinEnergy(100.0*MeV);                        108   SetMinEnergy(100.0*MeV);
100   SetMaxEnergy(500.0*GeV);                        109   SetMaxEnergy(500.0*GeV);
101                                                   110 
102   // Set the default verbose level to 0 - no o    111   // Set the default verbose level to 0 - no output.
103   verboseLevel = 0;                               112   verboseLevel = 0;
                                                   >> 113 }
104                                                   114 
105   // Creator model ID for the secondaries crea << 115 /*
106   secID_projectileDissociation = G4PhysicsMode << 116 G4EMDissociation::G4EMDissociation(const G4EMDissociation& emd)
107   secID_targetDissociation     = G4PhysicsMode << 117  : G4HadronicInteraction(emd)
                                                   >> 118 {
                                                   >> 119   if (emd.theExcitationHandler != 0) {
                                                   >> 120     theExcitationHandler = new G4ExcitationHandler;
                                                   >> 121     *theExcitationHandler = *emd.theExcitationHandler;
                                                   >> 122   }
                                                   >> 123 
                                                   >> 124   handlerDefinedInternally = emd.handlerDefinedInternally;
                                                   >> 125 
                                                   >> 126   if (emd.dissociationCrossSection != 0) {
                                                   >> 127     dissociationCrossSection = new G4EMDissociationCrossSection;
                                                   >> 128     *dissociationCrossSection = *emd.dissociationCrossSection;
                                                   >> 129   }
                                                   >> 130 
                                                   >> 131   if (emd.thePhotonSpectrum !- 0) {
                                                   >> 132     thePhotonSpectrum = new G4EMDissociationSpectrum;
                                                   >> 133     *thePhotonSpectrum = *emd.thePhotonSpectrum;
108 }                                                 134 }
                                                   >> 135 */
109                                                   136 
110 G4EMDissociation::G4EMDissociation (G4Excitati << 137 G4EMDissociation::G4EMDissociation (G4ExcitationHandler *aExcitationHandler)
111   G4HadronicInteraction("EMDissociation"),     << 
112   secID_projectileDissociation(-1), secID_targ << 
113 {                                                 138 {
114   // Send message to stdout to advise that the    139   // Send message to stdout to advise that the G4EMDissociation model is being
115   // used.                                        140   // used.
116   PrintWelcomeMessage();                          141   PrintWelcomeMessage();
117                                                   142   
118   theExcitationHandler     = aExcitationHandle    143   theExcitationHandler     = aExcitationHandler;
119   handlerDefinedInternally = false;               144   handlerDefinedInternally = false;
120                                                   145 
121   // This EM dissociation model needs access t    146   // This EM dissociation model needs access to the cross-sections held in
122   // G4EMDissociationCrossSection.                147   // G4EMDissociationCrossSection.
123   dissociationCrossSection = new G4EMDissociat    148   dissociationCrossSection = new G4EMDissociationCrossSection;
124   thePhotonSpectrum = new G4EMDissociationSpec    149   thePhotonSpectrum = new G4EMDissociationSpectrum;
125                                                   150 
126   // Set the minimum and maximum range for the    151   // Set the minimum and maximum range for the model (despite nomanclature, this
127   // is in energy per nucleon number)             152   // is in energy per nucleon number)    
128   SetMinEnergy(100.0*MeV);                        153   SetMinEnergy(100.0*MeV);
129   SetMaxEnergy(500.0*GeV);                        154   SetMaxEnergy(500.0*GeV);
130   verboseLevel = 0;                               155   verboseLevel = 0;
131                                                << 
132   // Creator model ID for the secondaries crea << 
133   secID_projectileDissociation = G4PhysicsMode << 
134   secID_targetDissociation     = G4PhysicsMode << 
135 }                                                 156 }
136                                                   157 
137                                                   158 
138 G4EMDissociation::~G4EMDissociation() {           159 G4EMDissociation::~G4EMDissociation() {
139   if (handlerDefinedInternally) delete theExci    160   if (handlerDefinedInternally) delete theExcitationHandler;
140   // delete dissociationCrossSection;             161   // delete dissociationCrossSection;
141   // Cross section deleted by G4CrossSectionRe    162   // Cross section deleted by G4CrossSectionRegistry; don't do it here
142   // Bug reported by Gong Ding in Bug Report #    163   // Bug reported by Gong Ding in Bug Report #1339
143   delete thePhotonSpectrum;                       164   delete thePhotonSpectrum;
144 }                                                 165 }
145                                                   166 
146                                                   167 
147 G4HadFinalState *G4EMDissociation::ApplyYourse    168 G4HadFinalState *G4EMDissociation::ApplyYourself
148   (const G4HadProjectile &theTrack, G4Nucleus     169   (const G4HadProjectile &theTrack, G4Nucleus &theTarget)
149 {                                                 170 {
150   // The secondaries will be returned in G4Had    171   // The secondaries will be returned in G4HadFinalState &theParticleChange -
151   // initialise this.                             172   // initialise this.
152                                                   173 
153   theParticleChange.Clear();                      174   theParticleChange.Clear();
154   theParticleChange.SetStatusChange(stopAndKil    175   theParticleChange.SetStatusChange(stopAndKill);
155                                                   176 
156   // Get relevant information about the projec    177   // Get relevant information about the projectile and target (A, Z) and
157   // energy/nuc, momentum, velocity, Lorentz f    178   // energy/nuc, momentum, velocity, Lorentz factor and rest-mass of the
158   // projectile.                                  179   // projectile.
159                                                   180 
160   const G4ParticleDefinition *definitionP = th    181   const G4ParticleDefinition *definitionP = theTrack.GetDefinition();
161   const G4double AP  = definitionP->GetBaryonN    182   const G4double AP  = definitionP->GetBaryonNumber();
162   const G4double ZP  = definitionP->GetPDGChar    183   const G4double ZP  = definitionP->GetPDGCharge();
163   G4LorentzVector pP = theTrack.Get4Momentum()    184   G4LorentzVector pP = theTrack.Get4Momentum();
164   G4double E         = theTrack.GetKineticEner    185   G4double E         = theTrack.GetKineticEnergy()/AP;
165   G4double MP        = theTrack.GetTotalEnergy    186   G4double MP        = theTrack.GetTotalEnergy() - E*AP;
166   G4double b         = pP.beta();                 187   G4double b         = pP.beta();
167   G4double AT        = theTarget.GetA_asInt();    188   G4double AT        = theTarget.GetA_asInt();
168   G4double ZT        = theTarget.GetZ_asInt();    189   G4double ZT        = theTarget.GetZ_asInt();
169   G4double MT        = G4NucleiProperties::Get    190   G4double MT        = G4NucleiProperties::GetNuclearMass(AT,ZT);
170                                                   191 
171   // Depending upon the verbosity level, outpu    192   // Depending upon the verbosity level, output the initial information on the
172   // projectile and target                        193   // projectile and target
173   if (verboseLevel >= 2) {                        194   if (verboseLevel >= 2) {
174     G4cout.precision(6);                          195     G4cout.precision(6);
175     G4cout <<"################################    196     G4cout <<"########################################"
176            <<"################################    197            <<"########################################"
177            <<G4endl;                              198            <<G4endl;
178     G4cout <<"IN G4EMDissociation" <<G4endl;      199     G4cout <<"IN G4EMDissociation" <<G4endl;
179     G4cout <<"Initial projectile A=" <<AP         200     G4cout <<"Initial projectile A=" <<AP 
180            <<", Z=" <<ZP                          201            <<", Z=" <<ZP
181            <<G4endl;                              202            <<G4endl; 
182     G4cout <<"Initial target     A=" <<AT         203     G4cout <<"Initial target     A=" <<AT
183            <<", Z=" <<ZT                          204            <<", Z=" <<ZT
184            <<G4endl;                              205            <<G4endl;
185     G4cout <<"Projectile momentum and Energy/n    206     G4cout <<"Projectile momentum and Energy/nuc = " <<pP <<" ," <<E <<G4endl;
186   }                                               207   }
187                                                   208 
188   // Initialise the variables which will be us    209   // Initialise the variables which will be used with the phase-space decay and
189   // to boost the secondaries from the interac    210   // to boost the secondaries from the interaction.
190                                                   211   
191   G4ParticleDefinition *typeNucleon  = NULL;      212   G4ParticleDefinition *typeNucleon  = NULL;
192   G4ParticleDefinition *typeDaughter = NULL;      213   G4ParticleDefinition *typeDaughter = NULL;
193   G4double Eg                        = 0.0;       214   G4double Eg                        = 0.0;
194   G4double mass                      = 0.0;       215   G4double mass                      = 0.0;
195   G4ThreeVector boost = G4ThreeVector(0.0, 0.0    216   G4ThreeVector boost = G4ThreeVector(0.0, 0.0, 0.0);
196                                                   217 
197   // Determine the cross-sections at the giant    218   // Determine the cross-sections at the giant dipole and giant quadrupole
198   // resonance energies for the projectile and    219   // resonance energies for the projectile and then target.  The information is
199   // initially provided in the G4PhysicsFreeVe    220   // initially provided in the G4PhysicsFreeVector individually for the E1
200   // and E2 fields. These are then summed.        221   // and E2 fields. These are then summed.
201                                                   222 
202   G4double bmin = thePhotonSpectrum->GetCloses    223   G4double bmin = thePhotonSpectrum->GetClosestApproach(AP, ZP, AT, ZT, b);
203   G4PhysicsFreeVector *crossSectionP = dissoci    224   G4PhysicsFreeVector *crossSectionP = dissociationCrossSection->
204     GetCrossSectionForProjectile(AP, ZP, AT, Z    225     GetCrossSectionForProjectile(AP, ZP, AT, ZT, b, bmin);
205   G4PhysicsFreeVector *crossSectionT = dissoci    226   G4PhysicsFreeVector *crossSectionT = dissociationCrossSection->
206     GetCrossSectionForTarget(AP, ZP, AT, ZT, b    227     GetCrossSectionForTarget(AP, ZP, AT, ZT, b, bmin);
207                                                   228 
208   G4double totCrossSectionP = (*crossSectionP)    229   G4double totCrossSectionP = (*crossSectionP)[0]+(*crossSectionP)[1];
209   G4double totCrossSectionT = (*crossSectionT)    230   G4double totCrossSectionT = (*crossSectionT)[0]+(*crossSectionT)[1];
210                                                   231 
211   // Now sample whether the interaction involv    232   // Now sample whether the interaction involved EM dissociation of the projectile
212   // or the target.                               233   // or the target.
213                                                << 234 
214   G4int secID = -1;  // Creator model ID for t << 
215   if (G4UniformRand() <                           235   if (G4UniformRand() <
216     totCrossSectionP / (totCrossSectionP + tot    236     totCrossSectionP / (totCrossSectionP + totCrossSectionT)) {
217                                                   237 
218     // It was the projectile which underwent E    238     // It was the projectile which underwent EM dissociation.  Define the Lorentz
219     // boost to be applied to the secondaries,    239     // boost to be applied to the secondaries, and sample whether a proton or a
220     // neutron was ejected.  Then determine th    240     // neutron was ejected.  Then determine the energy of the virtual gamma ray
221     // which passed from the target nucleus ..    241     // which passed from the target nucleus ... this will be used to define the
222     // excitation of the projectile.              242     // excitation of the projectile.
223                                                   243 
224     secID = secID_projectileDissociation;      << 
225     mass  = MP;                                   244     mass  = MP;
226     if (G4UniformRand() < dissociationCrossSec    245     if (G4UniformRand() < dissociationCrossSection->
227       GetWilsonProbabilityForProtonDissociatio    246       GetWilsonProbabilityForProtonDissociation (AP, ZP))
228     {                                             247     {
229       if (verboseLevel >= 2)                      248       if (verboseLevel >= 2)
230         G4cout <<"Projectile underwent EM diss    249         G4cout <<"Projectile underwent EM dissociation producing a proton"
231                <<G4endl;                          250                <<G4endl;
232       typeNucleon = G4Proton::ProtonDefinition    251       typeNucleon = G4Proton::ProtonDefinition();
233       typeDaughter = G4IonTable::GetIonTable() << 252       typeDaughter = G4ParticleTable::GetParticleTable()->
234       GetIon((G4int) ZP-1, (G4int) AP-1, 0.0);    253       GetIon((G4int) ZP-1, (G4int) AP-1, 0.0);
235     }                                             254     }
236     else                                          255     else
237     {                                             256     {
238       if (verboseLevel >= 2)                      257       if (verboseLevel >= 2)
239         G4cout <<"Projectile underwent EM diss    258         G4cout <<"Projectile underwent EM dissociation producing a neutron"
240                <<G4endl;                          259                <<G4endl;
241       typeNucleon = G4Neutron::NeutronDefiniti    260       typeNucleon = G4Neutron::NeutronDefinition();
242       typeDaughter = G4IonTable::GetIonTable() << 261       typeDaughter = G4ParticleTable::GetParticleTable()->
243       GetIon((G4int) ZP, (G4int) AP-1, 0.0);      262       GetIon((G4int) ZP, (G4int) AP-1, 0.0);
244     }                                             263     }
245     if (G4UniformRand() < (*crossSectionP)[0]/    264     if (G4UniformRand() < (*crossSectionP)[0]/totCrossSectionP)
246     {                                             265     {
247       Eg = crossSectionP->GetLowEdgeEnergy(0);    266       Eg = crossSectionP->GetLowEdgeEnergy(0);
248       if (verboseLevel >= 2)                      267       if (verboseLevel >= 2)
249         G4cout <<"Transition type was E1" <<G4    268         G4cout <<"Transition type was E1" <<G4endl;
250     }                                             269     }
251     else                                          270     else
252     {                                             271     {
253       Eg = crossSectionP->GetLowEdgeEnergy(1);    272       Eg = crossSectionP->GetLowEdgeEnergy(1);
254       if (verboseLevel >= 2)                      273       if (verboseLevel >= 2)
255         G4cout <<"Transition type was E2" <<G4    274         G4cout <<"Transition type was E2" <<G4endl;
256     }                                             275     }
257                                                   276 
258     // We need to define a Lorentz vector with    277     // We need to define a Lorentz vector with the original momentum, but total
259     // energy includes the projectile and virt    278     // energy includes the projectile and virtual gamma.  This is then used
260     // to calculate the boost required for the    279     // to calculate the boost required for the secondaries.
261                                                   280 
262     pP.setE( std::sqrt( pP.vect().mag2() + (ma << 281     pP.setE(pP.e()+Eg);
263     boost = pP.findBoostToCM();                   282     boost = pP.findBoostToCM();
264   }                                               283   }
265   else                                            284   else
266   {                                               285   {
267     // It was the target which underwent EM di    286     // It was the target which underwent EM dissociation.  Sample whether a
268     // proton or a neutron was ejected.  Then     287     // proton or a neutron was ejected.  Then determine the energy of the virtual 
269     // gamma ray which passed from the project    288     // gamma ray which passed from the projectile nucleus ... this will be used to
270     // define the excitation of the target.       289     // define the excitation of the target.
271                                                << 290 
272     secID = secID_targetDissociation;          << 
273     mass = MT;                                    291     mass = MT;
274     if (G4UniformRand() < dissociationCrossSec    292     if (G4UniformRand() < dissociationCrossSection->
275       GetWilsonProbabilityForProtonDissociatio    293       GetWilsonProbabilityForProtonDissociation (AT, ZT))
276     {                                             294     {
277       if (verboseLevel >= 2)                      295       if (verboseLevel >= 2)
278         G4cout <<"Target underwent EM dissocia    296         G4cout <<"Target underwent EM dissociation producing a proton"
279                <<G4endl;                          297                <<G4endl;
280       typeNucleon = G4Proton::ProtonDefinition    298       typeNucleon = G4Proton::ProtonDefinition();
281       typeDaughter = G4IonTable::GetIonTable() << 299       typeDaughter = G4ParticleTable::GetParticleTable()->
282       GetIon((G4int) ZT-1, (G4int) AT-1, 0.0);    300       GetIon((G4int) ZT-1, (G4int) AT-1, 0.0);
283     }                                             301     }
284     else                                          302     else
285     {                                             303     {
286       if (verboseLevel >= 2)                      304       if (verboseLevel >= 2)
287         G4cout <<"Target underwent EM dissocia    305         G4cout <<"Target underwent EM dissociation producing a neutron"
288                <<G4endl;                          306                <<G4endl;
289       typeNucleon = G4Neutron::NeutronDefiniti    307       typeNucleon = G4Neutron::NeutronDefinition();
290       typeDaughter = G4IonTable::GetIonTable() << 308       typeDaughter = G4ParticleTable::GetParticleTable()->
291       GetIon((G4int) ZT, (G4int) AT-1, 0.0);      309       GetIon((G4int) ZT, (G4int) AT-1, 0.0);
292     }                                             310     }
293     if (G4UniformRand() < (*crossSectionT)[0]/    311     if (G4UniformRand() < (*crossSectionT)[0]/totCrossSectionT)
294     {                                             312     {
295       Eg = crossSectionT->GetLowEdgeEnergy(0);    313       Eg = crossSectionT->GetLowEdgeEnergy(0);
296       if (verboseLevel >= 2)                      314       if (verboseLevel >= 2)
297         G4cout <<"Transition type was E1" <<G4    315         G4cout <<"Transition type was E1" <<G4endl;
298     }                                             316     }
299     else                                          317     else
300     {                                             318     {
301       Eg = crossSectionT->GetLowEdgeEnergy(1);    319       Eg = crossSectionT->GetLowEdgeEnergy(1);
302       if (verboseLevel >= 2)                      320       if (verboseLevel >= 2)
303         G4cout <<"Transition type was E2" <<G4    321         G4cout <<"Transition type was E2" <<G4endl;
304     }                                             322     }
305                                                   323 
306     // Add the projectile to theParticleChange    324     // Add the projectile to theParticleChange, less the energy of the
307     // not-so-virtual gamma-ray.  Not that at     325     // not-so-virtual gamma-ray.  Not that at the moment, no lateral momentum
308     // is transferred between the projectile a    326     // is transferred between the projectile and target nuclei.
309                                                   327 
310     G4ThreeVector v = pP.vect();                  328     G4ThreeVector v = pP.vect();
311     v.setMag(1.0);                                329     v.setMag(1.0);
312     G4DynamicParticle *changedP = new G4Dynami << 330     G4DynamicParticle *changedP = new G4DynamicParticle
313     theParticleChange.AddSecondary (changedP,  << 331       (const_cast<G4ParticleDefinition*>(definitionP), v, E*AP-Eg);
                                                   >> 332     theParticleChange.AddSecondary (changedP);
314     if (verboseLevel >= 2)                        333     if (verboseLevel >= 2)
315     {                                             334     {
316       G4cout <<"Projectile change:" <<G4endl;     335       G4cout <<"Projectile change:" <<G4endl;
317       changedP->DumpInfo();                       336       changedP->DumpInfo();
318     }                                             337     }
319   }                                               338   }
320                                                   339 
321   // Perform a two-body decay based on the res    340   // Perform a two-body decay based on the restmass energy of the parent and
322   // gamma-ray, and the masses of the daughter    341   // gamma-ray, and the masses of the daughters. In the frame of reference of
323   // the nucles, the angular distribution is s    342   // the nucles, the angular distribution is sampled isotropically, but the
324   // the nucleon and secondary nucleus are boo    343   // the nucleon and secondary nucleus are boosted if they've come from the
325   // projectile.                                  344   // projectile.
326                                                   345 
327   G4double e  = mass + Eg;                        346   G4double e  = mass + Eg;
328   G4double mass1 = typeNucleon->GetPDGMass();     347   G4double mass1 = typeNucleon->GetPDGMass();
329   G4double mass2 = typeDaughter->GetPDGMass();    348   G4double mass2 = typeDaughter->GetPDGMass();
330   G4double pp = (e+mass1+mass2)*(e+mass1-mass2    349   G4double pp = (e+mass1+mass2)*(e+mass1-mass2)*
331                 (e-mass1+mass2)*(e-mass1-mass2    350                 (e-mass1+mass2)*(e-mass1-mass2)/(4.0*e*e);
332   if (pp < 0.0) {                                 351   if (pp < 0.0) {
333     pp = 1.0*eV;                                  352     pp = 1.0*eV;
334 //    if (verboseLevel >`= 1)                     353 //    if (verboseLevel >`= 1)
335 //    {                                           354 //    {
336 //      G4cout <<"IN G4EMDissociation::ApplyYo    355 //      G4cout <<"IN G4EMDissociation::ApplyYoursef" <<G4endl;
337 //      G4cout <<"Error in mass of secondaries    356 //      G4cout <<"Error in mass of secondaries compared with primary:" <<G4endl;
338 //      G4cout <<"Rest mass of primary      =     357 //      G4cout <<"Rest mass of primary      = " <<mass <<" MeV" <<G4endl;
339 //      G4cout <<"Virtual gamma energy      =     358 //      G4cout <<"Virtual gamma energy      = " <<Eg   <<" MeV" <<G4endl;
340 //      G4cout <<"Rest mass of secondary #1 =     359 //      G4cout <<"Rest mass of secondary #1 = " <<mass1   <<" MeV" <<G4endl;
341 //      G4cout <<"Rest mass of secondary #2 =     360 //      G4cout <<"Rest mass of secondary #2 = " <<mass2   <<" MeV" <<G4endl;
342 //    }                                           361 //    }
343   }                                               362   }
344   else                                            363   else
345     pp = std::sqrt(pp);                           364     pp = std::sqrt(pp);
346   G4double costheta = 2.*G4UniformRand()-1.0;     365   G4double costheta = 2.*G4UniformRand()-1.0;
347   G4double sintheta = std::sqrt((1.0 - costhet    366   G4double sintheta = std::sqrt((1.0 - costheta)*(1.0 + costheta));
348   G4double phi      = 2.0*pi*G4UniformRand()*r    367   G4double phi      = 2.0*pi*G4UniformRand()*rad;
349   G4ThreeVector direction(sintheta*std::cos(ph    368   G4ThreeVector direction(sintheta*std::cos(phi),sintheta*std::sin(phi),costheta);
350   G4DynamicParticle *dynamicNucleon =             369   G4DynamicParticle *dynamicNucleon =
351     new G4DynamicParticle(typeNucleon, directi    370     new G4DynamicParticle(typeNucleon, direction*pp);
352   dynamicNucleon->Set4Momentum(dynamicNucleon-    371   dynamicNucleon->Set4Momentum(dynamicNucleon->Get4Momentum().boost(-boost));
353   G4DynamicParticle *dynamicDaughter =            372   G4DynamicParticle *dynamicDaughter =
354     new G4DynamicParticle(typeDaughter, -direc    373     new G4DynamicParticle(typeDaughter, -direction*pp);
355   dynamicDaughter->Set4Momentum(dynamicDaughte    374   dynamicDaughter->Set4Momentum(dynamicDaughter->Get4Momentum().boost(-boost));
356                                                   375 
357   // The "decay" products have to be transferr    376   // The "decay" products have to be transferred to the G4HadFinalState object.
358   // Furthermore, the residual nucleus should     377   // Furthermore, the residual nucleus should be de-excited.
359                                                   378 
360   theParticleChange.AddSecondary (dynamicNucle << 379   theParticleChange.AddSecondary (dynamicNucleon);
361   if (verboseLevel >= 2) {                        380   if (verboseLevel >= 2) {
362     G4cout <<"Nucleon from the EMD process:" <    381     G4cout <<"Nucleon from the EMD process:" <<G4endl;
363     dynamicNucleon->DumpInfo();                   382     dynamicNucleon->DumpInfo();
364   }                                               383   }
365                                                   384 
366   G4Fragment* theFragment = new                   385   G4Fragment* theFragment = new
367     G4Fragment(typeDaughter->GetBaryonNumber() << 386     G4Fragment((G4int) typeDaughter->GetBaryonNumber(),
368          G4lrint(typeDaughter->GetPDGCharge()/ << 387     (G4int) typeDaughter->GetPDGCharge(), dynamicDaughter->Get4Momentum());
369          dynamicDaughter->Get4Momentum());     << 
370                                                   388 
371   if (verboseLevel >= 2) {                        389   if (verboseLevel >= 2) {
372     G4cout <<"Dynamic properties of the prefra    390     G4cout <<"Dynamic properties of the prefragment:" <<G4endl;
373     G4cout.precision(6);                          391     G4cout.precision(6);
374     dynamicDaughter->DumpInfo();                  392     dynamicDaughter->DumpInfo();
375     G4cout <<"Nuclear properties of the prefra    393     G4cout <<"Nuclear properties of the prefragment:" <<G4endl;
376     G4cout <<theFragment <<G4endl;                394     G4cout <<theFragment <<G4endl;
377   }                                               395   }
378                                                   396 
379   G4ReactionProductVector* products =             397   G4ReactionProductVector* products =
380                       theExcitationHandler->Br    398                       theExcitationHandler->BreakItUp(*theFragment);
381   delete theFragment;                             399   delete theFragment;
382   theFragment = NULL;                             400   theFragment = NULL;
383                                                   401   
384   G4DynamicParticle* secondary = 0;               402   G4DynamicParticle* secondary = 0;
385   G4ReactionProductVector::iterator iter;         403   G4ReactionProductVector::iterator iter;
386   for (iter = products->begin(); iter != produ    404   for (iter = products->begin(); iter != products->end(); ++iter) {
387     secondary = new G4DynamicParticle((*iter)-    405     secondary = new G4DynamicParticle((*iter)->GetDefinition(),
388     (*iter)->GetTotalEnergy(), (*iter)->GetMom    406     (*iter)->GetTotalEnergy(), (*iter)->GetMomentum());
389     theParticleChange.AddSecondary (secondary, << 407     theParticleChange.AddSecondary (secondary);
390   }                                               408   }
391   delete products;                                409   delete products;
392                                                << 
393   delete crossSectionP;                        << 
394   delete crossSectionT;                        << 
395                                                   410 
396   if (verboseLevel >= 2)                          411   if (verboseLevel >= 2)
397     G4cout <<"################################    412     G4cout <<"########################################"
398            <<"################################    413            <<"########################################"
399            <<G4endl;                              414            <<G4endl;
400                                                   415  
401   return &theParticleChange;                      416   return &theParticleChange;
402 }                                                 417 }
403                                                   418 
404                                                   419 
405 void G4EMDissociation::PrintWelcomeMessage ()     420 void G4EMDissociation::PrintWelcomeMessage ()
406 {                                                 421 {
407   G4cout <<G4endl;                                422   G4cout <<G4endl;
408   G4cout <<" *********************************    423   G4cout <<" ****************************************************************"
409          <<G4endl;                                424          <<G4endl;
410   G4cout <<" EM dissociation model for nuclear    425   G4cout <<" EM dissociation model for nuclear-nuclear interactions activated"
411          <<G4endl;                                426          <<G4endl;
412   G4cout <<" (Written by QinetiQ Ltd for the E    427   G4cout <<" (Written by QinetiQ Ltd for the European Space Agency)"
413          <<G4endl;                                428          <<G4endl;
414   G4cout <<" *********************************    429   G4cout <<" ****************************************************************"
415          <<G4endl;                                430          <<G4endl;
416   G4cout << G4endl;                               431   G4cout << G4endl;
417                                                   432 
418   return;                                         433   return;
419 }                                                 434 }
420                                                   435 
421                                                   436