| 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$
>> 28 // GEANT4 tag $Name: $
27 // 29 //
28 30
29 // ------------------------------------------- 31 // ------------------------------------------------------------
30 // GEANT 4 class implementation file 32 // GEANT 4 class implementation file
31 // 33 //
32 // ---------------- G4FTFModel ---------- 34 // ---------------- G4FTFModel ----------------
33 // by Gunter Folger, May 1998. 35 // by Gunter Folger, May 1998.
34 // class implementing the excitation in 36 // class implementing the excitation in the FTF Parton String Model
35 // <<
36 // Vladimir Uzhinsky, November <<
37 // simulation of nucleus-nucleus interac <<
38 // ------------------------------------------- 37 // ------------------------------------------------------------
39 38
40 #include <utility> 39 #include <utility>
41 40
42 #include "G4FTFModel.hh" 41 #include "G4FTFModel.hh"
43 #include "G4ios.hh" 42 #include "G4ios.hh"
44 #include "G4PhysicalConstants.hh" 43 #include "G4PhysicalConstants.hh"
45 #include "G4SystemOfUnits.hh" 44 #include "G4SystemOfUnits.hh"
46 #include "G4FTFParameters.hh" 45 #include "G4FTFParameters.hh"
47 #include "G4FTFParticipants.hh" 46 #include "G4FTFParticipants.hh"
48 #include "G4DiffractiveSplitableHadron.hh" 47 #include "G4DiffractiveSplitableHadron.hh"
49 #include "G4InteractionContent.hh" 48 #include "G4InteractionContent.hh"
50 #include "G4LorentzRotation.hh" 49 #include "G4LorentzRotation.hh"
51 #include "G4ParticleDefinition.hh" 50 #include "G4ParticleDefinition.hh"
52 #include "G4ParticleTable.hh" 51 #include "G4ParticleTable.hh"
53 #include "G4IonTable.hh" 52 #include "G4IonTable.hh"
54 #include "G4KineticTrack.hh" <<
55 #include "G4HyperNucleiProperties.hh" <<
56 #include "G4Exp.hh" <<
57 #include "G4Log.hh" <<
58 <<
59 //============================================ <<
60 <<
61 //#define debugFTFmodel <<
62 //#define debugReggeonCascade <<
63 //#define debugPutOnMassShell <<
64 //#define debugAdjust <<
65 //#define debugBuildString <<
66 <<
67 <<
68 //============================================ <<
69 <<
70 G4FTFModel::G4FTFModel( const G4String& modelN <<
71 G4VPartonStringModel( modelName ), <<
72 theExcitation( new G4DiffractiveExcitation() <<
73 theElastic( new G4ElasticHNScattering() ), <<
74 theAnnihilation( new G4FTFAnnihilation() ) <<
75 { <<
76 // ---> JVY theParameters = 0; <<
77 theParameters = new G4FTFParameters(); <<
78 // <<
79 NumberOfInvolvedNucleonsOfTarget = 0; <<
80 NumberOfInvolvedNucleonsOfProjectile= 0; <<
81 for ( G4int i = 0; i < 250; ++i ) { <<
82 TheInvolvedNucleonsOfTarget[i] = 0; <<
83 TheInvolvedNucleonsOfProjectile[i] = 0; <<
84 } <<
85 <<
86 //LowEnergyLimit = 2000.0*MeV; <<
87 LowEnergyLimit = 1000.0*MeV; <<
88 <<
89 HighEnergyInter = true; <<
90 <<
91 G4LorentzVector tmp( 0.0, 0.0, 0.0, 0.0 ); <<
92 ProjectileResidual4Momentum = tmp; <<
93 ProjectileResidualMassNumber = 0; <<
94 ProjectileResidualCharge = 0; <<
95 ProjectileResidualLambdaNumber = 0; <<
96 ProjectileResidualExcitationEnergy = 0.0; <<
97 <<
98 TargetResidual4Momentum = tmp; <<
99 TargetResidualMassNumber = 0; <<
100 TargetResidualCharge = 0; <<
101 TargetResidualExcitationEnergy = 0.0; <<
102 <<
103 Bimpact = -1.0; <<
104 BinInterval = false; <<
105 Bmin = 0.0; <<
106 Bmax = 0.0; <<
107 NumberOfProjectileSpectatorNucleons = 0; <<
108 NumberOfTargetSpectatorNucleons = 0; <<
109 NumberOfNNcollisions = 0; <<
110 <<
111 SetEnergyMomentumCheckLevels( 2.0*perCent, 1 <<
112 } <<
113 <<
114 53
115 //============================================ << 54 // Class G4FTFModel
116 <<
117 struct DeleteVSplitableHadron { void operator( <<
118 55
>> 56 G4FTFModel::G4FTFModel(const G4String& modelName):G4VPartonStringModel(modelName),
>> 57 theExcitation(new G4DiffractiveExcitation()),
>> 58 theElastic(new G4ElasticHNScattering()),
>> 59 theAnnihilation(new G4FTFAnnihilation())
>> 60 {
>> 61 G4VPartonStringModel::SetThisPointer(this);
>> 62 theParameters=0;
>> 63 NumberOfInvolvedNucleon=0;
>> 64 NumberOfInvolvedNucleonOfProjectile=0;
>> 65 SetEnergyMomentumCheckLevels(2*perCent, 150*MeV);
>> 66 }
119 67
120 //============================================ << 68 struct DeleteVSplitableHadron { void operator()(G4VSplitableHadron * aH){ delete aH;} };
121 69
122 G4FTFModel::~G4FTFModel() { << 70 G4FTFModel::~G4FTFModel()
123 // Because FTF model can be called for vari << 71 {
124 // << 72 // Because FTF model can be called for various particles
125 // ---> NOTE (JVY): This statement below is << 73 // theParameters must be erased at the end of each call.
126 // theParameters must be erased at the end << 74 // Thus the delete is also in G4FTFModel::GetStrings() method
127 // Thus the delete is also in G4FTFModel::G << 75 if( theParameters != 0 ) delete theParameters;
128 // ---> JVY << 76 if( theExcitation != 0 ) delete theExcitation;
129 // << 77 if( theElastic != 0 ) delete theElastic;
130 if ( theParameters != 0 ) delete theParam << 78 if( theAnnihilation != 0 ) delete theAnnihilation;
131 if ( theExcitation != 0 ) delete theExcit << 79
132 if ( theElastic != 0 ) delete theElast << 80 // Erasing of strings created at annihilation
133 if ( theAnnihilation != 0 ) delete theAnnih << 81 if(theAdditionalString.size() != 0)
134 << 82 {
135 // Erasing of strings created at annihilati << 83 std::for_each(theAdditionalString.begin(), theAdditionalString.end(),
136 if ( theAdditionalString.size() != 0 ) { << 84 DeleteVSplitableHadron());
137 std::for_each( theAdditionalString.begin( <<
138 DeleteVSplitableHadron() ) <<
139 } 85 }
140 theAdditionalString.clear(); 86 theAdditionalString.clear();
141 87
142 // Erasing of target involved nucleons. << 88 // Erasing of target involved nucleons
143 if ( NumberOfInvolvedNucleonsOfTarget != 0 << 89 if( NumberOfInvolvedNucleon != 0)
144 for ( G4int i = 0; i < NumberOfInvolvedNu << 90 {
145 G4VSplitableHadron* aNucleon = TheInvol << 91 for(G4int i=0; i < NumberOfInvolvedNucleon; i++)
146 if ( aNucleon ) delete aNucleon; << 92 {
147 } << 93 G4VSplitableHadron * aNucleon = TheInvolvedNucleon[i]->GetSplitableHadron();
148 } << 94 if(aNucleon) delete aNucleon;
149 <<
150 // Erasing of projectile involved nucleons. <<
151 if ( NumberOfInvolvedNucleonsOfProjectile ! <<
152 for ( G4int i = 0; i < NumberOfInvolvedNu <<
153 G4VSplitableHadron* aNucleon = TheInvol <<
154 if ( aNucleon ) delete aNucleon; <<
155 } <<
156 } <<
157 } <<
158 <<
159 <<
160 //============================================ <<
161 <<
162 void G4FTFModel::Init( const G4Nucleus& aNucle <<
163 <<
164 theProjectile = aProjectile; <<
165 <<
166 G4double PlabPerParticle( 0.0 ); // Laborat <<
167 <<
168 #ifdef debugFTFmodel <<
169 G4cout << "FTF init Proj Name " << theProjec <<
170 << "FTF init Proj Mass " << theProjec <<
171 << " " << theProjectile.GetMomentum() <<
172 << "FTF init Proj B Q " << theProjec <<
173 << " " << (G4int) theProjectile.GetDe <<
174 << "FTF init Target A Z " << aNucleus <<
175 << " " << aNucleus.GetZ_asInt() << G4 <<
176 #endif <<
177 <<
178 theParticipants.Clean(); <<
179 <<
180 theParticipants.SetProjectileNucleus( 0 ); <<
181 <<
182 G4LorentzVector tmp( 0.0, 0.0, 0.0, 0.0 ); <<
183 ProjectileResidualMassNumber = 0; <<
184 ProjectileResidualCharge = 0; <<
185 ProjectileResidualLambdaNumber = 0; <<
186 ProjectileResidualExcitationEnergy = 0.0; <<
187 ProjectileResidual4Momentum = tmp; <<
188 <<
189 TargetResidualMassNumber = aNucleus.Ge <<
190 TargetResidualCharge = aNucleus.Ge <<
191 TargetResidualExcitationEnergy = 0.0; <<
192 TargetResidual4Momentum = tmp; <<
193 G4double TargetResidualMass = G4ParticleTabl <<
194 ->GetIonMass( <<
195 <<
196 TargetResidual4Momentum.setE( TargetResidual <<
197 <<
198 if ( std::abs( theProjectile.GetDefinition() <<
199 // Projectile is a hadron : meson or baryo <<
200 ProjectileResidualMassNumber = std::abs( t <<
201 ProjectileResidualCharge = G4int( theProje <<
202 PlabPerParticle = theProjectile.GetMomentu <<
203 ProjectileResidualExcitationEnergy = 0.0; <<
204 //G4double ProjectileResidualMass = thePro <<
205 ProjectileResidual4Momentum.setVect( thePr <<
206 ProjectileResidual4Momentum.setE( theProje <<
207 if ( PlabPerParticle < LowEnergyLimit ) { <<
208 HighEnergyInter = false; <<
209 } else { <<
210 HighEnergyInter = true; <<
211 } <<
212 } else { <<
213 if ( theProjectile.GetDefinition()->GetBar <<
214 // Projectile is a nucleus <<
215 ProjectileResidualMassNumber = theProjec <<
216 ProjectileResidualCharge = G4int( thePro <<
217 ProjectileResidualLambdaNumber = theProj <<
218 PlabPerParticle = theProjectile.GetMomen <<
219 if ( PlabPerParticle < LowEnergyLimit ) <<
220 HighEnergyInter = false; <<
221 } else { <<
222 HighEnergyInter = true; <<
223 } <<
224 theParticipants.InitProjectileNucleus( P <<
225 ProjectileResidualLambdaNumber <<
226 } else if ( theProjectile.GetDefinition()- <<
227 // Projectile is an anti-nucleus <<
228 ProjectileResidualMassNumber = std::abs( <<
229 ProjectileResidualCharge = std::abs( G4i <<
230 ProjectileResidualLambdaNumber = theProj <<
231 PlabPerParticle = theProjectile.GetMomen <<
232 if ( PlabPerParticle < LowEnergyLimit ) <<
233 HighEnergyInter = false; <<
234 } else { <<
235 HighEnergyInter = true; <<
236 } <<
237 theParticipants.InitProjectileNucleus( P <<
238 P <<
239 theParticipants.GetProjectileNucleus()-> <<
240 G4Nucleon* aNucleon; <<
241 while ( ( aNucleon = theParticipants.Get <<
242 if ( aNucleon->GetDefinition() == G4Pr <<
243 aNucleon->SetParticleType( G4AntiPro <<
244 } else if ( aNucleon->GetDefinition() <<
245 aNucleon->SetParticleType( G4AntiNeu <<
246 } else if ( aNucleon->GetDefinition() <<
247 aNucleon->SetParticleType( G4AntiLambda::D <<
248 } <<
249 } <<
250 } 95 }
>> 96 }
251 97
252 G4ThreeVector BoostVector = theProjectile. << 98 // Erasing of projectile involved nucleons
253 theParticipants.GetProjectileNucleus()->Do << 99 /*
254 theParticipants.GetProjectileNucleus()->Do << 100 if( NumberOfInvolvedNucleonOfProjectile != 0)
255 ProjectileResidualExcitationEnergy = 0.0; << 101 {
256 //G4double ProjectileResidualMass = thePro << 102 for(G4int i=0; i < NumberOfInvolvedNucleonOfProjectile; i++)
257 ProjectileResidual4Momentum.setVect( thePr << 103 {
258 ProjectileResidual4Momentum.setE( theProje << 104 G4VSplitableHadron * aNucleon = TheInvolvedNucleonOfProjectile[i]->GetSplitableHadron();
259 } << 105 if(aNucleon) delete aNucleon;
260 <<
261 // Init target nucleus (assumed to be never <<
262 theParticipants.Init( aNucleus.GetA_asInt(), <<
263 <<
264 NumberOfProjectileSpectatorNucleons = std::a <<
265 NumberOfTargetSpectatorNucleons = aNucleus.G <<
266 NumberOfNNcollisions = 0; <<
267 <<
268 // reset/recalculate everything for the new <<
269 theParameters->InitForInteraction( theProjec <<
270 aNucleus. <<
271 <<
272 if ( theAdditionalString.size() != 0 ) { <<
273 std::for_each( theAdditionalString.begin() <<
274 DeleteVSplitableHadron() ); <<
275 } <<
276 theAdditionalString.clear(); <<
277 <<
278 #ifdef debugFTFmodel <<
279 G4cout << "FTF end of Init" << G4endl << G4e <<
280 #endif <<
281 <<
282 // In the case of Hydrogen target, for non-i <<
283 // do NOT simulate quasi-elastic (by forcing <<
284 // elastic scatering in theParameters - whic <<
285 // This is necessary because in this case qu <<
286 // with only one nucleon would be identical <<
287 // and the latter is already included in the <<
288 // (i.e. G4HadronElasticProcess). <<
289 if ( std::abs( theProjectile.GetDefinition() <<
290 aNucleus.GetA_asInt() < 2 ) theParamete <<
291 <<
292 if ( SampleBinInterval() ) theParticipants.S <<
293 } <<
294 <<
295 <<
296 //============================================ <<
297 <<
298 G4ExcitedStringVector* G4FTFModel::GetStrings( <<
299 <<
300 #ifdef debugFTFmodel <<
301 G4cout << "G4FTFModel::GetStrings() " << G4e <<
302 #endif <<
303 <<
304 G4ExcitedStringVector* theStrings = new G4Ex <<
305 theParticipants.GetList( theProjectile, theP <<
306 <<
307 SetImpactParameter( theParticipants.GetImpac <<
308 <<
309 StoreInvolvedNucleon(); <<
310 <<
311 G4bool Success( true ); <<
312 <<
313 if ( HighEnergyInter ) { <<
314 ReggeonCascade(); <<
315 <<
316 #ifdef debugFTFmodel <<
317 G4cout << "FTF PutOnMassShell " << G4endl; <<
318 #endif <<
319 <<
320 Success = PutOnMassShell(); <<
321 <<
322 #ifdef debugFTFmodel <<
323 G4cout << "FTF PutOnMassShell Success? " <<
324 #endif <<
325 <<
326 } <<
327 <<
328 #ifdef debugFTFmodel <<
329 G4cout << "FTF ExciteParticipants " << G4end <<
330 #endif <<
331 <<
332 if ( Success ) Success = ExciteParticipants( <<
333 <<
334 #ifdef debugFTFmodel <<
335 G4cout << "FTF ExciteParticipants Success? " <<
336 #endif <<
337 <<
338 if ( Success ) { <<
339 <<
340 #ifdef debugFTFmodel <<
341 G4cout << "FTF BuildStrings "; <<
342 #endif <<
343 <<
344 BuildStrings( theStrings ); <<
345 <<
346 #ifdef debugFTFmodel <<
347 G4cout << "FTF BuildStrings " << theString <<
348 << "FTF GetResiduals of Nuclei " << <<
349 #endif <<
350 <<
351 GetResiduals(); <<
352 <<
353 /* <<
354 if ( theParameters != 0 ) { <<
355 delete theParameters; <<
356 theParameters = 0; <<
357 } <<
358 */ <<
359 } else if ( ! GetProjectileNucleus() ) { <<
360 // Erase the hadron projectile <<
361 std::vector< G4VSplitableHadron* > primari <<
362 theParticipants.StartLoop(); <<
363 while ( theParticipants.Next() ) { /* Loo <<
364 const G4InteractionContent& interaction <<
365 // Do not allow for duplicates <<
366 if ( primaries.end() == <<
367 std::find( primaries.begin(), prima <<
368 primaries.push_back( interaction.GetPr <<
369 } <<
370 } 106 }
371 std::for_each( primaries.begin(), primarie << 107 }
372 primaries.clear(); << 108 */
373 } <<
374 <<
375 // Cleaning of the memory <<
376 G4VSplitableHadron* aNucleon = 0; <<
377 <<
378 // Erase the projectile nucleons <<
379 for ( G4int i = 0; i < NumberOfInvolvedNucle <<
380 aNucleon = TheInvolvedNucleonsOfProjectile <<
381 if ( aNucleon ) delete aNucleon; <<
382 } <<
383 NumberOfInvolvedNucleonsOfProjectile = 0; <<
384 <<
385 // Erase the target nucleons <<
386 for ( G4int i = 0; i < NumberOfInvolvedNucle <<
387 aNucleon = TheInvolvedNucleonsOfTarget[i]- <<
388 if ( aNucleon ) delete aNucleon; <<
389 } <<
390 NumberOfInvolvedNucleonsOfTarget = 0; <<
391 <<
392 #ifdef debugFTFmodel <<
393 G4cout << "End of FTF. Go to fragmentation" <<
394 << "To continue - enter 1, to stop - <<
395 #endif <<
396 <<
397 theParticipants.Clean(); <<
398 <<
399 return theStrings; <<
400 } 109 }
401 110
>> 111 // ------------------------------------------------------------
>> 112 void G4FTFModel::Init(const G4Nucleus & aNucleus, const G4DynamicParticle & aProjectile)
>> 113 {
>> 114 theProjectile = aProjectile;
402 115
403 //============================================ << 116 G4double PlabPerParticle(0.); // Laboratory momentum Pz per particle/nucleon
404 <<
405 void G4FTFModel::StoreInvolvedNucleon() { <<
406 //To store nucleons involved in the interact <<
407 <<
408 NumberOfInvolvedNucleonsOfTarget = 0; <<
409 <<
410 G4V3DNucleus* theTargetNucleus = GetTargetNu <<
411 theTargetNucleus->StartLoop(); <<
412 <<
413 G4Nucleon* aNucleon; <<
414 while ( ( aNucleon = theTargetNucleus->GetNe <<
415 if ( aNucleon->AreYouHit() ) { <<
416 TheInvolvedNucleonsOfTarget[NumberOfInvo <<
417 NumberOfInvolvedNucleonsOfTarget++; <<
418 } <<
419 } <<
420 <<
421 #ifdef debugFTFmodel <<
422 G4cout << "G4FTFModel::StoreInvolvedNucleon <<
423 G4cout << "NumberOfInvolvedNucleonsOfTarget <<
424 << G4endl << G4endl; <<
425 #endif <<
426 <<
427 <<
428 if ( ! GetProjectileNucleus() ) return; // <<
429 <<
430 // The projectile is a nucleus or an anti-nu <<
431 <<
432 NumberOfInvolvedNucleonsOfProjectile = 0; <<
433 <<
434 G4V3DNucleus* theProjectileNucleus = GetProj <<
435 theProjectileNucleus->StartLoop(); <<
436 <<
437 G4Nucleon* aProjectileNucleon; <<
438 while ( ( aProjectileNucleon = theProjectile <<
439 if ( aProjectileNucleon->AreYouHit() ) { <<
440 // Projectile nucleon was involved in th <<
441 TheInvolvedNucleonsOfProjectile[NumberOf <<
442 NumberOfInvolvedNucleonsOfProjectile++; <<
443 } <<
444 } <<
445 <<
446 #ifdef debugFTFmodel <<
447 G4cout << "NumberOfInvolvedNucleonsOfProject <<
448 << G4endl << G4endl; <<
449 #endif <<
450 return; <<
451 } <<
452 117
>> 118 /*
>> 119 G4cout<<"FTF init Pro Name "<<theProjectile.GetDefinition()->GetParticleName()<<G4endl;
>> 120 G4cout<<"FTF init Pro Mass "<<theProjectile.GetMass()<<" "<<theProjectile.GetMomentum()<<G4endl;
>> 121 G4cout<<"FTF init Pro B Q "<<theProjectile.GetDefinition()->GetBaryonNumber()<<" "<<(G4int) theProjectile.GetDefinition()->GetPDGCharge()<<G4endl;
>> 122 G4cout<<"FTF init A Z "<<aNucleus.GetA_asInt()<<" "<<aNucleus.GetZ_asInt()<<G4endl;
>> 123 G4cout<<" "<<aNucleus.GetN()<<" "<<aNucleus.GetZ()<<G4endl;
>> 124 //G4int Uzhi; G4cin>>Uzhi;
>> 125 */
453 126
454 //============================================ << 127 theParticipants.SetProjectileNucleus(0);
>> 128 theParticipants.Init(aNucleus.GetA_asInt(),aNucleus.GetZ_asInt());
455 129
456 void G4FTFModel::ReggeonCascade() { << 130 if(std::abs(theProjectile.GetDefinition()->GetBaryonNumber()) <= 1)
457 // Implementation of the reggeon theory insp << 131 { // Projectile is a hadron
>> 132 PlabPerParticle=theProjectile.GetMomentum().z();
458 133
459 #ifdef debugReggeonCascade << 134 // S = sqr( theProjectile.GetMass() ) + sqr( ProtonMass ) +
460 G4cout << "G4FTFModel::ReggeonCascade ------ << 135 // 2*ProtonMass*theProjectile.GetTotalEnergy();
461 << "theProjectile.GetTotalMomentum() << 136 }
462 << "theProjectile.GetTotalEnergy() " <<
463 << "ExcitationE/WN " << theParameters <<
464 #endif <<
465 137
466 G4int InitNINt = NumberOfInvolvedNucleonsOfT <<
467 138
468 // Reggeon cascading in target nucleus << 139 if(theProjectile.GetDefinition()->GetBaryonNumber() > 1)
469 for ( G4int InvTN = 0; InvTN < InitNINt; Inv << 140 { // Projectile is a nucleus
470 G4Nucleon* aTargetNucleon = TheInvolvedNuc << 141 theParticipants.InitProjectileNucleus(
>> 142 theProjectile.GetDefinition()->GetBaryonNumber(),
>> 143 (G4int) theProjectile.GetDefinition()->GetPDGCharge() );
471 144
472 G4double CreationTime = aTargetNucleon->Ge << 145 G4ThreeVector BoostVector=theProjectile.GetMomentum()/theProjectile.GetTotalEnergy();
>> 146 theParticipants.theProjectileNucleus->DoLorentzBoost(BoostVector);
473 147
474 G4double XofWoundedNucleon = aTargetNucleo << 148 PlabPerParticle=theProjectile.GetMomentum().z()/
475 G4double YofWoundedNucleon = aTargetNucleo << 149 theProjectile.GetDefinition()->GetBaryonNumber();
476 <<
477 G4V3DNucleus* theTargetNucleus = GetTarget <<
478 theTargetNucleus->StartLoop(); <<
479 150
480 G4Nucleon* Neighbour(0); << 151 // S = 2.*sqr( ProtonMass ) + 2*ProtonMass*
481 while ( ( Neighbour = theTargetNucleus->Ge << 152 // theProjectile.GetTotalEnergy()/theProjectile.GetDefinition()->GetBaryonNumber();
482 if ( ! Neighbour->AreYouHit() ) { <<
483 G4double impact2 = sqr( XofWoundedNucl <<
484 sqr( YofWoundedNucl <<
485 <<
486 if ( G4UniformRand() < theParameters-> <<
487 G4Exp( -impact2 <<
488 ) { <<
489 // The neighbour nucleon is involved <<
490 TheInvolvedNucleonsOfTarget[ NumberO <<
491 NumberOfInvolvedNucleonsOfTarget++; <<
492 <<
493 G4VSplitableHadron* targetSplitable; <<
494 targetSplitable = new G4DiffractiveS <<
495 <<
496 Neighbour->Hit( targetSplitable ); <<
497 targetSplitable->SetTimeOfCreation( <<
498 targetSplitable->SetStatus( 3 ); // <<
499 } 153 }
500 } <<
501 } <<
502 } <<
503 154
504 #ifdef debugReggeonCascade << 155 if(theProjectile.GetDefinition()->GetBaryonNumber() < -1)
505 G4cout << "Final NumberOfInvolvedNucleonsOfT << 156 { // Projectile is an anti-nucleus
506 << NumberOfInvolvedNucleonsOfTarget < << 157 theParticipants.InitProjectileNucleus(
507 #endif << 158 std::abs( theProjectile.GetDefinition()->GetBaryonNumber()),
>> 159 std::abs((G4int) theProjectile.GetDefinition()->GetPDGCharge()) );
508 160
509 if ( ! GetProjectileNucleus() ) return; << 161 G4ThreeVector BoostVector=theProjectile.GetMomentum()/theProjectile.GetTotalEnergy();
510 162
511 // Nucleus-Nucleus Interaction : Destruction << 163 theParticipants.theProjectileNucleus->StartLoop();
512 G4int InitNINp = NumberOfInvolvedNucleonsOfP << 164 G4Nucleon * aNucleon;
>> 165 while ( (aNucleon = theParticipants.theProjectileNucleus->GetNextNucleon()) )
>> 166 {
>> 167 if(aNucleon->GetDefinition()->GetPDGEncoding() == 2212) // Proton
>> 168 {aNucleon->SetParticleType(G4AntiProton::AntiProton());}
513 169
514 //for ( G4int InvPN = 0; InvPN < NumberOfInv << 170 if(aNucleon->GetDefinition()->GetPDGEncoding() == 2112) // Neutron
515 for ( G4int InvPN = 0; InvPN < InitNINp; Inv << 171 {aNucleon->SetParticleType(G4AntiNeutron::AntiNeutron());}
516 G4Nucleon* aProjectileNucleon = TheInvolve << 172 } // end of while (theParticipant.theProjectileNucleus->GetNextNucleon())
517 173
518 G4double CreationTime = aProjectileNucleon << 174 theParticipants.theProjectileNucleus->DoLorentzBoost(BoostVector);
519 175
520 G4double XofWoundedNucleon = aProjectileNu << 176 PlabPerParticle= theProjectile.GetMomentum().z()/
521 G4double YofWoundedNucleon = aProjectileNu << 177 std::abs(theProjectile.GetDefinition()->GetBaryonNumber());
522 <<
523 G4V3DNucleus* theProjectileNucleus = GetPr <<
524 theProjectileNucleus->StartLoop(); <<
525 178
526 G4Nucleon* Neighbour( 0 ); << 179 // S = 2.*sqr( ProtonMass ) + 2*ProtonMass*
527 while ( ( Neighbour = theProjectileNucleus << 180 // theProjectile.GetTotalEnergy()/
528 if ( ! Neighbour->AreYouHit() ) { << 181 // std::abs(theProjectile.GetDefinition()->GetBaryonNumber());
529 G4double impact2= sqr( XofWoundedNucle <<
530 sqr( YofWoundedNucle <<
531 <<
532 if ( G4UniformRand() < theParameters-> <<
533 G4Exp( -impact2 <<
534 ) { <<
535 // The neighbour nucleon is involved <<
536 TheInvolvedNucleonsOfProjectile[ Num <<
537 NumberOfInvolvedNucleonsOfProjectile <<
538 <<
539 G4VSplitableHadron* projectileSplita <<
540 projectileSplitable = new G4Diffract <<
541 <<
542 Neighbour->Hit( projectileSplitable <<
543 projectileSplitable->SetTimeOfCreati <<
544 projectileSplitable->SetStatus( 3 ); <<
545 } 182 }
546 } <<
547 } <<
548 } <<
549 183
550 #ifdef debugReggeonCascade << 184 // ------------------------------------------------------------------------
551 G4cout << "NumberOfInvolvedNucleonsOfProject << 185 if( theParameters != 0 ) delete theParameters;
552 << NumberOfInvolvedNucleonsOfProjecti << 186 theParameters = new G4FTFParameters(theProjectile.GetDefinition(),
553 #endif << 187 aNucleus.GetA_asInt(),aNucleus.GetZ_asInt(),
554 } << 188 PlabPerParticle);
555 << 189 //G4cout<<" End Init "<<theProjectile.GetMomentum()<<G4endl;
556 << 190 // To turn on/off (1/0) elastic scattering close/open ...
557 //============================================ << 191 //theParameters->SetProbabilityOfElasticScatt(0.);
558 << 192 //G4cout<<" etProbabilityOfElasticScatt "<<theParameters->GetProbabilityOfElasticScatt()<<G4endl;
559 G4bool G4FTFModel::PutOnMassShell() { << 193 //G4cout<<" INIT ";
560 << 194 //G4int Uzhi; G4cin>>Uzhi;
561 G4bool isProjectileNucleus = false; <<
562 if ( GetProjectileNucleus() ) isProjectileNu <<
563 <<
564 #ifdef debugPutOnMassShell <<
565 G4cout << "PutOnMassShell start " << G4endl; <<
566 if ( isProjectileNucleus ) { <<
567 G4cout << "PutOnMassShell for Nucleus_Nucl <<
568 } <<
569 #endif <<
570 <<
571 G4LorentzVector Pprojectile( theProjectile.G <<
572 if ( Pprojectile.z() < 0.0 ) return false; <<
573 <<
574 G4bool isOk = true; <<
575 <<
576 G4LorentzVector Ptarget( 0.0, 0.0, 0.0, 0.0 <<
577 G4LorentzVector PtargetResidual( 0.0, 0.0, 0 <<
578 G4double SumMasses = 0.0; <<
579 G4V3DNucleus* theTargetNucleus = GetTargetNu <<
580 G4double TargetResidualMass = 0.0; <<
581 <<
582 #ifdef debugPutOnMassShell <<
583 G4cout << "Target : "; <<
584 #endif <<
585 isOk = ComputeNucleusProperties( theTargetNu <<
586 TargetResid <<
587 TargetResid <<
588 if ( ! isOk ) return false; <<
589 <<
590 G4double Mprojectile = 0.0; <<
591 G4double M2projectile = 0.0; <<
592 G4LorentzVector Pproj( 0.0, 0.0, 0.0, 0.0 ); <<
593 G4LorentzVector PprojResidual( 0.0, 0.0, 0.0 <<
594 G4V3DNucleus* thePrNucleus = GetProjectileNu <<
595 G4double PrResidualMass = 0.0; <<
596 <<
597 if ( ! isProjectileNucleus ) { // hadron-nu <<
598 Mprojectile = Pprojectile.mag(); <<
599 M2projectile = Pprojectile.mag2(); <<
600 SumMasses += Mprojectile + 20.0*MeV; <<
601 } else { // nucleus-nucleus or antinucleus- <<
602 #ifdef debugPutOnMassShell <<
603 G4cout << "Projectile : "; <<
604 #endif <<
605 isOk = ComputeNucleusProperties( thePrNucl <<
606 Projectil <<
607 Projectil <<
608 if ( ! isOk ) return false; <<
609 } <<
610 <<
611 G4LorentzVector Psum = Pprojectile + Ptarget <<
612 G4double SqrtS = Psum.mag(); <<
613 G4double S = Psum.mag2(); <<
614 <<
615 #ifdef debugPutOnMassShell <<
616 G4cout << "Psum " << Psum/GeV << " GeV" << G <<
617 << "SumMasses, PrResidualMass and Tar <<
618 << PrResidualMass/GeV << " " << Targe <<
619 #endif <<
620 <<
621 if ( SqrtS < SumMasses ) return false; // I <<
622 <<
623 // Try to consider also the excitation energ <<
624 // possible, with the available energy; othe <<
625 G4double savedSumMasses = SumMasses; <<
626 if ( isProjectileNucleus ) { <<
627 SumMasses -= std::sqrt( sqr( PrResidualMas <<
628 SumMasses += std::sqrt( sqr( PrResidualMas <<
629 + PprojResidual.pe <<
630 } <<
631 SumMasses -= std::sqrt( sqr( TargetResidualM <<
632 SumMasses += std::sqrt( sqr( TargetResidualM <<
633 + PtargetResidual.pe <<
634 <<
635 if ( SqrtS < SumMasses ) { <<
636 SumMasses = savedSumMasses; <<
637 if ( isProjectileNucleus ) ProjectileResid <<
638 TargetResidualExcitationEnergy = 0.0; <<
639 } <<
640 <<
641 TargetResidualMass += TargetResidualExcitati <<
642 if ( isProjectileNucleus ) PrResidualMass += <<
643 <<
644 #ifdef debugPutOnMassShell <<
645 if ( isProjectileNucleus ) { <<
646 G4cout << "PrResidualMass ProjResidualExci <<
647 << ProjectileResidualExcitationEnergy << <<
648 } <<
649 G4cout << "TargetResidualMass TargetResidual <<
650 << TargetResidualExcitationEnergy << <<
651 << "Sum masses " << SumMasses/GeV << <<
652 #endif <<
653 <<
654 // Sampling of nucleons what can transfer to <<
655 if ( isProjectileNucleus && thePrNucleus-> <<
656 isOk = GenerateDeltaIsobar( SqrtS, Numbe <<
657 TheInvolvedN <<
658 } <<
659 if ( theTargetNucleus->GetMassNumber() != 1 <<
660 isOk = isOk && GenerateDeltaIsobar( Sqrt <<
661 TheI <<
662 } <<
663 if ( ! isOk ) return false; <<
664 <<
665 // Now we know that it is kinematically poss <<
666 // of the involved nucleons (or correspondin <<
667 // We have to sample the kinematical variabl <<
668 // of the final state. The sampled kinematic <<
669 // Notice that the sampling of the transvers <<
670 // Fermi motion. <<
671 <<
672 G4LorentzRotation toCms( -1*Psum.boostVector <<
673 G4LorentzVector Ptmp = toCms*Pprojectile; <<
674 if ( Ptmp.pz() <= 0.0 ) return false; // "S <<
675 <<
676 G4LorentzRotation toLab( toCms.inverse() ); <<
677 <<
678 G4double YprojectileNucleus = 0.0; <<
679 if ( isProjectileNucleus ) { <<
680 Ptmp = toCms*Pproj; <<
681 YprojectileNucleus = Ptmp.rapidity(); <<
682 } <<
683 Ptmp = toCms*Ptarget; <<
684 G4double YtargetNucleus = Ptmp.rapidity(); <<
685 <<
686 // Ascribing of the involved nucleons Pt and <<
687 G4double DcorP = 0.0; <<
688 if ( isProjectileNucleus ) DcorP = theParame <<
689 G4double DcorT = theParameters->GetDof <<
690 G4double AveragePt2 = theParameters->GetPt2 <<
691 G4double maxPtSquare = theParameters->GetMax <<
692 <<
693 #ifdef debugPutOnMassShell <<
694 if ( isProjectileNucleus ) { <<
695 G4cout << "Y projectileNucleus " << Yproje <<
696 } <<
697 G4cout << "Y targetNucleus " << YtargetN <<
698 << "Dcor " << theParameters->GetDofNu <<
699 << " DcorP DcorT " << DcorP << " " << <<
700 #endif <<
701 <<
702 G4double M2proj = M2projectile; // Initiali <<
703 G4double WplusProjectile = 0.0; <<
704 G4double M2target = 0.0; <<
705 G4double WminusTarget = 0.0; <<
706 G4int NumberOfTries = 0; <<
707 G4double ScaleFactor = 2.0; <<
708 G4bool OuterSuccess = true; <<
709 <<
710 const G4int maxNumberOfLoops = 1000; <<
711 G4int loopCounter = 0; <<
712 do { // while ( ! OuterSuccess ) <<
713 OuterSuccess = true; <<
714 const G4int maxNumberOfInnerLoops = 10000; <<
715 do { // while ( SqrtS < Mprojectile + std <<
716 NumberOfTries++; <<
717 if ( NumberOfTries == 100*(NumberOfTries <<
718 // After many tries, it is convenient <<
719 // AveragePt2, so that the sampled mom <<
720 // involved nucleons (or corresponding delta <<
721 // it is more likely to satisfy the mo <<
722 ScaleFactor /= 2.0; <<
723 DcorP *= ScaleFactor; <<
724 DcorT *= ScaleFactor; <<
725 AveragePt2 *= ScaleFactor; <<
726 } <<
727 if ( isProjectileNucleus ) { <<
728 // Sampling of kinematical properties <<
729 isOk = SamplingNucleonKinematics( Aver <<
730 theP <<
731 PrRe <<
732 Numb <<
733 TheI <<
734 } <<
735 // Sampling of kinematical properties of <<
736 isOk = isOk && SamplingNucleonKinemati <<
737 <<
738 <<
739 <<
740 <<
741 #ifdef debugPutOnMassShell <<
742 G4cout << "SqrtS, Mp+Mt, Mp, Mt " << Sqr <<
743 << ( std::sqrt( M2proj ) + std::s <<
744 << std::sqrt( M2proj )/GeV << " " <<
745 #endif <<
746 if ( ! isOk ) return false; <<
747 } while ( ( SqrtS < std::sqrt( M2proj ) + <<
748 NumberOfTries < maxNumberOfInner <<
749 if ( NumberOfTries >= maxNumberOfInnerLoop <<
750 #ifdef debugPutOnMassShell <<
751 G4cout << "BAD situation: forced exit of <<
752 #endif <<
753 return false; <<
754 } <<
755 if ( isProjectileNucleus ) { <<
756 isOk = CheckKinematics( S, SqrtS, M2proj <<
757 NumberOfInvolved <<
758 TheInvolvedNucle <<
759 WminusTarget, Wp <<
760 } <<
761 isOk = isOk && CheckKinematics( S, SqrtS <<
762 NumberOf <<
763 WminusTa <<
764 if ( ! isOk ) return false; <<
765 } while ( ( ! OuterSuccess ) && <<
766 ++loopCounter < maxNumberOfLoops ) <<
767 if ( loopCounter >= maxNumberOfLoops ) { <<
768 #ifdef debugPutOnMassShell <<
769 G4cout << "BAD situation: forced exit of t <<
770 #endif <<
771 return false; <<
772 } <<
773 <<
774 // Now the sampling is completed, and we can <<
775 // whole system. This is done first in the c <<
776 // to the lab frame. The transverse momentum <<
777 // the recoil of each hadron (nucleon or del <<
778 // to conserve (by construction) the transve <<
779 <<
780 if ( ! isProjectileNucleus ) { // hadron-nu <<
781 <<
782 G4double Pzprojectile = WplusProjectile/2. <<
783 G4double Eprojectile = WplusProjectile/2. <<
784 Pprojectile.setPz( Pzprojectile ); <<
785 Pprojectile.setE( Eprojectile ); <<
786 <<
787 #ifdef debugPutOnMassShell <<
788 G4cout << "Proj after in CMS " << Pproject <<
789 #endif <<
790 <<
791 Pprojectile.transform( toLab ); <<
792 theProjectile.SetMomentum( Pprojectile.vec <<
793 theProjectile.SetTotalEnergy( Pprojectile. <<
794 <<
795 theParticipants.StartLoop(); <<
796 theParticipants.Next(); <<
797 G4VSplitableHadron* primary = theParticipa <<
798 primary->Set4Momentum( Pprojectile ); <<
799 <<
800 #ifdef debugPutOnMassShell <<
801 G4cout << "Final proj. mom in Lab. " << pr <<
802 #endif <<
803 <<
804 } else { // nucleus-nucleus or antinucleus- <<
805 <<
806 isOk = FinalizeKinematics( WplusProjectile <<
807 ProjectileResid <<
808 TheInvolvedNucl <<
809 <<
810 #ifdef debugPutOnMassShell <<
811 G4cout << "Projectile Residual4Momentum in <<
812 #endif <<
813 <<
814 if ( ! isOk ) return false; <<
815 <<
816 ProjectileResidual4Momentum.transform( toL <<
817 <<
818 #ifdef debugPutOnMassShell <<
819 G4cout << "Projectile Residual4Momentum in <<
820 #endif <<
821 <<
822 } <<
823 <<
824 isOk = FinalizeKinematics( WminusTarget, fal <<
825 TargetResidualMas <<
826 TheInvolvedNucleo <<
827 <<
828 #ifdef debugPutOnMassShell <<
829 G4cout << "Target Residual4Momentum in CMS " <<
830 #endif <<
831 <<
832 if ( ! isOk ) return false; <<
833 <<
834 TargetResidual4Momentum.transform( toLab ); <<
835 <<
836 #ifdef debugPutOnMassShell <<
837 G4cout << "Target Residual4Momentum in Lab " <<
838 #endif <<
839 <<
840 return true; <<
841 195
>> 196 if(theAdditionalString.size() != 0)
>> 197 {
>> 198 std::for_each(theAdditionalString.begin(), theAdditionalString.end(),
>> 199 DeleteVSplitableHadron());
>> 200 }
>> 201 theAdditionalString.clear();
>> 202 //G4cout<<" End Init theProjectile.GetMomentum()"<<theProjectile.GetMomentum()<<G4endl;
842 } 203 }
843 204
844 << 205 // ------------------------------------------------------------
845 //============================================ << 206 G4ExcitedStringVector * G4FTFModel::GetStrings()
846 << 207 {
847 G4bool G4FTFModel::ExciteParticipants() { << 208 G4ExcitedStringVector * theStrings(0);
848 << 209
849 #ifdef debugBuildString << 210 theParticipants.GetList(theProjectile,theParameters);
850 G4cout << "G4FTFModel::ExciteParticipants() << 211 // StoreInvolvedNucleon();
851 #endif << 212 //G4cout<<" GetList theProjectile.GetMomentum() GetBaryonNumber() "<<theProjectile.GetMomentum()<<" "<<theProjectile.GetDefinition()->GetBaryonNumber()<<G4endl;
852 << 213 G4bool Success(true);
853 G4bool Success( false ); << 214
854 G4int MaxNumOfInelCollisions = G4int( thePar << 215 if((std::abs(theProjectile.GetDefinition()->GetBaryonNumber()) <= 1) &&
855 if ( MaxNumOfInelCollisions > 0 ) { // Pla << 216 (theProjectile.GetDefinition()->GetBaryonNumber() != -1) )
856 G4double ProbMaxNumber = theParameters->Ge << 217 { // Standard variant of FTF for projectile hadron/nucleon
857 if ( G4UniformRand() < ProbMaxNumber ) Max << 218 //G4cout<<"Standard variant of FTF for projectile hadron/nucleon"<<G4endl;
858 } else { << 219 ReggeonCascade();
859 // Plab < Pbound, normal application of FT << 220 //G4cout<<"Success after Reggeon "<<Success<<" PutOnMasShell"<<G4endl;
860 MaxNumOfInelCollisions = 1; << 221 Success=PutOnMassShell();
861 } << 222 //G4cout<<"Success after PutOn "<<Success<<" GetResid"<<G4endl;
862 << 223 GetResidualNucleus();
863 #ifdef debugBuildString <<
864 G4cout << "MaxNumOfInelCollisions per hadron <<
865 #endif <<
866 <<
867 G4int CurrentInteraction( 0 ); <<
868 theParticipants.StartLoop(); <<
869 <<
870 G4bool InnerSuccess( true ); <<
871 while ( theParticipants.Next() ) { /* Loop <<
872 CurrentInteraction++; <<
873 const G4InteractionContent& collision = th <<
874 G4VSplitableHadron* projectile = collision <<
875 G4Nucleon* ProjectileNucleon = collision.G <<
876 G4VSplitableHadron* target = collision.Get <<
877 G4Nucleon* TargetNucleon = collision.GetTa <<
878 <<
879 #ifdef debugBuildString <<
880 G4cout << G4endl << "Interaction # Status <<
881 << collision.GetStatus() << G4endl <<
882 << target << G4endl << "projectile- <<
883 << projectile->GetStatus() << " " < <<
884 << "projectile->GetSoftC target->G <<
885 << " " << target->GetSoftCollisionC <<
886 #endif <<
887 <<
888 if ( collision.GetStatus() ) { <<
889 if ( G4UniformRand() < theParameters->Ge <<
890 // Elastic scattering <<
891 <<
892 #ifdef debugBuildString <<
893 G4cout << "Elastic scattering" << G4en <<
894 #endif <<
895 <<
896 if ( ! HighEnergyInter ) { <<
897 G4bool Annihilation = false; <<
898 G4bool Result = AdjustNucleons( proj <<
899 Targ <<
900 if ( ! Result ) continue; <<
901 } <<
902 InnerSuccess = theElastic->ElasticSca <<
903 } else if ( G4UniformRand() > theParamet <<
904 // Inelastic scattering <<
905 <<
906 #ifdef debugBuildString <<
907 G4cout << "Inelastic interaction" << G <<
908 << "MaxNumOfInelCollisions per <<
909 #endif <<
910 <<
911 if ( ! HighEnergyInter ) { <<
912 G4bool Annihilation = false; <<
913 G4bool Result = AdjustNucleons( proj <<
914 Targ <<
915 if ( ! Result ) continue; <<
916 } 224 }
917 if ( G4UniformRand() < << 225 //G4cout<<"Success after GetN "<<Success<<G4endl;
918 ( 1.0 - target->GetSoftCollisionC << 226 //G4int Uzhi; G4cin>>Uzhi;
919 ( 1.0 - projectile->GetSoftCollis << 227 if(theProjectile.GetDefinition()->GetBaryonNumber() > 1)
920 //if ( ! HighEnergyInter ) { << 228 { // Variant of FTF for projectile nuclei
921 // G4bool Annihilation = false; << 229 //G4cout<<"Variant of FTF for projectile nuclei"<<G4endl;
922 // G4bool Result = AdjustNucleons( << 230 StoreInvolvedNucleon();
923 // << 231 ReggeonCascade();
924 // if ( ! Result ) continue; << 232 Success=PutOnMassShell();
925 //} << 233 GetResidualNucleus();
926 if ( theExcitation->ExciteParticipan <<
927 InnerSuccess = true; <<
928 NumberOfNNcollisions++; <<
929 #ifdef debugBuildString <<
930 G4cout << "FTF excitation Successf <<
931 // G4cout << "After pro " << proj <<
932 // << projectile->Get4Momen <<
933 // << "After tar " << targ <<
934 // << target->Get4Momentum( <<
935 #endif <<
936 } else { <<
937 InnerSuccess = theElastic->Elastic <<
938 #ifdef debugBuildString <<
939 G4cout << "FTF excitation Non Inne <<
940 << InnerSuccess << G4endl; <<
941 #endif <<
942 } <<
943 } else { // The inelastic interactiti <<
944 #ifdef debugBuildString <<
945 G4cout << "Elastic scat. at rejectio <<
946 #endif <<
947 //if ( ! HighEnergyInter ) { <<
948 // G4bool Annihilation = false; <<
949 // G4bool Result = AdjustNucleons( <<
950 // <<
951 // if ( ! Result) continue; <<
952 //} <<
953 InnerSuccess = theElastic->ElasticSc <<
954 } 234 }
955 } else { // Annihilation <<
956 <<
957 #ifdef debugBuildString <<
958 G4cout << "Annihilation" << G4endl; <<
959 #endif <<
960 <<
961 // At last, annihilation <<
962 if ( ! HighEnergyInter ) { <<
963 G4bool Annihilation = true; <<
964 G4bool Result = AdjustNucleons( proj <<
965 Targ <<
966 if ( ! Result ) continue; <<
967 } <<
968 235
969 G4VSplitableHadron* AdditionalString = << 236 // G4bool LowE_Anti_Ion(false);
970 if ( theAnnihilation->Annihilate( proj << 237 if(theProjectile.GetDefinition()->GetBaryonNumber() <= -1)
971 InnerSuccess = true; << 238 { // Projectile is Anti-baryon or Anti-Nucleus
972 #ifdef debugBuildString << 239 //G4cout<<"Projectile is Anti-baryon or Anti-Nucleus "<<G4endl;
973 G4cout << "Annihilation successfull. << 240 //G4cout<<"Be4 Store"<<G4endl;
974 << AdditionalString << G4endl << 241 StoreInvolvedNucleon();
975 //G4cout << "After pro " << project << 242 if(theProjectile.GetTotalMomentum()/
976 //G4cout << "After tar " << target- << 243 std::abs(theProjectile.GetDefinition()->GetBaryonNumber()) > 5000.*MeV)
977 #endif << 244 {// High energy interaction
978 << 245 //G4cout<<"High energy interaction "<<G4endl;
979 if ( AdditionalString != 0 ) theAddi << 246 //G4cout<<"Regeon "<<G4endl;
980 << 247 ReggeonCascade();
981 NumberOfNNcollisions++; << 248 //G4cout<<"Put on mass "<<G4endl;
982 << 249 Success=PutOnMassShell();
983 // Skipping possible interactions of << 250 //G4cout<<"Residual "<<G4endl;
984 while ( theParticipants.Next() ) { << 251 GetResidualNucleus();
985 G4InteractionContent& acollision = << 252 }
986 G4VSplitableHadron* NextProjectile << 253 else
987 G4VSplitableHadron* NextTargetNucl << 254 {
988 if ( projectile == NextProjectileN << 255 //G4cout<<"Low energy interaction "<<G4endl;
989 acollision.SetStatus( 0 ); << 256 // LowE_Anti_Ion=true;
990 } << 257 Success=true;
>> 258 }
>> 259 }
>> 260 //G4cout<<"Before Excite Success "<<Success<<G4endl;
>> 261 Success=Success && ExciteParticipants();
>> 262 //G4cout<<"Success ExciteParticipants()? "<<Success<<G4endl;
>> 263 // if(LowE_Anti_Ion) Success=Success && GetResidualNucleusAfterAnnihilation();
>> 264
>> 265 if( Success )
>> 266 {
>> 267 theStrings = BuildStrings();
>> 268 //G4cout<<"BuildStrings OK"<<G4endl;
>> 269 if( theParameters != 0 )
>> 270 {
>> 271 delete theParameters;
>> 272 theParameters=0;
991 } 273 }
992 << 274 }
993 // Continue the interactions << 275 /*
994 theParticipants.StartLoop(); << 276 if( Success )
995 for ( G4int i = 0; i < CurrentIntera << 277 {
996 << 278 if( ExciteParticipants() )
997 /* << 279 {
998 if ( target->GetStatus() == 4 ) { << 280 //G4cout<<"Excite partic OK"<<G4endl;
999 // Skipping possible interactions << 281 theStrings = BuildStrings();
1000 while ( theParticipants.Next() ) << 282 //G4cout<<"Build String OK"<<G4endl;
1001 G4InteractionContent& acollisio << 283 if(LowE_Anti_Ion) GetResidualNucleusAfterAnnihilation();
1002 G4VSplitableHadron* NextProject << 284
1003 G4VSplitableHadron* NextTargetN << 285 if( theParameters != 0 )
1004 if ( target == NextTargetNucleo << 286 {
1005 } << 287 delete theParameters;
>> 288 theParameters=0;
1006 } 289 }
1007 theParticipants.StartLoop(); << 290 } else // if( ExciteParticipants() )
1008 for ( G4int I = 0; I < CurrentInter << 291 { Success=false;}
1009 */ << 292 } else // if( Success )
1010 << 293 { Success=false;}
1011 } << 294 */
1012 } << 295 if(!Success)
1013 } << 296 {
1014 << 297 // -------------- Erase the projectile ----------------
1015 if( InnerSuccess ) Success = true; << 298 //G4cout<<"Erase Proj"<<G4endl;
1016 << 299 std::vector<G4VSplitableHadron *> primaries;
1017 #ifdef debugBuildString << 300
1018 G4cout << "----------------------------- << 301 theParticipants.StartLoop(); // restart a loop
1019 << "projectile->GetStatus target-> << 302 while ( theParticipants.Next() )
1020 << " " << target->GetStatus() << G << 303 {
1021 << projectile->GetSoftCollisionCou << 304 const G4InteractionContent & interaction=theParticipants.GetInteraction();
1022 << G4endl << "ExciteParticipants() << 305
1023 #endif << 306 // do not allow for duplicates ...
1024 << 307 if ( primaries.end() == std::find(primaries.begin(), primaries.end(),
1025 } // end of while ( theParticipants.Next() << 308 interaction.GetProjectile()) )
1026 << 309 primaries.push_back(interaction.GetProjectile());
1027 return Success; << 310 }
1028 } << 311 std::for_each(primaries.begin(), primaries.end(), DeleteVSplitableHadron());
>> 312 primaries.clear();
>> 313 }
>> 314 // -------------- Cleaning of the memory --------------
>> 315 G4VSplitableHadron * aNucleon = 0;
>> 316 // -------------- Erase the projectile nucleon --------
>> 317 /*
>> 318 G4VSplitableHadron * aNucleon = 0;
>> 319 for(G4int i=0; i < NumberOfInvolvedNucleonOfProjectile; i++)
>> 320 {
>> 321 aNucleon = TheInvolvedNucleonOfProjectile[i]->GetSplitableHadron();
>> 322 if(aNucleon) delete aNucleon;
>> 323 }
1029 324
>> 325 NumberOfInvolvedNucleonOfProjectile=0;
>> 326 */ // Maybe it will be needed latter------------------
1030 327
1031 //=========================================== << 328 // -------------- Erase the target nucleons -----------
>> 329 //G4cout<<"Erase Target Ninv "<<NumberOfInvolvedNucleon<<G4endl;
>> 330 for(G4int i=0; i < NumberOfInvolvedNucleon; i++)
>> 331 {
>> 332 aNucleon = TheInvolvedNucleon[i]->GetSplitableHadron();
>> 333 if(aNucleon) delete aNucleon;
>> 334 }
1032 335
1033 G4bool G4FTFModel::AdjustNucleons( G4VSplitab << 336 NumberOfInvolvedNucleon=0;
1034 G4Nucleon* << 337 //G4cout<<"Go to fragmentation"<<G4endl;
1035 G4VSplitab << 338 return theStrings;
1036 G4Nucleon* <<
1037 G4bool Ann <<
1038 <<
1039 #ifdef debugAdjust <<
1040 G4cout << "AdjustNucleons ----------------- <<
1041 << "Proj is nucleus? " << GetProject <<
1042 << "Proj 4mom " << SelectedAntiBaryo <<
1043 << "Targ 4mom " << SelectedTargetNuc <<
1044 << "Pr ResidualMassNumber Pr Residua <<
1045 << ProjectileResidualMassNumber << " <<
1046 << ProjectileResidualExcitationEnerg <<
1047 << "Tr ResidualMassNumber Tr Residua <<
1048 << TargetResidualMassNumber << " " < <<
1049 << TargetResidualExcitationEnergy << <<
1050 << "Collis. pr tr " << SelectedAntiB <<
1051 << SelectedTargetNucleon->GetSoftCol <<
1052 #endif <<
1053 <<
1054 if ( SelectedAntiBaryon->GetSoftCollisionCo <<
1055 SelectedTargetNucleon->GetSoftCollisio <<
1056 return true; // Selected hadrons were adj <<
1057 } <<
1058 <<
1059 G4int interactionCase = 0; <<
1060 if ( ( ! GetProjectileNucleus() && <<
1061 SelectedAntiBaryon->GetSoftCollis <<
1062 SelectedTargetNucleon->GetSoftCol <<
1063 || <<
1064 ( SelectedAntiBaryon->GetSoftCollis <<
1065 SelectedTargetNucleon->GetSoftCol <<
1066 // The case of hadron-nucleus interaction <<
1067 // the case when projectile nuclear nucle <<
1068 // a collision, but target nucleon did no <<
1069 interactionCase = 1; <<
1070 #ifdef debugAdjust <<
1071 G4cout << "case 1, hA prcol=0 trcol=0, AA <<
1072 #endif <<
1073 if ( TargetResidualMassNumber < 1 ) { <<
1074 return false; <<
1075 } <<
1076 if ( SelectedAntiBaryon->Get4Momentum().r <<
1077 return false; <<
1078 } <<
1079 if ( TargetResidualMassNumber == 1 ) { <<
1080 TargetResidualMassNumber = 0; <<
1081 TargetResidualCharge = 0; <<
1082 TargetResidualExcitationEnergy = 0.0; <<
1083 SelectedTargetNucleon->Set4Momentum( Ta <<
1084 TargetResidual4Momentum = G4LorentzVect <<
1085 return true; <<
1086 } <<
1087 } else if ( SelectedAntiBaryon->GetSoftColl <<
1088 SelectedTargetNucleon->GetSoftC <<
1089 // It is assumed that in the case there i <<
1090 interactionCase = 2; <<
1091 #ifdef debugAdjust <<
1092 G4cout << "case 2, prcol=0 trcol#0" << G <<
1093 #endif <<
1094 if ( ProjectileResidualMassNumber < 1 ) { <<
1095 return false; <<
1096 } <<
1097 if ( ProjectileResidual4Momentum.rapidity <<
1098 SelectedTargetNucleon->Get4Momentum( <<
1099 return false; <<
1100 } <<
1101 if ( ProjectileResidualMassNumber == 1 ) <<
1102 ProjectileResidualMassNumber = 0; <<
1103 ProjectileResidualCharge = 0; <<
1104 ProjectileResidualExcitationEnergy = 0. <<
1105 SelectedAntiBaryon->Set4Momentum( Proje <<
1106 ProjectileResidual4Momentum = G4Lorentz <<
1107 return true; <<
1108 } <<
1109 } else { // It has to be a nucleus-nucleus <<
1110 interactionCase = 3; <<
1111 #ifdef debugAdjust <<
1112 G4cout << "case 3, prcol=0 trcol=0" << G <<
1113 #endif <<
1114 if ( ! GetProjectileNucleus() ) { <<
1115 return false; <<
1116 } <<
1117 #ifdef debugAdjust <<
1118 G4cout << "Proj res Init " << ProjectileR <<
1119 << "Targ res Init " << TargetResid <<
1120 << "ProjectileResidualMassNumber P <<
1121 << ProjectileResidualMassNumber << <<
1122 << " (" << ProjectileResidualLambd <<
1123 << "TargetResidualMassNumber Targe <<
1124 << " " << TargetResidualCharge << <<
1125 #endif <<
1126 } <<
1127 <<
1128 CommonVariables common; <<
1129 G4int returnCode = AdjustNucleonsAlgorithm_ <<
1130 <<
1131 <<
1132 G4bool returnResult = false; <<
1133 if ( returnCode == 0 ) { <<
1134 returnResult = true; // Successfully end <<
1135 } else if ( returnCode == 1 ) { <<
1136 // The part before sampling has been succ <<
1137 returnResult = AdjustNucleonsAlgorithm_Sa <<
1138 if ( returnResult ) { // The sampling ha <<
1139 AdjustNucleonsAlgorithm_afterSampling( <<
1140 <<
1141 } <<
1142 } <<
1143 339
1144 return returnResult; <<
1145 } 340 }
1146 341
1147 //------------------------------------------- 342 //-------------------------------------------------------------------
>> 343 void G4FTFModel::StoreInvolvedNucleon()
>> 344 { //--- To store nucleons involved in low energy interaction -------
>> 345 if(std::abs(theProjectile.GetDefinition()->GetBaryonNumber()) <= 1)
>> 346 { // the projectile is a hadron -----------
>> 347 //G4cout<<"the projectile is a hadron"<<G4endl;
>> 348 NumberOfInvolvedNucleon=0;
>> 349
>> 350 theParticipants.StartLoop();
>> 351
>> 352 while (theParticipants.Next())
>> 353 {
>> 354 //G4cout<<"theParticipants.Next()"<<G4endl;
>> 355 const G4InteractionContent & collision=theParticipants.GetInteraction();
>> 356 //G4cout<<"collision=theParticipants.GetInteraction()"<<G4endl;
>> 357 G4Nucleon * TargetNucleon=collision.GetTargetNucleon();
>> 358 //G4cout<<"TargetNucleon=collision.GetTargetNucleon()"<<G4endl;
>> 359
>> 360 TheInvolvedNucleon[NumberOfInvolvedNucleon]=TargetNucleon;
>> 361 NumberOfInvolvedNucleon++;
>> 362 //G4cout<<G4endl<<"Prim NumberOfInvolvedNucleon "<<NumberOfInvolvedNucleon<<G4endl;
>> 363 } // end of while (theParticipants.Next())
>> 364
>> 365 NumberOfInvolvedTargetNucleon=NumberOfInvolvedNucleon;
>> 366 // ---------------- Calculation of creation time for each target nucleon -----------
>> 367 //G4cout<<"theParticipants.StartLoop() "<<G4endl;
>> 368 theParticipants.StartLoop(); // restart a loop
>> 369 //G4cout<<"theParticipants.Next();"<<G4endl;
>> 370 theParticipants.Next();
>> 371 G4VSplitableHadron * primary = theParticipants.GetInteraction().GetProjectile();
>> 372 //G4cout<<"primary->Get4Momentum() "<<primary->Get4Momentum()<<G4endl;
>> 373 //G4cout<<"primary->Get4Momentum().pz() "<<primary->Get4Momentum().pz()<<G4endl;
>> 374 //G4cout<<"primary->Get4Momentum().e() "<<primary->Get4Momentum().e()<<G4endl;
>> 375
>> 376 G4double betta_z=primary->Get4Momentum().pz()/primary->Get4Momentum().e();
>> 377 //G4cout<<"betta_z "<<betta_z<<G4endl;
>> 378 primary->SetTimeOfCreation(0.);
>> 379
>> 380 G4double ZcoordinateOfPreviousCollision(0.);
>> 381 G4double ZcoordinateOfCurrentInteraction(0.);
>> 382 G4double TimeOfPreviousCollision(0.);
>> 383 G4double TimeOfCurrentCollision(0);
>> 384
>> 385 theParticipants.theNucleus->StartLoop();
>> 386 G4Nucleon * aNucleon;
>> 387 G4bool theFirstInvolvedNucleon(true);
>> 388 while ( (aNucleon = theParticipants.theNucleus->GetNextNucleon()) )
>> 389 {
>> 390 //G4cout<<"aNucleon->AreYouHit() "<<aNucleon->AreYouHit()<<G4endl;
>> 391 if(aNucleon->AreYouHit())
>> 392 {
>> 393 if(theFirstInvolvedNucleon)
>> 394 {
>> 395 ZcoordinateOfPreviousCollision=aNucleon->GetPosition().z();
>> 396 //G4cout<<"ZcoordinateOfPreviousCollision "<<ZcoordinateOfPreviousCollision/fermi<<G4endl;
>> 397 theFirstInvolvedNucleon=false;
>> 398 }
1148 399
1149 G4int G4FTFModel::AdjustNucleonsAlgorithm_bef << 400 ZcoordinateOfCurrentInteraction=aNucleon->GetPosition().z();
1150 << 401 //G4cout<<"ZcoordinateOfCurrentInteraction "<<ZcoordinateOfCurrentInteraction/fermi<<G4endl;
1151 << 402 //G4cout<<"TimeOfPreviousCollision "<<TimeOfPreviousCollision<<G4endl;
1152 << 403
1153 << 404 // A.R. 18-Oct-2011 : Protection needed for nuclear capture of
1154 << 405 // anti-proton at rest.
1155 << 406 if ( betta_z > 1.0e-10 ) {
1156 // First of the three utility methods used << 407 TimeOfCurrentCollision=TimeOfPreviousCollision+
1157 // This method returns an integer code - in << 408 (ZcoordinateOfCurrentInteraction-ZcoordinateOfPreviousCollision)/betta_z;
1158 // "0" : successfully ended and nothing e << 409 } else {
1159 // "1" : successfully completed, but the << 410 TimeOfCurrentCollision=TimeOfPreviousCollision;
1160 // "99" : unsuccessfully ended, nothing el << 411 }
1161 G4int returnCode = 99; << 412
1162 << 413 //G4cout<<"TimeOfCurrentCollision "<<TimeOfCurrentCollision<<G4endl;
1163 G4double ExcitationEnergyPerWoundedNucleon << 414 // It is assumed that the nucleons are ordered on increasing z-coordinate ------------
1164 << 415 aNucleon->GetSplitableHadron()->SetTimeOfCreation(TimeOfCurrentCollision);
1165 // some checks and initializations << 416
1166 if ( interactionCase == 1 ) { << 417 ZcoordinateOfPreviousCollision=ZcoordinateOfCurrentInteraction;
1167 common.Psum = SelectedAntiBaryon->Get4Mom << 418 TimeOfPreviousCollision=TimeOfCurrentCollision;
1168 #ifdef debugAdjust << 419 } // end of if(aNucleon->AreYouHit())
1169 G4cout << "Targ res Init " << TargetResid << 420 } // end of while (theParticipant.theNucleus->GetNextNucleon())
1170 #endif << 421 //
1171 common.Pprojectile = SelectedAntiBaryon-> << 422 return;
1172 } else if ( interactionCase == 2 ) { << 423 } // end of if(std::abs(theProjectile.GetDefinition()->GetBaryonNumber()) <= 1)
1173 common.Psum = ProjectileResidual4Momentum <<
1174 common.Pprojectile = ProjectileResidual4M <<
1175 } else if ( interactionCase == 3 ) { <<
1176 common.Psum = ProjectileResidual4Momentum <<
1177 common.Pprojectile = ProjectileResidual4M <<
1178 } <<
1179 <<
1180 // transform momenta to cms and then rotate <<
1181 common.toCms = G4LorentzRotation( -1*common <<
1182 common.Ptmp = common.toCms * common.Pprojec <<
1183 common.toCms.rotateZ( -1*common.Ptmp.phi() <<
1184 common.toCms.rotateY( -1*common.Ptmp.theta( <<
1185 common.Pprojectile.transform( common.toCms <<
1186 common.toLab = common.toCms.inverse(); <<
1187 common.SqrtS = common.Psum.mag(); <<
1188 common.S = sqr( common.SqrtS ); <<
1189 <<
1190 // get properties of the target residual an <<
1191 G4bool Stopping = false; <<
1192 if ( interactionCase == 1 ) { <<
1193 common.TResidualMassNumber = TargetResidu <<
1194 common.TResidualCharge = TargetResidual <<
1195 - G4int( TargetN <<
1196 common.TResidualExcitationEnergy = Targ <<
1197 - Exci <<
1198 if ( common.TResidualMassNumber <= 1 ) { <<
1199 common.TResidualExcitationEnergy = 0.0; <<
1200 } <<
1201 if ( common.TResidualMassNumber != 0 ) { <<
1202 common.TResidualMass = G4ParticleTable: <<
1203 ->GetIonMass( co <<
1204 } <<
1205 common.TNucleonMass = TargetNucleon->GetD <<
1206 common.SumMasses = SelectedAntiBaryon-> <<
1207 + common.TResidualMass <<
1208 #ifdef debugAdjust <<
1209 G4cout << "Annihilation " << Annihilation <<
1210 #endif <<
1211 } else if ( interactionCase == 2 ) { <<
1212 common.Ptarget = common.toCms * SelectedT <<
1213 common.TResidualMassNumber = ProjectileRe <<
1214 common.TResidualCharge = ProjectileResi <<
1215 - std::abs( G4in <<
1216 common.TResidualExcitationEnergy = Proj <<
1217 - Exci <<
1218 if ( common.TResidualMassNumber <= 1 ) { <<
1219 common.TResidualExcitationEnergy = 0.0; <<
1220 } <<
1221 if ( common.TResidualMassNumber != 0 ) { <<
1222 common.TResidualMass = G4ParticleTable: <<
1223 ->GetIonMass( co <<
1224 } <<
1225 common.TNucleonMass = ProjectileNucleon-> <<
1226 common.SumMasses = SelectedTargetNucleo <<
1227 + common.TResidualMass <<
1228 #ifdef debugAdjust <<
1229 G4cout << "SelectedTN.mag() PNMass + PRes <<
1230 << SelectedTargetNucleon->Get4Mome <<
1231 << common.TNucleonMass << " " << c <<
1232 #endif <<
1233 } else if ( interactionCase == 3 ) { <<
1234 common.Ptarget = common.toCms * TargetRes <<
1235 common.PResidualMassNumber = ProjectileRe <<
1236 common.PResidualCharge = ProjectileResi <<
1237 - std::abs( G4in <<
1238 common.PResidualLambdaNumber = Projectile <<
1239 if ( ProjectileNucleon->GetDefinition() = <<
1240 ProjectileNucleon->GetDefinition() == G4An <<
1241 --common.PResidualLambdaNumber; <<
1242 } <<
1243 common.PResidualExcitationEnergy = Proj <<
1244 - Exci <<
1245 if ( common.PResidualMassNumber <= 1 ) { <<
1246 common.PResidualExcitationEnergy = 0.0; <<
1247 } <<
1248 if ( common.PResidualMassNumber != 0 ) { <<
1249 if ( common.PResidualMassNumber == 1 ) <<
1250 if ( std::abs( common.PResidualCharge <<
1251 common.PResidualMass = G4Proton::De <<
1252 } else if ( common.PResidualLambdaNum <<
1253 common.PResidualMass = G4Lambda::De <<
1254 } else { <<
1255 common.PResidualMass = G4Neutron::D <<
1256 } <<
1257 } else { <<
1258 if ( common.PResidualLambdaNumber > 0 <<
1259 if ( common.PResidualMassNumber == <<
1260 common.PResidualMass = G4Lambda:: <<
1261 if ( std::abs( common.PResidualCharge ) <<
1262 common.PResidualMass += G4Proton::Def <<
1263 } else if ( common.PResidualLambdaNumbe <<
1264 common.PResidualMass += G4Neutron::De <<
1265 } else { <<
1266 common.PResidualMass += G4Lambda::Def <<
1267 } <<
1268 } else { <<
1269 common.PResidualMass = G4HyperNucleiPro <<
1270 std::abs( common.PRes <<
1271 common.PResidualLambdaN <<
1272 } <<
1273 } else { <<
1274 common.PResidualMass = G4ParticleTable::G <<
1275 GetIonMass( std::a <<
1276 } <<
1277 } <<
1278 } <<
1279 common.PNucleonMass = ProjectileNucleon-> <<
1280 common.TResidualMassNumber = TargetResidu <<
1281 common.TResidualCharge = TargetResidual <<
1282 - G4int( TargetN <<
1283 common.TResidualExcitationEnergy = Targ <<
1284 - Exci <<
1285 if ( common.TResidualMassNumber <= 1 ) { <<
1286 common.TResidualExcitationEnergy = 0.0; <<
1287 } <<
1288 if ( common.TResidualMassNumber != 0 ) { <<
1289 common.TResidualMass = G4ParticleTable: <<
1290 GetIonMass( comm <<
1291 } <<
1292 common.TNucleonMass = TargetNucleon->GetD <<
1293 common.SumMasses = common.PNucleonMass <<
1294 + common.TResidualMass <<
1295 #ifdef debugAdjust <<
1296 G4cout << "PNucleonMass PResidualMass TNu <<
1297 << " " << common.PResidualMass << <<
1298 << common.TResidualMass << G4endl <<
1299 << "PResidualExcitationEnergy " << <<
1300 << "TResidualExcitationEnergy " << <<
1301 #endif <<
1302 } // End-if on interactionCase <<
1303 <<
1304 if ( ! Annihilation ) { <<
1305 if ( common.SqrtS < common.SumMasses ) { <<
1306 #ifdef debugAdjust <<
1307 G4cout << "SqrtS < SumMasses " << commo <<
1308 #endif <<
1309 return returnCode; // Unsuccessfully e <<
1310 } <<
1311 if ( interactionCase == 1 || interactio <<
1312 if ( common.SqrtS < common.SumMasses + <<
1313 #ifdef debugAdjust <<
1314 G4cout << "TResidualExcitationEnergy <<
1315 #endif <<
1316 common.TResidualExcitationEnergy = co <<
1317 #ifdef debugAdjust <<
1318 G4cout << "TResidualExcitationEnergy <<
1319 #endif <<
1320 Stopping = true; <<
1321 return returnCode; // unsuccessfully <<
1322 } <<
1323 } else if ( interactionCase == 3 ) { <<
1324 #ifdef debugAdjust <<
1325 G4cout << "SqrtS < SumMasses + PResidua <<
1326 << common.SqrtS << " " << common <<
1327 << G4endl; <<
1328 #endif <<
1329 if ( common.SqrtS < common.SumMasses <<
1330 + common.TResidualE <<
1331 Stopping = true; <<
1332 if ( common.PResidualExcitationEnergy <<
1333 common.TResidualExcitationEnergy = <<
1334 } else if ( common.TResidualExcitatio <<
1335 common.PResidualExcitationEnergy = <<
1336 } else { <<
1337 G4double Fraction = ( common.SqrtS <<
1338 / ( common.PResidualExcitationEn <<
1339 common.PResidualExcitationEnergy *= <<
1340 common.TResidualExcitationEnergy *= <<
1341 } <<
1342 } <<
1343 } <<
1344 } // End-if on ! Annihilation <<
1345 <<
1346 if ( Annihilation ) { <<
1347 #ifdef debugAdjust <<
1348 G4cout << "SqrtS < SumMasses - TNucleonMa <<
1349 << common.SumMasses - common.TNucl <<
1350 #endif <<
1351 if ( common.SqrtS < common.SumMasses - co <<
1352 return returnCode; // unsuccessfully e <<
1353 } <<
1354 #ifdef debugAdjust <<
1355 G4cout << "SqrtS < SumMasses " << common. <<
1356 #endif <<
1357 if ( common.SqrtS < common.SumMasses ) { <<
1358 if ( interactionCase == 2 || interact <<
1359 common.TResidualExcitationEnergy = 0. <<
1360 } <<
1361 common.TNucleonMass = common.SqrtS - <<
1362 - common.TResidua <<
1363 #ifdef debugAdjust <<
1364 G4cout << "TNucleonMass " << common.TNu <<
1365 #endif <<
1366 common.SumMasses = common.SqrtS - commo <<
1367 Stopping = true; <<
1368 #ifdef debugAdjust <<
1369 G4cout << "SqrtS < SumMasses " << commo <<
1370 #endif <<
1371 } <<
1372 if ( interactionCase == 1 || interactio <<
1373 if ( common.SqrtS < common.SumMasses + <<
1374 common.TResidualExcitationEnergy = co <<
1375 Stopping = true; <<
1376 } <<
1377 } else if ( interactionCase == 3 ) { <<
1378 if ( common.SqrtS < common.SumMasses <<
1379 + common.TResidualE <<
1380 Stopping = true; <<
1381 if ( common.PResidualExcitationEnergy <<
1382 common.TResidualExcitationEnergy = <<
1383 } else if ( common.TResidualExcitatio <<
1384 common.PResidualExcitationEnergy = <<
1385 } else { <<
1386 G4double Fraction = ( common.SqrtS <<
1387 ( common.PResidualExcitationEnerg <<
1388 common.PResidualExcitationEnergy *= <<
1389 common.TResidualExcitationEnergy *= <<
1390 } <<
1391 } <<
1392 } <<
1393 } // End-if on Annihilation <<
1394 <<
1395 #ifdef debugAdjust <<
1396 G4cout << "Stopping " << Stopping << G4endl <<
1397 #endif <<
1398 <<
1399 if ( Stopping ) { <<
1400 // All 3-momenta of particles = 0 <<
1401 common.Ptmp.setPx( 0.0 ); common.Ptmp.set <<
1402 // New projectile <<
1403 if ( interactionCase == 1 ) { <<
1404 common.Ptmp.setE( SelectedAntiBaryon->G <<
1405 } else if ( interactionCase == 2 ) { <<
1406 common.Ptmp.setE( common.TNucleonMass ) <<
1407 } else if ( interactionCase == 3 ) { <<
1408 common.Ptmp.setE( common.PNucleonMass ) <<
1409 } <<
1410 #ifdef debugAdjust <<
1411 G4cout << "Proj stop " << common.Ptmp << <<
1412 #endif <<
1413 common.Pprojectile = common.Ptmp; <<
1414 common.Pprojectile.transform( common.toLa <<
1415 //---AR-Jul2019 : To avoid unphysical pro <<
1416 // original momentum of th <<
1417 G4LorentzVector saveSelectedAntiBaryon4Mo <<
1418 saveSelectedAntiBaryon4Momentum.transform <<
1419 //--- <<
1420 SelectedAntiBaryon->Set4Momentum( common. <<
1421 // New target nucleon <<
1422 if ( interactionCase == 1 || interactio <<
1423 common.Ptmp.setE( common.TNucleonMass ) <<
1424 } else if ( interactionCase == 2 ) { <<
1425 common.Ptmp.setE( SelectedTargetNucleon <<
1426 } <<
1427 #ifdef debugAdjust <<
1428 G4cout << "Targ stop " << common.Ptmp << <<
1429 #endif <<
1430 common.Ptarget = common.Ptmp; <<
1431 common.Ptarget.transform( common.toLab ); <<
1432 //---AR-Jul2019 : To avoid unphysical tar <<
1433 // momentum of the target <<
1434 G4LorentzVector saveSelectedTargetNucleon <<
1435 saveSelectedTargetNucleon4Momentum.transf <<
1436 //--- <<
1437 SelectedTargetNucleon->Set4Momentum( comm <<
1438 // New target residual <<
1439 if ( interactionCase == 1 || interactio <<
1440 common.Ptmp.setPx( 0.0 ); common.Ptmp.s <<
1441 TargetResidualMassNumber = common <<
1442 TargetResidualCharge = common <<
1443 TargetResidualExcitationEnergy = common <<
1444 //---AR-Jul2019 : To avoid unphysical t <<
1445 // original momentum of <<
1446 // This is a rough and s <<
1447 //common.Ptmp.setE( common.TResidualMas <<
1448 common.Ptmp.setPx( -saveSelectedTargetN <<
1449 common.Ptmp.setPy( -saveSelectedTargetN <<
1450 common.Ptmp.setPz( -saveSelectedTargetN <<
1451 common.Ptmp.setE( std::sqrt( sqr( commo <<
1452 //--- <<
1453 #ifdef debugAdjust <<
1454 G4cout << "Targ Resi stop " << common.P <<
1455 #endif <<
1456 common.Ptmp.transform( common.toLab ); <<
1457 TargetResidual4Momentum = common.Ptmp; <<
1458 } <<
1459 // New projectile residual <<
1460 if ( interactionCase == 2 || interactio <<
1461 common.Ptmp.setPx( 0.0 ); common.Ptmp.s <<
1462 if ( interactionCase == 2 ) { <<
1463 ProjectileResidualMassNumber = <<
1464 ProjectileResidualCharge = <<
1465 ProjectileResidualLambdaNumber = 0; // <<
1466 ProjectileResidualExcitationEnergy = <<
1467 common.Ptmp.setE( common.TResidualMas <<
1468 } else { <<
1469 ProjectileResidualMassNumber = <<
1470 ProjectileResidualCharge = <<
1471 ProjectileResidualLambdaNumber = common <<
1472 ProjectileResidualExcitationEnergy = <<
1473 //---AR-Jul2019 : To avoid unphysical <<
1474 // saved original mome <<
1475 // This is a rough and <<
1476 //common.Ptmp.setE( common.PResidualM <<
1477 common.Ptmp.setPx( -saveSelectedAntiB <<
1478 common.Ptmp.setPy( -saveSelectedAntiB <<
1479 common.Ptmp.setPz( -saveSelectedAntiB <<
1480 common.Ptmp.setE( std::sqrt( sqr( com <<
1481 //--- <<
1482 } <<
1483 #ifdef debugAdjust <<
1484 G4cout << "Proj Resi stop " << common.P <<
1485 #endif <<
1486 common.Ptmp.transform( common.toLab ); <<
1487 ProjectileResidual4Momentum = common.Pt <<
1488 } <<
1489 return returnCode = 0; // successfully e <<
1490 } // End-if on Stopping <<
1491 <<
1492 // Initializations before sampling <<
1493 if ( interactionCase == 1 ) { <<
1494 common.Mprojectile = common.Pprojectile. <<
1495 common.M2projectile = common.Pprojectile. <<
1496 common.TResidual4Momentum = common.toCms <<
1497 common.YtargetNucleus = common.TResidual4 <<
1498 common.TResidualMass += common.TResidualE <<
1499 } else if ( interactionCase == 2 ) { <<
1500 common.Mtarget = common.Ptarget.mag(); <<
1501 common.M2target = common.Ptarget.mag2(); <<
1502 common.TResidual4Momentum = common.toCms <<
1503 common.YprojectileNucleus = common.TResid <<
1504 common.TResidualMass += common.TResidualE <<
1505 } else if ( interactionCase == 3 ) { <<
1506 common.PResidual4Momentum = common.toCms <<
1507 common.YprojectileNucleus = common.PResid <<
1508 common.TResidual4Momentum = common.toCms* <<
1509 common.YtargetNucleus = common.TResidual4 <<
1510 common.PResidualMass += common.PResidualE <<
1511 common.TResidualMass += common.TResidualE <<
1512 } <<
1513 #ifdef debugAdjust <<
1514 G4cout << "YprojectileNucleus " << common.Y <<
1515 #endif <<
1516 424
1517 return returnCode = 1; // successfully com << 425 // The projectile is a nucleus or an anti-nucleus
1518 } << 426 //G4cout<<"FTF The projectile is a nucleus or an anti-nucleus"<<G4endl;
>> 427 NumberOfInvolvedNucleonOfProjectile=0;
>> 428
>> 429 G4V3DNucleus * ProjectileNucleus =theParticipants.GetProjectileNucleus();
>> 430 ProjectileNucleus->StartLoop();
>> 431
>> 432 G4Nucleon * ProjectileNucleon;
>> 433 while ( (ProjectileNucleon=ProjectileNucleus->GetNextNucleon()) )
>> 434 {
>> 435 if ( ProjectileNucleon->AreYouHit() )
>> 436 { // Projectile nucleon was involved in the interaction.
>> 437 TheInvolvedNucleonOfProjectile[NumberOfInvolvedNucleonOfProjectile]=
>> 438 ProjectileNucleon;
>> 439 NumberOfInvolvedNucleonOfProjectile++;
>> 440 }
>> 441 } // End of while ( (ProjectileNucleon=ProjectileNucleus->GetNextNucleon()) )
>> 442
>> 443 //------------------
>> 444 NumberOfInvolvedNucleon=0;
>> 445
>> 446 G4V3DNucleus * TargetNucleus =theParticipants.GetWoundedNucleus();
>> 447 TargetNucleus->StartLoop();
>> 448
>> 449 G4Nucleon * TargetNucleon;
>> 450 while ( (TargetNucleon=TargetNucleus->GetNextNucleon()) )
>> 451 {
>> 452 if ( TargetNucleon->AreYouHit() )
>> 453 { // Target nucleon was involved in the interaction.
>> 454 TheInvolvedNucleon[NumberOfInvolvedNucleon]=TargetNucleon;
>> 455 NumberOfInvolvedNucleon++;
>> 456 }
>> 457 } // End of while ( (TargetNucleon=TargetNucleus->GetNextNucleon()) )
>> 458 //G4cout<<"Store inv "<<NumberOfInvolvedNucleonOfProjectile<<" "<<NumberOfInvolvedNucleon<<G4endl;
1519 459
>> 460 NumberOfInvolvedTargetNucleon=NumberOfInvolvedNucleon;
>> 461 return;
>> 462 } // Uzhi 10 Feb. 2011
1520 463
1521 //------------------------------------------- 464 //-------------------------------------------------------------------
>> 465 void G4FTFModel::ReggeonCascade()
>> 466 { //--- Implementation of the reggeon theory inspired model-------
>> 467 //G4cout<<"In reggeon"<<G4endl;
>> 468
>> 469 if(std::abs(theProjectile.GetDefinition()->GetBaryonNumber()) > 1) return;
>> 470 // For Nucleus-nucleus or Anti-nucleus - nucleus interactions
>> 471 // the cascading will be implemented latter.
>> 472
>> 473 NumberOfInvolvedNucleon=0;
>> 474
>> 475 theParticipants.StartLoop();
>> 476 while (theParticipants.Next())
>> 477 {
>> 478 const G4InteractionContent & collision=theParticipants.GetInteraction();
>> 479
>> 480 G4Nucleon * TargetNucleon=collision.GetTargetNucleon();
>> 481
>> 482 TheInvolvedNucleon[NumberOfInvolvedNucleon]=TargetNucleon;
>> 483 NumberOfInvolvedNucleon++;
>> 484 //G4cout<<"Prim NumberOfInvolvedNucleon "<<NumberOfInvolvedNucleon<<G4endl;
>> 485 G4double XofWoundedNucleon = TargetNucleon->GetPosition().x();
>> 486 G4double YofWoundedNucleon = TargetNucleon->GetPosition().y();
>> 487
>> 488 theParticipants.theNucleus->StartLoop();
>> 489 G4Nucleon * Neighbour(0);
>> 490
>> 491 while ( (Neighbour = theParticipants.theNucleus->GetNextNucleon()) )
>> 492 {
>> 493 if(!Neighbour->AreYouHit())
>> 494 {
>> 495 G4double impact2= sqr(XofWoundedNucleon - Neighbour->GetPosition().x()) +
>> 496 sqr(YofWoundedNucleon - Neighbour->GetPosition().y());
>> 497
>> 498 if(G4UniformRand() < theParameters->GetCofNuclearDestruction()*
>> 499 std::exp(-impact2/theParameters->GetR2ofNuclearDestruction()))
>> 500 { // The neighbour nucleon is involved in the reggeon cascade
>> 501
>> 502 TheInvolvedNucleon[NumberOfInvolvedNucleon]=Neighbour;
>> 503 NumberOfInvolvedNucleon++;
>> 504 //G4cout<<"Seco NumberOfInvolvedNucleon "<<NumberOfInvolvedNucleon<<G4endl;
>> 505
>> 506 G4VSplitableHadron *targetSplitable;
>> 507 targetSplitable = new G4DiffractiveSplitableHadron(*Neighbour);
>> 508
>> 509 Neighbour->Hit(targetSplitable);
>> 510 targetSplitable->SetStatus(2);
>> 511 }
>> 512 } // end of if(!Neighbour->AreYouHit())
>> 513 } // end of while (theParticipant.theNucleus->GetNextNucleon())
>> 514 } // end of while (theParticipants.Next())
>> 515
>> 516 NumberOfInvolvedTargetNucleon=NumberOfInvolvedNucleon;
>> 517
>> 518 // ---------------- Calculation of creation time for each target nucleon -----------
>> 519 theParticipants.StartLoop(); // restart a loop
>> 520 theParticipants.Next();
>> 521 G4VSplitableHadron * primary = theParticipants.GetInteraction().GetProjectile();
>> 522 G4double betta_z=primary->Get4Momentum().pz()/primary->Get4Momentum().e();
>> 523 primary->SetTimeOfCreation(0.);
>> 524
>> 525 G4double ZcoordinateOfPreviousCollision(0.);
>> 526 G4double ZcoordinateOfCurrentInteraction(0.);
>> 527 G4double TimeOfPreviousCollision(0.);
>> 528 G4double TimeOfCurrentCollision(0);
>> 529
>> 530 theParticipants.theNucleus->StartLoop();
>> 531 G4Nucleon * aNucleon;
>> 532 G4bool theFirstInvolvedNucleon(true);
>> 533 while ( (aNucleon = theParticipants.theNucleus->GetNextNucleon()) )
>> 534 {
>> 535 if(aNucleon->AreYouHit())
>> 536 {
>> 537 if(theFirstInvolvedNucleon)
>> 538 {
>> 539 ZcoordinateOfPreviousCollision=aNucleon->GetPosition().z();
>> 540 theFirstInvolvedNucleon=false;
>> 541 }
1522 542
1523 G4bool G4FTFModel::AdjustNucleonsAlgorithm_Sa << 543 ZcoordinateOfCurrentInteraction=aNucleon->GetPosition().z();
1524 << 544 TimeOfCurrentCollision=TimeOfPreviousCollision+
1525 // Second of the three utility methods used << 545 (ZcoordinateOfCurrentInteraction-ZcoordinateOfPreviousCollision)/betta_z;
1526 // This method returns "false" if it fails << 546 // It is assumed that the nucleons are ordered on increasing z-coordinate ------------
1527 << 547 aNucleon->GetSplitableHadron()->SetTimeOfCreation(TimeOfCurrentCollision);
1528 // Ascribing of the involved nucleons Pt an << 548
1529 G4double Dcor = theParameters->GetDofNuclea << 549 ZcoordinateOfPreviousCollision=ZcoordinateOfCurrentInteraction;
1530 G4double DcorP = 0.0, DcorT = 0.0; << 550 TimeOfPreviousCollision=TimeOfCurrentCollision;
1531 if ( ProjectileResidualMassNumber != 0 ) Dc << 551 } // end of if(aNucleon->AreYouHit())
1532 if ( TargetResidualMassNumber != 0 ) Dc << 552 } // end of while (theParticipant.theNucleus->GetNextNucleon())
1533 G4double AveragePt2 = theParameters->GetPt2 << 553 //
1534 G4double maxPtSquare = theParameters->GetMa << 554 // The algorithm can be improved, but it will be more complicated, and will require
1535 << 555 // changes in G4DiffractiveExcitation.cc and G4ElasticHNScattering.cc
1536 G4double ScaleFactor = 1.0; << 556 } // Uzhi 26 July 2009
1537 G4bool OuterSuccess = true; <<
1538 const G4int maxNumberOfLoops = 1000; <<
1539 const G4int maxNumberOfTries = 10000; <<
1540 G4int loopCounter = 0; <<
1541 G4int NumberOfTries = 0; <<
1542 do { // Outmost do while loop <<
1543 OuterSuccess = true; <<
1544 G4bool loopCondition = false; <<
1545 do { // Intermediate do while loop <<
1546 if ( NumberOfTries == 100*(NumberOfTrie <<
1547 // At large number of tries it would <<
1548 ScaleFactor /= 2.0; <<
1549 DcorP *= ScaleFactor; <<
1550 DcorT *= ScaleFactor; <<
1551 AveragePt2 *= ScaleFactor; <<
1552 #ifdef debugAdjust <<
1553 //G4cout << "NumberOfTries ScaleFacto <<
1554 #endif <<
1555 } <<
1556 <<
1557 // Some kinematics <<
1558 if ( interactionCase == 1 ) { <<
1559 } else if ( interactionCase == 2 ) { <<
1560 #ifdef debugAdjust <<
1561 G4cout << "ProjectileResidualMassNumb <<
1562 #endif <<
1563 if ( ProjectileResidualMassNumber > 1 <<
1564 common.PtNucleon = GaussianPt( Aver <<
1565 } else { <<
1566 common.PtNucleon = G4ThreeVector( 0 <<
1567 } <<
1568 common.PtResidual = - common.PtNucleo <<
1569 common.Mprojectile = std::sqrt( sqr <<
1570 + std::sqrt( sqr <<
1571 #ifdef debugAdjust <<
1572 G4cout << "SqrtS < Mtarget + Mproject <<
1573 << " " << common.Mprojectile < <<
1574 #endif <<
1575 common.M2projectile = sqr( common.Mpr <<
1576 if ( common.SqrtS < common.Mtarget + <<
1577 OuterSuccess = false; <<
1578 loopCondition = true; <<
1579 continue; <<
1580 } <<
1581 } else if ( interactionCase == 3 ) { <<
1582 if ( ProjectileResidualMassNumber > 1 <<
1583 common.PtNucleonP = GaussianPt( Ave <<
1584 } else { <<
1585 common.PtNucleonP = G4ThreeVector( <<
1586 } <<
1587 common.PtResidualP = - common.PtNucle <<
1588 if ( TargetResidualMassNumber > 1 ) { <<
1589 common.PtNucleonT = GaussianPt( Ave <<
1590 } else { <<
1591 common.PtNucleonT = G4ThreeVector( <<
1592 } <<
1593 common.PtResidualT = - common.PtNucle <<
1594 common.Mprojectile = std::sqrt( sqr <<
1595 + std::sqrt( sqr <<
1596 common.M2projectile = sqr( common.Mpr <<
1597 common.Mtarget = std::sqrt( sqr( co <<
1598 + std::sqrt( sqr( co <<
1599 common.M2target = sqr( common.Mtarget <<
1600 if ( common.SqrtS < common.Mprojectil <<
1601 OuterSuccess = false; <<
1602 loopCondition = true; <<
1603 continue; <<
1604 } <<
1605 } // End-if on interactionCase <<
1606 557
1607 G4int numberOfTimesExecuteInnerLoop = 1 << 558 // ------------------------------------------------------------
1608 if ( interactionCase == 3 ) numberOfTim << 559 G4bool G4FTFModel::PutOnMassShell()
1609 for ( G4int iExecute = 0; iExecute < nu << 560 {
1610 << 561 //G4cout<<"PutOnMassShell start "<<G4endl;
1611 G4bool InnerSuccess = true; << 562 if(std::abs(theProjectile.GetDefinition()->GetBaryonNumber()) > 1) // !!!!
1612 G4bool isTargetToBeHandled = ( intera << 563 { // The projectile is a nucleus or an anti-nucleus
1613 ( inte << 564 //G4cout<<"PutOnMassShell AA "<<G4endl;
1614 G4bool condition = false; << 565 G4LorentzVector P_total(0.,0.,0.,0.);
1615 if ( isTargetToBeHandled ) { << 566
1616 condition = ( TargetResidualMassNum << 567 G4LorentzVector P_participants(0.,0.,0.,0.);
1617 } else { // Projectile to be handled << 568 G4int MultiplicityOfParticipants(0);
1618 condition = ( ProjectileResidualMas << 569
1619 } << 570 G4V3DNucleus * ProjectileNucleus =theParticipants.GetProjectileNucleus();
1620 if ( condition ) { << 571 ProjectileNucleus->StartLoop();
1621 const G4int maxNumberOfInnerLoops = << 572
1622 G4int innerLoopCounter = 0; << 573 G4Nucleon * ProjectileNucleon;
1623 do { // Inner do while loop << 574 while ( (ProjectileNucleon=ProjectileNucleus->GetNextNucleon()) )
1624 InnerSuccess = true; << 575 {
1625 if ( isTargetToBeHandled ) { << 576 if ( ProjectileNucleon->AreYouHit() )
1626 G4double Xcenter = 0.0; << 577 { // Projectile nucleon was involved in the interaction.
1627 if ( interactionCase == 1 ) { << 578 P_total+=ProjectileNucleon->Get4Momentum();
1628 common.PtNucleon = GaussianPt << 579 MultiplicityOfParticipants++;
1629 common.PtResidual = - common. << 580 P_participants+=ProjectileNucleon->Get4Momentum();
1630 common.Mtarget = std::sqrt( <<
1631 + std::sqrt( <<
1632 if ( common.SqrtS < common.Mp <<
1633 InnerSuccess = false; <<
1634 continue; <<
1635 } <<
1636 Xcenter = std::sqrt( sqr( com <<
1637 / common.Mtarget; <<
1638 } else { <<
1639 Xcenter = std::sqrt( sqr( com <<
1640 / common.Mtarget; <<
1641 } <<
1642 G4ThreeVector tmpX = GaussianPt <<
1643 common.XminusNucleon = Xcenter <<
1644 if ( common.XminusNucleon <= 0. <<
1645 InnerSuccess = false; <<
1646 continue; <<
1647 } <<
1648 common.XminusResidual = 1.0 - c <<
1649 } else { // Projectile to be han <<
1650 G4ThreeVector tmpX = GaussianPt <<
1651 G4double Xcenter = 0.0; <<
1652 if ( interactionCase == 2 ) { <<
1653 Xcenter = std::sqrt( sqr( com <<
1654 / common.Mprojectil <<
1655 } else { <<
1656 Xcenter = std::sqrt( sqr( com <<
1657 / common.Mprojectil <<
1658 } <<
1659 common.XplusNucleon = Xcenter + <<
1660 if ( common.XplusNucleon <= 0.0 <<
1661 InnerSuccess = false; <<
1662 continue; <<
1663 } <<
1664 common.XplusResidual = 1.0 - co <<
1665 } // End-if on isTargetToBeHandl <<
1666 } while ( ( ! InnerSuccess ) && <<
1667 ++innerLoopCounter < ma <<
1668 if ( innerLoopCounter >= maxNumberO <<
1669 #ifdef debugAdjust <<
1670 G4cout << "BAD situation: forced <<
1671 #endif <<
1672 return false; <<
1673 } 581 }
1674 } else { // condition is false << 582 } // End of while ( (ProjectileNucleon=ProjectileNucleus->GetNextNucleon()) )
1675 if ( isTargetToBeHandled ) { << 583 //G4cout<<"MultParts mom Pr "<<MultiplicityOfParticipants<<" "<<P_participants<<G4endl;
1676 common.XminusNucleon = 1.0; << 584 //------------------
1677 common.XminusResidual = 1.0; // << 585 G4int ResidualMassNumber(0);
1678 } else { // Projectile to be handl << 586 G4int ResidualCharge(0);
1679 common.XplusNucleon = 1.0; << 587 ResidualExcitationEnergy=0.;
1680 common.XplusResidual = 1.0; // << 588 G4LorentzVector PnuclearResidual(0.,0.,0.,0.);
1681 } << 589
1682 } // End-if on condition << 590 G4double ExcitationEnergyPerWoundedNucleon=
1683 << 591 theParameters->GetExcitationEnergyPerWoundedNucleon();
1684 } // End of for loop on iExecute << 592
1685 << 593 G4V3DNucleus * TargetNucleus =theParticipants.GetWoundedNucleus();
1686 if ( interactionCase == 1 ) { << 594 TargetNucleus->StartLoop();
1687 common.M2target = ( sqr( common.TN << 595
1688 / common.XminusN << 596 G4Nucleon * TargetNucleon;
1689 + ( sqr( common.TR << 597 while ( (TargetNucleon=TargetNucleus->GetNextNucleon()) )
1690 / common.XminusR << 598 {
1691 loopCondition = ( common.SqrtS < comm << 599 P_total+=TargetNucleon->Get4Momentum();
1692 } else if ( interactionCase == 2 ) { << 600 if ( TargetNucleon->AreYouHit() )
1693 #ifdef debugAdjust << 601 { // Target nucleon was involved in the interaction.
1694 G4cout << "TNucleonMass PtNucleon Xpl << 602 MultiplicityOfParticipants++;
1695 << common.PtNucleon << " " << << 603 P_participants+=TargetNucleon->Get4Momentum();
1696 << "TResidualMass PtResidual X << 604 ResidualExcitationEnergy+=ExcitationEnergyPerWoundedNucleon;
1697 << common.PtResidual << " " < << 605 }
1698 #endif << 606 else
1699 common.M2projectile = ( sqr( commo << 607 {
1700 / common.Xpl << 608 ResidualMassNumber+=1;
1701 + ( sqr( commo << 609 ResidualCharge+=(G4int) TargetNucleon->GetDefinition()->GetPDGCharge();
1702 / common.Xpl << 610 PnuclearResidual+=TargetNucleon->Get4Momentum();
1703 #ifdef debugAdjust << 611 }
1704 G4cout << "SqrtS < Mtarget + std::sqr << 612 } // End of while ( (TargetNucleon=TargetNucleus->GetNextNucleon()) )
1705 << common.Mtarget << " " << st <<
1706 << common.Mtarget + std::sqrt( <<
1707 #endif <<
1708 loopCondition = ( common.SqrtS < comm <<
1709 } else if ( interactionCase == 3 ) { <<
1710 #ifdef debugAdjust <<
1711 G4cout << "PtNucleonP " << common.PtN <<
1712 << "XplusNucleon XplusResidual <<
1713 << " " << common.XplusResidual <<
1714 << "PtNucleonT " << common.PtN <<
1715 << "XminusNucleon XminusResidu <<
1716 << " " << common.XminusResidua <<
1717 #endif <<
1718 common.M2projectile = ( sqr( common <<
1719 / common.Xplu <<
1720 + ( sqr( common <<
1721 / common.Xplu <<
1722 common.M2target = ( sqr( common.TN <<
1723 / common.XminusN <<
1724 + ( sqr( common.TR <<
1725 / common.XminusR <<
1726 loopCondition = ( common.SqrtS < ( <<
1727 + std::sqrt( common.M2target ) ) <<
1728 } // End-if on interactionCase <<
1729 <<
1730 } while ( loopCondition && <<
1731 ++NumberOfTries < maxNumberOfTr <<
1732 if ( NumberOfTries >= maxNumberOfTries ) <<
1733 #ifdef debugAdjust <<
1734 G4cout << "BAD situation: forced exit o <<
1735 #endif <<
1736 return false; <<
1737 } <<
1738 <<
1739 // kinematics <<
1740 G4double Yprojectile = 0.0, YprojectileNu <<
1741 G4double DecayMomentum2 = sqr( common.S ) <<
1742 - 2.0 * ( commo <<
1743 + com <<
1744 if ( interactionCase == 1 ) { <<
1745 common.WminusTarget = ( common.S - co <<
1746 + std::sqrt( De <<
1747 common.WplusProjectile = common.SqrtS - <<
1748 common.Pzprojectile = common.WplusPro <<
1749 - common.M2projec <<
1750 common.Eprojectile = common.WplusPro <<
1751 + common.M2projec <<
1752 Yprojectile = 0.5 * G4Log( ( common. <<
1753 / ( common. <<
1754 #ifdef debugAdjust <<
1755 G4cout << "DecayMomentum2 " << DecayMom <<
1756 << "WminusTarget WplusProjectile <<
1757 << " " << common.WplusProjectile <<
1758 << "Yprojectile " << Yprojectile <<
1759 #endif <<
1760 common.Mt2targetNucleon = sqr( common.T <<
1761 common.PztargetNucleon = - common.Wminu <<
1762 + common.Mt2ta <<
1763 / ( 2.0 * co <<
1764 common.EtargetNucleon = common.Wminus <<
1765 + common.Mt2tar <<
1766 / ( 2.0 * com <<
1767 YtargetNucleon = 0.5 * G4Log( ( commo <<
1768 / ( commo <<
1769 #ifdef debugAdjust <<
1770 G4cout << "YtN Ytr YtN-Ytr " << " " << <<
1771 << " " << YtargetNucleon - commo <<
1772 << "YtN Ypr YtN-Ypr " << " " << <<
1773 << " " << YtargetNucleon - Yproj <<
1774 #endif <<
1775 if ( std::abs( YtargetNucleon - common. <<
1776 Yprojectile < YtargetNucleon ) { <<
1777 OuterSuccess = false; <<
1778 continue; <<
1779 } <<
1780 } else if ( interactionCase == 2 ) { <<
1781 common.WplusProjectile = ( common.S + <<
1782 + std::sqrt( <<
1783 common.WminusTarget = common.SqrtS - co <<
1784 common.Pztarget = - common.WminusTarget <<
1785 common.Etarget = common.WminusTarget <<
1786 Ytarget = 0.5 * G4Log( ( common.Etarg <<
1787 / ( common.Etarg <<
1788 #ifdef debugAdjust <<
1789 G4cout << "DecayMomentum2 " << DecayMom <<
1790 << "WminusTarget WplusProjectile <<
1791 << " " << common.WplusProjectile <<
1792 << "Ytarget " << Ytarget << G4en <<
1793 #endif <<
1794 common.Mt2projectileNucleon = sqr( comm <<
1795 common.PzprojectileNucleon = common.W <<
1796 - common.M <<
1797 / ( 2.0 <<
1798 common.EprojectileNucleon = common.W <<
1799 + common.M <<
1800 / ( 2.0 <<
1801 YprojectileNucleon = 0.5 * G4Log( ( c <<
1802 / ( c <<
1803 #ifdef debugAdjust <<
1804 G4cout << "YpN Ypr YpN-Ypr " << " " << <<
1805 << " " << YprojectileNucleon - c <<
1806 << "YpN Ytr YpN-Ytr " << " " << <<
1807 << " " << YprojectileNucleon - Y <<
1808 #endif <<
1809 if ( std::abs( YprojectileNucleon - com <<
1810 Ytarget > YprojectileNucleon ) { <<
1811 OuterSuccess = false; <<
1812 continue; <<
1813 } <<
1814 } else if ( interactionCase == 3 ) { <<
1815 common.WplusProjectile = ( common.S + <<
1816 + std::sqrt( <<
1817 common.WminusTarget = common.SqrtS - co <<
1818 common.Mt2projectileNucleon = sqr( comm <<
1819 common.PzprojectileNucleon = common.W <<
1820 - common.M <<
1821 / ( 2.0 <<
1822 common.EprojectileNucleon = common.W <<
1823 + common.M <<
1824 / ( 2.0 <<
1825 YprojectileNucleon = 0.5 * G4Log( ( c <<
1826 / ( c <<
1827 common.Mt2targetNucleon = sqr( common.T <<
1828 common.PztargetNucleon = - common.Wminu <<
1829 + common.Mt2ta <<
1830 / ( 2.0 * co <<
1831 common.EtargetNucleon = common.Wminu <<
1832 + common.Mt2ta <<
1833 / ( 2.0 * co <<
1834 YtargetNucleon = 0.5 * G4Log( ( commo <<
1835 / ( commo <<
1836 if ( std::abs( YtargetNucleon - common. <<
1837 std::abs( YprojectileNucleon - com <<
1838 YprojectileNucleon < YtargetNucleo <<
1839 OuterSuccess = false; <<
1840 continue; <<
1841 } <<
1842 } // End-if on interactionCase <<
1843 <<
1844 } while ( ( ! OuterSuccess ) && <<
1845 ++loopCounter < maxNumberOfLoops <<
1846 if ( loopCounter >= maxNumberOfLoops ) { <<
1847 #ifdef debugAdjust <<
1848 G4cout << "BAD situation: forced exit of <<
1849 #endif <<
1850 return false; <<
1851 } <<
1852 <<
1853 return true; <<
1854 } <<
1855 <<
1856 //------------------------------------------- <<
1857 613
1858 void G4FTFModel::AdjustNucleonsAlgorithm_afte << 614 G4double ResidualMass(0.);
1859 << 615 if(ResidualMassNumber == 0)
1860 << 616 {
1861 << 617 ResidualMass=0.;
1862 // Third of the three utility methods used << 618 ResidualExcitationEnergy=0.;
1863 // and transform back. << 619 }
1864 << 620 else
1865 // New projectile << 621 {
1866 if ( interactionCase == 1 ) { << 622 ResidualMass=G4ParticleTable::GetParticleTable()->GetIonTable()->
1867 common.Pprojectile.setPz( common.Pzprojec << 623 GetIonMass(ResidualCharge ,ResidualMassNumber);
1868 common.Pprojectile.setE( common.Eprojecti << 624 if(ResidualMassNumber == 1) {ResidualExcitationEnergy=0.;}
1869 } else if ( interactionCase == 2 ) { << 625 }
1870 common.Pprojectile.setPx( common.PtNucleo << 626 // ResidualMass +=ResidualExcitationEnergy; // Will be given after checks
1871 common.Pprojectile.setPy( common.PtNucleo << 627
1872 common.Pprojectile.setPz( common.Pzprojec << 628 //G4cout<<"MultPars mom+Tr "<<MultiplicityOfParticipants<<" "<<P_participants<<G4endl;
1873 common.Pprojectile.setE( common.Eprojecti << 629 //G4cout<<"Res "<<ResidualMassNumber<<" "<<PnuclearResidual<<G4endl;
1874 } else if ( interactionCase == 3 ) { << 630 //G4cout<<"P_total "<<P_total<<G4endl;
1875 common.Pprojectile.setPx( common.PtNucleo << 631
1876 common.Pprojectile.setPy( common.PtNucleo << 632 //G4cout<<"Rez A Z M E* "<<ResidualMassNumber<<" "<<ResidualCharge<<" "<<ResidualMass<<" "<<ResidualExcitationEnergy<<G4endl;
1877 common.Pprojectile.setPz( common.Pzprojec << 633 //--------------
1878 common.Pprojectile.setE( common.Eprojecti << 634 G4double SqrtS=P_total.mag();
1879 } << 635 G4double S=P_total.mag2();
1880 #ifdef debugAdjust << 636
1881 G4cout << "Proj after in CMS " << common.Pp << 637 // if(theProjectile.GetDefinition()->GetBaryonNumber() > 1)
1882 #endif << 638 { // For nucleus-nucleus interactions
1883 common.Pprojectile.transform( common.toLab << 639 G4double MassOfParticipants=P_participants.mag();
1884 SelectedAntiBaryon->Set4Momentum( common.Pp << 640 G4double MassOfParticipants2=sqr(MassOfParticipants);
1885 #ifdef debugAdjust << 641
1886 G4cout << "Proj after in Lab " << common.Pp << 642 //G4cout<<"ResidualMass + MassOfParticipants "<<ResidualMass + MassOfParticipants<<G4endl;
1887 #endif << 643 if(SqrtS < ResidualMass + MassOfParticipants) {return false;}
1888 << 644
1889 // New target nucleon << 645 if(SqrtS < ResidualMass+ResidualExcitationEnergy + MassOfParticipants)
1890 if ( interactionCase == 1 ) { << 646 {ResidualExcitationEnergy=0.;}
1891 common.Ptarget.setPx( common.PtNucleon.x( << 647
1892 common.Ptarget.setPy( common.PtNucleon.y( << 648 ResidualMass +=ResidualExcitationEnergy;
1893 common.Ptarget.setPz( common.PztargetNucl << 649 //G4cout<<"Rez A Z M E* "<<ResidualMassNumber<<" "<<ResidualCharge<<" "<<ResidualMass<<" "<<ResidualExcitationEnergy<<G4endl;
1894 common.Ptarget.setE( common.EtargetNucleo << 650
1895 } else if ( interactionCase == 2 ) { << 651 G4double ResidualMass2=sqr(ResidualMass);
1896 common.Ptarget.setPz( common.Pztarget ); << 652
1897 common.Ptarget.setE( common.Etarget ); << 653 G4LorentzRotation toCms(-1*P_total.boostVector());
1898 } else if ( interactionCase == 3 ) { << 654
1899 common.Ptarget.setPx( common.PtNucleonT.x << 655 G4LorentzVector Pcms=toCms*P_participants;
1900 common.Ptarget.setPy( common.PtNucleonT.y << 656 //G4cout<<"Ppart in CMS "<<Ptmp<<G4endl;
1901 common.Ptarget.setPz( common.PztargetNucl << 657
1902 common.Ptarget.setE( common.EtargetNucleo << 658 if ( Pcms.pz() <= 0. )
1903 } << 659 { // "String" moving backwards in CMS, abort collision !!
1904 #ifdef debugAdjust << 660 //G4cout << " abort ColliDeleteVSplitableHadronsion!! " << G4endl;
1905 G4cout << "Targ after in CMS " << common.Pt << 661 return false;
1906 #endif << 662 }
1907 common.Ptarget.transform( common.toLab ); <<
1908 SelectedTargetNucleon->Set4Momentum( common <<
1909 #ifdef debugAdjust <<
1910 G4cout << "Targ after in Lab " << common.Pt <<
1911 #endif <<
1912 <<
1913 // New target residual <<
1914 if ( interactionCase == 1 || interactionC <<
1915 TargetResidualMassNumber = common.T <<
1916 TargetResidualCharge = common.T <<
1917 TargetResidualExcitationEnergy = common.T <<
1918 #ifdef debugAdjust <<
1919 G4cout << "TargetResidualMassNumber Targe <<
1920 << TargetResidualMassNumber << " " <<
1921 << TargetResidualExcitationEnergy <<
1922 #endif <<
1923 if ( TargetResidualMassNumber != 0 ) { <<
1924 G4double Mt2 = 0.0; <<
1925 if ( interactionCase == 1 ) { <<
1926 Mt2 = sqr( common.TResidualMass ) + c <<
1927 TargetResidual4Momentum.setPx( common <<
1928 TargetResidual4Momentum.setPy( common <<
1929 } else { // interactionCase == 3 <<
1930 Mt2 = sqr( common.TResidualMass ) + c <<
1931 TargetResidual4Momentum.setPx( common <<
1932 TargetResidual4Momentum.setPy( common <<
1933 } <<
1934 G4double Pz = - common.WminusTarget * c <<
1935 + Mt2 / ( 2.0 * common.Wm <<
1936 G4double E = common.WminusTarget * c <<
1937 + Mt2 / ( 2.0 * common.Wm <<
1938 TargetResidual4Momentum.setPz( Pz ); <<
1939 TargetResidual4Momentum.setE( E ) ; <<
1940 TargetResidual4Momentum.transform( comm <<
1941 } else { <<
1942 TargetResidual4Momentum = G4LorentzVect <<
1943 } <<
1944 #ifdef debugAdjust <<
1945 G4cout << "Tr N R " << common.Ptarget << <<
1946 #endif <<
1947 } <<
1948 <<
1949 // New projectile residual <<
1950 if ( interactionCase == 2 || interactionC <<
1951 if ( interactionCase == 2 ) { <<
1952 ProjectileResidualMassNumber = co <<
1953 ProjectileResidualCharge = co <<
1954 ProjectileResidualExcitationEnergy = co <<
1955 ProjectileResidualLambdaNumber = co <<
1956 } else { // interactionCase == 3 <<
1957 ProjectileResidualMassNumber = co <<
1958 ProjectileResidualCharge = co <<
1959 ProjectileResidualExcitationEnergy = co <<
1960 ProjectileResidualLambdaNumber = co <<
1961 } <<
1962 #ifdef debugAdjust <<
1963 G4cout << "ProjectileResidualMassNumber P <<
1964 << ProjectileResidualMassNumber << <<
1965 << ProjectileResidualLambdaNumber <<
1966 << ProjectileResidualExcitationEne <<
1967 #endif <<
1968 if ( ProjectileResidualMassNumber != 0 ) <<
1969 G4double Mt2 = 0.0; <<
1970 if ( interactionCase == 2 ) { <<
1971 Mt2 = sqr( common.TResidualMass ) + c <<
1972 ProjectileResidual4Momentum.setPx( co <<
1973 ProjectileResidual4Momentum.setPy( co <<
1974 } else { // interactionCase == 3 <<
1975 Mt2 = sqr( common.PResidualMass ) + c <<
1976 ProjectileResidual4Momentum.setPx( co <<
1977 ProjectileResidual4Momentum.setPy( co <<
1978 } <<
1979 G4double Pz = common.WplusProjectile <<
1980 - Mt2 / ( 2.0 * common.Wp <<
1981 G4double E = common.WplusProjectile <<
1982 + Mt2 / ( 2.0 * common.Wp <<
1983 ProjectileResidual4Momentum.setPz( Pz ) <<
1984 ProjectileResidual4Momentum.setE( E ); <<
1985 ProjectileResidual4Momentum.transform( <<
1986 } else { <<
1987 ProjectileResidual4Momentum = G4Lorentz <<
1988 } <<
1989 #ifdef debugAdjust <<
1990 G4cout << "Pr N R " << common.Pprojectile <<
1991 << " " << ProjectileResidual <<
1992 #endif <<
1993 } <<
1994 663
1995 } << 664 toCms.rotateZ(-1*Pcms.phi()); // Uzhi 5.12.09
>> 665 toCms.rotateY(-1*Pcms.theta()); // Uzhi 5.12.09
1996 666
>> 667 //G4cout<<"Mpa M res "<<ResidualMass<<" "<<MassOfParticipants<<" "<<SqrtS<<G4endl;
>> 668 G4LorentzRotation toLab(toCms.inverse());
1997 669
1998 //=========================================== << 670 G4double DecayMomentum2=
>> 671 sqr(S)+sqr(MassOfParticipants2)+ sqr(ResidualMass2) -
>> 672 2.*S*MassOfParticipants2 - 2.*S*ResidualMass2
>> 673 -2.*MassOfParticipants2*ResidualMass2;
>> 674
>> 675 if(DecayMomentum2 < 0.) return false;
>> 676
>> 677 DecayMomentum2/=(4.*S);
>> 678 G4double DecayMomentum = std::sqrt(DecayMomentum2);
>> 679 //G4cout<<"DecayMomentum "<<DecayMomentum<<G4endl;
>> 680
>> 681 G4LorentzVector New_P_participants(0.,0.,DecayMomentum,
>> 682 std::sqrt(DecayMomentum2+MassOfParticipants2));
>> 683 G4LorentzVector New_PnuclearResidual(0.,0.,-DecayMomentum,
>> 684 std::sqrt(DecayMomentum2+ResidualMass2));
>> 685
>> 686 //G4cout<<"MultPars mom "<<MultiplicityOfParticipants<<" "<<New_P_participants<<G4endl;
>> 687 //G4cout<<"Res "<<ResidualMassNumber<<" "<<New_PnuclearResidual<<G4endl;
>> 688 New_P_participants.transform(toLab);
>> 689 New_PnuclearResidual.transform(toLab);
>> 690
>> 691 //G4cout<<"MultPars mom "<<MultiplicityOfParticipants<<" "<<New_P_participants<<G4endl;
>> 692 //G4cout<<"Res "<<ResidualMassNumber<<" "<<New_PnuclearResidual<<G4endl;
>> 693
>> 694 G4LorentzVector DeltaP_participants=(New_P_participants - P_participants)/
>> 695 ((G4double) MultiplicityOfParticipants);
>> 696
>> 697 //G4cout<<"DeltaP_participants "<<DeltaP_participants<<G4endl;
>> 698 //-------------
>> 699 ProjectileNucleus->StartLoop();
>> 700 while ( (ProjectileNucleon=ProjectileNucleus->GetNextNucleon()) )
>> 701 {
>> 702 if ( ProjectileNucleon->AreYouHit() )
>> 703 { // Projectile nucleon was involved in the interaction.
>> 704 G4LorentzVector Ptmp=ProjectileNucleon->Get4Momentum() + DeltaP_participants;
>> 705 ProjectileNucleon->GetSplitableHadron()->Set4Momentum(Ptmp);
>> 706 ProjectileNucleon->SetMomentum(Ptmp);
>> 707 }
>> 708 } // End of while ( (ProjectileNucleon=ProjectileNucleus->GetNextNucleon()) )
>> 709
>> 710 //-------------
>> 711 TargetNucleus->StartLoop();
>> 712 while ( (TargetNucleon=TargetNucleus->GetNextNucleon()) )
>> 713 {
>> 714 if ( TargetNucleon->AreYouHit() )
>> 715 { // Target nucleon was involved in the interaction.
>> 716 G4LorentzVector Ptmp=TargetNucleon->Get4Momentum() + DeltaP_participants;
>> 717 TargetNucleon->GetSplitableHadron()->Set4Momentum(Ptmp);
>> 718 }
>> 719 } // End of while ( (TargetNucleon=TargetNucleus->GetNextNucleon()) )
>> 720
>> 721 Residual4Momentum = New_PnuclearResidual;
>> 722 // return true;
>> 723 } // End of if(theProjectile.GetDefinition()->GetBaryonNumber() > 1)
>> 724
>> 725 return true;
>> 726 }
>> 727 //---------------------------------------------------------------------
>> 728 // -------- The projectile is hadron, or baryon, or anti-baryon -------
>> 729 // -------------- Properties of the projectile ------------------------
>> 730 theParticipants.StartLoop(); // restart a loop
>> 731 theParticipants.Next();
>> 732 G4VSplitableHadron * primary = theParticipants.GetInteraction().GetProjectile();
>> 733 G4LorentzVector Pprojectile=primary->Get4Momentum();
>> 734
>> 735 // 13.06.2012 G4bool ProjectileIsAntiBaryon = primary->GetDefinition()->GetBaryonNumber() < 0;
>> 736
>> 737 //G4cout<<"PutOnMass Pprojectile "<<Pprojectile<<G4endl;
>> 738 // To get original projectile particle
>> 739
>> 740 if(Pprojectile.z() < 0.){return false;}
>> 741
>> 742 G4double Mprojectile = Pprojectile.mag();
>> 743 G4double M2projectile = Pprojectile.mag2();
>> 744 //-------------------------------------------------------------
>> 745 G4LorentzVector Psum = Pprojectile;
>> 746
>> 747 G4double SumMasses = Mprojectile + 20.*MeV; // 13.12.09
>> 748 // Separation energy for projectile
>> 749 // if(ProjectileIsAntiBaryon) {SumMasses = Mprojectile;}
>> 750 //G4cout<<"SumMasses Pr "<<SumMasses<<G4endl;
>> 751 //--------------- Target nucleus ------------------------------
>> 752 G4V3DNucleus *theNucleus = GetWoundedNucleus();
>> 753 G4int ResidualMassNumber=theNucleus->GetMassNumber();
>> 754 G4int ResidualCharge =theNucleus->GetCharge();
>> 755
>> 756 ResidualExcitationEnergy=0.;
>> 757 G4LorentzVector Ptarget(0.,0.,0.,0.);
>> 758 G4LorentzVector PnuclearResidual(0.,0.,0.,0.); // Uzhi 12.06.2012
>> 759
>> 760 G4double ExcitationEnergyPerWoundedNucleon=
>> 761 theParameters->GetExcitationEnergyPerWoundedNucleon();
>> 762
>> 763 //G4cout<<"ExcitationEnergyPerWoundedNucleon "<<ExcitationEnergyPerWoundedNucleon<<G4endl;
>> 764
>> 765 theNucleus->StartLoop();
>> 766
>> 767 while (G4Nucleon * aNucleon = theNucleus->GetNextNucleon())
>> 768 {
>> 769 Ptarget+=aNucleon->Get4Momentum();
>> 770
>> 771 if(aNucleon->AreYouHit())
>> 772 { // Involved nucleons
>> 773 //G4cout<<"PutOn Tr "<<aNucleon->Get4Momentum()<<G4endl;
>> 774 // Psum += aNucleon->Get4Momentum(); // Uzhi 20 Sept.
>> 775 // if(!ProjectileIsAntiBaryon) // Uzhi 13.06.2012
>> 776 // {
>> 777 SumMasses += std::sqrt(sqr(aNucleon->GetDefinition()->GetPDGMass()) //Uzhi 12.06.2012
>> 778 + aNucleon->Get4Momentum().perp2()); //Uzhi 12.06.2012
>> 779 SumMasses += 20.*MeV; // 13.12.09 Separation energy for a nucleon
>> 780 ResidualExcitationEnergy+=ExcitationEnergyPerWoundedNucleon;
>> 781 /* //Uzhi 13.06.2012
>> 782 } else
>> 783 {
>> 784 SumMasses += aNucleon->Get4Momentum().mag(); // 4.12.2010
>> 785 G4LorentzVector tmp=aNucleon->Get4Momentum();
>> 786 tmp.setE(aNucleon->Get4Momentum().mag()); // It is need to save mass 6.12.2011
>> 787 aNucleon->SetMomentum(tmp);
>> 788 }
>> 789 */ //Uzhi 13.06.2012
1999 790
2000 void G4FTFModel::BuildStrings( G4ExcitedStrin << 791 ResidualMassNumber--;
2001 // Loop over all collisions; find all prima << 792 ResidualCharge-=(G4int) aNucleon->GetDefinition()->GetPDGCharge();
2002 // (targets may be duplicate in the List (t << 793 }
2003 << 794 else
2004 G4ExcitedString* FirstString( 0 ); // I << 795 { // Spectator nucleons
2005 G4ExcitedString* SecondString( 0 ); // t << 796 PnuclearResidual += aNucleon->Get4Momentum(); // Uzhi 12.06.2012
2006 << 797 } // end of if(!aNucleon->AreYouHit())
2007 if ( ! GetProjectileNucleus() ) { << 798 } // end of while (theNucleus->GetNextNucleon())
2008 << 799
2009 std::vector< G4VSplitableHadron* > primar << 800 Psum += Ptarget;
2010 theParticipants.StartLoop(); << 801 PnuclearResidual.setPz(0.); PnuclearResidual.setE(0.); // Uzhi 12.06.2012
2011 while ( theParticipants.Next() ) { /* Lo << 802 //G4cout<<"ResidualCharge ,ResidualMassNumber "<<ResidualCharge<<" "<<ResidualMassNumber<<" "<<Ptarget<<" "<<PnuclearResidual<<G4endl;
2012 const G4InteractionContent& interaction << 803
2013 // do not allow for duplicates ... << 804 G4double ResidualMass(0.);
2014 if ( interaction.GetStatus() ) { << 805 if(ResidualMassNumber == 0)
2015 if ( primaries.end() == std::find( pr << 806 {
2016 in << 807 ResidualMass=0.;
2017 primaries.push_back( interaction.Ge << 808 ResidualExcitationEnergy=0.;
>> 809 }
>> 810 else
>> 811 {
>> 812 ResidualMass=G4ParticleTable::GetParticleTable()->GetIonTable()->
>> 813 GetIonMass(ResidualCharge ,ResidualMassNumber);
>> 814 if(ResidualMassNumber == 1) {ResidualExcitationEnergy=0.;}
2018 } 815 }
2019 } << 816
2020 } << 817 // ResidualMass +=ResidualExcitationEnergy; // Will be given after checks
2021 << 818 //G4cout<<"SumMasses And ResidualMass "<<SumMasses<<" "<<ResidualMass<<G4endl;
2022 #ifdef debugBuildString << 819 SumMasses += std::sqrt(sqr(ResidualMass)+PnuclearResidual.mag2()); // Uzhi 12.06.2012
2023 G4cout << "G4FTFModel::BuildStrings()" << << 820 //G4cout<<"SumMasses + ResM "<<SumMasses<<G4endl;
2024 << "Number of projectile strings " << 821 //G4cout<<"Psum "<<Psum<<G4endl;
2025 #endif << 822 //-------------------------------------------------------------
2026 << 823
2027 for ( unsigned int ahadron = 0; ahadron < << 824 G4double SqrtS=Psum.mag();
2028 G4bool isProjectile( true ); << 825 G4double S=Psum.mag2();
2029 //G4cout << "primaries[ ahadron ] " << << 826
2030 //if ( primaries[ ahadron ]->GetStatus( << 827 //G4cout<<"SqrtS < SumMasses "<<SqrtS<<" "<<SumMasses<<G4endl;
2031 FirstString = 0; SecondString = 0; << 828
2032 if ( primaries[ahadron]->GetStatus() == << 829 if(SqrtS < SumMasses) {return false;} // It is impossible to simulate
2033 theExcitation->CreateStrings( primarie << 830 // after putting nuclear nucleons
2034 FirstStr << 831 // on mass-shell
2035 NumberOfProjectileSpectatorNucleons--; << 832
2036 } else if ( primaries[ahadron]->GetStat << 833 SumMasses -= std::sqrt(sqr(ResidualMass)+PnuclearResidual.mag2()); // Uzhi 12.06.2012
2037 && primaries[ahadron]->GetSoft << 834 SumMasses += std::sqrt(sqr(ResidualMass+ResidualExcitationEnergy)
2038 theExcitation->CreateStrings( primarie << 835 +PnuclearResidual.mag2()); // Uzhi 12.06.2012
2039 FirstStr << 836 if(SqrtS < SumMasses) // Uzhi 12.06.2012
2040 NumberOfProjectileSpectatorNucleons--; << 837 {
2041 } else if ( primaries[ahadron]->GetStat << 838 SumMasses -= std::sqrt(sqr(ResidualMass+ResidualExcitationEnergy)
2042 && primaries[ahadron]->GetSoft << 839 +PnuclearResidual.mag2()); // Uzhi 12.06.2012
2043 G4LorentzVector ParticleMomentum=prima << 840 SumMasses += std::sqrt(sqr(ResidualMass)+PnuclearResidual.mag2());// Uzhi 12.06.2012
2044 G4KineticTrack* aTrack = new G4Kinetic << 841 ResidualExcitationEnergy=0.;
2045 << 842 }
2046 << 843
2047 << 844 ResidualMass +=ResidualExcitationEnergy;
2048 FirstString = new G4ExcitedString( aTr << 845 // SumMasses +=ResidualExcitationEnergy; // Uzhi 12.06.2012
2049 } else if (primaries[ahadron]->GetStatu << 846 //G4cout<<"ResidualMass SumMasses ResidualExcitationEnergy "<<ResidualMass<<" "<<SumMasses<<" "<<ResidualExcitationEnergy<<G4endl;
2050 G4LorentzVector ParticleMomentum=prima << 847 //-------------------------------------------------------------
2051 G4KineticTrack* aTrack = new G4Kinetic << 848 // Sampling of nucleons what are transfered to delta-isobars --
2052 << 849 G4int MaxNumberOfDeltas = (int)((SqrtS - SumMasses)/(400.*MeV));
2053 << 850 G4int NumberOfDeltas(0);
2054 << 851 //G4cout<<"MaxNumberOfDeltas "<<MaxNumberOfDeltas<<G4endl;
2055 FirstString = new G4ExcitedString( aTr << 852 if(theNucleus->GetMassNumber() != 1)
2056 NumberOfProjectileSpectatorNucleons--; << 853 {
2057 } else { << 854 //G4cout<<"NumberOfInvolvedNucleon "<<NumberOfInvolvedNucleon<<G4endl;
2058 G4cout << "Something wrong in FTF Mod << 855 G4double ProbDeltaIsobar(0.05); // Uzhi 6.07.2012
2059 } << 856 for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
2060 << 857 {
2061 if ( FirstString != 0 ) strings->push_ << 858 if((G4UniformRand() < ProbDeltaIsobar)&&(NumberOfDeltas < MaxNumberOfDeltas))
2062 if ( SecondString != 0 ) strings->push_ << 859 {
2063 << 860 NumberOfDeltas++;
2064 #ifdef debugBuildString << 861 G4VSplitableHadron * targetSplitable=TheInvolvedNucleon[i]->GetSplitableHadron();
2065 G4cout << "FirstString & SecondString? << 862 G4double MassDec=std::sqrt(sqr(targetSplitable->GetDefinition()->GetPDGMass())
2066 if ( FirstString->IsExcited() ) { << 863 + targetSplitable->Get4Momentum().perp2());
2067 G4cout << "Quarks on the FirstString << 864
2068 << " " << FirstString->GetLeft << 865 G4int PDGcode = targetSplitable->GetDefinition()->GetPDGEncoding();
2069 } else { << 866 G4ParticleDefinition* Old_def = targetSplitable->GetDefinition();
2070 G4cout << "Kinetic track is stored" < << 867
2071 } << 868 G4int newPDGcode = PDGcode/10; newPDGcode=newPDGcode*10+4; // Delta
2072 #endif << 869 G4ParticleDefinition* ptr =
2073 << 870 G4ParticleTable::GetParticleTable()->FindParticle(newPDGcode);
2074 } << 871 targetSplitable->SetDefinition(ptr);
2075 << 872 G4double MassInc=std::sqrt(sqr(targetSplitable->GetDefinition()->GetPDGMass())
2076 #ifdef debugBuildString << 873 + targetSplitable->Get4Momentum().perp2());
2077 if ( FirstString->IsExcited() ) { << 874 if(SqrtS < SumMasses + MassInc - MassDec) // Uzhi 12.06.2012
2078 G4cout << "Check 1 string " << strings- << 875 { // Change cannot be acsepted!
2079 << " " << strings->operator[](0) << 876 targetSplitable->SetDefinition(Old_def);
2080 } << 877 ProbDeltaIsobar = 0.;
2081 #endif << 878 } else
2082 << 879 { // Change is acsepted.
2083 std::for_each( primaries.begin(), primari << 880 SumMasses += (MassInc - MassDec);
2084 primaries.clear(); << 881 }
2085 << 882 }
2086 } else { // Projectile is a nucleus << 883 } // end of for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
2087 << 884 } // end of if(theNucleus.GetMassNumber() != 1)
2088 #ifdef debugBuildString << 885
2089 G4cout << "Building of projectile-like st << 886 //-------------------------------------------------------------
2090 #endif << 887
2091 << 888 G4LorentzRotation toCms(-1*Psum.boostVector());
2092 G4bool isProjectile = true; << 889 G4LorentzVector Ptmp=toCms*Pprojectile;
2093 for ( G4int ahadron = 0; ahadron < Number << 890 if ( Ptmp.pz() <= 0. )
2094 << 891 { // "String" moving backwards in CMS, abort collision !!
2095 #ifdef debugBuildString << 892 //G4cout << " abort ColliDeleteVSplitableHadronsion!! " << G4endl;
2096 G4cout << "Nucleon #, status, intCount << 893 return false;
2097 << TheInvolvedNucleonsOfProjecti << 894 }
2098 << " " << TheInvolvedNucleonsOfP << 895
2099 ->GetSoftCollisionCoun << 896 G4LorentzRotation toLab(toCms.inverse());
2100 #endif << 897
2101 << 898 Ptmp=toCms*Ptarget;
2102 G4VSplitableHadron* aProjectile = << 899 G4double YtargetNucleus=Ptmp.rapidity();
2103 TheInvolvedNucleonsOfProjectile[ ah << 900 //G4cout<<"YtargetNucleus "<<YtargetNucleus<<G4endl;
2104 << 901 //-------------------------------------------------------------
2105 #ifdef debugBuildString << 902 //------- Ascribing of the involved nucleons Pt and Xminus ----
2106 G4cout << G4endl << "ahadron aProjectil << 903 G4double Dcor = theParameters->GetDofNuclearDestruction()/
2107 << " " << aProjectile->GetStatus << 904 theNucleus->GetMassNumber();
2108 #endif << 905
2109 << 906 G4double AveragePt2 = theParameters->GetPt2ofNuclearDestruction();
2110 FirstString = 0; SecondString = 0; << 907 G4double maxPtSquare = theParameters->GetMaxPt2ofNuclearDestruction();
2111 if ( aProjectile->GetStatus() == 0 ) { << 908 //G4cout<<"Dcor "<<theParameters->GetDofNuclearDestruction()<<" "<<Dcor<<" AveragePt2 "<<AveragePt2<<G4endl;
2112 << 909 G4double M2target(0.);
2113 #ifdef debugBuildString << 910 G4double WminusTarget(0.);
2114 G4cout << "Case1 aProjectile->GetStat << 911 G4double WplusProjectile(0.);
2115 #endif << 912
2116 << 913 G4int NumberOfTries(0);
2117 theExcitation->CreateStrings( << 914 G4double ScaleFactor(1.);
2118 TheInvolvedNucleon << 915 G4bool OuterSuccess(true);
2119 isProjectile, Firs << 916 do // while (!OuterSuccess)
2120 NumberOfProjectileSpectatorNucleons-- << 917 {
2121 } else if ( aProjectile->GetStatus() == << 918 OuterSuccess=true;
2122 // Nucleon took part in diffractive i << 919
2123 << 920 do // while (SqrtS < Mprojectile + std::sqrt(M2target))
2124 #ifdef debugBuildString << 921 { // while (DecayMomentum < 0.)
2125 G4cout << "Case2 aProjectile->GetStat << 922
2126 #endif << 923 NumberOfTries++;
2127 << 924
2128 theExcitation->CreateStrings( << 925 if(NumberOfTries == 100*(NumberOfTries/100)) // 100
2129 TheInvolvedNucleon << 926 { // At large number of tries it would be better to reduce the values
2130 isProjectile, Firs << 927 ScaleFactor/=2.;
2131 NumberOfProjectileSpectatorNucleons-- << 928 Dcor *=ScaleFactor;
2132 } else if ( aProjectile->GetStatus() == << 929 AveragePt2 *=ScaleFactor;
2133 HighEnergyInter ) { << 930 }
2134 // Nucleon was considered as a parici <<
2135 // but the interaction was skipped du <<
2136 // It is now considered as an involve <<
2137 <<
2138 #ifdef debugBuildString <<
2139 G4cout << "Case3 aProjectile->GetStat <<
2140 #endif <<
2141 <<
2142 G4LorentzVector ParticleMomentum = aP <<
2143 G4KineticTrack* aTrack = new G4Kineti <<
2144 <<
2145 <<
2146 <<
2147 FirstString = new G4ExcitedString( aT <<
2148 <<
2149 #ifdef debugBuildString <<
2150 G4cout << " Strings are built for nuc <<
2151 << " the interaction was skipp <<
2152 #endif <<
2153 <<
2154 } else if ( aProjectile->GetStatus() == <<
2155 // Nucleon which was involved in the <<
2156 <<
2157 #ifdef debugBuildString <<
2158 G4cout << "Case4 aProjectile->GetStat <<
2159 #endif <<
2160 <<
2161 G4LorentzVector ParticleMomentum = aP <<
2162 G4KineticTrack* aTrack = new G4Kineti <<
2163 <<
2164 <<
2165 <<
2166 FirstString = new G4ExcitedString( aT <<
2167 <<
2168 #ifdef debugBuildString <<
2169 G4cout << " Strings are build for inv <<
2170 #endif <<
2171 <<
2172 if ( aProjectile->GetStatus() == 2 ) <<
2173 } else { <<
2174 <<
2175 #ifdef debugBuildString <<
2176 G4cout << "Case5 " << G4endl; <<
2177 #endif <<
2178 <<
2179 //TheInvolvedNucleonsOfProjectile[ ah <<
2180 //G4cout << TheInvolvedNucleonsOfProj <<
2181 <<
2182 #ifdef debugBuildString <<
2183 G4cout << " No string" << G4endl; <<
2184 #endif <<
2185 931
2186 } << 932 G4ThreeVector PtSum(0.,0.,0.);
>> 933 G4double XminusSum(0.);
>> 934 G4double Xminus(0.);
>> 935 G4bool InerSuccess=true;
>> 936
>> 937 do // while(!InerSuccess);
>> 938 {
>> 939 InerSuccess=true;
>> 940
>> 941 PtSum =G4ThreeVector(0.,0.,0.);
>> 942 XminusSum=0.;
>> 943
>> 944 for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
>> 945 {
>> 946 G4Nucleon * aNucleon = TheInvolvedNucleon[i];
>> 947
>> 948 G4ThreeVector tmpPt = GaussianPt(AveragePt2, maxPtSquare);
>> 949 PtSum += tmpPt;
>> 950
>> 951 G4ThreeVector tmpX=GaussianPt(Dcor*Dcor, 1.);
>> 952 Xminus=tmpX.x();
>> 953 XminusSum+=Xminus;
>> 954
>> 955 G4LorentzVector tmp(tmpPt.x(),tmpPt.y(),Xminus, // 6 Dec.2010
>> 956 aNucleon->Get4Momentum().e()); // 6 Dec.2010
>> 957
>> 958 aNucleon->SetMomentum(tmp);
>> 959 } // end of for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
>> 960
>> 961 //---------------------------------------------------------------------------
>> 962 G4double DeltaX = (PtSum.x()-PnuclearResidual.x())/NumberOfInvolvedNucleon;
>> 963 G4double DeltaY = (PtSum.y()-PnuclearResidual.y())/NumberOfInvolvedNucleon;
>> 964 G4double DeltaXminus(0.);
>> 965
>> 966 if(ResidualMassNumber == 0)
>> 967 {
>> 968 DeltaXminus = (XminusSum-1.)/NumberOfInvolvedNucleon;
>> 969 }
>> 970 else
>> 971 {
>> 972 DeltaXminus = -1./theNucleus->GetMassNumber();
>> 973 }
>> 974
>> 975 XminusSum=1.;
>> 976 M2target =0.;
>> 977
>> 978 for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
>> 979 {
>> 980 G4Nucleon * aNucleon = TheInvolvedNucleon[i];
>> 981
>> 982 Xminus = aNucleon->Get4Momentum().pz() - DeltaXminus;
>> 983 XminusSum-=Xminus;
>> 984
>> 985 if(ResidualMassNumber == 0)
>> 986 {
>> 987 if((Xminus <= 0.) || (Xminus > 1.)) {InerSuccess=false; break;}
>> 988 } else
>> 989 {
>> 990 if((Xminus <= 0.) || (Xminus > 1.) ||
>> 991 (XminusSum <=0.) || (XminusSum > 1.)) {InerSuccess=false; break;}
>> 992 }
>> 993
>> 994 G4double Px=aNucleon->Get4Momentum().px() - DeltaX;
>> 995 G4double Py=aNucleon->Get4Momentum().py() - DeltaY;
>> 996
>> 997 // if(!ProjectileIsAntiBaryon) // 4.12.2010
>> 998 // {
>> 999 M2target +=(aNucleon->GetSplitableHadron()->GetDefinition()->GetPDGMass()*
>> 1000 aNucleon->GetSplitableHadron()->GetDefinition()->GetPDGMass() +
>> 1001 Px*Px + Py*Py)/Xminus;
>> 1002 /*
>> 1003 } else
>> 1004 {
>> 1005 M2target +=(aNucleon->Get4Momentum().e() *
>> 1006 aNucleon->Get4Momentum().e() + // 6.12.2010
>> 1007 Px*Px + Py*Py)/Xminus;
>> 1008 }
>> 1009 */
>> 1010 G4LorentzVector tmp(Px,Py,Xminus,aNucleon->Get4Momentum().e()); // 6.12.2010
>> 1011 aNucleon->SetMomentum(tmp);
>> 1012 } // end of for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
>> 1013
>> 1014 if(InerSuccess && (ResidualMassNumber != 0))
>> 1015 {
>> 1016 M2target +=(sqr(ResidualMass) + PnuclearResidual.mag2())/XminusSum;
>> 1017 }
>> 1018 //G4cout<<"InerSuccess "<<InerSuccess<<" "<<SqrtS<<" "<<Mprojectile + std::sqrt(M2target)<<" "<<std::sqrt(M2target)<<G4endl;
>> 1019 //G4int Uzhi;G4cin>>Uzhi;
>> 1020 } while(!InerSuccess);
>> 1021 } while (SqrtS < Mprojectile + std::sqrt(M2target));
>> 1022 //-------------------------------------------------------------
>> 1023 G4double DecayMomentum2= S*S+M2projectile*M2projectile+M2target*M2target
>> 1024 -2.*S*M2projectile - 2.*S*M2target
>> 1025 -2.*M2projectile*M2target;
>> 1026
>> 1027 WminusTarget=(S-M2projectile+M2target+std::sqrt(DecayMomentum2))/2./SqrtS;
>> 1028 WplusProjectile=SqrtS - M2target/WminusTarget;
>> 1029
>> 1030 G4double Pzprojectile=WplusProjectile/2. - M2projectile/2./WplusProjectile;// 8.06.11
>> 1031 G4double Eprojectile =WplusProjectile/2. + M2projectile/2./WplusProjectile;// 8.06.11
>> 1032 G4double Yprojectile=0.5*std::log((Eprojectile+Pzprojectile)/
>> 1033 (Eprojectile-Pzprojectile)); // 1.07.11
>> 1034
>> 1035 //G4cout<<"Yprojectile "<<Yprojectile<<G4endl;
>> 1036 //G4LorentzVector TestPprojectile=Pprojectile;
>> 1037 //TestPprojectile.setPz(Pzprojectile); TestPprojectile.setE(Eprojectile);
>> 1038
>> 1039 //G4cout<<"DecayMomentum2 "<<DecayMomentum2<<G4endl;
>> 1040 //G4cout<<"WminusTarget WplusProjectile "<<WminusTarget<<" "<<WplusProjectile<<G4endl;
>> 1041 //G4int Uzhi;G4cin>>Uzhi;
>> 1042 //-------------------------------------------------------------
>> 1043 for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
>> 1044 {
>> 1045 G4Nucleon * aNucleon = TheInvolvedNucleon[i];
>> 1046 G4LorentzVector tmp=aNucleon->Get4Momentum();
>> 1047
>> 1048 G4double Mt2(0.);
>> 1049
>> 1050 // if(!ProjectileIsAntiBaryon) // 4.12.2010
>> 1051 // {
>> 1052 Mt2 = sqr(tmp.x())+sqr(tmp.y())+
>> 1053 sqr(aNucleon->GetSplitableHadron()->GetDefinition()->GetPDGMass());
>> 1054 /*
>> 1055 } else
>> 1056 {
>> 1057 Mt2 = sqr(tmp.x())+sqr(tmp.y())+ // 4.12.2010
>> 1058 sqr(aNucleon->Get4Momentum().e()); // sqr
>> 1059 }
>> 1060 */
>> 1061 G4double Xminus=tmp.z();
>> 1062
>> 1063 G4double Pz=-WminusTarget*Xminus/2. + Mt2/(2.*WminusTarget*Xminus);
>> 1064 G4double E = WminusTarget*Xminus/2. + Mt2/(2.*WminusTarget*Xminus);
>> 1065 G4double YtargetNucleon=0.5*std::log((E+Pz)/(E-Pz)); // 1.07.11 //Uzhi 20 Sept.
>> 1066
>> 1067 //G4cout<<"YtargetNucleon "<<YtargetNucleon<<G4endl;
>> 1068 //G4cout<<"YtargetNucleon-YtargetNucleus "<<YtargetNucleon-YtargetNucleus<<G4endl;
>> 1069 //G4cout<<"Yprojectile YtargetNucleon "<<Yprojectile<<" "<<YtargetNucleon<<G4endl;
>> 1070 if((std::abs(YtargetNucleon-YtargetNucleus) > 2) ||
>> 1071 (Yprojectile < YtargetNucleon)) {OuterSuccess=false; break;} // 1.07.11
>> 1072
>> 1073 } // end of for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
>> 1074 //if(ProjectileIsAntiBaryon) {G4int Uzhi;G4cin>>Uzhi;}
>> 1075 } while(!OuterSuccess);
>> 1076
>> 1077 //-------------------------------------------------------------
>> 1078 G4double Pzprojectile=WplusProjectile/2. - M2projectile/2./WplusProjectile;
>> 1079 G4double Eprojectile =WplusProjectile/2. + M2projectile/2./WplusProjectile;
>> 1080 Pprojectile.setPz(Pzprojectile); Pprojectile.setE(Eprojectile);
>> 1081 //G4cout<<"Proj after in CMS "<<Pprojectile<<G4endl;
>> 1082
>> 1083 Pprojectile.transform(toLab); // The work with the projectile
>> 1084 primary->Set4Momentum(Pprojectile); // is finished at the moment.
>> 1085 //G4cout<<"Final proj mom "<<primary->Get4Momentum()<<G4endl;
>> 1086
>> 1087 //-------------------------------------------------------------
>> 1088 G4ThreeVector Residual3Momentum(0.,0.,1.);
>> 1089
>> 1090 for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
>> 1091 {
>> 1092 G4Nucleon * aNucleon = TheInvolvedNucleon[i];
>> 1093 G4LorentzVector tmp=aNucleon->Get4Momentum();
>> 1094 //G4cout<<"trg "<<aNucleon->Get4Momentum()<<G4endl;
>> 1095 Residual3Momentum-=tmp.vect();
>> 1096
>> 1097 G4double Mt2(0.);
>> 1098
>> 1099 // if(!ProjectileIsAntiBaryon) // 4.12.2010
>> 1100 // {
>> 1101 Mt2 = sqr(tmp.x())+sqr(tmp.y())+
>> 1102 sqr(aNucleon->GetSplitableHadron()->GetDefinition()->GetPDGMass());
>> 1103 /*
>> 1104 } else
>> 1105 {
>> 1106 Mt2 = sqr(tmp.x())+sqr(tmp.y())+ // 4.12.2010
>> 1107 aNucleon->Get4Momentum().e()*aNucleon->Get4Momentum().e();
>> 1108 }
>> 1109 */
>> 1110 G4double Xminus=tmp.z();
>> 1111
>> 1112 G4double Pz=-WminusTarget*Xminus/2. + Mt2/(2.*WminusTarget*Xminus);
>> 1113 G4double E = WminusTarget*Xminus/2. + Mt2/(2.*WminusTarget*Xminus);
>> 1114
>> 1115 tmp.setPz(Pz);
>> 1116 tmp.setE(E);
>> 1117 //G4cout<<"Targ after in CMS "<<tmp<<G4endl;
>> 1118 tmp.transform(toLab);
>> 1119
>> 1120 aNucleon->SetMomentum(tmp);
>> 1121 //G4cout<<"Targ after in LAB "<<aNucleon->Get4Momentum()<<G4endl;
>> 1122 G4VSplitableHadron * targetSplitable=aNucleon->GetSplitableHadron();
>> 1123 targetSplitable->Set4Momentum(tmp);
>> 1124
>> 1125 } // end of for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
2187 1126
2188 if ( FirstString != 0 ) strings->push_ << 1127 //G4cout<<"ResidualMassNumber and Mom "<<ResidualMassNumber<<" "<<Residual3Momentum<<G4endl;
2189 if ( SecondString != 0 ) strings->push_ << 1128 G4double Mt2Residual=sqr(ResidualMass) +
2190 } << 1129 sqr(Residual3Momentum.x())+sqr(Residual3Momentum.y());
2191 } << 1130 //==========================
>> 1131 //G4cout<<"WminusTarget Residual3Momentum.z() "<<WminusTarget<<" "<<Residual3Momentum.z()<<G4endl;
>> 1132 G4double PzResidual=0.;
>> 1133 G4double EResidual =0.;
>> 1134 if(ResidualMassNumber != 0)
>> 1135 {
>> 1136 PzResidual=-WminusTarget*Residual3Momentum.z()/2. +
>> 1137 Mt2Residual/(2.*WminusTarget*Residual3Momentum.z());
>> 1138 EResidual = WminusTarget*Residual3Momentum.z()/2. +
>> 1139 Mt2Residual/(2.*WminusTarget*Residual3Momentum.z());
>> 1140 }
>> 1141 //==========================
>> 1142 Residual4Momentum.setPx(Residual3Momentum.x());
>> 1143 Residual4Momentum.setPy(Residual3Momentum.y());
>> 1144 Residual4Momentum.setPz(PzResidual);
>> 1145 Residual4Momentum.setE(EResidual);
>> 1146 //G4cout<<"Residual4Momentum in CMS Y "<<Residual4Momentum.beta()<<G4endl;
>> 1147 //G4int Uzhi; G4cin>>Uzhi;
>> 1148 Residual4Momentum.transform(toLab);
>> 1149 //G4cout<<"Residual4Momentum in Lab "<<Residual4Momentum<<G4endl;
>> 1150 //G4int Uzhi; G4cin>>Uzhi;
>> 1151 //-------------------------------------------------------------
>> 1152 return true;
>> 1153 }
2192 1154
2193 #ifdef debugBuildString << 1155 // ------------------------------------------------------------
2194 G4cout << "Building of target-like strings" << 1156 G4bool G4FTFModel::ExciteParticipants()
2195 #endif << 1157 {
2196 << 1158 //G4cout<<"G4FTFModel::ExciteParticipants() "<<G4endl;
2197 G4bool isProjectile = false; << 1159 G4bool Successfull(true); //(false); // 1.07.11
2198 for ( G4int ahadron = 0; ahadron < NumberOf <<
2199 G4VSplitableHadron* aNucleon = TheInvolve <<
2200 <<
2201 #ifdef debugBuildString <<
2202 G4cout << "Nucleon #, status, intCount " <<
2203 << aNucleon->GetStatus() << " " << <<
2204 #endif <<
2205 <<
2206 FirstString = 0 ; SecondString = 0; <<
2207 <<
2208 if ( aNucleon->GetStatus() == 0 ) { // A <<
2209 theExcitation->CreateStrings( aNucleon, <<
2210 FirstStri <<
2211 NumberOfTargetSpectatorNucleons--; <<
2212 <<
2213 #ifdef debugBuildString <<
2214 G4cout << " 1 case A string is build" < <<
2215 #endif <<
2216 <<
2217 } else if ( aNucleon->GetStatus() == 1 & <<
2218 // A nucleon took part in diffractive i <<
2219 theExcitation->CreateStrings( aNucleon, <<
2220 FirstStri <<
2221 <<
2222 #ifdef debugBuildString <<
2223 G4cout << " 2 case A string is build, n <<
2224 #endif <<
2225 <<
2226 NumberOfTargetSpectatorNucleons--; <<
2227 <<
2228 } else if ( aNucleon->GetStatus() == 1 & <<
2229 HighEnergyInter ) { <<
2230 // A nucleon was considered as a partic <<
2231 // its interactions were skipped. It wi <<
2232 // at high energies. <<
2233 <<
2234 G4LorentzVector ParticleMomentum = aNuc <<
2235 G4KineticTrack* aTrack = new G4KineticT <<
2236 <<
2237 <<
2238 <<
2239 <<
2240 FirstString = new G4ExcitedString( aTra <<
2241 <<
2242 #ifdef debugBuildString <<
2243 G4cout << "3 case A string is build" << <<
2244 #endif <<
2245 <<
2246 } else if ( aNucleon->GetStatus() == 1 & <<
2247 ! HighEnergyInter ) { <<
2248 // A nucleon was considered as a partic <<
2249 // its interactions were skipped. It wi <<
2250 // at low energies energies. <<
2251 aNucleon->SetStatus( 5 ); // 4->5 <<
2252 // ??? delete aNucleon; <<
2253 <<
2254 #ifdef debugBuildString <<
2255 G4cout << "4 case A string is not build <<
2256 #endif <<
2257 <<
2258 } else if ( aNucleon->GetStatus() == 2 | <<
2259 aNucleon->GetStatus() == 3 ) <<
2260 G4LorentzVector ParticleMomentum = aNuc <<
2261 G4KineticTrack* aTrack = new G4KineticT <<
2262 <<
2263 <<
2264 <<
2265 FirstString = new G4ExcitedString( aTra <<
2266 <<
2267 #ifdef debugBuildString <<
2268 G4cout << "5 case A string is build" << <<
2269 #endif <<
2270 <<
2271 if ( aNucleon->GetStatus() == 2 ) Numbe <<
2272 <<
2273 } else { <<
2274 <<
2275 #ifdef debugBuildString <<
2276 G4cout << "6 case No string" << G4endl; <<
2277 #endif <<
2278 1160
2279 } << 1161 theParticipants.StartLoop();
>> 1162 G4int CurrentInteraction(0); // Uzhi Feb26
2280 1163
2281 if ( FirstString != 0 ) strings->push_ba << 1164 G4int MaxNumOfInelCollisions=G4int(theParameters->GetMaxNumberOfCollisions());
2282 if ( SecondString != 0 ) strings->push_ba << 1165 //G4cout<<"MaxNumOfInelCollisions "<<MaxNumOfInelCollisions<<G4endl;
>> 1166 G4double NumberOfInel(0.);
>> 1167 //
>> 1168 if(MaxNumOfInelCollisions > 0)
>> 1169 { // Plab > Pbound, Normal application of FTF is possible
>> 1170 G4double ProbMaxNumber=theParameters->GetMaxNumberOfCollisions()-
>> 1171 MaxNumOfInelCollisions;
>> 1172 if(G4UniformRand() < ProbMaxNumber) {MaxNumOfInelCollisions++;}
>> 1173 NumberOfInel=MaxNumOfInelCollisions;
>> 1174 //G4cout<<"Plab > Pbound MaxNumOfInelCollisions "<<MaxNumOfInelCollisions<<G4endl;
>> 1175 } else
>> 1176 { // Plab < Pbound, Normal application of FTF is impossible,
>> 1177 // low energy corrections applied.
>> 1178 if(theParticipants.theNucleus->GetMassNumber() > 1)
>> 1179 {
>> 1180 NumberOfInel = theParameters->GetProbOfInteraction();
>> 1181 MaxNumOfInelCollisions = 1;
>> 1182 } else
>> 1183 { // Special case for hadron-nucleon interactions
>> 1184 NumberOfInel = 1.;
>> 1185 MaxNumOfInelCollisions = 1;
>> 1186 //G4cout<<"Plab < Pbound MaxNumOfInelCollisions "<<MaxNumOfInelCollisions<<G4endl;
>> 1187 }
>> 1188 } // end of if(MaxNumOfInelCollisions > 0)
>> 1189 //
>> 1190 //G4cout<<"MaxNumOfInelCollisions MaxNumOfInelCollisions "<<MaxNumOfInelCollisions<<" "<<MaxNumOfInelCollisions<<G4endl;
2283 1191
2284 } << 1192 while (theParticipants.Next())
>> 1193 {
>> 1194 CurrentInteraction++;
>> 1195 const G4InteractionContent & collision=theParticipants.GetInteraction();
2285 1196
2286 #ifdef debugBuildString << 1197 G4VSplitableHadron * projectile=collision.GetProjectile();
2287 G4cout << G4endl << "theAdditionalString.si << 1198 G4VSplitableHadron * target=collision.GetTarget();
2288 << G4endl << G4endl; << 1199 //
2289 #endif << 1200 //G4cout<<"Ppr tr "<<projectile<<" "<<target<<G4endl;
2290 << 1201 //G4cout<<"theInter Time "<<collision.GetInteractionTime()/fermi<<G4endl;
2291 isProjectile = true; << 1202 //G4cout<<"Int Status "<<collision.GetStatus()<<" "<<CurrentInteraction<<G4endl;
2292 if ( theAdditionalString.size() != 0 ) { << 1203 //G4cout<<"Proj M "<<projectile->Get4Momentum()<<" "<<projectile->Get4Momentum().mag()<<G4endl;
2293 for ( unsigned int ahadron = 0; ahadron << 1204 //G4cout<<"Targ M "<<target->Get4Momentum()<<" "<<target->Get4Momentum().mag()<<G4endl;
2294 //if ( theAdditionalString[ ahadron ]-> << 1205 //G4cout<<"ProbabilityOfElasticScatt "<<theParameters->GetProbabilityOfElasticScatt()<<G4endl;
2295 FirstString = 0; SecondString = 0; << 1206 //G4cout<<"ProbabilityOfAnnihilation "<<theParameters->GetProbabilityOfAnnihilation()<<G4endl;
2296 theExcitation->CreateStrings( theAdditi << 1207 //G4cout<<"projectile->GetStatus target->GetStatus "<<projectile->GetStatus()<<" "<<target->GetStatus()<<G4endl;
2297 FirstStri << 1208 //if((projectile->GetStatus() == 1) && (target->GetStatus() ==1))
2298 if ( FirstString != 0 ) strings->push_ << 1209 //
2299 if ( SecondString != 0 ) strings->push_ << 1210 if(collision.GetStatus()) // Uzhi Feb26
2300 } << 1211 {
2301 } <<
2302 1212
2303 //for ( unsigned int ahadron = 0; ahadron < << 1213 //theParameters->SetProbabilityOfElasticScatt(1.);
2304 // G4cout << ahadron << " " << strings->op << 1214 //G4cout<<"before pro "<<projectile->Get4Momentum()<<" "<<projectile->Get4Momentum().mag()<<G4endl;
2305 // << " " << strings->operator[]( a << 1215 //G4cout<<"before tar "<<target->Get4Momentum()<<" "<<target->Get4Momentum().mag()<<G4endl;
2306 //} << 1216 //G4cout<<"Prob el "<<theParameters->GetProbabilityOfElasticScatt()<<G4endl;
2307 //G4cout << "------------------------" << G << 1217 //G4cout<<"Prob an "<<theParameters->GetProbabilityOfAnnihilation()<<G4endl;
>> 1218
>> 1219 //G4cout<<"Pr Tr "<<projectile->GetDefinition()->GetPDGEncoding()<<" "<<target->GetDefinition()->GetPDGEncoding()<<G4endl;
>> 1220 //G4int Uzhi; G4cin>>Uzhi;
>> 1221
>> 1222 if(G4UniformRand()< theParameters->GetProbabilityOfElasticScatt())
>> 1223 { // Elastic scattering -------------------------
>> 1224 //G4cout<<"Elastic FTF"<<G4endl;
>> 1225 Successfull = theElastic->ElasticScattering(projectile, target, theParameters)
>> 1226 || Successfull; // 9.06.2012
>> 1227 }
>> 1228 else if(G4UniformRand() > theParameters->GetProbabilityOfAnnihilation())
>> 1229 { // Inelastic scattering ----------------------
>> 1230 //G4cout<<"Inelastic FTF"<<G4endl;
>> 1231 //G4cout<<"NumberOfInel MaxNumOfInelCollisions "<<NumberOfInel<<" "<<MaxNumOfInelCollisions<<G4endl;
>> 1232 if(G4UniformRand()< NumberOfInel/MaxNumOfInelCollisions)
>> 1233 {
>> 1234 if(theExcitation->ExciteParticipants(projectile, target,
>> 1235 theParameters, theElastic))
>> 1236 {
>> 1237 NumberOfInel--;
>> 1238 //G4cout<<"Excitation FTF Successfull "<<G4endl;
>> 1239 //G4cout<<"After pro "<<projectile->Get4Momentum()<<" "<<projectile->Get4Momentum().mag()<<G4endl;
>> 1240 //G4cout<<"After tar "<<target->Get4Momentum()<<" "<<target->Get4Momentum().mag()<<G4endl;
>> 1241 } else
>> 1242 {
>> 1243 //G4cout<<"Excitation FTF Non Successfull -> Elastic scattering "<<Successfull<<G4endl;
>> 1244
>> 1245 Successfull = theElastic->ElasticScattering(projectile, target, theParameters)
>> 1246 || Successfull; // 9.06.2012
>> 1247
>> 1248 /*
>> 1249 if(NumberOfInvolvedTargetNucleon > 1)
>> 1250 {
>> 1251 NumberOfInvolvedTargetNucleon--;
>> 1252 target->SetStatus(0); // 1->0 return nucleon to the target VU 10.04.2012
>> 1253 }
>> 1254 */
>> 1255 }
>> 1256 } else // If NumOfInel
>> 1257 {
>> 1258 //G4cout<<"Elastic at rejected inelastic scattering"<<G4endl;
>> 1259 Successfull = theElastic->ElasticScattering(projectile, target, theParameters)
>> 1260 || Successfull; // 9.06.2012
>> 1261 /*
>> 1262 if(theElastic->ElasticScattering(projectile, target, theParameters))
>> 1263 {
>> 1264 // Successfull = Successfull || true;
>> 1265 } else
>> 1266 {
>> 1267 // Successfull = Successfull || false;
>> 1268 //Successfull = Successfull && false;
>> 1269 Successfull = false; break; // 1.07.11
>> 1270
>> 1271 if(NumberOfInvolvedTargetNucleon > 1)
>> 1272 {
>> 1273 NumberOfInvolvedTargetNucleon--;
>> 1274 target->SetStatus(4); // 1->0 return nucleon to the target VU 10.04.2012
>> 1275 }
>> 1276 }
>> 1277 */
>> 1278 } // end if NumOfInel
>> 1279 }
>> 1280 else // Annihilation
>> 1281 {
>> 1282 //G4cout<<"Annihilation"<<G4endl;
>> 1283 //G4cout<<"Before pro "<<projectile->Get4Momentum()<<G4endl;
>> 1284 //G4cout<<"Before tar "<<target->Get4Momentum()<<G4endl;
>> 1285 //G4cout<<"Mom pro "<<theProjectile.GetTotalMomentum()<<G4endl;
>> 1286 //if(theProjectile.GetTotalMomentum() < 2000.*MeV)
>> 1287 {
>> 1288 while (theParticipants.Next())
>> 1289 {
>> 1290 G4InteractionContent & acollision=theParticipants.GetInteraction();//Uzhi Feb26
>> 1291
>> 1292 G4VSplitableHadron * NextProjectileNucleon=acollision.GetProjectile(); // Uzhi Feb26
>> 1293 G4VSplitableHadron * NextTargetNucleon =acollision.GetTarget();
>> 1294 if((projectile == NextProjectileNucleon) ||
>> 1295 (target == NextTargetNucleon )) acollision.SetStatus(0);
>> 1296 // if(target != NextTargetNucleon) NextTargetNucleon->SetStatus(0); // Uzhi Feb26
>> 1297 }
2308 1298
2309 return; << 1299 theParticipants.StartLoop();
>> 1300 for(G4int I=0; I < CurrentInteraction; I++) theParticipants.Next();
>> 1301
>> 1302 //-----------------------------------------
>> 1303 // 1Nov2011 AjustTargetNucleonForAnnihilation(projectile, target);
>> 1304 //-----------------------------------------
>> 1305 //G4cout<<"After Ajsd pro "<<projectile->Get4Momentum()<<G4endl;
>> 1306 //G4cout<<"After Ajst tar "<<target->Get4Momentum()<<G4endl;
2310 } 1307 }
>> 1308 G4VSplitableHadron *AdditionalString=0;
>> 1309 if(theAnnihilation->Annihilate(projectile, target, AdditionalString, theParameters))
>> 1310 {
>> 1311 Successfull = Successfull || true;
>> 1312 //G4cout<<G4endl<<"*AdditionalString "<<AdditionalString<<G4endl;
>> 1313 //G4cout<<"After pro "<<projectile->Get4Momentum()<<G4endl;
>> 1314 //G4cout<<"After tar "<<target->Get4Momentum()<<G4endl;
>> 1315
>> 1316 if(AdditionalString != 0) theAdditionalString.push_back(AdditionalString);
>> 1317
>> 1318 // break;
>> 1319
>> 1320 } else
>> 1321 {
>> 1322 //A.R. 25-Jul-2012 : commenting the next line to fix a Coverity
>> 1323 // "logically dead code".
>> 1324 //Successfull = Successfull || false;
2311 1325
>> 1326 // target->SetStatus(2);
>> 1327 }
>> 1328 }
>> 1329 //
>> 1330 } // End of if((projectile->GetStatus() == 1) && (target->GetStatus() ==1))
2312 1331
2313 //=========================================== << 1332 } // end of while (theParticipants.Next())
2314 <<
2315 void G4FTFModel::GetResiduals() { <<
2316 // This method is needed for the correct ap <<
2317 <<
2318 #ifdef debugFTFmodel <<
2319 G4cout << "GetResiduals(): HighEnergyInter? <<
2320 << HighEnergyInter << " " << GetProj <<
2321 #endif <<
2322 <<
2323 if ( HighEnergyInter ) { <<
2324 <<
2325 #ifdef debugFTFmodel <<
2326 G4cout << "NumberOfInvolvedNucleonsOfTarg <<
2327 #endif <<
2328 <<
2329 G4double DeltaExcitationE = TargetResidua <<
2330 G4double( Num <<
2331 G4LorentzVector DeltaPResidualNucleus = T <<
2332 G <<
2333 <<
2334 for ( G4int i = 0; i < NumberOfInvolvedNu <<
2335 G4Nucleon* aNucleon = TheInvolvedNucleo <<
2336 <<
2337 #ifdef debugFTFmodel <<
2338 G4VSplitableHadron* targetSplitable = a <<
2339 G4cout << i << " Hit? " << aNucleon->Ar <<
2340 if ( targetSplitable ) G4cout << i << " <<
2341 #endif <<
2342 <<
2343 G4LorentzVector tmp = -DeltaPResidualNu <<
2344 aNucleon->SetMomentum( tmp ); <<
2345 aNucleon->SetBindingEnergy( DeltaExcita <<
2346 } <<
2347 <<
2348 if ( TargetResidualMassNumber != 0 ) { <<
2349 G4ThreeVector bstToCM = TargetResidual4 <<
2350 1333
2351 G4V3DNucleus* theTargetNucleus = GetTar << 1334 //Successfull=true;
2352 G4LorentzVector residualMomentum( 0.0, << 1335 //G4cout<<"G4FTFModel::ExciteParticipants() Successfull "<<Successfull<<G4endl;
2353 G4Nucleon* aNucleon = 0; << 1336 //G4int Uzhi; G4cin>>Uzhi;
2354 theTargetNucleus->StartLoop(); << 1337 return Successfull;
2355 while ( ( aNucleon = theTargetNucleus-> << 1338 }
2356 if ( ! aNucleon->AreYouHit() ) { <<
2357 G4LorentzVector tmp = aNucleon->Get <<
2358 aNucleon->SetMomentum( tmp ); <<
2359 residualMomentum += tmp; <<
2360 } <<
2361 } <<
2362 1339
2363 residualMomentum /= TargetResidualMassN << 1340 //-------------------------------------------------------------------
>> 1341 void G4FTFModel::AjustTargetNucleonForAnnihilation(G4VSplitableHadron *SelectedAntiBaryon,
>> 1342 G4VSplitableHadron *SelectedTargetNucleon)
>> 1343 {
>> 1344 G4LorentzVector Pparticipants=SelectedAntiBaryon->Get4Momentum()+
>> 1345 SelectedTargetNucleon->Get4Momentum();
2364 1346
2365 G4double Mass = TargetResidual4Momentum << 1347 G4V3DNucleus *theNucleus = GetWoundedNucleus();
2366 G4double SumMasses = 0.0; << 1348 //G4cout<<"Init A mass "<<theNucleus->GetMass()<<" "<<theNucleus->GetMassNumber()<<" "<<theNucleus->GetCharge()<<G4endl;
2367 <<
2368 aNucleon = 0; <<
2369 theTargetNucleus->StartLoop(); <<
2370 while ( ( aNucleon = theTargetNucleus-> <<
2371 if ( ! aNucleon->AreYouHit() ) { <<
2372 G4LorentzVector tmp = aNucleon->Get <<
2373 G4double E = std::sqrt( tmp.vect(). <<
2374 sqr( aNucle <<
2375 tmp.setE( E ); aNucleon->SetMoment <<
2376 SumMasses += E; <<
2377 } <<
2378 } <<
2379 1349
2380 G4double Chigh = Mass / SumMasses; G4do << 1350 ResidualExcitationEnergy=0.;
2381 const G4int maxNumberOfLoops = 1000; << 1351 G4int ResidualCharge =theNucleus->GetCharge();
2382 G4int loopCounter = 0; << 1352 G4int ResidualMassNumber=theNucleus->GetMassNumber();
2383 do { << 1353
2384 C = ( Chigh + Clow ) / 2.0; << 1354 G4ThreeVector P3nuclearResidual(0.,0.,0.);
2385 SumMasses = 0.0; << 1355 G4LorentzVector PnuclearResidual(0.,0.,0.,0.);
2386 aNucleon = 0; << 1356
2387 theTargetNucleus->StartLoop(); << 1357
2388 while ( ( aNucleon = theTargetNucleus << 1358 G4double ExcitationEnergyPerWoundedNucleon=
2389 if ( ! aNucleon->AreYouHit() ) { << 1359 theParameters->GetExcitationEnergyPerWoundedNucleon();
2390 G4LorentzVector tmp = aNucleon->G << 1360 //-------
2391 G4double E = std::sqrt( tmp.vect( << 1361 G4Nucleon * aNucleon;
2392 sqr( aNuc << 1362 theNucleus->StartLoop();
2393 SumMasses += E; << 1363 G4int NumberOfHoles(0);
>> 1364 //G4cout<<"Start loop"<<G4endl;
>> 1365 while ((aNucleon = theNucleus->GetNextNucleon()))
>> 1366 {
>> 1367 G4int CurrentStatus=0;
>> 1368 if(aNucleon->AreYouHit())
>> 1369 {
>> 1370 if(aNucleon->GetSplitableHadron() == SelectedTargetNucleon)
>> 1371 {CurrentStatus=1;}
>> 1372 else
>> 1373 {
>> 1374 if(aNucleon->GetSplitableHadron()->GetSoftCollisionCount() == 0)
>> 1375 {CurrentStatus=0;}
>> 1376 else {CurrentStatus=1;}
2394 } 1377 }
>> 1378 }
>> 1379 //G4cout<<"CurrentStatus "<<CurrentStatus<<G4endl;
>> 1380 if(CurrentStatus != 0)
>> 1381 { // Participating nucleons
>> 1382 //G4cout<<" Partic "<<aNucleon->GetSplitableHadron()->GetStatus()<<" "<<aNucleon->GetSplitableHadron()->GetSoftCollisionCount()<<G4endl;
>> 1383 NumberOfHoles++;
>> 1384 ResidualExcitationEnergy+=ExcitationEnergyPerWoundedNucleon;
>> 1385 ResidualCharge-=(G4int) aNucleon->GetDefinition()->GetPDGCharge();
>> 1386 ResidualMassNumber--;
>> 1387 }
>> 1388 else
>> 1389 { // Spectator nucleon
>> 1390 PnuclearResidual+=aNucleon->Get4Momentum();
>> 1391 }
>> 1392 } // end of while (theNucleus->GetNextNucleon())
>> 1393
>> 1394 //G4cout<<"Res Z M "<<ResidualCharge<<" "<<ResidualMassNumber<<G4endl;
>> 1395 //-------------------------------
>> 1396 G4LorentzVector Psum=Pparticipants + PnuclearResidual; // 4-momentum in CMS
>> 1397
>> 1398 // Transform momenta to cms and then rotate parallel to z axis;
>> 1399 G4LorentzRotation toCms(-1*Psum.boostVector());
>> 1400
>> 1401 G4LorentzVector Ptmp=toCms*Psum;
>> 1402
>> 1403 toCms.rotateZ(-1*Ptmp.phi());
>> 1404 toCms.rotateY(-1*Ptmp.theta());
>> 1405
>> 1406 G4LorentzRotation toLab(toCms.inverse());
>> 1407
>> 1408 //-------------------------------
>> 1409 G4double SqrtS=Psum.mag();
>> 1410 G4double S =sqr(SqrtS);
>> 1411
>> 1412 G4double ResidualMass(0.);
>> 1413 if(ResidualMassNumber != 0)
>> 1414 {
>> 1415 ResidualMass=G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(
>> 1416 ResidualCharge,ResidualMassNumber);
>> 1417 } else {return;}
>> 1418
>> 1419 //G4cout<<"Res Mass E* "<<ResidualMass<<" "<<ResidualExcitationEnergy<<G4endl;
>> 1420
>> 1421 if(ResidualMass > SqrtS) {return;}
>> 1422 else
>> 1423 {
>> 1424 if(ResidualMass+ResidualExcitationEnergy > SqrtS)
>> 1425 ResidualExcitationEnergy = SqrtS-ResidualMass;
>> 1426 }
>> 1427
>> 1428 ResidualMass+=ResidualExcitationEnergy;
>> 1429 G4double ResidualMass2=sqr(ResidualMass);
>> 1430 //G4cout<<"New Res Mass E* "<<ResidualMass<<" "<<ResidualExcitationEnergy<<G4endl;
>> 1431
>> 1432 //-------
>> 1433 G4double ParticipantMass=Pparticipants.mag();
>> 1434 G4double ParticipantMass2=sqr(ParticipantMass);
>> 1435
>> 1436 if(ResidualMass + ParticipantMass > SqrtS) ParticipantMass=SqrtS-ResidualMass;
>> 1437
>> 1438 //G4cout<<"Parts P "<<Pparticipants<<G4endl;
>> 1439 //G4cout<<"Res Nuc "<<PnuclearResidual<<G4endl;
>> 1440
>> 1441 G4double DecayMomentum2=
>> 1442 sqr(S)+sqr(ParticipantMass2)+ sqr(ResidualMass2) -
>> 1443 2.*S*ParticipantMass2 - 2.*S*ResidualMass2
>> 1444 -2.*ParticipantMass2*ResidualMass2;
>> 1445
>> 1446 if(DecayMomentum2 < 0.) return;
>> 1447
>> 1448 DecayMomentum2/=(4.*S);
>> 1449 G4double DecayMomentum = std::sqrt(DecayMomentum2);
>> 1450 //G4cout<<"DecayMomentum "<<DecayMomentum<<G4endl;
>> 1451
>> 1452 G4LorentzVector New_Pparticipants(0.,0.,DecayMomentum,
>> 1453 std::sqrt(DecayMomentum2+ParticipantMass2));
>> 1454
>> 1455 G4LorentzVector New_PnuclearResidual(0.,0.,-DecayMomentum,
>> 1456 std::sqrt(DecayMomentum2+ResidualMass2));
>> 1457
>> 1458 //G4cout<<"New part P "<<New_Pparticipants<<" "<<New_Pparticipants.mag()<<G4endl;
>> 1459 //G4cout<<"New resd P "<<New_PnuclearResidual<<" "<<New_PnuclearResidual.mag()<<G4endl;
>> 1460
>> 1461 New_Pparticipants.transform(toLab);
>> 1462 New_PnuclearResidual.transform(toLab);
>> 1463 //G4cout<<"New part P "<<New_Pparticipants<<" "<<New_Pparticipants.mag()<<G4endl;
>> 1464 //G4cout<<"New resd P "<<New_PnuclearResidual<<" "<<New_PnuclearResidual.mag()<<G4endl;
>> 1465
>> 1466 G4LorentzVector DeltaP_participants=(Pparticipants - New_Pparticipants)/2.;
>> 1467 G4LorentzVector DeltaP_nuclearResidual=(PnuclearResidual - New_PnuclearResidual)/
>> 1468 (G4double) ResidualMassNumber;
>> 1469 //------------------
>> 1470
>> 1471 Ptmp=SelectedAntiBaryon->Get4Momentum() - DeltaP_participants;
>> 1472 SelectedAntiBaryon->Set4Momentum(Ptmp);
>> 1473
>> 1474 Ptmp=SelectedTargetNucleon->Get4Momentum() - DeltaP_participants;
>> 1475 SelectedTargetNucleon->Set4Momentum(Ptmp);
>> 1476 //-----------
>> 1477
>> 1478 //A.R. 25-Jul-2012 : Fix to Covery "Division by zero"
>> 1479 //G4double DeltaExcitationEnergy=ResidualExcitationEnergy/((G4double) NumberOfHoles);
>> 1480 G4double DeltaExcitationEnergy = 0.0;
>> 1481 if ( NumberOfHoles != 0 ) {
>> 1482 DeltaExcitationEnergy = ResidualExcitationEnergy / ((G4double) NumberOfHoles);
2395 } 1483 }
2396 << 1484
2397 if ( SumMasses > Mass ) Chigh = C; << 1485 // Re-definition of the wounded nucleon momenta
2398 else Clow = C; << 1486 theNucleus->StartLoop();
2399 << 1487 while ((aNucleon = theNucleus->GetNextNucleon()))
2400 } while ( Chigh - Clow > 0.01 && << 1488 {
2401 ++loopCounter < maxNumberOfLo << 1489 G4int CurrentStatus=0;
2402 if ( loopCounter >= maxNumberOfLoops ) << 1490 if(aNucleon->AreYouHit())
2403 #ifdef debugFTFmodel << 1491 {
2404 G4cout << "BAD situation: forced exit << 1492 if(aNucleon->GetSplitableHadron() == SelectedTargetNucleon)
2405 << "\t return immediately from << 1493 {CurrentStatus=1;}
2406 #endif << 1494 else
2407 return; << 1495 {
2408 } << 1496 if(aNucleon->GetSplitableHadron()->GetSoftCollisionCount() == 0)
2409 << 1497 {CurrentStatus=0;}
2410 aNucleon = 0; << 1498 else {CurrentStatus=1;}
2411 theTargetNucleus->StartLoop(); <<
2412 while ( ( aNucleon = theTargetNucleus-> <<
2413 if ( !aNucleon->AreYouHit() ) { <<
2414 G4LorentzVector tmp = aNucleon->Get <<
2415 G4double E = std::sqrt( tmp.vect(). <<
2416 sqr( aNucle <<
2417 tmp.setE( E ); tmp.boost( -bstToCM <<
2418 aNucleon->SetMomentum( tmp ); <<
2419 } <<
2420 } <<
2421 } <<
2422 <<
2423 if ( ! GetProjectileNucleus() ) return; <<
2424 <<
2425 #ifdef debugFTFmodel <<
2426 G4cout << "NumberOfInvolvedNucleonsOfProj <<
2427 << G4endl << "ProjectileResidualEx <<
2428 << ProjectileResidualExcitationEne <<
2429 #endif <<
2430 <<
2431 DeltaExcitationE = ProjectileResidualExci <<
2432 G4double( NumberOfInvo <<
2433 DeltaPResidualNucleus = ProjectileResidua <<
2434 G4double( NumberO <<
2435 <<
2436 for ( G4int i = 0; i < NumberOfInvolvedNu <<
2437 G4Nucleon* aNucleon = TheInvolvedNucleo <<
2438 <<
2439 #ifdef debugFTFmodel <<
2440 G4VSplitableHadron* projSplitable = aNu <<
2441 G4cout << i << " Hit? " << aNucleon->Ar <<
2442 if ( projSplitable ) G4cout << i << "St <<
2443 #endif <<
2444 <<
2445 G4LorentzVector tmp = -DeltaPResidualNu <<
2446 aNucleon->SetMomentum( tmp ); <<
2447 aNucleon->SetBindingEnergy( DeltaExcita <<
2448 } <<
2449 <<
2450 if ( ProjectileResidualMassNumber != 0 ) <<
2451 G4ThreeVector bstToCM = ProjectileResid <<
2452 <<
2453 G4V3DNucleus* theProjectileNucleus = Ge <<
2454 G4LorentzVector residualMomentum( 0.0, <<
2455 G4Nucleon* aNucleon = 0; <<
2456 theProjectileNucleus->StartLoop(); <<
2457 while ( ( aNucleon = theProjectileNucle <<
2458 if ( ! aNucleon->AreYouHit() ) { <<
2459 G4LorentzVector tmp = aNucleon->Get <<
2460 aNucleon->SetMomentum( tmp ); <<
2461 residualMomentum += tmp; <<
2462 } <<
2463 } <<
2464 <<
2465 residualMomentum /= ProjectileResidualM <<
2466 <<
2467 G4double Mass = ProjectileResidual4Mome <<
2468 G4double SumMasses= 0.0; <<
2469 <<
2470 aNucleon = 0; <<
2471 theProjectileNucleus->StartLoop(); <<
2472 while ( ( aNucleon = theProjectileNucle <<
2473 if ( ! aNucleon->AreYouHit() ) { <<
2474 G4LorentzVector tmp = aNucleon->Get <<
2475 G4double E=std::sqrt( tmp.vect().ma <<
2476 sqr(aNucleon- <<
2477 tmp.setE( E ); aNucleon->SetMoment <<
2478 SumMasses += E; <<
2479 } <<
2480 } <<
2481 <<
2482 G4double Chigh = Mass / SumMasses; G4do <<
2483 const G4int maxNumberOfLoops = 1000; <<
2484 G4int loopCounter = 0; <<
2485 do { <<
2486 C = ( Chigh + Clow ) / 2.0; <<
2487 <<
2488 SumMasses = 0.0; <<
2489 aNucleon = 0; <<
2490 theProjectileNucleus->StartLoop(); <<
2491 while ( ( aNucleon = theProjectileNuc <<
2492 if ( ! aNucleon->AreYouHit() ) { <<
2493 G4LorentzVector tmp = aNucleon->G <<
2494 G4double E = std::sqrt( tmp.vect( <<
2495 sqr( aNuc <<
2496 SumMasses += E; <<
2497 } 1499 }
2498 } << 1500 }
>> 1501 //G4cout<<"CurrentStatus "<<CurrentStatus<<G4endl;
>> 1502 if(CurrentStatus != 0)
>> 1503 { // Participating nucleons
>> 1504 //G4cout<<" Partic "<<aNucleon->GetSplitableHadron()->GetStatus()<<" "<<aNucleon->GetSplitableHadron()->GetSoftCollisionCount()<<G4endl;
>> 1505 aNucleon->SetBindingEnergy(DeltaExcitationEnergy);
>> 1506 }
>> 1507 else
>> 1508 { // Spectator nucleon of nuclear residual
>> 1509 Ptmp=aNucleon->Get4Momentum() - DeltaP_nuclearResidual;
>> 1510 aNucleon->SetMomentum(Ptmp);
>> 1511 }
>> 1512 } // end of while (theNucleus->GetNextNucleon())
2499 1513
2500 if ( SumMasses > Mass) Chigh = C; << 1514 //-------------------------------
2501 else Clow = C; <<
2502 <<
2503 } while ( Chigh - Clow > 0.01 && <<
2504 ++loopCounter < maxNumberOfLo <<
2505 if ( loopCounter >= maxNumberOfLoops ) <<
2506 #ifdef debugFTFmodel <<
2507 G4cout << "BAD situation: forced exit <<
2508 << "\t return immediately from <<
2509 #endif <<
2510 return; 1515 return;
2511 } <<
2512 <<
2513 aNucleon = 0; <<
2514 theProjectileNucleus->StartLoop(); <<
2515 while ( ( aNucleon = theProjectileNucle <<
2516 if ( ! aNucleon->AreYouHit() ) { <<
2517 G4LorentzVector tmp = aNucleon->Get <<
2518 G4double E = std::sqrt( tmp.vect(). <<
2519 sqr( aNucle <<
2520 tmp.setE( E ); tmp.boost( -bstToCM <<
2521 aNucleon->SetMomentum( tmp ); <<
2522 } <<
2523 } <<
2524 } // End of if ( ProjectileResidualMass <<
2525 <<
2526 #ifdef debugFTFmodel <<
2527 G4cout << "End projectile" << G4endl; <<
2528 #endif <<
2529 <<
2530 } else { // Related to the condition: if ( <<
2531 <<
2532 #ifdef debugFTFmodel <<
2533 G4cout << "Low energy interaction: Target <<
2534 << "Tr ResidualMassNumber Tr Resid <<
2535 << TargetResidualMassNumber << " " <<
2536 << TargetResidualExcitationEnergy <<
2537 #endif <<
2538 <<
2539 G4int NumberOfTargetParticipant( 0 ); <<
2540 for ( G4int i = 0; i < NumberOfInvolvedNu <<
2541 G4Nucleon* aNucleon = TheInvolvedNucleo <<
2542 G4VSplitableHadron* targetSplitable = a <<
2543 if ( targetSplitable->GetSoftCollisionC <<
2544 } <<
2545 <<
2546 G4double DeltaExcitationE( 0.0 ); <<
2547 G4LorentzVector DeltaPResidualNucleus = G <<
2548 <<
2549 if ( NumberOfTargetParticipant != 0 ) { <<
2550 DeltaExcitationE = TargetResidualExcita <<
2551 DeltaPResidualNucleus = TargetResidual4 <<
2552 } <<
2553 <<
2554 for ( G4int i = 0; i < NumberOfInvolvedNu <<
2555 G4Nucleon* aNucleon = TheInvolvedNucleo <<
2556 G4VSplitableHadron* targetSplitable = a <<
2557 if ( targetSplitable->GetSoftCollisionC <<
2558 G4LorentzVector tmp = -DeltaPResidual <<
2559 aNucleon->SetMomentum( tmp ); <<
2560 aNucleon->SetBindingEnergy( DeltaExci <<
2561 } else { <<
2562 delete targetSplitable; <<
2563 targetSplitable = 0; <<
2564 aNucleon->Hit( targetSplitable ); <<
2565 aNucleon->SetBindingEnergy( 0.0 ); <<
2566 } <<
2567 } <<
2568 <<
2569 #ifdef debugFTFmodel <<
2570 G4cout << "NumberOfTargetParticipant " << <<
2571 << "TargetResidual4Momentum " << <<
2572 #endif <<
2573 <<
2574 if ( ! GetProjectileNucleus() ) return; <<
2575 <<
2576 #ifdef debugFTFmodel <<
2577 G4cout << "Low energy interaction: Projec <<
2578 << "Pr ResidualMassNumber Pr Resid <<
2579 << ProjectileResidualMassNumber << <<
2580 << ProjectileResidualExcitationEne <<
2581 #endif <<
2582 <<
2583 G4int NumberOfProjectileParticipant( 0 ); <<
2584 for ( G4int i = 0; i < NumberOfInvolvedNu <<
2585 G4Nucleon* aNucleon = TheInvolvedNucleo <<
2586 G4VSplitableHadron* projectileSplitable <<
2587 if ( projectileSplitable->GetSoftCollis <<
2588 } <<
2589 <<
2590 #ifdef debugFTFmodel <<
2591 G4cout << "NumberOfProjectileParticipant" <<
2592 #endif <<
2593 <<
2594 DeltaExcitationE = 0.0; <<
2595 DeltaPResidualNucleus = G4LorentzVector( <<
2596 <<
2597 if ( NumberOfProjectileParticipant != 0 ) <<
2598 DeltaExcitationE = ProjectileResidualEx <<
2599 DeltaPResidualNucleus = ProjectileResid <<
2600 } <<
2601 //G4cout << "DeltaExcitationE DeltaPResid <<
2602 // << " " << DeltaPResidualNucleus <<
2603 for ( G4int i = 0; i < NumberOfInvolvedNu <<
2604 G4Nucleon* aNucleon = TheInvolvedNucleo <<
2605 G4VSplitableHadron* projectileSplitable <<
2606 if ( projectileSplitable->GetSoftCollis <<
2607 G4LorentzVector tmp = -DeltaPResidual <<
2608 aNucleon->SetMomentum( tmp ); <<
2609 aNucleon->SetBindingEnergy( DeltaExci <<
2610 } else { <<
2611 delete projectileSplitable; <<
2612 projectileSplitable = 0; <<
2613 aNucleon->Hit( projectileSplitable ); <<
2614 aNucleon->SetBindingEnergy( 0.0 ); <<
2615 } <<
2616 } <<
2617 <<
2618 #ifdef debugFTFmodel <<
2619 G4cout << "NumberOfProjectileParticipant <<
2620 << "ProjectileResidual4Momentum " <<
2621 #endif <<
2622 <<
2623 } // End of the condition: if ( HighEnergy <<
2624 <<
2625 #ifdef debugFTFmodel <<
2626 G4cout << "End GetResiduals --------------- <<
2627 #endif <<
2628 <<
2629 } 1516 }
2630 1517
>> 1518 // ------------------------------------------------------------
>> 1519 G4ExcitedStringVector * G4FTFModel::BuildStrings()
>> 1520 {
>> 1521 // Loop over all collisions; find all primaries, and all target ( targets may
>> 1522 // be duplicate in the List ( to unique G4VSplitableHadrons)
>> 1523
>> 1524 G4ExcitedStringVector * strings;
>> 1525 strings = new G4ExcitedStringVector();
>> 1526
>> 1527 std::vector<G4VSplitableHadron *> primaries;
>> 1528
>> 1529 G4ExcitedString * FirstString(0); // If there will be a kink,
>> 1530 G4ExcitedString * SecondString(0); // two strings will be produced.
2631 1531
2632 //=========================================== << 1532 theParticipants.StartLoop(); // restart a loop
2633 << 1533 //
2634 G4ThreeVector G4FTFModel::GaussianPt( G4doubl << 1534 while ( theParticipants.Next() )
2635 << 1535 {
2636 G4double Pt2( 0.0 ), Pt( 0.0 ); << 1536 const G4InteractionContent & interaction=theParticipants.GetInteraction();
>> 1537 // if(interaction.GetStatus() != 0) // Uzhi Feb26
>> 1538 {
>> 1539 // do not allow for duplicates ...
>> 1540 if ( primaries.end() == std::find(primaries.begin(), primaries.end(),
>> 1541 interaction.GetProjectile()) )
>> 1542 primaries.push_back(interaction.GetProjectile());
>> 1543 } // Uzhi Feb26
>> 1544 }
2637 1545
2638 if (AveragePt2 > 0.0) { << 1546 //G4cout<<G4endl<<"primaries.size() "<<primaries.size()<<G4endl;
2639 const G4double ymax = maxPtSquare/Average << 1547 for (unsigned int ahadron=0; ahadron < primaries.size() ; ahadron++)
2640 if ( ymax < 200. ) { << 1548 {
2641 Pt2 = -AveragePt2 * G4Log( 1.0 + G4Unif << 1549 G4bool isProjectile(0);
2642 } else { << 1550
2643 Pt2 = -AveragePt2 * G4Log( 1.0 - G4Unif << 1551 if(primaries[ahadron]->GetStatus() <= 1) {isProjectile=true; } // VU 10.04.2012
2644 } << 1552 // if(primaries[ahadron]->GetStatus() == 3) {isProjectile=false;}
2645 Pt = std::sqrt( Pt2 ); << 1553
2646 } << 1554 FirstString=0; SecondString=0;
>> 1555 theExcitation->CreateStrings(primaries[ahadron], isProjectile,
>> 1556 FirstString, SecondString,
>> 1557 theParameters);
>> 1558
>> 1559 if(FirstString != 0) strings->push_back(FirstString);
>> 1560 if(SecondString != 0) strings->push_back(SecondString);
>> 1561 //G4cout<<"Quarks in the string in FTF"<<FirstString->GetRightParton()->GetPDGcode()<<" "<<FirstString->GetLeftParton()->GetPDGcode()<<G4endl;
2647 1562
2648 G4double phi = G4UniformRand() * twopi; << 1563 //G4cout<<FirstString<<" "<<SecondString<<G4endl;
2649 << 1564 }
2650 return G4ThreeVector( Pt*std::cos(phi), Pt* <<
2651 } <<
2652 1565
2653 //=========================================== << 1566 //G4cout<<"Check "<<strings->operator[](0)->GetRightParton()->GetPDGcode()<<" "<<strings->operator[](0)->GetLeftParton()->GetPDGcode()<<G4endl;
>> 1567 //
>> 1568 // Looking for spectator nucleons of the projectile-----------
>> 1569 G4V3DNucleus * ProjectileNucleus =theParticipants.GetProjectileNucleus();
>> 1570 if(ProjectileNucleus)
>> 1571 {
>> 1572 ProjectileNucleus->StartLoop();
>> 1573
>> 1574 G4Nucleon * ProjectileNucleon;
>> 1575 while ( (ProjectileNucleon=ProjectileNucleus->GetNextNucleon()) )
>> 1576 {
>> 1577 if ( !ProjectileNucleon->AreYouHit() )
>> 1578 { // Projectile nucleon was involved in the interaction.
>> 1579
>> 1580 G4VSplitableHadron * ProjectileSplitable=0;
>> 1581 ProjectileSplitable= new G4DiffractiveSplitableHadron(*ProjectileNucleon);
>> 1582 ProjectileNucleon->Hit(0);
>> 1583
>> 1584 G4bool isProjectile=true;
>> 1585 FirstString=0; SecondString=0;
>> 1586 theExcitation->CreateStrings(ProjectileSplitable,
>> 1587 isProjectile,
>> 1588 FirstString, SecondString,
>> 1589 theParameters);
>> 1590 if(FirstString != 0) strings->push_back(FirstString);
>> 1591 if(SecondString != 0) strings->push_back(SecondString);
2654 1592
2655 G4bool G4FTFModel:: << 1593 delete ProjectileSplitable;
2656 ComputeNucleusProperties( G4V3DNucleus* nucle <<
2657 G4LorentzVector& nu <<
2658 G4LorentzVector& re <<
2659 G4double& sumMasses <<
2660 G4double& residualE <<
2661 G4double& residualM <<
2662 G4int& residualMass <<
2663 G4int& residualChar <<
2664 <<
2665 // This method, which is called only by Put <<
2666 // - either the target nucleus (which is n <<
2667 // of hadronic interaction (hadron-nucle <<
2668 // - or the projectile nucleus or antinucl <<
2669 // or antinucleus-nucleus interaction. <<
2670 // This method assumes that the all the par <<
2671 // the action of this method consists in mo <<
2672 // first one. The return value is "false" o <<
2673 // is null. <<
2674 <<
2675 if ( ! nucleus ) return false; <<
2676 <<
2677 G4double ExcitationEnergyPerWoundedNucleon <<
2678 theParameters->GetExcitationEnergyPerWoun <<
2679 <<
2680 // Loop over the nucleons of the nucleus. <<
2681 // The nucleons that have been involved in <<
2682 // Reggeon Cascading) will be candidate to <<
2683 // All the remaining nucleons will be the n <<
2684 // The variable sumMasses is the amount of <<
2685 // 1. transverse mass of each involved <<
2686 // 2. 20.0*MeV separation energy for ea <<
2687 // 3. transverse mass of the residual n <<
2688 // In this first evaluation of sumMasses, t <<
2689 // (residualExcitationEnergy, estimated by <<
2690 // nucleon) is not taken into account. <<
2691 G4int residualNumberOfLambdas = 0; // Proj <<
2692 G4Nucleon* aNucleon = 0; <<
2693 nucleus->StartLoop(); <<
2694 while ( ( aNucleon = nucleus->GetNextNucleo <<
2695 nucleusMomentum += aNucleon->Get4Momentum <<
2696 if ( aNucleon->AreYouHit() ) { // Involv <<
2697 // Consider in sumMasses the nominal, i <<
2698 // (not the current masses, which could <<
2699 sumMasses += std::sqrt( sqr( aNucleon-> <<
2700 + aNucleon->Ge <<
2701 sumMasses += 20.0*MeV; // Separation e <<
2702 <<
2703 //residualExcitationEnergy += Excitatio <<
2704 residualExcitationEnergy += -Excitation <<
2705 <<
2706 residualMassNumber--; <<
2707 // The absolute value below is needed o <<
2708 residualCharge -= std::abs( G4int( aNuc <<
2709 } else { // Spectator nucleons <<
2710 residualMomentum += aNucleon->Get4Momen <<
2711 if ( aNucleon->GetDefinition() == G4Lam <<
2712 aNucleon->GetDefinition() == G4AntiLambd <<
2713 ++residualNumberOfLambdas; <<
2714 } <<
2715 } <<
2716 } <<
2717 #ifdef debugPutOnMassShell <<
2718 G4cout << "ExcitationEnergyPerWoundedNucleo <<
2719 << "\t Residual Charge, MassNumber ( <<
2720 << residualMassNumber << " (" << residualN <<
2721 << G4endl << "\t Initial Momentum " <<
2722 << G4endl << "\t Residual Momentum <<
2723 #endif <<
2724 residualMomentum.setPz( 0.0 ); <<
2725 residualMomentum.setE( 0.0 ); <<
2726 if ( residualMassNumber == 0 ) { <<
2727 residualMass = 0.0; <<
2728 residualExcitationEnergy = 0.0; <<
2729 } else { <<
2730 if ( residualMassNumber == 1 ) { <<
2731 if ( std::abs( residualCharge ) == 1 ) <<
2732 residualMass = G4Proton::Definition() <<
2733 } else if ( residualNumberOfLambdas == <<
2734 residualMass = G4Lambda::Definition() <<
2735 } else { <<
2736 residualMass = G4Neutron::Definition( <<
2737 } <<
2738 residualExcitationEnergy = 0.0; <<
2739 } else { <<
2740 if ( residualNumberOfLambdas > 0 ) { <<
2741 if ( residualMassNumber == 2 ) { <<
2742 residualMass = G4Lambda::Definition()->Ge <<
2743 if ( std::abs( residualCharge ) == <<
2744 residualMass += G4Proton::Definit <<
2745 } else if ( residualNumberOfLambdas == 1 <<
2746 residualMass += G4Neutron::Definition() <<
2747 } else { <<
2748 residualMass += G4Lambda::Definition()- <<
2749 } 1594 }
2750 } else { << 1595 } // End of while ( (ProjectileNucleon=ProjectileNucleus->GetNextNucleon()) )
2751 residualMass = G4HyperNucleiProperties::G << 1596 } // End of if(ProjectileNucleus)
2752 residualNumberOfLambdas ) <<
2753 } <<
2754 } else { <<
2755 residualMass = G4ParticleTable::GetPa <<
2756 GetIonMass( std::abs( residu <<
2757 } <<
2758 } <<
2759 residualMass += residualExcitationEnergy; <<
2760 } <<
2761 sumMasses += std::sqrt( sqr( residualMass ) <<
2762 return true; <<
2763 } <<
2764 <<
2765 <<
2766 //=========================================== <<
2767 <<
2768 G4bool G4FTFModel:: <<
2769 GenerateDeltaIsobar( const G4double sqrtS, <<
2770 const G4int numberOfInvo <<
2771 G4Nucleon* involvedNucle <<
2772 G4double& sumMasses ) { <<
2773 <<
2774 // This method, which is called only by Put <<
2775 // re-interpret some of the involved nucleo <<
2776 // - either by replacing a proton (2212) wi <<
2777 // - or by replacing a neutron (2112) with <<
2778 // The on-shell mass of these delta-isobars <<
2779 // the corresponding nucleon on-shell mass. <<
2780 // the max number of deltas compatible with <<
2781 // The delta-isobars are considered with th <<
2782 // corresponding nucleons. <<
2783 // This method assumes that all the paramet <<
2784 // the action of this method consists in mo <<
2785 // sumMasses. The return value is "false" o <<
2786 // have unphysical values. <<
2787 <<
2788 if ( sqrtS < 0.0 || numberOfInvolvedNucle <<
2789 <<
2790 const G4double probDeltaIsobar = 0.05; <<
2791 <<
2792 G4int maxNumberOfDeltas = G4int( (sqrtS - s <<
2793 G4int numberOfDeltas = 0; <<
2794 <<
2795 for ( G4int i = 0; i < numberOfInvolvedNucl <<
2796 <<
2797 if ( G4UniformRand() < probDeltaIsobar & <<
2798 numberOfDeltas++; <<
2799 if ( ! involvedNucleons[i] ) continue; <<
2800 // Skip any eventual lambda (that can b <<
2801 if ( involvedNucleons[i]->GetDefinition <<
2802 involvedNucleons[i]->GetDefinition() == <<
2803 G4VSplitableHadron* splitableHadron = i <<
2804 G4double massNuc = std::sqrt( sqr( spli <<
2805 + splitab <<
2806 // The absolute value below is needed i <<
2807 G4int pdgCode = std::abs( splitableHadr <<
2808 const G4ParticleDefinition* old_def = s <<
2809 G4int newPdgCode = pdgCode/10; newPdgCo <<
2810 if ( splitableHadron->GetDefinition()-> <<
2811 const G4ParticleDefinition* ptr = <<
2812 G4ParticleTable::GetParticleTable()-> <<
2813 splitableHadron->SetDefinition( ptr ); <<
2814 G4double massDelta = std::sqrt( sqr( sp <<
2815 + split <<
2816 //G4cout << i << " " << sqrtS/GeV << " <<
2817 // << " " << massNuc << G4endl; <<
2818 if ( sqrtS < sumMasses + massDelta - ma <<
2819 splitableHadron->SetDefinition( old_d <<
2820 break; <<
2821 } else { // Change is accepted <<
2822 sumMasses += ( massDelta - massNuc ); <<
2823 } <<
2824 } <<
2825 } <<
2826 <<
2827 return true; <<
2828 } <<
2829 <<
2830 <<
2831 //=========================================== <<
2832 1597
2833 G4bool G4FTFModel:: << 1598 //G4cout<<G4endl<<"theAdditionalString.size() "<<theAdditionalString.size()<<G4endl;
2834 SamplingNucleonKinematics( G4double averagePt << 1599 if(theAdditionalString.size() != 0)
2835 const G4double max << 1600 {
2836 G4double dCor, << 1601 for (unsigned int ahadron=0; ahadron < theAdditionalString.size() ; ahadron++)
2837 G4V3DNucleus* nucl << 1602 {
2838 const G4LorentzVec << 1603 G4bool isProjectile(0);
2839 const G4double res << 1604
2840 const G4int residu << 1605 if(theAdditionalString[ahadron]->GetStatus() <= 1) {isProjectile=true; } // VU 10.04.2012
2841 const G4int number << 1606 // if(theAdditionalString[ahadron]->GetStatus() == 3) {isProjectile=false;}
2842 G4Nucleon* involve << 1607
2843 G4double& mass2 ) << 1608 FirstString=0; SecondString=0;
2844 << 1609 theExcitation->CreateStrings(theAdditionalString[ahadron], isProjectile,
2845 // This method, which is called only by Put << 1610 FirstString, SecondString,
2846 // - either the target nucleons: this for << 1611 theParameters);
2847 // (hadron-nucleus, nucleus-nucleus, ant << 1612
2848 // - or the projectile nucleons or antinuc << 1613 if(FirstString != 0) strings->push_back(FirstString);
2849 // nucleus-nucleus or antinucleus-nucleu << 1614 if(SecondString != 0) strings->push_back(SecondString);
2850 // This method assumes that all the paramet << 1615 //G4cout<<"Quarks in the string in FTF"<<FirstString->GetRightParton()->GetPDGcode()<<" "<<FirstString->GetLeftParton()->GetPDGcode()<<G4endl;
2851 // the action of this method consists in ch << 1616 //G4cout<<FirstString<<" "<<SecondString<<G4endl;
2852 // whose pointers are in the vector involve << 1617 }
2853 // variable mass2. <<
2854 #ifdef debugPutOnMassShell <<
2855 G4cout << "G4FTFModel::SamplingNucleonKinem <<
2856 G4cout << " averagePt2= " << averagePt2 << <<
2857 << " dCor= " << dCor << " resMass(GeV)= " <<
2858 << " resMassN= " << residualMassNumber <<
2859 << " nNuc= " << numberOfInvolvedNucl <<
2860 << " lv= " << pResidual << G4endl; <<
2861 #endif <<
2862 <<
2863 if ( ! nucleus || numberOfInvolvedNucleon <<
2864 <<
2865 if ( residualMassNumber == 0 && numberOfI <<
2866 dCor = 0.0; <<
2867 averagePt2 = 0.0; <<
2868 } <<
2869 <<
2870 G4bool success = true; <<
2871 <<
2872 G4double SumMasses = residualMass; <<
2873 G4double invN = 1.0 / (G4double)numberOfInv <<
2874 <<
2875 // to avoid problems due to precision lost <<
2876 const G4double eps = 1.e-10; <<
2877 const G4int maxNumberOfLoops = 1000; <<
2878 G4int loopCounter = 0; <<
2879 do { <<
2880 <<
2881 success = true; <<
2882 <<
2883 // Sampling of nucleon Pt <<
2884 G4ThreeVector ptSum( 0.0, 0.0, 0.0 ); <<
2885 if( averagePt2 > 0.0 ) { <<
2886 for ( G4int i = 0; i < numberOfInvolved <<
2887 G4Nucleon* aNucleon = involvedNucleons[i]; <<
2888 if ( ! aNucleon ) continue; <<
2889 G4ThreeVector tmpPt = GaussianPt( averagePt <<
2890 ptSum += tmpPt; <<
2891 G4LorentzVector tmp( tmpPt.x(), tmpPt.y(), <<
2892 aNucleon->SetMomentum( tmp ); <<
2893 } <<
2894 } <<
2895 <<
2896 G4double deltaPx = ( ptSum.x() - pResidua <<
2897 G4double deltaPy = ( ptSum.y() - pResidua <<
2898 <<
2899 SumMasses = residualMass; <<
2900 for ( G4int i = 0; i < numberOfInvolvedNu <<
2901 G4Nucleon* aNucleon = involvedNucleons[ <<
2902 if ( ! aNucleon ) continue; <<
2903 G4double px = aNucleon->Get4Momentum(). <<
2904 G4double py = aNucleon->Get4Momentum(). <<
2905 G4double MtN = std::sqrt( sqr( aNucleon <<
2906 + sqr( px ) + <<
2907 SumMasses += MtN; <<
2908 G4LorentzVector tmp( px, py, 0.0, MtN); <<
2909 aNucleon->SetMomentum( tmp ); <<
2910 } <<
2911 <<
2912 // Sampling X of nucleon <<
2913 G4double xSum = 0.0; <<
2914 <<
2915 for ( G4int i = 0; i < numberOfInvolvedNu <<
2916 G4Nucleon* aNucleon = involvedNucleons[ <<
2917 if ( ! aNucleon ) continue; <<
2918 <<
2919 G4double x = 0.0; <<
2920 if( 0.0 != dCor ) { <<
2921 G4ThreeVector tmpX = GaussianPt( dCor*dCor, <<
2922 x = tmpX.x(); <<
2923 } <<
2924 x += aNucleon->Get4Momentum().e()/SumMa <<
2925 if ( x < -eps || x > 1.0 + eps ) { <<
2926 success = false; <<
2927 break; <<
2928 } <<
2929 x = std::min(1.0, std::max(x, 0.0)); <<
2930 xSum += x; <<
2931 // The energy is in the lab (instead of <<
2932 <<
2933 G4LorentzVector tmp( aNucleon->Get4Mome <<
2934 aNucleon->Get4Mome <<
2935 x, aNucleon->Get4M <<
2936 aNucleon->SetMomentum( tmp ); <<
2937 } <<
2938 <<
2939 if ( xSum < -eps || xSum > 1.0 + eps ) su <<
2940 if ( ! success ) continue; <<
2941 <<
2942 G4double delta = ( residualMassNumber == <<
2943 <<
2944 xSum = 1.0; <<
2945 mass2 = 0.0; <<
2946 for ( G4int i = 0; i < numberOfInvolvedNu <<
2947 G4Nucleon* aNucleon = involvedNucleons[ <<
2948 if ( ! aNucleon ) continue; <<
2949 G4double x = aNucleon->Get4Momentum().p <<
2950 xSum -= x; <<
2951 <<
2952 if ( residualMassNumber == 0 ) { <<
2953 if ( x <= -eps || x > 1.0 + eps ) { <<
2954 success = false; <<
2955 break; <<
2956 } <<
2957 } else { <<
2958 if ( x <= -eps || x > 1.0 + eps || xS <<
2959 success = false; <<
2960 break; <<
2961 } 1618 }
2962 } << 1619 //G4cout<<"Check "<<strings->operator[](0)->GetRightParton()->GetPDGcode()<<" "<<strings->operator[](0)->GetLeftParton()->GetPDGcode()<<G4endl;
2963 x = std::min( 1.0, std::max(x, eps) ); << 1620 //G4cout<<"Check "<<strings->operator[](1)->GetRightParton()->GetPDGcode()<<" "<<strings->operator[](1)->GetLeftParton()->GetPDGcode()<<G4endl;
2964 << 1621 //
2965 mass2 += sqr( aNucleon->Get4Momentum(). << 1622 //G4cout<<G4endl<<"NumberOfInvolvedNucleon "<<NumberOfInvolvedNucleon<<G4endl;
>> 1623 for (G4int ahadron=0; ahadron < NumberOfInvolvedNucleon ; ahadron++)
>> 1624 {
>> 1625 //G4cout<<"Nucleon status & int# "<<ahadron<<" "<<TheInvolvedNucleon[ahadron]->GetSplitableHadron()->GetStatus()<<" "<<TheInvolvedNucleon[ahadron]->GetSplitableHadron()->GetSoftCollisionCount()<<G4endl;
>> 1626 if(TheInvolvedNucleon[ahadron]->GetSplitableHadron()->GetStatus() ==4)
>> 1627 { // A nucleon is returned back to the nucleus after annihilation act for example
>> 1628 //G4cout<<" Delete 0"<<G4endl;
>> 1629 delete TheInvolvedNucleon[ahadron]->GetSplitableHadron();
>> 1630 G4VSplitableHadron *aHit=0;
>> 1631 TheInvolvedNucleon[ahadron]->Hit(aHit);
>> 1632 }
>> 1633 else if((TheInvolvedNucleon[ahadron]->GetSplitableHadron()->GetStatus() <=1) && // VU 10.04.2012
>> 1634 (TheInvolvedNucleon[ahadron]->GetSplitableHadron()->GetSoftCollisionCount() ==0))
>> 1635 { // A nucleon is returned back to the nucleus after rejected interactions
>> 1636 // due to an annihilation before
>> 1637 //G4cout<<" Delete int# 0"<<G4endl;
>> 1638 delete TheInvolvedNucleon[ahadron]->GetSplitableHadron();
>> 1639 G4VSplitableHadron *aHit=0;
>> 1640 TheInvolvedNucleon[ahadron]->Hit(aHit);
>> 1641 }
>> 1642 else if((TheInvolvedNucleon[ahadron]->GetSplitableHadron()->GetStatus() <=1) && // VU 10.04.2012
>> 1643 (TheInvolvedNucleon[ahadron]->GetSplitableHadron()->GetSoftCollisionCount() !=0))
>> 1644 { // Nucleon which participate in the interactions,
>> 1645 //G4cout<<"Taken 1 !=0"<<G4endl;
>> 1646 G4bool isProjectile=false;
>> 1647 FirstString=0; SecondString=0;
>> 1648 theExcitation->CreateStrings(
>> 1649 TheInvolvedNucleon[ahadron]->GetSplitableHadron(),
>> 1650 isProjectile,
>> 1651 FirstString, SecondString,
>> 1652 theParameters);
>> 1653 if(FirstString != 0) strings->push_back(FirstString);
>> 1654 if(SecondString != 0) strings->push_back(SecondString);
>> 1655 }
>> 1656 else if(TheInvolvedNucleon[ahadron]->GetSplitableHadron()->GetStatus() ==2)
>> 1657 { // Nucleon which was involved in the Reggeon cascading
>> 1658 //G4cout<<"Taken st 2"<<G4endl;
>> 1659 G4bool isProjectile=false;
>> 1660 FirstString=0; SecondString=0;
>> 1661 theExcitation->CreateStrings(
>> 1662 TheInvolvedNucleon[ahadron]->GetSplitableHadron(),
>> 1663 isProjectile,
>> 1664 FirstString, SecondString,
>> 1665 theParameters);
>> 1666 if(FirstString != 0) strings->push_back(FirstString);
>> 1667 if(SecondString != 0) strings->push_back(SecondString);
>> 1668 //G4cout<<FirstString<<" "<<SecondString<<G4endl;
>> 1669 }
>> 1670 else if(TheInvolvedNucleon[ahadron]->GetSplitableHadron()->GetStatus() ==3)
>> 1671 { // Nucleon which has participated in annihilation and disappiered
>> 1672 //G4cout<<"Status 3 "<<G4endl;
>> 1673 TheInvolvedNucleon[ahadron]->SetBindingEnergy(theParameters->GetExcitationEnergyPerWoundedNucleon());
>> 1674 }
>> 1675 else {}
2966 1676
2967 G4LorentzVector tmp( aNucleon->Get4Mome << 1677 }
2968 x, aNucleon->Get4M << 1678 /*
2969 aNucleon->SetMomentum( tmp ); << 1679 G4cout<<"Check "<<strings->operator[](0)->GetRightParton()->GetPDGcode()<<" "<<strings->operator[](0)->GetLeftParton()->GetPDGcode()<<G4endl;
2970 } << 1680 G4cout<<"Check "<<strings->operator[](1)->GetRightParton()->GetPDGcode()<<" "<<strings->operator[](1)->GetLeftParton()->GetPDGcode()<<G4endl;
2971 if ( ! success ) continue; << 1681 //G4cout<<"Check "<<strings->operator[](2)->GetRightParton()->GetPDGcode()<<" "<<strings->operator[](2)->GetLeftParton()->GetPDGcode()<<G4endl;
2972 xSum = std::min( 1.0, std::max(xSum, eps) << 1682
>> 1683 G4cout<<"*** "<<strings->operator[](0)->GetRightParton()<<" "<<strings->operator[](0)->GetLeftParton()<<G4endl;
>> 1684 G4cout<<"*** "<<strings->operator[](1)->GetRightParton()<<" "<<strings->operator[](1)->GetLeftParton()<<G4endl;
>> 1685 //G4cout<<"*** "<<strings->operator[](2)->GetRightParton()<<" "<<strings->operator[](2)->GetLeftParton()<<G4endl;
>> 1686 */
>> 1687 std::for_each(primaries.begin(), primaries.end(), DeleteVSplitableHadron());
>> 1688 primaries.clear();
>> 1689 /*
>> 1690 G4cout<<"*** "<<strings->operator[](0)->GetRightParton()<<" "<<strings->operator[](0)->GetLeftParton()<<G4endl;
>> 1691 G4cout<<"*** "<<strings->operator[](1)->GetRightParton()<<" "<<strings->operator[](1)->GetLeftParton()<<G4endl;
>> 1692 G4cout<<"*** "<<strings->operator[](2)->GetRightParton()<<" "<<strings->operator[](2)->GetLeftParton()<<G4endl;
>> 1693
>> 1694 G4cout<<"Check "<<strings->operator[](0)->GetRightParton()->GetPDGcode()<<" "<<strings->operator[](0)->GetLeftParton()->GetPDGcode()<<G4endl;
>> 1695 G4cout<<"Check "<<strings->operator[](1)->GetRightParton()->GetPDGcode()<<" "<<strings->operator[](1)->GetLeftParton()->GetPDGcode()<<G4endl;
>> 1696 G4cout<<"Check "<<strings->operator[](2)->GetRightParton()->GetPDGcode()<<" "<<strings->operator[](2)->GetLeftParton()->GetPDGcode()<<G4endl;
>> 1697 */
2973 1698
2974 if ( residualMassNumber > 0 ) mass2 += ( << 1699 /*
2975 << 1700 for (unsigned int ahadron=0; ahadron < strings->size() ; ahadron++)
2976 #ifdef debugPutOnMassShell << 1701 {
2977 G4cout << "success: " << success << " Mt( << 1702 G4cout<<ahadron<<" "<<strings->operator[](ahadron)->GetRightParton()->GetPDGcode()<<" "<<strings->operator[](ahadron)->GetLeftParton()->GetPDGcode()<<G4endl;
2978 << std::sqrt( mass2 )/GeV << G4endl; <<
2979 #endif <<
2980 <<
2981 } while ( ( ! success ) && <<
2982 ++loopCounter < maxNumberOfLoops <<
2983 return ( loopCounter < maxNumberOfLoops ); <<
2984 } 1703 }
2985 << 1704 G4cout<<"------------------------"<<G4endl;
2986 << 1705 */
2987 //=========================================== << 1706 //G4int Uzhi; G4cin >> Uzhi;
2988 << 1707 return strings;
2989 G4bool G4FTFModel:: <<
2990 CheckKinematics( const G4double sValue, <<
2991 const G4double sqrtS, <<
2992 const G4double projectileMas <<
2993 const G4double targetMass2, <<
2994 const G4double nucleusY, <<
2995 const G4bool isProjectileNuc <<
2996 const G4int numberOfInvolved <<
2997 G4Nucleon* involvedNucleons[ <<
2998 G4double& targetWminus, <<
2999 G4double& projectileWplus, <<
3000 G4bool& success ) { <<
3001 <<
3002 // This method, which is called only by Put <<
3003 // kinematics is acceptable or not. <<
3004 // This method assumes that all the paramet <<
3005 // notice that the input boolean parameter <<
3006 // only in the case of nucleus or antinucle <<
3007 // The action of this method consists in co <<
3008 // and setting the parameter success to fal <<
3009 // be rejeted. <<
3010 <<
3011 G4double decayMomentum2 = sqr( sValue ) + s <<
3012 - 2.0*( sValue*( <<
3013 + project <<
3014 targetWminus = ( sValue - projectileMass2 + <<
3015 / 2.0 / sqrtS; <<
3016 projectileWplus = sqrtS - targetMass2/targe <<
3017 G4double projectilePz = projectileWplus/2.0 <<
3018 G4double projectileE = projectileWplus/2.0 <<
3019 G4double projectileY = 0.5 * G4Log( (proje <<
3020 (proje <<
3021 G4double targetPz = -targetWminus/2.0 + tar <<
3022 G4double targetE = targetWminus/2.0 + tar <<
3023 G4double targetY = 0.5 * G4Log( (targetE + <<
3024 <<
3025 #ifdef debugPutOnMassShell <<
3026 G4cout << "decayMomentum2 " << decayMomentu <<
3027 << "\t targetWminus projectileWplus <<
3028 << "\t projectileY targetY " << proj <<
3029 if ( isProjectileNucleus ) { <<
3030 G4cout << "Order# of Wounded nucleon i, n <<
3031 } else { <<
3032 G4cout << "Order# of Wounded nucleon i, n <<
3033 } <<
3034 G4cout << G4endl; <<
3035 #endif <<
3036 <<
3037 for ( G4int i = 0; i < numberOfInvolvedNucl <<
3038 G4Nucleon* aNucleon = involvedNucleons[i] <<
3039 if ( ! aNucleon ) continue; <<
3040 G4LorentzVector tmp = aNucleon->Get4Momen <<
3041 G4double mt2 = sqr( tmp.x() ) + sqr( tmp. <<
3042 sqr( aNucleon->GetSplitabl <<
3043 G4double x = tmp.z(); <<
3044 G4double pz = -targetWminus*x/2.0 + mt2/( <<
3045 G4double e = targetWminus*x/2.0 + mt2/( <<
3046 if ( isProjectileNucleus ) { <<
3047 pz = projectileWplus*x/2.0 - mt2/(2.0*p <<
3048 e = projectileWplus*x/2.0 + mt2/(2.0*p <<
3049 } <<
3050 G4double nucleonY = 0.5 * G4Log( (e + pz) <<
3051 <<
3052 #ifdef debugPutOnMassShell <<
3053 if( isProjectileNucleus ) { <<
3054 G4cout << " " << i << " " << nucleonY < <<
3055 } else { <<
3056 G4cout << " " << i << " " << nucleonY < <<
3057 } <<
3058 G4cout << G4endl; <<
3059 #endif <<
3060 <<
3061 if ( std::abs( nucleonY - nucleusY ) > 2 <<
3062 ( isProjectileNucleus && targetY > <<
3063 ( ! isProjectileNucleus && project <<
3064 success = false; <<
3065 break; <<
3066 } <<
3067 } <<
3068 return true; <<
3069 } <<
3070 <<
3071 <<
3072 //=========================================== <<
3073 <<
3074 G4bool G4FTFModel:: <<
3075 FinalizeKinematics( const G4double w, <<
3076 const G4bool isProjectile <<
3077 const G4LorentzRotation& <<
3078 const G4double residualMa <<
3079 const G4int residualMassN <<
3080 const G4int numberOfInvol <<
3081 G4Nucleon* involvedNucleo <<
3082 G4LorentzVector& residual4Momen <<
3083 <<
3084 // This method, which is called only by Put <<
3085 // this method is called when we are sure t <<
3086 // acceptable. <<
3087 // This method assumes that all the paramet <<
3088 // notice that the input boolean parameter <<
3089 // only in the case of nucleus or antinucle <<
3090 // because the sign of pz (in the center-of <<
3091 // with respect to the case of a normal had <<
3092 // The action of this method consists in mo <<
3093 // (in the lab frame) and computing the res <<
3094 // frame). <<
3095 <<
3096 G4ThreeVector residual3Momentum( 0.0, 0.0, <<
3097 <<
3098 for ( G4int i = 0; i < numberOfInvolvedNucl <<
3099 G4Nucleon* aNucleon = involvedNucleons[i] <<
3100 if ( ! aNucleon ) continue; <<
3101 G4LorentzVector tmp = aNucleon->Get4Momen <<
3102 residual3Momentum -= tmp.vect(); <<
3103 G4double mt2 = sqr( tmp.x() ) + sqr( tmp. <<
3104 sqr( aNucleon->GetSplitabl <<
3105 G4double x = tmp.z(); <<
3106 G4double pz = -w * x / 2.0 + mt2 / ( 2. <<
3107 G4double e = w * x / 2.0 + mt2 / ( 2. <<
3108 // Reverse the sign of pz in the case of <<
3109 if ( isProjectileNucleus ) pz *= -1.0; <<
3110 tmp.setPz( pz ); <<
3111 tmp.setE( e ); <<
3112 tmp.transform( boostFromCmsToLab ); <<
3113 aNucleon->SetMomentum( tmp ); <<
3114 G4VSplitableHadron* splitableHadron = aNu <<
3115 splitableHadron->Set4Momentum( tmp ); <<
3116 } <<
3117 <<
3118 G4double residualMt2 = sqr( residualMass ) <<
3119 + sqr( residual3Moment <<
3120 <<
3121 #ifdef debugPutOnMassShell <<
3122 if ( isProjectileNucleus ) { <<
3123 G4cout << "Wminus Proj and residual3Momen <<
3124 } else { <<
3125 G4cout << "Wplus Targ and residual3Momen <<
3126 } <<
3127 #endif <<
3128 <<
3129 G4double residualPz = 0.0; <<
3130 G4double residualE = 0.0; <<
3131 if ( residualMassNumber != 0 ) { <<
3132 residualPz = -w * residual3Momentum.z() / <<
3133 residualMt2 / ( 2.0 * w * r <<
3134 residualE = w * residual3Momentum.z() / <<
3135 residualMt2 / ( 2.0 * w * r <<
3136 // Reverse the sign of residualPz in the <<
3137 if ( isProjectileNucleus ) residualPz *= <<
3138 } <<
3139 <<
3140 residual4Momentum.setPx( residual3Momentum. <<
3141 residual4Momentum.setPy( residual3Momentum. <<
3142 residual4Momentum.setPz( residualPz ); <<
3143 residual4Momentum.setE( residualE ); <<
3144 <<
3145 return true; <<
3146 } 1708 }
>> 1709 // ------------------------------------------------------------
>> 1710 void G4FTFModel::GetResidualNucleus()
>> 1711 { // This method is needed for the correct application of G4PrecompoundModelInterface
>> 1712 G4double DeltaExcitationE=ResidualExcitationEnergy/
>> 1713 (G4double) NumberOfInvolvedNucleon;
>> 1714 G4LorentzVector DeltaPResidualNucleus = Residual4Momentum/
>> 1715 (G4double) NumberOfInvolvedNucleon;
>> 1716
>> 1717 for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
>> 1718 {
>> 1719 G4Nucleon * aNucleon = TheInvolvedNucleon[i];
>> 1720 // G4LorentzVector tmp=aNucleon->Get4Momentum()-DeltaPResidualNucleus;
>> 1721 G4LorentzVector tmp=-DeltaPResidualNucleus;
>> 1722 aNucleon->SetMomentum(tmp);
>> 1723 aNucleon->SetBindingEnergy(DeltaExcitationE);
>> 1724 } // end of for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
3147 1725
3148 << 1726 }
3149 //=========================================== <<
3150 1727
3151 void G4FTFModel::ModelDescription( std::ostre << 1728 // ------------------------------------------------------------
3152 desc << " FTF (Fritiof) Mod << 1729 G4ThreeVector G4FTFModel::GaussianPt(G4double AveragePt2, G4double maxPtSquare) const
3153 << "The FTF model is based on the well << 1730 { // @@ this method is used in FTFModel as well. Should go somewhere common!
3154 << "model (B. Andersson et al., Nucl. << 1731
3155 << "(1987)). Its first program impleme << 1732 G4double Pt2(0.);
3156 << "by B. Nilsson-Almquist and E. Sten << 1733 if(AveragePt2 <= 0.) {Pt2=0.;}
3157 << "Comm. 43, 387 (1987)). The Fritiof << 1734 else
3158 << "that all hadron-hadron interaction << 1735 {
3159 << "reactions, h_1+h_2->h_1'+h_2' wher << 1736 Pt2 = -AveragePt2 * std::log(1. + G4UniformRand() *
3160 << "are excited states of the hadrons << 1737 (std::exp(-maxPtSquare/AveragePt2)-1.));
3161 << "mass spectra. The excited hadrons << 1738 }
3162 << "QCD-strings, and the corresponding << 1739 G4double Pt=std::sqrt(Pt2);
3163 << "fragmentation model is applied for << 1740 G4double phi=G4UniformRand() * twopi;
3164 << "their decays. << 1741
3165 << " The Fritiof model assumes that << 1742 return G4ThreeVector (Pt*std::cos(phi), Pt*std::sin(phi), 0.);
3166 << "a hadron-nucleus interaction a str <<
3167 << "from the projectile can interact w <<
3168 << "nuclear nucleons and becomes into <<
3169 << "states. The probability of multipl <<
3170 << "calculated in the Glauber approxim <<
3171 << "of secondary particles was neglect <<
3172 << "to these, the original Fritiof mod <<
3173 << "cribe a nuclear destruction and sl <<
3174 << " In order to overcome the diffic <<
3175 << "the model by the reggeon theory in <<
3176 << "nuclear desctruction (Kh. Abdel-Wa <<
3177 << "nsky, Phys. Atom. Nucl. 60, 828 (1 <<
3178 << "(1997)). Momenta of the nucleons e <<
3179 << "leus in the reggeon cascading are <<
3180 << "to a Fermi motion algorithm presen <<
3181 << "Collaboration (M.I. Adamovich et a <<
3182 << "A358, 337 (1997)). <<
3183 << " New features were also added to <<
3184 << "implemented in Geant4: a simulatio <<
3185 << "ron-nucleon scatterings, a simulat <<
3186 << "reactions like NN>NN* in hadron-nu <<
3187 << "a separate simulation of single di <<
3188 << " diffractive events. These allowed <<
3189 << "model parameter tuning a wide set <<
3190 << "data. <<
3191 } 1743 }
3192 1744
>> 1745 void G4FTFModel::ModelDescription(std::ostream& desc) const
>> 1746 {
>> 1747 desc << "please add description here" << G4endl;
>> 1748 }
3193 1749