-
Notifications
You must be signed in to change notification settings - Fork 3
/
Copy pathgixos.cpp
122 lines (113 loc) · 4.16 KB
/
gixos.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
//scattering in pass-through geometry
//pass through geometry, beam pass through Si crystal
#include<iostream>
#include<fstream>
#include<sstream>
#include <vector>
#include <grace_np.h>
#include <sys/time.h>
#include <sys/resource.h>
#include "compton.h"
#include "block.h"
#include "randomInit.h" //init random number generator
using namespace std;
int main(int argc, char *argv[])
{
char randomBuffer[257];
randomInit(randomBuffer); //init random number generator
double pixelHeight=2./900.// angular now
,s5h=20/900.;// slit size
vector<double> dist_theta;
int theta_steps=(int)(1/pixelHeight+0.5);
double itheta_steps=theta_steps/(1.);
dist_theta.resize(theta_steps+200);
for (unsigned int i=0;i<dist_theta.size();i++) dist_theta.at(i)=0.;
int j;
/*
unsigned int xaxisDiv=8;
for (unsigned int ii=0;ii<=xaxisDiv;ii++){
if ( ii & (unsigned int)1) {
}else {
if (!ii){
continue;
}
int ii1=mcd(ii,xaxisDiv);
//cout<<ii<<' '<<xaxisDiv<<' '<<ii1<<endl;
if (xaxisDiv==ii){
continue;
}
}
}
*/
double diameter=0.;
if (argc>=2) { // read sample block diameter from command line
istringstream is0(argv[1]);
is0>>diameter;
if ( !( diameter>=0.01 && diameter < 1000.)) diameter=0.25;
} else diameter=0.25; // all units in inches
ostringstream os0;
os0<<"si-gixos-"<<diameter<<".txt";
string fn(os0.str().c_str());
cout<<"Deleting output file "<<fn<<endl;
unlink(fn.c_str());
int l1=100,l2=50000000;
ofstream out1;
double l3=0.;
int ii=1;
double muC= 0.15 * 2.33; // in cm^-1, compton
double muTotal=1.31 * 2.33; // in cm^-1, total
double compton_ratio_factor=muC/muTotal; //weight factor for each scattering
double psum;
double pEn=30.; //KeV
// double pLambda=E_to_l(pEn);
// double pk0=2*M_PI/pLambda;
//double fac1=1./(2.*pk0*l1); // q_z range up to 1 angstrom^-1
thetaDistribution tP0(pEn);
photon p0(diameter);
struct rusage r_start,r_end;
getrusage(RUSAGE_SELF, &r_start);
psum=0;
while (ii<l1){
//double t1= -ii*fac1;
//#pragma omp parallel for private( j,t0 ) reduction(+:sum) schedule(static)
for (int i=0;i<l2;i++){
//if (p0.initRefBuried(t1)) continue;
p0.initGixos();
while ((j=p0.propagateGixos(tP0.theta()))==1);
/*{
if(p0.scattered>1000){
cout<<p0.r.x<<' '<<p0.r.y<<' '<<p0.r.z<<endl;
cout<<p0.o.get_theta()<<' '<<p0.o.get_phi()<<endl;
break;
}
};
*/
// if( i && ((i>>18 ) <<18) == i) cout<<"i="<<i<<endl;
if (j==-1 ) {//scattered out of sample
// cout<<i<<' '<<p0.o.theta<<endl;
// cout<<i<<' '<<tP0.theta()<<' '<<p0.o.st<<' '<<p0.o.sp<<endl;
if (p0.scattered && p0.o.sp>=0. && fabs(p0.o.st)<s5h) dist_theta.at((unsigned int) (p0.o.sp*itheta_steps)) += pow(compton_ratio_factor,(int) p0.scattered);
}
//if( (i>>14)<<14 == i) cout<<i<<endl;
}
getrusage(RUSAGE_SELF, &r_end); //get running time
l3+=l2;
cout<<ii<<' '<<dist_theta.at(1)<<' '<<l3<<' '<<l2/(r_end.ru_utime.tv_sec -r_start.ru_utime.tv_sec+ double(1e-6)*(r_end.ru_utime.tv_usec -r_start.ru_utime.tv_usec ) )<<" P/s\n";
r_start=r_end;
/*
*/
out1.open(fn.c_str());
double fac2=1./(s5h*pixelHeight*l3);
double dx=1./itheta_steps;
double x0=0.5*dx;
for (unsigned int i2=0;x0<1.;i2++){
out1<<x0*(2.*M_PI/0.4132)<<' '<<dist_theta.at(i2)*fac2*(1-x0*x0)<<' '<<dist_theta.at(i2)<<' '<<l3<<endl;
x0+=dx;
}
//out1<<dist_theta.size()/itheta_steps<<' '<<dist_theta.at((dist_theta.size()-1)>>1)*fac2<<' '<<s5h<<' '<<pixelHeight<<endl;
//cout<<1./itheta_steps<<' '<<dist_theta.at((dist_theta.size()-1)>>1)*fac2<<"s5h(angular)= "<<s5h<<"pixelSize(angular)="<<pixelHeight<<endl;
out1.close();
ii++;
}
return(0);
}