Geant4.10
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4NeutronHPInelasticCompFS.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 // neutron_hp -- source file
27 // J.P. Wellisch, Nov-1996
28 // A prototype of the low energy neutron transport model.
29 //
30 // 070523 bug fix for G4FPE_DEBUG on by A. Howard ( and T. Koi)
31 // 070606 bug fix and migrate to enable to Partial cases by T. Koi
32 // 080603 bug fix for Hadron Hyper News #932 by T. Koi
33 // 080612 bug fix contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #4,6
34 // 080717 bug fix of calculation of residual momentum by T. Koi
35 // 080801 protect negative avalable energy by T. Koi
36 // introduce theNDLDataA,Z which has A and Z of NDL data by T. Koi
37 // 081024 G4NucleiPropertiesTable:: to G4NucleiProperties::
38 // 090514 Fix bug in IC electron emission case
39 // Contribution from Chao Zhang (Chao.Zhang@usd.edu) and Dongming Mei(Dongming.Mei@usd.edu)
40 // 100406 "nothingWasKnownOnHadron=1" then sample mu isotropic in CM
41 // add two_body_reaction
42 // 100909 add safty
43 // 101111 add safty for _nat_ data case in Binary reaction, but break conservation
44 // 110430 add Reaction Q value and break up flag (MF3::QI and LR)
45 //
47 #include "G4NeutronHPManager.hh"
48 #include "G4PhysicalConstants.hh"
49 #include "G4SystemOfUnits.hh"
50 #include "G4Nucleus.hh"
51 #include "G4NucleiProperties.hh"
52 #include "G4He3.hh"
53 #include "G4Alpha.hh"
54 #include "G4Electron.hh"
55 #include "G4NeutronHPDataUsed.hh"
56 #include "G4IonTable.hh"
57 
59 {
60  // char the[100] = {""};
61  // std::ostrstream ost(the, 100, std::ios::out);
62  // ost <<gammaPath<<"z"<<ZR<<".a"<<AR;
63  // G4String * aName = new G4String(the);
64  // std::ifstream from(*aName, std::ios::in);
65 
66  std::ostringstream ost;
67  ost <<gammaPath<<"z"<<ZR<<".a"<<AR;
68  G4String aName = ost.str();
69  std::ifstream from(aName, std::ios::in);
70 
71  if(!from) return; // no data found for this isotope
72  // std::ifstream theGammaData(*aName, std::ios::in);
73  std::ifstream theGammaData(aName, std::ios::in);
74 
75  theGammas.Init(theGammaData);
76  // delete aName;
77 }
78 
80 {
81  gammaPath = "/Inelastic/Gammas/";
82  if(!getenv("G4NEUTRONHPDATA"))
83  throw G4HadronicException(__FILE__, __LINE__, "Please setenv G4NEUTRONHPDATA to point to the neutron cross-section files.");
84  G4String tBase = getenv("G4NEUTRONHPDATA");
85  gammaPath = tBase+gammaPath;
86  G4String tString = dirName;
87  G4bool dbool;
88  G4NeutronHPDataUsed aFile = theNames.GetName(static_cast<G4int>(A), static_cast<G4int>(Z), M, tString, aFSType, dbool);
89  G4String filename = aFile.GetName();
90  SetAZMs( A, Z, M, aFile );
91  //theBaseA = aFile.GetA();
92  //theBaseZ = aFile.GetZ();
93  //theNDLDataA = (int)aFile.GetA();
94  //theNDLDataZ = aFile.GetZ();
95  //if(!dbool || ( Z<2.5 && ( std::abs(theBaseZ - Z)>0.0001 || std::abs(theBaseA - A)>0.0001)))
96  if ( !dbool || ( Z<2.5 && ( std::abs(theNDLDataZ - Z)>0.0001 || std::abs(theNDLDataA - A)>0.0001)) )
97  {
98  if(getenv("NeutronHPNamesLogging")) G4cout << "Skipped = "<< filename <<" "<<A<<" "<<Z<<G4endl;
99  hasAnyData = false;
100  hasFSData = false;
101  hasXsec = false;
102  return;
103  }
104  //theBaseA = A;
105  //theBaseZ = G4int(Z+.5);
106  //std::ifstream theData(filename, std::ios::in);
107  std::istringstream theData(std::ios::in);
108  G4NeutronHPManager::GetInstance()->GetDataStream(filename,theData);
109  if(!theData) //"!" is a operator of ios
110  {
111  hasAnyData = false;
112  hasFSData = false;
113  hasXsec = false;
114  //theData.close();
115  return;
116  }
117  // here we go
118  G4int infoType, dataType, dummy;
119  G4int sfType, it;
120  hasFSData = false;
121  while (theData >> infoType)
122  {
123  hasFSData = true;
124  theData >> dataType;
125  theData >> sfType >> dummy;
126  it = 50;
127  if(sfType>=600||(sfType<100&&sfType>=50)) it = sfType%50;
128  if(dataType==3)
129  {
130  //theData >> dummy >> dummy;
131  //TK110430
132  // QI and LR introudced since G4NDL3.15
133  G4double dqi;
134  G4int ilr;
135  theData >> dqi >> ilr;
136 
137  QI[ it ] = dqi*eV;
138  LR[ it ] = ilr;
139  theXsection[it] = new G4NeutronHPVector;
140  G4int total;
141  theData >> total;
142  theXsection[it]->Init(theData, total, eV);
143  //std::cout << theXsection[it]->GetXsec(1*MeV) << std::endl;
144  }
145  else if(dataType==4)
146  {
148  theAngularDistribution[it]->Init(theData);
149  }
150  else if(dataType==5)
151  {
153  theEnergyDistribution[it]->Init(theData);
154  }
155  else if(dataType==6)
156  {
158  theEnergyAngData[it]->Init(theData);
159  }
160  else if(dataType==12)
161  {
163  theFinalStatePhotons[it]->InitMean(theData);
164  }
165  else if(dataType==13)
166  {
168  theFinalStatePhotons[it]->InitPartials(theData);
169  }
170  else if(dataType==14)
171  {
172  theFinalStatePhotons[it]->InitAngular(theData);
173  }
174  else if(dataType==15)
175  {
176  theFinalStatePhotons[it]->InitEnergies(theData);
177  }
178  else
179  {
180  throw G4HadronicException(__FILE__, __LINE__, "Data-type unknown to G4NeutronHPInelasticCompFS");
181  }
182  }
183  //theData.close();
184 }
185 
187 {
188 
189 // it = 0 has without Photon
190  G4double running[50];
191  running[0] = 0;
192  unsigned int i;
193  for(i=0; i<50; i++)
194  {
195  if(i!=0) running[i]=running[i-1];
196  if(theXsection[i] != 0)
197  {
198  running[i] += std::max(0., theXsection[i]->GetXsec(eKinetic));
199  }
200  }
201  G4double random = G4UniformRand();
202  G4double sum = running[49];
203  G4int it = 50;
204  if(0!=sum)
205  {
206  G4int i0;
207  for(i0=0; i0<50; i0++)
208  {
209  it = i0;
210  if(random < running[i0]/sum) break;
211  }
212  }
213 //debug: it = 1;
214  return it;
215 }
216 
217 
218  //n,p,d,t,he3,a
220 {
221 
222 // prepare neutron
223  theResult.Clear();
224  G4double eKinetic = theTrack.GetKineticEnergy();
225  const G4HadProjectile *incidentParticle = &theTrack;
226  G4ReactionProduct theNeutron( const_cast<G4ParticleDefinition *>(incidentParticle->GetDefinition()) );
227  theNeutron.SetMomentum( incidentParticle->Get4Momentum().vect() );
228  theNeutron.SetKineticEnergy( eKinetic );
229 
230 // prepare target
231  G4int i;
232  for(i=0; i<50; i++)
233  { if(theXsection[i] != 0) { break; } }
234 
235  G4double targetMass=0;
236  G4double eps = 0.0001;
237  targetMass = ( G4NucleiProperties::GetNuclearMass(static_cast<G4int>(theBaseA+eps), static_cast<G4int>(theBaseZ+eps))) /
239 // if(theEnergyAngData[i]!=0)
240 // targetMass = theEnergyAngData[i]->GetTargetMass();
241 // else if(theAngularDistribution[i]!=0)
242 // targetMass = theAngularDistribution[i]->GetTargetMass();
243 // else if(theFinalStatePhotons[50]!=0)
244 // targetMass = theFinalStatePhotons[50]->GetTargetMass();
245  G4Nucleus aNucleus;
247  G4ThreeVector neuVelo = (1./incidentParticle->GetDefinition()->GetPDGMass())*theNeutron.GetMomentum();
248  theTarget = aNucleus.GetBiasedThermalNucleus( targetMass, neuVelo, theTrack.GetMaterial()->GetTemperature());
249 
250 // prepare the residual mass
251  G4double residualMass=0;
252  G4double residualZ = theBaseZ - aDefinition->GetPDGCharge();
253  G4double residualA = theBaseA - aDefinition->GetBaryonNumber()+1;
254  residualMass = ( G4NucleiProperties::GetNuclearMass(static_cast<G4int>(residualA+eps), static_cast<G4int>(residualZ+eps)) ) /
256 
257 // prepare energy in target rest frame
258  G4ReactionProduct boosted;
259  boosted.Lorentz(theNeutron, theTarget);
260  eKinetic = boosted.GetKineticEnergy();
261 // G4double momentumInCMS = boosted.GetTotalMomentum();
262 
263 // select exit channel for composite FS class.
264  G4int it = SelectExitChannel( eKinetic );
265 
266 // set target and neutron in the relevant exit channel
267  InitDistributionInitialState(theNeutron, theTarget, it);
268 
269  G4ReactionProductVector * thePhotons = 0;
270  G4ReactionProductVector * theParticles = 0;
271  G4ReactionProduct aHadron;
272  aHadron.SetDefinition(aDefinition); // what if only cross-sections exist ==> Na 23 11 @@@@
273  G4double availableEnergy = theNeutron.GetKineticEnergy() + theNeutron.GetMass() - aHadron.GetMass() +
274  (targetMass - residualMass)*G4Neutron::Neutron()->GetPDGMass();
275 //080730c
276  if ( availableEnergy < 0 )
277  {
278  //G4cout << "080730c Adjust availavleEnergy " << G4endl;
279  availableEnergy = 0;
280  }
281  G4int nothingWasKnownOnHadron = 0;
282  G4int dummy;
283  G4double eGamm = 0;
284  G4int iLevel=it-1;
285 
286 // TK without photon has it = 0
287  if( 50 == it )
288  {
289 
290 // TK Excitation level is not determined
291  iLevel=-1;
292  aHadron.SetKineticEnergy(availableEnergy*residualMass*G4Neutron::Neutron()->GetPDGMass()/
293  (aHadron.GetMass()+residualMass*G4Neutron::Neutron()->GetPDGMass()));
294 
295  //aHadron.SetMomentum(theNeutron.GetMomentum()*(1./theNeutron.GetTotalMomentum())*
296  // std::sqrt(aHadron.GetTotalEnergy()*aHadron.GetTotalEnergy()-
297  // aHadron.GetMass()*aHadron.GetMass()));
298 
299  //TK add safty 100909
300  G4double p2 = ( aHadron.GetTotalEnergy()*aHadron.GetTotalEnergy() - aHadron.GetMass()*aHadron.GetMass() );
301  G4double p = 0.0;
302  if ( p2 > 0.0 ) p = std::sqrt( p );
303 
304  aHadron.SetMomentum(theNeutron.GetMomentum()*(1./theNeutron.GetTotalMomentum())*p );
305 
306  }
307  else
308  {
309  while( iLevel!=-1 && theGammas.GetLevel(iLevel) == 0 ) { iLevel--; }
310  }
311 
312 
313  if ( theAngularDistribution[it] != 0 ) // MF4
314  {
315  if(theEnergyDistribution[it]!=0) // MF5
316  {
317  //************************************************************
318  /*
319  aHadron.SetKineticEnergy(theEnergyDistribution[it]->Sample(eKinetic, dummy));
320  G4double eSecN = aHadron.GetKineticEnergy();
321  */
322  //************************************************************
323  //EMendoza --> maximum allowable energy should be taken into account.
324  G4double dqi = 0.0;
325  if ( QI[it] < 0 || 849 < QI[it] ) dqi = QI[it]; //For backword compatibility QI introduced since G4NDL3.15
326  G4double MaxEne=eKinetic+dqi;
327  G4double eSecN;
328  do{
329  eSecN=theEnergyDistribution[it]->Sample(eKinetic, dummy);
330  }while(eSecN>MaxEne);
331  aHadron.SetKineticEnergy(eSecN);
332  //************************************************************
333  eGamm = eKinetic-eSecN;
334  for(iLevel=theGammas.GetNumberOfLevels()-1; iLevel>=0; iLevel--)
335  {
336  if(theGammas.GetLevelEnergy(iLevel)<eGamm) break;
337  }
338  G4double random = 2*G4UniformRand();
339  iLevel+=G4int(random);
340  if(iLevel>theGammas.GetNumberOfLevels()-1)iLevel = theGammas.GetNumberOfLevels()-1;
341  }
342  else
343  {
344  G4double eExcitation = 0;
345  if(iLevel>=0) eExcitation = theGammas.GetLevel(iLevel)->GetLevelEnergy();
346  while (eKinetic-eExcitation < 0 && iLevel>0)
347  {
348  iLevel--;
349  eExcitation = theGammas.GetLevel(iLevel)->GetLevelEnergy();
350  }
351  //110610TK BEGIN
352  //Use QI value for calculating excitation energy of residual.
353  G4bool useQI=false;
354  G4double dqi = QI[it];
355  if ( dqi < 0 || 849 < dqi ) useQI = true; //Former libraies does not have values of this range
356 
357  if ( useQI )
358  {
359  // QI introudced since G4NDL3.15
360  eExcitation = -QI[it];
361  //Re-evluate iLevel based on this eExcitation
362  iLevel = 0;
363  G4bool find = false;
364  G4int imaxEx = 0;
365  while( !theGammas.GetLevel(iLevel+1) == 0 )
366  {
367  G4double maxEx = 0.0;
368  if ( maxEx < theGammas.GetLevel(iLevel)->GetLevelEnergy() )
369  {
370  maxEx = theGammas.GetLevel(iLevel)->GetLevelEnergy();
371  imaxEx = iLevel;
372  }
373  if ( eExcitation < theGammas.GetLevel(iLevel)->GetLevelEnergy() )
374  {
375  find = true;
376  iLevel--;
377  // very small eExcitation, iLevel becomes -1, this is protected below.
378  if ( iLevel == -1 ) iLevel = 0; // But cause energy trouble.
379  break;
380  }
381  iLevel++;
382  }
383  // In case, cannot find proper level, then use the maximum level.
384  if ( !find ) iLevel = imaxEx;
385  }
386  //110610TK END
387 
388  if(getenv("InelasticCompFSLogging") && eKinetic-eExcitation < 0)
389  {
390  throw G4HadronicException(__FILE__, __LINE__, "SEVERE: InelasticCompFS: Consistency of data not good enough, please file report");
391  }
392  if(eKinetic-eExcitation < 0) eExcitation = 0;
393  if(iLevel!= -1) aHadron.SetKineticEnergy(eKinetic - eExcitation);
394 
395  }
397 
398  if( theFinalStatePhotons[it] == 0 )
399  {
400  //G4cout << "110610 USE Gamma Level" << G4endl;
401 // TK comment Most n,n* eneter to this
402  thePhotons = theGammas.GetDecayGammas(iLevel);
403  eGamm -= theGammas.GetLevelEnergy(iLevel);
404  if(eGamm>0) // @ ok for now, but really needs an efficient way of correllated sampling @
405  {
406  G4ReactionProduct * theRestEnergy = new G4ReactionProduct;
407  theRestEnergy->SetDefinition(G4Gamma::Gamma());
408  theRestEnergy->SetKineticEnergy(eGamm);
409  G4double costh = 2.*G4UniformRand()-1.;
410  G4double phi = twopi*G4UniformRand();
411  theRestEnergy->SetMomentum(eGamm*std::sin(std::acos(costh))*std::cos(phi),
412  eGamm*std::sin(std::acos(costh))*std::sin(phi),
413  eGamm*costh);
414  if(thePhotons == 0) { thePhotons = new G4ReactionProductVector; }
415  thePhotons->push_back(theRestEnergy);
416  }
417  }
418  }
419  else if(theEnergyAngData[it] != 0) // MF6
420  {
421  theParticles = theEnergyAngData[it]->Sample(eKinetic);
422  }
423  else
424  {
425  // @@@ what to do, if we have photon data, but no info on the hadron itself
426  nothingWasKnownOnHadron = 1;
427  }
428 
429  //G4cout << "theFinalStatePhotons it " << it << G4endl;
430  //G4cout << "theFinalStatePhotons[it] " << theFinalStatePhotons[it] << G4endl;
431  //G4cout << "theFinalStatePhotons it " << it << G4endl;
432  //G4cout << "theFinalStatePhotons[it] " << theFinalStatePhotons[it] << G4endl;
433  //G4cout << "thePhotons " << thePhotons << G4endl;
434 
435  if ( theFinalStatePhotons[it] != 0 )
436  {
437  // the photon distributions are in the Nucleus rest frame.
438  // TK residual rest frame
439  G4ReactionProduct boosted_tmp;
440  boosted_tmp.Lorentz(theNeutron, theTarget);
441  G4double anEnergy = boosted_tmp.GetKineticEnergy();
442  thePhotons = theFinalStatePhotons[it]->GetPhotons(anEnergy);
443  G4double aBaseEnergy = theFinalStatePhotons[it]->GetLevelEnergy();
444  G4double testEnergy = 0;
445  if(thePhotons!=0 && thePhotons->size()!=0)
446  { aBaseEnergy-=thePhotons->operator[](0)->GetTotalEnergy(); }
447  if(theFinalStatePhotons[it]->NeedsCascade())
448  {
449  while(aBaseEnergy>0.01*keV)
450  {
451  // cascade down the levels
452  G4bool foundMatchingLevel = false;
453  G4int closest = 2;
454  G4double deltaEold = -1;
455  for(G4int j=1; j<it; j++)
456  {
457  if(theFinalStatePhotons[j]!=0)
458  {
459  testEnergy = theFinalStatePhotons[j]->GetLevelEnergy();
460  }
461  else
462  {
463  testEnergy = 0;
464  }
465  G4double deltaE = std::abs(testEnergy-aBaseEnergy);
466  if(deltaE<0.1*keV)
467  {
468  G4ReactionProductVector * theNext =
469  theFinalStatePhotons[j]->GetPhotons(anEnergy);
470  thePhotons->push_back(theNext->operator[](0));
471  aBaseEnergy = testEnergy-theNext->operator[](0)->GetTotalEnergy();
472  delete theNext;
473  foundMatchingLevel = true;
474  break; // ===>
475  }
476  if(theFinalStatePhotons[j]!=0 && ( deltaE<deltaEold||deltaEold<0.) )
477  {
478  closest = j;
479  deltaEold = deltaE;
480  }
481  } // <=== the break goes here.
482  if(!foundMatchingLevel)
483  {
484  G4ReactionProductVector * theNext =
485  theFinalStatePhotons[closest]->GetPhotons(anEnergy);
486  thePhotons->push_back(theNext->operator[](0));
487  aBaseEnergy = aBaseEnergy-theNext->operator[](0)->GetTotalEnergy();
488  delete theNext;
489  }
490  }
491  }
492  }
493  unsigned int i0;
494  if(thePhotons!=0)
495  {
496  for(i0=0; i0<thePhotons->size(); i0++)
497  {
498  // back to lab
499  thePhotons->operator[](i0)->Lorentz(*(thePhotons->operator[](i0)), -1.*theTarget);
500  }
501  }
502  //G4cout << "nothingWasKnownOnHadron " << nothingWasKnownOnHadron << G4endl;
503  if(nothingWasKnownOnHadron)
504  {
505 // TKDB 100405
506 // In this case, hadron should be isotropic in CM
507 // mu and p should be correlated
508 //
509  G4double totalPhotonEnergy = 0.0;
510  if ( thePhotons != 0 )
511  {
512  unsigned int nPhotons = thePhotons->size();
513  unsigned int ii0;
514  for ( ii0=0; ii0<nPhotons; ii0++)
515  {
516  //thePhotons has energies at LAB system
517  totalPhotonEnergy += thePhotons->operator[](ii0)->GetTotalEnergy();
518  }
519  }
520 
521  //isotropic distribution in CM
522  G4double mu = 1.0 - 2 * G4UniformRand();
523 
524  // need momentums in target rest frame;
525  G4LorentzVector target_in_LAB ( theTarget.GetMomentum() , theTarget.GetTotalEnergy() );
526  G4ThreeVector boostToTargetRest = -target_in_LAB.boostVector();
527  G4LorentzVector proj_in_LAB = incidentParticle->Get4Momentum();
528 
529  G4DynamicParticle* proj = new G4DynamicParticle( G4Neutron::Neutron() , proj_in_LAB.boost( boostToTargetRest ) );
530  G4DynamicParticle* targ = new G4DynamicParticle( G4IonTable::GetIonTable()->GetIon ( (G4int)theBaseZ , (G4int)theBaseA , totalPhotonEnergy ) , G4ThreeVector(0) );
531  G4DynamicParticle* hadron = new G4DynamicParticle( aHadron.GetDefinition() , G4ThreeVector(0) ); // will be fill momentum
532 
533  two_body_reaction ( proj , targ , hadron , mu );
534 
535  G4LorentzVector hadron_in_trag_rest = hadron->Get4Momentum();
536  G4LorentzVector hadron_in_LAB = hadron_in_trag_rest.boost ( -boostToTargetRest );
537  aHadron.SetMomentum( hadron_in_LAB.v() );
538  aHadron.SetKineticEnergy ( hadron_in_LAB.e() - hadron_in_LAB.m() );
539 
540  delete proj;
541  delete targ;
542  delete hadron;
543 
544 //TKDB 100405
545 /*
546  G4double totalPhotonEnergy = 0;
547  if(thePhotons!=0)
548  {
549  unsigned int nPhotons = thePhotons->size();
550  unsigned int i0;
551  for(i0=0; i0<nPhotons; i0++)
552  {
553  totalPhotonEnergy += thePhotons->operator[](i0)->GetTotalEnergy();
554  }
555  }
556  availableEnergy -= totalPhotonEnergy;
557  residualMass += totalPhotonEnergy/G4Neutron::Neutron()->GetPDGMass();
558  aHadron.SetKineticEnergy(availableEnergy*residualMass*G4Neutron::Neutron()->GetPDGMass()/
559  (aHadron.GetMass()+residualMass*G4Neutron::Neutron()->GetPDGMass()));
560  G4double CosTheta = 1.0 - 2.0*G4UniformRand();
561  G4double SinTheta = std::sqrt(1.0 - CosTheta*CosTheta);
562  G4double Phi = twopi*G4UniformRand();
563  G4ThreeVector Vector(std::cos(Phi)*SinTheta, std::sin(Phi)*SinTheta, CosTheta);
564  //aHadron.SetMomentum(Vector* std::sqrt(aHadron.GetTotalEnergy()*aHadron.GetTotalEnergy()-
565  // aHadron.GetMass()*aHadron.GetMass()));
566  G4double p2 = aHadron.GetTotalEnergy()*aHadron.GetTotalEnergy()- aHadron.GetMass()*aHadron.GetMass();
567 
568  G4double p = 0.0;
569  if ( p2 > 0.0 )
570  p = std::sqrt ( p2 );
571 
572  aHadron.SetMomentum( Vector*p );
573 */
574 
575  }
576 
577 // fill the result
578 // Beware - the recoil is not necessarily in the particles...
579 // Can be calculated from momentum conservation?
580 // The idea is that the particles ar emitted forst, and the gammas only once the
581 // recoil is on the residual; assumption is that gammas do not contribute to
582 // the recoil.
583 // This needs more design @@@
584 
585  G4int nSecondaries = 2; // the hadron and the recoil
586  G4bool needsSeparateRecoil = false;
587  G4int totalBaryonNumber = 0;
588  G4int totalCharge = 0;
589  G4ThreeVector totalMomentum(0);
590  if(theParticles != 0)
591  {
592  nSecondaries = theParticles->size();
593  G4ParticleDefinition * aDef;
594  unsigned int ii0;
595  for(ii0=0; ii0<theParticles->size(); ii0++)
596  {
597  aDef = theParticles->operator[](ii0)->GetDefinition();
598  totalBaryonNumber+=aDef->GetBaryonNumber();
599  totalCharge+=G4int(aDef->GetPDGCharge()+eps);
600  totalMomentum += theParticles->operator[](ii0)->GetMomentum();
601  }
602  if(totalBaryonNumber!=G4int(theBaseA+eps+incidentParticle->GetDefinition()->GetBaryonNumber()))
603  {
604  needsSeparateRecoil = true;
605  nSecondaries++;
606  residualA = G4int(theBaseA+eps+incidentParticle->GetDefinition()->GetBaryonNumber()
607  -totalBaryonNumber);
608  residualZ = G4int(theBaseZ+eps+incidentParticle->GetDefinition()->GetPDGCharge()
609  -totalCharge);
610  }
611  }
612 
613  G4int nPhotons = 0;
614  if(thePhotons!=0) { nPhotons = thePhotons->size(); }
615  nSecondaries += nPhotons;
616 
617  G4DynamicParticle * theSec;
618 
619  if( theParticles==0 )
620  {
621  theSec = new G4DynamicParticle;
622  theSec->SetDefinition(aHadron.GetDefinition());
623  theSec->SetMomentum(aHadron.GetMomentum());
624  theResult.AddSecondary(theSec);
625 
626  aHadron.Lorentz(aHadron, theTarget);
627  G4ReactionProduct theResidual;
629  ->GetIon(static_cast<G4int>(residualZ), static_cast<G4int>(residualA), 0));
630  theResidual.SetKineticEnergy(aHadron.GetKineticEnergy()*aHadron.GetMass()/theResidual.GetMass());
631 
632  //080612TK contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #6
633  //theResidual.SetMomentum(-1.*aHadron.GetMomentum());
634  G4ThreeVector incidentNeutronMomentum = theNeutron.GetMomentum();
635  theResidual.SetMomentum(incidentNeutronMomentum - aHadron.GetMomentum());
636 
637  theResidual.Lorentz(theResidual, -1.*theTarget);
638  G4ThreeVector totalPhotonMomentum(0,0,0);
639  if(thePhotons!=0)
640  {
641  for(i=0; i<nPhotons; i++)
642  {
643  totalPhotonMomentum += thePhotons->operator[](i)->GetMomentum();
644  }
645  }
646  theSec = new G4DynamicParticle;
647  theSec->SetDefinition(theResidual.GetDefinition());
648  theSec->SetMomentum(theResidual.GetMomentum()-totalPhotonMomentum);
649  theResult.AddSecondary(theSec);
650  }
651  else
652  {
653  for(i0=0; i0<theParticles->size(); i0++)
654  {
655  theSec = new G4DynamicParticle;
656  theSec->SetDefinition(theParticles->operator[](i0)->GetDefinition());
657  theSec->SetMomentum(theParticles->operator[](i0)->GetMomentum());
658  theResult.AddSecondary(theSec);
659  delete theParticles->operator[](i0);
660  }
661  delete theParticles;
662  if(needsSeparateRecoil && residualZ!=0)
663  {
664  G4ReactionProduct theResidual;
666  ->GetIon(static_cast<G4int>(residualZ), static_cast<G4int>(residualA), 0));
667  G4double resiualKineticEnergy = theResidual.GetMass()*theResidual.GetMass();
668  resiualKineticEnergy += totalMomentum*totalMomentum;
669  resiualKineticEnergy = std::sqrt(resiualKineticEnergy) - theResidual.GetMass();
670 // cout << "Kinetic energy of the residual = "<<resiualKineticEnergy<<endl;
671  theResidual.SetKineticEnergy(resiualKineticEnergy);
672 
673  //080612TK contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #4
674  //theResidual.SetMomentum(-1.*totalMomentum);
675  //G4ThreeVector incidentNeutronMomentum = theNeutron.GetMomentum();
676  //theResidual.SetMomentum(incidentNeutronMomentum - aHadron.GetMomentum());
677 //080717 TK Comment still do NOT include photon's mometum which produce by thePhotons
678  theResidual.SetMomentum( theNeutron.GetMomentum() + theTarget.GetMomentum() - totalMomentum );
679 
680  theSec = new G4DynamicParticle;
681  theSec->SetDefinition(theResidual.GetDefinition());
682  theSec->SetMomentum(theResidual.GetMomentum());
683  theResult.AddSecondary(theSec);
684  }
685  }
686  if(thePhotons!=0)
687  {
688  for(i=0; i<nPhotons; i++)
689  {
690  theSec = new G4DynamicParticle;
691  //Bug reported Chao Zhang (Chao.Zhang@usd.edu), Dongming Mei(Dongming.Mei@usd.edu) Feb. 25, 2009
692  //theSec->SetDefinition(G4Gamma::Gamma());
693  theSec->SetDefinition( thePhotons->operator[](i)->GetDefinition() );
694  //But never cause real effect at least with G4NDL3.13 TK
695  theSec->SetMomentum(thePhotons->operator[](i)->GetMomentum());
696  theResult.AddSecondary(theSec);
697  delete thePhotons->operator[](i);
698  }
699 // some garbage collection
700  delete thePhotons;
701  }
702 
703 //080721
705  G4LorentzVector targ_4p_lab ( theTarget.GetMomentum() , std::sqrt( targ_pd->GetPDGMass()*targ_pd->GetPDGMass() + theTarget.GetMomentum().mag2() ) );
706  G4LorentzVector proj_4p_lab = theTrack.Get4Momentum();
707  G4LorentzVector init_4p_lab = proj_4p_lab + targ_4p_lab;
708  adjust_final_state ( init_4p_lab );
709 
710 // clean up the primary neutron
712 }
713 
714 
715 
716 #include "G4RotationMatrix.hh"
717 void G4NeutronHPInelasticCompFS::two_body_reaction ( G4DynamicParticle* proj, G4DynamicParticle* targ, G4DynamicParticle* hadron, G4double mu )
718 {
719 
720 // Target rest flame
721 // 4vector in targ rest frame;
722 // targ could have excitation energy (photon energy will be emiited) tricky but,,,
723 
724  G4LorentzVector before = proj->Get4Momentum() + targ->Get4Momentum();
725 
726  G4ThreeVector p3_proj = proj->GetMomentum();
727  G4ThreeVector d = p3_proj.unit();
728  G4RotationMatrix rot;
729  G4RotationMatrix rot1;
730  rot1.setPhi( pi/2 + d.phi() );
731  G4RotationMatrix rot2;
732  rot2.setTheta( d.theta() );
733  rot=rot2*rot1;
734  proj->SetMomentum( rot*p3_proj );
735 
736 // Now proj only has pz component;
737 
738 // mu in CM system
739 
740  //Valid only for neutron incidence
741  G4DynamicParticle* residual = new G4DynamicParticle ( G4IonTable::GetIonTable()->GetIon ( (G4int)( targ->GetDefinition()->GetPDGCharge() - hadron->GetDefinition()->GetPDGCharge() ) , (G4int)(targ->GetDefinition()->GetBaryonNumber() - hadron->GetDefinition()->GetBaryonNumber()+1) , 0 ) , G4ThreeVector(0) );
742 
743  G4double Q = proj->GetDefinition()->GetPDGMass() + targ->GetDefinition()->GetPDGMass()
744  - ( hadron->GetDefinition()->GetPDGMass() + residual->GetDefinition()->GetPDGMass() );
745 
746  // Non Relativistic Case
747  G4double A = targ->GetDefinition()->GetPDGMass() / proj->GetDefinition()->GetPDGMass();
748  G4double AA = hadron->GetDefinition()->GetPDGMass() / proj->GetDefinition()->GetPDGMass();
749  G4double E1 = proj->GetKineticEnergy();
750 
751 // 101111
752 // In _nat_ data (Q+E1) could become negative value, following line is safty for this case.
753  //if ( (Q+E1) < 0 )
754  if ( ( 1 + (1+A)/A*Q/E1 ) < 0 )
755  {
756 // 1.0e-6 eV is additional safty for numeric precision
757  Q = -( A/(1+A)*E1 ) + 1.0e-6*eV;
758  }
759 
760  G4double beta = std::sqrt ( A*(A+1-AA)/AA*( 1 + (1+A)/A*Q/E1 ) );
761  G4double gamma = AA/(A+1-AA)*beta;
762  G4double E3 = AA/std::pow((1+A),2)*(beta*beta+1+2*beta*mu)*E1;
763  G4double omega3 = (1+beta*mu)/std::sqrt(beta*beta+1+2*beta*mu);
764  if ( omega3 > 1.0 ) omega3 = 1.0;
765 
766  G4double E4 = (A+1-AA)/std::pow((1+A),2)*(gamma*gamma+1-2*gamma*mu)*E1;
767  G4double omega4 = (1-gamma*mu)/std::sqrt(gamma*gamma+1-2*gamma*mu);
768  if ( omega4 > 1.0 ) omega4 = 1.0;
769 
770  hadron->SetKineticEnergy ( E3 );
771 
772  G4double M = hadron->GetDefinition()->GetPDGMass();
773  G4double pmag = std::sqrt ((E3+M)*(E3+M)-M*M) ;
774  G4ThreeVector p ( 0 , pmag*std::sqrt(1-omega3*omega3), pmag*omega3 );
775 
776  G4double M4 = residual->GetDefinition()->GetPDGMass();
777  G4double pmag4 = std::sqrt ((E4+M4)*(E4+M4)-M4*M4) ;
778  G4ThreeVector p4 ( 0 , -pmag4*std::sqrt(1-omega4*omega4), pmag4*omega4 );
779 
780 // Rotate to orginal target rest flame.
781  p *= rot.inverse();
782  hadron->SetMomentum( p );
783 // Now hadron had 4 momentum in target rest flame
784 
785 // TypeA
786  p4 *= rot.inverse();
787  residual->SetMomentum ( p4 );
788 
789 //TypeB1
790  //residual->Set4Momentum ( p4_residual );
791 //TypeB2
792  //residual->SetMomentum ( p4_residual.v() );
793 
794 // Type A make difference in Momenutum
795 // Type B1 make difference in Mass of residual
796 // Type B2 make difference in total energy.
797 
798  delete residual;
799 
800 }
G4ReactionProductVector * GetPhotons(G4double anEnergy)
static G4double GetNuclearMass(const G4double A, const G4double Z)
void SetMomentum(const G4ThreeVector &momentum)
G4double GetTotalMomentum() const
G4NeutronHPLevel * GetLevel(G4int i)
void Lorentz(const G4ReactionProduct &p1, const G4ReactionProduct &p2)
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
void SetKineticEnergy(const G4double en)
static G4NeutronHPManager * GetInstance()
void SetMomentum(const G4double x, const G4double y, const G4double z)
const char * p
Definition: xmltok.h:285
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
Definition: G4IonTable.cc:449
void Init(std::istream &aDataFile)
Hep3Vector v() const
G4double Sample(G4double anEnergy, G4int &it)
virtual G4NeutronHPVector * GetXsec()
void Init(std::istream &aDataFile, G4int total, G4double ux=1., G4double uy=1.)
G4ParticleDefinition * GetDefinition() const
void GetDataStream(G4String, std::istringstream &iss)
G4NeutronHPDataUsed GetName(G4int A, G4int Z, G4String base, G4String rest, G4bool &active)
int G4int
Definition: G4Types.hh:78
HepRotation inverse() const
G4ParticleDefinition * GetDefinition() const
void SetStatusChange(G4HadFinalStateStatus aS)
std::vector< G4ReactionProduct * > G4ReactionProductVector
void InitAngular(std::istream &aDataFile)
G4ReactionProduct GetBiasedThermalNucleus(G4double aMass, G4ThreeVector aVelocity, G4double temp=-1) const
Definition: G4Nucleus.cc:108
Hep3Vector vect() const
#define G4UniformRand()
Definition: Randomize.hh:87
G4GLOB_DLL std::ostream G4cout
G4NeutronHPEnergyDistribution * theEnergyDistribution[51]
const G4ParticleDefinition * GetDefinition() const
void CompositeApply(const G4HadProjectile &theTrack, G4ParticleDefinition *aHadron)
bool G4bool
Definition: G4Types.hh:79
HepLorentzVector & boost(double, double, double)
G4double GetKineticEnergy() const
void setTheta(double theta)
Definition: RotationE.cc:266
G4int SelectExitChannel(G4double eKinetic)
G4ErrorTarget * theTarget
Definition: errprop.cc:59
void InitEnergies(std::istream &aDataFile)
void Init(std::istream &aDataFile)
double phi() const
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
void setPhi(double phi)
Definition: RotationE.cc:262
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
const G4LorentzVector & Get4Momentum() const
void InitGammas(G4double AR, G4double ZR)
void SetKineticEnergy(G4double aEnergy)
G4NeutronHPEnAngCorrelation * theEnergyAngData[51]
G4LorentzVector Get4Momentum() const
static G4IonTable * GetIonTable()
Definition: G4IonTable.hh:80
double theta() const
G4double GetKineticEnergy() const
void SampleAndUpdate(G4ReactionProduct &aNeutron)
G4NeutronHPAngular * theAngularDistribution[51]
G4double GetTotalEnergy() const
G4bool InitMean(std::istream &aDataFile)
G4double total(Particle const *const p1, Particle const *const p2)
G4double GetPDGMass() const
G4double GetLevelEnergy(G4int aLevel)
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4NeutronHPPhotonDist * theFinalStatePhotons[51]
Hep3Vector unit() const
void SetDefinition(G4ParticleDefinition *aParticleDefinition)
double mag2() const
G4ThreeVector GetMomentum() const
G4double GetTemperature() const
Definition: G4Material.hh:180
#define G4endl
Definition: G4ios.hh:61
const G4Material * GetMaterial() const
G4ReactionProductVector * Sample(G4double anEnergy)
void Init(G4double A, G4double Z, G4int M, G4String &dirName, G4String &aSFType)
G4ReactionProductVector * GetDecayGammas(G4int aLevel)
void SetAZMs(G4double anA, G4double aZ, G4int aM, G4NeutronHPDataUsed used)
double G4double
Definition: G4Types.hh:76
void InitDistributionInitialState(G4ReactionProduct &aNeutron, G4ReactionProduct &aTarget, G4int it)
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
void adjust_final_state(G4LorentzVector)
G4double GetPDGCharge() const
G4double GetLevelEnergy()
G4double GetMass() const
void InitPartials(std::istream &aDataFile)
void AddSecondary(G4DynamicParticle *aP)
void Init(std::istream &aDataFile)
G4ThreeVector GetMomentum() const