ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ExcitedStringDecay.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4ExcitedStringDecay.cc
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 #include "G4SampleResonance.hh"
33 
34 //#define debug_G4ExcitedStringDecay
35 //#define debug_G4ExcitedStringCorr
36 
37 
38 G4ExcitedStringDecay::G4ExcitedStringDecay() : G4VStringFragmentation(),theStringDecay(0)
39 {}
40 
41 
42 G4ExcitedStringDecay::G4ExcitedStringDecay(G4VLongitudinalStringDecay * aStringDecay)
43 : G4VStringFragmentation(),
44  theStringDecay(aStringDecay)
45 {}
46 
47 
48 G4ExcitedStringDecay::G4ExcitedStringDecay(const G4ExcitedStringDecay &)
49 : G4VStringFragmentation(),
50  theStringDecay(0)
51 {
52  throw G4HadronicException(__FILE__, __LINE__, "G4ExcitedStringDecay::copy ctor not accessible");
53 }
54 
55 
56 G4ExcitedStringDecay::~G4ExcitedStringDecay()
57 {
58 }
59 
60 
61 const G4ExcitedStringDecay & G4ExcitedStringDecay::operator=(const G4ExcitedStringDecay &)
62 {
63  throw G4HadronicException(__FILE__, __LINE__, "G4ExcitedStringDecay::operator= meant to not be accessable");
64  return *this;
65 }
66 
67 
68 G4bool G4ExcitedStringDecay::operator==(const G4ExcitedStringDecay &) const
69 {
70  return false;
71 }
72 
73 
74 G4bool G4ExcitedStringDecay::operator!=(const G4ExcitedStringDecay &) const
75 {
76  return true;
77 }
78 
79 
80 G4KineticTrackVector *G4ExcitedStringDecay::FragmentString(const G4ExcitedString &theString)
81 {
82  if ( theStringDecay == NULL ) theStringDecay=new G4LundStringFragmentation();
83  return theStringDecay->FragmentString(theString);
84 }
85 
86 
87 G4KineticTrackVector *G4ExcitedStringDecay::FragmentStrings(const G4ExcitedStringVector * theStrings)
88 {
89  G4LorentzVector KTsum(0.,0.,0.,0.);
90 
91  #ifdef debug_G4ExcitedStringDecay
92  G4cout<<G4endl;
93  G4cout<<"--------------------------- G4ExcitedStringDecay ----------------------"<<G4endl;
94  G4cout<<"Hadronization of Excited Strings: theStrings->size() "<<theStrings->size()<<G4endl;
95  #endif
96 
97  for ( unsigned int astring=0; astring < theStrings->size(); astring++)
98  // for ( unsigned int astring=0; astring < 1; astring++)
99  {
100  // G4cout<<"theStrings->operator[](astring)->IsExcited() "<<" "<<astring<<" "<<theStrings->operator[](astring)->IsExcited()<<G4endl;
101  if ( theStrings->operator[](astring)->IsExcited() )
102  {KTsum+= theStrings->operator[](astring)->Get4Momentum();}
103  else {KTsum+=theStrings->operator[](astring)->GetKineticTrack()->Get4Momentum();}
104  }
105 
106  G4LorentzRotation toCms( -1 * KTsum.boostVector() );
107  G4LorentzRotation toLab(toCms.inverse());
108  G4LorentzVector Ptmp;
109  KTsum=G4LorentzVector(0.,0.,0.,0.);
110 
111  for ( unsigned int astring=0; astring < theStrings->size(); astring++)
112  // for ( unsigned int astring=0; astring < 1; astring++)
113  {
114  if ( theStrings->operator[](astring)->IsExcited() )
115  {
116  Ptmp=toCms * theStrings->operator[](astring)->GetLeftParton()->Get4Momentum();
117  theStrings->operator[](astring)->GetLeftParton()->Set4Momentum(Ptmp);
118 
119  Ptmp=toCms * theStrings->operator[](astring)->GetRightParton()->Get4Momentum();
120  theStrings->operator[](astring)->GetRightParton()->Set4Momentum(Ptmp);
121 
122  KTsum+= theStrings->operator[](astring)->Get4Momentum();
123  }
124  else
125  {
126  Ptmp=toCms * theStrings->operator[](astring)->GetKineticTrack()->Get4Momentum();
127  theStrings->operator[](astring)->GetKineticTrack()->Set4Momentum(Ptmp);
128  KTsum+= theStrings->operator[](astring)->GetKineticTrack()->Get4Momentum();
129  }
130  }
131 
132  G4SampleResonance BrW;
133  const G4ParticleDefinition* TrackDefinition=0;
134 
135  G4KineticTrackVector * theResult = new G4KineticTrackVector;
136  G4int attempts(0);
137  G4bool success=false;
138  G4bool NeedEnergyCorrector=false;
139  do {
140  #ifdef debug_G4ExcitedStringDecay
141  G4cout<<"New try No "<<attempts<<" to hadronize strings"<<G4endl;
142  #endif
143 
144  std::for_each(theResult->begin() , theResult->end() , DeleteKineticTrack());
145  theResult->clear();
146 
147  attempts++;
148 
149  G4LorentzVector KTsecondaries(0.,0.,0.,0.);
150  NeedEnergyCorrector=false;
151 
152  for ( unsigned int astring=0; astring < theStrings->size(); astring++)
153  // for ( unsigned int astring=0; astring < 1; astring++) // Proj. Last Str. Decay for FTF
154  // for ( unsigned int astring=1; astring < 2; astring++) // Proj. Last Str. Decay for QGS
155  // for ( unsigned int astring=0; astring < 1; astring++) // For testing purposes
156  {
157  #ifdef debug_G4ExcitedStringDecay
158  G4cout<<"String No "<<astring+1<<" Excited? "<<theStrings->operator[](astring)->IsExcited()<<G4endl;
159 
160  G4cout<<"String No "<<astring+1<<" 4Momentum "<<theStrings->operator[](astring)->Get4Momentum()
161  <<" "<<theStrings->operator[](astring)->Get4Momentum().mag()<<G4endl;
162  #endif
163 
164  G4KineticTrackVector * generatedKineticTracks = NULL;
165  if ( theStrings->operator[](astring)->IsExcited() )
166  {
167  #ifdef debug_G4ExcitedStringDecay
168  G4cout<<"Fragment String with partons: "
169  <<theStrings->operator[](astring)->GetLeftParton()->GetPDGcode() <<" "
170  <<theStrings->operator[](astring)->GetRightParton()->GetPDGcode()<<" "
171  <<"Direction "<<theStrings->operator[](astring)->GetDirection()<<G4endl;
172  #endif
173  generatedKineticTracks=FragmentString(*theStrings->operator[](astring));
174  #ifdef debug_G4ExcitedStringDecay
175  G4cout<<"(G4ExcitedStringDecay) Number of produced hadrons = "
176  <<generatedKineticTracks->size()<<G4endl;
177  #endif
178  } else {
179  #ifdef debug_G4ExcitedStringDecay
180  G4cout<<" GetTrack from the String"<<G4endl;
181  #endif
182  G4LorentzVector Mom=theStrings->operator[](astring)->GetKineticTrack()->Get4Momentum();
183  G4KineticTrack * aTrack= new G4KineticTrack(
184  theStrings->operator[](astring)->GetKineticTrack()->GetDefinition(),
185  theStrings->operator[](astring)->GetKineticTrack()->GetFormationTime(),
186  G4ThreeVector(0), Mom);
187 
188  aTrack->SetPosition(theStrings->operator[](astring)->GetKineticTrack()->GetPosition());
189 
190  #ifdef debug_G4ExcitedStringDecay
191  G4cout<<" A particle stored in the track is "<<aTrack->GetDefinition()->GetParticleName()<<G4endl;
192  #endif
193 
194  generatedKineticTracks = new G4KineticTrackVector;
195  generatedKineticTracks->push_back(aTrack);
196  }
197 
198  if (generatedKineticTracks->size() == 0)
199  {
200  // G4cerr << "G4VPartonStringModel:No KineticTracks produced" << G4endl;
201  // continue;
202  success=false; NeedEnergyCorrector=false; break;
203  }
204 
205  G4LorentzVector KTsum1(0.,0.,0.,0.);
206  for ( unsigned int aTrack=0; aTrack<generatedKineticTracks->size();aTrack++)
207  {
208  #ifdef debug_G4ExcitedStringDecay
209  G4cout<<"Prod part No. "<<aTrack+1<<" "
210  <<(*generatedKineticTracks)[aTrack]->GetDefinition()->GetParticleName()<<" "
211  <<(*generatedKineticTracks)[aTrack]->Get4Momentum()
212  <<(*generatedKineticTracks)[aTrack]->Get4Momentum().mag()<<G4endl;
213  #endif
214  // --------------- Sampling mass of unstable hadronic resonances ----------------
215  TrackDefinition = (*generatedKineticTracks)[aTrack]->GetDefinition();
216 
217  if (TrackDefinition->IsShortLived())
218  {
219  G4double NewTrackMass =
220  BrW.SampleMass( TrackDefinition->GetPDGMass(), TrackDefinition->GetPDGWidth(),
221  BrW.GetMinimumMass( TrackDefinition ) + 10.0*MeV,
222  TrackDefinition->GetPDGMass() + 5.0*TrackDefinition->GetPDGWidth() );
223  G4LorentzVector Tmp=G4LorentzVector((*generatedKineticTracks)[aTrack]->Get4Momentum());
224  Tmp.setE(std::sqrt(sqr(NewTrackMass) + Tmp.vect().mag2()));
225 
226  (*generatedKineticTracks)[aTrack]->Set4Momentum(Tmp);
227 
228  #ifdef debug_G4ExcitedStringDecay
229  G4cout<<"Resonance *** "<<aTrack+1<<" "
230  <<(*generatedKineticTracks)[aTrack]->GetDefinition()->GetParticleName()<<" "
231  <<(*generatedKineticTracks)[aTrack]->Get4Momentum()
232  <<(*generatedKineticTracks)[aTrack]->Get4Momentum().mag()<<G4endl;
233  #endif
234  }
235  //-------------------------------------------------------------------------------
236 
237  theResult->push_back(generatedKineticTracks->operator[](aTrack));
238  KTsum1+= (*generatedKineticTracks)[aTrack]->Get4Momentum();
239  }
240  KTsecondaries+=KTsum1;
241 
242  #ifdef debug_G4ExcitedStringDecay
243  G4cout << "String secondaries(" <<generatedKineticTracks->size()<< ")"<<G4endl
244  <<"Init string momentum: "<< theStrings->operator[](astring)->Get4Momentum()<<G4endl
245  <<"Final hadrons momentum: "<< KTsum1 << G4endl;
246  #endif
247 
248  if ( KTsum1.e() > 0 &&
249  std::abs((KTsum1.e()-theStrings->operator[](astring)->Get4Momentum().e()) / KTsum1.e()) > perMillion )
250  {
251  NeedEnergyCorrector=true;
252  }
253 
254  #ifdef debug_G4ExcitedStringDecay
255  G4cout<<"NeedEnergyCorrection yes/no "<<NeedEnergyCorrector<<G4endl;
256  #endif
257 
258  // clean up
259  delete generatedKineticTracks;
260  success=true;
261  }
262 
263  if ( NeedEnergyCorrector ) success=EnergyAndMomentumCorrector(theResult, KTsum);
264  } while (!success && (attempts < 100)); /* Loop checking, 07.08.2015, A.Ribon */
265 
266  for ( unsigned int aTrack=0; aTrack<theResult->size();aTrack++)
267  {
268  Ptmp=(*theResult)[aTrack]->Get4Momentum();
269  Ptmp.transform( toLab);
270  (*theResult)[aTrack]->Set4Momentum(Ptmp);
271  }
272 
273  #ifdef debug_G4ExcitedStringDecay
274  G4cout<<"End of the strings fragmentation (G4ExcitedStringDecay)"<<G4endl;
275 
276  G4LorentzVector KTsum1(0.,0.,0.,0.);
277 
278  for ( unsigned int aTrack=0; aTrack<theResult->size();aTrack++)
279  {
280  G4cout << " corrected tracks .. " << (*theResult)[aTrack]->GetDefinition()->GetParticleName()
281  <<" " << (*theResult)[aTrack]->Get4Momentum()
282  <<" " << (*theResult)[aTrack]->Get4Momentum().mag()<< G4endl;
283  KTsum1+= (*theResult)[aTrack]->Get4Momentum();
284  }
285 
286  G4cout << "Needcorrector/success " << NeedEnergyCorrector << "/" << success
287  << ", Corrected total 4 momentum " << KTsum1 << G4endl;
288  if ( ! success ) G4cout << "failed to correct E/p" << G4endl;
289 
290  G4cout<<"End of the Hadronization (G4ExcitedStringDecay)"<<G4endl;
291  #endif
292 
293  if (!success)
294  {
295  if (theResult->size() != 0)
296  {
297  std::for_each(theResult->begin() , theResult->end() , DeleteKineticTrack());
298  theResult->clear();
299  delete theResult; theResult=0;
300  }
301  for ( unsigned int astring=0; astring < theStrings->size(); astring++)
302  // for ( unsigned int astring=0; astring < 1; astring++) // Need more correct. For testing purposes.
303  {
304  if ( theStrings->operator[](astring)->IsExcited() )
305  {
306  Ptmp=theStrings->operator[](astring)->GetLeftParton()->Get4Momentum();
307  Ptmp.transform( toLab);
308  theStrings->operator[](astring)->GetLeftParton()->Set4Momentum(Ptmp);
309 
310  Ptmp=theStrings->operator[](astring)->GetRightParton()->Get4Momentum();
311  Ptmp.transform( toLab);
312  theStrings->operator[](astring)->GetRightParton()->Set4Momentum(Ptmp);
313  }
314  else
315  {
316  Ptmp=theStrings->operator[](astring)->GetKineticTrack()->Get4Momentum();
317  Ptmp.transform( toLab);
318  theStrings->operator[](astring)->GetKineticTrack()->Set4Momentum(Ptmp);
319  }
320  }
321  }
322  return theResult;
323 }
324 
325 
326 G4bool G4ExcitedStringDecay::
327 EnergyAndMomentumCorrector(G4KineticTrackVector* Output, G4LorentzVector& TotalCollisionMom)
328 {
329  const int nAttemptScale = 500;
330  const double ErrLimit = 1.E-5;
331  if (Output->empty()) return TRUE;
332  G4LorentzVector SumMom;
333  G4double SumMass = 0;
334  G4double TotalCollisionMass = TotalCollisionMom.m();
335 
336  std::vector<G4double> HadronMass; G4double HadronM(0.);
337 
338  #ifdef debug_G4ExcitedStringCorr
339  G4cout<<G4endl<<"EnergyAndMomentumCorrector. Number of particles: "<<Output->size()<<G4endl;
340  #endif
341  // Calculate sum hadron 4-momenta and summing hadron mass
342  unsigned int cHadron;
343  for (cHadron = 0; cHadron < Output->size(); cHadron++)
344  {
345  SumMom += Output->operator[](cHadron)->Get4Momentum();
346  HadronM=Output->operator[](cHadron)->Get4Momentum().mag(); HadronMass.push_back(HadronM);
347  SumMass += Output->operator[](cHadron)->Get4Momentum().mag(); //GetDefinition()->GetPDGMass();
348  }
349 
350  #ifdef debug_G4ExcitedStringCorr
351  G4cout<<"Sum part mom "<<SumMom<<" "<<SumMom.mag()<<G4endl
352  <<"Sum str mom "<<TotalCollisionMom<<" "<<TotalCollisionMom.mag()<<G4endl;
353  G4cout<<"SumMass TotalCollisionMass "<<SumMass<<" "<<TotalCollisionMass<<G4endl;
354  #endif
355 
356  // Cannot correct a single particle
357  if (Output->size() < 2) return FALSE;
358 
359  if (SumMass > TotalCollisionMass) return FALSE;
360  SumMass = SumMom.m2();
361  if (SumMass < 0) return FALSE;
362  SumMass = std::sqrt(SumMass);
363 
364  // Compute c.m.s. hadron velocity and boost KTV to hadron c.m.s.
365  // G4ThreeVector Beta = -SumMom.boostVector();
366  G4ThreeVector Beta = -TotalCollisionMom.boostVector();
367  Output->Boost(Beta);
368 
369  // Scale total c.m.s. hadron energy (hadron system mass).
370  // It should be equal interaction mass
371  G4double Scale = 1;
372  G4int cAttempt = 0;
373  G4double Sum = 0;
374  G4bool success = false;
375  for (cAttempt = 0; cAttempt < nAttemptScale; cAttempt++)
376  {
377  Sum = 0;
378  for (cHadron = 0; cHadron < Output->size(); cHadron++)
379  {
380  HadronM = HadronMass.at(cHadron);
381  G4LorentzVector HadronMom = Output->operator[](cHadron)->Get4Momentum();
382  HadronMom.setVect(Scale*HadronMom.vect());
383  G4double E = std::sqrt(HadronMom.vect().mag2() + sqr(HadronM));
384  //sqr(Output->operator[](cHadron)->GetDefinition()->GetPDGMass()));
385  HadronMom.setE(E);
386  Output->operator[](cHadron)->Set4Momentum(HadronMom);
387  Sum += E;
388  }
389  Scale = TotalCollisionMass/Sum;
390  #ifdef debug_G4ExcitedStringCorr
391  G4cout << "Scale-1=" << Scale -1
392  << ", TotalCollisionMass=" << TotalCollisionMass
393  << ", Sum=" << Sum
394  << G4endl;
395  #endif
396  if (std::fabs(Scale - 1) <= ErrLimit)
397  {
398  success = true;
399  break;
400  }
401  }
402 
403  #ifdef debug_G4ExcitedStringCorr
404  if (!success)
405  {
406  G4cout << "G4ExcitedStringDecay::EnergyAndMomentumCorrector - Warning"<<G4endl;
407  G4cout << " Scale not unity at end of iteration loop: "<<TotalCollisionMass<<" "<<Sum<<" "<<Scale<<G4endl;
408  G4cout << " Number of secondaries: " << Output->size() << G4endl;
409  G4cout << " Wanted total energy: " << TotalCollisionMom.e() << G4endl;
410  G4cout << " Increase number of attempts or increase ERRLIMIT"<<G4endl;
411  // throw G4HadronicException(__FILE__, __LINE__, "G4ExcitedStringDecay failed to correct...");
412  }
413  #endif
414 
415  // Compute c.m.s. interaction velocity and KTV back boost
416  Beta = TotalCollisionMom.boostVector();
417  Output->Boost(Beta);
418 
419  return success;
420 }
421