Ignore:
Timestamp:
Nov 5, 2010, 3:45:55 PM (14 years ago)
Author:
garnier
Message:

update ti head

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/source/processes/hadronic/models/cascade/cascade/src/G4NonEquilibriumEvaporator.cc

    r1337 r1340  
    2323// * acceptance of all terms of the Geant4 Software license.          *
    2424// ********************************************************************
    25 //
    26 // $Id: G4NonEquilibriumEvaporator.cc,v 1.29.2.1 2010/06/25 09:44:52 gunter Exp $
    27 // Geant4 tag: $Name: geant4-09-04-beta-01 $
     25// $Id: G4NonEquilibriumEvaporator.cc,v 1.40 2010/09/24 20:51:05 mkelsey Exp $
     26// Geant4 tag: $Name: hadr-casc-V09-03-85 $
    2827//
    2928// 20100114  M. Kelsey -- Remove G4CascadeMomentum, use G4LorentzVector directly
    … …  
    3332// 20100413  M. Kelsey -- Pass buffers to paraMaker[Truncated]
    3433// 20100517  M. Kelsey -- Inherit from common base class
    35 
    36 #define RUN
    37 
    38 #include <cmath>
     34// 20100617  M. Kelsey -- Remove "RUN" preprocessor flag and all "#else" code
     35// 20100622  M. Kelsey -- Use local "bindingEnergy()" function to call through.
     36// 20100701  M. Kelsey -- Don't need to add excitation to nuclear mass; compute
     37//              new excitation energies properly (mass differences)
     38// 20100713  M. Kelsey -- Add conservation checking, diagnostic messages.
     39// 20100714  M. Kelsey -- Move conservation checking to base class
     40// 20100719  M. Kelsey -- Simplify EEXS calculations with particle evaporation.
     41// 20100724  M. Kelsey -- Replace std::vector<> D with trivial D[3] array.
     42// 20100914  M. Kelsey -- Migrate to integer A and Z: this involves replacing
     43//              a number of G4double terms with G4int, with consequent casts.
     44
    3945#include "G4NonEquilibriumEvaporator.hh"
    4046#include "G4CollisionOutput.hh"
    … …  
    4349#include "G4InuclSpecialFunctions.hh"
    4450#include "G4LorentzConvertor.hh"
    45 #include "G4NucleiProperties.hh"
    46 #include "G4HadTmpUtil.hh"
     51#include <cmath>
    4752
    4853using namespace G4InuclSpecialFunctions;
    … …  
    5055
    5156G4NonEquilibriumEvaporator::G4NonEquilibriumEvaporator()
    52   : G4VCascadeCollider("G4NonEquilibriumEvaporator") {}
     57  : G4CascadeColliderBase("G4NonEquilibriumEvaporator") {}
    5358
    5459
    … …  
    5762                                         G4CollisionOutput& output) {
    5863
    59   if (verboseLevel > 3) {
     64  if (verboseLevel) {
    6065    G4cout << " >>> G4NonEquilibriumEvaporator::collide" << G4endl;
    6166  }
    … …  
    6873  }
    6974
    70   const G4double a_cut = 5.0;
    71   const G4double z_cut = 3.0;
    72 
    73 #ifdef RUN
     75  if (verboseLevel > 2) {
     76    G4cout << " evaporating target: " << G4endl;
     77    target->printParticle();
     78  }
     79 
     80  const G4int a_cut = 5;
     81  const G4int z_cut = 3;
     82
    7483  const G4double eexs_cut = 0.1;
    75 #else
    76   const G4double eexs_cut = 100000.0;
    77 #endif
    7884
    7985  const G4double coul_coeff = 1.4;
    … …  
    8288  const G4double width_cut = 0.005;
    8389
    84     G4double A = nuclei_target->getA();
    85     G4double Z = nuclei_target->getZ();
    86     G4LorentzVector PEX = nuclei_target->getMomentum();
    87     G4LorentzVector pin = PEX;
    88     G4double EEXS = nuclei_target->getExitationEnergy();
    89     pin.setE(pin.e() + 0.001 * EEXS);
    90     G4InuclNuclei dummy_nuc;
    91     G4ExitonConfiguration config = nuclei_target->getExitonConfiguration(); 
    92 
    93     G4double QPP = config.protonQuasiParticles;
    94 
    95     G4double QNP = config.neutronQuasiParticles;
    96 
    97     G4double QPH = config.protonHoles;
    98 
    99     G4double QNH = config.neutronHoles;
    100 
    101     G4double QP = QPP + QNP;
    102     G4double QH = QPH + QNH;
    103     G4double QEX = QP + QH;
    104     G4InuclElementaryParticle dummy(small_ekin, 1);
    105     G4LorentzConvertor toTheExitonSystemRestFrame;
    106 
    107     toTheExitonSystemRestFrame.setBullet(dummy);
    108 
    109     G4double EFN = FermiEnergy(A, Z, 0);
    110     G4double EFP = FermiEnergy(A, Z, 1);
    111 
    112     G4double AR = A - QP;
    113     G4double ZR = Z - QPP; 
    114     G4int NEX = G4int(QEX + 0.5);
    115     G4LorentzVector ppout;
    116     G4bool try_again = (NEX > 0);
    117  
    118     // Buffer for parameter sets
    119     std::pair<G4double, G4double> parms;
    120 
    121     while (try_again) {
    122 
    123       if (A >= a_cut && Z >= z_cut && EEXS > eexs_cut) { // ok
    124 
    125         if (verboseLevel > 2) {
    126           G4cout << " A " << A << " Z " << Z << " EEXS " << EEXS << G4endl;
    127         }
    128 
    129         // update exiton system
    130         G4double nuc_mass = dummy_nuc.getNucleiMass(A, Z);
    131         PEX.setVectM(PEX.vect(), nuc_mass);
    132         toTheExitonSystemRestFrame.setTarget(PEX, nuc_mass);
    133         toTheExitonSystemRestFrame.toTheTargetRestFrame();
    134         G4double MEL = getMatrixElement(A);
    135         G4double E0 = getE0(A);
    136         G4double PL = getParLev(A, Z);
    137         G4double parlev = PL / A;
    138         G4double EG = PL * EEXS;
    139 
    140         if (QEX < std::sqrt(2.0 * EG)) { // ok
    141 
    142           paraMakerTruncated(Z, parms);
    143           const G4double& AK1 = parms.first;
    144           const G4double& CPA1 = parms.second;
    145 
    146           G4double VP = coul_coeff * Z * AK1 / (G4cbrt(A - 1.0) + 1.0) /
    147             (1.0 + EEXS / E0);
    148           //      G4double DM1 = bindingEnergy(A, Z);
    149           //      G4double BN = DM1 - bindingEnergy(A - 1.0, Z);
    150           //      G4double BP = DM1 - bindingEnergy(A - 1.0, Z - 1.0);
    151           G4double DM1 = G4NucleiProperties::GetBindingEnergy(G4lrint(A), G4lrint(Z));
    152           G4double BN = DM1 - G4NucleiProperties::GetBindingEnergy(G4lrint(A-1.0), G4lrint(Z));
    153           G4double BP = DM1 - G4NucleiProperties::GetBindingEnergy(G4lrint(A-1.0), G4lrint(Z-1.0));
    154           G4double EMN = EEXS - BN;
    155           G4double EMP = EEXS - BP - VP * A / (A - 1.0);
    156           G4double ESP = 0.0;
     90  G4int A = nuclei_target->getA();
     91  G4int Z = nuclei_target->getZ();
     92 
     93  G4LorentzVector PEX = nuclei_target->getMomentum();
     94  G4LorentzVector pin = PEX;            // Save original four-vector for later
     95 
     96  G4double EEXS = nuclei_target->getExitationEnergy();
     97 
     98  G4ExitonConfiguration config = nuclei_target->getExitonConfiguration(); 
     99 
     100  G4int QPP = config.protonQuasiParticles;
     101  G4int QNP = config.neutronQuasiParticles;
     102  G4int QPH = config.protonHoles;
     103  G4int QNH = config.neutronHoles;
     104 
     105  G4int QP = QPP + QNP;
     106  G4int QH = QPH + QNH;
     107  G4int QEX = QP + QH;
     108 
     109  G4InuclElementaryParticle dummy(small_ekin, 1);
     110  G4LorentzConvertor toTheExitonSystemRestFrame;
     111  //*** toTheExitonSystemRestFrame.setVerbose(verboseLevel);
     112  toTheExitonSystemRestFrame.setBullet(dummy);
     113 
     114  G4double EFN = FermiEnergy(A, Z, 0);
     115  G4double EFP = FermiEnergy(A, Z, 1);
     116 
     117  G4int AR = A - QP;
     118  G4int ZR = Z - QPP; 
     119  G4int NEX = QEX;
     120  G4LorentzVector ppout;
     121  G4bool try_again = (NEX > 0);
     122 
     123  // Buffer for parameter sets
     124  std::pair<G4double, G4double> parms;
     125 
     126  while (try_again) {
     127    if (A >= a_cut && Z >= z_cut && EEXS > eexs_cut) { // ok
     128      // update exiton system (include excitation energy!)
     129      G4double nuc_mass = G4InuclNuclei::getNucleiMass(A, Z, EEXS);
     130      PEX.setVectM(PEX.vect(), nuc_mass);
     131      toTheExitonSystemRestFrame.setTarget(PEX);
     132      toTheExitonSystemRestFrame.toTheTargetRestFrame();
     133     
     134      if (verboseLevel > 2) {
     135        G4cout << " A " << A << " Z " << Z << " mass " << nuc_mass
     136               << " EEXS " << EEXS << G4endl;
     137      }
     138     
     139      G4double MEL = getMatrixElement(A);
     140      G4double E0 = getE0(A);
     141      G4double PL = getParLev(A, Z);
     142      G4double parlev = PL / A;
     143      G4double EG = PL * EEXS;
     144     
     145      if (QEX < std::sqrt(2.0 * EG)) { // ok
     146        if (verboseLevel > 3)
     147          G4cout << " QEX " << QEX << " < sqrt(2*EG) " << std::sqrt(2.*EG)
     148                 << G4endl;
    157149       
    158           if (EMN > eexs_cut) { // ok
    159             G4int icase = 0;
     150        paraMakerTruncated(Z, parms);
     151        const G4double& AK1 = parms.first;
     152        const G4double& CPA1 = parms.second;
     153       
     154        G4double VP = coul_coeff * Z * AK1 / (G4cbrt(A-1) + 1.0) /
     155          (1.0 + EEXS / E0);
     156        G4double DM1 = bindingEnergy(A,Z);
     157        G4double BN = DM1 - bindingEnergy(A-1,Z);
     158        G4double BP = DM1 - bindingEnergy(A-1,Z-1);
     159        G4double EMN = EEXS - BN;
     160        G4double EMP = EEXS - BP - VP * A / (A-1);
     161        G4double ESP = 0.0;
     162       
     163        if (EMN > eexs_cut) { // ok
     164          G4int icase = 0;
    160165         
    161             if (NEX > 1) {
    162               G4double APH = 0.25 * (QP * QP + QH * QH + QP - 3.0 * QH);
    163               G4double APH1 = APH + 0.5 * (QP + QH);
    164               ESP = EEXS / QEX;
    165               G4double MELE = MEL / ESP / (A*A*A);
    166 
    167               if (ESP > 15.0) {
    168                 MELE *= std::sqrt(15.0 / ESP);
    169 
    170               } else if(ESP < 7.0) {
    171                 MELE *= std::sqrt(ESP / 7.0);
    172 
    173                 if (ESP < 2.0) MELE *= std::sqrt(ESP / 2.0);
    174               };   
     166          if (NEX > 1) {
     167            G4double APH = 0.25 * (QP * QP + QH * QH + QP - 3 * QH);
     168            G4double APH1 = APH + 0.5 * (QP + QH);
     169            ESP = EEXS / QEX;
     170            G4double MELE = MEL / ESP / (A*A*A);
     171           
     172            if (ESP > 15.0) {
     173              MELE *= std::sqrt(15.0 / ESP);
     174             
     175            } else if(ESP < 7.0) {
     176              MELE *= std::sqrt(ESP / 7.0);
     177             
     178              if (ESP < 2.0) MELE *= std::sqrt(ESP / 2.0);
     179            };   
    175180           
    176               G4double F1 = EG - APH;
    177               G4double F2 = EG - APH1;
     181            G4double F1 = EG - APH;
     182            G4double F2 = EG - APH1;
    178183           
    179               if (F1 > 0.0 && F2 > 0.0) {
    180                 G4double F = F2 / F1;
    181                 G4double M1 = 2.77 * MELE * PL;
    182                 std::vector<G4double> D(3, 0.0);
    183                 D[0] = M1 * F2 * F2 * std::pow(F, NEX - 1) / (QEX + 1.0);
    184 
    185                 if (D[0] > 0.0) {
    186 
    187                   if (NEX >= 2) {
    188                     D[1] = 0.0462 / parlev / G4cbrt(A) * QP * EEXS / QEX;
    189 
    190                     if (EMP > eexs_cut)
    191                       D[2] = D[1] * std::pow(EMP / EEXS, NEX) * (1.0 + CPA1);
    192                     D[1] *= std::pow(EMN / EEXS, NEX) * getAL(A);   
    193 
    194                     if (QNP < 1.0) D[1] = 0.0;
    195 
    196                     if (QPP < 1.0) D[2] = 0.0;
    197 
    198                     try_again = NEX > 1 && (D[1] > width_cut * D[0] ||
    199                                             D[2] > width_cut * D[0]);
    200 
    201                     if (try_again) {
    202                       G4double D5 = D[0] + D[1] + D[2];
    203                       G4double SL = D5 * inuclRndm();
    204                       G4double S1 = 0.;
    205 
    206                       for (G4int i = 0; i < 3; i++) {
    207                         S1 += D[i];     
    208 
    209                         if (SL <= S1) {
    210                           icase = i;
    211 
    212                           break;
    213                         };
     184            if (F1 > 0.0 && F2 > 0.0) {
     185              G4double F = F2 / F1;
     186              G4double M1 = 2.77 * MELE * PL;
     187              G4double D[3] = { 0., 0., 0. };
     188              D[0] = M1 * F2 * F2 * std::pow(F, NEX-1) / (QEX+1);
     189             
     190              if (D[0] > 0.0) {
     191               
     192                if (NEX >= 2) {
     193                  D[1] = 0.0462 / parlev / G4cbrt(A) * QP * EEXS / QEX;
     194                 
     195                  if (EMP > eexs_cut)
     196                    D[2] = D[1] * std::pow(EMP / EEXS, NEX) * (1.0 + CPA1);
     197                  D[1] *= std::pow(EMN / EEXS, NEX) * getAL(A);   
     198                 
     199                  if (QNP < 1) D[1] = 0.0;
     200                  if (QPP < 1) D[2] = 0.0;
     201                 
     202                  try_again = NEX > 1 && (D[1] > width_cut * D[0] ||
     203                                          D[2] > width_cut * D[0]);
     204                 
     205                  if (try_again) {
     206                    G4double D5 = D[0] + D[1] + D[2];
     207                    G4double SL = D5 * inuclRndm();
     208                    G4double S1 = 0.;
     209                   
     210                    for (G4int i = 0; i < 3; i++) {
     211                      S1 += D[i];       
     212                      if (SL <= S1) {
     213                        icase = i;
     214                        break;
    214215                      };
    215                     };
    216                   };
    217 
     216                    };
     217                  };
     218                };
     219              } else try_again = false;
     220            } else try_again = false;
     221          }
     222         
     223          if (try_again) {
     224            if (icase > 0) { // N -> N - 1 with particle escape
     225              if (verboseLevel > 3)
     226                G4cout << " try_again icase " << icase << G4endl;
     227             
     228              G4double V = 0.0;
     229              G4int ptype = 0;
     230              G4double B = 0.0;
     231             
     232              if (A < 3.0) try_again = false;
     233             
     234              if (try_again) {
     235               
     236                if (icase == 1) { // neutron escape
     237                 
     238                  if (QNP < 1) {
     239                    icase = 0;
     240                   
     241                  } else {
     242                    B = BN;
     243                    V = 0.0;
     244                    ptype = 2;           
     245                  };   
     246                 
     247                } else { // proton esape
     248                  if (QPP < 1) {
     249                    icase = 0;
     250                   
     251                  } else {
     252                    B = BP;
     253                    V = VP;
     254                    ptype = 1;
     255                   
     256                    if (Z-1 < 1) try_again = false;
     257                  };   
     258                };
     259               
     260                if (try_again && icase != 0) {
     261                  G4double EB = EEXS - B;
     262                  G4double E = EB - V * A / (A-1);
     263                 
     264                  if (E < 0.0) icase = 0;
     265                  else {
     266                    G4double E1 = EB - V;
     267                    G4double EEXS_new = -1.;
     268                    G4double EPART = 0.0;
     269                    G4int itry1 = 0;
     270                    G4bool bad = true;
     271                   
     272                    while (itry1 < itry_max && icase > 0 && bad) {
     273                      itry1++;
     274                      G4int itry = 0;               
     275                     
     276                      while (EEXS_new < 0.0 && itry < itry_max) {
     277                        itry++;
     278                        G4double R = inuclRndm();
     279                        G4double X;
     280                       
     281                        if (NEX == 2) {
     282                          X = 1.0 - std::sqrt(R);
     283                         
     284                        } else {                         
     285                          G4double QEX2 = 1.0 / QEX;
     286                          G4double QEX1 = 1.0 / (QEX-1);
     287                          X = std::pow(0.5 * R, QEX2);
     288                         
     289                          for (G4int i = 0; i < 1000; i++) {
     290                            G4double DX = X * QEX1 *
     291                              (1.0 + QEX2 * X * (1.0 - R / std::pow(X, NEX)) / (1.0 - X));
     292                            X -= DX;
     293                           
     294                            if (std::fabs(DX / X) < 0.01) break; 
     295                           
     296                          };
     297                        };
     298                        EPART = EB - X * E1;
     299                        EEXS_new = EB - EPART * A / (A-1);
     300                      } // while (EEXS_new < 0.0...
     301                     
     302                      if (itry == itry_max || EEXS_new < 0.0) {
     303                        icase = 0;
     304                        continue;
     305                      }
     306                     
     307                      if (verboseLevel > 2)
     308                        G4cout << " particle " << ptype << " escape " << G4endl;
     309                     
     310                      EPART /= GeV;             // From MeV to GeV
     311                     
     312                      G4InuclElementaryParticle particle(ptype);
     313                      particle.setModel(5);
     314                     
     315                      // generate particle momentum
     316                      G4double mass = particle.getMass();
     317                      G4double pmod = std::sqrt(EPART * (2.0 * mass + EPART));
     318                      G4LorentzVector mom = generateWithRandomAngles(pmod,mass);
     319                     
     320                      // Push evaporated paricle into current rest frame
     321                      mom = toTheExitonSystemRestFrame.backToTheLab(mom);
     322                     
     323                      // Adjust quasiparticle and nucleon counts
     324                      G4int QPP_new = QPP;
     325                      G4int QNP_new = QNP;
     326                     
     327                      G4int A_new = A-1;
     328                      G4int Z_new = Z;
     329                     
     330                      if (ptype == 1) {
     331                        QPP_new--;
     332                        Z_new--;
     333                      };
     334                     
     335                      if (ptype == 2) QNP_new--;
     336                     
     337                if (verboseLevel > 3) {
     338                  G4cout << " nucleus   px " << PEX.px() << " py " << PEX.py()
     339                         << " pz " << PEX.pz() << " E " << PEX.e() << G4endl
     340                         << " evaporate px " << mom.px() << " py " << mom.py()
     341                         << " pz " << mom.pz() << " E " << mom.e() << G4endl;
     342                }
     343               
     344                      // New excitation energy depends on residual nuclear state
     345                      G4double mass_new = G4InuclNuclei::getNucleiMass(A_new, Z_new);
     346                     
     347                      G4double EEXS_new = ((PEX-mom).m() - mass_new)*GeV;
     348                      if (EEXS_new < 0.) continue;      // Sanity check for new nucleus
     349                     
     350                      if (verboseLevel > 3)
     351                        G4cout << " EEXS_new " << EEXS_new << G4endl;
     352                     
     353                      PEX -= mom;
     354                      EEXS = EEXS_new;
     355                     
     356                      A = A_new;
     357                      Z = Z_new;
     358                     
     359                      NEX--;
     360                      QEX--;
     361                      QP--;
     362                      QPP = QPP_new;
     363                      QNP = QNP_new;
     364                     
     365                      particle.setMomentum(mom);
     366                      output.addOutgoingParticle(particle);
     367                      if (verboseLevel > 3) particle.printParticle();
     368                     
     369                      ppout += mom;
     370                      if (verboseLevel > 3) {
     371                        G4cout << " ppout px " << ppout.px() << " py " << ppout.py()
     372                               << " pz " << ppout.pz() << " E " << ppout.e() << G4endl;
     373                      }
     374
     375                      bad = false;
     376                    }           // while (itry1<itry_max && icase>0
     377                   
     378                    if (itry1 == itry_max) icase = 0;
     379                  }     // if (E < 0.) [else]
     380                }       // if (try_again && icase != 0)
     381              }         // if (try_again)
     382            }           // if (icase > 0)
     383           
     384            if (icase == 0 && try_again) { // N -> N + 2
     385              G4double TNN = 1.6 * EFN + ESP;
     386              G4double TNP = 1.6 * EFP + ESP;
     387              G4double XNUN = 1.0 / (1.6 + ESP / EFN);
     388              G4double XNUP = 1.0 / (1.6 + ESP / EFP);
     389              G4double SNN1 = csNN(TNP) * XNUP;
     390              G4double SNN2 = csNN(TNN) * XNUN;
     391              G4double SPN1 = csPN(TNP) * XNUP;
     392              G4double SPN2 = csPN(TNN) * XNUN;
     393              G4double PP = (QPP * SNN1 + QNP * SPN1) * ZR;
     394              G4double PN = (QPP * SPN2 + QNP * SNN2) * (AR - ZR);
     395              G4double PW = PP + PN;
     396              NEX += 2;
     397              QEX += 2;
     398              QP++;
     399              QH++;
     400              AR--;
     401             
     402              if (AR > 1) {
     403                G4double SL = PW * inuclRndm();
     404               
     405                if (SL > PP) {
     406                  QNP++;
     407                  QNH++;
    218408                } else {
    219                   try_again = false;
    220                 };
    221 
    222               } else {
    223                 try_again = false;
    224               };
    225             };           
    226 
    227             if (try_again) {
    228 
    229               if (icase > 0) { // N -> N - 1 with particle escape           
    230                 G4double V = 0.0;
    231                 G4int ptype = 0;
    232                 G4double B = 0.0;
    233 
    234                 if (A < 3.0) try_again = false;
    235 
    236                 if (try_again) {
    237 
    238                   if (icase == 1) { // neutron escape
    239 
    240                     if (QNP < 1.0) {
    241                       icase = 0;
    242 
    243                     } else {
    244                       B = BN;
    245                       V = 0.0;
    246                       ptype = 2;                 
    247                     };   
    248 
    249                   } else { // proton esape
    250                     if (QPP < 1.0) {
    251                       icase = 0;
    252 
    253                     } else {
    254                       B = BP;
    255                       V = VP;
    256                       ptype = 1;
    257 
    258                       if (Z - 1.0 < 1.0) try_again = false;
    259                     };   
    260                   };
    261                
    262                   if (try_again && icase != 0) {
    263                     G4double EB = EEXS - B;
    264                     G4double E = EB - V * A / (A - 1.0);
    265 
    266                     if (E < 0.0) {
    267                       icase = 0;
    268 
    269                     } else {
    270                       G4double E1 = EB - V;
    271                       G4double EEXS_new = -1.;
    272                       G4double EPART = 0.0;
    273                       G4int itry1 = 0;
    274                       G4bool bad = true;
    275 
    276                       while (itry1 < itry_max && icase > 0 && bad) {
    277                         itry1++;
    278                         G4int itry = 0;             
    279 
    280                         while (EEXS_new < 0.0 && itry < itry_max) {
    281                           itry++;
    282                           G4double R = inuclRndm();
    283                           G4double X;
    284 
    285                           if (NEX == 2) {
    286                             X = 1.0 - std::sqrt(R);
    287 
    288                           } else {                       
    289                             G4double QEX2 = 1.0 / QEX;
    290                             G4double QEX1 = 1.0 / (QEX - 1.0);
    291                             X = std::pow(0.5 * R, QEX2);
    292 
    293                             for (G4int i = 0; i < 1000; i++) {
    294                               G4double DX = X * QEX1 *
    295                                 (1.0 + QEX2 * X * (1.0 - R / std::pow(X, NEX)) / (1.0 - X));
    296                               X -= DX;
    297 
    298                               if (std::fabs(DX / X) < 0.01) break; 
    299 
    300                             };
    301                           };
    302                           EPART = EB - X * E1;
    303                           EEXS_new = EB - EPART * A / (A - 1.0);
    304                         };
    305                    
    306                         if (itry == itry_max || EEXS_new < 0.0) {
    307                           icase = 0;
    308 
    309                         } else { // real escape                 
    310                           G4InuclElementaryParticle particle(ptype);
    311 
    312                           particle.setModel(5);
    313                           G4double mass = particle.getMass();
    314                           EPART *= 0.001; // to the GeV
    315                           // generate particle momentum
    316                           G4double pmod = std::sqrt(EPART * (2.0 * mass + EPART));
    317                           G4LorentzVector mom = generateWithRandomAngles(pmod,mass);
    318                           G4LorentzVector mom_at_rest;
    319 
    320                           G4double QPP_new = QPP;
    321                           G4double Z_new = Z;
    322                        
    323                           if (ptype == 1) {
    324                             QPP_new -= 1.;
    325                             Z_new -= 1.0;
    326                           };
    327 
    328                           G4double QNP_new = QNP;
    329 
    330                           if (ptype == 2) QNP_new -= 1.0;
    331 
    332                           G4double A_new = A - 1.0;
    333                           G4double new_exiton_mass =
    334                             dummy_nuc.getNucleiMass(A_new, Z_new);
    335                           mom_at_rest.setVectM(-mom.vect(), new_exiton_mass);
    336 
    337                           G4LorentzVector part_mom =
    338                             toTheExitonSystemRestFrame.backToTheLab(mom);
    339                           part_mom.setVectM(part_mom.vect(), mass);
    340 
    341                           G4LorentzVector ex_mom =
    342                             toTheExitonSystemRestFrame.backToTheLab(mom_at_rest);
    343                           ex_mom.setVectM(ex_mom.vect(), new_exiton_mass);   
    344 
    345                           //             check energy conservation and set new exitation energy
    346                           EEXS_new = 1000.0 * (PEX.e() + 0.001 * EEXS -
    347                                                part_mom.e() - ex_mom.e());
    348 
    349                           if (EEXS_new > 0.0) { // everything ok
    350                             particle.setMomentum(part_mom);
    351                             output.addOutgoingParticle(particle);
    352                             ppout += part_mom;
    353 
    354                             A = A_new;
    355                             Z = Z_new;
    356                             PEX = ex_mom;
    357                             EEXS = EEXS_new;
    358                             NEX -= 1;
    359                             QEX -= 1;
    360                             QP -= 1.0;
    361                             QPP = QPP_new;
    362                             QNP = QNP_new;
    363                             bad = false;
    364                           };
    365                         };     
    366                       };   
    367 
    368                       if (itry1 == itry_max) icase = 0;
    369                     };   
    370                   };
    371                 };
    372               };
    373 
    374               if (icase == 0 && try_again) { // N -> N + 2
    375                 G4double TNN = 1.6 * EFN + ESP;
    376                 G4double TNP = 1.6 * EFP + ESP;
    377                 G4double XNUN = 1.0 / (1.6 + ESP / EFN);
    378                 G4double XNUP = 1.0 / (1.6 + ESP / EFP);
    379                 G4double SNN1 = csNN(TNP) * XNUP;
    380                 G4double SNN2 = csNN(TNN) * XNUN;
    381                 G4double SPN1 = csPN(TNP) * XNUP;
    382                 G4double SPN2 = csPN(TNN) * XNUN;
    383                 G4double PP = (QPP * SNN1 + QNP * SPN1) * ZR;
    384                 G4double PN = (QPP * SPN2 + QNP * SNN2) * (AR - ZR);
    385                 G4double PW = PP + PN;
    386                 NEX += 2;
    387                 QEX += 2.0;
    388                 QP += 1.0;
    389                 QH += 1.0;
    390                 AR -= 1.0;
    391 
    392                 if (AR > 1.0) {
    393                   G4double SL = PW * inuclRndm();
    394 
    395                   if (SL > PP) {
    396                     QNP += 1.0;
    397                     QNH += 1.0;
    398 
    399                   } else {
    400                     QPP += 1.0;
    401                     QPH += 1.0;
    402                     ZR -= 1.0;
    403 
    404                     if (ZR < 2.0) try_again = false;
    405 
    406                   };   
    407 
    408                 } else {
    409                   try_again = false;
    410                 };
    411               };
    412             };
    413 
    414           } else {
    415             try_again = false;
    416           };
    417 
    418         } else {
    419           try_again = false;
    420         };
    421 
    422       } else {
    423         try_again = false;
    424       }; 
    425     };
    426     // everything finished, set output nuclei
    427     // the exitation energy has to be re-set properly for the energy
    428     // conservation
    429 
     409                  QPP++;
     410                  QPH++;
     411                  ZR--;
     412                  if (ZR < 2) try_again = false;
     413                } 
     414              } else try_again = false;
     415            }   // if (icase==0 && try_again)
     416          }     // if (try_again)
     417        } else try_again = false;       // if (EMN > eexs_cut)
     418      } else try_again = false;         // if (QEX < sqrg(2*EG)
     419    } else try_again = false;           // if (A > a_cut ...
     420  }             // while (try_again)
     421 
     422  // everything finished, set output nuclei
     423
     424  if (output.numberOfOutgoingParticles() == 0) {
     425    output.addOutgoingNucleus(*nuclei_target);
     426  } else {
    430427    G4LorentzVector pnuc = pin - ppout;
    431     G4InuclNuclei nuclei(pnuc, A, Z);
    432 
    433     nuclei.setModel(5);
    434     nuclei.setEnergy();
    435 
    436     pnuc = nuclei.getMomentum();
    437     G4double eout = pnuc.e() + ppout.e(); 
    438     G4double eex_real = 1000.0 * (pin.e() - eout);       
    439 
    440     nuclei.setExitationEnergy(eex_real);
    441     output.addTargetFragment(nuclei);
    442 
     428    G4InuclNuclei nuclei(pnuc, A, Z, EEXS, 5);
     429   
     430    if (verboseLevel > 3) {
     431      G4cout << " remaining nucleus " << G4endl;
     432      nuclei.printParticle();
     433    }
     434    output.addOutgoingNucleus(nuclei);
     435  }
     436
     437  validateOutput(0, target, output);    // Check energy conservation, etc.
    443438  return;
    444439}
    445440
    446 G4double G4NonEquilibriumEvaporator::getMatrixElement(G4double A) const {
     441G4double G4NonEquilibriumEvaporator::getMatrixElement(G4int A) const {
    447442
    448443  if (verboseLevel > 3) {
    … …  
    452447  G4double me;
    453448
    454   if (A > 150.0) {
    455     me = 100.0;
    456 
    457   } else if(A > 20.0) {
    458     me = 140.0;
    459 
    460   } else {
    461     me = 70.0;
    462   };
     449  if (A > 150) me = 100.0;
     450  else if (A > 20) me = 140.0;
     451  else me = 70.0;
    463452 
    464453  return me;
    465454}
    466455
    467 G4double G4NonEquilibriumEvaporator::getE0(G4double ) const {
    468 
     456G4double G4NonEquilibriumEvaporator::getE0(G4int ) const {
    469457  if (verboseLevel > 3) {
    470458    G4cout << " >>> G4NonEquilibriumEvaporator::getEO" << G4endl;
    … …  
    476464}
    477465
    478 G4double G4NonEquilibriumEvaporator::getParLev(G4double A,
    479                                                G4double ) const {
     466G4double G4NonEquilibriumEvaporator::getParLev(G4int A,
     467                                               G4int ) const {
    480468
    481469  if (verboseLevel > 3) {
Note: See TracChangeset for help on using the changeset viewer.