80 #ifdef debugDoubleDiffraction
81 G4cout<<
G4endl<<
"G4QGSDiffractiveExcitation::ExciteParticipants - Double diffraction."<<
G4endl;
94 if(M0projectile < projectile->GetDefinition()->GetPDGMass())
105 if(M0target < target->GetDefinition()->GetPDGMass())
115 if (SqrtS < M0projectile + M0target) {
return false;}
118 G4double Mprojectile2 = M0projectile * M0projectile;
119 G4double Mtarget2 = M0target * M0target;
127 if ( Ptmp.
pz() <= 0. )
142 G4double PZcms2=(
S*
S+Mprojectile2*Mprojectile2+Mtarget2*Mtarget2-
143 2*
S*Mprojectile2-2*
S*Mtarget2-2*Mprojectile2*Mtarget2)/4./
S;
145 if (PZcms2 < 0) {
return false;}
151 if (Pprojectile.
z() > 0.)
153 Pprojectile.
setPz( PZcms);
154 Ptarget.
setPz( -PZcms);
158 Pprojectile.
setPz(-PZcms);
159 Ptarget.
setPz( PZcms);
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));
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;
174 G4int absPrPDGcode=std::abs(PrPDGcode);
178 if (M0projectile <= projectile->GetDefinition()->GetPDGMass())
180 if( absPrPDGcode > 1000 )
182 if ( absPrPDGcode > 4000 && absPrPDGcode < 6000 )
189 MinPrDiffMass = 1.16;
193 else if( absPrPDGcode == 211 || PrPDGcode == 111)
198 else if( absPrPDGcode == 321 || absPrPDGcode == 130 || absPrPDGcode == 310)
203 else if( absPrPDGcode > 400 && absPrPDGcode < 600)
210 MinPrDiffMass = 1.16;
216 MinPrDiffMass = M0projectile + 220.0*
MeV;
220 MinPrDiffMass = MinPrDiffMass *
GeV;
221 AveragePt2 = AveragePt2 *
GeV*
GeV;
225 if (SqrtS < MinPrDiffMass + MinTrDiffMass) {
return false;}
227 G4double MinPrDiffMass2 = MinPrDiffMass * MinPrDiffMass;
228 G4double MinTrDiffMass2 = MinTrDiffMass * MinTrDiffMass;
240 if (whilecount++ >= 500 && (whilecount%100)==0)
241 if (whilecount > 1000 ) {
250 ProjMassT2=MinPrDiffMass2+Pt2;
251 ProjMassT =std::sqrt(ProjMassT2);
253 TargMassT2=MinTrDiffMass2+Pt2;
254 TargMassT =std::sqrt(TargMassT2);
256 if (SqrtS < ProjMassT + TargMassT)
continue;
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);
263 G4double PMinusMin=std::sqrt(ProjMassT2+PZcms2)-PZcms;
266 PMinusNew=
ChooseP(PMinusMin,PMinusMax);
267 Qminus=PMinusNew-Pprojectile.
minus();
269 G4double TPlusMin=std::sqrt(TargMassT2+PZcms2)-PZcms;
272 TPlusNew=
ChooseP(TPlusMin, TPlusMax);
273 Qplus=-(TPlusNew-Ptarget.
plus());
275 Qmomentum.
setPz( (Qplus-Qminus)/2 );
276 Qmomentum.
setE( (Qplus+Qminus)/2 );
278 }
while ( (Pprojectile+Qmomentum).mag2() < MinPrDiffMass2 ||
279 (Ptarget -Qmomentum).mag2() < MinTrDiffMass2 );
281 Pprojectile += Qmomentum;
282 Ptarget -= Qmomentum;
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;
306 {
G4cout <<
" G4QGSDiffractiveExcitation::String() Error:No start parton found"<<
G4endl;
311 {
G4cout <<
" G4QGSDiffractiveExcitation::String() Error:No end parton found"<<
G4endl;
344 G4int Sign= isProjectile ? -1 : 1;
346 G4double endMinus = 0.5 * (tm1 + Sign*tm2);
352 Pstart.
setPz(0.5*(startPlus - startMinus));
353 Pstart.
setE(0.5*(startPlus + startMinus));
355 Pend.
setPz(0.5*(endPlus - endMinus));
356 Pend.
setE(0.5*(endPlus + endMinus));
361 #ifdef debugQGSdiffExictation
382 if (
Pmin <= 0. || range <=0. )
385 throw G4HadronicException(__FILE__, __LINE__,
"G4QGSDiffractiveExcitation::ChooseP : Invalid arguments ");
405 return G4ThreeVector (Pt*std::cos(phi), Pt*std::sin(phi), 0.);
G4double S(G4double temp)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
G4double G4Log(G4double x)
CLHEP::HepLorentzVector G4LorentzVector
static constexpr double twopi
static constexpr double GeV
static constexpr double MeV
CLHEP::Hep3Vector G4ThreeVector
G4GLOB_DLL std::ostream G4cout
HepLorentzRotation & rotateY(double delta)
HepLorentzRotation & rotateZ(double delta)
HepLorentzRotation inverse() const
Hep3Vector boostVector() const
HepLorentzVector & transform(const HepRotation &)
G4double GetPDGMass() const
G4int GetPDGEncoding() const
const G4String & GetParticleName() const
void Set4Momentum(const G4LorentzVector &aMomentum)
const G4LorentzVector & Get4Momentum() const
static G4Pow * GetInstance()
G4double powA(G4double A, G4double y) const
G4double ChooseP(G4double Pmin, G4double Pmax) const
G4QGSDiffractiveExcitation()
virtual ~G4QGSDiffractiveExcitation()
G4ThreeVector GaussianPt(G4double AveragePt2, G4double maxPtSquare) const
virtual G4bool ExciteParticipants(G4VSplitableHadron *aPartner, G4VSplitableHadron *bPartner, G4bool ProjectileDiffraction=TRUE) const
virtual G4ExcitedString * String(G4VSplitableHadron *aHadron, G4bool isProjectile) const
const G4ParticleDefinition * GetDefinition() const
void Set4Momentum(const G4LorentzVector &a4Momentum)
virtual G4Parton * GetNextParton()=0
const G4LorentzVector & Get4Momentum() const
const G4ThreeVector & GetPosition() const
static constexpr double GeV
T max(const T t1, const T t2)
brief Return the largest of the two arguments