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

Last change on this file was 1347, checked in by garnier, 14 years ago

geant4 tag 9.4
File size: 9.6 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//
27//080721 Create adjust_final_state method by T. Koi
28//080801 Residual reconstruction with theNDLDataA,Z (A, Z, and momentum are adjusted) by T. Koi
29//101110 Set lower limit for gamma energy(1keV) by T. Koi
30
31#include "G4NeutronHPFinalState.hh"
32
33#include "G4ParticleTable.hh"
34#include "G4Gamma.hh"
35
36
37
38void G4NeutronHPFinalState::adjust_final_state ( G4LorentzVector init_4p_lab )
39{
40
41   G4double minimum_energy = 1*keV;
42
43   if ( adjustResult != true ) return;
44
45   G4int nSecondaries = theResult.GetNumberOfSecondaries();
46
47   G4int sum_Z = 0;
48   G4int sum_A = 0;
49   G4int max_SecZ = 0;
50   G4int max_SecA = 0;
51   for ( int i = 0 ; i < nSecondaries ; i++ )
52   {
53      sum_Z += theResult.GetSecondary( i )->GetParticle()->GetDefinition()->GetAtomicNumber();
54      max_SecZ = std::max ( max_SecZ , theResult.GetSecondary( i )->GetParticle()->GetDefinition()->GetAtomicNumber() );
55      sum_A += theResult.GetSecondary( i )->GetParticle()->GetDefinition()->GetAtomicMass();
56      max_SecA = std::max ( max_SecA , theResult.GetSecondary( i )->GetParticle()->GetDefinition()->GetAtomicMass() );
57   }
58
59   G4ParticleDefinition* resi_pd = NULL;
60   if ( (int)(theBaseZ - sum_Z) == 0 && (int)(theBaseA + 1 - sum_A) == 0 )
61      resi_pd = G4ParticleTable::GetParticleTable()->GetIon ( max_SecZ , max_SecA , 0.0 );
62   else
63   {
64      resi_pd = G4ParticleTable::GetParticleTable()->GetIon ( int(theBaseZ - sum_Z) , (int)(theBaseA + 1 - sum_A) , 0.0 );
65      if ( resi_pd == NULL )
66      {
67        // theNDLDataZ,A has the Z and A of used NDL file
68        if ( (int)(theNDLDataZ - sum_Z) == 0 && (int)(theNDLDataA + 1 - sum_A) == 0 )
69        {
70           G4int dif_Z = ( int ) ( theNDLDataZ - theBaseZ );
71           G4int dif_A = ( int ) ( theNDLDataA - theBaseA );
72           resi_pd = G4ParticleTable::GetParticleTable()->GetIon ( max_SecZ - dif_Z , max_SecA - dif_A , 0.0 );
73           for ( int i = 0 ; i < nSecondaries ; i++ )
74           {
75              if ( theResult.GetSecondary( i )->GetParticle()->GetDefinition()->GetAtomicNumber() == max_SecZ 
76                && theResult.GetSecondary( i )->GetParticle()->GetDefinition()->GetAtomicMass() == max_SecA )
77              {
78                 G4ThreeVector p = theResult.GetSecondary( i )->GetParticle()->GetMomentum();
79                 p = p * resi_pd->GetPDGMass()/ G4ParticleTable::GetParticleTable()->GetIon ( max_SecZ , max_SecA , 0.0 )->GetPDGMass(); 
80                 theResult.GetSecondary( i )->GetParticle()->SetDefinition( resi_pd );
81                 theResult.GetSecondary( i )->GetParticle()->SetMomentum( p );
82              } 
83           }
84        }
85      }
86   }
87
88
89   G4LorentzVector secs_4p_lab( 0.0 );
90
91   G4int n_sec = theResult.GetNumberOfSecondaries();
92   G4double fast = 0;
93   G4double slow = 1;
94   G4int ifast = 0;
95   G4int islow = 0;
96   G4int ires = -1;
97
98   for ( G4int i = 0 ; i < n_sec ; i++ )
99   {
100
101      //G4cout << "HP_DB " << i
102      //       << " " << theResult.GetSecondary( i )->GetParticle()->GetDefinition()->GetParticleName()
103      //       << " 4p " << theResult.GetSecondary( i )->GetParticle()->Get4Momentum()
104      //       << " ke " << theResult.GetSecondary( i )->GetParticle()->Get4Momentum().e() - theResult.GetSecondary( i )->GetParticle()->GetDefinition()->GetPDGMass()
105      //       << G4endl;
106
107      secs_4p_lab += theResult.GetSecondary( i )->GetParticle()->Get4Momentum();
108
109      G4double beta = 0;
110      if ( theResult.GetSecondary( i )->GetParticle()->GetDefinition() != G4Gamma::Gamma() )
111      {
112         beta = theResult.GetSecondary( i )->GetParticle()->Get4Momentum().beta(); 
113      }
114      else
115      {
116         beta = 1;
117      }
118
119      if ( theResult.GetSecondary( i )->GetParticle()->GetDefinition() == resi_pd ) ires = i; 
120
121      if ( slow > beta && beta != 0 ) 
122      {
123         slow = beta; 
124         islow = i;
125      }
126
127      if ( fast <= beta ) 
128      {
129         if ( fast != 1 )
130         {
131            fast = beta; 
132            ifast = i;
133         } 
134         else
135         {
136//          fast is already photon then check E
137            G4double e = theResult.GetSecondary( i )->GetParticle()->Get4Momentum().e();
138            if ( e > theResult.GetSecondary( ifast )->GetParticle()->Get4Momentum().e() )   
139            {
140//             among photons, the highest E becomes the fastest
141               ifast = i; 
142            }
143         } 
144      }
145   }
146
147
148   G4LorentzVector dif_4p = init_4p_lab - secs_4p_lab;
149
150   //G4cout << "HP_DB dif_4p " << init_4p_lab - secs_4p_lab << G4endl;
151   //G4cout << "HP_DB dif_3p mag " << ( dif_4p.v() ).mag() << G4endl;
152   //G4cout << "HP_DB dif_e " << dif_4p.e() - ( dif_4p.v() ).mag()<< G4endl;
153
154   G4LorentzVector p4(0);
155   if ( ires == -1 ) 
156   {
157//    Create and Add Residual Nucleus
158      ires = nSecondaries;
159      nSecondaries += 1;
160
161      G4DynamicParticle* res = new G4DynamicParticle ( resi_pd , dif_4p.v() );   
162      theResult.AddSecondary ( res );   
163
164      p4 = res->Get4Momentum(); 
165      if ( slow > p4.beta() ) 
166      {
167         slow = p4.beta(); 
168         islow = ires;
169      }
170
171      dif_4p = init_4p_lab - ( secs_4p_lab + p4 ); 
172   }
173
174   //Which is bigger dif_p or dif_e
175
176   if ( dif_4p.v().mag() < std::abs( dif_4p.e() ) )
177   {
178
179      // Adjust p
180      //if ( dif_4p.v().mag() < 1*MeV )
181      if ( minimum_energy < dif_4p.v().mag() && dif_4p.v().mag() < 1*MeV )
182      {
183
184         nSecondaries += 1;
185         theResult.AddSecondary ( new G4DynamicParticle ( G4Gamma::Gamma() , dif_4p.v() ) );   
186
187      }
188      else
189      {
190         //G4cout << "HP_DB Difference in dif_p is too large (>1MeV) or too small(<1keV) to adjust, so that give up tuning" << G4endl;
191      }
192
193   }
194   else
195   {
196
197      // dif_p > dif_e
198      // at first momentum
199      // Move residual momentum
200
201      p4 = theResult.GetSecondary( ires )->GetParticle()->Get4Momentum();
202      theResult.GetSecondary( ires )->GetParticle()->SetMomentum( p4.v() + dif_4p.v() ); 
203      dif_4p = init_4p_lab - ( secs_4p_lab - p4 + theResult.GetSecondary( ires )->GetParticle()->Get4Momentum()  ); 
204
205      //G4cout << "HP_DB new residual kinetic energy " << theResult.GetSecondary( ires )->GetParticle()->GetKineticEnergy() << G4endl;
206
207   }
208
209   G4double dif_e = dif_4p.e() - ( dif_4p.v() ).mag();
210   //G4cout << "HP_DB dif_e " << dif_e << G4endl;
211
212   if ( dif_e > 0 )
213   {
214
215//    create 2 gamma
216
217      nSecondaries += 2;
218      G4double e1 = ( dif_4p.e() -dif_4p.v().mag() ) / 2;
219     
220      if ( minimum_energy < e1 ) 
221      {
222         G4double costh = 2.*G4UniformRand()-1.;
223         G4double phi = twopi*G4UniformRand();
224         G4ThreeVector dir( std::sin(std::acos(costh))*std::cos(phi), 
225                            std::sin(std::acos(costh))*std::sin(phi),
226                            costh);
227         theResult.AddSecondary ( new G4DynamicParticle ( G4Gamma::Gamma() , e1*dir ) );   
228         theResult.AddSecondary ( new G4DynamicParticle ( G4Gamma::Gamma() , -e1*dir ) );   
229      }
230      else
231      {
232         //G4cout << "HP_DB Difference is too small(<1keV) to adjust, so that neglect it" << G4endl;
233      }
234
235   }
236   else    //dif_e < 0
237   {
238
239//    At first reduce KE of the fastest secondary;
240      G4double ke0 = theResult.GetSecondary( ifast )->GetParticle()->GetKineticEnergy();
241      G4ThreeVector p0 = theResult.GetSecondary( ifast )->GetParticle()->GetMomentum();
242      G4ThreeVector dir = ( theResult.GetSecondary( ifast )->GetParticle()->GetMomentum() ).unit();
243
244      //G4cout << "HP_DB ifast " << ifast << " ke0 " << ke0 << G4endl;
245
246      if ( ke0 + dif_e > 0 )
247      {
248         theResult.GetSecondary( ifast )->GetParticle()->SetKineticEnergy( ke0 + dif_e ); 
249         G4ThreeVector dp = p0 - theResult.GetSecondary( ifast )->GetParticle()->GetMomentum(); 
250
251         G4ThreeVector p = theResult.GetSecondary( islow )->GetParticle()->GetMomentum();
252         //theResult.GetSecondary( islow )->GetParticle()->SetMomentum( p - dif_e*dir );
253         theResult.GetSecondary( islow )->GetParticle()->SetMomentum( p + dp ); 
254      }
255      else
256      {
257         //G4cout << "HP_DB Difference in dif_e too large ( <0MeV ) to adjust, so that give up tuning" << G4endl;
258      }
259
260   }
261
262}
Note: See TracBrowser for help on using the repository browser.