Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/parton_string/diffraction/src/G4FTFModel.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /processes/hadronic/models/parton_string/diffraction/src/G4FTFModel.cc (Version 11.3.0) and /processes/hadronic/models/parton_string/diffraction/src/G4FTFModel.cc (Version 9.4)


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