Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
Public Member Functions | List of all members
G4QGSDiffractiveExcitation Class Reference

#include <G4QGSDiffractiveExcitation.hh>

+ Inheritance diagram for G4QGSDiffractiveExcitation:

Public Member Functions

 G4QGSDiffractiveExcitation ()
 
virtual ~G4QGSDiffractiveExcitation ()
 
virtual G4bool ExciteParticipants (G4VSplitableHadron *aPartner, G4VSplitableHadron *bPartner) const
 
virtual G4ExcitedString * String (G4VSplitableHadron *aHadron, G4bool isProjectile) const
 

Detailed Description

Definition at line 51 of file G4QGSDiffractiveExcitation.hh.

Constructor & Destructor Documentation

◆ G4QGSDiffractiveExcitation()

G4QGSDiffractiveExcitation::G4QGSDiffractiveExcitation ( )

Definition at line 59 of file G4QGSDiffractiveExcitation.cc.

60{
61}

◆ ~G4QGSDiffractiveExcitation()

G4QGSDiffractiveExcitation::~G4QGSDiffractiveExcitation ( )
virtual

Definition at line 63 of file G4QGSDiffractiveExcitation.cc.

64{
65}

Member Function Documentation

◆ ExciteParticipants()

G4bool G4QGSDiffractiveExcitation::ExciteParticipants ( G4VSplitableHadron *  aPartner,
G4VSplitableHadron *  bPartner 
) const
virtual

Reimplemented in G4SingleDiffractiveExcitation.

Definition at line 68 of file G4QGSDiffractiveExcitation.cc.

70{
71
72 G4LorentzVector Pprojectile=projectile->Get4Momentum();
73
74 // -------------------- Projectile parameters -----------------------------------
75 G4bool PutOnMassShell=0;
76
77 // G4double M0projectile=projectile->GetDefinition()->GetPDGMass(); // With de-excitation
78 G4double M0projectile = Pprojectile.mag(); // Without de-excitation
79
80 if(M0projectile < projectile->GetDefinition()->GetPDGMass())
81 {
82 PutOnMassShell=1;
83 M0projectile=projectile->GetDefinition()->GetPDGMass();
84 }
85
86 G4double Mprojectile2 = M0projectile * M0projectile;
87
88 G4int PDGcode=projectile->GetDefinition()->GetPDGEncoding();
89 G4int absPDGcode=std::abs(PDGcode);
90 G4double ProjectileDiffCut;
91 G4double AveragePt2;
92
93 if( absPDGcode > 1000 ) //------Projectile is baryon --------
94 {
95 ProjectileDiffCut = 1.1; // GeV
96 AveragePt2 = 0.3; // GeV^2
97 }
98 else if( absPDGcode == 211 || PDGcode == 111) //------Projectile is Pion -----------
99 {
100 ProjectileDiffCut = 1.0; // GeV
101 AveragePt2 = 0.3; // GeV^2
102 }
103 else if( absPDGcode == 321 || PDGcode == -311) //------Projectile is Kaon -----------
104 {
105 ProjectileDiffCut = 1.1; // GeV
106 AveragePt2 = 0.3; // GeV^2
107 }
108 else //------Projectile is undefined, Nucleon assumed
109 {
110 ProjectileDiffCut = 1.1; // GeV
111 AveragePt2 = 0.3; // GeV^2
112 };
113
114 ProjectileDiffCut = ProjectileDiffCut * GeV;
115 AveragePt2 = AveragePt2 * GeV*GeV;
116
117 // -------------------- Target parameters ----------------------------------------------
118 G4LorentzVector Ptarget=target->Get4Momentum();
119
120 G4double M0target = Ptarget.mag();
121
122 if(M0target < target->GetDefinition()->GetPDGMass())
123 {
124 PutOnMassShell=1;
125 M0target=target->GetDefinition()->GetPDGMass();
126 }
127
128 G4double Mtarget2 = M0target * M0target; //Ptarget.mag2(); // for AA-inter.
129
130 G4double NuclearNucleonDiffCut = 1.1*GeV;
131
132 G4double ProjectileDiffCut2 = ProjectileDiffCut * ProjectileDiffCut;
133 G4double NuclearNucleonDiffCut2 = NuclearNucleonDiffCut * NuclearNucleonDiffCut;
134
135 // Transform momenta to cms and then rotate parallel to z axis;
136
137 G4LorentzVector Psum;
138 Psum=Pprojectile+Ptarget;
139
140 G4LorentzRotation toCms(-1*Psum.boostVector());
141
142 G4LorentzVector Ptmp=toCms*Pprojectile;
143
144 if ( Ptmp.pz() <= 0. )
145 {
146 // "String" moving backwards in CMS, abort collision !!
147 //G4cout << " abort Collision!! " << G4endl;
148 return false;
149 }
150
151 toCms.rotateZ(-1*Ptmp.phi());
152 toCms.rotateY(-1*Ptmp.theta());
153
154 G4LorentzRotation toLab(toCms.inverse());
155
156 Pprojectile.transform(toCms);
157 Ptarget.transform(toCms);
158
159 G4double Pt2;
160 G4double ProjMassT2, ProjMassT;
161 G4double TargMassT2, TargMassT;
162 G4double PZcms2, PZcms;
163 G4double PMinusNew, TPlusNew;
164
165 G4double S=Psum.mag2();
166 G4double SqrtS=std::sqrt(S);
167
168 if(SqrtS < 2200*MeV) {return false;} // The model cannot work for pp-interactions
169 // at Plab < 1.3 GeV/c. Uzhi
170
171 PZcms2=(S*S+Mprojectile2*Mprojectile2+Mtarget2*Mtarget2-
172 2*S*Mprojectile2-2*S*Mtarget2-2*Mprojectile2*Mtarget2)/4./S;
173 if(PZcms2 < 0)
174 {return false;} // It can be in an interaction with off-shell nuclear nucleon
175
176 PZcms = std::sqrt(PZcms2);
177
178 if(PutOnMassShell)
179 {
180 if(Pprojectile.z() > 0.)
181 {
182 Pprojectile.setPz( PZcms);
183 Ptarget.setPz( -PZcms);
184 }
185 else
186 {
187 Pprojectile.setPz(-PZcms);
188 Ptarget.setPz( PZcms);
189 };
190
191 Pprojectile.setE(std::sqrt(Mprojectile2+
192 Pprojectile.x()*Pprojectile.x()+
193 Pprojectile.y()*Pprojectile.y()+
194 PZcms2));
195 Ptarget.setE(std::sqrt( Mtarget2 +
196 Ptarget.x()*Ptarget.x()+
197 Ptarget.y()*Ptarget.y()+
198 PZcms2));
199 }
200
201 G4double maxPtSquare = PZcms2;
202
203 //G4cout << "Pprojectile aft boost : " << Pprojectile << G4endl;
204 //G4cout << "Ptarget aft boost : " << Ptarget << G4endl;
205 // G4cout << "cms aft boost : " << (Pprojectile+ Ptarget) << G4endl;
206
207 // G4cout << " Projectile Xplus / Xminus : " <<
208 // Pprojectile.plus() << " / " << Pprojectile.minus() << G4endl;
209 // G4cout << " Target Xplus / Xminus : " <<
210 // Ptarget.plus() << " / " << Ptarget.minus() << G4endl;
211
212 G4LorentzVector Qmomentum;
213 G4double Qminus, Qplus;
214
215 // /* Vova
216 G4int whilecount=0;
217 do {
218 // Generate pt
219
220 if (whilecount++ >= 500 && (whilecount%100)==0)
221 // G4cout << "G4QGSDiffractiveExcitation::ExciteParticipants possibly looping"
222 // << ", loop count/ maxPtSquare : "
223 // << whilecount << " / " << maxPtSquare << G4endl;
224 if (whilecount > 1000 )
225 {
226 Qmomentum=G4LorentzVector(0.,0.,0.,0.);
227 return false; // Ignore this interaction
228 }
229
230 Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0);
231
232 //G4cout << "generated Pt " << Qmomentum << G4endl;
233 //G4cout << "Pprojectile with pt : " << Pprojectile+Qmomentum << G4endl;
234 //G4cout << "Ptarget with pt : " << Ptarget-Qmomentum << G4endl;
235
236 // Momentum transfer
237 /* // Uzhi
238 G4double Xmin = minmass / ( Pprojectile.e() + Ptarget.e() );
239 G4double Xmax=1.;
240 G4double Xplus =ChooseX(Xmin,Xmax);
241 G4double Xminus=ChooseX(Xmin,Xmax);
242
243// G4cout << " X-plus " << Xplus << G4endl;
244// G4cout << " X-minus " << Xminus << G4endl;
245
246 G4double pt2=G4ThreeVector(Qmomentum.vect()).mag2();
247 G4double Qplus =-1 * pt2 / Xminus/Ptarget.minus();
248 G4double Qminus= pt2 / Xplus /Pprojectile.plus();
249 */ // Uzhi *
250
251 Pt2=G4ThreeVector(Qmomentum.vect()).mag2();
252 ProjMassT2=Mprojectile2+Pt2;
253 ProjMassT =std::sqrt(ProjMassT2);
254
255 TargMassT2=Mtarget2+Pt2;
256 TargMassT =std::sqrt(TargMassT2);
257
258 PZcms2=(S*S+ProjMassT2*ProjMassT2+
259 TargMassT2*TargMassT2-
260 2.*S*ProjMassT2-2.*S*TargMassT2-
261 2.*ProjMassT2*TargMassT2)/4./S;
262 if(PZcms2 < 0 ) {PZcms2=0;};
263 PZcms =std::sqrt(PZcms2);
264
265 G4double PMinusMin=std::sqrt(ProjMassT2+PZcms2)-PZcms;
266 G4double PMinusMax=SqrtS-TargMassT;
267
268 PMinusNew=ChooseP(PMinusMin,PMinusMax);
269 Qminus=PMinusNew-Pprojectile.minus();
270
271 G4double TPlusMin=std::sqrt(TargMassT2+PZcms2)-PZcms;
272 G4double TPlusMax=SqrtS-ProjMassT;
273
274 TPlusNew=ChooseP(TPlusMin, TPlusMax);
275 Qplus=-(TPlusNew-Ptarget.plus());
276
277 Qmomentum.setPz( (Qplus-Qminus)/2 );
278 Qmomentum.setE( (Qplus+Qminus)/2 );
279
280 //G4cout << "Qplus / Qminus " << Qplus << " / " << Qminus<<G4endl;
281 // G4cout << "pt2" << pt2 << G4endl;
282 // G4cout << "Qmomentum " << Qmomentum << G4endl;
283 // G4cout << " Masses (P/T) : " << (Pprojectile+Qmomentum).mag() <<
284 // " / " << (Ptarget-Qmomentum).mag() << G4endl;
285 /* // Uzhi
286 } while ( (Pprojectile+Qmomentum).mag2() <= Mprojectile2 ||
287 (Ptarget-Qmomentum).mag2() <= Mtarget2 );
288 */ // Uzhi *
289
290
291 } while (( (Pprojectile+Qmomentum).mag2() < Mprojectile2 || // Uzhi No without excitation
292 (Ptarget -Qmomentum).mag2() < Mtarget2 ) || // Uzhi
293 ( (Pprojectile+Qmomentum).mag2() < ProjectileDiffCut2 && // Uzhi No double Diffraction
294 (Ptarget -Qmomentum).mag2() < NuclearNucleonDiffCut2) );// Uzhi
295
296 if((Ptarget-Qmomentum).mag2() < NuclearNucleonDiffCut2) // Uzhi Projectile diffraction
297 {
298 G4double TMinusNew=SqrtS-PMinusNew;
299 Qminus=Ptarget.minus()-TMinusNew;
300 TPlusNew=TargMassT2/TMinusNew;
301 Qplus=Ptarget.plus()-TPlusNew;
302
303 Qmomentum.setPz( (Qplus-Qminus)/2 );
304 Qmomentum.setE( (Qplus+Qminus)/2 );
305 }
306 else if((Pprojectile+Qmomentum).mag2() < ProjectileDiffCut2) // Uzhi Target diffraction
307 {
308 G4double PPlusNew=SqrtS-TPlusNew;
309 Qplus=PPlusNew-Pprojectile.plus();
310 PMinusNew=ProjMassT2/PPlusNew;
311 Qminus=PMinusNew-Pprojectile.minus();
312
313 Qmomentum.setPz( (Qplus-Qminus)/2 );
314 Qmomentum.setE( (Qplus+Qminus)/2 );
315 };
316
317 Pprojectile += Qmomentum;
318 Ptarget -= Qmomentum;
319
320 // Vova
321
322 /*
323 Pprojectile.setPz(0.);
324 Pprojectile.setE(SqrtS-M0target);
325
326 Ptarget.setPz(0.);
327 Ptarget.setE(M0target);
328 */
329
330 //G4cout << "Pprojectile with Q : " << Pprojectile << G4endl;
331 //G4cout << "Ptarget with Q : " << Ptarget << G4endl;
332
333 // G4cout << "Projectile back: " << toLab * Pprojectile << G4endl;
334 // G4cout << "Target back: " << toLab * Ptarget << G4endl;
335
336 // Transform back and update SplitableHadron Participant.
337 Pprojectile.transform(toLab);
338 Ptarget.transform(toLab);
339
340 //G4cout << "Pprojectile with Q M: " << Pprojectile<<" "<< Pprojectile.mag() << G4endl;
341 //G4cout << "Ptarget with Q M: " << Ptarget <<" "<< Ptarget.mag() << G4endl;
342
343 //G4cout << "Target mass " << Ptarget.mag() << G4endl;
344
345 target->Set4Momentum(Ptarget);
346
347 //G4cout << "Projectile mass " << Pprojectile.mag() << G4endl;
348
349 projectile->Set4Momentum(Pprojectile);
350
351 return true;
352}
G4LorentzVector
CLHEP::HepLorentzVector G4LorentzVector
Definition: G4LorentzVector.hh:35
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
G4double
double G4double
Definition: G4Types.hh:64
G4int
int G4int
Definition: G4Types.hh:66
G4bool
bool G4bool
Definition: G4Types.hh:67
CLHEP::Hep3Vector::mag2
double mag2() const
CLHEP::HepLorentzVector::theta
double theta() const
CLHEP::HepLorentzVector::boostVector
Hep3Vector boostVector() const
Definition: LorentzVector.cc:177
CLHEP::HepLorentzVector::vect
Hep3Vector vect() const
CLHEP::HepLorentzVector::rotateZ
HepLorentzVector & rotateZ(double)
Definition: LorentzVector.cc:80
CLHEP::HepLorentzVector::rotateY
HepLorentzVector & rotateY(double)
Definition: LorentzVector.cc:76
CLHEP::HepLorentzVector::transform
HepLorentzVector & transform(const HepRotation &)

Referenced by G4QGSParticipants::SelectInteractions().

◆ String()

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

Definition at line 355 of file G4QGSDiffractiveExcitation.cc.

357{
358 hadron->SplitUp();
359 G4Parton *start= hadron->GetNextParton();
360 if ( start==NULL)
361 { G4cout << " G4FTFModel::String() Error:No start parton found"<< G4endl;
362 return NULL;
363 }
364 G4Parton *end = hadron->GetNextParton();
365 if ( end==NULL)
366 { G4cout << " G4FTFModel::String() Error:No end parton found"<< G4endl;
367 return NULL;
368 }
369
370 G4ExcitedString * string;
371 if ( isProjectile )
372 {
373 string= new G4ExcitedString(end,start, +1);
374 } else {
375 string= new G4ExcitedString(start,end, -1);
376 }
377
378 string->SetPosition(hadron->GetPosition());
379
380 // momenta of string ends
381 G4double ptSquared= hadron->Get4Momentum().perp2();
382 G4double transverseMassSquared= hadron->Get4Momentum().plus()
383 * hadron->Get4Momentum().minus();
384
385
386 G4double maxAvailMomentumSquared=
387 sqr( std::sqrt(transverseMassSquared) - std::sqrt(ptSquared) );
388
389 G4double widthOfPtSquare = 0.25; // Uzhi <Pt^2>=0.25 ??????????????????
390 G4ThreeVector pt=GaussianPt(widthOfPtSquare,maxAvailMomentumSquared);
391
392 G4LorentzVector Pstart(G4LorentzVector(pt,0.));
393 G4LorentzVector Pend;
394 Pend.setPx(hadron->Get4Momentum().px() - pt.x());
395 Pend.setPy(hadron->Get4Momentum().py() - pt.y());
396
397 G4double tm1=hadron->Get4Momentum().minus() +
398 ( Pend.perp2()-Pstart.perp2() ) / hadron->Get4Momentum().plus();
399
400 G4double tm2= std::sqrt( std::max(0., sqr(tm1) -
401 4. * Pend.perp2() * hadron->Get4Momentum().minus()
402 / hadron->Get4Momentum().plus() ));
403
404 G4int Sign= isProjectile ? -1 : 1;
405
406 G4double endMinus = 0.5 * (tm1 + Sign*tm2);
407 G4double startMinus= hadron->Get4Momentum().minus() - endMinus;
408
409 G4double startPlus= Pstart.perp2() / startMinus;
410 G4double endPlus = hadron->Get4Momentum().plus() - startPlus;
411
412 Pstart.setPz(0.5*(startPlus - startMinus));
413 Pstart.setE(0.5*(startPlus + startMinus));
414
415 Pend.setPz(0.5*(endPlus - endMinus));
416 Pend.setE(0.5*(endPlus + endMinus));
417
418 start->Set4Momentum(Pstart);
419 end->Set4Momentum(Pend);
420
421 #ifdef G4_FTFDEBUG
422 G4cout << " generated string flavors " << start->GetPDGcode() << " / " << end->GetPDGcode() << G4endl;
423 G4cout << " generated string momenta: quark " << start->Get4Momentum() << "mass : " <<start->Get4Momentum().mag()<< G4endl;
424 G4cout << " generated string momenta: Diquark " << end ->Get4Momentum() << "mass : " <<end->Get4Momentum().mag()<< G4endl;
425 G4cout << " sum of ends " << Pstart+Pend << G4endl;
426 G4cout << " Original " << hadron->Get4Momentum() << G4endl;
427 #endif
428
429 return string;
430}
G4endl
#define G4endl
Definition: G4ios.hh:52
G4cout
G4DLLIMPORT std::ostream G4cout
CLHEP::Hep3Vector::x
double x() const
CLHEP::Hep3Vector::y
double y() const
CLHEP::HepLorentzVector::perp2
double perp2() const
G4Parton::GetPDGcode
G4int GetPDGcode() const
Definition: G4Parton.hh:124
G4Parton::Set4Momentum
void Set4Momentum(const G4LorentzVector &aMomentum)
Definition: G4Parton.hh:145
G4Parton::Get4Momentum
const G4LorentzVector & Get4Momentum() const
Definition: G4Parton.hh:140
sqr
T sqr(const T &x)
Definition: templates.hh:145
The documentation for this class was generated from the following files: