Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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"
43 #include "G4PhysicalConstants.hh"
44 #include "G4SystemOfUnits.hh"
45 #include "G4FTFParameters.hh"
46 #include "G4FTFParticipants.hh"
47 #include "G4DiffractiveSplitableHadron.hh"
48 #include "G4InteractionContent.hh"
49 #include "G4LorentzRotation.hh"
50 #include "G4ParticleDefinition.hh"
51 #include "G4ParticleTable.hh"
52 #include "G4IonTable.hh"
53 
54 // Class G4FTFModel
55 
56 G4FTFModel::G4FTFModel(const G4String& modelName):G4VPartonStringModel(modelName),
57  theExcitation(new G4DiffractiveExcitation()),
58  theElastic(new G4ElasticHNScattering()),
59  theAnnihilation(new G4FTFAnnihilation())
60 {
61  G4VPartonStringModel::SetThisPointer(this);
62  theParameters=0;
63  NumberOfInvolvedNucleon=0;
64  NumberOfInvolvedNucleonOfProjectile=0;
65  SetEnergyMomentumCheckLevels(2*perCent, 150*MeV);
66 }
67 
68 struct DeleteVSplitableHadron { void operator()(G4VSplitableHadron * aH){ delete aH;} };
69 
70 G4FTFModel::~G4FTFModel()
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(),
84  DeleteVSplitableHadron());
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 // ------------------------------------------------------------
112 void 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 /*
119 G4cout<<"FTF init Pro Name "<<theProjectile.GetDefinition()->GetParticleName()<<G4endl;
120 G4cout<<"FTF init Pro Mass "<<theProjectile.GetMass()<<" "<<theProjectile.GetMomentum()<<G4endl;
121 G4cout<<"FTF init Pro B Q "<<theProjectile.GetDefinition()->GetBaryonNumber()<<" "<<(G4int) theProjectile.GetDefinition()->GetPDGCharge()<<G4endl;
122 G4cout<<"FTF init A Z "<<aNucleus.GetA_asInt()<<" "<<aNucleus.GetZ_asInt()<<G4endl;
123 G4cout<<" "<<aNucleus.GetN()<<" "<<aNucleus.GetZ()<<G4endl;
124 //G4int Uzhi; G4cin>>Uzhi;
125 */
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
168  {aNucleon->SetParticleType(G4AntiProton::AntiProton());}
169 
170  if(aNucleon->GetDefinition()->GetPDGEncoding() == 2112) // Neutron
171  {aNucleon->SetParticleType(G4AntiNeutron::AntiNeutron());}
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(),
199  DeleteVSplitableHadron());
200  }
201  theAdditionalString.clear();
202 //G4cout<<" End Init theProjectile.GetMomentum()"<<theProjectile.GetMomentum()<<G4endl;
203 }
204 
205 // ------------------------------------------------------------
206 G4ExcitedStringVector * G4FTFModel::GetStrings()
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 //-------------------------------------------------------------------
343 void G4FTFModel::StoreInvolvedNucleon()
344 { //--- To store nucleons involved in low energy interaction -------
345 if(std::abs(theProjectile.GetDefinition()->GetBaryonNumber()) <= 1)
346 { // the projectile is a hadron -----------
347 //G4cout<<"the projectile is a hadron"<<G4endl;
348  NumberOfInvolvedNucleon=0;
349 
350  theParticipants.StartLoop();
351 
352  while (theParticipants.Next())
353  {
354 //G4cout<<"theParticipants.Next()"<<G4endl;
355  const G4InteractionContent & collision=theParticipants.GetInteraction();
356 //G4cout<<"collision=theParticipants.GetInteraction()"<<G4endl;
357  G4Nucleon * TargetNucleon=collision.GetTargetNucleon();
358 //G4cout<<"TargetNucleon=collision.GetTargetNucleon()"<<G4endl;
359 
360  TheInvolvedNucleon[NumberOfInvolvedNucleon]=TargetNucleon;
361  NumberOfInvolvedNucleon++;
362 //G4cout<<G4endl<<"Prim NumberOfInvolvedNucleon "<<NumberOfInvolvedNucleon<<G4endl;
363  } // end of while (theParticipants.Next())
364 
365  NumberOfInvolvedTargetNucleon=NumberOfInvolvedNucleon;
366 // ---------------- Calculation of creation time for each target nucleon -----------
367 //G4cout<<"theParticipants.StartLoop() "<<G4endl;
368  theParticipants.StartLoop(); // restart a loop
369 //G4cout<<"theParticipants.Next();"<<G4endl;
370  theParticipants.Next();
371  G4VSplitableHadron * primary = theParticipants.GetInteraction().GetProjectile();
372 //G4cout<<"primary->Get4Momentum() "<<primary->Get4Momentum()<<G4endl;
373 //G4cout<<"primary->Get4Momentum().pz() "<<primary->Get4Momentum().pz()<<G4endl;
374 //G4cout<<"primary->Get4Momentum().e() "<<primary->Get4Momentum().e()<<G4endl;
375 
376  G4double betta_z=primary->Get4Momentum().pz()/primary->Get4Momentum().e();
377 //G4cout<<"betta_z "<<betta_z<<G4endl;
378  primary->SetTimeOfCreation(0.);
379 
380  G4double ZcoordinateOfPreviousCollision(0.);
381  G4double ZcoordinateOfCurrentInteraction(0.);
382  G4double TimeOfPreviousCollision(0.);
383  G4double TimeOfCurrentCollision(0);
384 
385  theParticipants.theNucleus->StartLoop();
386  G4Nucleon * aNucleon;
387  G4bool theFirstInvolvedNucleon(true);
388  while ( (aNucleon = theParticipants.theNucleus->GetNextNucleon()) )
389  {
390 //G4cout<<"aNucleon->AreYouHit() "<<aNucleon->AreYouHit()<<G4endl;
391  if(aNucleon->AreYouHit())
392  {
393  if(theFirstInvolvedNucleon)
394  {
395  ZcoordinateOfPreviousCollision=aNucleon->GetPosition().z();
396 //G4cout<<"ZcoordinateOfPreviousCollision "<<ZcoordinateOfPreviousCollision/fermi<<G4endl;
397  theFirstInvolvedNucleon=false;
398  }
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 //-------------------------------------------------------------------
465 void G4FTFModel::ReggeonCascade()
466 { //--- Implementation of the reggeon theory inspired model-------
467 //G4cout<<"In reggeon"<<G4endl;
468 
469  if(std::abs(theProjectile.GetDefinition()->GetBaryonNumber()) > 1) return;
470 // For Nucleus-nucleus or Anti-nucleus - nucleus interactions
471 // the cascading will be implemented latter.
472 
473  NumberOfInvolvedNucleon=0;
474 
475  theParticipants.StartLoop();
476  while (theParticipants.Next())
477  {
478  const G4InteractionContent & collision=theParticipants.GetInteraction();
479 
480  G4Nucleon * TargetNucleon=collision.GetTargetNucleon();
481 
482  TheInvolvedNucleon[NumberOfInvolvedNucleon]=TargetNucleon;
483  NumberOfInvolvedNucleon++;
484 //G4cout<<"Prim NumberOfInvolvedNucleon "<<NumberOfInvolvedNucleon<<G4endl;
485  G4double XofWoundedNucleon = TargetNucleon->GetPosition().x();
486  G4double YofWoundedNucleon = TargetNucleon->GetPosition().y();
487 
488  theParticipants.theNucleus->StartLoop();
489  G4Nucleon * Neighbour(0);
490 
491  while ( (Neighbour = theParticipants.theNucleus->GetNextNucleon()) )
492  {
493  if(!Neighbour->AreYouHit())
494  {
495  G4double impact2= sqr(XofWoundedNucleon - Neighbour->GetPosition().x()) +
496  sqr(YofWoundedNucleon - Neighbour->GetPosition().y());
497 
498  if(G4UniformRand() < theParameters->GetCofNuclearDestruction()*
499  std::exp(-impact2/theParameters->GetR2ofNuclearDestruction()))
500  { // The neighbour nucleon is involved in the reggeon cascade
501 
502  TheInvolvedNucleon[NumberOfInvolvedNucleon]=Neighbour;
503  NumberOfInvolvedNucleon++;
504 //G4cout<<"Seco NumberOfInvolvedNucleon "<<NumberOfInvolvedNucleon<<G4endl;
505 
506  G4VSplitableHadron *targetSplitable;
507  targetSplitable = new G4DiffractiveSplitableHadron(*Neighbour);
508 
509  Neighbour->Hit(targetSplitable);
510  targetSplitable->SetStatus(2);
511  }
512  } // end of if(!Neighbour->AreYouHit())
513  } // end of while (theParticipant.theNucleus->GetNextNucleon())
514  } // end of while (theParticipants.Next())
515 
516  NumberOfInvolvedTargetNucleon=NumberOfInvolvedNucleon;
517 
518 // ---------------- Calculation of creation time for each target nucleon -----------
519  theParticipants.StartLoop(); // restart a loop
520  theParticipants.Next();
521  G4VSplitableHadron * primary = theParticipants.GetInteraction().GetProjectile();
522  G4double betta_z=primary->Get4Momentum().pz()/primary->Get4Momentum().e();
523  primary->SetTimeOfCreation(0.);
524 
525  G4double ZcoordinateOfPreviousCollision(0.);
526  G4double ZcoordinateOfCurrentInteraction(0.);
527  G4double TimeOfPreviousCollision(0.);
528  G4double TimeOfCurrentCollision(0);
529 
530  theParticipants.theNucleus->StartLoop();
531  G4Nucleon * aNucleon;
532  G4bool theFirstInvolvedNucleon(true);
533  while ( (aNucleon = theParticipants.theNucleus->GetNextNucleon()) )
534  {
535  if(aNucleon->AreYouHit())
536  {
537  if(theFirstInvolvedNucleon)
538  {
539  ZcoordinateOfPreviousCollision=aNucleon->GetPosition().z();
540  theFirstInvolvedNucleon=false;
541  }
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 // ------------------------------------------------------------
559 G4bool 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=
591  theParameters->GetExcitationEnergyPerWoundedNucleon();
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  {
622  ResidualMass=G4ParticleTable::GetParticleTable()->GetIonTable()->
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=
761  theParameters->GetExcitationEnergyPerWoundedNucleon();
762 
763 //G4cout<<"ExcitationEnergyPerWoundedNucleon "<<ExcitationEnergyPerWoundedNucleon<<G4endl;
764 
765  theNucleus->StartLoop();
766 
767  while (G4Nucleon * aNucleon = theNucleus->GetNextNucleon())
768  {
769  Ptarget+=aNucleon->Get4Momentum();
770 
771  if(aNucleon->AreYouHit())
772  { // Involved nucleons
773 //G4cout<<"PutOn Tr "<<aNucleon->Get4Momentum()<<G4endl;
774 // Psum += aNucleon->Get4Momentum(); // Uzhi 20 Sept.
775 // if(!ProjectileIsAntiBaryon) // Uzhi 13.06.2012
776 // {
777  SumMasses += std::sqrt(sqr(aNucleon->GetDefinition()->GetPDGMass()) //Uzhi 12.06.2012
778  + aNucleon->Get4Momentum().perp2()); //Uzhi 12.06.2012
779  SumMasses += 20.*MeV; // 13.12.09 Separation energy for a nucleon
780  ResidualExcitationEnergy+=ExcitationEnergyPerWoundedNucleon;
781 /* //Uzhi 13.06.2012
782  } else
783  {
784  SumMasses += aNucleon->Get4Momentum().mag(); // 4.12.2010
785  G4LorentzVector tmp=aNucleon->Get4Momentum();
786  tmp.setE(aNucleon->Get4Momentum().mag()); // It is need to save mass 6.12.2011
787  aNucleon->SetMomentum(tmp);
788  }
789 */ //Uzhi 13.06.2012
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  {
812  ResidualMass=G4ParticleTable::GetParticleTable()->GetIonTable()->
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.perp2()); // Uzhi 16.05.2013
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.perp2()); // Uzhi 16.05.2013
834  SumMasses += std::sqrt(sqr(ResidualMass+ResidualExcitationEnergy)
835  +PnuclearResidual.perp2()); // Uzhi 16.05.2013
836  if(SqrtS < SumMasses) // Uzhi 12.06.2012
837  {
838  SumMasses -= std::sqrt(sqr(ResidualMass+ResidualExcitationEnergy)
839  +PnuclearResidual.perp2()); // Uzhi 16.05.2013
840  SumMasses += std::sqrt(sqr(ResidualMass)+PnuclearResidual.perp2());// Uzhi 16.05.2013
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
869  G4ParticleDefinition* ptr =
870  G4ParticleTable::GetParticleTable()->FindParticle(newPDGcode);
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.perp2())/XminusSum; // Uzhi 16.05.2013
1017  }
1018 //G4cout<<"InerSuccess "<<InerSuccess<<" "<<SqrtS<<" "<<Mprojectile + std::sqrt(M2target)<<" "<<std::sqrt(M2target)<<G4endl;
1019 //G4int Uzhi;G4cin>>Uzhi;
1020  } while(!InerSuccess);
1021  } while (SqrtS < Mprojectile + std::sqrt(M2target));
1022 //-------------------------------------------------------------
1023  G4double DecayMomentum2= S*S+M2projectile*M2projectile+M2target*M2target
1024  -2.*S*M2projectile - 2.*S*M2target
1025  -2.*M2projectile*M2target;
1026 
1027  WminusTarget=(S-M2projectile+M2target+std::sqrt(DecayMomentum2))/2./SqrtS;
1028  WplusProjectile=SqrtS - M2target/WminusTarget;
1029 
1030  G4double Pzprojectile=WplusProjectile/2. - M2projectile/2./WplusProjectile;// 8.06.11
1031  G4double Eprojectile =WplusProjectile/2. + M2projectile/2./WplusProjectile;// 8.06.11
1032  G4double Yprojectile=0.5*std::log((Eprojectile+Pzprojectile)/
1033  (Eprojectile-Pzprojectile)); // 1.07.11
1034 
1035 //G4cout<<"Yprojectile "<<Yprojectile<<G4endl;
1036 //G4LorentzVector TestPprojectile=Pprojectile;
1037 //TestPprojectile.setPz(Pzprojectile); TestPprojectile.setE(Eprojectile);
1038 
1039 //G4cout<<"DecayMomentum2 "<<DecayMomentum2<<G4endl;
1040 //G4cout<<"WminusTarget WplusProjectile "<<WminusTarget<<" "<<WplusProjectile<<G4endl;
1041 //G4int Uzhi;G4cin>>Uzhi;
1042 //-------------------------------------------------------------
1043  for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
1044  {
1045  G4Nucleon * aNucleon = TheInvolvedNucleon[i];
1046  G4LorentzVector tmp=aNucleon->Get4Momentum();
1047 
1048  G4double Mt2(0.);
1049 
1050 // if(!ProjectileIsAntiBaryon) // 4.12.2010
1051 // {
1052  Mt2 = sqr(tmp.x())+sqr(tmp.y())+
1053  sqr(aNucleon->GetSplitableHadron()->GetDefinition()->GetPDGMass());
1054 /*
1055  } else
1056  {
1057  Mt2 = sqr(tmp.x())+sqr(tmp.y())+ // 4.12.2010
1058  sqr(aNucleon->Get4Momentum().e()); // sqr
1059  }
1060 */
1061  G4double Xminus=tmp.z();
1062 
1063  G4double Pz=-WminusTarget*Xminus/2. + Mt2/(2.*WminusTarget*Xminus);
1064  G4double E = WminusTarget*Xminus/2. + Mt2/(2.*WminusTarget*Xminus);
1065  G4double YtargetNucleon=0.5*std::log((E+Pz)/(E-Pz)); // 1.07.11 //Uzhi 20 Sept.
1066 
1067 //G4cout<<"YtargetNucleon "<<YtargetNucleon<<G4endl;
1068 //G4cout<<"YtargetNucleon-YtargetNucleus "<<YtargetNucleon-YtargetNucleus<<G4endl;
1069 //G4cout<<"Yprojectile YtargetNucleon "<<Yprojectile<<" "<<YtargetNucleon<<G4endl;
1070 if((std::abs(YtargetNucleon-YtargetNucleus) > 2) ||
1071  (Yprojectile < YtargetNucleon)) {OuterSuccess=false; break;} // 1.07.11
1072 
1073  } // end of for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
1074 //if(ProjectileIsAntiBaryon) {G4int Uzhi;G4cin>>Uzhi;}
1075  } while(!OuterSuccess);
1076 
1077 //-------------------------------------------------------------
1078  G4double Pzprojectile=WplusProjectile/2. - M2projectile/2./WplusProjectile;
1079  G4double Eprojectile =WplusProjectile/2. + M2projectile/2./WplusProjectile;
1080  Pprojectile.setPz(Pzprojectile); Pprojectile.setE(Eprojectile);
1081 //G4cout<<"Proj after in CMS "<<Pprojectile<<G4endl;
1082 
1083  Pprojectile.transform(toLab); // The work with the projectile
1084  primary->Set4Momentum(Pprojectile); // is finished at the moment.
1085 //G4cout<<"Final proj mom "<<primary->Get4Momentum()<<G4endl;
1086 
1087 //-------------------------------------------------------------
1088  G4ThreeVector Residual3Momentum(0.,0.,1.);
1089 
1090  for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
1091  {
1092  G4Nucleon * aNucleon = TheInvolvedNucleon[i];
1093  G4LorentzVector tmp=aNucleon->Get4Momentum();
1094 //G4cout<<"trg "<<aNucleon->Get4Momentum()<<G4endl;
1095  Residual3Momentum-=tmp.vect();
1096 
1097  G4double Mt2(0.);
1098 
1099 // if(!ProjectileIsAntiBaryon) // 4.12.2010
1100 // {
1101  Mt2 = sqr(tmp.x())+sqr(tmp.y())+
1102  sqr(aNucleon->GetSplitableHadron()->GetDefinition()->GetPDGMass());
1103 /*
1104  } else
1105  {
1106  Mt2 = sqr(tmp.x())+sqr(tmp.y())+ // 4.12.2010
1107  aNucleon->Get4Momentum().e()*aNucleon->Get4Momentum().e();
1108  }
1109 */
1110  G4double Xminus=tmp.z();
1111 
1112  G4double Pz=-WminusTarget*Xminus/2. + Mt2/(2.*WminusTarget*Xminus);
1113  G4double E = WminusTarget*Xminus/2. + Mt2/(2.*WminusTarget*Xminus);
1114 
1115  tmp.setPz(Pz);
1116  tmp.setE(E);
1117 //G4cout<<"Targ after in CMS "<<tmp<<G4endl;
1118  tmp.transform(toLab);
1119 
1120  aNucleon->SetMomentum(tmp);
1121 //G4cout<<"Targ after in LAB "<<aNucleon->Get4Momentum()<<G4endl;
1122  G4VSplitableHadron * targetSplitable=aNucleon->GetSplitableHadron();
1123  targetSplitable->Set4Momentum(tmp);
1124 
1125  } // end of for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
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 // ------------------------------------------------------------
1156 G4bool 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 //
1210 if(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;
1269 Successfull = false; break; // 1.07.11
1270 
1271  if(NumberOfInvolvedTargetNucleon > 1)
1272  {
1273  NumberOfInvolvedTargetNucleon--;
1274  target->SetStatus(4); // 1->0 return nucleon to the target VU 10.04.2012
1275  }
1276  }
1277 */
1278  } // end if NumOfInel
1279  }
1280  else // Annihilation
1281  {
1282 //G4cout<<"Annihilation"<<G4endl;
1283 //G4cout<<"Before pro "<<projectile->Get4Momentum()<<G4endl;
1284 //G4cout<<"Before tar "<<target->Get4Momentum()<<G4endl;
1285 //G4cout<<"Mom pro "<<theProjectile.GetTotalMomentum()<<G4endl;
1286 //if(theProjectile.GetTotalMomentum() < 2000.*MeV)
1287 {
1288  while (theParticipants.Next())
1289  {
1290  G4InteractionContent & acollision=theParticipants.GetInteraction();//Uzhi Feb26
1291 
1292  G4VSplitableHadron * NextProjectileNucleon=acollision.GetProjectile(); // Uzhi Feb26
1293  G4VSplitableHadron * NextTargetNucleon =acollision.GetTarget();
1294  if((projectile == NextProjectileNucleon) ||
1295  (target == NextTargetNucleon )) acollision.SetStatus(0);
1296 // if(target != NextTargetNucleon) NextTargetNucleon->SetStatus(0); // Uzhi Feb26
1297  }
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 //-------------------------------------------------------------------
1341 void 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  {
1415  ResidualMass=G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(
1416  ResidualCharge,ResidualMassNumber);
1417  } else {return;}
1418 
1419 //G4cout<<"Res Mass E* "<<ResidualMass<<" "<<ResidualExcitationEnergy<<G4endl;
1420 
1421  if(ResidualMass > SqrtS) {return;}
1422  else
1423  {
1424  if(ResidualMass+ResidualExcitationEnergy > SqrtS)
1425  ResidualExcitationEnergy = SqrtS-ResidualMass;
1426  }
1427 
1428  ResidualMass+=ResidualExcitationEnergy;
1429  G4double ResidualMass2=sqr(ResidualMass);
1430 //G4cout<<"New Res Mass E* "<<ResidualMass<<" "<<ResidualExcitationEnergy<<G4endl;
1431 
1432 //-------
1433  G4double ParticipantMass=Pparticipants.mag();
1434  G4double ParticipantMass2=sqr(ParticipantMass);
1435 
1436  if(ResidualMass + ParticipantMass > SqrtS) ParticipantMass=SqrtS-ResidualMass;
1437 
1438 //G4cout<<"Parts P "<<Pparticipants<<G4endl;
1439 //G4cout<<"Res Nuc "<<PnuclearResidual<<G4endl;
1440 
1441  G4double DecayMomentum2=
1442  sqr(S)+sqr(ParticipantMass2)+ sqr(ResidualMass2) -
1443  2.*S*ParticipantMass2 - 2.*S*ResidualMass2
1444  -2.*ParticipantMass2*ResidualMass2;
1445 
1446  if(DecayMomentum2 < 0.) return;
1447 
1448  DecayMomentum2/=(4.*S);
1449  G4double DecayMomentum = std::sqrt(DecayMomentum2);
1450 //G4cout<<"DecayMomentum "<<DecayMomentum<<G4endl;
1451 
1452  G4LorentzVector New_Pparticipants(0.,0.,DecayMomentum,
1453  std::sqrt(DecayMomentum2+ParticipantMass2));
1454 
1455  G4LorentzVector New_PnuclearResidual(0.,0.,-DecayMomentum,
1456  std::sqrt(DecayMomentum2+ResidualMass2));
1457 
1458 //G4cout<<"New part P "<<New_Pparticipants<<" "<<New_Pparticipants.mag()<<G4endl;
1459 //G4cout<<"New resd P "<<New_PnuclearResidual<<" "<<New_PnuclearResidual.mag()<<G4endl;
1460 
1461  New_Pparticipants.transform(toLab);
1462  New_PnuclearResidual.transform(toLab);
1463 //G4cout<<"New part P "<<New_Pparticipants<<" "<<New_Pparticipants.mag()<<G4endl;
1464 //G4cout<<"New resd P "<<New_PnuclearResidual<<" "<<New_PnuclearResidual.mag()<<G4endl;
1465 
1466  G4LorentzVector DeltaP_participants=(Pparticipants - New_Pparticipants)/2.;
1467  G4LorentzVector DeltaP_nuclearResidual=(PnuclearResidual - New_PnuclearResidual)/
1468  (G4double) ResidualMassNumber;
1469 //------------------
1470 
1471  Ptmp=SelectedAntiBaryon->Get4Momentum() - DeltaP_participants;
1472  SelectedAntiBaryon->Set4Momentum(Ptmp);
1473 
1474  Ptmp=SelectedTargetNucleon->Get4Momentum() - DeltaP_participants;
1475  SelectedTargetNucleon->Set4Momentum(Ptmp);
1476 //-----------
1477 
1478  //A.R. 25-Jul-2012 : Fix to Covery "Division by zero"
1479  //G4double DeltaExcitationEnergy=ResidualExcitationEnergy/((G4double) NumberOfHoles);
1480  G4double DeltaExcitationEnergy = 0.0;
1481  if ( NumberOfHoles != 0 ) {
1482  DeltaExcitationEnergy = ResidualExcitationEnergy / ((G4double) NumberOfHoles);
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 // ------------------------------------------------------------
1519 G4ExcitedStringVector * G4FTFModel::BuildStrings()
1520 {
1521 // Loop over all collisions; find all primaries, and all target ( targets may
1522 // be duplicate in the List ( to unique G4VSplitableHadrons)
1523 
1524  G4ExcitedStringVector * strings;
1525  strings = new G4ExcitedStringVector();
1526 
1527  std::vector<G4VSplitableHadron *> primaries;
1528 
1529  G4ExcitedString * FirstString(0); // If there will be a kink,
1530  G4ExcitedString * SecondString(0); // two strings will be produced.
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 /*
1679 G4cout<<"Check "<<strings->operator[](0)->GetRightParton()->GetPDGcode()<<" "<<strings->operator[](0)->GetLeftParton()->GetPDGcode()<<G4endl;
1680 G4cout<<"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 
1683 G4cout<<"*** "<<strings->operator[](0)->GetRightParton()<<" "<<strings->operator[](0)->GetLeftParton()<<G4endl;
1684 G4cout<<"*** "<<strings->operator[](1)->GetRightParton()<<" "<<strings->operator[](1)->GetLeftParton()<<G4endl;
1685 //G4cout<<"*** "<<strings->operator[](2)->GetRightParton()<<" "<<strings->operator[](2)->GetLeftParton()<<G4endl;
1686 */
1687  std::for_each(primaries.begin(), primaries.end(), DeleteVSplitableHadron());
1688  primaries.clear();
1689 /*
1690 G4cout<<"*** "<<strings->operator[](0)->GetRightParton()<<" "<<strings->operator[](0)->GetLeftParton()<<G4endl;
1691 G4cout<<"*** "<<strings->operator[](1)->GetRightParton()<<" "<<strings->operator[](1)->GetLeftParton()<<G4endl;
1692 G4cout<<"*** "<<strings->operator[](2)->GetRightParton()<<" "<<strings->operator[](2)->GetLeftParton()<<G4endl;
1693 
1694 G4cout<<"Check "<<strings->operator[](0)->GetRightParton()->GetPDGcode()<<" "<<strings->operator[](0)->GetLeftParton()->GetPDGcode()<<G4endl;
1695 G4cout<<"Check "<<strings->operator[](1)->GetRightParton()->GetPDGcode()<<" "<<strings->operator[](1)->GetLeftParton()->GetPDGcode()<<G4endl;
1696 G4cout<<"Check "<<strings->operator[](2)->GetRightParton()->GetPDGcode()<<" "<<strings->operator[](2)->GetLeftParton()->GetPDGcode()<<G4endl;
1697 */
1698 
1699 /*
1700 for (unsigned int ahadron=0; ahadron < strings->size() ; ahadron++)
1701 {
1702 G4cout<<ahadron<<" "<<strings->operator[](ahadron)->GetRightParton()->GetPDGcode()<<" "<<strings->operator[](ahadron)->GetLeftParton()->GetPDGcode()<<G4endl;
1703 }
1704 G4cout<<"------------------------"<<G4endl;
1705 */
1706 //G4int Uzhi; G4cin >> Uzhi;
1707  return strings;
1708 }
1709 // ------------------------------------------------------------
1710 void G4FTFModel::GetResidualNucleus()
1711 { // This method is needed for the correct application of G4PrecompoundModelInterface
1712  G4double DeltaExcitationE=ResidualExcitationEnergy/
1713  (G4double) NumberOfInvolvedNucleon;
1714  G4LorentzVector DeltaPResidualNucleus = Residual4Momentum/
1715  (G4double) NumberOfInvolvedNucleon;
1716 
1717  for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
1718  {
1719  G4Nucleon * aNucleon = TheInvolvedNucleon[i];
1720 // G4LorentzVector tmp=aNucleon->Get4Momentum()-DeltaPResidualNucleus;
1721  G4LorentzVector tmp=-DeltaPResidualNucleus;
1722  aNucleon->SetMomentum(tmp);
1723  aNucleon->SetBindingEnergy(DeltaExcitationE);
1724  } // end of for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
1725 
1726 }
1727 
1728 // ------------------------------------------------------------
1729 G4ThreeVector 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 
1745 void G4FTFModel::ModelDescription(std::ostream& desc) const
1746 {
1747  desc << "please add description here" << G4endl;
1748 }
Generated on Sat May 25 2013 14:34:05 for Geant4 by   doxygen 1.8.4