G4PionMinusAbsorptionAtRest.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 //   G4PionMinusAbsorptionAtRest physics process
00027 //   Larry Felawka (TRIUMF), April 1998
00028 //---------------------------------------------------------------------
00029 
00030 #include <string.h>
00031 #include <cmath>
00032 #include <stdio.h>
00033 
00034 #include "G4PionMinusAbsorptionAtRest.hh"
00035 #include "G4DynamicParticle.hh"
00036 #include "G4ParticleTypes.hh"
00037 #include "G4SystemOfUnits.hh"
00038 #include "Randomize.hh" 
00039 #include "G4HadronicProcessStore.hh"
00040  #include "G4HadronicDeprecate.hh"
00041 
00042 #define MAX_SECONDARIES 100
00043 
00044 // constructor
00045  
00046 G4PionMinusAbsorptionAtRest::G4PionMinusAbsorptionAtRest(const G4String& processName,
00047                                       G4ProcessType   aType ) :
00048   G4VRestProcess (processName, aType),       // initialization
00049   massPionMinus(G4PionMinus::PionMinus()->GetPDGMass()/GeV),
00050   pdefGamma(G4Gamma::Gamma()),
00051   pdefPionZero(G4PionZero::PionZero()),
00052   pdefPionMinus(G4PionMinus::PionMinus()),
00053   pdefProton(G4Proton::Proton()),
00054   pdefNeutron(G4Neutron::Neutron()),
00055   pdefDeuteron(G4Deuteron::Deuteron()),
00056   pdefTriton(G4Triton::Triton()),
00057   pdefAlpha(G4Alpha::Alpha())
00058 {
00059   G4HadronicDeprecate("G4PiMinusAbsorptionAtRest");
00060 
00061   if (verboseLevel>0) {
00062     G4cout << GetProcessName() << " is created "<< G4endl;
00063   }
00064   SetProcessSubType(fHadronAtRest);
00065   pv   = new G4GHEKinematicsVector [MAX_SECONDARIES+1];
00066   eve  = new G4GHEKinematicsVector [MAX_SECONDARIES];
00067   gkin = new G4GHEKinematicsVector [MAX_SECONDARIES];
00068 
00069   G4HadronicProcessStore::Instance()->RegisterExtraProcess(this);
00070 }
00071  
00072 // destructor
00073  
00074 G4PionMinusAbsorptionAtRest::~G4PionMinusAbsorptionAtRest()
00075 {
00076   G4HadronicProcessStore::Instance()->DeRegisterExtraProcess(this);
00077   delete [] pv;
00078   delete [] eve;
00079   delete [] gkin;
00080 }
00081  
00082 void G4PionMinusAbsorptionAtRest::PreparePhysicsTable(const G4ParticleDefinition& p) 
00083 {
00084   G4HadronicProcessStore::Instance()->RegisterParticleForExtraProcess(this, &p);
00085 }
00086 
00087 void G4PionMinusAbsorptionAtRest::BuildPhysicsTable(const G4ParticleDefinition& p) 
00088 {
00089   G4HadronicProcessStore::Instance()->PrintInfo(&p);
00090 }
00091  
00092 // methods.............................................................................
00093  
00094 G4bool G4PionMinusAbsorptionAtRest::IsApplicable(
00095                                  const G4ParticleDefinition& particle
00096                                  )
00097 {
00098    return ( &particle == pdefPionMinus );
00099 
00100 }
00101  
00102 // Warning - this method may be optimized away if made "inline"
00103 G4int G4PionMinusAbsorptionAtRest::GetNumberOfSecondaries()
00104 {
00105   return ( ngkine );
00106 
00107 }
00108 
00109 // Warning - this method may be optimized away if made "inline"
00110 G4GHEKinematicsVector* G4PionMinusAbsorptionAtRest::GetSecondaryKinematics()
00111 {
00112   return ( &gkin[0] );
00113 
00114 }
00115 
00116 G4double G4PionMinusAbsorptionAtRest::AtRestGetPhysicalInteractionLength(
00117                                    const G4Track& track,
00118                                    G4ForceCondition* condition
00119                                    )
00120 {
00121   // beggining of tracking 
00122   ResetNumberOfInteractionLengthLeft();
00123 
00124   // condition is set to "Not Forced"
00125   *condition = NotForced;
00126 
00127   // get mean life time
00128   currentInteractionLength = GetMeanLifeTime(track, condition);
00129 
00130   if ((currentInteractionLength <0.0) || (verboseLevel>2)){
00131     G4cout << "G4PionMinusAbsorptionAtRestProcess::AtRestGetPhysicalInteractionLength ";
00132     G4cout << "[ " << GetProcessName() << "]" <<G4endl;
00133     track.GetDynamicParticle()->DumpInfo();
00134     G4cout << " in Material  " << track.GetMaterial()->GetName() <<G4endl;
00135     G4cout << "MeanLifeTime = " << currentInteractionLength/ns << "[ns]" <<G4endl;
00136   }
00137 
00138   return theNumberOfInteractionLengthLeft * currentInteractionLength;
00139 
00140 }
00141 
00142 G4VParticleChange* G4PionMinusAbsorptionAtRest::AtRestDoIt(
00143                                             const G4Track& track,
00144                                             const G4Step& 
00145                                             )
00146 //
00147 // Handles PionMinuss at rest; a PionMinus can either create secondaries or
00148 // do nothing (in which case it should be sent back to decay-handling
00149 // section
00150 //
00151 {
00152 
00153 //   Initialize ParticleChange
00154 //     all members of G4VParticleChange are set to equal to 
00155 //     corresponding member in G4Track
00156 
00157   aParticleChange.Initialize(track);
00158 
00159 //   Store some global quantities that depend on current material and particle
00160 
00161   globalTime = track.GetGlobalTime()/s;
00162   G4Material * aMaterial = track.GetMaterial();
00163   const G4int numberOfElements = aMaterial->GetNumberOfElements();
00164   const G4ElementVector* theElementVector = aMaterial->GetElementVector();
00165 
00166   const G4double* theAtomicNumberDensity = aMaterial->GetAtomicNumDensityVector();
00167   G4double normalization = 0;
00168   for ( G4int i1=0; i1 < numberOfElements; i1++ )
00169   {
00170     normalization += theAtomicNumberDensity[i1] ; // change when nucleon specific
00171                                                   // probabilities are included.
00172   }
00173   G4double runningSum= 0.;
00174   G4double random = G4UniformRand()*normalization;
00175   for ( G4int i2=0; i2 < numberOfElements; i2++ )
00176   {
00177     runningSum += theAtomicNumberDensity[i2]; // change when nucleon specific
00178                                               // probabilities are included.
00179     if (random<=runningSum)
00180     {
00181       targetCharge = G4double((*theElementVector)[i2]->GetZ());
00182       targetAtomicMass = (*theElementVector)[i2]->GetN();
00183     }
00184   }
00185   if (random>runningSum)
00186   {
00187     targetCharge = G4double((*theElementVector)[numberOfElements-1]->GetZ());
00188     targetAtomicMass = (*theElementVector)[numberOfElements-1]->GetN();
00189 
00190   }
00191 
00192   if (verboseLevel>1) {
00193     G4cout << "G4PionMinusAbsorptionAtRest::AtRestDoIt is invoked " <<G4endl;
00194     }
00195 
00196   G4ParticleMomentum momentum;
00197   G4float localtime;
00198 
00199   G4ThreeVector   position = track.GetPosition();
00200 
00201   GenerateSecondaries(); // Generate secondaries
00202 
00203   aParticleChange.SetNumberOfSecondaries( ngkine ); 
00204 
00205   for ( G4int isec = 0; isec < ngkine; isec++ ) {
00206     G4DynamicParticle* aNewParticle = new G4DynamicParticle;
00207     aNewParticle->SetDefinition( gkin[isec].GetParticleDef() );
00208     aNewParticle->SetMomentum( gkin[isec].GetMomentum() * GeV );
00209 
00210     localtime = globalTime + gkin[isec].GetTOF();
00211 
00212     G4Track* aNewTrack = new G4Track( aNewParticle, localtime*s, position );
00213                 aNewTrack->SetTouchableHandle(track.GetTouchableHandle());
00214     aParticleChange.AddSecondary( aNewTrack );
00215 
00216   }
00217 
00218   aParticleChange.ProposeLocalEnergyDeposit( 0.0*GeV );
00219 
00220   aParticleChange.ProposeTrackStatus(fStopAndKill); // Kill the incident PionMinus
00221 
00222 //   clear InteractionLengthLeft
00223 
00224   ResetNumberOfInteractionLengthLeft();
00225 
00226   return &aParticleChange;
00227 
00228 }
00229 
00230 
00231 void G4PionMinusAbsorptionAtRest::GenerateSecondaries()
00232 {
00233   static G4int index;
00234   static G4int l;
00235   static G4int nopt;
00236   static G4int i;
00237   // DHW 15 May 2011: unused: static G4ParticleDefinition* jnd;
00238 
00239   for (i = 1; i <= MAX_SECONDARIES; ++i) {
00240     pv[i].SetZero();
00241   }
00242 
00243   ngkine = 0;            // number of generated secondary particles
00244   ntot = 0;
00245   result.SetZero();
00246   result.SetMass( massPionMinus );
00247   result.SetKineticEnergyAndUpdate( 0. );
00248   result.SetTOF( 0. );
00249   result.SetParticleDef( pdefPionMinus );
00250 
00251   PionMinusAbsorption(&nopt);
00252 
00253   // *** CHECK WHETHER THERE ARE NEW PARTICLES GENERATED ***
00254   if (ntot != 0 || result.GetParticleDef() != pdefPionMinus) {
00255     // *** CURRENT PARTICLE IS NOT THE SAME AS IN THE BEGINNING OR/AND ***
00256     // *** ONE OR MORE SECONDARIES HAVE BEEN GENERATED ***
00257 
00258     // --- INITIAL PARTICLE TYPE HAS BEEN CHANGED ==> PUT NEW TYPE ON ---
00259     // --- THE GEANT TEMPORARY STACK ---
00260 
00261     // --- PUT PARTICLE ON THE STACK ---
00262     gkin[0] = result;
00263     gkin[0].SetTOF( result.GetTOF() * 5e-11 );
00264     ngkine = 1;
00265 
00266     // --- ALL QUANTITIES ARE TAKEN FROM THE GHEISHA STACK WHERE THE ---
00267     // --- CONVENTION IS THE FOLLOWING ---
00268 
00269     // --- ONE OR MORE SECONDARIES HAVE BEEN GENERATED ---
00270     for (l = 1; l <= ntot; ++l) {
00271       index = l - 1;
00272       // DHW 15 May 2011: unused: jnd = eve[index].GetParticleDef();
00273 
00274       // --- ADD PARTICLE TO THE STACK IF STACK NOT YET FULL ---
00275       if (ngkine < MAX_SECONDARIES) {
00276         gkin[ngkine] = eve[index];
00277         gkin[ngkine].SetTOF( eve[index].GetTOF() * 5e-11 );
00278         ++ngkine;
00279       }
00280     }
00281   }
00282   else {
00283     // --- NO SECONDARIES GENERATED AND PARTICLE IS STILL THE SAME ---
00284     // --- ==> COPY EVERYTHING BACK IN THE CURRENT GEANT STACK ---
00285     ngkine = 0;
00286     ntot = 0;
00287     globalTime += result.GetTOF() * G4float(5e-11);
00288   }
00289 
00290   // --- LIMIT THE VALUE OF NGKINE IN CASE OF OVERFLOW ---
00291   ngkine = G4int(std::min(ngkine,G4int(MAX_SECONDARIES)));
00292 
00293 } // GenerateSecondaries
00294 
00295 
00296 void G4PionMinusAbsorptionAtRest::PionMinusAbsorption(G4int *nopt)
00297 {
00298   static G4int i;
00299   static G4int nt, nbl;
00300   static G4float ran, tex;
00301   static G4int isw;
00302   static G4float ran2, tof1, ekin;
00303   static G4float ekin1, ekin2, black;
00304   static G4float pnrat;
00305   static G4ParticleDefinition* ipa1;
00306   static G4ParticleDefinition* inve;
00307 
00308   // *** CHARGED PION ABSORPTION BY A NUCLEUS ***
00309   // *** NVE 04-MAR-1988 CERN GENEVA ***
00310 
00311   // ORIGIN : H.FESEFELDT (09-JULY-1987)
00312 
00313   // PANOFSKY RATIO (PI- P --> N PI0/PI- P --> N GAMMA) = 3/2
00314   // FOR CAPTURE ON PROTON (HYDROGEN),
00315   // STAR PRODUCTION FOR HEAVIER ELEMENTS
00316 
00317   pv[1].SetZero();
00318   pv[1].SetMass( massPionMinus );
00319   pv[1].SetKineticEnergyAndUpdate( 0. );
00320   pv[1].SetTOF( result.GetTOF() );
00321   pv[1].SetParticleDef( result.GetParticleDef() );
00322   if (targetAtomicMass <= G4float(1.5)) {
00323     ran = G4UniformRand();
00324     isw = 1;
00325     if (ran < G4float(.33)) {
00326       isw = 2;
00327     }
00328     *nopt = isw;
00329     ran = G4UniformRand();
00330     tof1 = std::log(ran) * G4float(-25.);
00331     tof1 *= G4float(20.);
00332     if (isw != 1) {
00333       pv[2].SetZero();
00334       pv[2].SetMass( 0. );
00335       pv[2].SetKineticEnergyAndUpdate( .02 );
00336       pv[2].SetTOF( result.GetTOF() + tof1 );
00337       pv[2].SetParticleDef( pdefGamma );
00338     }
00339     else {
00340       pv[2] = pv[1];
00341       pv[2].SetTOF( result.GetTOF() + tof1 );
00342       pv[2].SetParticleDef( pdefPionZero );
00343     }
00344     result = pv[2];
00345   }
00346   else {
00347     // **
00348     // ** STAR PRODUCTION FOR PION ABSORPTION IN HEAVY ELEMENTS
00349     // **
00350     evapEnergy1 = G4float(.0135);
00351     evapEnergy3 = G4float(.0058);
00352     nt = 1;
00353     tex = evapEnergy1;
00354     black = std::log(targetAtomicMass) * G4float(.5);
00355     Poisso(black, &nbl);
00356     if (nbl <= 0) {
00357       nbl = 1;
00358     }
00359     if (nt + nbl > (MAX_SECONDARIES - 2)) {
00360       nbl = (MAX_SECONDARIES - 2) - nt;
00361     }
00362     ekin = tex / nbl;
00363     ekin2 = G4float(0.);
00364     for (i = 1; i <= nbl; ++i) {
00365       if (nt == (MAX_SECONDARIES - 2)) {
00366         continue;
00367       }
00368       ran2 = G4UniformRand();
00369       ekin1 = -G4double(ekin) * std::log(ran2);
00370       ekin2 += ekin1;
00371       ipa1 = pdefNeutron;
00372       pnrat = G4float(1.) - targetCharge / targetAtomicMass;
00373       if (G4UniformRand() > pnrat) {
00374         ipa1 = pdefProton;
00375       }
00376       ++nt;
00377       pv[nt].SetZero();
00378       pv[nt].SetMass( ipa1->GetPDGMass()/GeV );
00379       pv[nt].SetKineticEnergyAndUpdate( ekin1 );
00380       pv[nt].SetTOF( 2. );
00381       pv[nt].SetParticleDef( ipa1 );
00382       if (ekin2 > tex) {
00383         break;
00384       }
00385     }
00386     tex = evapEnergy3;
00387     black = std::log(targetAtomicMass) * G4float(.5);
00388     Poisso(black, &nbl);
00389     if (nt + nbl > (MAX_SECONDARIES - 2)) {
00390       nbl = (MAX_SECONDARIES - 2) - nt;
00391     }
00392     if (nbl <= 0) {
00393       nbl = 1;
00394     }
00395     ekin = tex / nbl;
00396     ekin2 = G4float(0.);
00397     for (i = 1; i <= nbl; ++i) {
00398       if (nt == (MAX_SECONDARIES - 2)) {
00399         continue;
00400       }
00401       ran2 = G4UniformRand();
00402       ekin1 = -G4double(ekin) * std::log(ran2);
00403       ekin2 += ekin1;
00404       ++nt;
00405       ran = G4UniformRand();
00406       inve= pdefDeuteron;
00407       if (ran > G4float(.6)) {
00408         inve = pdefTriton;
00409       }
00410       if (ran > G4float(.9)) {
00411         inve = pdefAlpha;
00412       }
00413       pv[nt].SetZero();
00414       pv[nt].SetMass( inve->GetPDGMass()/GeV );
00415       pv[nt].SetKineticEnergyAndUpdate( ekin1 );
00416       pv[nt].SetTOF( 2. );
00417       pv[nt].SetParticleDef( inve );
00418       if (ekin2 > tex) {
00419         break;
00420       }
00421     }
00422     // **
00423     // ** STORE ON EVENT COMMON
00424     // **
00425     ran = G4UniformRand();
00426     tof1 = std::log(ran) * G4float(-25.);
00427     tof1 *= G4float(20.);
00428     for (i = 2; i <= nt; ++i) {
00429       pv[i].SetTOF( result.GetTOF() + tof1 );
00430     }
00431     result = pv[2];
00432     for (i = 3; i <= nt; ++i) {
00433       if (ntot >= MAX_SECONDARIES) {
00434         break;
00435       }
00436       eve[ntot++] = pv[i];
00437     }
00438   }
00439 
00440 } // PionMinusAbsorption
00441 
00442 
00443 void G4PionMinusAbsorptionAtRest::Poisso(G4float xav, G4int *iran)
00444 {
00445   static G4int i;
00446   static G4float r, p1, p2, p3;
00447   static G4int fivex;
00448   static G4float rr, ran, rrr, ran1;
00449 
00450   // *** GENERATION OF POISSON DISTRIBUTION ***
00451   // *** NVE 16-MAR-1988 CERN GENEVA ***
00452   // ORIGIN : H.FESEFELDT (27-OCT-1983)
00453 
00454   // --- USE NORMAL DISTRIBUTION FOR <X> > 9.9 ---
00455   if (xav > G4float(9.9)) {
00456     // ** NORMAL DISTRIBUTION WITH SIGMA**2 = <X>
00457     Normal(&ran1);
00458     ran1 = xav + ran1 * std::sqrt(xav);
00459     *iran = G4int(ran1);
00460     if (*iran < 0) {
00461       *iran = 0;
00462     }
00463   }
00464   else {
00465     fivex = G4int(xav * G4float(5.));
00466     *iran = 0;
00467     if (fivex > 0) {
00468       r = std::exp(-G4double(xav));
00469       ran1 = G4UniformRand();
00470       if (ran1 > r) {
00471         rr = r;
00472         for (i = 1; i <= fivex; ++i) {
00473           ++(*iran);
00474           if (i <= 5) {
00475             rrr = std::pow(xav, G4float(i)) / NFac(i);
00476           }
00477           // ** STIRLING' S FORMULA FOR LARGE NUMBERS
00478           if (i > 5) {
00479             rrr = std::exp(i * std::log(xav) -
00480                       (i + G4float(.5)) * std::log(i * G4float(1.)) +
00481                       i - G4float(.9189385));
00482           }
00483           rr += r * rrr;
00484           if (ran1 <= rr) {
00485             break;
00486           }
00487         }
00488       }
00489     }
00490     else {
00491       // ** FOR VERY SMALL XAV TRY IRAN=1,2,3
00492       p1 = xav * std::exp(-G4double(xav));
00493       p2 = xav * p1 / G4float(2.);
00494       p3 = xav * p2 / G4float(3.);
00495       ran = G4UniformRand();
00496       if (ran >= p3) {
00497         if (ran >= p2) {
00498           if (ran >= p1) {
00499             *iran = 0;
00500           }
00501           else {
00502             *iran = 1;
00503           }
00504         }
00505         else {
00506           *iran = 2;
00507         }
00508       }
00509       else {
00510         *iran = 3;
00511       }
00512     }
00513   }
00514 
00515 } // Poisso
00516 
00517 
00518 G4int G4PionMinusAbsorptionAtRest::NFac(G4int n)
00519 {
00520   G4int ret_val;
00521   static G4int i, j;
00522 
00523   // *** NVE 16-MAR-1988 CERN GENEVA ***
00524   // ORIGIN : H.FESEFELDT (27-OCT-1983)
00525 
00526   ret_val = 1;
00527   j = n;
00528   if (j > 1) {
00529     if (j > 10) {
00530       j = 10;
00531     }
00532     for (i = 2; i <= j; ++i) {
00533       ret_val *= i;
00534     }
00535   }
00536   return ret_val;
00537 
00538 } // NFac
00539 
00540 
00541 void G4PionMinusAbsorptionAtRest::Normal(G4float *ran)
00542 {
00543   static G4int i;
00544 
00545   // *** NVE 14-APR-1988 CERN GENEVA ***
00546   // ORIGIN : H.FESEFELDT (27-OCT-1983)
00547 
00548   *ran = G4float(-6.);
00549   for (i = 1; i <= 12; ++i) {
00550     *ran += G4UniformRand();
00551   }
00552 
00553 } // Normal

Generated on Mon May 27 17:49:21 2013 for Geant4 by  doxygen 1.4.7