source: trunk/source/processes/hadronic/models/theo_high_energy/src/G4QuasiElasticChannel.cc@ 1199

Last change on this file since 1199 was 1196, checked in by garnier, 16 years ago

update CVS release candidate geant4.9.3.01
File size: 8.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//
27// $Id: G4QuasiElasticChannel.cc,v 1.7 2009/04/09 08:28:42 mkossov Exp $
28// GEANT4 tag $Name: geant4-09-03-cand-01 $
29//
30
31// Author : Gunter Folger March 2007
32// Modified by Mikhail Kossov. Apr2009, E/M conservation: ResidualNucleus is added (ResNuc)
33// Class Description
34// Final state production model for theoretical models of hadron inelastic
35// quasi elastic scattering in geant4;
36// Class Description - End
37//
38// Modified:
39
40#include "G4QuasiElasticChannel.hh"
41
42#include "G4Fancy3DNucleus.hh"
43
44
45#include "G4HadTmpUtil.hh" //lrint
46
47//#define debug_scatter
48
49G4QuasiElasticChannel::G4QuasiElasticChannel()
50{
51 theQuasiElastic=G4QuasiFreeRatios::GetPointer();
52}
53
54G4QuasiElasticChannel::~G4QuasiElasticChannel()
55{}
56
57G4double G4QuasiElasticChannel::GetFraction(G4Nucleus &theNucleus,
58 const G4DynamicParticle & thePrimary)
59{
60 std::pair<G4double,G4double> ratios;
61#ifdef debug_scatter
62 G4cout << "G4QuasiElasticChannel:: P=" << thePrimary.GetTotalMomentum()
63 << ", pPDG=" << thePrimary.GetDefinition()->GetPDGEncoding()
64 << ", Z = " << G4lrint(theNucleus.GetZ())
65 << ", N = " << G4lrint(theNucleus.GetN()-theNucleus.GetZ()) << G4endl;
66#endif
67 ratios=theQuasiElastic->GetRatios(thePrimary.GetTotalMomentum(),
68 thePrimary.GetDefinition()->GetPDGEncoding(),
69 G4lrint(theNucleus.GetZ()),
70 G4lrint(theNucleus.GetN()-theNucleus.GetZ()));
71#ifdef debug_scatter
72 G4cout << "G4QuasiElasticChannel::ratios " << ratios.first << " x " <<ratios.second
73 << " = " << ratios.first*ratios.second << G4endl;
74#endif
75
76 return ratios.first*ratios.second;
77 //return 0.; // Switch off quasielastic (temporary) M.K.
78 //return 1.; // Only quasielastic (temporary) M.K. (DANGEROSE! Crashes at A=1!)
79}
80
81G4KineticTrackVector * G4QuasiElasticChannel::Scatter(G4Nucleus &theNucleus,
82 const G4DynamicParticle & thePrimary)
83{
84 G4int A=G4lrint(theNucleus.GetN());
85 G4int Z=G4lrint(theNucleus.GetZ()); // M.K. ResNuc
86 // build Nucleus and choose random nucleon to scatter with
87 the3DNucleus.Init(theNucleus.GetN(),theNucleus.GetZ());
88#ifdef debug_scatter
89 G4cout<<"G4QElC::Scat: Before GetNucleons " << G4endl;
90#endif
91 const std::vector<G4Nucleon *> nucleons=the3DNucleus.GetNucleons();
92 G4double targetNucleusMass=the3DNucleus.GetMass(); // M.K. ResNuc
93#ifdef debug_scatter
94 G4cout<<"G4QElC::Scat: targetMass = " << targetNucleusMass << G4endl;
95#endif
96 G4LorentzVector targetNucleus4Mom(0.,0.,0.,targetNucleusMass); // M.K. ResNuc
97 G4int index;
98 do
99 {
100 index=G4lrint((A-1)*G4UniformRand());
101 } while (index < 0 || index >= static_cast<G4int>(nucleons.size()));
102#ifdef debug_scatter
103 G4cout<<"G4QElC::Scat: index = " << index << G4endl;
104#endif
105 G4ParticleDefinition * pDef= nucleons[index]->GetDefinition();
106 G4int resA=A-1; // M.K. ResNuc
107 G4int resZ=Z-static_cast<int>(pDef->GetPDGCharge()); // M.K. ResNuc
108 G4ParticleDefinition* resDef=G4Neutron::Neutron(); // Resolve t-p=nn problem M.K. ResNuc
109 G4double residualNucleusMass=resDef->GetPDGMass(); // M.K. ResNuc
110 if(resZ)
111 {
112 resDef=G4ParticleTable::GetParticleTable()->FindIon(resZ,resA,0,resZ);// M.K. ResNuc
113 residualNucleusMass=resDef->GetPDGMass(); // M.K. ResNuc
114 }
115 else residualNucleusMass*=resA; // resA=resN M.K. ResNuc
116#ifdef debug_scatter
117 G4cout<<"G4QElChan::Scatter: neutron - proton? A ="<<A<<", Z="<<Z<<", projName="
118 <<pDef->GetParticleName()<<G4endl;
119#endif
120 // G4LorentzVector pNucleon(G4ThreeVector(0,0,0),pDef->GetPDGMass());
121 G4LorentzVector pNucleon=nucleons[index]->Get4Momentum();
122 G4double residualNucleusEnergy=std::sqrt(residualNucleusMass*residualNucleusMass+
123 pNucleon.vect().mag2()); // M.K. ResNuc
124 pNucleon.setE(targetNucleusMass-residualNucleusEnergy); // M.K. ResNuc
125 G4LorentzVector residualNucleus4Mom=targetNucleus4Mom-pNucleon; // M.K. ResNuc
126
127 std::pair<G4LorentzVector,G4LorentzVector> result;
128
129 result=theQuasiElastic->Scatter(pDef->GetPDGEncoding(),pNucleon,
130 thePrimary.GetDefinition()->GetPDGEncoding(),
131 thePrimary.Get4Momentum());
132 G4LorentzVector scatteredHadron4Mom=result.second; // M.K. ResNuc
133 if (result.first.e() <= 0.)
134 {
135 //G4cout << "Warning - G4QuasiElasticChannel::Scatter no scattering" << G4endl;
136 //return 0; //no scatter
137 G4LorentzVector scatteredHadron4Mom=thePrimary.Get4Momentum(); // M.K. ResNuc
138 residualNucleus4Mom=G4LorentzVector(0.,0.,0.,targetNucleusMass); // M.K. ResNuc
139 resDef=G4ParticleTable::GetParticleTable()->FindIon(Z,A,0,Z); // M.K. ResNuc
140 }
141
142#ifdef debug_scatter
143 G4LorentzVector EpConservation=pNucleon+thePrimary.Get4Momentum()
144 - result.first - result.second;
145 if ( (EpConservation.vect().mag2() > .01*MeV*MeV )
146 || (std::abs(EpConservation.e()) > 0.1 * MeV ) )
147 {
148 G4cout << "Warning - G4QuasiElasticChannel::Scatter E-p non conservation : "
149 << EpConservation << G4endl;
150 }
151#endif
152
153 G4KineticTrackVector * ktv;
154 ktv=new G4KineticTrackVector();
155 G4KineticTrack * sPrim=new G4KineticTrack(thePrimary.GetDefinition(),
156 0.,G4ThreeVector(0), scatteredHadron4Mom);
157 ktv->push_back(sPrim);
158 if (result.first.e() > 0.)
159 {
160 G4KineticTrack * sNuc=new G4KineticTrack(pDef, 0.,G4ThreeVector(0), result.first);
161 ktv->push_back(sNuc);
162 }
163 if(resZ || resA==1) // For the only neutron or for tnuclei with Z>0 M.K. ResNuc
164 {
165 G4KineticTrack * rNuc=new G4KineticTrack(resDef,
166 0.,G4ThreeVector(0), residualNucleus4Mom); // M.K. ResNuc
167 ktv->push_back(rNuc); // M.K. ResNuc
168 }
169 else // The residual nucleus consists of only neutrons M.K. ResNuc
170 {
171 residualNucleus4Mom/=resA; // Split 4-mom of A*n system equally M.K. ResNuc
172 for(G4int in=0; in<resA; in++) // Loop over neutrons in A*n system. M.K. ResNuc
173 {
174 G4KineticTrack* rNuc=new G4KineticTrack(resDef,
175 0.,G4ThreeVector(0), residualNucleus4Mom); // M.K. ResNuc
176 ktv->push_back(rNuc); // M.K. ResNuc
177 }
178 }
179#ifdef debug_scatter
180 G4cout<<"G4QElC::Scat: Nucleon: "<<result.first <<" mass "<<result.first.mag() << G4endl;
181 G4cout<<"G4QElC::Scat: Project: "<<result.second<<" mass "<<result.second.mag()<< G4endl;
182#endif
183 return ktv;
184}
Note: See TracBrowser for help on using the repository browser.