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

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

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