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

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