Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4GIDI_target.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 /*
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.
59 # <<END-copyright>>
60 */
61 
62 #include <iostream>
63 #include <stdlib.h>
64 
65 #include "G4GIDI_target.hh"
66 #include "G4GIDI_mass.hh"
67 #include "G4GIDI_Misc.hh"
68 
69 using namespace std;
70 using namespace GIDI;
71 
72 /*
73 ***************************************************************
74 */
75 G4GIDI_target::G4GIDI_target( const char *fileName ) {
76 
77  init( fileName );
78 }
79 /*
80 ***************************************************************
81 */
82 G4GIDI_target::G4GIDI_target( string &fileName ) {
83 
84  init( fileName.c_str( ) );
85 }
86 /*
87 ***************************************************************
88 */
89 void G4GIDI_target::init( const char *fileName ) {
90 
91  int i, j, n, *p, *q;
92  tpia_channel *channel;
93 
94  smr_initialize( &smr );
95  sourceFilename = fileName;
96  target = tpia_target_createRead( &smr, fileName );
97  if( !smr_isOk( &smr ) ) {
98  smr_print( &smr, stderr, 1 );
99  throw 1;
100  }
101  name = target->targetID->name;
102  mass = G4GIDI_targetMass( target->targetID->name );
103  equalProbableBinSampleMethod = "constant";
104  elasticIndices = NULL;
105  nElasticIndices = nCaptureIndices = nFissionIndices = nOthersIndices = 0;
106 
107  if( ( n = tpia_target_numberOfChannels( &smr, target ) ) > 0 ) {
108  p = elasticIndices = (int *) xData_malloc2( NULL, n * sizeof( double ), 1, "elasticIndices" );
109  for( i = 0; i < n; i++ ) { /* Find elastic channel(s). */
110  channel = tpia_target_heated_getChannelAtIndex( target->baseHeatedTarget, i );
111  if( channel->ENDL_C == 10 ) {
112  *(p++) = i;
113  nElasticIndices++;
114  }
115  }
116  captureIndices = p;
117  for( i = 0; i < n; i++ ) { /* Find capture channel(s). */
118  channel = tpia_target_heated_getChannelAtIndex( target->baseHeatedTarget, i );
119  if( channel->ENDL_C == 46 ) {
120  *(p++) = i;
121  nCaptureIndices++;
122  }
123  }
124 
125  fissionIndices = p;
126  for( i = 0; i < n; i++ ) { /* Find fission channel(s). */
127  channel = tpia_target_heated_getChannelAtIndex( target->baseHeatedTarget, i );
128  if( channel->fission != NULL ) {
129  *(p++) = i;
130  nFissionIndices++;
131  }
132  }
133  othersIndices = p;
134  for( i = 0; i < n; i++ ) { /* Find other channel(s). */
135  for( j = 0, q = elasticIndices; j < nElasticIndices; j++, q++ ) if( *q == i ) break;
136  if( j < nElasticIndices ) continue;
137  for( j = 0, q = captureIndices; j < nCaptureIndices; j++, q++ ) if( *q == i ) break;
138  if( j < nCaptureIndices ) continue;
139  for( j = 0, q = fissionIndices; j < nFissionIndices; j++, q++ ) if( *q == i ) break;
140  if( j < nFissionIndices ) continue;
141  *p = i;
142  p++;
143  nOthersIndices++;
144  }
145  }
146 }
147 /*
148 ***************************************************************
149 */
151 
152  tpia_target_free( &smr, target );
153  xData_free( &smr, elasticIndices );
154  smr_release( &smr );
155 }
156 /*
157 ***************************************************************
158 */
159 string *G4GIDI_target::getName( void ) { return( &name ); }
160 /*
161 ***************************************************************
162 */
163 string *G4GIDI_target::getFilename( void ) { return( &sourceFilename ); }
164 /*
165 ***************************************************************
166 */
167 int G4GIDI_target::getZ( void ) {
168 
169  return( target->targetID->Z );
170 }
171 /*
172 ***************************************************************
173 */
174 int G4GIDI_target::getA( void ) {
175 
176  return( target->targetID->A );
177 }
178 /*
179 ***************************************************************
180 */
181 int G4GIDI_target::getM( void ) {
182 
183  return( target->targetID->m );
184 }
185 /*
186 ***************************************************************
187 */
188 double G4GIDI_target::getMass( void ) {
189 
190  return( mass );
191 }
192 /*
193 ***************************************************************
194 */
195 int G4GIDI_target::getTemperatures( double *temperatures ) {
196 
197  return( tpia_target_getTemperatures( &smr, target, temperatures ) );
198 }
199 /*
200 ***************************************************************
201 */
203 
204  return( tpia_target_readHeatedTarget( &smr, target, index, 0 ) );
205 }
206 /*
207 ***************************************************************
208 */
210 
211  return( equalProbableBinSampleMethod );
212 }
213 /*
214 ***************************************************************
215 */
217 
218  if( method == "constant" ) {
219  equalProbableBinSampleMethod = "constant"; }
220  if( method == "linear" ) {
221  equalProbableBinSampleMethod = "linear"; }
222  else {
223  return( 1 );
224  }
225  return( 0 );
226 }
227 /*
228 ***************************************************************
229 */
231 
232  return( tpia_target_numberOfChannels( &smr, target ) );
233 }
234 /*
235 ***************************************************************
236 */
238 
240 }
241 /*
242 ***************************************************************
243 */
244 vector<channelID> *G4GIDI_target::getChannelIDs( void ) {
245 
246  return( getChannelIDs2( target->baseHeatedTarget->channels, tpia_target_numberOfChannels( &smr, target ) ) );
247 }
248 /*
249 ***************************************************************
250 */
251 vector<channelID> *G4GIDI_target::getProductionChannelIDs( void ) {
252 
253  return( getChannelIDs2( target->baseHeatedTarget->productionChannels, tpia_target_numberOfProductionChannels( &smr, target ) ) );
254 }
255 /*
256 ***************************************************************
257 */
258 vector<channelID> *G4GIDI_target::getChannelIDs2( tpia_channel **channels, int n ) {
259 
260  int i = 0;
261  vector<channelID> *listOfChannels;
262 
263  listOfChannels = new vector<channelID>( n );
264  for( i = 0; i < n; i++ ) (*listOfChannels)[i].ID = channels[i]->outputChannel;
265  return( listOfChannels );
266 }
267 /*
268 ***************************************************************
269 */
271 
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 );
282 }
283 /*
284 ***************************************************************
285 */
286 double G4GIDI_target::getTotalCrossSectionAtE( double e_in, double temperature ) {
287 
288  return( tpia_target_getTotalCrossSectionAtTAndE( NULL, target, temperature, -1, e_in, tpia_crossSectionType_pointwise ) );
289 }
290 /*
291 ***************************************************************
292 */
293 double G4GIDI_target::getElasticCrossSectionAtE( double e_in, double temperature ) {
294 
295  return( sumChannelCrossSectionAtE( nElasticIndices, elasticIndices, e_in, temperature ) );
296 }
297 /*
298 ***************************************************************
299 */
300 double G4GIDI_target::getCaptureCrossSectionAtE( double e_in, double temperature ) {
301 
302  return( sumChannelCrossSectionAtE( nCaptureIndices, captureIndices, e_in, temperature ) );
303 }
304 /*
305 ***************************************************************
306 */
307 double G4GIDI_target::getFissionCrossSectionAtE( double e_in, double temperature ) {
308 
309  return( sumChannelCrossSectionAtE( nFissionIndices, fissionIndices, e_in, temperature ) );
310 }
311 /*
312 ***************************************************************
313 */
314 double G4GIDI_target::getOthersCrossSectionAtE( double e_in, double temperature ) {
315 
316  return( sumChannelCrossSectionAtE( nOthersIndices, othersIndices, e_in, temperature ) );
317 }
318 /*
319 ***************************************************************
320 */
321 double G4GIDI_target::sumChannelCrossSectionAtE( int nIndices, int *indices, double e_in, double temperature ) {
322 
323  int i;
324  double xsec = 0.;
325 
326  for( i = 0; i < nIndices; i++ )
327  xsec += tpia_target_getIndexChannelCrossSectionAtE( &smr, target, indices[i], temperature, -1, e_in, tpia_crossSectionType_pointwise );
328  return( xsec );
329 }
330 /*
331 ***************************************************************
332 */
333 int G4GIDI_target::sampleChannelCrossSectionAtE( int nIndices, int *indices, double e_in, double temperature,
334  double (*rng)( void * ), void *rngState ) {
335 
336  int i;
337  double xsec = 0., rxsec = sumChannelCrossSectionAtE( nIndices, indices, e_in, temperature ) * tpia_misc_drng( rng, rngState );
338 
339  for( i = 0; i < nIndices - 1; i++ ) {
340  xsec += tpia_target_getIndexChannelCrossSectionAtE( &smr, target, indices[i], temperature, -1, e_in, tpia_crossSectionType_pointwise );
341  if( xsec >= rxsec ) break;
342  }
343  return( indices[i] );
344 }
345 /*
346 ***************************************************************
347 */
348 //double G4GIDI_target::getElasticFinalState( double e_in, double temperature, double (*rng)( void * ), void *rngState ) {
349 double G4GIDI_target::getElasticFinalState( double e_in, double , double (*rng)( void * ), void *rngState ) {
350 
351  tpia_decaySamplingInfo decaySamplingInfo;
352  tpia_channel *channel = tpia_target_heated_getChannelAtIndex_smr( &smr, target->baseHeatedTarget, elasticIndices[0] );
353  tpia_product *product;
354 
355  decaySamplingInfo.e_in = e_in;
356  decaySamplingInfo.isVelocity = 0;
357  tpia_frame_setColumn( &smr, &(decaySamplingInfo.frame), 0, tpia_referenceFrame_lab );
358  decaySamplingInfo.samplingMethods = &(target->samplingMethods);
359  decaySamplingInfo.rng = rng;
360  decaySamplingInfo.rngState = rngState;
361  product = tpia_decayChannel_getFirstProduct( &(channel->decayChannel) );
362  tpia_angular_SampleMu( &smr, &(product->angular), &decaySamplingInfo );
363 
364  return( decaySamplingInfo.mu );
365 }
366 /*
367 ***************************************************************
368 */
369 vector<G4GIDI_Product> *G4GIDI_target::getCaptureFinalState( double e_in, double temperature, double (*rng)( void * ), void *rngState ) {
370 
371  return( getFinalState( nCaptureIndices, captureIndices, e_in, temperature, rng, rngState ) );
372 }
373 /*
374 ***************************************************************
375 */
376 vector<G4GIDI_Product> *G4GIDI_target::getFissionFinalState( double e_in, double temperature, double (*rng)( void * ), void *rngState ) {
377 
378  return( getFinalState( nFissionIndices, fissionIndices, e_in, temperature, rng, rngState ) );
379 }
380 /*
381 ***************************************************************
382 */
383 vector<G4GIDI_Product> *G4GIDI_target::getOthersFinalState( double e_in, double temperature, double (*rng)( void * ), void *rngState ) {
384 
385  return( getFinalState( nOthersIndices, othersIndices, e_in, temperature, rng, rngState ) );
386 }
387 /*
388 ***************************************************************
389 */
390 vector<G4GIDI_Product> *G4GIDI_target::getFinalState( int nIndices, int *indices, double e_in, double temperature, double (*rng)( void * ), void *rngState ) {
391 
392 #define nProductsMax 50
393  int index = 0, i, n;
394  vector<G4GIDI_Product> *products = NULL;
395  tpia_decaySamplingInfo decaySamplingInfo;
396  tpia_productOutgoingData productDatas[nProductsMax], *productData;
397 
398  decaySamplingInfo.e_in = e_in;
399  decaySamplingInfo.samplingMethods = &(target->samplingMethods);
400  decaySamplingInfo.isVelocity = 0;
401  tpia_frame_setColumn( &smr, &(decaySamplingInfo.frame), 0, tpia_referenceFrame_lab );
402  decaySamplingInfo.rng = rng;
403  decaySamplingInfo.rngState = rngState;
404 
405  if( nIndices == 0 ) {
406  return( NULL ); }
407  else {
408  if( nIndices == 1 ) {
409  index = indices[0]; }
410  else {
411  index = sampleChannelCrossSectionAtE( nIndices, indices, e_in, temperature, rng, rngState );
412  }
413  }
414  n = tpia_target_heated_sampleIndexChannelProductsAtE( &smr, target->baseHeatedTarget, index, &decaySamplingInfo, nProductsMax, productDatas );
415  if( n > 0 ) {
416  if( ( products = new vector<G4GIDI_Product>( n ) ) != NULL ) {
417  for( i = 0; i < n; i++ ) {
418  productData = &(productDatas[i]);
419  (*products)[i].A = productData->productID->A;
420  (*products)[i].Z = productData->productID->Z;
421  (*products)[i].m = productData->productID->m;
422  (*products)[i].kineticEnergy = productData->kineticEnergy;
423  (*products)[i].px = productData->px_vx;
424  (*products)[i].py = productData->py_vy;
425  (*products)[i].pz = productData->pz_vz;
426  }
427  }
428  }
429 
430  return( products );
431 #undef nProductsMax
432 }