Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/parton_string/diffraction/src/G4DiffractiveExcitation.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/G4DiffractiveExcitation.cc (Version 11.3.0) and /processes/hadronic/models/parton_string/diffraction/src/G4DiffractiveExcitation.cc (Version 9.3.p2)


  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 //                                             <<  27 // $Id: G4DiffractiveExcitation.cc,v 1.21 2009/12/15 19:14:31 vuzhinsk Exp $
 28                                                << 
 29 // -------------------------------------------     28 // ------------------------------------------------------------
 30 //      GEANT 4 class implemetation file           29 //      GEANT 4 class implemetation file
 31 //                                                 30 //
 32 //      ---------------- G4DiffractiveExcitati     31 //      ---------------- G4DiffractiveExcitation --------------
 33 //             by Gunter Folger, October 1998.     32 //             by Gunter Folger, October 1998.
 34 //        diffractive Excitation used by strin     33 //        diffractive Excitation used by strings models
 35 //             Take a projectile and a target      34 //             Take a projectile and a target
 36 //             excite the projectile and targe     35 //             excite the projectile and target
 37 //  Essential changed by V. Uzhinsky in Novemb     36 //  Essential changed by V. Uzhinsky in November - December 2006
 38 //  in order to put it in a correspondence wit     37 //  in order to put it in a correspondence with original FRITIOF
 39 //  model. Variant of FRITIOF with nucleon de-     38 //  model. Variant of FRITIOF with nucleon de-excitation is implemented.
 40 //  Other changes by V.Uzhinsky in May 2007 we     39 //  Other changes by V.Uzhinsky in May 2007 were introduced to fit
 41 //  meson-nucleon interactions. Additional cha     40 //  meson-nucleon interactions. Additional changes by V. Uzhinsky
 42 //  were introduced in December 2006. They tre     41 //  were introduced in December 2006. They treat diffraction dissociation
 43 //  processes more exactly.                        42 //  processes more exactly.
 44 //  Correct treatment of the diffraction disso << 
 45 //  Mass distributions for resonances and uu-d << 
 46 //  and dd-diquarks suppression in neutrons we << 
 47 // -------------------------------------------     43 // ---------------------------------------------------------------------
 48                                                    44 
                                                   >>  45 
 49 #include "globals.hh"                              46 #include "globals.hh"
 50 #include "Randomize.hh"                            47 #include "Randomize.hh"
 51 #include "G4PhysicalConstants.hh"              << 
 52 #include "G4SystemOfUnits.hh"                  << 
 53                                                    48 
 54 #include "G4DiffractiveExcitation.hh"              49 #include "G4DiffractiveExcitation.hh"
 55 #include "G4FTFParameters.hh"                      50 #include "G4FTFParameters.hh"
 56 #include "G4ElasticHNScattering.hh"                51 #include "G4ElasticHNScattering.hh"
 57                                                    52 
                                                   >>  53 #include "G4LorentzRotation.hh"
 58 #include "G4RotationMatrix.hh"                     54 #include "G4RotationMatrix.hh"
                                                   >>  55 #include "G4ThreeVector.hh"
 59 #include "G4ParticleDefinition.hh"                 56 #include "G4ParticleDefinition.hh" 
 60 #include "G4ParticleTable.hh"                  << 
 61 #include "G4SampleResonance.hh"                << 
 62 #include "G4VSplitableHadron.hh"                   57 #include "G4VSplitableHadron.hh"
 63 #include "G4ExcitedString.hh"                      58 #include "G4ExcitedString.hh"
                                                   >>  59 #include "G4ParticleTable.hh"
 64 #include "G4Neutron.hh"                            60 #include "G4Neutron.hh"
 65                                                <<  61 #include "G4ParticleDefinition.hh"
 66 #include "G4Exp.hh"                            << 
 67 #include "G4Log.hh"                            << 
 68 #include "G4Pow.hh"                            << 
 69                                                    62 
 70 //#include "G4ios.hh"                              63 //#include "G4ios.hh"
                                                   >>  64 //#include "UZHI_diffraction.hh"
 71                                                    65 
 72 //============================================ <<  66 G4DiffractiveExcitation::G4DiffractiveExcitation()
 73                                                <<  67 {
 74 //#define debugFTFexcitation                   << 
 75 //#define debug_heavyHadrons                   << 
 76                                                << 
 77 //============================================ << 
 78                                                << 
 79 G4DiffractiveExcitation::G4DiffractiveExcitati << 
 80                                                << 
 81                                                << 
 82 //============================================ << 
 83                                                << 
 84 G4DiffractiveExcitation::~G4DiffractiveExcitat << 
 85                                                << 
 86                                                << 
 87 //============================================ << 
 88                                                << 
 89 G4bool G4DiffractiveExcitation::ExciteParticip << 
 90                                                << 
 91                                                << 
 92                                                << 
 93                                                << 
 94   #ifdef debugFTFexcitation                    << 
 95   G4cout << G4endl << "FTF ExciteParticipants  << 
 96   #endif                                       << 
 97                                                << 
 98   CommonVariables common;                      << 
 99                                                << 
100   // Projectile parameters                     << 
101   common.Pprojectile = projectile->Get4Momentu << 
102   if ( common.Pprojectile.z() < 0.0 ) return f << 
103   common.ProjectilePDGcode = projectile->GetDe << 
104   common.absProjectilePDGcode = std::abs( comm << 
105   common.M0projectile = projectile->GetDefinit << 
106   G4double ProjectileRapidity = common.Pprojec << 
107                                                << 
108   // Target parameters                         << 
109   common.Ptarget = target->Get4Momentum();     << 
110   common.TargetPDGcode = target->GetDefinition << 
111   common.absTargetPDGcode = std::abs( common.T << 
112   common.M0target = target->GetDefinition()->G << 
113   G4double TargetRapidity = common.Ptarget.rap << 
114                                                << 
115   // Kinematical properties of the interaction << 
116   G4LorentzVector Psum = common.Pprojectile +  << 
117   common.S = Psum.mag2();                      << 
118   common.SqrtS = std::sqrt( common.S );        << 
119                                                << 
120   // Check off-shellness of the participants   << 
121   G4bool toBePutOnMassShell = true;  //Uzhi Au << 
122   common.MminProjectile = common.BrW.GetMinimu << 
123   /* Uzhi Aug.2019                             << 
124   if ( common.M0projectile < common.MminProjec << 
125     toBePutOnMassShell = true;                 << 
126     common.M0projectile = common.BrW.SampleMas << 
127                                                << 
128                                                << 
129   }                                            << 
130   */                                           << 
131   common.M0projectile2 = common.M0projectile * << 
132   common.ProjectileDiffStateMinMass    = thePa << 
133   common.ProjectileNonDiffStateMinMass = thePa << 
134   if ( common.M0projectile > common.Projectile << 
135     common.ProjectileDiffStateMinMass    = com << 
136     common.ProjectileNonDiffStateMinMass = com << 
137     if ( common.absProjectilePDGcode > 3000 )  << 
138       common.ProjectileDiffStateMinMass    +=  << 
139       common.ProjectileNonDiffStateMinMass +=  << 
140     }                                          << 
141   }                                            << 
142   common.MminTarget = common.BrW.GetMinimumMas << 
143   /* Uzhi Aug.2019                             << 
144   if ( common.M0target < common.MminTarget ) { << 
145     toBePutOnMassShell = true;                 << 
146     common.M0target = common.BrW.SampleMass( t << 
147                                              t << 
148                                              + << 
149   }                                            << 
150   */                                           << 
151   common.M0target2 = common.M0target * common. << 
152   common.TargetDiffStateMinMass    = theParame << 
153   common.TargetNonDiffStateMinMass = theParame << 
154   if ( common.M0target > common.TargetDiffStat << 
155     common.TargetDiffStateMinMass    = common. << 
156     common.TargetNonDiffStateMinMass = common. << 
157     if ( common.absTargetPDGcode > 3000 ) {  / << 
158       common.TargetDiffStateMinMass    += 140. << 
159       common.TargetNonDiffStateMinMass += 140. << 
160     }                                          << 
161   };                                           << 
162   #ifdef debugFTFexcitation                    << 
163   G4cout << "Proj Targ PDGcodes " << common.Pr << 
164          << "Mprojectile  Y " << common.Pproje << 
165          << "M0projectile Y " << common.M0proj << 
166   G4cout << "Mtarget      Y " << common.Ptarge << 
167          << "M0target     Y " << common.M0targ << 
168   G4cout << "Pproj " << common.Pprojectile <<  << 
169   #endif                                       << 
170                                                << 
171   // Transform momenta to cms and then rotate  << 
172   common.toCms = G4LorentzRotation( -1 * Psum. << 
173   G4LorentzVector Ptmp = common.toCms * common << 
174   if ( Ptmp.pz() <= 0.0 ) return false;  // "S << 
175   common.toCms.rotateZ( -1*Ptmp.phi() );       << 
176   common.toCms.rotateY( -1*Ptmp.theta() );     << 
177   common.toLab = common.toCms.inverse();       << 
178   common.Pprojectile.transform( common.toCms ) << 
179   common.Ptarget.transform( common.toCms );    << 
180                                                << 
181   G4double SumMasses = common.M0projectile + c << 
182   #ifdef debugFTFexcitation                    << 
183   G4cout << "SqrtS     " << common.SqrtS << G4 << 
184          << " " << common.M0target << " " << S << 
185   #endif                                       << 
186   if ( common.SqrtS < SumMasses ) return false << 
187                                                << 
188   common.PZcms2 = ( sqr( common.S ) + sqr( com << 
189                     - 2.0 * ( common.S * ( com << 
190                               + common.M0proje << 
191   #ifdef debugFTFexcitation                    << 
192   G4cout << "PZcms2 after toBePutOnMassShell " << 
193   #endif                                       << 
194   if ( common.PZcms2 < 0.0 ) return false;  // << 
195                                                << 
196   common.PZcms = std::sqrt( common.PZcms2 );   << 
197   if ( toBePutOnMassShell ) {                  << 
198     if ( common.Pprojectile.z() > 0.0 ) {      << 
199       common.Pprojectile.setPz(  common.PZcms  << 
200       common.Ptarget.setPz(     -common.PZcms  << 
201     } else {                                   << 
202       common.Pprojectile.setPz( -common.PZcms  << 
203       common.Ptarget.setPz(      common.PZcms  << 
204     }                                          << 
205     common.Pprojectile.setE( std::sqrt( common << 
206                                         + comm << 
207                                         + comm << 
208                                         + comm << 
209     common.Ptarget.setE( std::sqrt( common.M0t << 
210                                     + common.P << 
211                                     + common.P << 
212                                     + common.P << 
213   }                                            << 
214   #ifdef debugFTFexcitation                    << 
215   G4cout << "Start --------------------" << G4 << 
216          << " " << common.ProjectileDiffStateM << 
217          << G4endl                             << 
218          << "Targ M0 Mdif Mndif " << common.M0 << 
219          << " " << common.TargetNonDiffStateMi << 
220          << "Proj CMS " << common.Pprojectile  << 
221   #endif                                       << 
222                                                << 
223   // Check for possible quark exchange         << 
224   ProjectileRapidity = common.Pprojectile.rapi << 
225   TargetRapidity = common.Ptarget.rapidity();  << 
226   G4double QeNoExc = theParameters->GetProcPro << 
227   G4double QeExc   = theParameters->GetProcPro << 
228                      theParameters->GetProcPro << 
229   common.ProbProjectileDiffraction =           << 
230     theParameters->GetProcProb( 2, ProjectileR << 
231   common.ProbTargetDiffraction =               << 
232     theParameters->GetProcProb( 3, ProjectileR << 
233   common.ProbOfDiffraction = common.ProbProjec << 
234   #ifdef debugFTFexcitation                    << 
235   G4cout << "Proc Probs " << QeNoExc << " " << << 
236          << common.ProbProjectileDiffraction < << 
237          << "ProjectileRapidity " << Projectil << 
238   #endif                                       << 
239                                                << 
240   if ( QeNoExc + QeExc + common.ProbProjectile << 
241     QeNoExc = 1.0 - QeExc - common.ProbProject << 
242   }                                            << 
243   if ( QeExc + QeNoExc != 0.0 ) {              << 
244     common.ProbExc = QeExc / ( QeExc + QeNoExc << 
245   }                                            << 
246   if ( 1.0 - QeExc - QeNoExc > 0.0 ) {         << 
247     common.ProbProjectileDiffraction /= ( 1.0  << 
248     common.ProbTargetDiffraction     /= ( 1.0  << 
249   }                                            << 
250   #ifdef debugFTFexcitation                    << 
251   G4cout << "Proc Probs " << QeNoExc << " " << << 
252          << common.ProbProjectileDiffraction < << 
253          << "ProjectileRapidity " << Projectil << 
254   #endif                                       << 
255                                                << 
256   // Try out charge exchange                   << 
257   G4int returnCode = 1;                        << 
258   if ( G4UniformRand() < QeExc + QeNoExc ) {   << 
259     returnCode = ExciteParticipants_doChargeEx << 
260                                                << 
261   }                                            << 
262                                                << 
263   G4bool returnResult = false;                 << 
264   if ( returnCode == 0 ) {                     << 
265     returnResult = true;  // Successfully ende << 
266   } else if ( returnCode == 1 ) {              << 
267                                                << 
268     common.ProbOfDiffraction = common.ProbProj << 
269     #ifdef debugFTFexcitation                  << 
270     G4cout << "Excitation -------------------- << 
271            << "Proj M0 MdMin MndMin " << commo << 
272            << common.ProjectileDiffStateMinMas << 
273            << G4endl                           << 
274            << "Targ M0 MdMin MndMin " << commo << 
275            << " " << common.TargetNonDiffState << 
276            << G4endl                           << 
277            << "Prob: ProjDiff TargDiff + Sum " << 
278            << common.ProbTargetDiffraction <<  << 
279     #endif                                     << 
280     if ( common.ProbOfDiffraction != 0.0 ) {   << 
281       common.ProbProjectileDiffraction /= comm << 
282     } else {                                   << 
283       common.ProbProjectileDiffraction = 0.0;  << 
284     }                                          << 
285     #ifdef debugFTFexcitation                  << 
286     G4cout << "Prob: ProjDiff TargDiff + Sum " << 
287            << common.ProbTargetDiffraction <<  << 
288     #endif                                     << 
289     common.ProjectileDiffStateMinMass2    = sq << 
290     common.ProjectileNonDiffStateMinMass2 = sq << 
291     common.TargetDiffStateMinMass2        = sq << 
292     common.TargetNonDiffStateMinMass2     = sq << 
293     // Choose between diffraction and non-diff << 
294     if ( G4UniformRand() < common.ProbOfDiffra << 
295                                                << 
296       returnResult = ExciteParticipants_doDiff << 
297                                                << 
298     } else {                                   << 
299                                                << 
300       returnResult = ExciteParticipants_doNonD << 
301                                                << 
302     }                                          << 
303     if ( returnResult ) {                      << 
304       common.Pprojectile += common.Qmomentum;  << 
305       common.Ptarget     -= common.Qmomentum;  << 
306       // Transform back and update SplitableHa << 
307       common.Pprojectile.transform( common.toL << 
308       common.Ptarget.transform( common.toLab ) << 
309       #ifdef debugFTFexcitation                << 
310       G4cout << "Mproj " << common.Pprojectile << 
311              << G4endl;                        << 
312       #endif                                   << 
313       projectile->Set4Momentum( common.Pprojec << 
314       target->Set4Momentum( common.Ptarget );  << 
315       projectile->IncrementCollisionCount( 1 ) << 
316       target->IncrementCollisionCount( 1 );    << 
317     }                                          << 
318   }                                            << 
319                                                << 
320   return returnResult;                         << 
321 }                                                  68 }
322                                                    69 
323 //-------------------------------------------- <<  70 // ---------------------------------------------------------------------
                                                   >>  71 G4bool G4DiffractiveExcitation::
                                                   >>  72   ExciteParticipants(G4VSplitableHadron *projectile, 
                                                   >>  73                      G4VSplitableHadron *target,
                                                   >>  74                      G4FTFParameters    *theParameters,
                                                   >>  75                      G4ElasticHNScattering *theElastic) const  
                                                   >>  76 {
                                                   >>  77 // -------------------- Projectile parameters -----------------------
                                                   >>  78      G4LorentzVector Pprojectile=projectile->Get4Momentum();
                                                   >>  79 
                                                   >>  80      if(Pprojectile.z() < 0.)
                                                   >>  81      {
                                                   >>  82        target->SetStatus(2);
                                                   >>  83        return false;
                                                   >>  84      } 
                                                   >>  85 
                                                   >>  86      G4double ProjectileRapidity = Pprojectile.rapidity();
                                                   >>  87 
                                                   >>  88      G4int    ProjectilePDGcode=projectile->GetDefinition()->GetPDGEncoding();
                                                   >>  89      G4int    absProjectilePDGcode=std::abs(ProjectilePDGcode);
                                                   >>  90 
                                                   >>  91      G4bool PutOnMassShell(false);
                                                   >>  92 //   G4double M0projectile=projectile->GetDefinition()->GetPDGMass(); // With de-excitation
                                                   >>  93      G4double M0projectile = Pprojectile.mag();                       // Without de-excitation
                                                   >>  94 
                                                   >>  95      if(M0projectile < projectile->GetDefinition()->GetPDGMass())
                                                   >>  96      {
                                                   >>  97       PutOnMassShell=true;
                                                   >>  98       M0projectile=projectile->GetDefinition()->GetPDGMass();
                                                   >>  99      }
                                                   >> 100 
                                                   >> 101      G4double M0projectile2 = M0projectile * M0projectile;
                                                   >> 102 
                                                   >> 103      G4double ProjectileDiffStateMinMass=theParameters->GetProjMinDiffMass();
                                                   >> 104      G4double ProjectileNonDiffStateMinMass=theParameters->GetProjMinNonDiffMass();
                                                   >> 105      G4double ProbProjectileDiffraction=theParameters->GetProbabilityOfProjDiff();
                                                   >> 106 
                                                   >> 107 // -------------------- Target parameters -------------------------
                                                   >> 108      G4int    TargetPDGcode=target->GetDefinition()->GetPDGEncoding();
                                                   >> 109      G4int    absTargetPDGcode=std::abs(TargetPDGcode);
                                                   >> 110 
                                                   >> 111      G4LorentzVector Ptarget=target->Get4Momentum();
                                                   >> 112 
                                                   >> 113      G4double M0target = Ptarget.mag();
                                                   >> 114 
                                                   >> 115      G4double TargetRapidity = Ptarget.rapidity();
                                                   >> 116 
                                                   >> 117      if(M0target < target->GetDefinition()->GetPDGMass())
                                                   >> 118      {
                                                   >> 119       PutOnMassShell=true;
                                                   >> 120       M0target=target->GetDefinition()->GetPDGMass();
                                                   >> 121      }
324                                                   122 
325 G4int G4DiffractiveExcitation::                << 123      G4double M0target2 = M0target * M0target; 
326 ExciteParticipants_doChargeExchange( G4VSplita << 124  
327                                      G4VSplita << 125      G4double TargetDiffStateMinMass=theParameters->GetTarMinDiffMass();    
328                                      G4FTFPara << 126      G4double TargetNonDiffStateMinMass=theParameters->GetTarMinNonDiffMass();    
329                                      G4Elastic << 127      G4double ProbTargetDiffraction=theParameters->GetProbabilityOfTarDiff();
330                                      G4Diffrac << 128 
331   // First of the three utility methods used o << 129      G4double AveragePt2=theParameters->GetAveragePt2();
332   // it does the sampling for the charge-excha << 130 
333   // This method returns an integer code - ins << 131      G4double ProbOfDiffraction=ProbProjectileDiffraction +
334   //   "0" : successfully ended and nothing el << 132                                 ProbTargetDiffraction;
335   //   "1" : successfully completed, but the w << 133 
336   //  "99" : unsuccessfully ended, nothing els << 134      G4double SumMasses=M0projectile+M0target+200.*MeV;
337   G4int returnCode = 99;                       << 135 
338                                                << 136 // Kinematical properties of the interactions --------------
339   G4double DeltaProbAtQuarkExchange = theParam << 137      G4LorentzVector Psum;      // 4-momentum in CMS
340   G4ParticleDefinition* TestParticle = 0;      << 138      Psum=Pprojectile+Ptarget;
341   G4double MtestPr = 0.0, MtestTr = 0.0;       << 139      G4double S=Psum.mag2(); 
342                                                << 140 
343   #ifdef debugFTFexcitation                    << 141 // Transform momenta to cms and then rotate parallel to z axis;
344   G4cout << "Q exchange ---------------------- << 142      G4LorentzRotation toCms(-1*Psum.boostVector());
345   #endif                                       << 143 
346                                                << 144      G4LorentzVector Ptmp=toCms*Pprojectile;
347   // The target hadron is always a nucleon (i. << 145 
348   // never an antinucleon), therefore only a q << 146      if ( Ptmp.pz() <= 0. )
349   // exchanged between the projectile hadron a << 147      {
350   // we could get a quark-quark-antiquark syst << 148        target->SetStatus(2); 
351   // This implies that any projectile meson or << 149    // "String" moving backwards in  CMS, abort collision !!
352   // a constituent quark in all cases - can ha << 150        return false;
353   // hadron. Instead, any projectile anti-bary << 151      }
354   // with a target hadron (because it has only << 152 
355   // projectile baryons, instead can have char << 153      toCms.rotateZ(-1*Ptmp.phi());
356                                                << 154      toCms.rotateY(-1*Ptmp.theta());
357   G4int NewProjCode = 0, NewTargCode = 0, Proj << 155 
358   //  Projectile unpacking                     << 156      G4LorentzRotation toLab(toCms.inverse());
359   if ( common.absProjectilePDGcode < 1000 ) {  << 157 
360     UnpackMeson(  common.ProjectilePDGcode, Pr << 158      Pprojectile.transform(toCms);
361   } else {                                     << 159      Ptarget.transform(toCms);
362     UnpackBaryon( common.ProjectilePDGcode, Pr << 160 
363   }                                            << 161      G4double PZcms2, PZcms;
364   //  Target unpacking                         << 162 
365   G4int TargQ1 = 0, TargQ2 = 0, TargQ3 = 0;    << 163      G4double SqrtS=std::sqrt(S);
366   UnpackBaryon( common.TargetPDGcode, TargQ1,  << 164                
367   #ifdef debugFTFexcitation                    << 165      if(absProjectilePDGcode > 1000 && (SqrtS < 2300*MeV || SqrtS < SumMasses))
368   G4cout << "Proj Quarks " << ProjQ1 << " " << << 166      {target->SetStatus(2);  return false;}  // The model cannot work for
369          << "Targ Quarks " << TargQ1 << " " << << 167                                              // p+p-interactions
370   #endif                                       << 168                                              // at Plab < 1.62 GeV/c.
371                                                << 169 
372   // Sampling of exchanged quarks              << 170      if(( absProjectilePDGcode == 211 || ProjectilePDGcode ==  111) && 
373   G4int ProjExchangeQ = 0, TargExchangeQ = 0;  << 171         ((SqrtS < 1600*MeV) || (SqrtS < SumMasses)))
374   if ( common.absProjectilePDGcode < 1000 ) {  << 172      {target->SetStatus(2);  return false;}  // The model cannot work for
375                                                << 173                                              // Pi+p-interactions
376     G4bool isProjQ1Quark = false;              << 174                                              // at Plab < 1. GeV/c.
377     ProjExchangeQ = ProjQ2;                    << 175 
378     if ( ProjQ1 > 0 ) {  // ProjQ1 is a quark  << 176      if(( (absProjectilePDGcode == 321) || (ProjectilePDGcode == -311)   ||
379       isProjQ1Quark = true;                    << 177           (absProjectilePDGcode == 311) || (absProjectilePDGcode == 130) ||
380       ProjExchangeQ = ProjQ1;                  << 178           (absProjectilePDGcode == 310)) && ((SqrtS < 1600*MeV) || (SqrtS < SumMasses))) 
381     }                                          << 179      {target->SetStatus(2);  return false;}  // The model cannot work for
382     // Exchange of non-identical quarks is all << 180                                              // K+p-interactions
383     G4int NpossibleStates = 0;                 << 181                                              // at Plab < ??? GeV/c.  ???
384     if ( ProjExchangeQ != TargQ1 ) NpossibleSt << 182 
385     if ( ProjExchangeQ != TargQ2 ) NpossibleSt << 183      PZcms2=(S*S+M0projectile2*M0projectile2+M0target2*M0target2-
386     if ( ProjExchangeQ != TargQ3 ) NpossibleSt << 184              2*S*M0projectile2 - 2*S*M0target2 - 2*M0projectile2*M0target2)
387     G4int Nsampled = (G4int)G4RandFlat::shootI << 185              /4./S;
388     NpossibleStates = 0;                       << 186 
389     if ( ProjExchangeQ != TargQ1 ) {           << 187      if(PZcms2 < 0)
390       if ( ++NpossibleStates == Nsampled ) {   << 188      {target->SetStatus(2);  return false;}   // It can be in an interaction with 
391         TargExchangeQ = TargQ1; TargQ1 = ProjE << 189                                               // off-shell nuclear nucleon
392         isProjQ1Quark ? ProjQ1 = TargExchangeQ << 190 
393       }                                        << 191      PZcms = std::sqrt(PZcms2);
394     }                                          << 192 
395     if ( ProjExchangeQ != TargQ2 ) {           << 193      if(PutOnMassShell)
396       if ( ++NpossibleStates == Nsampled ) {   << 194      {
397         TargExchangeQ = TargQ2; TargQ2 = ProjE << 195       if(Pprojectile.z() > 0.)
398         isProjQ1Quark ? ProjQ1 = TargExchangeQ << 196       {
399       }                                        << 197        Pprojectile.setPz( PZcms);
400     }                                          << 198        Ptarget.setPz(    -PZcms);
401     if ( ProjExchangeQ != TargQ3 ) {           << 199       }
402       if ( ++NpossibleStates == Nsampled ) {   << 200       else
403         TargExchangeQ = TargQ3; TargQ3 = ProjE << 201       {
404         isProjQ1Quark ? ProjQ1 = TargExchangeQ << 202        Pprojectile.setPz(-PZcms);
405       }                                        << 203        Ptarget.setPz(     PZcms);
406     }                                          << 204       };
407     #ifdef debugFTFexcitation                  << 205 
408     G4cout << "Exchanged Qs in Pr Tr " << Proj << 206       Pprojectile.setE(std::sqrt(M0projectile2                  +
409     #endif                                     << 207                                  Pprojectile.x()*Pprojectile.x()+
410                                                << 208                                  Pprojectile.y()*Pprojectile.y()+
411     G4int aProjQ1 = std::abs( ProjQ1 ), aProjQ << 209                                  PZcms2));
412     G4bool ProjExcited = false;                << 210       Ptarget.setE(std::sqrt(M0target2              +
413     const G4int maxNumberOfAttempts = 50;      << 211                              Ptarget.x()*Ptarget.x()+
414     G4int attempts = 0;                        << 212                              Ptarget.y()*Ptarget.y()+
415     while ( attempts++ < maxNumberOfAttempts ) << 213                              PZcms2));
416                                                << 214      }
417       // Determination of a new projectile ID  << 215 
418       G4double ProbSpin0 = 0.5;                << 216      G4double maxPtSquare; // = PZcms2;
419       G4double Ksi = G4UniformRand();          << 217 
420       if ( aProjQ1 == aProjQ2 ) {              << 218 // Charge exchange can be possible for baryons -----------------
421         if ( G4UniformRand() < ProbSpin0 ) {   << 219 
422           if ( aProjQ1 < 3 ) {                 << 220 // Getting the values needed for exchange ----------------------
423             NewProjCode = 111;                 << 221      G4double MagQuarkExchange        =theParameters->GetMagQuarkExchange();
424             if ( Ksi < 0.5 ) {                 << 222      G4double SlopeQuarkExchange      =theParameters->GetSlopeQuarkExchange();
425               NewProjCode = 221;               << 223      G4double DeltaProbAtQuarkExchange=theParameters->GetDeltaProbAtQuarkExchange();
426               if ( Ksi < 0.25 ) {              << 224 
427                 NewProjCode = 331;             << 225 //     G4double NucleonMass=
428               }                                << 226 //              (G4ParticleTable::GetParticleTable()->FindParticle(2112))->GetPDGMass();     
429             }                                  << 227      G4double DeltaMass=
430           } else if ( aProjQ1 == 3 ) {         << 228               (G4ParticleTable::GetParticleTable()->FindParticle(2224))->GetPDGMass();
431               NewProjCode = 221;               << 229 
432               if ( Ksi < 0.5 ) {               << 230 // Check for possible quark excjane -----------------------------------
433                 NewProjCode = 331;             << 231      if(G4UniformRand() < MagQuarkExchange*
434               }                                << 232         std::exp(-SlopeQuarkExchange*(ProjectileRapidity - TargetRapidity)))
435           } else if ( aProjQ1 == 4 ) {         << 233      {    
436       NewProjCode = 441;                // eta << 234       G4int NewProjCode(0), NewTargCode(0);
437     } else if ( aProjQ1 == 5 ) {               << 235 
438       NewProjCode = 551;                // eta << 236       G4int ProjQ1(0), ProjQ2(0), ProjQ3(0);
439     }                                          << 237 
440         } else {                               << 238 //  Projectile unpacking --------------------------
441           if ( aProjQ1 < 3 ) {                 << 239       if(absProjectilePDGcode < 1000 )
442             NewProjCode = 113;                 << 240       {    // projectile is meson ----------------- 
443             if ( Ksi < 0.5 ) {                 << 241        UnpackMeson(ProjectilePDGcode, ProjQ1, ProjQ2);  
444               NewProjCode = 223;               << 242       } else 
445             }                                  << 243       {    // projectile is baryon ----------------
446           } else if ( aProjQ1 == 3 ) {         << 244        UnpackBaryon(ProjectilePDGcode, ProjQ1, ProjQ2, ProjQ3);
447             NewProjCode = 333;                 << 245       } // End of the hadron's unpacking ----------
448           } else if ( aProjQ1 == 4 ) {         << 246 
449       NewProjCode = 443;                // J/p << 247 //  Target unpacking ------------------------------
450     } else if ( aProjQ1 == 5 ) {               << 248       G4int TargQ1(0), TargQ2(0), TargQ3(0);
451       NewProjCode = 553;                // Ups << 249       UnpackBaryon(TargetPDGcode, TargQ1, TargQ2, TargQ3); 
452     }                                          << 250 
                                                   >> 251 // Sampling of exchanged quarks -------------------
                                                   >> 252       G4int ProjExchangeQ(0);
                                                   >> 253       G4int TargExchangeQ(0);
                                                   >> 254 
                                                   >> 255       if(absProjectilePDGcode < 1000 )
                                                   >> 256       {    // projectile is meson ----------------- 
                                                   >> 257        if(ProjQ1 > 0 ) // ProjQ1 is quark
                                                   >> 258        {  
                                                   >> 259         ProjExchangeQ = ProjQ1;
                                                   >> 260         if(ProjExchangeQ != TargQ1)
                                                   >> 261         {
                                                   >> 262          TargExchangeQ = TargQ1; TargQ1=ProjExchangeQ; ProjQ1=TargExchangeQ;
                                                   >> 263         } else
                                                   >> 264         if(ProjExchangeQ != TargQ2)
                                                   >> 265         {
                                                   >> 266          TargExchangeQ = TargQ2; TargQ2=ProjExchangeQ; ProjQ1=TargExchangeQ;
                                                   >> 267         } else          
                                                   >> 268         {
                                                   >> 269          TargExchangeQ = TargQ3;  TargQ3=ProjExchangeQ; ProjQ1=TargExchangeQ;
453         }                                         270         }
454       } else {                                 << 271        } else          // ProjQ2 is quark
455         if ( aProjQ1 > aProjQ2 ) {             << 272        {  
456           NewProjCode = aProjQ1*100 + aProjQ2* << 273         ProjExchangeQ = ProjQ2;
457         } else {                               << 274         if(ProjExchangeQ != TargQ1)
458           NewProjCode = aProjQ2*100 + aProjQ1* << 275         {
                                                   >> 276          TargExchangeQ = TargQ1; TargQ1=ProjExchangeQ; ProjQ2=TargExchangeQ;
                                                   >> 277         } else
                                                   >> 278         if(ProjExchangeQ != TargQ2)
                                                   >> 279         {
                                                   >> 280          TargExchangeQ = TargQ2; TargQ2=ProjExchangeQ; ProjQ2=TargExchangeQ;
                                                   >> 281         } else 
                                                   >> 282         {
                                                   >> 283          TargExchangeQ = TargQ3;  TargQ3=ProjExchangeQ; ProjQ2=TargExchangeQ;
459         }                                         284         }
460       }                                        << 
461       #ifdef debugFTFexcitation                << 
462       G4cout << "NewProjCode " << NewProjCode  << 
463       #endif                                   << 
464       // Decide (with 50% probability) whether << 
465       // but not in the case of charmed and bo << 
466       // there are no excited charmed and bott << 
467       ProjExcited = false;                     << 
468       if ( aProjQ1 <= 3  &&  aProjQ2 <= 3  &&  << 
469         NewProjCode += 2;  // Excited meson (J << 
470         ProjExcited = true;                    << 
471       }                                        << 
472                                                   285 
473       // So far we have used the absolute valu << 286        } // End of if(ProjQ1 > 0 ) // ProjQ1 is quark
474       // now look at their signed values to se << 
475       G4int value = ProjQ1, absValue = aProjQ1 << 
476       for ( G4int iQ = 0; iQ < 2; ++iQ ) {  // << 
477         if ( iQ == 1 ) {                       << 
478           value = ProjQ2; absValue = aProjQ2;  << 
479         }                                      << 
480         if ( absValue == 2  ||  absValue == 4  << 
481           Qquarks += 2*value/absValue;  // u,  << 
482         } else {                               << 
483           Qquarks -= value/absValue;    // d,  << 
484         }                                      << 
485       }                                        << 
486       // If Qquarks is positive, the sign of N << 
487       // If Qquarks is negative, then the sign << 
488       // If Qquarks is zero, then we need to d << 
489       // 1. If aProjQ1 and aProjQ2 are the sam << 
490       //    (because the antiparticle is the s << 
491       // 2. If aProjQ1 and aProjQ2 are not the << 
492       //    we have only two possibilities:    << 
493       //    a. aProjQ1 and aProjQ2 are two dif << 
494       //       (s,d) or (b,d), or (b,s). In th << 
495       //       is fine (because the heaviest o << 
496       //       to be anti-quark belonging to t << 
497       //       implies a meson with positive P << 
498       //    b. aProjQ1 and aProjQ2 are two dif << 
499       //       The heaviest of the two (c) has << 
500       //       in the projectile, therefore th << 
501       //       reverse: 421 -> -421 anti_D0 (c << 
502       if ( Qquarks < 0  || ( Qquarks == 0  &&  << 
503   NewProjCode *= -1;                           << 
504       }                                        << 
505       #ifdef debugFTFexcitation                << 
506       G4cout << "NewProjCode +2 or 0 " << NewP << 
507       G4cout<<"+++++++++++++++++++++++++++++++ << 
508       G4cout<<ProjQ1<<" "<<ProjQ2<<" "<<Qquark << 
509       G4cout<<"+++++++++++++++++++++++++++++++ << 
510       #endif                                   << 
511                                                << 
512       // Proj                                  << 
513       TestParticle = G4ParticleTable::GetParti << 
514       if ( ! TestParticle ) continue;          << 
515       common.MminProjectile = common.BrW.GetMi << 
516       if ( common.SqrtS - common.M0target < co << 
517       MtestPr = common.BrW.SampleMass( TestPar << 
518                                                << 
519       #ifdef debugFTFexcitation                << 
520       G4cout << "TestParticle Name " << NewPro << 
521              << G4endl                         << 
522              << "MtestPart MtestPart0 "<<Mtest << 
523              << "M0projectile projectile PDGMa << 
524              << projectile->GetDefinition()->G << 
525       #endif                                   << 
526                                                << 
527       // Targ                                  << 
528       NewTargCode = NewNucleonId( TargQ1, Targ << 
529       #ifdef debugFTFexcitation                << 
530       G4cout << "New TrQ " << TargQ1 << " " << << 
531              << "NewTargCode " << NewTargCode  << 
532       #endif                                   << 
533                                                << 
534       // If the target has not heavy (charm or << 
535       // see whether a Delta isobar can be cre << 
536       if ( TargQ1 <= 3  &&  TargQ2 <= 3  &&  T << 
537         if ( TargQ1 != TargQ2  &&  TargQ1 != T << 
538           if ( G4UniformRand() < 0.5 ) {       << 
539             NewTargCode += 2;                  << 
540           } else if ( G4UniformRand() < 0.75 ) << 
541             NewTargCode = 3122;  // Lambda     << 
542           }                                    << 
543         } else if ( TargQ1 == TargQ2  &&  Targ << 
544           NewTargCode += 2; ProjExcited = true << 
545         } else if ( target->GetDefinition()->G << 
546           if ( G4UniformRand() > DeltaProbAtQu << 
547             NewTargCode += 2; ProjExcited = tr << 
548           }                                    << 
549         } else if ( ! ProjExcited  &&          << 
550                     G4UniformRand() < DeltaPro << 
551                     common.SqrtS > common.M0pr << 
552                          G4ParticleTable::GetP << 
553           NewTargCode += 2;  // Create Delta i << 
554         }                                      << 
555       }                                        << 
556                                                   287 
557       // Protection against:                   << 288        G4int aProjQ1=std::abs(ProjQ1);
558       // - Lambda* (i.e. excited Lambda state) << 289        G4int aProjQ2=std::abs(ProjQ2);
559       // - Sigma*  (i.e. excited Sigma hyperon << 290        if(aProjQ1 == aProjQ2)          {NewProjCode = 111;} // Pi0-meson
560       // - Xi*     (i.e. excited Xi hyperon st << 291        else  // |ProjQ1| # |ProjQ2|
561       if ( NewTargCode == 3124  ||   // Lambda << 292        {
562      NewTargCode == 3224  ||   // Sigma*+ NOT  << 293         if(aProjQ1 > aProjQ2)          {NewProjCode = aProjQ1*100+aProjQ2*10+1;}
563      NewTargCode == 3214  ||   // Sigma*0 NOT  << 294         else                           {NewProjCode = aProjQ2*100+aProjQ1*10+1;}
564      NewTargCode == 3114  ||   // Sigma*- NOT  << 295        }
565      NewTargCode == 3324  ||   // Xi*0    NOT  << 296 
566      NewTargCode == 3314  ) {  // Xi*-    NOT  << 297        NewTargCode = NewNucleonId(TargQ1, TargQ2, TargQ3);
567   //G4cout << "G4DiffractiveExcitation::Excite << 298 
568   //       << NewTargCode << G4endl;           << 299        if( (TargQ1 == TargQ2) && (TargQ1 == TargQ3) &&
569   NewTargCode -= 2;  // Corresponding ground-s << 300            (SqrtS > M0projectile+DeltaMass))           {NewTargCode +=2;} //Create Delta isobar
570       }                                        << 301 
571                                                << 302        else if( target->GetDefinition()->GetPDGiIsospin() == 3 )         //Delta was the target
572       // Special treatment for charmed and bot << 303        { if(G4UniformRand() > DeltaProbAtQuarkExchange){NewTargCode +=2;} //Save   Delta isobar
573       // so we need to transform them by hand  << 304          else                                          {}     // De-excite initial Delta isobar
574       #ifdef debug_heavyHadrons                << 305        }
575       G4int initialNewTargCode = NewTargCode;  << 306 
576       #endif                                   << 307        else if((G4UniformRand() < DeltaProbAtQuarkExchange) &&         //Nucleon was the target
577       if      ( NewTargCode == 4322 ) NewTargC << 308                (SqrtS > M0projectile+DeltaMass))      {NewTargCode +=2;}  //Create Delta isobar
578       else if ( NewTargCode == 4312 ) NewTargC << 309        else                                           {}                 //Save initial nucleon
579       else if ( NewTargCode == 5312 ) NewTargC << 310 //       target->SetDefinition(                                          // Fix 15.12.09
580       else if ( NewTargCode == 5322 ) NewTargC << 311 //       G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode));// Fix 15.12.09 
581       #ifdef debug_heavyHadrons                << 312       } else 
582       if ( NewTargCode != initialNewTargCode ) << 313       {    // projectile is baryon ----------------
583   G4cout << "G4DiffractiveExcitation::ExcitePa << 314        G4double Same=0.; //0.3; //0.5;
584          << "\t target heavy baryon with pdgCo << 315        G4bool ProjDeltaHasCreated(false);
585          << " into pdgCode=" << NewTargCode << << 316        G4bool TargDeltaHasCreated(false);
586       }                                        << 317  
587       #endif                                   << 318        G4double Ksi=G4UniformRand();
588                                                << 319        if(G4UniformRand() < 0.5)     // Sampling exchange quark from proj. or targ.
589       TestParticle = G4ParticleTable::GetParti << 320        {   // Sampling exchanged quark from the projectile ---
590       if ( ! TestParticle ) continue;          << 321         if( Ksi < 0.333333 ) 
591       #ifdef debugFTFexcitation                << 322         {ProjExchangeQ = ProjQ1;}
592       G4cout << "New targ " << NewTargCode <<  << 323         else if( (0.333333 <= Ksi) && (Ksi < 0.666667))
593       #endif                                   << 324         {ProjExchangeQ = ProjQ2;}
594       common.MminTarget = common.BrW.GetMinimu << 325         else
595       if ( common.SqrtS - MtestPr < common.Mmi << 326         {ProjExchangeQ = ProjQ3;}
596       MtestTr = common.BrW.SampleMass( TestPar << 327 
597                                                << 328         if((ProjExchangeQ != TargQ1)||(G4UniformRand()<Same)) 
598       if ( common.SqrtS > MtestPr + MtestTr )  << 329         {
599                                                << 330          TargExchangeQ = TargQ1; TargQ1=ProjExchangeQ; ProjExchangeQ=TargExchangeQ;
600     }  // End of while loop                    << 331         } else
601     if ( attempts >= maxNumberOfAttempts ) ret << 332         if((ProjExchangeQ != TargQ2)||(G4UniformRand()<Same)) 
602                                                << 333         {
603     if ( MtestPr >= common.Pprojectile.mag()   << 334          TargExchangeQ = TargQ2; TargQ2=ProjExchangeQ; ProjExchangeQ=TargExchangeQ;
604       common.M0projectile = MtestPr;           << 335         } else 
605     }                                          << 336         {
606     #ifdef debugFTFexcitation                  << 337          TargExchangeQ = TargQ3;  TargQ3=ProjExchangeQ; ProjExchangeQ=TargExchangeQ;
607     G4cout << "M0projectile After check " << c << 
608     #endif                                     << 
609     common.M0projectile2 = common.M0projectile << 
610     common.ProjectileDiffStateMinMass    = com << 
611     common.ProjectileNonDiffStateMinMass = com << 
612     if ( MtestTr >= common.Ptarget.mag()  ||   << 
613       common.M0target = MtestTr;               << 
614     }                                          << 
615     common.M0target2 = common.M0target * commo << 
616     #ifdef debugFTFexcitation                  << 
617     G4cout << "New targ M0 M0^2 " << common.M0 << 
618     #endif                                     << 
619     common.TargetDiffStateMinMass    = common. << 
620     common.TargetNonDiffStateMinMass = common. << 
621                                                << 
622   } else {  // of the if ( common.absProjectil << 
623                                                << 
624     // If it is a projectile anti-baryon, no q << 
625     // therefore returns immediately 1 (which  << 
626     // needs to be continued").                << 
627     if ( common.ProjectilePDGcode < 0 ) return << 
628                                                << 
629     // Choose randomly, with equal probability << 
630     // projectile or target hadron for selecti << 
631     G4bool isProjectileExchangedQ = false;     << 
632     G4int firstQ      = TargQ1, secondQ      = << 
633     G4int otherFirstQ = ProjQ1, otherSecondQ = << 
634     if ( G4UniformRand() < 0.5 ) {             << 
635       isProjectileExchangedQ = true;           << 
636       firstQ      = ProjQ1; secondQ      = Pro << 
637       otherFirstQ = TargQ1; otherSecondQ = Tar << 
638     }                                          << 
639     // Choose randomly, with equal probability << 
640     // selected (projectile or target) hadron  << 
641     G4int exchangedQ = 0;                      << 
642     G4double Ksi = G4UniformRand();            << 
643     if ( Ksi < 0.333333 ) {                    << 
644       exchangedQ = firstQ;                     << 
645     } else if ( 0.333333 <= Ksi  &&  Ksi < 0.6 << 
646       exchangedQ = secondQ;                    << 
647     } else {                                   << 
648       exchangedQ = thirdQ;                     << 
649     }                                          << 
650     #ifdef debugFTFexcitation                  << 
651     G4cout << "Exchange Qs isProjectile Q " << << 
652     #endif                                     << 
653     // The exchanged quarks (one of the projec << 
654     // are always accepted if they have differ << 
655     // same flavour) they are accepted only wi << 
656     G4double probSame = theParameters->GetProb << 
657     const G4int MaxCount = 100;                << 
658     G4int count = 0, otherExchangedQ = 0;      << 
659     do {                                       << 
660       if ( exchangedQ != otherFirstQ  ||  G4Un << 
661         otherExchangedQ = otherFirstQ; otherFi << 
662       } else {                                 << 
663         if ( exchangedQ != otherSecondQ  ||  G << 
664           otherExchangedQ = otherSecondQ; othe << 
665         } else {                               << 
666           if ( exchangedQ != otherThirdQ  ||   << 
667             otherExchangedQ = otherThirdQ; oth << 
668           }                                    << 
669         }                                      << 
670       }                                        << 
671     } while ( otherExchangedQ == 0  &&  ++coun << 
672     if ( count >= MaxCount ) return returnCode << 
673     // Swap (between projectile and target had << 
674     // as "exchanged" quarks.                  << 
675     if ( Ksi < 0.333333 ) {                    << 
676       firstQ = exchangedQ;                     << 
677     } else if ( 0.333333 <= Ksi  &&  Ksi < 0.6 << 
678       secondQ = exchangedQ;                    << 
679     } else {                                   << 
680       thirdQ = exchangedQ;                     << 
681     }                                          << 
682     if ( isProjectileExchangedQ ) {            << 
683       ProjQ1 = firstQ;      ProjQ2 = secondQ;  << 
684       TargQ1 = otherFirstQ; TargQ2 = otherSeco << 
685     } else {                                   << 
686       TargQ1 = firstQ;      TargQ2 = secondQ;  << 
687       ProjQ1 = otherFirstQ; ProjQ2 = otherSeco << 
688     }                                          << 
689     #ifdef debugFTFexcitation                  << 
690     G4cout << "Exchange Qs Pr  Tr " << ( isPro << 
691            << " " << ( isProjectileExchangedQ  << 
692     #endif                                     << 
693                                                << 
694     NewProjCode = NewNucleonId( ProjQ1, ProjQ2 << 
695     NewTargCode = NewNucleonId( TargQ1, TargQ2 << 
696     // Decide whether the new projectile hadro << 
697     // then decide whether the new target hadr << 
698     // Notice that a Delta particle has the la << 
699     // whereas a nucleon has "2" (because its  << 
700     // If the new projectile hadron or the new << 
701     // constituent quark, then skip this part  << 
702     // excited charm and bottom hadrons).      << 
703     for ( G4int iHadron = 0; iHadron < 2; iHad << 
704       // First projectile hadron, then target  << 
705       G4int codeQ1 = ProjQ1, codeQ2 = ProjQ2,  << 
706       G4double massConstraint = common.M0targe << 
707       G4bool isHadronADelta = ( projectile->Ge << 
708       if ( iHadron == 1 ) {  // Target hadron  << 
709         codeQ1 = TargQ1, codeQ2 = TargQ2, code << 
710         massConstraint = common.M0projectile;  << 
711         isHadronADelta = ( target->GetDefiniti << 
712       }                                        << 
713       if ( codeQ1 > 3  ||  codeQ2 > 3  ||  cod << 
714       if ( codeQ1 == codeQ2  &&  codeQ1 == cod << 
715         newHadCode += 2;  // Delta++ (uuu) or  << 
716       } else if ( isHadronADelta ) {  // Hadro << 
717         if ( G4UniformRand() > DeltaProbAtQuar << 
718           newHadCode += 2;  // Delta+ (uud) or << 
719         } else {                               << 
720           newHadCode += 0;  // No delta (so th << 
721         }                                      << 
722       } else {  // Hadron (projectile or targe << 
723         if ( G4UniformRand() < DeltaProbAtQuar << 
724              common.SqrtS > G4ParticleTable::G << 
725                             + massConstraint ) << 
726           newHadCode += 2;  // Delta+ (uud) or << 
727         } else {                               << 
728           newHadCode += 0;  // No delta (so th << 
729         }                                         338         }
730       }                                        << 
731       if ( iHadron == 0 ) {  // Projectile had << 
732         NewProjCode = newHadCode;              << 
733       } else {               // Target hadron  << 
734         NewTargCode = newHadCode;              << 
735       }                                        << 
736     }                                          << 
737     #ifdef debugFTFexcitation                  << 
738     G4cout << "NewProjCode NewTargCode " << Ne << 
739     #endif                                     << 
740                                                << 
741     // Protection against:                     << 
742     // - Lambda* (i.e. excited Lambda state)   << 
743     // - Sigma*  (i.e. excited Sigma hyperon s << 
744     // - Xi*     (i.e. excited Xi hyperon stat << 
745     if ( NewProjCode == 3124  ||   // Lambda*  << 
746    NewProjCode == 3224  ||   // Sigma*+ NOT ex << 
747    NewProjCode == 3214  ||   // Sigma*0 NOT ex << 
748    NewProjCode == 3114  ||   // Sigma*- NOT ex << 
749    NewProjCode == 3324  ||   // Xi*0    NOT ex << 
750    NewProjCode == 3314  ) {  // Xi*-    NOT ex << 
751       //G4cout << "G4DiffractiveExcitation::Ex << 
752       //       << NewProjCode << G4endl;       << 
753       NewProjCode -= 2;  // Corresponding grou << 
754     }                                          << 
755     if ( NewTargCode == 3124  ||   // Lambda*  << 
756    NewTargCode == 3224  ||   // Sigma*+ NOT ex << 
757    NewTargCode == 3214  ||   // Sigma*0 NOT ex << 
758    NewTargCode == 3114  ||   // Sigma*- NOT ex << 
759    NewTargCode == 3324  ||   // Xi*0    NOT ex << 
760    NewTargCode == 3314  ) {  // Xi*-    NOT ex << 
761       //G4cout << "G4DiffractiveExcitation::Ex << 
762       //       << NewTargCode << G4endl;       << 
763       NewTargCode -= 2;  // Corresponding grou << 
764     }                                          << 
765                                                   339 
766     // Special treatment for charmed and botto << 340         if( Ksi < 0.333333 ) 
767     // so we need to transform them by hand to << 341         {ProjQ1=ProjExchangeQ;}
768     #ifdef debug_heavyHadrons                  << 342         else if( (0.333333 <= Ksi) && (Ksi < 0.666667))
769     G4int initialNewProjCode = NewProjCode, in << 343         {ProjQ2=ProjExchangeQ;}
770     #endif                                     << 344         else
771     if      ( NewProjCode == 4322 ) NewProjCod << 345         {ProjQ3=ProjExchangeQ;}
772     else if ( NewProjCode == 4312 ) NewProjCod << 346 
773     else if ( NewProjCode == 5312 ) NewProjCod << 347        } else
774     else if ( NewProjCode == 5322 ) NewProjCod << 348        {   // Sampling exchanged quark from the target -------
775     if      ( NewTargCode == 4322 ) NewTargCod << 349         if( Ksi < 0.333333 ) 
776     else if ( NewTargCode == 4312 ) NewTargCod << 350         {TargExchangeQ = TargQ1;}
777     else if ( NewTargCode == 5312 ) NewTargCod << 351         else if( (0.333333 <= Ksi) && (Ksi < 0.666667)) 
778     else if ( NewTargCode == 5322 ) NewTargCod << 352         {TargExchangeQ = TargQ2;}
779     #ifdef debug_heavyHadrons                  << 353         else
780     if ( NewProjCode != initialNewProjCode  || << 354         {TargExchangeQ = TargQ3;}
781       G4cout << "G4DiffractiveExcitation::Exci << 355 
782              << "\t heavy baryon into an exist << 356         if((TargExchangeQ != ProjQ1)||(G4UniformRand()<Same)) 
783       if ( NewProjCode != initialNewProjCode ) << 357         {
784   G4cout << "\t \t projectile baryon with pdgC << 358          ProjExchangeQ = ProjQ1; ProjQ1=TargExchangeQ; TargExchangeQ=ProjExchangeQ;
785          << " into pdgCode=" << NewProjCode << << 359         } else
786       }                                        << 360         if((TargExchangeQ != ProjQ2)||(G4UniformRand()<Same)) 
787       if ( NewTargCode != initialNewTargCode ) << 361         {
788         G4cout << "\t \t target baryon with pd << 362          ProjExchangeQ = ProjQ2; ProjQ2=TargExchangeQ; TargExchangeQ=ProjExchangeQ;
789          << " into pdgCode=" << NewTargCode << << 363         } else 
790       }                                        << 364         {
791     }                                          << 365          ProjExchangeQ = ProjQ3;  ProjQ3=TargExchangeQ; TargExchangeQ=ProjExchangeQ;
792     #endif                                     << 
793                                                << 
794     // Sampling of the masses of the projectil << 
795     // Because of energy conservation, the ord << 
796     // randomly, half of the time we sample fi << 
797     // then the projectile nucleon mass, and t << 
798     // sample first the projectile nucleon mas << 
799     G4VSplitableHadron* firstHadron  = target; << 
800     G4VSplitableHadron* secondHadron = project << 
801     G4int firstHadronCode = NewTargCode, secon << 
802     G4double massConstraint = common.M0project << 
803     G4bool isFirstTarget = true;               << 
804     if ( G4UniformRand() < 0.5 ) {             << 
805       // Sample first the projectile nucleon m << 
806       firstHadron  = projectile;      secondHa << 
807       firstHadronCode = NewProjCode;  secondHa << 
808       massConstraint = common.M0target;        << 
809       isFirstTarget = false;                   << 
810     }                                          << 
811     G4double Mtest1st = 0.0, Mtest2nd = 0.0, M << 
812     for ( int iSamplingCase = 0; iSamplingCase << 
813       G4VSplitableHadron* aHadron = firstHadro << 
814       G4int aHadronCode = firstHadronCode;     << 
815       if ( iSamplingCase == 1 ) {  // Second n << 
816         aHadron = secondHadron; aHadronCode =  << 
817       }                                        << 
818       G4double MtestHadron = 0.0, MminHadron = << 
819       if ( aHadron->GetStatus() == 1  ||  aHad << 
820         TestParticle = G4ParticleTable::GetPar << 
821         if ( ! TestParticle ) return returnCod << 
822         MminHadron = common.BrW.GetMinimumMass << 
823         if ( common.SqrtS - massConstraint < M << 
824         if ( TestParticle->GetPDGWidth() == 0. << 
825           MtestHadron = common.BrW.SampleMass( << 
826         } else {                               << 
827           const G4int maxNumberOfAttempts = 50 << 
828           G4int attempts = 0;                  << 
829           while ( attempts < maxNumberOfAttemp << 
830             attempts++;                        << 
831             MtestHadron = common.BrW.SampleMas << 
832                                                << 
833             if ( common.SqrtS < MtestHadron +  << 
834               continue;  // Kinematically unac << 
835             } else {                           << 
836               break;     // Kinematically acce << 
837             }                                  << 
838           }                                    << 
839           if ( attempts >= maxNumberOfAttempts << 
840         }                                         366         }
841       }                                        << 
842       if ( iSamplingCase == 0 ) {              << 
843         Mtest1st = MtestHadron;  Mmin1st = Mmi << 
844       } else {                                 << 
845         Mtest2nd = MtestHadron;  Mmin2nd = Mmi << 
846       }                                        << 
847     }  // End for loop on the two sampling cas << 
848     if ( isFirstTarget ) {                     << 
849       MtestTr = Mtest1st;    MtestPr = Mtest2n << 
850       common.MminTarget = Mmin1st;  common.Mmi << 
851     } else {                                   << 
852       MtestTr = Mtest2nd;    MtestPr = Mtest1s << 
853       common.MminTarget = Mmin2nd;  common.Mmi << 
854     }                                          << 
855                                                << 
856     if ( MtestPr != 0.0 ) {                    << 
857       common.M0projectile = MtestPr;           << 
858       common.M0projectile2 = common.M0projecti << 
859       common.ProjectileDiffStateMinMass =    c << 
860       common.ProjectileNonDiffStateMinMass = c << 
861     }                                          << 
862     if ( MtestTr != 0.0 ) {                    << 
863       common.M0target = MtestTr;               << 
864       common.M0target2 = common.M0target * com << 
865       common.TargetDiffStateMinMass    = commo << 
866       common.TargetNonDiffStateMinMass = commo << 
867     }                                          << 
868                                                << 
869   }  // End of if ( common.absProjectilePDGcod << 
870                                                << 
871   // If we assume that final state hadrons aft << 
872   // in the ground states, we have to put      << 
873   if ( common.SqrtS < common.M0projectile + co << 
874                                                << 
875   common.PZcms2 = ( sqr( common.S ) + sqr( com << 
876                     - 2.0 * ( common.S * ( com << 
877                               + common.M0proje << 
878   #ifdef debugFTFexcitation                    << 
879   G4cout << "At the end// NewProjCode " << New << 
880          << "At the end// NewTargCode " << New << 
881          << "M0pr  M0tr  SqS " << common.M0pro << 
882          << common.SqrtS << G4endl             << 
883          << "M0pr2 M0tr2 SqS " << common.M0pro << 
884          << common.SqrtS << G4endl             << 
885          << "PZcms2 after the change " << comm << 
886   #endif                                       << 
887   if ( common.PZcms2 < 0.0 ) return returnCode << 
888                                                << 
889   projectile->SetDefinition( G4ParticleTable:: << 
890   target->SetDefinition( G4ParticleTable::GetP << 
891   common.PZcms = std::sqrt( common.PZcms2 );   << 
892   common.Pprojectile.setPz( common.PZcms );    << 
893   common.Pprojectile.setE( std::sqrt( common.M << 
894   common.Ptarget.setPz(    -common.PZcms );    << 
895   common.Ptarget.setE( std::sqrt( common.M0tar << 
896   #ifdef debugFTFexcitation                    << 
897   G4cout << "Proj Targ and Proj+Targ in CMS" < << 
898          << common.Ptarget << G4endl << common << 
899   #endif                                       << 
900                                                << 
901   if ( projectile->GetStatus() != 0 ) projecti << 
902   if ( target->GetStatus()     != 0 ) target-> << 
903                                                << 
904   // Check for possible excitation of the part << 
905   if ( common.SqrtS < common.M0projectile + co << 
906        common.SqrtS < common.ProjectileDiffSta << 
907        common.ProbOfDiffraction == 0.0 ) commo << 
908                                                << 
909   if ( G4UniformRand() > common.ProbExc ) {  / << 
910     #ifdef debugFTFexcitation                  << 
911     G4cout << "Make elastic scattering of new  << 
912     #endif                                     << 
913     common.Pprojectile.transform( common.toLab << 
914     common.Ptarget.transform( common.toLab );  << 
915     projectile->Set4Momentum( common.Pprojecti << 
916     target->Set4Momentum( common.Ptarget );    << 
917     G4bool Result = theElastic->ElasticScatter << 
918     #ifdef debugFTFexcitation                  << 
919     G4cout << "Result of el. scatt " << Result << 
920            << G4endl << projectile->Get4Moment << 
921            << projectile->Get4Momentum() + tar << 
922            << (projectile->Get4Momentum() + ta << 
923     #endif                                     << 
924     if ( Result ) returnCode = 0;  // successf << 
925     return returnCode;                         << 
926   }                                            << 
927                                                << 
928   #ifdef debugFTFexcitation                    << 
929   G4cout << "Make excitation of new hadrons" < << 
930   #endif                                       << 
931                                                << 
932   // Redefinition of ProbOfDiffraction because << 
933   common.ProbOfDiffraction = common.ProbProjec << 
934   if ( common.ProbOfDiffraction != 0.0 ) {     << 
935     common.ProbProjectileDiffraction /= common << 
936     common.ProbTargetDiffraction     /= common << 
937   }                                            << 
938                                                << 
939   return returnCode = 1;  // successfully comp << 
940 }                                              << 
941                                                << 
942 //-------------------------------------------- << 
943                                                << 
944 G4bool G4DiffractiveExcitation::               << 
945 ExciteParticipants_doDiffraction( G4VSplitable << 
946                                   G4FTFParamet << 
947                                   G4Diffractiv << 
948   // Second of the three utility methods used  << 
949   // it does the sampling for the diffraction  << 
950                                                << 
951   G4bool isProjectileDiffraction = false;      << 
952   if ( G4UniformRand() < common.ProbProjectile << 
953     isProjectileDiffraction = true;            << 
954     #ifdef debugFTFexcitation                  << 
955     G4cout << "projectile diffraction" << G4en << 
956     #endif                                     << 
957     common.ProjMassT2 = common.ProjectileDiffS << 
958     common.ProjMassT  = common.ProjectileDiffS << 
959     common.TargMassT2 = common.M0target2;      << 
960     common.TargMassT  = common.M0target;       << 
961   } else {                                     << 
962     #ifdef debugFTFexcitation                  << 
963     G4cout << "Target diffraction" << G4endl;  << 
964     #endif                                     << 
965     common.ProjMassT2 = common.M0projectile2;  << 
966     common.ProjMassT  = common.M0projectile;   << 
967     common.TargMassT2 = common.TargetDiffState << 
968     common.TargMassT  = common.TargetDiffState << 
969   }                                            << 
970                                                << 
971   // Check whether the interaction is possible << 
972   if ( common.SqrtS < common.ProjMassT + commo << 
973                                                << 
974   common.PZcms2 = ( sqr( common.S ) + sqr( com << 
975                     - 2.0 * ( common.S * ( com << 
976                           + common.ProjMassT2  << 
977   if ( common.PZcms2 < 0.0 ) return false;     << 
978   common.maxPtSquare = common.PZcms2;          << 
979                                                << 
980   G4double DiffrAveragePt2 = theParameters->Ge << 
981   G4bool loopCondition = true;                 << 
982   G4int whilecount = 0;                        << 
983   do {  // Generate pt and mass of projectile  << 
984                                                << 
985     whilecount++;                              << 
986     if ( whilecount > 1000 ) {                 << 
987       common.Qmomentum = G4LorentzVector( 0.0, << 
988       return false;  //  Ignore this interacti << 
989     };                                         << 
990                                                << 
991     common.Qmomentum = G4LorentzVector( Gaussi << 
992     common.Pt2 = G4ThreeVector( common.Qmoment << 
993     if ( isProjectileDiffraction ) {  // proje << 
994       common.ProjMassT2 = common.ProjectileDif << 
995       common.TargMassT2 = common.M0target2 + c << 
996     } else {                          // targe << 
997       common.ProjMassT2 = common.M0projectile2 << 
998       common.TargMassT2 = common.TargetDiffSta << 
999     }                                          << 
1000     common.ProjMassT = std::sqrt( common.Proj << 
1001     common.TargMassT = std::sqrt( common.Targ << 
1002     if ( common.SqrtS < common.ProjMassT + co << 
1003                                               << 
1004     common.PZcms2 = ( sqr( common.S ) + sqr(  << 
1005                       - 2.0 * ( common.S * (  << 
1006                                 + common.Proj << 
1007     if ( common.PZcms2 < 0.0 ) continue;      << 
1008                                               << 
1009     common.PZcms = std::sqrt( common.PZcms2 ) << 
1010     if ( isProjectileDiffraction ) {  // proj << 
1011       common.PMinusMin = std::sqrt( common.Pr << 
1012       common.PMinusMax = common.SqrtS - commo << 
1013       common.PMinusNew = ChooseP( common.PMin << 
1014       common.TMinusNew = common.SqrtS - commo << 
1015       common.Qminus = common.Ptarget.minus()  << 
1016       common.TPlusNew = common.TargMassT2 / c << 
1017       common.Qplus = common.Ptarget.plus() -  << 
1018       common.Qmomentum.setPz( (common.Qplus - << 
1019       common.Qmomentum.setE(  (common.Qplus + << 
1020       loopCondition = ( ( common.Pprojectile  << 
1021                         common.ProjectileDiff << 
1022     } else {                          // targ << 
1023       common.TPlusMin = std::sqrt( common.Tar << 
1024       common.TPlusMax = common.SqrtS - common << 
1025       common.TPlusNew = ChooseP( common.TPlus << 
1026       common.PPlusNew = common.SqrtS - common << 
1027       common.Qplus = common.PPlusNew - common << 
1028       common.PMinusNew = common.ProjMassT2 /  << 
1029       common.Qminus = common.PMinusNew - comm << 
1030       common.Qmomentum.setPz( (common.Qplus - << 
1031       common.Qmomentum.setE(  (common.Qplus + << 
1032       loopCondition =  ( ( common.Ptarget - c << 
1033                          common.TargetDiffSta << 
1034     }                                         << 
1035                                                  367 
1036   } while ( loopCondition );  /* Loop checkin << 368         if( Ksi < 0.333333 ) 
1037           // Repeat the sampling because ther << 369         {TargQ1=TargExchangeQ;}
                                                   >> 370         else if( (0.333333 <= Ksi) && (Ksi < 0.666667)) 
                                                   >> 371         {TargQ2=TargExchangeQ;}
                                                   >> 372         else
                                                   >> 373         {TargQ3=TargExchangeQ;}
                                                   >> 374 
                                                   >> 375        } // End of sampling baryon
                                                   >> 376 
                                                   >> 377        NewProjCode = NewNucleonId(ProjQ1, ProjQ2, ProjQ3); // *****************************
                                                   >> 378 
                                                   >> 379        if((ProjQ1==ProjQ2) && (ProjQ1==ProjQ3)) {NewProjCode +=2; ProjDeltaHasCreated=true;}
                                                   >> 380        else if(projectile->GetDefinition()->GetPDGiIsospin() == 3)// Projectile was Delta
                                                   >> 381        { if(G4UniformRand() > DeltaProbAtQuarkExchange)
                                                   >> 382                                                 {NewProjCode +=2; ProjDeltaHasCreated=true;} 
                                                   >> 383          else                                   {NewProjCode +=0; ProjDeltaHasCreated=false;}
                                                   >> 384        }
                                                   >> 385        else                                                       // Projectile was Nucleon
                                                   >> 386        {
                                                   >> 387         if((G4UniformRand() < DeltaProbAtQuarkExchange) && (SqrtS > DeltaMass+M0target)) 
                                                   >> 388                                                 {NewProjCode +=2; ProjDeltaHasCreated=true;}
                                                   >> 389         else                                    {NewProjCode +=0; ProjDeltaHasCreated=false;}
                                                   >> 390        } 
                                                   >> 391 
                                                   >> 392        NewTargCode = NewNucleonId(TargQ1, TargQ2, TargQ3); // *****************************
                                                   >> 393 
                                                   >> 394        if((TargQ1==TargQ2) && (TargQ1==TargQ3)) {NewTargCode +=2; TargDeltaHasCreated=true;}  
                                                   >> 395        else if(target->GetDefinition()->GetPDGiIsospin() == 3)    // Target was Delta
                                                   >> 396        { if(G4UniformRand() > DeltaProbAtQuarkExchange)
                                                   >> 397                                                 {NewTargCode +=2; TargDeltaHasCreated=true;} 
                                                   >> 398          else                                   {NewTargCode +=0; TargDeltaHasCreated=false;}
                                                   >> 399        }
                                                   >> 400        else                                                       // Target was Nucleon
                                                   >> 401        {
                                                   >> 402         if((G4UniformRand() < DeltaProbAtQuarkExchange) && (SqrtS > M0projectile+DeltaMass)) 
                                                   >> 403                                                 {NewTargCode +=2; TargDeltaHasCreated=true;}
                                                   >> 404         else                                    {NewTargCode +=0; TargDeltaHasCreated=false;}
                                                   >> 405        }         
                                                   >> 406 
                                                   >> 407        if((absProjectilePDGcode == NewProjCode) && (absTargetPDGcode == NewTargCode))
                                                   >> 408        { // Nothing was changed! It is not right!?
                                                   >> 409        }
                                                   >> 410 // Forming baryons --------------------------------------------------
                                                   >> 411 
                                                   >> 412       } // End of if projectile is baryon ---------------------------
                                                   >> 413 
                                                   >> 414 
                                                   >> 415 // If we assume that final state hadrons after the charge exchange will be
                                                   >> 416 // in the ground states, we have to put ----------------------------------
                                                   >> 417 
                                                   >> 418       if(M0projectile < 
                                                   >> 419          (G4ParticleTable::GetParticleTable()->FindParticle(NewProjCode))->GetPDGMass())
                                                   >> 420       {
                                                   >> 421        M0projectile=
                                                   >> 422          (G4ParticleTable::GetParticleTable()->FindParticle(NewProjCode))->GetPDGMass();
                                                   >> 423        M0projectile2 = M0projectile * M0projectile;
                                                   >> 424       }
                                                   >> 425 
                                                   >> 426       if(M0target < 
                                                   >> 427          (G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode))->GetPDGMass())
                                                   >> 428       {
                                                   >> 429        M0target=
                                                   >> 430          (G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode))->GetPDGMass();
                                                   >> 431        M0target2 = M0target * M0target;
                                                   >> 432       }
                                                   >> 433 
                                                   >> 434       PZcms2=(S*S+M0projectile2*M0projectile2+M0target2*M0target2-
                                                   >> 435              2*S*M0projectile2 - 2*S*M0target2 - 2*M0projectile2*M0target2)
                                                   >> 436              /4./S;
                                                   >> 437 
                                                   >> 438       if(PZcms2 < 0) {return false;}  // It can be if energy is not sufficient for Delta
                                                   >> 439 //----------------------------------------------------------
                                                   >> 440       projectile->SetDefinition(
                                                   >> 441                   G4ParticleTable::GetParticleTable()->FindParticle(NewProjCode)); 
                                                   >> 442 
                                                   >> 443       target->SetDefinition(
                                                   >> 444                   G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode)); 
                                                   >> 445 //----------------------------------------------------------
                                                   >> 446       PZcms = std::sqrt(PZcms2);
                                                   >> 447 
                                                   >> 448       Pprojectile.setPz( PZcms);
                                                   >> 449       Pprojectile.setE(std::sqrt(M0projectile2+PZcms2));
                                                   >> 450 
                                                   >> 451       Ptarget.setPz(    -PZcms);
                                                   >> 452       Ptarget.setE(std::sqrt(M0target2+PZcms2));
                                                   >> 453 
                                                   >> 454       {
                                                   >> 455        Pprojectile.transform(toLab);
                                                   >> 456        Ptarget.transform(toLab);
                                                   >> 457 
                                                   >> 458        projectile->SetTimeOfCreation(target->GetTimeOfCreation());
                                                   >> 459        projectile->SetPosition(target->GetPosition());
                                                   >> 460 
                                                   >> 461        projectile->Set4Momentum(Pprojectile);
                                                   >> 462        target->Set4Momentum(Ptarget);
                                                   >> 463 
                                                   >> 464        G4bool Result= theElastic->ElasticScattering (projectile,target,theParameters);
                                                   >> 465 
                                                   >> 466        return Result;
                                                   >> 467       } 
                                                   >> 468      }  // End of charge exchange part ------------------------------
                                                   >> 469 
                                                   >> 470 // ------------------------------------------------------------------
                                                   >> 471      if(ProbOfDiffraction!=0.)
                                                   >> 472      {
                                                   >> 473       ProbProjectileDiffraction/=ProbOfDiffraction;
                                                   >> 474      }
                                                   >> 475      else
                                                   >> 476      {
                                                   >> 477       ProbProjectileDiffraction=0.;
                                                   >> 478      }
                                                   >> 479 
                                                   >> 480      G4double ProjectileDiffStateMinMass2    = ProjectileDiffStateMinMass    *
                                                   >> 481                                                ProjectileDiffStateMinMass;
                                                   >> 482      G4double ProjectileNonDiffStateMinMass2 = ProjectileNonDiffStateMinMass *
                                                   >> 483                                                ProjectileNonDiffStateMinMass;
                                                   >> 484 
                                                   >> 485      G4double TargetDiffStateMinMass2        = TargetDiffStateMinMass        *
                                                   >> 486                                                TargetDiffStateMinMass;
                                                   >> 487      G4double TargetNonDiffStateMinMass2     = TargetNonDiffStateMinMass     *
                                                   >> 488                                                TargetNonDiffStateMinMass;
                                                   >> 489 
                                                   >> 490      G4double Pt2;
                                                   >> 491      G4double ProjMassT2, ProjMassT;
                                                   >> 492      G4double TargMassT2, TargMassT;
                                                   >> 493      G4double PMinusMin, PMinusMax;
                                                   >> 494 //   G4double PPlusMin , PPlusMax;
                                                   >> 495      G4double TPlusMin , TPlusMax;
                                                   >> 496      G4double PMinusNew, PPlusNew, TPlusNew, TMinusNew;
                                                   >> 497 
                                                   >> 498      G4LorentzVector Qmomentum;
                                                   >> 499      G4double Qminus, Qplus;
                                                   >> 500 
                                                   >> 501      G4int whilecount=0;
                                                   >> 502 //   Choose a process ---------------------------
                                                   >> 503 
                                                   >> 504      if(G4UniformRand() < ProbOfDiffraction)
                                                   >> 505        {
                                                   >> 506         if(G4UniformRand() < ProbProjectileDiffraction)
                                                   >> 507         { //-------- projectile diffraction ---------------
                                                   >> 508          do {
                                                   >> 509 //             Generate pt
                                                   >> 510 //             if (whilecount++ >= 500 && (whilecount%100)==0)
                                                   >> 511 //       G4cout << "G4DiffractiveExcitation::ExciteParticipants possibly looping"
                                                   >> 512 //       << ", loop count/ maxPtSquare : "
                                                   >> 513 //             << whilecount << " / " << maxPtSquare << G4endl;
                                                   >> 514              if (whilecount > 1000 )
                                                   >> 515              {
                                                   >> 516               Qmomentum=G4LorentzVector(0.,0.,0.,0.);
                                                   >> 517               target->SetStatus(2);  return false;    //  Ignore this interaction
                                                   >> 518              };
                                                   >> 519 // --------------- Check that the interaction is possible -----------
                                                   >> 520              ProjMassT2=ProjectileDiffStateMinMass2;
                                                   >> 521              ProjMassT =ProjectileDiffStateMinMass;
                                                   >> 522 
                                                   >> 523              TargMassT2=M0target2;
                                                   >> 524              TargMassT =M0target;
                                                   >> 525 
                                                   >> 526              PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2-
                                                   >> 527                      2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
                                                   >> 528                     /4./S;
                                                   >> 529 
                                                   >> 530              if(PZcms2 < 0 ) 
                                                   >> 531              {
                                                   >> 532                target->SetStatus(2);  
                                                   >> 533                return false;
                                                   >> 534              }
                                                   >> 535              maxPtSquare=PZcms2;
                                                   >> 536 
                                                   >> 537              Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0);
                                                   >> 538              Pt2=G4ThreeVector(Qmomentum.vect()).mag2();
                                                   >> 539 
                                                   >> 540              ProjMassT2=ProjectileDiffStateMinMass2+Pt2;
                                                   >> 541              ProjMassT =std::sqrt(ProjMassT2);
                                                   >> 542 
                                                   >> 543              TargMassT2=M0target2+Pt2;
                                                   >> 544              TargMassT =std::sqrt(TargMassT2);
                                                   >> 545 
                                                   >> 546              PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2-
                                                   >> 547                      2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
                                                   >> 548                     /4./S;
                                                   >> 549 
                                                   >> 550              if(PZcms2 < 0 ) continue;
                                                   >> 551              PZcms =std::sqrt(PZcms2);
                                                   >> 552 
                                                   >> 553              PMinusMin=std::sqrt(ProjMassT2+PZcms2)-PZcms;
                                                   >> 554              PMinusMax=SqrtS-TargMassT;
                                                   >> 555 
                                                   >> 556              PMinusNew=ChooseP(PMinusMin, PMinusMax);
                                                   >> 557 // PMinusNew=1./sqrt(1./PMinusMin-G4UniformRand()*(1./PMinusMin-1./PMinusMax));
                                                   >> 558 
                                                   >> 559              TMinusNew=SqrtS-PMinusNew;
                                                   >> 560              Qminus=Ptarget.minus()-TMinusNew;
                                                   >> 561              TPlusNew=TargMassT2/TMinusNew;
                                                   >> 562              Qplus=Ptarget.plus()-TPlusNew;
                                                   >> 563 
                                                   >> 564              Qmomentum.setPz( (Qplus-Qminus)/2 );
                                                   >> 565              Qmomentum.setE(  (Qplus+Qminus)/2 );
                                                   >> 566           } while (
                                                   >> 567 ((Pprojectile+Qmomentum).mag2() <  ProjectileDiffStateMinMass2) ||  //No without excitation
                                                   >> 568 ((Ptarget    -Qmomentum).mag2() <  M0target2                  ));
                                                   >> 569         }
                                                   >> 570         else
                                                   >> 571         { // -------------- Target diffraction ----------------
                                                   >> 572          do {
                                                   >> 573 //             Generate pt
                                                   >> 574 //             if (whilecount++ >= 500 && (whilecount%100)==0)
                                                   >> 575 //       G4cout << "G4DiffractiveExcitation::ExciteParticipants possibly looping"
                                                   >> 576 //       << ", loop count/ maxPtSquare : "
                                                   >> 577 //             << whilecount << " / " << maxPtSquare << G4endl;
                                                   >> 578              if (whilecount > 1000 )
                                                   >> 579              {
                                                   >> 580               Qmomentum=G4LorentzVector(0.,0.,0.,0.);
                                                   >> 581               target->SetStatus(2);  return false;    //  Ignore this interaction
                                                   >> 582              };
                                                   >> 583 // --------------- Check that the interaction is possible -----------
                                                   >> 584              ProjMassT2=M0projectile2;
                                                   >> 585              ProjMassT =M0projectile;
                                                   >> 586 
                                                   >> 587              TargMassT2=TargetDiffStateMinMass2;
                                                   >> 588              TargMassT =TargetDiffStateMinMass;
                                                   >> 589 
                                                   >> 590              PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2-
                                                   >> 591                      2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
                                                   >> 592                     /4./S;
                                                   >> 593 
                                                   >> 594              if(PZcms2 < 0 ) 
                                                   >> 595              {
                                                   >> 596                target->SetStatus(2);  
                                                   >> 597                return false;
                                                   >> 598              }
                                                   >> 599              maxPtSquare=PZcms2;
                                                   >> 600 
                                                   >> 601              Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0);
                                                   >> 602              Pt2=G4ThreeVector(Qmomentum.vect()).mag2();
                                                   >> 603 
                                                   >> 604              ProjMassT2=M0projectile2+Pt2;
                                                   >> 605              ProjMassT =std::sqrt(ProjMassT2);
                                                   >> 606 
                                                   >> 607              TargMassT2=TargetDiffStateMinMass2+Pt2;
                                                   >> 608              TargMassT =std::sqrt(TargMassT2);
                                                   >> 609 
                                                   >> 610              PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2-
                                                   >> 611                      2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
                                                   >> 612                     /4./S;
                                                   >> 613 
                                                   >> 614              if(PZcms2 < 0 ) continue;
                                                   >> 615              PZcms =std::sqrt(PZcms2);
                                                   >> 616 
                                                   >> 617              TPlusMin=std::sqrt(TargMassT2+PZcms2)-PZcms;
                                                   >> 618              TPlusMax=SqrtS-ProjMassT;
                                                   >> 619 
                                                   >> 620              TPlusNew=ChooseP(TPlusMin, TPlusMax);
                                                   >> 621 
                                                   >> 622 //TPlusNew=TPlusMax;
                                                   >> 623 
                                                   >> 624              PPlusNew=SqrtS-TPlusNew;
                                                   >> 625              Qplus=PPlusNew-Pprojectile.plus();
                                                   >> 626              PMinusNew=ProjMassT2/PPlusNew;
                                                   >> 627              Qminus=PMinusNew-Pprojectile.minus();
                                                   >> 628 
                                                   >> 629              Qmomentum.setPz( (Qplus-Qminus)/2 );
                                                   >> 630              Qmomentum.setE(  (Qplus+Qminus)/2 );
                                                   >> 631 
                                                   >> 632           } while (
                                                   >> 633  ((Pprojectile+Qmomentum).mag2() <  M0projectile2          ) ||  //No without excitation
                                                   >> 634  ((Ptarget    -Qmomentum).mag2() <  TargetDiffStateMinMass2));
                                                   >> 635          }
                                                   >> 636         }
                                                   >> 637         else  //----------- Non-diffraction process ------------
                                                   >> 638         {
                                                   >> 639          do {
                                                   >> 640 //             Generate pt
                                                   >> 641 //             if (whilecount++ >= 500 && (whilecount%100)==0)
                                                   >> 642 //       G4cout << "G4DiffractiveExcitation::ExciteParticipants possibly looping"
                                                   >> 643 //       << ", loop count/ maxPtSquare : "
                                                   >> 644 //             << whilecount << " / " << maxPtSquare << G4endl;
                                                   >> 645              if (whilecount > 1000 )
                                                   >> 646              {
                                                   >> 647               Qmomentum=G4LorentzVector(0.,0.,0.,0.);
                                                   >> 648               target->SetStatus(2);  return false;    //  Ignore this interaction
                                                   >> 649              };
                                                   >> 650 // --------------- Check that the interaction is possible -----------
                                                   >> 651              ProjMassT2=ProjectileNonDiffStateMinMass2;
                                                   >> 652              ProjMassT =ProjectileNonDiffStateMinMass;
                                                   >> 653 
                                                   >> 654              TargMassT2=TargetNonDiffStateMinMass2;
                                                   >> 655              TargMassT =TargetNonDiffStateMinMass;
                                                   >> 656 
                                                   >> 657              PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2-
                                                   >> 658                     2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
                                                   >> 659                    /4./S;
                                                   >> 660 
                                                   >> 661              if(PZcms2 < 0 ) 
                                                   >> 662              {
                                                   >> 663                target->SetStatus(2);  
                                                   >> 664                return false;
                                                   >> 665              }
                                                   >> 666              maxPtSquare=PZcms2;
                                                   >> 667              Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0);
                                                   >> 668              Pt2=G4ThreeVector(Qmomentum.vect()).mag2();
                                                   >> 669 
                                                   >> 670              ProjMassT2=ProjectileNonDiffStateMinMass2+Pt2;
                                                   >> 671              ProjMassT =std::sqrt(ProjMassT2);
                                                   >> 672 
                                                   >> 673              TargMassT2=TargetNonDiffStateMinMass2+Pt2;
                                                   >> 674              TargMassT =std::sqrt(TargMassT2);
                                                   >> 675 
                                                   >> 676              PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2-
                                                   >> 677                     2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
                                                   >> 678                    /4./S;
                                                   >> 679 
                                                   >> 680              if(PZcms2 < 0 ) continue;
                                                   >> 681              PZcms =std::sqrt(PZcms2);
                                                   >> 682 
                                                   >> 683              PMinusMin=std::sqrt(ProjMassT2+PZcms2)-PZcms;
                                                   >> 684              PMinusMax=SqrtS-TargMassT;
                                                   >> 685 
                                                   >> 686              PMinusNew=ChooseP(PMinusMin, PMinusMax);
                                                   >> 687 
                                                   >> 688              Qminus=PMinusNew-Pprojectile.minus();
                                                   >> 689 
                                                   >> 690              TPlusMin=std::sqrt(TargMassT2+PZcms2)-PZcms;
                                                   >> 691 //           TPlusMax=SqrtS-PMinusNew;                      
                                                   >> 692              TPlusMax=SqrtS-ProjMassT;      
                                                   >> 693 
                                                   >> 694              TPlusNew=ChooseP(TPlusMin, TPlusMax);
                                                   >> 695 
                                                   >> 696              Qplus=-(TPlusNew-Ptarget.plus());
                                                   >> 697 
                                                   >> 698              Qmomentum.setPz( (Qplus-Qminus)/2 );
                                                   >> 699              Qmomentum.setE(  (Qplus+Qminus)/2 );
                                                   >> 700 
                                                   >> 701        } while (
                                                   >> 702  ((Pprojectile+Qmomentum).mag2() <  ProjectileNonDiffStateMinMass2) || //No double Diffraction
                                                   >> 703  ((Ptarget    -Qmomentum).mag2() <  TargetNonDiffStateMinMass2    ));
                                                   >> 704          }
                                                   >> 705 
                                                   >> 706      Pprojectile += Qmomentum;
                                                   >> 707      Ptarget     -= Qmomentum;
                                                   >> 708 
                                                   >> 709 // Transform back and update SplitableHadron Participant.
                                                   >> 710      Pprojectile.transform(toLab);
                                                   >> 711      Ptarget.transform(toLab);
                                                   >> 712 
                                                   >> 713 // Calculation of the creation time ---------------------
                                                   >> 714       projectile->SetTimeOfCreation(target->GetTimeOfCreation());
                                                   >> 715       projectile->SetPosition(target->GetPosition());
                                                   >> 716 // Creation time and position of target nucleon were determined at
                                                   >> 717 // ReggeonCascade() of G4FTFModel
                                                   >> 718 // ------------------------------------------------------
1038                                                  719 
1039   if ( isProjectileDiffraction ) {  // projec << 720      projectile->Set4Momentum(Pprojectile);
1040     projectile->SetStatus( 0 );               << 721      target->Set4Momentum(Ptarget);
1041     if ( projectile->GetStatus() == 2 ) proje << 
1042     if ( target->GetStatus() == 1  &&  target << 
1043   } else {                          // target << 
1044     target->SetStatus( 0 );                   << 
1045   }                                           << 
1046                                                  722 
1047   return true;                                << 723            projectile->IncrementCollisionCount(1);
1048 }                                             << 724            target->IncrementCollisionCount(1);
1049                                               << 
1050 //------------------------------------------- << 
1051                                               << 
1052 G4bool G4DiffractiveExcitation::              << 
1053 ExciteParticipants_doNonDiffraction( G4VSplit << 
1054                                      G4VSplit << 
1055                                      G4FTFPar << 
1056                                      G4Diffra << 
1057   // Third of the three utility methods used  << 
1058   // it does the sampling for the non-diffrac << 
1059                                               << 
1060   #ifdef debugFTFexcitation                   << 
1061   G4cout << "Non-diffraction process" << G4en << 
1062   #endif                                      << 
1063                                               << 
1064   // Check whether the interaction is possibl << 
1065   common.ProjMassT2 = common.ProjectileNonDif << 
1066   common.ProjMassT  = common.ProjectileNonDif << 
1067   common.TargMassT2 = common.TargetNonDiffSta << 
1068   common.TargMassT  = common.TargetNonDiffSta << 
1069   if ( common.SqrtS < common.ProjMassT + comm << 
1070                                               << 
1071   common.PZcms2 = ( sqr( common.S ) + sqr( co << 
1072                   - 2.0 * ( common.S * ( comm << 
1073                         + common.ProjMassT2 * << 
1074   common.maxPtSquare = common.PZcms2;         << 
1075                                               << 
1076   G4int whilecount = 0;                       << 
1077   do {  // Generate pt and masses             << 
1078                                               << 
1079     whilecount++;                             << 
1080     if ( whilecount > 1000 ) {                << 
1081       common.Qmomentum = G4LorentzVector( 0.0 << 
1082       return false;  // Ignore this interacti << 
1083     };                                        << 
1084                                               << 
1085     common.Qmomentum = G4LorentzVector( Gauss << 
1086                                               << 
1087     common.Pt2 = G4ThreeVector( common.Qmomen << 
1088     common.ProjMassT2 = common.ProjectileNonD << 
1089     common.ProjMassT  = std::sqrt( common.Pro << 
1090     common.TargMassT2 = common.TargetNonDiffS << 
1091     common.TargMassT  = std::sqrt( common.Tar << 
1092     if ( common.SqrtS < common.ProjMassT + co << 
1093                                               << 
1094     common.PZcms2 =( sqr( common.S ) + sqr( c << 
1095                      - 2.0 * ( common.S * ( c << 
1096                                + common.ProjM << 
1097     if ( common.PZcms2 < 0.0 ) continue;      << 
1098                                               << 
1099     common.PZcms = std::sqrt( common.PZcms2 ) << 
1100     common.PMinusMin = std::sqrt( common.Proj << 
1101     common.PMinusMax = common.SqrtS - common. << 
1102     common.TPlusMin = std::sqrt( common.TargM << 
1103     common.TPlusMax = common.SqrtS - common.P << 
1104                                               << 
1105     if ( G4UniformRand() <= 0.5 ) {  // Rando << 
1106       if ( G4UniformRand() < theParameters->G << 
1107         common.PMinusNew = ChooseP( common.PM << 
1108       } else {                                << 
1109         common.PMinusNew = ( common.PMinusMax << 
1110       }                                       << 
1111       if ( G4UniformRand() < theParameters->G << 
1112         common.TPlusNew = ChooseP( common.TPl << 
1113       } else {                                << 
1114         common.TPlusNew = ( common.TPlusMax - << 
1115       }                                       << 
1116     } else {                                  << 
1117       if ( G4UniformRand() < theParameters->G << 
1118         common.TPlusNew = ChooseP( common.TPl << 
1119       } else {                                << 
1120         common.TPlusNew = ( common.TPlusMax - << 
1121       }                                       << 
1122       if ( G4UniformRand() < theParameters->G << 
1123         common.PMinusNew = ChooseP( common.PM << 
1124       } else {                                << 
1125         common.PMinusNew = ( common.PMinusMax << 
1126       }                                       << 
1127     }                                         << 
1128                                                  725 
1129     common.Qminus = common.PMinusNew - common << 726      return true;
1130     common.Qplus = -( common.TPlusNew - commo << 
1131     common.Qmomentum.setPz( (common.Qplus - c << 
1132     common.Qmomentum.setE(  (common.Qplus + c << 
1133     #ifdef debugFTFexcitation                 << 
1134     G4cout <<"Sampled: Mpr, MdifPr, Mtr, Mdif << 
1135            << ( common.Pprojectile + common.Q << 
1136            << common.ProjectileNonDiffStateMi << 
1137            << ( common.Ptarget - common.Qmome << 
1138            << common.TargetNonDiffStateMinMas << 
1139     #endif                                    << 
1140                                               << 
1141   } while ( ( common.Pprojectile + common.Qmo << 
1142             common.ProjectileNonDiffStateMinM << 
1143             ( common.Ptarget - common.Qmoment << 
1144             common.TargetNonDiffStateMinMass2 << 
1145                                               << 
1146   projectile->SetStatus( 0 );                 << 
1147   target->SetStatus( 0 );                     << 
1148   return true;                                << 
1149 }                                                727 }
1150                                                  728 
                                                   >> 729 // ---------------------------------------------------------------------
                                                   >> 730 void G4DiffractiveExcitation::CreateStrings(G4VSplitableHadron * hadron, 
                                                   >> 731                                             G4bool isProjectile,
                                                   >> 732                                             G4ExcitedString * &FirstString, 
                                                   >> 733                                             G4ExcitedString * &SecondString,
                                                   >> 734                                             G4FTFParameters *theParameters) const
                                                   >> 735 {
                                                   >> 736   hadron->SplitUp();
                                                   >> 737   G4Parton *start= hadron->GetNextParton();
                                                   >> 738   if ( start==NULL)
                                                   >> 739   { G4cout << " G4FTFModel::String() Error:No start parton found"<< G4endl;
                                                   >> 740           FirstString=0; SecondString=0;
                                                   >> 741     return;
                                                   >> 742   }
                                                   >> 743   G4Parton *end  = hadron->GetNextParton();
                                                   >> 744   if ( end==NULL)
                                                   >> 745   { G4cout << " G4FTFModel::String() Error:No end parton found"<< G4endl;
                                                   >> 746           FirstString=0; SecondString=0;
                                                   >> 747     return;
                                                   >> 748   }
                                                   >> 749 
                                                   >> 750         G4LorentzVector Phadron=hadron->Get4Momentum();
                                                   >> 751 
                                                   >> 752         G4LorentzVector Pstart(0.,0.,0.,0.);
                                                   >> 753         G4LorentzVector Pend(0.,0.,0.,0.);
                                                   >> 754         G4LorentzVector Pkink(0.,0.,0.,0.);
                                                   >> 755         G4LorentzVector PkinkQ1(0.,0.,0.,0.);
                                                   >> 756         G4LorentzVector PkinkQ2(0.,0.,0.,0.);
                                                   >> 757 
                                                   >> 758         G4int PDGcode_startQ = std::abs(start->GetDefinition()->GetPDGEncoding());
                                                   >> 759         G4int PDGcode_endQ   = std::abs(  end->GetDefinition()->GetPDGEncoding());
                                                   >> 760 
                                                   >> 761 //--------------------------------------------------------------------------------        
                                                   >> 762         G4double Wmin(0.);
                                                   >> 763         if(isProjectile)
                                                   >> 764         {
                                                   >> 765           Wmin=theParameters->GetProjMinDiffMass();
                                                   >> 766         } else
                                                   >> 767         {
                                                   >> 768           Wmin=theParameters->GetTarMinDiffMass();
                                                   >> 769         } // end of if(isProjectile)
                                                   >> 770 
                                                   >> 771         G4double W = hadron->Get4Momentum().mag();
                                                   >> 772         G4double W2=W*W;
                                                   >> 773 
                                                   >> 774         G4double Pt(0.), x1(0.), x2(0.), x3(0.);
                                                   >> 775         G4bool Kink=false;
                                                   >> 776 
                                                   >> 777         if(W > Wmin)
                                                   >> 778         {                                        // Kink is possible
                                                   >> 779           G4double Pt2kink=theParameters->GetPt2Kink();
                                                   >> 780           Pt = std::sqrt(Pt2kink*(std::pow(W2/16./Pt2kink+1.,G4UniformRand()) - 1.));
                                                   >> 781 
                                                   >> 782           if(Pt > 500.*MeV)
                                                   >> 783           {
                                                   >> 784              G4double Ymax = std::log(W/2./Pt + std::sqrt(W2/4./Pt/Pt - 1.));
                                                   >> 785              G4double Y=Ymax*(1.- 2.*G4UniformRand());
                                                   >> 786 
                                                   >> 787              x1=1.-Pt/W*std::exp( Y);
                                                   >> 788              x3=1.-Pt/W*std::exp(-Y);
                                                   >> 789              x2=2.-x1-x3;
                                                   >> 790 
                                                   >> 791              G4double Mass_startQ = 650.*MeV;
                                                   >> 792              if(PDGcode_startQ <  3) Mass_startQ =  325.*MeV;
                                                   >> 793              if(PDGcode_startQ == 3) Mass_startQ =  500.*MeV;
                                                   >> 794              if(PDGcode_startQ == 4) Mass_startQ = 1600.*MeV;
                                                   >> 795 
                                                   >> 796              G4double Mass_endQ = 650.*MeV;
                                                   >> 797              if(PDGcode_endQ <  3) Mass_endQ =  325.*MeV;
                                                   >> 798              if(PDGcode_endQ == 3) Mass_endQ =  500.*MeV;
                                                   >> 799              if(PDGcode_endQ == 4) Mass_endQ = 1600.*MeV;
                                                   >> 800 
                                                   >> 801              G4double P2_1=W2*x1*x1/4.-Mass_endQ  *Mass_endQ;
                                                   >> 802              G4double P2_3=W2*x3*x3/4.-Mass_startQ*Mass_startQ;
                                                   >> 803      
                                                   >> 804              G4double P2_2=sqr((2.-x1-x3)*W/2.);
                                                   >> 805 
                                                   >> 806              if((P2_1 <= 0.) || (P2_3 <= 0.))
                                                   >> 807              { Kink=false;}
                                                   >> 808              else
                                                   >> 809              {
                                                   >> 810                G4double P_1=std::sqrt(P2_1);
                                                   >> 811                G4double P_2=std::sqrt(P2_2);
                                                   >> 812                G4double P_3=std::sqrt(P2_3);
                                                   >> 813 
                                                   >> 814                G4double CosT12=(P2_3-P2_1-P2_2)/(2.*P_1*P_2);
                                                   >> 815                G4double CosT13=(P2_2-P2_1-P2_3)/(2.*P_1*P_3);
                                                   >> 816 //             Pt=P_2*std::sqrt(1.-CosT12*CosT12);  // because system was rotated 11.12.09
                                                   >> 817 
                                                   >> 818                if((std::abs(CosT12) >1.) || (std::abs(CosT13) > 1.)) 
                                                   >> 819                { Kink=false;}
                                                   >> 820                else
                                                   >> 821                { 
                                                   >> 822                  Kink=true; 
                                                   >> 823                  Pt=P_2*std::sqrt(1.-CosT12*CosT12);  // because system was rotated 11.12.09
                                                   >> 824                  Pstart.setPx(-Pt); Pstart.setPy(0.); Pstart.setPz(P_3*CosT13); 
                                                   >> 825                  Pend.setPx(   0.); Pend.setPy(  0.); Pend.setPz(          P_1); 
                                                   >> 826                  Pkink.setPx(  Pt); Pkink.setPy( 0.); Pkink.setPz(  P_2*CosT12);
                                                   >> 827                  Pstart.setE(x3*W/2.);                
                                                   >> 828                  Pkink.setE(Pkink.vect().mag());
                                                   >> 829                  Pend.setE(x1*W/2.);
                                                   >> 830 
                                                   >> 831                  G4double XkQ=GetQuarkFractionOfKink(0.,1.);
                                                   >> 832                  if(Pkink.getZ() > 0.) 
                                                   >> 833                  {
                                                   >> 834                    if(XkQ > 0.5) {PkinkQ1=XkQ*Pkink;} else {PkinkQ1=(1.-XkQ)*Pkink;}
                                                   >> 835                  } else {
                                                   >> 836                    if(XkQ > 0.5) {PkinkQ1=(1.-XkQ)*Pkink;} else {PkinkQ1=XkQ*Pkink;}
                                                   >> 837                  }
                                                   >> 838 
                                                   >> 839                  PkinkQ2=Pkink - PkinkQ1;
                                                   >> 840 //------------------------- Minimizing Pt1^2+Pt3^2 ------------------------------
                                                   >> 841 
                                                   >> 842                  G4double Cos2Psi=(sqr(x1) -sqr(x3)+2.*sqr(x3*CosT13))/
                                                   >> 843                           std::sqrt(sqr(sqr(x1)-sqr(x3)) + sqr(2.*x1*x3*CosT13));
                                                   >> 844                  G4double Psi=std::acos(Cos2Psi);
                                                   >> 845 
                                                   >> 846 G4LorentzRotation Rotate;
                                                   >> 847 if(isProjectile) {Rotate.rotateY(Psi);}
                                                   >> 848 else             {Rotate.rotateY(pi-Psi);}                   
                                                   >> 849 Rotate.rotateZ(twopi*G4UniformRand());
                                                   >> 850 
                                                   >> 851 Pstart*=Rotate;
                                                   >> 852 Pkink*=Rotate;
                                                   >> 853 PkinkQ1*=Rotate;
                                                   >> 854 PkinkQ2*=Rotate;
                                                   >> 855 Pend*=Rotate;
                                                   >> 856 
                                                   >> 857                }
                                                   >> 858              }      // end of if((P2_1 < 0.) || (P2_3 <0.))
                                                   >> 859           }         // end of if(Pt > 500.*MeV)
                                                   >> 860         }           // end of if(W > Wmin) Check for a kink
                                                   >> 861 
                                                   >> 862 //--------------------------------------------------------------------------------
                                                   >> 863 
                                                   >> 864         if(Kink)
                                                   >> 865         {                                        // Kink is possible
                                                   >> 866           std::vector<G4double> QuarkProbabilitiesAtGluonSplitUp =
                                                   >> 867               theParameters->GetQuarkProbabilitiesAtGluonSplitUp();
                                                   >> 868 
                                                   >> 869           G4int QuarkInGluon(1); G4double Ksi=G4UniformRand();
                                                   >> 870           for(unsigned int Iq=0; Iq <3; Iq++)
                                                   >> 871           {
                                                   >> 872 
                                                   >> 873 if(Ksi > QuarkProbabilitiesAtGluonSplitUp[Iq]) QuarkInGluon++;}
                                                   >> 874 
                                                   >> 875           G4Parton * Gquark = new G4Parton(QuarkInGluon);
                                                   >> 876           G4Parton * Ganti_quark = new G4Parton(-QuarkInGluon);
                                                   >> 877 
                                                   >> 878 //-------------------------------------------------------------------------------
                                                   >> 879           G4LorentzRotation toCMS(-1*Phadron.boostVector());
                                                   >> 880 
                                                   >> 881           G4LorentzRotation toLab(toCMS.inverse());
                                                   >> 882 
                                                   >> 883           Pstart.transform(toLab);  start->Set4Momentum(Pstart);
                                                   >> 884           PkinkQ1.transform(toLab);
                                                   >> 885           PkinkQ2.transform(toLab);
                                                   >> 886           Pend.transform(toLab);    end->Set4Momentum(Pend);
                                                   >> 887 
                                                   >> 888           G4int absPDGcode=std::abs(hadron->GetDefinition()->GetPDGEncoding());
                                                   >> 889 
                                                   >> 890           if(absPDGcode < 1000)
                                                   >> 891           {                                // meson
                                                   >> 892       if ( isProjectile )
                                                   >> 893       {                              // Projectile
                                                   >> 894               if(end->GetDefinition()->GetPDGEncoding() > 0 )  // A quark on the end
                                                   >> 895               {                            // Quark on the end
                                                   >> 896                 FirstString = new G4ExcitedString(end   ,Ganti_quark, +1);
                                                   >> 897                 SecondString= new G4ExcitedString(Gquark,start      ,+1);
                                                   >> 898                 Ganti_quark->Set4Momentum(PkinkQ1);
                                                   >> 899                 Gquark->Set4Momentum(PkinkQ2);
                                                   >> 900 
                                                   >> 901               } else
                                                   >> 902               {                            // Anti_Quark on the end
                                                   >> 903                 FirstString = new G4ExcitedString(end        ,Gquark, +1);
                                                   >> 904                 SecondString= new G4ExcitedString(Ganti_quark,start ,+1);
                                                   >> 905                 Gquark->Set4Momentum(PkinkQ1);
                                                   >> 906                 Ganti_quark->Set4Momentum(PkinkQ2);
                                                   >> 907 
                                                   >> 908               }   // end of if(end->GetPDGcode() > 0)
                                                   >> 909             } else {                      // Target
                                                   >> 910               if(end->GetDefinition()->GetPDGEncoding() > 0 )  // A quark on the end
                                                   >> 911               {                           // Quark on the end
                                                   >> 912                 FirstString = new G4ExcitedString(Ganti_quark,end   ,-1);
                                                   >> 913                 SecondString= new G4ExcitedString(start      ,Gquark,-1);
                                                   >> 914                 Ganti_quark->Set4Momentum(PkinkQ2);
                                                   >> 915                 Gquark->Set4Momentum(PkinkQ1);
                                                   >> 916 
                                                   >> 917               } else
                                                   >> 918               {                            // Anti_Quark on the end
                                                   >> 919                 FirstString = new G4ExcitedString(Gquark,end        ,-1);
                                                   >> 920                 SecondString= new G4ExcitedString(start ,Ganti_quark,-1);
                                                   >> 921                 Gquark->Set4Momentum(PkinkQ2);
                                                   >> 922                 Ganti_quark->Set4Momentum(PkinkQ1);
                                                   >> 923 
                                                   >> 924               }   // end of if(end->GetPDGcode() > 0)
                                                   >> 925       }     // end of if ( isProjectile )
                                                   >> 926           } else  // if(absPDGCode < 1000)
                                                   >> 927           {                             // Baryon/AntiBaryon
                                                   >> 928       if ( isProjectile )
                                                   >> 929       {                              // Projectile
                                                   >> 930               if((end->GetDefinition()->GetParticleType() == "diquarks") &&
                                                   >> 931                  (end->GetDefinition()->GetPDGEncoding() > 0           )   ) 
                                                   >> 932               {                            // DiQuark on the end
                                                   >> 933                 FirstString = new G4ExcitedString(end        ,Gquark, +1);
                                                   >> 934                 SecondString= new G4ExcitedString(Ganti_quark,start ,+1);
                                                   >> 935                 Gquark->Set4Momentum(PkinkQ1);
                                                   >> 936                 Ganti_quark->Set4Momentum(PkinkQ2);
                                                   >> 937 
                                                   >> 938               } else
                                                   >> 939               {                            // Anti_DiQuark on the end or quark
                                                   >> 940                 FirstString = new G4ExcitedString(end   ,Ganti_quark, +1);
                                                   >> 941                 SecondString= new G4ExcitedString(Gquark,start      ,+1);
                                                   >> 942                 Ganti_quark->Set4Momentum(PkinkQ1);
                                                   >> 943                 Gquark->Set4Momentum(PkinkQ2);
                                                   >> 944 
                                                   >> 945               }   // end of if(end->GetPDGcode() > 0)
                                                   >> 946             } else {                      // Target
                                                   >> 947 
                                                   >> 948               if((end->GetDefinition()->GetParticleType() == "diquarks") &&
                                                   >> 949                  (end->GetDefinition()->GetPDGEncoding() > 0           )   ) 
                                                   >> 950               {                            // DiQuark on the end
                                                   >> 951                 FirstString = new G4ExcitedString(Gquark,end        ,-1);
                                                   >> 952 
                                                   >> 953                 SecondString= new G4ExcitedString(start ,Ganti_quark,-1);
                                                   >> 954                 Gquark->Set4Momentum(PkinkQ1);
                                                   >> 955                 Ganti_quark->Set4Momentum(PkinkQ2);
                                                   >> 956 
                                                   >> 957               } else
                                                   >> 958               {                            // Anti_DiQuark on the end or Q
                                                   >> 959                 FirstString = new G4ExcitedString(Ganti_quark,end   ,-1);
                                                   >> 960                 SecondString= new G4ExcitedString(start      ,Gquark,-1);
                                                   >> 961                 Gquark->Set4Momentum(PkinkQ2);
                                                   >> 962                 Ganti_quark->Set4Momentum(PkinkQ1);
                                                   >> 963 
                                                   >> 964               }   // end of if(end->GetPDGcode() > 0)
                                                   >> 965       }     // end of if ( isProjectile )
                                                   >> 966           }  // end of if(absPDGcode < 1000)
                                                   >> 967 
                                                   >> 968     FirstString->SetTimeOfCreation(hadron->GetTimeOfCreation());
                                                   >> 969     FirstString->SetPosition(hadron->GetPosition());
                                                   >> 970 
                                                   >> 971     SecondString->SetTimeOfCreation(hadron->GetTimeOfCreation());
                                                   >> 972     SecondString->SetPosition(hadron->GetPosition());
                                                   >> 973 
                                                   >> 974 // -------------------------------------------------------------------------  
                                                   >> 975         } else                                   // End of kink is possible
                                                   >> 976         {                                        // Kink is impossible
                                                   >> 977     if ( isProjectile )
                                                   >> 978     {
                                                   >> 979     FirstString= new G4ExcitedString(end,start, +1);
                                                   >> 980     } else {
                                                   >> 981     FirstString= new G4ExcitedString(start,end, -1);
                                                   >> 982     }
1151                                                  983 
1152 //=========================================== << 984           SecondString=0;
1153                                                  985 
1154 void G4DiffractiveExcitation::CreateStrings(  << 986     FirstString->SetTimeOfCreation(hadron->GetTimeOfCreation());
1155                                               << 987     FirstString->SetPosition(hadron->GetPosition());
1156                                               << 
1157                                               << 
1158                                               << 
1159                                               << 
1160   //G4cout << "Create Strings SplitUp " << ha << 
1161   //       << "Defin " << hadron->GetDefiniti << 
1162   //       << "Defin " << hadron->GetDefiniti << 
1163                                               << 
1164   G4bool HadronIsString = hadron->IsSplit();  << 
1165   if( ! HadronIsString ) hadron->SplitUp();   << 
1166                                               << 
1167   G4Parton* start = hadron->GetNextParton();  << 
1168   if ( start == nullptr ) {                   << 
1169     G4cout << " G4FTFModel::String() Error: N << 
1170     FirstString = 0; SecondString = 0;        << 
1171     return;                                   << 
1172   }                                           << 
1173                                               << 
1174   G4Parton* end = hadron->GetNextParton();    << 
1175   if ( end == nullptr ) {                     << 
1176     G4cout << " G4FTFModel::String() Error: N << 
1177     FirstString = 0; SecondString = 0;        << 
1178     return;                                   << 
1179   }                                           << 
1180                                               << 
1181   //G4cout << start << " " << start->GetPDGco << 
1182   //       << G4endl                          << 
1183   //       << "Create string " << start->GetP << 
1184                                               << 
1185   if ( HadronIsString ) {                     << 
1186     if ( isProjectile ) {                     << 
1187       FirstString = new G4ExcitedString( end, << 
1188     } else {                                  << 
1189       FirstString = new G4ExcitedString( end, << 
1190     }                                         << 
1191     FirstString->SetTimeOfCreation( hadron->G << 
1192     FirstString->SetPosition( hadron->GetPosi << 
1193     SecondString = 0;                         << 
1194     return;                                   << 
1195   }                                           << 
1196                                               << 
1197   G4LorentzVector Phadron = hadron->Get4Momen << 
1198   //G4cout << "String mom " << Phadron << G4e << 
1199   G4LorentzVector Pstart( 0.0, 0.0, 0.0, 0.0  << 
1200   G4LorentzVector Pend( 0.0, 0.0, 0.0, 0.0 ); << 
1201   G4LorentzVector Pkink( 0.0, 0.0, 0.0, 0.0 ) << 
1202   G4LorentzVector PkinkQ1( 0.0, 0.0, 0.0, 0.0 << 
1203   G4LorentzVector PkinkQ2( 0.0, 0.0, 0.0, 0.0 << 
1204                                               << 
1205   G4int PDGcode_startQ = std::abs( start->Get << 
1206   G4int PDGcode_endQ   = std::abs(   end->Get << 
1207   //G4cout << "PDGcode_startQ " << PDGcode_st << 
1208                                               << 
1209   G4double Wmin( 0.0 );                       << 
1210   if ( isProjectile ) {                       << 
1211     Wmin = theParameters->GetProjMinDiffMass( << 
1212   } else {                                    << 
1213     Wmin = theParameters->GetTarMinDiffMass() << 
1214   }                                           << 
1215                                               << 
1216   G4double W = hadron->Get4Momentum().mag();  << 
1217   G4double W2 = W*W;                          << 
1218   G4double Pt( 0.0 ), x1( 0.0 ), x3( 0.0 );   << 
1219   G4bool Kink = false;                        << 
1220                                               << 
1221   if ( ! ( ( start->GetDefinition()->GetParti << 
1222                end->GetDefinition()->GetParti << 
1223            ( start->GetDefinition()->GetParti << 
1224                end->GetDefinition()->GetParti << 
1225     // Kinky strings are allowed only for qq- << 
1226     // Kinky strings are impossible for other << 
1227     // according to the analysis of Pbar P in << 
1228                                               << 
1229     if ( W > Wmin ) {  // Kink is possible    << 
1230       if ( hadron->GetStatus() == 0 ) {       << 
1231         G4double Pt2kink = theParameters->Get << 
1232         if ( Pt2kink ) {                      << 
1233           Pt = std::sqrt( Pt2kink * ( G4Pow:: << 
1234         } else {                              << 
1235           Pt = 0.0;                           << 
1236         }                                     << 
1237       } else {                                << 
1238         Pt = 0.0;                             << 
1239       }                                       << 
1240                                                  988 
1241       if ( Pt > 500.0*MeV ) {                 << 989 // momenta of string ends
1242         G4double Ymax = G4Log( W/2.0/Pt + std << 990 //
1243         G4double Y = Ymax*( 1.0 - 2.0*G4Unifo << 991           G4double Momentum=hadron->Get4Momentum().vect().mag();
1244         x1 = 1.0 - Pt/W * G4Exp( Y );         << 992           G4double  Plus=hadron->Get4Momentum().e() + Momentum;
1245         x3 = 1.0 - Pt/W * G4Exp(-Y );         << 993           G4double Minus=hadron->Get4Momentum().e() - Momentum;
1246         //x2 = 2.0 - x1 - x3;                 << 994 
1247                                               << 995           G4ThreeVector tmp;
1248         G4double Mass_startQ = 650.0*MeV;     << 996           if(Momentum > 0.) 
1249         if ( PDGcode_startQ <  3 ) Mass_start << 997           {
1250         if ( PDGcode_startQ == 3 ) Mass_start << 998            tmp.set(hadron->Get4Momentum().px(),
1251         if ( PDGcode_startQ == 4 ) Mass_start << 999                    hadron->Get4Momentum().py(),
1252         if ( PDGcode_startQ == 5 ) Mass_start << 1000                    hadron->Get4Momentum().pz());
1253         G4double Mass_endQ = 650.0*MeV;       << 1001            tmp/=Momentum;
1254         if ( PDGcode_endQ <  3 ) Mass_endQ =  << 1002           }
1255         if ( PDGcode_endQ == 3 ) Mass_endQ =  << 1003           else
1256         if ( PDGcode_endQ == 4 ) Mass_endQ =  << 1004           {
1257         if ( PDGcode_endQ == 5 ) Mass_endQ =  << 1005            tmp.set(0.,0.,1.);
1258                                               << 
1259         G4double P2_1 = W2*x1*x1/4.0 - Mass_e << 
1260         G4double P2_3 = W2*x3*x3/4.0 - Mass_s << 
1261         G4double P2_2 = sqr( (2.0 - x1 - x3)* << 
1262         if ( P2_1 <= 0.0  ||  P2_3 <= 0.0 ) { << 
1263           Kink = false;                       << 
1264         } else {                              << 
1265           G4double P_1 = std::sqrt( P2_1 );   << 
1266           G4double P_2 = std::sqrt( P2_2 );   << 
1267           G4double P_3 = std::sqrt( P2_3 );   << 
1268           G4double CosT12 = ( P2_3 - P2_1 - P << 
1269           G4double CosT13 = ( P2_2 - P2_1 - P << 
1270                                               << 
1271           if ( std::abs( CosT12 ) > 1.0  ||   << 
1272             Kink = false;                     << 
1273           } else {                            << 
1274             Kink = true;                      << 
1275             Pt = P_2 * std::sqrt( 1.0 - CosT1 << 
1276             Pstart.setPx( -Pt ); Pstart.setPy << 
1277             Pend.setPx(   0.0 ); Pend.setPy(  << 
1278             Pkink.setPx(   Pt ); Pkink.setPy( << 
1279             Pstart.setE( x3*W/2.0 );          << 
1280             Pkink.setE( Pkink.vect().mag() ); << 
1281             Pend.setE( x1*W/2.0 );            << 
1282                                               << 
1283             G4double XkQ = GetQuarkFractionOf << 
1284             if ( Pkink.getZ() > 0.0 ) {       << 
1285               if ( XkQ > 0.5 ) {              << 
1286                 PkinkQ1 = XkQ*Pkink;          << 
1287               } else {                        << 
1288                 PkinkQ1 = (1.0 - XkQ)*Pkink;  << 
1289               }                               << 
1290             } else {                          << 
1291               if ( XkQ > 0.5 ) {              << 
1292                 PkinkQ1 = (1.0 - XkQ)*Pkink;  << 
1293               } else {                        << 
1294                 PkinkQ1 = XkQ*Pkink;          << 
1295               }                               << 
1296             }                                 << 
1297                                               << 
1298             PkinkQ2 = Pkink - PkinkQ1;        << 
1299             // Minimizing Pt1^2+Pt3^2         << 
1300             G4double Cos2Psi = ( sqr(x1) - sq << 
1301                                std::sqrt( sqr << 
1302             G4double Psi = std::acos( Cos2Psi << 
1303                                               << 
1304             G4LorentzRotation Rotate;         << 
1305             if ( isProjectile ) {             << 
1306               Rotate.rotateY( Psi );          << 
1307             } else {                          << 
1308               Rotate.rotateY( pi + Psi );     << 
1309             }                                 << 
1310             Rotate.rotateZ( twopi * G4Uniform << 
1311             Pstart *= Rotate;                 << 
1312             Pkink *= Rotate;                  << 
1313             PkinkQ1 *= Rotate;                << 
1314             PkinkQ2 *= Rotate;                << 
1315             Pend *= Rotate;                   << 
1316           }                                      1006           }
1317         }  // End of if ( P2_1 <= 0.0  ||  P2 << 
1318       }  // End of if ( Pt > 500.0*MeV )      << 
1319     } // End of if ( W > Wmin ) : check for a << 
1320   }  // end of qq-q string selection          << 
1321                                               << 
1322   if ( Kink ) {  // Kink is possible          << 
1323                                               << 
1324     //G4cout << "Kink is sampled!" << G4endl; << 
1325     std::vector< G4double > QuarkProbabilitie << 
1326         theParameters->GetQuarkProbabilitiesA << 
1327                                               << 
1328     G4int QuarkInGluon( 1 ); G4double Ksi = G << 
1329     for ( unsigned int Iq = 0; Iq < 3; Iq++ ) << 
1330       //G4cout << "Iq " << Iq << G4endl;      << 
1331       if ( Ksi > QuarkProbabilitiesAtGluonSpl << 
1332     }                                         << 
1333     //G4cout << "Last Iq " << QuarkInGluon << << 
1334     G4Parton* Gquark = new G4Parton( QuarkInG << 
1335     G4Parton* Ganti_quark = new G4Parton( -Qu << 
1336     //G4cout << "Lorentz " << G4endl;         << 
1337                                               << 
1338     G4LorentzRotation toCMS( -1 * Phadron.boo << 
1339     G4LorentzRotation toLab( toCMS.inverse()  << 
1340     //G4cout << "Pstart " << Pstart << G4endl << 
1341     //G4cout << "Pend   " << Pend << G4endl;  << 
1342     //G4cout << "Kink1  " <<PkinkQ1<<G4endl;  << 
1343     //G4cout << "Kink2  " <<PkinkQ2<<G4endl;  << 
1344     //G4cout << "Pstart " << Pstart << G4endl << 
1345                                               << 
1346     Pstart.transform( toLab );  start->Set4Mo << 
1347     PkinkQ1.transform( toLab );               << 
1348     PkinkQ2.transform( toLab );               << 
1349     Pend.transform( toLab );    end->Set4Mome << 
1350     //G4cout << "Pstart " << Pstart << G4endl << 
1351     //G4cout << "Pend   " << Pend << G4endl;  << 
1352     //G4cout << "Defin " << hadron->GetDefini << 
1353     //G4cout << "Defin " << hadron->GetDefini << 
1354                                               << 
1355     //G4int absPDGcode = std::abs( hadron->Ge << 
1356     G4int absPDGcode = 1500;                  << 
1357     if ( start->GetDefinition()->GetParticleS << 
1358          end->GetDefinition()->GetParticleSub << 
1359       absPDGcode = 110;                       << 
1360     }                                         << 
1361     //G4cout << "absPDGcode " << absPDGcode < << 
1362                                               << 
1363     if ( absPDGcode < 1000 ) {  // meson      << 
1364       if ( isProjectile ) { // Projectile     << 
1365         if ( end->GetDefinition()->GetPDGEnco << 
1366           FirstString  = new G4ExcitedString( << 
1367           SecondString = new G4ExcitedString( << 
1368           Ganti_quark->Set4Momentum( PkinkQ1  << 
1369           Gquark->Set4Momentum( PkinkQ2 );    << 
1370         } else {  // Anti_Quark on the end    << 
1371           FirstString  = new G4ExcitedString( << 
1372           SecondString = new G4ExcitedString( << 
1373           Gquark->Set4Momentum( PkinkQ1 );    << 
1374           Ganti_quark->Set4Momentum( PkinkQ2  << 
1375         }                                     << 
1376       } else {  // Target                     << 
1377         if ( end->GetDefinition()->GetPDGEnco << 
1378           FirstString  = new G4ExcitedString( << 
1379           SecondString = new G4ExcitedString( << 
1380           Ganti_quark->Set4Momentum( PkinkQ2  << 
1381           Gquark->Set4Momentum( PkinkQ1 );    << 
1382         } else {  // Anti_Quark on the end    << 
1383           FirstString  = new G4ExcitedString( << 
1384           SecondString = new G4ExcitedString( << 
1385           Gquark->Set4Momentum( PkinkQ2 );    << 
1386           Ganti_quark->Set4Momentum( PkinkQ1  << 
1387         }                                     << 
1388       }                                       << 
1389     } else {  // Baryon/AntiBaryon            << 
1390       if ( isProjectile ) {  // Projectile    << 
1391         if ( end->GetDefinition()->GetParticl << 
1392              end->GetDefinition()->GetPDGEnco << 
1393           FirstString  = new G4ExcitedString( << 
1394           SecondString = new G4ExcitedString( << 
1395           Gquark->Set4Momentum( PkinkQ1 );    << 
1396           Ganti_quark->Set4Momentum( PkinkQ2  << 
1397         } else {                            / << 
1398           FirstString  = new G4ExcitedString( << 
1399           SecondString = new G4ExcitedString( << 
1400           Ganti_quark->Set4Momentum( PkinkQ1  << 
1401           Gquark->Set4Momentum( PkinkQ2 );    << 
1402         }                                     << 
1403       } else {  // Target                     << 
1404         if ( end->GetDefinition()->GetParticl << 
1405              end->GetDefinition()->GetPDGEnco << 
1406           Gquark->Set4Momentum( PkinkQ1 );    << 
1407           Ganti_quark->Set4Momentum( PkinkQ2  << 
1408           FirstString  = new G4ExcitedString( << 
1409           SecondString = new G4ExcitedString( << 
1410         } else {  // Anti_DiQuark on the end  << 
1411           FirstString  = new G4ExcitedString( << 
1412           SecondString = new G4ExcitedString( << 
1413           Gquark->Set4Momentum( PkinkQ2 );    << 
1414           Ganti_quark->Set4Momentum( PkinkQ1  << 
1415         }                                     << 
1416       }                                       << 
1417     }                                         << 
1418                                               << 
1419     FirstString->SetTimeOfCreation( hadron->G << 
1420     FirstString->SetPosition( hadron->GetPosi << 
1421     SecondString->SetTimeOfCreation( hadron-> << 
1422     SecondString->SetPosition( hadron->GetPos << 
1423                                               << 
1424   } else {  // Kink is impossible             << 
1425                                               << 
1426     if ( isProjectile ) {                     << 
1427       FirstString = new G4ExcitedString( end, << 
1428     } else {                                  << 
1429       FirstString = new G4ExcitedString( end, << 
1430     }                                         << 
1431     FirstString->SetTimeOfCreation( hadron->G << 
1432     FirstString->SetPosition( hadron->GetPosi << 
1433     SecondString = 0;                         << 
1434     // momenta of string ends                 << 
1435     G4LorentzVector HadronMom = hadron->Get4M << 
1436     G4LorentzVector Pstart1 = G4LorentzVector << 
1437     G4LorentzVector   Pend1 = G4LorentzVector << 
1438     G4double Pz  = HadronMom.pz();            << 
1439     G4double Eh  = HadronMom.e();             << 
1440     G4double Pt2 = sqr( HadronMom.px() ) + sq << 
1441     G4double Mt2 = HadronMom.mt2();           << 
1442     G4double Exp = std::sqrt( sqr(Pz) + ( sqr << 
1443     G4double Pzq  = Pz/2.0 - Exp;             << 
1444     G4double Eq   = std::sqrt( sqr(Pzq) + Pt2 << 
1445     G4double Pzqq = Pz/2.0 + Exp;             << 
1446     G4double Eqq  = std::sqrt( sqr(Pzqq) + Pt << 
1447     start->Set4Momentum( Pstart1 );           << 
1448     end->Set4Momentum( Pend1 );               << 
1449     Pstart = Pstart1;  Pend = Pend1;          << 
1450                                               << 
1451   }  // End of "if (Kink)"                    << 
1452                                               << 
1453   //G4cout << "Quarks in the string at creati << 
1454   //       << " " << FirstString->GetLeftPart << 
1455   //       << FirstString << " " << SecondStr << 
1456                                               << 
1457   #ifdef G4_FTFDEBUG                          << 
1458   G4cout << " generated string flavors        << 
1459          << end->GetPDGcode() << G4endl << "  << 
1460          << start->Get4Momentum() << "mass :  << 
1461          << " generated string momenta: Diqua << 
1462          << end->Get4Momentum().mag() << G4en << 
1463          << Pstart + Pend << G4endl << " Orig << 
1464          << hadron->Get4Momentum() << " "<<ha << 
1465   #endif                                      << 
1466                                               << 
1467   return;                                     << 
1468 }                                             << 
1469                                                  1007 
                                                   >> 1008           G4LorentzVector Pstart(tmp,0.);
                                                   >> 1009           G4LorentzVector   Pend(tmp,0.);
1470                                                  1010 
1471 //=========================================== << 1011           if(isProjectile)
                                                   >> 1012           {
                                                   >> 1013            Pstart*=(-1.)*Minus/2.;
                                                   >> 1014            Pend  *=(+1.)*Plus /2.;
                                                   >> 1015           } 
                                                   >> 1016           else
                                                   >> 1017           {
                                                   >> 1018            Pstart*=(+1.)*Plus/2.;
                                                   >> 1019            Pend  *=(-1.)*Minus/2.;
                                                   >> 1020           }
1472                                                  1021 
1473 G4double G4DiffractiveExcitation::ChooseP( G4 << 1022           Momentum=-Pstart.mag();
1474   // Choose an x between Xmin and Xmax with P << 1023           Pstart.setT(Momentum);  // It is assumed that quark has m=0.
1475   // To be improved...                        << 
1476   G4double range = Pmax - Pmin;               << 
1477   if ( Pmin <= 0.0 || range <= 0.0 ) {        << 
1478     G4cout << " Pmin, range : " << Pmin << "  << 
1479     throw G4HadronicException( __FILE__, __LI << 
1480                                "G4Diffractive << 
1481   }                                           << 
1482   G4double P = Pmin * G4Pow::GetInstance()->p << 
1483   //G4double P = (Pmax - Pmin) * G4UniformRan << 
1484   return P;                                   << 
1485 }                                             << 
1486                                                  1024 
                                                   >> 1025           Momentum=-Pend.mag();
                                                   >> 1026           Pend.setT(Momentum);    // It is assumed that di-quark has m=0.
1487                                                  1027 
1488 //=========================================== << 1028     start->Set4Momentum(Pstart);
                                                   >> 1029     end->Set4Momentum(Pend);
                                                   >> 1030           SecondString=0;
                                                   >> 1031         }            // End of kink is impossible 
                                                   >> 1032 
                                                   >> 1033 #ifdef G4_FTFDEBUG
                                                   >> 1034     G4cout << " generated string flavors          " 
                                                   >> 1035                  << start->GetPDGcode() << " / " 
                                                   >> 1036                  << end->GetPDGcode()                    << G4endl;
                                                   >> 1037     G4cout << " generated string momenta:   quark " 
                                                   >> 1038                  << start->Get4Momentum() << "mass : " 
                                                   >> 1039                  <<start->Get4Momentum().mag()           << G4endl;
                                                   >> 1040     G4cout << " generated string momenta: Diquark " 
                                                   >> 1041                  << end ->Get4Momentum() 
                                                   >> 1042                  << "mass : " <<end->Get4Momentum().mag()<< G4endl;
                                                   >> 1043     G4cout << " sum of ends                       " << Pstart+Pend << G4endl;
                                                   >> 1044     G4cout << " Original                          " << hadron->Get4Momentum() << G4endl;
                                                   >> 1045 #endif
1489                                                  1046 
1490 G4ThreeVector G4DiffractiveExcitation::Gaussi << 1047         return;
1491   //  @@ this method is used in FTFModel as w << 1048    
1492   G4double Pt2( 0.0 );                        << 
1493   if ( AveragePt2 <= 0.0 ) {                  << 
1494     Pt2 = 0.0;                                << 
1495   } else {                                    << 
1496     Pt2 = -AveragePt2 * G4Log( 1.0 + G4Unifor << 
1497                                        ( G4Ex << 
1498   }                                           << 
1499   G4double Pt = std::sqrt( Pt2 );             << 
1500   G4double phi = G4UniformRand() * twopi;     << 
1501   return G4ThreeVector( Pt * std::cos( phi ), << 
1502 }                                                1049 }
1503                                                  1050 
1504                                                  1051 
1505 //=========================================== << 1052 // --------- private methods ----------------------
1506                                                  1053 
1507 G4double G4DiffractiveExcitation::GetQuarkFra << 1054 // ---------------------------------------------------------------------
1508   G4double z, yf;                             << 1055 G4double G4DiffractiveExcitation::ChooseP(G4double Pmin, G4double Pmax) const
1509   const G4int maxNumberOfLoops = 10000;       << 1056 {
1510   G4int loopCounter = 0;                      << 1057 // choose an x between Xmin and Xmax with P(x) ~ 1/x
1511   do {                                        << 1058 //  to be improved...
1512     z = zmin + G4UniformRand() * (zmax - zmin << 1059 
1513     yf = z*z + sqr(1.0 - z);                  << 1060   G4double range=Pmax-Pmin;                    
1514   } while ( ( G4UniformRand() > yf ) &&       << 1061 
1515             ++loopCounter < maxNumberOfLoops  << 1062   if ( Pmin <= 0. || range <=0. )
1516   if ( loopCounter >= maxNumberOfLoops ) {    << 1063   {
1517     z = 0.5*(zmin + zmax);  // Just something << 1064     G4cout << " Pmin, range : " << Pmin << " , " << range << G4endl;
1518   }                                           << 1065     throw G4HadronicException(__FILE__, __LINE__, "G4DiffractiveExcitation::ChooseP : Invalid arguments ");
1519   return z;                                   << 1066   }
1520 }                                             << 
1521                                               << 
1522                                               << 
1523 //=========================================== << 
1524                                                  1067 
1525 void G4DiffractiveExcitation::UnpackMeson( co << 1068   G4double P=Pmin * std::pow(Pmax/Pmin,G4UniformRand()); 
1526   G4int absIdPDG = std::abs( IdPDG );         << 1069   return P;
1527   if ( ! ( absIdPDG == 111 || absIdPDG == 221 << 
1528            absIdPDG == 441 || absIdPDG == 443 << 
1529     // All other projectile mesons (including << 
1530     Q1 =  absIdPDG / 100;                     << 
1531     Q2 = (absIdPDG % 100) / 10;               << 
1532     G4int anti = 1 - 2 * ( std::max( Q1, Q2 ) << 
1533     if ( IdPDG < 0 ) anti *= -1;              << 
1534     Q1 *= anti;                               << 
1535     Q2 *= -1 * anti;                          << 
1536   } else {                                    << 
1537     if ( absIdPDG == 441 || absIdPDG == 443 ) << 
1538       Q1 =  4; Q2 = -4;                       << 
1539     } else if ( absIdPDG == 553 ) {           << 
1540       Q1 =  5; Q2 = -5;                       << 
1541     } else {                                  << 
1542       if ( G4UniformRand() < 0.5 ) {          << 
1543   Q1 = 1; Q2 = -1;                            << 
1544       } else {                                << 
1545   Q1 = 2; Q2 = -2;                            << 
1546       }                                       << 
1547     }                                         << 
1548   }                                           << 
1549   return;                                     << 
1550 }                                                1070 }
1551                                                  1071 
1552                                               << 1072 // ---------------------------------------------------------------------
1553 //=========================================== << 1073 G4ThreeVector G4DiffractiveExcitation::GaussianPt(G4double AveragePt2, 
1554                                               << 1074                                                   G4double maxPtSquare) const
1555 void G4DiffractiveExcitation::UnpackBaryon( G << 1075 {            //  @@ this method is used in FTFModel as well. Should go somewhere common!
1556                                             G << 1076 
1557   Q1 = IdPDG          / 1000;                 << 1077   G4double Pt2(0.);
1558   Q2 = (IdPDG % 1000) / 100;                  << 1078         if(AveragePt2 <= 0.) {Pt2=0.;}
1559   Q3 = (IdPDG % 100)  / 10;                   << 1079         else
1560   return;                                     << 1080         {
                                                   >> 1081          Pt2 = -AveragePt2 * std::log(1. + G4UniformRand() * 
                                                   >> 1082                             (std::exp(-maxPtSquare/AveragePt2)-1.));
                                                   >> 1083         }
                                                   >> 1084   G4double Pt=std::sqrt(Pt2);
                                                   >> 1085   G4double phi=G4UniformRand() * twopi;
                                                   >> 1086   return G4ThreeVector (Pt*std::cos(phi), Pt*std::sin(phi), 0.);
1561 }                                                1087 }
1562                                                  1088 
                                                   >> 1089 // ---------------------------------------------------------------------
                                                   >> 1090 G4double G4DiffractiveExcitation::GetQuarkFractionOfKink(G4double zmin, G4double zmax) const
                                                   >> 1091     {
                                                   >> 1092        G4double z, yf;
                                                   >> 1093        do {
                                                   >> 1094            z = zmin + G4UniformRand()*(zmax-zmin);
                                                   >> 1095            yf = z*z +sqr(1 - z);  
                                                   >> 1096            } 
                                                   >> 1097        while (G4UniformRand() > yf); 
                                                   >> 1098        return z;
                                                   >> 1099     }
                                                   >> 1100 // ---------------------------------------------------------------------
                                                   >> 1101 void G4DiffractiveExcitation::UnpackMeson(const G4int IdPDG, G4int &Q1, G4int &Q2) const // Uzhi 7.09.09
                                                   >> 1102     {
                                                   >> 1103        G4int absIdPDG = std::abs(IdPDG);
                                                   >> 1104        Q1=  absIdPDG/ 100;
                                                   >> 1105        Q2= (absIdPDG %100)/10;
                                                   >> 1106      
                                                   >> 1107        G4int anti= 1 -2 * ( std::max( Q1, Q2 ) % 2 );
                                                   >> 1108 
                                                   >> 1109        if (IdPDG < 0 ) anti *=-1;
                                                   >> 1110        Q1 *= anti;
                                                   >> 1111        Q2 *= -1 * anti;
                                                   >> 1112        return;
                                                   >> 1113     }
                                                   >> 1114 // ---------------------------------------------------------------------
                                                   >> 1115 void G4DiffractiveExcitation::UnpackBaryon(G4int IdPDG, 
                                                   >> 1116                               G4int &Q1, G4int &Q2, G4int &Q3) const // Uzhi 7.09.09
                                                   >> 1117     {
                                                   >> 1118        Q1 =  IdPDG           / 1000;
                                                   >> 1119        Q2 = (IdPDG % 1000)  / 100;
                                                   >> 1120        Q3 = (IdPDG % 100)   / 10;
                                                   >> 1121        return;
                                                   >> 1122     }
                                                   >> 1123 // ---------------------------------------------------------------------
                                                   >> 1124 G4int G4DiffractiveExcitation::NewNucleonId(G4int Q1, G4int Q2, G4int Q3) const // Uzhi 7.09.09
                                                   >> 1125     {
                                                   >> 1126        G4int TmpQ(0);
                                                   >> 1127        if( Q3 > Q2 ) 
                                                   >> 1128        {
                                                   >> 1129         TmpQ = Q2;
                                                   >> 1130         Q2 = Q3;
                                                   >> 1131         Q3 = TmpQ;
                                                   >> 1132        } else if( Q3 > Q1 )
                                                   >> 1133        {
                                                   >> 1134         TmpQ = Q1;
                                                   >> 1135         Q1 = Q3;
                                                   >> 1136         Q3 = TmpQ;
                                                   >> 1137        }
                                                   >> 1138 
                                                   >> 1139        if( Q2 > Q1 ) 
                                                   >> 1140        {
                                                   >> 1141         TmpQ = Q1;
                                                   >> 1142         Q1 = Q2;
                                                   >> 1143         Q2 = TmpQ;
                                                   >> 1144        }
1563                                                  1145 
1564 //=========================================== << 1146        G4int NewCode = Q1*1000 + Q2* 100 + Q3*  10 + 2; 
1565                                               << 1147        return NewCode;
1566 G4int G4DiffractiveExcitation::NewNucleonId(  << 1148     }
1567   // Order the three integers in such a way t << 1149 // ---------------------------------------------------------------------
1568   G4int TmpQ( 0 );                            << 1150 G4DiffractiveExcitation::G4DiffractiveExcitation(const G4DiffractiveExcitation &)
1569   if ( Q3 > Q2 ) {                            << 1151 {
1570     TmpQ = Q2;                                << 1152   throw G4HadronicException(__FILE__, __LINE__, "G4DiffractiveExcitation copy contructor not meant to be called");
1571     Q2 = Q3;                                  << 
1572     Q3 = TmpQ;                                << 
1573   } else if ( Q3 > Q1 ) {                     << 
1574     TmpQ = Q1;                                << 
1575     Q1 = Q3;                                  << 
1576     Q3 = TmpQ;                                << 
1577   }                                           << 
1578   if ( Q2 > Q1 ) {                            << 
1579     TmpQ = Q1;                                << 
1580     Q1 = Q2;                                  << 
1581     Q2 = TmpQ;                                << 
1582   }                                           << 
1583   // By now Q1 >= Q2 >= Q3                    << 
1584   G4int NewCode = Q1*1000 + Q2*100 + Q3*10 +  << 
1585   return NewCode;                             << 
1586 }                                                1153 }
1587                                                  1154 
1588                                                  1155 
1589 //=========================================== << 1156 G4DiffractiveExcitation::~G4DiffractiveExcitation()
1590                                               << 1157 {
1591 G4DiffractiveExcitation::G4DiffractiveExcitat << 
1592   throw G4HadronicException( __FILE__, __LINE << 
1593                              "G4DiffractiveEx << 
1594 }                                                1158 }
1595                                                  1159 
1596                                                  1160 
1597 //=========================================== << 1161 const G4DiffractiveExcitation & G4DiffractiveExcitation::operator=(const G4DiffractiveExcitation &)
1598                                               << 1162 {
1599 const G4DiffractiveExcitation & G4Diffractive << 1163   throw G4HadronicException(__FILE__, __LINE__, "G4DiffractiveExcitation = operator meant to be called");
1600   throw G4HadronicException( __FILE__, __LINE << 1164   return *this;
1601                              "G4DiffractiveEx << 
1602   return *this;                               << 
1603 }                                                1165 }
1604                                                  1166 
1605                                                  1167 
1606 //=========================================== << 1168 int G4DiffractiveExcitation::operator==(const G4DiffractiveExcitation &) const
1607                                               << 1169 {
1608 G4bool G4DiffractiveExcitation::operator==( c << 1170   throw G4HadronicException(__FILE__, __LINE__, "G4DiffractiveExcitation == operator meant to be called");
1609   throw G4HadronicException( __FILE__, __LINE << 1171   return false;
1610                              "G4DiffractiveEx << 
1611 }                                                1172 }
1612                                                  1173 
1613                                               << 1174 int G4DiffractiveExcitation::operator!=(const G4DiffractiveExcitation &) const
1614 //=========================================== << 1175 {
1615                                               << 1176   throw G4HadronicException(__FILE__, __LINE__, "G4DiffractiveExcitation != operator meant to be called");
1616 G4bool G4DiffractiveExcitation::operator!= (  << 1177   return true;
1617   throw G4HadronicException( __FILE__, __LINE << 
1618                              "G4DiffractiveEx << 
1619 }                                                1178 }
1620                                               << 
1621                                                  1179