ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MCGIDI_outputChannel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file 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_outputChannel* outputChannel,
211  MCGIDI_decaySamplingInfo* decaySamplingInfo,
212  MCGIDI_sampledProductsDatas* productDatas,
213  double *masses_ )
214 {
215  int i1;
216  int multiplicity(0);
217  int secondTwoBody = 0, isDecayChannel = ( outputChannel->reaction == NULL );
218  double e_in = modes.getProjectileEnergy( );
219  MCGIDI_product *product;
220  double phi, p, masses[3];
221  MCGIDI_distribution *distribution;
222  MCGIDI_sampledProductsData productData[2];
223 
224  if (isDecayChannel) {
225  masses[0] = masses_[0]; /* More work may be needed here. */
226  masses[1] = masses_[1];
227  } else {
228  masses[0] = MCGIDI_reaction_getProjectileMass_MeV( smr, outputChannel->reaction );
229  masses[1] = MCGIDI_reaction_getTargetMass_MeV( smr, outputChannel->reaction );
230  }
231 
232  // Loop over all possible final state particles reachable from initial state
233  // List of these particles (products) was read in from GIDI
234  // Note: all particles satifying the sampling criteria are included in the
235  // final state, regardless of charge, energy or baryon number conservation
236 
237  for (i1 = 0; i1 < outputChannel->numberOfProducts; i1++) {
238  product = &(outputChannel->products[i1]);
241  modes, decaySamplingInfo,
242  productDatas, masses ) < 0 ) return( -1 );
243  } else {
244  distribution = &(product->distribution);
245  if( distribution->type == MCGIDI_distributionType_none_e ) continue;
246 
247  if (!secondTwoBody) {
248  // Sample multiplicity of final state particle at kinetic energy of projectile
249  // The multiplicity stored in GIDI is a real number whose fractional part is
250  // compared to a random number to decide what integer value is returned
251  if ((multiplicity = product->multiplicity) == 0) multiplicity =
252  MCGIDI_product_sampleMultiplicity(smr, product, e_in,
253  decaySamplingInfo->rng( decaySamplingInfo->rngState ) );
254  while (multiplicity > 0) {
255 
256  multiplicity--;
257  decaySamplingInfo->pop = product->pop;
258  decaySamplingInfo->mu = 0;
259  decaySamplingInfo->Ep = 0;
260  productData[0].isVelocity = decaySamplingInfo->isVelocity;
261  productData[0].pop = product->pop;
262  productData[0].delayedNeutronIndex = product->delayedNeutronIndex;
263  productData[0].delayedNeutronRate = product->delayedNeutronRate;
264  productData[0].birthTimeSec = 0;
265  if (product->delayedNeutronRate > 0) {
266  productData[0].birthTimeSec =
267  -G4Log( decaySamplingInfo->rng( decaySamplingInfo->rngState ) ) / product->delayedNeutronRate;
268  }
269 
270  switch( outputChannel->genre ) {
271 
273  secondTwoBody = 1;
274  MCGIDI_angular_sampleMu( smr, distribution->angular, modes, decaySamplingInfo );
275  if (smr_isOk(smr) ) {
276  phi = 2. * M_PI * decaySamplingInfo->rng( decaySamplingInfo->rngState );
277  MCGIDI_kinetics_2BodyReaction( smr, distribution->angular, e_in, decaySamplingInfo->mu, phi, productData );
278  if (!smr_isOk(smr) ) return( -1 );
279  productData[1].pop = product[1].pop;
280  productData[1].delayedNeutronIndex = product[1].delayedNeutronIndex;
281  productData[1].delayedNeutronRate = product->delayedNeutronRate;
282  productData[1].birthTimeSec = 0;
283  MCGIDI_sampledProducts_addProduct( smr, productDatas, productData );
284  if( !smr_isOk( smr ) ) return( -1 );
285  MCGIDI_sampledProducts_addProduct( smr, productDatas, &(productData[1]) );
286  if( !smr_isOk( smr ) ) return( -1 );
287  }
288  break;
289 
291 
293  // Get mass of final state particle, then get its distribution
294  // masses[0] and masses[1] are incident and target masses
295  masses[2] = MCGIDI_product_getMass_MeV( smr, product );
296  switch( distribution->type ) {
298  MCGIDI_uncorrelated_sampleDistribution( smr, distribution, modes, decaySamplingInfo );
299  break;
301  MCGIDI_energyAngular_sampleDistribution( smr, distribution, modes, decaySamplingInfo );
302  break;
304  MCGIDI_KalbachMann_sampleEp( smr, distribution->KalbachMann, modes, decaySamplingInfo );
305  break;
307  MCGIDI_angularEnergy_sampleDistribution( smr, distribution->angularEnergy, modes, decaySamplingInfo );
308  break;
309  default :
310  printf( "Unknown spectral data form product name = %s, channel genre = %d\n", product->pop->name, outputChannel->genre );
311  break;
312  }
313  break;
314 
316  printf( "Channel is undefined\n" );
317  break;
318 
320  printf( "Channel is twoBodyDecay\n" );
321  break;
322 
324  printf( "Channel is uncorrelatedDecay\n" );
325  break;
326 
327  default :
328  printf( "Unsupported channel genre = %d\n", outputChannel->genre );
329  break;
330  }
331 
332  if (!smr_isOk(smr) ) return( -1 );
333  if (!secondTwoBody) {
334  if (decaySamplingInfo->frame == xDataTOM_frame_centerOfMass) {
335  if (MCGIDI_kinetics_COM2Lab( smr, modes, decaySamplingInfo, masses) != 0 ) return( -1 );
336  }
337 
338  // Assign kinematics to final state product
339  productData[0].kineticEnergy = decaySamplingInfo->Ep;
340  p = std::sqrt( decaySamplingInfo->Ep * ( decaySamplingInfo->Ep + 2. * product->pop->mass_MeV ) );
341  if (productData[0].isVelocity) p *= MCGIDI_speedOfLight_cm_sec / std::sqrt( p * p + product->pop->mass_MeV * product->pop->mass_MeV );
342  productData[0].pz_vz = p * decaySamplingInfo->mu;
343  p = std::sqrt( 1. - decaySamplingInfo->mu * decaySamplingInfo->mu ) * p;
344  phi = 2. * M_PI * decaySamplingInfo->rng( decaySamplingInfo->rngState );
345  productData[0].px_vx = p * std::sin( phi );
346  productData[0].py_vy = p * std::cos( phi );
347  MCGIDI_sampledProducts_addProduct( smr, productDatas, productData );
348  if (!smr_isOk(smr) ) return( -1 );
349  }
350  } // while multiplicity
351 
352  } // if !secondTwoBody
353  } // if decay channel genre
354 
355  } // loop over possible final state products
356  return( productDatas->numberOfProducts );
357 
358 }
359 
360 #if defined __cplusplus
361 }
362 #endif
363