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

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

update CVS release candidate geant4.9.3.01
File size: 11.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.11 2009/10/06 10:10:36 vuzhinsk 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//G4cout<<"G4ElasticHNScattering::ElasticScattering"<<G4endl;
61// -------------------- Projectile parameters -----------------------------------
62 G4LorentzVector Pprojectile=projectile->Get4Momentum();
63//G4cout<<"Elastic P "<<Pprojectile<<G4endl;
64 if(Pprojectile.z() < 0.)
65 {
66 target->SetStatus(2);
67 return false;
68 }
69
70// G4double ProjectileRapidity = Pprojectile.rapidity();
71// G4int ProjectilePDGcode=projectile->GetDefinition()->GetPDGEncoding();
72// G4ParticleDefinition * ProjectileDefinition = projectile->GetDefinition();
73
74 G4bool PutOnMassShell(false);
75
76 G4double M0projectile = Pprojectile.mag();
77 if(M0projectile < projectile->GetDefinition()->GetPDGMass())
78 {
79 PutOnMassShell=true;
80 M0projectile=projectile->GetDefinition()->GetPDGMass();
81 }
82
83 G4double Mprojectile2 = M0projectile * M0projectile;
84
85 G4double AveragePt2=theParameters->GetAvaragePt2ofElasticScattering();
86
87// -------------------- Target parameters ----------------------------------------------
88// G4int TargetPDGcode=target->GetDefinition()->GetPDGEncoding();
89
90 G4LorentzVector Ptarget=target->Get4Momentum();
91//G4cout<<"Elastic T "<<Ptarget<<G4endl;
92// G4double TargetRapidity = Ptarget.rapidity();
93 G4double M0target = Ptarget.mag();
94
95 if(M0target < target->GetDefinition()->GetPDGMass())
96 {
97 PutOnMassShell=true;
98 M0target=target->GetDefinition()->GetPDGMass();
99 }
100
101 G4double Mtarget2 = M0target * M0target;
102
103// Transform momenta to cms and then rotate parallel to z axis;
104
105 G4LorentzVector Psum;
106 Psum=Pprojectile+Ptarget;
107
108 G4LorentzRotation toCms(-1*Psum.boostVector());
109
110 G4LorentzVector Ptmp=toCms*Pprojectile;
111
112 if ( Ptmp.pz() <= 0. )
113 {
114 // "String" moving backwards in CMS, abort collision !!
115 //G4cout << " abort Collision!! " << G4endl;
116 target->SetStatus(2);
117 return false;
118 }
119
120 toCms.rotateZ(-1*Ptmp.phi());
121 toCms.rotateY(-1*Ptmp.theta());
122
123 G4LorentzRotation toLab(toCms.inverse());
124
125 Pprojectile.transform(toCms);
126 Ptarget.transform(toCms);
127
128// ---------------------- Putting on mass-on-shell, if needed ------------------------
129 G4double PZcms2, PZcms;
130
131 G4double S=Psum.mag2();
132// G4double SqrtS=std::sqrt(S);
133
134 PZcms2=(S*S+Mprojectile2*Mprojectile2+Mtarget2*Mtarget2-
135 2*S*Mprojectile2-2*S*Mtarget2-2*Mprojectile2*Mtarget2)/4./S;
136
137 if(PZcms2 < 0.)
138 { // It can be in an interaction with off-shell nuclear nucleon
139 if(M0projectile > projectile->GetDefinition()->GetPDGMass())
140 { // An attempt to de-excite the projectile
141 // It is assumed that the target is in the ground state
142 M0projectile = projectile->GetDefinition()->GetPDGMass();
143 Mprojectile2=M0projectile*M0projectile;
144 PZcms2=(S*S+Mprojectile2*Mprojectile2+Mtarget2*Mtarget2-
145 2*S*Mprojectile2 - 2*S*Mtarget2 - 2*Mprojectile2*Mtarget2)
146 /4./S;
147 if(PZcms2 < 0.){ return false;} // Non succesful attempt after the de-excitation
148 }
149 else // if(M0projectile > projectile->GetDefinition()->GetPDGMass())
150 {
151 target->SetStatus(2); // Uzhi 17.07.09
152 return false; // The projectile was not excited,
153 // but the energy was too low to put
154 // the target nucleon on mass-shell
155 } // end of if(M0projectile > projectile->GetDefinition()->GetPDGMass())
156 } // end of if(PZcms2 < 0.)
157
158 PZcms = std::sqrt(PZcms2);
159
160 if(PutOnMassShell)
161 {
162 if(Pprojectile.z() > 0.)
163 {
164 Pprojectile.setPz( PZcms);
165 Ptarget.setPz( -PZcms);
166 }
167 else // if(Pprojectile.z() > 0.)
168 {
169 Pprojectile.setPz(-PZcms);
170 Ptarget.setPz( PZcms);
171 };
172
173 Pprojectile.setE(std::sqrt(Mprojectile2+
174 Pprojectile.x()*Pprojectile.x()+
175 Pprojectile.y()*Pprojectile.y()+
176 PZcms2));
177 Ptarget.setE(std::sqrt( Mtarget2 +
178 Ptarget.x()*Ptarget.x()+
179 Ptarget.y()*Ptarget.y()+
180 PZcms2));
181 } // end of if(PutOnMassShell)
182
183 G4double maxPtSquare = PZcms2;
184
185//----- Charge exchange between the projectile and the target, if possible
186// (G4UniformRand() < 0.5*std::exp(-0.5*(ProjectileRapidity - TargetRapidity))))
187/*
188 if((ProjectilePDGcode != TargetPDGcode) &&
189 ((ProjectilePDGcode > 1000) && (TargetPDGcode > 1000)) &&
190 (G4UniformRand() < 1.0*std::exp(-0.5*(ProjectileRapidity - TargetRapidity))))
191 {
192 projectile->SetDefinition(target->GetDefinition());
193 target->SetDefinition(ProjectileDefinition);
194 }
195*/
196// ------ Now we can calculate the transfered Pt --------------------------
197 G4double Pt2;
198 G4double ProjMassT2, ProjMassT;
199 G4double TargMassT2, TargMassT;
200
201 G4LorentzVector Qmomentum;
202
203 Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0);
204
205 Pt2=G4ThreeVector(Qmomentum.vect()).mag2();
206
207 ProjMassT2=Mprojectile2+Pt2;
208 ProjMassT =std::sqrt(ProjMassT2);
209
210 TargMassT2=Mtarget2+Pt2;
211 TargMassT =std::sqrt(TargMassT2);
212
213 PZcms2=(S*S+ProjMassT2*ProjMassT2+
214 TargMassT2*TargMassT2-
215 2.*S*ProjMassT2-2.*S*TargMassT2-
216 2.*ProjMassT2*TargMassT2)/4./S;
217 if(PZcms2 < 0 ) {PZcms2=0;};// to avoid the exactness problem
218 PZcms =std::sqrt(PZcms2);
219
220 Pprojectile.setPz( PZcms);
221 Ptarget.setPz( -PZcms);
222
223 Pprojectile += Qmomentum;
224 Ptarget -= Qmomentum;
225
226// Transform back and update SplitableHadron Participant.
227 Pprojectile.transform(toLab);
228 Ptarget.transform(toLab);
229/* // Maybe it will be needed for an exact calculations--------------------
230 G4double TargetMomentum=std::sqrt(Ptarget.x()*Ptarget.x()+
231 Ptarget.y()*Ptarget.y()+
232 Ptarget.z()*Ptarget.z());
233*/
234
235// Calculation of the creation time ---------------------
236 projectile->SetTimeOfCreation(target->GetTimeOfCreation());
237 projectile->SetPosition(target->GetPosition());
238// Creation time and position of target nucleon were determined at
239// ReggeonCascade() of G4FTFModel
240// ------------------------------------------------------
241
242 projectile->Set4Momentum(Pprojectile);
243 target->Set4Momentum(Ptarget);
244
245 projectile->IncrementCollisionCount(1);
246 target->IncrementCollisionCount(1);
247
248 return true;
249}
250
251
252// --------- private methods ----------------------
253
254G4ThreeVector G4ElasticHNScattering::GaussianPt(G4double AveragePt2, G4double maxPtSquare) const
255{ // @@ this method is used in FTFModel as well. Should go somewhere common!
256
257 G4double Pt2;
258 Pt2 = -AveragePt2 * std::log(1. + G4UniformRand() *
259 (std::exp(-maxPtSquare/AveragePt2)-1.));
260
261 G4double Pt=std::sqrt(Pt2);
262 G4double phi=G4UniformRand() * twopi;
263
264 return G4ThreeVector (Pt*std::cos(phi), Pt*std::sin(phi), 0.);
265}
266
267G4ElasticHNScattering::G4ElasticHNScattering(const G4ElasticHNScattering &)
268{
269 throw G4HadronicException(__FILE__, __LINE__, "G4ElasticHNScattering copy contructor not meant to be called");
270}
271
272
273G4ElasticHNScattering::~G4ElasticHNScattering()
274{
275}
276
277
278const G4ElasticHNScattering & G4ElasticHNScattering::operator=(const G4ElasticHNScattering &)
279{
280 throw G4HadronicException(__FILE__, __LINE__, "G4ElasticHNScattering = operator meant to be called");
281 return *this;
282}
283
284
285int G4ElasticHNScattering::operator==(const G4ElasticHNScattering &) const
286{
287 throw G4HadronicException(__FILE__, __LINE__, "G4ElasticHNScattering == operator meant to be called");
288 return false;
289}
290
291int G4ElasticHNScattering::operator!=(const G4ElasticHNScattering &) const
292{
293 throw G4HadronicException(__FILE__, __LINE__, "G4ElasticHNScattering != operator meant to be called");
294 return true;
295}
Note: See TracBrowser for help on using the repository browser.