source: trunk/source/processes/hadronic/models/parton_string/diffraction/src/G4DiffractiveExcitation.cc @ 1198

Last change on this file since 1198 was 1196, checked in by garnier, 15 years ago

update CVS release candidate geant4.9.3.01
File size: 58.0 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer                                           *
4// *                                                                  *
5// * The  Geant4 software  is  copyright of the Copyright Holders  of *
6// * the Geant4 Collaboration.  It is provided  under  the terms  and *
7// * conditions of the Geant4 Software License,  included in the file *
8// * LICENSE and available at  http://cern.ch/geant4/license .  These *
9// * include a list of copyright holders.                             *
10// *                                                                  *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work  make  any representation or  warranty, express or implied, *
14// * regarding  this  software system or assume any liability for its *
15// * use.  Please see the license in the file  LICENSE  and URL above *
16// * for the full disclaimer and the limitation of liability.         *
17// *                                                                  *
18// * This  code  implementation is the result of  the  scientific and *
19// * technical work of the GEANT4 collaboration.                      *
20// * By using,  copying,  modifying or  distributing the software (or *
21// * any work based  on the software)  you  agree  to acknowledge its *
22// * use  in  resulting  scientific  publications,  and indicate your *
23// * acceptance of all terms of the Geant4 Software license.          *
24// ********************************************************************
25//
26//
27// $Id: G4DiffractiveExcitation.cc,v 1.16 2009/10/06 10:10:36 vuzhinsk Exp $
28// ------------------------------------------------------------
29//      GEANT 4 class implemetation file
30//
31//      ---------------- G4DiffractiveExcitation --------------
32//             by Gunter Folger, October 1998.
33//        diffractive Excitation used by strings models
34//             Take a projectile and a target
35//             excite the projectile and target
36//  Essential changed by V. Uzhinsky in November - December 2006
37//  in order to put it in a correspondence with original FRITIOF
38//  model. Variant of FRITIOF with nucleon de-excitation is implemented.
39//  Other changes by V.Uzhinsky in May 2007 were introduced to fit
40//  meson-nucleon interactions. Additional changes by V. Uzhinsky
41//  were introduced in December 2006. They treat diffraction dissociation
42//  processes more exactly.
43// ---------------------------------------------------------------------
44
45
46#include "globals.hh"
47#include "Randomize.hh"
48
49#include "G4DiffractiveExcitation.hh"
50#include "G4FTFParameters.hh"
51#include "G4ElasticHNScattering.hh"
52
53#include "G4LorentzRotation.hh"
54#include "G4RotationMatrix.hh"
55#include "G4ThreeVector.hh"
56#include "G4ParticleDefinition.hh"
57#include "G4VSplitableHadron.hh"
58#include "G4ExcitedString.hh"
59#include "G4ParticleTable.hh"
60#include "G4Neutron.hh"
61#include "G4ParticleDefinition.hh"
62
63//#include "G4ios.hh"
64//#include "UZHI_diffraction.hh"
65
66G4DiffractiveExcitation::G4DiffractiveExcitation()
67{
68}
69
70// ---------------------------------------------------------------------
71G4bool G4DiffractiveExcitation::
72  ExciteParticipants(G4VSplitableHadron *projectile, 
73                     G4VSplitableHadron *target,
74                     G4FTFParameters    *theParameters,
75                     G4ElasticHNScattering *theElastic) const  // Uzhi 03.09.09
76{
77// -------------------- Projectile parameters -----------------------
78     G4LorentzVector Pprojectile=projectile->Get4Momentum();
79//G4cout<<"Excite P "<<Pprojectile<<G4endl;
80
81     if(Pprojectile.z() < 0.)
82     {
83       target->SetStatus(2);
84       return false;
85     } 
86
87     G4double ProjectileRapidity = Pprojectile.rapidity();
88
89     G4int    ProjectilePDGcode=projectile->GetDefinition()->GetPDGEncoding();
90//     G4ParticleDefinition * ProjectileDefinition = projectile->GetDefinition();
91     G4int    absProjectilePDGcode=std::abs(ProjectilePDGcode);
92
93     G4bool PutOnMassShell(false);
94//   G4double M0projectile=projectile->GetDefinition()->GetPDGMass(); // With de-excitation
95     G4double M0projectile = Pprojectile.mag();                       // Without de-excitation
96
97     if(M0projectile < projectile->GetDefinition()->GetPDGMass())
98     {
99      PutOnMassShell=true;
100      M0projectile=projectile->GetDefinition()->GetPDGMass();
101     }
102
103     G4double M0projectile2 = M0projectile * M0projectile;
104
105     G4double ProjectileDiffStateMinMass=theParameters->GetProjMinDiffMass();
106     G4double ProjectileNonDiffStateMinMass=theParameters->GetProjMinNonDiffMass();
107     G4double ProbProjectileDiffraction=theParameters->GetProbabilityOfProjDiff();
108
109// -------------------- Target parameters -------------------------
110     G4int    TargetPDGcode=target->GetDefinition()->GetPDGEncoding();
111     G4int    absTargetPDGcode=std::abs(TargetPDGcode);
112
113     G4LorentzVector Ptarget=target->Get4Momentum();
114
115     G4double M0target = Ptarget.mag();
116
117//G4cout<<"Excite T "<<Ptarget<<G4endl;
118//G4cout<<"PDGs "<<ProjectilePDGcode<<" "<<TargetPDGcode<<G4endl;
119//G4cout<<"Masses "<<M0projectile<<" "<<M0target<<G4endl;
120
121     G4double TargetRapidity = Ptarget.rapidity();
122
123     if(M0target < target->GetDefinition()->GetPDGMass())
124     {
125      PutOnMassShell=true;
126      M0target=target->GetDefinition()->GetPDGMass();
127     }
128
129     G4double M0target2 = M0target * M0target; 
130 
131     G4double TargetDiffStateMinMass=theParameters->GetTarMinDiffMass();   
132     G4double TargetNonDiffStateMinMass=theParameters->GetTarMinNonDiffMass();   
133     G4double ProbTargetDiffraction=theParameters->GetProbabilityOfTarDiff();
134
135     G4double AveragePt2=theParameters->GetAveragePt2();
136
137     G4double ProbOfDiffraction=ProbProjectileDiffraction +
138                                ProbTargetDiffraction;
139
140     G4double SumMasses=M0projectile+M0target+200.*MeV;
141
142// Kinematical properties of the interactions --------------
143     G4LorentzVector Psum;      // 4-momentum in CMS
144     Psum=Pprojectile+Ptarget;
145     G4double S=Psum.mag2(); 
146
147// Transform momenta to cms and then rotate parallel to z axis;
148     G4LorentzRotation toCms(-1*Psum.boostVector());
149
150     G4LorentzVector Ptmp=toCms*Pprojectile;
151     if ( Ptmp.pz() <= 0. )
152     {
153       target->SetStatus(2); 
154   // "String" moving backwards in  CMS, abort collision !!
155       return false;
156     }
157
158     toCms.rotateZ(-1*Ptmp.phi());
159     toCms.rotateY(-1*Ptmp.theta());
160
161     G4LorentzRotation toLab(toCms.inverse());
162
163     Pprojectile.transform(toCms);
164     Ptarget.transform(toCms);
165
166     G4double PZcms2, PZcms;
167
168     G4double SqrtS=std::sqrt(S);
169//G4cout<<" SqrtS "<<SqrtS<<" 2200 "<<G4endl;
170     if(absProjectilePDGcode > 1000 && SqrtS < 2300*MeV && SqrtS < SumMasses)
171     {target->SetStatus(2);  return false;}  // The model cannot work for
172                                             // p+p-interactions
173                                             // at Plab < 1.62 GeV/c.
174
175     if(( absProjectilePDGcode == 211 || ProjectilePDGcode ==  111) && 
176        (SqrtS < 1600*MeV) && (SqrtS < SumMasses))
177     {target->SetStatus(2);  return false;}  // The model cannot work for
178                                             // Pi+p-interactions
179                                             // at Plab < 1. GeV/c.
180
181     if(( (absProjectilePDGcode == 321) || (ProjectilePDGcode == -311)   ||
182          (absProjectilePDGcode == 311) || (absProjectilePDGcode == 130) ||
183          (absProjectilePDGcode == 310)) && (SqrtS < 1600*MeV) && (SqrtS < SumMasses)) 
184     {target->SetStatus(2);  return false;}  // The model cannot work for
185                                             // K+p-interactions
186                                             // at Plab < ??? GeV/c.  ???
187
188     PZcms2=(S*S+M0projectile2*M0projectile2+M0target2*M0target2-
189             2*S*M0projectile2 - 2*S*M0target2 - 2*M0projectile2*M0target2)
190             /4./S;
191
192     if(PZcms2 < 0)
193     {target->SetStatus(2);  return false;}   // It can be in an interaction with
194                                              // off-shell nuclear nucleon
195
196     PZcms = std::sqrt(PZcms2);
197
198     if(PutOnMassShell)
199     {
200      if(Pprojectile.z() > 0.)
201      {
202       Pprojectile.setPz( PZcms);
203       Ptarget.setPz(    -PZcms);
204      }
205      else
206      {
207       Pprojectile.setPz(-PZcms);
208       Ptarget.setPz(     PZcms);
209      };
210
211      Pprojectile.setE(std::sqrt(M0projectile2                  +
212                                 Pprojectile.x()*Pprojectile.x()+
213                                 Pprojectile.y()*Pprojectile.y()+
214                                 PZcms2));
215      Ptarget.setE(std::sqrt(M0target2              +
216                             Ptarget.x()*Ptarget.x()+
217                             Ptarget.y()*Ptarget.y()+
218                             PZcms2));
219     }
220
221     G4double maxPtSquare; // = PZcms2;
222
223// Charge exchange can be possible for baryons -----------------
224
225//G4cout<<1./std::exp(-2.0*(ProjectileRapidity - TargetRapidity))<<G4endl;
226//G4int Uzhi; G4cin>>Uzhi;
227// Getting the values needed for exchange ----------------------
228     G4double MagQuarkExchange        =theParameters->GetMagQuarkExchange();
229     G4double SlopeQuarkExchange      =theParameters->GetSlopeQuarkExchange();
230     G4double DeltaProbAtQuarkExchange=theParameters->GetDeltaProbAtQuarkExchange();
231
232//     G4double NucleonMass=
233              (G4ParticleTable::GetParticleTable()->FindParticle(2112))->GetPDGMass();     
234     G4double DeltaMass=
235              (G4ParticleTable::GetParticleTable()->FindParticle(2224))->GetPDGMass();
236
237// Check for possible quark excjane -----------------------------------
238     if(G4UniformRand() < MagQuarkExchange*
239        std::exp(-SlopeQuarkExchange*(ProjectileRapidity - TargetRapidity)))
240     {   
241
242//G4cout<<"Exchange "<<G4endl;
243
244      G4int NewProjCode(0), NewTargCode(0);
245
246//G4bool StandardExcitation(false); // =================================
247      G4int ProjQ1(0), ProjQ2(0), ProjQ3(0);
248
249//  Projectile unpacking --------------------------
250      if(absProjectilePDGcode < 1000 )
251      {    // projectile is meson -----------------
252       UnpackMeson(ProjectilePDGcode, ProjQ1, ProjQ2); 
253      } else 
254      {    // projectile is baryon ----------------
255       UnpackBaryon(ProjectilePDGcode, ProjQ1, ProjQ2, ProjQ3);
256      } // End of the hadron's unpacking ----------
257
258//  Target unpacking ------------------------------
259      G4int TargQ1(0), TargQ2(0), TargQ3(0);
260      UnpackBaryon(TargetPDGcode, TargQ1, TargQ2, TargQ3); 
261
262// Sampling of exchanged quarks -------------------
263      G4int ProjExchangeQ(0);
264      G4int TargExchangeQ(0);
265
266//G4cout<<"Targ "<<TargQ1<<" "<<TargQ2<<" "<<TargQ3<<G4endl;
267      if(absProjectilePDGcode < 1000 )
268      {    // projectile is meson -----------------
269//G4cout<<"Proj Q1 Q2 "<<ProjQ1<<" "<<ProjQ2<<G4endl;
270
271       if(ProjQ1 > 0 ) // ProjQ1 is quark
272       { 
273        ProjExchangeQ = ProjQ1;
274        if(ProjExchangeQ != TargQ1)
275        {
276         TargExchangeQ = TargQ1; TargQ1=ProjExchangeQ; ProjQ1=TargExchangeQ;
277        } else
278        if(ProjExchangeQ != TargQ2)
279        {
280         TargExchangeQ = TargQ2; TargQ2=ProjExchangeQ; ProjQ1=TargExchangeQ;
281        } else         
282        {
283         TargExchangeQ = TargQ3;  TargQ3=ProjExchangeQ; ProjQ1=TargExchangeQ;
284        }
285//G4cout<<" Pr Tr "<<ProjQ1<<" "<<TargQ1<<" "<<TargQ2<<" "<<TargQ3<<G4endl;
286       } else          // ProjQ2 is quark
287       { 
288        ProjExchangeQ = ProjQ2;
289        if(ProjExchangeQ != TargQ1)
290        {
291         TargExchangeQ = TargQ1; TargQ1=ProjExchangeQ; ProjQ2=TargExchangeQ;
292        } else
293        if(ProjExchangeQ != TargQ2)
294        {
295         TargExchangeQ = TargQ2; TargQ2=ProjExchangeQ; ProjQ2=TargExchangeQ;
296        } else 
297        {
298         TargExchangeQ = TargQ3;  TargQ3=ProjExchangeQ; ProjQ2=TargExchangeQ;
299        }
300
301       } // End of if(ProjQ1 > 0 ) // ProjQ1 is quark
302
303       G4int aProjQ1=std::abs(ProjQ1);
304       G4int aProjQ2=std::abs(ProjQ2);
305       if(aProjQ1 == aProjQ2)          {NewProjCode = 111;} // Pi0-meson
306       else  // |ProjQ1| # |ProjQ2|
307       {
308        if(aProjQ1 > aProjQ2)          {NewProjCode = aProjQ1*100+aProjQ2*10+1;}
309        else                           {NewProjCode = aProjQ2*100+aProjQ1*10+1;}
310       }
311
312//       projectile->SetDefinition(
313//       G4ParticleTable::GetParticleTable()->FindParticle(NewProjCode));
314
315       NewTargCode = NewNucleonId(TargQ1, TargQ2, TargQ3);
316 
317       if( (TargQ1 == TargQ2) && (TargQ1 == TargQ3) &&
318           (SqrtS > M0projectile+DeltaMass))           {NewTargCode +=2;} //Create Delta isobar
319
320       else if( target->GetDefinition()->GetPDGiIsospin() == 3 )         //Delta was the target
321       { if(G4UniformRand() > DeltaProbAtQuarkExchange){NewTargCode +=2;} //Save   Delta isobar
322         else                                          {}     // De-excite initial Delta isobar
323       }
324
325       else if((G4UniformRand() < DeltaProbAtQuarkExchange) &&         //Nucleon was the target
326               (SqrtS > M0projectile+DeltaMass))      {NewTargCode +=2;}  //Create Delta isobar
327       else                                           {}                 //Save initial nucleon
328
329       target->SetDefinition(
330       G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode)); 
331//G4cout<<"New PDGs "<<NewProjCode<<" "<<NewTargCode<<G4endl;
332//G4int Uzhi; G4cin>>Uzhi;
333       //}
334      } else 
335      {    // projectile is baryon ----------------
336G4double Same=0.; //0.3; //0.5;
337       G4bool ProjDeltaHasCreated(false);
338       G4bool TargDeltaHasCreated(false);
339 
340       G4double Ksi=G4UniformRand();
341       if(G4UniformRand() < 0.5)     // Sampling exchange quark from proj. or targ.
342       {   // Sampling exchanged quark from the projectile ---
343//G4cout<<"Proj Exc"<<G4endl;
344        if( Ksi < 0.333333 ) 
345        {ProjExchangeQ = ProjQ1;}
346        else if( (0.333333 <= Ksi) && (Ksi < 0.666667))
347        {ProjExchangeQ = ProjQ2;}
348        else
349        {ProjExchangeQ = ProjQ3;}
350//G4cout<<"ProjExchangeQ "<<ProjExchangeQ<<G4endl;
351
352        if((ProjExchangeQ != TargQ1)||(G4UniformRand()<Same)) // Vova
353        {
354         TargExchangeQ = TargQ1; TargQ1=ProjExchangeQ; ProjExchangeQ=TargExchangeQ;
355        } else
356        if((ProjExchangeQ != TargQ2)||(G4UniformRand()<Same)) // Vova
357        {
358         TargExchangeQ = TargQ2; TargQ2=ProjExchangeQ; ProjExchangeQ=TargExchangeQ;
359        } else 
360        {
361         TargExchangeQ = TargQ3;  TargQ3=ProjExchangeQ; ProjExchangeQ=TargExchangeQ;
362        }
363//G4cout<<"ProjExchangeQ "<<ProjExchangeQ<<G4endl;
364        if( Ksi < 0.333333 ) 
365        {ProjQ1=ProjExchangeQ;}
366        else if( (0.333333 <= Ksi) && (Ksi < 0.666667))
367        {ProjQ2=ProjExchangeQ;}
368        else
369        {ProjQ3=ProjExchangeQ;}
370
371/*
372        NewProjCode = NewNucleonId(ProjQ1, ProjQ2, ProjQ3); // *****************************
373
374        if((ProjQ1==ProjQ2) && (ProjQ1==ProjQ3)) {NewProjCode +=2; ProjDeltaHasCreated=true;}
375        else if(projectile->GetDefinition()->GetPDGiIsospin() == 3)// Projectile was Delta
376        { if(G4UniformRand() > DeltaProbAtQuarkExchange)
377                                                 {NewProjCode +=2; ProjDeltaHasCreated=true;}
378          else                                   {NewProjCode +=0; ProjDeltaHasCreated=false;}
379        }
380        else                                                       // Projectile was Nucleon
381        {
382         if((G4UniformRand() < DeltaProbAtQuarkExchange) && (SqrtS > DeltaMass+M0target))
383                                                 {NewProjCode +=2; ProjDeltaHasCreated=true;}
384         else                                    {NewProjCode +=0; ProjDeltaHasCreated=false;}
385        }
386
387        NewTargCode = NewNucleonId(TargQ1, TargQ2, TargQ3); // *****************************
388
389        if((TargQ1==TargQ2) && (TargQ1==TargQ3)) {NewTargCode +=2; TargDeltaHasCreated=true;} 
390        else if(target->GetDefinition()->GetPDGiIsospin() == 3)    // Target was Delta
391        { if(G4UniformRand() > DeltaProbAtQuarkExchange)
392                                                 {NewTargCode +=2; TargDeltaHasCreated=true;}
393          else                                   {NewTargCode +=0; TargDeltaHasCreated=false;}
394        }
395        else                                                       // Target was Nucleon
396        {
397         if((G4UniformRand() < DeltaProbAtQuarkExchange) && (SqrtS > M0projectile+DeltaMass))
398                                                 {NewTargCode +=2; TargDeltaHasCreated=true;}
399         else                                    {NewTargCode +=0; TargDeltaHasCreated=false;}
400        }         
401
402        if((absProjectilePDGcode == NewProjCode) && (absTargetPDGcode == NewTargCode))
403        { // Nothing was changed! It is not right!?
404        }
405*/
406/*-----------------------------------------------------------------------------------------
407        if(!ProjDeltaHasCreated && !TargDeltaHasCreated) // No Deltas were created then
408        { if( G4UniformRand() < 0.5) {NewProjCode +=2; ProjDeltaHasCreated = true;}
409          else                       {NewTargCode +=2; TargDeltaHasCreated = true;}
410        }
411
412        if(!(ProjDeltaHasCreated && TargDeltaHasCreated))
413        {
414         if((G4UniformRand() < DeltaProbAtQuarkExchange) && (SqrtS >2600.))
415         {
416G4cout<<"2 Deltas "<<G4endl;
417          if(!ProjDeltaHasCreated)
418               {NewProjCode +=2; ProjDeltaHasCreated = true;}
419          else {NewTargCode +=2; TargDeltaHasCreated = true;} // For Delta isobar
420         }
421        }
422*/
423       } else
424       {   // Sampling exchanged quark from the target -------
425//G4cout<<"Targ Exc"<<G4endl;
426        if( Ksi < 0.333333 ) 
427        {TargExchangeQ = TargQ1;}
428        else if( (0.333333 <= Ksi) && (Ksi < 0.666667)) 
429        {TargExchangeQ = TargQ2;}
430        else
431        {TargExchangeQ = TargQ3;}
432//G4cout<<"TargExchangeQ "<<TargExchangeQ<<G4endl;
433        if((TargExchangeQ != ProjQ1)||(G4UniformRand()<Same)) // Vova
434        {
435         ProjExchangeQ = ProjQ1; ProjQ1=TargExchangeQ; TargExchangeQ=ProjExchangeQ;
436        } else
437        if((TargExchangeQ != ProjQ2)||(G4UniformRand()<Same)) // Vova
438        {
439         ProjExchangeQ = ProjQ2; ProjQ2=TargExchangeQ; TargExchangeQ=ProjExchangeQ;
440        } else 
441        {
442         ProjExchangeQ = ProjQ3;  ProjQ3=TargExchangeQ; TargExchangeQ=ProjExchangeQ;
443        }
444//G4cout<<"TargExchangeQ "<<TargExchangeQ<<G4endl;
445        if( Ksi < 0.333333 ) 
446        {TargQ1=TargExchangeQ;}
447        else if( (0.333333 <= Ksi) && (Ksi < 0.666667)) 
448        {TargQ2=TargExchangeQ;}
449        else
450        {TargQ3=TargExchangeQ;}
451/*
452        NewProjCode = NewNucleonId(ProjQ1, ProjQ2, ProjQ3);
453        if( (ProjQ1 == ProjQ2) && (ProjQ1 == ProjQ3) ) 
454        {NewProjCode +=2; ProjDeltaHasCreated = true;}
455
456        NewTargCode = NewNucleonId(TargQ1, TargQ2, TargQ3);
457        if( (TargQ1 == TargQ2) && (TargQ1 == TargQ3) ) 
458        {NewTargCode +=2; TargDeltaHasCreated = true;}
459
460        if(!ProjDeltaHasCreated && !TargDeltaHasCreated)
461        {NewTargCode +=2; TargDeltaHasCreated = true;}
462
463        if(!(ProjDeltaHasCreated && TargDeltaHasCreated))
464        {
465         if((G4UniformRand() < DeltaProbAtQuarkExchange) && (SqrtS >2600.))
466         {
467G4cout<<"2 Deltas "<<G4endl;
468          if(!ProjDeltaHasCreated)
469               {NewProjCode +=2; ProjDeltaHasCreated = true;}
470          else {NewTargCode +=2; TargDeltaHasCreated = true;}
471         }
472        }
473*/
474       } // End of sampling baryon
475
476       NewProjCode = NewNucleonId(ProjQ1, ProjQ2, ProjQ3); // *****************************
477
478       if((ProjQ1==ProjQ2) && (ProjQ1==ProjQ3)) {NewProjCode +=2; ProjDeltaHasCreated=true;}
479       else if(projectile->GetDefinition()->GetPDGiIsospin() == 3)// Projectile was Delta
480       { if(G4UniformRand() > DeltaProbAtQuarkExchange)
481                                                {NewProjCode +=2; ProjDeltaHasCreated=true;} 
482         else                                   {NewProjCode +=0; ProjDeltaHasCreated=false;}
483       }
484       else                                                       // Projectile was Nucleon
485       {
486        if((G4UniformRand() < DeltaProbAtQuarkExchange) && (SqrtS > DeltaMass+M0target)) 
487                                                {NewProjCode +=2; ProjDeltaHasCreated=true;}
488        else                                    {NewProjCode +=0; ProjDeltaHasCreated=false;}
489       } 
490
491       NewTargCode = NewNucleonId(TargQ1, TargQ2, TargQ3); // *****************************
492
493       if((TargQ1==TargQ2) && (TargQ1==TargQ3)) {NewTargCode +=2; TargDeltaHasCreated=true;} 
494       else if(target->GetDefinition()->GetPDGiIsospin() == 3)    // Target was Delta
495       { if(G4UniformRand() > DeltaProbAtQuarkExchange)
496                                                {NewTargCode +=2; TargDeltaHasCreated=true;} 
497         else                                   {NewTargCode +=0; TargDeltaHasCreated=false;}
498       }
499       else                                                       // Target was Nucleon
500       {
501        if((G4UniformRand() < DeltaProbAtQuarkExchange) && (SqrtS > M0projectile+DeltaMass)) 
502                                                {NewTargCode +=2; TargDeltaHasCreated=true;}
503        else                                    {NewTargCode +=0; TargDeltaHasCreated=false;}
504       }         
505
506       if((absProjectilePDGcode == NewProjCode) && (absTargetPDGcode == NewTargCode))
507       { // Nothing was changed! It is not right!?
508       }
509// Forming baryons --------------------------------------------------
510
511//       if( ProjDeltaHasCreated && TargDeltaHasCreated ) {StandardExcitation = true;}
512
513//G4cout<<"New PDG "<<NewProjCode<<" "<<NewTargCode<<G4endl;
514//G4int Uzhi; G4cin>>Uzhi;
515      } // End of if projectile is baryon ---------------------------
516
517
518// If we assume that final state hadrons after the charge exchange will be
519// in the ground states, we have to put ----------------------------------
520      if(M0projectile < 
521         (G4ParticleTable::GetParticleTable()->FindParticle(NewProjCode))->GetPDGMass())
522      {
523       M0projectile=
524         (G4ParticleTable::GetParticleTable()->FindParticle(NewProjCode))->GetPDGMass();
525       M0projectile2 = M0projectile * M0projectile;
526      }
527
528      if(M0target < 
529         (G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode))->GetPDGMass())
530      {
531       M0target=
532         (G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode))->GetPDGMass();
533       M0target2 = M0target * M0target;
534      }
535
536      PZcms2=(S*S+M0projectile2*M0projectile2+M0target2*M0target2-
537             2*S*M0projectile2 - 2*S*M0target2 - 2*M0projectile2*M0target2)
538             /4./S;
539//G4cout<<"New Masses "<<M0projectile<<" "<<M0target<<G4endl;
540//G4cout<<"PZcms2 "<<PZcms2<<G4endl;
541
542      if(PZcms2 < 0) {return false;}  // It can be if energy is not sufficient for Delta
543//----------------------------------------------------------
544      projectile->SetDefinition(
545                  G4ParticleTable::GetParticleTable()->FindParticle(NewProjCode)); 
546
547      target->SetDefinition(
548                  G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode)); 
549//----------------------------------------------------------
550      PZcms = std::sqrt(PZcms2);
551
552      Pprojectile.setPz( PZcms);
553      Pprojectile.setE(std::sqrt(M0projectile2+PZcms2));
554
555      Ptarget.setPz(    -PZcms);
556      Ptarget.setE(std::sqrt(M0target2+PZcms2));
557
558//G4cout<<"Proj "<<Pprojectile<<" "<<Pprojectile.mag()<<G4endl;
559//G4cout<<"Targ "<<Ptarget<<" "<<Ptarget.mag()<<G4endl;
560
561//      if(!StandardExcitation)
562      {
563       Pprojectile.transform(toLab);
564       Ptarget.transform(toLab);
565
566       projectile->SetTimeOfCreation(target->GetTimeOfCreation());
567       projectile->SetPosition(target->GetPosition());
568
569       projectile->Set4Momentum(Pprojectile);
570       target->Set4Momentum(Ptarget);
571
572       G4bool Result= theElastic->ElasticScattering (projectile,target,theParameters);
573//G4cout<<"1 Delta result "<<Result<<"**********"<<G4endl;
574       return Result;
575      } 
576/*
577else
578      {
579       if(M0projectile > ProjectileDiffStateMinMass)
580       { ProjectileDiffStateMinMass +=200.*MeV; ProjectileNonDiffStateMinMass +=200.*MeV;}
581
582       if(M0target > TargetDiffStateMinMass)
583       { TargetDiffStateMinMass +=200.*MeV; TargetNonDiffStateMinMass +=200.*MeV;}
584
585       ProbOfDiffraction         = 1.;
586       ProbProjectileDiffraction =0.5;
587       ProbTargetDiffraction     =0.5;
588G4cout<<"Exc DD "<<M0projectile<<" "<<M0target<<" --------------------"<<G4endl;
589//       After that standard FTF excitation
590      }
591*/
592     }  // End of charge exchange part ------------------------------
593
594// ------------------------------------------------------------------
595//     G4double ProbOfDiffraction=ProbProjectileDiffraction +
596//                                ProbTargetDiffraction;
597//G4cout<<ProbOfDiffraction<<" "<<ProbProjectileDiffraction<<" "<<ProbTargetDiffraction<<G4endl;
598//G4int Uzhi; G4cin>>Uzhi;
599     if(ProbOfDiffraction!=0.)
600     {
601      ProbProjectileDiffraction/=ProbOfDiffraction;
602     }
603     else
604     {
605      ProbProjectileDiffraction=0.;
606     }
607
608//ProbOfDiffraction=1.; //0.5; // Vova
609
610     G4double ProjectileDiffStateMinMass2    = ProjectileDiffStateMinMass    *
611                                               ProjectileDiffStateMinMass;
612     G4double ProjectileNonDiffStateMinMass2 = ProjectileNonDiffStateMinMass *
613                                               ProjectileNonDiffStateMinMass;
614
615     G4double TargetDiffStateMinMass2        = TargetDiffStateMinMass        *
616                                               TargetDiffStateMinMass;
617     G4double TargetNonDiffStateMinMass2     = TargetNonDiffStateMinMass     *
618                                               TargetNonDiffStateMinMass;
619
620
621     G4double Pt2;
622     G4double ProjMassT2, ProjMassT;
623     G4double TargMassT2, TargMassT;
624     G4double PMinusMin, PMinusMax;
625//   G4double PPlusMin , PPlusMax;
626     G4double TPlusMin , TPlusMax;
627     G4double PMinusNew, PPlusNew, TPlusNew, TMinusNew;
628
629     G4LorentzVector Qmomentum;
630     G4double Qminus, Qplus;
631
632     G4int whilecount=0;
633//   Choose a process ---------------------------
634
635     if(G4UniformRand() < ProbOfDiffraction)
636       {
637        if(G4UniformRand() < ProbProjectileDiffraction)
638        { //-------- projectile diffraction ---------------
639//G4cout<<"Proj difr"<<G4endl;
640         do {
641//             Generate pt
642//             if (whilecount++ >= 500 && (whilecount%100)==0)
643//               G4cout << "G4DiffractiveExcitation::ExciteParticipants possibly looping"
644//               << ", loop count/ maxPtSquare : "
645//               << whilecount << " / " << maxPtSquare << G4endl;
646             if (whilecount > 1000 )
647             {
648              Qmomentum=G4LorentzVector(0.,0.,0.,0.);
649              target->SetStatus(2);  return false;    //  Ignore this interaction
650             };
651// --------------- Check that the interaction is possible -----------
652             ProjMassT2=ProjectileDiffStateMinMass2;
653             ProjMassT =ProjectileDiffStateMinMass;
654
655             TargMassT2=M0target2;
656             TargMassT =M0target;
657
658             PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2-
659                     2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
660                    /4./S;
661//G4cout<<"PZcms2 < 0 false "<<PZcms2<<G4endl;
662             if(PZcms2 < 0 ) 
663             {
664               target->SetStatus(2); 
665               return false;
666             }
667             maxPtSquare=PZcms2;
668
669             Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0);
670             Pt2=G4ThreeVector(Qmomentum.vect()).mag2();
671
672             ProjMassT2=ProjectileDiffStateMinMass2+Pt2;
673             ProjMassT =std::sqrt(ProjMassT2);
674
675             TargMassT2=M0target2+Pt2;
676             TargMassT =std::sqrt(TargMassT2);
677
678             PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2-
679                     2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
680                    /4./S;
681
682             if(PZcms2 < 0 ) continue;
683             PZcms =std::sqrt(PZcms2);
684
685             PMinusMin=std::sqrt(ProjMassT2+PZcms2)-PZcms;
686             PMinusMax=SqrtS-TargMassT;
687
688             PMinusNew=ChooseP(PMinusMin, PMinusMax);
689// PMinusNew=1./sqrt(1./PMinusMin-G4UniformRand()*(1./PMinusMin-1./PMinusMax));
690
691             TMinusNew=SqrtS-PMinusNew;
692             Qminus=Ptarget.minus()-TMinusNew;
693             TPlusNew=TargMassT2/TMinusNew;
694             Qplus=Ptarget.plus()-TPlusNew;
695
696             Qmomentum.setPz( (Qplus-Qminus)/2 );
697             Qmomentum.setE(  (Qplus+Qminus)/2 );
698          } while (
699((Pprojectile+Qmomentum).mag2() <  ProjectileDiffStateMinMass2) ||  //No without excitation
700((Ptarget    -Qmomentum).mag2() <  M0target2                  ));
701        }
702        else
703        { // -------------- Target diffraction ----------------
704         do {
705//G4cout<<"Targ difr"<<G4endl;
706//             Generate pt
707//             if (whilecount++ >= 500 && (whilecount%100)==0)
708//               G4cout << "G4DiffractiveExcitation::ExciteParticipants possibly looping"
709//               << ", loop count/ maxPtSquare : "
710//               << whilecount << " / " << maxPtSquare << G4endl;
711             if (whilecount > 1000 )
712             {
713              Qmomentum=G4LorentzVector(0.,0.,0.,0.);
714              target->SetStatus(2);  return false;    //  Ignore this interaction
715             };
716// --------------- Check that the interaction is possible -----------
717             ProjMassT2=M0projectile2;
718             ProjMassT =M0projectile;
719
720             TargMassT2=TargetDiffStateMinMass2;
721             TargMassT =TargetDiffStateMinMass;
722
723             PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2-
724                     2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
725                    /4./S;
726
727             if(PZcms2 < 0 ) 
728             {
729               target->SetStatus(2); 
730               return false;
731             }
732             maxPtSquare=PZcms2;
733
734             Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0);
735             Pt2=G4ThreeVector(Qmomentum.vect()).mag2();
736
737             ProjMassT2=M0projectile2+Pt2;
738             ProjMassT =std::sqrt(ProjMassT2);
739
740             TargMassT2=TargetDiffStateMinMass2+Pt2;
741             TargMassT =std::sqrt(TargMassT2);
742
743             PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2-
744                     2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
745                    /4./S;
746
747             if(PZcms2 < 0 ) continue;
748             PZcms =std::sqrt(PZcms2);
749
750             TPlusMin=std::sqrt(TargMassT2+PZcms2)-PZcms;
751             TPlusMax=SqrtS-ProjMassT;
752
753             TPlusNew=ChooseP(TPlusMin, TPlusMax);
754
755//TPlusNew=TPlusMax;
756
757             PPlusNew=SqrtS-TPlusNew;
758             Qplus=PPlusNew-Pprojectile.plus();
759             PMinusNew=ProjMassT2/PPlusNew;
760             Qminus=PMinusNew-Pprojectile.minus();
761
762             Qmomentum.setPz( (Qplus-Qminus)/2 );
763             Qmomentum.setE(  (Qplus+Qminus)/2 );
764
765          } while (
766 ((Pprojectile+Qmomentum).mag2() <  M0projectile2          ) ||  //No without excitation
767 ((Ptarget    -Qmomentum).mag2() <  TargetDiffStateMinMass2));
768         }
769        }
770        else  //----------- Non-diffraction process ------------
771        {
772//G4cout<<"Non difr"<<G4endl;
773         do {
774//             Generate pt
775//             if (whilecount++ >= 500 && (whilecount%100)==0)
776//               G4cout << "G4DiffractiveExcitation::ExciteParticipants possibly looping"
777//               << ", loop count/ maxPtSquare : "
778//               << whilecount << " / " << maxPtSquare << G4endl;
779             if (whilecount > 1000 )
780             {
781              Qmomentum=G4LorentzVector(0.,0.,0.,0.);
782              target->SetStatus(2);  return false;    //  Ignore this interaction
783             };
784// --------------- Check that the interaction is possible -----------
785             ProjMassT2=ProjectileNonDiffStateMinMass2;
786             ProjMassT =ProjectileNonDiffStateMinMass;
787
788             TargMassT2=TargetNonDiffStateMinMass2;
789             TargMassT =TargetNonDiffStateMinMass;
790
791             PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2-
792                    2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
793                   /4./S;
794
795             if(PZcms2 < 0 ) 
796             {
797               target->SetStatus(2); 
798               return false;
799             }
800             maxPtSquare=PZcms2;
801             Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0);
802             Pt2=G4ThreeVector(Qmomentum.vect()).mag2();
803
804             ProjMassT2=ProjectileNonDiffStateMinMass2+Pt2;
805             ProjMassT =std::sqrt(ProjMassT2);
806
807             TargMassT2=TargetNonDiffStateMinMass2+Pt2;
808             TargMassT =std::sqrt(TargMassT2);
809
810             PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2-
811                    2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
812                   /4./S;
813
814             if(PZcms2 < 0 ) continue;
815             PZcms =std::sqrt(PZcms2);
816
817             PMinusMin=std::sqrt(ProjMassT2+PZcms2)-PZcms;
818             PMinusMax=SqrtS-TargMassT;
819
820             PMinusNew=ChooseP(PMinusMin, PMinusMax);
821
822             Qminus=PMinusNew-Pprojectile.minus();
823
824             TPlusMin=std::sqrt(TargMassT2+PZcms2)-PZcms;
825//           TPlusMax=SqrtS-PMinusNew;                     
826             TPlusMax=SqrtS-ProjMassT;     
827
828             TPlusNew=ChooseP(TPlusMin, TPlusMax);
829
830             Qplus=-(TPlusNew-Ptarget.plus());
831
832             Qmomentum.setPz( (Qplus-Qminus)/2 );
833             Qmomentum.setE(  (Qplus+Qminus)/2 );
834
835       } while (
836 ((Pprojectile+Qmomentum).mag2() <  ProjectileNonDiffStateMinMass2) || //No double Diffraction
837 ((Ptarget    -Qmomentum).mag2() <  TargetNonDiffStateMinMass2    ));
838         }
839
840           Pprojectile += Qmomentum;
841           Ptarget     -= Qmomentum;
842
843//G4cout<<"Masses "<<Pprojectile.mag()<<" "<<Ptarget.mag()<<G4endl;
844//G4int Uzhi; G4cin>>Uzhi;
845
846// Transform back and update SplitableHadron Participant.
847           Pprojectile.transform(toLab);
848           Ptarget.transform(toLab);
849
850// Calculation of the creation time ---------------------
851      projectile->SetTimeOfCreation(target->GetTimeOfCreation());
852      projectile->SetPosition(target->GetPosition());
853// Creation time and position of target nucleon were determined at
854// ReggeonCascade() of G4FTFModel
855// ------------------------------------------------------
856
857           projectile->Set4Momentum(Pprojectile);
858           target->Set4Momentum(Ptarget);
859
860           projectile->IncrementCollisionCount(1);
861           target->IncrementCollisionCount(1);
862
863           return true;
864}
865
866// ---------------------------------------------------------------------
867//G4ExcitedString * G4DiffractiveExcitation::
868//           String(G4VSplitableHadron * hadron, G4bool isProjectile) const
869void G4DiffractiveExcitation::CreateStrings(G4VSplitableHadron * hadron, 
870                                            G4bool isProjectile,
871                                            G4ExcitedString * &FirstString, 
872                                            G4ExcitedString * &SecondString,
873                                            G4FTFParameters *theParameters) const
874{
875        hadron->SplitUp();
876        G4Parton *start= hadron->GetNextParton();
877        if ( start==NULL)
878        { G4cout << " G4FTFModel::String() Error:No start parton found"<< G4endl;
879          FirstString=0; SecondString=0;
880          return;
881        }
882        G4Parton *end  = hadron->GetNextParton();
883        if ( end==NULL)
884        { G4cout << " G4FTFModel::String() Error:No end parton found"<< G4endl;
885          FirstString=0; SecondString=0;
886          return;
887        }
888        G4LorentzVector Phadron=hadron->Get4Momentum();
889/*
890G4cout<<"Create strings had "<<hadron->GetDefinition()->GetParticleName()<<" "<<Phadron<<" "<<Phadron.mag()<<G4endl;
891G4cout<<"isProjectile "<<isProjectile<<G4endl;
892G4cout<<"start Q "<<start->GetDefinition()->GetPDGEncoding()<<G4endl;
893G4cout<<"end   Q "<<end->GetDefinition()->GetPDGEncoding()<<G4endl;
894*/
895        G4LorentzVector Pstart(0.,0.,0.,0.);
896        G4LorentzVector Pend(0.,0.,0.,0.);
897        G4LorentzVector Pkink(0.,0.,0.,0.);
898        G4LorentzVector PkinkQ1(0.,0.,0.,0.);
899        G4LorentzVector PkinkQ2(0.,0.,0.,0.);
900
901        G4int PDGcode_startQ = std::abs(start->GetDefinition()->GetPDGEncoding());
902        G4int PDGcode_endQ   = std::abs(  end->GetDefinition()->GetPDGEncoding());
903
904//--------------------------------------------------------------------------------       
905        G4double Wmin(0.);
906        if(isProjectile)
907        {
908          Wmin=theParameters->GetProjMinDiffMass();
909        } else
910        {
911          Wmin=theParameters->GetTarMinDiffMass();
912        } // end of if(isProjectile)
913
914        G4double W = hadron->Get4Momentum().mag();
915        G4double W2=W*W;
916//G4cout<<"W Wmin "<<W<<" "<<Wmin<<G4endl;
917        G4double Pt(0.), x1(0.), x2(0.), x3(0.);
918        G4bool Kink=false;
919
920        if(W > Wmin)
921        {                                        // Kink is possible
922          G4double Pt2kink=theParameters->GetPt2Kink();
923          Pt = std::sqrt(Pt2kink*(std::pow(W2/16./Pt2kink+1.,G4UniformRand()) - 1.));
924// Pt=0.;
925//G4cout<<"Pt2kink Pt Pt2 "<<Pt2kink<<" "<<Pt<<" "<<Pt*Pt<<G4endl;
926          if(Pt > 500.*MeV)
927          {
928             G4double Ymax = std::log(W/2./Pt + std::sqrt(W2/4./Pt/Pt - 1.));
929             G4double Y=Ymax*(1.- 2.*G4UniformRand());
930//G4cout<<"Ymax Y "<<Ymax<<" "<<Y<<G4endl;
931             x1=1.-Pt/W*std::exp( Y);
932             x3=1.-Pt/W*std::exp(-Y);
933             x2=2.-x1-x3;
934//G4cout<<"X1 X2 X3 "<<x1<<" "<<x2<<" "<<x3<<G4endl;
935             G4double Mass_startQ = 650.*MeV;
936             if(PDGcode_startQ <  3) Mass_startQ =  325.*MeV;
937             if(PDGcode_startQ == 3) Mass_startQ =  500.*MeV;
938             if(PDGcode_startQ == 4) Mass_startQ = 1600.*MeV;
939
940             G4double Mass_endQ = 650.*MeV;
941             if(PDGcode_endQ <  3) Mass_endQ =  325.*MeV;
942             if(PDGcode_endQ == 3) Mass_endQ =  500.*MeV;
943             if(PDGcode_endQ == 4) Mass_endQ = 1600.*MeV;
944
945             G4double P2_1=W2*x1*x1/4.-Mass_endQ  *Mass_endQ;
946             G4double P2_3=W2*x3*x3/4.-Mass_startQ*Mass_startQ;
947     
948             G4double P2_2=sqr((2.-x1-x3)*W/2.);
949//G4cout<<"P2_1 P2_2 P2_3 "<<P2_1<<" "<<P2_2<<" "<<P2_3<<G4endl;
950             if((P2_1 < 0.) || (P2_3 <0.))
951             { Kink=false;}
952             else
953             {
954               G4double P_1=std::sqrt(P2_1);
955               G4double P_2=std::sqrt(P2_2);
956               G4double P_3=std::sqrt(P2_3);
957//G4cout<<"P_1 P_2 P_3 "<<P_1<<" "<<P_2<<" "<<P_3<<G4endl;
958               G4double CosT12=(P2_3-P2_1-P2_2)/(2.*P_1*P_2);
959               G4double CosT13=(P2_2-P2_1-P2_3)/(2.*P_1*P_3);
960//G4cout<<"CosT12 CosT13 "<<CosT12<<" "<<CosT13<<" "<<std::acos(CosT12)*180./3.14159<<" "<<std::acos(CosT13)*180./3.14159<<G4endl;
961               if((std::abs(CosT12) >1.) || (std::abs(CosT13) > 1.)) 
962               { Kink=false;}
963               else
964               { 
965                 Kink=true; 
966                 Pstart.setPx(-Pt); Pstart.setPy(0.); Pstart.setPz(P_3*CosT13); 
967                 Pend.setPx(   0.); Pend.setPy(  0.); Pend.setPz(          P_1); 
968                 Pkink.setPx(  Pt); Pkink.setPy( 0.); Pkink.setPz(  P_2*CosT12);
969                 Pstart.setE(x3*W/2.);               
970                 Pkink.setE(Pkink.vect().mag());
971                 Pend.setE(x1*W/2.);
972
973                 G4double XkQ=GetQuarkFractionOfKink(0.,1.);
974                 if(Pkink.getZ() > 0.) 
975                 {
976                   if(XkQ > 0.5) {PkinkQ1=XkQ*Pkink;} else {PkinkQ1=(1.-XkQ)*Pkink;}
977                 } else {
978                   if(XkQ > 0.5) {PkinkQ1=(1.-XkQ)*Pkink;} else {PkinkQ1=XkQ*Pkink;}
979                 }
980
981                 PkinkQ2=Pkink - PkinkQ1;
982//------------------------- Minimizing Pt1^2+Pt3^2 ------------------------------
983/*
984G4cout<<"Pstart "<<Pstart<<G4endl;
985G4cout<<"Pkink  "<<Pkink <<G4endl;
986G4cout<<"Pkink1 "<<PkinkQ1<<G4endl;
987G4cout<<"Pkink2 "<<PkinkQ2<<G4endl;
988G4cout<<"Pend   "<<Pend  <<G4endl;
989*/
990                 G4double Cos2Psi=(sqr(x1) -sqr(x3)+2.*sqr(x3*CosT13))/
991                          std::sqrt(sqr(sqr(x1)-sqr(x3)) + sqr(2.*x1*x3*CosT13));
992                 G4double Psi=std::acos(Cos2Psi);
993//G4cout<<"Psi "<<Psi<<" "<<Psi*180./twopi<<G4endl;
994
995G4LorentzRotation Rotate;
996if(isProjectile) {Rotate.rotateY(Psi);}
997else             {Rotate.rotateY(pi-Psi);}                   
998Rotate.rotateZ(twopi*G4UniformRand());
999
1000//G4cout<<"Rotate"<<G4endl;
1001
1002Pstart*=Rotate;
1003Pkink*=Rotate;
1004PkinkQ1*=Rotate;
1005PkinkQ2*=Rotate;
1006Pend*=Rotate;
1007/*
1008G4cout<<"Pstart "<<Pstart<<G4endl;
1009G4cout<<"Pkink1 "<<PkinkQ1 <<G4endl;
1010G4cout<<"Pkink2 "<<PkinkQ2 <<G4endl;
1011G4cout<<"Pend   "<<Pend  <<G4endl;
1012*/
1013               }
1014             }      // end of if((P2_1 < 0.) || (P2_3 <0.))
1015          }         // end of if(Pt > 500.*MeV)
1016        }           // end of if(W > Wmin) Check for a kink
1017
1018//--------------------------------------------------------------------------------
1019
1020//G4cout<<"Kink "<<Kink<<G4endl;
1021
1022        if(Kink)
1023        {                                        // Kink is possible
1024          std::vector<G4double> QuarkProbabilitiesAtGluonSplitUp =
1025              theParameters->GetQuarkProbabilitiesAtGluonSplitUp();
1026
1027          G4int QuarkInGluon(1); G4double Ksi=G4UniformRand();
1028          for(unsigned int Iq=0; Iq <3; Iq++)
1029          {
1030//G4cout<<Iq<<" "<<QuarkProbabilitiesAtGluonSplitUp[Iq]<<G4endl;
1031if(Ksi > QuarkProbabilitiesAtGluonSplitUp[Iq]) QuarkInGluon++;}
1032
1033//G4cout<<"Gquark "<<QuarkInGluon<<G4endl;
1034          G4Parton * Gquark = new G4Parton(QuarkInGluon);
1035          G4Parton * Ganti_quark = new G4Parton(-QuarkInGluon);
1036/*
1037G4cout<<Gquark->GetDefinition()->GetParticleName()<<" "<<Gquark->GetDefinition()->GetPDGEncoding()<<" "<<Gquark->GetDefinition()->GetPDGMass()<<G4endl;
1038G4cout<<Ganti_quark->GetDefinition()->GetParticleName()<<" "<<Ganti_quark->GetDefinition()->GetPDGEncoding()<<" "<<Ganti_quark->GetDefinition()->GetPDGMass()<<G4endl;
1039*/
1040
1041//G4int Uzhi; G4cin>>Uzhi;
1042
1043//-------------------------------------------------------------------------------
1044/*
1045DefineMomentumInZ(G4double aLightConeMomentum, G4bool aDirection);
1046      void Set4Momentum(const G4LorentzVector & aMomentum);
1047      G4int PDGencoding;
1048      G4ParticleDefinition * theDefinition;
1049      G4LorentzVector theMomentum;
1050      G4ThreeVector   thePosition;
1051     
1052      G4int theColour;
1053      G4double theIsoSpinZ;
1054      G4double theSpinZ;
1055     
1056      G4double theX;
1057*/
1058//-------------------------------------------------------------------------------
1059
1060//G4cout<<"Phadron "<<Phadron<<" mass "<<Phadron.mag()<<G4endl;
1061          G4LorentzRotation toCMS(-1*Phadron.boostVector());
1062//G4cout<<"To lab"<<G4endl;
1063          G4LorentzRotation toLab(toCMS.inverse());
1064
1065          Pstart.transform(toLab);  start->Set4Momentum(Pstart);
1066          PkinkQ1.transform(toLab);
1067          PkinkQ2.transform(toLab);
1068          Pend.transform(toLab);    end->Set4Momentum(Pend);
1069/*
1070G4cout<<"Pstart "<<start->GetDefinition()->GetPDGEncoding()<<Pstart<<G4endl;
1071G4cout<<"Pkink1 "<<PkinkQ1 <<G4endl;
1072G4cout<<"Pkink2 "<<PkinkQ2 <<G4endl;
1073G4cout<<"Pend   "<<end->GetDefinition()->GetPDGEncoding()<<Pend  <<G4endl;
1074*/
1075// !!!    G4ExcitedString * FirstString(0); G4ExcitedString * SecondString(0);
1076          G4int absPDGcode=std::abs(hadron->GetDefinition()->GetPDGEncoding());
1077
1078//G4cout<<"absPDGcode "<<absPDGcode<<G4endl;
1079
1080          if(absPDGcode < 1000)
1081          {                                // meson
1082            if ( isProjectile )
1083            {                              // Projectile
1084              if(end->GetDefinition()->GetPDGEncoding() > 0 )  // A quark on the end
1085              {                            // Quark on the end
1086                FirstString = new G4ExcitedString(end   ,Ganti_quark, +1);
1087                SecondString= new G4ExcitedString(Gquark,start      ,+1);
1088                Ganti_quark->Set4Momentum(PkinkQ1);
1089                Gquark->Set4Momentum(PkinkQ2);
1090/*
1091G4cout<<" Proj Meson end Q"<<G4endl;
1092G4cout<<"First string  ============ "<<G4endl;
1093G4cout<<"end  Q "<<end->GetDefinition()->GetPDGEncoding()<<" "<<end->Get4Momentum()<<G4endl;
1094G4cout<<"G antiQ"<<Ganti_quark->GetDefinition()->GetPDGEncoding()<<" "<<Ganti_quark->Get4Momentum()<<G4endl;
1095G4cout<<"Sum P  "<<(Ganti_quark->Get4Momentum()+end->Get4Momentum())<<G4endl;
1096G4cout<<"Secondstring   ============ "<<G4endl;
1097G4cout<<"G    Q "<<Gquark->GetDefinition()->GetPDGEncoding()<<" "<<Gquark->Get4Momentum()<<G4endl;
1098G4cout<<"startQ "<<start->GetDefinition()->GetPDGEncoding()<<" "<<start->Get4Momentum()<<G4endl;
1099
1100G4cout<<"Sum P  "<<(Gquark->Get4Momentum()+start->Get4Momentum())<<G4endl;
1101*/
1102              } else
1103              {                            // Anti_Quark on the end
1104                FirstString = new G4ExcitedString(end        ,Gquark, +1);
1105                SecondString= new G4ExcitedString(Ganti_quark,start ,+1);
1106                Gquark->Set4Momentum(PkinkQ1);
1107                Ganti_quark->Set4Momentum(PkinkQ2);
1108/*
1109G4cout<<" Proj Meson end Qbar"<<G4endl;
1110G4cout<<"First string  ============ "<<G4endl;
1111G4cout<<"end  Q "<<end->GetDefinition()->GetPDGEncoding()<<" "<<end->Get4Momentum()<<G4endl;
1112G4cout<<"G     Q"<<Gquark->GetDefinition()->GetPDGEncoding()<<" "<<Gquark->Get4Momentum()<<G4endl;
1113G4cout<<"Sum P  "<<(Gquark->Get4Momentum()+end->Get4Momentum())<<G4endl;
1114G4cout<<"Secondstring   ============ "<<G4endl;
1115G4cout<<"G antQ "<<Ganti_quark->GetDefinition()->GetPDGEncoding()<<" "<<Ganti_quark->Get4Momentum()<<G4endl;
1116G4cout<<"startQ "<<start->GetDefinition()->GetPDGEncoding()<<" "<<start->Get4Momentum()<<G4endl;
1117G4cout<<"Sum P  "<<(Ganti_quark->Get4Momentum()+start->Get4Momentum())<<G4endl;
1118*/
1119              }   // end of if(end->GetPDGcode() > 0)
1120            } else {                      // Target
1121              if(end->GetDefinition()->GetPDGEncoding() > 0 )  // A quark on the end
1122              {                           // Quark on the end
1123                FirstString = new G4ExcitedString(Ganti_quark,end   ,-1);
1124                SecondString= new G4ExcitedString(start      ,Gquark,-1);
1125                Ganti_quark->Set4Momentum(PkinkQ2);
1126                Gquark->Set4Momentum(PkinkQ1);
1127/*
1128G4cout<<" Targ Meson end Q"<<G4endl;
1129G4cout<<"First string   ============ "<<G4endl;
1130G4cout<<"G antiQ"<<Ganti_quark->GetDefinition()->GetPDGEncoding()<<" "<<Ganti_quark->Get4Momentum()<<G4endl;
1131G4cout<<"end  Q "<<end->GetDefinition()->GetPDGEncoding()<<" "<<end->Get4Momentum()<<G4endl;
1132G4cout<<"Sum P  "<<(Ganti_quark->Get4Momentum()+end->Get4Momentum())<<G4endl;
1133G4cout<<"Secondstring   ============ "<<G4endl;
1134G4cout<<"startQ "<<start->GetDefinition()->GetPDGEncoding()<<" "<<start->Get4Momentum()<<G4endl;
1135G4cout<<"G    Q "<<Gquark->GetDefinition()->GetPDGEncoding()<<" "<<Gquark->Get4Momentum()<<G4endl;
1136G4cout<<"Sum P  "<<(Gquark->Get4Momentum()+start->Get4Momentum())<<G4endl;
1137*/
1138              } else
1139              {                            // Anti_Quark on the end
1140                FirstString = new G4ExcitedString(Gquark,end        ,-1);
1141                SecondString= new G4ExcitedString(start ,Ganti_quark,-1);
1142                Gquark->Set4Momentum(PkinkQ2);
1143                Ganti_quark->Set4Momentum(PkinkQ1);
1144/*
1145G4cout<<" Targ Meson end Qbar"<<G4endl;
1146G4cout<<"First string   ============ "<<G4endl;
1147G4cout<<"G     Q"<<Gquark->GetDefinition()->GetPDGEncoding()<<" "<<Gquark->Get4Momentum()<<G4endl;
1148G4cout<<"end  Q "<<end->GetDefinition()->GetPDGEncoding()<<" "<<end->Get4Momentum()<<G4endl;
1149G4cout<<"Sum P  "<<(Gquark->Get4Momentum()+end->Get4Momentum())<<G4endl;
1150G4cout<<"Secondstring   ============ "<<G4endl;
1151G4cout<<"startQ "<<start->GetDefinition()->GetPDGEncoding()<<" "<<start->Get4Momentum()<<G4endl;
1152G4cout<<"G antQ "<<Ganti_quark->GetDefinition()->GetPDGEncoding()<<" "<<Ganti_quark->Get4Momentum()<<G4endl;
1153G4cout<<"Sum P  "<<(Ganti_quark->Get4Momentum()+start->Get4Momentum())<<G4endl;
1154*/
1155              }   // end of if(end->GetPDGcode() > 0)
1156            }     // end of if ( isProjectile )
1157          } else  // if(absPDGCode < 1000)
1158          {                             // Baryon/AntiBaryon
1159            if ( isProjectile )
1160            {                              // Projectile
1161              if((end->GetDefinition()->GetParticleType() == "diquarks") &&
1162                 (end->GetDefinition()->GetPDGEncoding() > 0           )   ) 
1163              {                            // DiQuark on the end
1164                FirstString = new G4ExcitedString(end        ,Gquark, +1);
1165                SecondString= new G4ExcitedString(Ganti_quark,start ,+1);
1166                Gquark->Set4Momentum(PkinkQ1);
1167                Ganti_quark->Set4Momentum(PkinkQ2);
1168/*
1169G4cout<<" Proj baryon end QQ"<<G4endl;
1170G4cout<<"First string   ============ "<<G4endl;
1171G4cout<<"end OQ "<<end->GetDefinition()->GetPDGEncoding()<<" "<<end->Get4Momentum()<<G4endl;
1172G4cout<<"G    Q"<<Gquark->GetDefinition()->GetPDGEncoding()<<" "<<Gquark->Get4Momentum()<<G4endl;
1173G4cout<<"Sum P  "<<(Gquark->Get4Momentum()+end->Get4Momentum())<<G4endl;
1174G4cout<<"Secondstring   ============ "<<G4endl;
1175G4cout<<"G  Qbar"<<Ganti_quark->GetDefinition()->GetPDGEncoding()<<" "<<Ganti_quark->Get4Momentum()<<G4endl;
1176G4cout<<"startQ "<<start->GetDefinition()->GetPDGEncoding()<<" "<<start->Get4Momentum()<<G4endl;
1177G4cout<<"Sum P  "<<(Ganti_quark->Get4Momentum()+start->Get4Momentum())<<G4endl;
1178*/
1179              } else
1180              {                            // Anti_DiQuark on the end or quark
1181                FirstString = new G4ExcitedString(end   ,Ganti_quark, +1);
1182                SecondString= new G4ExcitedString(Gquark,start      ,+1);
1183                Ganti_quark->Set4Momentum(PkinkQ1);
1184                Gquark->Set4Momentum(PkinkQ2);
1185/*
1186G4cout<<" Proj baryon end Q"<<G4endl;
1187G4cout<<"First string   ============ "<<G4endl;
1188G4cout<<"end OQ "<<end->GetDefinition()->GetPDGEncoding()<<" "<<end->Get4Momentum()<<G4endl;
1189G4cout<<"G antQ"<<Ganti_quark->GetDefinition()->GetPDGEncoding()<<" "<<Ganti_quark->Get4Momentum()<<G4endl;
1190G4cout<<"Sum P  "<<(Ganti_quark->Get4Momentum()+end->Get4Momentum())<<G4endl;
1191G4cout<<"Secondstring   ============ "<<G4endl;
1192G4cout<<"G  Q   "<<Gquark->GetDefinition()->GetPDGEncoding()<<" "<<Gquark->Get4Momentum()<<G4endl;
1193G4cout<<"startQ "<<start->GetDefinition()->GetPDGEncoding()<<" "<<start->Get4Momentum()<<G4endl;
1194G4cout<<"Sum P  "<<(Gquark->Get4Momentum()+start->Get4Momentum())<<G4endl;
1195*/
1196              }   // end of if(end->GetPDGcode() > 0)
1197            } else {                      // Target
1198
1199              if((end->GetDefinition()->GetParticleType() == "diquarks") &&
1200                 (end->GetDefinition()->GetPDGEncoding() > 0           )   ) 
1201              {                            // DiQuark on the end
1202                FirstString = new G4ExcitedString(Gquark,end        ,-1);
1203
1204                SecondString= new G4ExcitedString(start ,Ganti_quark,-1);
1205                Gquark->Set4Momentum(PkinkQ1);
1206                Ganti_quark->Set4Momentum(PkinkQ2);
1207/*
1208G4cout<<" Targ baryon end QQ"<<G4endl;
1209G4cout<<"First string   ============ "<<G4endl;
1210G4cout<<"G  Q   "<<Gquark->GetDefinition()->GetPDGEncoding()<<" "<<Gquark->Get4Momentum()<<G4endl;
1211G4cout<<"end OQ "<<end->GetDefinition()->GetPDGEncoding()<<" "<<end->Get4Momentum()<<G4endl;
1212G4cout<<"Sum P  "<<(Gquark->Get4Momentum()+end->Get4Momentum())<<G4endl;
1213G4cout<<"Secondstring   ============ "<<G4endl;
1214G4cout<<"startQ "<<start->GetDefinition()->GetPDGEncoding()<<" "<<start->Get4Momentum()<<G4endl;
1215G4cout<<"G  Qbar"<<Ganti_quark->GetDefinition()->GetPDGEncoding()<<" "<<Ganti_quark->Get4Momentum()<<G4endl;
1216G4cout<<"Sum P  "<<(Ganti_quark->Get4Momentum()+start->Get4Momentum())<<G4endl;
1217*/
1218              } else
1219              {                            // Anti_DiQuark on the end or Q
1220                FirstString = new G4ExcitedString(Ganti_quark,end   ,-1);
1221                SecondString= new G4ExcitedString(start      ,Gquark,-1);
1222                Gquark->Set4Momentum(PkinkQ2);
1223                Ganti_quark->Set4Momentum(PkinkQ1);
1224/*
1225G4cout<<" Targ baryon end Q"<<G4endl;
1226G4cout<<"First string   ============ "<<G4endl;
1227G4cout<<"G  Qbar"<<Ganti_quark->GetDefinition()->GetPDGEncoding()<<" "<<Ganti_quark->Get4Momentum()<<G4endl;
1228G4cout<<"end O "<<end->GetDefinition()->GetPDGEncoding()<<" "<<end->Get4Momentum()<<G4endl;
1229G4cout<<"Sum P  "<<(Ganti_quark->Get4Momentum()+end->Get4Momentum())<<G4endl;
1230G4cout<<"Secondstring   ============ "<<G4endl;
1231G4cout<<"startQ "<<start->GetDefinition()->GetPDGEncoding()<<" "<<start->Get4Momentum()<<G4endl;
1232G4cout<<"G  Q   "<<Gquark->GetDefinition()->GetPDGEncoding()<<" "<<Gquark->Get4Momentum()<<G4endl;
1233G4cout<<"Sum P  "<<(Gquark->Get4Momentum()+start->Get4Momentum())<<G4endl;
1234*/
1235              }   // end of if(end->GetPDGcode() > 0)
1236            }     // end of if ( isProjectile )
1237          }  // end of if(absPDGcode < 1000)
1238
1239          FirstString->SetTimeOfCreation(hadron->GetTimeOfCreation());
1240          FirstString->SetPosition(hadron->GetPosition());
1241
1242          SecondString->SetTimeOfCreation(hadron->GetTimeOfCreation());
1243          SecondString->SetPosition(hadron->GetPosition());
1244
1245// ------------------------------------------------------------------------- 
1246        } else                                   // End of kink is possible
1247        {                                        // Kink is impossible
1248          if ( isProjectile )
1249          {
1250                FirstString= new G4ExcitedString(end,start, +1);
1251          } else {
1252                FirstString= new G4ExcitedString(start,end, -1);
1253          }
1254
1255          SecondString=0;
1256
1257          FirstString->SetTimeOfCreation(hadron->GetTimeOfCreation());
1258          FirstString->SetPosition(hadron->GetPosition());
1259
1260// momenta of string ends
1261//
1262          G4double Momentum=hadron->Get4Momentum().vect().mag();
1263          G4double  Plus=hadron->Get4Momentum().e() + Momentum;
1264          G4double Minus=hadron->Get4Momentum().e() - Momentum;
1265
1266          G4ThreeVector tmp;
1267          if(Momentum > 0.) 
1268          {
1269           tmp.set(hadron->Get4Momentum().px(),
1270                   hadron->Get4Momentum().py(),
1271                   hadron->Get4Momentum().pz());
1272           tmp/=Momentum;
1273          }
1274          else
1275          {
1276           tmp.set(0.,0.,1.);
1277          }
1278
1279          G4LorentzVector Pstart(tmp,0.);
1280          G4LorentzVector   Pend(tmp,0.);
1281
1282          if(isProjectile)
1283          {
1284           Pstart*=(-1.)*Minus/2.;
1285           Pend  *=(+1.)*Plus /2.;
1286          } 
1287          else
1288          {
1289           Pstart*=(+1.)*Plus/2.;
1290           Pend  *=(-1.)*Minus/2.;
1291          }
1292
1293          Momentum=-Pstart.mag();
1294          Pstart.setT(Momentum);  // It is assumed that quark has m=0.
1295
1296          Momentum=-Pend.mag();
1297          Pend.setT(Momentum);    // It is assumed that di-quark has m=0.
1298
1299          start->Set4Momentum(Pstart);
1300          end->Set4Momentum(Pend);
1301          SecondString=0;
1302        }            // End of kink is impossible
1303
1304#ifdef G4_FTFDEBUG
1305          G4cout << " generated string flavors          " 
1306                 << start->GetPDGcode() << " / " 
1307                 << end->GetPDGcode()                    << G4endl;
1308          G4cout << " generated string momenta:   quark " 
1309                 << start->Get4Momentum() << "mass : " 
1310                 <<start->Get4Momentum().mag()           << G4endl;
1311          G4cout << " generated string momenta: Diquark " 
1312                 << end ->Get4Momentum() 
1313                 << "mass : " <<end->Get4Momentum().mag()<< G4endl;
1314          G4cout << " sum of ends                       " << Pstart+Pend << G4endl;
1315          G4cout << " Original                          " << hadron->Get4Momentum() << G4endl;
1316#endif
1317
1318        return;
1319   
1320}
1321
1322
1323// --------- private methods ----------------------
1324
1325// ---------------------------------------------------------------------
1326G4double G4DiffractiveExcitation::ChooseP(G4double Pmin, G4double Pmax) const
1327{
1328// choose an x between Xmin and Xmax with P(x) ~ 1/x
1329//  to be improved...
1330
1331        G4double range=Pmax-Pmin;                   
1332
1333        if ( Pmin <= 0. || range <=0. )
1334        {
1335                G4cout << " Pmin, range : " << Pmin << " , " << range << G4endl;
1336                throw G4HadronicException(__FILE__, __LINE__, "G4DiffractiveExcitation::ChooseP : Invalid arguments ");
1337        }
1338
1339        G4double P=Pmin * std::pow(Pmax/Pmin,G4UniformRand()); 
1340        return P;
1341}
1342
1343// ---------------------------------------------------------------------
1344G4ThreeVector G4DiffractiveExcitation::GaussianPt(G4double AveragePt2, 
1345                                                  G4double maxPtSquare) const
1346{            //  @@ this method is used in FTFModel as well. Should go somewhere common!
1347
1348        G4double Pt2;
1349        Pt2 = -AveragePt2 * std::log(1. + G4UniformRand() * 
1350                           (std::exp(-maxPtSquare/AveragePt2)-1.));
1351
1352        G4double Pt=std::sqrt(Pt2);
1353        G4double phi=G4UniformRand() * twopi;
1354        return G4ThreeVector (Pt*std::cos(phi), Pt*std::sin(phi), 0.);
1355}
1356
1357// ---------------------------------------------------------------------
1358G4double G4DiffractiveExcitation::GetQuarkFractionOfKink(G4double zmin, G4double zmax) const
1359    {
1360       G4double z, yf;
1361       do {
1362           z = zmin + G4UniformRand()*(zmax-zmin);
1363           yf = z*z +sqr(1 - z);       
1364           } 
1365       while (G4UniformRand() > yf); 
1366       return z;
1367    }
1368// ---------------------------------------------------------------------
1369void G4DiffractiveExcitation::UnpackMeson(const G4int IdPDG, G4int &Q1, G4int &Q2) const // Uzhi 7.09.09
1370    {
1371       G4int absIdPDG = std::abs(IdPDG);
1372       Q1=  absIdPDG/ 100;
1373       Q2= (absIdPDG %100)/10;
1374           
1375       G4int anti= 1 -2 * ( std::max( Q1, Q2 ) % 2 );
1376
1377       if (IdPDG < 0 ) anti *=-1;
1378       Q1 *= anti;
1379       Q2 *= -1 * anti;
1380       return;
1381    }
1382// ---------------------------------------------------------------------
1383void G4DiffractiveExcitation::UnpackBaryon(G4int IdPDG, 
1384                              G4int &Q1, G4int &Q2, G4int &Q3) const // Uzhi 7.09.09
1385    {
1386       Q1 =  IdPDG           / 1000;
1387       Q2 = (IdPDG % 1000)  / 100;
1388       Q3 = (IdPDG % 100)   / 10;
1389       return;
1390    }
1391// ---------------------------------------------------------------------
1392G4int G4DiffractiveExcitation::NewNucleonId(G4int Q1, G4int Q2, G4int Q3) const // Uzhi 7.09.09
1393    {
1394       G4int TmpQ(0);
1395       if( Q3 > Q2 ) 
1396       {
1397        TmpQ = Q2;
1398        Q2 = Q3;
1399        Q3 = TmpQ;
1400       } else if( Q3 > Q1 )
1401       {
1402        TmpQ = Q1;
1403        Q1 = Q3;
1404        Q3 = TmpQ;
1405       }
1406
1407       if( Q2 > Q1 ) 
1408       {
1409        TmpQ = Q1;
1410        Q1 = Q2;
1411        Q2 = TmpQ;
1412       }
1413
1414       G4int NewCode = Q1*1000 + Q2* 100 + Q3*  10 + 2; 
1415       return NewCode;
1416    }
1417// ---------------------------------------------------------------------
1418G4DiffractiveExcitation::G4DiffractiveExcitation(const G4DiffractiveExcitation &)
1419{
1420        throw G4HadronicException(__FILE__, __LINE__, "G4DiffractiveExcitation copy contructor not meant to be called");
1421}
1422
1423
1424G4DiffractiveExcitation::~G4DiffractiveExcitation()
1425{
1426}
1427
1428
1429const G4DiffractiveExcitation & G4DiffractiveExcitation::operator=(const G4DiffractiveExcitation &)
1430{
1431        throw G4HadronicException(__FILE__, __LINE__, "G4DiffractiveExcitation = operator meant to be called");
1432        return *this;
1433}
1434
1435
1436int G4DiffractiveExcitation::operator==(const G4DiffractiveExcitation &) const
1437{
1438        throw G4HadronicException(__FILE__, __LINE__, "G4DiffractiveExcitation == operator meant to be called");
1439        return false;
1440}
1441
1442int G4DiffractiveExcitation::operator!=(const G4DiffractiveExcitation &) const
1443{
1444        throw G4HadronicException(__FILE__, __LINE__, "G4DiffractiveExcitation != operator meant to be called");
1445        return true;
1446}
Note: See TracBrowser for help on using the repository browser.