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.3)


  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"
 70 #include "G4IonTable.hh"                           73 #include "G4IonTable.hh"
 71 #include "G4DecayProducts.hh"                      74 #include "G4DecayProducts.hh"
 72 #include "G4DynamicParticle.hh"                    75 #include "G4DynamicParticle.hh"
 73 #include "G4Fragment.hh"                           76 #include "G4Fragment.hh"
 74 #include "G4ReactionProductVector.hh"              77 #include "G4ReactionProductVector.hh"
 75 #include "Randomize.hh"                            78 #include "Randomize.hh"
 76 #include "globals.hh"                              79 #include "globals.hh"
 77 #include "G4PhysicsModelCatalog.hh"            << 
 78                                                    80 
 79 G4EMDissociation::G4EMDissociation() :         <<  81 G4EMDissociation::G4EMDissociation():G4HadronicInteraction("EMDissociation") {
 80   G4HadronicInteraction("EMDissociation"),     <<  82 
 81   secID_projectileDissociation(-1), secID_targ << 
 82 {                                              << 
 83   // Send message to stdout to advise that the     83   // Send message to stdout to advise that the G4EMDissociation model is being
 84   // used.                                         84   // used.
 85   PrintWelcomeMessage();                           85   PrintWelcomeMessage();
 86                                                    86 
 87   // No de-excitation handler has been supplie     87   // No de-excitation handler has been supplied - define the default handler.
 88   theExcitationHandler            = new G4Exci     88   theExcitationHandler            = new G4ExcitationHandler;
                                                   >>  89   G4Evaporation* theEvaporation   = new G4Evaporation;
                                                   >>  90   G4FermiBreakUp* theFermiBreakUp = new G4FermiBreakUp;
                                                   >>  91   G4StatMF* theMF                 = new G4StatMF;
                                                   >>  92   theExcitationHandler->SetEvaporation(theEvaporation);
                                                   >>  93   theExcitationHandler->SetFermiModel(theFermiBreakUp);
                                                   >>  94   theExcitationHandler->SetMultiFragmentation(theMF);
                                                   >>  95   theExcitationHandler->SetMaxAandZForFermiBreakUp(12, 6);
 89   theExcitationHandler->SetMinEForMultiFrag(5.     96   theExcitationHandler->SetMinEForMultiFrag(5.0*MeV);
 90   handlerDefinedInternally = true;                 97   handlerDefinedInternally = true;
 91                                                    98 
 92   // This EM dissociation model needs access t     99   // This EM dissociation model needs access to the cross-sections held in
 93   // G4EMDissociationCrossSection.                100   // G4EMDissociationCrossSection.
 94   dissociationCrossSection = new G4EMDissociat    101   dissociationCrossSection = new G4EMDissociationCrossSection;
 95   thePhotonSpectrum = new G4EMDissociationSpec    102   thePhotonSpectrum = new G4EMDissociationSpectrum;
 96                                                   103 
 97   // Set the minimum and maximum range for the    104   // Set the minimum and maximum range for the model (despite nomanclature, this
 98   // is in energy per nucleon number).            105   // is in energy per nucleon number).    
 99   SetMinEnergy(100.0*MeV);                        106   SetMinEnergy(100.0*MeV);
100   SetMaxEnergy(500.0*GeV);                        107   SetMaxEnergy(500.0*GeV);
101                                                   108 
102   // Set the default verbose level to 0 - no o    109   // Set the default verbose level to 0 - no output.
103   verboseLevel = 0;                               110   verboseLevel = 0;
                                                   >> 111 }
                                                   >> 112 
                                                   >> 113 /*
                                                   >> 114 G4EMDissociation::G4EMDissociation(const G4EMDissociation& emd)
                                                   >> 115  : G4HadronicInteraction(emd)
                                                   >> 116 {
                                                   >> 117   if (emd.theExcitationHandler != 0) {
                                                   >> 118     theExcitationHandler = new G4ExcitationHandler;
                                                   >> 119     *theExcitationHandler = *emd.theExcitationHandler;
                                                   >> 120   }
                                                   >> 121 
                                                   >> 122   handlerDefinedInternally = emd.handlerDefinedInternally;
104                                                   123 
105   // Creator model ID for the secondaries crea << 124   if (emd.dissociationCrossSection != 0) {
106   secID_projectileDissociation = G4PhysicsMode << 125     dissociationCrossSection = new G4EMDissociationCrossSection;
107   secID_targetDissociation     = G4PhysicsMode << 126     *dissociationCrossSection = *emd.dissociationCrossSection;
                                                   >> 127   }
                                                   >> 128 
                                                   >> 129   if (emd.thePhotonSpectrum !- 0) {
                                                   >> 130     thePhotonSpectrum = new G4EMDissociationSpectrum;
                                                   >> 131     *thePhotonSpectrum = *emd.thePhotonSpectrum;
108 }                                                 132 }
                                                   >> 133 */
109                                                   134 
110 G4EMDissociation::G4EMDissociation (G4Excitati << 135 G4EMDissociation::G4EMDissociation (G4ExcitationHandler *aExcitationHandler)
111   G4HadronicInteraction("EMDissociation"),     << 
112   secID_projectileDissociation(-1), secID_targ << 
113 {                                                 136 {
114   // Send message to stdout to advise that the    137   // Send message to stdout to advise that the G4EMDissociation model is being
115   // used.                                        138   // used.
116   PrintWelcomeMessage();                          139   PrintWelcomeMessage();
117                                                   140   
118   theExcitationHandler     = aExcitationHandle    141   theExcitationHandler     = aExcitationHandler;
119   handlerDefinedInternally = false;               142   handlerDefinedInternally = false;
120                                                   143 
121   // This EM dissociation model needs access t    144   // This EM dissociation model needs access to the cross-sections held in
122   // G4EMDissociationCrossSection.                145   // G4EMDissociationCrossSection.
123   dissociationCrossSection = new G4EMDissociat    146   dissociationCrossSection = new G4EMDissociationCrossSection;
124   thePhotonSpectrum = new G4EMDissociationSpec    147   thePhotonSpectrum = new G4EMDissociationSpectrum;
125                                                   148 
126   // Set the minimum and maximum range for the    149   // Set the minimum and maximum range for the model (despite nomanclature, this
127   // is in energy per nucleon number)             150   // is in energy per nucleon number)    
128   SetMinEnergy(100.0*MeV);                        151   SetMinEnergy(100.0*MeV);
129   SetMaxEnergy(500.0*GeV);                        152   SetMaxEnergy(500.0*GeV);
130   verboseLevel = 0;                               153   verboseLevel = 0;
131                                                << 
132   // Creator model ID for the secondaries crea << 
133   secID_projectileDissociation = G4PhysicsMode << 
134   secID_targetDissociation     = G4PhysicsMode << 
135 }                                                 154 }
136                                                   155 
137                                                   156 
138 G4EMDissociation::~G4EMDissociation() {           157 G4EMDissociation::~G4EMDissociation() {
139   if (handlerDefinedInternally) delete theExci    158   if (handlerDefinedInternally) delete theExcitationHandler;
140   // delete dissociationCrossSection;             159   // delete dissociationCrossSection;
141   // Cross section deleted by G4CrossSectionRe    160   // Cross section deleted by G4CrossSectionRegistry; don't do it here
142   // Bug reported by Gong Ding in Bug Report #    161   // Bug reported by Gong Ding in Bug Report #1339
143   delete thePhotonSpectrum;                       162   delete thePhotonSpectrum;
144 }                                                 163 }
145                                                   164 
146                                                   165 
147 G4HadFinalState *G4EMDissociation::ApplyYourse    166 G4HadFinalState *G4EMDissociation::ApplyYourself
148   (const G4HadProjectile &theTrack, G4Nucleus     167   (const G4HadProjectile &theTrack, G4Nucleus &theTarget)
149 {                                                 168 {
150   // The secondaries will be returned in G4Had    169   // The secondaries will be returned in G4HadFinalState &theParticleChange -
151   // initialise this.                             170   // initialise this.
152                                                   171 
153   theParticleChange.Clear();                      172   theParticleChange.Clear();
154   theParticleChange.SetStatusChange(stopAndKil    173   theParticleChange.SetStatusChange(stopAndKill);
155                                                   174 
156   // Get relevant information about the projec    175   // Get relevant information about the projectile and target (A, Z) and
157   // energy/nuc, momentum, velocity, Lorentz f    176   // energy/nuc, momentum, velocity, Lorentz factor and rest-mass of the
158   // projectile.                                  177   // projectile.
159                                                   178 
160   const G4ParticleDefinition *definitionP = th    179   const G4ParticleDefinition *definitionP = theTrack.GetDefinition();
161   const G4double AP  = definitionP->GetBaryonN    180   const G4double AP  = definitionP->GetBaryonNumber();
162   const G4double ZP  = definitionP->GetPDGChar    181   const G4double ZP  = definitionP->GetPDGCharge();
163   G4LorentzVector pP = theTrack.Get4Momentum()    182   G4LorentzVector pP = theTrack.Get4Momentum();
164   G4double E         = theTrack.GetKineticEner    183   G4double E         = theTrack.GetKineticEnergy()/AP;
165   G4double MP        = theTrack.GetTotalEnergy    184   G4double MP        = theTrack.GetTotalEnergy() - E*AP;
166   G4double b         = pP.beta();                 185   G4double b         = pP.beta();
167   G4double AT        = theTarget.GetA_asInt();    186   G4double AT        = theTarget.GetA_asInt();
168   G4double ZT        = theTarget.GetZ_asInt();    187   G4double ZT        = theTarget.GetZ_asInt();
169   G4double MT        = G4NucleiProperties::Get    188   G4double MT        = G4NucleiProperties::GetNuclearMass(AT,ZT);
170                                                   189 
171   // Depending upon the verbosity level, outpu    190   // Depending upon the verbosity level, output the initial information on the
172   // projectile and target                        191   // projectile and target
173   if (verboseLevel >= 2) {                        192   if (verboseLevel >= 2) {
174     G4cout.precision(6);                          193     G4cout.precision(6);
175     G4cout <<"################################    194     G4cout <<"########################################"
176            <<"################################    195            <<"########################################"
177            <<G4endl;                              196            <<G4endl;
178     G4cout <<"IN G4EMDissociation" <<G4endl;      197     G4cout <<"IN G4EMDissociation" <<G4endl;
179     G4cout <<"Initial projectile A=" <<AP         198     G4cout <<"Initial projectile A=" <<AP 
180            <<", Z=" <<ZP                          199            <<", Z=" <<ZP
181            <<G4endl;                              200            <<G4endl; 
182     G4cout <<"Initial target     A=" <<AT         201     G4cout <<"Initial target     A=" <<AT
183            <<", Z=" <<ZT                          202            <<", Z=" <<ZT
184            <<G4endl;                              203            <<G4endl;
185     G4cout <<"Projectile momentum and Energy/n    204     G4cout <<"Projectile momentum and Energy/nuc = " <<pP <<" ," <<E <<G4endl;
186   }                                               205   }
187                                                   206 
188   // Initialise the variables which will be us    207   // Initialise the variables which will be used with the phase-space decay and
189   // to boost the secondaries from the interac    208   // to boost the secondaries from the interaction.
190                                                   209   
191   G4ParticleDefinition *typeNucleon  = NULL;      210   G4ParticleDefinition *typeNucleon  = NULL;
192   G4ParticleDefinition *typeDaughter = NULL;      211   G4ParticleDefinition *typeDaughter = NULL;
193   G4double Eg                        = 0.0;       212   G4double Eg                        = 0.0;
194   G4double mass                      = 0.0;       213   G4double mass                      = 0.0;
195   G4ThreeVector boost = G4ThreeVector(0.0, 0.0    214   G4ThreeVector boost = G4ThreeVector(0.0, 0.0, 0.0);
196                                                   215 
197   // Determine the cross-sections at the giant    216   // Determine the cross-sections at the giant dipole and giant quadrupole
198   // resonance energies for the projectile and    217   // resonance energies for the projectile and then target.  The information is
199   // initially provided in the G4PhysicsFreeVe    218   // initially provided in the G4PhysicsFreeVector individually for the E1
200   // and E2 fields. These are then summed.        219   // and E2 fields. These are then summed.
201                                                   220 
202   G4double bmin = thePhotonSpectrum->GetCloses    221   G4double bmin = thePhotonSpectrum->GetClosestApproach(AP, ZP, AT, ZT, b);
203   G4PhysicsFreeVector *crossSectionP = dissoci    222   G4PhysicsFreeVector *crossSectionP = dissociationCrossSection->
204     GetCrossSectionForProjectile(AP, ZP, AT, Z    223     GetCrossSectionForProjectile(AP, ZP, AT, ZT, b, bmin);
205   G4PhysicsFreeVector *crossSectionT = dissoci    224   G4PhysicsFreeVector *crossSectionT = dissociationCrossSection->
206     GetCrossSectionForTarget(AP, ZP, AT, ZT, b    225     GetCrossSectionForTarget(AP, ZP, AT, ZT, b, bmin);
207                                                   226 
208   G4double totCrossSectionP = (*crossSectionP)    227   G4double totCrossSectionP = (*crossSectionP)[0]+(*crossSectionP)[1];
209   G4double totCrossSectionT = (*crossSectionT)    228   G4double totCrossSectionT = (*crossSectionT)[0]+(*crossSectionT)[1];
210                                                   229 
211   // Now sample whether the interaction involv    230   // Now sample whether the interaction involved EM dissociation of the projectile
212   // or the target.                               231   // or the target.
213                                                << 232 
214   G4int secID = -1;  // Creator model ID for t << 
215   if (G4UniformRand() <                           233   if (G4UniformRand() <
216     totCrossSectionP / (totCrossSectionP + tot    234     totCrossSectionP / (totCrossSectionP + totCrossSectionT)) {
217                                                   235 
218     // It was the projectile which underwent E    236     // It was the projectile which underwent EM dissociation.  Define the Lorentz
219     // boost to be applied to the secondaries,    237     // boost to be applied to the secondaries, and sample whether a proton or a
220     // neutron was ejected.  Then determine th    238     // neutron was ejected.  Then determine the energy of the virtual gamma ray
221     // which passed from the target nucleus ..    239     // which passed from the target nucleus ... this will be used to define the
222     // excitation of the projectile.              240     // excitation of the projectile.
223                                                   241 
224     secID = secID_projectileDissociation;      << 
225     mass  = MP;                                   242     mass  = MP;
226     if (G4UniformRand() < dissociationCrossSec    243     if (G4UniformRand() < dissociationCrossSection->
227       GetWilsonProbabilityForProtonDissociatio    244       GetWilsonProbabilityForProtonDissociation (AP, ZP))
228     {                                             245     {
229       if (verboseLevel >= 2)                      246       if (verboseLevel >= 2)
230         G4cout <<"Projectile underwent EM diss    247         G4cout <<"Projectile underwent EM dissociation producing a proton"
231                <<G4endl;                          248                <<G4endl;
232       typeNucleon = G4Proton::ProtonDefinition    249       typeNucleon = G4Proton::ProtonDefinition();
233       typeDaughter = G4IonTable::GetIonTable()    250       typeDaughter = G4IonTable::GetIonTable()->
234       GetIon((G4int) ZP-1, (G4int) AP-1, 0.0);    251       GetIon((G4int) ZP-1, (G4int) AP-1, 0.0);
235     }                                             252     }
236     else                                          253     else
237     {                                             254     {
238       if (verboseLevel >= 2)                      255       if (verboseLevel >= 2)
239         G4cout <<"Projectile underwent EM diss    256         G4cout <<"Projectile underwent EM dissociation producing a neutron"
240                <<G4endl;                          257                <<G4endl;
241       typeNucleon = G4Neutron::NeutronDefiniti    258       typeNucleon = G4Neutron::NeutronDefinition();
242       typeDaughter = G4IonTable::GetIonTable()    259       typeDaughter = G4IonTable::GetIonTable()->
243       GetIon((G4int) ZP, (G4int) AP-1, 0.0);      260       GetIon((G4int) ZP, (G4int) AP-1, 0.0);
244     }                                             261     }
245     if (G4UniformRand() < (*crossSectionP)[0]/    262     if (G4UniformRand() < (*crossSectionP)[0]/totCrossSectionP)
246     {                                             263     {
247       Eg = crossSectionP->GetLowEdgeEnergy(0);    264       Eg = crossSectionP->GetLowEdgeEnergy(0);
248       if (verboseLevel >= 2)                      265       if (verboseLevel >= 2)
249         G4cout <<"Transition type was E1" <<G4    266         G4cout <<"Transition type was E1" <<G4endl;
250     }                                             267     }
251     else                                          268     else
252     {                                             269     {
253       Eg = crossSectionP->GetLowEdgeEnergy(1);    270       Eg = crossSectionP->GetLowEdgeEnergy(1);
254       if (verboseLevel >= 2)                      271       if (verboseLevel >= 2)
255         G4cout <<"Transition type was E2" <<G4    272         G4cout <<"Transition type was E2" <<G4endl;
256     }                                             273     }
257                                                   274 
258     // We need to define a Lorentz vector with    275     // We need to define a Lorentz vector with the original momentum, but total
259     // energy includes the projectile and virt    276     // energy includes the projectile and virtual gamma.  This is then used
260     // to calculate the boost required for the    277     // to calculate the boost required for the secondaries.
261                                                   278 
262     pP.setE( std::sqrt( pP.vect().mag2() + (ma << 279     pP.setE(pP.e()+Eg);
263     boost = pP.findBoostToCM();                   280     boost = pP.findBoostToCM();
264   }                                               281   }
265   else                                            282   else
266   {                                               283   {
267     // It was the target which underwent EM di    284     // It was the target which underwent EM dissociation.  Sample whether a
268     // proton or a neutron was ejected.  Then     285     // proton or a neutron was ejected.  Then determine the energy of the virtual 
269     // gamma ray which passed from the project    286     // gamma ray which passed from the projectile nucleus ... this will be used to
270     // define the excitation of the target.       287     // define the excitation of the target.
271                                                << 288 
272     secID = secID_targetDissociation;          << 
273     mass = MT;                                    289     mass = MT;
274     if (G4UniformRand() < dissociationCrossSec    290     if (G4UniformRand() < dissociationCrossSection->
275       GetWilsonProbabilityForProtonDissociatio    291       GetWilsonProbabilityForProtonDissociation (AT, ZT))
276     {                                             292     {
277       if (verboseLevel >= 2)                      293       if (verboseLevel >= 2)
278         G4cout <<"Target underwent EM dissocia    294         G4cout <<"Target underwent EM dissociation producing a proton"
279                <<G4endl;                          295                <<G4endl;
280       typeNucleon = G4Proton::ProtonDefinition    296       typeNucleon = G4Proton::ProtonDefinition();
281       typeDaughter = G4IonTable::GetIonTable()    297       typeDaughter = G4IonTable::GetIonTable()->
282       GetIon((G4int) ZT-1, (G4int) AT-1, 0.0);    298       GetIon((G4int) ZT-1, (G4int) AT-1, 0.0);
283     }                                             299     }
284     else                                          300     else
285     {                                             301     {
286       if (verboseLevel >= 2)                      302       if (verboseLevel >= 2)
287         G4cout <<"Target underwent EM dissocia    303         G4cout <<"Target underwent EM dissociation producing a neutron"
288                <<G4endl;                          304                <<G4endl;
289       typeNucleon = G4Neutron::NeutronDefiniti    305       typeNucleon = G4Neutron::NeutronDefinition();
290       typeDaughter = G4IonTable::GetIonTable()    306       typeDaughter = G4IonTable::GetIonTable()->
291       GetIon((G4int) ZT, (G4int) AT-1, 0.0);      307       GetIon((G4int) ZT, (G4int) AT-1, 0.0);
292     }                                             308     }
293     if (G4UniformRand() < (*crossSectionT)[0]/    309     if (G4UniformRand() < (*crossSectionT)[0]/totCrossSectionT)
294     {                                             310     {
295       Eg = crossSectionT->GetLowEdgeEnergy(0);    311       Eg = crossSectionT->GetLowEdgeEnergy(0);
296       if (verboseLevel >= 2)                      312       if (verboseLevel >= 2)
297         G4cout <<"Transition type was E1" <<G4    313         G4cout <<"Transition type was E1" <<G4endl;
298     }                                             314     }
299     else                                          315     else
300     {                                             316     {
301       Eg = crossSectionT->GetLowEdgeEnergy(1);    317       Eg = crossSectionT->GetLowEdgeEnergy(1);
302       if (verboseLevel >= 2)                      318       if (verboseLevel >= 2)
303         G4cout <<"Transition type was E2" <<G4    319         G4cout <<"Transition type was E2" <<G4endl;
304     }                                             320     }
305                                                   321 
306     // Add the projectile to theParticleChange    322     // Add the projectile to theParticleChange, less the energy of the
307     // not-so-virtual gamma-ray.  Not that at     323     // not-so-virtual gamma-ray.  Not that at the moment, no lateral momentum
308     // is transferred between the projectile a    324     // is transferred between the projectile and target nuclei.
309                                                   325 
310     G4ThreeVector v = pP.vect();                  326     G4ThreeVector v = pP.vect();
311     v.setMag(1.0);                                327     v.setMag(1.0);
312     G4DynamicParticle *changedP = new G4Dynami    328     G4DynamicParticle *changedP = new G4DynamicParticle (definitionP, v, E*AP-Eg);
313     theParticleChange.AddSecondary (changedP,  << 329     theParticleChange.AddSecondary (changedP);
314     if (verboseLevel >= 2)                        330     if (verboseLevel >= 2)
315     {                                             331     {
316       G4cout <<"Projectile change:" <<G4endl;     332       G4cout <<"Projectile change:" <<G4endl;
317       changedP->DumpInfo();                       333       changedP->DumpInfo();
318     }                                             334     }
319   }                                               335   }
320                                                   336 
321   // Perform a two-body decay based on the res    337   // Perform a two-body decay based on the restmass energy of the parent and
322   // gamma-ray, and the masses of the daughter    338   // gamma-ray, and the masses of the daughters. In the frame of reference of
323   // the nucles, the angular distribution is s    339   // the nucles, the angular distribution is sampled isotropically, but the
324   // the nucleon and secondary nucleus are boo    340   // the nucleon and secondary nucleus are boosted if they've come from the
325   // projectile.                                  341   // projectile.
326                                                   342 
327   G4double e  = mass + Eg;                        343   G4double e  = mass + Eg;
328   G4double mass1 = typeNucleon->GetPDGMass();     344   G4double mass1 = typeNucleon->GetPDGMass();
329   G4double mass2 = typeDaughter->GetPDGMass();    345   G4double mass2 = typeDaughter->GetPDGMass();
330   G4double pp = (e+mass1+mass2)*(e+mass1-mass2    346   G4double pp = (e+mass1+mass2)*(e+mass1-mass2)*
331                 (e-mass1+mass2)*(e-mass1-mass2    347                 (e-mass1+mass2)*(e-mass1-mass2)/(4.0*e*e);
332   if (pp < 0.0) {                                 348   if (pp < 0.0) {
333     pp = 1.0*eV;                                  349     pp = 1.0*eV;
334 //    if (verboseLevel >`= 1)                     350 //    if (verboseLevel >`= 1)
335 //    {                                           351 //    {
336 //      G4cout <<"IN G4EMDissociation::ApplyYo    352 //      G4cout <<"IN G4EMDissociation::ApplyYoursef" <<G4endl;
337 //      G4cout <<"Error in mass of secondaries    353 //      G4cout <<"Error in mass of secondaries compared with primary:" <<G4endl;
338 //      G4cout <<"Rest mass of primary      =     354 //      G4cout <<"Rest mass of primary      = " <<mass <<" MeV" <<G4endl;
339 //      G4cout <<"Virtual gamma energy      =     355 //      G4cout <<"Virtual gamma energy      = " <<Eg   <<" MeV" <<G4endl;
340 //      G4cout <<"Rest mass of secondary #1 =     356 //      G4cout <<"Rest mass of secondary #1 = " <<mass1   <<" MeV" <<G4endl;
341 //      G4cout <<"Rest mass of secondary #2 =     357 //      G4cout <<"Rest mass of secondary #2 = " <<mass2   <<" MeV" <<G4endl;
342 //    }                                           358 //    }
343   }                                               359   }
344   else                                            360   else
345     pp = std::sqrt(pp);                           361     pp = std::sqrt(pp);
346   G4double costheta = 2.*G4UniformRand()-1.0;     362   G4double costheta = 2.*G4UniformRand()-1.0;
347   G4double sintheta = std::sqrt((1.0 - costhet    363   G4double sintheta = std::sqrt((1.0 - costheta)*(1.0 + costheta));
348   G4double phi      = 2.0*pi*G4UniformRand()*r    364   G4double phi      = 2.0*pi*G4UniformRand()*rad;
349   G4ThreeVector direction(sintheta*std::cos(ph    365   G4ThreeVector direction(sintheta*std::cos(phi),sintheta*std::sin(phi),costheta);
350   G4DynamicParticle *dynamicNucleon =             366   G4DynamicParticle *dynamicNucleon =
351     new G4DynamicParticle(typeNucleon, directi    367     new G4DynamicParticle(typeNucleon, direction*pp);
352   dynamicNucleon->Set4Momentum(dynamicNucleon-    368   dynamicNucleon->Set4Momentum(dynamicNucleon->Get4Momentum().boost(-boost));
353   G4DynamicParticle *dynamicDaughter =            369   G4DynamicParticle *dynamicDaughter =
354     new G4DynamicParticle(typeDaughter, -direc    370     new G4DynamicParticle(typeDaughter, -direction*pp);
355   dynamicDaughter->Set4Momentum(dynamicDaughte    371   dynamicDaughter->Set4Momentum(dynamicDaughter->Get4Momentum().boost(-boost));
356                                                   372 
357   // The "decay" products have to be transferr    373   // The "decay" products have to be transferred to the G4HadFinalState object.
358   // Furthermore, the residual nucleus should     374   // Furthermore, the residual nucleus should be de-excited.
359                                                   375 
360   theParticleChange.AddSecondary (dynamicNucle << 376   theParticleChange.AddSecondary (dynamicNucleon);
361   if (verboseLevel >= 2) {                        377   if (verboseLevel >= 2) {
362     G4cout <<"Nucleon from the EMD process:" <    378     G4cout <<"Nucleon from the EMD process:" <<G4endl;
363     dynamicNucleon->DumpInfo();                   379     dynamicNucleon->DumpInfo();
364   }                                               380   }
365                                                   381 
366   G4Fragment* theFragment = new                   382   G4Fragment* theFragment = new
367     G4Fragment(typeDaughter->GetBaryonNumber() << 383     G4Fragment((G4int) typeDaughter->GetBaryonNumber(),
368          G4lrint(typeDaughter->GetPDGCharge()/ << 384     (G4int) typeDaughter->GetPDGCharge(), dynamicDaughter->Get4Momentum());
369          dynamicDaughter->Get4Momentum());     << 
370                                                   385 
371   if (verboseLevel >= 2) {                        386   if (verboseLevel >= 2) {
372     G4cout <<"Dynamic properties of the prefra    387     G4cout <<"Dynamic properties of the prefragment:" <<G4endl;
373     G4cout.precision(6);                          388     G4cout.precision(6);
374     dynamicDaughter->DumpInfo();                  389     dynamicDaughter->DumpInfo();
375     G4cout <<"Nuclear properties of the prefra    390     G4cout <<"Nuclear properties of the prefragment:" <<G4endl;
376     G4cout <<theFragment <<G4endl;                391     G4cout <<theFragment <<G4endl;
377   }                                               392   }
378                                                   393 
379   G4ReactionProductVector* products =             394   G4ReactionProductVector* products =
380                       theExcitationHandler->Br    395                       theExcitationHandler->BreakItUp(*theFragment);
381   delete theFragment;                             396   delete theFragment;
382   theFragment = NULL;                             397   theFragment = NULL;
383                                                   398   
384   G4DynamicParticle* secondary = 0;               399   G4DynamicParticle* secondary = 0;
385   G4ReactionProductVector::iterator iter;         400   G4ReactionProductVector::iterator iter;
386   for (iter = products->begin(); iter != produ    401   for (iter = products->begin(); iter != products->end(); ++iter) {
387     secondary = new G4DynamicParticle((*iter)-    402     secondary = new G4DynamicParticle((*iter)->GetDefinition(),
388     (*iter)->GetTotalEnergy(), (*iter)->GetMom    403     (*iter)->GetTotalEnergy(), (*iter)->GetMomentum());
389     theParticleChange.AddSecondary (secondary, << 404     theParticleChange.AddSecondary (secondary);
390   }                                               405   }
391   delete products;                                406   delete products;
392                                                   407 
393   delete crossSectionP;                           408   delete crossSectionP;
394   delete crossSectionT;                           409   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