/*******************************************************************************
* Instrument: ESS_butterfly_MCPL_test
*
* %I
* Written by: Peter Willendrup <pkwi@fysik.dtu.dk>
* Date: 2016-09-26
* Origin: ESS
* %INSTRUMENT_SITE: ESS
*
* Test instrument for the updated BF1 butterfly moderator design using MCPL input
*
* %D
* Test instrument for the updated BF1 butterfly moderator design using MCPL input files.
* 
* The below example gives a 50-50 (statistics-wise) cold/thermal beam at beamline N10.
* Example: mcrun -n1e6 ESS_butterfly_MCPL_test.instr sector=N beamline=10 cold=0.5
*
* Assumes access to binary MCPL datasets in . named [sector][beamline].mcpl.gz, i.e. W8.mcpl.gz.
*
* %P
* sector:        [str]  Defines the 'sector' of your instrument position. Valid values are "N","S","E" and "W"
* beamline:      [1]    Defines the 'beamline number' of your instrument position. Valid values are 1..10 or 1..11 depending on sector
* Lmin:          [AA]   Minimum wavelength simulated
* Lmax:          [AA]   Maximum wavelength simulated
* c_performance: [1]    Cold brilliance scalar performance multiplicator c_performance > 0
* t_performance: [1]    Thermal brilliance scalar performance multiplicator t_performance > 0
* index:         [1]    Target index for source focusing. Defaults to illuminate the "cold collimated" brilliance monitor, thereby suppressing "dist"
* dist:          [m]    Distance from origin to focusing rectangle; at (0,0,dist) - alternatively use target_index
* cold:          [1]    Defines the statistical fraction of events emitted from the cold part of the moderator
* Yheight:       [m]    Defines the moderator height. Valid values are 0.03 m and 0.06 m
* delta:         [m]    Parameter that allows to scan "collimator" position
*
* %L
* <reference/HTML link>
* Benchmarking website available at <a href="http://ess_butterfly.mcstas.org">http://ess_butterfly.mcstas.org</a>
* %E
*******************************************************************************/
DEFINE INSTRUMENT ESS_butterfly_test(string sector="N",
  int beamline=1,Lmin=0.2,Lmax=20,c_performance=1,t_performance=1,
  int index=8,dist=0,cold=1,Yheight=0.03,delta=0,thres=0.003,filter=0,repeat=1,
  E_smear=0.1,pos_smear=0.01,dir_smear=0.01)

DECLARE
%{
  int IsCold;
  double SrcX, SrcY, SrcZ;
  double Theta;
  double MonTransl;
  double XW, YH;
  char options1[256],options2[256],options3[256],options4[256];
  char srcdef[128];
  double WidthC=0.072,WidthT=0.108;
  double WL;
  double lambdamin, lambdamax;
  double CCold,CThermal;
  double SurfSign;
  double TCollmin;
  double TCollmax;
  double Emin,Emax;
  double EminTh=20, EmaxTh=100, EminC=0, EmaxC=20;
  double Eneutron;
  /* 10 beamlines in sector N and E  - plus one location added for drawing */
  double iBeamlinesN[] = { 30.0,  36.0,  42.0,  48.0,  54.0,  60.0,  66.0,  72.0,  78.0,  84.0,  90.0};
  double iBeamlinesE[] = {-30.0, -36.0, -42.0, -48.0, -54.0, -60.0, -66.0, -72.0, -78.0, -84.0, -90.0};
  /* 11 beamlines in sector S and W - plus one location added for drawing */
  double iBeamlinesW[] = { 150.0,  144.7,  138.0,  132.7,  126.0,  120.7,  114.0,  108.7,  102.0,  96.7,  90.0,  84.0};
  double iBeamlinesS[] = {-150.0, -144.7, -138.0, -132.7, -126.0, -120.7, -114.0, -108.7, -102.0, -96.7, -90.0, -84.0};
  double* iBeamlines;
  double ANGLE;
  double DeltaX,DeltaZ;
  char MCPLfile[128];
  double T0,L0;

  double calcAlpha(double length, double radius) {
    // calculate angle of arm after curved guide
    return RAD2DEG * length/radius;
  }	
  
  double calcX(double length, double radius) {
    // calculate position and angle of arm after curved guide
    double alpha = DEG2RAD * calcAlpha(length, radius);
    return radius*(1.0-cos(alpha));
  }	
  
  double calcZ(double length, double radius) {
    // calculate position and angle of arm after curved guide
    double alpha = DEG2RAD * calcAlpha(length, radius);
    return radius*sin(alpha);
  }
%}

INITIALIZE
%{
  lambdamin=Lmin; 
  lambdamax=Lmax;
  XW=1.05*(WidthC+2*WidthT);
  YH=1.05*Yheight;
  sprintf(options1,"user1 bins=201 limits=[-%g,%g]",XW/2,XW/2);
  sprintf(options4,"user1 bins=201 limits=[-%g,%g]",YH/2,YH/2);
  sprintf(options2,"user1 bins=201 limits=[-%g,%g], user2 bins=201 limits=[-%g,%g]",XW/2,XW/2,YH/2,YH/2);
  sprintf(options3,"user1 bins=201 limits=[-%g,%g], user2 bins=201 limits=[-%g,%g]",1.05*(WidthC/2),1.05*(WidthC/2),1.05*Yheight/2,1.05*Yheight/2);
  sprintf(srcdef,"2015");
  if (beamline==1) {
    TCollmin=0;
    TCollmax=0.058;
  } else if (beamline==2) {
    TCollmin=0;
    TCollmax=0.06;
  }
  else {
    TCollmin=0.011;
    TCollmax=0.071;
  }
  if (strcasestr(sector,"N")) {
    iBeamlines=iBeamlinesN;
    DeltaX=-0.0585; DeltaZ=0.0925;
  } else if (strcasestr(sector,"W")) {
    iBeamlines=iBeamlinesW;
    DeltaX=0.0585; DeltaZ=0.0925;
  } else if (strcasestr(sector,"S")) {
    iBeamlines=iBeamlinesS;
    DeltaX=0.0585; DeltaZ=-0.0925;
  } else if (strcasestr(sector,"E")) {
    iBeamlines=iBeamlinesE;
    DeltaX=-0.0585; DeltaZ=-0.0925;
  }
  ANGLE=iBeamlines[beamline-1]-90;
  if (filter==0)
    sprintf(MCPLfile,"%s%i.mcpl.gz",sector,beamline);
  else
    sprintf(MCPLfile,"%s%i_filtered.mcpl.gz",sector,beamline);
  printf("MCPLfile is %s\n",MCPLfile);
%}

TRACE

COMPONENT Origin = Progress_bar()
AT (0, 0, 0) ABSOLUTE

/* read neutrons from an mcpl file*/

COMPONENT vinROT2 = Arm()
AT(0,0,0) RELATIVE PREVIOUS
  ROTATED (0,-90,0) RELATIVE PREVIOUS

COMPONENT vinROT1 = Arm()
AT(0,0,0) RELATIVE PREVIOUS
  ROTATED (-90,0,0) RELATIVE PREVIOUS

  COMPONENT vin = MCPL_input(filename=MCPLfile,verbose=1,repeat_count=repeat,E_smear=E_smear,pos_smear=pos_smear,dir_smear=dir_smear)
AT(0,0,0) RELATIVE PREVIOUS
EXTEND %{
  if(p>mcipthres) {ABSORB;}
  else {SCATTER;}
  p*=1.56e16;
  p/=1e5;
  z=z-0.137;
%}
 
COMPONENT Sphere1 = PSD_monitor_4PI(filename="nonrotated", radius=2.2,restore_neutron=1)
AT (0,0,0) RELATIVE PREVIOUS

/* Focusing for this use of the source is a little unphysical: 1x1cm @ 1m ~ 1e-4 steradian. To be useful in a "proper" instrument, you should of course illuminate your beamport fully!*/
COMPONENT Source = ESS_butterfly(sector=sector,beamline=beamline,Lmin=Lmin,Lmax=Lmax,c_performance=c_performance,t_performance=t_performance,dist=dist,target_index=index,cold_frac=cold, yheight=Yheight,
				   focus_xw=0.12, focus_yh=0.12)
  WHEN (0==1) AT (DeltaX,0,DeltaZ) ABSOLUTE
  ROTATED (0, ANGLE, 0) ABSOLUTE

COMPONENT Sphere0 = PSD_monitor_4PI(filename="rotated", radius=2.2,restore_neutron=1)
AT (0,0,0) RELATIVE Source

COMPONENT Focus_cut=Shape(xwidth=0.01,yheight=0.01)
  AT(0,0,2) RELATIVE Source
EXTEND %{
  double xtmp,ytmp,ztmp,vxtmp,vytmp,vztmp;
  xtmp=x;ytmp=y;ztmp=z;
  vxtmp=vx;vytmp=vy;vztmp=vz;
  ALLOW_BACKPROP;
  PROP_Z0;
  SCATTER;
  
  if (fabs(x)>0.06 || fabs(y)>0.06) {
    ABSORB;
  } else {
    x=xtmp;y=ytmp;z=ztmp;
    vx=vxtmp;vy=vytmp;vz=vztmp;
  }
%}

COMPONENT BackTrace = Shape(xwidth=0.3,yheight=0.3)
  AT (0,0,0.08) RELATIVE Source
EXTEND %{
  /* Propagate back to a small rectangle in front of moderators */

  ALLOW_BACKPROP;
  PROP_Z0;
  SCATTER;
  
  /* Remove neutrons that are not from around the moderators */
  if (fabs(x)>0.12 || fabs(y)>0.03) {
   ABSORB;
  }
  double myL = (2*PI/V2K)/sqrt(vx*vx + vy*vy + vz*vz);
  if (myL<mcipLmin || myL>mcipLmax) {
    ABSORB;
  }
  t+=rand01()*ESS_SOURCE_DURATION;
  if (t>3*ESS_SOURCE_DURATION) {
    ABSORB;
  }
  /* Measure location and energy for later use */
  SrcX=x;SrcY=y;SrcZ=z;
  Eneutron=VS2E*(vx*vx + vy*vy + vz*vz);
  if (Eneutron>EminTh) {
    Emin=EminC;Emax=EmaxC;
    IsCold=0;
  } else {
    Emin=EminTh;Emax=EmaxTh;
    IsCold=1;
  }
  T0=t;
  L0=myL;
%}


/* These arms are just to ensure we get a good view of the monolith */
COMPONENT Arm1 = Arm()
  AT (0,0,2) RELATIVE Source

COMPONENT Arm2 = Arm()
  AT (0,0,3.5) RELATIVE Source



COMPONENT AutoTOFL0 = Monitor_nD(xwidth=XW, yheight=YH, user1=T0, user2=L0, options="user1 limits=[0 5e-3] bins=51, user2 limits=[0.1 20] bins=41", restore_neutron=1)
 AT (0, 0, 0.001) RELATIVE BackTrace

COMPONENT AutoTOF0 = Monitor_nD(xwidth=XW, yheight=YH, user1=T0, options="user1 limits=[0 5e-3] bins=51", restore_neutron=1)
 AT (0, 0, 0.001) RELATIVE PREVIOUS

  COMPONENT AutoL0 = Monitor_nD(xwidth=XW, yheight=YH, user1=L0, options="user1 limits=[0.1 20] bins=41", restore_neutron=1)
 AT (0, 0, 0.001) RELATIVE PREVIOUS

  COMPONENT PSD0= Monitor_nD(filename="flat",xwidth=0.2,yheight=0.2,user1=SrcX,user2=SrcY,options="user1 limits=[-0.1 0.1] bins=90, user2 limits=[-0.1 0.1] bins=90,", restore_neutron=1)
  AT (0,0,0.001) RELATIVE PREVIOUS

  COMPONENT PSD1=Monitor_nD(filename="flatC",xwidth=0.2,yheight=0.2,user1=SrcX,user2=SrcY,options="user1 limits=[-0.1 0.1] bins=90, user2 limits=[-0.1 0.1] bins=90,", restore_neutron=1)
WHEN (Eneutron<EminTh) AT (0,0,0.001) RELATIVE PREVIOUS

  COMPONENT PSD2=Monitor_nD(filename="flatT",xwidth=0.16,yheight=0.16,user1=SrcX,user2=SrcY,options="user1 limits=[-0.1 0.1] bins=90, user2 limits=[-0.1 0.1] bins=90,", restore_neutron=1)
WHEN (Eneutron>=EminTh) AT (0,0,0.001) RELATIVE PREVIOUS


/* Measures the horizontal emmision coordinate of all neutrons - gives the "apparent width" of the moderators as seen from the beamline */
COMPONENT MonND1 = Monitor_nD(xwidth=XW, yheight=YH, user1=SrcX, username1="Horizontal position / [m]", options=options1, restore_neutron=1)
AT (0, 0, 1) RELATIVE Source

/* Measures the horizontal emmision coordinate of all neutrons - gives the "apparent width" of the moderators as seen from the beamline */
COMPONENT CWidth = Monitor_nD(xwidth=XW, yheight=YH, user1=SrcX, username1="Horizontal position / [m]", options=options1, restore_neutron=1)
  WHEN(Eneutron<=EmaxC && Eneutron>=EminC) AT (0, 0, 1) RELATIVE Source

/* Measures the horizontal emmision coordinate of all neutrons - gives the "apparent width" of the moderators as seen from the beamline */
COMPONENT TWidth = Monitor_nD(xwidth=XW, yheight=YH, user1=SrcX, username1="Horizontal position / [m]", options=options1, restore_neutron=1)
  WHEN(Eneutron<=EmaxTh && Eneutron>=EminTh) AT (0, 0, 1) RELATIVE Source

/* Measures the vertical emmision coordinate of cold neutrons */
COMPONENT MonND2 = Monitor_nD(xwidth=XW, yheight=YH, user1=SrcY, username1="Vertical position COLD / [m]", options=options4, restore_neutron=1)
  WHEN(IsCold) AT (0, 0, 1) RELATIVE Source

/* Measures the vertical emmision coordinate of thermal neutrons */
COMPONENT MonND2_2 = Monitor_nD(xwidth=XW, yheight=YH, user1=SrcY, username1="Vertical position THERMAL/ [m]", options=options4, restore_neutron=1)
  WHEN(!IsCold) AT (0, 0, 1) RELATIVE Source

/* 2D-plot of emmision coordinates for all neutrons */
COMPONENT MonND3 = Monitor_nD(xwidth=XW, yheight=YH, user1=SrcX, username1="Horizontal position / [m]", user2=SrcY,username2="Vertical position / [m]", options=options2, restore_neutron=1)
 AT (0, 0, 1) RELATIVE Source

/* 2D-plot of (x,z) emmision coordinates for all neutrons */
COMPONENT MonND4 = Monitor_nD(xwidth=XW, yheight=YH, user1=SrcX, username1="Emission position / [m]", user2=SrcZ, username2="Z-component of position / [m]", options="user1 bins=201 limits=[-0.3,0.3], user2 bins=201 limits=[-0.3,0.3]", restore_neutron=1)
 AT (0, 0, 1) RELATIVE Source
  EXTEND %{
    //printf("Time is %g\n",t);
%}



COMPONENT AutoTOFL = Monitor_nD(xwidth=0.1, yheight=0.1, options="tof limits=[0 15e-3] bins=51, lambda limits=[0.1 20] bins=41", restore_neutron=1)
 AT (0, 0, 2) RELATIVE Source



/* Measures brilliance of the "full" cold source */
COMPONENT BrillmonCOLD = Brilliance_monitor(
    nlam = 101, nt = 101, filename = "brillCOLD", t_0 = -1000,
    t_1 =4e4, lambda_0 = lambdamin, lambda_1 = lambdamax,
    Freq =14, toflambda=1 ,tofcuts=0, srcarea=(100*0.072*100*Yheight), restore_neutron=1,source_dist=2,xwidth=0.1,yheight=0.1)
WHEN(IsCold)  AT (0, 0, 2) RELATIVE Source

/* Measures "collimated" brilliance of the cold source over fixed 6 cm wide area x central part vertically. */
/* Used for calibration of performance wrt. MCNP BF1 output, see http://ess_butterfly.mcstas.org */
COMPONENT BrillmonCOLD_COLL = Brilliance_monitor(
    nlam = 101, nt = 101, filename = "brillCOLD_COLL", t_0 = 0,
    t_1 = 4e4, lambda_0 = lambdamin, lambda_1 = lambdamax,
    Freq =14, toflambda=1,tofcuts=0, srcarea=(100*0.06*100*2*Yheight/2.5), restore_neutron=1,source_dist=2,xwidth=0.1,yheight=0.1)
  WHEN(fabs(SrcY)<Yheight/2.5 && SrcX < (0.071+delta) && SrcX > (0.011+delta))  AT (0, 0, 2) RELATIVE Source

/* Measures brilliance of the "full" thermal source */
COMPONENT BrillmonTHRM = Brilliance_monitor(
    nlam = 101, nt = 101, filename = "brillTHRM", t_0 = -1000,
    t_1 =4e4, lambda_0 = lambdamin, lambda_1 = lambdamax,
    Freq =14, toflambda=1,tofcuts=0, srcarea=(100*0.108*100*Yheight), restore_neutron=1,source_dist=2,xwidth=0.1,yheight=0.1)
  WHEN (!IsCold) AT (0, 0, 2) RELATIVE Source

/* Measures "collimated" brilliance of the thermal source over fixed 6 cm wide area (or smaller at beamlines no. 1,2) x central part vertically. */
/* Used for calibration of performance wrt. MCNP BF1 output, see http://ess_butterfly.mcstas.org */
COMPONENT BrillmonTHRM_COLL = Brilliance_monitor(
    nlam = 101, nt = 101, filename = "brillTHRM_COLL", t_0 = -1000,
    t_1 =4e4, lambda_0 = lambdamin, lambda_1 = lambdamax,
      Freq =14, toflambda=1,tofcuts=0, srcarea=(100*0.06*100*2*Yheight/2.5), restore_neutron=1,source_dist=2,xwidth=0.1,yheight=0.1)
  WHEN (fabs(SrcY)<Yheight/2.5 && -SrcX>(TCollmin+delta) && -SrcX<(TCollmax+delta))  AT (0, 0, 2) RELATIVE Source

  COMPONENT PSD0x= Monitor_nD(filename="flat_x",xwidth=0.1,yheight=0.1,user1=SrcX,user2=SrcY,options="user1 limits=[-0.08 0.08] bins=90, user2 limits=[-0.08 0.08] bins=90,", restore_neutron=1)
  AT (0,0,0.001) RELATIVE PREVIOUS

  COMPONENT PSD1x=Monitor_nD(filename="flatC_x",xwidth=0.1,yheight=0.1,user1=SrcX,user2=SrcY,options="user1 limits=[-0.08 0.08] bins=90, user2 limits=[-0.08 0.08] bins=90,", restore_neutron=1)
WHEN (Eneutron<EminTh) AT (0,0,0.001) RELATIVE PREVIOUS

  COMPONENT PSD2x=Monitor_nD(filename="flatT_x",xwidth=0.1,yheight=0.1,user1=SrcX,user2=SrcY,options="user1 limits=[-0.08 0.08] bins=90, user2 limits=[-0.08 0.08] bins=90,", restore_neutron=1)
WHEN (Eneutron>=EminTh) AT (0,0,0.001) RELATIVE PREVIOUS

COMPONENT GuideR = Guide_curved(
    w1 = 0.1, h1 = 0.1, l = 50, curvature=3000)
  AT (0,0,0.001) RELATIVE PREVIOUS

  COMPONENT RArm=Arm()
  AT (0.416657,0,50.002) RELATIVE GuideR
  ROTATED (0, calcAlpha(50,3000),0) RELATIVE GuideR
  
  COMPONENT RArm2=Arm()
  AT (calcX(50, 3000),0,calcZ(50, 3000)) RELATIVE GuideR
  ROTATED (0, calcAlpha(50,3000),0) RELATIVE GuideR
  
COMPONENT Monitor2_xy1 = Monitor_nD(
    options = "x limits=[-0.06 0.06] bins=51, y limits=[-0.06 0.06] bins=51,", xwidth = 0.12, yheight = 0.12)
 AT (0, 0, 0.05) RELATIVE RArm

/* /\* Uncomment these helper-arms to view "full" monolith *\/ */

COMPONENT DummyArm1 = Arm()
  AT (6,0,6) ABSOLUTE

COMPONENT DummyArm2 = Arm()
  AT (-6,0,6) ABSOLUTE

COMPONENT DummyArm3 = Arm()
  AT (-6,0,-6) ABSOLUTE

COMPONENT DummyArm4 = Arm()
  AT (6,0,-6) ABSOLUTE

COMPONENT DummyArm5 = Arm()
  AT (6,0,6) ABSOLUTE


FINALLY
%{
%}

END