| 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 // $Id$
26 // 27 //
27 // G4LMsdGenerator 28 // G4LMsdGenerator
28 // 29 //
29 // 30 //
30 31
31 #include "G4DynamicParticle.hh" 32 #include "G4DynamicParticle.hh"
32 #include "G4LMsdGenerator.hh" 33 #include "G4LMsdGenerator.hh"
33 #include "G4ReactionProductVector.hh" 34 #include "G4ReactionProductVector.hh"
34 #include "G4ReactionProduct.hh" 35 #include "G4ReactionProduct.hh"
35 #include "G4IonTable.hh" 36 #include "G4IonTable.hh"
36 #include "G4NucleiProperties.hh" 37 #include "G4NucleiProperties.hh"
37 #include "G4ParticleDefinition.hh" 38 #include "G4ParticleDefinition.hh"
38 #include "G4HadFinalState.hh" 39 #include "G4HadFinalState.hh"
39 #include "G4KineticTrack.hh" 40 #include "G4KineticTrack.hh"
40 #include "G4DecayKineticTracks.hh" 41 #include "G4DecayKineticTracks.hh"
41 #include "G4KineticTrackVector.hh" 42 #include "G4KineticTrackVector.hh"
42 #include "G4Log.hh" 43 #include "G4Log.hh"
43 #include "G4PhysicsModelCatalog.hh" <<
44 44
45 45
46 G4LMsdGenerator::G4LMsdGenerator(const G4Strin 46 G4LMsdGenerator::G4LMsdGenerator(const G4String& name)
47 : G4HadronicInteraction(name), secID(-1) << 47 : G4HadronicInteraction(name)
48 48
49 { 49 {
50 fPDGencoding = 0; 50 fPDGencoding = 0;
51 secID = G4PhysicsModelCatalog::GetModelID( " << 51
52 <<
53 // theParticleChange = new G4HadFinalState; 52 // theParticleChange = new G4HadFinalState;
54 } 53 }
55 54
56 G4LMsdGenerator::~G4LMsdGenerator() 55 G4LMsdGenerator::~G4LMsdGenerator()
57 { 56 {
58 // delete theParticleChange; 57 // delete theParticleChange;
59 } 58 }
60 59
61 void G4LMsdGenerator::ModelDescription(std::os 60 void G4LMsdGenerator::ModelDescription(std::ostream& outFile) const
62 { 61 {
63 outFile << GetModelName() <<" consists of a 62 outFile << GetModelName() <<" consists of a "
64 << " string model and a stage to de-excite 63 << " string model and a stage to de-excite the excited nuclear fragment."
65 << "\n<p>" 64 << "\n<p>"
66 << "The string model simulates the interac 65 << "The string model simulates the interaction of\n"
67 << "an incident hadron with a nucleu 66 << "an incident hadron with a nucleus, forming \n"
68 << "excited strings, decays these st 67 << "excited strings, decays these strings into hadrons,\n"
69 << "and leaves an excited nucleus. \ 68 << "and leaves an excited nucleus. \n"
70 << "<p>The string model:\n"; 69 << "<p>The string model:\n";
71 } 70 }
72 71
73 ////////////////////////////////////////////// 72 /////////////////////////////////////////////////////////////////
74 // 73 //
75 // Particle and kinematical limitation od diff 74 // Particle and kinematical limitation od diffraction dissociation
76 75
77 G4bool 76 G4bool
78 G4LMsdGenerator::IsApplicable( const G4HadProj 77 G4LMsdGenerator::IsApplicable( const G4HadProjectile& aTrack,
79 G4Nucleu 78 G4Nucleus& targetNucleus )
80 { 79 {
81 G4bool applied = false; 80 G4bool applied = false;
82 81
83 if( ( aTrack.GetDefinition() == G4Proton::Pr 82 if( ( aTrack.GetDefinition() == G4Proton::Proton() ||
84 aTrack.GetDefinition() == G4Neutron::Neutron 83 aTrack.GetDefinition() == G4Neutron::Neutron() ) &&
85 targetNucleus.GetA_asInt() >= 1 && 84 targetNucleus.GetA_asInt() >= 1 &&
86 // aTrack.GetKineticEnergy() > 1800*CLHE 85 // aTrack.GetKineticEnergy() > 1800*CLHEP::MeV
87 aTrack.GetKineticEnergy() > 300*CLHEP: 86 aTrack.GetKineticEnergy() > 300*CLHEP::MeV
88 ) // 750*CLHEP::MeV ) 87 ) // 750*CLHEP::MeV )
89 { 88 {
90 applied = true; 89 applied = true;
91 } 90 }
92 else if( ( aTrack.GetDefinition() == G4PionP 91 else if( ( aTrack.GetDefinition() == G4PionPlus::PionPlus() ||
93 aTrack.GetDefinition() == G4PionMinus:: 92 aTrack.GetDefinition() == G4PionMinus::PionMinus() ) &&
94 targetNucleus.GetA_asInt() >= 1 & 93 targetNucleus.GetA_asInt() >= 1 &&
95 aTrack.GetKineticEnergy() > 2340* 94 aTrack.GetKineticEnergy() > 2340*CLHEP::MeV )
96 { 95 {
97 applied = true; 96 applied = true;
98 } 97 }
99 else if( ( aTrack.GetDefinition() == G4KaonP 98 else if( ( aTrack.GetDefinition() == G4KaonPlus::KaonPlus() ||
100 aTrack.GetDefinition() == G4KaonMinus:: 99 aTrack.GetDefinition() == G4KaonMinus::KaonMinus() ) &&
101 targetNucleus.GetA_asInt() >= 1 & 100 targetNucleus.GetA_asInt() >= 1 &&
102 aTrack.GetKineticEnergy() > 1980* 101 aTrack.GetKineticEnergy() > 1980*CLHEP::MeV )
103 { 102 {
104 applied = true; 103 applied = true;
105 } 104 }
106 return applied; 105 return applied;
107 } 106 }
108 107
109 ////////////////////////////////////////////// 108 /////////////////////////////////////////////////////////////////
110 // 109 //
111 // Return dissociated particle products and re 110 // Return dissociated particle products and recoil nucleus
112 111
113 G4HadFinalState* 112 G4HadFinalState*
114 G4LMsdGenerator::ApplyYourself( const G4HadPro 113 G4LMsdGenerator::ApplyYourself( const G4HadProjectile& aTrack,
115 G4Nucleu 114 G4Nucleus& targetNucleus )
116 { 115 {
117 theParticleChange.Clear(); 116 theParticleChange.Clear();
118 117
119 const G4HadProjectile* aParticle = &aTrack; 118 const G4HadProjectile* aParticle = &aTrack;
120 G4double eTkin = aParticle->GetKineticEnergy 119 G4double eTkin = aParticle->GetKineticEnergy();
121 120
122 if( eTkin <= 1.*CLHEP::GeV && aTrack.GetDefi 121 if( eTkin <= 1.*CLHEP::GeV && aTrack.GetDefinition() != G4Proton::Proton())
123 { 122 {
124 theParticleChange.SetEnergyChange(eTkin); 123 theParticleChange.SetEnergyChange(eTkin);
125 theParticleChange.SetMomentumChange(aTrack 124 theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
126 return &theParticleChange; 125 return &theParticleChange;
127 } 126 }
128 127
129 G4int A = targetNucleus.GetA_asInt(); 128 G4int A = targetNucleus.GetA_asInt();
130 G4int Z = targetNucleus.GetZ_asInt(); 129 G4int Z = targetNucleus.GetZ_asInt();
131 130
132 G4double plab = aParticle->GetTotalMomentum( 131 G4double plab = aParticle->GetTotalMomentum();
133 G4double plab2 = plab*plab; 132 G4double plab2 = plab*plab;
134 133
135 const G4ParticleDefinition* theParticle = aP 134 const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
136 G4double partMass = theParticle->GetPDGMass( 135 G4double partMass = theParticle->GetPDGMass();
137 136
138 G4double oldE = partMass + eTkin; 137 G4double oldE = partMass + eTkin;
139 138
140 G4double targMass = G4NucleiProperties::Ge 139 G4double targMass = G4NucleiProperties::GetNuclearMass(A, Z);
141 G4double targMass2 = targMass*targMass; 140 G4double targMass2 = targMass*targMass;
142 141
143 G4LorentzVector partLV = aParticle->Get4Mome 142 G4LorentzVector partLV = aParticle->Get4Momentum();
144 143
145 G4double sumE = oldE + targMass; 144 G4double sumE = oldE + targMass;
146 G4double sumE2 = sumE*sumE; 145 G4double sumE2 = sumE*sumE;
147 146
148 G4ThreeVector p1 = partLV.vect(); 147 G4ThreeVector p1 = partLV.vect();
149 // G4cout<<"p1 = "<<p1<<G4endl; 148 // G4cout<<"p1 = "<<p1<<G4endl;
150 G4ParticleMomentum p1unit = p1.unit(); 149 G4ParticleMomentum p1unit = p1.unit();
151 150
152 G4double Mx = SampleMx(aParticle); // in GeV 151 G4double Mx = SampleMx(aParticle); // in GeV
153 G4double t = SampleT( aParticle, Mx); // in 152 G4double t = SampleT( aParticle, Mx); // in GeV
154 153
155 Mx *= CLHEP::GeV; 154 Mx *= CLHEP::GeV;
156 155
157 G4double Mx2 = Mx*Mx; 156 G4double Mx2 = Mx*Mx;
158 157
159 // equation for q|| based on sum-E-P and new 158 // equation for q|| based on sum-E-P and new invariant mass
160 159
161 G4double B = sumE2 + targMass2 - Mx2 - plab2 160 G4double B = sumE2 + targMass2 - Mx2 - plab2;
162 161
163 G4double a = 4*(plab2 - sumE2); 162 G4double a = 4*(plab2 - sumE2);
164 G4double b = 4*plab*B; 163 G4double b = 4*plab*B;
165 G4double c = B*B - 4*sumE2*targMass2; 164 G4double c = B*B - 4*sumE2*targMass2;
166 G4double det2 = b*b - 4*a*c; 165 G4double det2 = b*b - 4*a*c;
167 G4double qLong, det, eRetard; // , x2, x3, 166 G4double qLong, det, eRetard; // , x2, x3, e2;
168 167
169 if( det2 >= 0.) 168 if( det2 >= 0.)
170 { 169 {
171 det = std::sqrt(det2); 170 det = std::sqrt(det2);
172 qLong = (-b - det)/2./a; 171 qLong = (-b - det)/2./a;
173 eRetard = std::sqrt((plab-qLong)*(plab-qLo 172 eRetard = std::sqrt((plab-qLong)*(plab-qLong)+Mx2);
174 } 173 }
175 else 174 else
176 { 175 {
177 theParticleChange.SetEnergyChange(eTkin); 176 theParticleChange.SetEnergyChange(eTkin);
178 theParticleChange.SetMomentumChange(aTrack 177 theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
179 return &theParticleChange; 178 return &theParticleChange;
180 } 179 }
181 theParticleChange.SetStatusChange(stopAndKil 180 theParticleChange.SetStatusChange(stopAndKill);
182 181
183 plab -= qLong; 182 plab -= qLong;
184 183
185 G4ThreeVector pRetard = plab*p1unit; 184 G4ThreeVector pRetard = plab*p1unit;
186 185
187 G4ThreeVector pTarg = p1 - pRetard; 186 G4ThreeVector pTarg = p1 - pRetard;
188 187
189 G4double eTarg = std::sqrt( targMass2 + pTar 188 G4double eTarg = std::sqrt( targMass2 + pTarg.mag2()); // std::sqrt( targMass*targMass + pTarg.mag2() );
190 189
191 G4LorentzVector lvRetard(pRetard, eRetard); 190 G4LorentzVector lvRetard(pRetard, eRetard);
192 G4LorentzVector lvTarg(pTarg, eTarg); 191 G4LorentzVector lvTarg(pTarg, eTarg);
193 192
194 lvTarg += lvRetard; // sum LV 193 lvTarg += lvRetard; // sum LV
195 194
196 G4ThreeVector bst = lvTarg.boostVector(); 195 G4ThreeVector bst = lvTarg.boostVector();
197 196
198 lvRetard.boost(-bst); // to CNS 197 lvRetard.boost(-bst); // to CNS
199 198
200 G4ThreeVector pCMS = lvRetard.vect(); 199 G4ThreeVector pCMS = lvRetard.vect();
201 G4double momentumCMS = pCMS.mag(); 200 G4double momentumCMS = pCMS.mag();
202 G4double tMax = 4.0*momentumCMS*momen 201 G4double tMax = 4.0*momentumCMS*momentumCMS;
203 202
204 if( t > tMax ) t = tMax*G4UniformRand(); 203 if( t > tMax ) t = tMax*G4UniformRand();
205 204
206 G4double cost = 1. - 2.0*t/tMax; 205 G4double cost = 1. - 2.0*t/tMax;
207 206
208 207
209 G4double phi = G4UniformRand()*CLHEP::twopi 208 G4double phi = G4UniformRand()*CLHEP::twopi;
210 G4double sint; 209 G4double sint;
211 210
212 if( cost > 1.0 || cost < -1.0 ) // 211 if( cost > 1.0 || cost < -1.0 ) //
213 { 212 {
214 cost = 1.0; 213 cost = 1.0;
215 sint = 0.0; 214 sint = 0.0;
216 } 215 }
217 else // normal situation 216 else // normal situation
218 { 217 {
219 sint = std::sqrt( (1.0-cost)*(1.0+cost) ); 218 sint = std::sqrt( (1.0-cost)*(1.0+cost) );
220 } 219 }
221 G4ThreeVector v1( sint*std::cos(phi), sint*s 220 G4ThreeVector v1( sint*std::cos(phi), sint*std::sin(phi), cost);
222 221
223 v1 *= momentumCMS; 222 v1 *= momentumCMS;
224 223
225 G4LorentzVector lvRes( v1.x(),v1.y(),v1.z(), 224 G4LorentzVector lvRes( v1.x(),v1.y(),v1.z(), std::sqrt( momentumCMS*momentumCMS + Mx2));
226 225
227 lvRes.boost(bst); // to LS 226 lvRes.boost(bst); // to LS
228 227
229 lvTarg -= lvRes; 228 lvTarg -= lvRes;
230 229
231 G4double eRecoil = lvTarg.e() - targMass; 230 G4double eRecoil = lvTarg.e() - targMass;
232 231
233 if( eRecoil > 100.*CLHEP::MeV ) // add recoi 232 if( eRecoil > 100.*CLHEP::MeV ) // add recoil nucleus
234 { 233 {
235 G4ParticleDefinition * recoilDef = 0; 234 G4ParticleDefinition * recoilDef = 0;
236 235
237 if ( Z == 1 && A == 1 ) { recoilDef = 236 if ( Z == 1 && A == 1 ) { recoilDef = G4Proton::Proton(); }
238 else if ( Z == 1 && A == 2 ) { recoilDef = 237 else if ( Z == 1 && A == 2 ) { recoilDef = G4Deuteron::Deuteron(); }
239 else if ( Z == 1 && A == 3 ) { recoilDef = 238 else if ( Z == 1 && A == 3 ) { recoilDef = G4Triton::Triton(); }
240 else if ( Z == 2 && A == 3 ) { recoilDef = 239 else if ( Z == 2 && A == 3 ) { recoilDef = G4He3::He3(); }
241 else if ( Z == 2 && A == 4 ) { recoilDef = 240 else if ( Z == 2 && A == 4 ) { recoilDef = G4Alpha::Alpha(); }
242 else 241 else
243 { 242 {
244 recoilDef = 243 recoilDef =
245 G4ParticleTable::GetParticleTable()->GetIonT 244 G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon( Z, A, 0.0 );
246 } 245 }
247 G4DynamicParticle * aSec = new G4DynamicPa 246 G4DynamicParticle * aSec = new G4DynamicParticle( recoilDef, lvTarg);
248 theParticleChange.AddSecondary(aSec, secID << 247 theParticleChange.AddSecondary(aSec);
249 } 248 }
250 else if( eRecoil > 0.0 ) 249 else if( eRecoil > 0.0 )
251 { 250 {
252 theParticleChange.SetLocalEnergyDeposit( e 251 theParticleChange.SetLocalEnergyDeposit( eRecoil );
253 } 252 }
254 253
255 G4ParticleDefinition* ddPart = G4ParticleTab 254 G4ParticleDefinition* ddPart = G4ParticleTable::GetParticleTable()->
256 FindParticle( 255 FindParticle(fPDGencoding);
257 256
258 // G4cout<<fPDGencoding<<", "<<ddPart->GetPa 257 // G4cout<<fPDGencoding<<", "<<ddPart->GetParticleName()<<", "<<ddPart->GetPDGMass()<<" MeV; lvRes = "<<lvRes<<G4endl;
259 258
260 // G4DynamicParticle * aRes = new G4DynamicP 259 // G4DynamicParticle * aRes = new G4DynamicParticle( ddPart, lvRes);
261 // theParticleChange.AddSecondary(aRes); // 260 // theParticleChange.AddSecondary(aRes); // simply return resonance
262 261
263 262
264 263
265 // Recursive decay using methods of G4Kineti 264 // Recursive decay using methods of G4KineticTrack
266 265
267 G4KineticTrack ddkt( ddPart, 0., G4ThreeVect 266 G4KineticTrack ddkt( ddPart, 0., G4ThreeVector(0.,0.,0.), lvRes);
268 G4KineticTrackVector* ddktv = ddkt.Decay(); 267 G4KineticTrackVector* ddktv = ddkt.Decay();
269 268
270 G4DecayKineticTracks decay( ddktv ); 269 G4DecayKineticTracks decay( ddktv );
271 270
272 for( unsigned int i = 0; i < ddktv->size(); 271 for( unsigned int i = 0; i < ddktv->size(); i++ ) // add products to partchange
273 { 272 {
274 G4DynamicParticle * aNew = 273 G4DynamicParticle * aNew =
275 new G4DynamicParticle( ddktv->operator[] 274 new G4DynamicParticle( ddktv->operator[](i)->GetDefinition(),
276 ddktv->operator[] 275 ddktv->operator[](i)->Get4Momentum());
277 276
278 // G4cout<<" "<<i<<", "<<aNew->GetDe 277 // G4cout<<" "<<i<<", "<<aNew->GetDefinition()->GetParticleName()<<", "<<aNew->Get4Momentum()<<G4endl;
279 278
280 theParticleChange.AddSecondary(aNew, secID << 279 theParticleChange.AddSecondary(aNew);
281 delete ddktv->operator[](i); 280 delete ddktv->operator[](i);
282 } 281 }
283 delete ddktv; 282 delete ddktv;
284 283
285 return &theParticleChange; 284 return &theParticleChange;
286 } 285 }
287 286
288 ////////////////////////////////////// 287 //////////////////////////////////////
289 // 288 //
290 // Sample Mx as Roper resonances, set PDG enco 289 // Sample Mx as Roper resonances, set PDG encoding
291 290
292 G4double G4LMsdGenerator::SampleMx(const G4Had 291 G4double G4LMsdGenerator::SampleMx(const G4HadProjectile* aParticle)
293 { 292 {
294 G4double Mx = 0.; 293 G4double Mx = 0.;
295 G4int i; 294 G4int i;
296 G4double rand = G4UniformRand(); 295 G4double rand = G4UniformRand();
297 296
298 for( i = 0; i < 60; i++) 297 for( i = 0; i < 60; i++)
299 { 298 {
300 if( rand >= fProbMx[i][1] ) break; 299 if( rand >= fProbMx[i][1] ) break;
301 } 300 }
302 if(i <= 0) Mx = fProbMx[0][0]; 301 if(i <= 0) Mx = fProbMx[0][0];
303 else if(i >= 59) Mx = fProbMx[59][0]; 302 else if(i >= 59) Mx = fProbMx[59][0];
304 else Mx = fProbMx[i][0]; 303 else Mx = fProbMx[i][0];
305 304
306 fPDGencoding = 0; 305 fPDGencoding = 0;
307 306
308 if ( Mx <= 1.45 ) 307 if ( Mx <= 1.45 )
309 { 308 {
310 if( aParticle->GetDefinition() == G4Proton 309 if( aParticle->GetDefinition() == G4Proton::Proton() )
311 { 310 {
312 Mx = 1.44; 311 Mx = 1.44;
313 // fPDGencoding = 12212; 312 // fPDGencoding = 12212;
314 fPDGencoding = 2214; 313 fPDGencoding = 2214;
315 } 314 }
316 else if( aParticle->GetDefinition() == G4N 315 else if( aParticle->GetDefinition() == G4Neutron::Neutron() )
317 { 316 {
318 Mx = 1.44; 317 Mx = 1.44;
319 fPDGencoding = 12112; 318 fPDGencoding = 12112;
320 } 319 }
321 else if( aParticle->GetDefinition() == G4P 320 else if( aParticle->GetDefinition() == G4PionPlus::PionPlus() )
322 { 321 {
323 // Mx = 1.3; 322 // Mx = 1.3;
324 // fPDGencoding = 100211; 323 // fPDGencoding = 100211;
325 Mx = 1.26; 324 Mx = 1.26;
326 fPDGencoding = 20213; // a1(1260)+ 325 fPDGencoding = 20213; // a1(1260)+
327 } 326 }
328 else if( aParticle->GetDefinition() == G4P 327 else if( aParticle->GetDefinition() == G4PionMinus::PionMinus() )
329 { 328 {
330 // Mx = 1.3; 329 // Mx = 1.3;
331 // fPDGencoding = -100211; 330 // fPDGencoding = -100211;
332 Mx = 1.26; 331 Mx = 1.26;
333 fPDGencoding = -20213; // a1(1260)- 332 fPDGencoding = -20213; // a1(1260)-
334 } 333 }
335 else if( aParticle->GetDefinition() == G4K 334 else if( aParticle->GetDefinition() == G4KaonPlus::KaonPlus() )
336 { 335 {
337 Mx = 1.27; 336 Mx = 1.27;
338 fPDGencoding = 10323; 337 fPDGencoding = 10323;
339 } 338 }
340 else if( aParticle->GetDefinition() == G4K 339 else if( aParticle->GetDefinition() == G4KaonMinus::KaonMinus() )
341 { 340 {
342 Mx = 1.27; 341 Mx = 1.27;
343 fPDGencoding = -10323; 342 fPDGencoding = -10323;
344 } 343 }
345 } 344 }
346 else if ( Mx <= 1.55 ) 345 else if ( Mx <= 1.55 )
347 { 346 {
348 if( aParticle->GetDefinition() == G4Proton 347 if( aParticle->GetDefinition() == G4Proton::Proton() )
349 { 348 {
350 Mx = 1.52; 349 Mx = 1.52;
351 // fPDGencoding = 2124; 350 // fPDGencoding = 2124;
352 fPDGencoding = 2214; 351 fPDGencoding = 2214;
353 } 352 }
354 else if( aParticle->GetDefinition() == G4N 353 else if( aParticle->GetDefinition() == G4Neutron::Neutron() )
355 { 354 {
356 Mx = 1.52; 355 Mx = 1.52;
357 fPDGencoding = 1214; 356 fPDGencoding = 1214;
358 } 357 }
359 else if( aParticle->GetDefinition() == G4P 358 else if( aParticle->GetDefinition() == G4PionPlus::PionPlus() )
360 { 359 {
361 // Mx = 1.45; 360 // Mx = 1.45;
362 // fPDGencoding = 10211; 361 // fPDGencoding = 10211;
363 Mx = 1.32; 362 Mx = 1.32;
364 fPDGencoding = 215; // a2(1320)+ 363 fPDGencoding = 215; // a2(1320)+
365 } 364 }
366 else if( aParticle->GetDefinition() == G4P 365 else if( aParticle->GetDefinition() == G4PionMinus::PionMinus() )
367 { 366 {
368 // Mx = 1.45; 367 // Mx = 1.45;
369 // fPDGencoding = -10211; 368 // fPDGencoding = -10211;
370 Mx = 1.32; 369 Mx = 1.32;
371 fPDGencoding = -215; // a2(1320)- 370 fPDGencoding = -215; // a2(1320)-
372 } 371 }
373 else if( aParticle->GetDefinition() == G4K 372 else if( aParticle->GetDefinition() == G4KaonPlus::KaonPlus() )
374 { 373 {
375 Mx = 1.46; 374 Mx = 1.46;
376 fPDGencoding = 100321; 375 fPDGencoding = 100321;
377 } 376 }
378 else if( aParticle->GetDefinition() == G4K 377 else if( aParticle->GetDefinition() == G4KaonMinus::KaonMinus() )
379 { 378 {
380 Mx = 1.46; 379 Mx = 1.46;
381 fPDGencoding = -100321; 380 fPDGencoding = -100321;
382 } 381 }
383 } 382 }
384 else 383 else
385 { 384 {
386 if( aParticle->GetDefinition() == G4Proton 385 if( aParticle->GetDefinition() == G4Proton::Proton() )
387 { 386 {
388 Mx = 1.68; 387 Mx = 1.68;
389 // fPDGencoding = 12216; 388 // fPDGencoding = 12216;
390 fPDGencoding = 2214; 389 fPDGencoding = 2214;
391 } 390 }
392 else if( aParticle->GetDefinition() == G4N 391 else if( aParticle->GetDefinition() == G4Neutron::Neutron() )
393 { 392 {
394 Mx = 1.68; 393 Mx = 1.68;
395 fPDGencoding = 12116; 394 fPDGencoding = 12116;
396 } 395 }
397 else if( aParticle->GetDefinition() == G4P 396 else if( aParticle->GetDefinition() == G4PionPlus::PionPlus() )
398 { 397 {
399 Mx = 1.67; 398 Mx = 1.67;
400 fPDGencoding = 10215; // pi2(1670)+ 399 fPDGencoding = 10215; // pi2(1670)+
401 // Mx = 1.45; 400 // Mx = 1.45;
402 // fPDGencoding = 10211; 401 // fPDGencoding = 10211;
403 } 402 }
404 else if( aParticle->GetDefinition() == G4P 403 else if( aParticle->GetDefinition() == G4PionMinus::PionMinus() )
405 { 404 {
406 Mx = 1.67; // f0 problems->4pi vmg 2 405 Mx = 1.67; // f0 problems->4pi vmg 20.11.14
407 fPDGencoding = -10215; // pi2(1670)- 406 fPDGencoding = -10215; // pi2(1670)-
408 // Mx = 1.45; 407 // Mx = 1.45;
409 // fPDGencoding = -10211; 408 // fPDGencoding = -10211;
410 } 409 }
411 else if( aParticle->GetDefinition() == G4K 410 else if( aParticle->GetDefinition() == G4KaonPlus::KaonPlus() )
412 { 411 {
413 Mx = 1.68; 412 Mx = 1.68;
414 fPDGencoding = 30323; 413 fPDGencoding = 30323;
415 } 414 }
416 else if( aParticle->GetDefinition() == G4K 415 else if( aParticle->GetDefinition() == G4KaonMinus::KaonMinus() )
417 { 416 {
418 Mx = 1.68; 417 Mx = 1.68;
419 fPDGencoding = -30323; 418 fPDGencoding = -30323;
420 } 419 }
421 } 420 }
422 if(fPDGencoding == 0) 421 if(fPDGencoding == 0)
423 { 422 {
424 Mx = 1.44; 423 Mx = 1.44;
425 // fPDGencoding = 12212; 424 // fPDGencoding = 12212;
426 fPDGencoding = 2214; 425 fPDGencoding = 2214;
427 } 426 }
428 G4ParticleDefinition* myResonance = 427 G4ParticleDefinition* myResonance =
429 G4ParticleTable::GetParticleTable()->FindPar 428 G4ParticleTable::GetParticleTable()->FindParticle( fPDGencoding );
430 429
431 if ( myResonance ) Mx = myResonance->GetPDGM 430 if ( myResonance ) Mx = myResonance->GetPDGMass();
432 431
433 // G4cout<<"PDG-ID = "<<fPDGencoding<<"; wit 432 // G4cout<<"PDG-ID = "<<fPDGencoding<<"; with mass = "<<Mx/CLHEP::GeV<<" GeV"<<G4endl;
434 433
435 return Mx/CLHEP::GeV; 434 return Mx/CLHEP::GeV;
436 } 435 }
437 436
438 ////////////////////////////////////// 437 //////////////////////////////////////
439 // 438 //
440 // Sample t with kinematic limitations of Mx a 439 // Sample t with kinematic limitations of Mx and Tkin
441 440
442 G4double G4LMsdGenerator::SampleT( const G4Ha 441 G4double G4LMsdGenerator::SampleT( const G4HadProjectile* aParticle,
443 G4double Mx) 442 G4double Mx)
444 { 443 {
445 G4double t=0., b=0., rTkin = 50.*CLHEP::GeV, 444 G4double t=0., b=0., rTkin = 50.*CLHEP::GeV, eTkin = aParticle->GetKineticEnergy();
446 G4int i; 445 G4int i;
447 446
448 for( i = 0; i < 23; ++i) 447 for( i = 0; i < 23; ++i)
449 { 448 {
450 if( Mx <= fMxBdata[i][0] ) break; 449 if( Mx <= fMxBdata[i][0] ) break;
451 } 450 }
452 if( i <= 0 ) b = fMxBdata[0][1]; 451 if( i <= 0 ) b = fMxBdata[0][1];
453 else if( i >= 22 ) b = fMxBdata[22][1]; 452 else if( i >= 22 ) b = fMxBdata[22][1];
454 else b = fMxBdata[i][1]; 453 else b = fMxBdata[i][1];
455 454
456 if( eTkin > rTkin ) b *= 1. + G4Log(eTkin/rT 455 if( eTkin > rTkin ) b *= 1. + G4Log(eTkin/rTkin);
457 456
458 G4double rand = G4UniformRand(); 457 G4double rand = G4UniformRand();
459 458
460 t = -G4Log(rand)/b; 459 t = -G4Log(rand)/b;
461 460
462 t *= (CLHEP::GeV*CLHEP::GeV); // in G4 inter 461 t *= (CLHEP::GeV*CLHEP::GeV); // in G4 internal units
463 462
464 return t; 463 return t;
465 } 464 }
466 465
467 466
468 ////////////////////////////////////////////// 467 ////////////////////////////////////////////////
469 // 468 //
470 // Integral spectrum of Mx (GeV) 469 // Integral spectrum of Mx (GeV)
471 470
472 const G4double G4LMsdGenerator::fProbMx[60][2] 471 const G4double G4LMsdGenerator::fProbMx[60][2] =
473 { 472 {
474 {1.000000e+00, 1.000000e+00}, 473 {1.000000e+00, 1.000000e+00},
475 {1.025000e+00, 1.000000e+00}, 474 {1.025000e+00, 1.000000e+00},
476 {1.050000e+00, 1.000000e+00}, 475 {1.050000e+00, 1.000000e+00},
477 {1.075000e+00, 1.000000e+00}, 476 {1.075000e+00, 1.000000e+00},
478 {1.100000e+00, 9.975067e-01}, 477 {1.100000e+00, 9.975067e-01},
479 {1.125000e+00, 9.934020e-01}, 478 {1.125000e+00, 9.934020e-01},
480 {1.150000e+00, 9.878333e-01}, 479 {1.150000e+00, 9.878333e-01},
481 {1.175000e+00, 9.805002e-01}, 480 {1.175000e+00, 9.805002e-01},
482 {1.200000e+00, 9.716846e-01}, 481 {1.200000e+00, 9.716846e-01},
483 {1.225000e+00, 9.604761e-01}, 482 {1.225000e+00, 9.604761e-01},
484 {1.250000e+00, 9.452960e-01}, 483 {1.250000e+00, 9.452960e-01},
485 {1.275000e+00, 9.265278e-01}, 484 {1.275000e+00, 9.265278e-01},
486 {1.300000e+00, 9.053632e-01}, 485 {1.300000e+00, 9.053632e-01},
487 {1.325000e+00, 8.775566e-01}, 486 {1.325000e+00, 8.775566e-01},
488 {1.350000e+00, 8.441969e-01}, 487 {1.350000e+00, 8.441969e-01},
489 {1.375000e+00, 8.076336e-01}, 488 {1.375000e+00, 8.076336e-01},
490 {1.400000e+00, 7.682520e-01}, 489 {1.400000e+00, 7.682520e-01},
491 {1.425000e+00, 7.238306e-01}, 490 {1.425000e+00, 7.238306e-01},
492 {1.450000e+00, 6.769306e-01}, 491 {1.450000e+00, 6.769306e-01},
493 {1.475000e+00, 6.303898e-01}, 492 {1.475000e+00, 6.303898e-01},
494 {1.500000e+00, 5.824632e-01}, 493 {1.500000e+00, 5.824632e-01},
495 {1.525000e+00, 5.340696e-01}, 494 {1.525000e+00, 5.340696e-01},
496 {1.550000e+00, 4.873736e-01}, 495 {1.550000e+00, 4.873736e-01},
497 {1.575000e+00, 4.422901e-01}, 496 {1.575000e+00, 4.422901e-01},
498 {1.600000e+00, 3.988443e-01}, 497 {1.600000e+00, 3.988443e-01},
499 {1.625000e+00, 3.583727e-01}, 498 {1.625000e+00, 3.583727e-01},
500 {1.650000e+00, 3.205405e-01}, 499 {1.650000e+00, 3.205405e-01},
501 {1.675000e+00, 2.856655e-01}, 500 {1.675000e+00, 2.856655e-01},
502 {1.700000e+00, 2.537508e-01}, 501 {1.700000e+00, 2.537508e-01},
503 {1.725000e+00, 2.247863e-01}, 502 {1.725000e+00, 2.247863e-01},
504 {1.750000e+00, 1.985798e-01}, 503 {1.750000e+00, 1.985798e-01},
505 {1.775000e+00, 1.750252e-01}, 504 {1.775000e+00, 1.750252e-01},
506 {1.800000e+00, 1.539777e-01}, 505 {1.800000e+00, 1.539777e-01},
507 {1.825000e+00, 1.352741e-01}, 506 {1.825000e+00, 1.352741e-01},
508 {1.850000e+00, 1.187157e-01}, 507 {1.850000e+00, 1.187157e-01},
509 {1.875000e+00, 1.040918e-01}, 508 {1.875000e+00, 1.040918e-01},
510 {1.900000e+00, 9.118422e-02}, 509 {1.900000e+00, 9.118422e-02},
511 {1.925000e+00, 7.980909e-02}, 510 {1.925000e+00, 7.980909e-02},
512 {1.950000e+00, 6.979378e-02}, 511 {1.950000e+00, 6.979378e-02},
513 {1.975000e+00, 6.097771e-02}, 512 {1.975000e+00, 6.097771e-02},
514 {2.000000e+00, 5.322122e-02}, 513 {2.000000e+00, 5.322122e-02},
515 {2.025000e+00, 4.639628e-02}, 514 {2.025000e+00, 4.639628e-02},
516 {2.050000e+00, 4.039012e-02}, 515 {2.050000e+00, 4.039012e-02},
517 {2.075000e+00, 3.510275e-02}, 516 {2.075000e+00, 3.510275e-02},
518 {2.100000e+00, 3.044533e-02}, 517 {2.100000e+00, 3.044533e-02},
519 {2.125000e+00, 2.633929e-02}, 518 {2.125000e+00, 2.633929e-02},
520 {2.150000e+00, 2.271542e-02}, 519 {2.150000e+00, 2.271542e-02},
521 {2.175000e+00, 1.951295e-02}, 520 {2.175000e+00, 1.951295e-02},
522 {2.200000e+00, 1.667873e-02}, 521 {2.200000e+00, 1.667873e-02},
523 {2.225000e+00, 1.416633e-02}, 522 {2.225000e+00, 1.416633e-02},
524 {2.250000e+00, 1.193533e-02}, 523 {2.250000e+00, 1.193533e-02},
525 {2.275000e+00, 9.950570e-03}, 524 {2.275000e+00, 9.950570e-03},
526 {2.300000e+00, 8.181515e-03}, 525 {2.300000e+00, 8.181515e-03},
527 {2.325000e+00, 6.601664e-03}, 526 {2.325000e+00, 6.601664e-03},
528 {2.350000e+00, 5.188025e-03}, 527 {2.350000e+00, 5.188025e-03},
529 {2.375000e+00, 3.920655e-03}, 528 {2.375000e+00, 3.920655e-03},
530 {2.400000e+00, 2.782246e-03}, 529 {2.400000e+00, 2.782246e-03},
531 {2.425000e+00, 1.757765e-03}, 530 {2.425000e+00, 1.757765e-03},
532 {2.450000e+00, 8.341435e-04}, 531 {2.450000e+00, 8.341435e-04},
533 {2.475000e+00, 0.000000e+00} 532 {2.475000e+00, 0.000000e+00}
534 }; 533 };
535 534
536 ////////////////////////////////////////////// 535 //////////////////////////////////////////////
537 // 536 //
538 // Slope b (1/GeV/GeV) vs Mx (GeV) for t-sampl 537 // Slope b (1/GeV/GeV) vs Mx (GeV) for t-sampling over exp(-b*t)
539 538
540 const G4double G4LMsdGenerator::fMxBdata[23][2 539 const G4double G4LMsdGenerator::fMxBdata[23][2] =
541 { 540 {
542 {1.09014, 17.8620}, 541 {1.09014, 17.8620},
543 {1.12590, 19.2831}, 542 {1.12590, 19.2831},
544 {1.18549, 17.6907}, 543 {1.18549, 17.6907},
545 {1.21693, 16.4760}, 544 {1.21693, 16.4760},
546 {1.25194, 15.3867}, 545 {1.25194, 15.3867},
547 {1.26932, 14.4236}, 546 {1.26932, 14.4236},
548 {1.29019, 13.2931}, 547 {1.29019, 13.2931},
549 {1.30755, 12.2882}, 548 {1.30755, 12.2882},
550 {1.31790, 11.4509}, 549 {1.31790, 11.4509},
551 {1.33888, 10.6969}, 550 {1.33888, 10.6969},
552 {1.34911, 9.44130}, 551 {1.34911, 9.44130},
553 {1.37711, 8.56148}, 552 {1.37711, 8.56148},
554 {1.39101, 7.76593}, 553 {1.39101, 7.76593},
555 {1.42608, 6.88582}, 554 {1.42608, 6.88582},
556 {1.48593, 6.13019}, 555 {1.48593, 6.13019},
557 {1.53179, 5.87723}, 556 {1.53179, 5.87723},
558 {1.58111, 5.37308}, 557 {1.58111, 5.37308},
559 {1.64105, 4.95217}, 558 {1.64105, 4.95217},
560 {1.69037, 4.44803}, 559 {1.69037, 4.44803},
561 {1.81742, 3.89879}, 560 {1.81742, 3.89879},
562 {1.88096, 3.68693}, 561 {1.88096, 3.68693},
563 {1.95509, 3.43278}, 562 {1.95509, 3.43278},
564 {2.02219, 3.30445} 563 {2.02219, 3.30445}
565 }; 564 };
566 565
567 566
568 567
569 // 568 //
570 // 569 //
571 ///////////////////////////////////////////// 570 /////////////////////////////////////////////
572 571