| 1 | # flatgen.py |
|---|
| 2 | # S.King |
|---|
| 3 | # |
|---|
| 4 | # Generates a Flat Cell correction file for LOQ or SANS2D Mantid SANS data reduction |
|---|
| 5 | # NB: this is different to a Flat Cell file generated for COLETTE data reduction because |
|---|
| 6 | # the two programs handle the solid angle normalisation differently. |
|---|
| 7 | # |
|---|
| 8 | # THIS VERSION USES SetDetectorOffsets(bank, x, y, z, rot, radius, side) INSTEAD OF MoveInstrumentComponent |
|---|
| 9 | # BECAUSE THE LATTER DOES NOT ALLOW FOR THE FRONT DETECTOR ROT, RADIUS OR SIDE OFFSETS |
|---|
| 10 | # |
|---|
| 11 | # Make the reduction module available |
|---|
| 12 | #from SANSReduction import * |
|---|
| 13 | from ISISCommandInterface import * |
|---|
| 14 | # |
|---|
| 15 | # And need this to use getInstrument(), etc |
|---|
| 16 | from mantidsimple import * |
|---|
| 17 | # |
|---|
| 18 | # Pull in the external fortran algorithm |
|---|
| 19 | # IMPORTANT NOTE: This .pyd file must be located in a directory covered by the environment variable PATH or this script will crash! |
|---|
| 20 | import mantid_flatgen |
|---|
| 21 | |
|---|
| 22 | def FlatGen(verbose,instr,datain,padding,format,flood,quiet,workto,wires,radmin,radmax,xcent,ycent,xpixel,ypixel,spos,rxoff,ryoff,rzoff,fxoff,fyoff,fzoff,froff,fqoff,fsoff): |
|---|
| 23 | |
|---|
| 24 | maxspectra = 73800 |
|---|
| 25 | |
|---|
| 26 | # Convert shifts/offsets to metres for MoveInstrumentComponent later |
|---|
| 27 | spos = spos / 1000 |
|---|
| 28 | # NB: SetDetectorOffsets uses mm not m so don't need the following conversions |
|---|
| 29 | # rxoff = rxoff / 1000 |
|---|
| 30 | # ryoff = ryoff / 1000 |
|---|
| 31 | # rzoff = rzoff / 1000 |
|---|
| 32 | # fxoff = fxoff / 1000 |
|---|
| 33 | # fyoff = fyoff / 1000 |
|---|
| 34 | # fzoff = fzoff / 1000 |
|---|
| 35 | # NB: froff is not currently implemented |
|---|
| 36 | # froff = froff / 1000 |
|---|
| 37 | # fqoff = fqoff / 1000 |
|---|
| 38 | # fsoff = fsoff / 1000 |
|---|
| 39 | |
|---|
| 40 | if (quiet != -1): |
|---|
| 41 | sum_file = workto + 'MAN_' + str(flood) + '-' + str(quiet) + '_SUM_RAW.TXT' |
|---|
| 42 | col_type_flat_file = workto + 'COL_FLAT_CELL_' + str(flood) + '-' + str(quiet) + '.TXT' |
|---|
| 43 | flat_file = workto + 'MAN_FLAT_CELL_' + str(flood) + '-' + str(quiet) + '.TXT' |
|---|
| 44 | else: |
|---|
| 45 | sum_file = workto + 'MAN_' + str(flood) + '_SUM_RAW.TXT' |
|---|
| 46 | col_type_flat_file = workto + 'COL_FLAT_CELL_' + str(flood) + '.TXT' |
|---|
| 47 | flat_file = workto + 'MAN_FLAT_CELL_' + str(flood) + '.TXT' |
|---|
| 48 | |
|---|
| 49 | solid_angle_file = workto + 'MAN_' + str(flood) + '_SOLID_ANGLES.TXT' |
|---|
| 50 | |
|---|
| 51 | workspace0 = 'temporary_data' |
|---|
| 52 | workspace1 = 'TOF_Flood' |
|---|
| 53 | workspace2 = 'RawSum_Flood' |
|---|
| 54 | workspace3 = 'TOF_Quiet' |
|---|
| 55 | workspace4 = 'RawSum_Quiet' |
|---|
| 56 | workspace5 = 'RawSum_FloodMinusQuiet' |
|---|
| 57 | workspace6 = 'RawSumBySpectrum_Flood' |
|---|
| 58 | workspace7 = 'RawSumBySpectrum_FloodMinusQuiet' |
|---|
| 59 | workspace8 = 'COLETTE-type_Flat_Cell' |
|---|
| 60 | workspace9 = 'SolidAnglesBySpectrum' |
|---|
| 61 | workspace10 = 'MANTID-type_Flat_Cell' |
|---|
| 62 | |
|---|
| 63 | # Set instrument |
|---|
| 64 | if instr == 'SANS2D': |
|---|
| 65 | SANS2D() |
|---|
| 66 | instr_code = 2 |
|---|
| 67 | if verbose == 1: |
|---|
| 68 | print 'SANS2D data is expected' |
|---|
| 69 | elif instr == 'LOQ': |
|---|
| 70 | LOQ() |
|---|
| 71 | instr_code = 1 |
|---|
| 72 | if verbose == 1: |
|---|
| 73 | print 'LOQ data is expected' |
|---|
| 74 | else: |
|---|
| 75 | SANS2D() |
|---|
| 76 | instr_code = 2 |
|---|
| 77 | print 'WARNING! The instrument specified does not make sense! Will default to SANS2D' |
|---|
| 78 | |
|---|
| 79 | # Set data directory |
|---|
| 80 | DataPath(datain) |
|---|
| 81 | if verbose == 1: |
|---|
| 82 | print 'Will look for data in',datain |
|---|
| 83 | |
|---|
| 84 | if format == '.raw': |
|---|
| 85 | if verbose == 1: |
|---|
| 86 | print 'RAW files are expected' |
|---|
| 87 | elif format == '.nxs': |
|---|
| 88 | if verbose == 1: |
|---|
| 89 | print 'NXS files are expected' |
|---|
| 90 | else: |
|---|
| 91 | format = ".nxs" |
|---|
| 92 | print 'WARNING! The file format specified does not make sense! Will assume NXS' |
|---|
| 93 | |
|---|
| 94 | print ' ' |
|---|
| 95 | |
|---|
| 96 | # RAW ============================================================================================================================= |
|---|
| 97 | |
|---|
| 98 | if format == '.raw': |
|---|
| 99 | data_file = datain + instr + padding + str(flood) + format |
|---|
| 100 | print 'Processing run',data_file |
|---|
| 101 | print ' ' |
|---|
| 102 | |
|---|
| 103 | LoadRaw(Filename=data_file,OutputWorkspace=workspace1) |
|---|
| 104 | print ' sample position is at:',mtd[workspace1].getInstrument().getSample().getPos() |
|---|
| 105 | if instr == 'LOQ': |
|---|
| 106 | print ' main detector position:',mtd[workspace1].getInstrument().getComponentByName('main-detector-bank').getPos() |
|---|
| 107 | print ' HAB detector position:',mtd[workspace1].getInstrument().getComponentByName('HAB').getPos() |
|---|
| 108 | elif instr == 'SANS2D': |
|---|
| 109 | print ' front detector position:',mtd[workspace1].getInstrument().getComponentByName('front-detector').getPos() |
|---|
| 110 | print ' rear detector position:',mtd[workspace1].getInstrument().getComponentByName('rear-detector').getPos() |
|---|
| 111 | |
|---|
| 112 | print ' ' |
|---|
| 113 | print 'Applying following shifts:' |
|---|
| 114 | if instr == 'LOQ': |
|---|
| 115 | print ' to sample position: 0.0 0.0',spos |
|---|
| 116 | print ' to main detector position:',rxoff,ryoff,rzoff |
|---|
| 117 | # MoveInstrumentComponent(workspace1, 'main-detector-bank', X = rxoff, Y = ryoff, Z = rzoff, RelativePosition="1") |
|---|
| 118 | SetDetectorOffsets('main-detector-bank', rxoff, ryoff, rzoff, 0.0, 0.0, 0.0) |
|---|
| 119 | print ' to HAB detector position:',fxoff,fyoff,fzoff |
|---|
| 120 | # MoveInstrumentComponent(workspace1, 'HAB', X = fxoff, Y = fyoff, Z = fzoff, RelativePosition="1") |
|---|
| 121 | SetDetectorOffsets('HAB', fxoff, fyoff, fzoff, 0.0, 0.0, 0.0) |
|---|
| 122 | |
|---|
| 123 | elif instr == 'SANS2D': |
|---|
| 124 | print ' to sample position: 0.0 0.0',spos |
|---|
| 125 | print ' to front detector position:',fxoff,fyoff,fzoff |
|---|
| 126 | # MoveInstrumentComponent(workspace1, 'front-detector', X = fxoff, Y = fyoff, Z = fzoff, RelativePosition="1") |
|---|
| 127 | SetDetectorOffsets('front', fxoff, fyoff, fzoff, froff, fqoff, fsoff) |
|---|
| 128 | print ' to rear detector position:',rxoff,ryoff,rzoff |
|---|
| 129 | # MoveInstrumentComponent(workspace1, 'rear-detector', X = rxoff, Y = ryoff, Z = rzoff, RelativePosition="1") |
|---|
| 130 | SetDetectorOffsets('rear', rxoff, ryoff, rzoff, 0.0, 0.0, 0.0) |
|---|
| 131 | |
|---|
| 132 | # Re-load file to force offsets to be applied |
|---|
| 133 | # LoadRaw(Filename=data_file,OutputWorkspace=workspace1) |
|---|
| 134 | # MoveInstrumentComponent(workspace1, 'some-sample-holder', X = 0.0, Y = 0.0, Z = spos, RelativePosition="1") |
|---|
| 135 | |
|---|
| 136 | print ' ' |
|---|
| 137 | print ' sample position is now:',mtd[workspace1].getInstrument().getSample().getPos() |
|---|
| 138 | if instr == 'LOQ': |
|---|
| 139 | print ' main detector position now:',mtd[workspace1].getInstrument().getComponentByName('main-detector-bank').getPos() |
|---|
| 140 | print ' HAB detector position now:',mtd[workspace1].getInstrument().getComponentByName('HAB').getPos() |
|---|
| 141 | elif instr == 'SANS2D': |
|---|
| 142 | print ' front detector position now:',mtd[workspace1].getInstrument().getComponentByName('front-detector').getPos() |
|---|
| 143 | print ' rear detector position now:',mtd[workspace1].getInstrument().getComponentByName('rear-detector').getPos() |
|---|
| 144 | |
|---|
| 145 | Integration(workspace1,workspace2) |
|---|
| 146 | |
|---|
| 147 | if (quiet != -1): |
|---|
| 148 | data_file = datain + instr + padding + str(quiet) + format |
|---|
| 149 | print ' ' |
|---|
| 150 | print 'Processing run',data_file |
|---|
| 151 | # Apply position offsets/shifts to the quiet count run for good measure |
|---|
| 152 | if instr == 'LOQ': |
|---|
| 153 | # MoveInstrumentComponent(workspace1, 'main-detector-bank', X = rxoff, Y = ryoff, Z = rzoff, RelativePosition="1") |
|---|
| 154 | SetDetectorOffsets('main-detector-bank', rxoff, ryoff, rzoff, 0.0, 0.0, 0.0) |
|---|
| 155 | # MoveInstrumentComponent(workspace1, 'HAB', X = fxoff, Y = fyoff, Z = fzoff, RelativePosition="1") |
|---|
| 156 | SetDetectorOffsets('HAB', fxoff, fyoff, fzoff, 0.0, 0.0, 0.0) |
|---|
| 157 | elif instr == 'SANS2D': |
|---|
| 158 | # MoveInstrumentComponent(workspace3, 'front-detector', X = fxoff, Y = fyoff, Z = fzoff, RelativePosition="1") |
|---|
| 159 | SetDetectorOffsets('front', fxoff, fyoff, fzoff, froff, fqoff, fsoff) |
|---|
| 160 | # MoveInstrumentComponent(workspace3, 'rear-detector', X = rxoff, Y = ryoff, Z = rzoff, RelativePosition="1") |
|---|
| 161 | SetDetectorOffsets('rear', rxoff, ryoff, rzoff, 0.0, 0.0, 0.0) |
|---|
| 162 | # Tidy this up when SetDetectorOffsets is working properly |
|---|
| 163 | # LoadRaw(Filename=data_file,OutputWorkspace=workspace3) |
|---|
| 164 | # MoveInstrumentComponent(workspace3, 'some-sample-holder', X = 0.0, Y = 0.0, Z = spos, RelativePosition="1") |
|---|
| 165 | |
|---|
| 166 | Integration(workspace3,workspace4) |
|---|
| 167 | Minus(workspace2,workspace4,workspace5) |
|---|
| 168 | SaveRKH(workspace5,sum_file,0) |
|---|
| 169 | else: |
|---|
| 170 | SaveRKH(workspace2,sum_file,0) |
|---|
| 171 | |
|---|
| 172 | print ' ' |
|---|
| 173 | print 'Writing',sum_file |
|---|
| 174 | |
|---|
| 175 | # NEXUS =========================================================================================================================== |
|---|
| 176 | |
|---|
| 177 | elif format == '.nxs': |
|---|
| 178 | data_file = datain + instr + padding + str(flood) + format |
|---|
| 179 | print 'Processing run',data_file |
|---|
| 180 | print ' ' |
|---|
| 181 | |
|---|
| 182 | LoadNexus(Filename=data_file,OutputWorkspace=workspace1) |
|---|
| 183 | print ' sample position is at:',mtd[workspace1].getInstrument().getSample().getPos() |
|---|
| 184 | if instr == 'LOQ': |
|---|
| 185 | print ' main detector position:',mtd[workspace1].getInstrument().getComponentByName('main-detector-bank').getPos() |
|---|
| 186 | print ' HAB detector position:',mtd[workspace1].getInstrument().getComponentByName('HAB').getPos() |
|---|
| 187 | elif instr == 'SANS2D': |
|---|
| 188 | print ' front detector position:',mtd[workspace1].getInstrument().getComponentByName('front-detector').getPos() |
|---|
| 189 | print ' rear detector position:',mtd[workspace1].getInstrument().getComponentByName('rear-detector').getPos() |
|---|
| 190 | |
|---|
| 191 | print ' ' |
|---|
| 192 | print 'Applying following shifts:' |
|---|
| 193 | if instr == 'LOQ': |
|---|
| 194 | print ' to sample position: 0.0 0.0',spos |
|---|
| 195 | print ' to main detector position:',rxoff,ryoff,rzoff |
|---|
| 196 | # MoveInstrumentComponent(workspace1, 'main-detector-bank', X = rxoff, Y = ryoff, Z = rzoff, RelativePosition="1") |
|---|
| 197 | SetDetectorOffsets('main-detector-bank', rxoff, ryoff, rzoff, 0.0, 0.0, 0.0) |
|---|
| 198 | print ' to HAB detector position:',fxoff,fyoff,fzoff |
|---|
| 199 | # MoveInstrumentComponent(workspace1, 'HAB', X = fxoff, Y = fyoff, Z = fzoff, RelativePosition="1") |
|---|
| 200 | SetDetectorOffsets('HAB', fxoff, fyoff, fzoff, 0.0, 0.0, 0.0) |
|---|
| 201 | |
|---|
| 202 | elif instr == 'SANS2D': |
|---|
| 203 | print ' to sample position: 0.0 0.0',spos |
|---|
| 204 | print ' to front detector position:',fxoff,fyoff,fzoff |
|---|
| 205 | # MoveInstrumentComponent(workspace1, 'front-detector', X = fxoff, Y = fyoff, Z = fzoff, RelativePosition="1") |
|---|
| 206 | SetDetectorOffsets('front', fxoff, fyoff, fzoff, froff, fqoff, fsoff) |
|---|
| 207 | print ' to rear detector position:',rxoff,ryoff,rzoff |
|---|
| 208 | # MoveInstrumentComponent(workspace1, 'rear-detector', X = rxoff, Y = ryoff, Z = rzoff, RelativePosition="1") |
|---|
| 209 | SetDetectorOffsets('rear', rxoff, ryoff, rzoff, 0.0, 0.0, 0.0) |
|---|
| 210 | |
|---|
| 211 | # Re-load file to force offsets to be applied |
|---|
| 212 | # LoadNexus(Filename=data_file,OutputWorkspace=workspace1) |
|---|
| 213 | # MoveInstrumentComponent(workspace1, 'some-sample-holder', X = 0.0, Y = 0.0, Z = spos, RelativePosition="1") |
|---|
| 214 | |
|---|
| 215 | print ' ' |
|---|
| 216 | print ' sample position is now:',mtd[workspace1].getInstrument().getSample().getPos() |
|---|
| 217 | if instr == 'LOQ': |
|---|
| 218 | print ' main detector position now:',mtd[workspace1].getInstrument().getComponentByName('main-detector-bank').getPos() |
|---|
| 219 | print ' HAB detector position now:',mtd[workspace1].getInstrument().getComponentByName('HAB').getPos() |
|---|
| 220 | elif instr == 'SANS2D': |
|---|
| 221 | print ' front detector position now:',mtd[workspace1].getInstrument().getComponentByName('front-detector').getPos() |
|---|
| 222 | print ' rear detector position now:',mtd[workspace1].getInstrument().getComponentByName('rear-detector').getPos() |
|---|
| 223 | |
|---|
| 224 | Integration(workspace1,workspace2) |
|---|
| 225 | |
|---|
| 226 | if (quiet != -1): |
|---|
| 227 | data_file = datain + instr + padding + str(quiet) + format |
|---|
| 228 | print ' ' |
|---|
| 229 | print 'Processing run',data_file |
|---|
| 230 | # Apply position offsets/shifts to the quiet count run for good measure |
|---|
| 231 | if instr == 'LOQ': |
|---|
| 232 | # MoveInstrumentComponent(workspace1, 'main-detector-bank', X = rxoff, Y = ryoff, Z = rzoff, RelativePosition="1") |
|---|
| 233 | SetDetectorOffsets('main-detector-bank', rxoff, ryoff, rzoff, 0.0, 0.0, 0.0) |
|---|
| 234 | # MoveInstrumentComponent(workspace1, 'HAB', X = fxoff, Y = fyoff, Z = fzoff, RelativePosition="1") |
|---|
| 235 | SetDetectorOffsets('HAB', fxoff, fyoff, fzoff, 0.0, 0.0, 0.0) |
|---|
| 236 | elif instr == 'SANS2D': |
|---|
| 237 | # MoveInstrumentComponent(workspace3, 'front-detector', X = fxoff, Y = fyoff, Z = fzoff, RelativePosition="1") |
|---|
| 238 | SetDetectorOffsets('front', fxoff, fyoff, fzoff, froff, fqoff, fsoff) |
|---|
| 239 | # MoveInstrumentComponent(workspace3, 'rear-detector', X = rxoff, Y = ryoff, Z = rzoff, RelativePosition="1") |
|---|
| 240 | SetDetectorOffsets('rear', rxoff, ryoff, rzoff, 0.0, 0.0, 0.0) |
|---|
| 241 | # Tidy this up when SetDetectorOffsets is working properly |
|---|
| 242 | # LoadNexus(Filename=data_file,OutputWorkspace=workspace3) |
|---|
| 243 | # MoveInstrumentComponent(workspace3, 'some-sample-holder', X = 0.0, Y = 0.0, Z = spos, RelativePosition="1") |
|---|
| 244 | |
|---|
| 245 | Integration(workspace3,workspace4) |
|---|
| 246 | Minus(workspace2,workspace4,workspace5) |
|---|
| 247 | SaveRKH(workspace5,sum_file,0) |
|---|
| 248 | else: |
|---|
| 249 | SaveRKH(workspace2,sum_file,0) |
|---|
| 250 | |
|---|
| 251 | print ' ' |
|---|
| 252 | print 'Writing',sum_file |
|---|
| 253 | |
|---|
| 254 | # IMPORTANT NOTE: sum_file contains the integrated number of counts per spectrum, summed over all TOF's. If a Quiet Count run was specified it is |
|---|
| 255 | # the DIFFERENCE between the integral for the Flood Source run and the integral for the Quiet Count run. HOWEVER the corresponding Mantid |
|---|
| 256 | # workspaces - RawSum_Flood, RawSum_Quiet or RawSum_FloodMinusQuiet - contain the first and last TOF's as the X data AND NOT the spectrum number! |
|---|
| 257 | # |
|---|
| 258 | # To replace the X data **in Mantid** with a spectral index we need to do one of two things, EITHER: |
|---|
| 259 | # |
|---|
| 260 | # 1) call upon a Mantid user algorithm (ie, not one distributed with Mantid), or |
|---|
| 261 | # 2) use the Transpose algorithm (which was created more recently!) |
|---|
| 262 | # |
|---|
| 263 | # ================================================================================================================================= |
|---|
| 264 | # METHOD 1) |
|---|
| 265 | # |
|---|
| 266 | # if (quiet != -1): |
|---|
| 267 | # ConvertXToSpectrumNumber('RawSum_FloodMinusQuiet', 'RawSumBySpectrum_FloodMinusQuiet') |
|---|
| 268 | # else: |
|---|
| 269 | # ConvertXToSpectrumNumber('RawSum_Flood', 'RawSumBySpectrum_Flood') |
|---|
| 270 | # |
|---|
| 271 | # ================================================================================================================================= |
|---|
| 272 | # Algorithm courtesy of M Gigg |
|---|
| 273 | # IMPORTANT NOTE: This algorithm must be present in C:\MantidInstall\plugins\PythonAlgs or this script will crash! |
|---|
| 274 | # |
|---|
| 275 | # from MantidFramework import * |
|---|
| 276 | # |
|---|
| 277 | # class ConvertXToSpectrumNumber(PythonAlgorithm): |
|---|
| 278 | # |
|---|
| 279 | # def PyInit(self): |
|---|
| 280 | # self.declareWorkspaceProperty("InputWorkspace", "", Direction = Direction.Input) |
|---|
| 281 | # self.declareWorkspaceProperty("OutputWorkspace", "", Direction = Direction.Output) |
|---|
| 282 | # |
|---|
| 283 | # def PyExec(self): |
|---|
| 284 | # input_ws = self.getProperty("InputWorkspace") |
|---|
| 285 | # |
|---|
| 286 | # nhists = input_ws.getNumberHistograms() |
|---|
| 287 | # nbins = 1 |
|---|
| 288 | # output_ws = WorkspaceFactory.createMatrixWorkspaceFromCopy(input_ws, NVectors = nhists, YLength = nbins, XLength = 1) |
|---|
| 289 | # for i in range(nhists): |
|---|
| 290 | # # Set the data in the new workspace |
|---|
| 291 | # output_ws.dataX(i)[0] = i + 1 |
|---|
| 292 | # output_ws.dataY(i)[0] = input_ws.readY(i)[0] |
|---|
| 293 | # output_ws.dataE(i)[0] = input_ws.readE(i)[0] |
|---|
| 294 | # |
|---|
| 295 | # self.setProperty("OutputWorkspace", output_ws) |
|---|
| 296 | # |
|---|
| 297 | # mtd.registerPyAlgorithm(ConvertXToSpectrumNumber()) |
|---|
| 298 | # ================================================================================================================================= |
|---|
| 299 | # METHOD 2) |
|---|
| 300 | |
|---|
| 301 | if (quiet != -1): |
|---|
| 302 | Transpose(workspace5,workspace7) |
|---|
| 303 | else: |
|---|
| 304 | Transpose(workspace2,workspace6) |
|---|
| 305 | |
|---|
| 306 | # ================================================================================================================================= |
|---|
| 307 | |
|---|
| 308 | # Then extract the spectrum number values, the corresponding integrals and the number of bins from the converted workspaces... |
|---|
| 309 | if (quiet != -1): |
|---|
| 310 | RawSumBySpectrum_Workspace = mtd[workspace7] |
|---|
| 311 | else: |
|---|
| 312 | RawSumBySpectrum_Workspace = mtd[workspace6] |
|---|
| 313 | |
|---|
| 314 | xdata = list(RawSumBySpectrum_Workspace.readX(0)) |
|---|
| 315 | ndata = len(xdata) |
|---|
| 316 | ydata = list(RawSumBySpectrum_Workspace.readY(0)) |
|---|
| 317 | edata = list(RawSumBySpectrum_Workspace.readE(0)) |
|---|
| 318 | |
|---|
| 319 | # A small complication is that string and array lengths passed from Python to Fortran must match EXACTLY. Because LOQ and SANS2D have wildy different |
|---|
| 320 | # numbers of spectra (normally 17792 and 73736, respectively) it is not possible to dimension arrays in the Fortran subroutine to suit one case and expect |
|---|
| 321 | # Python to stillbe happy in the other. Thus it is necessary to either pad out a smaller array, or have two seperate Fortran routines doing the same thing. |
|---|
| 322 | # Here we shall pad (with thanks to W S Howells!)... |
|---|
| 323 | |
|---|
| 324 | # An Aside... Yes, we could have combined the previous two operations in: |
|---|
| 325 | # RawSumBySpectrum_Workspace = ConvertXToSpectrum('RawSum_...', 'RawSumBySpectrum_...').workspace() |
|---|
| 326 | |
|---|
| 327 | xarray = PadArray(xdata,maxspectra) |
|---|
| 328 | yarray = PadArray(ydata,maxspectra) |
|---|
| 329 | earray = PadArray(edata,maxspectra) |
|---|
| 330 | |
|---|
| 331 | # ...and pass the data to a wrapped version of RKH's Genie-II function routine FLATGEN to do whatever it does to generate |
|---|
| 332 | # a COLETTE-type Flat Cell file. The wrapping has been achieved with the F2Py utility within NumPy and it creates a Python .pyd file. |
|---|
| 333 | # |
|---|
| 334 | # IMPORTANT NOTE: The .pyd file must be located in a directory covered by the environment variable PATH (eg, drive:/Python25 or drive/Python25/scripts) |
|---|
| 335 | # or this script will crash! |
|---|
| 336 | # |
|---|
| 337 | print ' ' |
|---|
| 338 | print 'Passing', ndata, 'spectra to Mantid_Flatgen.Flatgen...' |
|---|
| 339 | |
|---|
| 340 | retcode,xout,yout,eout,totcnt,navgd,avgds,avgde = mantid_flatgen.flatgen(instr_code,wires,radmin,radmax,xcent,ycent,xpixel,ypixel,ndata,xarray,yarray,earray) |
|---|
| 341 | |
|---|
| 342 | print ' ' |
|---|
| 343 | print 'Retcode is', retcode |
|---|
| 344 | print 'See flatgenlog.txt for details' |
|---|
| 345 | |
|---|
| 346 | print ' ' |
|---|
| 347 | print 'total counts (all spectra) was', "%.3e" % totcnt |
|---|
| 348 | print ' ' |
|---|
| 349 | print navgd, 'spectra were averaged; result was', "%.3e" % avgds, '+/-', "%.3e" % avgde |
|---|
| 350 | |
|---|
| 351 | print ' ' |
|---|
| 352 | print 'Creating workspace...' |
|---|
| 353 | # xout/yout/eout are returned as float32 numpy.ndarray variables but CreateWorkspace requires list arguments |
|---|
| 354 | # CreateWorkspace creates one spectrum with all spectra numbers on its x-axis... |
|---|
| 355 | CreateWorkspace(workspace0,list(xout),list(yout),list(eout),NSpec=1) |
|---|
| 356 | # ...so use Transpose to create many spectra with the spectra axis copied from the original x-axis... |
|---|
| 357 | Transpose(workspace0,workspace0) |
|---|
| 358 | |
|---|
| 359 | # ...and then remove the trailing zeroes from the array padding |
|---|
| 360 | # NB: Cropping on the X value (rather than the index) does not currently (07/12/10) work properly when the X data is not a continuously increasing data set |
|---|
| 361 | CropWorkspace(workspace0,workspace8,StartWorkspaceIndex=0,EndWorkspaceIndex=ndata-1) |
|---|
| 362 | |
|---|
| 363 | print ' ' |
|---|
| 364 | print 'Writing COLETTE-type flat cell file:' |
|---|
| 365 | print col_type_flat_file |
|---|
| 366 | SaveRKH(workspace8,col_type_flat_file,0) |
|---|
| 367 | |
|---|
| 368 | # For Mantid SANS (SANS2D or LOQ) reduction ONLY we now need to 'normalise' the data by the solid angle for each pixel |
|---|
| 369 | # SolidAngle needs to know where it can get the detector information; best place is in a .raw or .nxs file! |
|---|
| 370 | SolidAngle(workspace1,workspace9,StartWorkspaceIndex=0,EndWorkspaceIndex=ndata-1) |
|---|
| 371 | SaveRKH(workspace9,solid_angle_file,0) |
|---|
| 372 | |
|---|
| 373 | print ' ' |
|---|
| 374 | print 'Dividing by solid angle...' |
|---|
| 375 | mtd.deleteWorkspace(workspace0) |
|---|
| 376 | Divide(workspace8,workspace9,workspace0) |
|---|
| 377 | |
|---|
| 378 | print ' ' |
|---|
| 379 | print 'Cleaning up workspace...' |
|---|
| 380 | # Need to tidy up the end of the workspace because 0's will have become infinity! |
|---|
| 381 | ReplaceSpecialValues(workspace0,workspace0,InfinityValue=0.0) |
|---|
| 382 | # SaveRKH(workspace0,workto+'MAN_Unclean_Flat_Cell.txt',0) |
|---|
| 383 | |
|---|
| 384 | # But also what were -1's on return from mantid_flatgen.flatgen have become large negative numbers so need to reset them |
|---|
| 385 | # No easy way to do this within Mantid (because data is in a workspace not an array or list), so will call another external Fortran routine |
|---|
| 386 | # So first package the data... |
|---|
| 387 | Transpose(workspace0,workspace0) |
|---|
| 388 | xdata = list(mtd[workspace0].readX(0)) |
|---|
| 389 | ndata = len(xdata) |
|---|
| 390 | ydata = list(mtd[workspace0].readY(0)) |
|---|
| 391 | edata = list(mtd[workspace0].readE(0)) |
|---|
| 392 | xarray = PadArray(xdata,maxspectra) |
|---|
| 393 | yarray = PadArray(ydata,maxspectra) |
|---|
| 394 | earray = PadArray(edata,maxspectra) |
|---|
| 395 | |
|---|
| 396 | print ' ' |
|---|
| 397 | print 'Passing', ndata, 'spectra to Mantid_Flatgen.Cleanup...' |
|---|
| 398 | # ...then feed it to the external subroutine... |
|---|
| 399 | retcode,xout,yout,eout = mantid_flatgen.cleanup(instr_code,wires,ndata,xarray,yarray,earray) |
|---|
| 400 | |
|---|
| 401 | print ' ' |
|---|
| 402 | print 'Retcode is', retcode |
|---|
| 403 | print 'See flatgenlog.txt for details' |
|---|
| 404 | |
|---|
| 405 | # ...then re-create a new Mantid workspace... |
|---|
| 406 | print ' ' |
|---|
| 407 | print 'Creating workspace...' |
|---|
| 408 | mtd.deleteWorkspace(workspace0) |
|---|
| 409 | CreateWorkspace(workspace0,list(xout),list(yout),list(eout),NSpec=1) |
|---|
| 410 | Transpose(workspace0,workspace0) |
|---|
| 411 | CropWorkspace(workspace0,workspace10,StartWorkspaceIndex=0,EndWorkspaceIndex=ndata-1) |
|---|
| 412 | |
|---|
| 413 | # ...and write it out! |
|---|
| 414 | print ' ' |
|---|
| 415 | print 'Writing Mantid-type flat cell file:' |
|---|
| 416 | print flat_file |
|---|
| 417 | SaveRKH(workspace10,flat_file,0) |
|---|
| 418 | |
|---|
| 419 | print ' ' |
|---|
| 420 | print 'Done' |
|---|
| 421 | |
|---|
| 422 | |
|---|
| 423 | |
|---|
| 424 | |
|---|
| 425 | def PadArray(inarray,nfixed): |
|---|
| 426 | npt=len(inarray) |
|---|
| 427 | padding = nfixed-npt |
|---|
| 428 | outarray=[] |
|---|
| 429 | outarray.extend(inarray) |
|---|
| 430 | outarray +=[0]*padding |
|---|
| 431 | return outarray |
|---|
| 432 | |
|---|
| 433 | |
|---|
| 434 | |
|---|
