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