87 G4ThreeVector the3IncidentPart = fCache.
Get().theProjectileRP->GetMomentum();
88 G4double nEnergy = fCache.
Get().theProjectileRP->GetTotalEnergy();
90 G4double tEnergy = fCache.
Get().theTarget->GetTotalEnergy();
94 G4double cmsMom = std::sqrt(the3CMS*the3CMS);
95 G4double sqrts = std::sqrt((totE-cmsMom)*(totE+cmsMom));
99 aIncidentPart.
Lorentz(*fCache.
Get().theProjectileRP, theCMS);
105 anEnergy = fCache.
Get().theProjectileRP->GetKineticEnergy();
112 fCache.
Get().theTotalMeanEnergy=0;
115 std::vector<int> nParticles;
128 for(i=0; i<nProducts; i++)
130 G4int massCode =
G4int(theProducts[i].GetMassCode());
133 sumZ += massCode/1000 * nPart;
134 sumA += massCode % 1000 * nPart;
136 if( getenv(
"G4ParticleHPDebug") )
G4cout << i <<
" G4ParticleHPEnAngCorrelation::MULTIPLICITY " << massCode <<
" sumZ " << sumZ <<
" sumA " << sumA <<
" NPART " << nPart <<
G4endl;
138 nParticles.push_back( nPart );
141 double targetZ = fCache.
Get().theTarget->GetDefinition()->GetAtomicNumber();
142 double targetA = fCache.
Get().theTarget->GetDefinition()->GetAtomicMass();
143 targetZ += fCache.
Get().theProjectileRP->GetDefinition()->GetAtomicNumber();
144 targetA += fCache.
Get().theProjectileRP->GetDefinition()->GetAtomicMass();
145 if ( bAdjustFinalState ) {
161 if ( ( sumZ != targetZ || sumA != targetA )
162 && ( sumZ > targetZ || sumA > targetA || (targetZ-sumZ) >= (targetA-sumA) ) ) {
167 if ( getenv(
"G4ParticleHPDebug") )
168 G4cerr <<
" WRONG MULTIPLICITY Z= " << sumZ
171 <<
" > " << targetA <<
G4endl;
181 G4Exception(
"G4ParticleHPEnAngCorrelation::Sample",
184 "Too many trials were done. Exiting current loop by force. You may have Probably, the result violating (baryon number) conservation law will be obtained.");
192 for(i=0; i<nProducts; i++)
195 it = theProducts[i].
Sample(anEnergy,nParticles[i]);
203 fCache.
Get().theTotalMeanEnergy += aMeanEnergy;
207 fCache.
Get().theTotalMeanEnergy = anEnergy/nProducts+theProducts[i].
GetQValue();
211 for(
unsigned int ii=0; ii<it->size(); ii++)
215 it->operator[](ii)->GetTotalEnergy());
217 if( getenv(
"G4PHPTEST") )
G4cout <<
" G4particleHPEnAngCorrelation COS THETA " << std::cos(it->operator[](ii)->GetMomentum().theta()) << G4endl;
218 it->operator[](ii)->SetMomentum(pTmp1.vect());
219 it->operator[](ii)->SetTotalEnergy(pTmp1.e());
220 if( getenv(
"G4PHPTEST") )
G4cout <<
" G4particleHPEnAngCorrelation COS THETA after toLab " << std::cos(it->operator[](ii)->GetMomentum().theta()) << G4endl;
224 it->operator[](ii)->Lorentz(*(it->operator[](ii)), -1.*(*fCache.
Get().theTarget));
226 else if(frameFlag==2 )
229 if( getenv(
"G4ParticleHPDebug") )
230 G4cout <<
"G4ParticleHPEnAngCorrelation: before Lorentz boost "<<
231 it->at(ii)->GetKineticEnergy()<<
" "<<
232 it->at(ii)->GetMomentum()<<
G4endl;
234 it->operator[](ii)->Lorentz(*(it->operator[](ii)), -1.*theCMS);
236 if( getenv(
"G4ParticleHPDebug") )
237 G4cout <<
"G4ParticleHPEnAngCorrelation: after Lorentz boost "<<
238 it->at(ii)->GetKineticEnergy()<<
" "<<
239 it->at(ii)->GetMomentum()<<
G4endl;
243 else if(frameFlag==3)
245 if ( theProducts[i].GetMassCode() > 4 )
248 it->operator[](ii)->Lorentz(*(it->operator[](ii)), -1.*(*fCache.
Get().theTarget));
250 if( getenv(
"G4ParticleHPDebug") )
251 G4cout <<
"G4ParticleHPEnAngCorrelation: after Lorentz boost "<<
252 it->at(ii)->GetKineticEnergy()<<
" "<<
253 it->at(ii)->GetMomentum()<<
G4endl;
259 it->operator[](ii)->Lorentz(*(it->operator[](ii)), -1.*theCMS);
261 if( getenv(
"G4ParticleHPDebug") )
262 G4cout <<
"G4ParticleHPEnAngCorrelation: after Lorentz boost "<<
263 it->at(ii)->GetKineticEnergy()<<
" "<<
264 it->at(ii)->GetMomentum()<<
G4endl;
270 throw G4HadronicException(__FILE__, __LINE__,
"G4ParticleHPEnAngCorrelation::Sample: The frame of the finalstate is not specified");
272 if( getenv(
"G4PHPTEST") )
G4cout << frameFlag <<
" G4particleHPEnAngCorrelation COS THETA after Lorentz " << std::cos(it->operator[](ii)->GetMomentum().theta()) << G4endl;
276 result->push_back(it->operator[](ii));
G4double G4ParticleHPJENDLHEData::G4double result
void Lorentz(const G4ReactionProduct &p1, const G4ReactionProduct &p2)
void SetMomentum(const G4double x, const G4double y, const G4double z)
G4double MeanEnergyOfThisInteraction()
HepLorentzRotation & rotateY(double delta)
std::vector< G4ReactionProduct * > G4ReactionProductVector
void SetMass(const G4double mas)
G4GLOB_DLL std::ostream G4cout
void SetTotalEnergy(const G4double en)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
HepLorentzRotation & rotateZ(double delta)
G4double GetTotalEnergy() const
G4int GetMultiplicity(G4double anEnergy)
G4ReactionProductVector * Sample(G4double anEnergy, G4int nParticles)
G4ThreeVector GetMomentum() const
HepLorentzRotation inverse() const
G4GLOB_DLL std::ostream G4cerr