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

Last change on this file since 967 was 962, checked in by garnier, 15 years ago

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