| Geant4 Cross Reference |
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 //
27 //
28 //// ------------------------------------------------------------
29 // GEANT 4 class implementation file
30 //
31 // ---------------- G4VPartonStringModel ----------------
32 // by Gunter Folger, May 1998.
33 // abstract class for all Parton String Models
34 // ------------------------------------------------------------
35 // debug switch
36 //#define debug_PartonStringModel
37 //#define debug_heavyHadrons
38 // ------------------------------------------------------------
39
40 #include "G4VPartonStringModel.hh"
41 #include "G4ios.hh"
42 #include "G4ShortLivedConstructor.hh"
43
44 #include "G4ParticleTable.hh"
45 #include "G4IonTable.hh"
46
47 G4VPartonStringModel::G4VPartonStringModel(const G4String& modelName)
48 : G4VHighEnergyGenerator(modelName),
49 stringFragmentationModel(nullptr)
50 {
51 // Make shure Shotrylived particles are constructed.
52 // VI: should not instantiate particles by any model
53 /*
54 G4ShortLivedConstructor ShortLived;
55 ShortLived.ConstructParticle();
56 */
57 }
58
59 G4VPartonStringModel::~G4VPartonStringModel()
60 {
61 }
62
63 G4KineticTrackVector * G4VPartonStringModel::Scatter(const G4Nucleus &theNucleus,
64 const G4DynamicParticle &aPrimary)
65 {
66 G4ExcitedStringVector * strings = nullptr;
67 G4DynamicParticle thePrimary=aPrimary;
68 G4LorentzVector SumStringMom(0.,0.,0.,0.);
69 G4KineticTrackVector * theResult = 0;
70 G4Nucleon * theNuclNucleon(nullptr);
71
72 #ifdef debug_PartonStringModel
73 G4cout<<G4endl;
74 G4cout<<"-----------------------Parton-String model is runnung ------------"<<G4endl;
75 G4cout<<"Projectile Name Mass "<<thePrimary.GetDefinition()->GetParticleName()<<" "
76 <<thePrimary.GetMass()<<G4endl;
77 G4cout<<" Momentum "<<thePrimary.Get4Momentum()<<G4endl;
78 G4cout<<"Target nucleus A Z "<<theNucleus.GetA_asInt()<<" "
79 <<theNucleus.GetZ_asInt()<<G4endl<<G4endl;
80 G4int Bsum=thePrimary.GetDefinition()->GetBaryonNumber() + theNucleus.GetA_asInt();
81 G4int Qsum=thePrimary.GetDefinition()->GetPDGCharge() + theNucleus.GetZ_asInt();
82 G4cout<<"Initial baryon number "<<Bsum<<G4endl;
83 G4cout<<"Initial charge "<<Qsum<<G4endl;
84 G4cout<<"-------------- Parton-String model: Generation of strings -------"<<G4endl<<G4endl;
85 Bsum -= theNucleus.GetA_asInt(); Qsum -= theNucleus.GetZ_asInt();
86 if(GetProjectileNucleus()) {
87 Bsum -= thePrimary.GetDefinition()->GetBaryonNumber();
88 Qsum -= thePrimary.GetDefinition()->GetPDGCharge();
89 }
90 G4int QsumSec(0), BsumSec(0);
91 G4LorentzVector SumPsecondr(0.,0.,0.,0.);
92 #endif
93
94 G4LorentzRotation toZ;
95 G4LorentzVector Ptmp=thePrimary.Get4Momentum();
96 toZ.rotateZ(-1*Ptmp.phi());
97 toZ.rotateY(-1*Ptmp.theta());
98 thePrimary.Set4Momentum(toZ*Ptmp);
99 G4LorentzRotation toLab(toZ.inverse());
100
101 G4bool Success=true;
102 G4int attempts = 0, maxAttempts=1000;
103 do
104 {
105 if (attempts++ > maxAttempts )
106 {
107 Init(theNucleus,thePrimary); // To put a nucleus into ground state
108 // But marks of hitted nucleons are left. They must be erased.
109 G4V3DNucleus * ResNucleus = GetWoundedNucleus();
110 theNuclNucleon = ResNucleus->StartLoop() ? ResNucleus->GetNextNucleon() : nullptr;
111 while( theNuclNucleon )
112 {
113 if(theNuclNucleon->AreYouHit()) theNuclNucleon->Hit(nullptr);
114 theNuclNucleon = ResNucleus->GetNextNucleon();
115 }
116
117 G4V3DNucleus * ProjResNucleus = GetProjectileNucleus();
118 if(ProjResNucleus != 0)
119 {
120 theNuclNucleon = ProjResNucleus->StartLoop() ? ProjResNucleus->GetNextNucleon() : nullptr;
121 while( theNuclNucleon )
122 {
123 if(theNuclNucleon->AreYouHit()) theNuclNucleon->Hit(nullptr);
124 theNuclNucleon = ProjResNucleus->GetNextNucleon();
125 }
126 }
127
128 G4ExceptionDescription ed;
129 ed << "Projectile Name Mass " <<thePrimary.GetDefinition()->GetParticleName()
130 << " " << thePrimary.GetMass()<< G4endl;
131 ed << " Momentum " << thePrimary.Get4Momentum() <<G4endl;
132 ed << "Target nucleus A Z " << theNucleus.GetA_asInt() << " "
133 << theNucleus.GetZ_asInt() <<G4endl;
134 ed << "Initial states of projectile and target nucleus will be returned!"<<G4endl;
135 G4Exception( "G4VPartonStringModel::Scatter(): fails to generate or fragment strings ",
136 "HAD_PARTON_STRING_001", JustWarning, ed );
137
138 G4ThreeVector Position(0.,0.,2*ResNucleus->GetOuterRadius());
139 G4KineticTrack* Hadron = new G4KineticTrack(aPrimary.GetParticleDefinition(), 0.,
140 Position, aPrimary.Get4Momentum());
141 if(theResult == nullptr) theResult = new G4KineticTrackVector();
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; continue; }
153
154 // G4double stringEnergy(0);
155 SumStringMom=G4LorentzVector(0.,0.,0.,0.);
156
157 #ifdef debug_PartonStringModel
158 G4cout<<"------------ Parton-String model: Number of produced strings ---- "<<strings->size()<<G4endl;
159 #endif
160
161 #ifdef debug_heavyHadrons
162 // Check charm and bottom numbers of the projectile:
163 G4int count_charm_projectile = thePrimary.GetDefinition()->GetQuarkContent( 4 ) -
164 thePrimary.GetDefinition()->GetAntiQuarkContent( 4 );
165 G4int count_bottom_projectile = thePrimary.GetDefinition()->GetQuarkContent( 5 ) -
166 thePrimary.GetDefinition()->GetAntiQuarkContent( 5 );
167 G4int count_charm_strings = 0, count_bottom_strings = 0;
168 G4int count_charm_hadrons = 0, count_bottom_hadrons = 0;
169 #endif
170
171 for ( unsigned int astring=0; astring < strings->size(); astring++)
172 {
173 // rotate string to lab frame, models have it aligned to z
174 if((*strings)[astring]->IsExcited())
175 {
176 // stringEnergy += (*strings)[astring]->GetLeftParton()->Get4Momentum().t();
177 // stringEnergy += (*strings)[astring]->GetRightParton()->Get4Momentum().t();
178 (*strings)[astring]->LorentzRotate(toLab);
179 SumStringMom+=(*strings)[astring]->Get4Momentum();
180 #ifdef debug_PartonStringModel
181 G4cout<<"String No "<<astring+1<<" "<<(*strings)[astring]->Get4Momentum()<<" "
182 <<(*strings)[astring]->Get4Momentum().mag()
183 <<" Partons "<<(*strings)[astring]->GetLeftParton()->GetDefinition()->GetPDGEncoding()
184 <<" "<<(*strings)[astring]->GetRightParton()->GetDefinition()->GetPDGEncoding()<<G4endl;
185 #endif
186
187 #ifdef debug_heavyHadrons
188 G4int left_charm = (*strings)[astring]->GetLeftParton()->GetDefinition()->GetQuarkContent( 4 );
189 G4int left_anticharm = (*strings)[astring]->GetLeftParton()->GetDefinition()->GetAntiQuarkContent( 4 );
190 G4int right_charm = (*strings)[astring]->GetRightParton()->GetDefinition()->GetQuarkContent( 4 );
191 G4int right_anticharm = (*strings)[astring]->GetRightParton()->GetDefinition()->GetAntiQuarkContent( 4 );
192 G4int left_bottom = (*strings)[astring]->GetLeftParton()->GetDefinition()->GetQuarkContent( 5 );
193 G4int left_antibottom = (*strings)[astring]->GetLeftParton()->GetDefinition()->GetAntiQuarkContent( 5 );
194 G4int right_bottom = (*strings)[astring]->GetRightParton()->GetDefinition()->GetQuarkContent( 5 );
195 G4int right_antibottom = (*strings)[astring]->GetRightParton()->GetDefinition()->GetAntiQuarkContent( 5 );
196 if ( left_charm != 0 || left_anticharm != 0 || right_charm != 0 || right_anticharm != 0 ||
197 left_bottom != 0 || left_antibottom != 0 || right_bottom != 0 || right_antibottom != 0 ) {
198 count_charm_strings += left_charm - left_anticharm + right_charm - right_anticharm;
199 count_bottom_strings += left_bottom - left_antibottom + right_bottom - right_antibottom;
200 G4cout << "G4VPartonStringModel::Scatter : string #" << astring << " ("
201 << (*strings)[astring]->GetLeftParton()->GetDefinition()->GetParticleName() << " , "
202 << (*strings)[astring]->GetRightParton()->GetDefinition()->GetParticleName() << ")" << G4endl;
203 }
204 #endif
205 }
206 else
207 {
208 // stringEnergy += (*strings)[astring]->GetKineticTrack()->Get4Momentum().t();
209 (*strings)[astring]->LorentzRotate(toLab);
210 SumStringMom+=(*strings)[astring]->GetKineticTrack()->Get4Momentum();
211 #ifdef debug_PartonStringModel
212 G4cout<<"A track No "<<astring+1<<" "
213 <<(*strings)[astring]->GetKineticTrack()->Get4Momentum()<<" "
214 <<(*strings)[astring]->GetKineticTrack()->Get4Momentum().mag()<<" "
215 <<(*strings)[astring]->GetKineticTrack()->GetDefinition()->GetParticleName()<<G4endl;
216 #endif
217
218 #ifdef debug_heavyHadrons
219 G4int charm = (*strings)[astring]->GetKineticTrack()->GetDefinition()->GetQuarkContent( 4 );
220 G4int anticharm = (*strings)[astring]->GetKineticTrack()->GetDefinition()->GetAntiQuarkContent( 4 );
221 G4int bottom = (*strings)[astring]->GetKineticTrack()->GetDefinition()->GetQuarkContent( 5 );
222 G4int antibottom = (*strings)[astring]->GetKineticTrack()->GetDefinition()->GetAntiQuarkContent( 5 );
223 if ( charm != 0 || anticharm != 0 || bottom != 0 || antibottom != 0 ) {
224 count_charm_strings += charm - anticharm;
225 count_bottom_strings += bottom - antibottom;
226 G4cout << "G4VPartonStringModel::Scatter : track #" << astring << "\t"
227 << (*strings)[astring]->GetKineticTrack()->GetDefinition()->GetParticleName() << G4endl;
228 }
229 #endif
230 }
231 }
232
233 #ifdef debug_heavyHadrons
234 if ( count_charm_projectile != count_charm_strings ) {
235 G4cout << "G4VPartonStringModel::Scatter : CHARM VIOLATION in String formation ! #projectile="
236 << count_charm_projectile << " ; #strings=" << count_charm_strings << G4endl;
237 }
238 if ( count_bottom_projectile != count_bottom_strings ) {
239 G4cout << "G4VPartonStringModel::Scatter : BOTTOM VIOLATION in String formation ! #projectile="
240 << count_bottom_projectile << " ; #strings=" << count_bottom_strings << G4endl;
241 }
242 #endif
243
244 #ifdef debug_PartonStringModel
245 G4cout<<G4endl<<"SumString4Mom "<<SumStringMom<<G4endl;
246 G4LorentzVector TargetResidual4Momentum(0.,0.,0.,0.);
247 G4LorentzVector ProjectileResidual4Momentum(0.,0.,0.,0.);
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 never a hypernucleus, whereas
254 // the projectile nucleus can be a light hypernucleus or anti-hypernucleus.
255
256 G4V3DNucleus * ProjResNucleus = GetProjectileNucleus();
257
258 G4int numberProtonProjectileResidual( 0 ), numberNeutronProjectileResidual( 0 );
259 G4int numberLambdaProjectileResidual( 0 );
260 if(ProjResNucleus != 0)
261 {
262 theNuclNucleon = ProjResNucleus->StartLoop() ? ProjResNucleus->GetNextNucleon() : nullptr;
263 G4int numberProtonProjectileHits( 0 ), numberNeutronProjectileHits( 0 );
264 G4int numberLambdaProjectileHits( 0 );
265 while( theNuclNucleon )
266 {
267 if(theNuclNucleon->AreYouHit())
268 {
269 G4LorentzVector tmp=toLab*theNuclNucleon->Get4Momentum();
270 const G4ParticleDefinition* def = theNuclNucleon->GetDefinition();
271 #ifdef debug_PartonStringModel
272 ProjectileResidual4Momentum += tmp;
273 hitsP++;
274 if ( def == G4Proton::Definition() || def == G4AntiProton::Definition() ) ++charged_hitsP;
275 ExcitationEp +=theNuclNucleon->GetBindingEnergy();
276 #endif
277 theNuclNucleon->SetMomentum(tmp);
278 if ( def == G4Proton::Definition() || def == G4AntiProton::Definition() ) ++numberProtonProjectileHits;
279 if ( def == G4Neutron::Definition() || def == G4AntiNeutron::Definition() ) ++numberNeutronProjectileHits;
280 if ( def == G4Lambda::Definition() || def == G4AntiLambda::Definition() ) ++numberLambdaProjectileHits;
281 }
282 theNuclNucleon = ProjResNucleus->GetNextNucleon();
283 }
284 G4int numberLambdaProjectile = 0;
285 if ( thePrimary.GetDefinition()->IsHypernucleus() ) {
286 numberLambdaProjectile = thePrimary.GetDefinition()->GetNumberOfLambdasInHypernucleus();
287 } else if ( thePrimary.GetDefinition()->IsAntiHypernucleus() ) {
288 numberLambdaProjectile = thePrimary.GetDefinition()->GetNumberOfAntiLambdasInAntiHypernucleus();
289 }
290 #ifdef debug_PartonStringModel
291 G4cout<<"Projectile residual A, Z (numberOfLambdasOrAntiLambdas) and E* "
292 <<thePrimary.GetDefinition()->GetBaryonNumber() - hitsP<<" "
293 <<thePrimary.GetDefinition()->GetPDGCharge() - charged_hitsP<<" ("
294 << numberLambdaProjectile - numberLambdaProjectileHits << ") "
295 <<ExcitationEp<<G4endl;
296 G4cout<<"Projectile residual 4 momentum "<<ProjectileResidual4Momentum<<G4endl;
297 #endif
298 numberProtonProjectileResidual = std::max( std::abs( G4int( thePrimary.GetDefinition()->GetPDGCharge() ) ) -
299 numberProtonProjectileHits, 0 );
300 numberLambdaProjectileResidual = std::max( numberLambdaProjectile - numberLambdaProjectileHits, 0 );
301 numberNeutronProjectileResidual = std::max( std::abs( thePrimary.GetDefinition()->GetBaryonNumber() ) -
302 std::abs( G4int( thePrimary.GetDefinition()->GetPDGCharge() ) ) -
303 numberLambdaProjectile - numberNeutronProjectileHits, 0 );
304 }
305
306 G4V3DNucleus * ResNucleus = GetWoundedNucleus();
307
308 // loop over wounded nucleus
309 theNuclNucleon = ResNucleus->StartLoop() ? ResNucleus->GetNextNucleon() : nullptr;
310 G4int numberProtonTargetHits( 0 ), numberNeutronTargetHits( 0 );
311 while( theNuclNucleon )
312 {
313 if(theNuclNucleon->AreYouHit())
314 {
315 G4LorentzVector tmp=toLab*theNuclNucleon->Get4Momentum();
316 #ifdef debug_PartonStringModel
317 TargetResidual4Momentum += tmp;
318 hitsT++;
319 if ( theNuclNucleon->GetDefinition() == G4Proton::Proton() ) ++charged_hitsT;
320 ExcitationEt +=theNuclNucleon->GetBindingEnergy();
321 #endif
322 theNuclNucleon->SetMomentum(tmp);
323 if ( theNuclNucleon->GetDefinition() == G4Proton::Proton() ) ++numberProtonTargetHits;
324 if ( theNuclNucleon->GetDefinition() == G4Neutron::Neutron() ) ++numberNeutronTargetHits;
325 }
326 theNuclNucleon = ResNucleus->GetNextNucleon();
327 }
328
329 #ifdef debug_PartonStringModel
330 G4cout<<"Target residual A, Z and E* "
331 <<theNucleus.GetA_asInt() - hitsT<<" "
332 <<theNucleus.GetZ_asInt() - charged_hitsT<<" "
333 <<ExcitationEt<<G4endl;
334 G4cout<<"Target residual 4 momentum "<<TargetResidual4Momentum<<G4endl;
335 Bsum+=( hitsT + hitsP);
336 Qsum+=(charged_hitsT + charged_hitsP);
337 G4cout<<"Hitted # of nucleons of projectile and target "<<hitsP<<" "<<hitsT<<G4endl;
338 G4cout<<"Hitted # of protons of projectile and target "
339 <<charged_hitsP<<" "<<charged_hitsT<<G4endl<<G4endl;
340 G4cout<<"Bsum Qsum "<<Bsum<<" "<<Qsum<<G4endl<<G4endl;
341 #endif
342
343 // Re-sample in the case of unphysical nuclear residual:
344 // 1 (H), 2 (2He), and 3 (3Li) protons alone without neutrons can exist, but not more;
345 // no bound states of 2 or more neutrons without protons can exist.
346 G4int numberProtonTargetResidual = theNucleus.GetZ_asInt() - numberProtonTargetHits;
347 G4int numberNeutronTargetResidual = theNucleus.GetA_asInt() - theNucleus.GetZ_asInt() - numberNeutronTargetHits;
348 G4bool unphysicalResidual = false;
349 if ( ( numberProtonTargetResidual > 3 && numberNeutronTargetResidual == 0 ) ||
350 ( numberProtonTargetResidual == 0 && numberNeutronTargetResidual > 1 ) ) {
351 unphysicalResidual = true;
352 //G4cout << "***UNPHYSICAL TARGET RESIDUAL*** Z=" << numberProtonTargetResidual
353 // << " ; N=" << numberNeutronTargetResidual;
354 }
355 // The projectile residual can be a hypernucleus or anti-hypernucleus:
356 // only the following combinations are currently allowed in Geant4:
357 // p-n-lambda (hypertriton), p-n-n-lambda (hyperH4), p-p-n-lambda (hyperAlpha),
358 // p-p-n-n-lambda (hyperHe5), n-n-lambda-lambda (doublehyperdoubleneutron),
359 // p-n-lambda-lambda (doubleHyperH4)
360 if ( ( numberProtonProjectileResidual > 3 && numberNeutronProjectileResidual == 0 ) ||
361 ( numberProtonProjectileResidual == 0 && numberNeutronProjectileResidual > 1 &&
362 numberLambdaProjectileResidual == 0 ) ||
363 ( numberProtonProjectileResidual == 0 && numberNeutronProjectileResidual <= 1 &&
364 numberLambdaProjectileResidual > 0 ) ||
365 ( numberProtonProjectileResidual == 0 && numberNeutronProjectileResidual > 2 &&
366 numberLambdaProjectileResidual > 0 ) ||
367 ( numberLambdaProjectileResidual > 2 ) ||
368 ( numberProtonProjectileResidual > 0 && numberNeutronProjectileResidual == 0 &&
369 numberLambdaProjectileResidual > 0 ) ||
370 ( numberProtonProjectileResidual > 1 && numberNeutronProjectileResidual > 1 &&
371 numberLambdaProjectileResidual > 1 )
372 ) {
373 unphysicalResidual = true;
374 //G4cout << "***UNPHYSICAL PROJECTILE RESIDUAL*** Z=" << numberProtonProjectileResidual
375 // << " ; N=" << numberNeutronProjectileResidual
376 // << " ; L=" << numberLambdaProjectileResidual;
377 }
378 if ( unphysicalResidual ) {
379 //G4cout << " -> REJECTING COLLISION because of unphysical residual !" << G4endl;
380 Success = false;
381 continue;
382 }
383
384 //=========================================================================================
385 // Fragment strings
386 #ifdef debug_PartonStringModel
387 G4cout<<"---------------- Attempt to fragment strings ------------- "<<attempts<<G4endl;
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(), theResult->end(), DeleteKineticTrack());
401 delete theResult;
402 }
403
404 theResult = stringFragmentationModel->FragmentStrings(strings);
405
406 #ifdef debug_PartonStringModel
407 G4cout<<"String fragmentation success (OK - #0, Not OK - 0) : "<<theResult<<G4endl;
408 #endif
409
410 if(theResult == 0) {Success=false; continue;}
411
412 #ifdef debug_PartonStringModel
413 G4cout<<"Parton-String model: Number of produced particles "<<theResult->size()<<G4endl;
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->size(); i++)
420 {
421 SumMass+=(*theResult)[i]->Get4Momentum().mag();
422 #ifdef debug_PartonStringModel
423 G4cout<<i<<" : "<<(*theResult)[i]->GetDefinition()->GetParticleName()<<" "
424 <<(*theResult)[i]->Get4Momentum()<<" "
425 <<(*theResult)[i]->Get4Momentum().mag()<<" "
426 <<(*theResult)[i]->GetDefinition()->GetPDGMass()<<G4endl;
427 SumPsecondr+=(*theResult)[i]->Get4Momentum();
428 BsumSec += (*theResult)[i]->GetDefinition()->GetBaryonNumber();
429 QsumSec += (*theResult)[i]->GetDefinition()->GetPDGCharge();
430 #endif
431
432 #ifdef debug_heavyHadrons
433 G4int charm = (*theResult)[i]->GetDefinition()->GetQuarkContent( 4 );
434 G4int anticharm = (*theResult)[i]->GetDefinition()->GetAntiQuarkContent( 4 );
435 G4int bottom = (*theResult)[i]->GetDefinition()->GetQuarkContent( 5 );
436 G4int antibottom = (*theResult)[i]->GetDefinition()->GetAntiQuarkContent( 5 );
437 if ( charm != 0 || anticharm != 0 || bottom != 0 || antibottom != 0 ) {
438 count_charm_hadrons += charm - anticharm;
439 count_bottom_hadrons += bottom - antibottom;
440 G4cout << "G4VPartonStringModel::Scatter : hadron #" << i << "\t"
441 << (*theResult)[i]->GetDefinition()->GetParticleName() << G4endl;
442 }
443 #endif
444 }
445
446 #ifdef debug_heavyHadrons
447 if ( count_charm_projectile != count_charm_hadrons ) {
448 G4cout << "G4VPartonStringModel::Scatter : CHARM VIOLATION in String hadronization ! #projectile="
449 << count_charm_projectile << " ; #hadrons=" << count_charm_hadrons << G4endl;
450 }
451 if ( count_bottom_projectile != count_bottom_hadrons ) {
452 G4cout << "G4VPartonStringModel::Scatter : BOTTOM VIOLATION in String hadronization ! #projectile="
453 << count_bottom_projectile << " ; #hadrons=" << count_bottom_hadrons << G4endl;
454 }
455 #endif
456
457 #ifdef debug_PartonStringModel
458 G4cout<<G4endl<<"-----------------------Parton-String model: balances -------------"<<G4endl;
459 if(Qsum != QsumSec) {
460 G4cout<<"Charge is not conserved!!! ----"<<G4endl;
461 G4cout<<" Qsum != QsumSec "<<Qsum<<" "<<QsumSec<<G4endl;
462 }
463 if(Bsum != BsumSec) {
464 G4cout<<"Baryon number is not conserved!!!"<<G4endl;
465 G4cout<<" Bsum != BsumSec "<<Bsum<<" "<<BsumSec<<G4endl;
466 }
467 #endif
468
469 if((SumMass > InvMass)||(SumMass == 0.)) {Success=false;}
470 std::for_each(strings->begin(), strings->end(), DeleteString() );
471 delete strings;
472
473 } while(!Success);
474
475 #ifdef debug_PartonStringModel
476 G4cout <<"Baryon number balance "<<Bsum-BsumSec<<G4endl;
477 G4cout <<"Charge balance "<<Qsum-QsumSec<<G4endl;
478 G4cout <<"4 momentum balance "<<SumStringMom-SumPsecondr<<G4endl;
479 G4cout<<"---------------------End of Parton-String model work -------------"<<G4endl<<G4endl;
480 #endif
481
482 return theResult;
483 }
484
485 void G4VPartonStringModel::ModelDescription(std::ostream& outFile) const
486 {
487 outFile << GetModelName() << " has no description yet.\n";
488 }
489
490 G4V3DNucleus * G4VPartonStringModel::GetProjectileNucleus() const
491 { return nullptr; }
492
493