81 #ifdef debugDoubleDiffraction
82 G4cout<<
G4endl<<
"G4QGSDiffractiveExcitation::ExciteParticipants - Double diffraction."<<
G4endl;
97 if(M0projectile < projectile->GetDefinition()->GetPDGMass())
108 if(M0target < target->GetDefinition()->GetPDGMass())
118 if(SqrtS < M0projectile + M0target) {
return false;}
121 G4double Mprojectile2 = M0projectile * M0projectile;
122 G4double Mtarget2 = M0target * M0target;
130 if ( Ptmp.pz() <= 0. )
138 toCms.
rotateY(-1*Ptmp.theta());
145 G4double PZcms2=(S*S+Mprojectile2*Mprojectile2+Mtarget2*Mtarget2-
146 2*S*Mprojectile2-2*S*Mtarget2-2*Mprojectile2*Mtarget2)/4./S;
148 if(PZcms2 < 0) {
return false;}
154 if(Pprojectile.z() > 0.)
156 Pprojectile.setPz( PZcms);
157 Ptarget.setPz( -PZcms);
161 Pprojectile.setPz(-PZcms);
162 Ptarget.setPz( PZcms);
165 Pprojectile.setE(std::sqrt(Mprojectile2+
sqr(Pprojectile.x())+
sqr(Pprojectile.y())+PZcms2));
166 Ptarget.setE( std::sqrt( Mtarget2+
sqr( Ptarget.x())+
sqr( Ptarget.y())+PZcms2));
171 #ifdef debugDoubleDiffraction
172 G4cout <<
"Pprojectile after boost to CMS: " << Pprojectile <<
" "<<Pprojectile.mag()<<
G4endl;
173 G4cout <<
"Ptarget after boost to CMS: " << Ptarget <<
" "<<Ptarget.mag() <<
G4endl;
177 G4int absPrPDGcode=std::abs(PrPDGcode);
181 if(M0projectile <= projectile->GetDefinition()->GetPDGMass())
183 if( absPrPDGcode > 1000 )
185 MinPrDiffMass = 1.16;
188 else if( absPrPDGcode == 211 || PrPDGcode == 111)
193 else if( absPrPDGcode == 321 || absPrPDGcode == 130 || absPrPDGcode == 310)
200 MinPrDiffMass = 1.16;
206 MinPrDiffMass = M0projectile + 220.0*
MeV;
210 MinPrDiffMass = MinPrDiffMass *
GeV;
211 AveragePt2 = AveragePt2 * GeV*
GeV;
215 if(SqrtS < MinPrDiffMass + MinTrDiffMass) {
return false;}
217 G4double MinPrDiffMass2 = MinPrDiffMass * MinPrDiffMass;
218 G4double MinTrDiffMass2 = MinTrDiffMass * MinTrDiffMass;
230 if (whilecount++ >= 500 && (whilecount%100)==0)
234 if (whilecount > 1000 )
244 ProjMassT2=MinPrDiffMass2+Pt2;
245 ProjMassT =std::sqrt(ProjMassT2);
247 TargMassT2=MinTrDiffMass2+Pt2;
248 TargMassT =std::sqrt(TargMassT2);
250 if(SqrtS < ProjMassT + TargMassT)
continue;
252 PZcms2=(S*S+ProjMassT2*ProjMassT2+
253 TargMassT2*TargMassT2-
254 2.*S*ProjMassT2-2.*S*TargMassT2-
255 2.*ProjMassT2*TargMassT2)/4./S;
256 if(PZcms2 < 0 ) {PZcms2=0;};
257 PZcms =std::sqrt(PZcms2);
259 G4double PMinusMin=std::sqrt(ProjMassT2+PZcms2)-PZcms;
262 PMinusNew=
ChooseP(PMinusMin,PMinusMax);
263 Qminus=PMinusNew-Pprojectile.minus();
265 G4double TPlusMin=std::sqrt(TargMassT2+PZcms2)-PZcms;
268 TPlusNew=
ChooseP(TPlusMin, TPlusMax);
269 Qplus=-(TPlusNew-Ptarget.plus());
271 Qmomentum.
setPz( (Qplus-Qminus)/2 );
272 Qmomentum.
setE( (Qplus+Qminus)/2 );
274 }
while ( (Pprojectile+Qmomentum).mag2() < MinPrDiffMass2 ||
275 (Ptarget -Qmomentum).mag2() < MinTrDiffMass2 );
277 Pprojectile += Qmomentum;
278 Ptarget -= Qmomentum;
284 #ifdef debugDoubleDiffraction
285 G4cout <<
"Pprojectile after boost to Lab: " << Pprojectile <<
" "<<Pprojectile.mag()<<
G4endl;
286 G4cout <<
"Ptarget after boost to Lab: " << Ptarget <<
" "<<Ptarget.mag() <<
G4endl;
302 {
G4cout <<
" G4QGSDiffractiveExcitation::String() Error:No start parton found"<<
G4endl;
307 {
G4cout <<
" G4QGSDiffractiveExcitation::String() Error:No end parton found"<<
G4endl;
348 G4int Sign= isProjectile ? -1 : 1;
350 G4double endMinus = 0.5 * (tm1 + Sign*tm2);
356 Pstart.
setPz(0.5*(startPlus - startMinus));
357 Pstart.
setE(0.5*(startPlus + startMinus));
359 Pend.
setPz(0.5*(endPlus - endMinus));
360 Pend.
setE(0.5*(endPlus + endMinus));
365 #ifdef debugQGSdiffExictation
386 if ( Pmin <= 0. || range <=0. )
388 G4cout <<
" Pmin, range : " << Pmin <<
" , " << range <<
G4endl;
389 throw G4HadronicException(__FILE__, __LINE__,
"G4QGSDiffractiveExcitation::ChooseP : Invalid arguments ");
409 return G4ThreeVector (Pt*std::cos(phi), Pt*std::sin(phi), 0.);
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
T max(const T t1, const T t2)
brief Return the largest of the two arguments
const G4ParticleDefinition * GetDefinition() const
HepLorentzVector & transform(const HepRotation &)
CLHEP::Hep3Vector G4ThreeVector
static constexpr double MeV
const G4String & GetParticleName() const
G4double ChooseP(G4double Pmin, G4double Pmax) const
virtual G4Parton * GetNextParton()=0
G4double GetPDGMass() const
G4double G4Log(G4double x)
G4QGSDiffractiveExcitation()
static G4Pow * GetInstance()
virtual ~G4QGSDiffractiveExcitation()
const G4ThreeVector & GetPosition() const
G4double powA(G4double A, G4double y) const
const G4ParticleDefinition const G4Material *G4double range
static constexpr double twopi
HepLorentzVector & rotateY(double)
const G4ThreeVector const G4double const
virtual G4bool ExciteParticipants(G4VSplitableHadron *aPartner, G4VSplitableHadron *bPartner, G4bool ProjectileDiffraction=TRUE) const
G4ThreeVector GaussianPt(G4double AveragePt2, G4double maxPtSquare) const
HepLorentzVector & rotateZ(double)
G4int GetPDGEncoding() const
G4GLOB_DLL std::ostream G4cout
Hep3Vector boostVector() const
const G4LorentzVector & Get4Momentum() const
CLHEP::HepLorentzVector G4LorentzVector
const G4LorentzVector & Get4Momentum() const
static constexpr double GeV
virtual G4ExcitedString * String(G4VSplitableHadron *aHadron, G4bool isProjectile) const
HepLorentzRotation & transform(const HepBoost &b)
void Set4Momentum(const G4LorentzVector &a4Momentum)
void Set4Momentum(const G4LorentzVector &aMomentum)