79{
80 #ifdef debugDoubleDiffraction
81 G4cout<<
G4endl<<
"G4QGSDiffractiveExcitation::ExciteParticipants - Double diffraction."<<
G4endl;
82 G4cout<<
"Proj Targ "<<projectile->GetDefinition()->GetParticleName()<<
" "<<target->GetDefinition()->GetParticleName()<<
G4endl;
83 G4cout<<
"Proj 4 Mom "<<projectile->Get4Momentum()<<
" "<<projectile->Get4Momentum().mag()<<
G4endl;
84 G4cout<<
"Targ 4 Mom "<<target->Get4Momentum() <<
" "<<target->Get4Momentum().mag() <<
G4endl;
85 #endif
86
88
89
91
93
94 if(M0projectile < projectile->GetDefinition()->GetPDGMass())
95 {
96 PutOnMassShell=1;
97 M0projectile=projectile->GetDefinition()->GetPDGMass();
98 }
99
100
102
104
105 if(M0target < target->GetDefinition()->GetPDGMass())
106 {
107 PutOnMassShell=1;
108 M0target=target->GetDefinition()->GetPDGMass();
109 }
110
114
115 if (SqrtS < M0projectile + M0target) {return false;}
116
117
118 G4double Mprojectile2 = M0projectile * M0projectile;
119 G4double Mtarget2 = M0target * M0target;
120
121
122
124
126
127 if ( Ptmp.
pz() <= 0. )
128 {
129
130
131 return false;
132 }
133
136
138
139 Pprojectile.transform(toCms);
140 Ptarget.transform(toCms);
141
142 G4double PZcms2=(
S*
S+Mprojectile2*Mprojectile2+Mtarget2*Mtarget2-
143 2*
S*Mprojectile2-2*
S*Mtarget2-2*Mprojectile2*Mtarget2)/4./
S;
144
145 if (PZcms2 < 0) {return false;}
146
148
149 if (PutOnMassShell)
150 {
151 if (Pprojectile.z() > 0.)
152 {
153 Pprojectile.setPz( PZcms);
154 Ptarget.setPz( -PZcms);
155 }
156 else
157 {
158 Pprojectile.setPz(-PZcms);
159 Ptarget.setPz( PZcms);
160 };
161
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));
164 }
165
167
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;
171 #endif
172
173 G4int PrPDGcode=projectile->GetDefinition()->GetPDGEncoding();
174 G4int absPrPDGcode=std::abs(PrPDGcode);
177
178 if (M0projectile <= projectile->GetDefinition()->GetPDGMass())
179 {
180 if( absPrPDGcode > 1000 )
181 {
182 if ( absPrPDGcode > 4000 && absPrPDGcode < 6000 )
183 {
184 MinPrDiffMass = projectile->GetDefinition()->GetPDGMass()/CLHEP::GeV + 0.25;
185 AveragePt2 = 0.3;
186 }
187 else
188 {
189 MinPrDiffMass = 1.16;
190 AveragePt2 = 0.3;
191 }
192 }
193 else if( absPrPDGcode == 211 || PrPDGcode == 111)
194 {
195 MinPrDiffMass = 1.0;
196 AveragePt2 = 0.3;
197 }
198 else if( absPrPDGcode == 321 || absPrPDGcode == 130 || absPrPDGcode == 310)
199 {
200 MinPrDiffMass = 1.1;
201 AveragePt2 = 0.3;
202 }
203 else if( absPrPDGcode > 400 && absPrPDGcode < 600)
204 {
205 MinPrDiffMass = projectile->GetDefinition()->GetPDGMass()/CLHEP::GeV + 0.25;
206 AveragePt2 = 0.3;
207 }
208 else
209 {
210 MinPrDiffMass = 1.16;
211 AveragePt2 = 0.3;
212 }
213 }
214 else
215 {
216 MinPrDiffMass = M0projectile + 220.0*MeV;
217 AveragePt2 = 0.3;
218 }
219
220 MinPrDiffMass = MinPrDiffMass * GeV;
221 AveragePt2 = AveragePt2 * GeV*GeV;
222
224
225 if (SqrtS < MinPrDiffMass + MinTrDiffMass) {return false;}
226
227 G4double MinPrDiffMass2 = MinPrDiffMass * MinPrDiffMass;
228 G4double MinTrDiffMass2 = MinTrDiffMass * MinTrDiffMass;
229
234
237
239 do {
240 if (whilecount++ >= 500 && (whilecount%100)==0)
241 if (whilecount > 1000 ) {
243 return false;
244 }
245
246
248
250 ProjMassT2=MinPrDiffMass2+Pt2;
251 ProjMassT =std::sqrt(ProjMassT2);
252
253 TargMassT2=MinTrDiffMass2+Pt2;
254 TargMassT =std::sqrt(TargMassT2);
255
256 if (SqrtS < ProjMassT + TargMassT) continue;
257
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);
262
263 G4double PMinusMin=std::sqrt(ProjMassT2+PZcms2)-PZcms;
265
266 PMinusNew=ChooseP(PMinusMin,PMinusMax);
267 Qminus=PMinusNew-Pprojectile.minus();
268
269 G4double TPlusMin=std::sqrt(TargMassT2+PZcms2)-PZcms;
271
272 TPlusNew=ChooseP(TPlusMin, TPlusMax);
273 Qplus=-(TPlusNew-Ptarget.plus());
274
275 Qmomentum.
setPz( (Qplus-Qminus)/2 );
276 Qmomentum.
setE( (Qplus+Qminus)/2 );
277
278 } while ( (Pprojectile+Qmomentum).mag2() < MinPrDiffMass2 ||
279 (Ptarget -Qmomentum).mag2() < MinTrDiffMass2 );
280
281 Pprojectile += Qmomentum;
282 Ptarget -= Qmomentum;
283
284
287
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;
291 #endif
292
293 target->Set4Momentum(Ptarget);
294 projectile->Set4Momentum(Pprojectile);
295
296 return true;
297}
G4double S(G4double temp)
CLHEP::HepLorentzVector G4LorentzVector
CLHEP::Hep3Vector G4ThreeVector
G4GLOB_DLL std::ostream G4cout
Hep3Vector boostVector() const
HepLorentzVector & rotateZ(double)
HepLorentzVector & rotateY(double)
HepLorentzVector & transform(const HepRotation &)