ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
GIDI_settings_particle.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file GIDI_settings_particle.cc
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 */
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 )
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 */
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 */
69 
70  if( mGroupX != NULL ) ptwX_free( mGroupX );
71 }
72 /*
73 =========================================================
74 */
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 */
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 */
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 */
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 */
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 */
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 */
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 }