Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

eedl.cc @ 1350

Last change on this file since 1350 was 1350, checked in by garnier, 13 years ago
update to last version 4.9.4
File size: 23.4 KB

Line
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	//
27	//
28	// -------------------------------------------------------------------
29	// GEANT 4 class file --- Copyright CERN 1998
30	// CERN Geneva Switzerland
31	//
32	//
33	// File name: eedl
34	// This code transform EEDL data library into
35	// Geant4 format
36	//
37	// Author: V.Ivanchenko
38	//
39	// Creation date: 7 November 2001
40	//
41	// Modifications:
42	//
43	// -------------------------------------------------------------------
44
45	#include "globals.hh"
46	#include "CLHEP/Matrix/Matrix.h"
47	#include "CLHEP/Matrix/SymMatrix.h"
48	#include "G4DataVector.hh"
49	#include "G4VDataSetAlgorithm.hh"
50	#include "G4LinInterpolation.hh"
51	#include "G4LogLogInterpolation.hh"
52	#include "G4LinLogInterpolation.hh"
53	#include "G4SemiLogInterpolation.hh"
54	#include "G4AtomicTransitionManager.hh"
55	#include "G4AtomicShell.hh"
56	#include "G4VEMDataSet.hh"
57	#include "G4EMDataSet.hh"
58	#include "G4ios.hh"
59	#include <fstream>
60	#include <iomanip>
61	#include <string>
62	#include <iostream.h>
63	#include <stdlib.h>
64	#include <strstream.h>
65	//#include "CLHEP/Hist/HBookFile.h"
66	//#include "CLHEP/Hist/Histogram.h"
67
68
69	int main(int argc,char** argv)
70	{
71
72	// -------------------------------------------------------------------
73	// Setup
74	G4int verbose = 0;
75	/*
76	HepHistogram* h1[100];
77	HepHistogram* h2[100];
78	hbookManager = new HBookFile(hFile, 58);
79	h[0] = hbookManager->histogram("Kinetic Energy (MeV)", 50,0.,gEnergy/MeV);
80	*/
81	G4cout.setf( ios::scientific, ios::floatfield );
82
83	// -------------------------------------------------------------------
84	// Control on input
85
86	if(argc < 2) {
87	cout << "Input file is not specified! Exit" << endl;
88	exit(1);
89	}
90
91	ifstream* fin = new ifstream();
92	string fname = argv[1];
93	fin->open(fname.c_str());
94	if( !fin->is_open()) {
95	cout << "Input file <" << fname << "> does not exist! Exit" << endl;
96	exit(1);
97	}
98
99	ofstream* fout_a = new ofstream();
100	string fname1 = "eedl/br-sp.dat";
101	fout_a->open(fname1.c_str(), std::ios::out\|std::ios::trunc);
102
103	ofstream* fout_b = new ofstream();
104	const string fname2 = "eedl/ion-sp-";
105	const string fname4 = "eedl/ion-ex-sig.dat";
106	const string fname5 = "eedl/ion-ex-av.dat";
107
108	ofstream* fout_c = new ofstream();
109	fout_c->open(fname4.c_str(), std::ios::out\|std::ios::trunc);
110	ofstream* fout_d = new ofstream();
111	fout_d->open(fname5.c_str(), std::ios::out\|std::ios::trunc);
112
113	// flags
114
115	G4bool end = false;
116	G4bool first = true;
117	G4bool secon = false;
118	G4bool infor = false;
119	G4bool bremspec= false;
120	G4bool ionis = false;
121	G4bool exsig = false;
122	G4bool exav = false;
123
124	//there can't be lines longer than nmax characters
125	const size_t nmax = 73;
126	char line[nmax];
127	std::string line1, line2;
128
129
130	// parameters
131
132	G4double ein = 0.0;
133	G4double eout = 0.0;
134	// G4double pb,fac,x,s0,sx,sx2,sx3,sx4,sx5,sx6,sy,syx,syx2,syx3,z;
135	// G4double xx,w,sz2,szx,szx2,szx3,syz;
136	G4double x, z, xx, pb;
137	G4double einold = 0.0;
138	G4double eioold = 0.0;
139	G4double eesold = 0.0;
140	G4double eeaold = 0.0;
141	G4DataVector e;
142	G4DataVector p;
143	G4DataVector dp;
144	e.clear();
145	p.clear();
146	dp.clear();
147	G4DataVector xion;
148	G4DataVector pion;
149	G4DataVector tion;
150	xion.clear();
151	pion.clear();
152	tion.clear();
153	G4int nibest = 0;
154	G4int nigood = 0;
155	G4int nibad = 0;
156	G4double zbest = 0.01;
157	G4double zgood = 0.10;
158
159	G4DataVector dv1;
160	G4DataVector dv2;
161	G4DataVector dv3;
162	G4DataVector dv4;
163	dv1.resize(25);
164	dv2.resize(25);
165	dv3.resize(25);
166	dv4.resize(25);
167
168	size_t counter = 0;
169	G4int Zold = -1;
170	G4int shell = 0;
171	G4int Z = 0;
172	G4int finalp, reactionId, dataId;
173	G4double delmax = 0.0;
174	G4double deltav = 0.0;
175	G4double nz = 0.0;
176	G4double ymax = 1.0;
177	G4double zdelmax = 0.0;
178	G4double dsumav = 0.0;
179	G4double dsumax = 0.0;
180	G4double zdeltav = 0.0;
181	G4double x4max = 0.0;
182	G4double znz = 0.0;
183	G4double be = 0.0;
184
185	G4double xdelmax = 0.0;
186	G4double zbad = 0;
187
188	// main loop
189
190	do {
191
192	// read next line
193
194	counter++;
195	for(size_t ii = 0; ii < nmax; ii++) {line[ii] = ' ';}
196	fin->getline( line, nmax);
197	line1 = std::string("");
198	line1 = std::string(line, nmax);
199	if(1 < verbose) {
200	cout << "Next line # " << counter << ": " << line1 << endl;
201	}
202
203	// analize line contence
204
205	// end of file
206	if(fin->eof()) {
207	end = true;
208	(*fout_a) << "-2 -2" << endl;
209	(*fout_b) << "-2 -2" << endl;
210	(*fout_c) << "-2 -2" << endl;
211	(*fout_d) << "-2 -2" << endl;
212
213	// end of data set
214	} else if(line[71] == '1') {
215
216	first = true;
217	infor = false;
218	ein = 0.0;
219
220	// first line in a data set
221	} else if(first) {
222
223	line2 = std::string(&line[0], 3);
224	Z = atoi(line2.c_str());
225	line2 = std::string(&line[10], 2);
226	finalp = atoi(line2.c_str());
227	if(0 < verbose) {
228	cout << "New data for Z= " << Z
229	<< " finalParticleId= " << finalp
230	<< endl;
231	}
232	first = false;
233	secon = true;
234
235	// second line in a data set
236	} else if(secon) {
237
238	line2 = std::string(&line[0], 2);
239	reactionId = atoi(line2.c_str());
240	line2 = std::string(&line[2], 3);
241	dataId = atoi(line2.c_str());
242	if(-1 < verbose) {
243	cout << " reactionId= " << reactionId
244	<< " dataId= " << dataId
245	<< endl;
246	}
247
248	// Bremsstrahlung data
249
250	if(reactionId == 82 && (dataId == 0 \|\| dataId == 21)) {
251	verbose = 1;
252	if(dataId == 21) bremspec= true;
253	einold = 0.0;
254	ein = 0.0;
255
256	// Delta-electron data
257
258	} else if(reactionId == 81 && (dataId == 0 \|\| dataId == 21)) {
259	verbose = 1;
260	eioold = 0.0;
261	ein = 0.0;
262	if(dataId == 21) {
263	ionis= true;
264	if(Z > Zold) {
265	if(Z > 1) {
266	(*fout_b) << "-1 -1" << endl;
267	(*fout_b) << "-2 -2" << endl;
268	fout_b->close();
269	}
270	shell = 0;
271	Zold = Z;
272	char nameChar[20] = {""};
273	std::ostrstream ost(nameChar, 20, std::ios::out);
274	ost << fname2 << Z << ".dat";
275	string fname3(nameChar);
276	fout_b->open(fname3.c_str(), std::ios::out\|std::ios::trunc);
277	} else {
278	shell++;
279	(*fout_b) << "-1 -1" << endl;
280	}
281	if(Z <= 99) {
282	be =((G4AtomicTransitionManager::Instance())->
283	Shell(Z, shell)->BindingEnergy());
284	if(0 < verbose) {
285	G4cout << "New ioni data Z= " << Z
286	<< " shell= " << shell
287	<< " be= " << be << G4endl;
288	}
289	} else {
290	ionis = false;
291	}
292	}
293	// Excitation
294
295	} else if(reactionId == 83 && (dataId == 0 \|\| dataId == 11)) {
296	verbose = 1;
297	eesold = 0.0;
298	eeaold = 0.0;
299	ein = 0.0;
300	if(dataId == 0) {
301	exsig = true;
302	(*fout_c) << "0 0" << endl;
303	}
304	if(dataId == 11) {
305	exav = true;
306	(*fout_d) << "0 0" << endl;
307	}
308
309	} else {
310	verbose = 0;
311	}
312
313	secon = false;
314	infor = true;
315
316	// line with information in a data set
317
318	} else if(infor) {
319
320	if(bremspec \|\| ionis) {
321	line2 = std::string(&line[0], 15);
322	ein = atof(line2.c_str());
323	line2 = std::string(&line[15], 13);
324	eout = atof(line2.c_str());
325	line2 = std::string(&line[28], 13);
326	pb = atof(line2.c_str());
327	if(1 < verbose) {
328	G4cout << ein << " " << eout << " " << pb << G4endl;
329	}
330	} else if(exsig \|\| exav) {
331	line2 = std::string(&line[0], 15);
332	ein = atof(line2.c_str());
333	line2 = std::string(&line[15], 13);
334	pb = atof(line2.c_str());
335	if(1 < verbose) {
336	G4cout << ein << " " << pb << G4endl;
337	}
338	}
339	}
340
341	// handle with Bremsstrahlung data
342
343	if(bremspec) {
344
345	// end of sample for given energy
346	if(einold > 1.0*eV && (first \|\| ein != einold)) {
347
348	size_t nn = e.size();
349	size_t i;
350	G4double x0 = 0.01;
351	size_t imax = 15;
352	G4double dx = 0.1;
353
354
355	G4DataVector* eee = new G4DataVector();
356	G4DataVector* ppp = new G4DataVector();
357	for (i=0; i<nn; i++) {
358	eee->push_back(e[i]);
359	ppp->push_back(p[i]);
360	}
361
362	G4VDataSetAlgorithm* inpl = new G4LinInterpolation();
363	G4VEMDataSet* dataSet = new G4EMDataSet(1, eee, ppp, inpl, 1., 1.);
364
365	G4DataVector* eeee = new G4DataVector();
366	G4DataVector* pppp = new G4DataVector();
367	G4VDataSetAlgorithm* inpll = new G4LinInterpolation();
368	G4double a = dataSet->FindValue(x0);
369	for(i=0; i<imax; i++) {
370	G4double en = dx*((G4double)i);
371	if(i == 0) en = 0.01;
372	if(i == 10) en = 0.95;
373	if(i == 11) en = 0.97;
374	if(i == 12) en = 0.99;
375	if(i == 13) en = 0.995;
376	if(i == 14) en = 1.00;
377	G4double pn = (dataSet->FindValue(en))/a;
378	eeee->push_back(en);
379	pppp->push_back(pn);
380	}
381	G4VEMDataSet* ndataSet = new G4EMDataSet(1, eeee, pppp, inpll, 1., 1.);
382
383
384	G4double delm = 0.0;
385	G4double dy, yy;
386	ymax = 1.0;
387	for (i=0; i<nn; i++) {
388	x = e[i];
389	if(x <= x0) {
390	yy = 1.0;
391	} else {
392	yy = ndataSet->FindValue(x);
393	}
394
395	if(yy > ymax) ymax = yy;
396
397	dy = std::abs(p[i] - yy*a)/p[i];
398
399	if(1 < verbose) {
400	cout << "x= " << x
401	<< " p[i]= " << p[i]
402	<< " dp= " << dy
403	<< endl;
404	}
405	if(dy > delm) delm = dy;
406	}
407
408
409	if(einold < 100.*MeV) {
410	if(delm > delmax) delmax = delm;
411	nz += 1.0;
412	deltav = (deltav*(nz - 1.0) + delm)/nz;
413	}
414
415	(*fout_a).precision(4);
416
417	(*fout_a) << einold/MeV << " "
418	<< ymax
419	<< endl;
420	for(i=0; i<imax; i++) {
421	(fout_a) << (pppp)[i] << " ";
422	}
423
424	(*fout_a) << endl;
425	if(0 < verbose) {
426	cout << "!!!!! e= " << einold/MeV << " x0= " << x0
427	<< " dymax= " << delm
428	<< " dm= " << delmax
429	<< " da= " << deltav
430	<< " ymax= " << ymax
431	<< endl;
432	}
433	einold = 0.0;
434	delete dataSet;
435	delete ndataSet;
436	}
437
438
439	// new sample
440	if(ein != einold) {
441	e.clear();
442	p.clear();
443	dp.clear();
444	einold = ein;
445	}
446
447	if(einold > 0.1*eV) {
448	x = eout/einold;
449	e.push_back(x);
450	p.push_back(pb*x);
451	dp.push_back(pb);
452
453	if(1 < verbose) {
454	cout << "x= " << x
455	<< " p= " << pb << " pp= " << x*pb << endl;
456	}
457	}
458	if(first) {
459	(*fout_a) << "-1 " << Z << endl;
460	bremspec = false;
461	}
462	}
463
464	if(ionis) {
465
466	// end of sample for given energy
467	size_t nn = xion.size();
468	if(first \|\| ein != eioold && (Z == 82 \|\| Z == 27 && Z == 1)) {
469	cout << "??? ein= " << ein << "; nn= " << nn << endl;
470	}
471	if(eioold > 1.0*eV && nn > 4 && (first \|\| ein != eioold)) {
472
473	G4double gam = (eioold + be)/electron_mass_c2 + 1.;
474	G4double g = (2.gam - 1.)/(gamgam);
475
476	G4double x1 = (0.1*eV + be)/(eioold + be);
477	if(x1 > 0.5) x1 = 0.5;
478	G4double x2 = 0.5;
479
480	x = xion[nn-1];
481	z = 1. - gx + (1. - g)xx + xx*(1./(1. - x) - g)/(1. - x);
482
483	G4double f1 = tion[nn-1]/z;
484	G4double f2 = 0.0;
485	G4double f3 = 0.0;
486
487	// aria interpolation
488	verbose = 1;
489	size_t jmax1 = 3;
490	size_t jmax2 = 16;
491	size_t jmax = jmax1 + jmax2 + 1;
492	G4int zpr = -54;
493	G4double range = (G4double)jmax;
494	G4double xmax = xion[0];
495	G4double tmax = tion[0];
496	G4double y = 0.0;
497
498	for(size_t i=0; i<nn; i++) {
499	tion[i] /= f1;
500	if(tion[i] > tmax && i < nn-1) {
501	tmax = tion[i];
502	xmax = xion[i];
503	}
504	}
505
506	G4double x3 = xmax * 10.0;
507	if(x3 > 0.5) x3 = 0.5;
508	G4double dx = (x3-x1)/range;
509	G4double x4 = x3 + dx;
510
511	for(i=0; i<nn; i++) {
512	x = xion[i];
513	if(x < 0.499 && x >= x3 ) {
514	z = 1. - gx + (1. - g)xx + xx*(1./(1. - x) - g)/(1. - x);
515	f2 += (tion[i] - z)/(0.5 - x);
516	f3 += 1.0/x;
517	}
518	}
519	if(f3 > 0.0) f2 /= f3;
520
521	G4double zdelm = 0.0;
522	G4double zdelm0= 0.0;
523	G4double xdelmax0 = 0.0;
524	G4double xdelmax1 = x3;
525	G4bool im = true;
526	G4bool out = false;
527	size_t iiiimax = 2;
528
529	for(size_t iiii=0; iiii<iiiimax; iiii++) {
530
531	G4double range1 = (G4double)jmax1;
532	G4double range2 = (G4double)jmax2;
533	jmax = jmax1 + jmax2 + 1;
534	// lin
535	G4double dx1 = (xmax - x1)/range1;
536
537	// log
538	G4double dx2 = std::log(x3/xmax)/range2;
539	x4 = x3*std::exp(dx2);
540
541	G4DataVector* eee = new G4DataVector();
542	G4DataVector* ppp = new G4DataVector();
543	G4VDataSetAlgorithm* inpl = new G4LinLogInterpolation();
544
545	eee->push_back(x1);
546	ppp->push_back(tion[0]);
547
548	for(i=0; i<nn; i++) {
549	if(xion[i] < x4) {
550	eee->push_back(xion[i]);
551	ppp->push_back(tion[i]);
552	}
553	}
554
555	G4VEMDataSet* dataSet = new G4EMDataSet(1,eee,ppp,inpl,1.,1.);
556
557	G4DataVector* eeee = new G4DataVector();
558	G4DataVector* pppp = new G4DataVector();
559	G4VDataSetAlgorithm* inpll = new G4LinInterpolation();
560
561	for(i=0; i<jmax; i++) {
562
563	if(i <= jmax1) {
564	x = x1 + dx1 * (G4double)i;
565
566	} else {
567
568	x = xmaxstd::exp(dx2 (G4double)(i - jmax1));
569	}
570
571	eeee->push_back(x);
572	if(i < jmax-1) {
573	y = dataSet->FindValue(x);
574	} else {
575	z = 1. - gx + (1. - g)xx + xx*(1./(1. - x) - g)/(1. - x);
576	y = z + f2*(0.5 - x)/x;
577	}
578	pppp->push_back(y);
579	}
580
581	G4VEMDataSet* ndataSet = new G4EMDataSet(1,eeee,pppp,inpll,1.,1.);
582
583
584	G4double dy = 0.0;
585	zdelm = 0.0;
586
587	x = x2;
588	z = 1. - gx + (1. - g)xx + xx*(1./(1. - x) - g)/(1. - x);
589	G4double zmax2 = z;
590	G4double xmax2 = x2;
591	G4double zmax = 0.0;
592	G4bool poin = true;
593	G4double yy, yyy;
594	G4double dsum = 0.0;
595	G4double sum = 0.0;
596	G4double x0 = 0.0;
597	G4double y0 = 0.0;
598	G4double yy0= 0.0;
599
600	for (i=0; i<nn; i++) {
601	x = xion[nn - i - 1];
602	y = tion[nn - i - 1];
603	z = 1. - gx + (1. - g)xx + xx*(1./(1. - x) - g)/(1. - x);
604	yyy = z + f2*(0.5 - x)/x;
605
606	if(x <= x3) {
607	yy = ndataSet->FindValue(x);
608
609	} else {
610	yy = yyy;
611	}
612
613
614	if(poin && yy > zmax2) {
615	poin = false;
616	xmax2 = xion[i+1];
617	} else if(!poin && yy*x > zmax && x >= xmax) {
618	zmax = yy*x;
619	}
620
621	dy = 1.0 - yy/y;
622
623	if((1 < verbose \|\| (Z == zpr)) && iiii == 1 ) {
624	cout << "x= " << x
625	<< " tion[i]= " << y
626	<< " delta= " << dy
627	<< " x1= " << x1
628	<< " x2= " << xmax
629	<< " x3= " << x3
630	<< endl;
631	}
632	dy = std::abs(dy);
633	if(dy > zdelm) {
634	zdelm = dy;
635	xdelmax0 = x;
636	}
637	if(im && dy > zgood) {
638	im = false;
639	xdelmax1 = x;
640	}
641	if(i > 0) {
642	G4double xsq = x*x;
643	G4double xsq0= x0*x0;
644	sum += (y/xsq + y0/xsq0)(x0 - x)0.5;
645	dsum += std::abs((y/xsq + y0/xsq0) - (yy/xsq + yy0/xsq0))(x0 - x)0.5;
646	}
647	yy0 = yy;
648	y0 = y;
649	x0 = x;
650	}
651
652
653	if(iiii == iiiimax-1 \|\| zdelm <= zgood ) {
654
655	out = true;
656	//if(Z == zpr) ndataSet->PrintData();
657	G4String qu;
658	if(dsum/sum <= zbest) {
659	nibest++;
660	qu = " ##BEST";
661	} else if(dsum/sum <= zgood) {
662	nigood++;
663	qu = " ##GOOD";
664	} else {
665	nibad++;
666	qu = " ##BAD";
667	}
668	if(zdelm > zdelmax) {
669	zdelmax = zdelm;
670	xdelmax = xdelmax0;
671	zbad = Z;
672	}
673
674	znz += 1.0;
675	zdeltav = (zdeltav*(znz - 1.0) + zdelm)/znz;
676
677	dsum /= sum;
678	if(dsum > dsumax) dsumax = dsum;
679	dsumav = (dsumav*(znz - 1.0) + dsum)/znz;
680
681	yy0 = be + eioold;
682	f2 *= yy0;
683	x1 *= yy0;
684	xmax *= yy0;
685	x3 *= yy0;
686
687	dv1 = dv2;
688	dv2 = dv3;
689	dv3[0] = eioold/MeV;
690	dv3[1] = f2;
691	dv3[2] = x1;
692	dv3[3] = xmax;
693	dv3[4] = x3;
694	for (i=5; i<25; i++) {
695	dv3[i] = (*pppp)[i-5];
696	}
697
698	for(G4double edu=1.0; edu<11.; edu *=3.0) {
699
700	if(edu > dv2[0] && edu < dv3[0]){
701
702	HepMatrix* m = new HepMatrix(3,3);
703	G4double u = std::log10(edu);
704	G4double u1 = std::log10(dv1[0]);
705	G4double u2 = std::log10(dv2[0]);
706	G4double u3 = std::log10(dv3[0]);
707	dv4[0] = edu;
708
709	(*m)(1,1) = 1.0;
710	(*m)(2,1) = u1;
711	(m)(3,1) = u1u1;
712	(*m)(1,2) = 1.0;
713	(*m)(2,2) = u2;
714	(m)(3,2) = u2u2;
715	(*m)(1,3) = 1.0;
716	(*m)(2,3) = u3;
717	(m)(3,3) = u3u3;
718	G4int j = 0;
719	HepMatrix invm = m->inverse(j);
720	G4double a1, a2, a3, q1, q2, q3;
721
722	for(i=1; i<25; i++) {
723	q1 = dv1[i];
724	q2 = dv2[i];
725	q3 = dv3[i];
726	a1 = q1invm(1,1) + q2invm(2,1) + q3*invm(3,1);
727	a2 = q1invm(1,2) + q2invm(2,2) + q3*invm(3,2);
728	a3 = q1invm(1,3) + q2invm(2,3) + q3*invm(3,3);
729	dv4[i] = a1 + a2u + a3u*u;
730	}
731	(*fout_b).precision(5);
732	(*fout_b) << dv4[0] << " " << dv4[1] << endl;
733	for (i=2; i<25; i++) {
734	if(i == 5) (*fout_b).precision(4);
735	(*fout_b) << dv4[i] << " ";
736	}
737	(*fout_b) << endl;
738	delete m;
739	}
740	}
741
742	(*fout_b).precision(5);
743	(*fout_b) << dv3[0] << " " << dv3[1] << endl;
744	for (i=2; i<25; i++) {
745	if(i == 5) (*fout_b).precision(4);
746	(*fout_b) << dv3[i] << " ";
747	}
748	(*fout_b) << endl;
749
750	if(0 < verbose \|\| (Z == zpr)) {
751	cout << "!! e= " << eioold/MeV
752	<< " dMax= " << zdelm
753	<< " dMaxAll= " << zdelmax
754	<< " xMaxAll= " << xdelmax
755	<< " dAv= " << zdeltav
756	<< " A= " << f2
757	<< " x1= " << x1
758	<< " x2= " << xmax
759	<< " x3= " << x3
760	<< endl;
761	cout << " dsum= " << dsum
762	<< " dSumMax= " << dsumax
763	<< " dSumAv= " << dsumav
764	<< " zmax= " << zmax
765	<< " zmax2= " << zmax2
766	<< " xmax2= " << xmax2
767	<< " jmax= " << jmax
768	<< " " << qu
769	<< endl;
770
771	for (i=0; i<pppp->size(); i++) {
772	cout << (*pppp)[i] << " ";
773	}
774	cout << " Zbad= " << zbad << endl;
775
776	}
777	} else if(iiii == 0) {
778	zdelm0 = zdelm;
779	/*
780	if(xdelmax1 > x3) {
781	//jmax3 = 12;
782	x3 = xdelmax1;
783	}
784	*/
785	}
786	delete dataSet;
787	delete ndataSet;
788	if(out) break;
789
790	}
791	einold = 0.0;
792	}
793
794
795	// new sample
796	if(ein != eioold) {
797	xion.clear();
798	pion.clear();
799	tion.clear();
800	eioold = ein;
801	}
802
803	if(eioold > 0.1*eV) {
804	x = (eout + be)/(eioold + be);
805	xion.push_back(x);
806	pion.push_back(pb);
807	xx = pbxx;
808	tion.push_back(xx);
809	if(1 < verbose) {
810	cout << "x= " << x
811	<< " p= " << pb
812	<< " px2= " << xx
813	<< " b= " << be
814	<< " ein= " << ein << endl;
815	}
816	}
817	if(first) {
818	ionis = false;
819	}
820	}
821	if(exsig) {
822
823	(*fout_c).precision(4);
824	if(ein != eesold) eesold = ein;
825	if(eesold >= 1.0eV && eesold <= 100.MeV) {
826	(*fout_c) << eesold/MeV << " "
827	<< pb*1.e-6
828	<< endl;
829	if(1 < verbose) {
830	cout << "!!!!! e= " << eesold/MeV << " sig= " << pb
831	<< endl;
832	}
833	}
834	if(first) {
835	(*fout_c) << "-1 " << Z <<endl;
836	exsig = false;
837	}
838	}
839	if(exav) {
840	(*fout_d).precision(4);
841	if(ein != eeaold) eeaold = ein;
842	if(eeaold > 1.0eV && eeaold <= 100.MeV) {
843	(*fout_d) << eeaold/MeV << " "
844	<< pb*1.e+6
845	<< endl;
846	if(1 < verbose) {
847	cout << "!!!!! e= " << eeaold/MeV << " Eav= " << pb/MeV
848	<< endl;
849	}
850	}
851	if(first) {
852	(*fout_d) << "-1 " << Z << endl;
853	exav = false;
854	}
855	}
856
857	} while(!end && counter < 1000000);
858
859	G4cout << "nibest = " << nibest << G4endl;
860	G4cout << "nigood = " << nigood << G4endl;
861	G4cout << "nibad = " << nibad << G4endl;
862	G4cout << "dMaxAll= " << zdelmax << G4endl;
863	G4cout << "xMaxAll= " << xdelmax << G4endl;
864	G4cout << "dAv = " << zdeltav << G4endl;
865	G4cout << "dSumMax= " << dsumax << G4endl;
866	G4cout << "dSumAv = " << dsumav << G4endl;
867	G4cout << "x4max = " << x4max << G4endl;
868	G4cout << "Zbad = " << zbad << G4endl;
869
870	G4cout << "###### End of test #####" << G4endl;
871	}
872

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: trunk/source/processes/electromagnetic/lowenergy/test/eedl.cc @ 1350

Download in other formats: