/******************************************************************************* * * McStas, neutron ray-tracing package * Copyright 1997-2002, All rights reserved * Risoe National Laboratory, Roskilde, Denmark * Institut Laue Langevin, Grenoble, France * * Component: PSD_Detector * * %I * Written by: Thorwald van Vuure * Date: Feb 21st, 2005 * Origin: ILL * Modified by: E. Farhi Aug 2008, allow event detector mode to save TOF signal * * Position-sensitive gas-filled detector with gaseous thermal-neutron * converter (box, cylinder or 'banana'). * * %D * An n times m pixel Position Sensitive Detector (PSD), box, cylinder or banana * filled with a mixture of thermal-neutron converter gas and stopping gas. * This model implements wires detection, including charge drift, parallax, etc. * The default is a simple 1D/2D position detector. However, setting * the parameter 'type' to "events" will save the signal as an event file which * contains all neutron parameters at detection points, including time. This is * especially suited for TOF detectors. Please use Histogrammer afterwards. * The gas creation event (neutron absorption in gas) is flagged as SCATTERED==1, * and the charge detection (signal) is flagged as SCATTERED==2. * * GEOMETRY: * A flat rectangular box of given depth is simulated, or alternatively * a cylinder (when giving radius), or a horizontal cylindrical 'banana' * detector (when giving the width along the arc). * Box position is at its center, the 'banana' and cylinder have their * position at focus point, symetric in angular range for the 'banana'. * * GAS SPECIFICATIONS: * The converter gas can be specified as He-3, BF3 or N2. * The filling pressure of converter gas must be specified along * with the zdepth: this gives an absorption probability (e.g. ~90% for * 0.18 nm neutrons at zdepth 3 cm and He-3 pressure 5 bar). * The stopping gas can be specified as one of the following gases: * CF4, C3H8, CO2, Xe/TMA or He, or any for which a table of stopping * powers is available. NB: each stopping gas table has specific * data on a pair of particles (from thermal-neutron absorptions in * either He-3, BF3 or N2) in a specific gas. Unusual choices like N2 * as a converter and He as a stopping gas probably require the * calculation of a new table. * Tip: to do a simulation in pure He-3, specify 0 bar for the stopping gas. * The pressures of the two gases together put a lower limit on the * achievable spatial resolution, no matter what pixel size is * specified. * Example: 1 bar CF4 gives ~2.5 mm spatial resolution, 3 bar CF4 gives * ~1 mm. C3H8 requires more pressure for the same resolution, then * Xe/TMA, then CO2 and finally He. * * IMPLEMENTED FEATURES: * Taken into account are the following effects: * - (rectangularly distributed) error in the recorded position due to * the range of the neutron capture products in the gas * - angular dependence of absorption efficiency * - border effect (events outside the sensitive volume tend to be * counted in the border pixels) * - parallax effect (with, if desired, an electrostatic lens to * partially correct for this - use "LensOn=1") * - wall effect * To correctly take the wall effect into account, a value must be * specified for the energy threshold for acceptance of an event as a * valid neutron event. This normally serves to discriminate from * gammas. Because these are not simulated in McStas, a value of 0 keV * can be chosen for this threshold if one is interested in seeing the * maximum number of neutrons; however, a more realistic value would be * 100 keV. This implies seeing the maximum number of neutrons as well, * except in geometries with dominating wall effect. * The border effect is a considerably complex phenomenon: it depends * on the precise electric field configuration near the border. It is * simulated here simply by introducing a dead zone, gas-filled, * bordering the sensitive volume. Events in there can still cause * signals in the detector. The dead zone can cause a shadow in the * sensitive volume, just as in reality. * The algorithm simply adds all events on the first 'pixel' in the * dead zone to the border pixels. This is a coarse approximation, * because the physical effect depends on the orientation of the * tracks and the energy threshold setting. The bottom line is that * tracking the position of the Center Of Gravity of the charge * deposition in the detector does not allow for accurate modeling of * the border effect, and therefore the border pixels should be * ignored, just as in a real detector. * Note that specifying 0 width of the dead zone next to the * sensitive volume, apart from being unrealistic, does not allow * accurate simulation of the border effect: in this case, there will * be a relative lack of events on the border pixels due to the wall * effect. * This component determines the position of the Center Of Gravity of * the charge cloud in the detector. It does not simulate independent * channels. * Manufacturing imperfections such as wire diameter variations and * spread in electronic amplifier component properties are not simulated. * This should allow an experimenter to identify these manufacturing * imperfections in his physical machine and correct for them. * * Example: PSD_Detector(xwidth=0.192, yheight=0.192, nx=64, ny=64, * zdepth=0.03, threshold=100, borderx=-1, bordery=-1, * PressureConv=5, PressureStop=1, * FN_Conv="Gas_tables/He3inHe.table", FN_Stop="Gas_tables/He3inCF4.table", * xChDivRelSigma=0, yChDivRelSigma=0, * filename="BIDIM19.psd") * * Example: PSD_Detector(xwidth=0.26, yheight=0.26, nx=128, ny=128, * zdepth=0.03, threshold=100, borderx=-1, bordery=-1, * PressureConv=5, PressureStop=1, * FN_Conv="Gas_tables/He3inHe.table", FN_Stop="Gas_tables/He3inCF4.table", * xChDivRelSigma=0, yChDivRelSigma=0, * filename="BIDIM26.psd") * * Example: PSD_Detector(radius=0.0127, yheight=0.20, nx=1, ny=128, * threshold=100, borderx=-1, * PressureConv=9.9, PressureStop=0.1, * FN_Conv="Gas_tables/He3inHe.table", FN_Stop="Gas_tables/He3inCO2.table", * yChDivRelSigma=0.006, * filename="ReuterStokes.psd") * * Example: PSD_Detector(awidth=1.533, yheight=0.4, nx=640, ny=256, * radius=0.73, zdepth=0.032, dc=0.026, threshold=100, * borderx=-1, bordery=-1, * PressureConv=5, PressureStop=1, * FN_Conv="Gas_tables/He3inHe.table", FN_Stop="Gas_tables/He3inCF4.table", * xChDivRelSigma=0, yChDivRelSigma=0.0037, * LensOn=1, filename="D19.psd") * * %P * INPUT PARAMETERS: * * xwidth: [m] Width of detector opening, in the case of a flat box * radius: [m] Radius of detector for cylindrical/banana shape * awidth: [m] Width of detector opening along the horizontal arc of the detector, measured along the inside arc of the sensitive volume in the case of a 'banana' detector * angle: [deg] Angular opening of detector in the case of a 'banana' detector, then overrides awidth parameter * yheight: [m] Height of detector opening * zdepth: [m] Depth of the sensitive volume of the detector * nx: [1] Number of pixel columns * ny: [1] Number of pixel rows * borderx: [m] Width of (gas-filled) border on the side of the detector * bordery: [m] Height of (gas-filled) border on the top/bottom of the detector giving rise to various border effects. Set to -2 to obtain the height of one pixel, -1 for the depth of the detector (default), or specify a non-negativedistance * FN_Conv: [text] Filename of the stopping power table of the converter gas * PressureConv: [bar] Pressure of gaseous thermal-neutron converter (e.g. He-3) * FN_Stop: [text] Filename of the stopping power table of the stopping gas * PressureStop: [bar] Pressure of stopping gas * threshold: [keV] Energy threshold for acceptance of an event * filename: [text] Name of file in which to store the detector image. The name of the component is used if file name is not specified * interpolate: [1] Performs energy interpolation if true * verbose: [1] Gives more information and spectrum * p_interact: [1] Fraction of statistical events to be treated by component. Set it to 0 to use real absorption probability * LensOn: [1] set to 1 to simulate an electrostatic lens according to P. Van Esch, NIM A 540 (2005) pp. 361-367. * dc: [m] Distance from the entrance window to the cathode plane, where the correction zone of the lens ends for banana shape and LensOn=1 * xChDivRelSigma: [1] Sigma of the error in position determination along the x axis (relative to xwidth), uniquely due to charge division readout. Set to 0 in case of independent channel readout along the x dimension * yChDivRelSigma: [1] Relative sigma due to charge division in y direction, relative to yheight. * type: [string] If set to 'events' detector signal saves signal into an event file to be read e.g. by Histogrammer. This is to be used for TOF detectors. * bufsize: [1] Size of neutron output buffer default is 0, i.e. save all - recommended. * restore_neutron: [1] If set, the monitor does not influence the neutron state * * OUTPUT PARAMETERS: * * PSD_N: [] Array of neutron counts * PSD_p: [] Array of neutron weight counts * PSD_p2: [] Array of second moments * SCATERED: [] is 1 when neutron is absorbed by gas, and reaches 2 when charges are detected * * %L * The test/example instrument Test_PSD_Detector.instr. * * %E *******************************************************************************/ DEFINE COMPONENT PSD_Detector DEFINITION PARAMETERS () SETTING PARAMETERS (int nx=128, int ny=128, xwidth=0, radius=0, awidth=0, yheight=0, zdepth=0.03, threshold=100, PressureConv=5, PressureStop=1, interpolate=1, p_interact=0, verbose=0, LensOn=1, dc=0, borderx=-2, bordery=-2, xChDivRelSigma=0, yChDivRelSigma=0, bufsize=0, restore_neutron=0, angle=0, string type=0, string filename=0, string FN_Conv="He3inHe.table", string FN_Stop="He3inCF4.table") OUTPUT PARAMETERS () /* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */ SHARE %{ %include "read_table-lib" %include "monitor_nd-lib" #ifndef PSD_Detector_SHARE #define PSD_Detector_SHARE #pragma acc routine seq double PSD_He_interp1value(double *array_x,double *array_y,long n_elements,double interp_x) { long cnt1 = 0; double a; // while (interp_x > array_x[++cnt1]); while (cnt1 < n_elements-1 && interp_x >= array_x[++cnt1]); a = (interp_x - array_x[cnt1 - 1]) / (array_x[cnt1] - array_x[cnt1 - 1]); return array_y[cnt1 - 1] + a * (array_y[cnt1] - array_y[cnt1 - 1]); } #endif %} DECLARE %{ MonitornD_Defines_type DEFS; /* for event files */ MonitornD_Variables_type Vars; DArray2d PSD_N; DArray2d PSD_p; DArray2d PSD_p2; double *EAP; // Array of energies for the proton from the list, with zero added on top double *EAT; // Array of energies for the triton from the list, with zero added on top double *M1P1; // Centers Of Gravity of transfer of energy by proton double *M1T1; // Centers Of Gravity of transfer of energy by triton double *PosAP; // positions of each packet of energy from the list, in m double *PosAT; // positions of each packet of energy from the list, in m double *PHSpectrum0; // 1 for each neutron double *PHSpectrum; // Pulse height spectrum procuced in the detector, 1 keV/channel double *PHSpectrum2; // 2nd order moments long PHSpectrum_n; double CrossSectionHe; // cross section for 0.18 nm neutrons in the converter double CountNeutrons; // counter over all neutrons to pass the component double GeomCumul; /* 'geometric' cumulative probability that simulated neutrons will encounter the detector on their path */ double AbsCumul; /* cumulative probability that simulated neutrons that encounter the detector on their path will be absorbed */ double SensVolCumul; /* cumulative probability that simulated neutrons that encounter the detector on their path will produce a signal (perhaps small) with COG in the sensitive volume of the detector */ double DetCumul; /* cumulative probability that simulated neutrons will be absorbed *and* produce a signal that will cross the threshold */ long nH_p; // number of table elements for the proton in Helium-3 long nH_t; // number of table elements for the triton in Helium-3 double FullEnergyP; double FullEnergyT; long VariousErrors; /* container for the occurrence of various errors, so that the algorithm gives the matching error messages only once each */ long DetectorType; // 1 for box, 2 for tube, 4 for cylinder ('banana'). Anything else gives an error. double rb; // radius of the back plane of the detector %} INITIALIZE %{ CountNeutrons=0; GeomCumul=0; AbsCumul=0; SensVolCumul=0; DetCumul=0; VariousErrors=0; DetectorType=0; long i,j; t_Table Part_He_3_p; t_Table Part_Stop_p; long nS_p; // number of table elements for the proton in the stopping gas t_Table Part_He_3_t; t_Table Part_Stop_t; long nS_t; // number of table elements for the triton in the stopping gas t_Table Part_He_3_n,Part_Stop_n; // third entry in the file, containing cross section only double *E_p; // energies from the table for protons, for which stopping powers are given double *dEdx_p; // resulting stopping powers of the counting gas mixture double *DAP; // Delta Array of energies for the proton from the list double *MuAP; // Distances in mu traversed for each delta energy in the list double *TempVar34; double *ETP1; // Energy transferred by proton double *PTP1; // Positions of transfer of energy by proton double *E_t; // energies from the table for tritons, for which stopping powers are given double *dEdx_t; // resulting stopping powers of the counting gas mixture double *DAT; // Delta Array of energies for the triton from the list double *MuAT; // Distances in mu traversed for each delta energy in the list double *TempVar35; double *ETT1; // Energy transferred by triton double *PTT1; // Positions of transfer of energy by triton /* GEOMETRY stuff ********************************************************* */ if (borderx==-2) borderx = xwidth/nx; else if (borderx==-1) borderx = zdepth; else if (borderx<0) { fprintf(stderr,"PSD_Detector: %s: Negative x border zone specified. Exit.\n", NAME_CURRENT_COMP); exit(-1); } if (bordery==-2) bordery = yheight/ny; else if (bordery==-1) bordery = zdepth; else if (bordery<0) { fprintf(stderr,"PSD_Detector: %s: Negative y border zone specified. Exit.\n", NAME_CURRENT_COMP); exit(-1); } if (xwidth>0) { /* panel */ DetectorType+=1; if (zdepth<0) { fprintf(stderr,"PSD_Detector: %s: Detector box has no zdepth. Exit.\n", NAME_CURRENT_COMP); exit(-1); } } if (radius && angle>0) awidth=angle/360*2*PI*radius; else if (radius && awidth>0) angle = awidth*360/2/PI/radius; if (radius>0 && !awidth) { DetectorType+=2; /* cylinder */ } if (awidth>0 && radius) { DetectorType+=4; /* banana */ if (zdepth<=0) { fprintf(stderr,"PSD_Detector: %s: Detector 'banana' has no zdepth. Exit.\n", NAME_CURRENT_COMP); exit(-1); } if (radius<=0) { fprintf(stderr,"PSD_Detector: %s: Non-positive radial distance to the sample was specified. Exit.\n", NAME_CURRENT_COMP); exit(-1); } if (!dc && zdepth && LensOn) dc = zdepth/2; if (((dc>zdepth) || (dc<0)) && (LensOn==1)) { fprintf(stderr,"PSD_Detector: %s: Electrostatic lens was turned on, but\n" " no valid zdepth of the cathode plane was specified. Exit.\n", NAME_CURRENT_COMP); exit(-1); } if (awidth+2*borderx>2*PI*radius) { fprintf(stderr,"PSD_Detector: %s: Detector comes perilously close to encompassing 2pi. Exit.\n", NAME_CURRENT_COMP); exit(-1); } } if (yheight<=0) { fprintf(stderr,"PSD_Detector: %s: Detector has no height (yheight). Exit.\n", NAME_CURRENT_COMP); exit(-1); } if ( (DetectorType!=1) && (DetectorType!=2) &&(DetectorType!=4) ) { fprintf(stderr,"PSD_Detector: %s: Detector has conflicting size\n" " specifications, i.e. combinations of xwidth, radius\n" " or awidth, or none at all. Exit.\n", NAME_CURRENT_COMP); exit(-1); } if (verbose) { printf("PSD_Detector: %s: Geometry is ", NAME_CURRENT_COMP); switch (DetectorType) { case 1: printf("box\n"); break; case 2: printf("tube\n"); break; case 4: printf("cylinder ('banana') opening angle=%g [deg] arc length=%g [m]\n", angle, awidth); break; default: printf("not defined\n"); } } /* Gas tables ************************************************************** */ PSD_N = create_darr2d(nx, ny); PSD_p = create_darr2d(nx, ny); PSD_p2 = create_darr2d(nx, ny); Table_Read(&Part_He_3_p,FN_Conv,1); nH_p = Part_He_3_p.rows; Table_Read(&Part_Stop_p,FN_Stop,1); nS_p = Part_Stop_p.rows; Table_Read(&Part_He_3_t,FN_Conv,2); nH_t = Part_He_3_t.rows; Table_Read(&Part_Stop_t,FN_Stop,2); nS_t = Part_Stop_t.rows; /* if Gas tables can not be found, try by pre-pending Gas_tables */ char tmp[256]; if (!nH_p) { sprintf(tmp, "Gas_tables%s%s", MC_PATHSEP_S, FN_Conv); Table_Read(&Part_He_3_p,tmp,1); nH_p = Part_He_3_p.rows; } if (!nS_p) { sprintf(tmp, "Gas_tables%s%s", MC_PATHSEP_S, FN_Stop); Table_Read(&Part_Stop_p,tmp,1); nS_p = Part_Stop_p.rows; } if (!nH_t) { sprintf(tmp, "Gas_tables%s%s", MC_PATHSEP_S, FN_Conv); Table_Read(&Part_He_3_t,tmp,2); nH_t = Part_He_3_t.rows; } if (!nS_t) { sprintf(tmp, "Gas_tables%s%s", MC_PATHSEP_S, FN_Stop); Table_Read(&Part_Stop_t,tmp,2); nS_t = Part_Stop_t.rows; } if (nH_p != nS_p || nH_t != nS_t) { fprintf(stderr,"PSD_Detector: %s: Data files for helium %s and stopping gas %s\n" " have different number of entries. Exit.\n", NAME_CURRENT_COMP, FN_Conv, FN_Stop); exit(-1); } E_p = (double *) malloc(nH_p*sizeof(double)); dEdx_p = (double *) malloc(nH_p*sizeof(double)); for (i=0; i1) p_interact=1; /* handle event files as in Virtual_output, with bufsize=0 ****************** */ long element_size; if (type && strlen(type) && strcmp(type, "NULL") && strcmp(type, "0")) { printf("PSD_Detector: %s: saving detector signal as events\n", NAME_CURRENT_COMP); strcpy(Vars.compcurname, NAME_CURRENT_COMP); if (bufsize > 0) sprintf(Vars.option, "list=%15g", bufsize); else strcpy(Vars.option, "list all"); strcat(Vars.option,", borders, x y z vx vy vz t sx sy sz"); Monitor_nD_Init(&DEFS, &Vars, 0.1, 0.1, 0, 0,0,0,0,0,0,0); /* dims for mcdisplay */ if (filename && strlen(filename) && strcmp(filename, "0") && strcmp(filename, "NULL")) strncpy(Vars.Mon_File, filename, 128); if (bufsize > 0) printf("Warning: PSD_Detector: %s: buffer size=%g not recommended\n", NAME_CURRENT_COMP, bufsize); if (bufsize > 0) printf( "PSD_Detector: %s: Beware virtual output generated file size (max %g Mo)\n" "WARNING Memory required is %g Mo\n", NAME_CURRENT_COMP, bufsize*element_size/1e6, bufsize*sizeof(double)/1e6); } rb=radius+zdepth; // radius of the back plane of the detector %} TRACE %{ long i,j; double RangeProton,RangeTriton; double NAvogadro,MolVolume,v,vt,mu; double length,GM,ZB,pa,RandomNumber,la,ta; double phi,theta; // direction of emission of the particles double vxp,vyp,vzp; // normalized speed of the particles double mlp,mlt; // 1st moment (COG) of energy deposition of the particles double xCOG,yCOG,zCOG,aCOG,rCOG,Energy; double tin,tout,tine,toute,tinb,toutb; // times in and out of cylindrical entrance window and back plane boxes long ISe,ISb; // type of intersection with a cylindrical box double EnergyP,EnergyT,SigmaEnergy; double lp,lt; // length proton- and triton tracks long intersect; long Impossibility; // used to check for impossible collisions double alimit; // borders of the sensitive volume + edges double aa; // alpha coordinate (along the arc of the detector) of the absorption double t1,t2; // time for the first and second charged particles to the next wall #define HUGE_VAL DBL_MAX #pragma acc atomic CountNeutrons = CountNeutrons + p; RangeProton = PosAP[nH_p]; RangeTriton = PosAT[nH_t]; v = sqrt(vx*vx + vy*vy + vz*vz) ; // speed of the neutron if (v>440000 && !(VariousErrors&1) ) { fprintf(stderr,"PSD_Detector: %s: Component cannot deal with\n" " high-energy neutrons. Absorbing.\n", NAME_CURRENT_COMP); #pragma acc atomic VariousErrors = VariousErrors + 1; ABSORB; } vt = sqrt( 2*0.025243*1.60218e-19 / ( 1.0086649 * 1.66053886e-27 ) ); // speed of 25.243 meV neutron (0.18 nm) NAvogadro = 6.022045e23; // Number of atoms per mol MolVolume = 24.7796e-3; // in m3, 24.7796 liter per mol at 1 bar at T=298 K mu = (NAvogadro*PressureConv*CrossSectionHe* vt ) / (MolVolume* v ); // in inverse m switch(DetectorType) { case 1: intersect = box_intersect(&tin,&tout,x,y,z,vx,vy,vz,xwidth+2*borderx,yheight+2*bordery,zdepth); if (tin==tout) intersect=0; break; case 2: intersect = cylinder_intersect(&tin,&tout,x,y,z,vx,vy,vz,radius,yheight+2*bordery); if (tin==tout) intersect=0; break; case 4: if (cylinder_intersect(&tine,&toute,x,y,z,vx,vy,vz,radius,yheight+2*bordery)==0 ) { ISe=4; } // case for no intersection with the entrance-window-cylinder else { if (tine<=0 && toute<=0) { ISe=1; // case for neutron outside the volume and moving away } if (tine<=0 && toute>0) { ISe=2; // case for neutron inside the volume and moving (obviously) out } if (tine>0 && toute>0) { ISe=3; // case for neutron outside the volume and moving towards it } if (tine>toute && !(VariousErrors&2) ) { fprintf(stderr,"PSD_Detector: %s: Something is seriously wrong.\n" " cylinder_intersect reported a later time\n" " for entering the cylinder than for leaving it.\n", NAME_CURRENT_COMP); #pragma acc atomic VariousErrors = VariousErrors + 2; ABSORB; } } if (cylinder_intersect(&tinb,&toutb,x,y,z,vx,vy,vz,rb,yheight+2*bordery)==0 || tinb==toutb) { ISb=4; } // case for no intersection with the back-plane-cylinder else { if (tinb<=0 && toutb<=0) { ISb=1; // case for neutron outside the volume and moving away } if (tinb<=0 && toutb>0) { ISb=2; // case for neutron inside the volume and moving (obviously) out } if (tinb>0 && toutb>0) { ISb=3; // case for neutron outside the volume and moving towards it } if (tinb>toutb && !(VariousErrors&2) ) { fprintf(stderr,"PSD_Detector: %s: Something is seriously wrong.\n" " cylinder_intersect reported a later time\n" " for entering the cylinder than for leaving it.\n", NAME_CURRENT_COMP); #pragma acc atomic VariousErrors = VariousErrors + 2; ABSORB; } } /******************************************************************** * Schematic representations of the possibilities for the path of the neutron * ----------------------------------------------------------------- * | ISb=1 ISb=2 ISb=3 ISb=4 * ----------------------------------------------------------------- * ISe=1 | Gone already in Impossible Impossible * | out thr. back * ----------------------------------------------------------------- * ISe=2 | Impossible in thr. entr. Impossible Impossible * | out thr. back * ----------------------------------------------------------------- * ISe=3 | Impossible already in in thr. back Impossible * | out thr. entr. out thr. entr. * ----------------------------------------------------------------- * ISe=4 | Gone already in in thr. back Gone * | out thr. back out thr. back * ----------------------------------------------------------------- * All we have to do is choose the appropriate entrance and exit times * from the tin and tout of both the entrance window cylinder and the * back plane cylinder. * There is a problem on the left and right sides of the detector. * We cannot check entrance or exit using cylinder_intersect there, so * we simply find interaction points in 2pi cylinder wall, and then * reject events that are too far from the sensitive volume. *********************************************************************/ Impossibility=1; // set this to 1: if it is not subsequently set to 0, an impossible neutron track has occurred if ( ((ISb==1) && ((ISe==1)||(ISe==4))) || ((ISb==4)&&(ISe==4)) ) { Impossibility=0; // these represent the cases that the neutron misses the detector intersect=0; } if ( (ISe==1)&&(ISb==2) ) { Impossibility = 0; intersect=1; tin=0; // neutron is already located within the sensitive volume tout=toutb; // neutron flies out the sensitive volume into the back plane } if ( (ISe==2)&&(ISb==2) ) { Impossibility = 0; intersect=1; tin=toute; // neutron comes in through the entrance window tout=toutb; // and flies out through the back plane } if ( (ISe==3)&&(ISb==2) ) { Impossibility = 0; intersect=1; tin=0; // neutron is already located within the sensitive volume tout=tine; // neutron flies out the sensitive volume into the entrance window } if ( (ISe==3)&&(ISb==3) ) { Impossibility = 0; intersect=1; tin=tinb; // neutron comes in through the back plane tout=tine; // and flies out through the entrance window } if ( (ISe==4)&&(ISb==2) ) { Impossibility = 0; intersect=1; tin=0; // neutron is already located within the sensitive volume tout=toutb; // neutron flies out the sensitive volume into the back plane } if ( (ISe==4)&&(ISb==3) ) { Impossibility = 0; intersect=1; tin=tinb; // neutron comes in through the back plane tout=toutb; // and flies out through the back plane } if ( tin==tout ) { intersect=0; /* This happens for instance when the neutron flies along the axis of the cylinders. The two cylinders overlap there, so zero time is spent in the gas volume, therefore no absorption. */ } if ( tin>tout && !(VariousErrors&4) ) { intersect=0; fprintf(stderr,"PSD_Detector: %s: A serious error occurred.\n" " A later time for entering the detector was calculated\n" " than for leaving it.\n", NAME_CURRENT_COMP); #pragma acc atomic VariousErrors = VariousErrors + 4; } if ( Impossibility==1 && !(VariousErrors&8) ) { fprintf(stderr,"PSD_Detector: %s: Something strange happened. A neutron\n" " followed a path deemed impossible in the algorithm.\n", NAME_CURRENT_COMP); intersect=0; #pragma acc atomic VariousErrors = VariousErrors + 8; } /* NB: at this point the neutron does not yet have to be absorbed within the sensitive volume of the detector - we still have to check polar angle to see if the event falls within the part of the cylinder wall that is covered by the detector. */ break; /* end switch case DetectorType==4 */ default: exit(-1); } // end switch DetectorType (1) if (intersect && tout<=0) { intersect=0; /* This is the case that the intersection along the trajectory was in the past, which is probably a common occurrance in a simulation, e.g. backscattering off the entrance window, and which should not be interpreted as a new detector event.*/ } if (intersect && tin<=0 && !(VariousErrors&16) ) { tin=0; fprintf(stderr,"PSD_Detector: %s: Warning: a neutron has been found\n" " 'tunneled' into the detector, rather than just entering\n" " through the entrance window or one of the other sides.\n", NAME_CURRENT_COMP); #pragma acc atomic VariousErrors = VariousErrors + 16; } if (intersect) { /* there is intersection, not 'grazing' of the detector. What kind of trajectory is followed is not important, since we have in- and out-going time. */ length = v*(tout-tin); // length (in m) of the trajectory through the sensitive volume GM = mu*length ; // number of radiation lengths in the gap thickness ZB = exp( -GM ); // zone boundary for generating uniformly distributed random numbers [ZB,1] pa = 1 - ZB; // Probability of absorption somewhere along the length of the trajectory through the sensitive volume if ( (p_interact <= 0 && (pa >= 1 || rand01() < pa)) || (p_interact > 0 && (p_interact >= 1 || rand01() < p_interact)) ) { if (p_interact > 0 && p_interact < 1) pa /= p_interact; if (p_interact <= 0) pa=1; RandomNumber = (1-ZB)*rand01() + ZB; // uniformly distributed random number between ZB and 1 la = -log(RandomNumber) / mu; // (in m) Absorption location (along trajectory), distributed exponentially declining ta = tin + la/v; // Absorption time alimit= (0.5*awidth + borderx)/radius; /* we take the gas-filled volume a bit larger than the sensitive volume alone, which is realistic. Later we throw away events that do not end up in the sensitive volume. */ aa=atan2(x+vx*ta,z+vz*ta); if ( DetectorType==1 || DetectorType==2 || (DetectorType==4 && aa>-alimit && aa=0) { t1=tin; t2=tout; } else { t1=tout; t2=tin; } break; case 2: intersect=cylinder_intersect(&tin,&tout,x,y,z,vxp,vyp,vzp,radius,yheight+2*bordery); if (tin>=0) { t1=tin; t2=tout; } else { t1=tout; t2=tin; } break; case 4: /****************************************************************** * What we have now is an absorption point somewhere between the two * cylinder walls, with distance 'zdepth' of the sensitive volume, and * a direction of emission of the reaction products. * What we need is very simple: for the first reaction product (the * proton for absorption in He-3) we need the smallest positive time, * for the second reaction product (triton) we choose the smallest * negative time. ******************************************************************/ if (cylinder_intersect(&tinb,&toutb,x,y,z,vxp,vyp,vzp,rb,yheight+2*bordery)==0 && !(VariousErrors&32) ) { fprintf(stderr,"PSD_Detector: %s: Oops. A neutron absorbed in the detector" " failed to send its reaction products towards the detector walls.\n", NAME_CURRENT_COMP); #pragma acc atomic VariousErrors = VariousErrors + 32; ABSORB; } if (cylinder_intersect(&tine,&toute,x,y,z,vxp,vyp,vzp,radius,yheight+2*bordery) ==0) { tine=-HUGE_VAL; // ugly fix for cases where the inner cylinder is missed toute=HUGE_VAL; } t1=tine; if (t1<0) { t1=toute; /* We are looking for a positive number. If t1 is negative, then choose the next number, toute, no matter what it is. */ } else { if ((toute=0)) { t1=toute; /* If the next number, toute, is positive yet smaller than t1, we want it. */ } } if (t1<0) { t1=tinb; /* We are looking for a positive number. If t1 is negative, then choose the next number, tinb, no matter what it is. */ } else { if ((tinb=0)) { t1=tinb; /* If the next number, tinb, is positive yet smaller than t1, we want it. */ } } if (t1<0) { t1=toutb; /* We are looking for a positive number. If t1 is negative, then choose the next number, toutb, no matter what it is. */ } else { if ((toutb=0)) { t1=toutb; /* If the next number, toutb, is positive yet smaller than t1, we want it. */ } } t2=tine; if (t2>0) { t2=toute; /* We are looking for a negative number. If t2 is positive, then choose the next number, toute, no matter what it is. */ } else { if ((toute>t2)&&(toute<=0)) { t2=toute; /* If the next number, toute, is negative yet larger than t2, we want it. */ } } if (t2>0) { t2=tinb; /* We are looking for a negative number. If t2 is positive, then choose the next number, tinb, no matter what it is. */ } else { if ((tinb>t2)&&(tinb<=0)) { t2=tinb; /* If the next number, tinb, is negative yet larger than t2, we want it. */ } } if (t2>0) { t2=toutb; /* We are looking for a negative number. If t2 is positive, then choose the next number, toutb, no matter what it is. */ } else { if ((toutb>t2)&&(toutb<=0)) { t2=toutb; /* If the next number, toutb, is negative yet larger than t2, we want it. */ } } break; default: fprintf(stderr,"PSD_Detector: %s: Detector has conflicting size\n" " specifications, i.e. combinations of xwidth, radius\n" " or awidth, or none at all. Exit.\n", NAME_CURRENT_COMP); exit(-1); } // end switch detectortype (2) if (t1*t2>=0 && !(VariousErrors&64) ) { fprintf(stderr,"PSD_Detector: %s: t1 was %g and t2 was %g.\n" " One is supposed to be negative, the other positive,\n" " neither zero.\n", NAME_CURRENT_COMP,t1,t2); #pragma acc atomic VariousErrors = VariousErrors + 64; } if ( t1= 0 && i < nx && j >= 0 && j < ny ) { SensVolCumul += p*pa; /* probability that a signal was produced in the sensitive volume, though not necessarily above the threshold. Note that, AT THIS POINT, the p has not been multiplied by pa yet. */ if (Energy>threshold) { // if the particles have deposited sufficient energy to be // recognised as a neutron event (and not a gamma event) double P = p*pa; double P2 = P*P; x = xCOG; y = yCOG; z = zCOG; SCATTER; // show point of detection for 3D view #pragma acc atomic PSD_N[i][j] = PSD_N[i][j] + 1; // one neutron tallied in the appropriate pixel #pragma acc atomic PSD_p[i][j] = PSD_p[i][j] + P; // incremented by neutron weight #pragma acc atomic PSD_p2[i][j] = PSD_p2[i][j] + P2; // 2nd order moments DetCumul += P; /* probability that a signal was produced in the sensitive volume above the energy threshold */ i = (long)floor(Energy); if (i >= 0 && i < PHSpectrum_n ) { double P2 = P*P; #pragma acc atomic PHSpectrum0[i] = PHSpectrum0[i] + 1; // one neutron tallied per pixel #pragma acc atomic PHSpectrum[i] = PHSpectrum[i] + P; // energy bin in pulse height spectrum incremented by neutron weight #pragma acc atomic PHSpectrum2[i] = PHSpectrum2[i] + P2; // 2nd order moments } /* save event file if activated */ if (type && strlen(type) && strcmp(type, "NULL") && strcmp(type, "0")) { double pp,P2; /* Vars.cp = P; Vars.cx = x; Vars.cvx = vx; Vars.csx = sx; Vars.cy = y; Vars.cvy = vy; Vars.csy = sy; Vars.cz = z; Vars.cvz = vz; Vars.csz = sz; Vars.ct = t; */ P2=P*P; pp = Monitor_nD_Trace(&DEFS, &Vars, _particle); #pragma acc atomic Vars.Nsum = Vars.Nsum + 1; #pragma acc atomic Vars.psum = Vars.psum + P; #pragma acc atomic Vars.p2sum = Vars.p2sum + P2; } // initial version ABSORB detected neutrons. // This was removed if user wants to analyze the behaviour of the detector } } /* storage */ p *= 1-pa; } // end if absorption in case of border effects (left and right) in the banana } /* end if p_interact */ } /* end if (intersect) */ if (restore_neutron) { RESTORE_NEUTRON(INDEX_CURRENT_COMP, x, y, z, vx, vy, vz, t, sx, sy, sz, p); } %} SAVE %{ char file[128]; if (type && strlen(type) && strcmp(type, "NULL") && strcmp(type, "0")) { /* event file */ Monitor_nD_Save(&DEFS, &Vars); DETECTOR_OUT(Vars.Nsum, Vars.psum, Vars.p2sum); } else { double width; if (xwidth>0) width=xwidth; if (radius>0) width=2*radius; if (awidth>0) width=awidth; if (filename && strlen(filename) && strcmp(filename, "0") && strcmp(filename, "NULL")) strncpy(file, filename, 128); else sprintf(file, "%s.dat", NAME_CURRENT_COMP); if (nx > 1 && ny > 1) DETECTOR_OUT_2D( "PSD Detector", "X position [m]", "Y position [m]", -width/2, width/2, -yheight/2, yheight/2, (double)nx, (double)ny, &PSD_N[0][0],&PSD_p[0][0],&PSD_p2[0][0], file); else if (nx == 1) DETECTOR_OUT_1D( "PSD Detector", "Y position [m]","Counts","Y", -yheight/2, yheight/2, (double)ny, &PSD_N[0][0],&PSD_p[0][0],&PSD_p2[0][0], file); else if (ny == 1) DETECTOR_OUT_1D( "PSD Detector", "X position [m]","Counts","X", -width/2, width/2, (double)nx, &PSD_N[0][0],&PSD_p[0][0],&PSD_p2[0][0], file); } if (verbose) { char xlabelstr[64]; snprintf(file, 128, "%s.en", filename && strlen(filename) && strcmp(filename,"0") && strcmp(filename,"NULL") ? filename : NAME_CURRENT_COMP); fprintf(stdout,"PSD_Detector: %s: statistics\n", NAME_CURRENT_COMP); fprintf(stdout," %g neutrons in the simulation, %g encounter the detector.\n", CountNeutrons, GeomCumul ); fprintf(stdout," Probability for a neutron to be absorbed in the detector is %g percent.\n", 100*AbsCumul/CountNeutrons ); fprintf(stdout," Probability for a neutron to be detected is %g percent.\n", 100*DetCumul/CountNeutrons ); fprintf(stdout," Fraction of neutrons not counted because their COG is outside\n" "the sensitive volume is %g.\n", 1-SensVolCumul/AbsCumul ); fprintf(stdout," Fraction of neutrons not counted because their signal is reduced below\n" "%g keV due to the wall effect is %g.\n", threshold,1-DetCumul/SensVolCumul ); fprintf(stdout," Theoretical limit to the position resolution in this gas is %g m FWHM\n" " (of the rectangular distribution). This implies a sigma of %g m.\n", 2*(M1P1[nH_p]-M1T1[nH_t]),(M1P1[nH_p]-M1T1[nH_t])/sqrt(3) ); snprintf(xlabelstr,64,"Energy [keV], threshold set to %g keV",threshold); DETECTOR_OUT_1D( "Pulse Height Spectrum", xlabelstr, "Counts [a.u]", "E", 0.0, (double)(PHSpectrum_n-1), PHSpectrum_n, PHSpectrum0, PHSpectrum, PHSpectrum2, file); } %} FINALLY %{ /* free pointers */ if (type && strlen(type) && strcmp(type, "NULL") && strcmp(type, "0")) { Monitor_nD_Finally(&DEFS, &Vars); if (bufsize) { printf("PSD_Detector: %s: Saved %lld events (from buffer) in file %s\n", NAME_CURRENT_COMP, Vars.Nsum, Vars.Mon_File); if (bufsize < Vars.Nsum) printf("WARNING When using this source, intensities must be multiplied\n" " by a factor %g\n", (double)Vars.Nsum/bufsize); } else printf("PSD_Detector: %s: Saved %lld events (all) in file %s\n", NAME_CURRENT_COMP, Vars.Nsum, Vars.Mon_File); } destroy_darr2d(PSD_N); destroy_darr2d(PSD_p); destroy_darr2d(PSD_p2); %} MCDISPLAY %{ double h; h=yheight; if (xwidth>0) { /* box */ double xmin = -0.5*xwidth; double xmax = 0.5*xwidth; double ymin = -0.5*yheight; double ymax = 0.5*yheight; double zmin = -0.5*zdepth; double zmax = 0.5*zdepth; multiline(5, xmin, ymin, zmin, xmax, ymin, zmin, xmax, ymax, zmin, xmin, ymax, zmin, xmin, ymin, zmin); multiline(5, xmin, ymin, zmax, xmax, ymin, zmax, xmax, ymax, zmax, xmin, ymax, zmax, xmin, ymin, zmax); line(xmin, ymin, zmin, xmin, ymin, zmax); line(xmax, ymin, zmin, xmax, ymin, zmax); line(xmin, ymax, zmin, xmin, ymax, zmax); line(xmax, ymax, zmin, xmax, ymax, zmax); } if (radius>0) { if (!awidth) { /* cylinder */ circle("xz", 0, h/2.0, 0, radius); circle("xz", 0, -h/2.0, 0, radius); line(-radius, -h/2.0, 0, -radius, +h/2.0, 0); line(+radius, -h/2.0, 0, +radius, +h/2.0, 0); line(0, -h/2.0, -radius, 0, +h/2.0, -radius); line(0, -h/2.0, +radius, 0, +h/2.0, +radius); if (borderx>0){ circle("xz", 0, h/2.0, 0, radius+borderx); circle("xz", 0, -h/2.0, 0, radius+borderx); } } else { int NH=24; int ih; for(ih = 0; ih < NH; ih++) { double phi0, phi1; double x0,y0,z0,x1,y1,z1,x2,y2,z2,x3,y3,z3; phi0 = (-angle/2+(angle/NH)*ih) *DEG2RAD; /* in xz plane */ phi1 = (-angle/2+(angle/NH)*(ih+1))*DEG2RAD; z0 = radius*cos(phi0); x0 = radius*sin(phi0); y0 = -yheight/2; z1 = radius*cos(phi0); x1 = radius*sin(phi0); y1 = yheight/2; z2 = radius*cos(phi1); x2 = radius*sin(phi1); y2 = y1; z3 = radius*cos(phi1); x3 = radius*sin(phi1); y3 = y0; mcdis_multiline(5, x0,y0,z0, x1,y1,z1, x2,y2,z2, x3,y3,z3, x0,y0,z0); } } } %} END