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

Last change on this file was 1347, checked in by garnier, 14 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.