source: trunk/source/processes/hadronic/models/neutron_hp/src/G4NeutronHPInelasticCompFS.cc@ 1350

Last change on this file since 1350 was 1347, checked in by garnier, 15 years ago

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