Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
MCGIDI_outputChannel.cc
이 파일의 문서화 페이지로 가기
1 /*
2 # <<BEGIN-copyright>>
3 # <<END-copyright>>
4 */
5 #include <string.h>
6 #include <cmath>
7 
8 #include "MCGIDI.h"
9 #include "MCGIDI_misc.h"
10 
11 #if defined __cplusplus
12 #include "G4Log.hh"
13 namespace GIDI {
14 using namespace GIDI;
15 #endif
16 
17 /*
18 ************************************************************
19 */
21 
22  MCGIDI_outputChannel *outputChannel;
23 
24  if( ( outputChannel = (MCGIDI_outputChannel *) smr_malloc2( smr, sizeof( MCGIDI_outputChannel ), 0, "outputChannel" ) ) == NULL ) return( NULL );
25  if( MCGIDI_outputChannel_initialize( smr, outputChannel ) ) outputChannel = MCGIDI_outputChannel_free( smr, outputChannel );
26  return( outputChannel );
27 }
28 /*
29 ************************************************************
30 */
32 
33  memset( outputChannel, 0, sizeof( MCGIDI_outputChannel ) );
34  return( 0 );
35 }
36 /*
37 ************************************************************
38 */
40 
41  MCGIDI_outputChannel_release( smr, outputChannel );
42  smr_freeMemory( (void **) &outputChannel );
43  return( NULL );
44 }
45 /*
46 ************************************************************
47 */
49 
50  int i;
51 
52  for( i = 0; i < outputChannel->numberOfProducts; i++ ) MCGIDI_product_release( smr, &(outputChannel->products[i]) );
53  smr_freeMemory( (void **) &(outputChannel->products) );
54  MCGIDI_outputChannel_initialize( smr, outputChannel );
55 
56  return( 0 );
57 }
58 /*
59 ************************************************************
60 */
62  MCGIDI_reaction *reaction, MCGIDI_product *parent ) {
63 
64  int n, delayedNeutronIndex = 0;
65  char const *genre, *Q;
66  xDataTOM_element *child;
67 
68  MCGIDI_outputChannel_initialize( smr, outputChannel );
69 
70  outputChannel->reaction = reaction;
71  outputChannel->parent = parent;
72  if( ( genre = xDataTOM_getAttributesValueInElement( element, "genre" ) ) == NULL ) goto err;
73  if( ( parent != NULL ) && ( strcmp( genre, "NBody" ) ) ) {
74  smr_setReportError2( smr, smr_unknownID, 1, "decay channel's genre can only be 'uncorreclated' (a.k.a. 'NBody') and not '%s'", genre );
75  goto err;
76  }
77  if( strcmp( genre, "twoBody" ) == 0 ) {
78  outputChannel->genre = MCGIDI_channelGenre_twoBody_e; }
79  else if( strcmp( genre, "NBody" ) == 0 ) {
80  outputChannel->genre = MCGIDI_channelGenre_uncorrelated_e; }
81  else if( strcmp( genre, "sumOfRemainingOutputChannels" ) == 0 ) {
83  else {
84  smr_setReportError2( smr, smr_unknownID, 1, "unsupported genre = '%s'", genre );
85  goto err;
86  }
87  if( ( Q = xDataTOM_getAttributesValueInElement( element, "Q" ) ) == NULL ) goto err;
88  outputChannel->QIsFloat = !MCGIDI_misc_PQUStringToDoubleInUnitOf( smr, Q, "MeV", &(outputChannel->Q) );
89 
90  if( ( n = xDataTOM_numberOfElementsByName( smr, element, "product" ) ) == 0 ) {
91  smr_setReportError2p( smr, smr_unknownID, 1, "outputChannel does not have any products" );
92  goto err;
93  }
94  if( ( outputChannel->products = (MCGIDI_product *) smr_malloc2( smr, n * sizeof( MCGIDI_product ), 0, "outputChannel->products" ) ) == NULL ) goto err;
95 
96  for( child = xDataTOME_getFirstElement( element ); child != NULL; child = xDataTOME_getNextElement( child ) ) {
97  if( strcmp( child->name, "product" ) == 0 ) {
98  if( MCGIDI_product_parseFromTOM( smr, child, outputChannel, pops, &(outputChannel->products[outputChannel->numberOfProducts]),
99  &delayedNeutronIndex ) ) goto err;
100  outputChannel->numberOfProducts++; }
101  else if( strcmp( child->name, "fissionEnergyReleased" ) == 0 ) { /* ????????? Need to support. */
102  continue; }
103  else {
104  printf( "outputChannel child not currently supported = %s\n", child->name );
105  }
106  }
107  if( outputChannel->genre == MCGIDI_channelGenre_twoBody_e ) {
108  double projectileMass_MeV, targetMass_MeV, productMass_MeV, residualMass_MeV;
109 
110  projectileMass_MeV = MCGIDI_reaction_getProjectileMass_MeV( smr, reaction );
111  targetMass_MeV = MCGIDI_reaction_getTargetMass_MeV( smr, reaction );
112  productMass_MeV = MCGIDI_product_getMass_MeV( smr, &(outputChannel->products[0]) );
113  residualMass_MeV = MCGIDI_product_getMass_MeV( smr, &(outputChannel->products[1]) );
114 
115  //TK 17-11-10 for v1.3
116  //A temporary fix for emission of gamma(2.2MeV) from n captured by H
117  // capture gamma D
118  if ( reaction->ENDF_MT == 102 && productMass_MeV == 0 && ( outputChannel->products[1].pop->A == 2 && outputChannel->products[1].pop->Z == 1 ) ) {
119  //include/PoPs_data.h:#define e_Mass 5.4857990943e-4 /* electron mass in AMU */
120  residualMass_MeV += 5.4857990943e-4*MCGIDI_AMU2MeV;
121  }
122 
123  MCGIDI_product_setTwoBodyMasses( smr, &(outputChannel->products[0]), projectileMass_MeV, targetMass_MeV, productMass_MeV, residualMass_MeV );
124  }
125 
126  return( 0 );
127 
128 err:
129  MCGIDI_outputChannel_release( smr, outputChannel );
130  return( 1 );
131 }
132 /*
133 ************************************************************
134 */
136 
137  return( outputChannel->numberOfProducts );
138 }
139 /*
140 ************************************************************
141 */
143 
144  if( ( i < 0 ) || ( i >= outputChannel->numberOfProducts ) ) {
145  smr_setReportError2( smr, smr_unknownID, 1, "bad product index = %d: outputChannel as only %d products", i, outputChannel->numberOfProducts );
146  return( NULL );
147  }
148  return( &(outputChannel->products[i]) );
149 }
150 /*
151 ************************************************************
152 */
153 int MCGIDI_outputChannel_getDomain( statusMessageReporting *smr, MCGIDI_outputChannel *outputChannel, double *EMin, double *EMax ) {
154 
155  if( outputChannel->reaction != NULL ) return( MCGIDI_reaction_getDomain( smr, outputChannel->reaction, EMin, EMax ) );
156  return( MCGIDI_product_getDomain( smr, outputChannel->parent, EMin, EMax ) );
157 }
158 /*
159 ************************************************************
160 */
162 
163  if( outputChannel->reaction != NULL ) return( MCGIDI_reaction_getTargetHeated( smr, outputChannel->reaction ) );
164  return( MCGIDI_product_getTargetHeated( smr, outputChannel->parent ) );
165 }
166 /*
167 ************************************************************
168 */
170 
171  if( outputChannel->reaction != NULL ) return( MCGIDI_reaction_getProjectileMass_MeV( smr, outputChannel->reaction ) );
172  return( MCGIDI_product_getProjectileMass_MeV( smr, outputChannel->parent ) );
173 }
174 /*
175 ************************************************************
176 */
178 
179  if( outputChannel->reaction != NULL ) return( MCGIDI_reaction_getTargetMass_MeV( smr, outputChannel->reaction ) );
180  return( MCGIDI_product_getTargetMass_MeV( smr, outputChannel->parent ) );
181 }
182 /*
183 ************************************************************
184 */
185 double MCGIDI_outputChannel_getQ_MeV( statusMessageReporting * /*smr*/, MCGIDI_outputChannel *outputChannel, double /*e_in*/ ) {
186 
187  return( outputChannel->Q );
188 }
189 /*
190 ************************************************************
191 */
193 
194  int iProduct;
195  double Q = outputChannel->Q;
196  MCGIDI_product *product;
197 
198  for( iProduct = 0; iProduct < outputChannel->numberOfProducts; iProduct++ ) {
199  product = &(outputChannel->products[iProduct]);
200  if( product->decayChannel.genre != MCGIDI_channelGenre_undefined_e ) Q += MCGIDI_outputChannel_getFinalQ( smr, &(product->decayChannel), e_in );
201  if( !smr_isOk( smr ) ) break;
202  }
203  return( Q );
204 }
205 /*
206 ************************************************************
207 */
209  MCGIDI_decaySamplingInfo *decaySamplingInfo, MCGIDI_sampledProductsDatas *productDatas, double *masses_ ) {
210 
211  int i1, multiplicity, secondTwoBody = 0, isDecayChannel = ( outputChannel->reaction == NULL );
212  double e_in = modes.getProjectileEnergy( );
213  MCGIDI_product *product;
214  double phi, p, masses[3];
215  MCGIDI_distribution *distribution;
216  MCGIDI_sampledProductsData productData[2];
217 
218  if( isDecayChannel ) {
219  masses[0] = masses_[0]; /* More work may be needed here. */
220  masses[1] = masses_[1]; }
221  else {
222  masses[0] = MCGIDI_reaction_getProjectileMass_MeV( smr, outputChannel->reaction );
223  masses[1] = MCGIDI_reaction_getTargetMass_MeV( smr, outputChannel->reaction );
224  }
225 
226  for( i1 = 0; i1 < outputChannel->numberOfProducts; i1++ ) {
227  product = &(outputChannel->products[i1]);
229  if( MCGIDI_outputChannel_sampleProductsAtE( smr, &(product->decayChannel), modes, decaySamplingInfo, productDatas, masses ) < 0 ) return( -1 ); }
230  else {
231  distribution = &(product->distribution);
232  if( distribution->type == MCGIDI_distributionType_none_e ) continue;
233  if( !secondTwoBody ) {
234  if( ( multiplicity = product->multiplicity ) == 0 ) multiplicity = MCGIDI_product_sampleMultiplicity( smr, product, e_in,
235  decaySamplingInfo->rng( decaySamplingInfo->rngState ) );
236  while( multiplicity > 0 ) {
237 
238  multiplicity--;
239  decaySamplingInfo->pop = product->pop;
240  decaySamplingInfo->mu = 0;
241  decaySamplingInfo->Ep = 0;
242  productData[0].isVelocity = decaySamplingInfo->isVelocity;
243  productData[0].pop = product->pop;
244  productData[0].delayedNeutronIndex = product->delayedNeutronIndex;
245  productData[0].delayedNeutronRate = product->delayedNeutronRate;
246  productData[0].birthTimeSec = 0;
247  if( product->delayedNeutronRate > 0 ) {
248  productData[0].birthTimeSec = -G4Log( decaySamplingInfo->rng( decaySamplingInfo->rngState ) ) / product->delayedNeutronRate;
249  }
250 
251  switch( outputChannel->genre ) {
253  secondTwoBody = 1;
254  MCGIDI_angular_sampleMu( smr, distribution->angular, modes, decaySamplingInfo );
255  if( smr_isOk( smr ) ) {
256  phi = 2. * M_PI * decaySamplingInfo->rng( decaySamplingInfo->rngState );
257  MCGIDI_kinetics_2BodyReaction( smr, distribution->angular, e_in, decaySamplingInfo->mu, phi, productData );
258  if( !smr_isOk( smr ) ) return( -1 );
259  productData[1].pop = product[1].pop;
260  productData[1].delayedNeutronIndex = product[1].delayedNeutronIndex;
261  productData[1].delayedNeutronRate = product->delayedNeutronRate;
262  productData[1].birthTimeSec = 0;
263  MCGIDI_sampledProducts_addProduct( smr, productDatas, productData );
264  if( !smr_isOk( smr ) ) return( -1 );
265  MCGIDI_sampledProducts_addProduct( smr, productDatas, &(productData[1]) );
266  if( !smr_isOk( smr ) ) return( -1 );
267  }
268  break;
271  masses[2] = MCGIDI_product_getMass_MeV( smr, product );
272  switch( distribution->type ) {
274  MCGIDI_uncorrelated_sampleDistribution( smr, distribution, modes, decaySamplingInfo );
275  break;
277  MCGIDI_energyAngular_sampleDistribution( smr, distribution, modes, decaySamplingInfo );
278  break;
280  MCGIDI_KalbachMann_sampleEp( smr, distribution->KalbachMann, modes, decaySamplingInfo );
281  break;
283  MCGIDI_angularEnergy_sampleDistribution( smr, distribution->angularEnergy, modes, decaySamplingInfo );
284  break;
285  default :
286  printf( "Unknown spectral data form product name = %s, channel genre = %d\n", product->pop->name, outputChannel->genre );
287  break;
288  }
289  break;
291  printf( "Channel is undefined\n" );
292  break;
294  printf( "Channel is twoBodyDecay\n" );
295  break;
297  printf( "Channel is uncorrelatedDecay\n" );
298  break;
299  default :
300  printf( "Unsupported channel genre = %d\n", outputChannel->genre );
301  break;
302  }
303  if( !smr_isOk( smr ) ) return( -1 );
304  if( !secondTwoBody ) {
305  if( decaySamplingInfo->frame == xDataTOM_frame_centerOfMass ) {
306  if( MCGIDI_kinetics_COM2Lab( smr, modes, decaySamplingInfo, masses ) != 0 ) return( -1 );
307  }
308  productData[0].kineticEnergy = decaySamplingInfo->Ep;
309  p = std::sqrt( decaySamplingInfo->Ep * ( decaySamplingInfo->Ep + 2. * product->pop->mass_MeV ) );
310  if( productData[0].isVelocity ) p *= MCGIDI_speedOfLight_cm_sec / std::sqrt( p * p + product->pop->mass_MeV * product->pop->mass_MeV );
311  productData[0].pz_vz = p * decaySamplingInfo->mu;
312  p = std::sqrt( 1. - decaySamplingInfo->mu * decaySamplingInfo->mu ) * p;
313  phi = 2. * M_PI * decaySamplingInfo->rng( decaySamplingInfo->rngState );
314  productData[0].px_vx = p * std::sin( phi );
315  productData[0].py_vy = p * std::cos( phi );
316  MCGIDI_sampledProducts_addProduct( smr, productDatas, productData );
317  if( !smr_isOk( smr ) ) return( -1 );
318  }
319  } // Loop checking, 11.06.2015, T. Koi
320  }
321  }
322  }
323  return( productDatas->numberOfProducts );
324 }
325 
326 #if defined __cplusplus
327 }
328 #endif
329 
int MCGIDI_outputChannel_getDomain(statusMessageReporting *smr, MCGIDI_outputChannel *outputChannel, double *EMin, double *EMax)
double(* rng)(void *)
Definition: MCGIDI.h:258
xDataTOM_element * xDataTOME_getNextElement(xDataTOM_element *element)
Definition: xDataTOM.cc:238
double MCGIDI_reaction_getTargetMass_MeV(statusMessageReporting *smr, MCGIDI_reaction *reaction)
MCGIDI_reaction * reaction
Definition: MCGIDI.h:392
MCGIDI_target_heated * MCGIDI_product_getTargetHeated(statusMessageReporting *smr, MCGIDI_product *product)
int MCGIDI_sampledProducts_addProduct(statusMessageReporting *smr, MCGIDI_sampledProductsDatas *sampledProductsDatas, MCGIDI_sampledProductsData *sampledProductsData)
int MCGIDI_outputChannel_initialize(statusMessageReporting *smr, MCGIDI_outputChannel *outputChannel)
MCGIDI_product * parent
Definition: MCGIDI.h:393
double MCGIDI_product_getProjectileMass_MeV(statusMessageReporting *smr, MCGIDI_product *product)
#define smr_setReportError2(smr, libraryID, code, fmt,...)
int MCGIDI_product_getDomain(statusMessageReporting *smr, MCGIDI_product *product, double *EMin, double *EMax)
const char * p
Definition: xmltok.h:285
int xDataTOM_numberOfElementsByName(statusMessageReporting *smr, xDataTOM_element *element, char const *name)
Definition: xDataTOM.cc:268
int MCGIDI_kinetics_COM2Lab(statusMessageReporting *smr, MCGIDI_quantitiesLookupModes &modes, MCGIDI_decaySamplingInfo *decaySamplingInfo, double masses[3])
int MCGIDI_energyAngular_sampleDistribution(statusMessageReporting *smr, MCGIDI_distribution *distribution, MCGIDI_quantitiesLookupModes &modes, MCGIDI_decaySamplingInfo *decaySamplingInfo)
int MCGIDI_kinetics_2BodyReaction(statusMessageReporting *smr, MCGIDI_angular *angular, double K, double mu, double phi, MCGIDI_sampledProductsData *outgoingData)
int MCGIDI_angular_sampleMu(statusMessageReporting *smr, MCGIDI_angular *angular, MCGIDI_quantitiesLookupModes &modes, MCGIDI_decaySamplingInfo *decaySamplingInfo)
int multiplicity
Definition: MCGIDI.h:404
double MCGIDI_reaction_getProjectileMass_MeV(statusMessageReporting *smr, MCGIDI_reaction *reaction)
MCGIDI_POP * pop
Definition: MCGIDI.h:401
int MCGIDI_reaction_getDomain(statusMessageReporting *smr, MCGIDI_reaction *reaction, double *EMin, double *EMax)
xDataTOM_element * xDataTOME_getFirstElement(xDataTOM_element *element)
Definition: xDataTOM.cc:230
double MCGIDI_outputChannel_getFinalQ(statusMessageReporting *smr, MCGIDI_outputChannel *outputChannel, double e_in)
double MCGIDI_outputChannel_getProjectileMass_MeV(statusMessageReporting *smr, MCGIDI_outputChannel *outputChannel)
int MCGIDI_misc_PQUStringToDoubleInUnitOf(statusMessageReporting *smr, char const *str, char const *toUnit, double *value)
Definition: MCGIDI_misc.cc:356
enum MCGIDI_distributionType type
Definition: MCGIDI.h:382
#define M_PI
Definition: SbMath.h:34
G4double G4Log(G4double x)
Definition: G4Log.hh:230
enum MCGIDI_channelGenre genre
Definition: MCGIDI.h:391
MCGIDI_angularEnergy * angularEnergy
Definition: MCGIDI.h:386
int MCGIDI_outputChannel_numberOfProducts(MCGIDI_outputChannel *outputChannel)
MCGIDI_outputChannel * MCGIDI_outputChannel_free(statusMessageReporting *smr, MCGIDI_outputChannel *outputChannel)
MCGIDI_outputChannel decayChannel
Definition: MCGIDI.h:412
double MCGIDI_outputChannel_getQ_MeV(statusMessageReporting *smr, MCGIDI_outputChannel *outputChannel, double e_in)
#define smr_malloc2(smr, size, zero, forItem)
int MCGIDI_product_release(statusMessageReporting *smr, MCGIDI_product *product)
int MCGIDI_product_parseFromTOM(statusMessageReporting *smr, xDataTOM_element *element, MCGIDI_outputChannel *outputChannel, MCGIDI_POPs *pops, MCGIDI_product *product, int *delayedNeutronIndex)
char const * xDataTOM_getAttributesValueInElement(xDataTOM_element *element, char const *name)
Definition: xDataTOM.cc:286
int smr_isOk(statusMessageReporting *smr)
int MCGIDI_outputChannel_release(statusMessageReporting *smr, MCGIDI_outputChannel *outputChannel)
double delayedNeutronRate
Definition: MCGIDI.h:406
printf("%d Experimental points found\n", nlines)
char * name
Definition: MCGIDI.h:232
MCGIDI_KalbachMann * KalbachMann
Definition: MCGIDI.h:387
MCGIDI_product * MCGIDI_outputChannel_getProductAtIndex(statusMessageReporting *smr, MCGIDI_outputChannel *outputChannel, int i)
enum xDataTOM_frame frame
Definition: MCGIDI.h:256
#define smr_unknownID
int MCGIDI_outputChannel_parseFromTOM(statusMessageReporting *smr, xDataTOM_element *element, MCGIDI_POPs *pops, MCGIDI_outputChannel *outputChannel, MCGIDI_reaction *reaction, MCGIDI_product *parent)
double MCGIDI_product_getTargetMass_MeV(statusMessageReporting *smr, MCGIDI_product *product)
#define smr_setReportError2p(smr, libraryID, code, fmt)
int MCGIDI_KalbachMann_sampleEp(statusMessageReporting *smr, MCGIDI_KalbachMann *KalbachMann, MCGIDI_quantitiesLookupModes &modes, MCGIDI_decaySamplingInfo *decaySamplingInfo)
int delayedNeutronIndex
Definition: MCGIDI.h:405
void * smr_freeMemory(void **p)
int MCGIDI_product_setTwoBodyMasses(statusMessageReporting *smr, MCGIDI_product *product, double projectileMass_MeV, double targetMass_MeV, double productMass_MeV, double residualMass_MeV)
int MCGIDI_uncorrelated_sampleDistribution(statusMessageReporting *smr, MCGIDI_distribution *distribution, MCGIDI_quantitiesLookupModes &modes, MCGIDI_decaySamplingInfo *decaySamplingInfo)
int MCGIDI_product_sampleMultiplicity(statusMessageReporting *smr, MCGIDI_product *product, double e_in, double r)
MCGIDI_target_heated * MCGIDI_outputChannel_getTargetHeated(statusMessageReporting *smr, MCGIDI_outputChannel *outputChannel)
MCGIDI_product * products
Definition: MCGIDI.h:397
MCGIDI_outputChannel * MCGIDI_outputChannel_new(statusMessageReporting *smr)
double mass_MeV
Definition: MCGIDI.h:235
double getProjectileEnergy(void) const
Definition: MCGIDI.h:97
int MCGIDI_outputChannel_sampleProductsAtE(statusMessageReporting *smr, MCGIDI_outputChannel *outputChannel, MCGIDI_quantitiesLookupModes &modes, MCGIDI_decaySamplingInfo *decaySamplingInfo, MCGIDI_sampledProductsDatas *productDatas, double *masses)
MCGIDI_target_heated * MCGIDI_reaction_getTargetHeated(statusMessageReporting *smr, MCGIDI_reaction *reaction)
#define MCGIDI_speedOfLight_cm_sec
Definition: MCGIDI.h:183
Char_t n[5]
static double Q[]
MCGIDI_angular * angular
Definition: MCGIDI.h:383
double MCGIDI_outputChannel_getTargetMass_MeV(statusMessageReporting *smr, MCGIDI_outputChannel *outputChannel)
MCGIDI_distribution distribution
Definition: MCGIDI.h:411
double MCGIDI_product_getMass_MeV(statusMessageReporting *smr, MCGIDI_product *product)
#define MCGIDI_AMU2MeV
Definition: MCGIDI.h:184
int MCGIDI_angularEnergy_sampleDistribution(statusMessageReporting *smr, MCGIDI_angularEnergy *angularEnergy, MCGIDI_quantitiesLookupModes &modes, MCGIDI_decaySamplingInfo *decaySamplingInfo)