/*******************************************************************************
* 
* Instrument: ESS_BEER_MCPL
* 
* %Identification
* Written by: Jan Saroun, saroun@ujf.cas.cz
* Date: June 2017
* Origin: NPI Rez
* %INSTRUMENT_SITE: ESS
* Version: 1.1, 27/11/2018 
* ToF Diffractometer BEER@ESS with MCPL input at the sample slit.
*
* %Description 
*
* Secondary part of the ToF diffractometer BEER@ESS. It takes MCPL file at the input as the source 
* of primary beam in front of the sample.
* Results:
* 1) xy, divergence, lambda and ToF_lambda plots of the primary beam (monitors the MCPL content)
* 2) projections of the sampling volume defined by the primary slit and secondary radial collimator  
* 3) Intensity vs. dhkl plot, using NPI_tof_dhkl_detector.comp
* 4) 2D plot (ToF,2theta),  using NPI_tof_theta_monitor
*  
* Modulation mode: 
* If modul=1, a modulation mode of BEER is assumed: the primary beam is modulated by a chopper placed close to the source.
* (set the lc distance appropriately). The NPI_tof_dhkl_detector.comp then performs data reduction in event mode, using 
* the given modulation parameters: lam0, mod_frq, mod_twidth and a table with primary dhkl estimates (dhkl.dat should be 
* in the current directory). As a result, the diffractogram with recombined peaks is produced, together with a 2D map 
* of estimated valid, empty and overlap regions on the ToF-2theta diagram. 
*
* If module=0, modulation parameters are ignored and the NPI_tof_dhkl_detector.comp performs a standard event based 
* data reduction producing a diffractogram on the basis of given nominal values of lam0 and distances.
* 
* (See NPI_tof_dhkl_detector.comp for more details)
* 
* NOTE: some instrument parameters are hard coded and have to be edited in the INITIALIZE and RUN sections.
*
* To use, please copy the relevant mcpl files from your $MCSTAS/data folder to the curring working folder
*
* Example 1: 
*  Simulation in medium resolution mode with pulse shaping choppers, wavelength range centred at 2 AA - use:
*  ESS_BEER_MCPL input=BEER_MR.mcpl repetition=50 pwdfile=duplex.laz lc=6.65 lam0=2 dlam=1.8 omega=45 chi=90 colw=1 modul=0 mod_frq=2240 mod_twidth=0.0029 mod_shift=0 only_event=-1 pinc=0.1 ptra=0 strain=0 ustrain=0 
*  Detector: psdtof_I=276.842
*
* Example 2: 
*  Simulation in the modulation mode, using the modulation chopper MCB with 8 slits (4 deg wide) rotating at 280 Hz. Wavelength range centred at 2 AA - use:
*  ESS_BEER_MCPL input=BEER_MCB.mcpl repetition=50 pwdfile=duplex.laz lc=9.35 lam0=2 dlam=1.8 omega=45 chi=90 colw=1 modul=1 mod_frq=2240 mod_twidth=0.0029 mod_shift=0 only_event=-1 pinc=0.1 ptra=0 strain=0 ustrain=0 
*  Detector: psdtof_I=98.3885
*
* %Example: input=BEER_MR.mcpl repetition=50 lc=6.65 modul=0 mod_twidth=0.0029 Detector: psdtof_I=138.421
* Example: input=BEER_MCB.mcpl repetition=50 lc=9.35 modul=1 mod_twidth=0.0029 Detector: psdtof_I=49.1942
*
* %Parameters
*
* input: [str]      Input MCPL file	
* repetition: [1]    Number of loops through the MCPL file	
* pwdfile: [str]    Input sample file for PowderN.comp
* lc: [m]           distance of the pulse definition chopper 
* lam0: [AA]        nominal wavelength (centre of the frame, determines the chopper phase) 
* dlam: [AA]        wavelength band width (only for filtering MCPL input and plot ranges)
* omega: [deg]      sample orientation (y-axis)
* chi: [deg]        sample orientation (z-axis)
* colw: [mm]        collimator width (0.5, 1, 2, 3 or 4)
* modul: [0|1]      modulation mode switch
* mod_frq: [Hz]     modulation frequency (chopper frequency x number of slits)
* mod_twidth: [s]   modulation frame width (should be ~> ESS pulse width)
* mod_shift: []     assumed line shift introduced to NPI_tof_dhkl_detector (modulation mode only)
* only_event: [1]   if > -1, filters out events with line_info.itype<>only_event after PowderN sample
* ptra: [ ]         p_transmit value passed to PowderN.comp
* pinc: [ ]         pinc value passed to PowderN.comp
* strain: [ppm]     Macro-strain (peak shift)
* ustrain: [ppm]    Micro-strain (peak broadening)
*
* %End
*******************************************************************************/
DEFINE INSTRUMENT ESS_BEER_MCPL(string input="BEER_MR.mcpl", int repetition=50, string pwdfile="duplex.laz", 
lc=6.65, lam0=2.0, dlam=1.8, omega=45, chi=90, colw=1, modul=0, mod_frq=2240, mod_twidth=0.0029, mod_shift=0,
only_event=-1, pinc=0.1, ptra=0.0, strain=0, ustrain=0)
DEPENDENCY "-DMCPLPATH=GETPATH(data)"

DECLARE
%{

  char optSGV1[256];
  char optSGV2[256]; 
  double col_ang,col_angv,col_h1,col_h2,col_rad,col_len,col_d;
  int col_ns;
  double mask_w,mask_dist,mask_h; 
  double det_th1,det_Linst,det_th2,det_t0,det_d1,det_d2,det_Lc,det_rad;
  double det_mod_dt, det_mod_twidth; 
  int det_modulation;
  char dhklTable[256];    
  double Lmin, Lmax, Emin, Emax;
 
  // perform r = m*v
  void MXV(Rotation m, double v[3], double r[3]) 
  {
        int i,j,k;
		for (i=0;i<3;i++) {   
			r[i]=0;  
			for (j=0;j<3;j++) { 
				r[i] += m[i][j]*v[j];
			}
		}
  };
  Rotation smi;
  int filt;
  double tof;
  #pragma acc declare create(smi,filt,tof)
%}

USERVARS
%{
  double scX;
  double scY;
  double scZ;
  double scP;
%}

INITIALIZE
%{
double hm = 2*PI*K2V; // h/m_n
printf("Using the input file: %s\n", input);
printf("Using the sample file: %s\n", pwdfile);

/*
----------------------------
Set instrument parameters
----------------------------
*/ 
double Linst=158; // source - sample distance


/* Collimator (ENGIN-X system) 
--------------------------------*/
const double cold[5]={0.100, 0.160, 0.310, 0.410, 0.490};  
col_ang=30; // horizontal angle
col_angv=30; // vertical angle
col_d=0.0001;
col_ns=160;
col_len=0.35;
// select distance according to width, allowed values 0.5, 1, 2, 3, 4
if (colw==0.5) {
	col_rad=cold[0];
} else if (colw==1) { 
	col_rad=cold[1];
} else if (colw==2) {
	col_rad=cold[2];
} else if (colw==3) {
	col_rad=cold[3];
} else if (colw==4) {
	col_rad=cold[4];
} else { 
	printf("WARNING: invalid collmator width (%g), selecting 2 mm.\n", colw);
	col_rad=cold[2];
}

// calculate height (add max. sample height 2 cm)
col_h1=0.02+2*col_rad*tan(col_angv*DEG2RAD/2);
col_h2=0.02+2*(col_rad+col_len)*tan(col_angv*DEG2RAD/2);
printf("INFO: Collimator radius=%g, height=%g -> %g\n",col_rad,col_h1,col_h2);

// mask before collimator
mask_w=2*col_rad*fabs(tan(col_ang*0.5*DEG2RAD));
mask_h=col_h1; 
mask_dist=col_rad*cos(col_ang*0.5*DEG2RAD)-0.01;

/* Detector parameters 
--------------------------------*/
// angular range
det_th1=75;
det_th2=105;
// detection radius
det_rad=2;
// distance from the source to the detector
det_Linst=Linst+det_rad;
// distance moderator - pulse chopper
det_Lc=lc;
// chopper delay
det_t0=lc/(hm/lam0);
// dhkl range
det_d1=0.8;
det_d2=2.4;

/* Modulation parameters 
--------------------------------*/
det_modulation=modul;
det_mod_dt=1/mod_frq;
det_mod_twidth=mod_twidth;

/* Sample settings 
--------------------------------*/
// d0 table for modulation mode
sprintf(dhklTable,"dhkl.tab");
// set sample orientation
Rotation sm;
rot_set_rotation(sm, 0.0,omega*DEG2RAD,chi*DEG2RAD); 
rot_transpose(sm, smi);
tof=det_Linst/hm*lam0*1e6; // tof in [us]
filt=only_event;

#pragma acc update device(smi,tof,filt)

/* 2D plot of the sampling gauge volume, in mm*/
sprintf(optSGV2,"user1 bins=201 limits=[-5,5], user2 bins=201 limits=[-5,5]");
/* 2D plot of the sampling gauge volume, in mm*/
sprintf(optSGV1,"user1 bins=201 limits=[-5,5]");

/* Calculate MCPL_input parameters 
-----------------------------------*/
double L2E;
Lmin = lam0 - 0.5*dlam;
Lmax = lam0 + 0.5*dlam;
L2E = pow(2*PI*K2V,2)*VS2E;
Emin = L2E/pow(Lmax,2);
Emax = L2E/pow(Lmin,2);

/* Messages 
--------------------------------*/
printf("Chopper delay = %g, mod_period = %g [ms]\n", det_t0*1000, det_mod_dt*1000);
printf("Mean ToF [ms] = %g\n", tof/1000);
printf("Selected energy range: %g to %g meV\n", Emin, Emax);
%}

TRACE
 
COMPONENT Origin = Progress_bar()
  AT (0, 0, 0) ABSOLUTE
  
/* Read neutrons from an mcpl file. */
COMPONENT src = MCPL_input(filename=input, polarisationuse=0, verbose=0, Emin=Emin, Emax=Emax, repeat_count=repetition, v_smear=0, pos_smear=0.0001, dir_smear=0.001)
AT( 0,0,0) RELATIVE PREVIOUS

COMPONENT xymon = PSD_monitor(
    nx=100, 
    ny=100,  
    filename="xymon.dat", 
    xwidth=0.01,  
    yheight=0.01,  
    restore_neutron=1) 
AT (0, 0, 0.001) RELATIVE src 

COMPONENT lmon = L_monitor( 
    nL=60, 
    filename="lmon.dat", 
    xwidth=0.05, 
    yheight=0.05, 
    Lmin=Lmin,  
    Lmax=Lmax, 
    restore_neutron=1)
AT (0, 0, 0.002) RELATIVE src   

COMPONENT hdiv_mon = Div1D_monitor(
    ndiv = 100, filename = "hdiv.dat", xwidth = 0.05,
    yheight = 0.05, maxdiv = 0.5, restore_neutron = 1)
  AT (0, 0, 0.003) RELATIVE src


COMPONENT toflam = TOFLambda_monitor(
    nL = 400, nt = 400, tmin = 70e3, tmax = 90e3,
    filename = "toflam.dat", xwidth = 0.05, yheight = 0.05,
    Lmin = 1.9, Lmax = 2.1, restore_neutron = 1)
  AT (0, 0, 0.003) RELATIVE src


/* Place the sample axis to the correct distance after the MCPL_input.
Depends on the configuration of the primary beam simulation. Here we assume 40 mm.*/  
COMPONENT Sample_axis = Arm() 
AT (0, 0, 0.04) RELATIVE src
ROTATED (0, 0, 0) RELATIVE src

/* Detector arm - defines scattering angle at the detector centre. */ 
COMPONENT Detector_arm = Arm() 
AT (0, 0, 0) RELATIVE Sample_axis
ROTATED (0, 90, 0) RELATIVE Sample_axis

COMPONENT sample = PowderN(
    reflections=pwdfile, 
    yheight=0.05,
    radius=0.0035, 
	d_omega=col_ang,
    d_phi=col_angv, 
	tth_sign=1,
//	sigma_abs=2.56,
//	sigma_inc=0.4,
 //   density=7.87,
	p_transmit=ptra, 
	p_inc=pinc,
	Strain=strain*1e-6,
	delta_d_d=ustrain*1e-6,
	focus_flip=0, 
	target_index=1) 
AT (0, 0, 0) RELATIVE Sample_axis
ROTATED (0, omega, chi) RELATIVE Sample_axis
EXTEND %{
  /*if ((filt>-1) && (itype!=filt)) {
 ABSORB; 
 }*/
// gauge coord. in mm
double r0[3] = {x,y,z};
double r[3]; 
MXV(smi, r0, r);
scX = r[0]*1000; 
scY = r[1]*1000;
scZ =  r[2]*1000;
scP = p; 
%}

/* Sample focuses at this component. */
COMPONENT col_mask = Slit(
    xwidth = mask_w, yheight = mask_h)
  AT (0, 0, mask_dist) RELATIVE Detector_arm

COMPONENT rad2 = Exact_radial_coll(
    theta_min=-col_ang*0.5,  
    theta_max=col_ang*0.5,  
    nslit=col_ns, 
    radius=col_rad, 
    length=col_len, 
    h_in=col_h1, 
    h_out=col_h2,  
	d=col_d, 
    verbose=1)
AT (0, 0, 0) RELATIVE Detector_arm

/* Gauge volume, XZ plot */
COMPONENT MonNDXZ = Monitor_nD(xwidth=2, yheight=2, user1="scX",
  username1="X, mm", user2="scZ", username2="Z, mm", options=optSGV2,
  filename="SGV_xz.dat", 
  restore_neutron=1) 
  AT (0, 0, col_rad+col_len+0.01) RELATIVE Detector_arm
 
 /* Gauge volume, X plot */
COMPONENT MonNDX = Monitor_nD(xwidth=2, yheight=2, user1="scX",
  username1="X, mm", options=optSGV1,
  filename="SGV_x.dat",   
  restore_neutron=1) 
  AT (0, 0, col_rad+col_len+0.01) RELATIVE Detector_arm
  
   /* Gauge volume, Y plot */
COMPONENT MonNDY = Monitor_nD(xwidth=2, yheight=2, user1="scY",
  username1="Y, mm", options=optSGV1,
  filename="SGV_y.dat", 
  restore_neutron=1) 
  AT (0, 0, col_rad+col_len+0.01) RELATIVE Detector_arm

   /* Gauge volume, Z plot */
COMPONENT MonNDZ = Monitor_nD(xwidth=2, yheight=2, user1="scZ",
  username1="Z, mm", options=optSGV1,
  filename="SGV_z.dat", 
  restore_neutron=1)
  AT (0, 0, col_rad+col_len+0.01) RELATIVE Detector_arm
  
  /* Detector with event mode data reduction.
  Generates 1D diffractogram and a 2D overlap map (in modulation mode).
  Should be centered at the sample axis and aligned with the primary beam.
  */ 
  COMPONENT dhklmon = NPI_tof_dhkl_detector(  
	nd=3000, 
	filename="dhkl.dat",  
	yheight=1.0, 
	zdepth=0.01, 
	radius=2, 
	amin=det_th1,  
	amax=det_th2, 
	d_min=det_d1,  
	d_max=det_d2, 
	time0=det_t0,
	Linst=det_Linst, 
	Lc = det_Lc,
	res_x=0.002,    
	res_y=0.005, 
	res_t=1e-6,    
	mu=1.0,  
	modulation=det_modulation,    
	mod_dt=det_mod_dt,    
	mod_twidth=det_mod_twidth ,    
	mod_shift=mod_shift,
	mod_d0_table=dhklTable,  
	restore_neutron=1)  
  AT (0, 0, 0) RELATIVE Sample_axis  
 
 /* ToF vs. 2theta map */
  COMPONENT psdtof = NPI_tof_theta_monitor( 
    nt = 800, na = 600, filename = "tof_theta.dat",
    radius = 2, yheight = 1, tmin = tof-40e3, tmax = tof+40e3, 
    amin = det_th1, amax = det_th2, restore_neutron = 1,verbose=0)       
  AT (0, 0, 0) RELATIVE Sample_axis 
  
 /* ToF vs. 2theta map, detail */ 
  COMPONENT psdtofDetail= NPI_tof_theta_monitor(
    nt = 400, na = 400, filename = "tof_theta_detail.dat", 
    radius = 2, yheight = 1, tmin = 100e3, tmax = 110e3, 
    amin = 75, amax = 80, restore_neutron = 1,verbose=0)       
  AT (0, 0, 0) RELATIVE Sample_axis     
 

END