Geant4  10.00.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 
33  theStringDecay(0)
34 {}
35 
38  theStringDecay(aStringDecay)
39 {}
40 
43  theStringDecay(0)
44 {
45  throw G4HadronicException(__FILE__, __LINE__, "G4ExcitedStringDecay::copy ctor not accessible");
46 }
47 
49 {
50 }
51 
53 {
54  throw G4HadronicException(__FILE__, __LINE__, "G4ExcitedStringDecay::operator= meant to not be accessable");
55  return *this;
56 }
57 
59 {
60  return 0;
61 }
62 
64 {
65  return 1;
66 }
67 
69  (const G4ExcitedString &theString)
70 {
71  if ( theStringDecay == NULL )
72 
73  theStringDecay=new G4LundStringFragmentation();
74 
75  return theStringDecay->FragmentString(theString);
76 }
77 
79  (const G4ExcitedStringVector * theStrings)
80 {
81  G4LorentzVector KTsum(0.,0.,0.,0.);
82 
83 //G4cout<<"theStrings->size() "<<theStrings->size()<<G4endl;
84  for ( unsigned int astring=0; astring < theStrings->size(); astring++)
85  {
86  KTsum+= theStrings->operator[](astring)->Get4Momentum();
87  }
88 
90  G4int attempts(0);
91  G4bool success=false;
92  G4bool NeedEnergyCorrector=false;
93  do {
94  //G4cout<<"Check of momentum at string fragmentations. New try."<<G4endl;
95  std::for_each(theResult->begin() , theResult->end() , DeleteKineticTrack());
96  theResult->clear();
97 
98  attempts++;
99  //G4cout<<G4endl<<"attempts "<<attempts<<G4endl;
100  G4LorentzVector KTsecondaries(0.,0.,0.,0.);
101  NeedEnergyCorrector=false;
102 
103  for ( unsigned int astring=0; astring < theStrings->size(); astring++)
104  {
105  //G4cout<<"String No "<<astring+1<<" "<<theStrings->operator[](astring)->Get4Momentum().mag2()<<" "<<theStrings->operator[](astring)->GetRightParton()->GetPDGcode()<<" "<<theStrings->operator[](astring)->GetLeftParton()->GetPDGcode()<<" "<<theStrings->operator[](astring)->Get4Momentum()<<G4endl;
106  //G4int Uzhi; G4cin >>Uzhi;
107  G4KineticTrackVector * generatedKineticTracks = NULL;
108  if ( theStrings->operator[](astring)->IsExcited() )
109  {
110  //G4cout<<"Fragment String"<<G4endl;
111  generatedKineticTracks=FragmentString(*theStrings->operator[](astring));
112  } else {
113  generatedKineticTracks = new G4KineticTrackVector;
114  generatedKineticTracks->push_back(theStrings->operator[](astring)->GetKineticTrack());
115  }
116 
117  if (generatedKineticTracks == NULL)
118  {
119  G4cerr << "G4VPartonStringModel:No KineticTracks produced" << G4endl;
120  continue;
121  }
122 
123  G4LorentzVector KTsum1(0.,0.,0.,0.);
124  for ( unsigned int aTrack=0; aTrack<generatedKineticTracks->size();aTrack++)
125  {
126  //--debug-- G4cout<<"Prod part "<<(*generatedKineticTracks)[aTrack]->GetDefinition()->GetParticleName()<<" "<<(*generatedKineticTracks)[aTrack]->Get4Momentum()<<G4endl;
127  theResult->push_back(generatedKineticTracks->operator[](aTrack));
128  KTsum1+= (*generatedKineticTracks)[aTrack]->Get4Momentum();
129  }
130  KTsecondaries+=KTsum1;
131 
132  //--debug--G4cout << "String secondaries(" <<generatedKineticTracks->size()<< ") momentum: "
133  //--debug--<< theStrings->operator[](astring)->Get4Momentum() << " " << KTsum1 << G4endl;
134  if ( KTsum1.e() > 0 && std::abs((KTsum1.e()-theStrings->operator[](astring)->Get4Momentum().e()) / KTsum1.e()) > perMillion )
135  {
136  //--debug-- G4cout << "String secondaries(" <<generatedKineticTracks->size()<< ") momentum: "
137  //--debug-- << theStrings->operator[](astring)->Get4Momentum() << " " << KTsum1 << G4endl;
138  NeedEnergyCorrector=true;
139  }
140 
141 // clean up
142  delete generatedKineticTracks;
143  }
144  //--debug G4cout << "Initial Strings / secondaries total 4 momentum " << KTsum << " " <<KTsecondaries << G4endl;
145 
146  success=true;
147  //G4cout<<"success "<<success<<G4endl;
148  if ( NeedEnergyCorrector ) success=EnergyAndMomentumCorrector(theResult, KTsum);
149  //G4cout<<"success after Ecorr "<<success<<G4endl;
150  } while(!success && (attempts < 10)); // It was 100 !!! Uzhi
151  //G4cout<<"End frag string"<<G4endl;
152 
153 #ifdef debug_ExcitedStringDecay
154  G4LorentzVector KTsum1=0;
155  for ( unsigned int aTrack=0; aTrack<theResult->size();aTrack++)
156  {
157  G4cout << " corrected tracks .. " << (*theResult)[aTrack]->GetDefinition()->GetParticleName()
158  <<" " << (*theResult)[aTrack]->Get4Momentum() << G4endl;
159  KTsum1+= (*theResult)[aTrack]->Get4Momentum();
160  }
161 
162  G4cout << "Needcorrector/success " << NeedEnergyCorrector << "/" << success << ", Corrected total 4 momentum " << KTsum1 << G4endl;
163  if ( ! success ) G4cout << "failed to correct E/p" << G4endl;
164 #endif
165 
166  return theResult;
167 }
168 
170  (G4KineticTrackVector* Output, G4LorentzVector& TotalCollisionMom)
171  {
172  const int nAttemptScale = 500;
173  const double ErrLimit = 1.E-5;
174  if (Output->empty())
175  return TRUE;
176  G4LorentzVector SumMom;
177  G4double SumMass = 0;
178  G4double TotalCollisionMass = TotalCollisionMom.m();
179 
180 //G4cout<<G4endl<<"EnergyAndMomentumCorrector "<<Output->size()<<G4endl;
181  // Calculate sum hadron 4-momenta and summing hadron mass
182  unsigned int cHadron;
183  for(cHadron = 0; cHadron < Output->size(); cHadron++)
184  {
185  SumMom += Output->operator[](cHadron)->Get4Momentum();
186  SumMass += Output->operator[](cHadron)->GetDefinition()->GetPDGMass();
187  }
188 
189 //G4cout<<"SumMass TotalCollisionMass "<<SumMass<<" "<<TotalCollisionMass<<G4endl;
190 
191  // Cannot correct a single particle
192  if (Output->size() < 2) return FALSE;
193 
194  if (SumMass > TotalCollisionMass) return FALSE;
195  SumMass = SumMom.m2();
196  if (SumMass < 0) return FALSE;
197  SumMass = std::sqrt(SumMass);
198 
199  // Compute c.m.s. hadron velocity and boost KTV to hadron c.m.s.
200  G4ThreeVector Beta = -SumMom.boostVector();
201  Output->Boost(Beta);
202 
203  // Scale total c.m.s. hadron energy (hadron system mass).
204  // It should be equal interaction mass
205  G4double Scale = 1;
206  G4int cAttempt = 0;
207  G4double Sum = 0;
208  G4bool success = false;
209  for(cAttempt = 0; cAttempt < nAttemptScale; cAttempt++)
210  {
211  Sum = 0;
212  for(cHadron = 0; cHadron < Output->size(); cHadron++)
213  {
214  G4LorentzVector HadronMom = Output->operator[](cHadron)->Get4Momentum();
215  HadronMom.setVect(Scale*HadronMom.vect());
216  G4double E = std::sqrt(HadronMom.vect().mag2() + sqr(Output->operator[](cHadron)->GetDefinition()->GetPDGMass()));
217  HadronMom.setE(E);
218  Output->operator[](cHadron)->Set4Momentum(HadronMom);
219  Sum += E;
220  }
221  Scale = TotalCollisionMass/Sum;
222 #ifdef debug_G4ExcitedStringDecay
223  G4cout << "Scale-1=" << Scale -1
224  << ", TotalCollisionMass=" << TotalCollisionMass
225  << ", Sum=" << Sum
226  << G4endl;
227 #endif
228  if (std::fabs(Scale - 1) <= ErrLimit)
229  {
230  success = true;
231  break;
232  }
233  }
234 #ifdef debug_G4ExcitedStringDecay
235  if(!success)
236  {
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;
242 // throw G4HadronicException(__FILE__, __LINE__, "G4ExcitedStringDecay failed to correct...");
243  }
244 #endif
245  // Compute c.m.s. interaction velocity and KTV back boost
246  Beta = TotalCollisionMom.boostVector();
247  Output->Boost(Beta);
248 
249  return success;
250  }
G4bool EnergyAndMomentumCorrector(G4KineticTrackVector *Output, G4LorentzVector &TotalCollisionMom)
std::vector< G4ExcitedString * > G4ExcitedStringVector
CLHEP::Hep3Vector G4ThreeVector
virtual G4KineticTrackVector * FragmentStrings(const G4ExcitedStringVector *theStrings)
const G4ExcitedStringDecay & operator=(const G4ExcitedStringDecay &right)
int G4int
Definition: G4Types.hh:78
G4KineticTrackVector * FragmentString(const G4ExcitedString &theString)
G4GLOB_DLL std::ostream G4cout
bool G4bool
Definition: G4Types.hh:79
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:298
#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
G4GLOB_DLL std::ostream G4cerr
CLHEP::HepLorentzVector G4LorentzVector