%uses parameter for actual detector and Marinelli beaker geometery
%adding mass attenuation coefficient for energy = 1500eV for detector and
%soil sample
%assuming the soil sample to be similar composition to that of the
%silicondioxide
%assuming the density of the soil is half of that of silicondioxide
% theoritical calculation is in cm
% detector diameter is 3.89 cm and the height is 4.72cm
% the base of the cylinder is at (0,0,0)
%calculates the total hits and the solid angle
% %draw Simple Cylinder and checks for hits for an individual trajectory
r=1.9450; % inner radius of the beaker in cm
thickness=3; % thickness of the sample
Hits=0; % keeps track of number of gamma that did not scattered in the detector
i=1;
n=100000;
PotentialHits=0; % keeps track of how many times it fulfilled condition 3
k=0; % keeps track of number of gammas that did not scattered in the sample
ps=0;
cs=0;
pd=0;
cd=0;
DistSample=0;
DistDetector=0;
while i<=n
%getting random points
theta1=randi([0,359])+rand(); % problem theta is never going to be zero or 360
z1=(randi([0,649])+rand())/100; % in cm units. height of the sample
R1=4.25+thickness*(rand())^.5; %creates random uniform distribution between the range of the thickness of the sample
x1=R1*cosd(theta1);
y1=R1*sind(theta1);
%getting random trajectory
R2=15; % R2 is the length of the trajectory
theta= acosd(1-2*rand()); %uniform theta distribution in degrees. REMEBER this is in FLOATING POINT
Phi= randi([0,359])+rand(); %uniform phi distribution in degrees. REMEBER this is in FLOATING POINT. also Phi will never be zero or 360
x2p= R2*sind(theta)*cosd(Phi); % measured in the source coordinate
y2p= R2*sind(theta)*sind(Phi); %measured in the source coordinate
z2p= R2*cosd(theta); %measured in the source coordinate
x2=x2p+x1; % measured in the cylinder coordinate
y2=y2p+y1; % measured in the cylinder coordinate
z2=z2p+z1; % measured in the cylinder coordinate
%counting hits
Condition=0; %initialization
%checking for exception conditions where Phi blows up
if (ceil(Phi*10000)/10000== 90.0000) && (ceil(theta*10000)/10000 == 90.0000 ) % comparing up to four sig fig
if z1>=1.78 && z1 <=6.50
xa=x1;
xb=x1;
xSource=xa-x1; % always zero
%check for the absolute value of x
if abs(xa)== r
%tangent condition so does not matter which direction the
%trajectory goes. But we are not going to count this as a hit.
Condition=0;
elseif abs(xa)<r
% intersects the circle at two points
Condition=1;
%since xaa=x1 the y of the source and the cylinder coordinate is
%the same
ya1=(r^2-(xa)^2)^.5;
yb1=-(r^2-(xb)^2)^.5;
ySourceCoa1=ya1-y1; % converting into source coordinate
ySourceCoa2=yb1-y1;
%check for direction
% for this condition the aPhiPrime1 and aPhiPrime3 will be positive
% 90 no matter what.
PhiPrime1= atand(ySourceCoa1/xSource)*1.0000001; %finds the Phi in the source coordinate and multipling by 1.0000001 to make it a floating point number
PhiPrime3= atand(ySourceCoa2/xSource)*1.0000001;
%for PhiPrime1
if xa>=0 && ySourceCoa1>=0
PhiPrime1=PhiPrime1+0;
elseif xa<0 && ySourceCoa1>=0
PhiPrime1=PhiPrime1 +180;
elseif xa<0 && ySourceCoa1<0
PhiPrime1=PhiPrime1+180;
else
PhiPrime1=PhiPrime1+360;
end
%for PhiPrime3
if xa>=0 && ySourceCoa2>=0
PhiPrime3=PhiPrime3+0;
elseif xa<0 && ySourceCoa2>=0
PhiPrime3=PhiPrime3+180;
elseif xa<0 && ySourceCoa2<0
PhiPrime3=PhiPrime3+180;
else
PhiPrime3=PhiPrime3+360;
end
% rounding to four sig digits
Phi=ceil(Phi*10000)/10000;
PhiPrime1=ceil(PhiPrime1*10000)/10000;
PhiPrime3=ceil(PhiPrime3*10000)/10000;
%check if the PhiPrime == Phi
if (PhiPrime1==Phi) && (PhiPrime3==Phi)
Condition=3;
PhiPrime2=0;
PhiPrime4=0;
end
else
% it does not intersect the circle at all
Condition=0;
end
end
else
%getting x and y projection
%r= radius of the cylinder
a=1+(tand(Phi))^2;
b=2*tand(Phi)*(y1-x1*tand(Phi));
c=(x1^2)*(tand(Phi))^2-2*x1*y1*tand(Phi)+(y1^2)-r^2;
xa= (-b+((b^2)-4*a*c)^.5)/(2*a);
xb= (-b-((b^2)-4*a*c)^.5)/(2*a);
if (isreal(xa)&& isreal(xb))==1
Condition=1;
%Checking for condition 2 (if it hits the detector)
xSourceCo1=xa-x1; %finding new value of x prime in the source coordinate
xSourceCo2=xb-x1;
zSourceCo1=(xSourceCo1 * cosd(theta))/(sind(theta)*cosd(Phi)); % finding value of z prime in the source coordinate
zSourceCo2=(xSourceCo2 * cosd(theta))/(sind(theta)*cosd(Phi));
zCylinderCo1=zSourceCo1+z1; %getting the value of z prime in the cylinder coordinate
zCylinderCo2=zSourceCo2+z1;
if (zCylinderCo1>=1.78 && zCylinderCo1<= 6.5) || (zCylinderCo2>=1.78 && zCylinderCo2<=6.5) %checking if the z value found is greater then the height of the cylinder in cm
Condition=2;
%checking for condition 3 (direction of the trajectory)
%getting all the values of y
ya1=((r^2)-(xa^2))^.5; % finds the y when the trajectory intersects the cylinder
ya2=-ya1;
yb1=((r^2)-(xb^2))^.5;
yb2=-yb1;
ySourceCo11=ya1-y1; % converting into source coordinate
ySourceCo12=ya2-y1;
ySourceCo21=yb1-y1;
ySourceCo22=yb2-y1;
%getting all the values of PhiPrime
PhiPrime11= atand(ySourceCo11/xSourceCo1); %finds the Phi in the source coordinate
PhiPrime12= atand(ySourceCo12/xSourceCo1);
PhiPrime21= atand(ySourceCo21/xSourceCo2);
PhiPrime22= atand(ySourceCo22/xSourceCo2);
%checking for the diffrent quadrent condition
%for PhiPrime 1
if xSourceCo1>=0 && ySourceCo11>=0
PhiPrime1=PhiPrime11;
elseif xSourceCo1<0 && ySourceCo11>=0
PhiPrime1=PhiPrime11+180;
elseif xSourceCo1<0 && ySourceCo11<0
PhiPrime1=PhiPrime11+180;
else
PhiPrime1=PhiPrime11+360;
end
%for PhiPrime 2
if xSourceCo1>=0 && ySourceCo12>=0
PhiPrime2=PhiPrime12;
elseif xSourceCo1<0 && ySourceCo12>=0
PhiPrime2=PhiPrime12+180;
elseif xSourceCo1<0 && ySourceCo12<0
PhiPrime2=PhiPrime12+180;
else
PhiPrime2=PhiPrime12+360;
end
%for PhiPrime 3
if xSourceCo2>=0 && ySourceCo21>=0
PhiPrime3=PhiPrime21;
elseif xSourceCo2<0 && ySourceCo21>=0
PhiPrime3=PhiPrime21+180;
elseif xSourceCo2<0 && ySourceCo21<0
PhiPrime3=PhiPrime21+180;
else
PhiPrime3=PhiPrime21+360;
end
%for PhiPrime 4
if xSourceCo2>=0 && ySourceCo22>=0
PhiPrime4=PhiPrime22;
elseif xSourceCo2<0 && ySourceCo22>=0
PhiPrime4=PhiPrime22+180;
elseif xSourceCo2<0 && ySourceCo22<0
PhiPrime4=PhiPrime22+180;
else
PhiPrime4=PhiPrime22+360;
end
Phi=ceil(Phi*10000)/10000;
PhiPrime1=ceil(PhiPrime1*10000)/10000;
PhiPrime2=ceil(PhiPrime2*10000)/10000;
PhiPrime3=ceil(PhiPrime3*10000)/10000;
PhiPrime4=ceil(PhiPrime4*10000)/10000;
%checking if the trajectory is in the right direction
if (PhiPrime1==Phi) || (PhiPrime2==Phi) || (PhiPrime3==Phi) || (PhiPrime4==Phi)
Condition=3;
end
end
end
end
% finding the two right PhiPrime to get the values of y and
% x
% there is two correct value for PhiPrime=Phi. One for xa and
% other for xb
if Condition==3
PotentialHits=PotentialHits+1; % keeps track of how many times it fulfilled condition 3
% PhiPrime1==PhiPrime3 and PhiPrime2==PhiPrime4 so
% checking for anyone of them is fine
if Phi==PhiPrime1
yf1=ya1;
zf1=zCylinderCo1;
end
if Phi==PhiPrime2
yf1=yb1;
zf1=zCylinderCo1;
end
if Phi==PhiPrime3
yf2=yb1;
zf2=zCylinderCo2;
end
if Phi==PhiPrime4
yf2=yb2;
zf2=zCylinderCo2;
end
% finding the distance between the intersection points
% and the source point
distanceSample1 =((xa-x1)^2 + (yf1-y1)^2 + (zf1-z1)^2)^.5; %distance from the source to the first intersection point in the cylinder
distanceSample2 =((xb-x1)^2 + (yf2-y1)^2+ (zf2-z1)^2)^.5; %distance from the source to the second intersection point in the cylinder
% Finds which intersection point is the first by
% comparing the distance from the source to the two
% intersection points
if distanceSample2>distanceSample1
FirstIntersectionDistance= distanceSample1; %subustracting the empty space between the beaker and the detector, whicn is atleast 2.3
else
FirstIntersectionDistance= distanceSample2;
end
DistSample=DistSample+FirstIntersectionDistance;
%attenuation due to compton scattering
Scs=-(log(1-rand()))*14.8960; %14.56=1/attenuation length(units=cm). For sample
Sj=-(log(1-rand()))*86586.0838; % 86586.08 = 1/attenuation length (units=cm). For sample.
PhotoSample=0;
ComptonSample=0;
if Scs>Sj %Photoelectric effect
ps=ps+1;
if FirstIntersectionDistance < Scs && FirstIntersectionDistance >= Sj
PhotoSample=1;
ComptonSample=0;
end
if FirstIntersectionDistance < Scs && FirstIntersectionDistance <Sj
PhotoSample=0;
ComptonSample=0;
end
if FirstIntersectionDistance > Scs
PhotoSample=1;
ComptonSample=0;
end
end
if Sj>Scs %Compton Scatter
cs=cs+1;
if FirstIntersectionDistance < Sj && FirstIntersectionDistance > Scs
PhotoSample=0;
ComptonSample=1;
end
if FirstIntersectionDistance < Scs && FirstIntersectionDistance <Sj
PhotoSample=0;
ComptonSample=0;
end
if FirstIntersectionDistance > Sj
PhotoSample=0;
ComptonSample=1;
end
end
if PhotoSample==0 && ComptonSample==0 %survives the attenuation in the detector
k=k+1;
distance =((xa-xb)^2 + (yf1-yf2)^2 + (zf1-zf2)^2)^.5; % distance travelled inside the detector
DistDetector=DistDetector+distance;
Sc=-(log(1-rand()))*4.12; %4.12=1/attenuation length(units=cm). For detector
Si=-(log(1-rand()))*570.58; % 243.76 = 1/attenuation length(units=cm). For detector
PhotoDetector=0;
ComptonDetector=0;
if Sc>Si %Photoelectric effect
pd=pd+1;
if distance < Sc && distance > Si
PhotoDetector=1;
ComptonDetector=0;
end
if distance < Sc && distance <Si
PhotoDetector=0;
ComptonDetector=0;
end
if distance > Sc
PhotoDetector=1;
ComptonDetector=0;
end
end
if Si>Sc %Compton Scatter
cd=cd+1;
if distance < Si && distance > Sc
PhotoDetector=0;
ComptonDetector=1;
Si2=-(log(1-rand()))*13.06;
if (distance-Sc)>Si2
PhotoDetector=1;
end
end
if distance < Sc && distance <Si
PhotoDetector=0;
ComptonDetector=0;
end
if distance > Si
PhotoDetector=0;
ComptonDetector=1;
if (distance-Sc)>Si
PhotoDetector=1;
end
end
end
%Checking if the ray satisfied both conditions
if PhotoDetector==1
Hits= Hits+1; %number of gammas that did not scattered in the detector
end
end
end
i=i+1;
end
PotentialHits
AvgDistSample=DistSample/PotentialHits
ps
cs
k
AvgDistDetector=DistDetector/k
pd
cd
Hits
%AvgtotalSc=totalSc/Indetector
% SamplePerCompton=((PotentialHits-NotCompton)/PotentialHits)*100
% SamplePerPhotoelectric=((NotCompton-NotPhotoelectric)/PotentialHits)*100
% DetectorPerCompton=((NotPhotoelectric-k)/NotPhotoelectric)*100
% DetectorPerPhotoelectric=(Hits/NotPhotoelectric)*100
%meanSi=totalSi/k;
%meanSj=totalSj/j;
%meanSc=totalSc/c;
%AvegDistance=totalDistance/k;
Efficiency=Hits/n
SolidAngle=PotentialHits/n
%fprintf('____________ \ntotal successful hits = %d \nTotal efficiency = %f \nAverage distance travelled by the trajectory inside the detector = %f\nMean scatter length in detector = %f\nMean Scatter length in sample = %f\nPercentage of gamma rays absorbed in the sample= %d\nPercentage of gamma rays absorbed the detector = %d\nSolid angle = %f\n',Hits, Efficiency,AvegDistance,meanSi,meanSj,((j-NotPhotoelectric)/j)*100,(Hits/NotPhotoelectric)*100,j/n)