Geant4-11
Public Member Functions | Private Member Functions
G4QGSDiffractiveExcitation Class Reference

#include <G4QGSDiffractiveExcitation.hh>

Inheritance diagram for G4QGSDiffractiveExcitation:
G4SingleDiffractiveExcitation

Public Member Functions

virtual G4bool ExciteParticipants (G4VSplitableHadron *aPartner, G4VSplitableHadron *bPartner, G4bool ProjectileDiffraction=TRUE) const
 
 G4QGSDiffractiveExcitation ()
 
virtual G4ExcitedStringString (G4VSplitableHadron *aHadron, G4bool isProjectile) const
 
virtual ~G4QGSDiffractiveExcitation ()
 

Private Member Functions

G4double ChooseP (G4double Pmin, G4double Pmax) const
 
 G4QGSDiffractiveExcitation (const G4QGSDiffractiveExcitation &right)
 
G4ThreeVector GaussianPt (G4double AveragePt2, G4double maxPtSquare) const
 
G4bool operator!= (const G4QGSDiffractiveExcitation &right) const
 
const G4QGSDiffractiveExcitationoperator= (const G4QGSDiffractiveExcitation &right)
 
G4bool operator== (const G4QGSDiffractiveExcitation &right) const
 

Detailed Description

Definition at line 50 of file G4QGSDiffractiveExcitation.hh.

Constructor & Destructor Documentation

◆ G4QGSDiffractiveExcitation() [1/2]

G4QGSDiffractiveExcitation::G4QGSDiffractiveExcitation ( )

Definition at line 68 of file G4QGSDiffractiveExcitation.cc.

69{
70}

◆ ~G4QGSDiffractiveExcitation()

G4QGSDiffractiveExcitation::~G4QGSDiffractiveExcitation ( )
virtual

Definition at line 72 of file G4QGSDiffractiveExcitation.cc.

73{
74}

◆ G4QGSDiffractiveExcitation() [2/2]

G4QGSDiffractiveExcitation::G4QGSDiffractiveExcitation ( const G4QGSDiffractiveExcitation right)
private

Member Function Documentation

◆ ChooseP()

G4double G4QGSDiffractiveExcitation::ChooseP ( G4double  Pmin,
G4double  Pmax 
) const
private

Definition at line 375 of file G4QGSDiffractiveExcitation.cc.

376{
377 // choose an x between Xmin and Xmax with P(x) ~ 1/x
378 // to be improved...
379
380 G4double range=Pmax-Pmin;
381
382 if ( Pmin <= 0. || range <=0. )
383 {
384 G4cout << " Pmin, range : " << Pmin << " , " << range << G4endl;
385 throw G4HadronicException(__FILE__, __LINE__, "G4QGSDiffractiveExcitation::ChooseP : Invalid arguments ");
386 }
387
388 G4double P;
390 //debug-hpw cout << "DiffractiveX "<<x<<G4endl;
391 return P;
392}
double G4double
Definition: G4Types.hh:83
#define G4endl
Definition: G4ios.hh:57
G4GLOB_DLL std::ostream G4cout
#define G4UniformRand()
Definition: Randomize.hh:52
static G4Pow * GetInstance()
Definition: G4Pow.cc:41
G4double powA(G4double A, G4double y) const
Definition: G4Pow.hh:230
static double P[]

References G4cout, G4endl, G4UniformRand, G4Pow::GetInstance(), P, anonymous_namespace{G4ChipsKaonMinusInelasticXS.cc}::Pmax, anonymous_namespace{G4ChipsKaonMinusInelasticXS.cc}::Pmin, and G4Pow::powA().

Referenced by ExciteParticipants().

◆ ExciteParticipants()

G4bool G4QGSDiffractiveExcitation::ExciteParticipants ( G4VSplitableHadron aPartner,
G4VSplitableHadron bPartner,
G4bool  ProjectileDiffraction = TRUE 
) const
virtual

Reimplemented in G4SingleDiffractiveExcitation.

Definition at line 77 of file G4QGSDiffractiveExcitation.cc.

79{
80 #ifdef debugDoubleDiffraction
81 G4cout<<G4endl<<"G4QGSDiffractiveExcitation::ExciteParticipants - Double diffraction."<<G4endl;
82 G4cout<<"Proj Targ "<<projectile->GetDefinition()->GetParticleName()<<" "<<target->GetDefinition()->GetParticleName()<<G4endl;
83 G4cout<<"Proj 4 Mom "<<projectile->Get4Momentum()<<" "<<projectile->Get4Momentum().mag()<<G4endl;
84 G4cout<<"Targ 4 Mom "<<target->Get4Momentum() <<" "<<target->Get4Momentum().mag() <<G4endl;
85 #endif
86
87 G4LorentzVector Pprojectile=projectile->Get4Momentum();
88
89 // -------------------- Projectile parameters -----------------------------------
90 G4bool PutOnMassShell=0;
91
92 G4double M0projectile = Pprojectile.mag(); // Without de-excitation
93
94 if(M0projectile < projectile->GetDefinition()->GetPDGMass())
95 {
96 PutOnMassShell=1;
97 M0projectile=projectile->GetDefinition()->GetPDGMass();
98 }
99
100 // -------------------- Target parameters ----------------------------------------------
101 G4LorentzVector Ptarget=target->Get4Momentum();
102
103 G4double M0target = Ptarget.mag();
104
105 if(M0target < target->GetDefinition()->GetPDGMass())
106 {
107 PutOnMassShell=1;
108 M0target=target->GetDefinition()->GetPDGMass();
109 }
110
111 G4LorentzVector Psum=Pprojectile+Ptarget;
112 G4double S=Psum.mag2();
113 G4double SqrtS=std::sqrt(S);
114
115 if (SqrtS < M0projectile + M0target) {return false;} // The model cannot work for pp-interactions
116 // at Plab < 1.3 GeV/c.
117
118 G4double Mprojectile2 = M0projectile * M0projectile;
119 G4double Mtarget2 = M0target * M0target; //Ptarget.mag2(); // for AA-inter.
120
121 // Transform momenta to cms and then rotate parallel to z axis;
122
123 G4LorentzRotation toCms(-1*Psum.boostVector());
124
125 G4LorentzVector Ptmp=toCms*Pprojectile;
126
127 if ( Ptmp.pz() <= 0. )
128 {
129 // "String" moving backwards in CMS, abort collision !!
130 //G4cout << " abort Collision!! " << G4endl;
131 return false;
132 }
133
134 toCms.rotateZ(-1*Ptmp.phi());
135 toCms.rotateY(-1*Ptmp.theta());
136
137 G4LorentzRotation toLab(toCms.inverse());
138
139 Pprojectile.transform(toCms);
140 Ptarget.transform(toCms);
141
142 G4double PZcms2=(S*S+Mprojectile2*Mprojectile2+Mtarget2*Mtarget2-
143 2*S*Mprojectile2-2*S*Mtarget2-2*Mprojectile2*Mtarget2)/4./S;
144
145 if (PZcms2 < 0) {return false;} // It can be in an interaction with off-shell nuclear nucleon
146
147 G4double PZcms = std::sqrt(PZcms2);
148
149 if (PutOnMassShell)
150 {
151 if (Pprojectile.z() > 0.)
152 {
153 Pprojectile.setPz( PZcms);
154 Ptarget.setPz( -PZcms);
155 }
156 else
157 {
158 Pprojectile.setPz(-PZcms);
159 Ptarget.setPz( PZcms);
160 };
161
162 Pprojectile.setE(std::sqrt(Mprojectile2+sqr(Pprojectile.x())+sqr(Pprojectile.y())+PZcms2));
163 Ptarget.setE( std::sqrt( Mtarget2+sqr( Ptarget.x())+sqr( Ptarget.y())+PZcms2));
164 }
165
166 G4double maxPtSquare = PZcms2;
167
168 #ifdef debugDoubleDiffraction
169 G4cout << "Pprojectile after boost to CMS: " << Pprojectile <<" "<<Pprojectile.mag()<<G4endl;
170 G4cout << "Ptarget after boost to CMS: " << Ptarget <<" "<<Ptarget.mag() <<G4endl;
171 #endif
172
173 G4int PrPDGcode=projectile->GetDefinition()->GetPDGEncoding();
174 G4int absPrPDGcode=std::abs(PrPDGcode);
175 G4double MinPrDiffMass(0.);
176 G4double AveragePt2(0.);
177
178 if (M0projectile <= projectile->GetDefinition()->GetPDGMass())
179 { // Normal projectile
180 if( absPrPDGcode > 1000 ) //------Projectile is baryon --------
181 {
182 if ( absPrPDGcode > 4000 && absPrPDGcode < 6000 ) // Projectile is a charm or bottom baryon
183 {
184 MinPrDiffMass = projectile->GetDefinition()->GetPDGMass()/CLHEP::GeV + 0.25; // GeV
185 AveragePt2 = 0.3; // GeV^2
186 }
187 else
188 {
189 MinPrDiffMass = 1.16; // GeV
190 AveragePt2 = 0.3; // GeV^2
191 }
192 }
193 else if( absPrPDGcode == 211 || PrPDGcode == 111) //------Projectile is Pion -----------
194 {
195 MinPrDiffMass = 1.0; // GeV
196 AveragePt2 = 0.3; // GeV^2
197 }
198 else if( absPrPDGcode == 321 || absPrPDGcode == 130 || absPrPDGcode == 310) //-Projectile is Kaon-
199 {
200 MinPrDiffMass = 1.1; // GeV
201 AveragePt2 = 0.3; // GeV^2
202 }
203 else if( absPrPDGcode > 400 && absPrPDGcode < 600) // Projectile is a charm or bottom meson
204 {
205 MinPrDiffMass = projectile->GetDefinition()->GetPDGMass()/CLHEP::GeV + 0.25; // GeV
206 AveragePt2 = 0.3; // GeV^2
207 }
208 else //------Projectile is undefined, Nucleon assumed
209 {
210 MinPrDiffMass = 1.16; // GeV
211 AveragePt2 = 0.3; // GeV^2
212 }
213 }
214 else
215 { // Excited projectile
216 MinPrDiffMass = M0projectile + 220.0*MeV;
217 AveragePt2 = 0.3;
218 }
219
220 MinPrDiffMass = MinPrDiffMass * GeV;
221 AveragePt2 = AveragePt2 * GeV*GeV;
222 //---------------------------------------------
223 G4double MinTrDiffMass = 1.16*GeV;
224
225 if (SqrtS < MinPrDiffMass + MinTrDiffMass) {return false;} // The model cannot work at low energy
226
227 G4double MinPrDiffMass2 = MinPrDiffMass * MinPrDiffMass;
228 G4double MinTrDiffMass2 = MinTrDiffMass * MinTrDiffMass;
229
230 G4double Pt2;
231 G4double ProjMassT2, ProjMassT;
232 G4double TargMassT2, TargMassT;
233 G4double PMinusNew, TPlusNew;
234
235 G4LorentzVector Qmomentum;
236 G4double Qminus, Qplus;
237
238 G4int whilecount=0;
239 do {
240 if (whilecount++ >= 500 && (whilecount%100)==0)
241 if (whilecount > 1000 ) {
242 Qmomentum=G4LorentzVector(0.,0.,0.,0.);
243 return false; // Ignore this interaction
244 }
245
246 // Generate pt
247 Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0);
248
249 Pt2=G4ThreeVector(Qmomentum.vect()).mag2();
250 ProjMassT2=MinPrDiffMass2+Pt2;
251 ProjMassT =std::sqrt(ProjMassT2);
252
253 TargMassT2=MinTrDiffMass2+Pt2;
254 TargMassT =std::sqrt(TargMassT2);
255
256 if (SqrtS < ProjMassT + TargMassT) continue;
257
258 PZcms2=(S*S+ProjMassT2*ProjMassT2+TargMassT2*TargMassT2-
259 2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)/4./S;
260 if (PZcms2 < 0 ) {PZcms2=0;};
261 PZcms =std::sqrt(PZcms2);
262
263 G4double PMinusMin=std::sqrt(ProjMassT2+PZcms2)-PZcms;
264 G4double PMinusMax=SqrtS-TargMassT;
265
266 PMinusNew=ChooseP(PMinusMin,PMinusMax);
267 Qminus=PMinusNew-Pprojectile.minus();
268
269 G4double TPlusMin=std::sqrt(TargMassT2+PZcms2)-PZcms;
270 G4double TPlusMax=SqrtS-ProjMassT;
271
272 TPlusNew=ChooseP(TPlusMin, TPlusMax);
273 Qplus=-(TPlusNew-Ptarget.plus());
274
275 Qmomentum.setPz( (Qplus-Qminus)/2 );
276 Qmomentum.setE( (Qplus+Qminus)/2 );
277
278 } while ( (Pprojectile+Qmomentum).mag2() < MinPrDiffMass2 ||
279 (Ptarget -Qmomentum).mag2() < MinTrDiffMass2 );
280
281 Pprojectile += Qmomentum;
282 Ptarget -= Qmomentum;
283
284 // Transform back and update SplitableHadron Participant.
285 Pprojectile.transform(toLab);
286 Ptarget.transform(toLab);
287
288 #ifdef debugDoubleDiffraction
289 G4cout << "Pprojectile after boost to Lab: " << Pprojectile <<" "<<Pprojectile.mag()<<G4endl;
290 G4cout << "Ptarget after boost to Lab: " << Ptarget <<" "<<Ptarget.mag() <<G4endl;
291 #endif
292
293 target->Set4Momentum(Ptarget);
294 projectile->Set4Momentum(Pprojectile);
295
296 return true;
297}
G4double S(G4double temp)
CLHEP::HepLorentzVector G4LorentzVector
static constexpr double GeV
Definition: G4SIunits.hh:203
static constexpr double MeV
Definition: G4SIunits.hh:200
CLHEP::Hep3Vector G4ThreeVector
bool G4bool
Definition: G4Types.hh:86
int G4int
Definition: G4Types.hh:85
double mag2() const
double theta() const
Hep3Vector boostVector() const
Hep3Vector vect() const
HepLorentzVector & rotateZ(double)
HepLorentzVector & rotateY(double)
HepLorentzVector & transform(const HepRotation &)
G4double ChooseP(G4double Pmin, G4double Pmax) const
G4ThreeVector GaussianPt(G4double AveragePt2, G4double maxPtSquare) const
static constexpr double GeV
T sqr(const T &x)
Definition: templates.hh:128

References CLHEP::HepLorentzVector::boostVector(), ChooseP(), G4cout, G4endl, GaussianPt(), G4VSplitableHadron::Get4Momentum(), G4VSplitableHadron::GetDefinition(), G4ParticleDefinition::GetParticleName(), G4ParticleDefinition::GetPDGEncoding(), G4ParticleDefinition::GetPDGMass(), CLHEP::GeV, GeV, CLHEP::HepLorentzRotation::inverse(), CLHEP::HepLorentzVector::mag(), CLHEP::HepLorentzVector::mag2(), CLHEP::Hep3Vector::mag2(), MeV, CLHEP::HepLorentzVector::minus(), CLHEP::HepLorentzVector::phi(), CLHEP::HepLorentzVector::plus(), CLHEP::HepLorentzVector::pz(), CLHEP::HepLorentzRotation::rotateY(), CLHEP::HepLorentzRotation::rotateZ(), S(), G4VSplitableHadron::Set4Momentum(), CLHEP::HepLorentzVector::setE(), CLHEP::HepLorentzVector::setPz(), sqr(), CLHEP::HepLorentzVector::theta(), CLHEP::HepLorentzVector::transform(), CLHEP::HepLorentzVector::vect(), CLHEP::HepLorentzVector::x(), CLHEP::HepLorentzVector::y(), and CLHEP::HepLorentzVector::z().

Referenced by G4QGSParticipants::PerformDiffractiveCollisions().

◆ GaussianPt()

G4ThreeVector G4QGSDiffractiveExcitation::GaussianPt ( G4double  AveragePt2,
G4double  maxPtSquare 
) const
private

Definition at line 395 of file G4QGSDiffractiveExcitation.cc.

396{ // @@ this method is used in FTFModel as well. Should go somewhere common!
397 G4double Pt2;
398
399 Pt2 = -AveragePt2 * G4Log(1. + G4UniformRand() * (G4Exp(-maxPtSquare/AveragePt2)-1.));
400
401 G4double Pt=std::sqrt(Pt2);
402
404
405 return G4ThreeVector (Pt*std::cos(phi), Pt*std::sin(phi), 0.);
406}
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:179
G4double G4Log(G4double x)
Definition: G4Log.hh:226
static constexpr double twopi
Definition: G4SIunits.hh:56

References G4Exp(), G4Log(), G4UniformRand, and twopi.

Referenced by ExciteParticipants(), and String().

◆ operator!=()

G4bool G4QGSDiffractiveExcitation::operator!= ( const G4QGSDiffractiveExcitation right) const
private

◆ operator=()

const G4QGSDiffractiveExcitation & G4QGSDiffractiveExcitation::operator= ( const G4QGSDiffractiveExcitation right)
private

◆ operator==()

G4bool G4QGSDiffractiveExcitation::operator== ( const G4QGSDiffractiveExcitation right) const
private

◆ String()

G4ExcitedString * G4QGSDiffractiveExcitation::String ( G4VSplitableHadron aHadron,
G4bool  isProjectile 
) const
virtual

Definition at line 300 of file G4QGSDiffractiveExcitation.cc.

302{
303 hadron->SplitUp();
304 G4Parton *start= hadron->GetNextParton();
305 if ( start==NULL)
306 { G4cout << " G4QGSDiffractiveExcitation::String() Error:No start parton found"<< G4endl;
307 return NULL;
308 }
309 G4Parton *end = hadron->GetNextParton();
310 if ( end==NULL)
311 { G4cout << " G4QGSDiffractiveExcitation::String() Error:No end parton found"<< G4endl;
312 return NULL;
313 }
314
315 G4ExcitedString * string;
316 if ( isProjectile )
317 {
318 string= new G4ExcitedString(end,start, +1);
319 } else {
320 string= new G4ExcitedString(start,end, -1);
321 }
322
323 string->SetPosition(hadron->GetPosition());
324
325 // momenta of string ends
326
327 G4double maxAvailMomentumSquared=sqr(hadron->Get4Momentum().mag()/2.);
328
329 G4double widthOfPtSquare = 0.5*sqr(GeV);
330 G4ThreeVector pt=GaussianPt(widthOfPtSquare,maxAvailMomentumSquared);
331
332 G4LorentzVector Pstart(G4LorentzVector(pt,0.));
333 G4LorentzVector Pend;
334 Pend.setPx(hadron->Get4Momentum().px() - pt.x());
335 Pend.setPy(hadron->Get4Momentum().py() - pt.y());
336
337 G4double tm1=hadron->Get4Momentum().minus() +
338 ( Pend.perp2()-Pstart.perp2() ) / hadron->Get4Momentum().plus();
339
340 G4double tm2= std::sqrt( std::max(0., sqr(tm1) -
341 4. * Pend.perp2() * hadron->Get4Momentum().minus()
342 / hadron->Get4Momentum().plus() ));
343
344 G4int Sign= isProjectile ? -1 : 1;
345
346 G4double endMinus = 0.5 * (tm1 + Sign*tm2);
347 G4double startMinus= hadron->Get4Momentum().minus() - endMinus;
348
349 G4double startPlus= Pstart.perp2() / startMinus;
350 G4double endPlus = hadron->Get4Momentum().plus() - startPlus;
351
352 Pstart.setPz(0.5*(startPlus - startMinus));
353 Pstart.setE(0.5*(startPlus + startMinus));
354
355 Pend.setPz(0.5*(endPlus - endMinus));
356 Pend.setE(0.5*(endPlus + endMinus));
357
358 start->Set4Momentum(Pstart);
359 end->Set4Momentum(Pend);
360
361 #ifdef debugQGSdiffExictation
362 G4cout << " generated string flavors " << start->GetPDGcode() << " / " << end->GetPDGcode() << G4endl;
363 G4cout << " generated string momenta: quark " << start->Get4Momentum() << "mass : " <<start->Get4Momentum().mag()<< G4endl;
364 G4cout << " generated string momenta: Diquark " << end ->Get4Momentum() << "mass : " <<end->Get4Momentum().mag()<< G4endl;
365 G4cout << " sum of ends " << Pstart+Pend << G4endl;
366 G4cout << " Original " << hadron->Get4Momentum() << G4endl;
367 #endif
368
369 return string;
370}
double x() const
double y() const
double perp2() const
G4int GetPDGcode() const
Definition: G4Parton.hh:127
void Set4Momentum(const G4LorentzVector &aMomentum)
Definition: G4Parton.hh:148
const G4LorentzVector & Get4Momentum() const
Definition: G4Parton.hh:143
T max(const T t1, const T t2)
brief Return the largest of the two arguments

References G4cout, G4endl, GaussianPt(), G4VSplitableHadron::Get4Momentum(), G4Parton::Get4Momentum(), G4VSplitableHadron::GetNextParton(), G4Parton::GetPDGcode(), G4VSplitableHadron::GetPosition(), GeV, CLHEP::HepLorentzVector::mag(), G4INCL::Math::max(), CLHEP::HepLorentzVector::minus(), CLHEP::HepLorentzVector::perp2(), CLHEP::HepLorentzVector::plus(), CLHEP::HepLorentzVector::px(), CLHEP::HepLorentzVector::py(), G4Parton::Set4Momentum(), CLHEP::HepLorentzVector::setE(), CLHEP::HepLorentzVector::setPx(), CLHEP::HepLorentzVector::setPy(), CLHEP::HepLorentzVector::setPz(), G4VSplitableHadron::SplitUp(), sqr(), CLHEP::Hep3Vector::x(), and CLHEP::Hep3Vector::y().


The documentation for this class was generated from the following files: