00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036 #include <string.h>
00037 #ifdef WIN32
00038 #define _USE_MATH_DEFINES
00039 #endif
00040 #include <cmath>
00041 #include "tpia_target.h"
00042
00043 #if defined __cplusplus
00044 namespace GIDI {
00045 using namespace GIDI;
00046 #endif
00047
00048
00049
00050
00051 tpia_product *tpia_decayChannel_getFirstProduct( tpia_decayChannel *decayChannel ) {
00052
00053 return( decayChannel->products );
00054 }
00055
00056
00057
00058 tpia_product *tpia_decayChannel_getNextProduct( tpia_product *product ) {
00059
00060 return( product->next );
00061 }
00062
00063
00064
00065 int tpia_decayChannel_sampleProductsAtE( statusMessageReporting *smr, tpia_decayChannel *decayChannel, tpia_decaySamplingInfo *decaySamplingInfo,
00066 int nProductData, tpia_productOutgoingData *productDatas ) {
00067
00068 int i, n = 0, multiplicity, secondTwoBody = 0, labFrame = tpia_referenceFrame_lab;
00069 tpia_product *product, *nextProduct;
00070 double phi, p;
00071
00072 if( nProductData < decayChannel->numberOfProducts ) {
00073 smr_setMessageError( smr, NULL, __FILE__, __LINE__, 1, "nProductData = %d < decayChannel->numberOfProducts = %d", nProductData,
00074 decayChannel->numberOfProducts );
00075 return( -1 );
00076 }
00077 for( i = 0, product = tpia_decayChannel_getFirstProduct( decayChannel ); product != NULL; i++, product = tpia_decayChannel_getNextProduct( product ) ) {
00078 nextProduct = tpia_decayChannel_getNextProduct( product );
00079 if( !secondTwoBody ) {
00080 if( ( multiplicity = product->multiplicity ) == 0 ) multiplicity = tpia_product_sampleMultiplicity( smr, product, decaySamplingInfo->e_in,
00081 tpia_misc_drng( decaySamplingInfo->rng, decaySamplingInfo->rngState ) );
00082 while( multiplicity > 0 ) {
00083 if( n >= nProductData ) break;
00084 multiplicity--;
00085 decaySamplingInfo->genre = product->genre;
00086 decaySamplingInfo->productID = product->productID;
00087 decaySamplingInfo->mu = 0;
00088 decaySamplingInfo->Ep = 0;
00089 productDatas[n].genre = product->genre;
00090 productDatas[n].isVelocity = decaySamplingInfo->isVelocity;
00091 tpia_frame_setColumns( smr, &(productDatas[n].frame), 1, &labFrame );
00092 productDatas[n].productID = product->productID;
00093 productDatas[n].decayChannel = &(product->decayChannel);
00094 if( strcmp( product->genre, "twoBody_angular" ) == 0 ) {
00095 secondTwoBody = 1;
00096 productDatas[n+1].productID = nextProduct->productID;
00097 productDatas[n].genre = product->genre;
00098 tpia_angular_SampleMu( smr, &(product->angular), decaySamplingInfo );
00099 if( smr_isOk( smr ) ) {
00100 phi = 2. * M_PI * tpia_misc_drng( decaySamplingInfo->rng, decaySamplingInfo->rngState );
00101 productDatas[n].isVelocity = decaySamplingInfo->isVelocity;
00102 productDatas[n].frame = decaySamplingInfo->frame;
00103 tpia_kinetics_2BodyReaction( smr, decayChannel, decaySamplingInfo->e_in, decaySamplingInfo->mu, phi, &(productDatas[n]) );
00104 } }
00105 else if( strcmp( product->genre, "NBody_Legendre" ) == 0 ) {
00106 tpia_Legendre_SampleEp( smr, &(product->Legendre), 1, decaySamplingInfo ); }
00107 else if( strcmp( product->genre, "NBody_angular_energy" ) == 0 ) {
00108 tpia_angular_SampleMu( smr, &(product->angular), decaySamplingInfo );
00109 tpia_angularEnergy_SampleEp( smr, &(product->angularEnergy), decaySamplingInfo ); }
00110 else if( strcmp( product->genre, "NBody_uncorrelate_Legendre" ) == 0 ) {
00111 tpia_angular_SampleMu( smr, &(product->angular), decaySamplingInfo );
00112 tpia_Legendre_SampleEp( smr, &(product->Legendre), 0, decaySamplingInfo ); }
00113 else if( strcmp( product->genre, "unknown" ) == 0 ) {
00114 }
00115 else {
00116 printf( "Unknown spectral data form product name = %s, genre = %s\n", product->productID->name, product->genre );
00117 }
00118 if( !smr_isOk( smr ) ) return( -1 );
00119 if( secondTwoBody ) {
00120 n++;
00121 productDatas[n].productID = nextProduct->productID;
00122 productDatas[n].genre = nextProduct->genre; }
00123 else {
00124 productDatas[n].kineticEnergy = decaySamplingInfo->Ep;
00125 p = std::sqrt( decaySamplingInfo->Ep * ( decaySamplingInfo->Ep + 2. * product->productID->fullMass_MeV ) );
00126 productDatas[n].pz_vz = p * decaySamplingInfo->mu;
00127 p = std::sqrt( 1. - decaySamplingInfo->mu * decaySamplingInfo->mu ) * p;
00128 phi = 2. * M_PI * tpia_misc_drng( decaySamplingInfo->rng, decaySamplingInfo->rngState );
00129 productDatas[n].px_vx = p * std::sin( phi );
00130 productDatas[n].py_vy = p * std::cos( phi );
00131 }
00132 n++;
00133 }
00134 }
00135 }
00136 return( n );
00137 }
00138
00139 #if defined __cplusplus
00140 }
00141 #endif