| 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 // 26 //
>> 27 // $Id: G4VPartonStringModel.cc 67999 2013-03-13 11:14:32Z gcosmo $
27 // 28 //
28 //// ----------------------------------------- 29 //// ------------------------------------------------------------
29 // GEANT 4 class implementation file 30 // GEANT 4 class implementation file
30 // 31 //
31 // ---------------- G4VPartonStringModel 32 // ---------------- G4VPartonStringModel ----------------
32 // by Gunter Folger, May 1998. 33 // by Gunter Folger, May 1998.
33 // abstract class for all Parton String M 34 // abstract class for all Parton String Models
34 // ------------------------------------------- 35 // ------------------------------------------------------------
35 // debug switch 36 // debug switch
36 //#define debug_PartonStringModel 37 //#define debug_PartonStringModel
37 //#define debug_heavyHadrons << 38
38 // ------------------------------------------- <<
39 39
40 #include "G4VPartonStringModel.hh" 40 #include "G4VPartonStringModel.hh"
41 #include "G4ios.hh" 41 #include "G4ios.hh"
42 #include "G4ShortLivedConstructor.hh" 42 #include "G4ShortLivedConstructor.hh"
43 43
44 #include "G4ParticleTable.hh" 44 #include "G4ParticleTable.hh"
45 #include "G4IonTable.hh" 45 #include "G4IonTable.hh"
46 46
>> 47 //#define debug_PartonStringModel
>> 48
47 G4VPartonStringModel::G4VPartonStringModel(con 49 G4VPartonStringModel::G4VPartonStringModel(const G4String& modelName)
48 : G4VHighEnergyGenerator(modelName), 50 : G4VHighEnergyGenerator(modelName),
49 stringFragmentationModel(nullptr) << 51 stringFragmentationModel(0),
>> 52 theThis(0)
50 { 53 {
51 // Make shure Shotrylived particles are con << 54 // Make shure Shotrylived partyicles are constructed.
52 // VI: should not instantiate particles by << 55 G4ShortLivedConstructor ShortLived;
53 /* << 56 ShortLived.ConstructParticle();
54 G4ShortLivedConstructor ShortLived; <<
55 ShortLived.ConstructParticle(); <<
56 */ <<
57 } 57 }
58 58
59 G4VPartonStringModel::~G4VPartonStringModel() 59 G4VPartonStringModel::~G4VPartonStringModel()
60 { 60 {
61 } 61 }
62 62
63 G4KineticTrackVector * G4VPartonStringModel::S 63 G4KineticTrackVector * G4VPartonStringModel::Scatter(const G4Nucleus &theNucleus,
64 64 const G4DynamicParticle &aPrimary)
65 { 65 {
66 G4ExcitedStringVector * strings = nullptr; << 66 G4ExcitedStringVector * strings = NULL;
>> 67
67 G4DynamicParticle thePrimary=aPrimary; 68 G4DynamicParticle thePrimary=aPrimary;
68 G4LorentzVector SumStringMom(0.,0.,0.,0.); <<
69 G4KineticTrackVector * theResult = 0; <<
70 G4Nucleon * theNuclNucleon(nullptr); <<
71 69
72 #ifdef debug_PartonStringModel << 70 #ifdef debug_PartonStringModel
73 G4cout<<G4endl; << 71 G4cout<<"Parton-String model is runnung ------------"<<G4endl;
74 G4cout<<"-----------------------Parton-Strin <<
75 G4cout<<"Projectile Name Mass "<<thePrimary. 72 G4cout<<"Projectile Name Mass "<<thePrimary.GetDefinition()->GetParticleName()<<" "
76 <<thePrimary. << 73 <<thePrimary.GetMass()<<G4endl;
77 G4cout<<" Momentum "<<thePrimary. 74 G4cout<<" Momentum "<<thePrimary.Get4Momentum()<<G4endl;
78 G4cout<<"Target nucleus A Z "<<theNucleus. 75 G4cout<<"Target nucleus A Z "<<theNucleus.GetA_asInt()<<" "
79 <<theNucleus. 76 <<theNucleus.GetZ_asInt()<<G4endl<<G4endl;
80 G4int Bsum=thePrimary.GetDefinition()->GetBa << 77 #endif
81 G4int Qsum=thePrimary.GetDefinition()->GetPD << 78
82 G4cout<<"Initial baryon number "<<Bsum<<G4en <<
83 G4cout<<"Initial charge "<<Qsum<<G4en <<
84 G4cout<<"-------------- Parton-String model: <<
85 Bsum -= theNucleus.GetA_asInt(); Qsum -= th <<
86 if(GetProjectileNucleus()) { <<
87 Bsum -= thePrimary.GetDefinition()->GetBar <<
88 Qsum -= thePrimary.GetDefinition()->GetPDG <<
89 } <<
90 G4int QsumSec(0), BsumSec(0); <<
91 G4LorentzVector SumPsecondr(0.,0.,0.,0.); <<
92 #endif <<
93 <<
94 G4LorentzRotation toZ; 79 G4LorentzRotation toZ;
95 G4LorentzVector Ptmp=thePrimary.Get4Momentum 80 G4LorentzVector Ptmp=thePrimary.Get4Momentum();
96 toZ.rotateZ(-1*Ptmp.phi()); 81 toZ.rotateZ(-1*Ptmp.phi());
97 toZ.rotateY(-1*Ptmp.theta()); 82 toZ.rotateY(-1*Ptmp.theta());
98 thePrimary.Set4Momentum(toZ*Ptmp); 83 thePrimary.Set4Momentum(toZ*Ptmp);
99 G4LorentzRotation toLab(toZ.inverse()); 84 G4LorentzRotation toLab(toZ.inverse());
100 85
101 G4bool Success=true; << 86 G4int attempts = 0, maxAttempts=20;
102 G4int attempts = 0, maxAttempts=1000; << 87 while ( strings == NULL )
103 do <<
104 { 88 {
105 if (attempts++ > maxAttempts ) << 89 if (attempts++ > maxAttempts )
106 { << 90 {
107 Init(theNucleus,thePrimary); // To put << 91 throw G4HadronicException(__FILE__, __LINE__,
108 / << 92 "G4VPartonStringModel::Scatter(): fails to generate strings");
109 G4V3DNucleus * ResNucleus = GetWoundedNu << 93 }
110 theNuclNucleon = ResNucleus->StartLoop() << 94 theThis->Init(theNucleus,thePrimary);
111 while( theNuclNucleon ) <<
112 { <<
113 if(theNuclNucleon->AreYouHit()) theNuc <<
114 theNuclNucleon = ResNucleus->GetNextNu <<
115 } <<
116 <<
117 G4V3DNucleus * ProjResNucleus = GetProje <<
118 if(ProjResNucleus != 0) <<
119 { <<
120 theNuclNucleon = ProjResNucleus->Start <<
121 while( theNuclNucleon ) <<
122 { <<
123 if(theNuclNucleon->AreYouHit()) theN <<
124 theNuclNucleon = ProjResNucleus->Get <<
125 } <<
126 } <<
127 <<
128 G4ExceptionDescription ed; <<
129 ed << "Projectile Name Mass " <<thePrima <<
130 << " " << thePrimary.GetMass()<< G4en <<
131 ed << "          Momentum " <<
132 ed << "Target nucleus  A Z " << theNu <<
133 << theNu <<
134 ed << "Initial states of projectile and <<
135 G4Exception( "G4VPartonStringModel::Scat <<
136 "HAD_PARTON_STRING_001", Ju <<
137 <<
138 G4ThreeVector Position(0.,0.,2*ResNucleu <<
139 G4KineticTrack* Hadron = new G4KineticTr <<
140 <<
141 if(theResult == nullptr) theResult = new <<
142 theResult->push_back(Hadron); <<
143 return theResult; <<
144 } <<
145 <<
146 Success=true; <<
147 <<
148 Init(theNucleus,thePrimary); <<
149 <<
150 strings = GetStrings(); <<
151 <<
152 if (strings->empty()) { Success=false; con <<
153 <<
154 // G4double stringEnergy(0); <<
155 SumStringMom=G4LorentzVector(0.,0.,0.,0.); <<
156 <<
157 #ifdef debug_PartonStringModel <<
158 G4cout<<"------------ Parton-String model: <<
159 #endif <<
160 <<
161 #ifdef debug_heavyHadrons <<
162 // Check charm and bottom numbers of the p <<
163 G4int count_charm_projectile = thePrimary <<
164 thePrimary <<
165 G4int count_bottom_projectile = thePrimary <<
166 thePrimary <<
167 G4int count_charm_strings = 0, count_botto <<
168 G4int count_charm_hadrons = 0, count_botto <<
169 #endif <<
170 <<
171 for ( unsigned int astring=0; astring < st <<
172 { <<
173 // rotate string to lab frame, models <<
174 if((*strings)[astring]->IsExcited()) <<
175 { <<
176 // stringEnergy += (*strings)[astring] <<
177 // stringEnergy += (*strings)[astring] <<
178 (*strings)[astring]->LorentzRotate(toL <<
179 SumStringMom+=(*strings)[astring]->Get <<
180 #ifdef debug_PartonStringModel <<
181 G4cout<<"String No "<<astring+1<<" "<< <<
182 <<(*strings)[astri <<
183 <<" Partons "<<(*strings)[astr <<
184 <<" "<<(*strings)[astri <<
185 #endif <<
186 <<
187 #ifdef debug_heavyHadrons <<
188 G4int left_charm = (*strings)[astring]- <<
189 G4int left_anticharm = (*strings)[astring]- <<
190 G4int right_charm = (*strings)[astring]- <<
191 G4int right_anticharm = (*strings)[astring]- <<
192 G4int left_bottom = (*strings)[astring] <<
193 G4int left_antibottom = (*strings)[astring] <<
194 G4int right_bottom = (*strings)[astring] <<
195 G4int right_antibottom = (*strings)[astring] <<
196 if ( left_charm != 0 || left_anticharm ! <<
197 left_bottom != 0 || left_antibottom ! <<
198 count_charm_strings += left_charm - left <<
199 count_bottom_strings += left_bottom - left <<
200 G4cout << "G4VPartonStringModel::Scatter : <<
201 << (*strings)[astring]->GetLeftParton()-> <<
202 << (*strings)[astring]->GetRightParton()- <<
203 } <<
204 #endif <<
205 } <<
206 else <<
207 { <<
208 // stringEnergy += (*strings)[astring] <<
209 (*strings)[astring]->LorentzRotate(toL <<
210 SumStringMom+=(*strings)[astring]->Get <<
211 #ifdef debug_PartonStringModel <<
212 G4cout<<"A track No "<<astring+1<<" " <<
213 <<(*strings)[astring]->GetKineti <<
214 <<(*strings)[astring]->GetKineti <<
215 <<(*strings)[astring]->GetKineti <<
216 #endif <<
217 <<
218 #ifdef debug_heavyHadrons <<
219 G4int charm = (*strings)[astring]->GetK <<
220 G4int anticharm = (*strings)[astring]->GetK <<
221 G4int bottom = (*strings)[astring]->GetK <<
222 G4int antibottom = (*strings)[astring]->GetK <<
223 if ( charm != 0 || anticharm != 0 | <<
224 count_charm_strings += charm - anticharm <<
225 count_bottom_strings += bottom - antibotto <<
226 G4cout << "G4VPartonStringModel::Scatter : <<
227 << (*strings)[astring]->GetKi <<
228 } <<
229 #endif <<
230 } <<
231 } <<
232 <<
233 #ifdef debug_heavyHadrons <<
234 if ( count_charm_projectile != count_charm <<
235 G4cout << "G4VPartonStringModel::Scatter <<
236 << count_charm_projectile << " ; #strin <<
237 } <<
238 if ( count_bottom_projectile != count_bott <<
239 G4cout << "G4VPartonStringModel::Scatter <<
240 << count_bottom_projectile << " ; #stri <<
241 } <<
242 #endif <<
243 <<
244 #ifdef debug_PartonStringModel <<
245 G4cout<<G4endl<<"SumString4Mom "<<SumStrin <<
246 G4LorentzVector TargetResidual4Momentum(0. <<
247 G4LorentzVector ProjectileResidual4Momentu <<
248 G4int hitsT(0), charged_hitsT(0); <<
249 G4int hitsP(0), charged_hitsP(0); <<
250 G4double ExcitationEt(0.), ExcitationEp(0. <<
251 #endif <<
252 <<
253 // We assume that the target nucleus is ne <<
254 // the projectile nucleus can be a light h <<
255 <<
256 G4V3DNucleus * ProjResNucleus = GetProject <<
257 <<
258 G4int numberProtonProjectileResidual( 0 ), <<
259 G4int numberLambdaProjectileResidual( 0 ); <<
260 if(ProjResNucleus != 0) <<
261 { <<
262 theNuclNucleon = ProjResNucleus->StartLo <<
263 G4int numberProtonProjectileHits( 0 ), n <<
264 G4int numberLambdaProjectileHits( 0 ); <<
265 while( theNuclNucleon ) <<
266 { <<
267 if(theNuclNucleon->AreYouHit()) <<
268 { <<
269 G4LorentzVector tmp=toLab*theNuclNuc <<
270 const G4ParticleDefinition* def = theNuclN <<
271 #ifdef debug_PartonStringModel <<
272 ProjectileResidual4Momentum += tmp; <<
273 hitsP++; <<
274 if ( def == G4Proton::Definition() | <<
275 ExcitationEp +=theNuclNucleon->GetBi <<
276 #endif <<
277 theNuclNucleon->SetMomentum(tmp); <<
278 if ( def == G4Proton::Definition() <<
279 if ( def == G4Neutron::Definition() <<
280 if ( def == G4Lambda::Definition() <<
281 } <<
282 theNuclNucleon = ProjResNucleus->GetNe <<
283 } <<
284 G4int numberLambdaProjectile = 0; <<
285 if ( thePrimary.GetDefinition()->IsHyper <<
286 numberLambdaProjectile = thePrimary.GetDefin <<
287 } else if ( thePrimary.GetDefinition()-> <<
288 numberLambdaProjectile = thePrimary.GetDefin <<
289 } <<
290 #ifdef debug_PartonStringModel <<
291 G4cout<<"Projectile residual A, Z (numbe <<
292 <<thePrimary.GetDefinition()->GetB <<
293 <<thePrimary.GetDefinition()->GetP <<
294 << numberLambdaProjectile - numberLambda <<
295 <<ExcitationEp<<G4endl; <<
296 G4cout<<"Projectile residual 4 momentum <<
297 #endif <<
298 numberProtonProjectileResidual = std::ma <<
299 numberProtonProjectileHits, 0 ); <<
300 numberLambdaProjectileResidual = std::ma <<
301 numberNeutronProjectileResidual = std::m <<
302 <<
303 numberLambdaProjectile - numberN <<
304 } <<
305 <<
306 G4V3DNucleus * ResNucleus = GetWoundedNucl <<
307 <<
308 // loop over wounded nucleus <<
309 theNuclNucleon = ResNucleus->StartLoop() ? <<
310 G4int numberProtonTargetHits( 0 ), numberN <<
311 while( theNuclNucleon ) <<
312 { <<
313 if(theNuclNucleon->AreYouHit()) <<
314 { <<
315 G4LorentzVector tmp=toLab*theNuclNucle <<
316 #ifdef debug_PartonStringModel <<
317 TargetResidual4Momentum += tmp; <<
318 hitsT++; <<
319 if ( theNuclNucleon->GetDefinition() = <<
320 ExcitationEt +=theNuclNucleon->GetBind <<
321 #endif <<
322 theNuclNucleon->SetMomentum(tmp); <<
323 if ( theNuclNucleon->GetDefinition() = <<
324 if ( theNuclNucleon->GetDefinition() = <<
325 } <<
326 theNuclNucleon = ResNucleus->GetNextNucl <<
327 } <<
328 <<
329 #ifdef debug_PartonStringModel <<
330 G4cout<<"Target residual A, Z and E* " <<
331 <<theNucleus.GetA_asInt() - hitsT<<" <<
332 <<theNucleus.GetZ_asInt() - charged_ <<
333 <<ExcitationEt<<G4endl; <<
334 G4cout<<"Target residual 4 momentum "<<Tar <<
335 Bsum+=( hitsT + hitsP); <<
336 Qsum+=(charged_hitsT + charged_hitsP); <<
337 G4cout<<"Hitted # of nucleons of projectil <<
338 G4cout<<"Hitted # of protons of projectile <<
339 <<charged_hitsP<<" "<<charged_hitsT< <<
340 G4cout<<"Bsum Qsum "<<Bsum<<" "<<Qsum<<G4e <<
341 #endif <<
342 <<
343 // Re-sample in the case of unphysical nuc <<
344 // 1 (H), 2 (2He), and 3 (3Li) protons alo <<
345 // no bound states of 2 or more neutrons w <<
346 G4int numberProtonTargetResidual = theNucl <<
347 G4int numberNeutronTargetResidual = theNuc <<
348 G4bool unphysicalResidual = false; <<
349 if ( ( numberProtonTargetResidual > 3 && <<
350 ( numberProtonTargetResidual == 0 && <<
351 unphysicalResidual = true; <<
352 //G4cout << "***UNPHYSICAL TARGET RESIDU <<
353 // << " ; N=" << numberNeutronTarg <<
354 } <<
355 // The projectile residual can be a hypern <<
356 // only the following combinations are cur <<
357 // p-n-lambda (hypertriton), p-n-n-lambda <<
358 // p-p-n-n-lambda (hyperHe5), n-n-lambda-l <<
359 // p-n-lambda-lambda (doubleHyperH4) <<
360 if ( ( numberProtonProjectileResidual > 3 <<
361 ( numberProtonProjectileResidual == 0 <<
362 numberLambdaProjectileResidual == 0 ) || <<
363 ( numberProtonProjectileResidual == 0 && <<
364 numberLambdaProjectileResidual > 0 ) || <<
365 ( numberProtonProjectileResidual == 0 && <<
366 numberLambdaProjectileResidual > 0 ) || <<
367 ( numberLambdaProjectileResidual > 2 ) || <<
368 ( numberProtonProjectileResidual > 0 && n <<
369 numberLambdaProjectileResidual > 0 ) || <<
370 ( numberProtonProjectileResidual > 1 && n <<
371 numberLambdaProjectileResidual > 1 ) <<
372 ) { <<
373 unphysicalResidual = true; <<
374 //G4cout << "***UNPHYSICAL PROJECTILE RE <<
375 // << " ; N=" << numberNeutronProj <<
376 // << " ; L=" << numberLambdaProje <<
377 } <<
378 if ( unphysicalResidual ) { <<
379 //G4cout << " -> REJECTING COLLISION bec <<
380 Success = false; <<
381 continue; <<
382 } <<
383 <<
384 //======================================== <<
385 // Fragment s <<
386 #ifdef debug_PartonStringModel <<
387 G4cout<<"---------------- Attempt to fragm <<
388 #endif <<
389 <<
390 G4double InvMass=SumStringMom.mag(); <<
391 G4double SumMass(0.); <<
392 <<
393 #ifdef debug_PartonStringModel <<
394 QsumSec=0; BsumSec=0; <<
395 SumPsecondr=G4LorentzVector(0.,0.,0.,0.); <<
396 #endif <<
397 <<
398 if(theResult != nullptr) <<
399 { <<
400 std::for_each(theResult->begin(), theRes <<
401 delete theResult; <<
402 } <<
403 <<
404 theResult = stringFragmentationModel->Frag <<
405 <<
406 #ifdef debug_PartonStringModel <<
407 G4cout<<"String fragmentation success (OK <<
408 #endif <<
409 <<
410 if(theResult == 0) {Success=false; continu <<
411 <<
412 #ifdef debug_PartonStringModel <<
413 G4cout<<"Parton-String model: Number of pr <<
414 SumPsecondr = G4LorentzVector(0.,0.,0.,0.) <<
415 QsumSec = 0; BsumSec = 0; <<
416 #endif <<
417 <<
418 SumMass=0.; <<
419 for ( unsigned int i=0; i < theResult->siz <<
420 { <<
421 SumMass+=(*theResult)[i]->Get4Momentum() <<
422 #ifdef debug_PartonStringModel <<
423 G4cout<<i<<" : "<<(*theResult)[i]->GetDe <<
424 <<(*theResult)[i]->Get4M <<
425 <<(*theResult)[i]->Get4M <<
426 <<(*theResult)[i]->GetDe <<
427 SumPsecondr+=(*theResult)[i]->Get4Moment <<
428 BsumSec += (*theResult)[i]->GetDefinitio <<
429 QsumSec += (*theResult)[i]->GetDefinitio <<
430 #endif <<
431 <<
432 #ifdef debug_heavyHadrons <<
433 G4int charm = (*theResult)[i]->GetD <<
434 G4int anticharm = (*theResult)[i]->GetD <<
435 G4int bottom = (*theResult)[i]->GetD <<
436 G4int antibottom = (*theResult)[i]->GetD <<
437 if ( charm != 0 || anticharm != 0 || <<
438 count_charm_hadrons += charm - antich <<
439 count_bottom_hadrons += bottom - antib <<
440 G4cout << "G4VPartonStringModel::Scatter : h <<
441 << (*theResult)[i]->GetDefiniti <<
442 } <<
443 #endif <<
444 } <<
445 <<
446 #ifdef debug_heavyHadrons <<
447 if ( count_charm_projectile != count_charm <<
448 G4cout << "G4VPartonStringModel::Scatter <<
449 << count_charm_projectile << " ; #hadro <<
450 } <<
451 if ( count_bottom_projectile != count_bott <<
452 G4cout << "G4VPartonStringModel::Scatter <<
453 << count_bottom_projectile << " ; #hadr <<
454 } <<
455 #endif <<
456 <<
457 #ifdef debug_PartonStringModel <<
458 G4cout<<G4endl<<"-----------------------Pa <<
459 if(Qsum != QsumSec) { <<
460 G4cout<<"Charge is not conserved!!! ---- <<
461 G4cout<<" Qsum != QsumSec "<<Qsum<<" "<< <<
462 } <<
463 if(Bsum != BsumSec) { <<
464 G4cout<<"Baryon number is not conserved! <<
465 G4cout<<" Bsum != BsumSec "<<Bsum<<" "<< <<
466 } <<
467 #endif <<
468 <<
469 if((SumMass > InvMass)||(SumMass == 0.)) { <<
470 std::for_each(strings->begin(), strings->e <<
471 delete strings; <<
472 <<
473 } while(!Success); <<
474 <<
475 #ifdef debug_PartonStringModel <<
476 G4cout <<"Baryon number balance "<<Bs <<
477 G4cout <<"Charge balance "<<Qs <<
478 G4cout <<"4 momentum balance "<<Su <<
479 G4cout<<"---------------------End of Parton- <<
480 #endif <<
481 95
>> 96 strings = GetStrings();
>> 97 }
>> 98
>> 99 G4double stringEnergy(0);
>> 100 G4LorentzVector SumStringMom(0.,0.,0.,0.);
>> 101
>> 102 #ifdef debug_PartonStringModel
>> 103 G4cout<<"Parton-String model: Number of produced strings "<<strings->size()<<G4endl;
>> 104 #endif
>> 105
>> 106 for ( unsigned int astring=0; astring < strings->size(); astring++)
>> 107 {
>> 108 // rotate string to lab frame, models have it aligned to z
>> 109 stringEnergy += (*strings)[astring]->GetLeftParton()->Get4Momentum().t();
>> 110 stringEnergy += (*strings)[astring]->GetRightParton()->Get4Momentum().t();
>> 111 (*strings)[astring]->LorentzRotate(toLab);
>> 112 SumStringMom+=(*strings)[astring]->Get4Momentum();
>> 113 }
>> 114
>> 115 G4double InvMass=SumStringMom.mag();
>> 116
>> 117 #ifdef debug_PartonStringModel
>> 118 G4cout<<"Parton-String model: SumStringMom "<<SumStringMom<<G4endl;
>> 119
>> 120 G4V3DNucleus * fancynucleus=theThis->GetWoundedNucleus();
>> 121
>> 122 // loop over wounded nucleus
>> 123 G4int hits(0), charged_hits(0);
>> 124 G4ThreeVector hitNucleonMomentum(0.,0.,0.);
>> 125 G4Nucleon * theCurrentNucleon = fancynucleus->StartLoop() ?
>> 126 fancynucleus->GetNextNucleon() : NULL;
>> 127 while( theCurrentNucleon )
>> 128 {
>> 129 if(theCurrentNucleon->AreYouHit())
>> 130 {
>> 131 hits++;
>> 132 hitNucleonMomentum += theCurrentNucleon->Get4Momentum().vect();
>> 133 if ( theCurrentNucleon->GetDefinition() == G4Proton::Proton() ) ++charged_hits;
>> 134 }
>> 135 theCurrentNucleon = fancynucleus->GetNextNucleon();
>> 136 }
>> 137
>> 138 G4int initialZ=fancynucleus->GetCharge();
>> 139 G4int initialA=fancynucleus->GetMassNumber();
>> 140 G4double initial_mass=
>> 141 G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(initialZ,initialA);
>> 142 G4double final_mass =
>> 143 G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(
>> 144 initialZ-charged_hits, initialA-hits);
>> 145 G4cout << "G4VPSM: " <<G4endl
>> 146 << "strE "<<stringEnergy <<G4endl
>> 147 << "hit nucleons "<<hits <<G4endl
>> 148 << "Primary "<<Ptmp.e() <<G4endl
>> 149 << "SumStringE "<<SumStringMom.e() <<G4endl
>> 150 << "Nucleus intial "<<initial_mass <<G4endl
>> 151 << "Nucleus final "<<final_mass <<G4endl
>> 152 << "Excitation estimate "
>> 153 <<Ptmp.e() + initial_mass - final_mass - stringEnergy << G4endl;
>> 154 G4cout << "momentum balance "
>> 155 << thePrimary.GetMomentum() + hitNucleonMomentum -
>> 156 SumStringMom.vect()<< G4endl;
>> 157 #endif
>> 158
>> 159 // Fragment strings
>> 160
>> 161 G4KineticTrackVector * theResult = 0;
>> 162 G4double SumMass(0.);
>> 163 attempts = 0;
>> 164 maxAttempts=100;
>> 165 do
>> 166 {
>> 167 attempts++;
>> 168 if(theResult != 0)
>> 169 {
>> 170 std::for_each(theResult->begin(), theResult->end(), DeleteKineticTrack());
>> 171 delete theResult;
>> 172 }
>> 173 theResult = stringFragmentationModel->FragmentStrings(strings);
>> 174 if(attempts > maxAttempts ) break;
>> 175
>> 176 #ifdef debug_PartonStringModel
>> 177 G4cout<<"Parton-String model: Number of produced particles "<<theResult->size()<<G4endl;
>> 178 G4LorentzVector SumPsecondr(0.,0.,0.,0.);
>> 179 #endif
>> 180
>> 181 SumMass=0.;
>> 182
>> 183 for ( unsigned int i=0; i < theResult->size(); i++)
>> 184 {
>> 185 SumMass+=(*theResult)[i]->GetDefinition()->GetPDGMass();
>> 186 //SumP+=(*theResult)[i]->Get4Momentum();
>> 187 #ifdef debug_PartonStringModel
>> 188 G4cout<<i<<" : "<<(*theResult)[i]->GetDefinition()->GetParticleName()<<" "
>> 189 <<(*theResult)[i]->Get4Momentum()<<" "
>> 190 <<(*theResult)[i]->Get4Momentum().mag()<<G4endl;
>> 191 SumPsecondr+=(*theResult)[i]->Get4Momentum();
>> 192 #endif
>> 193 }
>> 194 #ifdef debug_PartonStringModel
>> 195 G4cout<<"SumP secondaries "<<SumPsecondr<<G4endl;
>> 196 #endif
>> 197 } while(SumMass > InvMass);
>> 198
>> 199 std::for_each(strings->begin(), strings->end(), DeleteString() );
>> 200 delete strings;
>> 201
>> 202 #ifdef debug_PartonStringModel
>> 203 G4cout<<"End of string model work ------------"<<G4endl<<G4endl;
>> 204 #endif
482 return theResult; 205 return theResult;
483 } 206 }
484 207
485 void G4VPartonStringModel::ModelDescription(st 208 void G4VPartonStringModel::ModelDescription(std::ostream& outFile) const
486 { 209 {
487 outFile << GetModelName() << " has no descri << 210 outFile << GetModelName() << " has no description yet.\n";
488 } 211 }
489 212
490 G4V3DNucleus * G4VPartonStringModel::GetProjec << 213 G4V3DNucleus * G4VPartonStringModel::GetProjectileNucleus() const // Uzhi Jan. 2013
491 { return nullptr; } << 214 { return 0;}
492 <<
493 215