G4AntiProtonAnnihilationAtRest.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 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 //   G4AntiProtonAnnihilationAtRest physics process
00027 //   Larry Felawka (TRIUMF), April 1998
00028 //---------------------------------------------------------------------
00029 
00030 #include <string.h>
00031 #include <cmath>
00032 #include <stdio.h>
00033 
00034 #include "G4AntiProtonAnnihilationAtRest.hh"
00035 #include "G4SystemOfUnits.hh"
00036 #include "G4DynamicParticle.hh"
00037 #include "G4ParticleTypes.hh"
00038 #include "Randomize.hh" 
00039 #include "G4HadronicProcessStore.hh"
00040 #include "G4HadronicDeprecate.hh"
00041  
00042 #define MAX_SECONDARIES 100
00043 
00044 // constructor
00045  
00046 G4AntiProtonAnnihilationAtRest::G4AntiProtonAnnihilationAtRest(const G4String& processName,
00047                                       G4ProcessType   aType ) :
00048   G4VRestProcess (processName, aType),       // initialization
00049   massPionMinus(G4PionMinus::PionMinus()->GetPDGMass()/GeV),
00050   massProton(G4Proton::Proton()->GetPDGMass()/GeV),
00051   massPionZero(G4PionZero::PionZero()->GetPDGMass()/GeV),
00052   massAntiProton(G4AntiProton::AntiProton()->GetPDGMass()/GeV),
00053   massPionPlus(G4PionPlus::PionPlus()->GetPDGMass()/GeV),
00054   massGamma(G4Gamma::Gamma()->GetPDGMass()/GeV),
00055   pdefGamma(G4Gamma::Gamma()),
00056   pdefPionPlus(G4PionPlus::PionPlus()),
00057   pdefPionZero(G4PionZero::PionZero()),
00058   pdefPionMinus(G4PionMinus::PionMinus()),
00059   pdefProton(G4Proton::Proton()),
00060   pdefAntiProton(G4AntiProton::AntiProton()),
00061   pdefNeutron(G4Neutron::Neutron()),
00062   pdefDeuteron(G4Deuteron::Deuteron()),
00063   pdefTriton(G4Triton::Triton()),
00064   pdefAlpha(G4Alpha::Alpha())
00065 {
00066   G4HadronicDeprecate("G4AntiProtonAnnihilationAtRest");
00067   if (verboseLevel>0) {
00068     G4cout << GetProcessName() << " is created "<< G4endl;
00069   }
00070   SetProcessSubType(fHadronAtRest);
00071   pv   = new G4GHEKinematicsVector [MAX_SECONDARIES+1];
00072   eve  = new G4GHEKinematicsVector [MAX_SECONDARIES];
00073   gkin = new G4GHEKinematicsVector [MAX_SECONDARIES];
00074 
00075   G4HadronicProcessStore::Instance()->RegisterExtraProcess(this);
00076 }
00077  
00078 // destructor
00079  
00080 G4AntiProtonAnnihilationAtRest::~G4AntiProtonAnnihilationAtRest()
00081 {
00082   G4HadronicProcessStore::Instance()->DeRegisterExtraProcess(this);
00083   delete [] pv;
00084   delete [] eve;
00085   delete [] gkin;
00086 }
00087  
00088 void G4AntiProtonAnnihilationAtRest::PreparePhysicsTable(const G4ParticleDefinition& p) 
00089 {
00090   G4HadronicProcessStore::Instance()->RegisterParticleForExtraProcess(this, &p);
00091 }
00092 
00093 void G4AntiProtonAnnihilationAtRest::BuildPhysicsTable(const G4ParticleDefinition& p) 
00094 {
00095   G4HadronicProcessStore::Instance()->PrintInfo(&p);
00096 }
00097  
00098 // methods.............................................................................
00099  
00100 G4bool G4AntiProtonAnnihilationAtRest::IsApplicable(
00101                                  const G4ParticleDefinition& particle
00102                                  )
00103 {
00104    return ( &particle == pdefAntiProton );
00105 
00106 }
00107  
00108 // Warning - this method may be optimized away if made "inline"
00109 G4int G4AntiProtonAnnihilationAtRest::GetNumberOfSecondaries()
00110 {
00111   return ( ngkine );
00112 
00113 }
00114 
00115 // Warning - this method may be optimized away if made "inline"
00116 G4GHEKinematicsVector* G4AntiProtonAnnihilationAtRest::GetSecondaryKinematics()
00117 {
00118   return ( &gkin[0] );
00119 
00120 }
00121 
00122 G4double G4AntiProtonAnnihilationAtRest::AtRestGetPhysicalInteractionLength(
00123                                    const G4Track& track,
00124                                    G4ForceCondition* condition
00125                                    )
00126 {
00127   // beggining of tracking 
00128   ResetNumberOfInteractionLengthLeft();
00129 
00130   // condition is set to "Not Forced"
00131   *condition = NotForced;
00132 
00133   // get mean life time
00134   currentInteractionLength = GetMeanLifeTime(track, condition);
00135 
00136   if ((currentInteractionLength <0.0) || (verboseLevel>2)){
00137     G4cout << "G4AntiProtonAnnihilationAtRestProcess::AtRestGetPhysicalInteractionLength ";
00138     G4cout << "[ " << GetProcessName() << "]" <<G4endl;
00139     track.GetDynamicParticle()->DumpInfo();
00140     G4cout << " in Material  " << track.GetMaterial()->GetName() <<G4endl;
00141     G4cout << "MeanLifeTime = " << currentInteractionLength/ns << "[ns]" <<G4endl;
00142   }
00143 
00144   return theNumberOfInteractionLengthLeft * currentInteractionLength;
00145 
00146 }
00147 
00148 G4VParticleChange* G4AntiProtonAnnihilationAtRest::AtRestDoIt(
00149                                             const G4Track& track,
00150                                             const G4Step& 
00151                                             )
00152 //
00153 // Handles AntiProtons at rest; a AntiProton can either create secondaries or
00154 // do nothing (in which case it should be sent back to decay-handling
00155 // section
00156 //
00157 {
00158 
00159 //   Initialize ParticleChange
00160 //     all members of G4VParticleChange are set to equal to 
00161 //     corresponding member in G4Track
00162 
00163   aParticleChange.Initialize(track);
00164 
00165 //   Store some global quantities that depend on current material and particle
00166 
00167   globalTime = track.GetGlobalTime()/s;
00168   G4Material * aMaterial = track.GetMaterial();
00169   const G4int numberOfElements = aMaterial->GetNumberOfElements();
00170   const G4ElementVector* theElementVector = aMaterial->GetElementVector();
00171 
00172   const G4double* theAtomicNumberDensity = aMaterial->GetAtomicNumDensityVector();
00173   G4double normalization = 0;
00174   for ( G4int i1=0; i1 < numberOfElements; i1++ )
00175   {
00176     normalization += theAtomicNumberDensity[i1] ; // change when nucleon specific
00177                                                   // probabilities are included.
00178   }
00179   G4double runningSum= 0.;
00180   G4double random = G4UniformRand()*normalization;
00181   for ( G4int i2=0; i2 < numberOfElements; i2++ )
00182   {
00183     runningSum += theAtomicNumberDensity[i2]; // change when nucleon specific
00184                                               // probabilities are included.
00185     if (random<=runningSum)
00186     {
00187       targetCharge = G4double((*theElementVector)[i2]->GetZ());
00188       targetAtomicMass = (*theElementVector)[i2]->GetN();
00189     }
00190   }
00191   if (random>runningSum)
00192   {
00193     targetCharge = G4double((*theElementVector)[numberOfElements-1]->GetZ());
00194     targetAtomicMass = (*theElementVector)[numberOfElements-1]->GetN();
00195 
00196   }
00197 
00198   if (verboseLevel>1) {
00199     G4cout << "G4AntiProtonAnnihilationAtRest::AtRestDoIt is invoked " <<G4endl;
00200     }
00201 
00202   G4ParticleMomentum momentum;
00203   G4float localtime;
00204 
00205   G4ThreeVector   position = track.GetPosition();
00206 
00207   GenerateSecondaries(); // Generate secondaries
00208 
00209   aParticleChange.SetNumberOfSecondaries( ngkine ); 
00210 
00211   for ( G4int isec = 0; isec < ngkine; isec++ ) {
00212     G4DynamicParticle* aNewParticle = new G4DynamicParticle;
00213     aNewParticle->SetDefinition( gkin[isec].GetParticleDef() );
00214     aNewParticle->SetMomentum( gkin[isec].GetMomentum() * GeV );
00215 
00216     localtime = globalTime + gkin[isec].GetTOF();
00217 
00218     G4Track* aNewTrack = new G4Track( aNewParticle, localtime*s, position );
00219                 aNewTrack->SetTouchableHandle(track.GetTouchableHandle());
00220     aParticleChange.AddSecondary( aNewTrack );
00221 
00222   }
00223 
00224   aParticleChange.ProposeLocalEnergyDeposit( 0.0*GeV );
00225 
00226   aParticleChange.ProposeTrackStatus(fStopAndKill); // Kill the incident AntiProton
00227 
00228 //   clear InteractionLengthLeft
00229 
00230   ResetNumberOfInteractionLengthLeft();
00231 
00232   return &aParticleChange;
00233 
00234 }
00235 
00236 
00237 void G4AntiProtonAnnihilationAtRest::GenerateSecondaries()
00238 {
00239   static G4int index;
00240   static G4int l;
00241   static G4int nopt;
00242   static G4int i;
00243   // DHW 15 May 2011: unused: static G4ParticleDefinition* jnd;
00244 
00245   for (i = 1; i <= MAX_SECONDARIES; ++i) {
00246     pv[i].SetZero();
00247   }
00248 
00249   ngkine = 0;            // number of generated secondary particles
00250   ntot = 0;
00251   result.SetZero();
00252   result.SetMass( massAntiProton );
00253   result.SetKineticEnergyAndUpdate( 0. );
00254   result.SetTOF( 0. );
00255   result.SetParticleDef( pdefAntiProton );
00256 
00257   AntiProtonAnnihilation(&nopt);
00258 
00259   // *** CHECK WHETHER THERE ARE NEW PARTICLES GENERATED ***
00260   if (ntot != 0 || result.GetParticleDef() != pdefAntiProton) {
00261     // *** CURRENT PARTICLE IS NOT THE SAME AS IN THE BEGINNING OR/AND ***
00262     // *** ONE OR MORE SECONDARIES HAVE BEEN GENERATED ***
00263 
00264     // --- INITIAL PARTICLE TYPE HAS BEEN CHANGED ==> PUT NEW TYPE ON ---
00265     // --- THE GEANT TEMPORARY STACK ---
00266 
00267     // --- PUT PARTICLE ON THE STACK ---
00268     gkin[0] = result;
00269     gkin[0].SetTOF( result.GetTOF() * 5e-11 );
00270     ngkine = 1;
00271 
00272     // --- ALL QUANTITIES ARE TAKEN FROM THE GHEISHA STACK WHERE THE ---
00273     // --- CONVENTION IS THE FOLLOWING ---
00274 
00275     // --- ONE OR MORE SECONDARIES HAVE BEEN GENERATED ---
00276     for (l = 1; l <= ntot; ++l) {
00277       index = l - 1;
00278       // DHW 15 May 2011: unused: jnd = eve[index].GetParticleDef();
00279 
00280       // --- ADD PARTICLE TO THE STACK IF STACK NOT YET FULL ---
00281       if (ngkine < MAX_SECONDARIES) {
00282         gkin[ngkine] = eve[index];
00283         gkin[ngkine].SetTOF( eve[index].GetTOF() * 5e-11 );
00284         ++ngkine;
00285       }
00286     }
00287   }
00288   else {
00289     // --- NO SECONDARIES GENERATED AND PARTICLE IS STILL THE SAME ---
00290     // --- ==> COPY EVERYTHING BACK IN THE CURRENT GEANT STACK ---
00291     ngkine = 0;
00292     ntot = 0;
00293     globalTime += result.GetTOF() * G4float(5e-11);
00294   }
00295 
00296   // --- LIMIT THE VALUE OF NGKINE IN CASE OF OVERFLOW ---
00297   ngkine = G4int(std::min(ngkine,G4int(MAX_SECONDARIES)));
00298 
00299 } // GenerateSecondaries
00300 
00301 
00302 void G4AntiProtonAnnihilationAtRest::Poisso(G4float xav, G4int *iran)
00303 {
00304   static G4int i;
00305   static G4float r, p1, p2, p3;
00306   static G4int fivex;
00307   static G4float rr, ran, rrr, ran1;
00308 
00309   // *** GENERATION OF POISSON DISTRIBUTION ***
00310   // *** NVE 16-MAR-1988 CERN GENEVA ***
00311   // ORIGIN : H.FESEFELDT (27-OCT-1983)
00312 
00313   // --- USE NORMAL DISTRIBUTION FOR <X> > 9.9 ---
00314   if (xav > G4float(9.9)) {
00315     // ** NORMAL DISTRIBUTION WITH SIGMA**2 = <X>
00316     Normal(&ran1);
00317     ran1 = xav + ran1 * std::sqrt(xav);
00318     *iran = G4int(ran1);
00319     if (*iran < 0) {
00320       *iran = 0;
00321     }
00322   }
00323   else {
00324     fivex = G4int(xav * G4float(5.));
00325     *iran = 0;
00326     if (fivex > 0) {
00327       r = std::exp(-G4double(xav));
00328       ran1 = G4UniformRand();
00329       if (ran1 > r) {
00330         rr = r;
00331         for (i = 1; i <= fivex; ++i) {
00332           ++(*iran);
00333           if (i <= 5) {
00334             rrr = std::pow(xav, G4float(i)) / NFac(i);
00335           }
00336           // ** STIRLING' S FORMULA FOR LARGE NUMBERS
00337           if (i > 5) {
00338             rrr = std::exp(i * std::log(xav) -
00339                       (i + G4float(.5)) * std::log(i * G4float(1.)) +
00340                       i - G4float(.9189385));
00341           }
00342           rr += r * rrr;
00343           if (ran1 <= rr) {
00344             break;
00345           }
00346         }
00347       }
00348     }
00349     else {
00350       // ** FOR VERY SMALL XAV TRY IRAN=1,2,3
00351       p1 = xav * std::exp(-G4double(xav));
00352       p2 = xav * p1 / G4float(2.);
00353       p3 = xav * p2 / G4float(3.);
00354       ran = G4UniformRand();
00355       if (ran >= p3) {
00356         if (ran >= p2) {
00357           if (ran >= p1) {
00358             *iran = 0;
00359           }
00360           else {
00361             *iran = 1;
00362           }
00363         }
00364         else {
00365           *iran = 2;
00366         }
00367       }
00368       else {
00369         *iran = 3;
00370       }
00371     }
00372   }
00373 
00374 } // Poisso
00375 
00376 
00377 G4int G4AntiProtonAnnihilationAtRest::NFac(G4int n)
00378 {
00379   G4int ret_val;
00380 
00381   static G4int i, j;
00382 
00383   // *** NVE 16-MAR-1988 CERN GENEVA ***
00384   // ORIGIN : H.FESEFELDT (27-OCT-1983)
00385 
00386   ret_val = 1;
00387   j = n;
00388   if (j > 1) {
00389     if (j > 10) {
00390       j = 10;
00391     }
00392     for (i = 2; i <= j; ++i) {
00393       ret_val *= i;
00394     }
00395   }
00396   return ret_val;
00397 
00398 } // NFac
00399 
00400 
00401 void G4AntiProtonAnnihilationAtRest::Normal(G4float *ran)
00402 {
00403   static G4int i;
00404 
00405   // *** NVE 14-APR-1988 CERN GENEVA ***
00406   // ORIGIN : H.FESEFELDT (27-OCT-1983)
00407 
00408   *ran = G4float(-6.);
00409   for (i = 1; i <= 12; ++i) {
00410     *ran += G4UniformRand();
00411   }
00412 
00413 } // Normal
00414 
00415 
00416 void G4AntiProtonAnnihilationAtRest::AntiProtonAnnihilation(G4int *nopt)
00417 {
00418   static G4float brr[3] = { G4float(.125),G4float(.25),G4float(.5) };
00419 
00420   G4float r__1;
00421 
00422   static G4int i, ii, kk;
00423   static G4int nt;
00424   static G4float cfa, eka;
00425   static G4int ika, nbl;
00426   static G4float ran, pcm;
00427   static G4int isw;
00428   static G4float tex;
00429   static G4ParticleDefinition* ipa1;
00430   static G4float ran1, ran2, ekin, tkin;
00431   static G4float targ;
00432   static G4ParticleDefinition* inve;
00433   static G4float ekin1, ekin2, black;
00434   static G4float pnrat, rmnve1, rmnve2;
00435   static G4float ek, en;
00436 
00437   // *** ANTI PROTON ANNIHILATION AT REST ***
00438   // *** NVE 04-MAR-1988 CERN GENEVA ***
00439   // ORIGIN : H.FESEFELDT (09-JULY-1987)
00440 
00441   // NOPT=0    NO ANNIHILATION
00442   // NOPT=1    ANNIH.IN PI+ PI-
00443   // NOPT=2    ANNIH.IN PI0 PI0
00444   // NOPT=3    ANNIH.IN PI- PI0
00445   // NOPT=4    ANNIH.IN GAMMA GAMMA
00446 
00447   pv[1].SetZero();
00448   pv[1].SetMass( massAntiProton );
00449   pv[1].SetKineticEnergyAndUpdate( 0. );
00450   pv[1].SetTOF( result.GetTOF() );
00451   pv[1].SetParticleDef( result.GetParticleDef() );
00452   isw = 1;
00453   ran = G4UniformRand();
00454   if (ran > brr[0]) {
00455     isw = 2;
00456   }
00457   if (ran > brr[1]) {
00458     isw = 3;
00459   }
00460   if (ran > brr[2]) {
00461     isw = 4;
00462   }
00463   *nopt = isw;
00464   // **
00465   // **  EVAPORATION
00466   // **
00467   if (isw == 1) {
00468     rmnve1 = massPionPlus;
00469     rmnve2 = massPionMinus;
00470   }
00471   else if (isw == 2) {
00472     rmnve1 = massPionZero;
00473     rmnve2 = massPionZero;
00474   }
00475   else if (isw == 3) {
00476     rmnve1 = massPionMinus;
00477     rmnve2 = massPionZero;
00478   }
00479   else if (isw == 4) {
00480     rmnve1 = massGamma;
00481     rmnve2 = massGamma;
00482   }
00483   ek = massProton + massAntiProton - rmnve1 - rmnve2;
00484   tkin = ExNu(ek);
00485   ek -= tkin;
00486   if (ek < G4float(1e-4)) {
00487     ek = G4float(1e-4);
00488   }
00489   ek *= G4float(.5);
00490   en = ek + (rmnve1 + rmnve2) * G4float(.5);
00491   r__1 = en * en - rmnve1 * rmnve2;
00492   pcm = r__1 > 0 ? std::sqrt(r__1) : 0;
00493   pv[2].SetZero();
00494   pv[2].SetMass( rmnve1 );
00495   pv[3].SetZero();
00496   pv[3].SetMass( rmnve2 );
00497   if (isw > 3) {
00498     pv[2].SetMass( 0. );
00499     pv[3].SetMass( 0. );
00500   }
00501   pv[2].SetEnergyAndUpdate( std::sqrt(pv[2].GetMass()*pv[2].GetMass()+pcm*pcm) );
00502   pv[2].SetTOF( result.GetTOF() );
00503   pv[3].SetEnergy( std::sqrt(pv[3].GetMass()*pv[3].GetMass()+pcm*pcm) );
00504   pv[3].SetMomentumAndUpdate( -pv[2].GetMomentum().x(), -pv[2].GetMomentum().y(), -pv[2].GetMomentum().z() );
00505   pv[3].SetTOF( result.GetTOF() );
00506   switch ((int)isw) {
00507     case 1:
00508       pv[2].SetParticleDef( pdefPionPlus );
00509       pv[3].SetParticleDef( pdefPionMinus );
00510       break;
00511     case 2:
00512       pv[2].SetParticleDef( pdefPionZero );
00513       pv[3].SetParticleDef( pdefPionZero );
00514       break;
00515     case 3:
00516       pv[2].SetParticleDef( pdefPionMinus );
00517       pv[3].SetParticleDef( pdefPionZero );
00518       break;
00519     case 4:
00520       pv[2].SetParticleDef( pdefGamma );
00521       pv[3].SetParticleDef( pdefGamma );
00522       break;
00523     default:
00524       break;
00525   }
00526   nt = 3;
00527   if (targetAtomicMass >= G4float(1.5)) {
00528     cfa = (targetAtomicMass - G4float(1.)) /
00529       G4float(120.) * G4float(.025) *
00530       std::exp(-G4double(targetAtomicMass - G4float(1.)) / G4float(120.));
00531     targ = G4float(1.);
00532     tex = evapEnergy1;
00533     if (tex >= G4float(.001)) {
00534       black = (targ * G4float(1.25) +
00535                G4float(1.5)) * evapEnergy1 / (evapEnergy1 + evapEnergy3);
00536       Poisso(black, &nbl);
00537       if (G4float(G4int(targ) + nbl) > targetAtomicMass) {
00538         nbl = G4int(targetAtomicMass - targ);
00539       }
00540       if (nt + nbl > (MAX_SECONDARIES - 2)) {
00541         nbl = (MAX_SECONDARIES - 2) - nt;
00542       }
00543       if (nbl > 0) {
00544         ekin = tex / nbl;
00545         ekin2 = G4float(0.);
00546         for (i = 1; i <= nbl; ++i) {
00547           if (nt == (MAX_SECONDARIES - 2)) {
00548             continue;
00549           }
00550           if (ekin2 > tex) {
00551             break;
00552           }
00553           ran1 = G4UniformRand();
00554           Normal(&ran2);
00555           ekin1 = -G4double(ekin) * std::log(ran1) -
00556             cfa * (ran2 * G4float(.5) + G4float(1.));
00557           if (ekin1 < G4float(0.)) {
00558             ekin1 = std::log(ran1) * G4float(-.01);
00559           }
00560           ekin1 *= G4float(1.);
00561           ekin2 += ekin1;
00562           if (ekin2 > tex) {
00563             ekin1 = tex - (ekin2 - ekin1);
00564           }
00565           if (ekin1 < G4float(0.)) {
00566             ekin1 = G4float(.001);
00567           }
00568           ipa1 = pdefNeutron;
00569           pnrat = G4float(1.) - targetCharge / targetAtomicMass;
00570           if (G4UniformRand() > pnrat) {
00571             ipa1 = pdefProton;
00572           }
00573           ++nt;
00574           pv[nt].SetZero();
00575           pv[nt].SetMass( ipa1->GetPDGMass()/GeV );
00576           pv[nt].SetKineticEnergyAndUpdate( ekin1 );
00577           pv[nt].SetTOF( result.GetTOF() );
00578           pv[nt].SetParticleDef( ipa1 );
00579         }
00580         if (targetAtomicMass >= G4float(230.) && ek <= G4float(2.)) {
00581           ii = nt + 1;
00582           kk = 0;
00583           eka = ek;
00584           if (eka > G4float(1.)) {
00585             eka *= eka;
00586           }
00587           if (eka < G4float(.1)) {
00588             eka = G4float(.1);
00589           }
00590           ika = G4int(G4float(3.6) / eka);
00591           for (i = 1; i <= nt; ++i) {
00592             --ii;
00593             if (pv[ii].GetParticleDef() != pdefProton) {
00594               continue;
00595             }
00596             ipa1 = pdefNeutron;
00597             pv[ii].SetMass( ipa1->GetPDGMass()/GeV );
00598             pv[ii].SetParticleDef( ipa1 );
00599             ++kk;
00600             if (kk > ika) {
00601               break;
00602             }
00603           }
00604         }
00605       }
00606     }
00607     // **
00608     // ** THEN ALSO DEUTERONS, TRITONS AND ALPHAS
00609     // **
00610     tex = evapEnergy3;
00611     if (tex >= G4float(.001)) {
00612       black = (targ * G4float(1.25) + G4float(1.5)) * evapEnergy3 /
00613         (evapEnergy1 + evapEnergy3);
00614       Poisso(black, &nbl);
00615       if (nt + nbl > (MAX_SECONDARIES - 2)) {
00616         nbl = (MAX_SECONDARIES - 2) - nt;
00617       }
00618       if (nbl > 0) {
00619         ekin = tex / nbl;
00620         ekin2 = G4float(0.);
00621         for (i = 1; i <= nbl; ++i) {
00622           if (nt == (MAX_SECONDARIES - 2)) {
00623             continue;
00624           }
00625           if (ekin2 > tex) {
00626             break;
00627           }
00628           ran1 = G4UniformRand();
00629           Normal(&ran2);
00630           ekin1 = -G4double(ekin) * std::log(ran1) -
00631             cfa * (ran2 * G4float(.5) + G4float(1.));
00632           if (ekin1 < G4float(0.)) {
00633             ekin1 = std::log(ran1) * G4float(-.01);
00634           }
00635           ekin1 *= G4float(1.);
00636           ekin2 += ekin1;
00637           if (ekin2 > tex) {
00638             ekin1 = tex - (ekin2 - ekin1);
00639           }
00640           if (ekin1 < G4float(0.)) {
00641             ekin1 = G4float(.001);
00642           }
00643           ran = G4UniformRand();
00644           inve = pdefDeuteron;
00645           if (ran > G4float(.6)) {
00646             inve = pdefTriton;
00647           }
00648           if (ran > G4float(.9)) {
00649             inve = pdefAlpha;
00650           }
00651           ++nt;
00652           pv[nt].SetZero();
00653           pv[nt].SetMass( inve->GetPDGMass()/GeV );
00654           pv[nt].SetKineticEnergyAndUpdate( ekin1 );
00655           pv[nt].SetTOF( result.GetTOF() );
00656           pv[nt].SetParticleDef( inve );
00657         }
00658       }
00659     }
00660   }
00661   result = pv[2];
00662   if (nt == 2) {
00663     return;
00664   }
00665   for (i = 3; i <= nt; ++i) {
00666     if (ntot >= MAX_SECONDARIES) {
00667       return;
00668     }
00669     eve[ntot++] = pv[i];
00670   }
00671 
00672 } // AntiProtonAnnihilation
00673 
00674 
00675 G4double G4AntiProtonAnnihilationAtRest::ExNu(G4float ek1)
00676 {
00677   G4float ret_val, r__1;
00678 
00679   static G4float cfa, gfa, ran1, ran2, ekin1, atno3;
00680   static G4int magic;
00681   static G4float fpdiv;
00682 
00683   // *** NUCLEAR EVAPORATION AS FUNCTION OF ATOMIC NUMBER ATNO ***
00684   // *** AND KINETIC ENERGY EKIN OF PRIMARY PARTICLE ***
00685   // *** NVE 04-MAR-1988 CERN GENEVA ***
00686   // ORIGIN : H.FESEFELDT (10-DEC-1986)
00687 
00688   ret_val = G4float(0.);
00689   if (targetAtomicMass >= G4float(1.5)) {
00690     magic = 0;
00691     if (G4int(targetCharge + G4float(.1)) == 82) {
00692       magic = 1;
00693     }
00694     ekin1 = ek1;
00695     if (ekin1 < G4float(.1)) {
00696       ekin1 = G4float(.1);
00697     }
00698     if (ekin1 > G4float(4.)) {
00699       ekin1 = G4float(4.);
00700     }
00701     // **   0.35 VALUE AT 1 GEV
00702     // **   0.05 VALUE AT 0.1 GEV
00703     cfa = G4float(.13043478260869565);
00704     cfa = cfa * std::log(ekin1) + G4float(.35);
00705     if (cfa < G4float(.15)) {
00706       cfa = G4float(.15);
00707     }
00708     ret_val = cfa * G4float(7.716) * std::exp(-G4double(cfa));
00709     atno3 = targetAtomicMass;
00710     if (atno3 > G4float(120.)) {
00711       atno3 = G4float(120.);
00712     }
00713     cfa = (atno3 - G4float(1.)) /
00714       G4float(120.) * std::exp(-G4double(atno3 - G4float(1.)) / G4float(120.));
00715     ret_val *= cfa;
00716     r__1 = ekin1;
00717     fpdiv = G4float(1.) - r__1 * r__1 * G4float(.25);
00718     if (fpdiv < G4float(.5)) {
00719       fpdiv = G4float(.5);
00720     }
00721     gfa = (targetAtomicMass - G4float(1.)) /
00722       G4float(70.) * G4float(2.) *
00723       std::exp(-G4double(targetAtomicMass - G4float(1.)) / G4float(70.));
00724     evapEnergy1 = ret_val * fpdiv;
00725     evapEnergy3 = ret_val - evapEnergy1;
00726     Normal(&ran1);
00727     Normal(&ran2);
00728     if (magic == 1) {
00729       ran1 = G4float(0.);
00730       ran2 = G4float(0.);
00731     }
00732     evapEnergy1 *= ran1 * gfa + G4float(1.);
00733     if (evapEnergy1 < G4float(0.)) {
00734       evapEnergy1 = G4float(0.);
00735     }
00736     evapEnergy3 *= ran2 * gfa + G4float(1.);
00737     if (evapEnergy3 < G4float(0.)) {
00738       evapEnergy3 = G4float(0.);
00739     }
00740     while ((ret_val = evapEnergy1 + evapEnergy3) >= ek1) {
00741       evapEnergy1 *= G4float(1.) - G4UniformRand() * G4float(.5);
00742       evapEnergy3 *= G4float(1.) - G4UniformRand() * G4float(.5);
00743     }
00744   }
00745   return ret_val;
00746 
00747 } // ExNu

Generated on Mon May 27 17:47:40 2013 for Geant4 by  doxygen 1.4.7