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