5 #include "FullModelHadronicProcess.hh"
6 #include "G4Version.hh"
7 #include "G4ProcessManager.hh"
8 #include "G4ProcessHelper.hh"
9 #include "G4ParticleTable.hh"
10 #include "FullModelReactionDynamics.hh"
12 #if G4VERSION_NUMBER < 1100
13 #include "G4HadReentrentException.hh"
15 #include "G4HadronicException.hh"
21 FullModelHadronicProcess::FullModelHadronicProcess(
const G4String& processName) :
22 G4VDiscreteProcess(processName,fUserDefined),m_theParticle(0),m_newParticle(0),m_cache(0)
25 m_theHelper = G4ProcessHelper::Instance();
30 G4bool FullModelHadronicProcess::IsApplicable(
const G4ParticleDefinition& aP)
33 return m_theHelper->ApplicabilityTester(aP);
37 G4double FullModelHadronicProcess::GetMicroscopicCrossSection(
const G4DynamicParticle *aParticle,
38 const G4Element *anElement,
42 G4double InclXsec = m_theHelper->GetInclusiveCrossSection(aParticle,anElement);
46 G4double HighestEnergyLimit = 10 *
CLHEP::TeV ;
47 G4double LowestEnergyLimit = 1 *
CLHEP::eV;
49 G4double ParticleEnergy = aParticle->GetKineticEnergy();
52 if (ParticleEnergy > HighestEnergyLimit || ParticleEnergy < LowestEnergyLimit){
61 G4double FullModelHadronicProcess::GetMeanFreePath(
const G4Track& aTrack, G4double, G4ForceCondition*)
63 G4Material *aMaterial = aTrack.GetMaterial();
64 const G4DynamicParticle *aParticle = aTrack.GetDynamicParticle();
67 G4int nElements = aMaterial->GetNumberOfElements();
69 const G4double *theAtomicNumDensityVector =
70 aMaterial->GetAtomicNumDensityVector();
71 G4double aTemp = aMaterial->GetTemperature();
73 for( G4int
i=0;
i<nElements; ++
i )
76 GetMicroscopicCrossSection( aParticle, (*aMaterial->GetElementVector())[
i], aTemp);
77 sigma += theAtomicNumDensityVector[
i] * xSection;
85 G4VParticleChange* FullModelHadronicProcess::PostStepDoIt(
const G4Track& aTrack,
89 aParticleChange.Initialize(aTrack);
91 const G4DynamicParticle* IncidentRhadron = aTrack.GetDynamicParticle();
93 const G4ThreeVector& aPosition = aTrack.GetPosition();
94 const G4int theIncidentPDG = IncidentRhadron->GetDefinition()->GetPDGEncoding();
95 G4ParticleTable* theParticleTable = G4ParticleTable::GetParticleTable();
96 std::vector<G4ParticleDefinition*> theParticleDefinitions;
97 G4bool IncidentSurvives =
false;
98 G4bool TargetSurvives =
false;
99 G4Nucleus targetNucleus(aTrack.GetMaterial());
100 G4ParticleDefinition* outgoingRhadron=0;
101 G4ParticleDefinition* outgoingCloud=0;
102 G4ParticleDefinition* outgoingTarget=0;
105 G4ThreeVector p_0 = IncidentRhadron->GetMomentum();
109 G4double e_kin_0 = IncidentRhadron->GetKineticEnergy();
111 G4DynamicParticle* cloudParticle =
new G4DynamicParticle();
113 cloudParticle->SetDefinition(CustomIncident->
GetCloud());
115 double scale=cloudParticle->GetDefinition()->GetPDGMass()/IncidentRhadron->GetDefinition()->GetPDGMass();
117 G4LorentzVector cloudMomentum;
118 cloudMomentum.setVectM(IncidentRhadron->GetMomentum()*
scale,cloudParticle->GetDefinition()->GetPDGMass());
119 G4LorentzVector gluinoMomentum;
121 gluinoMomentum.setVectM(IncidentRhadron->GetMomentum()*(1.-
scale),CustomIncident->
GetSpectator()->GetPDGMass());
123 G4LorentzVector FullRhadron4Momentum = IncidentRhadron->Get4Momentum();
124 G4LorentzVector Cloud4Momentum = cloudMomentum;
126 cloudParticle->Set4Momentum(cloudMomentum);
128 G4DynamicParticle* OrgPart = cloudParticle;
130 G4double ek = OrgPart->GetKineticEnergy();
131 G4double amas = OrgPart->GetDefinition()->GetPDGMass();
132 G4double tkin = targetNucleus.Cinema( ek );
134 OrgPart->SetKineticEnergy( ek );
135 G4double
et = ek + amas;
136 G4double
p = std::sqrt( std::abs((
et-amas)*(
et+amas)) );
137 G4double pp = OrgPart->GetMomentum().mag();
140 G4ThreeVector
momentum = OrgPart->GetMomentum();
141 OrgPart->SetMomentum(
momentum * (
p/pp) );
145 if(ek>0.) {tkin = targetNucleus.EvaporationEffects( ek ); ek -= tkin;}
146 if(ek+gluinoMomentum.e()-gluinoMomentum.m()<=0.1*
CLHEP::MeV||ek<=0.) {
147 G4cout<<
"Stopping particle...:"<<G4endl;
148 G4cout<<
"ek: "<<ek/
CLHEP::MeV<<
" MeV"<<G4endl;
149 aParticleChange.ProposeTrackStatus( fStopButAlive );
151 return &aParticleChange;
153 OrgPart->SetKineticEnergy( ek );
155 p = std::sqrt( std::abs((
et-amas)*(
et+amas)) );
156 pp = OrgPart->GetMomentum().mag();
160 G4ThreeVector
momentum = OrgPart->GetMomentum();
161 OrgPart->SetMomentum(
momentum * (
p/pp) );
168 G4ParticleDefinition* aTarget;
169 ReactionProduct
rp = m_theHelper->GetFinalState(aTrack,aTarget);
180 G4int NumberOfSecondaries =
rp.size();
184 G4LorentzVector Target4Momentum;
185 Target4Momentum.setVectM(CLHEP::Hep3Vector(0.,0.,0.),aTarget->GetPDGMass());
187 G4LorentzVector psum_full,psum_cloud;
188 psum_full = FullRhadron4Momentum + Target4Momentum;
189 psum_cloud = Cloud4Momentum + Target4Momentum;
190 G4ThreeVector trafo_full_cms = (-1)*psum_full.boostVector();
191 G4ThreeVector trafo_cloud_cms = (-1)*psum_cloud.boostVector();
200 G4ParticleDefinition* tempDef = theParticleTable->FindParticle(*
it);
202 if (tempDef==aTarget) TargetSurvives =
true;
206 outgoingRhadron = tempDef;
208 if(outgoingCloud == 0)
210 G4cout <<
"FullModelHadronicProcess::PostStepDoIt() Definition of outgoing particle cloud for " << tempDef->GetParticleName() <<
" not available!!" << G4endl;
215 if (tempCust==0&&
rp.size()==2) outgoingTarget = tempDef;
216 if (tempDef->GetPDGEncoding()==theIncidentPDG) {
217 IncidentSurvives =
true;
219 theParticleDefinitions.push_back(tempDef);
223 if (NULL==outgoingRhadron) {
228 description <<
"PostStepDoIt: outgoingRHadron is null";
229 G4Exception(
"FullModelHadronicProcess",
"NoOutGoingRHadron", FatalException,
description);
232 if (outgoingTarget==0) outgoingTarget = theParticleTable->FindParticle(
rp[1]);
236 if(IncidentSurvives) NumberOfSecondaries--;
237 aParticleChange.SetNumberOfSecondaries(NumberOfSecondaries);
244 const G4HadProjectile* originalIncident =
new G4HadProjectile(*OrgPart);
248 G4HadProjectile* org2 =
new G4HadProjectile(*OrgPart);
249 G4LorentzRotation trans = org2->GetTrafoToLab();
252 G4DynamicParticle *originalTarget =
new G4DynamicParticle;
253 originalTarget->SetDefinition(aTarget);
255 G4ReactionProduct targetParticle(aTarget);
258 G4ReactionProduct currentParticle(
const_cast<G4ParticleDefinition *
>(originalIncident->GetDefinition() ) );
259 currentParticle.SetMomentum( originalIncident->Get4Momentum().vect() );
260 currentParticle.SetKineticEnergy( originalIncident->GetKineticEnergy() );
263 G4ReactionProduct modifiedOriginal = currentParticle;
265 currentParticle.SetSide( 1 );
266 targetParticle.SetSide( -1 );
267 G4bool incidentHasChanged =
false;
268 if (!IncidentSurvives) incidentHasChanged =
true;
269 G4bool targetHasChanged =
false;
270 if (!TargetSurvives) targetHasChanged =
true;
271 G4bool quasiElastic =
false;
272 if (
rp.size()==2) quasiElastic =
true;
273 G4FastVector<G4ReactionProduct,GHADLISTSIZE>
vec;
282 for (G4int
i = 0;
i!=NumberOfSecondaries;
i++){
283 if(theParticleDefinitions[
i]!=aTarget
284 && theParticleDefinitions[
i]!=originalIncident->GetDefinition()
285 && theParticleDefinitions[
i]!=outgoingRhadron
286 && theParticleDefinitions[
i]!=outgoingTarget)
288 G4ReactionProduct* pa =
new G4ReactionProduct;
289 pa->SetDefinition( theParticleDefinitions[
i] );
290 (G4UniformRand() < 0.5) ? pa->SetSide( -1 ) : pa->SetSide( 1 );
291 vec.SetElement( vecLen++, pa );
297 if (incidentHasChanged) currentParticle.SetDefinitionAndUpdateE( outgoingCloud );
298 if (incidentHasChanged) modifiedOriginal.SetDefinition( outgoingCloud );
299 if (targetHasChanged) targetParticle.SetDefinitionAndUpdateE( outgoingTarget );
302 CalculateMomenta(
vec, vecLen,
303 originalIncident, originalTarget, modifiedOriginal,
304 targetNucleus, currentParticle, targetParticle,
305 incidentHasChanged, targetHasChanged, quasiElastic );
308 G4String cPname = currentParticle.GetDefinition()->GetParticleName();
311 aParticleChange.SetNumberOfSecondaries(vecLen+NumberOfSecondaries);
313 G4LorentzVector cloud_p4_new;
315 cloud_p4_new.setVectM(currentParticle.GetMomentum(),currentParticle.GetMass());
316 cloud_p4_new *= trans;
318 G4LorentzVector cloud_p4_old_full = Cloud4Momentum;
319 cloud_p4_old_full.boost(trafo_full_cms);
320 G4LorentzVector cloud_p4_old_cloud = Cloud4Momentum;
321 cloud_p4_old_cloud.boost(trafo_cloud_cms);
322 G4LorentzVector cloud_p4_full = cloud_p4_new;
323 cloud_p4_full.boost(trafo_full_cms);
324 G4LorentzVector cloud_p4_cloud = cloud_p4_new;
325 cloud_p4_cloud.boost(trafo_cloud_cms);
327 G4LorentzVector p_g_cms = gluinoMomentum;
328 p_g_cms.boost(trafo_full_cms);
330 G4LorentzVector p4_new;
331 p4_new.setVectM( cloud_p4_new.v() + gluinoMomentum.v(), outgoingRhadron->GetPDGMass() );
334 p_new = p4_new.vect();
336 aParticleChange.ProposeLocalEnergyDeposit((p4_new-cloud_p4_new-gluinoMomentum).
m());
338 if( incidentHasChanged )
340 G4DynamicParticle*
p0 =
new G4DynamicParticle;
341 p0->SetDefinition( outgoingRhadron );
342 p0->SetMomentum( p_new );
345 G4Track* Track0 =
new G4Track(
p0,
346 aTrack.GetGlobalTime(),
348 aParticleChange.AddSecondary(Track0);
353 if(
p0->GetKineticEnergy()>e_kin_0+1.e-9) {
354 G4cout<<
"ALAAAAARM!!! (incident changed from "
355 <<IncidentRhadron->GetDefinition()->GetParticleName()
356 <<
" to "<<
p0->GetDefinition()->GetParticleName()<<
")"<<G4endl;
357 G4cout<<
"Energy loss: "<<(e_kin_0-
p0->GetKineticEnergy())/
CLHEP::GeV<<
" GeV (should be positive)"<<G4endl;
359 if(
rp.size()!=3) G4cout<<
"DOUBLE ALAAAAARM!!!"<<G4endl;
363 if(std::abs(
p0->GetKineticEnergy()-e_kin_0)>100*
CLHEP::GeV) {
364 G4cout<<
"Diff. too big"<<G4endl;
367 aParticleChange.ProposeTrackStatus( fStopAndKill );
372 G4double
p = p_new.mag();
374 aParticleChange.ProposeMomentumDirection( p_new.x()/
p, p_new.y()/
p, p_new.z()/
p );
376 aParticleChange.ProposeMomentumDirection( 1.0, 0.0, 0.0 );
378 G4double aE = sqrt(
p*
p+(outgoingRhadron->GetPDGMass()*outgoingRhadron->GetPDGMass()) );
379 e_kin = aE - outgoingRhadron->GetPDGMass();
382 G4cout<<
"ALAAAAARM!!!"<<G4endl;
383 G4cout<<
"Energy loss: "<<(e_kin_0-e_kin)/
CLHEP::GeV<<
" GeV (should be positive)"<<G4endl;
384 if(
rp.size()!=3) G4cout<<
"DOUBLE ALAAAAARM!!!"<<G4endl;
387 G4cout<<
"Diff. too big"<<G4endl;
391 aParticleChange.ProposeEnergy( aE- outgoingRhadron->GetPDGMass() );
393 G4cout<<
"Diff. too big"<<G4endl;
398 if( targetParticle.GetMass() > 0.0 )
400 G4DynamicParticle *
p1 =
new G4DynamicParticle;
401 p1->SetDefinition( targetParticle.GetDefinition() );
403 G4ThreeVector
momentum = targetParticle.GetMomentum();
406 p1->SetMomentum( (trans*
p1->Get4Momentum()).vect());
407 G4Track* Track1 =
new G4Track(
p1,
408 aTrack.GetGlobalTime(),
410 aParticleChange.AddSecondary(Track1);
412 G4DynamicParticle *pa;
415 for(
int i=0;
i<vecLen; ++
i )
417 pa =
new G4DynamicParticle();
418 pa->SetDefinition(
vec[
i]->GetDefinition() );
419 pa->SetMomentum(
vec[
i]->GetMomentum() );
420 pa->Set4Momentum(trans*(pa->Get4Momentum()));
421 G4ThreeVector
pvec = pa->GetMomentum();
422 G4Track* Trackn =
new G4Track(pa,
423 aTrack.GetGlobalTime(),
425 aParticleChange.AddSecondary(Trackn);
429 delete originalIncident;
430 delete originalTarget;
433 ClearNumberOfInteractionLengthLeft();
436 return &aParticleChange;
441 void FullModelHadronicProcess::CalculateMomenta(
442 G4FastVector<G4ReactionProduct,GHADLISTSIZE> &
vec,
444 const G4HadProjectile *originalIncident,
445 const G4DynamicParticle *originalTarget,
446 G4ReactionProduct &modifiedOriginal,
447 G4Nucleus &targetNucleus,
448 G4ReactionProduct ¤tParticle,
449 G4ReactionProduct &targetParticle,
450 G4bool &incidentHasChanged,
451 G4bool &targetHasChanged,
452 G4bool quasiElastic )
454 FullModelReactionDynamics theReactionDynamics;
457 m_what = originalIncident->Get4Momentum().vect();
470 theReactionDynamics.TwoBody(
vec, vecLen,
471 modifiedOriginal, originalTarget,
472 currentParticle, targetParticle,
473 targetNucleus, targetHasChanged );
479 G4ReactionProduct leadingStrangeParticle;
480 G4bool leadFlag = MarkLeadingStrangeParticle( currentParticle,
482 leadingStrangeParticle );
489 G4bool finishedGenXPt =
false;
490 G4bool annihilation =
false;
491 if( originalIncident->GetDefinition()->GetPDGEncoding() < 0 &&
492 currentParticle.GetMass() == 0.0 && targetParticle.GetMass() == 0.0 )
496 G4double ekcor = 1.0;
497 G4double ek = originalIncident->GetKineticEnergy();
500 const G4double tarmas = originalTarget->GetDefinition()->GetPDGMass();
502 const G4double atomicWeight = targetNucleus.AtomicMass(targetNucleus.GetA_asInt(),targetNucleus.GetZ_asInt());
503 ek = 2*tarmas + ek*(1.+ekcor/atomicWeight);
505 G4double tkin = targetNucleus.Cinema( ek );
508 modifiedOriginal.SetKineticEnergy( ekOrg );
511 const G4double twsup[] = { 1.0, 0.7, 0.5, 0.3, 0.2, 0.1 };
512 G4double rand1 = G4UniformRand();
513 G4double rand2 = G4UniformRand();
514 if (annihilation || vecLen >= 6 ||
515 (modifiedOriginal.GetKineticEnergy()/
CLHEP::GeV >= 1.0 &&
517 (originalIncident->GetDefinition() == G4KaonPlus::KaonPlus() ||
518 originalIncident->GetDefinition() == G4KaonMinus::KaonMinus() ||
519 originalIncident->GetDefinition() == G4KaonZeroLong::KaonZeroLong() ||
520 originalIncident->GetDefinition() == G4KaonZeroShort::KaonZeroShort()))
521 || rand2 > twsup[vecLen])))
523 theReactionDynamics.GenerateXandPt(
vec, vecLen,
524 modifiedOriginal, originalIncident,
525 currentParticle, targetParticle,
526 targetNucleus, incidentHasChanged,
527 targetHasChanged, leadFlag,
528 leadingStrangeParticle);
535 G4bool finishedTwoClu =
false;
536 if( modifiedOriginal.GetTotalMomentum()/
CLHEP::MeV < 1.0 )
539 for(G4int
i=0;
i<vecLen;
i++)
delete vec[
i];
545 theReactionDynamics.SuppressChargedPions(
vec, vecLen,
546 modifiedOriginal, currentParticle,
548 incidentHasChanged );
552 finishedTwoClu = theReactionDynamics.TwoCluster(
vec, vecLen,
553 modifiedOriginal, originalIncident,
554 currentParticle, targetParticle,
555 targetNucleus, incidentHasChanged,
556 targetHasChanged, leadFlag,
557 leadingStrangeParticle );
559 #if G4VERSION_NUMBER < 1100
560 catch(G4HadReentrentException& aC)
563 throw G4HadReentrentException(__FILE__, __LINE__,
"Failing to calculate momenta");
566 catch(G4HadronicException& aC)
569 throw G4HadronicException(__FILE__, __LINE__,
"Failing to calculate momenta");
599 theReactionDynamics.TwoBody(
vec, vecLen,
600 modifiedOriginal, originalTarget,
601 currentParticle, targetParticle,
602 targetNucleus, targetHasChanged );
606 G4bool FullModelHadronicProcess::MarkLeadingStrangeParticle(
607 const G4ReactionProduct ¤tParticle,
608 const G4ReactionProduct &targetParticle,
609 G4ReactionProduct &leadParticle )
617 if( (currentParticle.GetMass() >= G4KaonPlus::KaonPlus()->GetPDGMass()) &&
622 leadParticle = currentParticle;
624 else if( (targetParticle.GetMass() >= G4KaonPlus::KaonPlus()->GetPDGMass()) &&
629 leadParticle = targetParticle;
634 void FullModelHadronicProcess::Rotate(G4FastVector<G4ReactionProduct,GHADLISTSIZE> &
vec, G4int &vecLen)
639 for(
i=0;
i<vecLen; ++
i )
647 const G4DynamicParticle* FullModelHadronicProcess::FindRhadron(G4ParticleChange* aParticleChange)
649 G4int nsec = aParticleChange->GetNumberOfSecondaries();
650 if (nsec==0)
return 0;
652 G4bool
found =
false;
653 while (
i!=nsec && !
found){
656 if (
dynamic_cast<CustomParticle*
>(aParticleChange->GetSecondary(
i)->GetDynamicParticle()->GetDefinition())!=0)
found =
true;
660 if(
found)
return aParticleChange->GetSecondary(
i)->GetDynamicParticle();