| Geant4 Cross Reference |
1 /* 1 /*
2 # <<BEGIN-copyright>> 2 # <<BEGIN-copyright>>
3 # <<END-copyright>> 3 # <<END-copyright>>
4 */ 4 */
5 5
6 #include <iostream> 6 #include <iostream>
7 #include <cmath> 7 #include <cmath>
8 8
9 #include "GIDI_settings.hh" 9 #include "GIDI_settings.hh"
10 10
11 using namespace GIDI; 11 using namespace GIDI;
12 12
13 /* 13 /*
14 ============================================== 14 =========================================================
15 */ 15 */
16 GIDI_settings_particle::GIDI_settings_particle 16 GIDI_settings_particle::GIDI_settings_particle( int PoPId, bool transporting, int energyMode ) : mGroup( ) {
17 17
18 initialize( PoPId, transporting, energyMod 18 initialize( PoPId, transporting, energyMode );
19 } 19 }
20 /* 20 /*
21 ============================================== 21 =========================================================
22 */ 22 */
23 GIDI_settings_particle::GIDI_settings_particle 23 GIDI_settings_particle::GIDI_settings_particle( GIDI_settings_particle const &particle ) {
24 24
25 initialize( particle.mPoPId, particle.mTra 25 initialize( particle.mPoPId, particle.mTransporting, particle.mEnergyMode );
26 setGroup( particle.mGroup ); 26 setGroup( particle.mGroup );
27 for( std::vector<GIDI_settings_processedFl 27 for( std::vector<GIDI_settings_processedFlux>::const_iterator iter = particle.mProcessedFluxes.begin( ); iter != particle.mProcessedFluxes.end( ); ++iter ) {
28 mProcessedFluxes.push_back( *iter ); 28 mProcessedFluxes.push_back( *iter );
29 } 29 }
30 } 30 }
31 /* 31 /*
32 ============================================== 32 =========================================================
33 */ 33 */
34 int GIDI_settings_particle::initialize( int Po 34 int GIDI_settings_particle::initialize( int PoPId, bool transporting, int energyMode ) {
35 35
36 mPoPId = PoPId; 36 mPoPId = PoPId;
37 mTransporting = transporting; 37 mTransporting = transporting;
38 int energyMode_ = ( energyMode & GIDI_se 38 int energyMode_ = ( energyMode & GIDI_settings_projectileEnergyMode_continuousEnergy )
39 + ( energyMode & GIDI_se 39 + ( energyMode & GIDI_settings_projectileEnergyMode_grouped )
40 // + ( energyMode & GIDI_se 40 // + ( energyMode & GIDI_settings_projectileEnergyMode_fixedGrid ) // Currently not supported.
41 ; 41 ;
42 42
43 if( energyMode_ != energyMode ) throw 1; 43 if( energyMode_ != energyMode ) throw 1;
44 mEnergyMode = energyMode; 44 mEnergyMode = energyMode;
45 45
46 mGroupX = NULL; 46 mGroupX = NULL;
47 setGroup( mGroup ); 47 setGroup( mGroup );
48 return( 0 ); 48 return( 0 );
49 } 49 }
50 /* 50 /*
51 ============================================== 51 =========================================================
52 */ 52 */
53 void GIDI_settings_particle::setGroup( GIDI_se 53 void GIDI_settings_particle::setGroup( GIDI_settings_group const &group ) {
54 54
55 nfu_status status_nf; 55 nfu_status status_nf;
56 56
57 mGroup = group; 57 mGroup = group;
58 58
59 if( mGroupX != NULL ) ptwX_free( mGroupX ) 59 if( mGroupX != NULL ) ptwX_free( mGroupX );
60 mGroupX = NULL; 60 mGroupX = NULL;
61 if( mGroup.size( ) > 0 ) { 61 if( mGroup.size( ) > 0 ) {
62 if( ( mGroupX = ptwX_create( (int) mGr 62 if( ( mGroupX = ptwX_create( (int) mGroup.size( ), (int) mGroup.size( ), mGroup.pointer( ), &status_nf ) ) == NULL ) throw 1;
63 } 63 }
64 } 64 }
65 /* 65 /*
66 ============================================== 66 =========================================================
67 */ 67 */
68 GIDI_settings_particle::~GIDI_settings_particl 68 GIDI_settings_particle::~GIDI_settings_particle( ) {
69 69
70 if( mGroupX != NULL ) ptwX_free( mGroupX ) 70 if( mGroupX != NULL ) ptwX_free( mGroupX );
71 } 71 }
72 /* 72 /*
73 ============================================== 73 =========================================================
74 */ 74 */
75 int GIDI_settings_particle::addFlux( statusMes 75 int GIDI_settings_particle::addFlux( statusMessageReporting* /* smr */, GIDI_settings_flux const &flux ) {
76 76
77 double temperature = flux.getTemperature( 77 double temperature = flux.getTemperature( );
78 std::vector<GIDI_settings_processedFlux>:: 78 std::vector<GIDI_settings_processedFlux>::iterator iter;
79 79
80 for( iter = mProcessedFluxes.begin( ); ite 80 for( iter = mProcessedFluxes.begin( ); iter != mProcessedFluxes.end( ); ++iter ) {
81 if( temperature <= iter->getTemperatur 81 if( temperature <= iter->getTemperature( ) ) break;
82 } 82 }
83 // BRB need to check if temperature is the sam 83 // BRB need to check if temperature is the same.
84 mProcessedFluxes.insert( iter, GIDI_settin 84 mProcessedFluxes.insert( iter, GIDI_settings_processedFlux( flux, mGroupX ) );
85 return( 0 ); 85 return( 0 );
86 } 86 }
87 /* 87 /*
88 ============================================== 88 =========================================================
89 */ 89 */
90 GIDI_settings_processedFlux const *GIDI_settin 90 GIDI_settings_processedFlux const *GIDI_settings_particle::nearestFluxToTemperature( double temperature ) const {
91 91
92 double priorTemperature, lastTemperature; 92 double priorTemperature, lastTemperature;
93 std::vector<GIDI_settings_processedFlux>:: 93 std::vector<GIDI_settings_processedFlux>::const_iterator iter;
94 94
95 if( mProcessedFluxes.size( ) == 0 ) return 95 if( mProcessedFluxes.size( ) == 0 ) return( NULL );
96 96
97 priorTemperature = mProcessedFluxes[0].get 97 priorTemperature = mProcessedFluxes[0].getTemperature( );
98 //TK adds next line 98 //TK adds next line
99 lastTemperature = mProcessedFluxes[0].get 99 lastTemperature = mProcessedFluxes[0].getTemperature( );
100 for( iter = mProcessedFluxes.begin( ); ite 100 for( iter = mProcessedFluxes.begin( ); iter != mProcessedFluxes.end( ); ++iter ) {
101 lastTemperature = iter->getTemperature 101 lastTemperature = iter->getTemperature( );
102 if( lastTemperature > temperature ) br 102 if( lastTemperature > temperature ) break;
103 //TK add next line 103 //TK add next line
104 priorTemperature = iter->getTemperatur 104 priorTemperature = iter->getTemperature( );
105 } 105 }
106 if( iter == mProcessedFluxes.end( ) ) { 106 if( iter == mProcessedFluxes.end( ) ) {
107 --iter; } 107 --iter; }
108 else { 108 else {
109 //if( fabs( lastTemperature - temperat 109 //if( fabs( lastTemperature - temperature ) < fabs( temperature - priorTemperature ) ) --iter;
110 //TK modified above line 110 //TK modified above line
111 if( std::fabs( lastTemperature - tempe 111 if( std::fabs( lastTemperature - temperature ) > std::fabs( temperature - priorTemperature ) ) --iter;
112 } 112 }
113 return( &(*iter) ); 113 return( &(*iter) );
114 } 114 }
115 /* 115 /*
116 ============================================== 116 =========================================================
117 */ 117 */
118 ptwXPoints *GIDI_settings_particle::groupFunct 118 ptwXPoints *GIDI_settings_particle::groupFunction( statusMessageReporting *smr, ptwXYPoints *ptwXY1, double temperature, int order ) const {
119 119
120 if( mGroupX == NULL ) return( NULL ); 120 if( mGroupX == NULL ) return( NULL );
121 GIDI_settings_processedFlux const *process 121 GIDI_settings_processedFlux const *processedFlux = nearestFluxToTemperature( temperature );
122 if( processedFlux == NULL ) return( NULL ) 122 if( processedFlux == NULL ) return( NULL );
123 return( processedFlux->groupFunction( smr, 123 return( processedFlux->groupFunction( smr, mGroupX, ptwXY1, order ) );
124 } 124 }
125 125
126 /* ---- GIDI_settings_processedFlux ---- */ 126 /* ---- GIDI_settings_processedFlux ---- */
127 /* 127 /*
128 ============================================== 128 =========================================================
129 */ 129 */
130 GIDI_settings_processedFlux::GIDI_settings_pro 130 GIDI_settings_processedFlux::GIDI_settings_processedFlux( GIDI_settings_flux const &flux, ptwXPoints *groupX ) : mFlux( flux ) {
131 131
132 nfu_status status_nf; 132 nfu_status status_nf;
133 ptwXYPoints *fluxXY = NULL; 133 ptwXYPoints *fluxXY = NULL;
134 ptwXPoints *groupedFluxX; 134 ptwXPoints *groupedFluxX;
135 GIDI_settings_flux_order const *fluxOrder; 135 GIDI_settings_flux_order const *fluxOrder;
136 double const *energies, *fluxes; 136 double const *energies, *fluxes;
137 137
138 for( int order = 0; order < (int) flux.siz 138 for( int order = 0; order < (int) flux.size( ); ++order ) {
139 fluxOrder = flux[order]; 139 fluxOrder = flux[order];
140 energies = fluxOrder->getEnergies( ); 140 energies = fluxOrder->getEnergies( );
141 fluxes = fluxOrder->getFluxes( ); 141 fluxes = fluxOrder->getFluxes( );
142 if( ( fluxXY = ptwXY_createFrom_Xs_Ys( 142 if( ( fluxXY = ptwXY_createFrom_Xs_Ys( ptwXY_interpolationLinLin, NULL, 12, 1e-3, fluxOrder->size( ), 10,
143 fluxOrder->size( ), energies, flux 143 fluxOrder->size( ), energies, fluxes, &status_nf, 0 ) ) == NULL ) goto err;
144 mFluxXY.push_back( fluxXY ); 144 mFluxXY.push_back( fluxXY );
145 if( ( groupedFluxX = ptwXY_groupOneFun 145 if( ( groupedFluxX = ptwXY_groupOneFunction( fluxXY, groupX, ptwXY_group_normType_none, NULL, &status_nf ) ) == NULL ) goto err;
146 mGroupedFlux.push_back( groupedFluxX ) 146 mGroupedFlux.push_back( groupedFluxX );
147 } 147 }
148 return; 148 return;
149 149
150 err: 150 err:
151 throw 1; 151 throw 1;
152 } 152 }
153 /* 153 /*
154 ============================================== 154 =========================================================
155 */ 155 */
156 GIDI_settings_processedFlux::GIDI_settings_pro 156 GIDI_settings_processedFlux::GIDI_settings_processedFlux( GIDI_settings_processedFlux const &flux ) : mFlux( flux.mFlux ) {
157 157
158 nfu_status status_nf; 158 nfu_status status_nf;
159 ptwXYPoints *fluxXY; 159 ptwXYPoints *fluxXY;
160 ptwXPoints *fluxX; 160 ptwXPoints *fluxX;
161 161
162 for( int order = 0; order < (int) mFlux.si 162 for( int order = 0; order < (int) mFlux.size( ); ++order ) {
163 if( ( fluxXY = ptwXY_clone( flux.mFlux 163 if( ( fluxXY = ptwXY_clone( flux.mFluxXY[order], &status_nf ) ) == NULL ) goto err;
164 mFluxXY.push_back( fluxXY ); 164 mFluxXY.push_back( fluxXY );
165 if( ( fluxX = ptwX_clone( flux.mGroupe 165 if( ( fluxX = ptwX_clone( flux.mGroupedFlux[order], &status_nf ) ) == NULL ) goto err;
166 mGroupedFlux.push_back( fluxX ); 166 mGroupedFlux.push_back( fluxX );
167 } 167 }
168 return; 168 return;
169 169
170 err: 170 err:
171 for( std::vector<ptwXYPoints *>::iterator 171 for( std::vector<ptwXYPoints *>::iterator iter = mFluxXY.begin( ); iter != mFluxXY.end( ); ++iter ) ptwXY_free( *iter );
172 for( std::vector<ptwXPoints *>::iterator i 172 for( std::vector<ptwXPoints *>::iterator iter = mGroupedFlux.begin( ); iter != mGroupedFlux.end( ); ++iter ) ptwX_free( *iter );
173 throw 1; 173 throw 1;
174 } 174 }
175 /* 175 /*
176 ============================================== 176 =========================================================
177 */ 177 */
178 GIDI_settings_processedFlux& GIDI_settings_pro <<
179 if ( this != &flux ) { <<
180 // Clean-up existing things <<
181 for( std::vector<ptwXYPoints *>::iterator <<
182 for( std::vector<ptwXPoints *>::iterator i <<
183 // Now assign <<
184 mFlux = flux.mFlux; <<
185 nfu_status status_nf; <<
186 ptwXYPoints *fluxXY; <<
187 ptwXPoints *fluxX; <<
188 for( int order = 0; order < (int) mFlux.si <<
189 if( ( fluxXY = ptwXY_clone( flux.mFlux <<
190 mFluxXY.push_back( fluxXY ); <<
191 if( ( fluxX = ptwX_clone( flux.mGroupe <<
192 mGroupedFlux.push_back( fluxX ); <<
193 } <<
194 } <<
195 return *this; <<
196 <<
197 err: <<
198 for( std::vector<ptwXYPoints *>::iterator <<
199 for( std::vector<ptwXPoints *>::iterator i <<
200 throw 1; <<
201 } <<
202 /* <<
203 ============================================== <<
204 */ <<
205 GIDI_settings_processedFlux::~GIDI_settings_pr 178 GIDI_settings_processedFlux::~GIDI_settings_processedFlux( ) {
206 179
207 for( std::vector<ptwXYPoints *>::iterator 180 for( std::vector<ptwXYPoints *>::iterator iter = mFluxXY.begin( ); iter != mFluxXY.end( ); ++iter ) ptwXY_free( *iter );
208 for( std::vector<ptwXPoints *>::iterator i 181 for( std::vector<ptwXPoints *>::iterator iter = mGroupedFlux.begin( ); iter != mGroupedFlux.end( ); ++iter ) ptwX_free( *iter );
209 } 182 }
210 /* 183 /*
211 ============================================== 184 =========================================================
212 */ 185 */
213 ptwXPoints *GIDI_settings_processedFlux::group 186 ptwXPoints *GIDI_settings_processedFlux::groupFunction( statusMessageReporting * /*smr*/, ptwXPoints *groupX, ptwXYPoints *ptwXY1, int order ) const {
214 187
215 nfu_status status_nf; 188 nfu_status status_nf;
216 ptwXYPoints *fluxXY; 189 ptwXYPoints *fluxXY;
217 ptwXPoints *groupedX; 190 ptwXPoints *groupedX;
218 191
219 if( groupX == NULL ) return( NULL ); 192 if( groupX == NULL ) return( NULL );
220 if( order < 0 ) order = 0; 193 if( order < 0 ) order = 0;
221 if( order >= (int) mFluxXY.size( ) ) order 194 if( order >= (int) mFluxXY.size( ) ) order = (int) mFluxXY.size( ) - 1;
222 195
223 fluxXY = ptwXY_xSlice( mFluxXY[order], ptw 196 fluxXY = ptwXY_xSlice( mFluxXY[order], ptwXY_getXMin( ptwXY1 ), ptwXY_getXMax( ptwXY1 ), 10, 1, &status_nf );
224 // if( fluxXY == NULL ) smr_setReportError2 197 // if( fluxXY == NULL ) smr_setReportError2( smr, smr_unknownID, 1, "ptwXY_xSlice error %d: %s", status_nf, nfu_statusMessage( status_nf ) );
225 198
226 groupedX = ptwXY_groupTwoFunctions( ptwXY1 199 groupedX = ptwXY_groupTwoFunctions( ptwXY1, fluxXY, groupX, ptwXY_group_normType_norm, mGroupedFlux[order], &status_nf );
227 ptwXY_free( fluxXY ); 200 ptwXY_free( fluxXY );
228 // if( groupedX == NULL ) smr_setReportErro 201 // if( groupedX == NULL ) smr_setReportError2( smr, smr_unknownID, 1, "ptwXY_groupTwoFunctions error %d: %s", status_nf, nfu_statusMessage( status_nf ) );
229 return( groupedX ); 202 return( groupedX );
230 } 203 }
231 204