| 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 83684 2014-09-09 12:37:39Z 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 // ------------------------------------------- <<
39 38
40 #include "G4VPartonStringModel.hh" 39 #include "G4VPartonStringModel.hh"
41 #include "G4ios.hh" 40 #include "G4ios.hh"
42 #include "G4ShortLivedConstructor.hh" 41 #include "G4ShortLivedConstructor.hh"
43 42
44 #include "G4ParticleTable.hh" 43 #include "G4ParticleTable.hh"
45 #include "G4IonTable.hh" 44 #include "G4IonTable.hh"
46 45
47 G4VPartonStringModel::G4VPartonStringModel(con 46 G4VPartonStringModel::G4VPartonStringModel(const G4String& modelName)
48 : G4VHighEnergyGenerator(modelName), 47 : G4VHighEnergyGenerator(modelName),
49 stringFragmentationModel(nullptr) << 48 stringFragmentationModel(0),
>> 49 theThis(0)
50 { 50 {
51 // Make shure Shotrylived particles are con << 51 // Make shure Shotrylived partyicles are constructed.
52 // VI: should not instantiate particles by << 52 G4ShortLivedConstructor ShortLived;
53 /* << 53 ShortLived.ConstructParticle();
54 G4ShortLivedConstructor ShortLived; <<
55 ShortLived.ConstructParticle(); <<
56 */ <<
57 } 54 }
58 55
59 G4VPartonStringModel::~G4VPartonStringModel() 56 G4VPartonStringModel::~G4VPartonStringModel()
60 { 57 {
61 } 58 }
62 59
63 G4KineticTrackVector * G4VPartonStringModel::S 60 G4KineticTrackVector * G4VPartonStringModel::Scatter(const G4Nucleus &theNucleus,
64 61 const G4DynamicParticle &aPrimary)
65 { 62 {
66 G4ExcitedStringVector * strings = nullptr; << 63 G4ExcitedStringVector * strings = NULL;
>> 64
67 G4DynamicParticle thePrimary=aPrimary; 65 G4DynamicParticle thePrimary=aPrimary;
68 G4LorentzVector SumStringMom(0.,0.,0.,0.); <<
69 G4KineticTrackVector * theResult = 0; <<
70 G4Nucleon * theNuclNucleon(nullptr); <<
71 66
72 #ifdef debug_PartonStringModel << 67 #ifdef debug_PartonStringModel
73 G4cout<<G4endl; 68 G4cout<<G4endl;
74 G4cout<<"-----------------------Parton-Strin 69 G4cout<<"-----------------------Parton-String model is runnung ------------"<<G4endl;
75 G4cout<<"Projectile Name Mass "<<thePrimary. 70 G4cout<<"Projectile Name Mass "<<thePrimary.GetDefinition()->GetParticleName()<<" "
76 <<thePrimary. 71 <<thePrimary.GetMass()<<G4endl;
77 G4cout<<" Momentum "<<thePrimary. 72 G4cout<<" Momentum "<<thePrimary.Get4Momentum()<<G4endl;
78 G4cout<<"Target nucleus A Z "<<theNucleus. 73 G4cout<<"Target nucleus A Z "<<theNucleus.GetA_asInt()<<" "
79 <<theNucleus. 74 <<theNucleus.GetZ_asInt()<<G4endl<<G4endl;
80 G4int Bsum=thePrimary.GetDefinition()->GetBa << 75 #endif
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 76
94 G4LorentzRotation toZ; 77 G4LorentzRotation toZ;
95 G4LorentzVector Ptmp=thePrimary.Get4Momentum 78 G4LorentzVector Ptmp=thePrimary.Get4Momentum();
96 toZ.rotateZ(-1*Ptmp.phi()); 79 toZ.rotateZ(-1*Ptmp.phi());
97 toZ.rotateY(-1*Ptmp.theta()); 80 toZ.rotateY(-1*Ptmp.theta());
98 thePrimary.Set4Momentum(toZ*Ptmp); 81 thePrimary.Set4Momentum(toZ*Ptmp);
99 G4LorentzRotation toLab(toZ.inverse()); 82 G4LorentzRotation toLab(toZ.inverse());
100 83
101 G4bool Success=true; << 84 G4int attempts = 0, maxAttempts=20;
102 G4int attempts = 0, maxAttempts=1000; << 85 while ( strings == NULL )
103 do <<
104 { 86 {
105 if (attempts++ > maxAttempts ) << 87 if (attempts++ > maxAttempts )
106 { << 88 {
107 Init(theNucleus,thePrimary); // To put << 89 throw G4HadronicException(__FILE__, __LINE__,
108 / << 90 "G4VPartonStringModel::Scatter(): fails to generate strings");
109 G4V3DNucleus * ResNucleus = GetWoundedNu << 91 }
110 theNuclNucleon = ResNucleus->StartLoop() << 92 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 93
154 // G4double stringEnergy(0); << 94 strings = GetStrings();
155 SumStringMom=G4LorentzVector(0.,0.,0.,0.); << 95 }
156 << 96
157 #ifdef debug_PartonStringModel << 97 G4double stringEnergy(0);
158 G4cout<<"------------ Parton-String model: << 98 G4LorentzVector SumStringMom(0.,0.,0.,0.);
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 99
306 G4V3DNucleus * ResNucleus = GetWoundedNucl << 100 #ifdef debug_PartonStringModel
>> 101 G4cout<<"Parton-String model: Number of produced strings "<<strings->size()<<G4endl;
>> 102 #endif
307 103
308 // loop over wounded nucleus << 104 for ( unsigned int astring=0; astring < strings->size(); astring++)
309 theNuclNucleon = ResNucleus->StartLoop() ? << 105 {
310 G4int numberProtonTargetHits( 0 ), numberN << 106 // rotate string to lab frame, models have it aligned to z
311 while( theNuclNucleon ) << 107 if((*strings)[astring]->IsExcited())
312 { 108 {
313 if(theNuclNucleon->AreYouHit()) << 109 stringEnergy += (*strings)[astring]->GetLeftParton()->Get4Momentum().t();
314 { << 110 stringEnergy += (*strings)[astring]->GetRightParton()->Get4Momentum().t();
315 G4LorentzVector tmp=toLab*theNuclNucle << 111 (*strings)[astring]->LorentzRotate(toLab);
316 #ifdef debug_PartonStringModel << 112 SumStringMom+=(*strings)[astring]->Get4Momentum();
317 TargetResidual4Momentum += tmp; << 113 #ifdef debug_PartonStringModel
318 hitsT++; << 114 G4cout<<"String No "<<astring+1<<" "<<(*strings)[astring]->Get4Momentum()<<" "
319 if ( theNuclNucleon->GetDefinition() = << 115 <<(*strings)[astring]->Get4Momentum().mag()<<G4endl;
320 ExcitationEt +=theNuclNucleon->GetBind << 116 #endif
321 #endif << 117 }
322 theNuclNucleon->SetMomentum(tmp); << 118 else
323 if ( theNuclNucleon->GetDefinition() = << 119 {
324 if ( theNuclNucleon->GetDefinition() = << 120 stringEnergy += (*strings)[astring]->GetKineticTrack()->Get4Momentum().t();
325 } << 121 (*strings)[astring]->LorentzRotate(toLab);
326 theNuclNucleon = ResNucleus->GetNextNucl << 122 SumStringMom+=(*strings)[astring]->GetKineticTrack()->Get4Momentum();
327 } << 123 #ifdef debug_PartonStringModel
328 << 124 G4cout<<"A track No "<<astring+1<<" "
329 #ifdef debug_PartonStringModel << 125 <<(*strings)[astring]->GetKineticTrack()->Get4Momentum()<<" "
330 G4cout<<"Target residual A, Z and E* " << 126 <<(*strings)[astring]->GetKineticTrack()->Get4Momentum().mag()<<G4endl;
331 <<theNucleus.GetA_asInt() - hitsT<<" << 127 #endif
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 } 128 }
>> 129 }
383 130
384 //======================================== << 131 G4double InvMass=SumStringMom.mag();
385 // Fragment s << 132 G4V3DNucleus * ResNucleus=theThis->GetWoundedNucleus();
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 133
398 if(theResult != nullptr) << 134 // loop over wounded nucleus
399 { << 135 G4Nucleon * theNuclNucleon = ResNucleus->StartLoop() ?
400 std::for_each(theResult->begin(), theRes << 136 ResNucleus->GetNextNucleon() : NULL;
401 delete theResult; << 137 while( theNuclNucleon )
402 } << 138 {
>> 139 if(theNuclNucleon->AreYouHit())
>> 140 {
>> 141 G4LorentzVector tmp=toLab*theNuclNucleon->Get4Momentum();
>> 142 theNuclNucleon->SetMomentum(tmp);
>> 143 }
>> 144 theNuclNucleon = ResNucleus->GetNextNucleon();
>> 145 }
403 146
404 theResult = stringFragmentationModel->Frag << 147 G4V3DNucleus * ProjResNucleus=theThis->GetProjectileNucleus();
405 148
406 #ifdef debug_PartonStringModel << 149 #ifdef debug_PartonStringModel
407 G4cout<<"String fragmentation success (OK << 150 G4ThreeVector hitNucleonMomentum(0.,0.,0.);
408 #endif << 151 #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 152
418 SumMass=0.; << 153 if(ProjResNucleus != 0)
419 for ( unsigned int i=0; i < theResult->siz << 154 {
>> 155 theNuclNucleon = ProjResNucleus->StartLoop() ?
>> 156 ProjResNucleus->GetNextNucleon() : NULL;
>> 157 while( theNuclNucleon )
420 { 158 {
421 SumMass+=(*theResult)[i]->Get4Momentum() << 159 if(theNuclNucleon->AreYouHit())
>> 160 {
>> 161 G4LorentzVector tmp=toLab*theNuclNucleon->Get4Momentum();
422 #ifdef debug_PartonStringModel 162 #ifdef debug_PartonStringModel
423 G4cout<<i<<" : "<<(*theResult)[i]->GetDe << 163 hitNucleonMomentum += tmp.vect();
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 164 #endif
431 << 165 theNuclNucleon->SetMomentum(tmp);
432 #ifdef debug_heavyHadrons << 166 }
433 G4int charm = (*theResult)[i]->GetD << 167 theNuclNucleon = ProjResNucleus->GetNextNucleon();
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 } 168 }
>> 169 }
445 170
446 #ifdef debug_heavyHadrons << 171 #ifdef debug_PartonStringModel
447 if ( count_charm_projectile != count_charm << 172 G4cout<<"Parton-String model: SumStringMom "<<SumStringMom<<G4endl;
448 G4cout << "G4VPartonStringModel::Scatter << 173
449 << count_charm_projectile << " ; #hadro << 174 G4V3DNucleus * fancynucleus=theThis->GetWoundedNucleus();
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 175
469 if((SumMass > InvMass)||(SumMass == 0.)) { << 176 // loop over wounded nucleus
470 std::for_each(strings->begin(), strings->e << 177 G4int hits(0), charged_hits(0);
471 delete strings; << 178 // G4ThreeVector hitNucleonMomentum(0.,0.,0.); // Uzhi Feb. 2014
472 << 179 G4Nucleon * theCurrentNucleon = fancynucleus->StartLoop() ?
473 } while(!Success); << 180 fancynucleus->GetNextNucleon() : NULL;
474 << 181 while( theCurrentNucleon )
475 #ifdef debug_PartonStringModel << 182 {
476 G4cout <<"Baryon number balance "<<Bs << 183 if(theCurrentNucleon->AreYouHit())
477 G4cout <<"Charge balance "<<Qs << 184 {
478 G4cout <<"4 momentum balance "<<Su << 185 hits++;
479 G4cout<<"---------------------End of Parton- << 186 hitNucleonMomentum += theCurrentNucleon->Get4Momentum().vect();
480 #endif << 187 if ( theCurrentNucleon->GetDefinition() == G4Proton::Proton() ) ++charged_hits;
>> 188 }
>> 189 theCurrentNucleon = fancynucleus->GetNextNucleon();
>> 190 }
>> 191
>> 192 G4int initialZ=fancynucleus->GetCharge();
>> 193 G4int initialA=fancynucleus->GetMassNumber();
>> 194 G4double initial_mass=
>> 195 G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(initialZ,initialA);
>> 196 G4double final_mass(0.);
>> 197 if(initialA-hits != 0) final_mass =
>> 198 G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(
>> 199 initialZ-charged_hits, initialA-hits);
>> 200 G4cout << "G4VPSM: " <<G4endl
>> 201 << "strE "<<stringEnergy <<G4endl
>> 202 << "Hit targeet nucleons "<<hits <<G4endl
>> 203 << "Primary "<<Ptmp.e() <<G4endl
>> 204 << "SumStringE "<<SumStringMom.e() <<G4endl
>> 205 << "Target nucleus intial "<<initial_mass <<G4endl
>> 206 << "Target nucleus final "<<final_mass <<G4endl
>> 207 << "Excitation estimate "
>> 208 <<Ptmp.e() + initial_mass - final_mass - stringEnergy << G4endl;
>> 209 G4cout << "momentum balance "
>> 210 << thePrimary.GetMomentum() + hitNucleonMomentum -
>> 211 SumStringMom.vect()<< G4endl;
>> 212 #endif
481 213
>> 214 // Fragment strings
>> 215
>> 216 G4KineticTrackVector * theResult = 0;
>> 217 G4double SumMass(0.);
>> 218 attempts = 0;
>> 219 maxAttempts=100;
>> 220 do
>> 221 {
>> 222 attempts++;
>> 223 if(theResult != 0)
>> 224 {
>> 225 std::for_each(theResult->begin(), theResult->end(), DeleteKineticTrack());
>> 226 delete theResult;
>> 227 }
>> 228
>> 229 theResult = stringFragmentationModel->FragmentStrings(strings);
>> 230
>> 231 if(attempts > maxAttempts ) break;
>> 232
>> 233 #ifdef debug_PartonStringModel
>> 234 G4cout<<"Parton-String model: Number of produced particles "<<theResult->size()<<G4endl;
>> 235 G4LorentzVector SumPsecondr(0.,0.,0.,0.);
>> 236 #endif
>> 237
>> 238 SumMass=0.;
>> 239
>> 240 for ( unsigned int i=0; i < theResult->size(); i++)
>> 241 {
>> 242 SumMass+=(*theResult)[i]->GetDefinition()->GetPDGMass();
>> 243 //SumP+=(*theResult)[i]->Get4Momentum();
>> 244 #ifdef debug_PartonStringModel
>> 245 G4cout<<i<<" : "<<(*theResult)[i]->GetDefinition()->GetParticleName()<<" "
>> 246 <<(*theResult)[i]->Get4Momentum()<<" "
>> 247 <<(*theResult)[i]->Get4Momentum().mag()<<G4endl;
>> 248 SumPsecondr+=(*theResult)[i]->Get4Momentum();
>> 249 #endif
>> 250 }
>> 251 #ifdef debug_PartonStringModel
>> 252 G4cout<<"SumP secondaries "<<SumPsecondr<<G4endl;
>> 253 #endif
>> 254 } while(SumMass > InvMass);
>> 255
>> 256 std::for_each(strings->begin(), strings->end(), DeleteString() );
>> 257 delete strings;
>> 258
>> 259 #ifdef debug_PartonStringModel
>> 260 G4cout<<"End of string model work ------------"<<G4endl<<G4endl;
>> 261 #endif
482 return theResult; 262 return theResult;
483 } 263 }
484 264
485 void G4VPartonStringModel::ModelDescription(st 265 void G4VPartonStringModel::ModelDescription(std::ostream& outFile) const
486 { 266 {
487 outFile << GetModelName() << " has no descri << 267 outFile << GetModelName() << " has no description yet.\n";
488 } 268 }
489 269
490 G4V3DNucleus * G4VPartonStringModel::GetProjec 270 G4V3DNucleus * G4VPartonStringModel::GetProjectileNucleus() const
491 { return nullptr; } << 271 { return 0;}
492 272
493 273