60{
61 #ifdef debugSingleDiffraction
63 #endif
64
66 G4double Mprojectile = projectile->GetDefinition()->GetPDGMass();
67 G4double Mprojectile2=
sqr(projectile->GetDefinition()->GetPDGMass());
68
70 G4double Mtarget = target->GetDefinition()->GetPDGMass();
71 G4double Mtarget2=
sqr(target->GetDefinition()->GetPDGMass());
72
73 #ifdef debugSingleDiffraction
74 G4cout<<
"Proj Targ "<<projectile->GetDefinition()->GetPDGEncoding()<<
" "<<target->GetDefinition()->GetPDGEncoding()<<
G4endl;
76 <<
" "<<Ptarget <<
" "<<Ptarget.
mag() <<
G4endl;
77 #endif
78
82
83 #ifdef debugSingleDiffraction
84 G4cout<<
"SqrtS-Mprojectile-Mtarget "<<SqrtS<<
" "<<Mprojectile<<
" "<<Mtarget
85 <<
" "<<SqrtS-Mprojectile-Mtarget<<
G4endl;
86 #endif
87 if (SqrtS-Mprojectile-Mtarget <= 250.0*MeV) {
88 #ifdef debugSingleDiffraction
89 G4cerr<<
"Projectile: "<<projectile->GetDefinition()->GetPDGEncoding()<<
" "
90 <<Pprojectile<<
" "<<Pprojectile.
mag()<<
G4endl;
91 G4cerr<<
"Target: "<<target->GetDefinition()->GetPDGEncoding()<<
" "
92 <<Ptarget<<
" "<<Ptarget.mag()<<
G4endl;
93 G4cerr<<
"sqrt(S) = "<<SqrtS<<
" Mp + Mt = "<<Pprojectile.
mag()+Ptarget.mag()<<
G4endl;
94 #endif
95 return true;
96 }
97
99
101
102 if ( Ptmp.
pz() <= 0. )
103 {
104
105
106 return false;
107 }
108
111
113
114 Pprojectile.transform(toCms);
115 Ptarget.transform(toCms);
116 #ifdef debugSingleDiffraction
117 G4cout <<
"Pprojectile in CMS " << Pprojectile <<
G4endl;
119 #endif
121
122 G4double ProjectileMinDiffrMass(0.), TargetMinDiffrMass(0.);
124 G4int absPDGcode=std::abs(projectile->GetDefinition()->GetPDGEncoding());
125
126 if ( ProjectileDiffraction ) {
127 if ( absPDGcode > 1000 )
128 {
129 if ( absPDGcode > 4000 && absPDGcode < 6000 )
130 {
131 ProjectileMinDiffrMass = projectile->GetDefinition()->GetPDGMass()/CLHEP::GeV + 0.25;
132 AveragePt2 = 0.3;
133 }
134 else
135 {
136 ProjectileMinDiffrMass = 1.16;
137 AveragePt2 = 0.3;
138 }
139 }
140 else if( absPDGcode == 211 || absPDGcode == 111)
141 {
142 ProjectileMinDiffrMass = 1.0;
143 AveragePt2 = 0.3;
144 }
145 else if( absPDGcode == 321 || absPDGcode == 130 || absPDGcode == 310)
146 {
147 ProjectileMinDiffrMass = 1.1;
148 AveragePt2 = 0.3;
149 }
150 else if( absPDGcode == 22)
151 {
152 ProjectileMinDiffrMass = 0.25;
153 AveragePt2 = 0.36;
154 }
155 else if( absPDGcode > 400 && absPDGcode < 600)
156 {
157 ProjectileMinDiffrMass = projectile->GetDefinition()->GetPDGMass()/CLHEP::GeV + 0.25;
158 AveragePt2 = 0.3;
159 }
160 else
161 {
162 ProjectileMinDiffrMass = 1.1;
163 AveragePt2 = 0.3;
164 };
165
166 ProjectileMinDiffrMass = ProjectileMinDiffrMass * GeV;
167 Mprojectile2=
sqr(ProjectileMinDiffrMass);
168 }
169 else
170 {
171 TargetMinDiffrMass = 1.16*GeV;
172 Mtarget2 =
sqr( TargetMinDiffrMass) ;
173 AveragePt2 = 0.3;
174 }
175
176 AveragePt2 = AveragePt2 * GeV*GeV;
177
183 G4double PMinusNew, PPlusNew, TPlusNew, TMinusNew;
184
187
189 do {
190 whilecount++;
191
192 if (whilecount > 1000 )
193 {
195 return false;
196 }
197
198
200
202
203 ProjMassT2 = Mprojectile2 + Pt2;
204 ProjMassT = std::sqrt( ProjMassT2 );
205 TargMassT2 = Mtarget2 + Pt2;
206 TargMassT = std::sqrt( TargMassT2 );
207
208 #ifdef debugSingleDiffraction
209 G4cout<<whilecount<<
" "<<Pt2<<
" "<<ProjMassT<<
" "<<TargMassT<<
" "<<SqrtS<<
" "<<
S<<
" "<<ProjectileDiffraction<<
G4endl;
210 #endif
211 if ( SqrtS < ProjMassT + TargMassT ) continue;
212
213 PZcms2 = (
S*
S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2
214 - 2.0*
S*ProjMassT2 - 2.0*
S*TargMassT2 - 2.0*ProjMassT2*TargMassT2 ) / 4.0 /
S;
215
216 if ( PZcms2 < 0 ) continue;
217
218 PZcms = std::sqrt( PZcms2 );
219
220 if ( ProjectileDiffraction )
221 {
222 PMinusMin = std::sqrt( ProjMassT2 + PZcms2 ) - PZcms;
223 PMinusMax = SqrtS - TargMassT;
224
225 PMinusNew = ChooseX( PMinusMin, PMinusMax );
226 TMinusNew = SqrtS - PMinusNew;
227
228 Qminus = Ptarget.minus() - TMinusNew;
229 TPlusNew = TargMassT2 / TMinusNew;
230 Qplus = Ptarget.plus() - TPlusNew;
231
232 } else {
233 TPlusMin = std::sqrt( TargMassT2 + PZcms2 ) - PZcms;
234 TPlusMax = SqrtS - ProjMassT;
235
236 TPlusNew = ChooseX( TPlusMin, TPlusMax );
237 PPlusNew = SqrtS - TPlusNew;
238
239 Qplus = PPlusNew - Pprojectile.plus();
240 PMinusNew = ProjMassT2 / PPlusNew;
241 Qminus = PMinusNew - Pprojectile.minus();
242 }
243
244 Qmomentum.
setPz( (Qplus - Qminus)/2 );
245 Qmomentum.
setE( (Qplus + Qminus)/2 );
246
247 #ifdef debugSingleDiffraction
248 G4cout<<ProjectileDiffraction<<
" "<<( Pprojectile + Qmomentum ).mag2()<<
" "<< Mprojectile2<<
G4endl;
249 G4cout<<!ProjectileDiffraction<<
" "<<( Ptarget - Qmomentum ).mag2()<<
" "<< Mtarget2<<
G4endl;
250 #endif
251
252 } while ( ( ProjectileDiffraction&&( Pprojectile + Qmomentum ).mag2() < Mprojectile2 ) ||
253 (!ProjectileDiffraction&&( Ptarget - Qmomentum ).mag2() < Mtarget2 ) );
254
255
256 Pprojectile += Qmomentum;
257
258 Ptarget -= Qmomentum;
259
260
263
264 #ifdef debugSingleDiffraction
265 G4cout <<
"Pprojectile in Lab. " << Pprojectile <<
G4endl;
267 G4cout <<
"G4SingleDiffractiveExcitation- Projectile mass " << Pprojectile.mag() <<
G4endl;
268 G4cout <<
"G4SingleDiffractiveExcitation- Target mass " << Ptarget.mag() <<
G4endl;
269 #endif
270
271 target->Set4Momentum(Ptarget);
272 projectile->Set4Momentum(Pprojectile);
273
274 return true;
275}
CLHEP::HepLorentzVector G4LorentzVector
CLHEP::Hep3Vector G4ThreeVector
G4GLOB_DLL std::ostream G4cerr
G4GLOB_DLL std::ostream G4cout
Hep3Vector boostVector() const
HepLorentzVector & rotateZ(double)
HepLorentzVector & rotateY(double)
HepLorentzVector & transform(const HepRotation &)