G4GIDI_target.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 /*
00027 # <<BEGIN-copyright>>
00028 # Copyright (c) 2010, Lawrence Livermore National Security, LLC. 
00029 # Produced at the Lawrence Livermore National Laboratory 
00030 # Written by Bret R. Beck, beck6@llnl.gov. 
00031 # CODE-461393
00032 # All rights reserved. 
00033 #  
00034 # This file is part of GIDI. For details, see nuclear.llnl.gov. 
00035 # Please also read the "Additional BSD Notice" at nuclear.llnl.gov. 
00036 # 
00037 # Redistribution and use in source and binary forms, with or without modification, 
00038 # are permitted provided that the following conditions are met: 
00039 #
00040 #      1) Redistributions of source code must retain the above copyright notice, 
00041 #         this list of conditions and the disclaimer below.
00042 #      2) Redistributions in binary form must reproduce the above copyright notice, 
00043 #         this list of conditions and the disclaimer (as noted below) in the 
00044 #          documentation and/or other materials provided with the distribution.
00045 #      3) Neither the name of the LLNS/LLNL nor the names of its contributors may be 
00046 #         used to endorse or promote products derived from this software without 
00047 #         specific prior written permission. 
00048 #
00049 # THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY 
00050 # EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES 
00051 # OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT 
00052 # SHALL LAWRENCE LIVERMORE NATIONAL SECURITY, LLC, THE U.S. DEPARTMENT OF ENERGY OR 
00053 # CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 
00054 # CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS 
00055 # OR SERVICES;  LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED 
00056 # AND ON  ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 
00057 # (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, 
00058 # EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 
00059 # <<END-copyright>>
00060 */
00061 
00062 #include <iostream>
00063 #include <stdlib.h>
00064 
00065 #include "G4GIDI_target.hh"
00066 #include "G4GIDI_mass.hh"
00067 #include "G4GIDI_Misc.hh"
00068 
00069 using namespace std;
00070 using namespace GIDI;
00071 
00072 /*
00073 ***************************************************************
00074 */
00075 G4GIDI_target::G4GIDI_target( const char *fileName ) {
00076 
00077     init( fileName );
00078 }
00079 /*
00080 ***************************************************************
00081 */
00082 G4GIDI_target::G4GIDI_target( string &fileName ) {
00083 
00084     init( fileName.c_str( ) );
00085 }
00086 /*
00087 ***************************************************************
00088 */
00089 void G4GIDI_target::init( const char *fileName ) {
00090 
00091     int i, j, n, *p, *q;
00092     tpia_channel *channel;
00093 
00094     smr_initialize( &smr );
00095     sourceFilename = fileName;
00096     target = tpia_target_createRead( &smr, fileName );
00097     if( !smr_isOk( &smr ) ) {
00098         smr_print( &smr, stderr, 1 );
00099         throw 1;
00100     }
00101     name = target->targetID->name;
00102     mass = G4GIDI_targetMass( target->targetID->name );
00103     equalProbableBinSampleMethod = "constant";
00104     elasticIndices = NULL;
00105     nElasticIndices = nCaptureIndices = nFissionIndices = nOthersIndices = 0;
00106 
00107     if( ( n = tpia_target_numberOfChannels( &smr, target ) ) > 0 ) {
00108         p = elasticIndices = (int *) xData_malloc2( NULL, n * sizeof( double ), 1, "elasticIndices" );
00109         for( i = 0; i < n; i++ ) {      /* Find elastic channel(s). */
00110             channel = tpia_target_heated_getChannelAtIndex( target->baseHeatedTarget, i );
00111             if( channel->ENDL_C == 10 ) {
00112                 *(p++) = i;
00113                 nElasticIndices++;
00114             }
00115         }
00116         captureIndices = p;
00117         for( i = 0; i < n; i++ ) {      /* Find capture channel(s). */
00118             channel = tpia_target_heated_getChannelAtIndex( target->baseHeatedTarget, i );
00119             if( channel->ENDL_C == 46 ) {
00120                 *(p++) = i;
00121                 nCaptureIndices++;
00122             }
00123         }
00124 
00125         fissionIndices = p;
00126         for( i = 0; i < n; i++ ) {      /* Find fission channel(s). */
00127             channel = tpia_target_heated_getChannelAtIndex( target->baseHeatedTarget, i );
00128             if( channel->fission != NULL ) {
00129                 *(p++) = i;
00130                 nFissionIndices++;
00131             }
00132         }
00133         othersIndices = p;
00134         for( i = 0; i < n; i++ ) {      /* Find other channel(s). */
00135             for( j = 0, q = elasticIndices; j < nElasticIndices; j++, q++ ) if( *q == i ) break;
00136             if( j < nElasticIndices ) continue;
00137             for( j = 0, q = captureIndices; j < nCaptureIndices; j++, q++ ) if( *q == i ) break;
00138             if( j < nCaptureIndices ) continue;
00139             for( j = 0, q = fissionIndices; j < nFissionIndices; j++, q++ ) if( *q == i ) break;
00140             if( j < nFissionIndices ) continue;
00141             *p = i;
00142             p++;
00143             nOthersIndices++;
00144         }
00145     }
00146 }
00147 /*
00148 ***************************************************************
00149 */
00150 G4GIDI_target::~G4GIDI_target( ) {
00151 
00152     tpia_target_free( &smr, target );
00153     xData_free( &smr, elasticIndices );
00154     smr_release( &smr );
00155 }
00156 /*
00157 ***************************************************************
00158 */
00159 string *G4GIDI_target::getName( void ) { return( &name ); }
00160 /*
00161 ***************************************************************
00162 */
00163 string *G4GIDI_target::getFilename( void ) { return( &sourceFilename ); }
00164 /*
00165 ***************************************************************
00166 */
00167 int G4GIDI_target::getZ( void ) {
00168    
00169     return( target->targetID->Z );
00170 }
00171 /*
00172 ***************************************************************
00173 */
00174 int G4GIDI_target::getA( void ) {
00175    
00176     return( target->targetID->A );
00177 }
00178 /*
00179 ***************************************************************
00180 */
00181 int G4GIDI_target::getM( void ) {
00182    
00183     return( target->targetID->m );
00184 }
00185 /*
00186 ***************************************************************
00187 */
00188 double G4GIDI_target::getMass( void ) {
00189 
00190     return( mass );
00191 }
00192 /*
00193 ***************************************************************
00194 */
00195 int G4GIDI_target::getTemperatures( double *temperatures ) {
00196 
00197     return( tpia_target_getTemperatures( &smr, target, temperatures ) );
00198 }
00199 /*
00200 ***************************************************************
00201 */
00202 int G4GIDI_target::readTemperature( int index ) {
00203 
00204     return( tpia_target_readHeatedTarget( &smr, target, index, 0 ) );
00205 }
00206 /*
00207 ***************************************************************
00208 */
00209 string G4GIDI_target::getEqualProbableBinSampleMethod( void ) {
00210 
00211     return( equalProbableBinSampleMethod );
00212 }
00213 /*
00214 ***************************************************************
00215 */
00216 int G4GIDI_target::setEqualProbableBinSampleMethod( string method ) {
00217 
00218     if( method == "constant" ) {
00219         equalProbableBinSampleMethod = "constant"; }
00220     if( method == "linear" ) {
00221         equalProbableBinSampleMethod = "linear"; }
00222     else {
00223         return( 1 );
00224     }
00225     return( 0 );
00226 }
00227 /*
00228 ***************************************************************
00229 */
00230 int G4GIDI_target::getNumberOfChannels( void ) {
00231 
00232     return( tpia_target_numberOfChannels( &smr, target ) );
00233 }
00234 /*
00235 ***************************************************************
00236 */
00237 int G4GIDI_target::getNumberOfProductionChannels( void ) {
00238 
00239     return( tpia_target_numberOfProductionChannels( &smr, target ) );
00240 }
00241 /*
00242 ***************************************************************
00243 */
00244 vector<channelID> *G4GIDI_target::getChannelIDs( void ) { 
00245 
00246     return( getChannelIDs2( target->baseHeatedTarget->channels, tpia_target_numberOfChannels( &smr, target ) ) );
00247 }
00248 /*
00249 ***************************************************************
00250 */
00251 vector<channelID> *G4GIDI_target::getProductionChannelIDs( void ) {
00252 
00253     return( getChannelIDs2( target->baseHeatedTarget->productionChannels, tpia_target_numberOfProductionChannels( &smr, target ) ) );
00254 }
00255 /*
00256 ***************************************************************
00257 */
00258 vector<channelID> *G4GIDI_target::getChannelIDs2( tpia_channel **channels, int n ) {
00259 
00260     int i = 0;
00261     vector<channelID> *listOfChannels;
00262 
00263     listOfChannels = new vector<channelID>( n );
00264     for( i = 0; i < n; i++ ) (*listOfChannels)[i].ID = channels[i]->outputChannel;
00265     return( listOfChannels );
00266 }
00267 /*
00268 ***************************************************************
00269 */
00270 vector<double> *G4GIDI_target::getEnergyGridAtTIndex( int index ) {
00271 
00272     xData_Int i, n;
00273     double *dEnergyGrid;
00274     vector<double> *energyGrid;
00275     vector<double>::iterator iter;
00276 
00277     n = tpia_target_getEnergyGridAtTIndex( &smr, target, index, &dEnergyGrid );
00278     if( n < 0 ) return( NULL );
00279     energyGrid = new vector<double>( n );
00280     for( i = 0, iter = energyGrid->begin( ); i < n; i++, iter++ ) *iter = dEnergyGrid[i];
00281     return( energyGrid );
00282 }
00283 /*
00284 ***************************************************************
00285 */
00286 double G4GIDI_target::getTotalCrossSectionAtE( double e_in, double temperature ) {
00287 
00288     return( tpia_target_getTotalCrossSectionAtTAndE( NULL, target, temperature, -1, e_in, tpia_crossSectionType_pointwise ) );
00289 }
00290 /*
00291 ***************************************************************
00292 */
00293 double G4GIDI_target::getElasticCrossSectionAtE( double e_in, double temperature ) {
00294 
00295     return( sumChannelCrossSectionAtE( nElasticIndices, elasticIndices, e_in, temperature ) );
00296 }
00297 /*
00298 ***************************************************************
00299 */
00300 double G4GIDI_target::getCaptureCrossSectionAtE( double e_in, double temperature ) {
00301 
00302     return( sumChannelCrossSectionAtE( nCaptureIndices, captureIndices, e_in, temperature ) );
00303 }
00304 /*
00305 ***************************************************************
00306 */
00307 double G4GIDI_target::getFissionCrossSectionAtE( double e_in, double temperature ) {
00308 
00309     return( sumChannelCrossSectionAtE( nFissionIndices, fissionIndices, e_in, temperature ) );
00310 }
00311 /*
00312 ***************************************************************
00313 */
00314 double G4GIDI_target::getOthersCrossSectionAtE( double e_in, double temperature ) {
00315 
00316     return( sumChannelCrossSectionAtE( nOthersIndices, othersIndices, e_in, temperature ) );
00317 }
00318 /*
00319 ***************************************************************
00320 */
00321 double G4GIDI_target::sumChannelCrossSectionAtE( int nIndices, int *indices, double e_in, double temperature ) {
00322 
00323     int i;
00324     double xsec = 0.;
00325 
00326     for( i = 0; i < nIndices; i++ ) 
00327         xsec += tpia_target_getIndexChannelCrossSectionAtE( &smr, target, indices[i], temperature, -1, e_in, tpia_crossSectionType_pointwise );
00328     return( xsec );
00329 }
00330 /*
00331 ***************************************************************
00332 */
00333 int G4GIDI_target::sampleChannelCrossSectionAtE( int nIndices, int *indices, double e_in, double temperature,
00334         double (*rng)( void * ), void *rngState ) {
00335 
00336     int i;
00337     double xsec = 0., rxsec = sumChannelCrossSectionAtE( nIndices, indices, e_in, temperature ) * tpia_misc_drng( rng, rngState );
00338 
00339     for( i = 0; i < nIndices - 1; i++ ) {
00340         xsec += tpia_target_getIndexChannelCrossSectionAtE( &smr, target, indices[i], temperature, -1, e_in, tpia_crossSectionType_pointwise );
00341         if( xsec >= rxsec ) break;
00342     }
00343     return( indices[i] );
00344 }
00345 /*
00346 ***************************************************************
00347 */
00348 //double G4GIDI_target::getElasticFinalState( double e_in, double temperature, double (*rng)( void * ), void *rngState ) {
00349 double G4GIDI_target::getElasticFinalState( double e_in, double , double (*rng)( void * ), void *rngState ) {
00350 
00351     tpia_decaySamplingInfo decaySamplingInfo;
00352     tpia_channel *channel = tpia_target_heated_getChannelAtIndex_smr( &smr, target->baseHeatedTarget, elasticIndices[0] );
00353     tpia_product *product;
00354 
00355     decaySamplingInfo.e_in = e_in;
00356     decaySamplingInfo.isVelocity = 0;
00357     tpia_frame_setColumn( &smr, &(decaySamplingInfo.frame), 0, tpia_referenceFrame_lab );
00358     decaySamplingInfo.samplingMethods = &(target->samplingMethods);
00359     decaySamplingInfo.rng = rng;
00360     decaySamplingInfo.rngState = rngState;
00361     product = tpia_decayChannel_getFirstProduct( &(channel->decayChannel) );
00362     tpia_angular_SampleMu( &smr, &(product->angular), &decaySamplingInfo );
00363 
00364     return( decaySamplingInfo.mu );
00365 }
00366 /*
00367 ***************************************************************
00368 */
00369 vector<G4GIDI_Product> *G4GIDI_target::getCaptureFinalState( double e_in, double temperature, double (*rng)( void * ), void *rngState ) {
00370 
00371     return( getFinalState( nCaptureIndices, captureIndices, e_in, temperature, rng, rngState ) );
00372 }
00373 /*
00374 ***************************************************************
00375 */
00376 vector<G4GIDI_Product> *G4GIDI_target::getFissionFinalState( double e_in, double temperature, double (*rng)( void * ), void *rngState ) {
00377 
00378     return( getFinalState( nFissionIndices, fissionIndices, e_in, temperature, rng, rngState ) );
00379 }
00380 /*
00381 ***************************************************************
00382 */
00383 vector<G4GIDI_Product> *G4GIDI_target::getOthersFinalState( double e_in, double temperature, double (*rng)( void * ), void *rngState ) {
00384 
00385     return( getFinalState( nOthersIndices, othersIndices, e_in, temperature, rng, rngState ) );
00386 }
00387 /*
00388 ***************************************************************
00389 */
00390 vector<G4GIDI_Product> *G4GIDI_target::getFinalState( int nIndices, int *indices, double e_in, double temperature, double (*rng)( void * ), void *rngState ) {
00391 
00392 #define nProductsMax 50
00393     int index = 0, i, n;
00394     vector<G4GIDI_Product> *products = NULL;
00395     tpia_decaySamplingInfo decaySamplingInfo;
00396     tpia_productOutgoingData productDatas[nProductsMax], *productData;
00397 
00398     decaySamplingInfo.e_in = e_in;
00399     decaySamplingInfo.samplingMethods = &(target->samplingMethods);
00400     decaySamplingInfo.isVelocity = 0;
00401     tpia_frame_setColumn( &smr, &(decaySamplingInfo.frame), 0, tpia_referenceFrame_lab );
00402     decaySamplingInfo.rng = rng;
00403     decaySamplingInfo.rngState = rngState;
00404 
00405     if( nIndices == 0 ) {
00406         return( NULL ); }
00407     else {
00408         if( nIndices == 1 ) {
00409             index = indices[0]; }
00410         else {
00411             index = sampleChannelCrossSectionAtE( nIndices, indices, e_in, temperature, rng, rngState );
00412         }
00413     }
00414     n = tpia_target_heated_sampleIndexChannelProductsAtE( &smr, target->baseHeatedTarget, index, &decaySamplingInfo, nProductsMax, productDatas );
00415     if( n > 0 ) {
00416         if( ( products = new vector<G4GIDI_Product>( n ) ) != NULL ) {
00417             for( i = 0; i < n; i++ ) {
00418                 productData = &(productDatas[i]);
00419                 (*products)[i].A = productData->productID->A;
00420                 (*products)[i].Z = productData->productID->Z;
00421                 (*products)[i].m = productData->productID->m;
00422                 (*products)[i].kineticEnergy = productData->kineticEnergy;
00423                 (*products)[i].px = productData->px_vx;
00424                 (*products)[i].py = productData->py_vy;
00425                 (*products)[i].pz = productData->pz_vz;
00426             }
00427         }
00428     }
00429 
00430     return( products );
00431 #undef nProductsMax
00432 }

Generated on Mon May 27 17:48:23 2013 for Geant4 by  doxygen 1.4.7