source: trunk/source/processes/hadronic/models/parton_string/diffraction/src/G4ElasticHNScattering.cc@ 1315

Last change on this file since 1315 was 1228, checked in by garnier, 16 years ago

update geant4.9.3 tag
File size: 10.8 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: G4ElasticHNScattering.cc,v 1.14 2009/12/16 17:51:13 gunter Exp $
28// ------------------------------------------------------------
29// GEANT 4 class implemetation file
30//
31// ---------------- G4ElasticHNScattering --------------
32// by V. Uzhinsky, March 2008.
33// elastic scattering used by Fritiof model
34// Take a projectile and a target
35// scatter the projectile and target
36// ---------------------------------------------------------------------
37
38
39#include "globals.hh"
40#include "Randomize.hh"
41
42#include "G4ElasticHNScattering.hh"
43#include "G4LorentzRotation.hh"
44#include "G4ThreeVector.hh"
45#include "G4ParticleDefinition.hh"
46#include "G4VSplitableHadron.hh"
47#include "G4ExcitedString.hh"
48#include "G4FTFParameters.hh"
49//#include "G4ios.hh"
50
51G4ElasticHNScattering::G4ElasticHNScattering()
52{
53}
54
55G4bool G4ElasticHNScattering::
56 ElasticScattering (G4VSplitableHadron *projectile,
57 G4VSplitableHadron *target,
58 G4FTFParameters *theParameters) const
59{
60// -------------------- Projectile parameters -----------------------------------
61 G4LorentzVector Pprojectile=projectile->Get4Momentum();
62
63 if(Pprojectile.z() < 0.)
64 {
65 target->SetStatus(2);
66 return false;
67 }
68
69 G4bool PutOnMassShell(false);
70
71 G4double M0projectile = Pprojectile.mag();
72 if(M0projectile < projectile->GetDefinition()->GetPDGMass())
73 {
74 PutOnMassShell=true;
75 M0projectile=projectile->GetDefinition()->GetPDGMass();
76 }
77
78 G4double Mprojectile2 = M0projectile * M0projectile;
79
80 G4double AveragePt2=theParameters->GetAvaragePt2ofElasticScattering();
81
82// -------------------- Target parameters ----------------------------------------------
83
84 G4LorentzVector Ptarget=target->Get4Momentum();
85
86 G4double M0target = Ptarget.mag();
87
88 if(M0target < target->GetDefinition()->GetPDGMass())
89 {
90 PutOnMassShell=true;
91 M0target=target->GetDefinition()->GetPDGMass();
92 }
93
94 G4double Mtarget2 = M0target * M0target;
95
96// Transform momenta to cms and then rotate parallel to z axis;
97
98 G4LorentzVector Psum;
99 Psum=Pprojectile+Ptarget;
100
101 G4LorentzRotation toCms(-1*Psum.boostVector());
102
103 G4LorentzVector Ptmp=toCms*Pprojectile;
104
105 if ( Ptmp.pz() <= 0. )
106 {
107 // "String" moving backwards in CMS, abort collision !!
108 //G4cout << " abort Collision!! " << G4endl;
109 target->SetStatus(2);
110 return false;
111 }
112
113 toCms.rotateZ(-1*Ptmp.phi());
114 toCms.rotateY(-1*Ptmp.theta());
115
116 G4LorentzRotation toLab(toCms.inverse());
117
118 Pprojectile.transform(toCms);
119 Ptarget.transform(toCms);
120
121// ---------------------- Putting on mass-on-shell, if needed ------------------------
122 G4double PZcms2, PZcms;
123
124 G4double S=Psum.mag2();
125// G4double SqrtS=std::sqrt(S);
126
127 PZcms2=(S*S+Mprojectile2*Mprojectile2+Mtarget2*Mtarget2-
128 2*S*Mprojectile2-2*S*Mtarget2-2*Mprojectile2*Mtarget2)/4./S;
129
130 if(PZcms2 < 0.)
131 { // It can be in an interaction with off-shell nuclear nucleon
132 if(M0projectile > projectile->GetDefinition()->GetPDGMass())
133 { // An attempt to de-excite the projectile
134 // It is assumed that the target is in the ground state
135 M0projectile = projectile->GetDefinition()->GetPDGMass();
136 Mprojectile2=M0projectile*M0projectile;
137 PZcms2=(S*S+Mprojectile2*Mprojectile2+Mtarget2*Mtarget2-
138 2*S*Mprojectile2 - 2*S*Mtarget2 - 2*Mprojectile2*Mtarget2)
139 /4./S;
140 if(PZcms2 < 0.){ return false;} // Non succesful attempt after the de-excitation
141 }
142 else // if(M0projectile > projectile->GetDefinition()->GetPDGMass())
143 {
144 target->SetStatus(2);
145 return false; // The projectile was not excited,
146 // but the energy was too low to put
147 // the target nucleon on mass-shell
148 } // end of if(M0projectile > projectile->GetDefinition()->GetPDGMass())
149 } // end of if(PZcms2 < 0.)
150
151 PZcms = std::sqrt(PZcms2);
152
153 if(PutOnMassShell)
154 {
155 if(Pprojectile.z() > 0.)
156 {
157 Pprojectile.setPz( PZcms);
158 Ptarget.setPz( -PZcms);
159 }
160 else // if(Pprojectile.z() > 0.)
161 {
162 Pprojectile.setPz(-PZcms);
163 Ptarget.setPz( PZcms);
164 };
165
166 Pprojectile.setE(std::sqrt(Mprojectile2+
167 Pprojectile.x()*Pprojectile.x()+
168 Pprojectile.y()*Pprojectile.y()+
169 PZcms2));
170 Ptarget.setE(std::sqrt( Mtarget2 +
171 Ptarget.x()*Ptarget.x()+
172 Ptarget.y()*Ptarget.y()+
173 PZcms2));
174 } // end of if(PutOnMassShell)
175
176 G4double maxPtSquare = PZcms2;
177
178// ------ Now we can calculate the transfered Pt --------------------------
179 G4double Pt2;
180 G4double ProjMassT2, ProjMassT;
181 G4double TargMassT2, TargMassT;
182
183 G4LorentzVector Qmomentum;
184
185 Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0);
186
187 Pt2=G4ThreeVector(Qmomentum.vect()).mag2();
188
189 ProjMassT2=Mprojectile2+Pt2;
190 ProjMassT =std::sqrt(ProjMassT2);
191
192 TargMassT2=Mtarget2+Pt2;
193 TargMassT =std::sqrt(TargMassT2);
194
195 PZcms2=(S*S+ProjMassT2*ProjMassT2+
196 TargMassT2*TargMassT2-
197 2.*S*ProjMassT2-2.*S*TargMassT2-
198 2.*ProjMassT2*TargMassT2)/4./S;
199 if(PZcms2 < 0 ) {PZcms2=0;};// to avoid the exactness problem
200 PZcms =std::sqrt(PZcms2);
201
202 Pprojectile.setPz( PZcms);
203 Ptarget.setPz( -PZcms);
204
205 Pprojectile += Qmomentum;
206 Ptarget -= Qmomentum;
207
208// Transform back and update SplitableHadron Participant.
209 Pprojectile.transform(toLab);
210 Ptarget.transform(toLab);
211/* // Maybe it will be needed for an exact calculations--------------------
212 G4double TargetMomentum=std::sqrt(Ptarget.x()*Ptarget.x()+
213 Ptarget.y()*Ptarget.y()+
214 Ptarget.z()*Ptarget.z());
215*/
216
217// Calculation of the creation time ---------------------
218 projectile->SetTimeOfCreation(target->GetTimeOfCreation());
219 projectile->SetPosition(target->GetPosition());
220// Creation time and position of target nucleon were determined at
221// ReggeonCascade() of G4FTFModel
222// ------------------------------------------------------
223
224 projectile->Set4Momentum(Pprojectile);
225 target->Set4Momentum(Ptarget);
226
227 projectile->IncrementCollisionCount(1);
228 target->IncrementCollisionCount(1);
229
230 return true;
231}
232
233
234// --------- private methods ----------------------
235
236G4ThreeVector G4ElasticHNScattering::GaussianPt(G4double AveragePt2, G4double maxPtSquare) const
237{ // @@ this method is used in FTFModel as well. Should go somewhere common!
238
239 G4double Pt2(0.);
240 if(AveragePt2 <= 0.) {Pt2=0.;}
241 else
242 {
243 Pt2 = -AveragePt2 * std::log(1. + G4UniformRand() *
244 (std::exp(-maxPtSquare/AveragePt2)-1.));
245 }
246 G4double Pt=std::sqrt(Pt2);
247 G4double phi=G4UniformRand() * twopi;
248
249 return G4ThreeVector (Pt*std::cos(phi), Pt*std::sin(phi), 0.);
250}
251
252G4ElasticHNScattering::G4ElasticHNScattering(const G4ElasticHNScattering &)
253{
254 throw G4HadronicException(__FILE__, __LINE__, "G4ElasticHNScattering copy contructor not meant to be called");
255}
256
257
258G4ElasticHNScattering::~G4ElasticHNScattering()
259{
260}
261
262
263const G4ElasticHNScattering & G4ElasticHNScattering::operator=(const G4ElasticHNScattering &)
264{
265 throw G4HadronicException(__FILE__, __LINE__, "G4ElasticHNScattering = operator meant to be called");
266 return *this;
267}
268
269
270int G4ElasticHNScattering::operator==(const G4ElasticHNScattering &) const
271{
272 throw G4HadronicException(__FILE__, __LINE__, "G4ElasticHNScattering == operator meant to be called");
273 return false;
274}
275
276int G4ElasticHNScattering::operator!=(const G4ElasticHNScattering &) const
277{
278 throw G4HadronicException(__FILE__, __LINE__, "G4ElasticHNScattering != operator meant to be called");
279 return true;
280}
Note: See TracBrowser for help on using the repository browser.