| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // Historic fragment from M.Komogorov; clean-u 26 // Historic fragment from M.Komogorov; clean-up still necessary @@@
27 27
28 #include "G4ExcitedStringDecay.hh" 28 #include "G4ExcitedStringDecay.hh"
29 #include "G4SystemOfUnits.hh" 29 #include "G4SystemOfUnits.hh"
30 #include "G4KineticTrack.hh" 30 #include "G4KineticTrack.hh"
31 #include "G4LundStringFragmentation.hh" 31 #include "G4LundStringFragmentation.hh"
32 #include "G4HadronicInteractionRegistry.hh" 32 #include "G4HadronicInteractionRegistry.hh"
33 33
34 #include "G4SampleResonance.hh" 34 #include "G4SampleResonance.hh"
35 35
36 //#define debug_G4ExcitedStringDecay 36 //#define debug_G4ExcitedStringDecay
37 //#define debug_G4ExcitedStringCorr 37 //#define debug_G4ExcitedStringCorr
38 38
39 G4ExcitedStringDecay::G4ExcitedStringDecay(G4V 39 G4ExcitedStringDecay::G4ExcitedStringDecay(G4VLongitudinalStringDecay* ptr)
40 : G4VStringFragmentation(), theStringDecay(p 40 : G4VStringFragmentation(), theStringDecay(ptr)
41 { 41 {
42 if(!ptr) { 42 if(!ptr) {
43 G4HadronicInteraction* p = 43 G4HadronicInteraction* p =
44 G4HadronicInteractionRegistry::Instance( 44 G4HadronicInteractionRegistry::Instance()->FindModel("LundStringFragmentation");
45 theStringDecay = static_cast<G4VLongitudin 45 theStringDecay = static_cast<G4VLongitudinalStringDecay*>(p);
46 if(!theStringDecay) { theStringDecay = new 46 if(!theStringDecay) { theStringDecay = new G4LundStringFragmentation(); }
47 } 47 }
48 SetModelName(theStringDecay->GetModelName()) 48 SetModelName(theStringDecay->GetModelName());
49 } 49 }
50 50
51 G4ExcitedStringDecay::~G4ExcitedStringDecay() 51 G4ExcitedStringDecay::~G4ExcitedStringDecay()
52 {} 52 {}
53 53
54 G4KineticTrackVector *G4ExcitedStringDecay::Fr 54 G4KineticTrackVector *G4ExcitedStringDecay::FragmentString(const G4ExcitedString &theString)
55 { 55 {
56 return theStringDecay->FragmentString(theStr 56 return theStringDecay->FragmentString(theString);
57 } 57 }
58 58
59 G4KineticTrackVector *G4ExcitedStringDecay::Fr 59 G4KineticTrackVector *G4ExcitedStringDecay::FragmentStrings(const G4ExcitedStringVector * theStrings)
60 { 60 {
61 G4LorentzVector KTsum(0.,0.,0.,0.); 61 G4LorentzVector KTsum(0.,0.,0.,0.);
62 62
63 #ifdef debug_G4ExcitedStringDecay 63 #ifdef debug_G4ExcitedStringDecay
64 G4cout<<G4endl; 64 G4cout<<G4endl;
65 G4cout<<"--------------------------- G4Excit 65 G4cout<<"--------------------------- G4ExcitedStringDecay ----------------------"<<G4endl;
66 G4cout<<"Hadronization of Excited Strings: t 66 G4cout<<"Hadronization of Excited Strings: theStrings->size() "<<theStrings->size()<<G4endl;
67 #endif 67 #endif
68 68
69 for ( unsigned int astring=0; astring < theS 69 for ( unsigned int astring=0; astring < theStrings->size(); astring++)
70 // for ( unsigned int astring=0; astring < 1 70 // for ( unsigned int astring=0; astring < 1; astring++)
71 { 71 {
72 // G4cout<<"theStrings->operator[](astring 72 // G4cout<<"theStrings->operator[](astring)->IsExcited() "<<" "<<astring<<" "<<theStrings->operator[](astring)->IsExcited()<<G4endl;
73 if ( theStrings->operator[](astring)->IsEx 73 if ( theStrings->operator[](astring)->IsExcited() )
74 {KTsum+= theStrings->operator[](astri 74 {KTsum+= theStrings->operator[](astring)->Get4Momentum();}
75 else {KTsum+=theStrings->operator[](astrin 75 else {KTsum+=theStrings->operator[](astring)->GetKineticTrack()->Get4Momentum();}
76 } 76 }
77 77
78 G4LorentzRotation toCms( -1 * KTsum.boostVec 78 G4LorentzRotation toCms( -1 * KTsum.boostVector() );
79 G4LorentzRotation toLab(toCms.inverse()); 79 G4LorentzRotation toLab(toCms.inverse());
80 G4LorentzVector Ptmp; 80 G4LorentzVector Ptmp;
81 KTsum=G4LorentzVector(0.,0.,0.,0.); 81 KTsum=G4LorentzVector(0.,0.,0.,0.);
82 82
83 for ( unsigned int astring=0; astring < theS 83 for ( unsigned int astring=0; astring < theStrings->size(); astring++)
84 // for ( unsigned int astring=0; astring < 1 84 // for ( unsigned int astring=0; astring < 1; astring++)
85 { 85 {
86 if ( theStrings->operator[](astring)->IsEx 86 if ( theStrings->operator[](astring)->IsExcited() )
87 { 87 {
88 Ptmp=toCms * theStrings->operator[](astr 88 Ptmp=toCms * theStrings->operator[](astring)->GetLeftParton()->Get4Momentum();
89 theStrings->operator[](astring)->GetLeft 89 theStrings->operator[](astring)->GetLeftParton()->Set4Momentum(Ptmp);
90 90
91 Ptmp=toCms * theStrings->operator[](astr 91 Ptmp=toCms * theStrings->operator[](astring)->GetRightParton()->Get4Momentum();
92 theStrings->operator[](astring)->GetRigh 92 theStrings->operator[](astring)->GetRightParton()->Set4Momentum(Ptmp);
93 93
94 KTsum+= theStrings->operator[](astring)- 94 KTsum+= theStrings->operator[](astring)->Get4Momentum();
95 } 95 }
96 else 96 else
97 { 97 {
98 Ptmp=toCms * theStrings->operator[](astr 98 Ptmp=toCms * theStrings->operator[](astring)->GetKineticTrack()->Get4Momentum();
99 theStrings->operator[](astring)->GetKine 99 theStrings->operator[](astring)->GetKineticTrack()->Set4Momentum(Ptmp);
100 KTsum+= theStrings->operator[](astring)- 100 KTsum+= theStrings->operator[](astring)->GetKineticTrack()->Get4Momentum();
101 } 101 }
102 } 102 }
103 103
104 G4SampleResonance BrW; 104 G4SampleResonance BrW;
105 const G4ParticleDefinition* TrackDefinition= 105 const G4ParticleDefinition* TrackDefinition=0;
106 106
107 G4KineticTrackVector * theResult = new G4Kin 107 G4KineticTrackVector * theResult = new G4KineticTrackVector;
108 G4int attempts(0); 108 G4int attempts(0);
109 G4bool success=false; 109 G4bool success=false;
110 G4bool NeedEnergyCorrector=false; 110 G4bool NeedEnergyCorrector=false;
111 do { 111 do {
112 #ifdef debug_G4ExcitedStringDecay 112 #ifdef debug_G4ExcitedStringDecay
113 G4cout<<"New try No "<<attempts<<" to 113 G4cout<<"New try No "<<attempts<<" to hadronize strings"<<G4endl;
114 #endif 114 #endif
115 115
116 std::for_each(theResult->begin() , theResult 116 std::for_each(theResult->begin() , theResult->end() , DeleteKineticTrack());
117 theResult->clear(); 117 theResult->clear();
118 118
119 attempts++; 119 attempts++;
120 120
121 G4LorentzVector KTsecondaries(0.,0.,0.,0.); 121 G4LorentzVector KTsecondaries(0.,0.,0.,0.);
122 NeedEnergyCorrector=false; 122 NeedEnergyCorrector=false;
123 123
124 for ( unsigned int astring=0; astring < theS 124 for ( unsigned int astring=0; astring < theStrings->size(); astring++)
125 // for ( unsigned int astring=0; astri 125 // for ( unsigned int astring=0; astring < 1; astring++) // Proj. Last Str. Decay for FTF
126 // for ( unsigned int astring=1; astri 126 // for ( unsigned int astring=1; astring < 2; astring++) // Proj. Last Str. Decay for QGS
127 // for ( unsigned int astring=0; astri 127 // for ( unsigned int astring=0; astring < 1; astring++) // For testing purposes
128 { 128 {
129 #ifdef debug_G4ExcitedStringDecay 129 #ifdef debug_G4ExcitedStringDecay
130 G4cout<<"String No "<<astring+1<<" E 130 G4cout<<"String No "<<astring+1<<" Excited? "<<theStrings->operator[](astring)->IsExcited()<<G4endl;
131 131
132 G4cout<<"String No "<<astring+1<<" 4 132 G4cout<<"String No "<<astring+1<<" 4Momentum "<<theStrings->operator[](astring)->Get4Momentum()
133 <<" "<<theStrings->operator[](astrin 133 <<" "<<theStrings->operator[](astring)->Get4Momentum().mag()<<G4endl;
134 #endif 134 #endif
135 135
136 G4KineticTrackVector * generatedKine 136 G4KineticTrackVector * generatedKineticTracks = NULL;
137 if ( theStrings->operator[](astring)->IsEx 137 if ( theStrings->operator[](astring)->IsExcited() )
138 { 138 {
139 #ifdef debug_G4ExcitedStringDecay 139 #ifdef debug_G4ExcitedStringDecay
140 G4cout<<"Fragment String with par 140 G4cout<<"Fragment String with partons: "
141 <<theStrings->operator[](as 141 <<theStrings->operator[](astring)->GetLeftParton()->GetPDGcode() <<" "
142 <<theStrings->operator[](as 142 <<theStrings->operator[](astring)->GetRightParton()->GetPDGcode()<<" "
143 <<"Direction "<<theStrings- 143 <<"Direction "<<theStrings->operator[](astring)->GetDirection()<<G4endl;
144 #endif 144 #endif
145 generatedKineticTracks=FragmentString 145 generatedKineticTracks=FragmentString(*theStrings->operator[](astring));
146 #ifdef debug_G4ExcitedStringDecay 146 #ifdef debug_G4ExcitedStringDecay
147 G4cout<<"(G4ExcitedStringDecay) N 147 G4cout<<"(G4ExcitedStringDecay) Number of produced hadrons = "
148 <<generatedKineticTracks->s 148 <<generatedKineticTracks->size()<<G4endl;
149 #endif 149 #endif
150 } else { 150 } else {
151 #ifdef debug_G4ExcitedStringDecay 151 #ifdef debug_G4ExcitedStringDecay
152 G4cout<<" GetTrack from the Str 152 G4cout<<" GetTrack from the String"<<G4endl;
153 #endif 153 #endif
154 G4LorentzVector Mom=theStrings->o 154 G4LorentzVector Mom=theStrings->operator[](astring)->GetKineticTrack()->Get4Momentum();
155 G4KineticTrack * aTrack= new G4Ki 155 G4KineticTrack * aTrack= new G4KineticTrack(
156 theStrings->operat 156 theStrings->operator[](astring)->GetKineticTrack()->GetDefinition(),
157 theStrings->operat 157 theStrings->operator[](astring)->GetKineticTrack()->GetFormationTime(),
158 G4ThreeVector(0), 158 G4ThreeVector(0), Mom);
159 159
160 aTrack->SetPosition(theStrings->o 160 aTrack->SetPosition(theStrings->operator[](astring)->GetKineticTrack()->GetPosition());
161 161
162 #ifdef debug_G4ExcitedStringDecay 162 #ifdef debug_G4ExcitedStringDecay
163 G4cout<<" A particle stored in 163 G4cout<<" A particle stored in the track is "<<aTrack->GetDefinition()->GetParticleName()<<G4endl;
164 #endif 164 #endif
165 165
166 generatedKineticTracks = new G4KineticT 166 generatedKineticTracks = new G4KineticTrackVector;
167 generatedKineticTracks->push_back(aTrac 167 generatedKineticTracks->push_back(aTrack);
168 } 168 }
169 169
170 if (generatedKineticTracks == nullptr || g 170 if (generatedKineticTracks == nullptr || generatedKineticTracks->size() == 0)
171 { 171 {
172 // G4cerr << "G4VPartonStringMode 172 // G4cerr << "G4VPartonStringModel:No KineticTracks produced" << G4endl;
173 delete generatedKineticTracks; << 173 // continue;
174 generatedKineticTracks = nullptr; << 174 success=false; NeedEnergyCorrector=false; break;
175 continue; <<
176 } 175 }
177 176
178 G4LorentzVector KTsum1(0.,0.,0.,0.); 177 G4LorentzVector KTsum1(0.,0.,0.,0.);
179 for ( unsigned int aTrack=0; aTrack< 178 for ( unsigned int aTrack=0; aTrack<generatedKineticTracks->size();aTrack++)
180 { 179 {
181 #ifdef debug_G4ExcitedStringDecay 180 #ifdef debug_G4ExcitedStringDecay
182 G4cout<<"Prod part No. "<<aTrack+ 181 G4cout<<"Prod part No. "<<aTrack+1<<" "
183 <<(*generatedKineticTracks) 182 <<(*generatedKineticTracks)[aTrack]->GetDefinition()->GetParticleName()<<" "
184 <<(*generatedKineticTracks) 183 <<(*generatedKineticTracks)[aTrack]->Get4Momentum()
185 <<(*generatedKineticTracks) 184 <<(*generatedKineticTracks)[aTrack]->Get4Momentum().mag()<<G4endl;
186 #endif 185 #endif
187 // --------------- Sampling mass 186 // --------------- Sampling mass of unstable hadronic resonances ----------------
188 TrackDefinition = (*generatedKine 187 TrackDefinition = (*generatedKineticTracks)[aTrack]->GetDefinition();
189 188
190 if (TrackDefinition->IsShortLived()) 189 if (TrackDefinition->IsShortLived())
191 { 190 {
192 G4double NewTrackMass = 191 G4double NewTrackMass =
193 BrW.SampleMass( TrackDefiniti 192 BrW.SampleMass( TrackDefinition->GetPDGMass(), TrackDefinition->GetPDGWidth(),
194 BrW.GetMinimu 193 BrW.GetMinimumMass( TrackDefinition ) + 10.0*MeV,
195 TrackDefiniti 194 TrackDefinition->GetPDGMass() + 5.0*TrackDefinition->GetPDGWidth() );
196 G4LorentzVector Tmp=G4LorentzVe 195 G4LorentzVector Tmp=G4LorentzVector((*generatedKineticTracks)[aTrack]->Get4Momentum());
197 Tmp.setE(std::sqrt(sqr(NewTrack 196 Tmp.setE(std::sqrt(sqr(NewTrackMass) + Tmp.vect().mag2()));
198 197
199 (*generatedKineticTracks)[aTrac 198 (*generatedKineticTracks)[aTrack]->Set4Momentum(Tmp);
200 199
201 #ifdef debug_G4ExcitedStringDec 200 #ifdef debug_G4ExcitedStringDecay
202 G4cout<<"Resonance *** "<<aTrac 201 G4cout<<"Resonance *** "<<aTrack+1<<" "
203 <<(*generatedKineticTrack 202 <<(*generatedKineticTracks)[aTrack]->GetDefinition()->GetParticleName()<<" "
204 <<(*generatedKineticTrack 203 <<(*generatedKineticTracks)[aTrack]->Get4Momentum()
205 <<(*generatedKineticTrack 204 <<(*generatedKineticTracks)[aTrack]->Get4Momentum().mag()<<G4endl;
206 #endif 205 #endif
207 } 206 }
208 //------------------------------- 207 //-------------------------------------------------------------------------------
209 208
210 theResult->push_back(generatedKin 209 theResult->push_back(generatedKineticTracks->operator[](aTrack));
211 KTsum1+= (*generatedKineticTracks 210 KTsum1+= (*generatedKineticTracks)[aTrack]->Get4Momentum();
212 } 211 }
213 KTsecondaries+=KTsum1; 212 KTsecondaries+=KTsum1;
214 213
215 #ifdef debug_G4ExcitedStringDecay 214 #ifdef debug_G4ExcitedStringDecay
216 G4cout << "String secondaries(" <<ge 215 G4cout << "String secondaries(" <<generatedKineticTracks->size()<< ")"<<G4endl
217 <<"Init string momentum: "< 216 <<"Init string momentum: "<< theStrings->operator[](astring)->Get4Momentum()<<G4endl
218 <<"Final hadrons momentum: "< 217 <<"Final hadrons momentum: "<< KTsum1 << G4endl;
219 #endif 218 #endif
220 219
221 if ( KTsum1.e() > 0 && 220 if ( KTsum1.e() > 0 &&
222 std::abs((KTsum1.e()-theString 221 std::abs((KTsum1.e()-theStrings->operator[](astring)->Get4Momentum().e()) / KTsum1.e()) > perMillion )
223 { 222 {
224 NeedEnergyCorrector=true; 223 NeedEnergyCorrector=true;
225 } 224 }
226 225
227 #ifdef debug_G4ExcitedStringDecay 226 #ifdef debug_G4ExcitedStringDecay
228 G4cout<<"NeedEnergyCorrection yes/no 227 G4cout<<"NeedEnergyCorrection yes/no "<<NeedEnergyCorrector<<G4endl;
229 #endif 228 #endif
230 229
231 // clean up 230 // clean up
232 delete generatedKineticTracks; 231 delete generatedKineticTracks;
233 success=true; 232 success=true;
234 } 233 }
235 234
236 if ( NeedEnergyCorrector ) success=EnergyAnd 235 if ( NeedEnergyCorrector ) success=EnergyAndMomentumCorrector(theResult, KTsum);
237 } while (!success && (attempts < 100)); /* 236 } while (!success && (attempts < 100)); /* Loop checking, 07.08.2015, A.Ribon */
238 237
239 for ( unsigned int aTrack=0; aTrack<theResul 238 for ( unsigned int aTrack=0; aTrack<theResult->size();aTrack++)
240 { 239 {
241 Ptmp=(*theResult)[aTrack]->Get4Momentum(); 240 Ptmp=(*theResult)[aTrack]->Get4Momentum();
242 Ptmp.transform( toLab); 241 Ptmp.transform( toLab);
243 (*theResult)[aTrack]->Set4Momentum(Ptmp); 242 (*theResult)[aTrack]->Set4Momentum(Ptmp);
244 } 243 }
245 244
246 #ifdef debug_G4ExcitedStringDecay 245 #ifdef debug_G4ExcitedStringDecay
247 G4cout<<"End of the strings fragmentation (G 246 G4cout<<"End of the strings fragmentation (G4ExcitedStringDecay)"<<G4endl;
248 247
249 G4LorentzVector KTsum1(0.,0.,0.,0.); 248 G4LorentzVector KTsum1(0.,0.,0.,0.);
250 249
251 for ( unsigned int aTrack=0; aTrack<theResul 250 for ( unsigned int aTrack=0; aTrack<theResult->size();aTrack++)
252 { 251 {
253 G4cout << " corrected tracks .. " << (*t 252 G4cout << " corrected tracks .. " << (*theResult)[aTrack]->GetDefinition()->GetParticleName()
254 <<" " << (*theResult)[aTrack]->Get4Mome 253 <<" " << (*theResult)[aTrack]->Get4Momentum()
255 <<" " << (*theResult)[aTrack]->Get4Mome 254 <<" " << (*theResult)[aTrack]->Get4Momentum().mag()<< G4endl;
256 KTsum1+= (*theResult)[aTrack]->Get4Momen 255 KTsum1+= (*theResult)[aTrack]->Get4Momentum();
257 } 256 }
258 257
259 G4cout << "Needcorrector/success " << NeedEn 258 G4cout << "Needcorrector/success " << NeedEnergyCorrector << "/" << success
260 << ", Corrected total 4 momentum " < 259 << ", Corrected total 4 momentum " << KTsum1 << G4endl;
261 if ( ! success ) G4cout << "failed to correc 260 if ( ! success ) G4cout << "failed to correct E/p" << G4endl;
262 261
263 G4cout<<"End of the Hadronization (G4Excited 262 G4cout<<"End of the Hadronization (G4ExcitedStringDecay)"<<G4endl;
264 #endif 263 #endif
265 264
266 if (!success) 265 if (!success)
267 { 266 {
268 if (theResult->size() != 0) 267 if (theResult->size() != 0)
269 { 268 {
270 std::for_each(theResult->begin() , theRe 269 std::for_each(theResult->begin() , theResult->end() , DeleteKineticTrack());
271 theResult->clear(); 270 theResult->clear();
272 delete theResult; theResult=0; 271 delete theResult; theResult=0;
273 } 272 }
274 for ( unsigned int astring=0; astring < th 273 for ( unsigned int astring=0; astring < theStrings->size(); astring++)
275 // for ( unsigned int astring=0; astring < 274 // for ( unsigned int astring=0; astring < 1; astring++) // Need more correct. For testing purposes.
276 { 275 {
277 if ( theStrings->operator[](astring)->Is 276 if ( theStrings->operator[](astring)->IsExcited() )
278 { 277 {
279 Ptmp=theStrings->operator[](astring)-> 278 Ptmp=theStrings->operator[](astring)->GetLeftParton()->Get4Momentum();
280 Ptmp.transform( toLab); 279 Ptmp.transform( toLab);
281 theStrings->operator[](astring)->GetLe 280 theStrings->operator[](astring)->GetLeftParton()->Set4Momentum(Ptmp);
282 281
283 Ptmp=theStrings->operator[](astring)-> 282 Ptmp=theStrings->operator[](astring)->GetRightParton()->Get4Momentum();
284 Ptmp.transform( toLab); 283 Ptmp.transform( toLab);
285 theStrings->operator[](astring)->GetRi 284 theStrings->operator[](astring)->GetRightParton()->Set4Momentum(Ptmp);
286 } 285 }
287 else 286 else
288 { 287 {
289 Ptmp=theStrings->operator[](astring)-> 288 Ptmp=theStrings->operator[](astring)->GetKineticTrack()->Get4Momentum();
290 Ptmp.transform( toLab); 289 Ptmp.transform( toLab);
291 theStrings->operator[](astring)->GetKi 290 theStrings->operator[](astring)->GetKineticTrack()->Set4Momentum(Ptmp);
292 } 291 }
293 } 292 }
294 } 293 }
295 return theResult; 294 return theResult;
296 } 295 }
297 296
298 297
299 G4bool G4ExcitedStringDecay:: 298 G4bool G4ExcitedStringDecay::
300 EnergyAndMomentumCorrector(G4KineticTrackVecto 299 EnergyAndMomentumCorrector(G4KineticTrackVector* Output, G4LorentzVector& TotalCollisionMom)
301 { 300 {
302 const int nAttemptScale = 500; 301 const int nAttemptScale = 500;
303 const double ErrLimit = 1.E-5; 302 const double ErrLimit = 1.E-5;
304 if (Output->empty()) return TRUE; 303 if (Output->empty()) return TRUE;
305 G4LorentzVector SumMom; 304 G4LorentzVector SumMom;
306 G4double SumMass = 0; 305 G4double SumMass = 0;
307 G4double TotalCollisionMass = Total 306 G4double TotalCollisionMass = TotalCollisionMom.m();
308 307
309 std::vector<G4double> HadronMass; G4double 308 std::vector<G4double> HadronMass; G4double HadronM(0.);
310 309
311 #ifdef debug_G4ExcitedStringCorr 310 #ifdef debug_G4ExcitedStringCorr
312 G4cout<<G4endl<<"EnergyAndMomentumCorrecto 311 G4cout<<G4endl<<"EnergyAndMomentumCorrector. Number of particles: "<<Output->size()<<G4endl;
313 #endif 312 #endif
314 // Calculate sum hadron 4-momenta and summ 313 // Calculate sum hadron 4-momenta and summing hadron mass
315 unsigned int cHadron; 314 unsigned int cHadron;
316 for (cHadron = 0; cHadron < Output->size() 315 for (cHadron = 0; cHadron < Output->size(); cHadron++)
317 { 316 {
318 SumMom += Output->operator[](cHadron) 317 SumMom += Output->operator[](cHadron)->Get4Momentum();
319 HadronM=Output->operator[](cHadron)->G 318 HadronM=Output->operator[](cHadron)->Get4Momentum().mag(); HadronMass.push_back(HadronM);
320 SumMass += Output->operator[](cHadron) 319 SumMass += Output->operator[](cHadron)->Get4Momentum().mag(); //GetDefinition()->GetPDGMass();
321 } 320 }
322 321
323 #ifdef debug_G4ExcitedStringCorr 322 #ifdef debug_G4ExcitedStringCorr
324 G4cout<<"Sum part mom "<<SumMom<<" "<<SumM 323 G4cout<<"Sum part mom "<<SumMom<<" "<<SumMom.mag()<<G4endl
325 <<"Sum str mom "<<TotalCollisionMom 324 <<"Sum str mom "<<TotalCollisionMom<<" "<<TotalCollisionMom.mag()<<G4endl;
326 G4cout<<"SumMass TotalCollisionMass "<<Sum 325 G4cout<<"SumMass TotalCollisionMass "<<SumMass<<" "<<TotalCollisionMass<<G4endl;
327 #endif 326 #endif
328 327
329 // Cannot correct a single particle 328 // Cannot correct a single particle
330 if (Output->size() < 2) return FALSE; 329 if (Output->size() < 2) return FALSE;
331 330
332 if (SumMass > TotalCollisionMass) return F 331 if (SumMass > TotalCollisionMass) return FALSE;
333 SumMass = SumMom.m2(); 332 SumMass = SumMom.m2();
334 if (SumMass < 0) return FALSE; 333 if (SumMass < 0) return FALSE;
335 SumMass = std::sqrt(SumMass); 334 SumMass = std::sqrt(SumMass);
336 335
337 // Compute c.m.s. hadron velocity and boos 336 // Compute c.m.s. hadron velocity and boost KTV to hadron c.m.s.
338 // G4ThreeVector Beta = -SumMom.boostVecto 337 // G4ThreeVector Beta = -SumMom.boostVector();
339 G4ThreeVector Beta = -TotalCollisionMom.bo 338 G4ThreeVector Beta = -TotalCollisionMom.boostVector();
340 Output->Boost(Beta); 339 Output->Boost(Beta);
341 340
342 // Scale total c.m.s. hadron energy (hadro 341 // Scale total c.m.s. hadron energy (hadron system mass).
343 // It should be equal interaction mass 342 // It should be equal interaction mass
344 G4double Scale = 1; 343 G4double Scale = 1;
345 G4int cAttempt = 0; 344 G4int cAttempt = 0;
346 G4double Sum = 0; 345 G4double Sum = 0;
347 G4bool success = false; 346 G4bool success = false;
348 for (cAttempt = 0; cAttempt < nAttemptScal 347 for (cAttempt = 0; cAttempt < nAttemptScale; cAttempt++)
349 { 348 {
350 Sum = 0; 349 Sum = 0;
351 for (cHadron = 0; cHadron < Output->size 350 for (cHadron = 0; cHadron < Output->size(); cHadron++)
352 { 351 {
353 HadronM = HadronMass.at(cHadron); 352 HadronM = HadronMass.at(cHadron);
354 G4LorentzVector HadronMom = Output->op 353 G4LorentzVector HadronMom = Output->operator[](cHadron)->Get4Momentum();
355 HadronMom.setVect(Scale*HadronMom.vect 354 HadronMom.setVect(Scale*HadronMom.vect());
356 G4double E = std::sqrt(HadronMom.vect( 355 G4double E = std::sqrt(HadronMom.vect().mag2() + sqr(HadronM));
357 356 //sqr(Output->operator[](cHadron)->GetDefinition()->GetPDGMass()));
358 HadronMom.setE(E); 357 HadronMom.setE(E);
359 Output->operator[](cHadron)->Set4Momen 358 Output->operator[](cHadron)->Set4Momentum(HadronMom);
360 Sum += E; 359 Sum += E;
361 } 360 }
362 Scale = TotalCollisionMass/Sum; 361 Scale = TotalCollisionMass/Sum;
363 #ifdef debug_G4ExcitedStringCorr 362 #ifdef debug_G4ExcitedStringCorr
364 G4cout << "Scale-1=" << Scale -1 363 G4cout << "Scale-1=" << Scale -1
365 << ", TotalCollisionMass=" << To 364 << ", TotalCollisionMass=" << TotalCollisionMass
366 << ", Sum=" << Sum 365 << ", Sum=" << Sum
367 << G4endl; 366 << G4endl;
368 #endif 367 #endif
369 if (std::fabs(Scale - 1) <= ErrLimit) 368 if (std::fabs(Scale - 1) <= ErrLimit)
370 { 369 {
371 success = true; 370 success = true;
372 break; 371 break;
373 } 372 }
374 } 373 }
375 374
376 #ifdef debug_G4ExcitedStringCorr 375 #ifdef debug_G4ExcitedStringCorr
377 if (!success) 376 if (!success)
378 { 377 {
379 G4cout << "G4ExcitedStringDecay::EnergyA 378 G4cout << "G4ExcitedStringDecay::EnergyAndMomentumCorrector - Warning"<<G4endl;
380 G4cout << " Scale not unity at end of 379 G4cout << " Scale not unity at end of iteration loop: "<<TotalCollisionMass<<" "<<Sum<<" "<<Scale<<G4endl;
381 G4cout << " Number of secondaries: " < 380 G4cout << " Number of secondaries: " << Output->size() << G4endl;
382 G4cout << " Wanted total energy: " << 381 G4cout << " Wanted total energy: " << TotalCollisionMom.e() << G4endl;
383 G4cout << " Increase number of attempt 382 G4cout << " Increase number of attempts or increase ERRLIMIT"<<G4endl;
384 // throw G4HadronicException(__FILE__, _ 383 // throw G4HadronicException(__FILE__, __LINE__, "G4ExcitedStringDecay failed to correct...");
385 } 384 }
386 #endif 385 #endif
387 386
388 // Compute c.m.s. interaction velocity and 387 // Compute c.m.s. interaction velocity and KTV back boost
389 Beta = TotalCollisionMom.boostVector(); 388 Beta = TotalCollisionMom.boostVector();
390 Output->Boost(Beta); 389 Output->Boost(Beta);
391 390
392 return success; 391 return success;
393 } 392 }
394 393
395 394