Geant4  10.02.p01
G4ExcitedStringDecay.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // Historic fragment from M.Komogorov; clean-up still necessary @@@
27 
28 #include "G4ExcitedStringDecay.hh"
29 #include "G4SystemOfUnits.hh"
30 #include "G4KineticTrack.hh"
31 
32 //#define debug_G4ExcitedStringDecay
33 //#define debug_G4ExcitedStringCorr
34 
36 {}
37 
40  theStringDecay(aStringDecay)
41 {}
42 
45  theStringDecay(0)
46 {
47  throw G4HadronicException(__FILE__, __LINE__, "G4ExcitedStringDecay::copy ctor not accessible");
48 }
49 
51 {
52 }
53 
55 {
56  throw G4HadronicException(__FILE__, __LINE__, "G4ExcitedStringDecay::operator= meant to not be accessable");
57  return *this;
58 }
59 
61 {
62  return 0;
63 }
64 
66 {
67  return 1;
68 }
69 
71 {
73  return theStringDecay->FragmentString(theString);
74 }
75 
77 {
78  G4LorentzVector KTsum(0.,0.,0.,0.);
79 
80 #ifdef debug_G4ExcitedStringDecay
81  G4cout<<G4endl;
82  G4cout<<"--------------------------- G4ExcitedStringDecay ----------------------"<<G4endl;
83  G4cout<<"Hadronization of Excited Strings: theStrings->size() "<<theStrings->size()<<G4endl;
84 #endif
85 
86  for ( unsigned int astring=0; astring < theStrings->size(); astring++)
87  {
88  if ( theStrings->operator[](astring)->IsExcited() )
89  {KTsum+= theStrings->operator[](astring)->Get4Momentum();}
90  else {KTsum+=theStrings->operator[](astring)->GetKineticTrack()->Get4Momentum();}
91  }
92 
93  G4LorentzRotation toCms( -1 * KTsum.boostVector() ); // Uzhi 22 June 2014
94  G4LorentzRotation toLab(toCms.inverse()); // Uzhi 22 June 2014
95  G4LorentzVector Ptmp;
96  KTsum=G4LorentzVector(0.,0.,0.,0.);
97  for ( unsigned int astring=0; astring < theStrings->size(); astring++) // Uzhi 22 June 2014
98  {
99  if ( theStrings->operator[](astring)->IsExcited() )
100  {
101  Ptmp=toCms * theStrings->operator[](astring)->GetLeftParton()->Get4Momentum();
102  theStrings->operator[](astring)->GetLeftParton()->Set4Momentum(Ptmp);
103 
104  Ptmp=toCms * theStrings->operator[](astring)->GetRightParton()->Get4Momentum();
105  theStrings->operator[](astring)->GetRightParton()->Set4Momentum(Ptmp);
106 
107  KTsum+= theStrings->operator[](astring)->Get4Momentum();
108  }
109  else
110  {
111  Ptmp=toCms * theStrings->operator[](astring)->GetKineticTrack()->Get4Momentum();
112  theStrings->operator[](astring)->GetKineticTrack()->Set4Momentum(Ptmp);
113  KTsum+= theStrings->operator[](astring)->GetKineticTrack()->Get4Momentum();
114  }
115  } // Uzhi 22 June 2014
116 
117  G4KineticTrackVector * theResult = new G4KineticTrackVector;
118  G4int attempts(0);
119  G4bool success=false;
120  G4bool NeedEnergyCorrector=false;
121  do {
122 #ifdef debug_G4ExcitedStringDecay
123  G4cout<<"New try No "<<attempts<<" to hadronize strings"<<G4endl;
124 #endif
125 
126  std::for_each(theResult->begin() , theResult->end() , DeleteKineticTrack());
127  theResult->clear();
128 
129  attempts++;
130 
131  G4LorentzVector KTsecondaries(0.,0.,0.,0.);
132  NeedEnergyCorrector=false;
133 
134  for ( unsigned int astring=0; astring < theStrings->size(); astring++)
135 // for ( unsigned int astring=0; astring < 1; astring++)
136  {
137 #ifdef debug_G4ExcitedStringDecay
138  G4cout<<"String No "<<astring+1<<" Excited? "<<theStrings->operator[](astring)->IsExcited()<<G4endl;
139 
140  G4cout<<"String No "<<astring+1<<" 4Momentum "<<theStrings->operator[](astring)->Get4Momentum()
141  <<" "<<theStrings->operator[](astring)->Get4Momentum().mag()<<G4endl;
142 #endif
143 
144  G4KineticTrackVector * generatedKineticTracks = NULL;
145  if ( theStrings->operator[](astring)->IsExcited() )
146  {
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;
152 #endif
153  generatedKineticTracks=FragmentString(*theStrings->operator[](astring));
154 #ifdef debug_G4ExcitedStringDecay
155  G4cout<<"(G4ExcitedStringDecay) Number of produced hadrons = "<<generatedKineticTracks->size()<<G4endl;
156 #endif
157  } else {
158 #ifdef debug_G4ExcitedStringDecay
159  G4cout<<" GetTrack from the String"<<G4endl;
160 #endif
161  G4LorentzVector Mom=theStrings->operator[](astring)->GetKineticTrack()->Get4Momentum();
162  G4KineticTrack * aTrack= new G4KineticTrack(
163  theStrings->operator[](astring)->GetKineticTrack()->GetDefinition(),
164  theStrings->operator[](astring)->GetKineticTrack()->GetFormationTime(),
165  G4ThreeVector(0), Mom);
166 
167  aTrack->SetPosition(theStrings->operator[](astring)->GetKineticTrack()->GetPosition());
168 
169 #ifdef debug_G4ExcitedStringDecay
170  G4cout<<" A particle stored in the track is "<<aTrack->GetDefinition()->GetParticleName()<<G4endl;
171 #endif
172 
173  generatedKineticTracks = new G4KineticTrackVector;
174  generatedKineticTracks->push_back(aTrack);
175  }
176 
177 // if (generatedKineticTracks == NULL) // Uzhi 02.06.2015
178  if (generatedKineticTracks->size() == 0) // Uzhi 02.06.2015
179  {
180 // G4cerr << "G4VPartonStringModel:No KineticTracks produced" << G4endl; // Uzhi 02.06.2015
181 // continue; // Uzhi 02.06.2015
182  success=false; NeedEnergyCorrector=false; break; // Uzhi 02.06.2015
183  }
184 
185  G4LorentzVector KTsum1(0.,0.,0.,0.);
186  for ( unsigned int aTrack=0; aTrack<generatedKineticTracks->size();aTrack++)
187  {
188 #ifdef debug_G4ExcitedStringDecay
189  G4cout<<"Prod part No. "<<aTrack+1<<" "
190  <<(*generatedKineticTracks)[aTrack]->GetDefinition()->GetParticleName()<<" "
191  <<(*generatedKineticTracks)[aTrack]->Get4Momentum()<<G4endl;
192 #endif
193  theResult->push_back(generatedKineticTracks->operator[](aTrack));
194  KTsum1+= (*generatedKineticTracks)[aTrack]->Get4Momentum();
195  }
196  KTsecondaries+=KTsum1;
197 
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;
202 #endif
203 
204  if ( KTsum1.e() > 0 && std::abs((KTsum1.e()-theStrings->operator[](astring)->Get4Momentum().e()) / KTsum1.e()) > perMillion )
205  {
206  NeedEnergyCorrector=true;
207  }
208 
209 #ifdef debug_G4ExcitedStringDecay
210  G4cout<<"NeedEnergyCorrection yes/no "<<NeedEnergyCorrector<<G4endl;
211 #endif
212 
213 // clean up
214  delete generatedKineticTracks;
215  success=true; // Uzhi 02.06.2015
216  }
217 
218 // success=true; // Uzhi 02.06.2015
219 
220  if ( NeedEnergyCorrector ) success=EnergyAndMomentumCorrector(theResult, KTsum);
221  } while(!success && (attempts < 10)); /* Loop checking, 07.08.2015, A.Ribon */
222 
223  for ( unsigned int aTrack=0; aTrack<theResult->size();aTrack++) // Uzhi 22 June 2014
224  {
225  Ptmp=(*theResult)[aTrack]->Get4Momentum();
226  Ptmp.transform( toLab);
227  (*theResult)[aTrack]->Set4Momentum(Ptmp);
228  } // Uzhi 22 June 2014
229 
230 #ifdef debug_G4ExcitedStringDecay
231  G4cout<<"End of the strings fragmentation (G4ExcitedStringDecay)"<<G4endl;
232 
233  G4LorentzVector KTsum1(0.,0.,0.,0.);
234 
235  for ( unsigned int aTrack=0; aTrack<theResult->size();aTrack++)
236  {
237  G4cout << " corrected tracks .. " << (*theResult)[aTrack]->GetDefinition()->GetParticleName()
238  <<" " << (*theResult)[aTrack]->Get4Momentum()
239  <<" " << (*theResult)[aTrack]->Get4Momentum().mag()<< G4endl;
240  KTsum1+= (*theResult)[aTrack]->Get4Momentum();
241  }
242 
243  G4cout << "Needcorrector/success " << NeedEnergyCorrector << "/" << success << ", Corrected total 4 momentum " << KTsum1 << G4endl;
244  if ( ! success ) G4cout << "failed to correct E/p" << G4endl;
245 
246  G4cout<<"End of the Hadronization (G4ExcitedStringDecay)"<<G4endl;
247 #endif
248 
249  for ( unsigned int astring=0; astring < theStrings->size(); astring++) // Uzhi 24 Oct. 2014
250  {
251  if ( theStrings->operator[](astring)->IsExcited() )
252  {
253  Ptmp=theStrings->operator[](astring)->GetLeftParton()->Get4Momentum();
254  Ptmp.transform( toLab);
255  theStrings->operator[](astring)->GetLeftParton()->Set4Momentum(Ptmp);
256 
257  Ptmp=theStrings->operator[](astring)->GetRightParton()->Get4Momentum();
258  Ptmp.transform( toLab);
259  theStrings->operator[](astring)->GetRightParton()->Set4Momentum(Ptmp);
260  }
261  else
262  {
263  Ptmp=theStrings->operator[](astring)->GetKineticTrack()->Get4Momentum();
264  Ptmp.transform( toLab);
265  theStrings->operator[](astring)->GetKineticTrack()->Set4Momentum(Ptmp);
266  }
267  } // Uzhi 24 Oct. 2014
268 
269  return theResult;
270 }
271 
273  (G4KineticTrackVector* Output, G4LorentzVector& TotalCollisionMom)
274  {
275  const int nAttemptScale = 500;
276  const double ErrLimit = 1.E-5;
277  if (Output->empty()) return TRUE;
278  G4LorentzVector SumMom;
279  G4double SumMass = 0;
280  G4double TotalCollisionMass = TotalCollisionMom.m();
281 
282 #ifdef debug_G4ExcitedStringCorr
283  G4cout<<G4endl<<"EnergyAndMomentumCorrector. Number of particles: "<<Output->size()<<G4endl;
284 #endif
285  // Calculate sum hadron 4-momenta and summing hadron mass
286  unsigned int cHadron;
287  for(cHadron = 0; cHadron < Output->size(); cHadron++)
288  {
289  SumMom += Output->operator[](cHadron)->Get4Momentum();
290  SumMass += Output->operator[](cHadron)->GetDefinition()->GetPDGMass();
291  }
292 
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;
297 #endif
298 
299  // Cannot correct a single particle
300  if (Output->size() < 2) return FALSE;
301 
302  if (SumMass > TotalCollisionMass) return FALSE;
303  SumMass = SumMom.m2();
304  if (SumMass < 0) return FALSE;
305  SumMass = std::sqrt(SumMass);
306 
307  // Compute c.m.s. hadron velocity and boost KTV to hadron c.m.s.
308 // G4ThreeVector Beta = -SumMom.boostVector(); // Uzhi 22 June 2014
309  G4ThreeVector Beta = -TotalCollisionMom.boostVector(); // Uzhi 22 June 2014
310  Output->Boost(Beta);
311 
312  // Scale total c.m.s. hadron energy (hadron system mass).
313  // It should be equal interaction mass
314  G4double Scale = 1;
315  G4int cAttempt = 0;
316  G4double Sum = 0;
317  G4bool success = false;
318  for(cAttempt = 0; cAttempt < nAttemptScale; cAttempt++)
319  {
320  Sum = 0;
321  for(cHadron = 0; cHadron < Output->size(); cHadron++)
322  {
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()));
326  HadronMom.setE(E);
327  Output->operator[](cHadron)->Set4Momentum(HadronMom);
328  Sum += E;
329  }
330  Scale = TotalCollisionMass/Sum;
331 #ifdef debug_G4ExcitedStringCorr
332  G4cout << "Scale-1=" << Scale -1
333  << ", TotalCollisionMass=" << TotalCollisionMass
334  << ", Sum=" << Sum
335  << G4endl;
336 #endif
337  if (std::fabs(Scale - 1) <= ErrLimit)
338  {
339  success = true;
340  break;
341  }
342  }
343 #ifdef debug_G4ExcitedStringCorr
344  if(!success)
345  {
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;
351 // throw G4HadronicException(__FILE__, __LINE__, "G4ExcitedStringDecay failed to correct...");
352  }
353 #endif
354  // Compute c.m.s. interaction velocity and KTV back boost
355  Beta = TotalCollisionMom.boostVector();
356  Output->Boost(Beta);
357 
358  return success;
359  }
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
int G4int
Definition: G4Types.hh:78
const G4String & GetParticleName() const
G4KineticTrackVector * FragmentString(const G4ExcitedString &theString)
G4GLOB_DLL std::ostream G4cout
bool G4bool
Definition: G4Types.hh:79
void SetPosition(const G4ThreeVector aPosition)
int operator!=(const G4ExcitedStringDecay &right) const
#define FALSE
Definition: globals.hh:52
#define TRUE
Definition: globals.hh:55
void Boost(G4ThreeVector &Velocity)
static const double perMillion
Definition: G4SIunits.hh:331
#define G4endl
Definition: G4ios.hh:61
T sqr(const T &x)
Definition: templates.hh:145
double G4double
Definition: G4Types.hh:76
int operator==(const G4ExcitedStringDecay &right) const
const G4ParticleDefinition * GetDefinition() const
CLHEP::HepLorentzVector G4LorentzVector