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

Last change on this file since 1347 was 1347, checked in by garnier, 14 years ago

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