Merge pull request e0404#536 from HomolkaN/dev_varRBErobOpt_Update

Extensive TOPAS workflow update:

.gitignore

@@ -1 +1,2 @@
 MCsquare/bin/Materials/
+topas/MCrun

4D/matRad_calc4dDose.m

@@ -84,12 +84,12 @@
     resultGUI.accSqrtBetaDose = matRad_doseAcc(ct,resultGUI.phaseSqrtBetaDose, cst, 'DDM');
 
     % only compute where we have biologically defined tissue
-    ix = dij.alphaX~=0; 
+    ix = dij.ax(:,1)~=0; 
 
     resultGUI.accEffect = resultGUI.accAlphaDose + resultGUI.accSqrtBetaDose.^2;
 
     resultGUI.accRBExD     = zeros(ct.cubeDim);
-    resultGUI.accRBExD(ix) = ((sqrt(dij.alphaX(ix).^2 + 4 .* dij.betaX(ix) .* resultGUI.accEffect(ix)) - dij.alphaX(ix))./(2.*dij.betaX(ix)));
+    resultGUI.accRBExD(ix) = ((sqrt(dij.ax(ix).^2 + 4 .* dij.bx(ix) .* resultGUI.accEffect(ix)) - dij.ax(ix))./(2.*dij.bx(ix)));
 
 end
 

4D/matRad_doseAcc.m

@@ -61,10 +61,10 @@
         error('dose accumulation via direct dose mapping (DDM) requires pull dvfs');
     end    
 
-    [Y,X,Z] = meshgrid(yGridVec,xGridVec,zGridVec);
+    [X,Y,Z] = meshgrid(xGridVec,yGridVec,zGridVec);
 
     % TODODODODODOD %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-    ix = [];1:prod(ct.cubeDim);%[];
+    ix = [];
     for i = 1:size(cst,1)
         ix = unique([ix; cst{i,4}{1}]);
     end
@@ -74,8 +74,8 @@
         dvf_x_i = squeeze(ct.dvf{1,i}(1,:,:,:))/ct.resolution.x;
         dvf_y_i = squeeze(ct.dvf{1,i}(2,:,:,:))/ct.resolution.y;
         dvf_z_i = squeeze(ct.dvf{1,i}(3,:,:,:))/ct.resolution.z;
-        
-        d_ref = matRad_interp3(yGridVec,xGridVec',zGridVec,phaseCubes{i}, ...
+
+        d_ref = matRad_interp3(yGridVec,xGridVec,zGridVec,permute(phaseCubes{i},[2 1 3]), ...
                          Y(ix) + dvf_y_i(ix), ...     
                          X(ix) + dvf_x_i(ix), ... 
                          Z(ix) + dvf_z_i(ix), ...

ChangeLog.txt

@@ -1,3 +1,37 @@
+Monte-Carlo Update
+    Workflow of the Monte Carlo pipeline including MCsquare and TOPAS have been completely overhauled
+    - Restructured the MCEmittanceBaseData class, fit and calculation pipeline
+        - Added example 13 for generating analytical data file by fitting to given machine emittance
+            - added function to plot particleBaseDataEntry
+            - edited function to fit base data
+            - added function to generate a single pencil beam
+        - MCemmittanceBaseData can calculate meanEnergy and spread also for carbon and helium ions
+    - Added support for 4D calculation
+        - added 4D accumulation in calcDoseDirectMC
+        - Fixed number of errors in 4D calculation workflow
+    - TOPAS Updates:
+        - Added comments to whole pipeline
+        - Implemented dij calculation
+        - Restructured resampling in calcParticleDoseMCtopas in separate function
+        - calcParticleDoseMCtopas now generates dij in matRad format for calcDoseDirect
+        - Merged support functions for TOPAS into topasConfig class
+        - modular Schneider Converter: Converter is generated on demand rather than read from file which allows a variety of different options
+        - modular TOPAS scorer which can be individually turned on and off
+        - Export feature for TOPAS to run externally (includes functions to read from external folders) 
+    - MCsquare Updates:
+        - Merged support functions for MCsquare into MCsquareConfig class
+        - added variable in pln to contain an existing BDL file if available
+        - renamed MCsquare property "Num_Primaries" -> "numHistories" to be in line with other Monte Carlo (this is written to the BDL file in original format)
+    - Implemented calculation of std of physicalDose
+
+    Changes to matRad workflow:
+    - Added a class constructor for pln in MatRad_config, which loads requested classes in pln and writes default values that were not set
+    - Changed matRad_calcCubes to accept a variety of different fields for Monte Carlo, without changing the current usage
+    - Enabled helium in calcdoseDirect
+    - Bugfix for coordinate system in resizeCstToGrid function 
+    - Added optional initial weights to fluence Optimization
+    - Added flag in stf to catch specific error where no energies could be found
+
 Development Changes
 	- Bugfixes:
 		- Fix SFUD wrapper for new pln & cst format

MCsquare/bin/Scanners/default/HU_Density_Conversion.txt

@@ -15,4 +15,4 @@
 588		1.335
 1046		1.560
 1565		1.825
-5000		1.825
+5000		1.825

MCsquare/bin/Scanners/default/HU_Material_Conversion.txt

@@ -25,5 +25,4 @@
 1200	59		# Schneider_Marrow_Bone12
 1300	60		# Schneider_Marrow_Bone13
 1400	61		# Schneider_Marrow_Bone14
-1500	62		# Schneider_Marrow_Bone15
-
+1500	62		# Schneider_Marrow_Bone15

MCsquare/MatRad_MCsquareBaseData.m → MCsquare/matRad_MCsquareBaseData.m

@@ -1,5 +1,5 @@
-classdef MatRad_MCsquareBaseData < MatRad_MCemittanceBaseData
-    % MatRad_MCsquareBaseData class for calculating MCsquare base data and
+classdef matRad_MCsquareBaseData < matRad_MCemittanceBaseData
+    % matRad_MCsquareBaseData class for calculating MCsquare base data and
     % writing it to a .txt file, for MCsquare to use
     % 
     % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -16,9 +16,13 @@
     % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
     methods (Access = public)
-        function obj = MatRad_MCsquareBaseData(machine,stf)
-            %Call MatRad_MCemmitanceBaseData constructor
-            obj = obj@MatRad_MCemittanceBaseData(machine,stf);             
+        function obj = matRad_MCsquareBaseData(machine,stf)
+            %Call matRad_MCemmitanceBaseData constructor
+            if nargin < 2
+                stf = [];
+            end
+
+            obj = obj@matRad_MCemittanceBaseData(machine, stf);
         end
 
         function obj = writeMCsquareData(obj,filepath)
@@ -35,7 +39,7 @@
             end
 
             %remove field not needed for MCsquare base data
-            selectedData = rmfield(selectedData, 'FWHMatIso');
+%             selectedData = rmfield(selectedData, 'FWHMatIso');
 
             %write MCsqaure data base file
             try
@@ -77,9 +81,11 @@
                 end
                 fprintf(fileID, '\n');
 
+                indices = obj.selectedFocus(obj.energyIndex);
                 for k = 1:size(selectedData,2)
                     for m = 1:numel(fn)
-                        fprintf(fileID, '%g', selectedData(k).(fn{m}));
+                        tmp = selectedData(k).(fn{m});
+                        fprintf(fileID, '%g ', tmp(indices(k)));
                         fprintf(fileID, '\t');
                     end
                     fprintf(fileID, '\n');