Ignore:
Timestamp:
Jun 18, 2010, 11:42:07 AM (14 years ago)
Author:
garnier
Message:

update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r962 r1315  
    2323// * acceptance of all terms of the Geant4 Software license.          *
    2424// ********************************************************************
     25// $Id: G4NonEquilibriumEvaporator.cc,v 1.29 2010/05/21 17:56:34 mkelsey Exp $
     26// Geant4 tag: $Name: geant4-09-04-beta-cand-01 $
    2527//
     28// 20100114  M. Kelsey -- Remove G4CascadeMomentum, use G4LorentzVector directly
     29// 20100309  M. Kelsey -- Use new generateWithRandomAngles for theta,phi stuff;
     30//              eliminate some unnecessary std::pow()
     31// 20100412  M. Kelsey -- Pass G4CollisionOutput by ref to ::collide()
     32// 20100413  M. Kelsey -- Pass buffers to paraMaker[Truncated]
     33// 20100517  M. Kelsey -- Inherit from common base class
     34
    2635#define RUN
    2736
    2837#include <cmath>
    2938#include "G4NonEquilibriumEvaporator.hh"
     39#include "G4CollisionOutput.hh"
    3040#include "G4InuclElementaryParticle.hh"
    3141#include "G4InuclNuclei.hh"
     42#include "G4InuclSpecialFunctions.hh"
    3243#include "G4LorentzConvertor.hh"
     44#include "G4NucleiProperties.hh"
     45#include "G4HadTmpUtil.hh"
     46
     47using namespace G4InuclSpecialFunctions;
     48
    3349
    3450G4NonEquilibriumEvaporator::G4NonEquilibriumEvaporator()
    35   : verboseLevel(1) {
    36 
    37   if (verboseLevel > 3) {
    38     G4cout << " >>> G4NonEquilibriumEvaporator::G4NonEquilibriumEvaporator" << G4endl;
    39   }
    40 }
    41 
    42 G4CollisionOutput G4NonEquilibriumEvaporator::collide(G4InuclParticle* /*bullet*/,
    43                                                       G4InuclParticle* target) {
     51  : G4VCascadeCollider("G4NonEquilibriumEvaporator") {}
     52
     53
     54void G4NonEquilibriumEvaporator::collide(G4InuclParticle* /*bullet*/,
     55                                         G4InuclParticle* target,
     56                                         G4CollisionOutput& output) {
    4457
    4558  if (verboseLevel > 3) {
    … …  
    4760  }
    4861
    49   const G4double one_third = 1.0/3.0;
     62  // Sanity check
     63  G4InuclNuclei* nuclei_target = dynamic_cast<G4InuclNuclei*>(target);
     64  if (!nuclei_target) {
     65    G4cerr << " NonEquilibriumEvaporator -> target is not nuclei " << G4endl;   
     66    return;
     67  }
     68
    5069  const G4double a_cut = 5.0;
    5170  const G4double z_cut = 3.0;
    … …  
    6180  const G4double small_ekin = 1.0e-6;
    6281  const G4double width_cut = 0.005;
    63   G4CollisionOutput output;
    64 
    65   if (G4InuclNuclei* nuclei_target = dynamic_cast<G4InuclNuclei*>(target)) {
    66     //  initialization
     82
    6783    G4double A = nuclei_target->getA();
    6884    G4double Z = nuclei_target->getZ();
    69     G4CascadeMomentum PEX = nuclei_target->getMomentum();
    70     G4CascadeMomentum pin = PEX;
     85    G4LorentzVector PEX = nuclei_target->getMomentum();
     86    G4LorentzVector pin = PEX;
    7187    G4double EEXS = nuclei_target->getExitationEnergy();
    72     pin[0] += 0.001 * EEXS;
     88    pin.setE(pin.e() + 0.001 * EEXS);
    7389    G4InuclNuclei dummy_nuc;
    7490    G4ExitonConfiguration config = nuclei_target->getExitonConfiguration(); 
    … …  
    88104    G4LorentzConvertor toTheExitonSystemRestFrame;
    89105
    90     toTheExitonSystemRestFrame.setBullet(dummy.getMomentum(), dummy.getMass());
     106    toTheExitonSystemRestFrame.setBullet(dummy);
    91107
    92108    G4double EFN = FermiEnergy(A, Z, 0);
    … …  
    96112    G4double ZR = Z - QPP; 
    97113    G4int NEX = G4int(QEX + 0.5);
    98     G4CascadeMomentum ppout;
    99     G4bool try_again = NEX > 0 ? true : false;
     114    G4LorentzVector ppout;
     115    G4bool try_again = (NEX > 0);
    100116 
     117    // Buffer for parameter sets
     118    std::pair<G4double, G4double> parms;
     119
    101120    while (try_again) {
    102121
    … …  
    109128        // update exiton system
    110129        G4double nuc_mass = dummy_nuc.getNucleiMass(A, Z);
    111         PEX[0] = std::sqrt(PEX[1] * PEX[1] + PEX[2] * PEX[2] + PEX[3] * PEX[3] +
    112                       nuc_mass * nuc_mass);     
     130        PEX.setVectM(PEX.vect(), nuc_mass);
    113131        toTheExitonSystemRestFrame.setTarget(PEX, nuc_mass);
    114132        toTheExitonSystemRestFrame.toTheTargetRestFrame();
    … …  
    121139        if (QEX < std::sqrt(2.0 * EG)) { // ok
    122140
    123           std::pair<G4double, G4double> parms = paraMakerTruncated(Z);
    124 
    125           G4double AK1 = parms.first;
    126           G4double CPA1 = parms.second;
    127           G4double VP = coul_coeff * Z * AK1 / (std::pow(A - 1.0, one_third) + 1.0) /
     141          paraMakerTruncated(Z, parms);
     142          const G4double& AK1 = parms.first;
     143          const G4double& CPA1 = parms.second;
     144
     145          G4double VP = coul_coeff * Z * AK1 / (G4cbrt(A - 1.0) + 1.0) /
    128146            (1.0 + EEXS / E0);
    129           G4double DM1 = bindingEnergy(A, Z);
    130           G4double BN = DM1 - bindingEnergy(A - 1.0, Z);
    131           G4double BP = DM1 - bindingEnergy(A - 1.0, Z - 1.0);
     147          //      G4double DM1 = bindingEnergy(A, Z);
     148          //      G4double BN = DM1 - bindingEnergy(A - 1.0, Z);
     149          //      G4double BP = DM1 - bindingEnergy(A - 1.0, Z - 1.0);
     150          G4double DM1 = G4NucleiProperties::GetBindingEnergy(G4lrint(A), G4lrint(Z));
     151          G4double BN = DM1 - G4NucleiProperties::GetBindingEnergy(G4lrint(A-1.0), G4lrint(Z));
     152          G4double BP = DM1 - G4NucleiProperties::GetBindingEnergy(G4lrint(A-1.0), G4lrint(Z-1.0));
    132153          G4double EMN = EEXS - BN;
    133154          G4double EMP = EEXS - BP - VP * A / (A - 1.0);
    … …  
    141162              G4double APH1 = APH + 0.5 * (QP + QH);
    142163              ESP = EEXS / QEX;
    143               G4double MELE = MEL / ESP / std::pow(A, 3);
     164              G4double MELE = MEL / ESP / (A*A*A);
    144165
    145166              if (ESP > 15.0) {
    … …  
    164185
    165186                  if (NEX >= 2) {
    166                     D[1] = 0.0462 / parlev / std::pow(A, one_third) * QP * EEXS / QEX;
     187                    D[1] = 0.0462 / parlev / G4cbrt(A) * QP * EEXS / QEX;
    167188
    168189                    if (EMP > eexs_cut)
    … …  
    293314                          // generate particle momentum
    294315                          G4double pmod = std::sqrt(EPART * (2.0 * mass + EPART));
    295                           G4CascadeMomentum mom;
    296                           std::pair<G4double, G4double> COS_SIN = randomCOS_SIN();
    297                           G4double FI = randomPHI();
    298                           G4double P1 = pmod * COS_SIN.second;
    299                           mom[1] = P1 * std::cos(FI);
    300                           mom[2] = P1 * std::sin(FI);
    301                           mom[3] = pmod * COS_SIN.first;
    302                           G4CascadeMomentum mom_at_rest;
    303 
    304                           for (G4int i = 1; i < 4; i++) mom_at_rest[i] = -mom[i];
     316                          G4LorentzVector mom = generateWithRandomAngles(pmod,mass);
     317                          G4LorentzVector mom_at_rest;
    305318
    306319                          G4double QPP_new = QPP;
    … …  
    319332                          G4double new_exiton_mass =
    320333                            dummy_nuc.getNucleiMass(A_new, Z_new);
    321                           mom_at_rest[0] = std::sqrt(mom_at_rest[1] * mom_at_rest[1] +
    322                                                 mom_at_rest[2] * mom_at_rest[2] +
    323                                                 mom_at_rest[3] * mom_at_rest[3] +
    324                                                 new_exiton_mass * new_exiton_mass);
    325                           mom[0] = std::sqrt(mom[1] * mom[1] + mom[2] * mom[2] +
    326                                         mom[3] * mom[3] + mass * mass);
    327 
    328                           G4CascadeMomentum part_mom =
     334                          mom_at_rest.setVectM(-mom.vect(), new_exiton_mass);
     335
     336                          G4LorentzVector part_mom =
    329337                            toTheExitonSystemRestFrame.backToTheLab(mom);
    330 
    331                           part_mom[0] = std::sqrt(part_mom[1] * part_mom[1] +
    332                                              part_mom[2] * part_mom[2] + part_mom[3] * part_mom[3] +
    333                                              mass * mass);
    334 
    335                           G4CascadeMomentum ex_mom =
     338                          part_mom.setVectM(part_mom.vect(), mass);
     339
     340                          G4LorentzVector ex_mom =
    336341                            toTheExitonSystemRestFrame.backToTheLab(mom_at_rest);
    337 
    338                           ex_mom[0] = std::sqrt(ex_mom[1] * ex_mom[1] + ex_mom[2] * ex_mom[2]
    339                                            + ex_mom[3] * ex_mom[3] + new_exiton_mass * new_exiton_mass);   
     342                          ex_mom.setVectM(ex_mom.vect(), new_exiton_mass);   
     343
    340344                          //             check energy conservation and set new exitation energy
    341                           EEXS_new = 1000.0 * (PEX[0] + 0.001 * EEXS -
    342                                                part_mom[0] - ex_mom[0]);
     345                          EEXS_new = 1000.0 * (PEX.e() + 0.001 * EEXS -
     346                                               part_mom.e() - ex_mom.e());
    343347
    344348                          if (EEXS_new > 0.0) { // everything ok
    345349                            particle.setMomentum(part_mom);
    346350                            output.addOutgoingParticle(particle);
    347 
    348                             for (G4int i = 0; i < 4; i++) ppout[i] += part_mom[i];
     351                            ppout += part_mom;
     352
    349353                            A = A_new;
    350354                            Z = Z_new;
    … …  
    423427    // conservation
    424428
    425     G4CascadeMomentum pnuc;
    426 
    427     for (G4int i = 1; i < 4; i++) pnuc[i] = pin[i] - ppout[i];
     429    G4LorentzVector pnuc = pin - ppout;
    428430    G4InuclNuclei nuclei(pnuc, A, Z);
    429431
    … …  
    432434
    433435    pnuc = nuclei.getMomentum();
    434     G4double eout = pnuc[0] + ppout[0]; 
    435     G4double eex_real = 1000.0 * (pin[0] - eout);       
     436    G4double eout = pnuc.e() + ppout.e(); 
     437    G4double eex_real = 1000.0 * (pin.e() - eout);       
    436438
    437439    nuclei.setExitationEnergy(eex_real);
    438440    output.addTargetFragment(nuclei);
    439441
    440   } else {
    441     G4cout << " NonEquilibriumEvaporator -> target is not nuclei " << G4endl;   
    442 
    443   };
    444 
    445   return output;
     442  return;
    446443}
    447444
Note: See TracChangeset for help on using the changeset viewer.