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

Last change on this file since 829 was 819, checked in by garnier, 16 years ago

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