Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/lend/src/GIDI_settings_particle.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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