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