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

Last change on this file since 1201 was 1196, checked in by garnier, 16 years ago

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