G4NonEquilibriumEvaporator.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 // $Id$
00027 //
00028 // 20100114  M. Kelsey -- Remove G4CascadeMomentum, use G4LorentzVector directly
00029 // 20100309  M. Kelsey -- Use new generateWithRandomAngles for theta,phi stuff;
00030 //              eliminate some unnecessary std::pow()
00031 // 20100412  M. Kelsey -- Pass G4CollisionOutput by ref to ::collide()
00032 // 20100413  M. Kelsey -- Pass buffers to paraMaker[Truncated]
00033 // 20100517  M. Kelsey -- Inherit from common base class
00034 // 20100617  M. Kelsey -- Remove "RUN" preprocessor flag and all "#else" code
00035 // 20100622  M. Kelsey -- Use local "bindingEnergy()" function to call through.
00036 // 20100701  M. Kelsey -- Don't need to add excitation to nuclear mass; compute
00037 //              new excitation energies properly (mass differences)
00038 // 20100713  M. Kelsey -- Add conservation checking, diagnostic messages.
00039 // 20100714  M. Kelsey -- Move conservation checking to base class
00040 // 20100719  M. Kelsey -- Simplify EEXS calculations with particle evaporation.
00041 // 20100724  M. Kelsey -- Replace std::vector<> D with trivial D[3] array.
00042 // 20100914  M. Kelsey -- Migrate to integer A and Z: this involves replacing
00043 //              a number of G4double terms with G4int, with consequent casts.
00044 // 20110214  M. Kelsey -- Follow G4InuclParticle::Model enumerator migration
00045 // 20110922  M. Kelsey -- Follow G4InuclParticle::print(ostream&) migration
00046 // 20120608  M. Kelsey -- Fix variable-name "shadowing" compiler warnings.
00047 // 20121009  M. Kelsey -- Add some high-verbosity debugging output
00048 
00049 #include <cmath>
00050 
00051 #include "G4NonEquilibriumEvaporator.hh"
00052 #include "G4SystemOfUnits.hh"
00053 #include "G4CollisionOutput.hh"
00054 #include "G4InuclElementaryParticle.hh"
00055 #include "G4InuclNuclei.hh"
00056 #include "G4InuclSpecialFunctions.hh"
00057 #include "G4LorentzConvertor.hh"
00058 
00059 using namespace G4InuclSpecialFunctions;
00060 
00061 
00062 G4NonEquilibriumEvaporator::G4NonEquilibriumEvaporator()
00063   : G4CascadeColliderBase("G4NonEquilibriumEvaporator") {}
00064 
00065 
00066 void G4NonEquilibriumEvaporator::collide(G4InuclParticle* /*bullet*/,
00067                                          G4InuclParticle* target,
00068                                          G4CollisionOutput& output) {
00069 
00070   if (verboseLevel) {
00071     G4cout << " >>> G4NonEquilibriumEvaporator::collide" << G4endl;
00072   }
00073 
00074   // Sanity check
00075   G4InuclNuclei* nuclei_target = dynamic_cast<G4InuclNuclei*>(target);
00076   if (!nuclei_target) {
00077     G4cerr << " NonEquilibriumEvaporator -> target is not nuclei " << G4endl;    
00078     return;
00079   }
00080 
00081   if (verboseLevel > 2) G4cout << " evaporating target:\n" << *target << G4endl;
00082   
00083   const G4int a_cut = 5;
00084   const G4int z_cut = 3;
00085 
00086   const G4double eexs_cut = 0.1;
00087 
00088   const G4double coul_coeff = 1.4;
00089   const G4int itry_max = 1000;
00090   const G4double small_ekin = 1.0e-6;
00091   const G4double width_cut = 0.005;
00092 
00093   G4int A = nuclei_target->getA();
00094   G4int Z = nuclei_target->getZ();
00095   
00096   G4LorentzVector PEX = nuclei_target->getMomentum();
00097   G4LorentzVector pin = PEX;            // Save original four-vector for later
00098   
00099   G4double EEXS = nuclei_target->getExitationEnergy();
00100   
00101   G4ExitonConfiguration config = nuclei_target->getExitonConfiguration();  
00102 
00103   G4int QPP = config.protonQuasiParticles;
00104   G4int QNP = config.neutronQuasiParticles; 
00105   G4int QPH = config.protonHoles;
00106   G4int QNH = config.neutronHoles; 
00107   
00108   G4int QP = QPP + QNP;
00109   G4int QH = QPH + QNH;
00110   G4int QEX = QP + QH;
00111   
00112   G4InuclElementaryParticle dummy(small_ekin, 1);
00113   G4LorentzConvertor toTheExitonSystemRestFrame;
00114   //*** toTheExitonSystemRestFrame.setVerbose(verboseLevel);
00115   toTheExitonSystemRestFrame.setBullet(dummy);
00116   
00117   G4double EFN = FermiEnergy(A, Z, 0);
00118   G4double EFP = FermiEnergy(A, Z, 1);
00119   
00120   G4int AR = A - QP;
00121   G4int ZR = Z - QPP;  
00122   G4int NEX = QEX;
00123   G4LorentzVector ppout;
00124   G4bool try_again = (NEX > 0);
00125   
00126   // Buffer for parameter sets
00127   std::pair<G4double, G4double> parms;
00128   
00129   while (try_again) {
00130     if (A >= a_cut && Z >= z_cut && EEXS > eexs_cut) { // ok
00131       // update exiton system (include excitation energy!)
00132       G4double nuc_mass = G4InuclNuclei::getNucleiMass(A, Z, EEXS); 
00133       PEX.setVectM(PEX.vect(), nuc_mass);
00134       toTheExitonSystemRestFrame.setTarget(PEX);
00135       toTheExitonSystemRestFrame.toTheTargetRestFrame();
00136       
00137       if (verboseLevel > 2) {
00138         G4cout << " A " << A << " Z " << Z << " mass " << nuc_mass
00139                << " EEXS " << EEXS << G4endl; 
00140       }
00141       
00142       G4double MEL = getMatrixElement(A);
00143       G4double E0 = getE0(A);
00144       G4double PL = getParLev(A, Z);
00145       G4double parlev = PL / A;
00146       G4double EG = PL * EEXS;
00147       
00148       if (QEX < std::sqrt(2.0 * EG)) { // ok
00149         if (verboseLevel > 3)
00150           G4cout << " QEX " << QEX << " < sqrt(2*EG) " << std::sqrt(2.*EG)
00151                  << " NEX " << NEX << G4endl;
00152         
00153         paraMakerTruncated(Z, parms);
00154         const G4double& AK1 = parms.first;
00155         const G4double& CPA1 = parms.second;
00156         
00157         G4double VP = coul_coeff * Z * AK1 / (G4cbrt(A-1) + 1.0) /
00158           (1.0 + EEXS / E0);
00159         G4double DM1 = bindingEnergy(A,Z);
00160         G4double BN = DM1 - bindingEnergy(A-1,Z);
00161         G4double BP = DM1 - bindingEnergy(A-1,Z-1);
00162         G4double EMN = EEXS - BN;
00163         G4double EMP = EEXS - BP - VP * A / (A-1);
00164         G4double ESP = 0.0;
00165 
00166         if (verboseLevel > 3) {
00167           G4cout << " AK1 " << AK1 << " CPA1 " << " VP " << VP
00168                  << "\n bind(A,Z) " << DM1 << " dBind(N) " << BN 
00169                  << " dBind(P) " << BP
00170                  << "\n EMN " << EMN << " EMP " << EMP << G4endl;
00171         }
00172 
00173         if (EMN > eexs_cut) { // ok
00174           G4int icase = 0;
00175           
00176           if (NEX > 1) {
00177             G4double APH = 0.25 * (QP * QP + QH * QH + QP - 3 * QH);
00178             G4double APH1 = APH + 0.5 * (QP + QH);
00179             ESP = EEXS / QEX;
00180             G4double MELE = MEL / ESP / (A*A*A);
00181 
00182             if (verboseLevel > 3)
00183               G4cout << " APH " << APH << " APH1 " << APH1 << " ESP " << ESP
00184                      << G4endl;
00185 
00186             if (ESP > 15.0) {
00187               MELE *= std::sqrt(15.0 / ESP);
00188             } else if(ESP < 7.0) {
00189               MELE *= std::sqrt(ESP / 7.0);
00190               if (ESP < 2.0) MELE *= std::sqrt(ESP / 2.0);
00191             };    
00192 
00193             G4double F1 = EG - APH;
00194             G4double F2 = EG - APH1;
00195 
00196             if (verboseLevel > 3)
00197               G4cout << " MELE " << MELE << " F1 " << F1 << " F2 " << F2
00198                      << G4endl;
00199             
00200             if (F1 > 0.0 && F2 > 0.0) {
00201               G4double F = F2 / F1;
00202               G4double M1 = 2.77 * MELE * PL;
00203               G4double D[3] = { 0., 0., 0. };
00204               D[0] = M1 * F2 * F2 * std::pow(F, NEX-1) / (QEX+1);
00205               
00206               if (D[0] > 0.0) {
00207                 
00208                 if (NEX >= 2) {
00209                   D[1] = 0.0462 / parlev / G4cbrt(A) * QP * EEXS / QEX;
00210                   
00211                   if (EMP > eexs_cut) 
00212                     D[2] = D[1] * std::pow(EMP / EEXS, NEX) * (1.0 + CPA1);
00213                   D[1] *= std::pow(EMN / EEXS, NEX) * getAL(A);   
00214                   
00215                   if (QNP < 1) D[1] = 0.0;
00216                   if (QPP < 1) D[2] = 0.0;
00217                   
00218                   try_again = NEX > 1 && (D[1] > width_cut * D[0] || 
00219                                           D[2] > width_cut * D[0]);
00220                   
00221                   if (try_again) {
00222                     G4double D5 = D[0] + D[1] + D[2];
00223                     G4double SL = D5 * inuclRndm();
00224                     G4double S1 = 0.;
00225 
00226                     if (verboseLevel > 3)
00227                       G4cout << " D5 " << D5 << " SL " << SL << G4endl;
00228 
00229                     for (G4int i = 0; i < 3; i++) {
00230                       S1 += D[i];       
00231                       if (SL <= S1) {
00232                         icase = i;
00233                         break;
00234                       }
00235                     }
00236 
00237                     if (verboseLevel > 3)
00238                       G4cout << " got icase " << icase << G4endl;
00239                   }                             // if (try_again)
00240                 }                               // if (NEX >= 2)
00241               } else try_again = false;         // if (D[0] > 0)
00242             } else try_again = false;           // if (F1>0 && F2>0)
00243           }                                     // if (NEX > 1)
00244           
00245           if (try_again) {
00246             if (icase > 0) {                    // N -> N-1 with particle escape
00247               if (verboseLevel > 3)
00248                 G4cout << " try_again icase " << icase << G4endl;
00249               
00250               G4double V = 0.0;
00251               G4int ptype = 0;
00252               G4double B = 0.0;
00253               
00254               if (A < 3.0) try_again = false;
00255               
00256               if (try_again) { 
00257                 
00258                 if (icase == 1) { // neutron escape
00259                   if (verboseLevel > 3)
00260                     G4cout << " trying neutron escape" << G4endl;
00261 
00262                   if (QNP < 1) icase = 0;
00263                   else {
00264                     B = BN;
00265                     V = 0.0;
00266                     ptype = 2;            
00267                   };    
00268                 } else { // proton esape
00269                   if (verboseLevel > 3)
00270                     G4cout << " trying proton escape" << G4endl;
00271 
00272                   if (QPP < 1) icase = 0;
00273                   else {
00274                     B = BP;
00275                     V = VP;
00276                     ptype = 1;
00277                     
00278                     if (Z-1 < 1) try_again = false;
00279                   };   
00280                 };
00281                 
00282                 if (try_again && icase != 0) {
00283                   if (verboseLevel > 3)
00284                     G4cout << " ptype " << ptype << " B " << B << " V " << V
00285                            << G4endl;
00286 
00287                   G4double EB = EEXS - B;
00288                   G4double E = EB - V * A / (A-1);
00289                   
00290                   if (E < 0.0) icase = 0;
00291                   else {
00292                     G4double E1 = EB - V;
00293                     G4double EEXS_new = -1.;
00294                     G4double EPART = 0.0;
00295                     G4int itry1 = 0;
00296                     G4bool bad = true;
00297                     
00298                     while (itry1 < itry_max && icase > 0 && bad) {
00299                       itry1++;
00300                       G4int itry = 0;               
00301                       
00302                       while (EEXS_new < 0.0 && itry < itry_max) {
00303                         itry++;
00304                         G4double R = inuclRndm();
00305                         G4double X;
00306                         
00307                         if (NEX == 2) {
00308                           X = 1.0 - std::sqrt(R);
00309                           
00310                         } else {                         
00311                           G4double QEX2 = 1.0 / QEX;
00312                           G4double QEX1 = 1.0 / (QEX-1);
00313                           X = std::pow(0.5 * R, QEX2);
00314                           
00315                           for (G4int i = 0; i < 1000; i++) {
00316                             G4double DX = X * QEX1 * 
00317                               (1.0 + QEX2 * X * (1.0 - R / std::pow(X, NEX)) / (1.0 - X));
00318                             X -= DX;
00319                             
00320                             if (std::fabs(DX / X) < 0.01) break;  
00321                             
00322                           };
00323                         }; 
00324                         EPART = EB - X * E1;
00325                         EEXS_new = EB - EPART * A / (A-1);
00326                       } // while (EEXS_new < 0.0...
00327                       
00328                       if (itry == itry_max || EEXS_new < 0.0) {
00329                         icase = 0;
00330                         continue;
00331                       }
00332                       
00333                       if (verboseLevel > 2)
00334                         G4cout << " particle " << ptype << " escape " << G4endl;
00335                       
00336                       EPART /= GeV;             // From MeV to GeV
00337                       
00338                       G4InuclElementaryParticle particle(ptype);
00339                       particle.setModel(G4InuclParticle::NonEquilib);
00340                       
00341                       // generate particle momentum
00342                       G4double mass = particle.getMass();
00343                       G4double pmod = std::sqrt(EPART * (2.0 * mass + EPART));
00344                       G4LorentzVector mom = generateWithRandomAngles(pmod,mass);
00345                       
00346                       // Push evaporated paricle into current rest frame
00347                       mom = toTheExitonSystemRestFrame.backToTheLab(mom);
00348                       
00349                       // Adjust quasiparticle and nucleon counts
00350                       G4int QPP_new = QPP;
00351                       G4int QNP_new = QNP;
00352                       
00353                       G4int A_new = A-1;
00354                       G4int Z_new = Z;
00355                       
00356                       if (ptype == 1) {
00357                         QPP_new--;
00358                         Z_new--;
00359                       };
00360                       
00361                       if (ptype == 2) QNP_new--;
00362                       
00363                       if (verboseLevel > 3) {
00364                         G4cout << " nucleus   px " << PEX.px() << " py " << PEX.py()
00365                                << " pz " << PEX.pz() << " E " << PEX.e() << G4endl
00366                                << " evaporate px " << mom.px() << " py " << mom.py()
00367                                << " pz " << mom.pz() << " E " << mom.e() << G4endl;
00368                       }
00369                 
00370                       // New excitation energy depends on residual nuclear state
00371                       G4double mass_new = G4InuclNuclei::getNucleiMass(A_new, Z_new);
00372                       
00373                       EEXS_new = ((PEX-mom).m() - mass_new)*GeV;
00374                       if (EEXS_new < 0.) continue;      // Sanity check for new nucleus
00375                       
00376                       if (verboseLevel > 3)
00377                         G4cout << " EEXS_new " << EEXS_new << G4endl;
00378                       
00379                       PEX -= mom;
00380                       EEXS = EEXS_new;
00381                       
00382                       A = A_new;
00383                       Z = Z_new;
00384                       
00385                       NEX--;
00386                       QEX--;
00387                       QP--;
00388                       QPP = QPP_new;
00389                       QNP = QNP_new;
00390                       
00391                       particle.setMomentum(mom);
00392                       output.addOutgoingParticle(particle);
00393                       ppout += mom;
00394                       if (verboseLevel > 3) {
00395                         G4cout << particle << G4endl
00396                                << " ppout px " << ppout.px() << " py " << ppout.py()
00397                                << " pz " << ppout.pz() << " E " << ppout.e() << G4endl;
00398                       }
00399 
00400                       bad = false;
00401                     }           // while (itry1<itry_max && icase>0
00402                     
00403                     if (itry1 == itry_max) icase = 0;
00404                   }     // if (E < 0.) [else]
00405                 }       // if (try_again && icase != 0)
00406               }         // if (try_again)
00407             }           // if (icase > 0)
00408             
00409             if (icase == 0 && try_again) { // N -> N + 2 
00410               if (verboseLevel > 3) G4cout << " adding excitons" << G4endl;
00411 
00412               G4double TNN = 1.6 * EFN + ESP;
00413               G4double TNP = 1.6 * EFP + ESP;
00414               G4double XNUN = 1.0 / (1.6 + ESP / EFN);
00415               G4double XNUP = 1.0 / (1.6 + ESP / EFP);
00416               G4double SNN1 = csNN(TNP) * XNUP;
00417               G4double SNN2 = csNN(TNN) * XNUN;
00418               G4double SPN1 = csPN(TNP) * XNUP;
00419               G4double SPN2 = csPN(TNN) * XNUN;
00420               G4double PP = (QPP * SNN1 + QNP * SPN1) * ZR;
00421               G4double PN = (QPP * SPN2 + QNP * SNN2) * (AR - ZR);
00422               G4double PW = PP + PN;
00423               NEX += 2;
00424               QEX += 2; 
00425               QP++;
00426               QH++;
00427               AR--;
00428               
00429               if (AR > 1) {
00430                 G4double SL = PW * inuclRndm();
00431                 
00432                 if (SL > PP) {
00433                   QNP++;
00434                   QNH++;
00435                 } else {
00436                   QPP++;
00437                   QPH++;
00438                   ZR--;
00439                   if (ZR < 2) try_again = false;
00440                 }  
00441               } else try_again = false;
00442             }   // if (icase==0 && try_again)
00443           }     // if (try_again)
00444         } else try_again = false;       // if (EMN > eexs_cut)
00445       } else try_again = false;         // if (QEX < sqrg(2*EG)
00446     } else try_again = false;           // if (A > a_cut ...
00447   }             // while (try_again)
00448   
00449   // everything finished, set output nuclei
00450 
00451   if (output.numberOfOutgoingParticles() == 0) {
00452     output.addOutgoingNucleus(*nuclei_target);
00453   } else {
00454     G4LorentzVector pnuc = pin - ppout;
00455     G4InuclNuclei nuclei(pnuc, A, Z, EEXS, G4InuclParticle::NonEquilib);
00456     
00457     if (verboseLevel > 3) G4cout << " remaining nucleus\n" << nuclei << G4endl;
00458     output.addOutgoingNucleus(nuclei);
00459   }
00460 
00461   validateOutput(0, target, output);    // Check energy conservation, etc.
00462   return;
00463 }
00464 
00465 G4double G4NonEquilibriumEvaporator::getMatrixElement(G4int A) const {
00466 
00467   if (verboseLevel > 3) {
00468     G4cout << " >>> G4NonEquilibriumEvaporator::getMatrixElement" << G4endl;
00469   }
00470 
00471   G4double me;
00472 
00473   if (A > 150) me = 100.0;
00474   else if (A > 20) me = 140.0;
00475   else me = 70.0;
00476  
00477   return me;
00478 }
00479 
00480 G4double G4NonEquilibriumEvaporator::getE0(G4int ) const {
00481   if (verboseLevel > 3) {
00482     G4cout << " >>> G4NonEquilibriumEvaporator::getEO" << G4endl;
00483   }
00484 
00485   const G4double e0 = 200.0;
00486 
00487   return e0;   
00488 }
00489 
00490 G4double G4NonEquilibriumEvaporator::getParLev(G4int A, 
00491                                                G4int ) const {
00492 
00493   if (verboseLevel > 3) {
00494     G4cout << " >>> G4NonEquilibriumEvaporator::getParLev" << G4endl;
00495   }
00496 
00497   //  const G4double par = 0.125;
00498   G4double pl = 0.125 * A;
00499 
00500   return pl; 
00501 }

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