G4EMDissociation.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // *                                                                  *
00021 // * Parts of this code which have been  developed by QinetiQ Ltd     *
00022 // * under contract to the European Space Agency (ESA) are the        *
00023 // * intellectual property of ESA. Rights to use, copy, modify and    *
00024 // * redistribute this software for general public use are granted    *
00025 // * in compliance with any licensing, distribution and development   *
00026 // * policy adopted by the Geant4 Collaboration. This code has been   *
00027 // * written by QinetiQ Ltd for the European Space Agency, under ESA  *
00028 // * contract 17191/03/NL/LvH (Aurora Programme).                     *
00029 // *                                                                  *
00030 // * By using,  copying,  modifying or  distributing the software (or *
00031 // * any work based  on the software)  you  agree  to acknowledge its *
00032 // * use  in  resulting  scientific  publications,  and indicate your *
00033 // * acceptance of all terms of the Geant4 Software license.          *
00034 // ********************************************************************
00035 //
00036 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00037 //
00038 // MODULE:              G4EMDissociation.cc
00039 //
00040 // Version:             B.1
00041 // Date:                15/04/04
00042 // Author:              P R Truscott
00043 // Organisation:        QinetiQ Ltd, UK
00044 // Customer:            ESA/ESTEC, NOORDWIJK
00045 // Contract:            17191/03/NL/LvH
00046 //
00047 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00048 //
00049 // CHANGE HISTORY
00050 // --------------
00051 //
00052 // 17 October 2003, P R Truscott, QinetiQ Ltd, UK
00053 // Created.
00054 //
00055 // 15 March 2004, P R Truscott, QinetiQ Ltd, UK
00056 // Beta release
00057 //
00058 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00060 //
00061 #include "G4EMDissociation.hh"
00062 #include "G4PhysicalConstants.hh"
00063 #include "G4SystemOfUnits.hh"
00064 #include "G4Evaporation.hh"
00065 #include "G4FermiBreakUp.hh"
00066 #include "G4StatMF.hh"
00067 #include "G4ParticleDefinition.hh"
00068 #include "G4LorentzVector.hh"
00069 #include "G4PhysicsFreeVector.hh"
00070 #include "G4EMDissociationCrossSection.hh"
00071 #include "G4Proton.hh"
00072 #include "G4Neutron.hh"
00073 #include "G4ParticleTable.hh"
00074 #include "G4IonTable.hh"
00075 #include "G4GeneralPhaseSpaceDecay.hh"
00076 #include "G4DecayProducts.hh"
00077 #include "G4DynamicParticle.hh"
00078 #include "G4Fragment.hh"
00079 #include "G4ReactionProductVector.hh"
00080 #include "Randomize.hh"
00081 #include "globals.hh"
00082 
00083 G4EMDissociation::G4EMDissociation():G4HadronicInteraction("EMDissociation") {
00084 
00085   // Send message to stdout to advise that the G4EMDissociation model is being
00086   // used.
00087   PrintWelcomeMessage();
00088 
00089   // No de-excitation handler has been supplied - define the default handler.
00090   theExcitationHandler            = new G4ExcitationHandler;
00091   G4Evaporation* theEvaporation   = new G4Evaporation;
00092   G4FermiBreakUp* theFermiBreakUp = new G4FermiBreakUp;
00093   G4StatMF* theMF                 = new G4StatMF;
00094   theExcitationHandler->SetEvaporation(theEvaporation);
00095   theExcitationHandler->SetFermiModel(theFermiBreakUp);
00096   theExcitationHandler->SetMultiFragmentation(theMF);
00097   theExcitationHandler->SetMaxAandZForFermiBreakUp(12, 6);
00098   theExcitationHandler->SetMinEForMultiFrag(5.0*MeV);
00099   handlerDefinedInternally = true;
00100 
00101   // This EM dissociation model needs access to the cross-sections held in
00102   // G4EMDissociationCrossSection.
00103   dissociationCrossSection = new G4EMDissociationCrossSection;
00104   thePhotonSpectrum = new G4EMDissociationSpectrum;
00105 
00106   // Set the minimum and maximum range for the model (despite nomanclature, this
00107   // is in energy per nucleon number).    
00108   SetMinEnergy(100.0*MeV);
00109   SetMaxEnergy(500.0*GeV);
00110 
00111   // Set the default verbose level to 0 - no output.
00112   verboseLevel = 0;
00113 }
00114 
00115 /*
00116 G4EMDissociation::G4EMDissociation(const G4EMDissociation& emd)
00117  : G4HadronicInteraction(emd)
00118 {
00119   if (emd.theExcitationHandler != 0) {
00120     theExcitationHandler = new G4ExcitationHandler;
00121     *theExcitationHandler = *emd.theExcitationHandler;
00122   }
00123 
00124   handlerDefinedInternally = emd.handlerDefinedInternally;
00125 
00126   if (emd.dissociationCrossSection != 0) {
00127     dissociationCrossSection = new G4EMDissociationCrossSection;
00128     *dissociationCrossSection = *emd.dissociationCrossSection;
00129   }
00130 
00131   if (emd.thePhotonSpectrum !- 0) {
00132     thePhotonSpectrum = new G4EMDissociationSpectrum;
00133     *thePhotonSpectrum = *emd.thePhotonSpectrum;
00134 }
00135 */
00136 
00137 G4EMDissociation::G4EMDissociation (G4ExcitationHandler *aExcitationHandler)
00138 {
00139 //
00140 //
00141 // Send message to stdout to advise that the G4EMDissociation model is being
00142 // used.
00143 //
00144   PrintWelcomeMessage();
00145   
00146   theExcitationHandler     = aExcitationHandler;
00147   handlerDefinedInternally = false;
00148 //
00149 //
00150 // This EM dissociation model needs access to the cross-sections held in
00151 // G4EMDissociationCrossSection.
00152 //  
00153   dissociationCrossSection = new G4EMDissociationCrossSection;
00154   thePhotonSpectrum = new G4EMDissociationSpectrum;
00155 //
00156 //
00157 // Set the minimum and maximum range for the model (despite nomanclature, this
00158 // is in energy per nucleon number).  
00159 //    
00160   SetMinEnergy(100.0*MeV);
00161   SetMaxEnergy(500.0*GeV);
00162 //
00163 //
00164 // Set the default verbose level to 0 - no output.
00165 //
00166   verboseLevel = 0;
00167 }
00168 
00169 
00170 G4EMDissociation::~G4EMDissociation() {
00171   if (handlerDefinedInternally) delete theExcitationHandler;
00172   // delete dissociationCrossSection;
00173   // Cross section deleted by G4CrossSectionRegistry; don't do it here
00174   // Bug reported by Gong Ding in Bug Report #1339
00175   delete thePhotonSpectrum;
00176 }
00177 
00178 
00179 G4HadFinalState *G4EMDissociation::ApplyYourself
00180   (const G4HadProjectile &theTrack, G4Nucleus &theTarget)
00181 {
00182 //
00183 //
00184 // The secondaries will be returned in G4HadFinalState &theParticleChange -
00185 // initialise this.
00186 //
00187   theParticleChange.Clear();
00188   theParticleChange.SetStatusChange(stopAndKill);
00189 //
00190 //
00191 // Get relevant information about the projectile and target (A, Z) and
00192 // energy/nuc, momentum, velocity, Lorentz factor and rest-mass of the
00193 // projectile.
00194 //
00195   const G4ParticleDefinition *definitionP = theTrack.GetDefinition();
00196   const G4double AP  = definitionP->GetBaryonNumber();
00197   const G4double ZP  = definitionP->GetPDGCharge();
00198   G4LorentzVector pP = theTrack.Get4Momentum();
00199   G4double E         = theTrack.GetKineticEnergy()/AP;
00200   G4double MP        = theTrack.GetTotalEnergy() - E*AP;
00201   G4double b         = pP.beta();
00202   G4double AT        = theTarget.GetA_asInt();
00203   G4double ZT        = theTarget.GetZ_asInt();
00204   G4double MT        = G4NucleiProperties::GetNuclearMass(AT,ZT);
00205 //
00206 //
00207 // Depending upon the verbosity level, output the initial information on the
00208 // projectile and target.
00209 //
00210   if (verboseLevel >= 2)
00211   {
00212     G4cout.precision(6);
00213     G4cout <<"########################################"
00214            <<"########################################"
00215            <<G4endl;
00216     G4cout <<"IN G4EMDissociation" <<G4endl;
00217     G4cout <<"Initial projectile A=" <<AP 
00218            <<", Z=" <<ZP
00219            <<G4endl; 
00220     G4cout <<"Initial target     A=" <<AT
00221            <<", Z=" <<ZT
00222            <<G4endl;
00223     G4cout <<"Projectile momentum and Energy/nuc = " <<pP <<" ," <<E <<G4endl;
00224   }
00225 //
00226 //
00227 // Initialise the variables which will be used with the phase-space decay and
00228 // to boost the secondaries from the interaction.
00229 //  
00230   G4ParticleDefinition *typeNucleon  = NULL;
00231   G4ParticleDefinition *typeDaughter = NULL;
00232   G4double Eg                        = 0.0;
00233   G4double mass                      = 0.0;
00234   G4ThreeVector boost = G4ThreeVector(0.0, 0.0, 0.0);
00235 //
00236 //
00237 // Determine the cross-sections at the giant dipole and giant quadrupole
00238 // resonance energies for the projectile and then target.  The information is
00239 // initially provided in the G4PhysicsFreeVector individually for the E1
00240 // and E2 fields. These are then summed.
00241 //
00242   G4double bmin = thePhotonSpectrum->GetClosestApproach(AP, ZP, AT, ZT, b);
00243   G4PhysicsFreeVector *crossSectionP = dissociationCrossSection->
00244     GetCrossSectionForProjectile(AP, ZP, AT, ZT, b, bmin);
00245   G4PhysicsFreeVector *crossSectionT = dissociationCrossSection->
00246     GetCrossSectionForTarget(AP, ZP, AT, ZT, b, bmin);
00247 
00248   G4double totCrossSectionP = (*crossSectionP)[0]+(*crossSectionP)[1];
00249   G4double totCrossSectionT = (*crossSectionT)[0]+(*crossSectionT)[1];
00250 //
00251 //
00252 // Now sample whether the interaction involved EM dissociation of the projectile
00253 // or the target.
00254 //
00255   if (G4UniformRand() <
00256     totCrossSectionP / (totCrossSectionP + totCrossSectionT))
00257   {
00258 //
00259 //
00260 // It was the projectile which underwent EM dissociation.  Define the Lorentz
00261 // boost to be applied to the secondaries, and sample whether a proton or a
00262 // neutron was ejected.  Then determine the energy of the virtual gamma ray
00263 // which passed from the target nucleus ... this will be used to define the
00264 // excitation of the projectile.
00265 //
00266     mass  = MP;
00267     if (G4UniformRand() < dissociationCrossSection->
00268       GetWilsonProbabilityForProtonDissociation (AP, ZP))
00269     {
00270       if (verboseLevel >= 2)
00271         G4cout <<"Projectile underwent EM dissociation producing a proton"
00272                <<G4endl;
00273       typeNucleon = G4Proton::ProtonDefinition();
00274       typeDaughter = G4ParticleTable::GetParticleTable()->
00275       GetIon((G4int) ZP-1, (G4int) AP-1, 0.0);
00276     }
00277     else
00278     {
00279       if (verboseLevel >= 2)
00280         G4cout <<"Projectile underwent EM dissociation producing a neutron"
00281                <<G4endl;
00282       typeNucleon = G4Neutron::NeutronDefinition();
00283       typeDaughter = G4ParticleTable::GetParticleTable()->
00284       GetIon((G4int) ZP, (G4int) AP-1, 0.0);
00285     }
00286     if (G4UniformRand() < (*crossSectionP)[0]/totCrossSectionP)
00287     {
00288       Eg = crossSectionP->GetLowEdgeEnergy(0);
00289       if (verboseLevel >= 2)
00290         G4cout <<"Transition type was E1" <<G4endl;
00291     }
00292     else
00293     {
00294       Eg = crossSectionP->GetLowEdgeEnergy(1);
00295       if (verboseLevel >= 2)
00296         G4cout <<"Transition type was E2" <<G4endl;
00297     }
00298 //
00299 //
00300 // We need to define a Lorentz vector with the original momentum, but total
00301 // energy includes the projectile and virtual gamma.  This is then used
00302 // to calculate the boost required for the secondaries.
00303 //
00304     pP.setE(pP.e()+Eg);
00305     boost = pP.findBoostToCM();
00306   }
00307   else
00308   {
00309 //
00310 //
00311 // It was the target which underwent EM dissociation.  Sample whether a
00312 // proton or a neutron was ejected.  Then determine the energy of the virtual 
00313 // gamma ray which passed from the projectile nucleus ... this will be used to
00314 // define the excitation of the target.
00315 //
00316     mass = MT;
00317     if (G4UniformRand() < dissociationCrossSection->
00318       GetWilsonProbabilityForProtonDissociation (AT, ZT))
00319     {
00320       if (verboseLevel >= 2)
00321         G4cout <<"Target underwent EM dissociation producing a proton"
00322                <<G4endl;
00323       typeNucleon = G4Proton::ProtonDefinition();
00324       typeDaughter = G4ParticleTable::GetParticleTable()->
00325       GetIon((G4int) ZT-1, (G4int) AT-1, 0.0);
00326     }
00327     else
00328     {
00329       if (verboseLevel >= 2)
00330         G4cout <<"Target underwent EM dissociation producing a neutron"
00331                <<G4endl;
00332       typeNucleon = G4Neutron::NeutronDefinition();
00333       typeDaughter = G4ParticleTable::GetParticleTable()->
00334       GetIon((G4int) ZT, (G4int) AT-1, 0.0);
00335     }
00336     if (G4UniformRand() < (*crossSectionT)[0]/totCrossSectionT)
00337     {
00338       Eg = crossSectionT->GetLowEdgeEnergy(0);
00339       if (verboseLevel >= 2)
00340         G4cout <<"Transition type was E1" <<G4endl;
00341     }
00342     else
00343     {
00344       Eg = crossSectionT->GetLowEdgeEnergy(1);
00345       if (verboseLevel >= 2)
00346         G4cout <<"Transition type was E2" <<G4endl;
00347     }
00348 //
00349 //
00350 // Add the projectile to theParticleChange, less the energy of the
00351 // not-so-virtual gamma-ray.  Not that at the moment, no lateral momentum
00352 // is transferred between the projectile and target nuclei.
00353 //
00354     G4ThreeVector v = pP.vect();
00355     v.setMag(1.0);
00356     G4DynamicParticle *changedP = new G4DynamicParticle
00357       (const_cast<G4ParticleDefinition*>(definitionP), v, E*AP-Eg);
00358     theParticleChange.AddSecondary (changedP);
00359     if (verboseLevel >= 2)
00360     {
00361       G4cout <<"Projectile change:" <<G4endl;
00362       changedP->DumpInfo();
00363     }
00364   }
00365 //
00366 //
00367 // Perform a two-body decay based on the restmass energy of the parent and
00368 // gamma-ray, and the masses of the daughters. In the frame of reference of
00369 // the nucles, the angular distribution is sampled isotropically, but the
00370 // the nucleon and secondary nucleus are boosted if they've come from the
00371 // projectile.
00372 //
00373   G4double e  = mass + Eg;
00374   G4double mass1 = typeNucleon->GetPDGMass();
00375   G4double mass2 = typeDaughter->GetPDGMass();
00376   G4double pp = (e+mass1+mass2)*(e+mass1-mass2)*
00377                 (e-mass1+mass2)*(e-mass1-mass2)/(4.0*e*e);
00378   if (pp < 0.0)
00379   {
00380     pp = 1.0*eV;
00381 //    if (verboseLevel >`= 1)
00382 //    {
00383 //      G4cout <<"IN G4EMDissociation::ApplyYoursef" <<G4endl;
00384 //      G4cout <<"Error in mass of secondaries compared with primary:" <<G4endl;
00385 //      G4cout <<"Rest mass of primary      = " <<mass <<" MeV" <<G4endl;
00386 //      G4cout <<"Virtual gamma energy      = " <<Eg   <<" MeV" <<G4endl;
00387 //      G4cout <<"Rest mass of secondary #1 = " <<mass1   <<" MeV" <<G4endl;
00388 //      G4cout <<"Rest mass of secondary #2 = " <<mass2   <<" MeV" <<G4endl;
00389 //    }
00390   }
00391   else
00392     pp = std::sqrt(pp);
00393   G4double costheta = 2.*G4UniformRand()-1.0;
00394   G4double sintheta = std::sqrt((1.0 - costheta)*(1.0 + costheta));
00395   G4double phi      = 2.0*pi*G4UniformRand()*rad;
00396   G4ThreeVector direction(sintheta*std::cos(phi),sintheta*std::sin(phi),costheta);
00397   G4DynamicParticle *dynamicNucleon =
00398     new G4DynamicParticle(typeNucleon, direction*pp);
00399   dynamicNucleon->Set4Momentum(dynamicNucleon->Get4Momentum().boost(-boost));
00400   G4DynamicParticle *dynamicDaughter =
00401     new G4DynamicParticle(typeDaughter, -direction*pp);
00402   dynamicDaughter->Set4Momentum(dynamicDaughter->Get4Momentum().boost(-boost));
00403 //
00404 //
00405 // The "decay" products have to be transferred to the G4HadFinalState object.
00406 // Furthermore, the residual nucleus should be de-excited.
00407 //
00408   theParticleChange.AddSecondary (dynamicNucleon);
00409   if (verboseLevel >= 2)
00410   {
00411     G4cout <<"Nucleon from the EMD process:" <<G4endl;
00412     dynamicNucleon->DumpInfo();
00413   }
00414 
00415   G4Fragment *theFragment = new
00416     G4Fragment((G4int) typeDaughter->GetBaryonNumber(),
00417     (G4int) typeDaughter->GetPDGCharge(), dynamicDaughter->Get4Momentum());
00418   if (verboseLevel >= 2)
00419   {
00420     G4cout <<"Dynamic properties of the prefragment:" <<G4endl;
00421     G4cout.precision(6);
00422     dynamicDaughter->DumpInfo();
00423     G4cout <<"Nuclear properties of the prefragment:" <<G4endl;
00424     G4cout <<theFragment <<G4endl;
00425   }
00426   G4ReactionProductVector *products =
00427     theExcitationHandler->BreakItUp(*theFragment);
00428   delete theFragment;
00429   theFragment = NULL;
00430   
00431   G4ReactionProductVector::iterator iter;
00432   for (iter = products->begin(); iter != products->end(); ++iter)
00433   {
00434     G4DynamicParticle *secondary =
00435       new G4DynamicParticle((*iter)->GetDefinition(),
00436       (*iter)->GetTotalEnergy(), (*iter)->GetMomentum());
00437     theParticleChange.AddSecondary (secondary);
00438   }
00439 
00440   if (verboseLevel >= 2)
00441     G4cout <<"########################################"
00442            <<"########################################"
00443            <<G4endl;
00444  
00445   return &theParticleChange;
00446 }
00448 //
00449 void G4EMDissociation::PrintWelcomeMessage ()
00450 {
00451   G4cout <<G4endl;
00452   G4cout <<" ****************************************************************"
00453          <<G4endl;
00454   G4cout <<" EM dissociation model for nuclear-nuclear interactions activated"
00455          <<G4endl;
00456   G4cout <<" (Written by QinetiQ Ltd for the European Space Agency)"
00457          <<G4endl;
00458   G4cout <<" ****************************************************************"
00459          <<G4endl;
00460   G4cout << G4endl;
00461 
00462   return;
00463 }
00465 //

Generated on Mon May 27 17:48:08 2013 for Geant4 by  doxygen 1.4.7