Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Bring ViewModISIS up to date with McStas tree? #1

Open
wants to merge 2 commits into
base: master
Choose a base branch
from
Open

Bring ViewModISIS up to date with McStas tree? #1

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
bfb2f46
Name change
willend Apr 26, 2018
2ca1ad4
Improvements on docs, input file handling etc. from McStas tree
willend Apr 26, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
152 changes: 82 additions & 70 deletions ViewModISISver1.comp → ViewModISIS.comp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@
* Risoe National Laboratory, Roskilde, Denmark
* Institut Laue Langevin, Grenoble, France
*
* Component: ViewModISISver1
* Component: ViewModISIS
*
*
* %I
Expand All @@ -24,56 +24,45 @@
* The Face argument determines which TS1 or TS2 beamline is to be sampled by using corresponding Etable file.
* Neutrons are created having a range of energies determined by the E0 and E1 arguments.
* Trajectories are produced such that they pass through the shutter front face (RECOMENDED) or moderator face (defined by
* modXsize and modZsize) and a focusing rectangle (defined by xw, yh and dist).
* xw and yh) and a focusing rectangle (defined by focus_xw, focus_yh and dist).
*
*
* %P
* INPUT PARAMETERS:
*
* Face: (word) Etable filename
* Face: [string] Etable filename
* E0: [meV] Lower edge of energy distribution
* E1: [meV] Upper edge of energy distribution
* modPosition: [int] Defines the initial surface for neutron distribution production. Possible values = 0 (at moderator) or 1 (at shutter insert)
* xw: [m] Moderator width
* yh: [m] Moderator height
* focus_xw: [m] Width of focusing rectangle
* focus_yh: [m] Height of focusing rectangle
* dist: [m] Distance from source surface to the focusing rectangle
* verbose: [int] Flag to output debugging information
* beamcurrent [uA] ISIS beam current
*
* E0: (meV) Lower edge of energy distribution
*
* E1: (meV) Upper edge of energy distribution
*
* tally: (m) Distance from shutter face to point tally in full simulation (TEMPORARY TESTING PARAMETER ONLY)
* - not used in this version of component -
*
* modPosition: Defines the initial surface for neutron distribution production
* Possible values = 0 or 1
* 0 = "moderator"; 1 (RECOMENDED) = front face of shutter insert
*
* modXsize: (m) Moderator width
*
* modZsize: (m) Moderator height
*
* xw: (m) Width of focusing rectangle
*
* yh: (m) Height of focusing rectangle
*
* dist: (m) Distance from source surface to the focusing rectangle
*
*
* Example 1: ViewModISISver1(Face="msMerlin_160.mcstas", E0 = E_min, E1 = E_max,
* tally=8.4, modPosition=0, modXsize=0.12, modZsize = 0.12, xw = 0.094, yh = 0.094, dist = 1.6)
* Example 1: ViewModISISver1(Face="TS1_S04_Merlin.mcstas", E0 = E_min, E1 = E_max,
* dist = 1.7, focus_xw = 0.094, focus_yh = 0.094, modPosition=0,
* xw=0.12,yh=0.12)
*
* MERLIN simulation.
* modPosition=0 so initial surface is at (near to) the moderator face.
* In this example, xw and yh are chosen to be identical to the shutter opening dimension.
* dist = 1.6 is the 'real' distance to the shutter front face => This is TimeOffset value (=160 [cm])
* from msMerlin_160.mcstas file.
* In this example, focus_xw and focus_yh are chosen to be identical to the shutter opening dimension.
* dist = 1.7 is the 'real' distance to the shutter front face => This is TimeOffset value (=170 [cm])
* from TS1_S04_Merlin.mcstas file.
*
*
* Example 2: ViewModISISver1(Face="Let_timeTest_155.mcstas", E0 = E_min, E1 = E_max,
* tally=8.4, modPosition=1, modXsize=0.196, modZsize = 0.12, xw = 0.04, yh = 0.094, dist = 0.5)
* modPosition=1, xw=0.196, yh = 0.12, focus_xw = 0.04, focus_yh = 0.094, dist = 0.5)
*
* LET simulation.
* modPosition=1 so initial surface is at front face of shutter insert.
*
* (IMPORTANT) If modPosition=1, the modXsize and modZsize values are obsolete. The dimensions of
* (IMPORTANT) If modPosition=1, the xw and yh values are obsolete. The dimensions of
* initial surface are automatically calculated using "RDUM" values in Let_timeTest_155.mcstas file.
*
* In this example, xw and yh are arbitrary chosen to be identical to the shutter opening dimension.
* In this example, focus_xw and focus_yh are arbitrary chosen to be identical to the shutter opening dimension.
*
*
*
Expand All @@ -83,8 +72,8 @@
* %E
*******************************************************************************/

DEFINE COMPONENT ViewModISISver1
SETTING PARAMETERS (string Face,E0, E1,tally,modPosition,modXsize,modZsize,xw,yh,dist)
DEFINE COMPONENT ViewModISIS
SETTING PARAMETERS (string Face="TS1_S04_Merlin.mcstas",E0, E1,modPosition=0,xw=0.12,yh=0.12,focus_xw=0.094,focus_yh=0.094,dist=1.7,int verbose=0, beamcurrent=1)
DECLARE
%{

Expand Down Expand Up @@ -589,8 +578,8 @@ processHeader(FILE* TFile)
TS.XAxis=sqrt(TS.XAxis)/100.0;
if (!modPosition)
{
TS.ZAxis=modZsize;
TS.XAxis=modXsize;
TS.ZAxis=yh;
TS.XAxis=xw;
}
fprintf(stderr,"Time off sec == %g \n",TS.timeOffset);
fprintf(stderr,"mod size z == %g \n",TS.ZAxis);
Expand Down Expand Up @@ -696,7 +685,7 @@ readHtable(FILE* TFile,const double Einit,const double Eend)
TS.EInt=(double*) malloc(TS.nEnergy*sizeof(double)); // Runing integral of energy
TS.TimeBin=(double*) malloc(TS.nTime*sizeof(double));
TS.EnergyBin=(double*) malloc(TS.nEnergy*sizeof(double));
fprintf(stderr,"NUMBER %d %d \n",TS.nEnergy,TS.nTime);
if (verbose) fprintf(stderr,"NUMBER %d %d \n",TS.nEnergy,TS.nTime);

Tsum=0.0;
Ea=0.0;
Expand Down Expand Up @@ -735,8 +724,8 @@ readHtable(FILE* TFile,const double Einit,const double Eend)
{
Ftime=0;
TS.EInt[eIndex+1]=Tsum;

fprintf(stderr,"Energy %g %g \n",EtoLambda(Ea),(TS.EInt[eIndex+1]-TS.EInt[eIndex])*3.744905847e14);
if (verbose) fprintf(stderr,"Energy %g %g \n",EtoLambda(Ea),(TS.EInt[eIndex+1]-TS.EInt[eIndex])*3.744905847e14);
}
else
{
Expand Down Expand Up @@ -938,43 +927,66 @@ openFile(char* FileName)
strcpy(extFileName,FileName);
strcpy(extFileName+fLen,".mcstas");

/* Now check for the requested moderator file. In terms of precedence,
1) Use MCTABLES location if available
2) Is a local file available in PWD?
3) Is there an ISIS_tables in PWD?
4) Is the file available from the MCSTAS/contrib/ISIS_tables folder?
- otherwise fail */

/* Is the file located in working dir? */
fprintf(stderr,"Searching for %s in multiple locations... -- \n",FileName);

/* 1) Is MCTABLES set and file located there? */
if (getenv("MCTABLES"))
{
strcpy(mct, getenv("MCTABLES"));
sprintf(testFileName, "%s%c%s", mct, MC_PATHSEP_C, FileName);
efile=fopen(testFileName,"r");
if (efile) {
fprintf(stderr," - Found in MCTABLES folder %s!\n",mct);
return efile;
}
}

/* 2) Is the file located in working dir? */
efile=fopen(FileName,"r");
if (efile) return efile;
if (efile) {
fprintf(stderr," - Found in current directory!\n");
return efile;
}

efile=fopen(extFileName,"r");
if (efile) return efile;

/* Is the file in a local 'tablefiles' folder? */
sprintf(testFileName,"%c%s%s", MC_PATHSEP_C,"ISIS_tables",FileName);
/* 3) Is the file in a local 'tablefiles' folder? */
sprintf(testFileName,"%s%c%s","ISIS_tables",MC_PATHSEP_C,FileName);
efile=fopen(testFileName,"r");
if (efile) return efile;

if (efile) {
fprintf(stderr," - Found in local ISIS_tables directory!\n");
return efile;
}

/* Is MCTABLES set, files located there? */
if (getenv("MCTABLES"))
{
strcpy(mct, getenv("MCTABLES"));
sprintf(testFileName, "%s%c%s", mct, MC_PATHSEP_C, FileName);
efile=fopen(testFileName,"r");
if (efile) return efile;
}
/* 4) Is the file available within the MCSTAS install dir? */
sprintf(testFileName,"%s%c%s%c%s%c%s",MCSTAS,MC_PATHSEP_C,"contrib",MC_PATHSEP_C,"ISIS_tables",MC_PATHSEP_C,FileName);
efile=fopen(testFileName,"r");
if (efile) {
fprintf(stderr," - Found in MCSTAS system dir: \n %s%c%s%c%s%\n",MCSTAS,MC_PATHSEP_C,"contrib",MC_PATHSEP_C,"ISIS_tables");
return efile;
}

fprintf(stderr,"Searching -- %s, MCTABLES directory\n",testFileName);
fprintf(stderr,"\nPlease check your McStas installation and/or MCTABLES environment variable!\n");
/* If we reach here, the file was not found - raise fatal error */
fprintf(stderr,"% - Not found! \nPlease check your McStas installation and/or MCTABLES environment variable!\n",FileName);
exit(1);
return efile;
}

double strArea()
{
/*
Returns the mean Str view of the viewport
This integrates over each point on the window xw to yh
This integrates over each point on the window focus_xw to focus_yh
View port is symmetric so use only 1/4 of the view
for the calcuation.
Control Values TS.XAxis TS.ZAxis xw yh
Control Values TS.XAxis TS.ZAxis focus_xw focus_yh
*/

double A;
Expand All @@ -989,7 +1001,7 @@ double strArea()
D2=dist*dist;
A=0.0;
fprintf(stderr,"TS axis == %g %g\n",TS.XAxis,TS.ZAxis);
fprintf(stderr,"AW axis == %g %g\n",xw,yh);
fprintf(stderr,"AW axis == %g %g\n",focus_xw,focus_yh);
for(i=0;i<NStep;i++) // Mod X
{
Mx=(i+0.5)*TS.XAxis/(NStep*2.0);
Expand All @@ -1000,27 +1012,27 @@ double strArea()
for(aa= -NStep;aa<=NStep;aa++) //view port
for(bb= -NStep;bb<=NStep;bb++)
{
Vx=aa*xw/(2.0*NStep+1.0);
Vy=bb*yh/(2.0*NStep+1.0);
Vx=aa*focus_xw/(2.0*NStep+1.0);
Vy=bb*focus_yh/(2.0*NStep+1.0);
A+=1.0/((Mx-Vx)*(Mx-Vx)+(My-Vy)*(My-Vy)+D2);
}
}
}
//change to Mx*My
A*= (TS.XAxis*TS.ZAxis)/(NStep*NStep*(2.0*NStep+1.0)*(2.0*NStep+1.0));
// Correct for the area of the viewport. (tables are per cm^2)
A*=xw*yh*10000;
A*=focus_xw*focus_yh*10000;

// testing only - Goran
// if (!modPosition)
// {
// projArea=xw*yh/tally/tally*(tally+dist)*(tally+dist);
// projArea=focus_xw*focus_yh/tally/tally*(tally+dist)*(tally+dist);
// A*=TS.XAxis*TS.ZAxis/(TS.XAxis*TS.ZAxis-projArea);
// }

fprintf(stderr,"Viewport == %g %g Moderator size == (%g * %g) m^2 \n",xw,yh,TS.XAxis,TS.ZAxis);
fprintf(stderr,"Viewport == %g %g Moderator size == (%g * %g) m^2 \n",focus_xw,focus_yh,TS.XAxis,TS.ZAxis);
fprintf(stderr,"Dist == %g (metres) \n",dist);
fprintf(stderr,"Viewport Solid angle == %g str\n",A/(xw*yh*10000.0));
fprintf(stderr,"Viewport Solid angle == %g str\n",A/(focus_xw*focus_yh*10000.0));
fprintf(stderr,"Solid angle used == %g str\n",A);
return A;
}
Expand Down Expand Up @@ -1076,8 +1088,8 @@ TRACE

// if (modPosition) z=TS.rdumMid; // Goran

xf = xw*(0.5-rand01()); /* Choose focusing position uniformly */
yf = yh*(0.5-rand01()); /* Choose focusing position uniformly */
xf = focus_xw*(0.5-rand01()); /* Choose focusing position uniformly */
yf = focus_yh*(0.5-rand01()); /* Choose focusing position uniformly */

getPoint(&Tval,&Eval,&rtE0,&rtE1);

Expand All @@ -1101,13 +1113,13 @@ TRACE
t+=TS.rdumMid/vz;
}

p=angleArea*Ival/Nsim;
p=beamcurrent*angleArea*Ival/Nsim;

// testing only - Goran
// if (!modPosition) p*=TS.timeOffset;


if (!(Ncount % 100000))
if (verbose && !(Ncount % 100000))
{
fprintf(stderr,"FPos[%d]=> %g %g %g \n",Ncount,x,y,z);
fprintf(stderr,"FDir[%d]=> %g %g %g \n",Ncount,vx,vy,vz);
Expand Down