61 #ifdef debugSingleDiffraction
73 #ifdef debugSingleDiffraction
75 G4cout<<
"Pr Tr 4-Mom "<<Pprojectile<<
" "<<Pprojectile.
mag()<<G4endl
76 <<
" "<<Ptarget <<
" "<<Ptarget.mag() <<
G4endl;
83 #ifdef debugSingleDiffraction
84 G4cout<<
"SqrtS-Mprojectile-Mtarget "<<SqrtS<<
" "<<Mprojectile<<
" "<<Mtarget
85 <<
" "<<SqrtS-Mprojectile-Mtarget<<
G4endl;
87 if(SqrtS-Mprojectile-Mtarget <= 250.0*
MeV) {
88 #ifdef debugSingleDiffraction
90 <<Pprojectile<<
" "<<Pprojectile.
mag()<<
G4endl;
92 <<Ptarget<<
" "<<Ptarget.mag()<<
G4endl;
93 G4cerr<<
"sqrt(S) = "<<SqrtS<<
" Mp + Mt = "<<Pprojectile.
mag()+Ptarget.mag()<<
G4endl;
102 if ( Ptmp.pz() <= 0. )
110 toCms.
rotateY(-1*Ptmp.theta());
116 #ifdef debugSingleDiffraction
117 G4cout <<
"Pprojectile in CMS " << Pprojectile <<
G4endl;
123 G4double ProjectileMinDiffrMass(0.), TargetMinDiffrMass(0.);
127 if ( ProjectileDiffraction ) {
128 if( absPDGcode > 1000 )
130 ProjectileMinDiffrMass = 1.16;
133 else if( absPDGcode == 211 || absPDGcode == 111)
135 ProjectileMinDiffrMass = 1.0;
138 else if( absPDGcode == 321 || absPDGcode == 130 || absPDGcode == 310)
140 ProjectileMinDiffrMass = 1.1;
143 else if( absPDGcode == 22)
145 ProjectileMinDiffrMass = 0.25;
150 ProjectileMinDiffrMass = 1.1;
154 ProjectileMinDiffrMass = ProjectileMinDiffrMass *
GeV;
155 Mprojectile2=
sqr(ProjectileMinDiffrMass);
159 TargetMinDiffrMass = 1.16*
GeV;
160 Mtarget2 =
sqr( TargetMinDiffrMass) ;
164 AveragePt2 = AveragePt2 *
GeV*
GeV;
173 G4double PMinusNew, PPlusNew, TPlusNew, TMinusNew;
182 if (whilecount > 1000 )
193 ProjMassT2 = Mprojectile2 + Pt2;
194 ProjMassT = std::sqrt( ProjMassT2 );
195 TargMassT2 = Mtarget2 + Pt2;
196 TargMassT = std::sqrt( TargMassT2 );
198 #ifdef debugSingleDiffraction
199 G4cout<<whilecount<<
" "<<Pt2<<
" "<<ProjMassT<<
" "<<TargMassT<<
" "<<SqrtS<<
" "<<S<<
" "<<ProjectileDiffraction<<
G4endl;
201 if ( SqrtS < ProjMassT + TargMassT )
continue;
203 PZcms2 = ( S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2
204 - 2.0*S*ProjMassT2 - 2.0*S*TargMassT2 - 2.0*ProjMassT2*TargMassT2 ) / 4.0 / S;
206 if ( PZcms2 < 0 )
continue;
208 PZcms = std::sqrt( PZcms2 );
210 if ( ProjectileDiffraction )
212 PMinusMin = std::sqrt( ProjMassT2 + PZcms2 ) - PZcms;
213 PMinusMax = SqrtS - TargMassT;
215 PMinusNew =
ChooseX( PMinusMin, PMinusMax );
216 TMinusNew = SqrtS - PMinusNew;
218 Qminus = Ptarget.minus() - TMinusNew;
219 TPlusNew = TargMassT2 / TMinusNew;
220 Qplus = Ptarget.plus() - TPlusNew;
223 TPlusMin = std::sqrt( TargMassT2 + PZcms2 ) - PZcms;
224 TPlusMax = SqrtS - ProjMassT;
226 TPlusNew =
ChooseX( TPlusMin, TPlusMax );
227 PPlusNew = SqrtS - TPlusNew;
229 Qplus = PPlusNew - Pprojectile.plus();
230 PMinusNew = ProjMassT2 / PPlusNew;
231 Qminus = PMinusNew - Pprojectile.minus();
235 Qmomentum.
setPz( (Qplus - Qminus)/2 );
236 Qmomentum.
setE( (Qplus + Qminus)/2 );
238 #ifdef debugSingleDiffraction
239 G4cout<<ProjectileDiffraction<<
" "<<( Pprojectile + Qmomentum ).mag2()<<
" "<< Mprojectile2<<
G4endl;
240 G4cout<<!ProjectileDiffraction<<
" "<<( Ptarget - Qmomentum ).mag2()<<
" "<< Mtarget2<<
G4endl;
243 }
while ( ( ProjectileDiffraction&&( Pprojectile + Qmomentum ).mag2() < Mprojectile2 ) ||
244 (!ProjectileDiffraction&&( Ptarget - Qmomentum ).mag2() < Mtarget2 ) );
248 Pprojectile += Qmomentum;
250 Ptarget -= Qmomentum;
256 #ifdef debugSingleDiffraction
257 G4cout <<
"Pprojectile in Lab. " << Pprojectile <<
G4endl;
259 G4cout <<
"G4SingleDiffractiveExcitation- Projectile mass " << Pprojectile.mag() <<
G4endl;
260 G4cout <<
"G4SingleDiffractiveExcitation- Target mass " << Ptarget.mag() <<
G4endl;
278 if ( Xmin <= 0. || range <=0. )
280 G4cout <<
" Xmin, range : " << Xmin <<
" , " << range <<
G4endl;
281 throw G4HadronicException(__FILE__, __LINE__,
"G4SingleDiffractiveExcitation::ChooseX : Invalid arguments ");
294 const G4int maxNumberOfLoops = 1000;
295 G4int loopCounter = 0;
298 }
while ( ( pt2 > maxPtSquare) && ++loopCounter < maxNumberOfLoops );
299 if ( loopCounter >= maxNumberOfLoops ) {
300 pt2 = 0.99*maxPtSquare;
307 return G4ThreeVector (pt2*std::cos(phi), pt2*std::sin(phi), 0.);
const G4ParticleDefinition * GetDefinition() const
HepLorentzVector & transform(const HepRotation &)
CLHEP::Hep3Vector G4ThreeVector
static constexpr double MeV
G4double GetPDGMass() const
G4double G4Log(G4double x)
~G4SingleDiffractiveExcitation()
G4double ChooseX(G4double Xmin, G4double Xmax) const
static G4Pow * GetInstance()
G4double powA(G4double A, G4double y) const
const G4ParticleDefinition const G4Material *G4double range
G4bool ExciteParticipants(G4VSplitableHadron *aPartner, G4VSplitableHadron *bPartner, G4bool ProjectileDiffraction) const
static constexpr double twopi
HepLorentzVector & rotateY(double)
G4GLOB_DLL std::ostream G4cerr
G4ThreeVector GaussianPt(G4double widthSquare, G4double maxPtSquare) const
HepLorentzVector & rotateZ(double)
G4int GetPDGEncoding() const
G4GLOB_DLL std::ostream G4cout
Hep3Vector boostVector() const
const G4LorentzVector & Get4Momentum() const
G4SingleDiffractiveExcitation()
CLHEP::HepLorentzVector G4LorentzVector
static constexpr double GeV
HepLorentzRotation & transform(const HepBoost &b)
void Set4Momentum(const G4LorentzVector &a4Momentum)