375 lines
9.8 KiB
Common Lisp
375 lines
9.8 KiB
Common Lisp
#pragma OPENCL EXTENSION cl_khr_fp64 : enable
|
|
|
|
/******Random numbers********/
|
|
|
|
/* struct for random number state */
|
|
typedef struct {
|
|
|
|
double s10;
|
|
double s11;
|
|
double s12;
|
|
double s20;
|
|
double s21;
|
|
double s22;
|
|
double z;
|
|
bool gen;
|
|
|
|
} RNDState;
|
|
|
|
#define NORM 2.328306549295728e-10
|
|
#define M1 4294967087.0
|
|
#define M2 4294944443.0
|
|
#define A12 1403580.0
|
|
#define A13N 810728.0
|
|
#define A21 527612.0
|
|
#define A23N 1370589.0
|
|
|
|
/* MRG32k3a uniform random number generator */
|
|
double rand_uniform(RNDState *s) {
|
|
long k;
|
|
double p1, p2;
|
|
|
|
/* Component 1 */
|
|
p1 = A12 * (*s).s11 - A13N * (*s).s10;
|
|
k = p1 / M1;
|
|
p1 -= k * M1;
|
|
if (p1 < 0.0)
|
|
p1 += M1;
|
|
(*s).s10 = (*s).s11;
|
|
(*s).s11 = (*s).s12;
|
|
(*s).s12 = p1;
|
|
|
|
/* Component 2 */
|
|
p2 = A21 * (*s).s22 - A23N * (*s).s20;
|
|
k = p2 / M2;
|
|
p2 -= k * M2;
|
|
if (p2 < 0.0)
|
|
p2 += M2;
|
|
(*s).s20 = (*s).s21;
|
|
(*s).s21 = (*s).s22;
|
|
(*s).s22 = p2;
|
|
|
|
/* Combination */
|
|
if (p1 <= p2)
|
|
return ((p1 - p2 + M1) * NORM);
|
|
else
|
|
return ((p1 - p2) * NORM);
|
|
}
|
|
|
|
/* get random variable with gaussian distribution */
|
|
double rand_normal(RNDState *s, double mu, double sigma) {
|
|
|
|
const double two_pi = 2.0 * 3.141592653589793223846;
|
|
double z0;
|
|
|
|
if (!(*s).gen) {
|
|
(*s).gen = true;
|
|
return (*s).z * sigma + mu;
|
|
}
|
|
|
|
double u1, u2;
|
|
u1 = rand_uniform(s);
|
|
u2 = rand_uniform(s);
|
|
|
|
z0 = sqrt(-2.0 * log(u1)) * cos(two_pi * u2);
|
|
(*s).z = sqrt(-2.0 * log(u1)) * sin(two_pi * u2);
|
|
(*s).gen = false;
|
|
|
|
return z0 * sigma + mu;
|
|
|
|
|
|
}
|
|
|
|
/* initialize random states */
|
|
__kernel void initRand(__global RNDState *s, unsigned int seed, int N) {
|
|
|
|
int id = get_global_id(0);
|
|
|
|
if (id < N) {
|
|
RNDState tmp;
|
|
int tmp_seed = 2*id;// * 0x100000000ULL;
|
|
tmp.s10 = 12345 + tmp_seed;
|
|
tmp.s11 = 12345 + tmp_seed;
|
|
tmp.s12 = 12345 + tmp_seed;
|
|
tmp.s20 = 12345 + tmp_seed;
|
|
tmp.s21 = 12345 + tmp_seed;
|
|
tmp.s22 = 12345 + tmp_seed;
|
|
|
|
|
|
tmp.z = 0;
|
|
tmp.gen = true;
|
|
|
|
s[id] = tmp;
|
|
}
|
|
|
|
}
|
|
|
|
/* create random numbers and fill an array */
|
|
__kernel void createRandoms(__global RNDState *states, __global double *data, int size) {
|
|
|
|
int idx = get_global_id(0);
|
|
|
|
if (idx < size) {
|
|
RNDState s = states[idx];
|
|
data[idx] = rand_uniform(&s);
|
|
states[idx] = s;
|
|
}
|
|
|
|
}
|
|
|
|
|
|
/**********Degrader**********/
|
|
enum PARAMS { POSITION,
|
|
ZSIZE,
|
|
M_P,
|
|
C,
|
|
RHO_M,
|
|
PI,
|
|
AVO,
|
|
R_E,
|
|
eM_E,
|
|
Z_M,
|
|
A_M,
|
|
A2_C,
|
|
A3_C,
|
|
A4_C,
|
|
A5_C,
|
|
Z_P,
|
|
X0_M,
|
|
I_M,
|
|
DT_M};
|
|
|
|
|
|
typedef struct {
|
|
int label;
|
|
unsigned localID;
|
|
double3 Rincol;
|
|
double3 Pincol;
|
|
} PART;
|
|
|
|
double Dot(double3 d1, double3 d2) {
|
|
return (d1.x * d2.x + d1.y * d2.y + d1.z * d2.z);
|
|
}
|
|
|
|
/* check if particle is in degrader material */
|
|
bool checkHit(double z, double position, double zsize) {
|
|
return ( ( z > position) && ( z <= position + zsize) );
|
|
}
|
|
|
|
/* calculate particles energy loss */
|
|
void energyLoss(double *Eng, bool *pdead, double deltat, RNDState *s, __local double *par) {
|
|
|
|
double dEdx = 0.0;
|
|
double gamma = ( (*Eng) + par[M_P]) / par[M_P];
|
|
|
|
double gamma2 = gamma * gamma;
|
|
|
|
double beta = sqrt(1.0 - 1.0 / gamma2);
|
|
double beta2 = beta * beta;
|
|
double deltas = deltat * beta * par[C];
|
|
double deltasrho = deltas * 100 * par[RHO_M];
|
|
double K = 4.0 * par[PI] * par[AVO] * par[R_E] * par[R_E] * par[eM_E] * 1E7;
|
|
double sigma_E = sqrt(K * par[eM_E] * par[RHO_M] * (par[Z_M]/par[A_M])* deltas * 1E5);
|
|
|
|
if (((*Eng) > 0.00001) && ((*Eng) < 0.0006)) {
|
|
double Ts = ((*Eng)*1E6)/1.0073;
|
|
double epsilon_low = par[A2_C]*pow(Ts,0.45);
|
|
double epsilon_high = (par[A3_C]/Ts)*log(1+(par[A4_C]/Ts)+(par[A5_C]*Ts));
|
|
double epsilon = (epsilon_low*epsilon_high)/(epsilon_low + epsilon_high);
|
|
dEdx = - epsilon /(1E21*(par[A_M]/par[AVO]));
|
|
double delta_Eave = deltasrho * dEdx;
|
|
double delta_E = delta_Eave + rand_normal(s, 0, sigma_E);
|
|
|
|
(*Eng) = (*Eng) + delta_E / 1E3;
|
|
}
|
|
|
|
if ((*Eng) >= 0.0006) {
|
|
double Tmax = 2.0 * par[eM_E] * 1e9 * beta2 * gamma2 /
|
|
(1.0 + 2.0 * gamma * par[eM_E] / par[M_P] +
|
|
(par[eM_E] / par[M_P]) * (par[eM_E] / par[M_P]));
|
|
dEdx = -K * par[Z_P] * par[Z_P] * par[Z_M] / (par[A_M] * beta2) *
|
|
(1.0 / 2.0 * log(2 * par[eM_E] * 1e9 * beta2 * gamma2 *
|
|
Tmax / par[I_M] / par[I_M]) - beta2);
|
|
|
|
double delta_Eave = deltasrho * dEdx;
|
|
double delta_E = delta_Eave + rand_normal(s, 0, sigma_E);
|
|
|
|
(*Eng) = (*Eng)+delta_E / 1E3;
|
|
}
|
|
|
|
(*pdead) = (((*Eng)<1E-4) || (dEdx>0));
|
|
|
|
}
|
|
|
|
/* rotate partocle */
|
|
void Rot(double3 *P, double3 *R, double xplane,
|
|
double normP, double thetacou, double deltas, int coord,
|
|
__local double *par)
|
|
{
|
|
double Psixz;
|
|
double pxz;
|
|
|
|
double px = (*P).x;
|
|
double pz = (*P).z;
|
|
double x = (*R).x;
|
|
double z = (*R).z;
|
|
|
|
if (px>=0 && pz>=0) Psixz = atan(px/pz);
|
|
else if (px>0 && pz<0)
|
|
Psixz = atan(px/pz) + par[PI];
|
|
else if (px<0 && pz>0)
|
|
Psixz = atan(px/pz) + 2*par[PI];
|
|
else
|
|
Psixz = atan(px/pz) + par[PI];
|
|
|
|
pxz = sqrt(px*px + pz*pz);
|
|
if(coord==1) {
|
|
(*R).x = x + deltas * px/normP + xplane*cos(Psixz);
|
|
(*R).z = z - xplane * sin(Psixz);
|
|
}
|
|
if(coord==2) {
|
|
(*R).x = x + deltas * px/normP + xplane*cos(Psixz);
|
|
(*R).z = z - xplane * sin(Psixz) + deltas * pz / normP;
|
|
}
|
|
(*P).x = pxz*cos(Psixz)*sin(thetacou) + pxz*sin(Psixz)*cos(thetacou);
|
|
(*P).z = -pxz*sin(Psixz)*sin(thetacou) + pxz*cos(Psixz)*cos(thetacou);
|
|
}
|
|
|
|
|
|
void coulombScat(double3 *R, double3 *P, double deltat,
|
|
RNDState *s, __local double* par) {
|
|
|
|
double dotP = Dot((*P), (*P));
|
|
|
|
double Eng = sqrt(dotP + 1.0) * par[M_P] - par[M_P];
|
|
double gamma = (Eng + par[M_P]) / par[M_P];
|
|
double beta = sqrt(1.0 - 1.0 / (gamma * gamma));
|
|
double normP = sqrt(dotP);
|
|
double deltas = deltat * beta * par[C];
|
|
double theta0 = 13.6e6 / (beta * sqrt(dotP) * par[M_P] * 1e9) *
|
|
par[Z_P] * sqrt(deltas / par[X0_M]) * (1.0 + 0.038 * log(deltas / par[X0_M]));
|
|
|
|
// x-direction: See Physical Review, "Multiple Scattering"
|
|
double z1 = rand_normal(s, 0, 1);//curand_normal_double(&state);//gsl_ran_gaussian(rGen_m,1.0);
|
|
double z2 = rand_normal(s, 0, 1);//curand_normal_double(&state);//gsl_ran_gaussian(rGen_m,1.0);
|
|
double thetacou = z2 * theta0;
|
|
|
|
while(fabs(thetacou) > 3.5 * theta0) {
|
|
z1 = rand_normal(s, 0, 1);//curand_normal_double(&state);//gsl_ran_gaussian(rGen_m,1.0);
|
|
z2 = rand_normal(s, 0, 1);//curand_normal_double(&state);//gsl_ran_gaussian(rGen_m,1.0);
|
|
thetacou = z2 * theta0;
|
|
}
|
|
|
|
double xplane = z1 * deltas * theta0 / sqrt(12.0) + z2 * deltas * theta0 / 2.0;
|
|
int coord = 1;
|
|
Rot(P, R, xplane, normP, thetacou, deltas, coord, par);
|
|
|
|
double P2 = rand_uniform(s);//curand_uniform_double(&state);//gsl_rng_uniform(rGen_m);
|
|
if(P2 < 0.0047) {
|
|
double P3 = rand_uniform(s);//curand_uniform_double(&state);//gsl_rng_uniform(rGen_m);
|
|
double thetaru = 2.5 * sqrt(1 / P3) * sqrt(2.0) * theta0;
|
|
double P4 = rand_uniform(s);//curand_uniform_double(&state);//gsl_rng_uniform(rGen_m);
|
|
if(P4 > 0.5)
|
|
thetaru = -thetaru;
|
|
coord = 0; // no change in coordinates but one in momenta-direction
|
|
Rot(P, R, xplane, normP, thetaru, deltas, coord, par);
|
|
}
|
|
|
|
// y-direction: See Physical Review, "Multiple Scattering"
|
|
z1 = rand_normal(s, 0, 1);//curand_normal_double(&state);//gsl_ran_gaussian(rGen_m,1.0);
|
|
z2 = rand_normal(s, 0, 1);//curand_normal_double(&state);//gsl_ran_gaussian(rGen_m,1.0);
|
|
thetacou = z2 * theta0;
|
|
|
|
while(fabs(thetacou) > 3.5 * theta0) {
|
|
z1 = rand_normal(s, 0, 1);//curand_normal_double(&state);//gsl_ran_gaussian(rGen_m,1.0);
|
|
z2 = rand_normal(s, 0, 1);//curand_normal_double(&state);//gsl_ran_gaussian(rGen_m,1.0);
|
|
thetacou = z2 * theta0;
|
|
}
|
|
|
|
double yplane = z1 * deltas * theta0 / sqrt(12.0) + z2 * deltas * theta0 / 2.0;
|
|
coord = 2;
|
|
Rot(P, R, yplane, normP, thetacou, deltas, coord, par);
|
|
|
|
P2 = rand_uniform(s);//curand_uniform_double(&state);//gsl_rng_uniform(rGen_m);
|
|
if(P2 < 0.0047) {
|
|
double P3 = rand_uniform(s);//curand_uniform_double(&state);//gsl_rng_uniform(rGen_m);
|
|
double thetaru = 2.5 * sqrt(1 / P3) * sqrt(2.0) * theta0;
|
|
double P4 = rand_uniform(s);//curand_uniform_double(&state);//gsl_rng_uniform(rGen_m);
|
|
if(P4 > 0.5)
|
|
thetaru = -thetaru;
|
|
coord = 0; // no change in coordinates but one in momenta-direction
|
|
Rot(P, R, yplane, normP, thetaru, deltas, coord, par);
|
|
}
|
|
|
|
}
|
|
|
|
#define NUMPARAMS 19
|
|
__kernel void kernelCollimatorPhysics(__global PART *data, __global double *par,
|
|
__global RNDState *state, int numparticles,
|
|
__local double *p)
|
|
{
|
|
|
|
//get global id
|
|
int tid = get_local_id(0);
|
|
int idx = get_global_id(0);
|
|
|
|
printf("idx:\n");//, idx);
|
|
|
|
//transfer params to local memory
|
|
if (tid < NUMPARAMS)
|
|
p[tid] = par[tid];
|
|
|
|
barrier(CLK_LOCAL_MEM_FENCE);
|
|
|
|
RNDState s;
|
|
double3 R, P;
|
|
int l = 0;
|
|
if (idx < numparticles) {
|
|
R = data[idx].Rincol;
|
|
P = data[idx].Pincol;
|
|
s = state[idx];
|
|
}
|
|
|
|
double sq = sqrt(1.0 + Dot(P, P));
|
|
bool pdead = false;
|
|
bool hit = checkHit(R.z, p[POSITION], p[ZSIZE]);
|
|
double Eng;
|
|
|
|
if (hit) {
|
|
Eng = (sq - 1) * p[M_P];
|
|
energyLoss(&Eng, &pdead, p[DT_M], &s, p);
|
|
} else {
|
|
R.x = R.x + p[DT_M] * p[C] * P.x / sq;
|
|
R.y = R.y + p[DT_M] * p[C] * P.y / sq;
|
|
R.z = R.z + p[DT_M] * p[C] * P.z / sq;
|
|
l = -2;
|
|
}
|
|
|
|
if (hit && !pdead) {
|
|
double ptot = sqrt((p[M_P] + Eng) * (p[M_P] + Eng) - (p[M_P] * p[M_P])) / p[M_P];
|
|
sq = sqrt(Dot(P, P));
|
|
P.x = P.x * ptot / sq;
|
|
P.y = P.y * ptot / sq;
|
|
P.z = P.z * ptot / sq;
|
|
coulombScat(&R, &P, p[DT_M], &s, p);
|
|
}
|
|
|
|
if (hit && pdead)
|
|
l = -1;
|
|
|
|
if (idx < numparticles) {
|
|
data[idx].Rincol = R;
|
|
data[idx].Pincol = P;
|
|
data[idx].label = l;
|
|
state[idx] = s;
|
|
}
|
|
|
|
}
|
|
|
|
|
|
/* count dead particles and particles leaving material - boost compute? */
|
|
|
|
/* sort particles so dead and leaving particles are at the end of PART array - boost compute */
|
|
|
|
|