| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 /* 26 /*
27 # <<BEGIN-copyright>> 27 # <<BEGIN-copyright>>
>> 28 # Copyright (c) 2010, Lawrence Livermore National Security, LLC.
>> 29 # Produced at the Lawrence Livermore National Laboratory
>> 30 # Written by Bret R. Beck, beck6@llnl.gov.
>> 31 # CODE-461393
>> 32 # All rights reserved.
>> 33 #
>> 34 # This file is part of GIDI. For details, see nuclear.llnl.gov.
>> 35 # Please also read the "Additional BSD Notice" at nuclear.llnl.gov.
>> 36 #
>> 37 # Redistribution and use in source and binary forms, with or without modification,
>> 38 # are permitted provided that the following conditions are met:
>> 39 #
>> 40 # 1) Redistributions of source code must retain the above copyright notice,
>> 41 # this list of conditions and the disclaimer below.
>> 42 # 2) Redistributions in binary form must reproduce the above copyright notice,
>> 43 # this list of conditions and the disclaimer (as noted below) in the
>> 44 # documentation and/or other materials provided with the distribution.
>> 45 # 3) Neither the name of the LLNS/LLNL nor the names of its contributors may be
>> 46 # used to endorse or promote products derived from this software without
>> 47 # specific prior written permission.
>> 48 #
>> 49 # THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
>> 50 # EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
>> 51 # OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT
>> 52 # SHALL LAWRENCE LIVERMORE NATIONAL SECURITY, LLC, THE U.S. DEPARTMENT OF ENERGY OR
>> 53 # CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
>> 54 # CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
>> 55 # OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED
>> 56 # AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
>> 57 # (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE,
>> 58 # EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
28 # <<END-copyright>> 59 # <<END-copyright>>
29 */ 60 */
>> 61
30 #include <iostream> 62 #include <iostream>
31 #include <stdlib.h> 63 #include <stdlib.h>
32 #include <cmath> <<
33 64
34 #include "G4GIDI_target.hh" 65 #include "G4GIDI_target.hh"
35 #include "G4GIDI_mass.hh" 66 #include "G4GIDI_mass.hh"
36 #include "G4GIDI_Misc.hh" 67 #include "G4GIDI_Misc.hh"
37 68
38 using namespace std; 69 using namespace std;
39 using namespace GIDI; 70 using namespace GIDI;
40 71
41 /* 72 /*
42 ********************************************** 73 ***************************************************************
43 */ 74 */
44 G4GIDI_target::G4GIDI_target( char const *file << 75 G4GIDI_target::G4GIDI_target( const char *fileName ) {
45 76
46 init( fileName ); 77 init( fileName );
47 } 78 }
48 /* 79 /*
49 ********************************************** 80 ***************************************************************
50 */ 81 */
51 G4GIDI_target::G4GIDI_target( string const &fi << 82 G4GIDI_target::G4GIDI_target( string &fileName ) {
52 83
53 init( fileName.c_str( ) ); 84 init( fileName.c_str( ) );
54 } 85 }
55 /* 86 /*
56 ********************************************** 87 ***************************************************************
57 */ 88 */
58 void G4GIDI_target::init( char const *fileName << 89 void G4GIDI_target::init( const char *fileName ) {
59 90
60 int i, j, n, *p, *q, ir; << 91 int i, j, n, *p, *q;
61 MCGIDI_reaction *reaction; << 92 tpia_channel *channel;
62 93
63 smr_initialize( &smr, smr_status_Ok, 1 ); << 94 smr_initialize( &smr );
64 sourceFilename = fileName; 95 sourceFilename = fileName;
65 target = MCGIDI_target_newRead( &smr, file << 96 target = tpia_target_createRead( &smr, fileName );
66 if( !smr_isOk( &smr ) ) { 97 if( !smr_isOk( &smr ) ) {
67 smr_print( &smr, 1 ); << 98 smr_print( &smr, stderr, 1 );
68 throw 1; 99 throw 1;
69 } 100 }
70 projectilesPOPID = target->projectilePOP-> << 101 name = target->targetID->name;
71 name = target->targetPOP->name; << 102 mass = G4GIDI_targetMass( target->targetID->name );
72 mass = G4GIDI_targetMass( target->targetPO <<
73 equalProbableBinSampleMethod = "constant"; 103 equalProbableBinSampleMethod = "constant";
74 elasticIndices = NULL; 104 elasticIndices = NULL;
75 nElasticIndices = nCaptureIndices = nFissi 105 nElasticIndices = nCaptureIndices = nFissionIndices = nOthersIndices = 0;
76 106
77 if( ( n = MCGIDI_target_numberOfReactions( << 107 if( ( n = tpia_target_numberOfChannels( &smr, target ) ) > 0 ) {
78 if( ( p = elasticIndices = (int *) smr << 108 p = elasticIndices = (int *) xData_malloc2( NULL, n * sizeof( double ), 1, "elasticIndices" );
79 smr_print( &smr, 1 ); <<
80 throw 1; <<
81 } <<
82 for( i = 0; i < n; i++ ) { /* Fin 109 for( i = 0; i < n; i++ ) { /* Find elastic channel(s). */
83 reaction = MCGIDI_target_heated_ge << 110 channel = tpia_target_heated_getChannelAtIndex( target->baseHeatedTarget, i );
84 if( MCGIDI_reaction_getENDF_MTNumb << 111 if( channel->ENDL_C == 10 ) {
85 *(p++) = i; 112 *(p++) = i;
86 nElasticIndices++; 113 nElasticIndices++;
87 } 114 }
88 } 115 }
89 captureIndices = p; 116 captureIndices = p;
90 for( i = 0; i < n; i++ ) { /* Fin 117 for( i = 0; i < n; i++ ) { /* Find capture channel(s). */
91 reaction = MCGIDI_target_heated_ge << 118 channel = tpia_target_heated_getChannelAtIndex( target->baseHeatedTarget, i );
92 if( MCGIDI_reaction_getENDF_MTNumb << 119 if( channel->ENDL_C == 46 ) {
93 *(p++) = i; 120 *(p++) = i;
94 nCaptureIndices++; 121 nCaptureIndices++;
95 } 122 }
96 } 123 }
97 124
98 fissionIndices = p; 125 fissionIndices = p;
99 for( i = 0; i < n; i++ ) { /* Fin 126 for( i = 0; i < n; i++ ) { /* Find fission channel(s). */
100 reaction = MCGIDI_target_heated_ge << 127 channel = tpia_target_heated_getChannelAtIndex( target->baseHeatedTarget, i );
101 ir = MCGIDI_reaction_getENDF_MTNum << 128 if( channel->fission != NULL ) {
102 if( ( ir != 18 ) && ( ir != 19 ) & << 129 *(p++) = i;
103 *(p++) = i; << 130 nFissionIndices++;
104 nFissionIndices++; << 131 }
105 } 132 }
106 othersIndices = p; 133 othersIndices = p;
107 for( i = 0; i < n; i++ ) { /* Fin 134 for( i = 0; i < n; i++ ) { /* Find other channel(s). */
108 for( j = 0, q = elasticIndices; j 135 for( j = 0, q = elasticIndices; j < nElasticIndices; j++, q++ ) if( *q == i ) break;
109 if( j < nElasticIndices ) continue 136 if( j < nElasticIndices ) continue;
110 for( j = 0, q = captureIndices; j 137 for( j = 0, q = captureIndices; j < nCaptureIndices; j++, q++ ) if( *q == i ) break;
111 if( j < nCaptureIndices ) continue 138 if( j < nCaptureIndices ) continue;
112 for( j = 0, q = fissionIndices; j 139 for( j = 0, q = fissionIndices; j < nFissionIndices; j++, q++ ) if( *q == i ) break;
113 if( j < nFissionIndices ) continue 140 if( j < nFissionIndices ) continue;
114 *p = i; 141 *p = i;
115 p++; 142 p++;
116 nOthersIndices++; 143 nOthersIndices++;
117 } 144 }
118 #if 0 <<
119 printf( "elastic %d: ", nElasticIndices ); <<
120 for( i = 0; i < nElasticIndices; i++ ) printf( <<
121 printf( "\ncapture %d: ", nCaptureIndices ); <<
122 for( i = 0; i < nCaptureIndices; i++ ) printf( <<
123 printf( "\nfission %d: ", nFissionIndices ); <<
124 for( i = 0; i < nFissionIndices; i++ ) printf( <<
125 printf( "\nothers %d: ", nOthersIndices ); <<
126 for( i = 0; i < nOthersIndices; i++ ) printf( <<
127 printf( "\n" ); <<
128 #endif <<
129 } 145 }
130 } 146 }
131 /* 147 /*
132 ********************************************** 148 ***************************************************************
133 */ 149 */
134 G4GIDI_target::~G4GIDI_target( ) { 150 G4GIDI_target::~G4GIDI_target( ) {
135 151
136 MCGIDI_target_free( &smr, target ); << 152 tpia_target_free( &smr, target );
137 smr_freeMemory( (void **) &elasticIndices << 153 xData_free( &smr, elasticIndices );
138 smr_release( &smr ); 154 smr_release( &smr );
139 } 155 }
140 /* 156 /*
141 ********************************************** 157 ***************************************************************
142 */ 158 */
143 string *G4GIDI_target::getName( void ) { retur 159 string *G4GIDI_target::getName( void ) { return( &name ); }
144 /* 160 /*
145 ********************************************** 161 ***************************************************************
146 */ 162 */
147 string *G4GIDI_target::getFilename( void ) { r 163 string *G4GIDI_target::getFilename( void ) { return( &sourceFilename ); }
148 /* 164 /*
149 ********************************************** 165 ***************************************************************
150 */ 166 */
151 int G4GIDI_target::getZ( void ) { 167 int G4GIDI_target::getZ( void ) {
152 168
153 return( target->targetPOP->Z ); << 169 return( target->targetID->Z );
154 } 170 }
155 /* 171 /*
156 ********************************************** 172 ***************************************************************
157 */ 173 */
158 int G4GIDI_target::getA( void ) { 174 int G4GIDI_target::getA( void ) {
159 175
160 return( target->targetPOP->A ); << 176 return( target->targetID->A );
161 } 177 }
162 /* 178 /*
163 ********************************************** 179 ***************************************************************
164 */ 180 */
165 int G4GIDI_target::getM( void ) { 181 int G4GIDI_target::getM( void ) {
166 182
167 return( target->targetPOP->m ); << 183 return( target->targetID->m );
168 } 184 }
169 /* 185 /*
170 ********************************************** 186 ***************************************************************
171 */ 187 */
172 double G4GIDI_target::getMass( void ) { 188 double G4GIDI_target::getMass( void ) {
173 189
174 return( mass ); 190 return( mass );
175 } 191 }
176 /* 192 /*
177 ********************************************** 193 ***************************************************************
178 */ 194 */
179 int G4GIDI_target::getTemperatures( double *te 195 int G4GIDI_target::getTemperatures( double *temperatures ) {
180 196
181 return( MCGIDI_target_getTemperatures( &sm << 197 return( tpia_target_getTemperatures( &smr, target, temperatures ) );
182 } 198 }
183 /* 199 /*
184 ********************************************** 200 ***************************************************************
185 */ 201 */
186 int G4GIDI_target::readTemperature( int index 202 int G4GIDI_target::readTemperature( int index ) {
187 203
188 return( MCGIDI_target_readHeatedTarget( &s << 204 return( tpia_target_readHeatedTarget( &smr, target, index, 0 ) );
189 } 205 }
190 /* 206 /*
191 ********************************************** 207 ***************************************************************
192 */ 208 */
193 string G4GIDI_target::getEqualProbableBinSampl 209 string G4GIDI_target::getEqualProbableBinSampleMethod( void ) {
194 210
195 return( equalProbableBinSampleMethod ); 211 return( equalProbableBinSampleMethod );
196 } 212 }
197 /* 213 /*
198 ********************************************** 214 ***************************************************************
199 */ 215 */
200 int G4GIDI_target::setEqualProbableBinSampleMe 216 int G4GIDI_target::setEqualProbableBinSampleMethod( string method ) {
201 217
202 if( method == "constant" ) { 218 if( method == "constant" ) {
203 equalProbableBinSampleMethod = "consta 219 equalProbableBinSampleMethod = "constant"; }
204 if( method == "linear" ) { 220 if( method == "linear" ) {
205 equalProbableBinSampleMethod = "linear 221 equalProbableBinSampleMethod = "linear"; }
206 else { 222 else {
207 return( 1 ); 223 return( 1 );
208 } 224 }
209 return( 0 ); 225 return( 0 );
210 } 226 }
211 /* 227 /*
212 ********************************************** 228 ***************************************************************
213 */ 229 */
214 int G4GIDI_target::getNumberOfChannels( void ) 230 int G4GIDI_target::getNumberOfChannels( void ) {
215 231
216 return( MCGIDI_target_numberOfReactions( & << 232 return( tpia_target_numberOfChannels( &smr, target ) );
217 } 233 }
218 /* 234 /*
219 ********************************************** 235 ***************************************************************
220 */ 236 */
221 int G4GIDI_target::getNumberOfProductionChanne 237 int G4GIDI_target::getNumberOfProductionChannels( void ) {
222 238
223 return( MCGIDI_target_numberOfProductionRe << 239 return( tpia_target_numberOfProductionChannels( &smr, target ) );
224 } 240 }
225 /* 241 /*
226 ********************************************** 242 ***************************************************************
227 */ 243 */
228 channelID G4GIDI_target::getChannelsID( int ch << 244 vector<channelID> *G4GIDI_target::getChannelIDs( void ) {
229 245
230 MCGIDI_reaction *reaction; << 246 return( getChannelIDs2( target->baseHeatedTarget->channels, tpia_target_numberOfChannels( &smr, target ) ) );
>> 247 }
>> 248 /*
>> 249 ***************************************************************
>> 250 */
>> 251 vector<channelID> *G4GIDI_target::getProductionChannelIDs( void ) {
231 252
232 if( ( reaction = MCGIDI_target_heated_getR << 253 return( getChannelIDs2( target->baseHeatedTarget->productionChannels, tpia_target_numberOfProductionChannels( &smr, target ) ) );
233 smr_print( &smr, 1 ); <<
234 throw 1; <<
235 } <<
236 return( string( reaction->outputChannelStr <<
237 } 254 }
238 /* 255 /*
239 ********************************************** 256 ***************************************************************
240 */ 257 */
241 vector<channelID> *G4GIDI_target::getChannelID << 258 vector<channelID> *G4GIDI_target::getChannelIDs2( tpia_channel **channels, int n ) {
242 259
243 int i, n = MCGIDI_target_numberOfReactions << 260 int i = 0;
244 MCGIDI_reaction *reaction; <<
245 vector<channelID> *listOfChannels; 261 vector<channelID> *listOfChannels;
246 262
247 listOfChannels = new vector<channelID>( n 263 listOfChannels = new vector<channelID>( n );
248 for( i = 0; i < n; i++ ) { << 264 for( i = 0; i < n; i++ ) (*listOfChannels)[i].ID = channels[i]->outputChannel;
249 reaction = MCGIDI_target_heated_getRea <<
250 (*listOfChannels)[i] = reaction->outpu <<
251 } <<
252 return( listOfChannels ); 265 return( listOfChannels );
253 } 266 }
254 /* 267 /*
255 ********************************************** 268 ***************************************************************
256 */ 269 */
257 vector<channelID> *G4GIDI_target::getProductio << 270 vector<double> *G4GIDI_target::getEnergyGridAtTIndex( int index ) {
258 271
259 return( NULL ); << 272 xData_Int i, n;
>> 273 double *dEnergyGrid;
>> 274 vector<double> *energyGrid;
>> 275 vector<double>::iterator iter;
>> 276
>> 277 n = tpia_target_getEnergyGridAtTIndex( &smr, target, index, &dEnergyGrid );
>> 278 if( n < 0 ) return( NULL );
>> 279 energyGrid = new vector<double>( n );
>> 280 for( i = 0, iter = energyGrid->begin( ); i < n; i++, iter++ ) *iter = dEnergyGrid[i];
>> 281 return( energyGrid );
260 } 282 }
261 /* 283 /*
262 ********************************************** 284 ***************************************************************
263 */ 285 */
264 double G4GIDI_target::getTotalCrossSectionAtE( 286 double G4GIDI_target::getTotalCrossSectionAtE( double e_in, double temperature ) {
265 287
266 MCGIDI_quantitiesLookupModes mode( project << 288 return( tpia_target_getTotalCrossSectionAtTAndE( NULL, target, temperature, -1, e_in, tpia_crossSectionType_pointwise ) );
267 <<
268 mode.setProjectileEnergy( e_in ); <<
269 mode.setCrossSectionMode( MCGIDI_quantityL <<
270 mode.setTemperature( temperature ); <<
271 <<
272 return( MCGIDI_target_getTotalCrossSection <<
273 } 289 }
274 /* 290 /*
275 ********************************************** 291 ***************************************************************
276 */ 292 */
277 double G4GIDI_target::getElasticCrossSectionAt 293 double G4GIDI_target::getElasticCrossSectionAtE( double e_in, double temperature ) {
278 294
279 return( sumChannelCrossSectionAtE( nElasti 295 return( sumChannelCrossSectionAtE( nElasticIndices, elasticIndices, e_in, temperature ) );
280 } 296 }
281 /* 297 /*
282 ********************************************** 298 ***************************************************************
283 */ 299 */
284 double G4GIDI_target::getCaptureCrossSectionAt 300 double G4GIDI_target::getCaptureCrossSectionAtE( double e_in, double temperature ) {
285 301
286 return( sumChannelCrossSectionAtE( nCaptur 302 return( sumChannelCrossSectionAtE( nCaptureIndices, captureIndices, e_in, temperature ) );
287 } 303 }
288 /* 304 /*
289 ********************************************** 305 ***************************************************************
290 */ 306 */
291 double G4GIDI_target::getFissionCrossSectionAt 307 double G4GIDI_target::getFissionCrossSectionAtE( double e_in, double temperature ) {
292 308
293 return( sumChannelCrossSectionAtE( nFissio 309 return( sumChannelCrossSectionAtE( nFissionIndices, fissionIndices, e_in, temperature ) );
294 } 310 }
295 /* 311 /*
296 ********************************************** 312 ***************************************************************
297 */ 313 */
298 double G4GIDI_target::getOthersCrossSectionAtE 314 double G4GIDI_target::getOthersCrossSectionAtE( double e_in, double temperature ) {
299 315
300 return( sumChannelCrossSectionAtE( nOthers 316 return( sumChannelCrossSectionAtE( nOthersIndices, othersIndices, e_in, temperature ) );
301 } 317 }
302 /* 318 /*
303 ********************************************** 319 ***************************************************************
304 */ 320 */
305 double G4GIDI_target::sumChannelCrossSectionAt 321 double G4GIDI_target::sumChannelCrossSectionAtE( int nIndices, int *indices, double e_in, double temperature ) {
306 322
307 int i; 323 int i;
308 double xsec = 0.; 324 double xsec = 0.;
309 MCGIDI_quantitiesLookupModes mode( project <<
310 <<
311 mode.setProjectileEnergy( e_in ); <<
312 mode.setCrossSectionMode( MCGIDI_quantityL <<
313 mode.setTemperature( temperature ); <<
314 325
315 for( i = 0; i < nIndices; i++ ) 326 for( i = 0; i < nIndices; i++ )
316 xsec += MCGIDI_target_getIndexReaction << 327 xsec += tpia_target_getIndexChannelCrossSectionAtE( &smr, target, indices[i], temperature, -1, e_in, tpia_crossSectionType_pointwise );
317 return( xsec ); 328 return( xsec );
318 } 329 }
319 /* 330 /*
320 ********************************************** 331 ***************************************************************
321 */ 332 */
322 int G4GIDI_target::sampleChannelCrossSectionAt 333 int G4GIDI_target::sampleChannelCrossSectionAtE( int nIndices, int *indices, double e_in, double temperature,
323 double (*rng)( void * ), void *rngStat 334 double (*rng)( void * ), void *rngState ) {
324 335
325 int i; 336 int i;
326 double xsec = 0., rxsec = sumChannelCrossS << 337 double xsec = 0., rxsec = sumChannelCrossSectionAtE( nIndices, indices, e_in, temperature ) * tpia_misc_drng( rng, rngState );
327 MCGIDI_quantitiesLookupModes mode( project <<
328 <<
329 mode.setProjectileEnergy( e_in ); <<
330 mode.setCrossSectionMode( MCGIDI_quantityL <<
331 mode.setTemperature( temperature ); <<
332 338
333 for( i = 0; i < nIndices - 1; i++ ) { 339 for( i = 0; i < nIndices - 1; i++ ) {
334 xsec += MCGIDI_target_getIndexReaction << 340 xsec += tpia_target_getIndexChannelCrossSectionAtE( &smr, target, indices[i], temperature, -1, e_in, tpia_crossSectionType_pointwise );
335 if( xsec >= rxsec ) break; 341 if( xsec >= rxsec ) break;
336 } 342 }
337 return( indices[i] ); 343 return( indices[i] );
338 } 344 }
339 /* 345 /*
340 ********************************************** 346 ***************************************************************
341 */ 347 */
342 double G4GIDI_target::getElasticFinalState( do << 348 //double G4GIDI_target::getElasticFinalState( double e_in, double temperature, double (*rng)( void * ), void *rngState ) {
343 << 349 double G4GIDI_target::getElasticFinalState( double e_in, double , double (*rng)( void * ), void *rngState ) {
344 MCGIDI_decaySamplingInfo decaySamplingInfo <<
345 MCGIDI_reaction *reaction = MCGIDI_target_ <<
346 MCGIDI_product *product; <<
347 MCGIDI_quantitiesLookupModes mode( project <<
348 <<
349 if( ( product = MCGIDI_outputChannel_getPr <<
350 smr_print( &smr, 1 ); <<
351 throw 1; <<
352 } <<
353 350
354 mode.setProjectileEnergy( e_in ); << 351 tpia_decaySamplingInfo decaySamplingInfo;
355 mode.setCrossSectionMode( MCGIDI_quantityL << 352 tpia_channel *channel = tpia_target_heated_getChannelAtIndex_smr( &smr, target->baseHeatedTarget, elasticIndices[0] );
356 mode.setTemperature( temperature ); << 353 tpia_product *product;
357 354
>> 355 decaySamplingInfo.e_in = e_in;
358 decaySamplingInfo.isVelocity = 0; 356 decaySamplingInfo.isVelocity = 0;
>> 357 tpia_frame_setColumn( &smr, &(decaySamplingInfo.frame), 0, tpia_referenceFrame_lab );
>> 358 decaySamplingInfo.samplingMethods = &(target->samplingMethods);
359 decaySamplingInfo.rng = rng; 359 decaySamplingInfo.rng = rng;
360 decaySamplingInfo.rngState = rngState; 360 decaySamplingInfo.rngState = rngState;
361 if( MCGIDI_product_sampleMu( &smr, product << 361 product = tpia_decayChannel_getFirstProduct( &(channel->decayChannel) );
362 smr_print( &smr, 1 ); << 362 tpia_angular_SampleMu( &smr, &(product->angular), &decaySamplingInfo );
363 throw 1; <<
364 } <<
365 363
366 return( decaySamplingInfo.mu ); 364 return( decaySamplingInfo.mu );
367 } 365 }
368 /* 366 /*
369 ********************************************** 367 ***************************************************************
370 */ 368 */
371 vector<G4GIDI_Product> *G4GIDI_target::getCapt 369 vector<G4GIDI_Product> *G4GIDI_target::getCaptureFinalState( double e_in, double temperature, double (*rng)( void * ), void *rngState ) {
372 370
373 return( getFinalState( nCaptureIndices, ca 371 return( getFinalState( nCaptureIndices, captureIndices, e_in, temperature, rng, rngState ) );
374 } 372 }
375 /* 373 /*
376 ********************************************** 374 ***************************************************************
377 */ 375 */
378 vector<G4GIDI_Product> *G4GIDI_target::getFiss 376 vector<G4GIDI_Product> *G4GIDI_target::getFissionFinalState( double e_in, double temperature, double (*rng)( void * ), void *rngState ) {
379 377
380 return( getFinalState( nFissionIndices, fi 378 return( getFinalState( nFissionIndices, fissionIndices, e_in, temperature, rng, rngState ) );
381 } 379 }
382 /* 380 /*
383 ********************************************** 381 ***************************************************************
384 */ 382 */
385 vector<G4GIDI_Product> *G4GIDI_target::getOthe 383 vector<G4GIDI_Product> *G4GIDI_target::getOthersFinalState( double e_in, double temperature, double (*rng)( void * ), void *rngState ) {
386 384
387 return( getFinalState( nOthersIndices, oth 385 return( getFinalState( nOthersIndices, othersIndices, e_in, temperature, rng, rngState ) );
388 } 386 }
389 /* 387 /*
390 ********************************************** 388 ***************************************************************
391 */ 389 */
392 vector<G4GIDI_Product> *G4GIDI_target::getFina << 390 vector<G4GIDI_Product> *G4GIDI_target::getFinalState( int nIndices, int *indices, double e_in, double temperature, double (*rng)( void * ), void *rngState ) {
393 double (*rng)( void * ), void *rngState ) <<
394 391
>> 392 #define nProductsMax 50
395 int index = 0, i, n; 393 int index = 0, i, n;
396 vector<G4GIDI_Product> *products = NULL; 394 vector<G4GIDI_Product> *products = NULL;
397 MCGIDI_decaySamplingInfo decaySamplingInfo << 395 tpia_decaySamplingInfo decaySamplingInfo;
398 MCGIDI_sampledProductsDatas sampledProduct << 396 tpia_productOutgoingData productDatas[nProductsMax], *productData;
399 MCGIDI_sampledProductsData *productData; <<
400 MCGIDI_quantitiesLookupModes mode( project <<
401 397
>> 398 decaySamplingInfo.e_in = e_in;
>> 399 decaySamplingInfo.samplingMethods = &(target->samplingMethods);
402 decaySamplingInfo.isVelocity = 0; 400 decaySamplingInfo.isVelocity = 0;
>> 401 tpia_frame_setColumn( &smr, &(decaySamplingInfo.frame), 0, tpia_referenceFrame_lab );
403 decaySamplingInfo.rng = rng; 402 decaySamplingInfo.rng = rng;
404 decaySamplingInfo.rngState = rngState; 403 decaySamplingInfo.rngState = rngState;
405 404
406 if( nIndices == 0 ) { 405 if( nIndices == 0 ) {
407 return( NULL ); } 406 return( NULL ); }
408 else { 407 else {
409 if( nIndices == 1 ) { 408 if( nIndices == 1 ) {
410 index = indices[0]; } 409 index = indices[0]; }
411 else { 410 else {
412 index = sampleChannelCrossSectionA 411 index = sampleChannelCrossSectionAtE( nIndices, indices, e_in, temperature, rng, rngState );
413 } 412 }
414 } 413 }
415 << 414 n = tpia_target_heated_sampleIndexChannelProductsAtE( &smr, target->baseHeatedTarget, index, &decaySamplingInfo, nProductsMax, productDatas );
416 MCGIDI_sampledProducts_initialize( &smr, & <<
417 if( !smr_isOk( &smr ) ) { <<
418 smr_print( &smr, 1 ); <<
419 throw 1; <<
420 } <<
421 <<
422 mode.setProjectileEnergy( e_in ); <<
423 mode.setCrossSectionMode( MCGIDI_quantityL <<
424 mode.setTemperature( temperature ); <<
425 <<
426 n = MCGIDI_target_heated_sampleIndexReacti <<
427 &decaySamplingInfo, &sampledProduc <<
428 if( !smr_isOk( &smr ) ) { <<
429 smr_print( &smr, 1 ); <<
430 throw 1; <<
431 } <<
432 if( n > 0 ) { 415 if( n > 0 ) {
433 if( ( products = new vector<G4GIDI_Pro 416 if( ( products = new vector<G4GIDI_Product>( n ) ) != NULL ) {
434 for( i = 0; i < n; i++ ) { 417 for( i = 0; i < n; i++ ) {
435 productData = &(sampledProduct << 418 productData = &(productDatas[i]);
436 (*products)[i].A = productData << 419 (*products)[i].A = productData->productID->A;
437 (*products)[i].Z = productData << 420 (*products)[i].Z = productData->productID->Z;
438 (*products)[i].m = productData << 421 (*products)[i].m = productData->productID->m;
439 (*products)[i].kineticEnergy = 422 (*products)[i].kineticEnergy = productData->kineticEnergy;
440 (*products)[i].px = productDat 423 (*products)[i].px = productData->px_vx;
441 (*products)[i].py = productDat 424 (*products)[i].py = productData->py_vy;
442 (*products)[i].pz = productDat 425 (*products)[i].pz = productData->pz_vz;
443 (*products)[i].birthTimeSec = produc <<
444 } 426 }
445 } 427 }
446 } 428 }
447 MCGIDI_sampledProducts_release( &smr, &sam <<
448 429
449 return( products ); 430 return( products );
450 } << 431 #undef nProductsMax
451 /* <<
452 ********************************************** <<
453 */ <<
454 double G4GIDI_target::getReactionsThreshold( i <<
455 <<
456 return( MCGIDI_target_heated_getReactionsT <<
457 } <<
458 /* <<
459 ********************************************** <<
460 */ <<
461 double G4GIDI_target::getReactionsDomain( int <<
462 <<
463 return( MCGIDI_target_heated_getReactionsD <<
464 } 432 }
465 433