| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // 26 //
>> 27 // $Id: G4QuasiElasticChannel.cc,v 1.9 2010-09-17 11:34:29 gunter Exp $
>> 28 // GEANT4 tag $Name: not supported by cvs2svn $
27 // 29 //
28 30
29 // Author : Gunter Folger March 2007 31 // Author : Gunter Folger March 2007
30 // Modified by Mikhail Kossov. Apr2009, E/M co 32 // Modified by Mikhail Kossov. Apr2009, E/M conservation: ResidualNucleus is added (ResNuc)
31 // Class Description 33 // Class Description
32 // Final state production model for theoretica 34 // Final state production model for theoretical models of hadron inelastic
33 // quasi elastic scattering in geant4; 35 // quasi elastic scattering in geant4;
34 // Class Description - End 36 // Class Description - End
35 // 37 //
36 // Modified: 38 // Modified:
37 // 20110805 M. Kelsey -- Follow change to G4V 39 // 20110805 M. Kelsey -- Follow change to G4V3DNucleus::GetNucleons()
38 // 20110808 M. Kelsey -- Move #includes from 40 // 20110808 M. Kelsey -- Move #includes from .hh, add many missing
39 41
40 #include "G4QuasiElasticChannel.hh" 42 #include "G4QuasiElasticChannel.hh"
41 43
42 #include "G4Fancy3DNucleus.hh" 44 #include "G4Fancy3DNucleus.hh"
43 #include "G4DynamicParticle.hh" 45 #include "G4DynamicParticle.hh"
44 #include "G4HadTmpUtil.hh" /* lrint */ 46 #include "G4HadTmpUtil.hh" /* lrint */
45 #include "G4KineticTrack.hh" 47 #include "G4KineticTrack.hh"
46 #include "G4KineticTrackVector.hh" 48 #include "G4KineticTrackVector.hh"
47 #include "G4LorentzVector.hh" 49 #include "G4LorentzVector.hh"
48 #include "G4Neutron.hh" 50 #include "G4Neutron.hh"
49 #include "G4Nucleon.hh" 51 #include "G4Nucleon.hh"
50 #include "G4Nucleus.hh" 52 #include "G4Nucleus.hh"
51 #include "G4ParticleDefinition.hh" 53 #include "G4ParticleDefinition.hh"
52 #include "G4ParticleTable.hh" 54 #include "G4ParticleTable.hh"
53 #include "G4IonTable.hh" <<
54 #include "G4QuasiElRatios.hh" 55 #include "G4QuasiElRatios.hh"
55 #include "globals.hh" 56 #include "globals.hh"
56 #include <vector> 57 #include <vector>
57 #include "G4PhysicsModelCatalog.hh" <<
58 58
59 //#define debug_scatter 59 //#define debug_scatter
60 60
61 61
62 G4QuasiElasticChannel::G4QuasiElasticChannel() 62 G4QuasiElasticChannel::G4QuasiElasticChannel()
63 : G4HadronicInteraction("QuasiElastic"), << 63 : theQuasiElastic(G4QuasiElRatios::GetPointer()),
64 theQuasiElastic(new G4QuasiElRatios), << 64 the3DNucleus(new G4Fancy3DNucleus) {}
65 the3DNucleus(new G4Fancy3DNucleus), <<
66 secID(-1) { <<
67 secID = G4PhysicsModelCatalog::GetModelID( " <<
68 } <<
69 65
70 G4QuasiElasticChannel::~G4QuasiElasticChannel( 66 G4QuasiElasticChannel::~G4QuasiElasticChannel()
71 { 67 {
72 delete the3DNucleus; 68 delete the3DNucleus;
73 delete theQuasiElastic; <<
74 } 69 }
75 70
76 G4double G4QuasiElasticChannel::GetFraction(G4 71 G4double G4QuasiElasticChannel::GetFraction(G4Nucleus &theNucleus,
77 const G4DynamicParticle & thePrimary) 72 const G4DynamicParticle & thePrimary)
78 { 73 {
79 #ifdef debug_scatter 74 #ifdef debug_scatter
80 G4cout << "G4QuasiElasticChannel:: P=" < 75 G4cout << "G4QuasiElasticChannel:: P=" << thePrimary.GetTotalMomentum()
81 << ", pPDG=" << thePrimary.GetDef 76 << ", pPDG=" << thePrimary.GetDefinition()->GetPDGEncoding()
82 << ", Z = " << theNucleus.GetZ_a 77 << ", Z = " << theNucleus.GetZ_asInt())
83 << ", N = " << theNucleus.GetN_a 78 << ", N = " << theNucleus.GetN_asInt())
84 << ", A = " << theNucleus.GetA_a 79 << ", A = " << theNucleus.GetA_asInt() << G4endl;
85 #endif 80 #endif
86 81
87 std::pair<G4double,G4double> ratios; 82 std::pair<G4double,G4double> ratios;
88 ratios=theQuasiElastic->GetRatios(thePrimary 83 ratios=theQuasiElastic->GetRatios(thePrimary.GetTotalMomentum(),
89 thePrimary 84 thePrimary.GetDefinition()->GetPDGEncoding(),
90 theNucleus 85 theNucleus.GetZ_asInt(),
91 theNucleus 86 theNucleus.GetN_asInt());
92 #ifdef debug_scatter 87 #ifdef debug_scatter
93 G4cout << "G4QuasiElasticChannel::ratios 88 G4cout << "G4QuasiElasticChannel::ratios " << ratios.first << " x " <<ratios.second
94 << " = " << ratios.first*ratios. 89 << " = " << ratios.first*ratios.second << G4endl;
95 #endif 90 #endif
96 91
97 return ratios.first*ratios.second; 92 return ratios.first*ratios.second;
98 } 93 }
99 94
100 G4KineticTrackVector * G4QuasiElasticChannel:: 95 G4KineticTrackVector * G4QuasiElasticChannel::Scatter(G4Nucleus &theNucleus,
101 96 const G4DynamicParticle & thePrimary)
102 { 97 {
103 G4int A=theNucleus.GetA_asInt(); 98 G4int A=theNucleus.GetA_asInt();
104 G4int Z=theNucleus.GetZ_asInt(); 99 G4int Z=theNucleus.GetZ_asInt();
105 // build Nucleus and choose random nucleon 100 // build Nucleus and choose random nucleon to scatter with
106 the3DNucleus->Init(theNucleus.GetA_asInt(),t 101 the3DNucleus->Init(theNucleus.GetA_asInt(),theNucleus.GetZ_asInt());
107 const std::vector<G4Nucleon>& nucleons=the3D 102 const std::vector<G4Nucleon>& nucleons=the3DNucleus->GetNucleons();
108 G4double targetNucleusMass=the3DNucleus->Get << 103 G4double targetNucleusMass=the3DNucleus->GetMass(); // M.K. ResNuc
109 G4LorentzVector targetNucleus4Mom(0.,0.,0.,t << 104 G4LorentzVector targetNucleus4Mom(0.,0.,0.,targetNucleusMass); // M.K. ResNuc
110 G4int index; 105 G4int index;
111 do { << 106 do
>> 107 {
112 index=G4lrint((A-1)*G4UniformRand()); 108 index=G4lrint((A-1)*G4UniformRand());
113 } while (index < 0 || index >= static_cast<G << 109 } while (index < 0 || index >= static_cast<G4int>(nucleons.size()));
114 << 110 G4ParticleDefinition * pDef= nucleons[index].GetDefinition();
115 const G4ParticleDefinition * pDef= nucleons[ <<
116 111
117 G4int resA=A - 1; << 112 G4int resA=A-1; // M.K. ResNuc
118 G4int resZ=Z - static_cast<int>(pDef->GetPDG << 113 G4int resZ=Z-static_cast<int>(pDef->GetPDGCharge()); // M.K. ResNuc
119 const G4ParticleDefinition* resDef; << 114 G4ParticleDefinition* resDef=G4Neutron::Neutron(); // Resolve t-p=nn problem M.K. ResNuc
120 G4double residualNucleusMass; << 115 G4double residualNucleusMass=resDef->GetPDGMass(); // M.K. ResNuc
121 if(resZ) 116 if(resZ)
122 { 117 {
123 resDef=G4ParticleTable::GetParticleTable() << 118 resDef=G4ParticleTable::GetParticleTable()->FindIon(resZ,resA,0,resZ);// M.K. ResNuc
124 residualNucleusMass=resDef->GetPDGMass(); << 119 residualNucleusMass=resDef->GetPDGMass(); // M.K. ResNuc
125 } 120 }
126 else { 121 else {
127 resDef=G4Neutron::Neutron(); << 122 residualNucleusMass*=resA;
128 residualNucleusMass=resA * G4Neutron::Neut <<
129 } 123 }
130 #ifdef debug_scatter 124 #ifdef debug_scatter
131 G4cout<<"G4QElChan::Scatter: neutron - pr 125 G4cout<<"G4QElChan::Scatter: neutron - proton? A ="<<A<<", Z="<<Z<<", projName="
132 <<pDef->GetParticleName()<<G4endl; 126 <<pDef->GetParticleName()<<G4endl;
133 #endif 127 #endif
134 128
135 G4LorentzVector pNucleon=nucleons[index].Get 129 G4LorentzVector pNucleon=nucleons[index].Get4Momentum();
136 G4double residualNucleusEnergy=std::sqrt(sqr << 130 G4double residualNucleusEnergy=std::sqrt(residualNucleusMass*residualNucleusMass+
137 pNu << 131 pNucleon.vect().mag2()); // M.K. ResNuc
138 pNucleon.setE(targetNucleusMass-residualNucl << 132 pNucleon.setE(targetNucleusMass-residualNucleusEnergy); // M.K. ResNuc
139 G4LorentzVector residualNucleus4Mom=targetNu << 133 G4LorentzVector residualNucleus4Mom=targetNucleus4Mom-pNucleon; // M.K. ResNuc
140 134
141 std::pair<G4LorentzVector,G4LorentzVector> r 135 std::pair<G4LorentzVector,G4LorentzVector> result;
142 136
143 result=theQuasiElastic->Scatter(pDef->GetPDG 137 result=theQuasiElastic->Scatter(pDef->GetPDGEncoding(),pNucleon,
144 thePrimary.G 138 thePrimary.GetDefinition()->GetPDGEncoding(),
145 thePrimary.G << 139 thePrimary.Get4Momentum());
146 G4LorentzVector scatteredHadron4Mom; << 140 G4LorentzVector scatteredHadron4Mom=result.second; // M.K. ResNuc
147 if (result.first.e() > 0.) << 141 if (result.first.e() <= 0.)
148 scatteredHadron4Mom=result.second; << 142 {
149 else { //scatter failed <<
150 //G4cout << "Warning - G4QuasiElasticChann 143 //G4cout << "Warning - G4QuasiElasticChannel::Scatter no scattering" << G4endl;
151 //return 0; //no scatter 144 //return 0; //no scatter
152 scatteredHadron4Mom=thePrimary.Get4Momentu << 145 G4LorentzVector scatteredHadron4Mom=thePrimary.Get4Momentum(); // M.K. ResNuc
153 residualNucleus4Mom=G4LorentzVector(0.,0., << 146 residualNucleus4Mom=G4LorentzVector(0.,0.,0.,targetNucleusMass); // M.K. ResNuc
154 resDef=G4ParticleTable::GetParticleTable() << 147 resDef=G4ParticleTable::GetParticleTable()->FindIon(Z,A,0,Z); // M.K. ResNuc
155 } 148 }
156 149
157 #ifdef debug_scatter 150 #ifdef debug_scatter
158 G4LorentzVector EpConservation=pNucleon+theP 151 G4LorentzVector EpConservation=pNucleon+thePrimary.Get4Momentum()
159 - result.firs 152 - result.first - result.second;
160 if ( (EpConservation.vect().mag2() > .01*M 153 if ( (EpConservation.vect().mag2() > .01*MeV*MeV )
161 || (std::abs(EpConservation.e()) > 0.1 * 154 || (std::abs(EpConservation.e()) > 0.1 * MeV ) )
162 { 155 {
163 G4cout << "Warning - G4QuasiElasticChannel 156 G4cout << "Warning - G4QuasiElasticChannel::Scatter E-p non conservation : "
164 << EpConservation << G4endl; 157 << EpConservation << G4endl;
165 } 158 }
166 #endif 159 #endif
167 160
168 G4KineticTrackVector * ktv = new G4KineticTr << 161 G4KineticTrackVector * ktv;
>> 162 ktv=new G4KineticTrackVector();
169 G4KineticTrack * sPrim=new G4KineticTrack(th 163 G4KineticTrack * sPrim=new G4KineticTrack(thePrimary.GetDefinition(),
170 0. 164 0.,G4ThreeVector(0), scatteredHadron4Mom);
171 sPrim->SetCreatorModelID( secID ); <<
172 ktv->push_back(sPrim); 165 ktv->push_back(sPrim);
173 if (result.first.e() > 0.) 166 if (result.first.e() > 0.)
174 { 167 {
175 G4KineticTrack * sNuc=new G4KineticTrack(p 168 G4KineticTrack * sNuc=new G4KineticTrack(pDef, 0.,G4ThreeVector(0), result.first);
176 sNuc->SetCreatorModelID( secID ); <<
177 ktv->push_back(sNuc); 169 ktv->push_back(sNuc);
178 } 170 }
179 if(resZ || resA==1) // For the only neutron << 171 if(resZ || resA==1) // For the only neutron or for tnuclei with Z>0 M.K. ResNuc
180 { 172 {
181 G4KineticTrack * rNuc=new G4KineticTrack(r 173 G4KineticTrack * rNuc=new G4KineticTrack(resDef,
182 0.,G4ThreeVector(0) << 174 0.,G4ThreeVector(0), residualNucleus4Mom); // M.K. ResNuc
183 rNuc->SetCreatorModelID( secID ); << 175 ktv->push_back(rNuc); // M.K. ResNuc
184 ktv->push_back(rNuc); <<
185 } 176 }
186 else // The residual nucleus consists of onl << 177 else // The residual nucleus consists of only neutrons M.K. ResNuc
187 { 178 {
188 residualNucleus4Mom/=resA; // Split 4- << 179 residualNucleus4Mom/=resA; // Split 4-mom of A*n system equally M.K. ResNuc
189 for(G4int in=0; in<resA; in++) // Loop ove << 180 for(G4int in=0; in<resA; in++) // Loop over neutrons in A*n system. M.K. ResNuc
190 { 181 {
191 G4KineticTrack* rNuc=new G4KineticTrack( 182 G4KineticTrack* rNuc=new G4KineticTrack(resDef,
192 0.,G4ThreeVector(0) << 183 0.,G4ThreeVector(0), residualNucleus4Mom); // M.K. ResNuc
193 rNuc->SetCreatorModelID( secID ); << 184 ktv->push_back(rNuc); // M.K. ResNuc
194 ktv->push_back(rNuc); <<
195 } 185 }
196 } 186 }
197 #ifdef debug_scatter 187 #ifdef debug_scatter
198 G4cout<<"G4QElC::Scat: Nucleon: "<<result.fi 188 G4cout<<"G4QElC::Scat: Nucleon: "<<result.first <<" mass "<<result.first.mag() << G4endl;
199 G4cout<<"G4QElC::Scat: Project: "<<result.se 189 G4cout<<"G4QElC::Scat: Project: "<<result.second<<" mass "<<result.second.mag()<< G4endl;
200 #endif 190 #endif
201 return ktv; 191 return ktv;
202 } 192 }
203 193