Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
G4FTFModel.cc
Go to the documentation of this file.
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26//
27// $Id$
28// GEANT4 tag $Name: $
29//
30
31// ------------------------------------------------------------
32// GEANT 4 class implementation file
33//
34// ---------------- G4FTFModel ----------------
35// by Gunter Folger, May 1998.
36// class implementing the excitation in the FTF Parton String Model
37// ------------------------------------------------------------
38
39#include <utility>
40
41#include "G4FTFModel.hh"
42#include "G4ios.hh"
44#include "G4SystemOfUnits.hh"
45#include "G4FTFParameters.hh"
46#include "G4FTFParticipants.hh"
49#include "G4LorentzRotation.hh"
51#include "G4ParticleTable.hh"
52#include "G4IonTable.hh"
53
54// Class G4FTFModel
55
57 theExcitation(new G4DiffractiveExcitation()),
58 theElastic(new G4ElasticHNScattering()),
59 theAnnihilation(new G4FTFAnnihilation())
60{
62 theParameters=0;
63 NumberOfInvolvedNucleon=0;
64 NumberOfInvolvedNucleonOfProjectile=0;
65 SetEnergyMomentumCheckLevels(2*perCent, 150*MeV);
66}
67
68struct DeleteVSplitableHadron { void operator()(G4VSplitableHadron * aH){ delete aH;} };
69
71{
72// Because FTF model can be called for various particles
73// theParameters must be erased at the end of each call.
74// Thus the delete is also in G4FTFModel::GetStrings() method
75 if( theParameters != 0 ) delete theParameters;
76 if( theExcitation != 0 ) delete theExcitation;
77 if( theElastic != 0 ) delete theElastic;
78 if( theAnnihilation != 0 ) delete theAnnihilation;
79
80// Erasing of strings created at annihilation
81 if(theAdditionalString.size() != 0)
82 {
83 std::for_each(theAdditionalString.begin(), theAdditionalString.end(),
85 }
86 theAdditionalString.clear();
87
88// Erasing of target involved nucleons
89 if( NumberOfInvolvedNucleon != 0)
90 {
91 for(G4int i=0; i < NumberOfInvolvedNucleon; i++)
92 {
93 G4VSplitableHadron * aNucleon = TheInvolvedNucleon[i]->GetSplitableHadron();
94 if(aNucleon) delete aNucleon;
95 }
96 }
97
98// Erasing of projectile involved nucleons
99/*
100 if( NumberOfInvolvedNucleonOfProjectile != 0)
101 {
102 for(G4int i=0; i < NumberOfInvolvedNucleonOfProjectile; i++)
103 {
104 G4VSplitableHadron * aNucleon = TheInvolvedNucleonOfProjectile[i]->GetSplitableHadron();
105 if(aNucleon) delete aNucleon;
106 }
107 }
108*/
109}
110
111// ------------------------------------------------------------
112void G4FTFModel::Init(const G4Nucleus & aNucleus, const G4DynamicParticle & aProjectile)
113{
114 theProjectile = aProjectile;
115
116 G4double PlabPerParticle(0.); // Laboratory momentum Pz per particle/nucleon
117
118/*
119G4cout<<"FTF init Pro Name "<<theProjectile.GetDefinition()->GetParticleName()<<G4endl;
120G4cout<<"FTF init Pro Mass "<<theProjectile.GetMass()<<" "<<theProjectile.GetMomentum()<<G4endl;
121G4cout<<"FTF init Pro B Q "<<theProjectile.GetDefinition()->GetBaryonNumber()<<" "<<(G4int) theProjectile.GetDefinition()->GetPDGCharge()<<G4endl;
122G4cout<<"FTF init A Z "<<aNucleus.GetA_asInt()<<" "<<aNucleus.GetZ_asInt()<<G4endl;
123G4cout<<" "<<aNucleus.GetN()<<" "<<aNucleus.GetZ()<<G4endl;
124//G4int Uzhi; G4cin>>Uzhi;
125*/
126
127 theParticipants.SetProjectileNucleus(0);
128 theParticipants.Init(aNucleus.GetA_asInt(),aNucleus.GetZ_asInt());
129
130 if(std::abs(theProjectile.GetDefinition()->GetBaryonNumber()) <= 1)
131 { // Projectile is a hadron
132 PlabPerParticle=theProjectile.GetMomentum().z();
133
134// S = sqr( theProjectile.GetMass() ) + sqr( ProtonMass ) +
135// 2*ProtonMass*theProjectile.GetTotalEnergy();
136 }
137
138
139 if(theProjectile.GetDefinition()->GetBaryonNumber() > 1)
140 { // Projectile is a nucleus
141 theParticipants.InitProjectileNucleus(
142 theProjectile.GetDefinition()->GetBaryonNumber(),
143 (G4int) theProjectile.GetDefinition()->GetPDGCharge() );
144
145 G4ThreeVector BoostVector=theProjectile.GetMomentum()/theProjectile.GetTotalEnergy();
146 theParticipants.theProjectileNucleus->DoLorentzBoost(BoostVector);
147
148 PlabPerParticle=theProjectile.GetMomentum().z()/
149 theProjectile.GetDefinition()->GetBaryonNumber();
150
151// S = 2.*sqr( ProtonMass ) + 2*ProtonMass*
152// theProjectile.GetTotalEnergy()/theProjectile.GetDefinition()->GetBaryonNumber();
153 }
154
155 if(theProjectile.GetDefinition()->GetBaryonNumber() < -1)
156 { // Projectile is an anti-nucleus
157 theParticipants.InitProjectileNucleus(
158 std::abs( theProjectile.GetDefinition()->GetBaryonNumber()),
159 std::abs((G4int) theProjectile.GetDefinition()->GetPDGCharge()) );
160
161 G4ThreeVector BoostVector=theProjectile.GetMomentum()/theProjectile.GetTotalEnergy();
162
163 theParticipants.theProjectileNucleus->StartLoop();
164 G4Nucleon * aNucleon;
165 while ( (aNucleon = theParticipants.theProjectileNucleus->GetNextNucleon()) )
166 {
167 if(aNucleon->GetDefinition()->GetPDGEncoding() == 2212) // Proton
169
170 if(aNucleon->GetDefinition()->GetPDGEncoding() == 2112) // Neutron
172 } // end of while (theParticipant.theProjectileNucleus->GetNextNucleon())
173
174 theParticipants.theProjectileNucleus->DoLorentzBoost(BoostVector);
175
176 PlabPerParticle= theProjectile.GetMomentum().z()/
177 std::abs(theProjectile.GetDefinition()->GetBaryonNumber());
178
179// S = 2.*sqr( ProtonMass ) + 2*ProtonMass*
180// theProjectile.GetTotalEnergy()/
181// std::abs(theProjectile.GetDefinition()->GetBaryonNumber());
182 }
183
184// ------------------------------------------------------------------------
185 if( theParameters != 0 ) delete theParameters;
186 theParameters = new G4FTFParameters(theProjectile.GetDefinition(),
187 aNucleus.GetA_asInt(),aNucleus.GetZ_asInt(),
188 PlabPerParticle);
189//G4cout<<" End Init "<<theProjectile.GetMomentum()<<G4endl;
190// To turn on/off (1/0) elastic scattering close/open ...
191//theParameters->SetProbabilityOfElasticScatt(0.);
192//G4cout<<" etProbabilityOfElasticScatt "<<theParameters->GetProbabilityOfElasticScatt()<<G4endl;
193//G4cout<<" INIT ";
194//G4int Uzhi; G4cin>>Uzhi;
195
196 if(theAdditionalString.size() != 0)
197 {
198 std::for_each(theAdditionalString.begin(), theAdditionalString.end(),
200 }
201 theAdditionalString.clear();
202//G4cout<<" End Init theProjectile.GetMomentum()"<<theProjectile.GetMomentum()<<G4endl;
203}
204
205// ------------------------------------------------------------
207{
208 G4ExcitedStringVector * theStrings(0);
209
210 theParticipants.GetList(theProjectile,theParameters);
211// StoreInvolvedNucleon();
212//G4cout<<" GetList theProjectile.GetMomentum() GetBaryonNumber() "<<theProjectile.GetMomentum()<<" "<<theProjectile.GetDefinition()->GetBaryonNumber()<<G4endl;
213 G4bool Success(true);
214
215 if((std::abs(theProjectile.GetDefinition()->GetBaryonNumber()) <= 1) &&
216 (theProjectile.GetDefinition()->GetBaryonNumber() != -1) )
217 { // Standard variant of FTF for projectile hadron/nucleon
218//G4cout<<"Standard variant of FTF for projectile hadron/nucleon"<<G4endl;
219 ReggeonCascade();
220//G4cout<<"Success after Reggeon "<<Success<<" PutOnMasShell"<<G4endl;
221 Success=PutOnMassShell();
222//G4cout<<"Success after PutOn "<<Success<<" GetResid"<<G4endl;
223 GetResidualNucleus();
224 }
225//G4cout<<"Success after GetN "<<Success<<G4endl;
226//G4int Uzhi; G4cin>>Uzhi;
227 if(theProjectile.GetDefinition()->GetBaryonNumber() > 1)
228 { // Variant of FTF for projectile nuclei
229//G4cout<<"Variant of FTF for projectile nuclei"<<G4endl;
230 StoreInvolvedNucleon();
231 ReggeonCascade();
232 Success=PutOnMassShell();
233 GetResidualNucleus();
234 }
235
236// G4bool LowE_Anti_Ion(false);
237 if(theProjectile.GetDefinition()->GetBaryonNumber() <= -1)
238 { // Projectile is Anti-baryon or Anti-Nucleus
239//G4cout<<"Projectile is Anti-baryon or Anti-Nucleus "<<G4endl;
240//G4cout<<"Be4 Store"<<G4endl;
241 StoreInvolvedNucleon();
242 if(theProjectile.GetTotalMomentum()/
243 std::abs(theProjectile.GetDefinition()->GetBaryonNumber()) > 5000.*MeV)
244 {// High energy interaction
245//G4cout<<"High energy interaction "<<G4endl;
246//G4cout<<"Regeon "<<G4endl;
247 ReggeonCascade();
248//G4cout<<"Put on mass "<<G4endl;
249 Success=PutOnMassShell();
250//G4cout<<"Residual "<<G4endl;
251 GetResidualNucleus();
252 }
253 else
254 {
255//G4cout<<"Low energy interaction "<<G4endl;
256// LowE_Anti_Ion=true;
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;
273 }
274 }
275/*
276 if( Success )
277 {
278 if( ExciteParticipants() )
279 {
280//G4cout<<"Excite partic OK"<<G4endl;
281 theStrings = BuildStrings();
282//G4cout<<"Build String OK"<<G4endl;
283 if(LowE_Anti_Ion) GetResidualNucleusAfterAnnihilation();
284
285 if( theParameters != 0 )
286 {
287 delete theParameters;
288 theParameters=0;
289 }
290 } else // if( ExciteParticipants() )
291 { Success=false;}
292 } else // if( Success )
293 { Success=false;}
294*/
295 if(!Success)
296 {
297 // -------------- Erase the projectile ----------------
298//G4cout<<"Erase Proj"<<G4endl;
299 std::vector<G4VSplitableHadron *> primaries;
300
301 theParticipants.StartLoop(); // restart a loop
302 while ( theParticipants.Next() )
303 {
304 const G4InteractionContent & interaction=theParticipants.GetInteraction();
305
306 // do not allow for duplicates ...
307 if ( primaries.end() == std::find(primaries.begin(), primaries.end(),
308 interaction.GetProjectile()) )
309 primaries.push_back(interaction.GetProjectile());
310 }
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 }
324
325 NumberOfInvolvedNucleonOfProjectile=0;
326*/ // Maybe it will be needed latter------------------
327
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 }
335
336 NumberOfInvolvedNucleon=0;
337//G4cout<<"Go to fragmentation"<<G4endl;
338 return theStrings;
339
340}
341
342//-------------------------------------------------------------------
343void G4FTFModel::StoreInvolvedNucleon()
344{ //--- To store nucleons involved in low energy interaction -------
345if(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 }
399
400 ZcoordinateOfCurrentInteraction=aNucleon->GetPosition().z();
401//G4cout<<"ZcoordinateOfCurrentInteraction "<<ZcoordinateOfCurrentInteraction/fermi<<G4endl;
402//G4cout<<"TimeOfPreviousCollision "<<TimeOfPreviousCollision<<G4endl;
403
404 // A.R. 18-Oct-2011 : Protection needed for nuclear capture of
405 // anti-proton at rest.
406 if ( betta_z > 1.0e-10 ) {
407 TimeOfCurrentCollision=TimeOfPreviousCollision+
408 (ZcoordinateOfCurrentInteraction-ZcoordinateOfPreviousCollision)/betta_z;
409 } else {
410 TimeOfCurrentCollision=TimeOfPreviousCollision;
411 }
412
413//G4cout<<"TimeOfCurrentCollision "<<TimeOfCurrentCollision<<G4endl;
414// It is assumed that the nucleons are ordered on increasing z-coordinate ------------
415 aNucleon->GetSplitableHadron()->SetTimeOfCreation(TimeOfCurrentCollision);
416
417 ZcoordinateOfPreviousCollision=ZcoordinateOfCurrentInteraction;
418 TimeOfPreviousCollision=TimeOfCurrentCollision;
419 } // end of if(aNucleon->AreYouHit())
420 } // end of while (theParticipant.theNucleus->GetNextNucleon())
421//
422 return;
423} // end of if(std::abs(theProjectile.GetDefinition()->GetBaryonNumber()) <= 1)
424
425// The projectile is a nucleus or an anti-nucleus
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;
459
460 NumberOfInvolvedTargetNucleon=NumberOfInvolvedNucleon;
461 return;
462} // Uzhi 10 Feb. 2011
463
464//-------------------------------------------------------------------
465void 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 }
542
543 ZcoordinateOfCurrentInteraction=aNucleon->GetPosition().z();
544 TimeOfCurrentCollision=TimeOfPreviousCollision+
545 (ZcoordinateOfCurrentInteraction-ZcoordinateOfPreviousCollision)/betta_z;
546// It is assumed that the nucleons are ordered on increasing z-coordinate ------------
547 aNucleon->GetSplitableHadron()->SetTimeOfCreation(TimeOfCurrentCollision);
548
549 ZcoordinateOfPreviousCollision=ZcoordinateOfCurrentInteraction;
550 TimeOfPreviousCollision=TimeOfCurrentCollision;
551 } // end of if(aNucleon->AreYouHit())
552 } // end of while (theParticipant.theNucleus->GetNextNucleon())
553//
554// The algorithm can be improved, but it will be more complicated, and will require
555// changes in G4DiffractiveExcitation.cc and G4ElasticHNScattering.cc
556} // Uzhi 26 July 2009
557
558// ------------------------------------------------------------
559G4bool G4FTFModel::PutOnMassShell()
560{
561//G4cout<<"PutOnMassShell start "<<G4endl;
562 if(std::abs(theProjectile.GetDefinition()->GetBaryonNumber()) > 1) // !!!!
563 { // The projectile is a nucleus or an anti-nucleus
564//G4cout<<"PutOnMassShell AA "<<G4endl;
565 G4LorentzVector P_total(0.,0.,0.,0.);
566
567 G4LorentzVector P_participants(0.,0.,0.,0.);
568 G4int MultiplicityOfParticipants(0);
569
570 G4V3DNucleus * ProjectileNucleus =theParticipants.GetProjectileNucleus();
571 ProjectileNucleus->StartLoop();
572
573 G4Nucleon * ProjectileNucleon;
574 while ( (ProjectileNucleon=ProjectileNucleus->GetNextNucleon()) )
575 {
576 if ( ProjectileNucleon->AreYouHit() )
577 { // Projectile nucleon was involved in the interaction.
578 P_total+=ProjectileNucleon->Get4Momentum();
579 MultiplicityOfParticipants++;
580 P_participants+=ProjectileNucleon->Get4Momentum();
581 }
582 } // End of while ( (ProjectileNucleon=ProjectileNucleus->GetNextNucleon()) )
583//G4cout<<"MultParts mom Pr "<<MultiplicityOfParticipants<<" "<<P_participants<<G4endl;
584//------------------
585 G4int ResidualMassNumber(0);
586 G4int ResidualCharge(0);
587 ResidualExcitationEnergy=0.;
588 G4LorentzVector PnuclearResidual(0.,0.,0.,0.);
589
590 G4double ExcitationEnergyPerWoundedNucleon=
592
593 G4V3DNucleus * TargetNucleus =theParticipants.GetWoundedNucleus();
594 TargetNucleus->StartLoop();
595
596 G4Nucleon * TargetNucleon;
597 while ( (TargetNucleon=TargetNucleus->GetNextNucleon()) )
598 {
599 P_total+=TargetNucleon->Get4Momentum();
600 if ( TargetNucleon->AreYouHit() )
601 { // Target nucleon was involved in the interaction.
602 MultiplicityOfParticipants++;
603 P_participants+=TargetNucleon->Get4Momentum();
604 ResidualExcitationEnergy+=ExcitationEnergyPerWoundedNucleon;
605 }
606 else
607 {
608 ResidualMassNumber+=1;
609 ResidualCharge+=(G4int) TargetNucleon->GetDefinition()->GetPDGCharge();
610 PnuclearResidual+=TargetNucleon->Get4Momentum();
611 }
612 } // End of while ( (TargetNucleon=TargetNucleus->GetNextNucleon()) )
613
614 G4double ResidualMass(0.);
615 if(ResidualMassNumber == 0)
616 {
617 ResidualMass=0.;
618 ResidualExcitationEnergy=0.;
619 }
620 else
621 {
623 GetIonMass(ResidualCharge ,ResidualMassNumber);
624 if(ResidualMassNumber == 1) {ResidualExcitationEnergy=0.;}
625 }
626// ResidualMass +=ResidualExcitationEnergy; // Will be given after checks
627
628//G4cout<<"MultPars mom+Tr "<<MultiplicityOfParticipants<<" "<<P_participants<<G4endl;
629//G4cout<<"Res "<<ResidualMassNumber<<" "<<PnuclearResidual<<G4endl;
630//G4cout<<"P_total "<<P_total<<G4endl;
631
632//G4cout<<"Rez A Z M E* "<<ResidualMassNumber<<" "<<ResidualCharge<<" "<<ResidualMass<<" "<<ResidualExcitationEnergy<<G4endl;
633//--------------
634 G4double SqrtS=P_total.mag();
635 G4double S=P_total.mag2();
636
637// if(theProjectile.GetDefinition()->GetBaryonNumber() > 1)
638 { // For nucleus-nucleus interactions
639 G4double MassOfParticipants=P_participants.mag();
640 G4double MassOfParticipants2=sqr(MassOfParticipants);
641
642//G4cout<<"ResidualMass + MassOfParticipants "<<ResidualMass + MassOfParticipants<<G4endl;
643 if(SqrtS < ResidualMass + MassOfParticipants) {return false;}
644
645 if(SqrtS < ResidualMass+ResidualExcitationEnergy + MassOfParticipants)
646 {ResidualExcitationEnergy=0.;}
647
648 ResidualMass +=ResidualExcitationEnergy;
649//G4cout<<"Rez A Z M E* "<<ResidualMassNumber<<" "<<ResidualCharge<<" "<<ResidualMass<<" "<<ResidualExcitationEnergy<<G4endl;
650
651 G4double ResidualMass2=sqr(ResidualMass);
652
653 G4LorentzRotation toCms(-1*P_total.boostVector());
654
655 G4LorentzVector Pcms=toCms*P_participants;
656//G4cout<<"Ppart in CMS "<<Ptmp<<G4endl;
657
658 if ( Pcms.pz() <= 0. )
659 { // "String" moving backwards in CMS, abort collision !!
660 //G4cout << " abort ColliDeleteVSplitableHadronsion!! " << G4endl;
661 return false;
662 }
663
664 toCms.rotateZ(-1*Pcms.phi()); // Uzhi 5.12.09
665 toCms.rotateY(-1*Pcms.theta()); // Uzhi 5.12.09
666
667//G4cout<<"Mpa M res "<<ResidualMass<<" "<<MassOfParticipants<<" "<<SqrtS<<G4endl;
668 G4LorentzRotation toLab(toCms.inverse());
669
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=
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
790
791 ResidualMassNumber--;
792 ResidualCharge-=(G4int) aNucleon->GetDefinition()->GetPDGCharge();
793 }
794 else
795 { // Spectator nucleons
796 PnuclearResidual += aNucleon->Get4Momentum(); // Uzhi 12.06.2012
797 } // end of if(!aNucleon->AreYouHit())
798 } // end of while (theNucleus->GetNextNucleon())
799
800 Psum += Ptarget;
801 PnuclearResidual.setPz(0.); PnuclearResidual.setE(0.); // Uzhi 12.06.2012
802//G4cout<<"ResidualCharge ,ResidualMassNumber "<<ResidualCharge<<" "<<ResidualMassNumber<<" "<<Ptarget<<" "<<PnuclearResidual<<G4endl;
803
804 G4double ResidualMass(0.);
805 if(ResidualMassNumber == 0)
806 {
807 ResidualMass=0.;
808 ResidualExcitationEnergy=0.;
809 }
810 else
811 {
813 GetIonMass(ResidualCharge ,ResidualMassNumber);
814 if(ResidualMassNumber == 1) {ResidualExcitationEnergy=0.;}
815 }
816
817// ResidualMass +=ResidualExcitationEnergy; // Will be given after checks
818//G4cout<<"SumMasses And ResidualMass "<<SumMasses<<" "<<ResidualMass<<G4endl;
819 SumMasses += std::sqrt(sqr(ResidualMass)+PnuclearResidual.mag2()); // Uzhi 12.06.2012
820//G4cout<<"SumMasses + ResM "<<SumMasses<<G4endl;
821//G4cout<<"Psum "<<Psum<<G4endl;
822//-------------------------------------------------------------
823
824 G4double SqrtS=Psum.mag();
825 G4double S=Psum.mag2();
826
827//G4cout<<"SqrtS < SumMasses "<<SqrtS<<" "<<SumMasses<<G4endl;
828
829 if(SqrtS < SumMasses) {return false;} // It is impossible to simulate
830 // after putting nuclear nucleons
831 // on mass-shell
832
833 SumMasses -= std::sqrt(sqr(ResidualMass)+PnuclearResidual.mag2()); // Uzhi 12.06.2012
834 SumMasses += std::sqrt(sqr(ResidualMass+ResidualExcitationEnergy)
835 +PnuclearResidual.mag2()); // Uzhi 12.06.2012
836 if(SqrtS < SumMasses) // Uzhi 12.06.2012
837 {
838 SumMasses -= std::sqrt(sqr(ResidualMass+ResidualExcitationEnergy)
839 +PnuclearResidual.mag2()); // Uzhi 12.06.2012
840 SumMasses += std::sqrt(sqr(ResidualMass)+PnuclearResidual.mag2());// Uzhi 12.06.2012
841 ResidualExcitationEnergy=0.;
842 }
843
844 ResidualMass +=ResidualExcitationEnergy;
845// SumMasses +=ResidualExcitationEnergy; // Uzhi 12.06.2012
846//G4cout<<"ResidualMass SumMasses ResidualExcitationEnergy "<<ResidualMass<<" "<<SumMasses<<" "<<ResidualExcitationEnergy<<G4endl;
847//-------------------------------------------------------------
848// Sampling of nucleons what are transfered to delta-isobars --
849 G4int MaxNumberOfDeltas = (int)((SqrtS - SumMasses)/(400.*MeV));
850 G4int NumberOfDeltas(0);
851//G4cout<<"MaxNumberOfDeltas "<<MaxNumberOfDeltas<<G4endl;
852 if(theNucleus->GetMassNumber() != 1)
853 {
854//G4cout<<"NumberOfInvolvedNucleon "<<NumberOfInvolvedNucleon<<G4endl;
855 G4double ProbDeltaIsobar(0.05); // Uzhi 6.07.2012
856 for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
857 {
858 if((G4UniformRand() < ProbDeltaIsobar)&&(NumberOfDeltas < MaxNumberOfDeltas))
859 {
860 NumberOfDeltas++;
861 G4VSplitableHadron * targetSplitable=TheInvolvedNucleon[i]->GetSplitableHadron();
862 G4double MassDec=std::sqrt(sqr(targetSplitable->GetDefinition()->GetPDGMass())
863 + targetSplitable->Get4Momentum().perp2());
864
865 G4int PDGcode = targetSplitable->GetDefinition()->GetPDGEncoding();
866 G4ParticleDefinition* Old_def = targetSplitable->GetDefinition();
867
868 G4int newPDGcode = PDGcode/10; newPDGcode=newPDGcode*10+4; // Delta
871 targetSplitable->SetDefinition(ptr);
872 G4double MassInc=std::sqrt(sqr(targetSplitable->GetDefinition()->GetPDGMass())
873 + targetSplitable->Get4Momentum().perp2());
874 if(SqrtS < SumMasses + MassInc - MassDec) // Uzhi 12.06.2012
875 { // Change cannot be acsepted!
876 targetSplitable->SetDefinition(Old_def);
877 ProbDeltaIsobar = 0.;
878 } else
879 { // Change is acsepted.
880 SumMasses += (MassInc - MassDec);
881 }
882 }
883 } // end of for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
884 } // end of if(theNucleus.GetMassNumber() != 1)
885
886//-------------------------------------------------------------
887
888 G4LorentzRotation toCms(-1*Psum.boostVector());
889 G4LorentzVector Ptmp=toCms*Pprojectile;
890 if ( Ptmp.pz() <= 0. )
891 { // "String" moving backwards in CMS, abort collision !!
892 //G4cout << " abort ColliDeleteVSplitableHadronsion!! " << G4endl;
893 return false;
894 }
895
896 G4LorentzRotation toLab(toCms.inverse());
897
898 Ptmp=toCms*Ptarget;
899 G4double YtargetNucleus=Ptmp.rapidity();
900//G4cout<<"YtargetNucleus "<<YtargetNucleus<<G4endl;
901//-------------------------------------------------------------
902//------- Ascribing of the involved nucleons Pt and Xminus ----
903 G4double Dcor = theParameters->GetDofNuclearDestruction()/
904 theNucleus->GetMassNumber();
905
906 G4double AveragePt2 = theParameters->GetPt2ofNuclearDestruction();
907 G4double maxPtSquare = theParameters->GetMaxPt2ofNuclearDestruction();
908//G4cout<<"Dcor "<<theParameters->GetDofNuclearDestruction()<<" "<<Dcor<<" AveragePt2 "<<AveragePt2<<G4endl;
909 G4double M2target(0.);
910 G4double WminusTarget(0.);
911 G4double WplusProjectile(0.);
912
913 G4int NumberOfTries(0);
914 G4double ScaleFactor(1.);
915 G4bool OuterSuccess(true);
916 do // while (!OuterSuccess)
917 {
918 OuterSuccess=true;
919
920 do // while (SqrtS < Mprojectile + std::sqrt(M2target))
921 { // while (DecayMomentum < 0.)
922
923 NumberOfTries++;
924
925 if(NumberOfTries == 100*(NumberOfTries/100)) // 100
926 { // At large number of tries it would be better to reduce the values
927 ScaleFactor/=2.;
928 Dcor *=ScaleFactor;
929 AveragePt2 *=ScaleFactor;
930 }
931
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())+
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;
1070if((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())+
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++ )
1126
1127//G4cout<<"ResidualMassNumber and Mom "<<ResidualMassNumber<<" "<<Residual3Momentum<<G4endl;
1128 G4double Mt2Residual=sqr(ResidualMass) +
1129 sqr(Residual3Momentum.x())+sqr(Residual3Momentum.y());
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}
1154
1155// ------------------------------------------------------------
1156G4bool G4FTFModel::ExciteParticipants()
1157{
1158//G4cout<<"G4FTFModel::ExciteParticipants() "<<G4endl;
1159 G4bool Successfull(true); //(false); // 1.07.11
1160
1161 theParticipants.StartLoop();
1162 G4int CurrentInteraction(0); // Uzhi Feb26
1163
1164 G4int MaxNumOfInelCollisions=G4int(theParameters->GetMaxNumberOfCollisions());
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;
1191
1192 while (theParticipants.Next())
1193 {
1194 CurrentInteraction++;
1195 const G4InteractionContent & collision=theParticipants.GetInteraction();
1196
1197 G4VSplitableHadron * projectile=collision.GetProjectile();
1198 G4VSplitableHadron * target=collision.GetTarget();
1199//
1200//G4cout<<"Ppr tr "<<projectile<<" "<<target<<G4endl;
1201//G4cout<<"theInter Time "<<collision.GetInteractionTime()/fermi<<G4endl;
1202//G4cout<<"Int Status "<<collision.GetStatus()<<" "<<CurrentInteraction<<G4endl;
1203//G4cout<<"Proj M "<<projectile->Get4Momentum()<<" "<<projectile->Get4Momentum().mag()<<G4endl;
1204//G4cout<<"Targ M "<<target->Get4Momentum()<<" "<<target->Get4Momentum().mag()<<G4endl;
1205//G4cout<<"ProbabilityOfElasticScatt "<<theParameters->GetProbabilityOfElasticScatt()<<G4endl;
1206//G4cout<<"ProbabilityOfAnnihilation "<<theParameters->GetProbabilityOfAnnihilation()<<G4endl;
1207//G4cout<<"projectile->GetStatus target->GetStatus "<<projectile->GetStatus()<<" "<<target->GetStatus()<<G4endl;
1208//if((projectile->GetStatus() == 1) && (target->GetStatus() ==1))
1209//
1210if(collision.GetStatus()) // Uzhi Feb26
1211{
1212
1213//theParameters->SetProbabilityOfElasticScatt(1.);
1214//G4cout<<"before pro "<<projectile->Get4Momentum()<<" "<<projectile->Get4Momentum().mag()<<G4endl;
1215//G4cout<<"before tar "<<target->Get4Momentum()<<" "<<target->Get4Momentum().mag()<<G4endl;
1216//G4cout<<"Prob el "<<theParameters->GetProbabilityOfElasticScatt()<<G4endl;
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;
1269Successfull = 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 }
1298
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;
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;
1325
1326// target->SetStatus(2);
1327 }
1328 }
1329//
1330} // End of if((projectile->GetStatus() == 1) && (target->GetStatus() ==1))
1331
1332 } // end of while (theParticipants.Next())
1333
1334//Successfull=true;
1335//G4cout<<"G4FTFModel::ExciteParticipants() Successfull "<<Successfull<<G4endl;
1336//G4int Uzhi; G4cin>>Uzhi;
1337 return Successfull;
1338}
1339
1340//-------------------------------------------------------------------
1341void G4FTFModel::AjustTargetNucleonForAnnihilation(G4VSplitableHadron *SelectedAntiBaryon,
1342 G4VSplitableHadron *SelectedTargetNucleon)
1343{
1344 G4LorentzVector Pparticipants=SelectedAntiBaryon->Get4Momentum()+
1345 SelectedTargetNucleon->Get4Momentum();
1346
1347 G4V3DNucleus *theNucleus = GetWoundedNucleus();
1348//G4cout<<"Init A mass "<<theNucleus->GetMass()<<" "<<theNucleus->GetMassNumber()<<" "<<theNucleus->GetCharge()<<G4endl;
1349
1350 ResidualExcitationEnergy=0.;
1351 G4int ResidualCharge =theNucleus->GetCharge();
1352 G4int ResidualMassNumber=theNucleus->GetMassNumber();
1353
1354 G4ThreeVector P3nuclearResidual(0.,0.,0.);
1355 G4LorentzVector PnuclearResidual(0.,0.,0.,0.);
1356
1357
1358 G4double ExcitationEnergyPerWoundedNucleon=
1359 theParameters->GetExcitationEnergyPerWoundedNucleon();
1360//-------
1361 G4Nucleon * aNucleon;
1362 theNucleus->StartLoop();
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;}
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 {
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);
1483 }
1484
1485// Re-definition of the wounded nucleon momenta
1486 theNucleus->StartLoop();
1487 while ((aNucleon = theNucleus->GetNextNucleon()))
1488 {
1489 G4int CurrentStatus=0;
1490 if(aNucleon->AreYouHit())
1491 {
1492 if(aNucleon->GetSplitableHadron() == SelectedTargetNucleon)
1493 {CurrentStatus=1;}
1494 else
1495 {
1496 if(aNucleon->GetSplitableHadron()->GetSoftCollisionCount() == 0)
1497 {CurrentStatus=0;}
1498 else {CurrentStatus=1;}
1499 }
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())
1513
1514//-------------------------------
1515 return;
1516}
1517
1518// ------------------------------------------------------------
1519G4ExcitedStringVector * 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.
1531
1532 theParticipants.StartLoop(); // restart a loop
1533//
1534 while ( theParticipants.Next() )
1535 {
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 }
1545
1546//G4cout<<G4endl<<"primaries.size() "<<primaries.size()<<G4endl;
1547 for (unsigned int ahadron=0; ahadron < primaries.size() ; ahadron++)
1548 {
1549 G4bool isProjectile(0);
1550
1551 if(primaries[ahadron]->GetStatus() <= 1) {isProjectile=true; } // VU 10.04.2012
1552// if(primaries[ahadron]->GetStatus() == 3) {isProjectile=false;}
1553
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;
1562
1563//G4cout<<FirstString<<" "<<SecondString<<G4endl;
1564 }
1565
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);
1592
1593 delete ProjectileSplitable;
1594 }
1595 } // End of while ( (ProjectileNucleon=ProjectileNucleus->GetNextNucleon()) )
1596 } // End of if(ProjectileNucleus)
1597
1598//G4cout<<G4endl<<"theAdditionalString.size() "<<theAdditionalString.size()<<G4endl;
1599 if(theAdditionalString.size() != 0)
1600 {
1601 for (unsigned int ahadron=0; ahadron < theAdditionalString.size() ; ahadron++)
1602 {
1603 G4bool isProjectile(0);
1604
1605 if(theAdditionalString[ahadron]->GetStatus() <= 1) {isProjectile=true; } // VU 10.04.2012
1606// if(theAdditionalString[ahadron]->GetStatus() == 3) {isProjectile=false;}
1607
1608 FirstString=0; SecondString=0;
1609 theExcitation->CreateStrings(theAdditionalString[ahadron], isProjectile,
1610 FirstString, SecondString,
1611 theParameters);
1612
1613 if(FirstString != 0) strings->push_back(FirstString);
1614 if(SecondString != 0) strings->push_back(SecondString);
1615//G4cout<<"Quarks in the string in FTF"<<FirstString->GetRightParton()->GetPDGcode()<<" "<<FirstString->GetLeftParton()->GetPDGcode()<<G4endl;
1616//G4cout<<FirstString<<" "<<SecondString<<G4endl;
1617 }
1618 }
1619//G4cout<<"Check "<<strings->operator[](0)->GetRightParton()->GetPDGcode()<<" "<<strings->operator[](0)->GetLeftParton()->GetPDGcode()<<G4endl;
1620//G4cout<<"Check "<<strings->operator[](1)->GetRightParton()->GetPDGcode()<<" "<<strings->operator[](1)->GetLeftParton()->GetPDGcode()<<G4endl;
1621//
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 {}
1676
1677 }
1678/*
1679G4cout<<"Check "<<strings->operator[](0)->GetRightParton()->GetPDGcode()<<" "<<strings->operator[](0)->GetLeftParton()->GetPDGcode()<<G4endl;
1680G4cout<<"Check "<<strings->operator[](1)->GetRightParton()->GetPDGcode()<<" "<<strings->operator[](1)->GetLeftParton()->GetPDGcode()<<G4endl;
1681//G4cout<<"Check "<<strings->operator[](2)->GetRightParton()->GetPDGcode()<<" "<<strings->operator[](2)->GetLeftParton()->GetPDGcode()<<G4endl;
1682
1683G4cout<<"*** "<<strings->operator[](0)->GetRightParton()<<" "<<strings->operator[](0)->GetLeftParton()<<G4endl;
1684G4cout<<"*** "<<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/*
1690G4cout<<"*** "<<strings->operator[](0)->GetRightParton()<<" "<<strings->operator[](0)->GetLeftParton()<<G4endl;
1691G4cout<<"*** "<<strings->operator[](1)->GetRightParton()<<" "<<strings->operator[](1)->GetLeftParton()<<G4endl;
1692G4cout<<"*** "<<strings->operator[](2)->GetRightParton()<<" "<<strings->operator[](2)->GetLeftParton()<<G4endl;
1693
1694G4cout<<"Check "<<strings->operator[](0)->GetRightParton()->GetPDGcode()<<" "<<strings->operator[](0)->GetLeftParton()->GetPDGcode()<<G4endl;
1695G4cout<<"Check "<<strings->operator[](1)->GetRightParton()->GetPDGcode()<<" "<<strings->operator[](1)->GetLeftParton()->GetPDGcode()<<G4endl;
1696G4cout<<"Check "<<strings->operator[](2)->GetRightParton()->GetPDGcode()<<" "<<strings->operator[](2)->GetLeftParton()->GetPDGcode()<<G4endl;
1697*/
1698
1699/*
1700for (unsigned int ahadron=0; ahadron < strings->size() ; ahadron++)
1701{
1702G4cout<<ahadron<<" "<<strings->operator[](ahadron)->GetRightParton()->GetPDGcode()<<" "<<strings->operator[](ahadron)->GetLeftParton()->GetPDGcode()<<G4endl;
1703}
1704G4cout<<"------------------------"<<G4endl;
1705*/
1706//G4int Uzhi; G4cin >> Uzhi;
1707 return strings;
1708}
1709// ------------------------------------------------------------
1710void 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++ )
1725
1726}
1727
1728// ------------------------------------------------------------
1729G4ThreeVector G4FTFModel::GaussianPt(G4double AveragePt2, G4double maxPtSquare) const
1730{ // @@ this method is used in FTFModel as well. Should go somewhere common!
1731
1732 G4double Pt2(0.);
1733 if(AveragePt2 <= 0.) {Pt2=0.;}
1734 else
1735 {
1736 Pt2 = -AveragePt2 * std::log(1. + G4UniformRand() *
1737 (std::exp(-maxPtSquare/AveragePt2)-1.));
1738 }
1739 G4double Pt=std::sqrt(Pt2);
1740 G4double phi=G4UniformRand() * twopi;
1741
1742 return G4ThreeVector (Pt*std::cos(phi), Pt*std::sin(phi), 0.);
1743}
1744
1745void G4FTFModel::ModelDescription(std::ostream& desc) const
1746{
1747 desc << "please add description here" << G4endl;
1748}
std::vector< G4ExcitedString * > G4ExcitedStringVector
CLHEP::Hep3Vector G4ThreeVector
double G4double
Definition: G4Types.hh:64
int G4int
Definition: G4Types.hh:66
bool G4bool
Definition: G4Types.hh:67
#define G4endl
Definition: G4ios.hh:52
#define G4UniformRand()
Definition: Randomize.hh:53
double z() const
double x() const
double y() const
double theta() const
Hep3Vector boostVector() const
Hep3Vector vect() const
HepLorentzVector & rotateZ(double)
double perp2() const
HepLorentzVector & rotateY(double)
HepLorentzVector & transform(const HepRotation &)
static G4AntiNeutron * AntiNeutron()
static G4AntiProton * AntiProton()
Definition: G4AntiProton.cc:93
virtual G4bool ExciteParticipants(G4VSplitableHadron *aPartner, G4VSplitableHadron *bPartner, G4FTFParameters *theParameters, G4ElasticHNScattering *theElastic) const
virtual void CreateStrings(G4VSplitableHadron *aHadron, G4bool isProjectile, G4ExcitedString *&FirstString, G4ExcitedString *&SecondString, G4FTFParameters *theParameters) const
virtual G4bool ElasticScattering(G4VSplitableHadron *aPartner, G4VSplitableHadron *bPartner, G4FTFParameters *theParameters) const
virtual G4bool Annihilate(G4VSplitableHadron *aPartner, G4VSplitableHadron *bPartner, G4VSplitableHadron *&AdditionalString, G4FTFParameters *theParameters) const
G4V3DNucleus * GetWoundedNucleus() const
Definition: G4FTFModel.hh:121
void Init(const G4Nucleus &aNucleus, const G4DynamicParticle &aProjectile)
Definition: G4FTFModel.cc:112
G4FTFModel(const G4String &modelName="FTF")
Definition: G4FTFModel.cc:56
virtual void ModelDescription(std::ostream &) const
Definition: G4FTFModel.cc:1745
G4ExcitedStringVector * GetStrings()
Definition: G4FTFModel.cc:206
G4double GetProbabilityOfAnnihilation()
G4double GetMaxNumberOfCollisions()
G4double GetPt2ofNuclearDestruction()
G4double GetMaxPt2ofNuclearDestruction()
G4double GetProbabilityOfElasticScatt()
G4double GetExcitationEnergyPerWoundedNucleon()
G4double GetProbOfInteraction()
G4double GetDofNuclearDestruction()
G4double GetR2ofNuclearDestruction()
G4double GetCofNuclearDestruction()
void SetProjectileNucleus(G4V3DNucleus *aNucleus)
G4V3DNucleus * GetProjectileNucleus()
void GetList(const G4ReactionProduct &thePrimary, G4FTFParameters *theParameters)
void InitProjectileNucleus(G4int theZ, G4int theA)
G4V3DNucleus * theProjectileNucleus
G4InteractionContent & GetInteraction()
G4Nucleon * GetTargetNucleon() const
G4VSplitableHadron * GetProjectile() const
void SetStatus(G4int aValue)
G4VSplitableHadron * GetTarget() const
G4double GetIonMass(G4int Z, G4int A, G4int L=0) const
!! Only ground states are supported now
Definition: G4IonTable.cc:774
G4bool AreYouHit() const
Definition: G4Nucleon.hh:97
void SetMomentum(G4LorentzVector &aMomentum)
Definition: G4Nucleon.hh:70
G4VSplitableHadron * GetSplitableHadron() const
Definition: G4Nucleon.hh:96
virtual const G4LorentzVector & Get4Momentum() const
Definition: G4Nucleon.hh:72
void Hit(G4VSplitableHadron *aHit)
Definition: G4Nucleon.hh:90
virtual const G4ThreeVector & GetPosition() const
Definition: G4Nucleon.hh:68
void SetParticleType(G4Proton *aProton)
Definition: G4Nucleon.hh:77
void SetBindingEnergy(G4double anEnergy)
Definition: G4Nucleon.hh:74
virtual G4ParticleDefinition * GetDefinition() const
Definition: G4Nucleon.hh:85
G4int GetA_asInt() const
Definition: G4Nucleus.hh:109
G4int GetZ_asInt() const
Definition: G4Nucleus.hh:115
G4double GetPDGCharge() const
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
static G4ParticleTable * GetParticleTable()
G4IonTable * GetIonTable()
G4double GetTotalMomentum() const
G4double GetTotalEnergy() const
G4ThreeVector GetMomentum() const
G4ParticleDefinition * GetDefinition() const
virtual G4Nucleon * GetNextNucleon()=0
virtual G4int GetCharge()=0
virtual G4bool StartLoop()=0
virtual void DoLorentzBoost(const G4LorentzVector &theBoost)=0
virtual G4int GetMassNumber()=0
void SetEnergyMomentumCheckLevels(G4double relativeLevel, G4double AbsoluteLevel)
virtual void Init(G4int theZ, G4int theA)
G4V3DNucleus * theNucleus
virtual G4V3DNucleus * GetWoundedNucleus() const
void SetThisPointer(G4VPartonStringModel *aPointer)
void SetTimeOfCreation(G4double aTime)
void SetStatus(const G4int aStatus)
G4ParticleDefinition * GetDefinition() const
void Set4Momentum(const G4LorentzVector &a4Momentum)
const G4LorentzVector & Get4Momentum() const
void SetDefinition(G4ParticleDefinition *aDefinition)
void operator()(G4VSplitableHadron *aH)
Definition: G4FTFModel.cc:68
T sqr(const T &x)
Definition: templates.hh:145