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

Last change on this file since 1196 was 1055, checked in by garnier, 15 years ago

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