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

Last change on this file since 1330 was 1315, checked in by garnier, 15 years ago

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