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

Last change on this file since 1357 was 1347, checked in by garnier, 15 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.