| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
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 M
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(con
48 : G4VHighEnergyGenerator(modelName),
49 stringFragmentationModel(nullptr)
50 {
51 // Make shure Shotrylived particles are con
52 // VI: should not instantiate particles by
53 /*
54 G4ShortLivedConstructor ShortLived;
55 ShortLived.ConstructParticle();
56 */
57 }
58
59 G4VPartonStringModel::~G4VPartonStringModel()
60 {
61 }
62
63 G4KineticTrackVector * G4VPartonStringModel::S
64
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-Strin
75 G4cout<<"Projectile Name Mass "<<thePrimary.
76 <<thePrimary.
77 G4cout<<" Momentum "<<thePrimary.
78 G4cout<<"Target nucleus A Z "<<theNucleus.
79 <<theNucleus.
80 G4int Bsum=thePrimary.GetDefinition()->GetBa
81 G4int Qsum=thePrimary.GetDefinition()->GetPD
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;
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
108 /
109 G4V3DNucleus * ResNucleus = GetWoundedNu
110 theNuclNucleon = ResNucleus->StartLoop()
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
482 return theResult;
483 }
484
485 void G4VPartonStringModel::ModelDescription(st
486 {
487 outFile << GetModelName() << " has no descri
488 }
489
490 G4V3DNucleus * G4VPartonStringModel::GetProjec
491 { return nullptr; }
492
493