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

Last change on this file since 1036 was 968, checked in by garnier, 17 years ago

fichier ajoutes
File size: 11.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: G4ElasticHNScattering.cc,v 1.3 2008/05/19 12:56: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" // Uzhi 29.03.08
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
62 G4LorentzVector Pprojectile=projectile->Get4Momentum();
63
64// -------------------- Projectile parameters -----------------------------------
65 G4bool PutOnMassShell=0;
66
67 G4double M0projectile = Pprojectile.mag();
68
69 if(M0projectile < projectile->GetDefinition()->GetPDGMass())
70 {
71 PutOnMassShell=1;
72 M0projectile=projectile->GetDefinition()->GetPDGMass();
73 }
74
75 G4double Mprojectile2 = M0projectile * M0projectile;
76
77// G4double AveragePt2=theParameters->GetSlope(); // Uzhi ???
78// AveragePt2 = AveragePt2 * GeV*GeV;
79
80 G4double AveragePt2=theParameters->GetAvaragePt2ofElasticScattering();
81
82// -------------------- Target parameters ----------------------------------------------
83 G4LorentzVector Ptarget=target->Get4Momentum();
84
85 G4double M0target = Ptarget.mag();
86//G4cout<<" Mp Mt Pt2 "<<M0projectile<<" "<<M0target<<" "<<AveragePt2/GeV/GeV<<G4endl;
87
88 if(M0target < target->GetDefinition()->GetPDGMass())
89 {
90 PutOnMassShell=1;
91 M0target=target->GetDefinition()->GetPDGMass();
92 }
93
94 G4double Mtarget2 = M0target * M0target; //Ptarget.mag2();
95 // for AA-inter.
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. ) // Uzhi ???
106 {
107 // "String" moving backwards in CMS, abort collision !!
108 //G4cout << " abort Collision!! " << G4endl;
109 return false;
110 }
111
112 toCms.rotateZ(-1*Ptmp.phi());
113 toCms.rotateY(-1*Ptmp.theta());
114
115 G4LorentzRotation toLab(toCms.inverse());
116
117 Pprojectile.transform(toCms);
118 Ptarget.transform(toCms);
119
120// ---------------------- Sampling of transfered Pt ------------------------
121 G4double Pt2;
122 G4double ProjMassT2, ProjMassT;
123 G4double TargMassT2, TargMassT;
124 G4double PZcms2, PZcms;
125
126 G4double S=Psum.mag2();
127// G4double SqrtS=std::sqrt(S);
128
129 PZcms2=(S*S+Mprojectile2*Mprojectile2+Mtarget2*Mtarget2-
130 2*S*Mprojectile2-2*S*Mtarget2-2*Mprojectile2*Mtarget2)/4./S;
131 if(PZcms2 < 0)
132 {return false;} // It can be in an interaction with off-shell nuclear nucleon
133
134 PZcms = std::sqrt(PZcms2);
135
136 if(PutOnMassShell)
137 {
138 if(Pprojectile.z() > 0.)
139 {
140 Pprojectile.setPz( PZcms);
141 Ptarget.setPz( -PZcms);
142 }
143 else
144 {
145 Pprojectile.setPz(-PZcms);
146 Ptarget.setPz( PZcms);
147 };
148
149 Pprojectile.setE(std::sqrt(Mprojectile2+
150 Pprojectile.x()*Pprojectile.x()+
151 Pprojectile.y()*Pprojectile.y()+
152 PZcms2));
153 Ptarget.setE(std::sqrt( Mtarget2 +
154 Ptarget.x()*Ptarget.x()+
155 Ptarget.y()*Ptarget.y()+
156 PZcms2));
157 }
158
159 G4double maxPtSquare = PZcms2;
160 G4LorentzVector Qmomentum;
161
162 Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0);
163
164 Pt2=G4ThreeVector(Qmomentum.vect()).mag2();
165//G4cout<<"Pt2 GeV^2 "<<(Pt2)/GeV/GeV<<G4endl;
166
167 ProjMassT2=Mprojectile2+Pt2;
168 ProjMassT =std::sqrt(ProjMassT2);
169
170 TargMassT2=Mtarget2+Pt2;
171 TargMassT =std::sqrt(TargMassT2);
172
173 PZcms2=(S*S+ProjMassT2*ProjMassT2+
174 TargMassT2*TargMassT2-
175 2.*S*ProjMassT2-2.*S*TargMassT2-
176 2.*ProjMassT2*TargMassT2)/4./S;
177 if(PZcms2 < 0 ) {PZcms2=0;};
178 PZcms =std::sqrt(PZcms2);
179
180 Pprojectile.setPz( PZcms); // Uzhi Proj can move backward
181 Ptarget.setPz( -PZcms); // Uzhi Proj can move backward
182
183//G4cout << "Qplus / Qminus " << Qplus << " / " << Qminus<<G4endl;
184// G4cout << "pt2" << pt2 << G4endl;
185// G4cout << "Qmomentum " << Qmomentum << G4endl;
186// G4cout << " Masses (P/T) : " << (Pprojectile+Qmomentum).mag() <<
187// " / " << (Ptarget-Qmomentum).mag() << G4endl;
188
189 Pprojectile += Qmomentum;
190 Ptarget -= Qmomentum;
191
192//G4cout << "Pprojectile with Q : " << Pprojectile << G4endl;
193//G4cout << "Ptarget with Q : " << Ptarget << G4endl;
194
195// G4cout << "Projectile back: " << toLab * Pprojectile << G4endl;
196// G4cout << "Target back: " << toLab * Ptarget << G4endl;
197
198// Transform back and update SplitableHadron Participant.
199 Pprojectile.transform(toLab);
200 Ptarget.transform(toLab);
201
202//G4cout << "Pprojectile with Q M: " << Pprojectile<<" "<< Pprojectile.mag() << G4endl;
203//G4cout << "Ptarget with Q M: " << Ptarget <<" "<< Ptarget.mag() << G4endl;
204
205//G4cout << "Target mass " << Ptarget.mag() << G4endl;
206
207//G4cout << "Projectile mass " << Pprojectile.mag() << G4endl;
208
209 G4double ZcoordinateOfCurrentInteraction = target->GetPosition().z();
210// It is assumed that nucleon z-coordinates are ordered on increasing -----------
211
212 G4double betta_z=projectile->Get4Momentum().pz()/projectile->Get4Momentum().e();
213
214 G4double ZcoordinateOfPreviousCollision=projectile->GetPosition().z();
215 if(projectile->GetSoftCollisionCount()==0) {
216 projectile->SetTimeOfCreation(0.);
217 target->SetTimeOfCreation(0.);
218 ZcoordinateOfPreviousCollision=ZcoordinateOfCurrentInteraction;
219 }
220
221 G4ThreeVector thePosition(projectile->GetPosition().x(),
222 projectile->GetPosition().y(),
223 ZcoordinateOfCurrentInteraction);
224 projectile->SetPosition(thePosition);
225
226 G4double TimeOfPreviousCollision=projectile->GetTimeOfCreation();
227 G4double TimeOfCurrentCollision=TimeOfPreviousCollision+
228 (ZcoordinateOfCurrentInteraction-ZcoordinateOfPreviousCollision)/betta_z;
229
230 projectile->SetTimeOfCreation(TimeOfCurrentCollision);
231 target->SetTimeOfCreation(TimeOfCurrentCollision);
232
233 projectile->Set4Momentum(Pprojectile);
234 target->Set4Momentum(Ptarget);
235
236 projectile->IncrementCollisionCount(1);
237 target->IncrementCollisionCount(1);
238
239 return true;
240}
241
242
243// --------- private methods ----------------------
244
245G4ThreeVector G4ElasticHNScattering::GaussianPt(G4double AveragePt2, G4double maxPtSquare) const
246{ // @@ this method is used in FTFModel as well. Should go somewhere common!
247
248 G4double Pt2;
249 Pt2 = -AveragePt2 * std::log(1. + G4UniformRand() *
250 (std::exp(-maxPtSquare/AveragePt2)-1.));
251
252 G4double Pt=std::sqrt(Pt2);
253 G4double phi=G4UniformRand() * twopi;
254
255 return G4ThreeVector (Pt*std::cos(phi), Pt*std::sin(phi), 0.);
256}
257
258G4ElasticHNScattering::G4ElasticHNScattering(const G4ElasticHNScattering &)
259{
260 throw G4HadronicException(__FILE__, __LINE__, "G4ElasticHNScattering copy contructor not meant to be called");
261}
262
263
264G4ElasticHNScattering::~G4ElasticHNScattering()
265{
266}
267
268
269const G4ElasticHNScattering & G4ElasticHNScattering::operator=(const G4ElasticHNScattering &)
270{
271 throw G4HadronicException(__FILE__, __LINE__, "G4ElasticHNScattering = operator meant to be called");
272 return *this;
273}
274
275
276int G4ElasticHNScattering::operator==(const G4ElasticHNScattering &) const
277{
278 throw G4HadronicException(__FILE__, __LINE__, "G4ElasticHNScattering == operator meant to be called");
279 return false;
280}
281
282int G4ElasticHNScattering::operator!=(const G4ElasticHNScattering &) const
283{
284 throw G4HadronicException(__FILE__, __LINE__, "G4ElasticHNScattering != operator meant to be called");
285 return true;
286}
Note: See TracBrowser for help on using the repository browser.