66 {
67 projectile->IncrementCollisionCount( 1 );
68 target->IncrementCollisionCount( 1 );
69
70 if ( projectile->Get4Momentum().z() < 0.0 ) return false;
71
72
75 G4double M0projectile2 = M0projectile * M0projectile;
76
77
80 G4double M0target2 = M0target * M0target;
81
83
84
86 Psum = Pprojectile + Ptarget;
89 if ( Ptmp.
pz() <= 0.0 )
return false;
90
91
95 Pprojectile.transform( toCms );
96 Ptarget.transform( toCms );
97
101 if ( SqrtS < M0projectile + M0target ) return false;
102
103 PZcms2 = (
S*
S +
sqr( M0projectile2 ) +
sqr( M0target2 )
104 - 2*
S*M0projectile2 - 2*
S*M0target2 - 2*M0projectile2*M0target2 ) / 4.0 /
S;
105
106 PZcms = ( PZcms2 > 0.0 ? std::sqrt( PZcms2 ) : 0.0 );
107
109
110
115
116 const G4int maxNumberOfLoops = 1000;
117 G4int loopCounter = 0;
118 do {
119 Qmomentum =
G4LorentzVector( GaussianPt( AveragePt2, maxPtSquare ), 0.0 );
121 ProjMassT2 = M0projectile2 + Pt2;
122 ProjMassT = std::sqrt( ProjMassT2 );
123 TargMassT2 = M0target2 + Pt2;
124 TargMassT = std::sqrt( TargMassT2 );
125 } while ( ( SqrtS < ProjMassT + TargMassT ) &&
126 ++loopCounter < maxNumberOfLoops );
127 if ( loopCounter >= maxNumberOfLoops ) {
128 return false;
129 }
130
131 PZcms2 = (
S*
S +
sqr( ProjMassT2 ) +
sqr( TargMassT2 )
132 - 2.0*
S*ProjMassT2 - 2.0*
S*TargMassT2 - 2.0*ProjMassT2*TargMassT2 ) / 4.0 /
S;
133
134 if ( PZcms2 < 0.0 ) { PZcms2 = 0.0; };
135 PZcms = std::sqrt( PZcms2 );
136 Pprojectile.setPz( PZcms );
137 Ptarget.setPz( -PZcms );
138 Pprojectile += Qmomentum;
139 Ptarget -= Qmomentum;
140
141
144
145
146 projectile->SetTimeOfCreation( target->GetTimeOfCreation() );
147 projectile->SetPosition( target->GetPosition() );
148
149
150
151
152 projectile->Set4Momentum( Pprojectile );
153 target->Set4Momentum( Ptarget );
154
155
156
157
158 return true;
159}
G4double S(G4double temp)
CLHEP::HepLorentzVector G4LorentzVector
CLHEP::Hep3Vector G4ThreeVector
Hep3Vector boostVector() const
HepLorentzVector & rotateZ(double)
HepLorentzVector & rotateY(double)
HepLorentzVector & transform(const HepRotation &)
G4double GetAvaragePt2ofElasticScattering()