Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
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