BL18_pxtal_BOTH.instr

/***************************************************************************
*
* Written by: A. Cereser and E. Knudsen
* Based on code by K. Oikawa
* September 2013 - February 2016
* Written at DTU and J-PARC
*
* Instrument: BL18 SENJU at J-PARC
*
* mcstas ./BL18_pxtal_event.instr
*
* %P
* L_min (AA) lower edge of wavelength band to be emitted from the Source
* L_max (AA) upper edge of wavelength band to be emitted from the Source
* NX ( ) Number of voxels along X
* NY ( ) Number of voxels along Y
* NZ ( ) Number of voxels along Z
* XW (m) Width of the polycrystal hull.
* YH (m) Height of the polycrystal hull
* ZD (m) Depth og the polycrystal hull
* map_fn ( ) name of file containing the orientation map
* orts_fn ( ) name of file containing orientations
* lau_fn ( ) name of file ocntaining the reflection list
* AX ( ) X-component of unit cell vector a
* ...
* CZ ( ) Z-component of unit cell vector c
* samplename ( ) name of the sample to produce SENJU style filenames
* runnumber ( ) runnumber to produce SENJU style filenames
*
* These are some lines from the original J-parc simulation
*	E_min=81.8049/lmax/lmax;
*	E_max=81.8049/lmin/lmin;
*	tof_min=dd*505.555*35.6*sin((TTT)/360*PI)-delt;
*	tof_max=dd*505.555*35.6*sin((TTT)/360*PI)+delt;
*	E_min=81.8049*(252.778*35.6/tof_max)*(252.778*35.6/tof_max);
*	E_max=81.8049*(252.778*35.6/tof_min)*(252.778*35.6/tof_min);
*
**************************************************************************/
DEFINE INSTRUMENT BL18_SENJU(L_min=0.4, L_max=4, NX=1, NY=1, NZ=1, XW=0.01,
    YH=0.01, ZD=0.01, string map_fn="default.map_101x101_rectified",
    string orts_fn="default.orts",string lau_fn="Fe.lau",
    AX = 2.886,AY = 0, AZ = 0.0, BX = 0.0, BY = 2.886, BZ = 0.0, CX = 0.0,
    CY = 0.0, CZ = 2.886,omega=0, string samplename="C0127", int runnumber=0,
    int side=0
)

DECLARE
%{
  double lmin, lmax, E_min, E_max, tof_min, tof_max;
  double TTS, UML;
  double dd, delt;
  double shtwin, shthin, shtwout, shthout, nchann;
  double ddd, mos, ndd, detx;
  double alpha,wsht, Pi, nD, m5;
  double TG, TG1, TG2, TG3, TG4;
  double mx_Vd, my_Vd;
  double dsp, lam, lw, TT, TTT, TTD;
  double tt0,tt1;
  double ux,uy,uz;
  double kux,kuy,kuz;
  double anx,any,anz;
  int xtal_scatter;
  char crystal_file[]="Al.lau";
  //the crystal files are in /usr/local/lib/mcstas-2.0/data

  void compute_global_vector(char *s,double kx,double ky, double kz,double *gkx,
  double *gky, double *gkz, Rotation R){
    Rotation TT;
    rot_transpose(R, TT);
    /* now make the coordinate system change */
    coords_get(rot_apply(TT,coords_set(kx,ky,kz)),gkx,gky,gkz);
  }

  void compute_global_coords(char *s,double x,double y,double z,double *gx,
  double *gy,double *gz,Coords A,Rotation R){
    Rotation TT;
    rot_transpose(R, TT);
    /* now make the coordinate system change */
    coords_get(coords_add(rot_apply(TT,coords_set(x,y,z)),A),gx,gy,gz);
  }
  int sx_scatter;
  int anton;
  int oidx;

  long long nid;

  /*The lines below have been added to deal with the polycrystalline sample*/
  %include "read_table-lib.h"
  int nx;                 /*Number of voronoi regions in each direction*/
  int ny;                 /*The number has to be the same in pixellation.m */
  int nz;
  int ix,iy,iz,crystal_miss;
  int xtals;

  double posx,posy,posz;

  t_Table *map_table,*orts_table;

  Coords ao,bo,co;

  /*heres where we store the orientation dependent hkl_info pointers*/
  struct hkl_info_struct *poly_hkl_info;
  int initd;

  int sx_background;
  /*flag to indicate whether or not this is a void in the crystal*/

  char det_fn[41][256];
  char *det_fn_001=det_fn[1];
  char *det_fn_002=det_fn[2];
  char *det_fn_003=det_fn[3];
  char *det_fn_004=det_fn[4];
  char *det_fn_005=det_fn[5];
  char *det_fn_006=det_fn[6];
  char *det_fn_007=det_fn[7];
  char *det_fn_008=det_fn[8];
  char *det_fn_009=det_fn[9];
  char *det_fn_010=det_fn[10];
  char *det_fn_011=det_fn[11];
  char *det_fn_012=det_fn[12];
  char *det_fn_013=det_fn[13];
  char *det_fn_014=det_fn[14];
  char *det_fn_015=det_fn[15];
  char *det_fn_016=det_fn[16];
  char *det_fn_017=det_fn[17];
  char *det_fn_018=det_fn[18];
  char *det_fn_019=det_fn[19];
  char *det_fn_020=det_fn[20];
  char *det_fn_021=det_fn[21];
  char *det_fn_022=det_fn[22];
  char *det_fn_023=det_fn[23];
  char *det_fn_024=det_fn[24];
  char *det_fn_025=det_fn[25];
  char *det_fn_026=det_fn[26];
  char *det_fn_027=det_fn[27];
  char *det_fn_028=det_fn[28];
  char *det_fn_029=det_fn[29];
  char *det_fn_030=det_fn[30];
  char *det_fn_031=det_fn[31];
  char *det_fn_032=det_fn[32];
  char *det_fn_033=det_fn[33];
  char *det_fn_034=det_fn[34];
  char *det_fn_035=det_fn[35];
  char *det_fn_036=det_fn[36];
  char *det_fn_037=det_fn[37];
  char *det_fn_038=det_fn[38];
  char *det_fn_039=det_fn[39];
  char *det_fn_040=det_fn[40];

%}

INITIALIZE
%{
  Pi=3.14159;
  alpha=3.0;
  wsht=3e-3;

  /*Section below added to work with the polycrystalline sample*/

  /*These should come from analyzing the map file instead*/
  nx=NX;
  ny=NY;
  nz=NZ;


  /*Force components initialization to avoid inifite recursion*/
  if (initd!=1){
    initd=1;
    mcinit();

    int status;
    map_table=malloc(sizeof(t_Table));
    orts_table=malloc(sizeof(t_Table));

    if ( (status=Table_Read(orts_table,orts_fn,0))==-1){
      fprintf(stderr,"Error (%s): Could not parse file \"%s\"\n",
      NAME_CURRENT_COMP,orts_fn);
      exit(-1);
    }
    if ( (status=Table_Read(map_table,map_fn,0))==-1){
      fprintf(stderr,"Error (%s): Could not parse file \"%s\"\n",
      NAME_CURRENT_COMP,map_fn);
      exit(-1);
    }

    ao=coords_set(MC_GETPAR2(SX,ax),MC_GETPAR2(SX,ay),MC_GETPAR2(SX,az));
    bo=coords_set(MC_GETPAR2(SX,bx),MC_GETPAR2(SX,by),MC_GETPAR2(SX,bz));
    co=coords_set(MC_GETPAR2(SX,cx),MC_GETPAR2(SX,cy),MC_GETPAR2(SX,cz));

    /*allocate memory for all orientation hkl_infos*/
    if( (poly_hkl_info=malloc(orts_table->rows*sizeof(struct hkl_info_struct)))==NULL ){
        fprintf(stderr,"Error (%s): Memory allocation error for orts hkl_infos,
        aborting.\n",NAME_CURRENT_COMP);
        exit(-1);
    }

    /*get important parameters from actual component*/
    char *laufn=MC_GETPAR2(SX,reflections);
    double SC_mosaic=MC_GETPAR2(SX,mosaic);
    double SC_mosaic_a=MC_GETPAR2(SX,mosaic_a);
    double SC_mosaic_b=MC_GETPAR2(SX,mosaic_b);
    double SC_mosaic_c=MC_GETPAR2(SX,mosaic_c);
    double *SC_mosaic_ABin=NULL;

    /*loop over all orientations and set up info-lists for each one.*/
    int i;
    for (i=0;i<orts_table->rows;i++){
      poly_hkl_info[i]=MC_GETPAR2(SX,hkl_info);
      /*by def. we are dealing with crystal prms. in real space*/
      poly_hkl_info[i].recip=0;

      /*apply rotation to crystal vectors*/
      Coords a,b,c;
      Rotation U;
      memcpy(*U,&(orts_table->data[i*9]),sizeof(U[0][0])*9);
      a=rot_apply(U,ao);
      b=rot_apply(U,bo);
      c=rot_apply(U,co);
      /*set the crystal parameters to the rotated ones*/
      coords_get(a,&(poly_hkl_info[i].m_ax), &(poly_hkl_info[i].m_ay),
      &(poly_hkl_info[i].m_az));
      coords_get(b,&(poly_hkl_info[i].m_bx), &(poly_hkl_info[i].m_by),
      &(poly_hkl_info[i].m_bz));
      coords_get(c,&(poly_hkl_info[i].m_cx), &(poly_hkl_info[i].m_cy),
      &(poly_hkl_info[i].m_cz));

      /*call read_hkl_list function*/
      if (!read_hkl_data(laufn, &(poly_hkl_info[i]), SC_mosaic, SC_mosaic_a,
      SC_mosaic_b, SC_mosaic_c, SC_mosaic_ABin))
        exit(-1);


    }
    /*This should only be done once*/
    /*adjust the tunnelling probability to avoid statistical extinction*/
    int nmax;
    if (NX>NY) nmax=NX;
    else nmax=NY;
    if(nmax<NZ) nmax=NZ;

    MC_GETPAR2(SX,p_transmit)=exp(log(MC_GETPAR2(SX,p_transmit))/(double)nmax);
    /*end of polycrystal*/

    return;
    /*do other initialization stuff*/
  }else{
    /*this block gets done before component initialization*/
    /*set up detector filenames*/

    int j;
    for (j=1;j<=40;j++){
      sprintf(det_fn[j],"%s_%06d_%03d",samplename,runnumber,j);
      if(j==1){
        printf("detector filenames:\n%s\n .\n .\n .\n",det_fn[j]);
      }else if(j==40){
        printf("%s\n",det_fn[j]);
      }
    }

    /*position SENJU far field detectors according to calibration*/

  }
%}

TRACE

COMPONENT Origin = Progress_bar()
  AT (0,0,0) ABSOLUTE
EXTEND
%{
 anton=0;
 xtals=0;
 nid=mcget_run_num();
 //printf("%ld\n",nid);

%}

COMPONENT a1 = Arm()
  AT (0,0,0) ABSOLUTE

COMPONENT mod = SNS_source(
  filename="source_BL18.txt",
  xwidth=0.1, yheight=0.1,
  dist = 15.19,
  focus_xw = 0.041, focus_yh = 0.041,
  Emax=pow((2*M_PI/L_min)*K2V,2.0)*VS2E, Emin=pow((2*M_PI/L_max)*K2V,2.0)*VS2E
) AT (0,0,0) RELATIVE a1

/*Uncomment the section below to get a pencil beam

COMPONENT mod=Modr_dcHPb_flux(
	xws=0.1, yhs=0.1, angle=0.0,
	Emin=E_min, Emax=E_max,
  dist = 15.19,
  xw = 0.041, yh = 0.041
) AT (0,0,0)  RELATIVE a1
*/

COMPONENT slt_in = Slit(
	xwidth = 0.070, yheight = 0.070)
	AT (0, 0, 2.271) RELATIVE a1

COMPONENT slt_out = Slit(
	xwidth = 0.044, yheight = 0.044)
	AT (0, 0, 15.140) RELATIVE a1

COMPONENT guide1 = Guide(w1=0.03951,
			 h1=0.03951,
			 w2=0.04146,
			 h2=0.04146,
			 l=0.8-0.0000001, R0=0.99, Qc=0.0217, alpha=4,m=4.0,W=0.0015)
			 AT (0, 0, 15.2) RELATIVE a1 ROTATED (0,0,0) RELATIVE a1

COMPONENT guide2 = Guide(w1=0.04146,
			 h1=0.04146,
			 w2=0.04347,
			 h2=0.04347,
			 l=1.0-0.0000001, R0=0.99, Qc=0.0217, alpha=3,m=3.2,W=0.001)
			 AT (0, 0, 16.0) RELATIVE a1 ROTATED (0,0,0) RELATIVE a1

COMPONENT guide3 = Guide(w1=0.04347,
			 h1=0.04347,
			 w2=0.04507,
			 h2=0.04507,
			 l=1.0-0.0000001, R0=0.99, Qc=0.0217,  alpha=3,m=3.2,W=0.001)
			 AT (0, 0, 17.0) RELATIVE a1 ROTATED (0,0,0) RELATIVE a1

COMPONENT guide4 = Guide(w1=0.04507,
			 h1=0.04507,
			 w2=0.04630,
			 h2=0.04630,
			 l=1.0-0.0000001, R0=0.99, Qc=0.0217, alpha=3,m=3.2,W=0.001)
			 AT (0, 0, 18.0) RELATIVE a1 ROTATED (0,0,0) RELATIVE a1

COMPONENT guide5 = Guide(w1=0.04630,
			 h1=0.04630,
			 w2=0.04718,
			 h2=0.04718,
			 l=1.0-0.0000001, R0=0.99, Qc=0.0217, alpha=3,m=3.2,W=0.001)
			 AT (0, 0, 19.0) RELATIVE a1 ROTATED (0,0,0) RELATIVE a1

COMPONENT guide6 = Guide(w1=0.04718,
			 h1=0.04718,
			 w2=0.04774,
			 h2=0.04774,
			 l=1.0-0.0000001, R0=0.99, Qc=0.0217,  alpha=3,m=3.2,W=0.001)
			 AT (0, 0, 20.0) RELATIVE a1 ROTATED (0,0,0) RELATIVE a1

COMPONENT guide7 = Guide(w1=0.04774,
			 h1=0.04774,
			 w2=0.04799,
			 h2=0.04799,
			 l=1.0-0.0000001, R0=0.99, Qc=0.0217,  alpha=3,m=3.2,W=0.001)
			 AT (0, 0, 21.0) RELATIVE a1 ROTATED (0,0,0) RELATIVE a1

COMPONENT guide8 = Guide(w1=0.04799,
			 h1=0.04799,
			 w2=0.04799,
			 h2=0.04799,
			 l=0.6-0.0000001, R0=0.99, Qc=0.0217,  alpha=3,m=3.2,W=0.001)
			 AT (0, 0, 22.0) RELATIVE a1 ROTATED (0,0,0) RELATIVE a1

COMPONENT guide9 = Guide(w1=0.04799,
			 h1=0.04799,
			 w2=0.04774,
			 h2=0.04774,
			 l=1.0-0.0000001, R0=0.99, Qc=0.0217,  alpha=3,m=3.2,W=0.001)
			 AT (0, 0, 22.6) RELATIVE a1 ROTATED (0,0,0) RELATIVE a1

COMPONENT guide10 = Guide(w1=0.04774,
			 h1=0.04774,
			 w2=0.04718,
			 h2=0.04718,
			 l=1.0-0.0000001, R0=0.99, Qc=0.0217,  alpha=3,m=3.2,W=0.001)
			 AT (0, 0, 23.6) RELATIVE a1 ROTATED (0,0,0) RELATIVE a1

COMPONENT guide11 = Guide(w1=0.04718,
			 h1=0.04718,
			 w2=0.04630,
			 h2=0.04630,
			 l=1.0-0.0000001, R0=0.99, Qc=0.0217,  alpha=3,m=3.2,W=0.001)
			 AT (0, 0, 24.6) RELATIVE a1 ROTATED (0,0,0) RELATIVE a1

COMPONENT guide12 = Guide(w1=0.04630,
			 h1=0.04630,
			 w2=0.04507,
			 h2=0.04507,
			 l=1.0-0.0000001, R0=0.99, Qc=0.0217,  alpha=3,m=3.2,W=0.001)
			 AT (0, 0, 25.6) RELATIVE a1 ROTATED (0,0,0) RELATIVE a1

COMPONENT guide13 = Guide(w1=0.04507,
			 h1=0.04507,
			 w2=0.04347,
			 h2=0.04347,
			 l=1.0-0.0000001, R0=0.99, Qc=0.0217,  alpha=3,m=3.2,W=0.001)
			 AT (0, 0, 26.6) RELATIVE a1 ROTATED (0,0,0) RELATIVE a1

COMPONENT guide14 = Guide(w1=0.04347,
			 h1=0.04347,
			 w2=0.04146,
			 h2=0.04146,
			 l=1.0-0.0000001, R0=0.99, Qc=0.0217,  alpha=3,m=3.2,W=0.001)
			 AT (0, 0, 27.6) RELATIVE a1 ROTATED (0,0,0) RELATIVE a1

COMPONENT guide15 = Guide(w1=0.04146,
			 h1=0.04146,
			 w2=0.03840,
			 h2=0.03840,
			 l=1.2-0.0000001, R0=0.99, Qc=0.0217, alpha=4,m=4.0,W=0.0015)
			 AT (0, 0, 28.6) RELATIVE a1 ROTATED (0,0,0) RELATIVE a1

COMPONENT guide16 = Guide(w1=0.03840,
			 h1=0.03840,
			 w2=0.03519,
			 h2=0.03519,
			 l=1.0-0.0000001, R0=0.99, Qc=0.0217, alpha=4,m=4.0,W=0.0015)
			 AT (0, 0, 29.8) RELATIVE a1 ROTATED (0,0,0) RELATIVE a1

COMPONENT guide17 = Guide(w1=0.03519,
			 h1=0.03519,
			 w2=0.03120,
			 h2=0.03120,
			 l=1.0-0.0000001, R0=0.99, Qc=0.0217, alpha=4,m=4.0,W=0.0015)
			 AT (0, 0, 30.8) RELATIVE a1 ROTATED (0,0,0) RELATIVE a1

COMPONENT slt1 = Slit(
	xwidth = 0.0292, yheight = 0.035)
AT (0, 0, 31.845) RELATIVE a1 /*31.8624 given by present J-PARC staff*/

COMPONENT slt2 = Slit(
    xwidth = 0.0115, yheight = 0.022)
AT (0, 0, 33.615) RELATIVE a1 /*33.6274 measured by J-PARC staff*/

/******************/
/*sample position*/
/****************/

COMPONENT psdsamp = PSD_monitor(
 xwidth = 0.0512, yheight = 0.0512,restore_neutron=1,
 nx = 128, ny = 128, filename = "PSD_samp.txt"
) AT (0,0,34.8025) RELATIVE a1

COMPONENT tofsampl = TOF_monitor(
  xwidth=0.0512, yheight=0.0512, restore_neutron=1, nt=2424, tmin=1, tmax=40000,
  filename="TOF_samp.txt"
) AT (0,0,0) RELATIVE PREVIOUS

COMPONENT lsamp = L_monitor(
  xwidth=0.0512, yheight=0.0512, restore_neutron=1,Lmin=0, Lmax=20, nL=100,
  filename="L_samp.txt"
)AT (0,0,0)RELATIVE PREVIOUS

SPLIT 20 COMPONENT splitpoint = Arm()
AT(0,0,0) RELATIVE PREVIOUS
EXTEND
%{
 ix=iy=iz=-1;
%}

COMPONENT sx_hull=Arm()
AT(0,0,0) RELATIVE PREVIOUS
ROTATED (0,omega,0) RELATIVE Origin
EXTEND
%{
  do {
    double xw,yh,zd;
    xw=nx*MC_GETPAR2(SX,xwidth);
    yh=ny*MC_GETPAR2(SX,yheight);
    zd=nz*MC_GETPAR2(SX,zdepth);
    crystal_miss=0;
    int s;//,oidx;
    double t0,t1;
    if ( (s=box_intersect(&t0,&t1,x,y,z,vx,vy,vz,xw,yh,zd) )){
      /*miss crystal hull or already past it*/
      if(s==0 || t1<=0){
        crystal_miss=1;
        break;
      }

      if(t0>0){
        PROP_DT(t0+DBL_EPSILON);
      }
      /*find out which subcrystal we hit. The dbl_eps are there to protect
      against rounding errors which may occur in the propagation routines so
      the neutron does not switch to the next subcrystal*/
      int iix,iiy,iiz;
      iix = floor((x + xw*0.5)*nx/xw);
      iiy = floor((y + yh*0.5)*ny/yh);
      iiz = floor((z + zd*0.5)*nz/zd);
      if (iiz==nz || iiz<0 || iix==nx || iix<0 || iiy==ny || iiy<0 ){
        /*Neutron is outside the crystal hull. This should never happen*/
        fprintf(stderr,"Warning (%s): Neutron has unpxectedly left the crystal
        hull: (ix,iy,iz) = (%d %d %d) (x,y,z)=(%g %g %g). Terminating\n",
        NAME_CURRENT_COMP,iix,iiy,iiz,x,y,z);
        ABSORB;
      }
      /*avoid getting stuck on a boundary*/
      int count=0;
      if(iix == ix && iiy==iy && iiz==iz){
        /*neutron is still in the same voxel as before - move it to enter the
        next one.*/
        //fprintf(stderr,"I am stuck...hopping to try to get unstuck.\n");
        //fprintf(stderr,"ix,iy,iz= %d %d %d nid=%ld\n",ix,iy,iz,nid);
        PROP_DT(DBL_EPSILON*10.0);
        iix = floor((x + xw*0.5)*nx/xw);
        iiy = floor((y + yh*0.5)*ny/yh);
        iiz = floor((z + zd*0.5)*nz/zd);
        if(iix == ix && iiy==iy && iiz==iz){
          fprintf(stderr,"I am still stuck...terminating neutron.\n");
          ABSORB;
        }
      }

      ix=iix;iy=iiy;iz=iiz;

      //SCATTER;

      /*inject the right hkl_info struct and carry over the "old" neutron ki*/
      /* to skip the orientation search.*/
      /*orientation index is in the 4th col of map_table. -1 since array is*/
      /*zero-indexed.*/
      oidx=Table_Index(*map_table,ix*(ny*nz)+iy*nz+iz,3)-1;
      double kix_carry=MC_GETPAR(SX,hkl_info).kix;
      double kiy_carry=MC_GETPAR(SX,hkl_info).kiy;
      double kiz_carry=MC_GETPAR(SX,hkl_info).kiz;
      MC_GETPAR(SX,hkl_info)=(poly_hkl_info[oidx]);
      MC_GETPAR(SX,hkl_info).kix=kix_carry;
      MC_GETPAR(SX,hkl_info).kiy=kiy_carry;
      MC_GETPAR(SX,hkl_info).kiz=kiz_carry;

      posx=(ix-nx*0.5+.5)*xw/nx;
      posy=(iy-ny*0.5+.5)*yh/ny;
      posz=(iz-nz*0.5+.5)*zd/nz;
      /*translate neutron to account for crystal centre position*/
      x-=posx;y-=posy;z-=posz;

      if(oidx!=-1){
        sx_background=0;
      }else{
        /*This is a void in the polycrystal*/
        sx_background=1;
      }
    }else{
      /*we miss the compound crystal completely so flag that*/
      crystal_miss=1;
    }
  }while(0);
%}
JUMP sx_exit WHEN (crystal_miss)

COMPONENT SX = Single_crystal(
    reflections = lau_fn,
    xwidth = XW/(double)NX, yheight = YH/(double)NY,
    zdepth = ZD/(double)NZ, delta_d_d = 0.01, mosaic = 0.01,sigma_inc=0.4,
    sigma_abs=2.56,barns=1,
    ax = AX, ay = AY, az = AZ, bx = BX, by = BY, bz = BZ, cx = CX,cy = CY,
    cz = CZ,
    p_transmit = 0.9)
WHEN (!sx_background) AT (0, 0, 0) RELATIVE sx_hull
EXTEND
%{
#ifdef MCSTAS_DEBUG
  printf("I have gone through a crystal: scatter=%d, oidx=%d\n",mcScattered,oidx);
#endif
  do {
    if (hkl_info.type='t'){
      //printf("COHSCATTORT: %d\n",oidx);
    }
    if(SCATTERED) xtals++;
    int s;
    double t0,t1;
    if ( (s=box_intersect(&t0,&t1,x,y,z,vx,vy,vz,xwidth,yheight,zdepth) ) &&
    t1>0 ){
      PROP_DT(t1+DBL_EPSILON);
    }
    /*transform back to regular coordinates*/
    x+=posx;y+=posy;z+=posz;
    SCATTER;
  }while(0);
%}

COMPONENT SX_background=Arm()
WHEN (sx_background) AT(0,0,0) RELATIVE sx_hull
EXTEND
%{
#ifdef MCSTAS_DEBUG
  printf("I have gone through background: scatter=%d\n",mcScattered);
#endif
  do {
    int s;
    double t0,t1;
    if ( (s=box_intersect(&t0,&t1,x,y,z,vx,vy,vz,MC_GETPAR2(SX,xwidth),
    MC_GETPAR2(SX,yheight),MC_GETPAR2(SX,zdepth) ) ) && t1>0 ){
      PROP_DT(t1+DBL_EPSILON);
    }
    /*transform back to regular coordinates*/
    x+=posx;y+=posy;z+=posz;
    SCATTER;
  }while(0);
%}

COMPONENT sx_exit = Arm()
AT(0,0,0) RELATIVE sx_hull
EXTEND
%{
  //printf("exit miss: %d\n",crystal_miss);
%}
JUMP sx_hull WHEN (crystal_miss==0)

COMPONENT psd_direct_beam0 = PSD_monitor(
      xwidth=XW,yheight=YH, nx=101, ny=101, filename="psd_direct_beam0",
      restore_neutron=1)
AT(0,0,ZD/2.0+1e-6) RELATIVE psdsamp
GROUP g0

COMPONENT psd_direct_beam_arm = Arm()
AT(0,0,0) RELATIVE PREVIOUS
GROUP g0
EXTEND
%{
    SCATTER;
%}

COMPONENT psd_anton = PSD_monitor(
    xwidth=0.028, yheight=0.028, nx=301, ny=301,filename="psd_anton.dat",
    restore_neutron=1)
AT(0,0,0.2) RELATIVE psdsamp
GROUP g1

COMPONENT psd_anton_arm = Arm()
AT(0,0,0) RELATIVE PREVIOUS
GROUP g1
EXTEND
%{
    SCATTER;
%}
/******************/
/* far field detectors. u = up, d = down */
/****************/

%include "detector_arms.instr"

/* Here's the 4pi detector*/

COMPONENT 4pi_psd=PSD_monitor_4PI(
    nx=301, ny=301,radius=0.2,filename="4pi_psd", restore_neutron=1
)
AT (0,0,0) RELATIVE psdsamp

COMPONENT det_014 = HDB_monitor_new(
        xwidth=0.256, yheight=0.256, filename=det_fn_014, noextrafiles=0)
WHEN (side & 1) AT(0,0,0) RELATIVE det_014a
GROUP det

COMPONENT det_015 = COPY(det_014)(filename=det_fn_015)
WHEN (side & 1) AT (0,0,0) RELATIVE det_015a
GROUP det

COMPONENT det_013 = COPY(det_014)(filename=det_fn_013)
WHEN (side & 1) AT (0,0,0) RELATIVE det_013a
GROUP det

COMPONENT det_011 = COPY(det_014)(filename=det_fn_011)
WHEN (side & 1) AT (0,0,0) RELATIVE det_011a
GROUP det

COMPONENT det_012 = COPY(det_014)(filename=det_fn_012)
WHEN (side & 1) AT (0,0,0) RELATIVE det_012a
GROUP det

COMPONENT det_010 = COPY(det_014)(filename=det_fn_010)
WHEN (side & 1) AT (0,0,0) RELATIVE det_010a
GROUP det

COMPONENT det_008 = COPY(det_014)(filename=det_fn_008)
WHEN (side & 1) AT (0,0,0) RELATIVE det_008a
GROUP det

COMPONENT det_009 = COPY(det_014)(filename=det_fn_009)
WHEN (side & 1) AT (0,0,0) RELATIVE det_009a
GROUP det

COMPONENT det_007 = COPY(det_014)(filename=det_fn_007)
WHEN (side & 1) AT (0,0,0) RELATIVE det_007a
GROUP det

COMPONENT det_005 = COPY(det_014)(filename=det_fn_005)
WHEN (side & 1) AT (0,0,0) RELATIVE det_005a
GROUP det

COMPONENT det_006 = COPY(det_014)(filename=det_fn_006)
WHEN (side & 1) AT (0,0,0) RELATIVE det_006a
GROUP det

COMPONENT det_004 = COPY(det_014)(filename=det_fn_004)
WHEN (side & 1) AT (0,0,0) RELATIVE det_004a
GROUP det

COMPONENT det_002 = COPY(det_014)(filename=det_fn_002)
WHEN (side & 1) AT (0,0,0) RELATIVE det_002a
GROUP det

COMPONENT det_003 = COPY(det_014)(filename=det_fn_003)
WHEN (side & 1) AT (0,0,0) RELATIVE det_003a
GROUP det

COMPONENT det_001 = COPY(det_014)(filename=det_fn_001)
WHEN (side & 1) AT (0,0,0) RELATIVE det_001a
GROUP det


/* N monitors */

COMPONENT det_017 = HDB_monitor_new(
        xwidth=0.256, yheight=0.256, filename=det_fn_017, noextrafiles=0)
WHEN (side & 2) AT(0,0,0) RELATIVE det_017a
GROUP det

COMPONENT det_018 = COPY(det_017)(filename=det_fn_018)
WHEN (side & 2) AT (0,0,0) RELATIVE det_018a
GROUP det

COMPONENT det_016 = COPY(det_017)(filename=det_fn_016)
WHEN (side & 2) AT (0,0,0) RELATIVE det_016a
GROUP det

COMPONENT det_020 = COPY(det_017)(filename=det_fn_020)
WHEN (side & 2) AT (0,0,0) RELATIVE det_020a
GROUP det

COMPONENT det_021 = COPY(det_017)(filename=det_fn_021)
WHEN (side & 2) AT (0,0,0) RELATIVE det_021a
GROUP det

COMPONENT det_019 = COPY(det_017)(filename=det_fn_019)
WHEN (side & 2) AT (0,0,0) RELATIVE det_019a
GROUP det

COMPONENT det_023 = COPY(det_017)(filename=det_fn_023)
WHEN (side & 2) AT (0,0,0) RELATIVE det_023a
GROUP det

COMPONENT det_024 = COPY(det_017)(filename=det_fn_024)
WHEN (side & 2) AT (0,0,0) RELATIVE det_024a
GROUP det

COMPONENT det_022 = COPY(det_017)(filename=det_fn_022)
WHEN (side & 2) AT (0,0,0) RELATIVE det_022a
GROUP det

COMPONENT det_026 = COPY(det_017)(filename=det_fn_026)
WHEN (side & 2) AT (0,0,0) RELATIVE det_026a
GROUP det

COMPONENT det_027 = COPY(det_017)(filename=det_fn_027)
WHEN (side & 2) AT (0,0,0) RELATIVE det_027a
GROUP det

COMPONENT det_025 = COPY(det_017)(filename=det_fn_025)
WHEN (side & 2) AT (0,0,0) RELATIVE det_025a
GROUP det

COMPONENT det_029 = COPY(det_017)(filename=det_fn_029)
WHEN (side & 2) AT (0,0,0) RELATIVE det_029a
GROUP det

COMPONENT det_030 = COPY(det_017)(filename=det_fn_030)
WHEN (side & 2) AT (0,0,0) RELATIVE det_030a
GROUP det

COMPONENT det_028 = COPY(det_017)(filename=det_fn_028)
WHEN (side & 2) AT (0,0,0) RELATIVE det_028a
GROUP det

COMPONENT det_032 = COPY(det_017)(filename=det_fn_032)
WHEN (side & 2) AT (0,0,0) RELATIVE det_032a
GROUP det

COMPONENT det_033 = COPY(det_017)(filename=det_fn_033)
WHEN (side & 2) AT (0,0,0) RELATIVE det_033a
GROUP det

COMPONENT det_031 = COPY(det_017)(filename=det_fn_031)
WHEN (side & 2) AT (0,0,0) RELATIVE det_031a
GROUP det

COMPONENT det_035 = COPY(det_017)(filename=det_fn_035)
WHEN (side & 2) AT (0,0,0) RELATIVE det_035a
GROUP det

COMPONENT det_036 = COPY(det_017)(filename=det_fn_036)
WHEN (side & 2) AT (0,0,0) RELATIVE det_036a
GROUP det

COMPONENT det_034 = COPY(det_017)(filename=det_fn_034)
WHEN (side & 2) AT (0,0,0) RELATIVE det_034a
GROUP det

COMPONENT det_037 = COPY(det_014)(filename=det_fn_037)
WHEN (side & 2) AT(0,0,0) RELATIVE det_037a
GROUP det

COMPONENT catch_all = Arm()
AT(0,0,0) RELATIVE PREVIOUS
GROUP det
EXTEND
%{
  SCATTER;
%}

COMPONENT mon_MCP = FITS_monitor_new(
        xwidth=0.028, yheight=0.028, filename="mon_MCP",short_floats=1,nx=258,
        ny=258,nt=400,tmin=1, tmax=40000,restore_neutron=1
)
WHEN (side & 4) AT(0,0,0.02) RELATIVE psdsamp
GROUP mcp

COMPONENT mon_MCP_not = Monitor_nD(
    radius=0.025, filename="mon_MCP_not",
    options="t bins=400 limits=[1e-6 40000e-6]", restore_neutron=1)
AT(0,0,0) RELATIVE psdsamp
GROUP mcp

COMPONENT mon_MCP_catchall = Arm()
AT(0,0,0) RELATIVE PREVIOUS
GROUP mcp
EXTEND
%{
  SCATTER;
%}


FINALLY
%{
%}

END