diff --git a/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/FilterEvents.h b/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/FilterEvents.h index 8e491e7b5179..53bfc7a93c5e 100644 --- a/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/FilterEvents.h +++ b/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/FilterEvents.h @@ -148,8 +148,6 @@ namespace Algorithms /// TOF detector/sample correction type TOFCorrectionType m_tofCorrType; - /// TOF detector/sample correcton operation - TOFCorrectionOp m_tofCorrOperation; }; diff --git a/Code/Mantid/Framework/Algorithms/src/FilterEvents.cpp b/Code/Mantid/Framework/Algorithms/src/FilterEvents.cpp index d2fd4b289a6e..98d8afa23c1e 100644 --- a/Code/Mantid/Framework/Algorithms/src/FilterEvents.cpp +++ b/Code/Mantid/Framework/Algorithms/src/FilterEvents.cpp @@ -152,14 +152,6 @@ namespace Algorithms setPropertySettings("IncidentEnergy", new VisibleWhenProperty("CorrectionToSample", IS_EQUAL_TO, "Direct")); - vector corrops; - corrops.push_back("Shift"); - corrops.push_back("Multiply"); - declareProperty("CustumCorrectionOp", "Multiply", boost::make_shared(corrops), - "Type of operation for TOF correction of customized type. "); - setPropertySettings("CustumCorrectionOp", - new VisibleWhenProperty("CorrectionToSample", IS_EQUAL_TO, "Customized")); - declareProperty("SplitSampleLogs", true, "If selected, all sample logs will be splitted by the " "event splitters. It is not recommended for fast event log splitters. "); @@ -309,14 +301,6 @@ namespace Algorithms m_detCorrectWorkspace = getProperty("DetectorTOFCorrectionWorkspace"); if (!m_detCorrectWorkspace) throw runtime_error("In case of customized TOF correction, correction workspace must be given!"); - - string corrop = getPropertyValue("CustumCorrectionOp"); - if (corrop.compare("Shift") == 0) - m_tofCorrOperation = ShiftOp; - else if (corrop.compare("Multiply") == 0) - m_tofCorrOperation = MultiplyOp; - else - throw runtime_error("Impossible situation for CustomCorrectionOp!"); } m_splitSampleLogs = getProperty("SplitSampleLogs"); @@ -688,7 +672,12 @@ namespace Algorithms if (par) { efix = par->value(); - g_log.debug() << "Detector: " << det->getID() << " EFixed: " << efix << "\n"; + g_log.debug() << "Detector: " << det->getID() << " of spectrum " << i << " EFixed: " << efix << "\n"; + } + else + { + g_log.warning() << "Detector: " << det->getID() << " of spectrum " << i << " does not have EFixed set up." + << "\n"; } } catch (std::runtime_error&) @@ -705,6 +694,8 @@ namespace Algorithms // Calculate shift shift = -1. * l2 / sqrt(efix * twomev_d_mass); + + g_log.notice() << "Detector " << i << ": " << "L2 = " << l2 << ", EFix = " << efix << ", Shift = " << shift << "\n"; } else { @@ -728,30 +719,64 @@ namespace Algorithms //---------------------------------------------------------------------------------------------- /** Set up corrections with customized TOF correction input + * The first column must be either DetectorID or Spectrum (from 0... as workspace index) + * The second column must be Correction or CorrectFactor, a number between 0 and 1, i.e, [0, 1] + * The third column is optional as shift in unit of second */ void FilterEvents::setupCustomizedTOFCorrection() { - // Check input workspace vector colnames = m_detCorrectWorkspace->getColumnNames(); + bool hasshift = false; if (colnames.size() < 2) throw runtime_error("Input table workspace is not valide."); - else if (colnames[0].compare("DetectorID") || colnames[1].compare("Correction")) - throw runtime_error("Input table workspace has wrong column definition."); + else if (colnames.size() >= 3) + hasshift = true; + + bool usedetid; + if (colnames[0].compare("DetectorID") == 0) usedetid = true; + else if (colnames[0].compare("Spectrum") == 0) usedetid = false; + else + { + usedetid = false; + stringstream errss; + errss << "First column must be either DetectorID or Spectrum. " << colnames[0] + << " is not supported. "; + throw runtime_error(errss.str()); + } // Parse detector and its TOF offset (i.e., correction) to a map map correctmap; + map shiftmap; size_t numrows = m_detCorrectWorkspace->rowCount(); for (size_t i = 0; i < numrows; ++i) { TableRow row = m_detCorrectWorkspace->getRow(i); + // Parse to map detid_t detid; - double offset; + double offset, shift; row >> detid >> offset; + if (offset >= 0 && offset <= 1) + { + // Valid offset (factor value) + correctmap.insert(make_pair(detid, offset)); + } + else + { + // Error, throw! + stringstream errss; + errss << "Correction (i.e., offset) equal to " << offset << " of row " << "is out of range [0, 1]."; + throw runtime_error(errss.str()); + } - correctmap.insert(make_pair(detid, offset)); - } + // Shift + if (hasshift) + { + row >> shift; + shiftmap.insert(make_pair(detid, shift)); + } + } // ENDFOR(row i) // Check size of TOF correction map size_t numhist = m_eventWS->getNumberHistograms(); @@ -767,54 +792,87 @@ namespace Algorithms errss << "Input correction table workspace has more detectors (" << correctmap.size() << ") than input workspace " << m_eventWS->name() << "'s spectra number (" << numhist << ".\n"; - g_log.error(errss.str()); throw runtime_error(errss.str()); } - // Get vector IDs - vector vecDetIDs; - vecDetIDs.resize(numhist, 0); - // Set up the detector IDs to vecDetIDs and set up the initial value - for (size_t i = 0; i < numhist; ++i) + // Apply to m_detTofOffsets and m_detTofShifts + if (usedetid) { - // It is assumed that there is one detector per spectra. - // If there are more than 1 spectrum, it is very likely to have problem with correction factor - const DataObjects::EventList events = m_eventWS->getEventList(i); - std::set detids = events.getDetectorIDs(); - std::set::iterator detit; - if (detids.size() != 1) + // Get vector IDs + vector vecDetIDs; + vecDetIDs.resize(numhist, 0); + // Set up the detector IDs to vecDetIDs and set up the initial value + for (size_t i = 0; i < numhist; ++i) { - // Check whether there are more than 1 detector per spectra. - stringstream errss; - errss << "The assumption is that one spectrum has one and only one detector. " - << "Error is found at spectrum " << i << ". It has " << detids.size() << " detectors."; - throw runtime_error(errss.str()); + // It is assumed that there is one detector per spectra. + // If there are more than 1 spectrum, it is very likely to have problem with correction factor + const DataObjects::EventList events = m_eventWS->getEventList(i); + std::set detids = events.getDetectorIDs(); + std::set::iterator detit; + if (detids.size() != 1) + { + // Check whether there are more than 1 detector per spectra. + stringstream errss; + errss << "The assumption is that one spectrum has one and only one detector. " + << "Error is found at spectrum " << i << ". It has " << detids.size() << " detectors."; + throw runtime_error(errss.str()); + } + detid_t detid = 0; + for (detit=detids.begin(); detit!=detids.end(); ++detit) + detid = *detit; + vecDetIDs[i] = detid; } - detid_t detid = 0; - for (detit=detids.begin(); detit!=detids.end(); ++detit) - detid = *detit; - vecDetIDs[i] = detid; - } - // Map correction map to list - map::iterator fiter; - for (size_t i = 0; i < numhist; ++i) + // Map correction map to list + map::iterator fiter; + for (size_t i = 0; i < numhist; ++i) + { + detid_t detid = vecDetIDs[i]; + // correction (factor) map + fiter = correctmap.find(detid); + if (fiter != correctmap.end()) m_detTofOffsets[i] = fiter->second; + else + { + stringstream errss; + errss << "Detector " << "w/ ID << " << detid << " of spectrum " << i << " in Eventworkspace " << m_eventWS->name() + << " cannot be found in input TOF calibration workspace. "; + throw runtime_error(errss.str()); + } + // correction shift map + fiter = shiftmap.find(detid); + if (fiter != shiftmap.end()) + m_detTofShifts[i] = fiter->second; + } // ENDFOR (each spectrum i) + } + else { - detid_t detid = vecDetIDs[i]; - fiter = correctmap.find(detid); - if (fiter != correctmap.end()) + // It is spectrum ID already + map::iterator fiter; + // correction factor + for (fiter = correctmap.begin(); fiter != correctmap.end(); ++fiter) { - m_detTofOffsets[i] = fiter->second; + size_t wsindex = static_cast(fiter->first); + if (wsindex < numhist) m_detTofOffsets[wsindex] = fiter->second; + else + { + stringstream errss; + errss << "Workspace index " << wsindex << " is out of range."; + throw runtime_error(errss.str()); + } } - else + // correction shift + for (fiter = shiftmap.begin(); fiter != shiftmap.end(); ++fiter) { - stringstream errss; - errss << "Detector " << "w/ ID << " << detid << " of spectrum " << i << " in Eventworkspace " << m_eventWS->name() - << " cannot be found in input TOF calibration workspace. "; - g_log.error(errss.str()); - throw runtime_error(errss.str()); + size_t wsindex = static_cast(fiter->first); + if (wsindex < numhist) m_detTofShifts[wsindex] = fiter->second; + else + { + stringstream errss; + errss << "Workspace index " << wsindex << " is out of range."; + throw runtime_error(errss.str()); + } } - } // ENDFOR (each spectrum i) + } return; } diff --git a/Code/Mantid/Framework/Algorithms/test/FilterEventsTest.h b/Code/Mantid/Framework/Algorithms/test/FilterEventsTest.h index 54a7ab23d3f5..8ba80c2c5292 100644 --- a/Code/Mantid/Framework/Algorithms/test/FilterEventsTest.h +++ b/Code/Mantid/Framework/Algorithms/test/FilterEventsTest.h @@ -275,7 +275,6 @@ class FilterEventsTest : public CxxTest::TestSuite TS_ASSERT_THROWS_NOTHING(filter.setProperty("OutputWorkspaceBaseName", "SplittedDataX")); TS_ASSERT_THROWS_NOTHING(filter.setProperty("CorrectionToSample", "Customized")); TS_ASSERT_THROWS_NOTHING(filter.setProperty("DetectorTOFCorrectionWorkspace", "TimeCorrectionTableX")); - TS_ASSERT_THROWS_NOTHING(filter.setProperty("CustumCorrectionOp", "Multiply")); TS_ASSERT_THROWS_NOTHING(filter.setProperty("SplitterWorkspace", splws)); // 3. Execute @@ -315,6 +314,118 @@ class FilterEventsTest : public CxxTest::TestSuite return; } + //---------------------------------------------------------------------------------------------- + /** Test filtering with correction of direct geometry + */ + void test_FilterDGCorrection() + { + DataObjects::EventWorkspace_sptr ws = createEventWorkspaceDirect(0, 1000000); + AnalysisDataService::Instance().addOrReplace("MockDirectEventWS", ws); + + MatrixWorkspace_sptr splws = createMatrixSplittersDG(); + AnalysisDataService::Instance().addOrReplace("SplitterTableX", splws); + + // Run the filtering + FilterEvents filter; + filter.initialize(); + + filter.setProperty("InputWorkspace", ws->name()); + filter.setProperty("OutputWorkspaceBaseName", "SplittedDataDG"); + filter.setProperty("CorrectionToSample", "Direct"); + filter.setProperty("SplitterWorkspace", "SplitterTableX"); + + TS_ASSERT_THROWS_NOTHING(filter.execute()); + TS_ASSERT(filter.isExecuted()); + + // Check + EventWorkspace_sptr ws5 = boost::dynamic_pointer_cast( + AnalysisDataService::Instance().retrieve("SplittedDataDG_5")); + TS_ASSERT(ws5); + if (ws5) + { + TS_ASSERT_EQUALS(ws5->getNumberEvents(), 0); + } + + EventWorkspace_sptr ws7 = boost::dynamic_pointer_cast( + AnalysisDataService::Instance().retrieve("SplittedDataDG_7")); + TS_ASSERT(ws7); + if (ws7) + { + TS_ASSERT_EQUALS(ws7->getNumberEvents(), ws7->getNumberHistograms()); + } + + // FIXME - Should find a way to delete all workspaces holding splitted events + + AnalysisDataService::Instance().remove("MockDirectEventWS"); + AnalysisDataService::Instance().remove("SplitterTableX"); + + return; + } + + //---------------------------------------------------------------------------------------------- + /** Test filtering with correction to indirect geometry inelastic instrument + */ + void test_FilterIndirectGeometryCorrection() + { + // Create workspaces for filtering + DataObjects::EventWorkspace_sptr ws = createEventWorkspaceInDirect(0, 1000000); + AnalysisDataService::Instance().addOrReplace("MockIndirectEventWS", ws); + + MatrixWorkspace_sptr splws = createMatrixSplittersDG(); + AnalysisDataService::Instance().addOrReplace("SplitterTableX", splws); + + // Run the filtering + FilterEvents filter; + filter.initialize(); + + filter.setProperty("InputWorkspace", "MockIndirectEventWS"); + filter.setProperty("OutputWorkspaceBaseName", "SplittedDataDG"); + filter.setProperty("CorrectionToSample", "Indirect"); + filter.setProperty("SplitterWorkspace", "SplitterTableX"); + filter.setProperty("OutputTOFCorrectionWorkspace", "MockIndGeoCorrWS"); + + TS_ASSERT_THROWS_NOTHING(filter.execute()); + TS_ASSERT(filter.isExecuted()); + + // Check + MatrixWorkspace_sptr outcorrws = boost::dynamic_pointer_cast( + AnalysisDataService::Instance().retrieve("MockIndGeoCorrWS")); + TS_ASSERT(outcorrws); + if (outcorrws) + { + TS_ASSERT_EQUALS(outcorrws->getNumberHistograms(), ws->getNumberHistograms()); + TS_ASSERT_EQUALS(outcorrws->readX(0).size(), 2); + + Kernel::V3D samplepos = ws->getInstrument()->getSample()->getPos(); + + for (size_t iws = 0; iws < outcorrws->getNumberHistograms(); ++iws) + { + const ParameterMap& pmap = ws->constInstrumentParameters(); + + IDetector_const_sptr det = ws->getDetector(iws); + Kernel::V3D detpos = det->getPos(); + Parameter_sptr par = pmap.getRecursive(det.get(),"Efixed"); + double efix = par->value(); + + double l2 = samplepos.distance(detpos); + + double shift = -l2/sqrt(efix*2.*PhysicalConstants::meV/PhysicalConstants::NeutronMass); + std::cout << "Detector " << iws << ": L2 = " << l2 << ", EFix = " << efix << "\n"; + + TS_ASSERT_DELTA(outcorrws->readY(iws)[0], 1., 1.0E-9); + TS_ASSERT_DELTA(outcorrws->readY(iws)[1], shift, 1.0E-9); + + } + } + + // Clean + // FIXME - Need to write the code to remove all the workspaces holding split events + AnalysisDataService::Instance().remove("MockIndirectEventWS"); + + return; + + } + //---------------------------------------------------------------------------------------------- /** Create an EventWorkspace. This workspace has @@ -356,6 +467,153 @@ class FilterEventsTest : public CxxTest::TestSuite } + //---------------------------------------------------------------------------------------------- + /** Create an EventWorkspace to mimic direct inelastic scattering insturment. + * This workspace will have the same neutron events as the test case in EventList + * + * @param runstart_i64 : absolute run start time in int64_t format with unit nanosecond + * @param pulsedt : pulse length in int64_t format with unit nanosecond + */ + DataObjects::EventWorkspace_sptr createEventWorkspaceDirect(int64_t runstart_i64, int64_t pulsedt) + { + // Create an EventWorkspace with 10 banks with 1 detector each. No events is generated + DataObjects::EventWorkspace_sptr eventWS = + WorkspaceCreationHelper::createEventWorkspaceWithFullInstrument(10, 1, true); + + // L1 = 10 + Kernel::V3D samplepos = eventWS->getInstrument()->getSample()->getPos(); + Kernel::V3D sourcepos = eventWS->getInstrument()->getSource()->getPos(); + std::cout << "sample position: " << samplepos.toString() << "\n"; + std::cout << "source position: " << sourcepos.toString() << "\n"; + double l1 = samplepos.distance(sourcepos); + + Kernel::DateAndTime runstart(runstart_i64); + + EventList fakeevlist = fake_uniform_time_sns_data(runstart_i64, pulsedt); + + // Set properties: (1) run_start time; (2) Ei + eventWS->mutableRun().addProperty("run_start", runstart.toISO8601String(), true); + + double shift = 2.E-4; + double ei = (l1*l1*PhysicalConstants::NeutronMass) / (shift*shift*2.*PhysicalConstants::meV); + + eventWS->mutableRun().addProperty("Ei", ei, true); + + // Add neutrons + for (size_t i = 0; i < eventWS->getNumberHistograms(); i ++) + { + DataObjects::EventList* elist = eventWS->getEventListPtr(i); + + for (size_t ievent = 0; ievent < fakeevlist.getNumberEvents(); ++ievent) + { + TofEvent tofevent = fakeevlist.getEvent(ievent); + elist->addEventQuickly(tofevent); + // std::cout << "Add event " << ievent << " as " << tofevent.tof() << ". Size of elist = " << elist->getNumberEvents() << "\n"; + } // FOR each pulse + } // For each bank + + double constshift = l1/sqrt(ei * 2. * PhysicalConstants::meV / PhysicalConstants::NeutronMass); + std::cout << "Fake direct inelastic scattering: L1 = " << l1 << ", Shift = " << shift << ": Ei = " + << ei << "; check shift = " << constshift << "\n" + << "- Number of events = " << eventWS->getNumberEvents() << ", Number of spectra = " + << eventWS->getNumberHistograms() << "\n" + << "- Number of faked events = " << fakeevlist.getNumberEvents() << "\n"; + + return eventWS; + } + + //---------------------------------------------------------------------------------------------- + /** Create an EventWorkspace to mimic direct inelastic scattering insturment. This workspace has + * @param runstart_i64 : absolute run start time in int64_t format with unit nanosecond + * @param pulsedt : pulse length in int64_t format with unit nanosecond + */ + DataObjects::EventWorkspace_sptr createEventWorkspaceInDirect(int64_t runstart_i64, int64_t pulsedt) + { + // Create an EventWorkspace with 10 banks with 1 detector each. No events is generated + DataObjects::EventWorkspace_sptr eventWS = + WorkspaceCreationHelper::createEventWorkspaceWithFullInstrument(10, 1, true); + + + // Add EFixed to each detector + const ParameterMap& pmap = eventWS->constInstrumentParameters(); + + for (size_t i = 0; i < 10; ++i) + { + Geometry::IDetector_const_sptr det = eventWS->getDetector(i); + if (det) + std::cout << "detector " << i << ": " << det->getPos().toString() << "\n"; + else + std::cout << "detector of workspace index " << i << " is not defined. " << "\n"; + + Parameter_sptr par = pmap.getRecursive(det.get(),"Efixed"); + if (par) + { + double efix = par->value(); + std::cout << "Detector: " << det->getID() << " EFixed: " << efix << "\n"; + } + else + { + std::cout << "Detector: " << det->getID() << " has no EFixed set up. Set to 2.08" << "\n"; + eventWS->setEFixed(det->getID(), 2.08); + } + } + + // Add neutrons + EventList fakeevlist = fake_uniform_time_sns_data(runstart_i64, pulsedt); + for (size_t i = 0; i < eventWS->getNumberHistograms(); i ++) + { + DataObjects::EventList* elist = eventWS->getEventListPtr(i); + + for (size_t ievent = 0; ievent < fakeevlist.getNumberEvents(); ++ievent) + { + TofEvent tofevent = fakeevlist.getEvent(ievent); + elist->addEventQuickly(tofevent); + } // FOR each pulse + } // For each bank + + return eventWS; + } + + //---------------------------------------------------------------------------------------------- + /** Create an EventWorkspace. This workspace has + * @param runstart_i64 : absolute run start time in int64_t format with unit nanosecond + * @param pulsedt : pulse length in int64_t format with unit nanosecond + * @param todft : time interval between 2 adjacent event in same pulse in int64_t format of unit nanosecond + * @param numpulses : number of pulses in the event workspace + */ + DataObjects::EventWorkspace_sptr createEventWorkspaceElastic(int64_t runstart_i64, int64_t pulsedt, int64_t tofdt, + size_t numpulses) + { + // 1. Create an EventWorkspace with 10 detectors + DataObjects::EventWorkspace_sptr eventWS = + WorkspaceCreationHelper::createEventWorkspaceWithFullInstrument(10, 1, true); + + Kernel::DateAndTime runstart(runstart_i64); + + // 2. Set run_start time + eventWS->mutableRun().addProperty("run_start", runstart.toISO8601String(), true); + + for (size_t i = 0; i < eventWS->getNumberHistograms(); i ++) + { + DataObjects::EventList* elist = eventWS->getEventListPtr(i); + + for (int64_t pid = 0; pid < static_cast(numpulses); pid ++) + { + int64_t pulsetime_i64 = pid*pulsedt+runstart.totalNanoseconds(); + Kernel::DateAndTime pulsetime(pulsetime_i64); + for (size_t e = 0; e < 10; e ++) + { + double tof = static_cast(e*tofdt/1000); + DataObjects::TofEvent event(tof, pulsetime); + elist->addEventQuickly(event); + } + } // FOR each pulse + } // For each bank + + return eventWS; + } + + //---------------------------------------------------------------------------------------------- /** Create a Splitter for output * Region: @@ -463,6 +721,62 @@ class FilterEventsTest : public CxxTest::TestSuite return corrtable; } + //---------------------------------------------------------------------------------------------- + /** Fake uniform time data more close to SNS case + * A list of 1000 events + * Pulse length: 1000000 * nano-second + */ + EventList fake_uniform_time_sns_data(int64_t runstart, int64_t pulselength) + { + //Clear the list + EventList el = EventList(); + + //Create some mostly-reasonable fake data. + srand(1234); //Fixed random seed + for (int time = 0; time < 1000; time++) + { + //All pulse times from 0 to 999 in seconds + Kernel::DateAndTime pulsetime(static_cast(time*pulselength + runstart)); + el += TofEvent( rand()%1000, pulsetime ); //Kernel::DateAndTime(time*1.0, 0.0) ); + } + + return el; + } + + API::MatrixWorkspace_sptr createMatrixSplittersDG() + { + MatrixWorkspace_sptr spws = boost::dynamic_pointer_cast( + WorkspaceFactory::Instance().create("Workspace2D", 1, 11, 10)); + + MantidVec& vec_splitTimes = spws->dataX(0); + MantidVec& vec_splitGroup = spws->dataY(0); + + vec_splitTimes[0] = 1000000; + vec_splitTimes[1] = 1300000; // Rule in 1,339,000 + vec_splitTimes[2] = 2000000; + vec_splitTimes[3] = 2190000; // Rule out 2,155,000 + vec_splitTimes[4] = 4000000; + vec_splitTimes[5] = 5000000; + vec_splitTimes[6] = 5500000; // Rule in 5,741,000 + vec_splitTimes[7] = 7000000; + vec_splitTimes[8] = 8000000; + vec_splitTimes[9] = 9000000; + vec_splitTimes[10] = 10000000; + + + vec_splitGroup[0] = 2; + vec_splitGroup[1] = 5; + vec_splitGroup[2] = 4; + vec_splitGroup[3] = -1; + vec_splitGroup[4] = 6; + vec_splitGroup[5] = 7; + vec_splitGroup[6] = 8; + vec_splitGroup[7] = -1; + vec_splitGroup[8] = 1; + vec_splitGroup[9] = 3; + + return spws; + } };