40 theStringDecay(aStringDecay)
47 throw G4HadronicException(__FILE__, __LINE__,
"G4ExcitedStringDecay::copy ctor not accessible");
56 throw G4HadronicException(__FILE__, __LINE__,
"G4ExcitedStringDecay::operator= meant to not be accessable");
80 #ifdef debug_G4ExcitedStringDecay
82 G4cout<<
"--------------------------- G4ExcitedStringDecay ----------------------"<<
G4endl;
83 G4cout<<
"Hadronization of Excited Strings: theStrings->size() "<<theStrings->size()<<
G4endl;
86 for (
unsigned int astring=0; astring < theStrings->size(); astring++)
88 if ( theStrings->operator[](astring)->IsExcited() )
89 {KTsum+= theStrings->operator[](astring)->Get4Momentum();}
90 else {KTsum+=theStrings->operator[](astring)->GetKineticTrack()->Get4Momentum();}
97 for (
unsigned int astring=0; astring < theStrings->size(); astring++)
99 if ( theStrings->operator[](astring)->IsExcited() )
101 Ptmp=toCms * theStrings->operator[](astring)->GetLeftParton()->Get4Momentum();
102 theStrings->operator[](astring)->GetLeftParton()->Set4Momentum(Ptmp);
104 Ptmp=toCms * theStrings->operator[](astring)->GetRightParton()->Get4Momentum();
105 theStrings->operator[](astring)->GetRightParton()->Set4Momentum(Ptmp);
107 KTsum+= theStrings->operator[](astring)->Get4Momentum();
111 Ptmp=toCms * theStrings->operator[](astring)->GetKineticTrack()->Get4Momentum();
112 theStrings->operator[](astring)->GetKineticTrack()->Set4Momentum(Ptmp);
113 KTsum+= theStrings->operator[](astring)->GetKineticTrack()->Get4Momentum();
120 G4bool NeedEnergyCorrector=
false;
122 #ifdef debug_G4ExcitedStringDecay
123 G4cout<<
"New try No "<<attempts<<
" to hadronize strings"<<
G4endl;
132 NeedEnergyCorrector=
false;
134 for (
unsigned int astring=0; astring < theStrings->size(); astring++)
137 #ifdef debug_G4ExcitedStringDecay
138 G4cout<<
"String No "<<astring+1<<
" Excited? "<<theStrings->operator[](astring)->IsExcited()<<
G4endl;
140 G4cout<<
"String No "<<astring+1<<
" 4Momentum "<<theStrings->operator[](astring)->Get4Momentum()
141 <<
" "<<theStrings->operator[](astring)->Get4Momentum().mag()<<
G4endl;
145 if ( theStrings->operator[](astring)->IsExcited() )
147 #ifdef debug_G4ExcitedStringDecay
148 G4cout<<
"Fragment String with partons: "
149 <<theStrings->operator[](astring)->GetLeftParton()->GetPDGcode() <<
" "
150 <<theStrings->operator[](astring)->GetRightParton()->GetPDGcode()<<
" "
151 <<
"Direction "<<theStrings->operator[](astring)->GetDirection()<<
G4endl;
153 generatedKineticTracks=
FragmentString(*theStrings->operator[](astring));
154 #ifdef debug_G4ExcitedStringDecay
155 G4cout<<
"(G4ExcitedStringDecay) Number of produced hadrons = "<<generatedKineticTracks->size()<<
G4endl;
158 #ifdef debug_G4ExcitedStringDecay
161 G4LorentzVector Mom=theStrings->operator[](astring)->GetKineticTrack()->Get4Momentum();
163 theStrings->operator[](astring)->GetKineticTrack()->GetDefinition(),
164 theStrings->operator[](astring)->GetKineticTrack()->GetFormationTime(),
167 aTrack->
SetPosition(theStrings->operator[](astring)->GetKineticTrack()->GetPosition());
169 #ifdef debug_G4ExcitedStringDecay
174 generatedKineticTracks->push_back(aTrack);
178 if (generatedKineticTracks->size() == 0)
182 success=
false; NeedEnergyCorrector=
false;
break;
186 for (
unsigned int aTrack=0; aTrack<generatedKineticTracks->size();aTrack++)
188 #ifdef debug_G4ExcitedStringDecay
189 G4cout<<
"Prod part No. "<<aTrack+1<<
" "
190 <<(*generatedKineticTracks)[aTrack]->GetDefinition()->GetParticleName()<<
" "
191 <<(*generatedKineticTracks)[aTrack]->Get4Momentum()<<
G4endl;
193 theResult->push_back(generatedKineticTracks->operator[](aTrack));
194 KTsum1+= (*generatedKineticTracks)[aTrack]->Get4Momentum();
196 KTsecondaries+=KTsum1;
198 #ifdef debug_G4ExcitedStringDecay
199 G4cout <<
"String secondaries(" <<generatedKineticTracks->size()<<
")"<<G4endl
200 <<
"Init string momentum: "<< theStrings->operator[](astring)->Get4Momentum()<<G4endl
201 <<
"Final hadrons momentum: "<< KTsum1 <<
G4endl;
204 if ( KTsum1.e() > 0 && std::abs((KTsum1.e()-theStrings->operator[](astring)->Get4Momentum().e()) / KTsum1.e()) >
perMillion )
206 NeedEnergyCorrector=
true;
209 #ifdef debug_G4ExcitedStringDecay
210 G4cout<<
"NeedEnergyCorrection yes/no "<<NeedEnergyCorrector<<
G4endl;
214 delete generatedKineticTracks;
221 }
while(!success && (attempts < 10));
223 for (
unsigned int aTrack=0; aTrack<theResult->size();aTrack++)
225 Ptmp=(*theResult)[aTrack]->Get4Momentum();
226 Ptmp.transform( toLab);
227 (*theResult)[aTrack]->Set4Momentum(Ptmp);
230 #ifdef debug_G4ExcitedStringDecay
231 G4cout<<
"End of the strings fragmentation (G4ExcitedStringDecay)"<<
G4endl;
235 for (
unsigned int aTrack=0; aTrack<theResult->size();aTrack++)
237 G4cout <<
" corrected tracks .. " << (*theResult)[aTrack]->GetDefinition()->GetParticleName()
238 <<
" " << (*theResult)[aTrack]->Get4Momentum()
239 <<
" " << (*theResult)[aTrack]->Get4Momentum().mag()<<
G4endl;
240 KTsum1+= (*theResult)[aTrack]->Get4Momentum();
243 G4cout <<
"Needcorrector/success " << NeedEnergyCorrector <<
"/" << success <<
", Corrected total 4 momentum " << KTsum1 <<
G4endl;
244 if ( ! success )
G4cout <<
"failed to correct E/p" <<
G4endl;
246 G4cout<<
"End of the Hadronization (G4ExcitedStringDecay)"<<
G4endl;
249 for (
unsigned int astring=0; astring < theStrings->size(); astring++)
251 if ( theStrings->operator[](astring)->IsExcited() )
253 Ptmp=theStrings->operator[](astring)->GetLeftParton()->Get4Momentum();
254 Ptmp.transform( toLab);
255 theStrings->operator[](astring)->GetLeftParton()->Set4Momentum(Ptmp);
257 Ptmp=theStrings->operator[](astring)->GetRightParton()->Get4Momentum();
258 Ptmp.transform( toLab);
259 theStrings->operator[](astring)->GetRightParton()->Set4Momentum(Ptmp);
263 Ptmp=theStrings->operator[](astring)->GetKineticTrack()->Get4Momentum();
264 Ptmp.transform( toLab);
265 theStrings->operator[](astring)->GetKineticTrack()->Set4Momentum(Ptmp);
275 const int nAttemptScale = 500;
276 const double ErrLimit = 1.E-5;
277 if (Output->empty())
return TRUE;
280 G4double TotalCollisionMass = TotalCollisionMom.m();
282 #ifdef debug_G4ExcitedStringCorr
283 G4cout<<
G4endl<<
"EnergyAndMomentumCorrector. Number of particles: "<<Output->size()<<
G4endl;
286 unsigned int cHadron;
287 for(cHadron = 0; cHadron < Output->size(); cHadron++)
289 SumMom += Output->operator[](cHadron)->Get4Momentum();
290 SumMass += Output->operator[](cHadron)->GetDefinition()->GetPDGMass();
293 #ifdef debug_G4ExcitedStringCorr
294 G4cout<<
"Sum part mom "<<SumMom<<
" "<<SumMom.mag()<<
G4endl
295 <<
"Sum str mom "<<TotalCollisionMom<<
" "<<TotalCollisionMom.mag()<<
G4endl;
296 G4cout<<
"SumMass TotalCollisionMass "<<SumMass<<
" "<<TotalCollisionMass<<
G4endl;
300 if (Output->size() < 2)
return FALSE;
302 if (SumMass > TotalCollisionMass)
return FALSE;
303 SumMass = SumMom.m2();
304 if (SumMass < 0)
return FALSE;
305 SumMass = std::sqrt(SumMass);
318 for(cAttempt = 0; cAttempt < nAttemptScale; cAttempt++)
321 for(cHadron = 0; cHadron < Output->size(); cHadron++)
323 G4LorentzVector HadronMom = Output->operator[](cHadron)->Get4Momentum();
324 HadronMom.setVect(Scale*HadronMom.vect());
325 G4double E = std::sqrt(HadronMom.vect().mag2() +
sqr(Output->operator[](cHadron)->GetDefinition()->GetPDGMass()));
327 Output->operator[](cHadron)->Set4Momentum(HadronMom);
330 Scale = TotalCollisionMass/Sum;
331 #ifdef debug_G4ExcitedStringCorr
332 G4cout <<
"Scale-1=" << Scale -1
333 <<
", TotalCollisionMass=" << TotalCollisionMass
337 if (std::fabs(Scale - 1) <= ErrLimit)
343 #ifdef debug_G4ExcitedStringCorr
346 G4cout <<
"G4ExcitedStringDecay::EnergyAndMomentumCorrector - Warning"<<
G4endl;
347 G4cout <<
" Scale not unity at end of iteration loop: "<<TotalCollisionMass<<
" "<<Sum<<
" "<<Scale<<
G4endl;
348 G4cout <<
" Number of secondaries: " << Output->size() <<
G4endl;
349 G4cout <<
" Wanted total energy: " << TotalCollisionMom.e() <<
G4endl;
350 G4cout <<
" Increase number of attempts or increase ERRLIMIT"<<
G4endl;
355 Beta = TotalCollisionMom.boostVector();
G4bool EnergyAndMomentumCorrector(G4KineticTrackVector *Output, G4LorentzVector &TotalCollisionMom)
std::vector< G4ExcitedString * > G4ExcitedStringVector
CLHEP::Hep3Vector G4ThreeVector
CLHEP::HepLorentzRotation G4LorentzRotation
virtual G4KineticTrackVector * FragmentString(const G4ExcitedString &theString)=0
virtual G4KineticTrackVector * FragmentStrings(const G4ExcitedStringVector *theStrings)
const G4ExcitedStringDecay & operator=(const G4ExcitedStringDecay &right)
G4VLongitudinalStringDecay * theStringDecay
const G4String & GetParticleName() const
G4KineticTrackVector * FragmentString(const G4ExcitedString &theString)
G4GLOB_DLL std::ostream G4cout
void SetPosition(const G4ThreeVector aPosition)
int operator!=(const G4ExcitedStringDecay &right) const
void Boost(G4ThreeVector &Velocity)
static const double perMillion
int operator==(const G4ExcitedStringDecay &right) const
const G4ParticleDefinition * GetDefinition() const
virtual ~G4ExcitedStringDecay()
CLHEP::HepLorentzVector G4LorentzVector