38 theStringDecay(aStringDecay)
45 throw G4HadronicException(__FILE__, __LINE__,
"G4ExcitedStringDecay::copy ctor not accessible");
54 throw G4HadronicException(__FILE__, __LINE__,
"G4ExcitedStringDecay::operator= meant to not be accessable");
71 if ( theStringDecay == NULL )
84 for (
unsigned int astring=0; astring < theStrings->size(); astring++)
86 KTsum+= theStrings->operator[](astring)->Get4Momentum();
92 G4bool NeedEnergyCorrector=
false;
101 NeedEnergyCorrector=
false;
103 for (
unsigned int astring=0; astring < theStrings->size(); astring++)
108 if ( theStrings->operator[](astring)->IsExcited() )
111 generatedKineticTracks=FragmentString(*theStrings->operator[](astring));
114 generatedKineticTracks->push_back(theStrings->operator[](astring)->GetKineticTrack());
117 if (generatedKineticTracks == NULL)
119 G4cerr <<
"G4VPartonStringModel:No KineticTracks produced" <<
G4endl;
124 for (
unsigned int aTrack=0; aTrack<generatedKineticTracks->size();aTrack++)
127 theResult->push_back(generatedKineticTracks->operator[](aTrack));
128 KTsum1+= (*generatedKineticTracks)[aTrack]->Get4Momentum();
130 KTsecondaries+=KTsum1;
134 if ( KTsum1.
e() > 0 && std::abs((KTsum1.
e()-theStrings->operator[](astring)->Get4Momentum().e()) / KTsum1.
e()) >
perMillion )
138 NeedEnergyCorrector=
true;
142 delete generatedKineticTracks;
148 if ( NeedEnergyCorrector ) success=EnergyAndMomentumCorrector(theResult, KTsum);
150 }
while(!success && (attempts < 10));
153 #ifdef debug_ExcitedStringDecay
155 for (
unsigned int aTrack=0; aTrack<theResult->size();aTrack++)
157 G4cout <<
" corrected tracks .. " << (*theResult)[aTrack]->GetDefinition()->GetParticleName()
158 <<
" " << (*theResult)[aTrack]->Get4Momentum() <<
G4endl;
159 KTsum1+= (*theResult)[aTrack]->Get4Momentum();
162 G4cout <<
"Needcorrector/success " << NeedEnergyCorrector <<
"/" << success <<
", Corrected total 4 momentum " << KTsum1 <<
G4endl;
163 if ( ! success )
G4cout <<
"failed to correct E/p" <<
G4endl;
169 G4bool G4ExcitedStringDecay::EnergyAndMomentumCorrector
172 const int nAttemptScale = 500;
173 const double ErrLimit = 1.E-5;
178 G4double TotalCollisionMass = TotalCollisionMom.
m();
182 unsigned int cHadron;
183 for(cHadron = 0; cHadron < Output->size(); cHadron++)
185 SumMom += Output->operator[](cHadron)->Get4Momentum();
186 SumMass += Output->operator[](cHadron)->GetDefinition()->GetPDGMass();
192 if (Output->size() < 2)
return FALSE;
194 if (SumMass > TotalCollisionMass)
return FALSE;
195 SumMass = SumMom.
m2();
196 if (SumMass < 0)
return FALSE;
197 SumMass = std::sqrt(SumMass);
209 for(cAttempt = 0; cAttempt < nAttemptScale; cAttempt++)
212 for(cHadron = 0; cHadron < Output->size(); cHadron++)
214 G4LorentzVector HadronMom = Output->operator[](cHadron)->Get4Momentum();
216 G4double E = std::sqrt(HadronMom.
vect().
mag2() +
sqr(Output->operator[](cHadron)->GetDefinition()->GetPDGMass()));
218 Output->operator[](cHadron)->Set4Momentum(HadronMom);
221 Scale = TotalCollisionMass/Sum;
222 #ifdef debug_G4ExcitedStringDecay
223 G4cout <<
"Scale-1=" << Scale -1
224 <<
", TotalCollisionMass=" << TotalCollisionMass
228 if (std::fabs(Scale - 1) <= ErrLimit)
234 #ifdef debug_G4ExcitedStringDecay
237 G4cout <<
"G4ExcitedStringDecay::EnergyAndMomentumCorrector - Warning"<<
G4endl;
238 G4cout <<
" Scale not unity at end of iteration loop: "<<TotalCollisionMass<<
" "<<Sum<<
" "<<Scale<<
G4endl;
239 G4cout <<
" Number of secondaries: " << Output->size() <<
G4endl;
240 G4cout <<
" Wanted total energy: " << TotalCollisionMom.
e() <<
G4endl;
241 G4cout <<
" Increase number of attempts or increase ERRLIMIT"<<
G4endl;
Hep3Vector boostVector() const
std::vector< G4ExcitedString * > G4ExcitedStringVector
virtual G4KineticTrackVector * FragmentString(const G4ExcitedString &theString)=0
virtual G4KineticTrackVector * FragmentStrings(const G4ExcitedStringVector *theStrings)
G4GLOB_DLL std::ostream G4cout
void Boost(G4ThreeVector &Velocity)
simulatedPeak Scale(1/simulationNormalisationFactor)
void setVect(const Hep3Vector &)
virtual ~G4ExcitedStringDecay()
G4GLOB_DLL std::ostream G4cerr