[819] | 1 | // |
---|
| 2 | // ******************************************************************** |
---|
| 3 | // * License and Disclaimer * |
---|
| 4 | // * * |
---|
| 5 | // * The Geant4 software is copyright of the Copyright Holders of * |
---|
| 6 | // * the Geant4 Collaboration. It is provided under the terms and * |
---|
| 7 | // * conditions of the Geant4 Software License, included in the file * |
---|
| 8 | // * LICENSE and available at http://cern.ch/geant4/license . These * |
---|
| 9 | // * include a list of copyright holders. * |
---|
| 10 | // * * |
---|
| 11 | // * Neither the authors of this software system, nor their employing * |
---|
| 12 | // * institutes,nor the agencies providing financial support for this * |
---|
| 13 | // * work make any representation or warranty, express or implied, * |
---|
| 14 | // * regarding this software system or assume any liability for its * |
---|
| 15 | // * use. Please see the license in the file LICENSE and URL above * |
---|
| 16 | // * for the full disclaimer and the limitation of liability. * |
---|
| 17 | // * * |
---|
| 18 | // * This code implementation is the result of the scientific and * |
---|
| 19 | // * technical work of the GEANT4 collaboration. * |
---|
| 20 | // * By using, copying, modifying or distributing the software (or * |
---|
| 21 | // * any work based on the software) you agree to acknowledge its * |
---|
| 22 | // * use in resulting scientific publications, and indicate your * |
---|
| 23 | // * acceptance of all terms of the Geant4 Software license. * |
---|
| 24 | // ******************************************************************** |
---|
| 25 | // |
---|
[1196] | 26 | // $Id: G4Abla.cc,v 1.20 2009/11/18 10:43:14 kaitanie Exp $ |
---|
[819] | 27 | // Translation of INCL4.2/ABLA V3 |
---|
| 28 | // Pekka Kaitaniemi, HIP (translation) |
---|
| 29 | // Christelle Schmidt, IPNL (fission code) |
---|
| 30 | // Alain Boudard, CEA (contact person INCL/ABLA) |
---|
| 31 | // Aatos Heikkinen, HIP (project coordination) |
---|
| 32 | |
---|
| 33 | #include <time.h> |
---|
| 34 | |
---|
| 35 | #include "G4Abla.hh" |
---|
| 36 | #include "G4InclAblaDataFile.hh" |
---|
| 37 | #include "Randomize.hh" |
---|
[962] | 38 | #include "G4InclRandomNumbers.hh" |
---|
| 39 | #include "G4Ranecu.hh" |
---|
| 40 | #include "G4AblaFissionSimfis18.hh" |
---|
| 41 | #include "G4AblaFission.hh" |
---|
[819] | 42 | |
---|
| 43 | G4Abla::G4Abla() |
---|
| 44 | { |
---|
[962] | 45 | verboseLevel = 0; |
---|
| 46 | ilast = 0; |
---|
[819] | 47 | } |
---|
| 48 | |
---|
| 49 | G4Abla::G4Abla(G4Hazard *hazard, G4Volant *volant) |
---|
| 50 | { |
---|
| 51 | verboseLevel = 0; |
---|
[962] | 52 | ilast = 0; |
---|
[819] | 53 | volant = volant; // ABLA internal particle data |
---|
| 54 | volant->iv = 0; |
---|
| 55 | hazard = hazard; // Random seeds |
---|
| 56 | |
---|
[962] | 57 | randomGenerator = new G4InclGeant4Random(); |
---|
| 58 | //randomGenerator = new G4Ranecu(); |
---|
[819] | 59 | varntp = new G4VarNtp(); |
---|
| 60 | pace = new G4Pace(); |
---|
| 61 | ald = new G4Ald(); |
---|
| 62 | eenuc = new G4Eenuc(); |
---|
| 63 | ec2sub = new G4Ec2sub(); |
---|
| 64 | ecld = new G4Ecld(); |
---|
| 65 | fb = new G4Fb(); |
---|
| 66 | fiss = new G4Fiss(); |
---|
| 67 | opt = new G4Opt(); |
---|
| 68 | } |
---|
| 69 | |
---|
| 70 | G4Abla::G4Abla(G4Hazard *aHazard, G4Volant *aVolant, G4VarNtp *aVarntp) |
---|
| 71 | { |
---|
| 72 | verboseLevel = 0; |
---|
[962] | 73 | ilast = 0; |
---|
[819] | 74 | volant = aVolant; // ABLA internal particle data |
---|
| 75 | volant->iv = 0; |
---|
| 76 | hazard = aHazard; // Random seeds |
---|
| 77 | varntp = aVarntp; // Output data structure |
---|
| 78 | varntp->ntrack = 0; |
---|
| 79 | |
---|
[962] | 80 | randomGenerator = new G4InclGeant4Random(); |
---|
| 81 | //randomGenerator = new G4Ranecu(); |
---|
| 82 | |
---|
| 83 | // ABLA fission |
---|
| 84 | // fissionModel = new G4AblaFissionSimfis18(hazard, randomGenerator); |
---|
| 85 | fissionModel = new G4AblaFission(hazard, randomGenerator); |
---|
| 86 | if(verboseLevel > 0) { |
---|
| 87 | fissionModel->about(); |
---|
| 88 | } |
---|
| 89 | verboseLevel = 0; |
---|
[819] | 90 | pace = new G4Pace(); |
---|
| 91 | ald = new G4Ald(); |
---|
| 92 | eenuc = new G4Eenuc(); |
---|
| 93 | ec2sub = new G4Ec2sub(); |
---|
| 94 | ecld = new G4Ecld(); |
---|
| 95 | fb = new G4Fb(); |
---|
| 96 | fiss = new G4Fiss(); |
---|
| 97 | opt = new G4Opt(); |
---|
| 98 | } |
---|
| 99 | |
---|
| 100 | G4Abla::~G4Abla() |
---|
| 101 | { |
---|
[962] | 102 | delete randomGenerator; |
---|
[819] | 103 | delete pace; |
---|
| 104 | delete ald; |
---|
| 105 | delete eenuc; |
---|
| 106 | delete ec2sub; |
---|
| 107 | delete ecld; |
---|
| 108 | delete fb; |
---|
| 109 | delete fiss; |
---|
| 110 | delete opt; |
---|
| 111 | } |
---|
| 112 | |
---|
| 113 | // Main interface to the evaporation |
---|
| 114 | |
---|
| 115 | // Possible problem with generic Geant4 interface: ABLA evaporation |
---|
| 116 | // needs angular momentum information (calculated by INCL) to |
---|
| 117 | // work. Maybe there is a way to obtain this information from |
---|
| 118 | // G4Fragment? |
---|
| 119 | |
---|
| 120 | void G4Abla::breakItUp(G4double nucleusA, G4double nucleusZ, G4double nucleusMass, G4double excitationEnergy, |
---|
| 121 | G4double angularMomentum, G4double recoilEnergy, G4double momX, G4double momY, G4double momZ, |
---|
| 122 | G4int eventnumber) |
---|
| 123 | { |
---|
| 124 | const G4double uma = 931.4942; |
---|
| 125 | const G4double melec = 0.511; |
---|
| 126 | const G4double fmp = 938.27231; |
---|
| 127 | const G4double fmn = 939.56563; |
---|
| 128 | |
---|
[962] | 129 | G4double alrem = 0.0, berem = 0.0, garem = 0.0; |
---|
| 130 | G4double R[4][4]; // Rotation matrix |
---|
[819] | 131 | G4double csdir1[4]; |
---|
| 132 | G4double csdir2[4]; |
---|
| 133 | G4double csrem[4]; |
---|
| 134 | G4double pfis_rem[4]; |
---|
| 135 | G4double pf1_rem[4]; |
---|
[962] | 136 | for(G4int init_i = 0; init_i < 4; init_i++) { |
---|
| 137 | csdir1[init_i] = 0.0; |
---|
| 138 | csdir2[init_i] = 0.0; |
---|
| 139 | csrem[init_i] = 0.0; |
---|
| 140 | pfis_rem[init_i] = 0.0; |
---|
| 141 | pf1_rem[init_i] = 0.0; |
---|
| 142 | for(G4int init_j = 0; init_j < 4; init_j++) { |
---|
| 143 | R[init_i][init_j] = 0.0; |
---|
| 144 | } |
---|
| 145 | } |
---|
[819] | 146 | |
---|
[962] | 147 | G4double plab1 = 0.0, gam1 = 0.0, eta1 = 0.0; |
---|
| 148 | G4double plab2 = 0.0, gam2 = 0.0, eta2 = 0.0; |
---|
| 149 | |
---|
| 150 | G4double sitet = 0.0; |
---|
| 151 | G4double stet1 = 0.0; |
---|
| 152 | G4double stet2 = 0.0; |
---|
| 153 | |
---|
| 154 | G4int nbpevap = 0; |
---|
| 155 | G4int mempaw = 0, memiv = 0; |
---|
| 156 | |
---|
[819] | 157 | G4double e_evapo = 0.0; |
---|
[962] | 158 | G4double el = 0.0; |
---|
| 159 | G4double fmcv = 0.0; |
---|
[819] | 160 | |
---|
[962] | 161 | G4double aff1 = 0.0; |
---|
| 162 | G4double zff1 = 0.0; |
---|
| 163 | G4double eff1 = 0.0; |
---|
| 164 | G4double aff2 = 0.0; |
---|
| 165 | G4double zff2 = 0.0; |
---|
| 166 | G4double eff2 = 0.0; |
---|
[819] | 167 | |
---|
[962] | 168 | G4double v1 = 0.0, v2 = 0.0; |
---|
[819] | 169 | |
---|
| 170 | G4double t2 = 0.0; |
---|
[962] | 171 | G4double ctet1 = 0.0; |
---|
[819] | 172 | G4double ctet2 = 0.0; |
---|
[962] | 173 | G4double phi1 = 0.0; |
---|
[819] | 174 | G4double phi2 = 0.0; |
---|
| 175 | G4double p2 = 0.0; |
---|
[962] | 176 | G4double epf2_out = 0.0 ; |
---|
[819] | 177 | G4int lma_pf1 = 0, lmi_pf1 = 0; |
---|
| 178 | G4int lma_pf2 = 0, lmi_pf2 = 0; |
---|
| 179 | G4int nopart = 0; |
---|
| 180 | |
---|
[962] | 181 | G4double cst = 0.0, sst = 0.0, csf = 0.0, ssf = 0.0; |
---|
[819] | 182 | |
---|
[962] | 183 | G4double zf = 0.0, af = 0.0, mtota = 0.0, pleva = 0.0, pxeva = 0.0, pyeva = 0.0; |
---|
[819] | 184 | G4int ff = 0; |
---|
| 185 | G4int inum = eventnumber; |
---|
[962] | 186 | G4int inttype = 0; |
---|
[819] | 187 | G4double esrem = excitationEnergy; |
---|
| 188 | |
---|
| 189 | G4double aprf = nucleusA; |
---|
| 190 | G4double zprf = nucleusZ; |
---|
| 191 | G4double mcorem = nucleusMass; |
---|
| 192 | G4double ee = excitationEnergy; |
---|
| 193 | G4double jprf = angularMomentum; // actually root-mean-squared |
---|
| 194 | |
---|
| 195 | G4double erecrem = recoilEnergy; |
---|
[962] | 196 | G4double trem = 0.0; |
---|
[819] | 197 | G4double pxrem = momX; |
---|
| 198 | G4double pyrem = momY; |
---|
| 199 | G4double pzrem = momZ; |
---|
| 200 | |
---|
[962] | 201 | G4double remmass = 0.0; |
---|
[819] | 202 | |
---|
[962] | 203 | volant->clear(); // Clean up an initialize ABLA output. |
---|
| 204 | varntp->clear(); // Clean up an initialize ABLA output. |
---|
[819] | 205 | varntp->ntrack = 0; |
---|
[962] | 206 | volant->iv = 0; |
---|
| 207 | //volant->iv = 1; |
---|
[819] | 208 | |
---|
| 209 | G4double pcorem = std::sqrt(erecrem*(erecrem +2.*938.2796*nucleusA)); |
---|
| 210 | // G4double pcorem = std::sqrt(std::pow(momX,2) + std::pow(momY,2) + std::pow(momZ,2)); |
---|
[962] | 211 | if(pcorem != 0) { // Guard against division by zero. |
---|
| 212 | alrem = pxrem/pcorem; |
---|
| 213 | berem = pyrem/pcorem; |
---|
| 214 | garem = pzrem/pcorem; |
---|
| 215 | } else { |
---|
| 216 | alrem = 0.0; |
---|
| 217 | berem = 0.0; |
---|
| 218 | garem = 0.0; |
---|
| 219 | } |
---|
| 220 | |
---|
| 221 | G4int idebug = 0; |
---|
| 222 | G4int itest = 0; |
---|
| 223 | if(idebug == 1) { |
---|
| 224 | zprf = 81.; |
---|
| 225 | aprf = 201.; |
---|
| 226 | // ee = 86.5877686; |
---|
| 227 | ee = 300.0; |
---|
| 228 | jprf = 32.; |
---|
| 229 | zf = 0.; |
---|
| 230 | af = 0.; |
---|
| 231 | mtota = 0.; |
---|
| 232 | pleva = 0.; |
---|
| 233 | pxeva = 0.; |
---|
| 234 | pyeva = 0.; |
---|
| 235 | ff = -1; |
---|
| 236 | inttype = 0; |
---|
| 237 | inum = 2; |
---|
| 238 | } |
---|
| 239 | if(itest == 1) { |
---|
| 240 | G4cout <<" PK::: EVAPORA event " << inum << G4endl; |
---|
| 241 | G4cout <<" PK::: zprf = " << zprf << G4endl; |
---|
| 242 | G4cout <<" PK::: aprf = " << aprf << G4endl; |
---|
| 243 | G4cout <<" PK::: ee = " << ee << G4endl; |
---|
| 244 | G4cout <<" PK::: eeDiff = " << ee - 86.5877686 << G4endl; |
---|
| 245 | G4cout <<" PK::: jprf = " << jprf << G4endl; |
---|
| 246 | G4cout <<" PK::: zf = " << zf << G4endl; |
---|
| 247 | G4cout <<" PK::: af = " << af << G4endl; |
---|
| 248 | G4cout <<" PK::: mtota = " << mtota << G4endl; |
---|
| 249 | G4cout <<" PK::: pleva = " << pleva << G4endl; |
---|
| 250 | G4cout <<" PK::: pxeva = " << pxeva << G4endl; |
---|
| 251 | G4cout <<" PK::: pyeva = " << pyeva << G4endl; |
---|
| 252 | G4cout <<" PK::: ff = " << ff << G4endl; |
---|
| 253 | G4cout <<" PK::: inttype = " << inttype << G4endl; |
---|
| 254 | G4cout <<" PK::: inum = " << inum << G4endl; |
---|
| 255 | } |
---|
| 256 | |
---|
| 257 | if(esrem >= 1.0e-3) { |
---|
[819] | 258 | evapora(zprf,aprf,ee,jprf, &zf, &af, &mtota, &pleva, &pxeva, &pyeva, &ff, &inttype, &inum); |
---|
| 259 | } |
---|
| 260 | else { |
---|
| 261 | ff = 0; |
---|
| 262 | zf = zprf; |
---|
| 263 | af = aprf; |
---|
| 264 | pxeva = pxrem; |
---|
| 265 | pyeva = pyrem; |
---|
| 266 | pleva = pzrem; |
---|
[962] | 267 | } |
---|
[819] | 268 | |
---|
[962] | 269 | if (ff == 1) { |
---|
| 270 | // Fission: |
---|
| 271 | // variable ee: Energy of fissioning nucleus above the fission barrier. |
---|
| 272 | // Calcul des impulsions des particules evaporees (avant fission) |
---|
| 273 | // dans le systeme labo. |
---|
[819] | 274 | |
---|
[962] | 275 | if(verboseLevel > 2) |
---|
| 276 | G4cout << __FILE__ << ":" << __LINE__ << " Entering fission code." << G4endl; |
---|
[819] | 277 | trem = double(erecrem); |
---|
| 278 | remmass = pace2(aprf,zprf) + aprf*uma - zprf*melec; // canonic |
---|
[962] | 279 | // remmass = mcorem + double(esrem); // ok |
---|
| 280 | // remmass = mcorem; //cugnon |
---|
[819] | 281 | varntp->kfis = 1; |
---|
| 282 | G4double gamrem = (remmass + trem)/remmass; |
---|
| 283 | G4double etrem = std::sqrt(trem*(trem + 2.0*remmass))/remmass; |
---|
[962] | 284 | |
---|
[819] | 285 | // This is not treated as accurately as for the non fission case for which |
---|
| 286 | // the remnant mass is computed to satisfy the energy conservation |
---|
| 287 | // of evaporated particles. But it is not bad and more canonical! |
---|
| 288 | remmass = pace2(aprf,zprf) + aprf*uma - zprf*melec+double(esrem); // !canonic |
---|
| 289 | // Essais avec la masse de KHS (9/2002): |
---|
| 290 | el = 0.0; |
---|
| 291 | mglms(aprf,zprf,0,&el); |
---|
| 292 | remmass = zprf*fmp + (aprf-zprf)*fmn + el + double(esrem); |
---|
| 293 | gamrem = std::sqrt(std::pow(pcorem,2) + std::pow(remmass,2))/remmass; |
---|
| 294 | etrem = pcorem/remmass; |
---|
| 295 | |
---|
| 296 | csrem[0] = 0.0; // Should not be used. |
---|
| 297 | csrem[1] = alrem; |
---|
| 298 | csrem[2] = berem; |
---|
| 299 | csrem[3] = garem; |
---|
[962] | 300 | if(verboseLevel > 1) { |
---|
| 301 | G4cout << __FILE__ << ":" << __LINE__ << " momRem = (" << pxrem << ", " << pyrem << ", " << pzrem << ")" << G4endl; |
---|
| 302 | G4cout << __FILE__ << ":" << __LINE__ << " csrem = (" << csrem[1] << ", " << csrem[2] << ", " << csrem[3] << ")" << G4endl; |
---|
| 303 | } |
---|
| 304 | |
---|
[819] | 305 | // C Pour Vérif Remnant = evapo(Pre fission) + Noyau_fissionant (systÚme Remnant) |
---|
| 306 | G4double bil_e = 0.0; |
---|
| 307 | G4double bil_px = 0.0; |
---|
| 308 | G4double bil_py = 0.0; |
---|
| 309 | G4double bil_pz = 0.0; |
---|
| 310 | G4double masse = 0.0; |
---|
| 311 | |
---|
| 312 | for(G4int iloc = 1; iloc <= volant->iv; iloc++) { //DO iloc=1,iv |
---|
| 313 | mglms(double(volant->acv[iloc]),double(volant->zpcv[iloc]),0,&el); |
---|
| 314 | masse = volant->zpcv[iloc]*fmp + (volant->acv[iloc] - volant->zpcv[iloc])*fmn + el; |
---|
| 315 | bil_e = bil_e + std::sqrt(std::pow(volant->pcv[iloc],2) + std::pow(masse,2)); |
---|
| 316 | bil_px = bil_px + volant->pcv[iloc]*(volant->xcv[iloc]); |
---|
| 317 | bil_py = bil_py + volant->pcv[iloc]*(volant->ycv[iloc]); |
---|
| 318 | bil_pz = bil_pz + volant->pcv[iloc]*(volant->zcv[iloc]); |
---|
| 319 | } // enddo |
---|
| 320 | // C Ce bilan (impulsion nulle) est parfait. (Bil_Px=Bil_Px+PXEVA....) |
---|
| 321 | |
---|
| 322 | G4int ndec = 1; |
---|
| 323 | |
---|
| 324 | if(volant->iv != 0) { //then |
---|
| 325 | if(verboseLevel > 2) { |
---|
| 326 | G4cout <<"varntp->ntrack = " << varntp->ntrack << G4endl; |
---|
| 327 | G4cout <<"1st Translab: Adding indices from " << ndec << " to " << volant->iv << G4endl; |
---|
| 328 | } |
---|
| 329 | nopart = varntp->ntrack - 1; |
---|
| 330 | translab(gamrem,etrem,csrem,nopart,ndec); |
---|
| 331 | if(verboseLevel > 2) { |
---|
| 332 | G4cout <<"Translab complete!" << G4endl; |
---|
| 333 | G4cout <<"varntp->ntrack = " << varntp->ntrack << G4endl; |
---|
| 334 | } |
---|
| 335 | } |
---|
| 336 | nbpevap = volant->iv; // nombre de particules d'evaporation traitees |
---|
| 337 | |
---|
| 338 | // C |
---|
| 339 | // C Now calculation of the fission fragment distribution including |
---|
| 340 | // C evaporation from the fragments. |
---|
| 341 | // C |
---|
| 342 | |
---|
| 343 | // C Distribution of the fission fragments: |
---|
| 344 | |
---|
| 345 | // void fissionDistri(G4double a,G4double z,G4double e, |
---|
| 346 | // G4double &a1,G4double &z1,G4double &e1,G4double &v1, |
---|
| 347 | // G4double &a2,G4double &z2,G4double &e2,G4double &v2); |
---|
| 348 | |
---|
[962] | 349 | //fissionModel->fissionDistri(af,zf,ee,aff1,zff1,eff1,v1,aff2,zff2,eff2,v2); |
---|
| 350 | fissionModel->doFission(af,zf,ee,aff1,zff1,eff1,v1,aff2,zff2,eff2,v2); |
---|
[819] | 351 | // C verif des A et Z decimaux: |
---|
| 352 | G4int na_f = int(std::floor(af + 0.5)); |
---|
| 353 | G4int nz_f = int(std::floor(zf + 0.5)); |
---|
| 354 | varntp->izfis = nz_f; // copie dans le ntuple |
---|
| 355 | varntp->iafis = na_f; |
---|
| 356 | G4int na_pf1 = int(std::floor(aff1 + 0.5)); |
---|
| 357 | G4int nz_pf1 = int(std::floor(zff1 + 0.5)); |
---|
| 358 | G4int na_pf2 = int(std::floor(aff2 + 0.5)); |
---|
| 359 | G4int nz_pf2 = int(std::floor(zff2 + 0.5)); |
---|
| 360 | |
---|
| 361 | if((na_f != (na_pf1+na_pf2)) || (nz_f != (nz_pf1+nz_pf2))) { |
---|
| 362 | if(verboseLevel > 2) { |
---|
| 363 | G4cout <<"problemes arrondis dans la fission " << G4endl; |
---|
| 364 | G4cout << "af,zf,aff1,zff1,aff2,zff2" << G4endl; |
---|
| 365 | G4cout << af <<" , " << zf <<" , " << aff1 <<" , " << zff1 <<" , " << aff2 <<" , " << zff2 << G4endl; |
---|
| 366 | G4cout << "a,z,a1,z1,a2,z2 integer" << G4endl; |
---|
| 367 | G4cout << na_f <<" , " << nz_f <<" , " << na_pf1 <<" , " << nz_pf1 <<" , " << na_pf2 <<" , " << nz_pf2 << G4endl; |
---|
| 368 | } |
---|
| 369 | } |
---|
| 370 | |
---|
| 371 | // Calcul de l'impulsion des PF dans le syteme noyau de fission: |
---|
| 372 | G4int kboud = idnint(zf); |
---|
| 373 | G4int jboud = idnint(af-zf); |
---|
| 374 | //G4double ef = fb->efa[kboud][jboud]; // barriere de fission |
---|
| 375 | G4double ef = fb->efa[jboud][kboud]; // barriere de fission |
---|
| 376 | varntp->estfis = ee + ef; // copie dans le ntuple |
---|
| 377 | |
---|
| 378 | // C MASSEF = pace2(AF,ZF) |
---|
| 379 | // C MASSEF = MASSEF + AF*UMA - ZF*MELEC + EE + EF |
---|
| 380 | // C MASSE1 = pace2(DBLE(AFF1),DBLE(ZFF1)) |
---|
| 381 | // C MASSE1 = MASSE1 + AFF1*UMA - ZFF1*MELEC + EFF1 |
---|
| 382 | // C MASSE2 = pace2(DBLE(AFF2),DBLE(ZFF2)) |
---|
| 383 | // C MASSE2 = MASSE2 + AFF2*UMA - ZFF2*MELEC + EFF2 |
---|
| 384 | // C WRITE(6,*) 'MASSEF,MASSE1,MASSE2',MASSEF,MASSE1,MASSE2 |
---|
| 385 | // C MGLMS est la fonction de masse cohérente avec KHS evapo-fis. |
---|
| 386 | // C Attention aux parametres, ici 0=OPTSHP, NO microscopic correct. |
---|
| 387 | mglms(af,zf,0,&el); |
---|
| 388 | G4double massef = zf*fmp + (af - zf)*fmn + el + ee + ef; |
---|
| 389 | mglms(double(aff1),double(zff1),0,&el); |
---|
| 390 | G4double masse1 = zff1*fmp + (aff1-zff1)*fmn + el + eff1; |
---|
| 391 | mglms(aff2,zff2,0,&el); |
---|
| 392 | G4double masse2 = zff2*fmp + (aff2 - zff2)*fmn + el + eff2; |
---|
| 393 | // C WRITE(6,*) 'MASSEF,MASSE1,MASSE2',MASSEF,MASSE1,MASSE2 |
---|
| 394 | G4double b = massef - masse1 - masse2; |
---|
| 395 | if(b < 0.0) { //then |
---|
| 396 | b=0.0; |
---|
| 397 | if(verboseLevel > 2) { |
---|
| 398 | G4cout <<"anomalie dans la fission: " << G4endl; |
---|
| 399 | G4cout << inum<< " , " << af<< " , " <<zf<< " , " <<massef<< " , " <<aff1<< " , " <<zff1<< " , " <<masse1<< " , " <<aff2<< " , " <<zff2<< " , " << masse2 << G4endl; |
---|
| 400 | } |
---|
| 401 | } //endif |
---|
| 402 | G4double t1 = b*(b + 2.0*masse2)/(2.0*massef); |
---|
| 403 | G4double p1 = std::sqrt(t1*(t1 + 2.0*masse1)); |
---|
| 404 | |
---|
| 405 | G4double rndm; |
---|
[962] | 406 | rndm = randomGenerator->getRandom(); |
---|
| 407 | // standardRandom(&rndm, &(hazard->igraine[13])); |
---|
[819] | 408 | ctet1 = 2.0*rndm - 1.0; |
---|
[962] | 409 | // standardRandom(&rndm,&(hazard->igraine[9])); |
---|
| 410 | rndm = randomGenerator->getRandom(); |
---|
[819] | 411 | phi1 = rndm*2.0*3.141592654; |
---|
| 412 | |
---|
| 413 | // C ----Coefs de la transformation de Lorentz (noyau de fission -> Remnant) |
---|
| 414 | G4double peva = std::pow(pxeva,2) + std::pow(pyeva,2) + std::pow(pleva,2); |
---|
| 415 | G4double gamfis = std::sqrt(std::pow(massef,2) + peva)/massef; |
---|
| 416 | peva = std::sqrt(peva); |
---|
| 417 | G4double etfis = peva/massef; |
---|
| 418 | |
---|
[962] | 419 | G4double epf1_in = 0.0; |
---|
| 420 | G4double epf1_out = 0.0; |
---|
[819] | 421 | |
---|
| 422 | // C ----Matrice de rotation (noyau de fission -> Remnant) |
---|
| 423 | if(peva >= 1.0e-4) { |
---|
| 424 | sitet = std::sqrt(std::pow(pxeva,2)+std::pow(pyeva,2))/peva; |
---|
| 425 | } |
---|
| 426 | if(sitet > 1.0e-5) { //then |
---|
| 427 | G4double cstet = pleva/peva; |
---|
| 428 | G4double siphi = pyeva/(sitet*peva); |
---|
| 429 | G4double csphi = pxeva/(sitet*peva); |
---|
| 430 | |
---|
| 431 | R[1][1] = cstet*csphi; |
---|
| 432 | R[1][2] = -siphi; |
---|
| 433 | R[1][3] = sitet*csphi; |
---|
| 434 | R[2][1] = cstet*siphi; |
---|
| 435 | R[2][2] = csphi; |
---|
| 436 | R[2][3] = sitet*siphi; |
---|
| 437 | R[3][1] = -sitet; |
---|
| 438 | R[3][2] = 0.0; |
---|
| 439 | R[3][3] = cstet; |
---|
| 440 | } |
---|
| 441 | else { |
---|
| 442 | R[1][1] = 1.0; |
---|
| 443 | R[1][2] = 0.0; |
---|
| 444 | R[1][3] = 0.0; |
---|
| 445 | R[2][1] = 0.0; |
---|
| 446 | R[2][2] = 1.0; |
---|
| 447 | R[2][3] = 0.0; |
---|
| 448 | R[3][1] = 0.0; |
---|
| 449 | R[3][2] = 0.0; |
---|
| 450 | R[3][3] = 1.0; |
---|
| 451 | } // endif |
---|
| 452 | // c test de verif: |
---|
| 453 | |
---|
| 454 | if((zff1 <= 0.0) || (aff1 <= 0.0) || (aff1 < zff1)) { //then |
---|
| 455 | if(verboseLevel > 2) { |
---|
| 456 | G4cout <<"zf = " << zf <<" af = " << af <<"ee = " << ee <<"zff1 = " << zff1 <<"aff1 = " << aff1 << G4endl; |
---|
| 457 | } |
---|
| 458 | } |
---|
| 459 | else { |
---|
| 460 | // C ---------------------- PF1 will evaporate |
---|
| 461 | epf1_in = double(eff1); |
---|
| 462 | epf1_out = epf1_in; |
---|
| 463 | // void evapora(G4double zprf, G4double aprf, G4double ee, G4double jprf, |
---|
| 464 | // G4double *zf_par, G4double *af_par, G4double *mtota_par, |
---|
| 465 | // G4double *pleva_par, G4double *pxeva_par, G4double *pyeva_par, |
---|
| 466 | // G4double *ff_par, G4int *inttype_par, G4int *inum_par); |
---|
| 467 | G4double zf1, af1, malpha1, ffpleva1, ffpxeva1, ffpyeva1; |
---|
| 468 | G4int ff1, ftype1; |
---|
| 469 | evapora(zff1, aff1, epf1_out, 0.0, &zf1, &af1, &malpha1, &ffpleva1, |
---|
| 470 | &ffpxeva1, &ffpyeva1, &ff1, &ftype1, &inum); |
---|
| 471 | // C On ajoute le fragment: |
---|
| 472 | volant->iv = volant->iv + 1; |
---|
| 473 | volant->acv[volant->iv] = af1; |
---|
| 474 | volant->zpcv[volant->iv] = zf1; |
---|
| 475 | if(verboseLevel > 2) { |
---|
[962] | 476 | G4cout << __FILE__ << ":" << __LINE__ << " Added: zf1 = " << zf1 << " af1 = " << af1 << " at index " << volant->iv << G4endl; |
---|
| 477 | volant->dump(); |
---|
| 478 | } |
---|
| 479 | if(verboseLevel > 2) { |
---|
[819] | 480 | G4cout <<"Added fission fragment: a = " << volant->acv[volant->iv] << " z = " << volant->zpcv[volant->iv] << G4endl; |
---|
| 481 | } |
---|
| 482 | peva = std::sqrt(std::pow(ffpxeva1,2) + std::pow(ffpyeva1,2) + std::pow(ffpleva1,2)); |
---|
| 483 | volant->pcv[volant->iv] = peva; |
---|
| 484 | if(peva > 0.001) { // then |
---|
| 485 | volant->xcv[volant->iv] = ffpxeva1/peva; |
---|
| 486 | volant->ycv[volant->iv] = ffpyeva1/peva; |
---|
| 487 | volant->zcv[volant->iv] = ffpleva1/peva; |
---|
| 488 | } |
---|
| 489 | else { |
---|
| 490 | volant->xcv[volant->iv] = 1.0; |
---|
| 491 | volant->ycv[volant->iv] = 0.0; |
---|
| 492 | volant->zcv[volant->iv] = 0.0; |
---|
| 493 | } // end if |
---|
| 494 | |
---|
| 495 | // C Pour Vérif evapo de PF1 dans le systeme du Noyau Fissionant |
---|
| 496 | G4double bil1_e = 0.0; |
---|
| 497 | G4double bil1_px = 0.0; |
---|
| 498 | G4double bil1_py=0.0; |
---|
| 499 | G4double bil1_pz=0.0; |
---|
| 500 | for(G4int iloc = nbpevap + 1; iloc <= volant->iv; iloc++) { //do iloc=nbpevap+1,iv |
---|
| 501 | // for(G4int iloc = nbpevap + 1; iloc <= volant->iv + 1; iloc++) { //do iloc=nbpevap+1,iv |
---|
| 502 | mglms(volant->acv[iloc], volant->zpcv[iloc],0,&el); |
---|
| 503 | masse = volant->zpcv[iloc]*fmp + (volant->acv[iloc] - volant->zpcv[iloc])*fmn + el; |
---|
| 504 | bil1_e = bil1_e + std::sqrt(std::pow(volant->pcv[iloc],2) + std::pow(masse,2)); |
---|
| 505 | bil1_px = bil1_px + volant->pcv[iloc]*(volant->xcv[iloc]); |
---|
| 506 | bil1_py = bil1_py + volant->pcv[iloc]*(volant->ycv[iloc]); |
---|
| 507 | bil1_pz = bil1_pz + volant->pcv[iloc]*(volant->zcv[iloc]); |
---|
| 508 | } // enddo |
---|
| 509 | |
---|
| 510 | // Calcul des cosinus directeurs de PF1 dans le Remnant et calcul |
---|
| 511 | // des coefs pour la transformation de Lorentz Systeme PF --> Systeme Remnant |
---|
| 512 | translabpf(masse1,t1,p1,ctet1,phi1,gamfis,etfis,R,&plab1,&gam1,&eta1,csdir1); |
---|
| 513 | |
---|
| 514 | // calcul des impulsions des particules evaporees dans le systeme Remnant: |
---|
[962] | 515 | if(verboseLevel > 2) |
---|
[819] | 516 | G4cout <<"2nd Translab (pf1 evap): Adding indices from " << nbpevap+1 << " to " << volant->iv << G4endl; |
---|
| 517 | nopart = varntp->ntrack - 1; |
---|
| 518 | translab(gam1,eta1,csdir1,nopart,nbpevap+1); |
---|
| 519 | if(verboseLevel > 2) { |
---|
| 520 | G4cout <<"After translab call... varntp->ntrack = " << varntp->ntrack << G4endl; |
---|
| 521 | } |
---|
| 522 | memiv = nbpevap + 1; // memoires pour la future transformation |
---|
| 523 | mempaw = nopart; // remnant->labo pour pf1 et pf2. |
---|
| 524 | lmi_pf1 = nopart + nbpevap + 1; // indices min et max dans /var_ntp/ |
---|
| 525 | lma_pf1 = nopart + volant->iv; // des particules issues de pf1 |
---|
| 526 | nbpevap = volant->iv; // nombre de particules d'evaporation traitees |
---|
| 527 | } // end if |
---|
| 528 | // C --------------------- End of PF1 calculation |
---|
| 529 | |
---|
| 530 | // c test de verif: |
---|
| 531 | if((zff2 <= 0.0) || (aff2 <= 0.0) || (aff2 <= zff2)) { //then |
---|
| 532 | if(verboseLevel > 2) { |
---|
| 533 | G4cout << zf << " " << af << " " << ee << " " << zff2 << " " << aff2 << G4endl; |
---|
| 534 | } |
---|
| 535 | } |
---|
| 536 | else { |
---|
| 537 | // C ---------------------- PF2 will evaporate |
---|
| 538 | G4double epf2_in = double(eff2); |
---|
| 539 | G4double epf2_out = epf2_in; |
---|
| 540 | // void evapora(G4double zprf, G4double aprf, G4double ee, G4double jprf, |
---|
| 541 | // G4double *zf_par, G4double *af_par, G4double *mtota_par, |
---|
| 542 | // G4double *pleva_par, G4double *pxeva_par, G4double *pyeva_par, |
---|
| 543 | // G4double *ff_par, G4int *inttype_par, G4int *inum_par); |
---|
| 544 | G4double zf2, af2, malpha2, ffpleva2, ffpxeva2, ffpyeva2; |
---|
| 545 | G4int ff2, ftype2; |
---|
| 546 | evapora(zff2,aff2,epf2_out,0.0,&zf2,&af2,&malpha2,&ffpleva2, |
---|
| 547 | &ffpxeva2,&ffpyeva2,&ff2,&ftype2,&inum); |
---|
| 548 | // C On ajoute le fragment: |
---|
| 549 | volant->iv = volant->iv + 1; |
---|
| 550 | volant->acv[volant->iv] = af2; |
---|
| 551 | volant->zpcv[volant->iv] = zf2; |
---|
[962] | 552 | if(verboseLevel > 2) |
---|
| 553 | G4cout << __FILE__ << ":" << __LINE__ << " Added: zf2 = " << zf2 << " af2 = " << af2 << " at index " << volant->iv << G4endl; |
---|
[819] | 554 | if(verboseLevel > 2) { |
---|
| 555 | G4cout <<"Added fission fragment: a = " << volant->acv[volant->iv] << " z = " << volant->zpcv[volant->iv] << G4endl; |
---|
| 556 | } |
---|
| 557 | peva = std::sqrt(std::pow(ffpxeva2,2) + std::pow(ffpyeva2,2) + std::pow(ffpleva2,2)); |
---|
| 558 | volant->pcv[volant->iv] = peva; |
---|
| 559 | // exit(0); |
---|
| 560 | if(peva > 0.001) { //then |
---|
| 561 | volant->xcv[volant->iv] = ffpxeva2/peva; |
---|
| 562 | volant->ycv[volant->iv] = ffpyeva2/peva; |
---|
| 563 | volant->zcv[volant->iv] = ffpleva2/peva; |
---|
| 564 | } |
---|
| 565 | else { |
---|
| 566 | volant->xcv[volant->iv] = 1.0; |
---|
| 567 | volant->ycv[volant->iv] = 0.0; |
---|
| 568 | volant->zcv[volant->iv] = 0.0; |
---|
| 569 | } //end if |
---|
| 570 | // C Pour Vérif evapo de PF1 dans le systeme du Noyau Fissionant |
---|
| 571 | G4double bil2_e = 0.0; |
---|
| 572 | G4double bil2_px = 0.0; |
---|
| 573 | G4double bil2_py = 0.0; |
---|
| 574 | G4double bil2_pz = 0.0; |
---|
| 575 | // for(G4int iloc = nbpevap + 1; iloc <= volant->iv; iloc++) { //do iloc=nbpevap+1,iv |
---|
| 576 | for(G4int iloc = nbpevap + 1; iloc <= volant->iv; iloc++) { //do iloc=nbpevap+1,iv |
---|
| 577 | mglms(volant->acv[iloc],volant->zpcv[iloc],0,&el); |
---|
| 578 | masse = volant->zpcv[iloc]*fmp + (volant->acv[iloc] - volant->zpcv[iloc])*fmn + el; |
---|
| 579 | bil2_e = bil2_e + std::sqrt(std::pow(volant->pcv[iloc],2) + std::pow(masse,2)); |
---|
| 580 | bil2_px = bil2_px + volant->pcv[iloc]*(volant->xcv[iloc]); |
---|
| 581 | bil2_py = bil2_py + volant->pcv[iloc]*(volant->ycv[iloc]); |
---|
| 582 | bil2_pz = bil2_pz + volant->pcv[iloc]*(volant->zcv[iloc]); |
---|
| 583 | } //enddo |
---|
| 584 | |
---|
| 585 | // C ----Calcul des cosinus directeurs de PF2 dans le Remnant et calcul |
---|
| 586 | // c des coefs pour la transformation de Lorentz Systeme PF --> Systeme Remnant |
---|
| 587 | G4double t2 = b - t1; |
---|
| 588 | // G4double ctet2 = -ctet1; |
---|
| 589 | ctet2 = -1.0*ctet1; |
---|
| 590 | phi2 = dmod(phi1+3.141592654,6.283185308); |
---|
| 591 | G4double p2 = std::sqrt(t2*(t2+2.0*masse2)); |
---|
| 592 | |
---|
| 593 | // void translabpf(G4double masse1, G4double t1, G4double p1, G4double ctet1, |
---|
| 594 | // G4double phi1, G4double gamrem, G4double etrem, G4double R[][4], |
---|
| 595 | // G4double *plab1, G4double *gam1, G4double *eta1, G4double csdir[]); |
---|
| 596 | translabpf(masse2,t2,p2,ctet2,phi2,gamfis,etfis,R,&plab2,&gam2,&eta2,csdir2); |
---|
| 597 | // C |
---|
| 598 | // C calcul des impulsions des particules evaporees dans le systeme Remnant: |
---|
| 599 | // c |
---|
[962] | 600 | if(verboseLevel > 2) |
---|
| 601 | G4cout <<"3rd Translab (pf2 evap): Adding indices from " << nbpevap+1 << " to " << volant->iv << G4endl; |
---|
[819] | 602 | nopart = varntp->ntrack - 1; |
---|
| 603 | translab(gam2,eta2,csdir2,nopart,nbpevap+1); |
---|
| 604 | lmi_pf2 = nopart + nbpevap + 1; // indices min et max dans /var_ntp/ |
---|
| 605 | lma_pf2 = nopart + volant->iv; // des particules issues de pf2 |
---|
| 606 | } // end if |
---|
| 607 | // C --------------------- End of PF2 calculation |
---|
| 608 | |
---|
| 609 | // C Pour vérifications: calculs du noyau fissionant et des PF dans |
---|
| 610 | // C le systeme du remnant. |
---|
| 611 | for(G4int iloc = 1; iloc <= 3; iloc++) { // do iloc=1,3 |
---|
| 612 | pfis_rem[iloc] = 0.0; |
---|
| 613 | } // enddo |
---|
| 614 | G4double efis_rem, pfis_trav[4]; |
---|
| 615 | lorab(gamfis,etfis,massef,pfis_rem,&efis_rem,pfis_trav); |
---|
| 616 | rotab(R,pfis_trav,pfis_rem); |
---|
| 617 | |
---|
| 618 | stet1 = std::sqrt(1.0 - std::pow(ctet1,2)); |
---|
| 619 | pf1_rem[1] = p1*stet1*std::cos(phi1); |
---|
| 620 | pf1_rem[2] = p1*stet1*std::sin(phi1); |
---|
| 621 | pf1_rem[3] = p1*ctet1; |
---|
| 622 | G4double e1_rem; |
---|
| 623 | lorab(gamfis,etfis,masse1+t1,pf1_rem,&e1_rem,pfis_trav); |
---|
| 624 | rotab(R,pfis_trav,pf1_rem); |
---|
| 625 | |
---|
| 626 | stet2 = std::sqrt(1.0 - std::pow(ctet2,2)); |
---|
| 627 | |
---|
| 628 | G4double pf2_rem[4]; |
---|
| 629 | G4double e2_rem; |
---|
| 630 | pf2_rem[1] = p2*stet2*std::cos(phi2); |
---|
| 631 | pf2_rem[2] = p2*stet2*std::sin(phi2); |
---|
| 632 | pf2_rem[3] = p2*ctet2; |
---|
| 633 | lorab(gamfis,etfis,masse2+t2,pf2_rem,&e2_rem,pfis_trav); |
---|
| 634 | rotab(R,pfis_trav,pf2_rem); |
---|
| 635 | // C Verif 0: Remnant = evapo_pre_fission + Noyau Fissionant |
---|
| 636 | bil_e = remmass - efis_rem - bil_e; |
---|
| 637 | bil_px = bil_px + pfis_rem[1]; |
---|
| 638 | bil_py = bil_py + pfis_rem[2]; |
---|
| 639 | bil_pz = bil_pz + pfis_rem[3]; |
---|
| 640 | // C Verif 1: noyau fissionant = PF1 + PF2 dans le systeme remnant |
---|
| 641 | // G4double bilan_e = efis_rem - e1_rem - e2_rem; |
---|
| 642 | // G4double bilan_px = pfis_rem[1] - pf1_rem[1] - pf2_rem[1]; |
---|
| 643 | // G4double bilan_py = pfis_rem[2] - pf1_rem[2] - pf2_rem[2]; |
---|
| 644 | // G4double bilan_pz = pfis_rem[3] - pf1_rem[3] - pf2_rem[3]; |
---|
| 645 | // C Verif 2: PF1 et PF2 egaux a toutes leurs particules evaporees |
---|
| 646 | // C (Systeme remnant) |
---|
| 647 | if((lma_pf1-lmi_pf1) != 0) { //then |
---|
| 648 | G4double bil_e_pf1 = e1_rem - epf1_out; |
---|
| 649 | G4double bil_px_pf1 = pf1_rem[1]; |
---|
| 650 | G4double bil_py_pf1 = pf1_rem[2]; |
---|
| 651 | G4double bil_pz_pf1 = pf1_rem[3]; |
---|
| 652 | for(G4int ipf1 = lmi_pf1; ipf1 <= lma_pf1; ipf1++) { //do ipf1=lmi_pf1,lma_pf1 |
---|
[1196] | 653 | if(varntp->enerj[ipf1] <= 0.0) continue; // Safeguard against a division by zero |
---|
[819] | 654 | bil_e_pf1 = bil_e_pf1 - (std::pow(varntp->plab[ipf1],2) + std::pow(varntp->enerj[ipf1],2))/(2.0*(varntp->enerj[ipf1])); |
---|
| 655 | cst = std::cos(varntp->tetlab[ipf1]/57.2957795); |
---|
| 656 | sst = std::sin(varntp->tetlab[ipf1]/57.2957795); |
---|
| 657 | csf = std::cos(varntp->philab[ipf1]/57.2957795); |
---|
| 658 | ssf = std::sin(varntp->philab[ipf1]/57.2957795); |
---|
| 659 | bil_px_pf1 = bil_px_pf1 - varntp->plab[ipf1]*sst*csf; |
---|
| 660 | bil_py_pf1 = bil_py_pf1 - varntp->plab[ipf1]*sst*ssf; |
---|
| 661 | bil_pz_pf1 = bil_pz_pf1 - varntp->plab[ipf1]*cst; |
---|
| 662 | } // enddo |
---|
| 663 | } //endif |
---|
| 664 | |
---|
| 665 | if((lma_pf2-lmi_pf2) != 0) { //then |
---|
| 666 | G4double bil_e_pf2 = e2_rem - epf2_out; |
---|
| 667 | G4double bil_px_pf2 = pf2_rem[1]; |
---|
| 668 | G4double bil_py_pf2 = pf2_rem[2]; |
---|
| 669 | G4double bil_pz_pf2 = pf2_rem[3]; |
---|
| 670 | for(G4int ipf2 = lmi_pf2; ipf2 <= lma_pf2; ipf2++) { //do ipf2=lmi_pf2,lma_pf2 |
---|
[1196] | 671 | if(varntp->enerj[ipf2] <= 0.0) continue; // Safeguard against a division by zero |
---|
[819] | 672 | bil_e_pf2 = bil_e_pf2 - (std::pow(varntp->plab[ipf2],2) + std::pow(varntp->enerj[ipf2],2))/(2.0*(varntp->enerj[ipf2])); |
---|
| 673 | G4double cst = std::cos(varntp->tetlab[ipf2]/57.2957795); |
---|
| 674 | G4double sst = std::sin(varntp->tetlab[ipf2]/57.2957795); |
---|
| 675 | G4double csf = std::cos(varntp->philab[ipf2]/57.2957795); |
---|
| 676 | G4double ssf = std::sin(varntp->philab[ipf2]/57.2957795); |
---|
| 677 | bil_px_pf2 = bil_px_pf2 - varntp->plab[ipf2]*sst*csf; |
---|
| 678 | bil_py_pf2 = bil_py_pf2 - varntp->plab[ipf2]*sst*ssf; |
---|
| 679 | bil_pz_pf2 = bil_pz_pf2 - varntp->plab[ipf2]*cst; |
---|
| 680 | } // enddo |
---|
| 681 | } //endif |
---|
| 682 | // C |
---|
| 683 | // C ---- Transformation systeme Remnant -> systeme labo. (evapo de PF1 ET PF2) |
---|
| 684 | // C |
---|
| 685 | // G4double mempaw, memiv; |
---|
[962] | 686 | if(verboseLevel > 2) |
---|
[819] | 687 | G4cout <<"4th Translab: Adding indices from " << memiv << " to " << volant->iv << G4endl; |
---|
| 688 | translab(gamrem,etrem,csrem,mempaw,memiv); |
---|
| 689 | // C ******************* END of fission calculations ************************ |
---|
[962] | 690 | if(verboseLevel > 2) { |
---|
| 691 | G4cout <<"Dump at the end of fission event " << G4endl; |
---|
| 692 | volant->dump(); |
---|
| 693 | G4cout <<"End of dump." << G4endl; |
---|
| 694 | } |
---|
[819] | 695 | } |
---|
| 696 | else { |
---|
| 697 | // C ************************ Evapo sans fission ***************************** |
---|
| 698 | // C Here, FF=0, --> Evapo sans fission, on ajoute le fragment: |
---|
| 699 | // C ************************************************************************* |
---|
| 700 | varntp->kfis = 0; |
---|
| 701 | if(verboseLevel > 2) { |
---|
| 702 | G4cout <<"Evaporation without fission" << G4endl; |
---|
| 703 | } |
---|
| 704 | volant->iv = volant->iv + 1; |
---|
| 705 | volant->acv[volant->iv] = af; |
---|
| 706 | volant->zpcv[volant->iv] = zf; |
---|
[962] | 707 | if(verboseLevel > 2) |
---|
| 708 | G4cout << __FILE__ << ":" << __LINE__ << " Added: zf = " << zf << " af = " << af << " at index " << volant->iv << G4endl; |
---|
[819] | 709 | G4double peva = std::sqrt(std::pow(pxeva,2)+std::pow(pyeva,2)+std::pow(pleva,2)); |
---|
| 710 | volant->pcv[volant->iv] = peva; |
---|
| 711 | if(peva > 0.001) { //then |
---|
| 712 | volant->xcv[volant->iv] = pxeva/peva; |
---|
| 713 | volant->ycv[volant->iv] = pyeva/peva; |
---|
| 714 | volant->zcv[volant->iv] = pleva/peva; |
---|
| 715 | } |
---|
| 716 | else { |
---|
| 717 | volant->xcv[volant->iv] = 1.0; |
---|
| 718 | volant->ycv[volant->iv] = 0.0; |
---|
| 719 | volant->zcv[volant->iv] = 0.0; |
---|
| 720 | } // end if |
---|
| 721 | |
---|
| 722 | // C |
---|
| 723 | // C calcul des impulsions des particules evaporees dans le systeme labo: |
---|
| 724 | // c |
---|
| 725 | trem = double(erecrem); |
---|
| 726 | // C REMMASS = pace2(APRF,ZPRF) + APRF*UMA - ZPRF*MELEC !Canonic |
---|
| 727 | // C REMMASS = MCOREM + DBLE(ESREM) !OK |
---|
| 728 | remmass = mcorem; //Cugnon |
---|
| 729 | // C GAMREM = (REMMASS + TREM)/REMMASS !OK |
---|
| 730 | // C ETREM = DSQRT(TREM*(TREM + 2.*REMMASS))/REMMASS !OK |
---|
| 731 | csrem[0] = 0.0; // Should not be used. |
---|
| 732 | csrem[1] = alrem; |
---|
| 733 | csrem[2] = berem; |
---|
| 734 | csrem[3] = garem; |
---|
| 735 | |
---|
| 736 | // for(G4int j = 1; j <= volant->iv; j++) { //do j=1,iv |
---|
| 737 | for(G4int j = 1; j <= volant->iv; j++) { //do j=1,iv |
---|
| 738 | if(volant->acv[j] == 0) { |
---|
| 739 | if(verboseLevel > 2) { |
---|
| 740 | G4cout <<"volant->acv[" << j << "] = 0" << G4endl; |
---|
| 741 | G4cout <<"volant->iv = " << volant->iv << G4endl; |
---|
| 742 | } |
---|
| 743 | } |
---|
| 744 | if(volant->acv[j] > 0) { |
---|
| 745 | mglms(volant->acv[j],volant->zpcv[j],0,&el); |
---|
| 746 | fmcv = volant->zpcv[j]*fmp + (volant->acv[j] - volant->zpcv[j])*fmn + el; |
---|
| 747 | e_evapo = e_evapo + std::sqrt(std::pow(volant->pcv[j],2) + std::pow(fmcv,2)); |
---|
| 748 | } |
---|
| 749 | } // enddo |
---|
| 750 | |
---|
| 751 | // C Redefinition pour conservation d'impulsion!!! |
---|
| 752 | // C this mass obtained by energy balance is very close to the |
---|
| 753 | // C mass of the remnant computed by pace2 + excitation energy (EE). (OK) |
---|
| 754 | remmass = e_evapo; |
---|
| 755 | |
---|
| 756 | G4double gamrem = std::sqrt(std::pow(pcorem,2)+std::pow(remmass,2))/remmass; |
---|
| 757 | G4double etrem = pcorem/remmass; |
---|
| 758 | |
---|
[962] | 759 | if(verboseLevel > 2) |
---|
[819] | 760 | G4cout <<"5th Translab (no fission): Adding indices from " << 1 << " to " << volant->iv << G4endl; |
---|
| 761 | nopart = varntp->ntrack - 1; |
---|
| 762 | translab(gamrem,etrem,csrem,nopart,1); |
---|
| 763 | |
---|
[962] | 764 | if(verboseLevel > 2) { |
---|
| 765 | G4cout <<"Dump at the end of evaporation event " << G4endl; |
---|
| 766 | volant->dump(); |
---|
| 767 | G4cout <<"End of dump." << G4endl; |
---|
| 768 | } |
---|
[819] | 769 | // C End of the (FISSION - NO FISSION) condition (FF=1 or 0) |
---|
| 770 | } //end if |
---|
| 771 | // C *********************** FIN de l'EVAPO KHS ******************** |
---|
| 772 | } |
---|
| 773 | |
---|
| 774 | // Evaporation code |
---|
| 775 | void G4Abla::initEvapora() |
---|
| 776 | { |
---|
| 777 | // 37 C PROJECTILE AND TARGET PARAMETERS + CROSS SECTIONS |
---|
| 778 | // 38 C COMMON /ABLAMAIN/ AP,ZP,AT,ZT,EAP,BETA,BMAXNUC,CRTOT,CRNUC, |
---|
| 779 | // 39 C R_0,R_P,R_T, IMAX,IRNDM,PI, |
---|
| 780 | // 40 C BFPRO,SNPRO,SPPRO,SHELL |
---|
| 781 | // 41 C |
---|
| 782 | // 42 C AP,ZP,AT,ZT - PROJECTILE AND TARGET MASSES |
---|
| 783 | // 43 C EAP,BETA - BEAM ENERGY PER NUCLEON, V/C |
---|
| 784 | // 44 C BMAXNUC - MAX. IMPACT PARAMETER FOR NUCL. REAC. |
---|
| 785 | // 45 C CRTOT,CRNUC - TOTAL AND NUCLEAR REACTION CROSS SECTION |
---|
| 786 | // 46 C R_0,R_P,R_T, - RADIUS PARAMETER, PROJECTILE+ TARGET RADII |
---|
| 787 | // 47 C IMAX,IRNDM,PI - MAXIMUM NUMBER OF EVENTS, DUMMY, 3.141... |
---|
| 788 | // 48 C BFPRO - FISSION BARRIER OF THE PROJECTILE |
---|
| 789 | // 49 C SNPRO - NEUTRON SEPARATION ENERGY OF THE PROJECTILE |
---|
| 790 | // 50 C SPPRO - PROTON " " " " " |
---|
| 791 | // 51 C SHELL - GROUND STATE SHELL CORRECTION |
---|
| 792 | // 52 C--------------------------------------------------------------------- |
---|
| 793 | // 53 C |
---|
| 794 | // 54 C ENERGIES WIDTHS AND CROSS SECTIONS FOR EM EXCITATION |
---|
| 795 | // 55 C COMMON /EMDPAR/ EGDR,EGQR,FWHMGDR,FWHMGQR,CREMDE1,CREMDE2, |
---|
| 796 | // 56 C AE1,BE1,CE1,AE2,BE2,CE2,SR1,SR2,XR |
---|
| 797 | // 57 C |
---|
| 798 | // 58 C EGDR,EGQR - MEAN ENERGY OF GDR AND GQR |
---|
| 799 | // 59 C FWHMGDR,FWHMGQR - FWHM OF GDR, GQR |
---|
| 800 | // 60 C CREMDE1,CREMDE2 - EM CROSS SECTION FOR E1 AND E2 |
---|
| 801 | // 61 C AE1,BE1,CE1 - ARRAYS TO CALCULATE |
---|
| 802 | // 62 C AE2,BE2,CE2 - THE EXCITATION ENERGY AFTER E.M. EXC. |
---|
| 803 | // 63 C SR1,SR2,XR - WITH MONTE CARLO |
---|
| 804 | // 64 C--------------------------------------------------------------------- |
---|
| 805 | // 65 C |
---|
| 806 | // 66 C DEFORMATIONS AND G.S. SHELL EFFECTS |
---|
| 807 | // 67 C COMMON /ECLD/ ECGNZ,ECFNZ,VGSLD,ALPHA |
---|
| 808 | // 68 C |
---|
| 809 | // 69 C ECGNZ - GROUND STATE SHELL CORR. FRLDM FOR A SPHERICAL G.S. |
---|
| 810 | // 70 C ECFNZ - SHELL CORRECTION FOR THE SADDLE POINT (NOW: == 0) |
---|
| 811 | // 71 C VGSLD - DIFFERENCE BETWEEN DEFORMED G.S. AND LDM VALUE |
---|
| 812 | // 72 C ALPHA - ALPHA GROUND STATE DEFORMATION (THIS IS NOT BETA2!) |
---|
| 813 | // 73 C BETA2 = SQRT(5/(4PI)) * ALPHA |
---|
| 814 | // 74 C--------------------------------------------------------------------- |
---|
| 815 | // 75 C |
---|
| 816 | // 76 C ARRAYS FOR EXCITATION ENERGY BY STATISTICAL HOLE ENERY MODEL |
---|
| 817 | // 77 C COMMON /EENUC/ SHE, XHE |
---|
| 818 | // 78 C |
---|
| 819 | // 79 C SHE, XHE - ARRAYS TO CALCULATE THE EXC. ENERGY AFTER |
---|
| 820 | // 80 C ABRASION BY THE STATISTICAL HOLE ENERGY MODEL |
---|
| 821 | // 81 C--------------------------------------------------------------------- |
---|
| 822 | // 82 C |
---|
| 823 | // 83 C G.S. SHELL EFFECT |
---|
| 824 | // 84 C COMMON /EC2SUB/ ECNZ |
---|
| 825 | // 85 C |
---|
| 826 | // 86 C ECNZ G.S. SHELL EFFECT FOR THE MASSES (IDENTICAL TO ECGNZ) |
---|
| 827 | // 87 C--------------------------------------------------------------------- |
---|
| 828 | // 88 C |
---|
| 829 | // 89 C OPTIONS AND PARAMETERS FOR FISSION CHANNEL |
---|
| 830 | // 90 C COMMON /FISS/ AKAP,BET,HOMEGA,KOEFF,IFIS, |
---|
| 831 | // 91 C OPTSHP,OPTXFIS,OPTLES,OPTCOL |
---|
| 832 | // 92 C |
---|
| 833 | // 93 C AKAP - HBAR**2/(2* MN * R_0**2) = 10 MEV |
---|
| 834 | // 94 C BET - REDUCED NUCLEAR FRICTION COEFFICIENT IN (10**21 S**-1) |
---|
| 835 | // 95 C HOMEGA - CURVATURE OF THE FISSION BARRIER = 1 MEV |
---|
| 836 | // 96 C KOEFF - COEFFICIENT FOR THE LD FISSION BARRIER == 1.0 |
---|
| 837 | // 97 C IFIS - 0/1 FISSION CHANNEL OFF/ON |
---|
| 838 | // 98 C OPTSHP - INTEGER SWITCH FOR SHELL CORRECTION IN MASSES/ENERGY |
---|
| 839 | // 99 C = 0 NO MICROSCOPIC CORRECTIONS IN MASSES AND ENERGY |
---|
| 840 | // 100 C = 1 SHELL , NO PAIRING |
---|
| 841 | // 101 C = 2 PAIRING, NO SHELL |
---|
| 842 | // 102 C = 3 SHELL AND PAIRING |
---|
| 843 | // 103 C OPTCOL - 0/1 COLLECTIVE ENHANCEMENT SWITCHED ON/OFF |
---|
| 844 | // 104 C OPTXFIS- 0,1,2 FOR MYERS & SWIATECKI, DAHLINGER, ANDREYEV |
---|
| 845 | // 105 C FISSILITY PARAMETER. |
---|
| 846 | // 106 C OPTLES - CONSTANT TEMPERATURE LEVEL DENSITY FOR A,Z > TH-224 |
---|
| 847 | // 107 C OPTCOL - 0/1 COLLECTIVE ENHANCEMENT OFF/ON |
---|
| 848 | // 108 C--------------------------------------------------------------------- |
---|
| 849 | // 109 C |
---|
| 850 | // 110 C OPTIONS |
---|
| 851 | // 111 C COMMON /OPT/ OPTEMD,OPTCHA,EEFAC |
---|
| 852 | // 112 C |
---|
| 853 | // 113 C OPTEMD - 0/1 NO EMD / INCL. EMD |
---|
| 854 | // 114 C OPTCHA - 0/1 0 GDR / 1 HYPERGEOMETRICAL PREFRAGMENT-CHARGE-DIST. |
---|
| 855 | // 115 C *** RECOMMENDED IS OPTCHA = 1 *** |
---|
| 856 | // 116 C EEFAC - EXCITATION ENERGY FACTOR, 2.0 RECOMMENDED |
---|
| 857 | // 117 C--------------------------------------------------------------------- |
---|
| 858 | // 118 C |
---|
| 859 | // 119 C FISSION BARRIERS |
---|
| 860 | // 120 C COMMON /FB/ EFA |
---|
| 861 | // 121 C EFA - ARRAY OF FISSION BARRIERS |
---|
| 862 | // 122 C--------------------------------------------------------------------- |
---|
| 863 | // 123 C |
---|
| 864 | // 124 C p LEVEL DENSITY PARAMETERS |
---|
| 865 | // 125 C COMMON /ALD/ AV,AS,AK,OPTAFAN |
---|
| 866 | // 126 C AV,AS,AK - VOLUME,SURFACE,CURVATURE DEPENDENCE OF THE |
---|
| 867 | // 127 C LEVEL DENSITY PARAMETER |
---|
| 868 | // 128 C OPTAFAN - 0/1 AF/AN >=1 OR AF/AN ==1 |
---|
| 869 | // 129 C RECOMMENDED IS OPTAFAN = 0 |
---|
| 870 | // 130 C--------------------------------------------------------------------- |
---|
| 871 | // 131 C ____________________________________________________________________ |
---|
| 872 | // 132 C / |
---|
| 873 | // 133 C / INITIALIZES PARAMETERS IN COMMON /ABRAMAIN/, /EMDPAR/, /ECLD/ ... |
---|
| 874 | // 134 C / PROJECTILE PARAMETERS, EMD PARAMETERS, SHELL CORRECTION TABLES. |
---|
| 875 | // 135 C / CALCULATES MAXIMUM IMPACT PARAMETER FOR NUCLEAR COLLISIONS AND |
---|
| 876 | // 136 C / TOTAL GEOMETRICAL CROSS SECTION + EMD CROSS SECTIONS |
---|
| 877 | // 137 C ____________________________________________________________________ |
---|
| 878 | // 138 C |
---|
| 879 | // 139 C |
---|
| 880 | // 201 C |
---|
| 881 | // 202 C---------- SET INPUT VALUES |
---|
| 882 | // 203 C |
---|
| 883 | // 204 C *** INPUT FROM UNIT 10 IN THE FOLLOWING SEQUENCE ! |
---|
| 884 | // 205 C AP1 = INTEGER ! |
---|
| 885 | // 206 C ZP1 = INTEGER ! |
---|
| 886 | // 207 C AT1 = INTEGER ! |
---|
| 887 | // 208 C ZT1 = INTEGER ! |
---|
| 888 | // 209 C EAP = REAL ! |
---|
| 889 | // 210 C IMAX = INTEGER ! |
---|
| 890 | // 211 C IFIS = INTEGER SWITCH FOR FISSION |
---|
| 891 | // 212 C OPTSHP = INTEGER SWITCH FOR SHELL CORRECTION IN MASSES/ENERGY |
---|
| 892 | // 213 C =0 NO MICROSCOPIC CORRECTIONS IN MASSES AND ENERGY |
---|
| 893 | // 214 C =1 SHELL , NO PAIRING CORRECTION |
---|
| 894 | // 215 C =2 PAIRING, NO SHELL CORRECTION |
---|
| 895 | // 216 C =3 SHELL AND PAIRING CORRECTION IN MASSES AND ENERGY |
---|
| 896 | // 217 C OPTEMD =0,1 0 NO EMD, 1 INCL. EMD |
---|
| 897 | // 218 C ELECTROMAGNETIC DISSOZIATION IS CALCULATED AS WELL. |
---|
| 898 | // 219 C OPTCHA =0,1 0 GDR- , 1 HYPERGEOMETRICAL PREFRAGMENT-CHARGE-DIST. |
---|
| 899 | // 220 C RECOMMENDED IS OPTCHA=1 |
---|
| 900 | // 221 C OPTCOL =0,1 COLLECTIVE ENHANCEMENT SWITCHED ON 1 OR OFF 0 IN DENSN |
---|
| 901 | // 222 C OPTAFAN=0,1 SWITCH FOR AF/AN = 1 IN DENSNIV 0 AF/AN>1 1 AF/AN=1 |
---|
| 902 | // 223 C AKAP = REAL ALWAYS EQUALS 10 |
---|
| 903 | // 224 C BET = REAL REDUCED FRICTION COEFFICIENT / 10**(+21) S**(-1) |
---|
| 904 | // 225 C HOMEGA = REAL CURVATURE / MEV RECOMMENDED = 1. MEV |
---|
| 905 | // 226 C KOEFF = REAL COEFFICIENT FOR FISSION BARRIER |
---|
| 906 | // 227 C OPTXFIS= INTEGER 0,1,2 FOR MYERS & SWIATECKI, DAHLINGER, ANDREYEV |
---|
| 907 | // 228 C FISSILITY PARAMETER. |
---|
| 908 | // 229 C EEFAC = REAL EMPIRICAL FACTOR FOR THE EXCITATION ENERGY |
---|
| 909 | // 230 C RECOMMENDED 2.D0, STATISTICAL ABRASION MODELL 1.D0 |
---|
| 910 | // 231 C AV = REAL KOEFFICIENTS FOR CALCULATION OF A(TILDE) |
---|
| 911 | // 232 C AS = REAL LEVEL DENSITY PARAMETER |
---|
| 912 | // 233 C AK = REAL |
---|
| 913 | // 234 C |
---|
| 914 | // 235 C This following inputs will be initialized in the main through the |
---|
| 915 | // 236 C common /ABLAMAIN/ (A.B.) |
---|
| 916 | // 237 |
---|
| 917 | |
---|
| 918 | // switch-fission.1=on.0=off |
---|
| 919 | fiss->ifis = 1; |
---|
| 920 | |
---|
| 921 | // shell+pairing.0-1-2-3 |
---|
| 922 | fiss->optshp = 0; |
---|
| 923 | |
---|
| 924 | // optemd =0,1 0 no emd, 1 incl. emd |
---|
| 925 | opt->optemd = 1; |
---|
| 926 | // read(10,*,iostat=io) dum(10),optcha |
---|
| 927 | opt->optcha = 1; |
---|
| 928 | |
---|
| 929 | // not.to.be.changed.(akap) |
---|
| 930 | fiss->akap = 10.0; |
---|
| 931 | |
---|
| 932 | // nuclear.viscosity.(beta) |
---|
| 933 | fiss->bet = 1.5; |
---|
| 934 | |
---|
| 935 | // potential-curvature |
---|
| 936 | fiss->homega = 1.0; |
---|
| 937 | |
---|
| 938 | // fission-barrier-coefficient |
---|
| 939 | fiss->koeff = 1.; |
---|
| 940 | |
---|
| 941 | //collective enhancement switched on 1 or off 0 in densn (qr=val or =1.) |
---|
| 942 | fiss->optcol = 0; |
---|
| 943 | |
---|
| 944 | // switch-for-low-energy-sys |
---|
| 945 | fiss->optles = 0; |
---|
| 946 | |
---|
| 947 | opt->eefac = 2.; |
---|
| 948 | |
---|
| 949 | ald->optafan = 0; |
---|
| 950 | |
---|
| 951 | ald->av = 0.073e0; |
---|
| 952 | ald->as = 0.095e0; |
---|
| 953 | ald->ak = 0.0e0; |
---|
| 954 | |
---|
| 955 | if(verboseLevel > 3) { |
---|
| 956 | G4cout <<"ifis " << fiss->ifis << G4endl; |
---|
| 957 | G4cout <<"optshp " << fiss->optshp << G4endl; |
---|
| 958 | G4cout <<"optemd " << opt->optemd << G4endl; |
---|
| 959 | G4cout <<"optcha " << opt->optcha << G4endl; |
---|
| 960 | G4cout <<"akap " << fiss->akap << G4endl; |
---|
| 961 | G4cout <<"bet " << fiss->bet << G4endl; |
---|
| 962 | G4cout <<"homega " << fiss->homega << G4endl; |
---|
| 963 | G4cout <<"koeff " << fiss->koeff << G4endl; |
---|
| 964 | G4cout <<"optcol " << fiss->optcol << G4endl; |
---|
| 965 | G4cout <<"optles " << fiss->optles << G4endl; |
---|
| 966 | G4cout <<"eefac " << opt->eefac << G4endl; |
---|
| 967 | G4cout <<"optafan " << ald->optafan << G4endl; |
---|
| 968 | G4cout <<"av " << ald->av << G4endl; |
---|
| 969 | G4cout <<"as " << ald->as << G4endl; |
---|
| 970 | G4cout <<"ak " << ald->ak << G4endl; |
---|
| 971 | } |
---|
| 972 | fiss->optxfis = 1; |
---|
| 973 | |
---|
| 974 | G4InclAblaDataFile *dataInterface = new G4InclAblaDataFile(); |
---|
| 975 | if(dataInterface->readData() == true) { |
---|
| 976 | if(verboseLevel > 0) { |
---|
| 977 | G4cout <<"G4Abla: Datafiles read successfully." << G4endl; |
---|
| 978 | } |
---|
| 979 | } |
---|
| 980 | else { |
---|
| 981 | G4Exception("ERROR: Failed to read datafiles."); |
---|
| 982 | } |
---|
| 983 | |
---|
[962] | 984 | for(int z = 0; z < 99; z++) { //do 30 z = 0,98,1 |
---|
[819] | 985 | for(int n = 0; n < 154; n++) { //do 31 n = 0,153,1 |
---|
| 986 | ecld->ecfnz[n][z] = 0.e0; |
---|
| 987 | ec2sub->ecnz[n][z] = dataInterface->getEcnz(n,z); |
---|
| 988 | ecld->ecgnz[n][z] = dataInterface->getEcnz(n,z); |
---|
| 989 | ecld->alpha[n][z] = dataInterface->getAlpha(n,z); |
---|
| 990 | ecld->vgsld[n][z] = dataInterface->getVgsld(n,z); |
---|
[962] | 991 | // if(ecld->ecgnz[n][z] != 0.0) G4cout <<"ecgnz[" << n << "][" << z << "] = " << ecld->ecgnz[n][z] << G4endl; |
---|
[819] | 992 | } |
---|
| 993 | } |
---|
| 994 | |
---|
| 995 | for(int z = 0; z < 500; z++) { |
---|
| 996 | for(int a = 0; a < 500; a++) { |
---|
| 997 | pace->dm[z][a] = dataInterface->getPace2(z,a); |
---|
| 998 | } |
---|
| 999 | } |
---|
| 1000 | |
---|
| 1001 | delete dataInterface; |
---|
| 1002 | } |
---|
| 1003 | |
---|
| 1004 | void G4Abla::qrot(G4double z, G4double a, G4double bet, G4double sig, G4double u, G4double *qr) |
---|
| 1005 | { |
---|
[962] | 1006 | G4double ucr = 10.0; // Critical energy for damping. |
---|
| 1007 | G4double dcr = 40.0; // Width of damping. |
---|
| 1008 | G4double ponq = 0.0, dn = 0.0, n = 0.0, dz = 0.0; |
---|
[819] | 1009 | |
---|
| 1010 | if(((std::fabs(bet)-1.15) < 0) || ((std::fabs(bet)-1.15) == 0)) { |
---|
| 1011 | goto qrot10; |
---|
| 1012 | } |
---|
| 1013 | |
---|
| 1014 | if((std::fabs(bet)-1.15) > 0) { |
---|
| 1015 | goto qrot11; |
---|
| 1016 | } |
---|
| 1017 | |
---|
| 1018 | qrot10: |
---|
| 1019 | n = a - z; |
---|
| 1020 | dz = std::fabs(z - 82.0); |
---|
| 1021 | if (n > 104) { |
---|
| 1022 | dn = std::fabs(n-126.e0); |
---|
| 1023 | } |
---|
| 1024 | else { |
---|
| 1025 | dn = std::fabs(n - 82.0); |
---|
| 1026 | } |
---|
| 1027 | |
---|
| 1028 | bet = 0.022 + 0.003*dn + 0.005*dz; |
---|
| 1029 | |
---|
| 1030 | sig = 25.0*std::pow(bet,2) * sig; |
---|
| 1031 | |
---|
| 1032 | qrot11: |
---|
| 1033 | ponq = (u - ucr)/dcr; |
---|
| 1034 | |
---|
| 1035 | if (ponq > 700.0) { |
---|
| 1036 | ponq = 700.0; |
---|
| 1037 | } |
---|
| 1038 | if (sig < 1.0) { |
---|
| 1039 | sig = 1.0; |
---|
| 1040 | } |
---|
| 1041 | (*qr) = 1.0/(1.0 + std::exp(ponq)) * (sig - 1.0) + 1.0; |
---|
| 1042 | |
---|
| 1043 | if ((*qr) < 1.0) { |
---|
| 1044 | (*qr) = 1.0; |
---|
| 1045 | } |
---|
| 1046 | |
---|
| 1047 | return; |
---|
| 1048 | } |
---|
| 1049 | |
---|
| 1050 | void G4Abla::mglw(G4double a, G4double z, G4double *el) |
---|
| 1051 | { |
---|
| 1052 | // MODEL DE LA GOUTTE LIQUIDE DE C. F. WEIZSACKER. |
---|
| 1053 | // USUALLY AN OBSOLETE OPTION |
---|
| 1054 | |
---|
[962] | 1055 | G4int a1 = 0, z1 = 0; |
---|
[819] | 1056 | G4double xv = 0.0, xs = 0.0, xc = 0.0, xa = 0.0; |
---|
| 1057 | |
---|
| 1058 | a1 = idnint(a); |
---|
| 1059 | z1 = idnint(z); |
---|
| 1060 | |
---|
| 1061 | if ((a <= 0.01) || (z < 0.01)) { |
---|
| 1062 | (*el) = 1.0e38; |
---|
| 1063 | } |
---|
| 1064 | else { |
---|
| 1065 | xv = -15.56*a; |
---|
| 1066 | xs = 17.23*std::pow(a,(2.0/3.0)); |
---|
| 1067 | |
---|
| 1068 | if (a > 1.0) { |
---|
| 1069 | xc = 0.7*z*(z-1.0)*std::pow((a-1.0),(-1.e0/3.e0)); |
---|
| 1070 | } |
---|
| 1071 | else { |
---|
| 1072 | xc = 0.0; |
---|
| 1073 | } |
---|
| 1074 | } |
---|
| 1075 | |
---|
| 1076 | xa = 23.6*(std::pow((a-2.0*z),2)/a); |
---|
| 1077 | (*el) = xv+xs+xc+xa; |
---|
| 1078 | return; |
---|
| 1079 | } |
---|
| 1080 | |
---|
| 1081 | void G4Abla::mglms(G4double a, G4double z, G4int refopt4, G4double *el) |
---|
| 1082 | { |
---|
| 1083 | // USING FUNCTION EFLMAC(IA,IZ,0) |
---|
| 1084 | // |
---|
| 1085 | // REFOPT4 = 0 : WITHOUT MICROSCOPIC CORRECTIONS |
---|
| 1086 | // REFOPT4 = 1 : WITH SHELL CORRECTION |
---|
| 1087 | // REFOPT4 = 2 : WITH PAIRING CORRECTION |
---|
| 1088 | // REFOPT4 = 3 : WITH SHELL- AND PAIRING CORRECTION |
---|
| 1089 | |
---|
| 1090 | // 1839 C----------------------------------------------------------------------- |
---|
| 1091 | // 1840 C A1 LOCAL MASS NUMBER (INTEGER VARIABLE OF A) |
---|
| 1092 | // 1841 C Z1 LOCAL NUCLEAR CHARGE (INTEGER VARIABLE OF Z) |
---|
| 1093 | // 1842 C REFOPT4 OPTION, SPECIFYING THE MASS FORMULA (SEE ABOVE) |
---|
| 1094 | // 1843 C A MASS NUMBER |
---|
| 1095 | // 1844 C Z NUCLEAR CHARGE |
---|
| 1096 | // 1845 C DEL PAIRING CORRECTION |
---|
| 1097 | // 1846 C EL BINDING ENERGY |
---|
| 1098 | // 1847 C ECNZ( , ) TABLE OF SHELL CORRECTIONS |
---|
| 1099 | // 1848 C----------------------------------------------------------------------- |
---|
| 1100 | // 1849 C |
---|
| 1101 | G4int a1 = idnint(a); |
---|
| 1102 | G4int z1 = idnint(z); |
---|
| 1103 | |
---|
| 1104 | if ( (a1 <= 0) || (z1 <= 0) || ((a1-z1) <= 0) ) { //then |
---|
| 1105 | // modif pour récupérer une masse p et n correcte: |
---|
| 1106 | (*el) = 0.0; |
---|
| 1107 | return; |
---|
| 1108 | // goto mglms50; |
---|
| 1109 | } |
---|
| 1110 | else { |
---|
| 1111 | // binding energy incl. pairing contr. is calculated from |
---|
| 1112 | // function eflmac |
---|
| 1113 | (*el) = eflmac(a1,z1,0,refopt4); |
---|
| 1114 | if (refopt4 > 0) { |
---|
| 1115 | if (refopt4 != 2) { |
---|
| 1116 | (*el) = (*el) + ec2sub->ecnz[a1-z1][z1]; |
---|
| 1117 | //(*el) = (*el) + ec2sub->ecnz[z1][a1-z1]; |
---|
| 1118 | } |
---|
| 1119 | } |
---|
| 1120 | } |
---|
| 1121 | return; |
---|
| 1122 | } |
---|
| 1123 | |
---|
| 1124 | G4double G4Abla::spdef(G4int a, G4int z, G4int optxfis) |
---|
| 1125 | { |
---|
| 1126 | |
---|
| 1127 | // INPUT: A,Z,OPTXFIS MASS AND CHARGE OF A NUCLEUS, |
---|
| 1128 | // OPTION FOR FISSILITY |
---|
| 1129 | // OUTPUT: SPDEF |
---|
| 1130 | |
---|
| 1131 | // ALPHA2 SADDLE POINT DEF. COHEN&SWIATECKI ANN.PHYS. 22 (1963) 406 |
---|
| 1132 | // RANGING FROM FISSILITY X=0.30 TO X=1.00 IN STEPS OF 0.02 |
---|
| 1133 | |
---|
[962] | 1134 | G4int index = 0; |
---|
| 1135 | G4double x = 0.0, v = 0.0, dx = 0.0; |
---|
[819] | 1136 | |
---|
| 1137 | const G4int alpha2Size = 37; |
---|
| 1138 | // The value 0.0 at alpha2[0] added by PK. |
---|
| 1139 | G4double alpha2[alpha2Size] = {0.0, 2.5464e0, 2.4944e0, 2.4410e0, 2.3915e0, 2.3482e0, |
---|
| 1140 | 2.3014e0, 2.2479e0, 2.1982e0, 2.1432e0, 2.0807e0, 2.0142e0, 1.9419e0, |
---|
| 1141 | 1.8714e0, 1.8010e0, 1.7272e0, 1.6473e0, 1.5601e0, 1.4526e0, 1.3164e0, |
---|
| 1142 | 1.1391e0, 0.9662e0, 0.8295e0, 0.7231e0, 0.6360e0, 0.5615e0, 0.4953e0, |
---|
| 1143 | 0.4354e0, 0.3799e0, 0.3274e0, 0.2779e0, 0.2298e0, 0.1827e0, 0.1373e0, |
---|
| 1144 | 0.0901e0, 0.0430e0, 0.0000e0}; |
---|
| 1145 | |
---|
| 1146 | dx = 0.02; |
---|
| 1147 | x = fissility(a,z,optxfis); |
---|
| 1148 | |
---|
| 1149 | if (x > 1.0) { |
---|
| 1150 | x = 1.0; |
---|
| 1151 | } |
---|
| 1152 | |
---|
| 1153 | if (x < 0.0) { |
---|
| 1154 | x = 0.0; |
---|
| 1155 | } |
---|
| 1156 | |
---|
| 1157 | v = (x - 0.3)/dx + 1.0; |
---|
| 1158 | index = idnint(v); |
---|
| 1159 | |
---|
| 1160 | if (index < 1) { |
---|
| 1161 | return(alpha2[1]); // alpha2[0] -> alpha2[1] |
---|
| 1162 | } |
---|
| 1163 | |
---|
| 1164 | if (index == 36) { //then // :::FIXME:: Possible off-by-one bug... |
---|
| 1165 | return(alpha2[36]); |
---|
| 1166 | } |
---|
| 1167 | else { |
---|
| 1168 | return(alpha2[index] + (alpha2[index+1] - alpha2[index]) / dx * ( x - (0.3e0 + dx*(index-1)))); //:::FIXME::: Possible off-by-one |
---|
| 1169 | } |
---|
| 1170 | |
---|
| 1171 | return alpha2[0]; // The algorithm is not supposed to reach this point. |
---|
| 1172 | } |
---|
| 1173 | |
---|
| 1174 | G4double G4Abla::fissility(int a,int z, int optxfis) |
---|
| 1175 | { |
---|
| 1176 | // CALCULATION OF FISSILITY PARAMETER |
---|
| 1177 | // |
---|
| 1178 | // INPUT: A,Z INTEGER MASS & CHARGE OF NUCLEUS |
---|
| 1179 | // OPTXFIS = 0 : MYERS, SWIATECKI |
---|
| 1180 | // 1 : DAHLINGER |
---|
| 1181 | // 2 : ANDREYEV |
---|
| 1182 | |
---|
[962] | 1183 | G4double aa = 0.0, zz = 0.0, i = 0.0; |
---|
[819] | 1184 | G4double fissilityResult = 0.0; |
---|
| 1185 | |
---|
| 1186 | aa = double(a); |
---|
| 1187 | zz = double(z); |
---|
| 1188 | i = double(a-2*z) / aa; |
---|
| 1189 | |
---|
| 1190 | // myers & swiatecki droplet modell |
---|
| 1191 | if (optxfis == 0) { //then |
---|
| 1192 | fissilityResult = std::pow(zz,2) / aa /50.8830e0 / (1.0e0 - 1.7826e0 * std::pow(i,2)); |
---|
| 1193 | } |
---|
| 1194 | |
---|
| 1195 | if (optxfis == 1) { |
---|
| 1196 | // dahlinger fit: |
---|
| 1197 | fissilityResult = std::pow(zz,2) / aa * std::pow((49.22e0*(1.e0 - 0.3803e0*std::pow(i,2) - 20.489e0*std::pow(i,4))),(-1)); |
---|
| 1198 | } |
---|
| 1199 | |
---|
| 1200 | if (optxfis == 2) { |
---|
| 1201 | // dubna fit: |
---|
| 1202 | fissilityResult = std::pow(zz,2) / aa /(48.e0*(1.e0 - 17.22e0*std::pow(i,4))); |
---|
| 1203 | } |
---|
| 1204 | |
---|
| 1205 | return fissilityResult; |
---|
| 1206 | } |
---|
| 1207 | |
---|
| 1208 | void G4Abla::evapora(G4double zprf, G4double aprf, G4double ee, G4double jprf, |
---|
| 1209 | G4double *zf_par, G4double *af_par, G4double *mtota_par, |
---|
| 1210 | G4double *pleva_par, G4double *pxeva_par, G4double *pyeva_par, |
---|
| 1211 | G4int *ff_par, G4int *inttype_par, G4int *inum_par) |
---|
| 1212 | { |
---|
| 1213 | G4double zf = (*zf_par); |
---|
| 1214 | G4double af = (*af_par); |
---|
| 1215 | G4double mtota = (*mtota_par); |
---|
| 1216 | G4double pleva = (*pleva_par); |
---|
| 1217 | G4double pxeva = (*pxeva_par); |
---|
| 1218 | G4double pyeva = (*pyeva_par); |
---|
| 1219 | G4int ff = (*ff_par); |
---|
| 1220 | G4int inttype = (*inttype_par); |
---|
| 1221 | G4int inum = (*inum_par); |
---|
| 1222 | |
---|
[962] | 1223 | G4int idebug = 0; |
---|
| 1224 | |
---|
| 1225 | if(idebug == 666) { |
---|
| 1226 | zprf = 81.; |
---|
| 1227 | aprf = 201.; |
---|
| 1228 | // ee = 86.5877686; |
---|
| 1229 | ee = 300.0; |
---|
| 1230 | jprf = 32.; |
---|
| 1231 | zf = 0.; |
---|
| 1232 | af = 0.; |
---|
| 1233 | mtota = 0.; |
---|
| 1234 | pleva = 0.; |
---|
| 1235 | pxeva = 0.; |
---|
| 1236 | pyeva = 0.; |
---|
| 1237 | ff = -1; |
---|
| 1238 | inttype = 0; |
---|
| 1239 | inum = 2; |
---|
| 1240 | G4cout <<" PK::: EVAPORA event " << inum << G4endl; |
---|
| 1241 | G4cout <<" PK::: zprf = " << zprf << G4endl; |
---|
| 1242 | G4cout <<" PK::: aprf = " << aprf << G4endl; |
---|
| 1243 | G4cout <<" PK::: ee = " << ee << G4endl; |
---|
| 1244 | G4cout <<" PK::: eeDiff = " << ee - 86.5877686 << G4endl; |
---|
| 1245 | G4cout <<" PK::: jprf = " << jprf << G4endl; |
---|
| 1246 | G4cout <<" PK::: zf = " << zf << G4endl; |
---|
| 1247 | G4cout <<" PK::: af = " << af << G4endl; |
---|
| 1248 | G4cout <<" PK::: mtota = " << mtota << G4endl; |
---|
| 1249 | G4cout <<" PK::: pleva = " << pleva << G4endl; |
---|
| 1250 | G4cout <<" PK::: pxeva = " << pxeva << G4endl; |
---|
| 1251 | G4cout <<" PK::: pyeva = " << pyeva << G4endl; |
---|
| 1252 | G4cout <<" PK::: ff = " << ff << G4endl; |
---|
| 1253 | G4cout <<" PK::: inttype = " << inttype << G4endl; |
---|
| 1254 | G4cout <<" PK::: inum = " << inum << G4endl; |
---|
| 1255 | } |
---|
[819] | 1256 | // 533 C |
---|
| 1257 | // 534 C INPUT: |
---|
| 1258 | // 535 C |
---|
| 1259 | // 536 C ZPRF, APRF, EE(EE IS MODIFIED!), JPRF |
---|
| 1260 | // 537 C |
---|
| 1261 | // 538 C PROJECTILE AND TARGET PARAMETERS + CROSS SECTIONS |
---|
| 1262 | // 539 C COMMON /ABRAMAIN/ AP,ZP,AT,ZT,EAP,BETA,BMAXNUC,CRTOT,CRNUC, |
---|
| 1263 | // 540 C R_0,R_P,R_T, IMAX,IRNDM,PI, |
---|
| 1264 | // 541 C BFPRO,SNPRO,SPPRO,SHELL |
---|
| 1265 | // 542 C |
---|
| 1266 | // 543 C AP,ZP,AT,ZT - PROJECTILE AND TARGET MASSES |
---|
| 1267 | // 544 C EAP,BETA - BEAM ENERGY PER NUCLEON, V/C |
---|
| 1268 | // 545 C BMAXNUC - MAX. IMPACT PARAMETER FOR NUCL. REAC. |
---|
| 1269 | // 546 C CRTOT,CRNUC - TOTAL AND NUCLEAR REACTION CROSS SECTION |
---|
| 1270 | // 547 C R_0,R_P,R_T, - RADIUS PARAMETER, PROJECTILE+ TARGET RADII |
---|
| 1271 | // 548 C IMAX,IRNDM,PI - MAXIMUM NUMBER OF EVENTS, DUMMY, 3.141... |
---|
| 1272 | // 549 C BFPRO - FISSION BARRIER OF THE PROJECTILE |
---|
| 1273 | // 550 C SNPRO - NEUTRON SEPARATION ENERGY OF THE PROJECTILE |
---|
| 1274 | // 551 C SPPRO - PROTON " " " " " |
---|
| 1275 | // 552 C SHELL - GROUND STATE SHELL CORRECTION |
---|
| 1276 | // 553 C |
---|
| 1277 | // 554 C--------------------------------------------------------------------- |
---|
| 1278 | // 555 C FISSION BARRIERS |
---|
| 1279 | // 556 C COMMON /FB/ EFA |
---|
| 1280 | // 557 C EFA - ARRAY OF FISSION BARRIERS |
---|
| 1281 | // 558 C--------------------------------------------------------------------- |
---|
| 1282 | // 559 C OUTPUT: |
---|
| 1283 | // 560 C ZF, AF, MTOTA, PLEVA, PTEVA, FF, INTTYPE, INUM |
---|
| 1284 | // 561 C |
---|
| 1285 | // 562 C ZF,AF - CHARGE AND MASS OF FINAL FRAGMENT AFTER EVAPORATION |
---|
| 1286 | // 563 C MTOTA _ NUMBER OF EVAPORATED ALPHAS |
---|
| 1287 | // 564 C PLEVA,PXEVA,PYEVA - MOMENTUM RECOIL BY EVAPORATION |
---|
| 1288 | // 565 C INTTYPE - TYPE OF REACTION 0/1 NUCLEAR OR ELECTROMAGNETIC |
---|
| 1289 | // 566 C FF - 0/1 NO FISSION / FISSION EVENT |
---|
| 1290 | // 567 C INUM - EVENTNUMBER |
---|
| 1291 | // 568 C ____________________________________________________________________ |
---|
| 1292 | // 569 C / |
---|
| 1293 | // 570 C / CALCUL DE LA MASSE ET CHARGE FINALES D'UNE CHAINE D'EVAPORATION |
---|
| 1294 | // 571 C / |
---|
| 1295 | // 572 C / PROCEDURE FOR CALCULATING THE FINAL MASS AND CHARGE VALUES OF A |
---|
| 1296 | // 573 C / SPECIFIC EVAPORATION CHAIN, STARTING POINT DEFINED BY (APRF, ZPRF, |
---|
| 1297 | // 574 C / EE) |
---|
| 1298 | // 575 C / On ajoute les 3 composantes de l'impulsion (PXEVA,PYEVA,PLEVA) |
---|
| 1299 | // 576 C / (actuellement PTEVA n'est pas correct; mauvaise norme...) |
---|
| 1300 | // 577 C /____________________________________________________________________ |
---|
| 1301 | // 578 C |
---|
| 1302 | // 612 C |
---|
| 1303 | // 613 C----------------------------------------------------------------------- |
---|
| 1304 | // 614 C IRNDM DUMMY ARGUMENT FOR RANDOM-NUMBER FUNCTION |
---|
| 1305 | // 615 C SORTIE LOCAL HELP VARIABLE TO END THE EVAPORATION CHAIN |
---|
| 1306 | // 616 C ZF NUCLEAR CHARGE OF THE FRAGMENT |
---|
| 1307 | // 617 C ZPRF NUCLEAR CHARGE OF THE PREFRAGMENT |
---|
| 1308 | // 618 C AF MASS NUMBER OF THE FRAGMENT |
---|
| 1309 | // 619 C APRF MASS NUMBER OF THE PREFRAGMENT |
---|
| 1310 | // 620 C EPSILN ENERGY BURNED IN EACH EVAPORATION STEP |
---|
| 1311 | // 621 C MALPHA LOCAL MASS CONTRIBUTION TO MTOTA IN EACH EVAPORATION |
---|
| 1312 | // 622 C STEP |
---|
| 1313 | // 623 C EE EXCITATION ENERGY (VARIABLE) |
---|
| 1314 | // 624 C PROBP PROTON EMISSION PROBABILITY |
---|
| 1315 | // 625 C PROBN NEUTRON EMISSION PROBABILITY |
---|
| 1316 | // 626 C PROBA ALPHA-PARTICLE EMISSION PROBABILITY |
---|
| 1317 | // 627 C PTOTL TOTAL EMISSION PROBABILITY |
---|
| 1318 | // 628 C E LOWEST PARTICLE-THRESHOLD ENERGY |
---|
| 1319 | // 629 C SN NEUTRON SEPARATION ENERGY |
---|
| 1320 | // 630 C SBP PROTON SEPARATION ENERGY PLUS EFFECTIVE COULOMB |
---|
| 1321 | // 631 C BARRIER |
---|
| 1322 | // 632 C SBA ALPHA-PARTICLE SEPARATION ENERGY PLUS EFFECTIVE |
---|
| 1323 | // 633 C COULOMB BARRIER |
---|
| 1324 | // 634 C BP EFFECTIVE PROTON COULOMB BARRIER |
---|
| 1325 | // 635 C BA EFFECTIVE ALPHA COULOMB BARRIER |
---|
| 1326 | // 636 C MTOTA TOTAL MASS OF THE EVAPORATED ALPHA PARTICLES |
---|
| 1327 | // 637 C X UNIFORM RANDOM NUMBER FOR NUCLEAR CHARGE |
---|
| 1328 | // 638 C AMOINS LOCAL MASS NUMBER OF EVAPORATED PARTICLE |
---|
| 1329 | // 639 C ZMOINS LOCAL NUCLEAR CHARGE OF EVAPORATED PARTICLE |
---|
| 1330 | // 640 C ECP KINETIC ENERGY OF PROTON WITHOUT COULOMB |
---|
| 1331 | // 641 C REPULSION |
---|
| 1332 | // 642 C ECN KINETIC ENERGY OF NEUTRON |
---|
| 1333 | // 643 C ECA KINETIC ENERGY OF ALPHA PARTICLE WITHOUT COULOMB |
---|
| 1334 | // 644 C REPULSION |
---|
| 1335 | // 645 C PLEVA TRANSVERSAL RECOIL MOMENTUM OF EVAPORATION |
---|
| 1336 | // 646 C PTEVA LONGITUDINAL RECOIL MOMENTUM OF EVAPORATION |
---|
| 1337 | // 647 C FF FISSION FLAG |
---|
| 1338 | // 648 C INTTYPE INTERACTION TYPE FLAG |
---|
| 1339 | // 649 C RNDX RECOIL MOMENTUM IN X-DIRECTION IN A SINGLE STEP |
---|
| 1340 | // 650 C RNDY RECOIL MOMENTUM IN Y-DIRECTION IN A SINGLE STEP |
---|
| 1341 | // 651 C RNDZ RECOIL MOMENTUM IN Z-DIRECTION IN A SINGLE STEP |
---|
| 1342 | // 652 C RNDN NORMALIZATION OF RECOIL MOMENTUM FOR EACH STEP |
---|
| 1343 | // 653 C----------------------------------------------------------------------- |
---|
| 1344 | // 654 C |
---|
| 1345 | // 655 SAVE |
---|
| 1346 | // SAVE -> static |
---|
| 1347 | |
---|
[962] | 1348 | static G4int sortie = 0; |
---|
| 1349 | static G4double epsiln = 0.0, probp = 0.0, probn = 0.0, proba = 0.0, ptotl = 0.0, e = 0.0; |
---|
| 1350 | static G4double sn = 0.0, sbp = 0.0, sba = 0.0, x = 0.0, amoins = 0.0, zmoins = 0.0; |
---|
| 1351 | G4double ecn = 0.0, ecp = 0.0,eca = 0.0, bp = 0.0, ba = 0.0; |
---|
| 1352 | static G4double pteva = 0.0; |
---|
[819] | 1353 | |
---|
[962] | 1354 | static G4int itest = 0; |
---|
| 1355 | static G4double probf = 0.0; |
---|
[819] | 1356 | |
---|
[962] | 1357 | static G4int k = 0, j = 0, il = 0; |
---|
[819] | 1358 | |
---|
[962] | 1359 | static G4double ctet1 = 0.0, stet1 = 0.0, phi1 = 0.0; |
---|
| 1360 | static G4double sbfis = 0.0, rnd = 0.0; |
---|
| 1361 | static G4double selmax = 0.0; |
---|
| 1362 | static G4double segs = 0.0; |
---|
| 1363 | static G4double ef = 0.0; |
---|
| 1364 | static G4int irndm = 0; |
---|
[819] | 1365 | |
---|
[962] | 1366 | static G4double pc = 0.0, malpha = 0.0; |
---|
[819] | 1367 | |
---|
| 1368 | zf = zprf; |
---|
| 1369 | af = aprf; |
---|
| 1370 | pleva = 0.0; |
---|
| 1371 | pteva = 0.0; |
---|
| 1372 | pxeva = 0.0; |
---|
| 1373 | pyeva = 0.0; |
---|
| 1374 | |
---|
| 1375 | sortie = 0; |
---|
| 1376 | ff = 0; |
---|
| 1377 | |
---|
| 1378 | itest = 0; |
---|
| 1379 | if (itest == 1) { |
---|
| 1380 | G4cout << "***************************" << G4endl; |
---|
| 1381 | } |
---|
| 1382 | |
---|
| 1383 | evapora10: |
---|
| 1384 | |
---|
| 1385 | if (itest == 1) { |
---|
| 1386 | G4cout <<"------zf,af,ee------" << idnint(zf) << "," << idnint(af) << "," << ee << G4endl; |
---|
| 1387 | } |
---|
| 1388 | |
---|
| 1389 | // calculation of the probabilities for the different decay channels |
---|
| 1390 | // plus separation energies and kinetic energies of the particles |
---|
| 1391 | direct(zf,af,ee,jprf,&probp,&probn,&proba,&probf,&ptotl, |
---|
| 1392 | &sn,&sbp,&sba,&ecn,&ecp,&eca,&bp,&ba,inttype,inum,itest); //:::FIXME::: Call |
---|
[962] | 1393 | if((eca+ba) < 0) { |
---|
| 1394 | eca = 0.0; |
---|
| 1395 | ba = 0.0; |
---|
| 1396 | } |
---|
[819] | 1397 | k = idnint(zf); |
---|
| 1398 | j = idnint(af-zf); |
---|
| 1399 | |
---|
| 1400 | // now ef is calculated from efa that depends on the subroutine |
---|
| 1401 | // barfit which takes into account the modification on the ang. mom. |
---|
| 1402 | // jb mvr 6-aug-1999 |
---|
| 1403 | // note *** shell correction! (ecgnz) jb mvr 20-7-1999 |
---|
| 1404 | il = idnint(jprf); |
---|
| 1405 | barfit(k,k+j,il,&sbfis,&segs,&selmax); |
---|
| 1406 | |
---|
| 1407 | if ((fiss->optshp == 1) || (fiss->optshp == 3)) { //then |
---|
| 1408 | // fb->efa[k][j] = G4double(sbfis) - ecld->ecgnz[j][k]; |
---|
| 1409 | fb->efa[j][k] = G4double(sbfis) - ecld->ecgnz[j][k]; |
---|
| 1410 | } |
---|
| 1411 | else { |
---|
| 1412 | fb->efa[j][k] = G4double(sbfis); |
---|
| 1413 | // fb->efa[j][k] = G4double(sbfis); |
---|
| 1414 | } //end if |
---|
| 1415 | ef = fb->efa[j][k]; |
---|
| 1416 | // ef = fb->efa[j][k]; |
---|
| 1417 | // here the final steps of the evaporation are calculated |
---|
| 1418 | if ((sortie == 1) || (ptotl == 0.e0)) { |
---|
| 1419 | e = dmin1(sn,sbp,sba); |
---|
| 1420 | if (e > 1.0e30) { |
---|
| 1421 | if(verboseLevel > 2) { |
---|
| 1422 | G4cout <<"erreur a la sortie evapora,e>1.e30,af=" << af <<" zf=" << zf << G4endl; |
---|
| 1423 | } |
---|
| 1424 | } |
---|
| 1425 | if (zf <= 6.0) { |
---|
| 1426 | goto evapora100; |
---|
| 1427 | } |
---|
| 1428 | if (e < 0.0) { |
---|
| 1429 | if (sn == e) { |
---|
| 1430 | af = af - 1.e0; |
---|
| 1431 | } |
---|
| 1432 | else if (sbp == e) { |
---|
| 1433 | af = af - 1.0; |
---|
| 1434 | zf = zf - 1.0; |
---|
| 1435 | } |
---|
| 1436 | else if (sba == e) { |
---|
| 1437 | af = af - 4.0; |
---|
| 1438 | zf = zf - 2.0; |
---|
| 1439 | } |
---|
| 1440 | if (af < 2.5) { |
---|
| 1441 | goto evapora100; |
---|
| 1442 | } |
---|
| 1443 | goto evapora10; |
---|
| 1444 | } |
---|
| 1445 | goto evapora100; |
---|
| 1446 | } |
---|
| 1447 | irndm = irndm + 1; |
---|
| 1448 | |
---|
| 1449 | // here the normal evaporation cascade starts |
---|
| 1450 | |
---|
| 1451 | // random number for the evaporation |
---|
| 1452 | // x = double(Rndm(irndm))*ptotl; |
---|
[962] | 1453 | // x = double(haz(1))*ptotl; |
---|
| 1454 | x = randomGenerator->getRandom() * ptotl; |
---|
[819] | 1455 | // G4cout <<"proba = " << proba << G4endl; |
---|
| 1456 | // G4cout <<"probp = " << probp << G4endl; |
---|
| 1457 | // G4cout <<"probn = " << probn << G4endl; |
---|
| 1458 | // G4cout <<"probf = " << probf << G4endl; |
---|
| 1459 | |
---|
| 1460 | itest = 0; |
---|
| 1461 | if (x < proba) { |
---|
| 1462 | // alpha evaporation |
---|
| 1463 | if (itest == 1) { |
---|
[962] | 1464 | G4cout <<"PK::: < alpha evaporation >" << G4endl; |
---|
[819] | 1465 | } |
---|
| 1466 | amoins = 4.0; |
---|
| 1467 | zmoins = 2.0; |
---|
| 1468 | epsiln = sba + eca; |
---|
| 1469 | pc = std::sqrt(std::pow((1.0 + (eca+ba)/3.72834e3),2) - 1.0) * 3.72834e3; |
---|
| 1470 | malpha = 4.0; |
---|
| 1471 | |
---|
| 1472 | // volant: |
---|
| 1473 | volant->iv = volant->iv + 1; |
---|
| 1474 | volant->acv[volant->iv] = 4.; |
---|
| 1475 | volant->zpcv[volant->iv] = 2.; |
---|
| 1476 | volant->pcv[volant->iv] = pc; |
---|
| 1477 | } |
---|
| 1478 | else if (x < proba+probp) { |
---|
| 1479 | // proton evaporation |
---|
| 1480 | if (itest == 1) { |
---|
[962] | 1481 | G4cout <<"PK::: < proton evaporation >" << G4endl; |
---|
[819] | 1482 | } |
---|
| 1483 | amoins = 1.0; |
---|
| 1484 | zmoins = 1.0; |
---|
| 1485 | epsiln = sbp + ecp; |
---|
| 1486 | pc = std::sqrt(std::pow((1.0 + (ecp + bp)/9.3827e2),2) - 1.0) * 9.3827e2; |
---|
| 1487 | malpha = 0.0; |
---|
| 1488 | // volant: |
---|
| 1489 | volant->iv = volant->iv + 1; |
---|
[962] | 1490 | volant->acv[volant->iv] = 1.0; |
---|
[819] | 1491 | volant->zpcv[volant->iv] = 1.; |
---|
| 1492 | volant->pcv[volant->iv] = pc; |
---|
| 1493 | } |
---|
| 1494 | else if (x < proba+probp+probn) { |
---|
| 1495 | // neutron evaporation |
---|
| 1496 | if (itest == 1) { |
---|
[962] | 1497 | G4cout <<"PK::: < neutron evaporation >" << G4endl; |
---|
[819] | 1498 | } |
---|
| 1499 | amoins = 1.0; |
---|
| 1500 | zmoins = 0.0; |
---|
| 1501 | epsiln = sn + ecn; |
---|
| 1502 | pc = std::sqrt(std::pow((1.0 + (ecn)/9.3956e2),2) - 1.0) * 9.3956e2; |
---|
[962] | 1503 | if(itest == 1) { |
---|
| 1504 | G4cout <<"PK::: pc " << pc << G4endl; |
---|
| 1505 | } |
---|
[819] | 1506 | malpha = 0.0; |
---|
| 1507 | |
---|
| 1508 | // volant: |
---|
| 1509 | volant->iv = volant->iv + 1; |
---|
| 1510 | volant->acv[volant->iv] = 1.; |
---|
| 1511 | volant->zpcv[volant->iv] = 0.; |
---|
| 1512 | volant->pcv[volant->iv] = pc; |
---|
[962] | 1513 | |
---|
| 1514 | if(volant->getTotalMass() > 209 && verboseLevel > 0) { |
---|
| 1515 | volant->dump(); |
---|
| 1516 | G4cout <<"DEBUGA Total = " << volant->getTotalMass() << G4endl; |
---|
| 1517 | } |
---|
[819] | 1518 | } |
---|
| 1519 | else { |
---|
| 1520 | // fission |
---|
| 1521 | // in case of fission-events the fragment nucleus is the mother nucleus |
---|
| 1522 | // before fission occurs with excitation energy above the fis.- barrier. |
---|
| 1523 | // fission fragment mass distribution is calulated in subroutine fisdis |
---|
| 1524 | if (itest == 1) { |
---|
[962] | 1525 | G4cout <<"PK::: < fission >" << G4endl; |
---|
[819] | 1526 | } |
---|
| 1527 | amoins = 0.0; |
---|
| 1528 | zmoins = 0.0; |
---|
| 1529 | epsiln = ef; |
---|
| 1530 | |
---|
| 1531 | malpha = 0.0; |
---|
| 1532 | pc = 0.0; |
---|
| 1533 | ff = 1; |
---|
[962] | 1534 | // ff = 0; // For testing, allows to disable fission! |
---|
[819] | 1535 | } |
---|
| 1536 | |
---|
| 1537 | if (itest == 1) { |
---|
[962] | 1538 | G4cout << std::setprecision(9) <<"PK::: SN,SBP,SBA,EF " << sn << " " << sbp << " " << sba <<" " << ef << G4endl; |
---|
| 1539 | G4cout << std::setprecision(9) <<"PK::: PROBN,PROBP,PROBA,PROBF,PTOTL " <<" "<< probn <<" "<< probp <<" "<< proba <<" "<< probf <<" "<< ptotl << G4endl; |
---|
[819] | 1540 | } |
---|
| 1541 | |
---|
| 1542 | // calculation of the daughter nucleus |
---|
| 1543 | af = af - amoins; |
---|
| 1544 | zf = zf - zmoins; |
---|
| 1545 | ee = ee - epsiln; |
---|
| 1546 | if (ee <= 0.01) { |
---|
| 1547 | ee = 0.01; |
---|
| 1548 | } |
---|
| 1549 | mtota = mtota + malpha; |
---|
| 1550 | |
---|
| 1551 | if(ff == 0) { |
---|
[962] | 1552 | rnd = randomGenerator->getRandom(); |
---|
| 1553 | // standardRandom(&rnd,&(hazard->igraine[8])); |
---|
[819] | 1554 | ctet1 = 2.0*rnd - 1.0; |
---|
[962] | 1555 | rnd = randomGenerator->getRandom(); |
---|
| 1556 | // standardRandom(&rnd,&(hazard->igraine[4])); |
---|
[819] | 1557 | phi1 = rnd*2.0*3.141592654; |
---|
| 1558 | stet1 = std::sqrt(1.0 - std::pow(ctet1,2)); |
---|
| 1559 | volant->xcv[volant->iv] = stet1*std::cos(phi1); |
---|
| 1560 | volant->ycv[volant->iv] = stet1*std::sin(phi1); |
---|
| 1561 | volant->zcv[volant->iv] = ctet1; |
---|
| 1562 | pxeva = pxeva - pc * volant->xcv[volant->iv]; |
---|
| 1563 | pyeva = pyeva - pc * volant->ycv[volant->iv]; |
---|
| 1564 | pleva = pleva - pc * ctet1; |
---|
| 1565 | } |
---|
| 1566 | |
---|
| 1567 | // condition for end of evaporation |
---|
| 1568 | if ((af < 2.5) || (ff == 1)) { |
---|
| 1569 | goto evapora100; |
---|
| 1570 | } |
---|
| 1571 | goto evapora10; |
---|
| 1572 | |
---|
| 1573 | evapora100: |
---|
| 1574 | (*zf_par) = zf; |
---|
| 1575 | (*af_par) = af; |
---|
| 1576 | (*mtota_par) = mtota; |
---|
| 1577 | (*pleva_par) = pleva; |
---|
| 1578 | (*pxeva_par) = pxeva; |
---|
| 1579 | (*pyeva_par) = pyeva; |
---|
| 1580 | (*ff_par) = ff; |
---|
| 1581 | (*inttype_par) = inttype; |
---|
| 1582 | (*inum_par) = inum; |
---|
| 1583 | |
---|
| 1584 | return; |
---|
| 1585 | } |
---|
| 1586 | |
---|
| 1587 | void G4Abla::direct(G4double zprf, G4double a, G4double ee, G4double jprf, |
---|
| 1588 | G4double *probp_par, G4double *probn_par, G4double *proba_par, |
---|
| 1589 | G4double *probf_par, G4double *ptotl_par, G4double *sn_par, |
---|
| 1590 | G4double *sbp_par, G4double *sba_par, G4double *ecn_par, |
---|
| 1591 | G4double *ecp_par,G4double *eca_par, G4double *bp_par, |
---|
| 1592 | G4double *ba_par, G4int inttype, G4int inum, G4int itest) |
---|
| 1593 | { |
---|
[962] | 1594 | G4int dummy0 = 0; |
---|
[819] | 1595 | |
---|
| 1596 | G4double probp = (*probp_par); |
---|
| 1597 | G4double probn = (*probn_par); |
---|
| 1598 | G4double proba = (*proba_par); |
---|
| 1599 | G4double probf = (*probf_par); |
---|
| 1600 | G4double ptotl = (*ptotl_par); |
---|
| 1601 | G4double sn = (*sn_par); |
---|
| 1602 | G4double sbp = (*sbp_par); |
---|
| 1603 | G4double sba = (*sba_par); |
---|
| 1604 | G4double ecn = (*ecn_par); |
---|
| 1605 | G4double ecp = (*ecp_par); |
---|
| 1606 | G4double eca = (*eca_par); |
---|
| 1607 | G4double bp = (*bp_par); |
---|
| 1608 | G4double ba = (*ba_par); |
---|
| 1609 | |
---|
| 1610 | // CALCULATION OF PARTICLE-EMISSION PROBABILITIES & FISSION / |
---|
| 1611 | // BASED ON THE SIMPLIFIED FORMULAS FOR THE DECAY WIDTH BY / |
---|
| 1612 | // MORETTO, ROCHESTER MEETING TO AVOID COMPUTING TIME / |
---|
| 1613 | // INTENSIVE INTEGRATION OF THE LEVEL DENSITIES / |
---|
| 1614 | // USES EFFECTIVE COULOMB BARRIERS AND AN AVERAGE KINETIC ENERGY/ |
---|
| 1615 | // OF THE EVAPORATED PARTICLES / |
---|
| 1616 | // COLLECTIVE ENHANCMENT OF THE LEVEL DENSITY IS INCLUDED / |
---|
| 1617 | // DYNAMICAL HINDRANCE OF FISSION IS INCLUDED BY A STEP FUNCTION/ |
---|
| 1618 | // APPROXIMATION. SEE A.R. JUNGHANS DIPLOMA THESIS / |
---|
| 1619 | // SHELL AND PAIRING STRUCTURES IN THE LEVEL DENSITY IS INCLUDED/ |
---|
| 1620 | |
---|
| 1621 | // INPUT: |
---|
| 1622 | // ZPRF,A,EE CHARGE, MASS, EXCITATION ENERGY OF COMPOUND |
---|
| 1623 | // NUCLEUS |
---|
| 1624 | // JPRF ROOT-MEAN-SQUARED ANGULAR MOMENTUM |
---|
| 1625 | |
---|
| 1626 | // DEFORMATIONS AND G.S. SHELL EFFECTS |
---|
| 1627 | // COMMON /ECLD/ ECGNZ,ECFNZ,VGSLD,ALPHA |
---|
| 1628 | |
---|
| 1629 | // ECGNZ - GROUND STATE SHELL CORR. FRLDM FOR A SPHERICAL G.S. |
---|
| 1630 | // ECFNZ - SHELL CORRECTION FOR THE SADDLE POINT (NOW: == 0) |
---|
| 1631 | // VGSLD - DIFFERENCE BETWEEN DEFORMED G.S. AND LDM VALUE |
---|
| 1632 | // ALPHA - ALPHA GROUND STATE DEFORMATION (THIS IS NOT BETA2!) |
---|
| 1633 | // BETA2 = SQRT((4PI)/5) * ALPHA |
---|
| 1634 | |
---|
| 1635 | //OPTIONS AND PARAMETERS FOR FISSION CHANNEL |
---|
| 1636 | //COMMON /FISS/ AKAP,BET,HOMEGA,KOEFF,IFIS, |
---|
| 1637 | // OPTSHP,OPTXFIS,OPTLES,OPTCOL |
---|
| 1638 | // |
---|
| 1639 | // AKAP - HBAR**2/(2* MN * R_0**2) = 10 MEV, R_0 = 1.4 FM |
---|
| 1640 | // BET - REDUCED NUCLEAR FRICTION COEFFICIENT IN (10**21 S**-1) |
---|
| 1641 | // HOMEGA - CURVATURE OF THE FISSION BARRIER = 1 MEV |
---|
| 1642 | // KOEFF - COEFFICIENT FOR THE LD FISSION BARRIER == 1.0 |
---|
| 1643 | // IFIS - 0/1 FISSION CHANNEL OFF/ON |
---|
| 1644 | // OPTSHP - INTEGER SWITCH FOR SHELL CORRECTION IN MASSES/ENERGY |
---|
| 1645 | // = 0 NO MICROSCOPIC CORRECTIONS IN MASSES AND ENERGY |
---|
| 1646 | // = 1 SHELL , NO PAIRING |
---|
| 1647 | // = 2 PAIRING, NO SHELL |
---|
| 1648 | // = 3 SHELL AND PAIRING |
---|
| 1649 | // OPTCOL - 0/1 COLLECTIVE ENHANCEMENT SWITCHED ON/OFF |
---|
| 1650 | // OPTXFIS- 0,1,2 FOR MYERS & SWIATECKI, DAHLINGER, ANDREYEV |
---|
| 1651 | // FISSILITY PARAMETER. |
---|
| 1652 | // OPTLES - CONSTANT TEMPERATURE LEVEL DENSITY FOR A,Z > TH-224 |
---|
| 1653 | // OPTCOL - 0/1 COLLECTIVE ENHANCEMENT OFF/ON |
---|
| 1654 | |
---|
| 1655 | // LEVEL DENSITY PARAMETERS |
---|
| 1656 | // COMMON /ALD/ AV,AS,AK,OPTAFAN |
---|
| 1657 | // AV,AS,AK - VOLUME,SURFACE,CURVATURE DEPENDENCE OF THE |
---|
| 1658 | // LEVEL DENSITY PARAMETER |
---|
| 1659 | // OPTAFAN - 0/1 AF/AN >=1 OR AF/AN ==1 |
---|
| 1660 | // RECOMMENDED IS OPTAFAN = 0 |
---|
| 1661 | |
---|
| 1662 | // FISSION BARRIERS |
---|
| 1663 | // COMMON /FB/ EFA |
---|
| 1664 | // EFA - ARRAY OF FISSION BARRIERS |
---|
| 1665 | |
---|
| 1666 | |
---|
| 1667 | // OUTPUT: PROBN,PROBP,PROBA,PROBF,PTOTL: |
---|
| 1668 | // - EMISSION PROBABILITIES FOR N EUTRON, P ROTON, A LPHA |
---|
| 1669 | // PARTICLES, F ISSION AND NORMALISATION |
---|
| 1670 | // SN,SBP,SBA: SEPARATION ENERGIES N P A |
---|
| 1671 | // INCLUDING EFFECTIVE BARRIERS |
---|
| 1672 | // ECN,ECP,ECA,BP,BA |
---|
| 1673 | // - AVERAGE KINETIC ENERGIES (2*T) AND EFFECTIVE BARRIERS |
---|
| 1674 | |
---|
[962] | 1675 | static G4double bk = 0.0; |
---|
| 1676 | static G4int afp = 0; |
---|
| 1677 | static G4double at = 0.0; |
---|
| 1678 | static G4double bs = 0.0; |
---|
| 1679 | static G4double bshell = 0.0; |
---|
| 1680 | static G4double cf = 0.0; |
---|
| 1681 | static G4double dconst = 0.0; |
---|
| 1682 | static G4double defbet = 0.0; |
---|
| 1683 | static G4double denomi = 0.0; |
---|
| 1684 | static G4double densa = 0.0; |
---|
| 1685 | static G4double densf = 0.0; |
---|
| 1686 | static G4double densg = 0.0; |
---|
| 1687 | static G4double densn = 0.0; |
---|
| 1688 | static G4double densp = 0.0; |
---|
| 1689 | static G4double edyn = 0.0; |
---|
| 1690 | static G4double eer = 0.0; |
---|
| 1691 | static G4double ef = 0.0; |
---|
| 1692 | static G4double ft = 0.0; |
---|
| 1693 | static G4double ga = 0.0; |
---|
| 1694 | static G4double gf = 0.0; |
---|
| 1695 | static G4double gn = 0.0; |
---|
| 1696 | static G4double gngf = 0.0; |
---|
| 1697 | static G4double gp = 0.0; |
---|
| 1698 | static G4double gsum = 0.0; |
---|
[819] | 1699 | static G4double hbar = 6.582122e-22; // = 0.0; |
---|
[962] | 1700 | static G4double iflag = 0.0; |
---|
| 1701 | static G4int il = 0; |
---|
| 1702 | static G4int imaxwell = 0; |
---|
| 1703 | static G4int in = 0; |
---|
| 1704 | static G4int iz = 0; |
---|
| 1705 | static G4int j = 0; |
---|
| 1706 | static G4int k = 0; |
---|
| 1707 | static G4double ma1z = 0.0; |
---|
| 1708 | static G4double ma1z1 = 0.0; |
---|
| 1709 | static G4double ma4z2 = 0.0; |
---|
| 1710 | static G4double maz = 0.0; |
---|
| 1711 | static G4double nprf = 0.0; |
---|
| 1712 | static G4double nt = 0.0; |
---|
| 1713 | static G4double parc = 0.0; |
---|
[819] | 1714 | static G4double pi = 3.14159265; |
---|
[962] | 1715 | static G4double pt = 0.0; |
---|
| 1716 | static G4double ra = 0.0; |
---|
| 1717 | static G4double rat = 0.0; |
---|
| 1718 | static G4double refmod = 0.0; |
---|
| 1719 | static G4double rf = 0.0; |
---|
| 1720 | static G4double rn = 0.0; |
---|
| 1721 | static G4double rnd = 0.0; |
---|
| 1722 | static G4double rnt = 0.0; |
---|
| 1723 | static G4double rp = 0.0; |
---|
| 1724 | static G4double rpt = 0.0; |
---|
| 1725 | static G4double sa = 0.0; |
---|
| 1726 | static G4double sbf = 0.0; |
---|
| 1727 | static G4double sbfis = 0.0; |
---|
| 1728 | static G4double segs = 0.0; |
---|
| 1729 | static G4double selmax = 0.0; |
---|
| 1730 | static G4double sp = 0.0; |
---|
| 1731 | static G4double tauc = 0.0; |
---|
| 1732 | static G4double tconst = 0.0; |
---|
| 1733 | static G4double temp = 0.0; |
---|
| 1734 | static G4double ts1 = 0.0; |
---|
| 1735 | static G4double tsum = 0.0; |
---|
| 1736 | static G4double wf = 0.0; |
---|
| 1737 | static G4double wfex = 0.0; |
---|
| 1738 | static G4double xx = 0.0; |
---|
| 1739 | static G4double y = 0.0; |
---|
[819] | 1740 | |
---|
| 1741 | imaxwell = 1; |
---|
| 1742 | inttype = 0; |
---|
| 1743 | |
---|
| 1744 | // limiting of excitation energy where fission occurs |
---|
| 1745 | // Note, this is not the dynamical hindrance (see end of routine) |
---|
| 1746 | edyn = 1000.0; |
---|
| 1747 | |
---|
| 1748 | // no limit if statistical model is calculated. |
---|
[962] | 1749 | if (fiss->bet <= 1.0e-16) { |
---|
[819] | 1750 | edyn = 10000.0; |
---|
| 1751 | } |
---|
| 1752 | |
---|
| 1753 | // just a change of name until the end of this subroutine |
---|
| 1754 | eer = ee; |
---|
[962] | 1755 | if(verboseLevel > 2) |
---|
| 1756 | G4cout << __FILE__ << ":" << __LINE__ << " eer = " << eer << G4endl; |
---|
[819] | 1757 | if (inum == 1) { |
---|
| 1758 | ilast = 1; |
---|
| 1759 | } |
---|
| 1760 | |
---|
| 1761 | // calculation of masses |
---|
| 1762 | // refmod = 1 ==> myers,swiatecki model |
---|
| 1763 | // refmod = 0 ==> weizsaecker model |
---|
| 1764 | refmod = 1; // Default = 1 |
---|
| 1765 | |
---|
| 1766 | if (refmod == 1) { |
---|
| 1767 | mglms(a,zprf,fiss->optshp,&maz); |
---|
| 1768 | mglms(a-1.0,zprf,fiss->optshp,&ma1z); |
---|
| 1769 | mglms(a-1.0,zprf-1.0,fiss->optshp,&ma1z1); |
---|
| 1770 | mglms(a-4.0,zprf-2.0,fiss->optshp,&ma4z2); |
---|
| 1771 | } |
---|
| 1772 | else { |
---|
| 1773 | mglw(a,zprf,&maz); |
---|
| 1774 | mglw(a-1.0,zprf,&ma1z); |
---|
| 1775 | mglw(a-1.0,zprf-1.0,&ma1z1); |
---|
| 1776 | mglw(a-4.0,zprf-2.0,&ma4z2); |
---|
| 1777 | } |
---|
| 1778 | |
---|
| 1779 | // separation energies and effective barriers |
---|
| 1780 | sn = ma1z - maz; |
---|
| 1781 | sp = ma1z1 - maz; |
---|
| 1782 | sa = ma4z2 - maz - 28.29688; |
---|
| 1783 | if (zprf < 1.0e0) { |
---|
| 1784 | sbp = 1.0e75; |
---|
| 1785 | goto direct30; |
---|
| 1786 | } |
---|
| 1787 | |
---|
| 1788 | // parameterisation gaimard: |
---|
| 1789 | // bp = 1.44*(zprf-1.d0)/(1.22*std::pow((a - 1.0),(1.0/3.0))+5.6) |
---|
| 1790 | // parameterisation khs (12-99) |
---|
| 1791 | bp = 1.44*(zprf - 1.0)/(2.1*std::pow((a - 1.0),(1.0/3.0)) + 0.0); |
---|
| 1792 | |
---|
| 1793 | sbp = sp + bp; |
---|
| 1794 | if (a-4.0 <= 0.0) { |
---|
| 1795 | sba = 1.0e+75; |
---|
| 1796 | goto direct30; |
---|
| 1797 | } |
---|
| 1798 | |
---|
| 1799 | // new effective barrier for alpha evaporation d=6.1: khs |
---|
| 1800 | // ba = 2.88d0*(zprf-2.d0)/(1.22d0*(a-4.d0)**(1.d0/3.d0)+6.1d0) |
---|
| 1801 | // parametrisation khs (12-99) |
---|
| 1802 | ba = 2.88*(zprf - 2.0)/(2.2*std::pow((a - 4.0),(1.0/3.0)) + 0.0); |
---|
| 1803 | |
---|
| 1804 | sba = sa + ba; |
---|
| 1805 | direct30: |
---|
| 1806 | |
---|
| 1807 | // calculation of surface and curvature integrals needed to |
---|
| 1808 | // to calculate the level density parameter (in densniv) |
---|
| 1809 | if (fiss->ifis > 0) { |
---|
| 1810 | k = idnint(zprf); |
---|
| 1811 | j = idnint(a - zprf); |
---|
| 1812 | |
---|
| 1813 | // now ef is calculated from efa that depends on the subroutine |
---|
| 1814 | // barfit which takes into account the modification on the ang. mom. |
---|
| 1815 | // jb mvr 6-aug-1999 |
---|
| 1816 | // note *** shell correction! (ecgnz) jb mvr 20-7-1999 |
---|
| 1817 | il = idnint(jprf); |
---|
| 1818 | barfit(k,k+j,il,&sbfis,&segs,&selmax); |
---|
| 1819 | if ((fiss->optshp == 1) || (fiss->optshp == 3)) { |
---|
| 1820 | // fb->efa[k][j] = G4double(sbfis) - ecld->ecgnz[j][k]; |
---|
| 1821 | // fb->efa[j][k] = G4double(sbfis) - ecld->ecgnz[j][k]; |
---|
| 1822 | fb->efa[j][k] = double(sbfis) - ecld->ecgnz[j][k]; |
---|
| 1823 | } |
---|
| 1824 | else { |
---|
| 1825 | // fb->efa[k][j] = G4double(sbfis); |
---|
| 1826 | fb->efa[j][k] = double(sbfis); |
---|
| 1827 | } |
---|
| 1828 | // ef = fb->efa[k][j]; |
---|
| 1829 | ef = fb->efa[j][k]; |
---|
| 1830 | |
---|
| 1831 | // to avoid negative values for impossible nuclei |
---|
| 1832 | // the fission barrier is set to zero if smaller than zero. |
---|
| 1833 | // if (fb->efa[k][j] < 0.0) { |
---|
| 1834 | // fb->efa[k][j] = 0.0; |
---|
| 1835 | // } |
---|
| 1836 | if (fb->efa[j][k] < 0.0) { |
---|
| 1837 | if(verboseLevel > 2) { |
---|
| 1838 | G4cout <<"Setting fission barrier to 0" << G4endl; |
---|
| 1839 | } |
---|
| 1840 | fb->efa[j][k] = 0.0; |
---|
| 1841 | } |
---|
| 1842 | |
---|
| 1843 | // factor with jprf should be 0.0025d0 - 0.01d0 for |
---|
| 1844 | // approximate influence of ang. momentum on bfis a.j. 22.07.96 |
---|
| 1845 | // 0.0 means no angular momentum |
---|
| 1846 | |
---|
| 1847 | if (ef < 0.0) { |
---|
| 1848 | ef = 0.0; |
---|
| 1849 | } |
---|
| 1850 | xx = fissility((k+j),k,fiss->optxfis); |
---|
| 1851 | |
---|
| 1852 | y = 1.00 - xx; |
---|
| 1853 | if (y < 0.0) { |
---|
| 1854 | y = 0.0; |
---|
| 1855 | } |
---|
| 1856 | if (y > 1.0) { |
---|
| 1857 | y = 1.0; |
---|
| 1858 | } |
---|
| 1859 | bs = bipol(1,y); |
---|
| 1860 | bk = bipol(2,y); |
---|
| 1861 | } |
---|
| 1862 | else { |
---|
| 1863 | ef = 1.0e40; |
---|
| 1864 | bs = 1.0; |
---|
| 1865 | bk = 1.0; |
---|
| 1866 | } |
---|
| 1867 | sbf = ee - ef; |
---|
| 1868 | |
---|
| 1869 | afp = idnint(a); |
---|
| 1870 | iz = idnint(zprf); |
---|
| 1871 | in = afp - iz; |
---|
| 1872 | bshell = ecld->ecfnz[in][iz]; |
---|
| 1873 | |
---|
| 1874 | // ld saddle point deformation |
---|
| 1875 | // here: beta2 = std::sqrt(5/(4pi)) * alpha2 |
---|
| 1876 | |
---|
| 1877 | // for the ground state def. 1.5d0 should be used |
---|
| 1878 | // because this was just the factor to produce the |
---|
| 1879 | // alpha-deformation table 1.5d0 should be used |
---|
| 1880 | // a.r.j. 6.8.97 |
---|
| 1881 | defbet = 1.58533e0 * spdef(idnint(a),idnint(zprf),fiss->optxfis); |
---|
| 1882 | |
---|
| 1883 | // level density and temperature at the saddle point |
---|
| 1884 | // G4cout <<"a = " << a << G4endl; |
---|
| 1885 | // G4cout <<"zprf = " << zprf << G4endl; |
---|
| 1886 | // G4cout <<"ee = " << ee << G4endl; |
---|
| 1887 | // G4cout <<"ef = " << ef << G4endl; |
---|
| 1888 | // G4cout <<"bshell = " << bshell << G4endl; |
---|
| 1889 | // G4cout <<"bs = " << bs << G4endl; |
---|
| 1890 | // G4cout <<"bk = " << bk << G4endl; |
---|
| 1891 | // G4cout <<"defbet = " << defbet << G4endl; |
---|
| 1892 | densniv(a,zprf,ee,ef,&densf,bshell,bs,bk,&temp,int(fiss->optshp),int(fiss->optcol),defbet); |
---|
| 1893 | // G4cout <<"densf = " << densf << G4endl; |
---|
| 1894 | // G4cout <<"temp = " << temp << G4endl; |
---|
| 1895 | ft = temp; |
---|
| 1896 | if (iz >= 2) { |
---|
| 1897 | bshell = ecld->ecgnz[in][iz-1] - ecld->vgsld[in][iz-1]; |
---|
| 1898 | defbet = 1.5 * (ecld->alpha[in][iz-1]); |
---|
| 1899 | |
---|
| 1900 | // level density and temperature in the proton daughter |
---|
| 1901 | densniv(a-1.0,zprf-1.0e0,ee,sbp,&densp, bshell,1.e0,1.e0,&temp,int(fiss->optshp),int(fiss->optcol),defbet); |
---|
| 1902 | pt = temp; |
---|
| 1903 | if (imaxwell == 1) { |
---|
| 1904 | // valentina - random kinetic energy in a maxwelliam distribution |
---|
| 1905 | // modif juin/2002 a.b. c.v. for light targets; limit on the energy |
---|
| 1906 | // from the maxwell distribution. |
---|
| 1907 | rpt = pt; |
---|
| 1908 | ecp = 2.0 * pt; |
---|
| 1909 | if(rpt <= 1.0e-3) { |
---|
| 1910 | goto direct2914; |
---|
| 1911 | } |
---|
| 1912 | iflag = 0; |
---|
| 1913 | direct1914: |
---|
| 1914 | ecp = fmaxhaz(rpt); |
---|
| 1915 | iflag = iflag + 1; |
---|
| 1916 | if(iflag >= 10) { |
---|
[962] | 1917 | rnd = randomGenerator->getRandom(); |
---|
| 1918 | // standardRandom(&rnd,&(hazard->igraine[5])); |
---|
[819] | 1919 | ecp=std::sqrt(rnd)*(eer-sbp); |
---|
| 1920 | goto direct2914; |
---|
| 1921 | } |
---|
| 1922 | if((ecp+sbp) > eer) { |
---|
| 1923 | goto direct1914; |
---|
| 1924 | } |
---|
| 1925 | } |
---|
| 1926 | else { |
---|
| 1927 | ecp = 2.0 * pt; |
---|
| 1928 | } |
---|
| 1929 | |
---|
| 1930 | direct2914: |
---|
| 1931 | dummy0 = 0; |
---|
| 1932 | // G4cout <<""<<G4endl; |
---|
| 1933 | } |
---|
| 1934 | else { |
---|
| 1935 | densp = 0.0; |
---|
| 1936 | ecp = 0.0; |
---|
| 1937 | pt = 0.0; |
---|
| 1938 | } |
---|
| 1939 | |
---|
| 1940 | if (in >= 2) { |
---|
| 1941 | bshell = ecld->ecgnz[in-1][iz] - ecld->vgsld[in-1][iz]; |
---|
| 1942 | defbet = 1.5e0 * (ecld->alpha[in-1][iz]); |
---|
| 1943 | |
---|
| 1944 | // level density and temperature in the neutron daughter |
---|
| 1945 | densniv(a-1.0,zprf,ee,sn,&densn,bshell, 1.e0,1.e0,&temp,int(fiss->optshp),int(fiss->optcol),defbet); |
---|
| 1946 | nt = temp; |
---|
[962] | 1947 | if(itest == 1) |
---|
| 1948 | G4cout <<"PK::: nt = " << nt <<G4endl; |
---|
[819] | 1949 | if (imaxwell == 1) { |
---|
| 1950 | // valentina - random kinetic energy in a maxwelliam distribution |
---|
| 1951 | // modif juin/2002 a.b. c.v. for light targets; limit on the energy |
---|
| 1952 | // from the maxwell distribution. |
---|
| 1953 | rnt = nt; |
---|
| 1954 | ecn = 2.0 * nt; |
---|
| 1955 | if(rnt <= 1.e-3) { |
---|
| 1956 | goto direct2915; |
---|
| 1957 | } |
---|
| 1958 | |
---|
| 1959 | iflag=0; |
---|
[962] | 1960 | direct1915: |
---|
[819] | 1961 | ecn = fmaxhaz(rnt); |
---|
[962] | 1962 | if(verboseLevel > 2) { |
---|
| 1963 | G4cout <<"rnt = " << rnt << G4endl; |
---|
| 1964 | G4cout << __FILE__ << ":" << __LINE__ << " ecn = " << ecn << G4endl; |
---|
| 1965 | } |
---|
[819] | 1966 | iflag=iflag+1; |
---|
| 1967 | if(iflag >= 10) { |
---|
[962] | 1968 | rnd = randomGenerator->getRandom(); |
---|
| 1969 | // standardRandom(&rnd,&(hazard->igraine[6])); |
---|
[819] | 1970 | ecn = std::sqrt(rnd)*(eer-sn); |
---|
[962] | 1971 | if(verboseLevel > 2) |
---|
| 1972 | G4cout << __FILE__ << ":" << __LINE__ << " ecn = " << ecn << G4endl; |
---|
[819] | 1973 | goto direct2915; |
---|
| 1974 | } |
---|
[962] | 1975 | if((ecn+sn) > eer) { |
---|
| 1976 | goto direct1915; |
---|
| 1977 | } |
---|
| 1978 | } |
---|
| 1979 | else { |
---|
| 1980 | ecn = 2.e0 * nt; |
---|
| 1981 | if(verboseLevel > 2) |
---|
| 1982 | G4cout << __FILE__ << ":" << __LINE__ << " ecn = " << ecn << G4endl; |
---|
| 1983 | } |
---|
| 1984 | // if((ecn + sn) <= eer) { |
---|
| 1985 | // ecn = 2.0 * nt; |
---|
| 1986 | // G4cout << __FILE__ << ":" << __LINE__ << " ecn = " << ecn << G4endl; |
---|
| 1987 | // } |
---|
[819] | 1988 | direct2915: |
---|
| 1989 | dummy0 = 0; |
---|
| 1990 | // G4cout <<"" <<G4endl; |
---|
[962] | 1991 | } |
---|
[819] | 1992 | else { |
---|
| 1993 | densn = 0.0; |
---|
| 1994 | ecn = 0.0; |
---|
| 1995 | nt = 0.0; |
---|
| 1996 | } |
---|
| 1997 | |
---|
| 1998 | if ((in >= 3) && (iz >= 3)) { |
---|
| 1999 | bshell = ecld->ecgnz[in-2][iz-2] - ecld->vgsld[in-2][iz-2]; |
---|
| 2000 | defbet = 1.5 * (ecld->alpha[in-2][iz-2]); |
---|
| 2001 | |
---|
[962] | 2002 | // For debugging: |
---|
| 2003 | // bshell = -10.7; |
---|
| 2004 | // defbet = -0.06105; |
---|
| 2005 | // G4cout <<"ecgnz N = " << in-2 << G4endl; |
---|
| 2006 | // G4cout <<"ecgnz Z = " << iz-2 << G4endl; |
---|
| 2007 | // G4cout <<"bshell = " << bshell << G4endl; |
---|
| 2008 | // G4cout <<"defbet = " << defbet << G4endl; |
---|
[819] | 2009 | // level density and temperature in the alpha daughter |
---|
| 2010 | densniv(a-4.0,zprf-2.0e0,ee,sba,&densa,bshell,1.e0,1.e0,&temp,int(fiss->optshp),int(fiss->optcol),defbet); |
---|
[962] | 2011 | // G4cout <<"densa = " << densa << G4endl; |
---|
| 2012 | // G4cout <<"temp = " << temp << G4endl; |
---|
[819] | 2013 | |
---|
| 2014 | // valentina - random kinetic energy in a maxwelliam distribution |
---|
| 2015 | at = temp; |
---|
| 2016 | if (imaxwell == 1) { |
---|
| 2017 | // modif juin/2002 a.b. c.v. for light targets; limit on the energy |
---|
| 2018 | // from the maxwell distribution. |
---|
| 2019 | rat = at; |
---|
| 2020 | eca= 2.e0 * at; |
---|
| 2021 | if(rat <= 1.e-3) { |
---|
| 2022 | goto direct2916; |
---|
| 2023 | } |
---|
| 2024 | iflag=0; |
---|
| 2025 | direct1916: |
---|
| 2026 | eca = fmaxhaz(rat); |
---|
| 2027 | iflag=iflag+1; |
---|
| 2028 | if(iflag >= 10) { |
---|
[962] | 2029 | rnd = randomGenerator->getRandom(); |
---|
| 2030 | // standardRandom(&rnd,&(hazard->igraine[7])); |
---|
[819] | 2031 | eca=std::sqrt(rnd)*(eer-sba); |
---|
| 2032 | goto direct2916; |
---|
| 2033 | } |
---|
| 2034 | if((eca+sba) > eer) { |
---|
| 2035 | goto direct1916; |
---|
| 2036 | } |
---|
[962] | 2037 | } |
---|
| 2038 | else { |
---|
| 2039 | eca = 2.0 * at; |
---|
| 2040 | } |
---|
[819] | 2041 | direct2916: |
---|
| 2042 | dummy0 = 0; |
---|
| 2043 | // G4cout <<"" << G4endl; |
---|
[962] | 2044 | } |
---|
| 2045 | else { |
---|
| 2046 | densa = 0.0; |
---|
| 2047 | eca = 0.0; |
---|
| 2048 | at = 0.0; |
---|
| 2049 | } |
---|
| 2050 | //} // PK |
---|
[819] | 2051 | |
---|
| 2052 | // special treatment for unbound nuclei |
---|
| 2053 | if (sn < 0.0) { |
---|
| 2054 | probn = 1.0; |
---|
| 2055 | probp = 0.0; |
---|
| 2056 | proba = 0.0; |
---|
| 2057 | probf = 0.0; |
---|
| 2058 | goto direct70; |
---|
| 2059 | } |
---|
| 2060 | if (sbp < 0.0) { |
---|
| 2061 | probp = 1.0; |
---|
| 2062 | probn = 0.0; |
---|
| 2063 | proba = 0.0; |
---|
| 2064 | probf = 0.0; |
---|
| 2065 | goto direct70; |
---|
| 2066 | } |
---|
| 2067 | |
---|
| 2068 | if ((a < 50.e0) || (ee > edyn)) { // no fission if e*> edyn or mass < 50 |
---|
| 2069 | // G4cout <<"densf = 0.0" << G4endl; |
---|
| 2070 | densf = 0.e0; |
---|
| 2071 | } |
---|
| 2072 | |
---|
| 2073 | bshell = ecld->ecgnz[in][iz] - ecld->vgsld[in][iz]; |
---|
| 2074 | defbet = 1.5e0 * (ecld->alpha[in][iz]); |
---|
| 2075 | |
---|
| 2076 | // compound nucleus level density |
---|
| 2077 | densniv(a,zprf,ee,0.0e0,&densg,bshell,1.e0,1.e0,&temp,int(fiss->optshp),int(fiss->optcol),defbet); |
---|
| 2078 | |
---|
| 2079 | if ( densg > 0.e0) { |
---|
| 2080 | // calculation of the partial decay width |
---|
| 2081 | // used for both the time scale and the evaporation decay width |
---|
| 2082 | gp = (std::pow(a,(2.0/3.0))/fiss->akap)*densp/densg/pi*std::pow(pt,2); |
---|
| 2083 | gn = (std::pow(a,(2.0/3.0))/fiss->akap)*densn/densg/pi*std::pow(nt,2); |
---|
| 2084 | ga = (std::pow(a,(2.0/3.0))/fiss->akap)*densa/densg/pi*2.0*std::pow(at,2); |
---|
| 2085 | gf = densf/densg/pi/2.0*ft; |
---|
| 2086 | |
---|
| 2087 | if(itest == 1) { |
---|
| 2088 | G4cout <<"gn,gp,ga,gf " << gn <<","<< gp <<","<< ga <<","<< gf << G4endl; |
---|
| 2089 | } |
---|
| 2090 | } |
---|
| 2091 | else { |
---|
| 2092 | if(verboseLevel > 2) { |
---|
| 2093 | G4cout <<"direct: densg <= 0.e0 " << a <<","<< zprf <<","<< ee << G4endl; |
---|
| 2094 | } |
---|
| 2095 | } |
---|
| 2096 | |
---|
| 2097 | gsum = ga + gp + gn; |
---|
| 2098 | if (gsum > 0.0) { |
---|
| 2099 | ts1 = hbar / gsum; |
---|
| 2100 | } |
---|
| 2101 | else { |
---|
| 2102 | ts1 = 1.0e99; |
---|
| 2103 | } |
---|
| 2104 | |
---|
| 2105 | if (inum > ilast) { // new event means reset the time scale |
---|
| 2106 | tsum = 0; |
---|
| 2107 | } |
---|
| 2108 | |
---|
| 2109 | // calculate the relative probabilities for all decay channels |
---|
| 2110 | if (densf == 0.0) { |
---|
| 2111 | if (densp == 0.0) { |
---|
| 2112 | if (densn == 0.0) { |
---|
| 2113 | if (densa == 0.0) { |
---|
| 2114 | // no reaction is possible |
---|
| 2115 | probf = 0.0; |
---|
| 2116 | probp = 0.0; |
---|
| 2117 | probn = 0.0; |
---|
| 2118 | proba = 0.0; |
---|
| 2119 | goto direct70; |
---|
| 2120 | } |
---|
| 2121 | |
---|
| 2122 | // alpha evaporation is the only open channel |
---|
| 2123 | rf = 0.0; |
---|
| 2124 | rp = 0.0; |
---|
| 2125 | rn = 0.0; |
---|
| 2126 | ra = 1.0; |
---|
| 2127 | goto direct50; |
---|
| 2128 | } |
---|
| 2129 | |
---|
| 2130 | // alpha emission and neutron emission |
---|
| 2131 | rf = 0.0; |
---|
| 2132 | rp = 0.0; |
---|
| 2133 | rn = 1.0; |
---|
| 2134 | ra = densa*2.0/densn*std::pow((at/nt),2); |
---|
| 2135 | goto direct50; |
---|
| 2136 | } |
---|
| 2137 | // alpha, proton and neutron emission |
---|
| 2138 | rf = 0.0; |
---|
| 2139 | rp = 1.0; |
---|
| 2140 | rn = densn/densp*std::pow((nt/pt),2); |
---|
| 2141 | ra = densa*2.0/densp*std::pow((at/pt),2); |
---|
| 2142 | goto direct50; |
---|
| 2143 | } |
---|
| 2144 | |
---|
| 2145 | // here fission has taken place |
---|
| 2146 | rf = 1.0; |
---|
| 2147 | |
---|
| 2148 | // cramers and weidenmueller factors for the dynamical hindrances of |
---|
| 2149 | // fission |
---|
[962] | 2150 | if (fiss->bet <= 1.0e-16) { |
---|
[819] | 2151 | cf = 1.0; |
---|
| 2152 | wf = 1.0; |
---|
| 2153 | } |
---|
| 2154 | else if (sbf > 0.0e0) { |
---|
[962] | 2155 | cf = cram(fiss->bet,fiss->homega); |
---|
[819] | 2156 | // if fission barrier ef=0.d0 then fission is the only possible |
---|
| 2157 | // channel. to avoid std::log(0) in function tau |
---|
| 2158 | // a.j. 7/28/93 |
---|
| 2159 | if (ef <= 0.0) { |
---|
| 2160 | rp = 0.0; |
---|
| 2161 | rn = 0.0; |
---|
| 2162 | ra = 0.0; |
---|
| 2163 | goto direct50; |
---|
| 2164 | } |
---|
| 2165 | else { |
---|
| 2166 | // transient time tau() |
---|
[962] | 2167 | tauc = tau(fiss->bet,fiss->homega,ef,ft); |
---|
[819] | 2168 | } |
---|
| 2169 | wfex = (tauc - tsum)/ts1; |
---|
| 2170 | |
---|
| 2171 | if (wfex < 0.0) { |
---|
| 2172 | wf = 1.0; |
---|
| 2173 | } |
---|
| 2174 | else { |
---|
| 2175 | wf = std::exp( -wfex); |
---|
| 2176 | } |
---|
| 2177 | } |
---|
| 2178 | else { |
---|
| 2179 | cf=1.0; |
---|
| 2180 | wf=1.0; |
---|
| 2181 | } |
---|
| 2182 | |
---|
| 2183 | if(verboseLevel > 2) { |
---|
| 2184 | G4cout <<"tsum,wf,cf " << tsum <<","<< wf <<","<< cf << G4endl; |
---|
| 2185 | } |
---|
| 2186 | |
---|
| 2187 | tsum = tsum + ts1; |
---|
| 2188 | |
---|
| 2189 | // change by g.k. and a.h. 5.9.95 |
---|
| 2190 | tconst = 0.7; |
---|
| 2191 | dconst = 12.0/std::sqrt(a); |
---|
| 2192 | nprf = a - zprf; |
---|
| 2193 | |
---|
| 2194 | if (fiss->optshp >= 2) { //then |
---|
| 2195 | parite(nprf,&parc); |
---|
| 2196 | dconst = dconst*parc; |
---|
| 2197 | } |
---|
| 2198 | else { |
---|
| 2199 | dconst= 0.0; |
---|
| 2200 | } |
---|
| 2201 | if ((ee <= 17.e0) && (fiss->optles == 1) && (iz >= 90) && (in >= 134)) { //then |
---|
| 2202 | // constant changed to 5.0 accord to moretto & vandenbosch a.j. 19.3.96 |
---|
| 2203 | gngf = std::pow(a,(2.0/3.0))*tconst/10.0*std::exp((ef-sn+dconst)/tconst); |
---|
| 2204 | |
---|
| 2205 | // if the excitation energy is so low that densn=0 ==> gn = 0 |
---|
| 2206 | // fission remains the only channel. |
---|
| 2207 | // a. j. 10.1.94 |
---|
| 2208 | if (gn == 0.0) { |
---|
| 2209 | rn = 0.0; |
---|
| 2210 | rp = 0.0; |
---|
| 2211 | ra = 0.0; |
---|
| 2212 | } |
---|
| 2213 | else { |
---|
| 2214 | rn=gngf; |
---|
| 2215 | rp=gngf*gp/gn; |
---|
| 2216 | ra=gngf*ga/gn; |
---|
| 2217 | } |
---|
| 2218 | } else { |
---|
| 2219 | rn = gn/(gf*cf); |
---|
| 2220 | // G4cout <<"rn = " << G4endl; |
---|
| 2221 | // G4cout <<"gn = " << gn << " gf = " << gf << " cf = " << cf << G4endl; |
---|
| 2222 | rp = gp/(gf*cf); |
---|
| 2223 | ra = ga/(gf*cf); |
---|
| 2224 | } |
---|
| 2225 | direct50: |
---|
| 2226 | // relative decay probabilities |
---|
| 2227 | denomi = rp+rn+ra+rf; |
---|
| 2228 | // decay probabilities after transient time |
---|
| 2229 | probf = rf/denomi; |
---|
| 2230 | probp = rp/denomi; |
---|
| 2231 | probn = rn/denomi; |
---|
| 2232 | proba = ra/denomi; |
---|
| 2233 | |
---|
| 2234 | // new treatment of grange-weidenmueller factor, 5.1.2000, khs !!! |
---|
| 2235 | |
---|
| 2236 | // decay probabilites with transient time included |
---|
| 2237 | probf = probf * wf; |
---|
| 2238 | if(probf == 1.0) { |
---|
| 2239 | probp = 0.0; |
---|
| 2240 | probn = 0.0; |
---|
| 2241 | proba = 0.0; |
---|
| 2242 | } |
---|
| 2243 | else { |
---|
| 2244 | probp = probp * (wf + (1.e0-wf)/(1.e0-probf)); |
---|
| 2245 | probn = probn * (wf + (1.e0-wf)/(1.e0-probf)); |
---|
| 2246 | proba = proba * (wf + (1.e0-wf)/(1.e0-probf)); |
---|
| 2247 | } |
---|
| 2248 | direct70: |
---|
| 2249 | ptotl = probp+probn+proba+probf; |
---|
| 2250 | |
---|
| 2251 | ee = eer; |
---|
| 2252 | ilast = inum; |
---|
| 2253 | |
---|
| 2254 | // Return values: |
---|
| 2255 | (*probp_par) = probp; |
---|
| 2256 | (*probn_par) = probn; |
---|
| 2257 | (*proba_par) = proba; |
---|
| 2258 | (*probf_par) = probf; |
---|
| 2259 | (*ptotl_par) = ptotl; |
---|
| 2260 | (*sn_par) = sn; |
---|
| 2261 | (*sbp_par) = sbp; |
---|
| 2262 | (*sba_par) = sba; |
---|
| 2263 | (*ecn_par) = ecn; |
---|
| 2264 | (*ecp_par) = ecp; |
---|
| 2265 | (*eca_par) = eca; |
---|
| 2266 | (*bp_par) = bp; |
---|
| 2267 | (*ba_par) = ba; |
---|
| 2268 | } |
---|
| 2269 | |
---|
| 2270 | void G4Abla::densniv(G4double a, G4double z, G4double ee, G4double esous, G4double *dens, G4double bshell, G4double bs, G4double bk, |
---|
| 2271 | G4double *temp, G4int optshp, G4int optcol, G4double defbet) |
---|
| 2272 | { |
---|
| 2273 | // 1498 C |
---|
| 2274 | // 1499 C INPUT: |
---|
| 2275 | // 1500 C A,EE,ESOUS,OPTSHP,BS,BK,BSHELL,DEFBET |
---|
| 2276 | // 1501 C |
---|
| 2277 | // 1502 C LEVEL DENSITY PARAMETERS |
---|
| 2278 | // 1503 C COMMON /ALD/ AV,AS,AK,OPTAFAN |
---|
| 2279 | // 1504 C AV,AS,AK - VOLUME,SURFACE,CURVATURE DEPENDENCE OF THE |
---|
| 2280 | // 1505 C LEVEL DENSITY PARAMETER |
---|
| 2281 | // 1506 C OPTAFAN - 0/1 AF/AN >=1 OR AF/AN ==1 |
---|
| 2282 | // 1507 C RECOMMENDED IS OPTAFAN = 0 |
---|
| 2283 | // 1508 C--------------------------------------------------------------------- |
---|
| 2284 | // 1509 C OUTPUT: DENS,TEMP |
---|
| 2285 | // 1510 C |
---|
| 2286 | // 1511 C ____________________________________________________________________ |
---|
| 2287 | // 1512 C / |
---|
| 2288 | // 1513 C / PROCEDURE FOR CALCULATING THE STATE DENSITY OF A COMPOUND NUCLEUS |
---|
| 2289 | // 1514 C /____________________________________________________________________ |
---|
| 2290 | // 1515 C |
---|
| 2291 | // 1516 INTEGER AFP,IZ,OPTSHP,OPTCOL,J,OPTAFAN |
---|
| 2292 | // 1517 REAL*8 A,EE,ESOUS,DENS,E,Y0,Y1,Y2,Y01,Y11,Y21,PA,BS,BK,TEMP |
---|
| 2293 | // 1518 C=====INSERTED BY KUDYAEV=============================================== |
---|
| 2294 | // 1519 COMMON /ALD/ AV,AS,AK,OPTAFAN |
---|
| 2295 | // 1520 REAL*8 ECR,ER,DELTAU,Z,DELTPP,PARA,PARZ,FE,HE,ECOR,ECOR1,Pi6 |
---|
| 2296 | // 1521 REAL*8 BSHELL,DELTA0,AV,AK,AS,PONNIV,PONFE,DEFBET,QR,SIG,FP |
---|
| 2297 | // 1522 C======================================================================= |
---|
| 2298 | // 1523 C |
---|
| 2299 | // 1524 C |
---|
| 2300 | // 1525 C----------------------------------------------------------------------- |
---|
| 2301 | // 1526 C A MASS NUMBER OF THE DAUGHTER NUCLEUS |
---|
| 2302 | // 1527 C EE EXCITATION ENERGY OF THE MOTHER NUCLEUS |
---|
| 2303 | // 1528 C ESOUS SEPARATION ENERGY PLUS EFFECTIVE COULOMB BARRIER |
---|
| 2304 | // 1529 C DENS STATE DENSITY OF DAUGHTER NUCLEUS AT EE-ESOUS-EC |
---|
| 2305 | // 1530 C BSHELL SHELL CORRECTION |
---|
| 2306 | // 1531 C TEMP NUCLEAR TEMPERATURE |
---|
| 2307 | // 1532 C E LOCAL EXCITATION ENERGY OF THE DAUGHTER NUCLEUS |
---|
| 2308 | // 1533 C E1 LOCAL HELP VARIABLE |
---|
| 2309 | // 1534 C Y0,Y1,Y2,Y01,Y11,Y21 |
---|
| 2310 | // 1535 C LOCAL HELP VARIABLES |
---|
| 2311 | // 1536 C PA LOCAL STATE-DENSITY PARAMETER |
---|
| 2312 | // 1537 C EC KINETIC ENERGY OF EMITTED PARTICLE WITHOUT |
---|
| 2313 | // 1538 C COULOMB REPULSION |
---|
| 2314 | // 1539 C IDEN FAKTOR FOR SUBSTRACTING KINETIC ENERGY IDEN*TEMP |
---|
| 2315 | // 1540 C DELTA0 PAIRING GAP 12 FOR GROUND STATE |
---|
| 2316 | // 1541 C 14 FOR SADDLE POINT |
---|
| 2317 | // 1542 C EITERA HELP VARIABLE FOR TEMPERATURE ITERATION |
---|
| 2318 | // 1543 C----------------------------------------------------------------------- |
---|
| 2319 | // 1544 C |
---|
| 2320 | // 1545 C |
---|
[962] | 2321 | G4double afp = 0.0; |
---|
| 2322 | G4double delta0 = 0.0; |
---|
| 2323 | G4double deltau = 0.0; |
---|
| 2324 | G4double deltpp = 0.0; |
---|
| 2325 | G4double e = 0.0; |
---|
[819] | 2326 | G4double ecor = 0.0; |
---|
[962] | 2327 | G4double ecor1 = 0.0; |
---|
| 2328 | G4double ecr = 0.0; |
---|
| 2329 | G4double er = 0.0; |
---|
| 2330 | G4double fe = 0.0; |
---|
| 2331 | G4double fp = 0.0; |
---|
| 2332 | G4double he = 0.0; |
---|
| 2333 | G4double iz = 0.0; |
---|
| 2334 | G4double pa = 0.0; |
---|
| 2335 | G4double para = 0.0; |
---|
| 2336 | G4double parz = 0.0; |
---|
| 2337 | G4double ponfe = 0.0; |
---|
| 2338 | G4double ponniv = 0.0; |
---|
| 2339 | G4double qr = 0.0; |
---|
| 2340 | G4double sig = 0.0; |
---|
| 2341 | G4double y01 = 0.0; |
---|
| 2342 | G4double y11 = 0.0; |
---|
| 2343 | G4double y2 = 0.0; |
---|
| 2344 | G4double y21 = 0.0; |
---|
| 2345 | G4double y1 = 0.0; |
---|
| 2346 | G4double y0 = 0.0; |
---|
[819] | 2347 | |
---|
| 2348 | G4double pi6 = std::pow(3.1415926535,2) / 6.0; |
---|
| 2349 | ecr=10.0; |
---|
| 2350 | er=28.0; |
---|
| 2351 | afp=idnint(a); |
---|
| 2352 | iz=idnint(z); |
---|
| 2353 | |
---|
| 2354 | // level density parameter |
---|
| 2355 | if((ald->optafan == 1)) { |
---|
| 2356 | pa = (ald->av)*a + (ald->as)*std::pow(a,(2.e0/3.e0)) + (ald->ak)*std::pow(a,(1.e0/3.e0)); |
---|
| 2357 | } |
---|
| 2358 | else { |
---|
| 2359 | pa = (ald->av)*a + (ald->as)*bs*std::pow(a,(2.e0/3.e0)) + (ald->ak)*bk*std::pow(a,(1.e0/3.e0)); |
---|
| 2360 | } |
---|
| 2361 | |
---|
| 2362 | fp = 0.01377937231e0 * std::pow(a,(5.e0/3.e0)) * (1.e0 + defbet/3.e0); |
---|
| 2363 | |
---|
| 2364 | // pairing corrections |
---|
| 2365 | if (bs > 1.0) { |
---|
| 2366 | delta0 = 14; |
---|
| 2367 | } |
---|
| 2368 | else { |
---|
| 2369 | delta0 = 12; |
---|
| 2370 | } |
---|
| 2371 | |
---|
| 2372 | if (esous > 1.0e30) { |
---|
| 2373 | (*dens) = 0.0; |
---|
| 2374 | (*temp) = 0.0; |
---|
| 2375 | goto densniv100; |
---|
| 2376 | } |
---|
| 2377 | |
---|
| 2378 | e = ee - esous; |
---|
| 2379 | |
---|
| 2380 | if (e < 0.0) { |
---|
| 2381 | (*dens) = 0.0; |
---|
| 2382 | (*temp) = 0.0; |
---|
| 2383 | goto densniv100; |
---|
| 2384 | } |
---|
| 2385 | |
---|
| 2386 | // shell corrections |
---|
| 2387 | if (optshp > 0) { |
---|
| 2388 | deltau = bshell; |
---|
| 2389 | if (optshp == 2) { |
---|
| 2390 | deltau = 0.0; |
---|
| 2391 | } |
---|
| 2392 | if (optshp >= 2) { |
---|
| 2393 | // pairing energy shift with condensation energy a.r.j. 10.03.97 |
---|
| 2394 | // deltpp = -0.25e0* (delta0/std::pow(std::sqrt(a),2)) * pa /pi6 + 2.e0*delta0/std::sqrt(a); |
---|
| 2395 | deltpp = -0.25e0* std::pow((delta0/std::sqrt(a)),2) * pa /pi6 + 2.e0*delta0/std::sqrt(a); |
---|
| 2396 | |
---|
| 2397 | parite(a,¶); |
---|
| 2398 | if (para < 0.0) { |
---|
| 2399 | e = e - delta0/std::sqrt(a); |
---|
| 2400 | } else { |
---|
| 2401 | parite(z,&parz); |
---|
| 2402 | if (parz > 0.e0) { |
---|
| 2403 | e = e - 2.0*delta0/std::sqrt(a); |
---|
| 2404 | } else { |
---|
| 2405 | e = e; |
---|
| 2406 | } |
---|
| 2407 | } |
---|
| 2408 | } else { |
---|
| 2409 | deltpp = 0.0; |
---|
| 2410 | } |
---|
| 2411 | } else { |
---|
| 2412 | deltau = 0.0; |
---|
| 2413 | deltpp = 0.0; |
---|
| 2414 | } |
---|
| 2415 | if (e < 0.0) { |
---|
| 2416 | e = 0.0; |
---|
| 2417 | (*temp) = 0.0; |
---|
| 2418 | } |
---|
| 2419 | |
---|
| 2420 | // washing out is made stronger ! g.k. 3.7.96 |
---|
| 2421 | ponfe = -2.5*pa*e*std::pow(a,(-4.0/3.0)); |
---|
| 2422 | |
---|
| 2423 | if (ponfe < -700.0) { |
---|
| 2424 | ponfe = -700.0; |
---|
| 2425 | } |
---|
| 2426 | fe = 1.0 - std::exp(ponfe); |
---|
| 2427 | if (e < ecr) { |
---|
| 2428 | // priv. comm. k.-h. schmidt |
---|
| 2429 | he = 1.0 - std::pow((1.0 - e/ecr),2); |
---|
| 2430 | } |
---|
| 2431 | else { |
---|
| 2432 | he = 1.0; |
---|
| 2433 | } |
---|
| 2434 | |
---|
| 2435 | // Excitation energy corrected for pairing and shell effects |
---|
| 2436 | // washing out with excitation energy is included. |
---|
| 2437 | ecor = e + deltau*fe + deltpp*he; |
---|
| 2438 | |
---|
| 2439 | if (ecor <= 0.1) { |
---|
| 2440 | ecor = 0.1; |
---|
| 2441 | } |
---|
| 2442 | |
---|
| 2443 | // statt 170.d0 a.r.j. 8.11.97 |
---|
| 2444 | |
---|
| 2445 | // iterative procedure according to grossjean and feldmeier |
---|
| 2446 | // to avoid the singularity e = 0 |
---|
| 2447 | if (ee < 5.0) { |
---|
| 2448 | y1 = std::sqrt(pa*ecor); |
---|
| 2449 | for(int j = 0; j < 5; j++) { |
---|
| 2450 | y2 = pa*ecor*(1.e0-std::exp(-y1)); |
---|
| 2451 | y1 = std::sqrt(y2); |
---|
| 2452 | } |
---|
| 2453 | |
---|
| 2454 | y0 = pa/y1; |
---|
| 2455 | (*temp)=1.0/y0; |
---|
| 2456 | (*dens) = std::exp(y0*ecor)/ (std::pow((std::pow(ecor,3)*y0),0.5)*std::pow((1.0-0.5*y0*ecor*std::exp(-y1)),0.5))* std::exp(y1)*(1.0-std::exp(-y1))*0.1477045; |
---|
| 2457 | if (ecor < 1.0) { |
---|
| 2458 | ecor1=1.0; |
---|
| 2459 | y11 = std::sqrt(pa*ecor1); |
---|
| 2460 | for(int j = 0; j < 7; j++) { |
---|
| 2461 | y21 = pa*ecor1*(1.0-std::exp(-y11)); |
---|
| 2462 | y11 = std::sqrt(y21); |
---|
| 2463 | } |
---|
| 2464 | |
---|
| 2465 | y01 = pa/y11; |
---|
| 2466 | (*dens) = (*dens)*std::pow((y01/y0),1.5); |
---|
| 2467 | (*temp) = (*temp)*std::pow((y01/y0),1.5); |
---|
| 2468 | } |
---|
| 2469 | } |
---|
| 2470 | else { |
---|
| 2471 | ponniv = 2.0*std::sqrt(pa*ecor); |
---|
| 2472 | if (ponniv > 700.0) { |
---|
| 2473 | ponniv = 700.0; |
---|
| 2474 | } |
---|
| 2475 | |
---|
| 2476 | // fermi gas state density |
---|
| 2477 | (*dens) = std::pow(pa,(-0.25e0))*std::pow(ecor,(-1.25e0))*std::exp(ponniv) * 0.1477045e0; |
---|
| 2478 | (*temp) = std::sqrt(ecor/pa); |
---|
| 2479 | } |
---|
| 2480 | densniv100: |
---|
| 2481 | |
---|
| 2482 | // spin cutoff parameter |
---|
| 2483 | sig = fp * (*temp); |
---|
| 2484 | |
---|
| 2485 | // collective enhancement |
---|
| 2486 | if (optcol == 1) { |
---|
| 2487 | qrot(z,a,defbet,sig,ecor,&qr); |
---|
| 2488 | } |
---|
| 2489 | else { |
---|
| 2490 | qr = 1.0; |
---|
| 2491 | } |
---|
| 2492 | |
---|
| 2493 | (*dens) = (*dens) * qr; |
---|
[962] | 2494 | if(verboseLevel > 2) { |
---|
| 2495 | G4cout <<"PK::: dens = " << (*dens) << G4endl; |
---|
| 2496 | G4cout <<"PK::: AFP, IZ, ECOR, ECOR1 " << afp << " " << iz << " " << ecor << " " << ecor1 << G4endl; |
---|
| 2497 | } |
---|
[819] | 2498 | } |
---|
| 2499 | |
---|
| 2500 | |
---|
| 2501 | G4double G4Abla::bfms67(G4double zms, G4double ams) |
---|
| 2502 | { |
---|
| 2503 | // This subroutine calculates the fission barriers |
---|
| 2504 | // of the liquid-drop model of Myers and Swiatecki (1967). |
---|
| 2505 | // Analytic parameterization of Dahlinger 1982 |
---|
| 2506 | // replaces tables. Barrier heights from Myers and Swiatecki !!! |
---|
| 2507 | |
---|
[962] | 2508 | G4double nms = 0.0, ims = 0.0, ksims = 0.0, xms = 0.0, ums = 0.0; |
---|
[819] | 2509 | |
---|
| 2510 | nms = ams - zms; |
---|
| 2511 | ims = (nms-zms)/ams; |
---|
| 2512 | ksims= 50.15e0 * (1.- 1.78 * std::pow(ims,2)); |
---|
| 2513 | xms = std::pow(zms,2) / (ams * ksims); |
---|
| 2514 | ums = 0.368e0-5.057e0*xms+8.93e0*std::pow(xms,2)-8.71*std::pow(xms,3); |
---|
| 2515 | return(0.7322e0*std::pow(zms,2)/std::pow(ams,(0.333333e0))*std::pow(10.e0,ums)); |
---|
| 2516 | } |
---|
| 2517 | |
---|
| 2518 | void G4Abla::lpoly(G4double x, G4int n, G4double pl[]) |
---|
| 2519 | { |
---|
| 2520 | // THIS SUBROUTINE CALCULATES THE ORDINARY LEGENDRE POLYNOMIALS OF |
---|
| 2521 | // ORDER 0 TO N-1 OF ARGUMENT X AND STORES THEM IN THE VECTOR PL. |
---|
| 2522 | // THEY ARE CALCULATED BY RECURSION RELATION FROM THE FIRST TWO |
---|
| 2523 | // POLYNOMIALS. |
---|
| 2524 | // WRITTEN BY A.J.SIERK LANL T-9 FEBRUARY, 1984 |
---|
| 2525 | |
---|
| 2526 | // NOTE: PL AND X MUST BE DOUBLE PRECISION ON 32-BIT COMPUTERS! |
---|
| 2527 | |
---|
| 2528 | pl[0] = 1.0; |
---|
| 2529 | pl[1] = x; |
---|
| 2530 | |
---|
| 2531 | for(int i = 2; i < n; i++) { |
---|
| 2532 | pl[i] = ((2*double(i+1) - 3.0)*x*pl[i-1] - (double(i+1) - 2.0)*pl[i-2])/(double(i+1)-1.0); |
---|
| 2533 | } |
---|
| 2534 | } |
---|
| 2535 | |
---|
| 2536 | G4double G4Abla::eflmac(G4int ia, G4int iz, G4int flag, G4int optshp) |
---|
| 2537 | { |
---|
| 2538 | // CHANGED TO CALCULATE TOTAL BINDING ENERGY INSTEAD OF MASS EXCESS. |
---|
| 2539 | // SWITCH FOR PAIRING INCLUDED AS WELL. |
---|
| 2540 | // BINDING = EFLMAC(IA,IZ,0,OPTSHP) |
---|
| 2541 | // FORTRAN TRANSCRIPT OF /U/GREWE/LANG/EEX/FRLDM.C |
---|
| 2542 | // A.J. 15.07.96 |
---|
| 2543 | |
---|
| 2544 | // this function will calculate the liquid-drop nuclear mass for spheri |
---|
| 2545 | // configuration according to the preprint NUCLEAR GROUND-STATE |
---|
| 2546 | // MASSES and DEFORMATIONS by P. M"oller et al. from August 16, 1993 p. |
---|
| 2547 | // All constants are taken from this publication for consistency. |
---|
| 2548 | |
---|
| 2549 | // Parameters: |
---|
| 2550 | // a: nuclear mass number |
---|
| 2551 | // z: nuclear charge |
---|
| 2552 | // flag: 0 - return mass excess |
---|
| 2553 | // otherwise - return pairing (= -1/2 dpn + 1/2 (Dp + Dn)) |
---|
| 2554 | |
---|
[962] | 2555 | G4double eflmacResult = 0.0; |
---|
[819] | 2556 | |
---|
[962] | 2557 | G4int in = 0; |
---|
| 2558 | G4double z = 0.0, n = 0.0, a = 0.0, av = 0.0, as = 0.0; |
---|
| 2559 | G4double a0 = 0.0, c1 = 0.0, c4 = 0.0, b1 = 0.0, b3 = 0.0; |
---|
| 2560 | G4double f = 0.0, ca = 0.0, w = 0.0, dp = 0.0, dn = 0.0, dpn = 0.0, efl = 0.0; |
---|
| 2561 | G4double rmac = 0.0, bs = 0.0, h = 0.0, r0 = 0.0, kf = 0.0, ks = 0.0; |
---|
| 2562 | G4double kv = 0.0, rp = 0.0, ay = 0.0, aden = 0.0, x0 = 0.0, y0 = 0.0; |
---|
| 2563 | G4double mh = 0.0, mn = 0.0, esq = 0.0, ael = 0.0, i = 0.0; |
---|
| 2564 | G4double pi = 3.141592653589793238e0; |
---|
[819] | 2565 | |
---|
| 2566 | // fundamental constants |
---|
| 2567 | // hydrogen-atom mass excess |
---|
| 2568 | mh = 7.289034; |
---|
| 2569 | |
---|
| 2570 | // neutron mass excess |
---|
| 2571 | mn = 8.071431; |
---|
| 2572 | |
---|
| 2573 | // electronic charge squared |
---|
| 2574 | esq = 1.4399764; |
---|
| 2575 | |
---|
| 2576 | // constants from considerations other than nucl. masses |
---|
| 2577 | // electronic binding |
---|
| 2578 | ael = 1.433e-5; |
---|
| 2579 | |
---|
| 2580 | // proton rms radius |
---|
| 2581 | rp = 0.8; |
---|
| 2582 | |
---|
| 2583 | // nuclear radius constant |
---|
| 2584 | r0 = 1.16; |
---|
| 2585 | |
---|
| 2586 | // range of yukawa-plus-expon. potential |
---|
| 2587 | ay = 0.68; |
---|
| 2588 | |
---|
| 2589 | // range of yukawa function used to generate |
---|
| 2590 | // nuclear charge distribution |
---|
| 2591 | aden= 0.70; |
---|
| 2592 | |
---|
| 2593 | // constants from considering odd-even mass differences |
---|
| 2594 | // average pairing gap |
---|
| 2595 | rmac= 4.80; |
---|
| 2596 | |
---|
| 2597 | // neutron-proton interaction |
---|
| 2598 | h = 6.6; |
---|
| 2599 | |
---|
| 2600 | // wigner constant |
---|
| 2601 | w = 30.0; |
---|
| 2602 | |
---|
| 2603 | // adjusted parameters |
---|
| 2604 | // volume energy |
---|
| 2605 | av = 16.00126; |
---|
| 2606 | |
---|
| 2607 | // volume asymmetry |
---|
| 2608 | kv = 1.92240; |
---|
| 2609 | |
---|
| 2610 | // surface energy |
---|
| 2611 | as = 21.18466; |
---|
| 2612 | |
---|
| 2613 | // surface asymmetry |
---|
| 2614 | ks = 2.345; |
---|
| 2615 | // a^0 constant |
---|
| 2616 | a0 = 2.615; |
---|
| 2617 | |
---|
| 2618 | // charge asymmetry |
---|
| 2619 | ca = 0.10289; |
---|
| 2620 | |
---|
| 2621 | // we will account for deformation by using the microscopic |
---|
| 2622 | // corrections tabulated from p. 68ff */ |
---|
| 2623 | bs = 1.0; |
---|
| 2624 | |
---|
| 2625 | z = double(iz); |
---|
| 2626 | a = double(ia); |
---|
| 2627 | in = ia - iz; |
---|
| 2628 | n = double(in); |
---|
| 2629 | dn = rmac*bs/std::pow(n,(1.0/3.0)); |
---|
| 2630 | dp = rmac*bs/std::pow(z,(1.0/3.0)); |
---|
| 2631 | dpn = h/bs/std::pow(a,(2.0/3.0)); |
---|
| 2632 | |
---|
| 2633 | c1 = 3.0/5.0*esq/r0; |
---|
| 2634 | c4 = 5.0/4.0*std::pow((3.0/(2.0*pi)),(2.0/3.0)) * c1; |
---|
| 2635 | kf = std::pow((9.0*pi*z/(4.0*a)),(1.0/3.0))/r0; |
---|
| 2636 | |
---|
| 2637 | f = -1.0/8.0*rp*rp*esq/std::pow(r0,3) * (145.0/48.0 - 327.0/2880.0*std::pow(kf,2) * std::pow(rp,2) + 1527.0/1209600.0*std::pow(kf,4) * std::pow(rp,4)); |
---|
| 2638 | i = (n-z)/a; |
---|
| 2639 | |
---|
| 2640 | x0 = r0 * std::pow(a,(1.0/3.0)) / ay; |
---|
| 2641 | y0 = r0 * std::pow(a,(1.0/3.0)) / aden; |
---|
| 2642 | |
---|
| 2643 | b1 = 1.0 - 3.0/(std::pow(x0,2)) + (1.0 + x0) * (2.0 + 3.0/x0 + 3.0/std::pow(x0,2)) * std::exp(-2.0*x0); |
---|
| 2644 | |
---|
| 2645 | b3 = 1.0 - 5.0/std::pow(y0,2) * (1.0 - 15.0/(8.0*y0) + 21.0/(8.0 * std::pow(y0,3)) |
---|
| 2646 | - 3.0/4.0 * (1.0 + 9.0/(2.0*y0) + 7.0/std::pow(y0,2) |
---|
| 2647 | + 7.0/(2.0 * std::pow(y0,3))) * std::exp(-2.0*y0)); |
---|
| 2648 | |
---|
| 2649 | // now calulation of total binding energy a.j. 16.7.96 |
---|
| 2650 | |
---|
| 2651 | efl = -1.0 * av*(1.0 - kv*i*i)*a + as*(1.0 - ks*i*i)*b1 * std::pow(a,(2.0/3.0)) + a0 |
---|
| 2652 | + c1*z*z*b3/std::pow(a,(1.0/3.0)) - c4*std::pow(z,(4.0/3.0))/std::pow(a,(1.e0/3.e0)) |
---|
| 2653 | + f*std::pow(z,2)/a -ca*(n-z) - ael * std::pow(z,(2.39e0)); |
---|
| 2654 | |
---|
| 2655 | if ((in == iz) && (mod(in,2) == 1) && (mod(iz,2) == 1)) { |
---|
| 2656 | // n and z odd and equal |
---|
| 2657 | efl = efl + w*(utilabs(i)+1.e0/a); |
---|
| 2658 | } |
---|
| 2659 | else { |
---|
| 2660 | efl= efl + w* utilabs(i); |
---|
| 2661 | } |
---|
| 2662 | |
---|
| 2663 | // pairing is made optional |
---|
| 2664 | if (optshp >= 2) { |
---|
| 2665 | // average pairing |
---|
| 2666 | if ((mod(in,2) == 1) && (mod(iz,2) == 1)) { |
---|
| 2667 | efl = efl - dpn; |
---|
| 2668 | } |
---|
| 2669 | if (mod(in,2) == 1) { |
---|
| 2670 | efl = efl + dn; |
---|
| 2671 | } |
---|
| 2672 | if (mod(iz,2) == 1) { |
---|
| 2673 | efl = efl + dp; |
---|
| 2674 | } |
---|
| 2675 | // end if for pairing term |
---|
| 2676 | } |
---|
| 2677 | |
---|
| 2678 | if (flag != 0) { |
---|
| 2679 | eflmacResult = (0.5*(dn + dp) - 0.5*dpn); |
---|
| 2680 | } |
---|
| 2681 | else { |
---|
| 2682 | eflmacResult = efl; |
---|
| 2683 | } |
---|
| 2684 | |
---|
| 2685 | return eflmacResult; |
---|
| 2686 | } |
---|
| 2687 | |
---|
| 2688 | void G4Abla::appariem(G4double a, G4double z, G4double *del) |
---|
| 2689 | { |
---|
| 2690 | // CALCUL DE LA CORRECTION, DUE A L'APPARIEMENT, DE L'ENERGIE DE |
---|
| 2691 | // LIAISON D'UN NOYAU |
---|
| 2692 | // PROCEDURE FOR CALCULATING THE PAIRING CORRECTION TO THE BINDING |
---|
| 2693 | // ENERGY OF A SPECIFIC NUCLEUS |
---|
| 2694 | |
---|
[962] | 2695 | double para = 0.0, parz = 0.0; |
---|
[819] | 2696 | // A MASS NUMBER |
---|
| 2697 | // Z NUCLEAR CHARGE |
---|
| 2698 | // PARA HELP VARIABLE FOR PARITY OF A |
---|
| 2699 | // PARZ HELP VARIABLE FOR PARITY OF Z |
---|
| 2700 | // DEL PAIRING CORRECTION |
---|
| 2701 | |
---|
| 2702 | parite(a, ¶); |
---|
| 2703 | |
---|
| 2704 | if (para < 0.0) { |
---|
| 2705 | (*del) = 0.0; |
---|
| 2706 | } |
---|
| 2707 | else { |
---|
| 2708 | parite(z, &parz); |
---|
| 2709 | if (parz > 0.0) { |
---|
| 2710 | (*del) = -12.0/std::sqrt(a); |
---|
| 2711 | } |
---|
| 2712 | else { |
---|
| 2713 | (*del) = 12.0/std::sqrt(a); |
---|
| 2714 | } |
---|
| 2715 | } |
---|
| 2716 | } |
---|
| 2717 | |
---|
| 2718 | void G4Abla::parite(G4double n, G4double *par) |
---|
| 2719 | { |
---|
| 2720 | // CALCUL DE LA PARITE DU NOMBRE N |
---|
| 2721 | // |
---|
| 2722 | // PROCEDURE FOR CALCULATING THE PARITY OF THE NUMBER N. |
---|
| 2723 | // RETURNS -1 IF N IS ODD AND +1 IF N IS EVEN |
---|
| 2724 | |
---|
[962] | 2725 | G4double n1 = 0.0, n2 = 0.0, n3 = 0.0; |
---|
[819] | 2726 | |
---|
| 2727 | // N NUMBER TO BE TESTED |
---|
| 2728 | // N1,N2 HELP VARIABLES |
---|
| 2729 | // PAR HELP VARIABLE FOR PARITY OF N |
---|
| 2730 | |
---|
| 2731 | n3 = double(idnint(n)); |
---|
| 2732 | n1 = n3/2.0; |
---|
| 2733 | n2 = n1 - dint(n1); |
---|
| 2734 | |
---|
| 2735 | if (n2 > 0.0) { |
---|
| 2736 | (*par) = -1.0; |
---|
| 2737 | } |
---|
| 2738 | else { |
---|
| 2739 | (*par) = 1.0; |
---|
| 2740 | } |
---|
| 2741 | } |
---|
| 2742 | |
---|
| 2743 | G4double G4Abla::tau(G4double bet, G4double homega, G4double ef, G4double t) |
---|
| 2744 | { |
---|
| 2745 | // INPUT : BET, HOMEGA, EF, T |
---|
| 2746 | // OUTPUT: TAU - RISE TIME IN WHICH THE FISSION WIDTH HAS REACHED |
---|
| 2747 | // 90 PERCENT OF ITS FINAL VALUE |
---|
| 2748 | // |
---|
| 2749 | // BETA - NUCLEAR VISCOSITY |
---|
| 2750 | // HOMEGA - CURVATURE OF POTENTIAL |
---|
| 2751 | // EF - FISSION BARRIER |
---|
| 2752 | // T - NUCLEAR TEMPERATURE |
---|
| 2753 | |
---|
[962] | 2754 | G4double tauResult = 0.0; |
---|
[819] | 2755 | |
---|
[962] | 2756 | G4double tlim = 8.e0 * ef; |
---|
[819] | 2757 | if (t > tlim) { |
---|
| 2758 | t = tlim; |
---|
| 2759 | } |
---|
| 2760 | |
---|
| 2761 | // modified bj and khs 6.1.2000: |
---|
| 2762 | if (bet/(std::sqrt(2.0)*10.0*(homega/6.582122)) <= 1.0) { |
---|
| 2763 | tauResult = std::log(10.0*ef/t)/(bet*1.0e21); |
---|
| 2764 | } |
---|
| 2765 | else { |
---|
| 2766 | tauResult = std::log(10.0*ef/t)/ (2.0*std::pow((10.0*homega/6.582122),2))*(bet*1.0e-21); |
---|
| 2767 | } //end if |
---|
| 2768 | |
---|
| 2769 | return tauResult; |
---|
| 2770 | } |
---|
| 2771 | |
---|
| 2772 | G4double G4Abla::cram(G4double bet, G4double homega) |
---|
| 2773 | { |
---|
| 2774 | // INPUT : BET, HOMEGA NUCLEAR VISCOSITY + CURVATURE OF POTENTIAL |
---|
| 2775 | // OUTPUT: KRAMERS FAKTOR - REDUCTION OF THE FISSION PROBABILITY |
---|
| 2776 | // INDEPENDENT OF EXCITATION ENERGY |
---|
| 2777 | |
---|
| 2778 | G4double rel = bet/(20.0*homega/6.582122); |
---|
| 2779 | G4double cramResult = std::sqrt(1.0 + std::pow(rel,2)) - rel; |
---|
| 2780 | // limitation introduced 6.1.2000 by khs |
---|
| 2781 | |
---|
| 2782 | if (cramResult > 1.0) { |
---|
| 2783 | cramResult = 1.0; |
---|
| 2784 | } |
---|
| 2785 | |
---|
| 2786 | return cramResult; |
---|
| 2787 | } |
---|
| 2788 | |
---|
| 2789 | G4double G4Abla::bipol(int iflag, G4double y) |
---|
| 2790 | { |
---|
| 2791 | // CALCULATION OF THE SURFACE BS OR CURVATURE BK OF A NUCLEUS |
---|
| 2792 | // RELATIVE TO THE SPHERICAL CONFIGURATION |
---|
| 2793 | // BASED ON MYERS, DROPLET MODEL FOR ARBITRARY SHAPES |
---|
| 2794 | |
---|
| 2795 | // INPUT: IFLAG - 0/1 BK/BS CALCULATION |
---|
| 2796 | // Y - (1 - X) COMPLEMENT OF THE FISSILITY |
---|
| 2797 | |
---|
| 2798 | // LINEAR INTERPOLATION OF BS BK TABLE |
---|
| 2799 | |
---|
[962] | 2800 | int i = 0; |
---|
[819] | 2801 | |
---|
[962] | 2802 | G4double bipolResult = 0.0; |
---|
[819] | 2803 | |
---|
| 2804 | const int bsbkSize = 54; |
---|
| 2805 | |
---|
| 2806 | G4double bk[bsbkSize] = {0.0, 1.00000,1.00087,1.00352,1.00799,1.01433,1.02265,1.03306, |
---|
| 2807 | 1.04576,1.06099,1.07910,1.10056,1.12603,1.15651,1.19348, |
---|
| 2808 | 1.23915,1.29590,1.35951,1.41013,1.44103,1.46026,1.47339, |
---|
| 2809 | 1.48308,1.49068,1.49692,1.50226,1.50694,1.51114,1.51502, |
---|
| 2810 | 1.51864,1.52208,1.52539,1.52861,1.53177,1.53490,1.53803, |
---|
| 2811 | 1.54117,1.54473,1.54762,1.55096,1.55440,1.55798,1.56173, |
---|
| 2812 | 1.56567,1.56980,1.57413,1.57860,1.58301,1.58688,1.58688, |
---|
| 2813 | 1.58688,1.58740,1.58740, 0.0}; //Zeroes at bk[0], and at the end added by PK |
---|
| 2814 | |
---|
| 2815 | G4double bs[bsbkSize] = {0.0, 1.00000,1.00086,1.00338,1.00750,1.01319, |
---|
| 2816 | 1.02044,1.02927,1.03974, |
---|
| 2817 | 1.05195,1.06604,1.08224,1.10085,1.12229,1.14717,1.17623,1.20963, |
---|
| 2818 | 1.24296,1.26532,1.27619,1.28126,1.28362,1.28458,1.28477,1.28450, |
---|
| 2819 | 1.28394,1.28320,1.28235,1.28141,1.28042,1.27941,1.27837,1.27732, |
---|
| 2820 | 1.27627,1.27522,1.27418,1.27314,1.27210,1.27108,1.27006,1.26906, |
---|
| 2821 | 1.26806,1.26707,1.26610,1.26514,1.26418,1.26325,1.26233,1.26147, |
---|
| 2822 | 1.26147,1.26147,1.25992,1.25992, 0.0}; |
---|
| 2823 | |
---|
| 2824 | i = idint(y/(2.0e-02)) + 1; |
---|
| 2825 | |
---|
| 2826 | if(i >= bsbkSize) { |
---|
| 2827 | if(verboseLevel > 2) { |
---|
| 2828 | G4cout <<"G4Abla error: index i = " << i << " is greater than array size permits." << G4endl; |
---|
| 2829 | } |
---|
| 2830 | bipolResult = 0.0; |
---|
| 2831 | } |
---|
| 2832 | else { |
---|
| 2833 | if (iflag == 1) { |
---|
| 2834 | bipolResult = bs[i] + (bs[i+1] - bs[i])/2.0e-02 * (y - 2.0e-02*(i - 1)); |
---|
| 2835 | } |
---|
| 2836 | else { |
---|
| 2837 | bipolResult = bk[i] + (bk[i+1] - bk[i])/2.0e-02 * (y - 2.0e-02*(i - 1)); |
---|
| 2838 | } |
---|
| 2839 | } |
---|
| 2840 | |
---|
| 2841 | return bipolResult; |
---|
| 2842 | } |
---|
| 2843 | |
---|
| 2844 | void G4Abla::barfit(G4int iz, G4int ia, G4int il, G4double *sbfis, G4double *segs, G4double *selmax) |
---|
| 2845 | { |
---|
| 2846 | // 2223 C VERSION FOR 32BIT COMPUTER |
---|
| 2847 | // 2224 C THIS SUBROUTINE RETURNS THE BARRIER HEIGHT BFIS, THE |
---|
| 2848 | // 2225 C GROUND-STATE ENERGY SEGS, IN MEV, AND THE ANGULAR MOMENTUM |
---|
| 2849 | // 2226 C AT WHICH THE FISSION BARRIER DISAPPEARS, LMAX, IN UNITS OF |
---|
| 2850 | // 2227 C H-BAR, WHEN CALLED WITH INTEGER AGUMENTS IZ, THE ATOMIC |
---|
| 2851 | // 2228 C NUMBER, IA, THE ATOMIC MASS NUMBER, AND IL, THE ANGULAR |
---|
| 2852 | // 2229 C MOMENTUM IN UNITS OF H-BAR. (PLANCK'S CONSTANT DIVIDED BY |
---|
| 2853 | // 2230 C 2*PI). |
---|
| 2854 | // 2231 C |
---|
| 2855 | // 2232 C THE FISSION BARRIER FO IL = 0 IS CALCULATED FROM A 7TH |
---|
| 2856 | // 2233 C ORDER FIT IN TWO VARIABLES TO 638 CALCULATED FISSION |
---|
| 2857 | // 2234 C BARRIERS FOR Z VALUES FROM 20 TO 110. THESE 638 BARRIERS ARE |
---|
| 2858 | // 2235 C FIT WITH AN RMS DEVIATION OF 0.10 MEV BY THIS 49-PARAMETER |
---|
| 2859 | // 2236 C FUNCTION. |
---|
| 2860 | // 2237 C IF BARFIT IS CALLED WITH (IZ,IA) VALUES OUTSIDE THE RANGE OF |
---|
| 2861 | // 2238 C THE BARRIER HEIGHT IS SET TO 0.0, AND A MESSAGE IS PRINTED |
---|
| 2862 | // 2239 C ON THE DEFAULT OUTPUT FILE. |
---|
| 2863 | // 2240 C |
---|
| 2864 | // 2241 C FOR IL VALUES NOT EQUAL TO ZERO, THE VALUES OF L AT WHICH |
---|
| 2865 | // 2242 C THE BARRIER IS 80% AND 20% OF THE L=0 VALUE ARE RESPECTIVELY |
---|
| 2866 | // 2243 C FIT TO 20-PARAMETER FUNCTIONS OF Z AND A, OVER A MORE |
---|
| 2867 | // 2244 C RESTRICTED RANGE OF A VALUES, THAN IS THE CASE FOR L = 0. |
---|
| 2868 | // 2245 C THE VALUE OF L WHERE THE BARRIER DISAPPEARS, LMAX IS FIT TO |
---|
| 2869 | // 2246 C A 24-PARAMETER FUNCTION OF Z AND A, WITH THE SAME RANGE OF |
---|
| 2870 | // 2247 C Z AND A VALUES AS L-80 AND L-20. |
---|
| 2871 | // 2248 C ONCE AGAIN, IF AN (IZ,IA) PAIR IS OUTSIDE OF THE RANGE OF |
---|
| 2872 | // 2249 C VALIDITY OF THE FIT, THE BARRIER VALUE IS SET TO 0.0 AND A |
---|
| 2873 | // 2250 C MESSAGE IS PRINTED. THESE THREE VALUES (BFIS(L=0),L-80, AND |
---|
| 2874 | // 2251 C L-20) AND THE CONSTRINTS OF BFIS = 0 AND D(BFIS)/DL = 0 AT |
---|
| 2875 | // 2252 C L = LMAX AND L=0 LEAD TO A FIFTH-ORDER FIT TO BFIS(L) FOR |
---|
| 2876 | // 2253 C L>L-20. THE FIRST THREE CONSTRAINTS LEAD TO A THIRD-ORDER FIT |
---|
| 2877 | // 2254 C FOR THE REGION L < L-20. |
---|
| 2878 | // 2255 C |
---|
| 2879 | // 2256 C THE GROUND STATE ENERGIES ARE CALCULATED FROM A |
---|
| 2880 | // 2257 C 120-PARAMETER FIT IN Z, A, AND L TO 214 GROUND-STATE ENERGIES |
---|
| 2881 | // 2258 C FOR 36 DIFFERENT Z AND A VALUES. |
---|
| 2882 | // 2259 C (THE RANGE OF Z AND A IS THE SAME AS FOR L-80, L-20, AND |
---|
| 2883 | // 2260 C L-MAX) |
---|
| 2884 | // 2261 C |
---|
| 2885 | // 2262 C THE CALCULATED BARRIERS FROM WHICH THE FITS WERE MADE WERE |
---|
| 2886 | // 2263 C CALCULATED IN 1983-1984 BY A. J. SIERK OF LOS ALAMOS |
---|
| 2887 | // 2264 C NATIONAL LABORATORY GROUP T-9, USING YUKAWA-PLUS-EXPONENTIAL |
---|
| 2888 | // 2265 C G4DOUBLE FOLDED NUCLEAR ENERGY, EXACT COULOMB DIFFUSENESS |
---|
| 2889 | // 2266 C CORRECTIONS, AND DIFFUSE-MATTER MOMENTS OF INERTIA. |
---|
| 2890 | // 2267 C THE PARAMETERS OF THE MODEL R-0 = 1.16 FM, AS 21.13 MEV, |
---|
| 2891 | // 2268 C KAPPA-S = 2.3, A = 0.68 FM. |
---|
| 2892 | // 2269 C THE DIFFUSENESS OF THE MATTER AND CHARGE DISTRIBUTIONS USED |
---|
| 2893 | // 2270 C CORRESPONDS TO A SURFACE DIFFUSENESS PARAMETER (DEFINED BY |
---|
| 2894 | // 2271 C MYERS) OF 0.99 FM. THE CALCULATED BARRIERS FOR L = 0 ARE |
---|
| 2895 | // 2272 C ACCURATE TO A LITTLE LESS THAN 0.1 MEV; THE OUTPUT FROM |
---|
| 2896 | // 2273 C THIS SUBROUTINE IS A LITTLE LESS ACCURATE. WORST ERRORS MAY BE |
---|
| 2897 | // 2274 C AS LARGE AS 0.5 MEV; CHARACTERISTIC UNCERTAINY IS IN THE RANGE |
---|
| 2898 | // 2275 C OF 0.1-0.2 MEV. THE RMS DEVIATION OF THE GROUND-STATE FIT |
---|
| 2899 | // 2276 C FROM THE 214 INPUT VALUES IS 0.20 MEV. THE MAXIMUM ERROR |
---|
| 2900 | // 2277 C OCCURS FOR LIGHT NUCLEI IN THE REGION WHERE THE GROUND STATE |
---|
| 2901 | // 2278 C IS PROLATE, AND MAY BE GREATER THAN 1.0 MEV FOR VERY NEUTRON |
---|
| 2902 | // 2279 C DEFICIENT NUCLEI, WITH L NEAR LMAX. FOR MOST NUCLEI LIKELY TO |
---|
| 2903 | // 2280 C BE ENCOUNTERED IN REAL EXPERIMENTS, THE MAXIMUM ERROR IS |
---|
| 2904 | // 2281 C CLOSER TO 0.5 MEV, AGAIN FOR LIGHT NUCLEI AND L NEAR LMAX. |
---|
| 2905 | // 2282 C |
---|
| 2906 | // 2283 C WRITTEN BY A. J. SIERK, LANL T-9 |
---|
| 2907 | // 2284 C VERSION 1.0 FEBRUARY, 1984 |
---|
| 2908 | // 2285 C |
---|
| 2909 | // 2286 C THE FOLLOWING IS NECESSARY FOR 32-BIT MACHINES LIKE DEC VAX, |
---|
| 2910 | // 2287 C IBM, ETC |
---|
| 2911 | |
---|
| 2912 | G4double pa[7],pz[7],pl[10]; |
---|
[962] | 2913 | for(G4int init_i = 0; init_i < 7; init_i++) { |
---|
| 2914 | pa[init_i] = 0.0; |
---|
| 2915 | pz[init_i] = 0.0; |
---|
| 2916 | } |
---|
| 2917 | for(G4int init_i = 0; init_i < 10; init_i++) { |
---|
| 2918 | pl[init_i] = 0.0; |
---|
| 2919 | } |
---|
[819] | 2920 | |
---|
[962] | 2921 | G4double a = 0.0, z = 0.0, amin = 0.0, amax = 0.0, amin2 = 0.0; |
---|
| 2922 | G4double amax2 = 0.0, aa = 0.0, zz = 0.0, bfis = 0.0; |
---|
| 2923 | G4double bfis0 = 0.0, ell = 0.0, el = 0.0, egs = 0.0, el80 = 0.0, el20 = 0.0; |
---|
| 2924 | G4double elmax = 0.0, sel80 = 0.0, sel20 = 0.0, x = 0.0, y = 0.0, q = 0.0, qa = 0.0, qb = 0.0; |
---|
| 2925 | G4double aj = 0.0, ak = 0.0, a1 = 0.0, a2 = 0.0; |
---|
[819] | 2926 | |
---|
[962] | 2927 | G4int i = 0, j = 0, k = 0, m = 0; |
---|
| 2928 | G4int l = 0; |
---|
| 2929 | |
---|
[819] | 2930 | G4double emncof[4][5] = {{-9.01100e+2,-1.40818e+3, 2.77000e+3,-7.06695e+2, 8.89867e+2}, |
---|
| 2931 | {1.35355e+4,-2.03847e+4, 1.09384e+4,-4.86297e+3,-6.18603e+2}, |
---|
| 2932 | {-3.26367e+3, 1.62447e+3, 1.36856e+3, 1.31731e+3, 1.53372e+2}, |
---|
| 2933 | {7.48863e+3,-1.21581e+4, 5.50281e+3,-1.33630e+3, 5.05367e-2}}; |
---|
| 2934 | |
---|
| 2935 | G4double elmcof[4][5] = {{1.84542e+3,-5.64002e+3, 5.66730e+3,-3.15150e+3, 9.54160e+2}, |
---|
| 2936 | {-2.24577e+3, 8.56133e+3,-9.67348e+3, 5.81744e+3,-1.86997e+3}, |
---|
| 2937 | {2.79772e+3,-8.73073e+3, 9.19706e+3,-4.91900e+3, 1.37283e+3}, |
---|
| 2938 | {-3.01866e+1, 1.41161e+3,-2.85919e+3, 2.13016e+3,-6.49072e+2}}; |
---|
| 2939 | |
---|
| 2940 | G4double emxcof[4][6] = {{9.43596e4,-2.241997e5,2.223237e5,-1.324408e5,4.68922e4,-8.83568e3}, |
---|
| 2941 | {-1.655827e5,4.062365e5,-4.236128e5,2.66837e5,-9.93242e4,1.90644e4}, |
---|
| 2942 | {1.705447e5,-4.032e5,3.970312e5,-2.313704e5,7.81147e4,-1.322775e4}, |
---|
| 2943 | {-9.274555e4,2.278093e5,-2.422225e5,1.55431e5,-5.78742e4,9.97505e3}}; |
---|
| 2944 | |
---|
| 2945 | G4double elzcof[7][7] = {{5.11819909e+5,-1.30303186e+6, 1.90119870e+6,-1.20628242e+6, 5.68208488e+5, 5.48346483e+4,-2.45883052e+4}, |
---|
| 2946 | {-1.13269453e+6, 2.97764590e+6,-4.54326326e+6, 3.00464870e+6, -1.44989274e+6,-1.02026610e+5, 6.27959815e+4}, |
---|
| 2947 | {1.37543304e+6,-3.65808988e+6, 5.47798999e+6,-3.78109283e+6, 1.84131765e+6, 1.53669695e+4,-6.96817834e+4}, |
---|
| 2948 | {-8.56559835e+5, 2.48872266e+6,-4.07349128e+6, 3.12835899e+6, -1.62394090e+6, 1.19797378e+5, 4.25737058e+4}, |
---|
| 2949 | {3.28723311e+5,-1.09892175e+6, 2.03997269e+6,-1.77185718e+6, 9.96051545e+5,-1.53305699e+5,-1.12982954e+4}, |
---|
| 2950 | {4.15850238e+4, 7.29653408e+4,-4.93776346e+5, 6.01254680e+5, -4.01308292e+5, 9.65968391e+4,-3.49596027e+3}, |
---|
| 2951 | {-1.82751044e+5, 3.91386300e+5,-3.03639248e+5, 1.15782417e+5, -4.24399280e+3,-6.11477247e+3, 3.66982647e+2}}; |
---|
| 2952 | |
---|
| 2953 | const G4int sizex = 5; |
---|
| 2954 | const G4int sizey = 6; |
---|
| 2955 | const G4int sizez = 4; |
---|
| 2956 | |
---|
| 2957 | G4double egscof[sizey][sizey][sizez]; |
---|
| 2958 | |
---|
| 2959 | G4double egs1[sizey][sizex] = {{1.927813e5, 7.666859e5, 6.628436e5, 1.586504e5,-7.786476e3}, |
---|
| 2960 | {-4.499687e5,-1.784644e6,-1.546968e6,-4.020658e5,-3.929522e3}, |
---|
| 2961 | {4.667741e5, 1.849838e6, 1.641313e6, 5.229787e5, 5.928137e4}, |
---|
| 2962 | {-3.017927e5,-1.206483e6,-1.124685e6,-4.478641e5,-8.682323e4}, |
---|
| 2963 | {1.226517e5, 5.015667e5, 5.032605e5, 2.404477e5, 5.603301e4}, |
---|
| 2964 | {-1.752824e4,-7.411621e4,-7.989019e4,-4.175486e4,-1.024194e4}}; |
---|
| 2965 | |
---|
| 2966 | G4double egs2[sizey][sizex] = {{-6.459162e5,-2.903581e6,-3.048551e6,-1.004411e6,-6.558220e4}, |
---|
| 2967 | {1.469853e6, 6.564615e6, 6.843078e6, 2.280839e6, 1.802023e5}, |
---|
| 2968 | {-1.435116e6,-6.322470e6,-6.531834e6,-2.298744e6,-2.639612e5}, |
---|
| 2969 | {8.665296e5, 3.769159e6, 3.899685e6, 1.520520e6, 2.498728e5}, |
---|
| 2970 | {-3.302885e5,-1.429313e6,-1.512075e6,-6.744828e5,-1.398771e5}, |
---|
| 2971 | {4.958167e4, 2.178202e5, 2.400617e5, 1.167815e5, 2.663901e4}}; |
---|
| 2972 | |
---|
| 2973 | G4double egs3[sizey][sizex] = {{3.117030e5, 1.195474e6, 9.036289e5, 6.876190e4,-6.814556e4}, |
---|
| 2974 | {-7.394913e5,-2.826468e6,-2.152757e6,-2.459553e5, 1.101414e5}, |
---|
| 2975 | {7.918994e5, 3.030439e6, 2.412611e6, 5.228065e5, 8.542465e3}, |
---|
| 2976 | {-5.421004e5,-2.102672e6,-1.813959e6,-6.251700e5,-1.184348e5}, |
---|
| 2977 | {2.370771e5, 9.459043e5, 9.026235e5, 4.116799e5, 1.001348e5}, |
---|
| 2978 | {-4.227664e4,-1.738756e5,-1.795906e5,-9.292141e4,-2.397528e4}}; |
---|
| 2979 | |
---|
| 2980 | G4double egs4[sizey][sizex] = {{-1.072763e5,-5.973532e5,-6.151814e5, 7.371898e4, 1.255490e5}, |
---|
| 2981 | {2.298769e5, 1.265001e6, 1.252798e6,-2.306276e5,-2.845824e5}, |
---|
| 2982 | {-2.093664e5,-1.100874e6,-1.009313e6, 2.705945e5, 2.506562e5}, |
---|
| 2983 | {1.274613e5, 6.190307e5, 5.262822e5,-1.336039e5,-1.115865e5}, |
---|
| 2984 | {-5.715764e4,-2.560989e5,-2.228781e5,-3.222789e3, 1.575670e4}, |
---|
| 2985 | {1.189447e4, 5.161815e4, 4.870290e4, 1.266808e4, 2.069603e3}}; |
---|
| 2986 | |
---|
| 2987 | for(i = 0; i < sizey; i++) { |
---|
| 2988 | for(j = 0; j < sizex; j++) { |
---|
| 2989 | // egscof[i][j][0] = egs1[i][j]; |
---|
| 2990 | // egscof[i][j][1] = egs2[i][j]; |
---|
| 2991 | // egscof[i][j][2] = egs3[i][j]; |
---|
| 2992 | // egscof[i][j][3] = egs4[i][j]; |
---|
| 2993 | egscof[i][j][0] = egs1[i][j]; |
---|
| 2994 | egscof[i][j][1] = egs2[i][j]; |
---|
| 2995 | egscof[i][j][2] = egs3[i][j]; |
---|
| 2996 | egscof[i][j][3] = egs4[i][j]; |
---|
| 2997 | } |
---|
| 2998 | } |
---|
| 2999 | |
---|
| 3000 | // the program starts here |
---|
| 3001 | if (iz < 19 || iz > 111) { |
---|
| 3002 | goto barfit900; |
---|
| 3003 | } |
---|
| 3004 | |
---|
| 3005 | if(iz > 102 && il > 0) { |
---|
| 3006 | goto barfit902; |
---|
| 3007 | } |
---|
| 3008 | |
---|
| 3009 | z=double(iz); |
---|
| 3010 | a=double(ia); |
---|
| 3011 | el=double(il); |
---|
| 3012 | amin= 1.2e0*z + 0.01e0*z*z; |
---|
| 3013 | amax= 5.8e0*z - 0.024e0*z*z; |
---|
| 3014 | |
---|
| 3015 | if(a < amin || a > amax) { |
---|
| 3016 | goto barfit910; |
---|
| 3017 | } |
---|
| 3018 | |
---|
| 3019 | // angul.mom.zero barrier |
---|
| 3020 | aa=2.5e-3*a; |
---|
| 3021 | zz=1.0e-2*z; |
---|
| 3022 | ell=1.0e-2*el; |
---|
| 3023 | bfis0 = 0.0; |
---|
| 3024 | lpoly(zz,7,pz); |
---|
| 3025 | lpoly(aa,7,pa); |
---|
| 3026 | |
---|
| 3027 | for(i = 0; i < 7; i++) { //do 10 i=1,7 |
---|
| 3028 | for(j = 0; j < 7; j++) { //do 10 j=1,7 |
---|
| 3029 | bfis0=bfis0+elzcof[j][i]*pz[i]*pa[j]; |
---|
| 3030 | //bfis0=bfis0+elzcof[i][j]*pz[j]*pa[i]; |
---|
| 3031 | } |
---|
| 3032 | } |
---|
| 3033 | |
---|
| 3034 | bfis=bfis0; |
---|
| 3035 | |
---|
| 3036 | (*sbfis)=bfis; |
---|
| 3037 | egs=0.0; |
---|
| 3038 | (*segs)=egs; |
---|
| 3039 | |
---|
| 3040 | // values of l at which the barrier |
---|
| 3041 | // is 20%(el20) and 80%(el80) of l=0 value |
---|
| 3042 | amin2 = 1.4e0*z + 0.009e0*z*z; |
---|
| 3043 | amax2 = 20.e0 + 3.0e0*z; |
---|
| 3044 | |
---|
| 3045 | if((a < amin2-5.e0 || a > amax2+10.e0) && il > 0) { |
---|
| 3046 | goto barfit920; |
---|
| 3047 | } |
---|
| 3048 | |
---|
| 3049 | lpoly(zz,5,pz); |
---|
| 3050 | lpoly(aa,4,pa); |
---|
| 3051 | el80=0.0; |
---|
| 3052 | el20=0.0; |
---|
| 3053 | elmax=0.0; |
---|
| 3054 | |
---|
| 3055 | for(i = 0; i < 4; i++) { |
---|
| 3056 | for(j = 0; j < 5; j++) { |
---|
| 3057 | // el80 = el80 + elmcof[j][i]*pz[j]*pa[i]; |
---|
| 3058 | // el20 = el20 + emncof[j][i]*pz[j]*pa[i]; |
---|
| 3059 | el80 = el80 + elmcof[i][j]*pz[j]*pa[i]; |
---|
| 3060 | el20 = el20 + emncof[i][j]*pz[j]*pa[i]; |
---|
| 3061 | } |
---|
| 3062 | } |
---|
| 3063 | |
---|
| 3064 | sel80 = el80; |
---|
| 3065 | sel20 = el20; |
---|
| 3066 | |
---|
| 3067 | // value of l (elmax) where barrier disapp. |
---|
| 3068 | lpoly(zz,6,pz); |
---|
| 3069 | lpoly(ell,9,pl); |
---|
| 3070 | |
---|
| 3071 | for(i = 0; i < 4; i++) { //do 30 i= 1,4 |
---|
| 3072 | for(j = 0; j < 6; j++) { //do 30 j=1,6 |
---|
| 3073 | //elmax = elmax + emxcof[j][i]*pz[j]*pa[i]; |
---|
| 3074 | // elmax = elmax + emxcof[j][i]*pz[i]*pa[j]; |
---|
| 3075 | elmax = elmax + emxcof[i][j]*pz[j]*pa[i]; |
---|
| 3076 | } |
---|
| 3077 | } |
---|
| 3078 | |
---|
| 3079 | (*selmax)=elmax; |
---|
| 3080 | |
---|
| 3081 | // value of barrier at ang.mom. l |
---|
| 3082 | if(il < 1){ |
---|
| 3083 | return; |
---|
| 3084 | } |
---|
| 3085 | |
---|
| 3086 | x = sel20/(*selmax); |
---|
| 3087 | y = sel80/(*selmax); |
---|
| 3088 | |
---|
| 3089 | if(el <= sel20) { |
---|
| 3090 | // low l |
---|
| 3091 | q = 0.2/(std::pow(sel20,2)*std::pow(sel80,2)*(sel20-sel80)); |
---|
| 3092 | qa = q*(4.0*std::pow(sel80,3) - std::pow(sel20,3)); |
---|
| 3093 | qb = -q*(4.0*std::pow(sel80,2) - std::pow(sel20,2)); |
---|
| 3094 | bfis = bfis*(1.0 + qa*std::pow(el,2) + qb*std::pow(el,3)); |
---|
| 3095 | } |
---|
| 3096 | else { |
---|
| 3097 | // high l |
---|
| 3098 | aj = (-20.0*std::pow(x,5) + 25.e0*std::pow(x,4) - 4.0)*std::pow((y-1.0),2)*y*y; |
---|
| 3099 | ak = (-20.0*std::pow(y,5) + 25.0*std::pow(y,4) - 1.0) * std::pow((x-1.0),2)*x*x; |
---|
| 3100 | q = 0.2/(std::pow((y-x)*((1.0-x)*(1.0-y)*x*y),2)); |
---|
| 3101 | qa = q*(aj*y - ak*x); |
---|
| 3102 | qb = -q*(aj*(2.0*y + 1.0) - ak*(2.0*x + 1.0)); |
---|
| 3103 | z = el/(*selmax); |
---|
| 3104 | a1 = 4.0*std::pow(z,5) - 5.0*std::pow(z,4) + 1.0; |
---|
| 3105 | a2 = qa*(2.e0*z + 1.e0); |
---|
| 3106 | bfis=bfis*(a1 + (z - 1.e0)*(a2 + qb*z)*z*z*(z - 1.e0)); |
---|
| 3107 | } |
---|
| 3108 | |
---|
| 3109 | if(bfis <= 0.0) { |
---|
| 3110 | bfis=0.0; |
---|
| 3111 | } |
---|
| 3112 | |
---|
| 3113 | if(el > (*selmax)) { |
---|
| 3114 | bfis=0.0; |
---|
| 3115 | } |
---|
| 3116 | (*sbfis)=bfis; |
---|
| 3117 | |
---|
| 3118 | // now calculate rotating ground state energy |
---|
| 3119 | if(el > (*selmax)) { |
---|
| 3120 | return; |
---|
| 3121 | } |
---|
| 3122 | |
---|
| 3123 | for(k = 0; k < 4; k++) { |
---|
| 3124 | for(l = 0; l < 6; l++) { |
---|
| 3125 | for(m = 0; m < 5; m++) { |
---|
| 3126 | //egs = egs + egscof[l][m][k]*pz[l]*pa[k]*pl[2*m-1]; |
---|
| 3127 | egs = egs + egscof[l][m][k]*pz[l]*pa[k]*pl[2*m]; |
---|
| 3128 | // egs = egs + egscof[m][l][k]*pz[l]*pa[k]*pl[2*m-1]; |
---|
| 3129 | } |
---|
| 3130 | } |
---|
| 3131 | } |
---|
| 3132 | |
---|
| 3133 | (*segs)=egs; |
---|
| 3134 | if((*segs) < 0.0) { |
---|
| 3135 | (*segs)=0.0; |
---|
| 3136 | } |
---|
| 3137 | |
---|
| 3138 | return; |
---|
| 3139 | |
---|
| 3140 | barfit900: //continue |
---|
| 3141 | (*sbfis)=0.0; |
---|
| 3142 | // for z<19 sbfis set to 1.0e3 |
---|
| 3143 | if (iz < 19) { |
---|
| 3144 | (*sbfis) = 1.0e3; |
---|
| 3145 | } |
---|
| 3146 | (*segs)=0.0; |
---|
| 3147 | (*selmax)=0.0; |
---|
| 3148 | return; |
---|
| 3149 | |
---|
| 3150 | barfit902: |
---|
| 3151 | (*sbfis)=0.0; |
---|
| 3152 | (*segs)=0.0; |
---|
| 3153 | (*selmax)=0.0; |
---|
| 3154 | return; |
---|
| 3155 | |
---|
| 3156 | barfit910: |
---|
| 3157 | (*sbfis)=0.0; |
---|
| 3158 | (*segs)=0.0; |
---|
| 3159 | (*selmax)=0.0; |
---|
| 3160 | return; |
---|
| 3161 | |
---|
| 3162 | barfit920: |
---|
| 3163 | (*sbfis)=0.0; |
---|
| 3164 | (*segs)=0.0; |
---|
| 3165 | (*selmax)=0.0; |
---|
| 3166 | return; |
---|
| 3167 | } |
---|
| 3168 | |
---|
| 3169 | G4double G4Abla::expohaz(G4int k, G4double T) |
---|
| 3170 | { |
---|
| 3171 | // TIRAGE ALEATOIRE DANS UNE EXPONENTIELLLE : Y=EXP(-X/T) |
---|
| 3172 | |
---|
| 3173 | return (-1.0*T*std::log(haz(k))); |
---|
| 3174 | } |
---|
| 3175 | |
---|
| 3176 | G4double G4Abla::fd(G4double E) |
---|
| 3177 | { |
---|
| 3178 | // DISTRIBUTION DE MAXWELL |
---|
| 3179 | |
---|
| 3180 | return (E*std::exp(-E)); |
---|
| 3181 | } |
---|
| 3182 | |
---|
| 3183 | G4double G4Abla::f(G4double E) |
---|
| 3184 | { |
---|
| 3185 | // FONCTION INTEGRALE DE FD(E) |
---|
| 3186 | return (1.0 - (E + 1.0) * std::exp(-E)); |
---|
| 3187 | } |
---|
| 3188 | |
---|
| 3189 | G4double G4Abla::fmaxhaz(G4double T) |
---|
| 3190 | { |
---|
| 3191 | // tirage aleatoire dans une maxwellienne |
---|
| 3192 | // t : temperature |
---|
| 3193 | // |
---|
| 3194 | // declaration des variables |
---|
| 3195 | // |
---|
| 3196 | |
---|
| 3197 | const int pSize = 101; |
---|
| 3198 | static G4double p[pSize]; |
---|
| 3199 | |
---|
| 3200 | // ial generateur pour le cascade (et les iy pour eviter les correlations) |
---|
[962] | 3201 | static G4int i = 0; |
---|
[819] | 3202 | static G4int itest = 0; |
---|
| 3203 | // programme principal |
---|
| 3204 | |
---|
| 3205 | // calcul des p(i) par approximation de newton |
---|
| 3206 | p[pSize-1] = 8.0; |
---|
| 3207 | G4double x = 0.1; |
---|
[962] | 3208 | G4double x1 = 0.0; |
---|
| 3209 | G4double y = 0.0; |
---|
[819] | 3210 | |
---|
| 3211 | if (itest == 1) { |
---|
| 3212 | goto fmaxhaz120; |
---|
| 3213 | } |
---|
| 3214 | |
---|
| 3215 | for(i = 1; i <= 99; i++) { |
---|
| 3216 | fmaxhaz20: |
---|
| 3217 | x1 = x - (f(x) - double(i)/100.0)/fd(x); |
---|
| 3218 | x = x1; |
---|
| 3219 | if (std::fabs(f(x) - double(i)/100.0) < 1e-5) { |
---|
| 3220 | goto fmaxhaz100; |
---|
| 3221 | } |
---|
| 3222 | goto fmaxhaz20; |
---|
| 3223 | fmaxhaz100: |
---|
| 3224 | p[i] = x; |
---|
| 3225 | } //end do |
---|
| 3226 | |
---|
| 3227 | // itest = 1; |
---|
| 3228 | itest = 0; |
---|
| 3229 | // tirage aleatoire et calcul du x correspondant |
---|
| 3230 | // par regression lineaire |
---|
| 3231 | fmaxhaz120: |
---|
[962] | 3232 | y = randomGenerator->getRandom(); |
---|
| 3233 | // standardRandom(&y, &(hazard->igraine[17])); |
---|
[819] | 3234 | i = nint(y*100); |
---|
| 3235 | |
---|
| 3236 | // 2590 c ici on evite froidement les depassements de tableaux....(a.b. 3/9/99) |
---|
| 3237 | if(i == 0) { |
---|
| 3238 | goto fmaxhaz120; |
---|
| 3239 | } |
---|
| 3240 | |
---|
| 3241 | if (i == 1) { |
---|
| 3242 | x = p[i]*y*100; |
---|
| 3243 | } |
---|
| 3244 | else { |
---|
| 3245 | x = (p[i] - p[i-1])*(y*100 - i) + p[i]; |
---|
| 3246 | } |
---|
| 3247 | |
---|
| 3248 | return(x*T); |
---|
| 3249 | } |
---|
| 3250 | |
---|
| 3251 | G4double G4Abla::pace2(G4double a, G4double z) |
---|
| 3252 | { |
---|
| 3253 | // PACE2 |
---|
| 3254 | // Cette fonction retourne le defaut de masse du noyau A,Z en MeV |
---|
| 3255 | // Révisée pour a, z flottants 25/4/2002 = |
---|
| 3256 | |
---|
[962] | 3257 | G4double pace2 = 0.0; |
---|
[819] | 3258 | |
---|
| 3259 | G4int ii = idint(a+0.5); |
---|
| 3260 | G4int jj = idint(z+0.5); |
---|
| 3261 | |
---|
| 3262 | if(ii <= 0 || jj < 0) { |
---|
| 3263 | pace2=0.; |
---|
| 3264 | return pace2; |
---|
| 3265 | } |
---|
| 3266 | |
---|
| 3267 | if(jj > 300) { |
---|
| 3268 | pace2=0.0; |
---|
| 3269 | } |
---|
| 3270 | else { |
---|
| 3271 | pace2=pace->dm[ii][jj]; |
---|
| 3272 | } |
---|
| 3273 | pace2=pace2/1000.; |
---|
| 3274 | |
---|
| 3275 | if(pace->dm[ii][jj] == 0.) { |
---|
| 3276 | if(ii < 12) { |
---|
| 3277 | pace2=-500.; |
---|
| 3278 | } |
---|
| 3279 | else { |
---|
| 3280 | guet(&a, &z, &pace2); |
---|
| 3281 | pace2=pace2-ii*931.5; |
---|
| 3282 | pace2=pace2/1000.; |
---|
| 3283 | } |
---|
| 3284 | } |
---|
| 3285 | |
---|
| 3286 | return pace2; |
---|
| 3287 | } |
---|
| 3288 | |
---|
| 3289 | void G4Abla::guet(G4double *x_par, G4double *z_par, G4double *find_par) |
---|
| 3290 | { |
---|
| 3291 | // TABLE DE MASSES ET FORMULE DE MASSE TIRE DU PAPIER DE BRACK-GUET |
---|
| 3292 | // Gives the theoritical value for mass excess... |
---|
| 3293 | // Révisée pour x, z flottants 25/4/2002 |
---|
| 3294 | |
---|
| 3295 | //real*8 x,z |
---|
| 3296 | // dimension q(0:50,0:70) |
---|
| 3297 | G4double x = (*x_par); |
---|
| 3298 | G4double z = (*z_par); |
---|
| 3299 | G4double find = (*find_par); |
---|
| 3300 | |
---|
| 3301 | const G4int qrows = 50; |
---|
| 3302 | const G4int qcols = 70; |
---|
| 3303 | G4double q[qrows][qcols]; |
---|
[962] | 3304 | for(G4int init_i = 0; init_i < qrows; init_i++) { |
---|
| 3305 | for(G4int init_j = 0; init_j < qcols; init_j++) { |
---|
| 3306 | q[init_i][init_j] = 0.0; |
---|
| 3307 | } |
---|
| 3308 | } |
---|
[819] | 3309 | |
---|
| 3310 | G4int ix=G4int(std::floor(x+0.5)); |
---|
| 3311 | G4int iz=G4int(std::floor(z+0.5)); |
---|
| 3312 | G4double zz = iz; |
---|
| 3313 | G4double xx = ix; |
---|
| 3314 | find = 0.0; |
---|
| 3315 | G4double avol = 15.776; |
---|
| 3316 | G4double asur = -17.22; |
---|
| 3317 | G4double ac = -10.24; |
---|
| 3318 | G4double azer = 8.0; |
---|
| 3319 | G4double xjj = -30.03; |
---|
| 3320 | G4double qq = -35.4; |
---|
| 3321 | G4double c1 = -0.737; |
---|
| 3322 | G4double c2 = 1.28; |
---|
| 3323 | |
---|
| 3324 | if(ix <= 7) { |
---|
| 3325 | q[0][1]=939.50; |
---|
| 3326 | q[1][1]=938.21; |
---|
| 3327 | q[1][2]=1876.1; |
---|
| 3328 | q[1][3]=2809.39; |
---|
| 3329 | q[2][4]=3728.34; |
---|
| 3330 | q[2][3]=2809.4; |
---|
| 3331 | q[2][5]=4668.8; |
---|
| 3332 | q[2][6]=5606.5; |
---|
| 3333 | q[3][5]=4669.1; |
---|
| 3334 | q[3][6]=5602.9; |
---|
| 3335 | q[3][7]=6535.27; |
---|
| 3336 | q[4][6]=5607.3; |
---|
| 3337 | q[4][7]=6536.1; |
---|
| 3338 | q[5][7]=6548.3; |
---|
| 3339 | find=q[iz][ix]; |
---|
| 3340 | } |
---|
| 3341 | else { |
---|
| 3342 | G4double xneu=xx-zz; |
---|
| 3343 | G4double si=(xneu-zz)/xx; |
---|
| 3344 | G4double x13=std::pow(xx,.333); |
---|
| 3345 | G4double ee1=c1*zz*zz/x13; |
---|
| 3346 | G4double ee2=c2*zz*zz/xx; |
---|
| 3347 | G4double aux=1.+(9.*xjj/4./qq/x13); |
---|
| 3348 | G4double ee3=xjj*xx*si*si/aux; |
---|
| 3349 | G4double ee4=avol*xx+asur*(std::pow(xx,.666))+ac*x13+azer; |
---|
| 3350 | G4double tota = ee1 + ee2 + ee3 + ee4; |
---|
| 3351 | find = 939.55*xneu+938.77*zz - tota; |
---|
| 3352 | } |
---|
| 3353 | |
---|
| 3354 | (*x_par) = x; |
---|
| 3355 | (*z_par) = z; |
---|
| 3356 | (*find_par) = find; |
---|
| 3357 | } |
---|
| 3358 | |
---|
| 3359 | // SUBROUTINE TRANSLAB(GAMREM,ETREM,CSREM,NOPART,NDEC) |
---|
| 3360 | void G4Abla::translab(G4double gamrem, G4double etrem, G4double csrem[4], G4int nopart, G4int ndec) |
---|
| 3361 | { |
---|
| 3362 | // c Ce subroutine transforme dans un repere 1 les impulsions pcv des |
---|
| 3363 | // c particules acv, zcv et de cosinus directeurs xcv, ycv, zcv calculees |
---|
| 3364 | // c dans un repere 2. |
---|
| 3365 | // c La transformation de lorentz est definie par GAMREM (gamma) et |
---|
| 3366 | // c ETREM (eta). La direction du repere 2 dans 1 est donnees par les |
---|
| 3367 | // c cosinus directeurs ALREM,BEREM,GAREM (axe oz du repere 2). |
---|
| 3368 | // c L'axe oy(2) est fixe par le produit vectoriel oz(1)*oz(2). |
---|
| 3369 | // c Le calcul est fait pour les particules de NDEC a iv du common volant. |
---|
| 3370 | // C Resultats dans le NTUPLE (common VAR_NTP) decale de NOPART (cascade). |
---|
| 3371 | |
---|
| 3372 | // REAL*8 GAMREM,ETREM,ER,PLABI(3),PLABF(3),R(3,3) |
---|
| 3373 | // real*8 MASSE,PTRAV2,CSREM(3),UMA,MELEC,EL |
---|
| 3374 | // real*4 acv,zpcv,pcv,xcv,ycv,zcv |
---|
| 3375 | // common/volant/acv(200),zpcv(200),pcv(200),xcv(200), |
---|
| 3376 | // s ycv(200),zcv(200),iv |
---|
| 3377 | |
---|
| 3378 | // parameter (max=250) |
---|
| 3379 | // real*4 EXINI,ENERJ,BIMPACT,PLAB,TETLAB,PHILAB,ESTFIS |
---|
| 3380 | // integer AVV,ZVV,JREMN,KFIS,IZFIS,IAFIS |
---|
| 3381 | // common/VAR_NTP/MASSINI,MZINI,EXINI,MULNCASC,MULNEVAP, |
---|
| 3382 | // +MULNTOT,BIMPACT,JREMN,KFIS,ESTFIS,IZFIS,IAFIS,NTRACK, |
---|
| 3383 | // +ITYPCASC(max),AVV(max),ZVV(max),ENERJ(max),PLAB(max), |
---|
| 3384 | // +TETLAB(max),PHILAB(max) |
---|
| 3385 | |
---|
| 3386 | // DATA UMA,MELEC/931.4942,0.511/ |
---|
| 3387 | |
---|
| 3388 | // C Matrice de rotation dans le labo: |
---|
[962] | 3389 | G4double avv = 0.0, zvv = 0.0, enerj = 0.0, plab = 0.0, tetlab = 0.0, philab = 0.0; |
---|
| 3390 | G4int itypcasc = 0; |
---|
[819] | 3391 | G4double sitet = std::sqrt(std::pow(csrem[1],2)+std::pow(csrem[2],2)); |
---|
[962] | 3392 | G4double cstet = 0.0, siphi = 0.0, csphi = 0.0; |
---|
[819] | 3393 | G4double R[4][4]; |
---|
[962] | 3394 | for(G4int init_i = 0; init_i < 4; init_i++) { |
---|
| 3395 | for(G4int init_j = 0; init_j < 4; init_j++) { |
---|
| 3396 | R[init_i][init_j] = 0.0; |
---|
| 3397 | } |
---|
| 3398 | } |
---|
| 3399 | if(verboseLevel > 1) |
---|
| 3400 | G4cout <<"csrem = " << csrem[1] << ", " << csrem[2] << ", " << csrem[3] << ")" << G4endl; |
---|
[819] | 3401 | if(sitet > 1.0e-6) { //then |
---|
| 3402 | cstet = csrem[3]; |
---|
| 3403 | siphi = csrem[2]/sitet; |
---|
| 3404 | csphi = csrem[1]/sitet; |
---|
| 3405 | |
---|
| 3406 | R[1][1] = cstet*csphi; |
---|
| 3407 | R[1][2] = -siphi; |
---|
| 3408 | R[1][3] = sitet*csphi; |
---|
| 3409 | R[2][1] = cstet*siphi; |
---|
| 3410 | R[2][2] = csphi; |
---|
| 3411 | R[2][3] = sitet*siphi; |
---|
| 3412 | R[3][1] = -sitet; |
---|
| 3413 | R[3][2] = 0.0; |
---|
| 3414 | R[3][3] = cstet; |
---|
| 3415 | } |
---|
| 3416 | else { |
---|
| 3417 | R[1][1] = 1.0; |
---|
| 3418 | R[1][2] = 0.0; |
---|
| 3419 | R[1][3] = 0.0; |
---|
| 3420 | R[2][1] = 0.0; |
---|
| 3421 | R[2][2] = 1.0; |
---|
| 3422 | R[2][3] = 0.0; |
---|
| 3423 | R[3][1] = 0.0; |
---|
| 3424 | R[3][2] = 0.0; |
---|
| 3425 | R[3][3] = 1.0; |
---|
| 3426 | } //endif |
---|
| 3427 | |
---|
| 3428 | G4int intp = 0; |
---|
| 3429 | G4double el = 0.0; |
---|
| 3430 | G4double masse = 0.0; |
---|
| 3431 | G4double er = 0.0; |
---|
| 3432 | G4double plabi[4]; |
---|
| 3433 | G4double ptrav2 = 0.0; |
---|
| 3434 | G4double plabf[4]; |
---|
| 3435 | G4double bidon = 0.0; |
---|
[962] | 3436 | for(G4int init_i = 0; init_i < 4; init_i++) { |
---|
| 3437 | plabi[init_i] = 0.0; |
---|
| 3438 | plabf[init_i] = 0.0; |
---|
| 3439 | } |
---|
| 3440 | ndec = 1; |
---|
[819] | 3441 | for(G4int i = ndec; i <= volant->iv; i++) { //do i=ndec,iv |
---|
| 3442 | intp = i + nopart; |
---|
[962] | 3443 | if(volant->copied[i]) continue; // Avoid double copying |
---|
| 3444 | #ifdef USE_LEGACY_CODE |
---|
[819] | 3445 | varntp->ntrack = varntp->ntrack + 1; |
---|
[962] | 3446 | #endif |
---|
[819] | 3447 | if(nint(volant->acv[i]) == 0 && nint(volant->zpcv[i]) == 0) { |
---|
[962] | 3448 | if(verboseLevel > -1) { |
---|
| 3449 | G4cout << __FILE__ << ":" << __LINE__ << " Error: Particles with A = 0 Z = 0 detected! " << G4endl; |
---|
[819] | 3450 | } |
---|
| 3451 | continue; |
---|
| 3452 | } |
---|
| 3453 | if(varntp->ntrack >= VARNTPSIZE) { |
---|
| 3454 | if(verboseLevel > 2) { |
---|
| 3455 | G4cout <<"Error! Output data structure not big enough!" << G4endl; |
---|
| 3456 | } |
---|
| 3457 | } |
---|
[962] | 3458 | #ifdef USE_LEGACY_CODE |
---|
[819] | 3459 | varntp->avv[intp] = nint(volant->acv[i]); |
---|
| 3460 | varntp->zvv[intp] = nint(volant->zpcv[i]); |
---|
| 3461 | varntp->itypcasc[intp] = 0; |
---|
[962] | 3462 | #else |
---|
| 3463 | avv = nint(volant->acv[i]); |
---|
| 3464 | zvv = nint(volant->zpcv[i]); |
---|
| 3465 | itypcasc = 0; |
---|
| 3466 | #endif |
---|
[819] | 3467 | // transformation de lorentz remnan --> labo: |
---|
[962] | 3468 | #ifdef USE_LEGACY_CODE |
---|
[819] | 3469 | if (varntp->avv[intp] == -1) { //then |
---|
[962] | 3470 | #else |
---|
| 3471 | if(avv == -1) { |
---|
| 3472 | #endif |
---|
[819] | 3473 | masse = 138.00; // cugnon |
---|
| 3474 | // c if (avv(intp).eq.1) masse=938.2796 !cugnon |
---|
| 3475 | // c if (avv(intp).eq.4) masse=3727.42 !ok |
---|
| 3476 | } |
---|
| 3477 | else { |
---|
| 3478 | mglms(double(volant->acv[i]),double(volant->zpcv[i]),0, &el); |
---|
| 3479 | masse = volant->zpcv[i]*938.27 + (volant->acv[i] - volant->zpcv[i])*939.56 + el; |
---|
| 3480 | } //end if |
---|
| 3481 | |
---|
| 3482 | er = std::sqrt(std::pow(volant->pcv[i],2) + std::pow(masse,2)); |
---|
| 3483 | plabi[1] = volant->pcv[i]*(volant->xcv[i]); |
---|
| 3484 | plabi[2] = volant->pcv[i]*(volant->ycv[i]); |
---|
| 3485 | plabi[3] = er*etrem + gamrem*(volant->pcv[i])*(volant->zcv[i]); |
---|
| 3486 | |
---|
| 3487 | ptrav2 = std::pow(plabi[1],2) + std::pow(plabi[2],2) + std::pow(plabi[3],2); |
---|
[962] | 3488 | #ifdef USE_LEGACY_CODE |
---|
[819] | 3489 | varntp->plab[intp] = std::sqrt(ptrav2); |
---|
[962] | 3490 | #else |
---|
| 3491 | plab = std::sqrt(ptrav2); |
---|
| 3492 | #endif |
---|
| 3493 | #ifdef USE_LEGACY_CODE |
---|
| 3494 | if(std::abs(varntp->plab[intp] - 122.009) < 1.0e-3) { |
---|
| 3495 | G4cout <<__FILE__ << ":" << __LINE__ << " Error: varntp->plab["<< intp <<"] = " << varntp->plab[intp] << G4endl; |
---|
| 3496 | |
---|
| 3497 | volant->dump(); |
---|
| 3498 | varntp->dump(); |
---|
| 3499 | # |
---|
| 3500 | G4cout <<"varntp->plab[intp] = " << varntp->plab[intp] << G4endl; |
---|
| 3501 | G4cout <<"ndec (starting index for loop) = " << ndec << G4endl; |
---|
| 3502 | G4cout <<"volant->iv (stopping index for loop) = " << volant->iv << G4endl; |
---|
| 3503 | G4cout <<"i (current position) = " << i << G4endl; |
---|
| 3504 | G4cout <<"intp (index for writing to varntp) = " << intp << G4endl; |
---|
| 3505 | // exit(0); |
---|
| 3506 | } |
---|
| 3507 | #endif |
---|
| 3508 | |
---|
| 3509 | #ifdef USE_LEGACY_CODE |
---|
[819] | 3510 | varntp->enerj[intp] = std::sqrt(ptrav2 + std::pow(masse,2)) - masse; |
---|
[962] | 3511 | #else |
---|
| 3512 | enerj = std::sqrt(ptrav2 + std::pow(masse,2)) - masse; |
---|
| 3513 | #endif |
---|
[819] | 3514 | // Rotation dans le labo: |
---|
| 3515 | for(G4int j = 1; j <= 3; j++) { //do j=1,3 |
---|
| 3516 | plabf[j] = 0.0; |
---|
| 3517 | for(G4int k = 1; k <= 3; k++) { //do k=1,3 |
---|
[962] | 3518 | //plabf[j] = plabf[j] + R[k][j]*plabi[k]; // :::Fixme::: (indices?) |
---|
| 3519 | plabf[j] = plabf[j] + R[j][k]*plabi[k]; // :::Fixme::: (indices?) |
---|
[819] | 3520 | } // end do |
---|
| 3521 | } // end do |
---|
| 3522 | // C impulsions dans le nouveau systeme copiees dans /volant/ |
---|
[962] | 3523 | #ifdef USE_LEGACY_CODE |
---|
[819] | 3524 | volant->pcv[i] = varntp->plab[intp]; |
---|
[962] | 3525 | #else |
---|
| 3526 | volant->pcv[i] = plab; |
---|
| 3527 | #endif |
---|
[819] | 3528 | ptrav2 = std::sqrt(std::pow(plabf[1],2) + std::pow(plabf[2],2) + std::pow(plabf[3],2)); |
---|
| 3529 | if(ptrav2 >= 1.0e-6) { //then |
---|
| 3530 | volant->xcv[i] = plabf[1]/ptrav2; |
---|
| 3531 | volant->ycv[i] = plabf[2]/ptrav2; |
---|
| 3532 | volant->zcv[i] = plabf[3]/ptrav2; |
---|
| 3533 | } |
---|
| 3534 | else { |
---|
| 3535 | volant->xcv[i] = 1.0; |
---|
| 3536 | volant->ycv[i] = 0.0; |
---|
| 3537 | volant->zcv[i] = 0.0; |
---|
| 3538 | } //endif |
---|
| 3539 | // impulsions dans le nouveau systeme copiees dans /VAR_NTP/ |
---|
[962] | 3540 | #ifdef USE_LEGACY_CODE |
---|
[819] | 3541 | if(varntp->plab[intp] >= 1.0e-6) { //then |
---|
[962] | 3542 | #else |
---|
| 3543 | if(plab >= 1.0e-6) { //then |
---|
| 3544 | #endif |
---|
| 3545 | #ifdef USE_LEGACY_CODE |
---|
[819] | 3546 | bidon = plabf[3]/(varntp->plab[intp]); |
---|
[962] | 3547 | #else |
---|
| 3548 | bidon = plabf[3]/plab; |
---|
| 3549 | #endif |
---|
[819] | 3550 | if(bidon > 1.0) { |
---|
| 3551 | bidon = 1.0; |
---|
| 3552 | } |
---|
| 3553 | if(bidon < -1.0) { |
---|
| 3554 | bidon = -1.0; |
---|
| 3555 | } |
---|
[962] | 3556 | #ifdef USE_LEGACY_CODE |
---|
[819] | 3557 | varntp->tetlab[intp] = std::acos(bidon); |
---|
| 3558 | sitet = std::sin(varntp->tetlab[intp]); |
---|
| 3559 | varntp->philab[intp] = std::atan2(plabf[2],plabf[1]); |
---|
| 3560 | varntp->tetlab[intp] = varntp->tetlab[intp]*57.2957795; |
---|
| 3561 | varntp->philab[intp] = varntp->philab[intp]*57.2957795; |
---|
[962] | 3562 | #else |
---|
| 3563 | tetlab = std::acos(bidon); |
---|
| 3564 | sitet = std::sin(tetlab); |
---|
| 3565 | philab = std::atan2(plabf[2],plabf[1]); |
---|
| 3566 | tetlab = tetlab*57.2957795; |
---|
| 3567 | philab = philab*57.2957795; |
---|
| 3568 | #endif |
---|
[819] | 3569 | } |
---|
| 3570 | else { |
---|
[962] | 3571 | #ifdef USE_LEGACY_CODE |
---|
[819] | 3572 | varntp->tetlab[intp] = 90.0; |
---|
| 3573 | varntp->philab[intp] = 0.0; |
---|
[962] | 3574 | #else |
---|
| 3575 | tetlab = 90.0; |
---|
| 3576 | philab = 0.0; |
---|
| 3577 | #endif |
---|
[819] | 3578 | } // endif |
---|
[962] | 3579 | volant->copied[i] = true; |
---|
| 3580 | #ifndef USE_LEGACY_CODE |
---|
| 3581 | varntp->addParticle(avv, zvv, enerj, plab, tetlab, philab); |
---|
| 3582 | #else |
---|
| 3583 | G4cout <<__FILE__ << ":" << __LINE__ << " volant -> varntp: " << " intp = " << intp << "Avv = " << varntp->avv[intp] << " Zvv = " << varntp->zvv[intp] << " Plab = " << varntp->plab[intp] << G4endl; |
---|
| 3584 | G4cout <<__FILE__ << ":" << __LINE__ << " volant index i = " << i << " varnpt index intp = " << intp << G4endl; |
---|
| 3585 | #endif |
---|
[819] | 3586 | } // end do |
---|
| 3587 | } |
---|
| 3588 | // C------------------------------------------------------------------------- |
---|
| 3589 | |
---|
| 3590 | // SUBROUTINE TRANSLABPF(MASSE1,T1,P1,CTET1,PHI1,GAMREM,ETREM,R, |
---|
| 3591 | // s PLAB1,GAM1,ETA1,CSDIR) |
---|
| 3592 | void G4Abla::translabpf(G4double masse1, G4double t1, G4double p1, G4double ctet1, |
---|
| 3593 | G4double phi1, G4double gamrem, G4double etrem, G4double R[][4], |
---|
| 3594 | G4double *plab1, G4double *gam1, G4double *eta1, G4double csdir[]) |
---|
| 3595 | { |
---|
| 3596 | // C Calcul de l'impulsion du PF (PLAB1, cos directeurs CSDIR(3)) dans le |
---|
| 3597 | // C systeme remnant et des coefs de Lorentz GAM1,ETA1 de passage |
---|
| 3598 | // c du systeme PF --> systeme remnant. |
---|
| 3599 | // c |
---|
| 3600 | // C Input: MASSE1, T1 (energie cinetique), CTET1,PHI1 (cosTHETA et PHI) |
---|
| 3601 | // C (le PF dans le systeme du Noyau de Fission (NF)). |
---|
| 3602 | // C GAMREM,ETREM les coefs de Lorentz systeme NF --> syst remnant, |
---|
| 3603 | // C R(3,3) la matrice de rotation systeme NF--> systeme remnant. |
---|
| 3604 | // C |
---|
| 3605 | // C |
---|
| 3606 | // REAL*8 MASSE1,T1,P1,CTET1,PHI1,GAMREM,ETREM,R(3,3), |
---|
| 3607 | // s PLAB1,GAM1,ETA1,CSDIR(3),ER,SITET,PLABI(3),PLABF(3) |
---|
| 3608 | |
---|
| 3609 | G4double er = t1 + masse1; |
---|
| 3610 | |
---|
| 3611 | G4double sitet = std::sqrt(1.0 - std::pow(ctet1,2)); |
---|
| 3612 | |
---|
| 3613 | G4double plabi[4]; |
---|
| 3614 | G4double plabf[4]; |
---|
[962] | 3615 | for(G4int init_i = 0; init_i < 4; init_i++) { |
---|
| 3616 | plabi[init_i] = 0.0; |
---|
| 3617 | plabf[init_i] = 0.0; |
---|
| 3618 | } |
---|
| 3619 | |
---|
[819] | 3620 | // C ----Transformation de Lorentz Noyau fissionnant --> Remnant: |
---|
| 3621 | plabi[1] = p1*sitet*std::cos(phi1); |
---|
| 3622 | plabi[2] = p1*sitet*std::sin(phi1); |
---|
| 3623 | plabi[3] = er*etrem + gamrem*p1*ctet1; |
---|
| 3624 | |
---|
| 3625 | // C ----Rotation du syst Noyaut Fissionant vers syst remnant: |
---|
| 3626 | for(G4int j = 1; j <= 3; j++) { // do j=1,3 |
---|
| 3627 | plabf[j] = 0.0; |
---|
| 3628 | for(G4int k = 1; k <= 3; k++) { //do k=1,3 |
---|
[962] | 3629 | plabf[j] = plabf[j] + R[j][k]*plabi[k]; |
---|
| 3630 | //plabf[j] = plabf[j] + R[k][j]*plabi[k]; |
---|
[819] | 3631 | } //end do |
---|
| 3632 | } //end do |
---|
| 3633 | // C ----Cosinus directeurs et coefs de la transf de Lorentz dans le |
---|
| 3634 | // c nouveau systeme: |
---|
| 3635 | (*plab1) = std::pow(plabf[1],2) + std::pow(plabf[2],2) + std::pow(plabf[3],2); |
---|
| 3636 | (*gam1) = std::sqrt(std::pow(masse1,2) + (*plab1))/masse1; |
---|
| 3637 | (*plab1) = std::sqrt((*plab1)); |
---|
| 3638 | (*eta1) = (*plab1)/masse1; |
---|
| 3639 | |
---|
| 3640 | if((*plab1) <= 1.0e-6) { //then |
---|
| 3641 | csdir[1] = 0.0; |
---|
| 3642 | csdir[2] = 0.0; |
---|
| 3643 | csdir[3] = 1.0; |
---|
| 3644 | } |
---|
| 3645 | else { |
---|
| 3646 | for(G4int i = 1; i <= 3; i++) { //do i=1,3 |
---|
| 3647 | csdir[i] = plabf[i]/(*plab1); |
---|
| 3648 | } // end do |
---|
| 3649 | } //endif |
---|
| 3650 | } |
---|
| 3651 | |
---|
| 3652 | // SUBROUTINE LOR_AB(GAM,ETA,Ein,Pin,Eout,Pout) |
---|
| 3653 | void G4Abla::lorab(G4double gam, G4double eta, G4double ein, G4double pin[], |
---|
| 3654 | G4double *eout, G4double pout[]) |
---|
| 3655 | { |
---|
| 3656 | // C Transformation de lorentz brute pour vérifs. |
---|
| 3657 | // C P(3) = P_longitudinal (transformé) |
---|
| 3658 | // C P(1) et P(2) = P_transvers (non transformés) |
---|
| 3659 | // DIMENSION Pin(3),Pout(3) |
---|
| 3660 | // REAL*8 GAM,ETA,Ein |
---|
| 3661 | |
---|
| 3662 | pout[1] = pin[1]; |
---|
| 3663 | pout[2] = pin[2]; |
---|
| 3664 | (*eout) = gam*ein + eta*pin[3]; |
---|
| 3665 | pout[3] = eta*ein + gam*pin[3]; |
---|
| 3666 | } |
---|
| 3667 | |
---|
| 3668 | // SUBROUTINE ROT_AB(R,Pin,Pout) |
---|
| 3669 | void G4Abla::rotab(G4double R[4][4], G4double pin[4], G4double pout[4]) |
---|
| 3670 | { |
---|
| 3671 | // C Rotation d'un vecteur |
---|
| 3672 | // DIMENSION Pin(3),Pout(3) |
---|
| 3673 | // REAL*8 R(3,3) |
---|
| 3674 | |
---|
| 3675 | for(G4int i = 1; i <= 3; i++) { // do i=1,3 |
---|
| 3676 | pout[i] = 0.0; |
---|
| 3677 | for(G4int j = 1; j <= 3; j++) { //do j=1,3 |
---|
[962] | 3678 | pout[i] = pout[i] + R[i][j]*pin[j]; |
---|
| 3679 | //pout[i] = pout[i] + R[j][i]*pin[j]; |
---|
[819] | 3680 | } // enddo |
---|
| 3681 | } //enddo |
---|
| 3682 | } |
---|
| 3683 | |
---|
| 3684 | // Methods related to the internal ABLA random number generator. In |
---|
| 3685 | // the future the random number generation must be factored into its |
---|
| 3686 | // own class |
---|
| 3687 | |
---|
| 3688 | void G4Abla::standardRandom(G4double *rndm, G4long *seed) |
---|
| 3689 | { |
---|
| 3690 | (*seed) = (*seed); // Avoid warning during compilation. |
---|
| 3691 | // Use Geant4 G4UniformRand |
---|
[962] | 3692 | (*rndm) = randomGenerator->getRandom(); |
---|
[819] | 3693 | } |
---|
| 3694 | |
---|
| 3695 | G4double G4Abla::haz(G4int k) |
---|
| 3696 | { |
---|
| 3697 | const G4int pSize = 110; |
---|
| 3698 | static G4double p[pSize]; |
---|
[962] | 3699 | static G4long ix = 0, i = 0; |
---|
| 3700 | static G4double x = 0.0, y = 0.0, a = 0.0, haz = 0.0; |
---|
[819] | 3701 | // k =< -1 on initialise |
---|
| 3702 | // k = -1 c'est reproductible |
---|
| 3703 | // k < -1 || k > -1 ce n'est pas reproductible |
---|
| 3704 | |
---|
| 3705 | // Zero is invalid random seed. Set proper value from our random seed collection: |
---|
| 3706 | if(ix == 0) { |
---|
| 3707 | ix = hazard->ial; |
---|
| 3708 | } |
---|
| 3709 | |
---|
| 3710 | if (k <= -1) { //then |
---|
| 3711 | if(k == -1) { //then |
---|
| 3712 | ix = 0; |
---|
| 3713 | } |
---|
| 3714 | else { |
---|
| 3715 | x = 0.0; |
---|
| 3716 | y = secnds(int(x)); |
---|
| 3717 | ix = int(y * 100 + 43543000); |
---|
| 3718 | if(mod(ix,2) == 0) { |
---|
| 3719 | ix = ix + 1; |
---|
| 3720 | } |
---|
| 3721 | } |
---|
| 3722 | |
---|
| 3723 | // Here we are using random number generator copied from INCL code |
---|
| 3724 | // instead of the CERNLIB one! This causes difficulties for |
---|
| 3725 | // automatic testing since the random number generators, and thus |
---|
| 3726 | // the behavior of the routines in C++ and FORTRAN versions is no |
---|
| 3727 | // longer exactly the same! |
---|
[962] | 3728 | x = randomGenerator->getRandom(); |
---|
| 3729 | // standardRandom(&x, &ix); |
---|
[819] | 3730 | for(G4int i = 0; i < pSize; i++) { //do i=1,110 |
---|
[962] | 3731 | p[i] = randomGenerator->getRandom(); |
---|
| 3732 | // standardRandom(&(p[i]), &ix); |
---|
[819] | 3733 | } |
---|
[962] | 3734 | a = randomGenerator->getRandom(); |
---|
[819] | 3735 | standardRandom(&a, &ix); |
---|
| 3736 | k = 0; |
---|
| 3737 | } |
---|
| 3738 | |
---|
| 3739 | i = nint(100*a)+1; |
---|
| 3740 | haz = p[i]; |
---|
[962] | 3741 | a = randomGenerator->getRandom(); |
---|
| 3742 | // standardRandom(&a, &ix); |
---|
[819] | 3743 | p[i] = a; |
---|
| 3744 | |
---|
| 3745 | hazard->ial = ix; |
---|
| 3746 | return haz; |
---|
| 3747 | } |
---|
| 3748 | |
---|
| 3749 | // Utilities |
---|
| 3750 | |
---|
| 3751 | G4double G4Abla::min(G4double a, G4double b) |
---|
| 3752 | { |
---|
| 3753 | if(a < b) { |
---|
| 3754 | return a; |
---|
| 3755 | } |
---|
| 3756 | else { |
---|
| 3757 | return b; |
---|
| 3758 | } |
---|
| 3759 | } |
---|
| 3760 | |
---|
| 3761 | G4int G4Abla::min(G4int a, G4int b) |
---|
| 3762 | { |
---|
| 3763 | if(a < b) { |
---|
| 3764 | return a; |
---|
| 3765 | } |
---|
| 3766 | else { |
---|
| 3767 | return b; |
---|
| 3768 | } |
---|
| 3769 | } |
---|
| 3770 | |
---|
| 3771 | G4double G4Abla::max(G4double a, G4double b) |
---|
| 3772 | { |
---|
| 3773 | if(a > b) { |
---|
| 3774 | return a; |
---|
| 3775 | } |
---|
| 3776 | else { |
---|
| 3777 | return b; |
---|
| 3778 | } |
---|
| 3779 | } |
---|
| 3780 | |
---|
| 3781 | G4int G4Abla::max(G4int a, G4int b) |
---|
| 3782 | { |
---|
| 3783 | if(a > b) { |
---|
| 3784 | return a; |
---|
| 3785 | } |
---|
| 3786 | else { |
---|
| 3787 | return b; |
---|
| 3788 | } |
---|
| 3789 | } |
---|
| 3790 | |
---|
| 3791 | G4int G4Abla::nint(G4double number) |
---|
| 3792 | { |
---|
[962] | 3793 | G4double intpart = 0.0; |
---|
| 3794 | G4double fractpart = 0.0; |
---|
[819] | 3795 | fractpart = std::modf(number, &intpart); |
---|
| 3796 | if(number == 0) { |
---|
| 3797 | return 0; |
---|
| 3798 | } |
---|
| 3799 | if(number > 0) { |
---|
| 3800 | if(fractpart < 0.5) { |
---|
| 3801 | return int(std::floor(number)); |
---|
| 3802 | } |
---|
| 3803 | else { |
---|
| 3804 | return int(std::ceil(number)); |
---|
| 3805 | } |
---|
| 3806 | } |
---|
| 3807 | if(number < 0) { |
---|
| 3808 | if(fractpart < -0.5) { |
---|
| 3809 | return int(std::floor(number)); |
---|
| 3810 | } |
---|
| 3811 | else { |
---|
| 3812 | return int(std::ceil(number)); |
---|
| 3813 | } |
---|
| 3814 | } |
---|
| 3815 | |
---|
| 3816 | return int(std::floor(number)); |
---|
| 3817 | } |
---|
| 3818 | |
---|
| 3819 | G4int G4Abla::secnds(G4int x) |
---|
| 3820 | { |
---|
| 3821 | time_t mytime; |
---|
| 3822 | tm *mylocaltime; |
---|
| 3823 | |
---|
| 3824 | time(&mytime); |
---|
| 3825 | mylocaltime = localtime(&mytime); |
---|
| 3826 | |
---|
| 3827 | if(x == 0) { |
---|
| 3828 | return(mylocaltime->tm_hour*60*60 + mylocaltime->tm_min*60 + mylocaltime->tm_sec); |
---|
| 3829 | } |
---|
| 3830 | else { |
---|
| 3831 | return(mytime - x); |
---|
| 3832 | } |
---|
| 3833 | } |
---|
| 3834 | |
---|
| 3835 | G4int G4Abla::mod(G4int a, G4int b) |
---|
| 3836 | { |
---|
| 3837 | if(b != 0) { |
---|
| 3838 | return (a - (a/b)*b); |
---|
| 3839 | } |
---|
| 3840 | else { |
---|
| 3841 | return 0; |
---|
| 3842 | } |
---|
| 3843 | } |
---|
| 3844 | |
---|
| 3845 | G4double G4Abla::dmod(G4double a, G4double b) |
---|
| 3846 | { |
---|
| 3847 | if(b != 0) { |
---|
| 3848 | return (a - (a/b)*b); |
---|
| 3849 | } |
---|
| 3850 | else { |
---|
| 3851 | return 0.0; |
---|
| 3852 | } |
---|
| 3853 | } |
---|
| 3854 | |
---|
| 3855 | G4double G4Abla::dint(G4double a) |
---|
| 3856 | { |
---|
| 3857 | G4double value = 0.0; |
---|
| 3858 | |
---|
| 3859 | if(a < 0.0) { |
---|
| 3860 | value = double(std::ceil(a)); |
---|
| 3861 | } |
---|
| 3862 | else { |
---|
| 3863 | value = double(std::floor(a)); |
---|
| 3864 | } |
---|
| 3865 | |
---|
| 3866 | return value; |
---|
| 3867 | } |
---|
| 3868 | |
---|
| 3869 | G4int G4Abla::idint(G4double a) |
---|
| 3870 | { |
---|
| 3871 | G4int value = 0; |
---|
| 3872 | |
---|
| 3873 | if(a < 0) { |
---|
| 3874 | value = int(std::ceil(a)); |
---|
| 3875 | } |
---|
| 3876 | else { |
---|
| 3877 | value = int(std::floor(a)); |
---|
| 3878 | } |
---|
| 3879 | |
---|
| 3880 | return value; |
---|
| 3881 | } |
---|
| 3882 | |
---|
| 3883 | G4int G4Abla::idnint(G4double value) |
---|
| 3884 | { |
---|
| 3885 | G4double valueCeil = int(std::ceil(value)); |
---|
| 3886 | G4double valueFloor = int(std::floor(value)); |
---|
| 3887 | |
---|
| 3888 | if(std::fabs(value - valueCeil) < std::fabs(value - valueFloor)) { |
---|
| 3889 | return int(valueCeil); |
---|
| 3890 | } |
---|
| 3891 | else { |
---|
| 3892 | return int(valueFloor); |
---|
| 3893 | } |
---|
| 3894 | } |
---|
| 3895 | |
---|
| 3896 | G4double G4Abla::dmin1(G4double a, G4double b, G4double c) |
---|
| 3897 | { |
---|
| 3898 | if(a < b && a < c) { |
---|
| 3899 | return a; |
---|
| 3900 | } |
---|
| 3901 | if(b < a && b < c) { |
---|
| 3902 | return b; |
---|
| 3903 | } |
---|
| 3904 | if(c < a && c < b) { |
---|
| 3905 | return c; |
---|
| 3906 | } |
---|
| 3907 | return a; |
---|
| 3908 | } |
---|
| 3909 | |
---|
| 3910 | G4double G4Abla::utilabs(G4double a) |
---|
| 3911 | { |
---|
| 3912 | if(a > 0) { |
---|
| 3913 | return a; |
---|
| 3914 | } |
---|
| 3915 | if(a < 0) { |
---|
| 3916 | return (-1*a); |
---|
| 3917 | } |
---|
| 3918 | if(a == 0) { |
---|
| 3919 | return a; |
---|
| 3920 | } |
---|
| 3921 | |
---|
| 3922 | return a; |
---|
| 3923 | } |
---|