| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 #include <utility>
27
28 #include "G4QGSParticipants.hh"
29 #include "globals.hh"
30 #include "G4SystemOfUnits.hh"
31 #include "G4LorentzVector.hh"
32 #include "G4V3DNucleus.hh"
33 #include "G4ParticleTable.hh"
34 #include "G4IonTable.hh"
35 #include "G4PhysicalConstants.hh"
36
37 #include "G4Exp.hh"
38 #include "G4Log.hh"
39 #include "G4Pow.hh"
40
41 //#define debugQGSParticipants
42 //#define debugPutOnMassShell
43
44 // Class G4QGSParticipants
45
46 // Promoting model parameters from local varia
47 G4ThreadLocal G4int G4QGSParticipants_NPart =
48
49 G4QGSParticipants::G4QGSParticipants()
50 : theDiffExcitaton(), ModelMode(SOFT), nCutM
51 ThresholdParameter(0.0*GeV), QGSMThreshold
52 theNucleonRadius(1.5*fermi), theCurrentVel
53 theProjectileSplitable(nullptr), Regge(nul
54 InteractionMode(ALL), alpha(-0.5), beta(2.
55 NumberOfInvolvedNucleonsOfTarget(0), Numbe
56 ProjectileResidual4Momentum(G4LorentzVecto
57 ProjectileResidualCharge(0), ProjectileRes
58 TargetResidual4Momentum(G4LorentzVector())
59 TargetResidualCharge(0), TargetResidualExc
60 CofNuclearDestruction(0.0), R2ofNuclearDes
61 ExcitationEnergyPerWoundedNucleon(0.0), Do
62 Pt2ofNuclearDestruction(0.0), MaxPt2ofNucl
63 {
64 for (size_t i=0; i < 250; ++i) {
65 TheInvolvedNucleonsOfTarget[i] = nullptr;
66 TheInvolvedNucleonsOfProjectile[i] = nullp
67 }
68 // Parameters setting
69 SetCofNuclearDestruction( 1.0 );
70 SetR2ofNuclearDestruction( 1.5*fermi*fermi )
71 SetDofNuclearDestruction( 0.3 );
72 SetPt2ofNuclearDestruction( 0.075*GeV*GeV );
73 SetMaxPt2ofNuclearDestruction( 1.0*GeV*GeV )
74 SetExcitationEnergyPerWoundedNucleon( 40.0*M
75
76 sigmaPt=0.25*sqr(GeV);
77 }
78
79 G4QGSParticipants::G4QGSParticipants(const G4Q
80 : G4VParticipants(), ModelMode(right.ModelMo
81 ThresholdParameter(right.ThresholdParamete
82 QGSMThreshold(right.QGSMThreshold),
83 theNucleonRadius(right.theNucleonRadius),
84 theCurrentVelocity(right.theCurrentVelocit
85 theProjectileSplitable(right.theProjectile
86 Regge(right.Regge), InteractionMode(right.
87 alpha(right.alpha), beta(right.beta), sigm
88 NumberOfInvolvedNucleonsOfTarget(right.Num
89 NumberOfInvolvedNucleonsOfProjectile(right
90 ProjectileResidual4Momentum(right.Projecti
91 ProjectileResidualMassNumber(right.Project
92 ProjectileResidualCharge(right.ProjectileR
93 ProjectileResidualExcitationEnergy(right.P
94 TargetResidual4Momentum(right.TargetResidu
95 TargetResidualMassNumber(right.TargetResid
96 TargetResidualCharge(right.TargetResidualC
97 TargetResidualExcitationEnergy(right.Targe
98 CofNuclearDestruction(right.CofNuclearDest
99 R2ofNuclearDestruction(right.R2ofNuclearDe
100 ExcitationEnergyPerWoundedNucleon(right.Ex
101 DofNuclearDestruction(right.DofNuclearDest
102 Pt2ofNuclearDestruction(right.Pt2ofNuclear
103 MaxPt2ofNuclearDestruction(right.MaxPt2ofN
104 {
105 for (size_t i=0; i < 250; ++i) {
106 TheInvolvedNucleonsOfTarget[i] = right.The
107 TheInvolvedNucleonsOfProjectile[i] = right
108 }
109 }
110
111 G4QGSParticipants::~G4QGSParticipants() {}
112
113 void G4QGSParticipants::BuildInteractions(cons
114 {
115 theProjectile = thePrimary;
116
117 Regge = new G4Reggeons(theProjectile.GetDefi
118
119 SetProjectileNucleus( 0 );
120
121 NumberOfInvolvedNucleonsOfProjectile= 0;
122 G4LorentzVector tmp( 0.0, 0.0, 0.0, 0.0 );
123 ProjectileResidualMassNumber = 0;
124 ProjectileResidualCharge = 0;
125 ProjectileResidualExcitationEnergy = 0.0;
126 ProjectileResidual4Momentum = tmp;
127
128 NumberOfInvolvedNucleonsOfTarget= 0;
129 TargetResidualMassNumber = theNucleus-
130 TargetResidualCharge = theNucleus-
131 TargetResidualExcitationEnergy = 0.0;
132
133 theNucleus->StartLoop();
134 G4Nucleon* NuclearNucleon;
135 while ( ( NuclearNucleon = theNucleus->GetNe
136 tmp+=NuclearNucleon->Get4Momentum();
137 }
138 TargetResidual4Momentum = tmp;
139
140 if ( std::abs( theProjectile.GetDefinition()
141 // Projectile is a hadron : meson or baryo
142 ProjectileResidualMassNumber = std::abs( t
143 ProjectileResidualCharge = G4int( theProje
144 ProjectileResidualExcitationEnergy = 0.0;
145 ProjectileResidual4Momentum.setVect( thePr
146 ProjectileResidual4Momentum.setE( theProje
147 }
148
149 #ifdef debugQGSParticipants
150 G4cout <<G4endl<< "G4QGSParticipants::Buil
151 << "thePrimary " << thePrimary.GetD
152 <<ProjectileResidual4Momentum<<G4en
153 G4cout << "Target A and Z " << theNucleus
154 << TargetResidual4Momentum<<G4endl;
155 #endif
156
157 G4bool Success( true );
158
159 const G4int maxNumberOfLoops = 1000;
160 G4int loopCounter = 0;
161 do
162 {
163 const G4int maxNumberOfInternalLoops = 100
164 G4int internalLoopCounter = 0;
165 do
166 {
167 if(std::abs(theProjectile.GetDefinition(
168 {
169 SelectInteractions(theProjectile); //
170 }
171 else
172 {
173 GetList( theProjectile ); // Get list
174 }
175
176 if ( theInteractions.size() == 0 ) retur
177
178 StoreInvolvedNucleon(); // Store p
179
180 ReggeonCascade(); // Make re
181
182 Success = PutOnMassShell(); // Ascribe
183
184 if(!Success) PrepareInitialState( thePri
185
186 } while( (!Success) && ++internalLoopCount
187
188 if ( internalLoopCounter >= maxNumberOfInt
189 Success = false;
190 }
191
192 if ( Success ) {
193 #ifdef debugQGSParticipants
194 G4cout<<G4endl<<"PerformDiffractiveCol
195 #endif
196
197 PerformDiffractiveCollisions();
198
199 #ifdef debugQGSParticipants
200 G4cout<<G4endl<<"SplitHadrons();" <<G4
201 #endif
202
203 SplitHadrons();
204
205 if( theProjectileSplitable && theProject
206 #ifdef debugQGSParticipants
207 G4cout<<"Perform non-Diffractive Col
208 #endif
209 Success = DeterminePartonMomenta();
210 }
211
212 if(!Success) PrepareInitialState( thePri
213 }
214 } while( (!Success) && ++loopCounter < maxNu
215
216 if ( loopCounter >= maxNumberOfLoops ) {
217 Success = false;
218 #ifdef debugQGSParticipants
219 G4cout<<"NOT Successful ======" <<G4endl
220 #endif
221 }
222
223 if ( Success ) {
224 CreateStrings(); // To create strings
225
226 GetResiduals(); // To calculate residual
227
228 #ifdef debugQGSParticipants
229 G4cout<<"All O.K. ======" <<G4endl;
230 #endif
231 }
232
233 // clean-up, if necessary
234 #ifdef debugQGSParticipants
235 G4cout<<"Clearing "<<G4endl;
236 G4cout<<"theInteractions.size() "<<theInte
237 #endif
238
239 if( Regge ) delete Regge;
240
241 std::for_each( theInteractions.begin(), theI
242 theInteractions.clear();
243
244 // Erasing of target involved nucleons.
245 #ifdef debugQGSParticipants
246 G4cout<<"Erasing of involved target nucleo
247 #endif
248
249 if ( NumberOfInvolvedNucleonsOfTarget != 0 )
250 for ( G4int i = 0; i < NumberOfInvolvedNu
251 G4VSplitableHadron* aNucleon = TheInvolv
252 if ( (aNucleon != 0 ) && (aNucleon->GetS
253 }
254 }
255
256 // Erasing of projectile involved nucleons.
257 if ( NumberOfInvolvedNucleonsOfProjectile !=
258 for ( G4int i = 0; i < NumberOfInvolvedNu
259 G4VSplitableHadron* aNucleon = TheInvol
260 if ( aNucleon ) delete aNucleon;
261 }
262 }
263
264 #ifdef debugQGSParticipants
265 G4cout<<"Delition of target nucleons from
266 <<G4endl<<G4endl;
267 #endif
268 std::for_each(theTargets.begin(), theTargets
269 theTargets.clear();
270
271 if ( theProjectileSplitable ) {
272 delete theProjectileSplitable;
273 theProjectileSplitable = 0;
274 }
275 }
276
277 //============================================
278 void G4QGSParticipants::PrepareInitialState( c
279 {
280 // Clearing of the arrays
281 // Erasing of the projectile
282 G4InteractionContent* anIniteraction = theIn
283 G4VSplitableHadron* pProjectile = anIniterac
284 if( pProjectile ) delete pProjectile;
285
286 std::for_each(theInteractions.begin(), theIn
287 theInteractions.clear();
288
289 // Erasing of the envolved nucleons and targ
290 theNucleus->StartLoop();
291 G4Nucleon* aNucleon;
292 while ( ( aNucleon = theNucleus->GetNextNucl
293 {
294 if ( aNucleon->AreYouHit() ) {
295 G4VSplitableHadron* splaNucleon = aNucle
296 if ( (splaNucleon != 0) && (splaNucleon-
297 aNucleon->Hit(nullptr);
298 NumberOfInvolvedNucleonsOfTarget--;
299 }
300 }
301
302 // Erasing of nuclear nucleons participated
303 std::for_each(theTargets.begin(), theTargets
304 theTargets.clear();
305
306 // Preparation to a new attempt
307 theProjectile = thePrimary;
308
309 theNucleus->Init(theNucleus->GetMassNumber()
310 theNucleus->SortNucleonsIncZ();
311 DoLorentzBoost(-theCurrentVelocity); // Lor
312
313 if (theNucleus->GetMassNumber() == 1)
314 {
315 G4ThreeVector aPos = G4ThreeVector(0.,0.,0
316 theNucleus->StartLoop();
317 G4Nucleon* tNucleon=theNucleus->GetNextNuc
318 tNucleon->SetPosition(aPos);
319 }
320
321 G4LorentzVector Tmp( 0.0, 0.0, 0.0, 0.0 );
322 NumberOfInvolvedNucleonsOfTarget= 0;
323 TargetResidualMassNumber = theNucleus-
324 TargetResidualCharge = theNucleus-
325 TargetResidualExcitationEnergy = 0.0;
326
327 G4Nucleon* NuclearNucleon;
328 while ( ( NuclearNucleon = theNucleus->GetNe
329 {Tmp+=NuclearNucleon->Get4Moment
330
331 TargetResidual4Momentum = Tmp;
332 }
333
334 //============================================
335 void G4QGSParticipants::GetList( const G4React
336 #ifdef debugQGSParticipants
337 G4cout<<G4endl<<"G4QGSParticipants::GetLis
338 #endif
339
340 // Direction: True - Proj, False - Target
341 theProjectileSplitable = new G4QGSMSplitable
342 theProjectileSplitable->SetStatus(1);
343
344 G4LorentzVector aPrimaryMomentum(thePrimary.
345 G4LorentzVector aNucleonMomentum(0.,0.,0., 9
346
347 G4double SS=(aPrimaryMomentum + aNucleonMome
348
349 Regge->SetS(SS);
350
351 //--------------------------------------
352 theNucleus->StartLoop();
353 G4Nucleon * tNucleon = theNucleus->GetNextNu
354
355 if ( ! tNucleon ) {
356 #ifdef debugQGSParticipants
357 G4cout << "QGSM - BAD situation: pNucleo
358 #endif
359 return;
360 }
361
362 G4double theNucleusOuterR = theNucleus->GetO
363
364 if (theNucleus->GetMassNumber() == 1)
365 {
366 G4ThreeVector aPos = G4ThreeVector(0.,0.,0
367 tNucleon->SetPosition(aPos);
368 theNucleusOuterR = 0.;
369 }
370
371 // Determination of participating nucleons o
372
373 std::for_each(theInteractions.begin(), theIn
374 theInteractions.clear();
375
376 G4int MaxPower=thePrimary.GetMomentum().mag(
377
378 const G4int maxNumberOfLoops = 1000;
379
380 G4int NumberOfTries = 0;
381 G4double Scale = 1.0;
382
383 G4int loopCounter = -1;
384 while( (theInteractions.size() == 0) && ++lo
385 {
386 InteractionMode = ALL; // Mode = ALL, WIT
387
388 // choose random impact parameter of a col
389 std::pair<G4double, G4double> theImpactPar
390
391 NumberOfTries++;
392 if( NumberOfTries == 100*(NumberOfTries/10
393
394 theImpactParameter = theNucleus->ChooseImp
395 G4double impactX = theImpactParameter.firs
396 G4double impactY = theImpactParameter.seco
397
398 #ifdef debugQGSParticipants
399 G4cout<<"InteractionMode "<<InteractionM
400 G4cout<<"Impact parameter (fm ) "<<std::
401 G4int nucleonCount = -1;
402 #endif
403
404 // loop over nucleons to find collisions
405 theNucleus->StartLoop();
406 G4QGSParticipants_NPart = 0;
407
408 G4double Power=MaxPower;
409
410 while( (tNucleon = theNucleus->GetNextNucl
411 {
412 if(Power <= 0.) break;
413
414 G4LorentzVector nucleonMomentum=tNucleon
415
416 G4double Distance2 = sqr(impactX - tNucl
417 sqr(impactY - tNucl
418
419 G4double Pint(0.); //
420 G4double Pprd(0.), Ptrd(0.), Pdd(0.); //
421 G4double Pnd (0.), Pnvr(0.); //
422 G4int NcutPomerons(0); //
423
424 Regge->GetProbabilities(std::sqrt(Distan
425 Pint, Pprd, Ptrd, Pdd, Pnd, Pnvr);
426 #ifdef debugQGSParticipants
427 nucleonCount++;
428 G4cout<<"Nucleon & its impact paramete
429 G4cout<<"Probability of interaction:
430 G4cout<<"Probability of PrD, TrD, DD: "<<
431 G4cout<<"Probability of NonDiff, QuarkExc.:
432 #endif
433
434 if (Pint > G4UniformRand())
435 { // An inte
436
437 G4double rndNumber = G4UniformRand();
438 G4int InteractionType(0);
439
440 if((InteractionMode==ALL)||(Interactio
441 {
442 if( rndNumber < Pprd ) {
443 else if( rndNumber < Pprd+Ptrd) {
444 else if( rndNumber < Pprd+Ptrd+Pdd) {
445 else if( rndNumber < Pprd+Ptrd+Pdd+Pnd ) {
446 NcutPomerons = Regge->ncPom
447 else {Intera
448 }
449 else // InteractionMode == NON_DIFF
450 {
451 InteractionMode = NON_DIFF;
452 if( rndNumber < Ptrd ) {Interact
453 else if( rndNumber < Ptrd + Pnd) {Interact
454 }
455
456 if( (InteractionType == NonD) && (Ncut
457
458 G4QGSParticipants_NPart ++;
459 G4QGSMSplitableHadron* aTargetSPB = ne
460 tNucleon->Hit(aTargetSPB);
461
462 #ifdef debugQGSParticipants
463 G4cout<<"An interaction is happend."
464 G4cout<<"Target nucleon - "<<nucleon
465 <<tNucleon->GetDefinition()->G
466 G4cout<<"Interaction type:"<<Interac
467 <<" (0 -PrD, 1 - TrD, 2 - DD,
468 G4cout<<"New Inter. mode:"<<Interac
469 <<" (0 -ALL, 1 - WITHOUT_R, 2
470 if( InteractionType == NonD )
471 G4cout<<"Number of cutted pomerons
472 #endif
473
474 if((InteractionType == PrD) || (Intera
475 (InteractionType == Qexc))
476 { //
477 #ifdef debugQGSParticipants
478 G4cout<<"Diffractive-like interact
479 #endif
480
481 G4InteractionContent * aInteraction
482 theProjectileSplitable->SetStatus(1*
483
484 aInteraction->SetTarget(aTargetSPB);
485 aInteraction->SetTargetNucleon(tNucl
486 aTargetSPB->SetCollisionCount(0);
487 aTargetSPB->SetStatus(1);
488
489 aInteraction->SetNumberOfDiffractive
490 aInteraction->SetNumberOfSoftCollisi
491 aInteraction->SetStatus(InteractionT
492 theInteractions.push_back(aInteracti
493 }
494 else
495 { // non
496 #ifdef debugQGSParticipants
497 G4cout<<"Non-diffractive interacti
498 #endif
499
500 G4int nCuts;
501
502 G4int Vncut=0;
503 for(nCuts = 0; nCuts < NcutPomerons; nCuts
504 {
505 if( G4UniformRand() < Power/MaxPower ){V
506 }
507 nCuts=Vncut;
508
509 if( nCuts == 0 ) {delete aTargetSPB; tNucl
510
511 #ifdef debugQGSParticipants
512 G4cout<<"Number of cuts in the int
513 #endif
514
515 aTargetSPB->IncrementCollisionCount(nCuts)
516 aTargetSPB->SetStatus(0);
517 theTargets.push_back(aTargetSPB);
518
519 theProjectileSplitable->IncrementCollision
520 theProjectileSplitable->SetStatus(0*
521
522 G4InteractionContent * aInteraction = new
523 aInteraction->SetTarget(aTargetSPB);
524 aInteraction->SetTargetNucleon(tNucl
525 aInteraction->SetNumberOfSoftCollisions(nC
526 aInteraction->SetStatus(InteractionT
527 theInteractions.push_back(aInteraction);
528 }
529 } // End of if (Pint > G4UniformRand(
530 } // End of while( (tNucleon = theNucl
531
532 #ifdef debugQGSParticipants
533 G4cout << G4endl<<"Number of wounded nuc
534 #endif
535
536 } // End of while( (theInteractions.size()
537
538 if ( loopCounter >= maxNumberOfLoops ) {
539 #ifdef debugQGSParticipants
540 G4cout <<"BAD situation: forced loop exi
541 #endif
542 // Perhaps there is something to set here.
543 // Decrease impact parameter ??
544 // Select collisions with only diffraction
545 // Selecy only non-diffractive interaction
546 }
547 //------------------------------------------
548 std::vector<G4InteractionContent*>::iterator
549
550 if( theInteractions.size() != 0)
551 {
552 if( InteractionMode == ALL ) // It can be
553 { // Only the
554 i = theInteractions.end()-1;
555
556 while ( theInteractions.size() != 1 )
557 {
558 G4InteractionContent* anInteraction =
559 G4Nucleon * pNucleon = anInteraction->
560 delete anInteraction->GetTarget();
561 delete *i;
562 i=theInteractions.erase(i);
563 i--;
564 }
565 }
566 else
567 { // All quark
568 i = theInteractions.begin();
569 while ( i != theInteractions.end() )
570 {
571 G4InteractionContent* anInteraction =
572
573 if( anInteraction->GetStatus() == Qexc
574 {
575 G4Nucleon* aTargetNucleon = a
576 aTargetNucleon->Hit(nullptr);
577
578 delete anInteraction->GetTarget();
579 delete *i;
580 i=theInteractions.erase(i);
581 }
582 else
583 {
584 i++;
585 }
586 }
587 }
588 }
589 }
590
591 //============================================
592 void G4QGSParticipants::StoreInvolvedNucleon()
593 { //To store nucleons involved in the interact
594
595 NumberOfInvolvedNucleonsOfTarget = 0;
596
597 theNucleus->StartLoop();
598
599 G4Nucleon* aNucleon;
600 while ( ( aNucleon = theNucleus->GetNextNucl
601 if ( aNucleon->AreYouHit() ) {
602 TheInvolvedNucleonsOfTarget[NumberOfInvo
603 NumberOfInvolvedNucleonsOfTarget++;
604 }
605 }
606
607 #ifdef debugQGSParticipants
608 G4cout << G4endl<<"G4QGSParticipants::Stor
609 <<"Stored # of wounded nucleons of
610 << NumberOfInvolvedNucleonsOfTarget
611 #endif
612 return;
613 }
614
615 //============================================
616
617 void G4QGSParticipants::ReggeonCascade()
618 { // Implementation of the reggeon theory insp
619 #ifdef debugQGSParticipants
620 G4cout << G4endl<<"Reggeon cascading .....
621 G4cout<<"C of nucl. desctruction "<<GetCof
622 <<" R2 "<<GetR2ofNuclearDestruction(
623 #endif
624
625 G4int InitNINt = NumberOfInvolvedNucleonsOfT
626
627 // Reggeon cascading in target nucleus
628 for ( G4int InvTN = 0; InvTN < InitNINt; Inv
629 G4Nucleon* aTargetNucleon = TheInvolvedNuc
630
631 G4double CreationTime = aTargetNucleon->Ge
632
633 G4double XofWoundedNucleon = aTargetNucleo
634 G4double YofWoundedNucleon = aTargetNucleo
635
636 G4V3DNucleus* theTargetNucleus = theNucleu
637 theTargetNucleus->StartLoop();
638
639 #ifdef debugQGSParticipants
640 G4int TrgNuc=0;
641 #endif
642 G4Nucleon* Neighbour(0);
643 while ( ( Neighbour = theTargetNucleus->Ge
644 #ifdef debugQGSParticipants
645 TrgNuc++;
646 #endif
647 if ( ! Neighbour->AreYouHit() ) {
648 G4double impact2 = sqr( XofWoundedNucl
649 sqr( YofWoundedNucl
650
651 if ( G4UniformRand() < GetCofNuclearDe
652 G4Exp( -impact2
653 ) {
654 // The neighbour nucleon is involved
655 #ifdef debugQGSParticipants
656 G4cout<<"Target nucleon involved i
657 <<Neighbour->GetDefinition()
658 #endif
659 TheInvolvedNucleonsOfTarget[ NumberO
660 NumberOfInvolvedNucleonsOfTarget++;
661
662 G4QGSMSplitableHadron* targetSplitab
663
664 Neighbour->Hit( targetSplitable );
665 targetSplitable->SetTimeOfCreation(
666 targetSplitable->SetStatus( 2 );
667 targetSplitable->SetCollisionCount(0
668
669 G4InteractionContent * anInteraction
670 anInteraction->SetTarget(targetSplit
671 anInteraction->SetTargetNucleon(Neig
672
673 anInteraction->SetNumberOfDiffractiv
674 anInteraction->SetNumberOfSoftCollis
675 anInteraction->SetStatus(3);
676 theInteractions.push_back(anInteract
677 }
678 }
679 }
680 }
681
682 #ifdef debugQGSParticipants
683 G4cout <<"Number of new involved nucleons
684 #endif
685 return;
686 }
687
688 //============================================
689
690 G4bool G4QGSParticipants::PutOnMassShell() {
691
692 G4bool isProjectileNucleus = false;
693 if ( GetProjectileNucleus() ) {
694 isProjectileNucleus = true;
695 }
696
697 #ifdef debugPutOnMassShell
698 G4cout <<G4endl<< "PutOnMassShell start ..
699 if ( isProjectileNucleus ) {G4cout << "Put
700 #endif
701
702 G4LorentzVector Pprojectile( theProjectile.G
703 if ( Pprojectile.z() < 0.0 ) {
704 return false;
705 }
706
707 G4bool isOk = true;
708
709 G4LorentzVector Ptarget( 0.0, 0.0, 0.0, 0.0
710 G4LorentzVector PtargetResidual( 0.0, 0.0, 0
711 G4double SumMasses = 0.0;
712 G4V3DNucleus* theTargetNucleus = GetTargetNu
713 G4double TargetResidualMass = 0.0;
714
715 #ifdef debugPutOnMassShell
716 G4cout << "Target : ";
717 #endif
718
719 isOk = ComputeNucleusProperties( theTargetNu
720 TargetResid
721 TargetResid
722
723 if ( ! isOk ) return false;
724
725 G4double Mprojectile = 0.0;
726 G4double M2projectile = 0.0;
727 G4LorentzVector Pproj( 0.0, 0.0, 0.0, 0.0 );
728 G4LorentzVector PprojResidual( 0.0, 0.0, 0.0
729 G4V3DNucleus* thePrNucleus = GetProjectileNu
730 G4double PrResidualMass = 0.0;
731
732 if ( ! isProjectileNucleus ) { // hadron-nu
733 Mprojectile = Pprojectile.mag();
734 M2projectile = Pprojectile.mag2();
735 SumMasses += Mprojectile + 20.0*MeV;
736 } else { // nucleus-nucleus or antinucleus-
737
738 #ifdef debugPutOnMassShell
739 G4cout << "Projectile : ";
740 #endif
741
742 isOk = ComputeNucleusProperties( thePrNucl
743 Projectil
744 Projectil
745 if ( ! isOk ) return false;
746 }
747
748 G4LorentzVector Psum = Pprojectile + Ptarget
749 G4double SqrtS = Psum.mag();
750 G4double S = Psum.mag2();
751
752 #ifdef debugPutOnMassShell
753 G4cout << "Pproj "<<Pprojectile<<G4endl;
754 G4cout << "Ptarg "<<Ptarget<<G4endl;
755 G4cout << "Psum " << Psum/GeV << " GeV" <<
756 << "SumMasses, PrResidualMass and T
757 << PrResidualMass/GeV << " " << Tar
758 G4cout << "Ptar res. "<<PtargetResidual<<G
759 #endif
760
761 if ( SqrtS < SumMasses ) {
762 return false; // It is impossible to simu
763 }
764
765 // Try to consider also the excitation energ
766 // possible, with the available energy; othe
767
768 G4double savedSumMasses = SumMasses;
769 if ( isProjectileNucleus ) {
770 SumMasses -= std::sqrt( sqr( PrResidualMas
771 SumMasses += std::sqrt( sqr( PrResidualMas
772 + PprojResidual.pe
773 }
774 SumMasses -= std::sqrt( sqr( TargetResidualM
775 SumMasses += std::sqrt( sqr( TargetResidualM
776 + PtargetResidual.pe
777
778 if ( SqrtS < SumMasses ) {
779 SumMasses = savedSumMasses;
780 if ( isProjectileNucleus ) {
781 ProjectileResidualExcitationEnergy = 0.0
782 }
783 TargetResidualExcitationEnergy = 0.0;
784 }
785
786 TargetResidualMass += TargetResidualExcitati
787 if ( isProjectileNucleus ) {
788 PrResidualMass += ProjectileResidualExcita
789 }
790
791 #ifdef debugPutOnMassShell
792 if ( isProjectileNucleus ) {
793 G4cout << "PrResidualMass ProjResidualEx
794 << ProjectileResidualExcitationEnergy <
795 }
796 G4cout << "TargetResidualMass TargetResidu
797 << TargetResidualExcitationEnergy <
798 << "Sum masses " << SumMasses/GeV <
799 #endif
800
801 // Sampling of nucleons what can transfer to
802 if ( isProjectileNucleus && thePrNucleus->
803 isOk = GenerateDeltaIsobar( SqrtS, NumberO
804 TheInvolvedNuc
805 }
806 if ( theTargetNucleus->GetMassNumber() != 1
807 isOk = isOk &&
808 GenerateDeltaIsobar( SqrtS, NumberO
809 TheInvolvedNuc
810 }
811 if ( ! isOk ) return false;
812
813 // Now we know that it is kinematically poss
814 // of the involved nucleons (or correspondin
815 // We have to sample the kinematical variabl
816 // of the final state. The sampled kinematic
817 // Notice that the sampling of the transvers
818 // Fermi motion.
819
820 // If target is nucleon - return ?
821
822 G4LorentzRotation toCms( -1*Psum.boostVector
823 G4LorentzVector Ptmp = toCms*Pprojectile;
824 if ( Ptmp.pz() <= 0.0 ) { // "String" movin
825 return false;
826 }
827
828 G4LorentzRotation toLab( toCms.inverse() );
829
830 G4double YprojectileNucleus = 0.0;
831 if ( isProjectileNucleus ) {
832 Ptmp = toCms*Pproj;
833 YprojectileNucleus = Ptmp.rapidity();
834 }
835 Ptmp = toCms*Ptarget;
836 G4double YtargetNucleus = Ptmp.rapidity();
837
838 // Ascribing of the involved nucleons Pt and
839 G4double DcorP = 0.0;
840 if ( isProjectileNucleus ) {
841 DcorP = GetDofNuclearDestruction() / thePr
842 }
843 G4double DcorT = GetDofNuclearDestruct
844 G4double AveragePt2 = GetPt2ofNuclearDestru
845 G4double maxPtSquare = GetMaxPt2ofNuclearDes
846
847 #ifdef debugPutOnMassShell
848 if ( isProjectileNucleus ) {
849 G4cout << "Y projectileNucleus " << Ypro
850 }
851 G4cout << "Y targetNucleus " << Ytarge
852 << "Dcor " << GetDofNuclearDestruct
853 << " DcorP DcorT " << DcorP << " "
854 #endif
855
856 G4double M2proj = M2projectile; // Initiali
857 G4double WplusProjectile = 0.0;
858 G4double M2target = 0.0;
859 G4double WminusTarget = 0.0;
860 G4int NumberOfTries = 0;
861 G4double ScaleFactor = 1.0;
862 G4bool OuterSuccess = true;
863
864 const G4int maxNumberOfLoops = 1000;
865 G4int loopCounter = 0;
866 do {
867 G4double sqrtM2proj = 0.0, sqrtM2target =
868 OuterSuccess = true;
869 const G4int maxNumberOfTries = 1000;
870 do {
871 NumberOfTries++;
872 if ( NumberOfTries == 100*(NumberOfTries
873 // After many tries, it is convenient
874 // AveragePt2, so that the sampled mom
875 // involved nucleons (or corresponding delta
876 // it is more likely to satisfy the mo
877 ScaleFactor /= 2.0;
878 DcorP *= ScaleFactor;
879 DcorT *= ScaleFactor;
880 AveragePt2 *= ScaleFactor;
881 }
882 if ( isProjectileNucleus ) {
883 // Sampling of kinematical properties
884 isOk = SamplingNucleonKinematics( Aver
885 theP
886 PrRe
887 Numb
888 TheI
889 }
890 // Sampling of kinematical properties of
891 isOk = isOk &&
892 SamplingNucleonKinematics( Averag
893 theTar
894 Target
895 Number
896 TheInv
897
898 if ( M2proj < 0.0 ) {
899 if( M2proj < -0.000001 ) {
900 G4ExceptionDescription ed;
901 ed << "Projectile " << theProjectile.GetDe
902 << " Target (Z,A)=(" << theTargetNucle
903 << ") M2proj=" << M2proj << " ->
904 G4Exception( "G4QGSParticipants::PutOnMass
905 "HAD_QGSPARTICIPANTS_002", JustWarn
906 }
907 M2proj = 0.0;
908 }
909 sqrtM2proj = std::sqrt( M2proj );
910 if ( M2target < 0.0 ) {
911 G4ExceptionDescription ed;
912 ed << "Projectile " << theProjectile.G
913 << " Target (Z,A)=(" << theTargetN
914 << ") M2target=" << M2target << "
915 G4Exception( "G4QGSParticipants::PutOn
916 "HAD_QGSPARTICIPANTS_003",
917 M2target = 0.0;
918 };
919 sqrtM2target = std::sqrt( M2target );
920
921 #ifdef debugPutOnMassShell
922 G4cout << "SqrtS, Mp+Mt, Mp, Mt " << S
923 << ( sqrtM2proj + sqrtM2target
924 << sqrtM2proj/GeV << " " << sqr
925 #endif
926
927 if ( ! isOk ) return false;
928 } while ( ( SqrtS < ( sqrtM2proj + sqrtM2t
929 ++NumberOfTries < maxNumberOfTri
930 if ( NumberOfTries >= maxNumberOfTries ) {
931 return false;
932 }
933 if ( isProjectileNucleus ) {
934 isOk = CheckKinematics( S, SqrtS, M2proj
935 NumberOfInvolved
936 TheInvolvedNucle
937 WminusTarget, Wp
938 }
939 isOk = isOk &&
940 CheckKinematics( S, SqrtS, M2proj,
941 NumberOfInvolvedNu
942 WminusTarget, Wplu
943 if ( ! isOk ) return false;
944 } while ( ( ! OuterSuccess ) &&
945 ++loopCounter < maxNumberOfLoops )
946 if ( loopCounter >= maxNumberOfLoops ) {
947 return false;
948 }
949
950 // Now the sampling is completed, and we can
951 // whole system. This is done first in the c
952 // to the lab frame. The transverse momentum
953 // the recoil of each hadron (nucleon or del
954 // to conserve (by construction) the transve
955
956 if ( ! isProjectileNucleus ) { // hadron-nu
957
958 G4double Pzprojectile = WplusProjectile/2.
959 G4double Eprojectile = WplusProjectile/2.
960 Pprojectile.setPz( Pzprojectile );
961 Pprojectile.setE( Eprojectile );
962
963 #ifdef debugPutOnMassShell
964 G4cout << "Proj after in CMS " << Pproje
965 #endif
966
967 Pprojectile.transform( toLab );
968 theProjectile.SetMomentum( Pprojectile.vec
969 theProjectile.SetTotalEnergy( Pprojectile.
970
971 if ( theProjectileSplitable ) theProjectil
972
973 #ifdef debugPutOnMassShell
974 G4cout << "Final proj. mom in Lab. " <<t
975 <<t
976 #endif
977
978 } else { // nucleus-nucleus or antinucleus-
979
980 isOk = FinalizeKinematics( WplusProjectile
981 ProjectileResid
982 TheInvolvedNucl
983
984 #ifdef debugPutOnMassShell
985 G4cout << "Projectile Residual4Momentum
986 #endif
987
988 if ( ! isOk ) return false;
989
990 ProjectileResidual4Momentum.transform( toL
991
992 #ifdef debugPutOnMassShell
993 G4cout << "Projectile Residual4Momentum
994 #endif
995
996 }
997
998 isOk = FinalizeKinematics( WminusTarget, fal
999 TargetResidualMas
1000 TheInvolvedNucle
1001
1002 #ifdef debugPutOnMassShell
1003 G4cout << "Target Residual4Momentum in CM
1004 #endif
1005
1006 if ( ! isOk ) return false;
1007
1008 TargetResidual4Momentum.transform( toLab );
1009
1010 #ifdef debugPutOnMassShell
1011 G4cout << "Target Residual4Momentum in La
1012 #endif
1013
1014 return true;
1015
1016 }
1017
1018 //===========================================
1019
1020 G4ThreeVector G4QGSParticipants::GaussianPt(
1021 // @@ this method is used in FTFModel as w
1022
1023 G4double Pt2( 0.0 ), Pt(0.0);
1024 if ( AveragePt2 > 0.0 ) {
1025 G4double x = maxPtSquare/AveragePt2;
1026 Pt2 = (x < 200) ?
1027 -AveragePt2 * G4Log( 1.0 + G4UniformRan
1028 : -AveragePt2 * G4Log( 1.0 - G4UniformR
1029 Pt = std::sqrt( Pt2 );
1030 }
1031 G4double phi = G4UniformRand() * twopi;
1032
1033 return G4ThreeVector( Pt*std::cos(phi), Pt*
1034 }
1035 //===========================================
1036
1037 G4bool G4QGSParticipants::
1038 ComputeNucleusProperties( G4V3DNucleus* nucle
1039 G4LorentzVector& nu
1040 G4LorentzVector& re
1041 G4double& sumMasses
1042 G4double& residualE
1043 G4double& residualM
1044 G4int& residualMass
1045 G4int& residualChar
1046
1047 // This method, which is called only by Put
1048 // - either the target nucleus (which is n
1049 // of hadronic interaction (hadron-nucle
1050 // - or the projectile nucleus or antinucl
1051 // or antinucleus-nucleus interaction.
1052 // This method assumes that the all the par
1053 // the action of this method consists in mo
1054 // first one. The return value is "false" o
1055 // is null.
1056
1057 if ( ! nucleus ) return false;
1058
1059 G4double ExcitationEPerWoundedNucleon = Get
1060
1061 // Loop over the nucleons of the nucleus.
1062 // The nucleons that have been involved in
1063 // Reggeon Cascading) will be candidate to
1064 // All the remaining nucleons will be the n
1065 // The variable sumMasses is the amount of
1066 // 1. transverse mass of each involved
1067 // 2. 20.0*MeV separation energy for ea
1068 // 3. transverse mass of the residual n
1069 // In this first evaluation of sumMasses, t
1070 // (residualExcitationEnergy, estimated by
1071 // nucleon) is not taken into account.
1072 G4Nucleon* aNucleon = 0;
1073 nucleus->StartLoop();
1074 while ( ( aNucleon = nucleus->GetNextNucleo
1075 nucleusMomentum += aNucleon->Get4Momentum
1076 if ( aNucleon->AreYouHit() ) { // Involv
1077 // Consider in sumMasses the nominal, i
1078 // (not the current masses, which could
1079
1080 sumMasses += std::sqrt( sqr( aNucleon->
1081 + aNucleon->Ge
1082 sumMasses += 20.0*MeV; // Separation e
1083
1084 //residualExcitationEnergy += Excitatio
1085 residualExcitationEnergy += -Excitation
1086 residualMassNumber--;
1087 // The absolute value below is needed o
1088 residualCharge -= std::abs( G4int( aNuc
1089 } else { // Spectator nucleons
1090 residualMomentum += aNucleon->Get4Momen
1091 }
1092 }
1093 #ifdef debugPutOnMassShell
1094 G4cout << "ExcitationEnergyPerWoundedNucl
1095 << "\t Residual Charge, MassNumber
1096 << G4endl << "\t Initial Momentum
1097 << G4endl << "\t Residual Momentum
1098 #endif
1099
1100 residualMomentum.setPz( 0.0 );
1101 residualMomentum.setE( 0.0 );
1102 if ( residualMassNumber == 0 ) {
1103 residualMass = 0.0;
1104 residualExcitationEnergy = 0.0;
1105 } else {
1106 residualMass = G4ParticleTable::GetPartic
1107 GetIonMass( residualChar
1108 if ( residualMassNumber == 1 ) {
1109 residualExcitationEnergy = 0.0;
1110 }
1111 }
1112 sumMasses += std::sqrt( sqr( residualMass )
1113 return true;
1114 }
1115
1116
1117 //===========================================
1118
1119 G4bool G4QGSParticipants::
1120 GenerateDeltaIsobar( const G4double sqrtS,
1121 const G4int numberOfInvo
1122 G4Nucleon* involvedNucle
1123 G4double& sumMasses ) {
1124
1125 // This method, which is called only by Put
1126 // re-interpret some of the involved nucleo
1127 // - either by replacing a proton (2212) wi
1128 // - or by replacing a neutron (2112) with
1129 // The on-shell mass of these delta-isobars
1130 // the corresponding nucleon on-shell mass.
1131 // the max number of deltas compatible with
1132 // The delta-isobars are considered with th
1133 // corresponding nucleons.
1134 // This method assumes that all the paramet
1135 // the action of this method consists in mo
1136 // sumMasses. The return value is "false" o
1137 // have unphysical values.
1138
1139 if ( sqrtS < 0.0 || numberOfInvolvedNucle
1140
1141 //const G4double ProbDeltaIsobar = 0.05; /
1142 //const G4double ProbDeltaIsobar = 0.25; /
1143 const G4double probDeltaIsobar = 0.10; //
1144
1145 G4int maxNumberOfDeltas = G4int( (sqrtS - s
1146 G4int numberOfDeltas = 0;
1147
1148 for ( G4int i = 0; i < numberOfInvolvedNucl
1149 //G4cout << "i maxNumberOfDeltas probDelt
1150 // << " " << probDeltaIsobar << G4e
1151 if ( G4UniformRand() < probDeltaIsobar &
1152 numberOfDeltas++;
1153 if ( ! involvedNucleons[i] ) continue;
1154 G4VSplitableHadron* splitableHadron = i
1155 G4double massNuc = std::sqrt( sqr( spli
1156 + splitab
1157 //AR The absolute value below is needed
1158 G4int pdgCode = std::abs( splitableHadr
1159 const G4ParticleDefinition* old_def = s
1160 G4int newPdgCode = pdgCode/10; newPdgCo
1161 if ( splitableHadron->GetDefinition()->
1162 const G4ParticleDefinition* ptr =
1163 G4ParticleTable::GetParticleTable()->
1164 splitableHadron->SetDefinition( ptr );
1165 G4double massDelta = std::sqrt( sqr( sp
1166 + split
1167 //G4cout << i << " " << sqrtS/GeV << "
1168 // << " " << massNuc << G4endl;
1169 if ( sqrtS < sumMasses + massDelta - ma
1170 splitableHadron->SetDefinition( old_d
1171 break;
1172 } else { // Change is accepted
1173 sumMasses += ( massDelta - massNuc );
1174 }
1175 }
1176 }
1177 //G4cout << "maxNumberOfDeltas numberOfDelt
1178 // << numberOfDeltas << G4endl;
1179 return true;
1180 }
1181
1182
1183 //===========================================
1184
1185 G4bool G4QGSParticipants::
1186 SamplingNucleonKinematics( G4double averagePt
1187 const G4double max
1188 G4double dCor,
1189 G4V3DNucleus* nucl
1190 const G4LorentzVec
1191 const G4double res
1192 const G4int residu
1193 const G4int number
1194 G4Nucleon* involve
1195 G4double& mass2 )
1196
1197 // This method, which is called only by Put
1198 // - either the target nucleons: this for
1199 // (hadron-nucleus, nucleus-nucleus, ant
1200 // - or the projectile nucleons or antinuc
1201 // nucleus-nucleus or antinucleus-nucleu
1202 // This method assumes that all the paramet
1203 // the action of this method consists in ch
1204 // whose pointers are in the vector involve
1205 // variable mass2.
1206
1207 if ( ! nucleus ) return false;
1208
1209 if ( residualMassNumber == 0 && numberOfI
1210 dCor = 0.0;
1211 averagePt2 = 0.0;
1212 }
1213
1214 G4bool success = true;
1215
1216 G4double SumMasses = residualMass;
1217 for ( G4int i = 0; i < numberOfInvolvedNucl
1218 G4Nucleon* aNucleon = involvedNucleons[i]
1219 if ( ! aNucleon ) continue;
1220 SumMasses += aNucleon->GetSplitableHadron
1221 }
1222
1223 const G4int maxNumberOfLoops = 1000;
1224 G4int loopCounter = 0;
1225 do {
1226
1227 success = true;
1228 G4ThreeVector ptSum( 0.0, 0.0, 0.0 );
1229 G4double xSum = 0.0;
1230
1231 for ( G4int i = 0; i < numberOfInvolvedNu
1232 G4Nucleon* aNucleon = involvedNucleons[
1233 if ( ! aNucleon ) continue;
1234 G4ThreeVector tmpPt = GaussianPt( avera
1235 ptSum += tmpPt;
1236 G4ThreeVector tmpX = GaussianPt( dCor*d
1237 G4double x = tmpX.x() +
1238 aNucleon->GetSplitableHadr
1239 if ( x < 0.0 || x > 1.0 ) {
1240 success = false;
1241 break;
1242 }
1243 xSum += x;
1244 //AR The energy is in the lab (instead
1245 G4LorentzVector tmp( tmpPt.x(), tmpPt.y
1246 aNucleon->SetMomentum( tmp );
1247 }
1248
1249 if ( xSum < 0.0 || xSum > 1.0 ) success
1250
1251 if ( ! success ) continue;
1252
1253 G4double deltaPx = ( ptSum.x() - pResidua
1254 G4double deltaPy = ( ptSum.y() - pResidua
1255 G4double delta = 0.0;
1256 if ( residualMassNumber == 0 ) {
1257 delta = ( xSum - 1.0 ) / numberOfInvolv
1258 } else {
1259 delta = 0.0;
1260 }
1261
1262 xSum = 1.0;
1263 mass2 = 0.0;
1264 for ( G4int i = 0; i < numberOfInvolvedNu
1265 G4Nucleon* aNucleon = involvedNucleons[
1266 if ( ! aNucleon ) continue;
1267 G4double x = aNucleon->Get4Momentum().p
1268 xSum -= x;
1269 if ( residualMassNumber == 0 ) {
1270 if ( x <= 0.0 || x > 1.0 ) {
1271 success = false;
1272 break;
1273 }
1274 } else {
1275 if ( x <= 0.0 || x > 1.0 || xSum
1276 success = false;
1277 break;
1278 }
1279 }
1280 G4double px = aNucleon->Get4Momentum().
1281 G4double py = aNucleon->Get4Momentum().
1282 mass2 += ( sqr( aNucleon->GetSplitableH
1283 + sqr( px ) + sqr( py ) )
1284 G4LorentzVector tmp( px, py, x, aNucleo
1285 aNucleon->SetMomentum( tmp );
1286 }
1287
1288 if ( success && residualMassNumber != 0
1289 mass2 += ( sqr( residualMass ) + pResid
1290 }
1291
1292 #ifdef debugPutOnMassShell
1293 G4cout << "success " << success << G4en
1294 #endif
1295
1296 } while ( ( ! success ) &&
1297 ++loopCounter < maxNumberOfLoops
1298 if ( loopCounter >= maxNumberOfLoops ) {
1299 success = false;
1300 }
1301
1302 return success;
1303 }
1304
1305
1306 //===========================================
1307
1308 G4bool G4QGSParticipants::
1309 CheckKinematics( const G4double sValue,
1310 const G4double sqrtS,
1311 const G4double projectileMas
1312 const G4double targetMass2,
1313 const G4double nucleusY,
1314 const G4bool isProjectileNuc
1315 const G4int numberOfInvolved
1316 G4Nucleon* involvedNucleons[
1317 G4double& targetWminus,
1318 G4double& projectileWplus,
1319 G4bool& success ) {
1320
1321 // This method, which is called only by Put
1322 // kinematics is acceptable or not.
1323 // This method assumes that all the paramet
1324 // notice that the input boolean parameter
1325 // only in the case of nucleus or antinucle
1326 // The action of this method consists in co
1327 // and setting the parameter success to fal
1328 // be rejeted.
1329
1330 G4double decayMomentum2 = sqr( sValue ) + s
1331 - 2.0*sValue*proj
1332 - 2.0*projectileM
1333 targetWminus = ( sValue - projectileMass2 +
1334 / 2.0 / sqrtS;
1335 projectileWplus = sqrtS - targetMass2/targe
1336 G4double projectilePz = projectileWplus/2.0
1337 G4double projectileE = projectileWplus/2.0
1338 G4double projectileY(1.0e5);
1339 if (projectileE - projectilePz > 0.) {
1340 projectileY = 0.5 * G4Log( (proje
1341 (proje
1342 }
1343 G4double targetPz = -targetWminus/2.0 + tar
1344 G4double targetE = targetWminus/2.0 + tar
1345 G4double targetY = 0.5 * G4Log( (targetE +
1346
1347 #ifdef debugPutOnMassShell
1348 G4cout << "decayMomentum2 " << decayMomen
1349 << "\t targetWminus projectileWplu
1350 << "\t projectileY targetY " << pr
1351 #endif
1352
1353 for ( G4int i = 0; i < numberOfInvolvedNucl
1354 G4Nucleon* aNucleon = involvedNucleons[i]
1355 if ( ! aNucleon ) continue;
1356 G4LorentzVector tmp = aNucleon->Get4Momen
1357 G4double mt2 = sqr( tmp.x() ) + sqr( tmp.
1358 sqr( aNucleon->GetSplitabl
1359 G4double x = tmp.z();
1360 G4double pz = -targetWminus*x/2.0 + mt2/(
1361 G4double e = targetWminus*x/2.0 + mt2/(
1362 if ( isProjectileNucleus ) {
1363 pz = projectileWplus*x/2.0 - mt2/(2.0*p
1364 e = projectileWplus*x/2.0 + mt2/(2.0*p
1365 }
1366 G4double nucleonY = 0.5 * G4Log( (e + pz)
1367
1368 #ifdef debugPutOnMassShell
1369 G4cout << "i nY pY nY-AY AY " << i << "
1370 #endif
1371
1372 if ( std::abs( nucleonY - nucleusY ) > 2
1373 ( isProjectileNucleus && targetY >
1374 ( ! isProjectileNucleus && project
1375 success = false;
1376 break;
1377 }
1378 }
1379 return true;
1380 }
1381
1382
1383 //===========================================
1384
1385 G4bool G4QGSParticipants::
1386 FinalizeKinematics( const G4double w,
1387 const G4bool isProjectile
1388 const G4LorentzRotation&
1389 const G4double residualMa
1390 const G4int residualMassN
1391 const G4int numberOfInvol
1392 G4Nucleon* involvedNucleo
1393 G4LorentzVector& residual4Momen
1394
1395 // This method, which is called only by Put
1396 // this method is called when we are sure t
1397 // acceptable.
1398 // This method assumes that all the paramet
1399 // notice that the input boolean parameter
1400 // only in the case of nucleus or antinucle
1401 // because the sign of pz (in the center-of
1402 // with respect to the case of a normal had
1403 // The action of this method consists in mo
1404 // (in the lab frame) and computing the res
1405 // frame).
1406
1407 G4ThreeVector residual3Momentum( 0.0, 0.0,
1408
1409 for ( G4int i = 0; i < numberOfInvolvedNucl
1410 G4Nucleon* aNucleon = involvedNucleons[i]
1411 if ( ! aNucleon ) continue;
1412 G4LorentzVector tmp = aNucleon->Get4Momen
1413 residual3Momentum -= tmp.vect();
1414 G4double mt2 = sqr( tmp.x() ) + sqr( tmp.
1415 sqr( aNucleon->GetSplitabl
1416 G4double x = tmp.z();
1417 G4double pz = -w * x / 2.0 + mt2 / ( 2.
1418 G4double e = w * x / 2.0 + mt2 / ( 2.
1419 // Reverse the sign of pz in the case of
1420 if ( isProjectileNucleus ) pz *= -1.0;
1421 tmp.setPz( pz );
1422 tmp.setE( e );
1423 tmp.transform( boostFromCmsToLab );
1424 aNucleon->SetMomentum( tmp );
1425 G4VSplitableHadron* splitableHadron = aNu
1426 splitableHadron->Set4Momentum( tmp );
1427 #ifdef debugPutOnMassShell
1428 G4cout << "Target involved nucleon No,
1429 << i<<" "<<aNucleon->GetDefinitio
1430 #endif
1431 }
1432
1433 G4double residualMt2 = sqr( residualMass )
1434 + sqr( residual3Moment
1435
1436 #ifdef debugPutOnMassShell
1437 G4cout <<G4endl<< "w residual3Momentum.z(
1438 #endif
1439
1440 G4double residualPz = 0.0;
1441 G4double residualE = 0.0;
1442 if ( residualMassNumber != 0 ) {
1443 residualPz = -w * residual3Momentum.z() /
1444 residualMt2 / ( 2.0 * w * r
1445 residualE = w * residual3Momentum.z() /
1446 residualMt2 / ( 2.0 * w * r
1447 // Reverse the sign of residualPz in the
1448 if ( isProjectileNucleus ) residualPz *=
1449 }
1450
1451 residual4Momentum.setPx( residual3Momentum.
1452 residual4Momentum.setPy( residual3Momentum.
1453 residual4Momentum.setPz( residualPz );
1454 residual4Momentum.setE( residualE );
1455
1456 return true;
1457 }
1458
1459 //===========================================
1460 void G4QGSParticipants::PerformDiffractiveCol
1461 {
1462 #ifdef debugQGSParticipants
1463 G4cout<<G4endl<<"PerformDiffractiveCollisi
1464 <<"theInteractions.size() "<<theInte
1465 #endif
1466
1467 unsigned int i;
1468 for (i = 0; i < theInteractions.size(); i++
1469 {
1470 G4InteractionContent* anIniteraction = th
1471 #ifdef debugQGSParticipants
1472 G4cout<<"Interaction # and its status "
1473 <<i<<" "<<theInteractions[i]->Get
1474 #endif
1475
1476 G4int InterStatus = theInteractions[i]->G
1477 if ( (InterStatus == PrD) || (InterStatus
1478 { // Selection of diffractive interactio
1479 #ifdef debugQGSParticipants
1480 G4cout<<"Simulation of diffractive in
1481 #endif
1482
1483 G4VSplitableHadron* aTarget = anInitera
1484
1485 #ifdef debugQGSParticipants
1486 G4cout<<"The proj. before inter "
1487 <<theProjectileSplitable->Get4M
1488 <<theProjectileSplitable->Get4M
1489 G4cout<<"The targ. before inter " <<a
1490 <<aTarget->Get4Momentum().mag()
1491 #endif
1492
1493 if ( InterStatus == PrD )
1494 theSingleDiffExcitation.ExcitePartici
1495
1496 if ( InterStatus == TrD )
1497 theSingleDiffExcitation.ExcitePartici
1498
1499 if ( InterStatus == DD )
1500 theDiffExcitaton.ExciteParticipants(t
1501
1502 #ifdef debugQGSParticipants
1503 G4cout<<"The proj. after inter " <<t
1504 <<theProjectileSplitable->Get4M
1505 G4cout<<"The targ. after inter " <<aTarget
1506 <<aTarget->Get4Momentum().mag()
1507 #endif
1508 }
1509
1510 if ( InterStatus == Qexc )
1511 { // Quark exchange process
1512 #ifdef debugQGSParticipants
1513 G4cout<<"Simulation of interaction wi
1514 #endif
1515 G4VSplitableHadron* aTarget = anInitera
1516
1517 #ifdef debugQGSParticipants
1518 G4cout<<"The proj. before inter " <<t
1519 <<theProjectileSplitable->Get4M
1520 G4cout<<"The targ. before inter "<<aT
1521 <<aTarget->Get4Momentum().mag()
1522 #endif
1523
1524 theQuarkExchange.ExciteParticipants(the
1525
1526 #ifdef debugQGSParticipants
1527 G4cout<<"The proj. after inter " <<t
1528 <<theProjectileSplitable->Get4M
1529 G4cout<<"The targ. after inter " <<aTarget
1530 <<aTarget->Get4Momentum().mag()
1531 #endif
1532 }
1533 }
1534 }
1535
1536 //===========================================
1537 G4bool G4QGSParticipants::DeterminePartonMome
1538 {
1539 if ( ! theProjectileSplitable ) return fals
1540
1541 const G4double aHugeValue = 1.0e+10;
1542
1543 #ifdef debugQGSParticipants
1544 G4cout<<G4endl<<"DeterminePartonMomenta()
1545 G4cout<<"theProjectile status (0 -nondiff
1546 #endif
1547
1548 if (theProjectileSplitable->GetStatus() !=
1549
1550 G4LorentzVector Projectile4Momentum = theP
1551 G4LorentzVector Psum = Projectile4Momentum;
1552
1553 G4double VqM_pr(0.), VaqM_pr(0.), VqM_tr(35
1554 if (std::abs(theProjectile.GetDefinition()-
1555
1556 #ifdef debugQGSParticipants
1557 G4cout<<"Projectile 4 momentum "<<Psum<<G
1558 <<"Target nucleon momenta at start"
1559 G4int NuclNo=0;
1560 #endif
1561
1562 std::vector<G4VSplitableHadron*>::iterator
1563
1564 for (i = theTargets.begin(); i != theTarget
1565 {
1566 Psum += (*i)->Get4Momentum();
1567 #ifdef debugQGSParticipants
1568 G4cout<<"Nusleus nucleon # and its 4Mom
1569 NuclNo++;
1570 #endif
1571 }
1572
1573 G4LorentzRotation toCms( -1*Psum.boostVecto
1574
1575 G4LorentzVector Ptmp = toCms*Projectile4Mom
1576
1577 toCms.rotateZ( -1*Ptmp.phi() );
1578 toCms.rotateY( -1*Ptmp.theta() );
1579 G4LorentzRotation toLab(toCms.inverse());
1580 Projectile4Momentum.transform( toCms );
1581 // Ptarget.transform( toCms );
1582
1583 #ifdef debugQGSParticipants
1584 G4cout<<G4endl<<"In CMS---------------"<<
1585 G4cout<<"Projectile 4 Mom "<<Projectile4M
1586 NuclNo=0;
1587 #endif
1588
1589 G4LorentzVector Target4Momentum(0.,0.,0.,0.
1590 for(i = theTargets.begin(); i != theTargets
1591 {
1592 G4LorentzVector tmp= (*i)->Get4Momentum()
1593 (*i)->Set4Momentum( tmp );
1594 #ifdef debugQGSParticipants
1595 G4cout<<"Target nucleon # and 4Mom "<<"
1596 NuclNo++;
1597 #endif
1598 Target4Momentum += tmp;
1599 }
1600
1601 G4double S = Psum.mag2();
1602 G4double SqrtS = std::sqrt(S);
1603
1604 #ifdef debugQGSParticipants
1605 G4cout<<"Sum of target nucleons 4 momentu
1606 G4cout<<"Target nucleons mom: px, py, z_1
1607 NuclNo=0;
1608 #endif
1609
1610 //G4double PplusProjectile = Projectile4Mom
1611 G4double PminusTarget = Target4Momentum.
1612
1613 for(i = theTargets.begin(); i != theTargets
1614 {
1615 G4LorentzVector tmp = (*i)->Get4Momentum(
1616
1617 //AR-19Jan2017 : the following line is ca
1618 // is built in optimized mo
1619 // To fix it, I get the mas
1620 // mass from the Lorentz ve
1621 // square root. If the mass
1622 // exception is thrown, and
1623 //G4double Mass = tmp.mag();
1624 G4double Mass2 = tmp.mag2();
1625 G4double Mass = 0.0;
1626 if ( Mass2 < 0.0 ) {
1627 G4ExceptionDescription ed;
1628 ed << "Projectile " << theProjectile.Ge
1629 << "Â 4-momentum " << Psum << G4end
1630 ed << "LorentzVector tmp " << tmp << "
1631 G4Exception( "G4QGSParticipants::Determ
1632 "HAD_QGSPARTICIPANTS_001",
1633 } else {
1634 Mass = std::sqrt( Mass2 );
1635 }
1636
1637 tmp.setPz(tmp.minus()/PminusTarget); tm
1638 (*i)->Set4Momentum(tmp);
1639 #ifdef debugQGSParticipants
1640 G4cout<<"Target nucleons # and mom: "<<
1641 NuclNo++;
1642 #endif
1643 }
1644
1645 //+++++++++++++++++++++++++++++++++++++++++
1646
1647 G4double SigPt = sigmaPt;
1648 G4Parton* aParton(0);
1649 G4ThreeVector aPtVector(0.,0.,0.);
1650 G4LorentzVector tmp(0.,0.,0.,0.);
1651
1652 G4double Mt(0.);
1653 G4double ProjSumMt(0.), ProjSumMt2perX(0.);
1654 G4double TargSumMt(0.), TargSumMt2perX(0.);
1655
1656
1657 G4double aBeta = beta; // Member of the c
1658
1659 const G4ParticleDefinition* theProjectileDe
1660 if (theProjectileDefinition == G4PionMinus:
1661 if (theProjectileDefinition == G4Gamma::Gam
1662 if (theProjectileDefinition == G4PionPlus::
1663 if (theProjectileDefinition == G4PionZero::
1664 if (theProjectileDefinition == G4KaonPlus::
1665 if (theProjectileDefinition == G4KaonMinus:
1666
1667 G4double Xmin = 0.;
1668
1669 G4bool Success = true; G4int attempt = 0;
1670 const G4int maxNumberOfAttempts = 1000;
1671 do
1672 {
1673 attempt++; if( attempt == 100*(attempt/
1674
1675 ProjSumMt=0.; ProjSumMt2perX=0.;
1676 TargSumMt=0.; TargSumMt2perX=0.;
1677
1678 Success = true;
1679 G4int nSeaPair = theProjectileSplitable->
1680 #ifdef debugQGSParticipants
1681 G4cout<<"attempt ----------------------
1682 G4cout<<"nSeaPair of proj "<<nSeaPair<<
1683 #endif
1684
1685 G4double SumPx = 0.;
1686 G4double SumPy = 0.;
1687 G4double SumZ = 0.;
1688 G4int NumberOfUnsampledSeaQuarks = 2*nSea
1689
1690 G4double Qmass=0.;
1691 for (G4int aSeaPair = 0; aSeaPair < nSeaP
1692 {
1693 aParton = theProjectileSplitable->GetNe
1694 #ifdef debugQGSParticipants
1695 G4cout<<"Sea quarks: "<<aSeaPair<<" "
1696 #endif
1697 aPtVector = GaussianPt(SigPt, aHugeValu
1698 tmp.setPx(aPtVector.x()); tmp.setPy(aPt
1699 SumPx += aPtVector.x(); SumPy += aPtV
1700 Mt = std::sqrt(aPtVector.mag2()+sqr(Qma
1701 ProjSumMt += Mt;
1702
1703 // Sampling of Z fraction
1704 tmp.setPz(SampleX(Xmin, NumberOfUnsampl
1705 SumZ += tmp.z();
1706
1707 NumberOfUnsampledSeaQuarks--;
1708 ProjSumMt2perX +=sqr(Mt)/tmp.pz();
1709 tmp.setE(sqr(Mt));
1710 aParton->Set4Momentum(tmp);
1711
1712 aParton = theProjectileSplitable->GetNe
1713 #ifdef debugQGSParticipants
1714 G4cout<<" "<<aParton->GetDefinition()
1715 G4cout<<" "<<tmp<<" "<<S
1716 #endif
1717 aPtVector = GaussianPt(SigPt, aHugeValu
1718 tmp.setPx(aPtVector.x()); tmp.setPy(aPt
1719 SumPx += aPtVector.x(); SumPy += aPtV
1720 Mt = std::sqrt(aPtVector.mag2()+sqr(Qma
1721 ProjSumMt += Mt;
1722
1723 // Sampling of Z fraction
1724 tmp.setPz(SampleX(Xmin, NumberOfUnsampl
1725 SumZ += tmp.z();
1726
1727 NumberOfUnsampledSeaQuarks--;
1728 ProjSumMt2perX +=sqr(Mt)/tmp.pz();
1729 tmp.setE(sqr(Mt));
1730 aParton->Set4Momentum(tmp);
1731 #ifdef debugQGSParticipants
1732 G4cout<<" "<<tmp<<" "<<S
1733 #endif
1734 }
1735
1736 // For valence quark
1737 aParton = theProjectileSplitable->GetNext
1738 #ifdef debugQGSParticipants
1739 G4cout<<"Val quark of Pr"<<" "<<aParton
1740 #endif
1741 aPtVector = GaussianPt(SigPt, aHugeValue)
1742 tmp.setPx(aPtVector.x()); tmp.setPy(aPtVe
1743 SumPx += aPtVector.x(); SumPy += aPtVec
1744 Mt = std::sqrt(aPtVector.mag2()+sqr(VqM_p
1745 ProjSumMt += Mt;
1746
1747 // Sampling of Z fraction
1748 tmp.setPz(SampleX(Xmin, NumberOfUnsampled
1749 SumZ += tmp.z();
1750
1751 ProjSumMt2perX +=sqr(Mt)/tmp.pz();
1752 tmp.setE(sqr(Mt));
1753 aParton->Set4Momentum(tmp);
1754
1755 // For valence di-quark
1756 aParton = theProjectileSplitable->GetNext
1757 #ifdef debugQGSParticipants
1758 G4cout<<" "<<aParton->GetDefinition()->
1759 G4cout<<" "<<tmp<<" "<<Sum
1760 #endif
1761 tmp.setPx(-SumPx); tmp.setPy(-SumPy);
1762 //Uzhi 2019 Mt = std::sqrt(aPtVector.mag
1763 Mt = std::sqrt( sqr(SumPx) + sqr(SumPy) +
1764 ProjSumMt += Mt;
1765 tmp.setPz(1.-SumZ);
1766
1767 ProjSumMt2perX +=sqr(Mt)/tmp.pz(); // QQ
1768 tmp.setE(sqr(Mt));
1769 aParton->Set4Momentum(tmp);
1770 #ifdef debugQGSParticipants
1771 G4cout<<" "<<tmp<<" "<<Sum
1772 NuclNo=0;
1773 #endif
1774
1775 // End of work with the projectile
1776
1777 // Work with target nucleons
1778
1779 for(i = theTargets.begin(); i != theTarge
1780 {
1781 nSeaPair = (*i)->GetSoftCollisionCount(
1782 #ifdef debugQGSParticipants
1783 G4cout<<"nSeaPair of target N "<<nSea
1784 <<"Target nucleon 4Mom "<<(*i)-
1785 #endif
1786
1787 SumPx = (*i)->Get4Momentum().px() * (-1
1788 SumPy = (*i)->Get4Momentum().py() * (-1
1789 SumZ = 0.;
1790
1791 NumberOfUnsampledSeaQuarks = 2*nSeaPair
1792
1793 Qmass=0;
1794 for (G4int aSeaPair = 0; aSeaPair < nSe
1795 {
1796 aParton = (*i)->GetNextParton(); //
1797 #ifdef debugQGSParticipants
1798 G4cout<<"Sea quarks: "<<aSeaPair<<"
1799 #endif
1800 aPtVector = GaussianPt(SigPt, aHugeVa
1801 tmp.setPx(aPtVector.x()); tmp.setPy(a
1802 SumPx += aPtVector.x(); SumPy += aP
1803 Mt=std::sqrt(aPtVector.mag2()+sqr(Qma
1804 TargSumMt += Mt;
1805
1806 // Sampling of Z fraction
1807 tmp.setPz(SampleX(Xmin, NumberOfUnsam
1808 SumZ += tmp.z();
1809 tmp.setPz((*i)->Get4Momentum().pz()*t
1810 NumberOfUnsampledSeaQuarks--;
1811 TargSumMt2perX +=sqr(Mt)/tmp.pz();
1812 tmp.setE(sqr(Mt));
1813 aParton->Set4Momentum(tmp);
1814
1815 aParton = (*i)->GetNextAntiParton();
1816 #ifdef debugQGSParticipants
1817 G4cout<<" "<<aParton->GetDefinition
1818 G4cout<<" "<<tmp<<" "<
1819 #endif
1820 aPtVector = GaussianPt(SigPt, aHugeVa
1821 tmp.setPx(aPtVector.x()); tmp.setPy(a
1822 SumPx += aPtVector.x(); SumPy += aP
1823 Mt=std::sqrt(aPtVector.mag2()+sqr(Qma
1824 TargSumMt += Mt;
1825
1826 // Sampling of Z fraction
1827 tmp.setPz(SampleX(Xmin, NumberOfUnsam
1828 SumZ += tmp.z();
1829 tmp.setPz((*i)->Get4Momentum().pz()*t
1830 NumberOfUnsampledSeaQuarks--;
1831 TargSumMt2perX +=sqr(Mt)/tmp.pz();
1832 tmp.setE(sqr(Mt));
1833 aParton->Set4Momentum(tmp);
1834 #ifdef debugQGSParticipants
1835 G4cout<<" "<<tmp<<" "<
1836 #endif
1837 }
1838
1839 // Valence quark
1840 aParton = (*i)->GetNextParton(); // f
1841 #ifdef debugQGSParticipants
1842 G4cout<<"Val quark of Tr"<<" "<<aPart
1843 #endif
1844 aPtVector = GaussianPt(SigPt, aHugeValu
1845 tmp.setPx(aPtVector.x()); tmp.setPy(aPt
1846 SumPx += aPtVector.x(); SumPy += aPtV
1847 Mt=std::sqrt(aPtVector.mag2()+sqr(VqM_t
1848 TargSumMt += Mt;
1849
1850 // Sampling of Z fraction
1851 tmp.setPz(SampleX(Xmin, NumberOfUnsampl
1852 SumZ += tmp.z();
1853 tmp.setPz((*i)->Get4Momentum().pz()*tmp
1854 TargSumMt2perX +=sqr(Mt)/tmp.pz();
1855 tmp.setE(sqr(Mt));
1856 aParton->Set4Momentum(tmp);
1857
1858 // Valence di-quark
1859 aParton = (*i)->GetNextAntiParton();
1860 #ifdef debugQGSParticipants
1861 G4cout<<" "<<aParton->GetDefinition()
1862 G4cout<<" "<<tmp<<" "<<S
1863 #endif
1864 tmp.setPx(-SumPx); tmp
1865 //Uzhi 2019 Mt=std::sqrt(aPtVector.mag
1866 Mt=std::sqrt( sqr(SumPx) + sqr(SumPy) +
1867 TargSumMt += Mt;
1868
1869 tmp.setPz((*i)->Get4Momentum().pz()*(1.
1870 TargSumMt2perX +=sqr(Mt)/tmp.pz();
1871 tmp.setE(sqr(Mt));
1872 aParton->Set4Momentum(tmp);
1873 #ifdef debugQGSParticipants
1874 G4cout<<" "<<tmp<<" "<<1
1875 #endif
1876
1877 } // End of for(i = theTargets.begin();
1878
1879 if( ProjSumMt + TargSumMt > Sqr
1880 Success = false; continue;}
1881 if( std::sqrt(ProjSumMt2perX) + std::sqrt
1882 Success = false; continue;}
1883
1884 } while( (!Success) &&
1885 attempt < maxNumberOfAttempts );
1886
1887 if ( attempt >= maxNumberOfAttempts ) {
1888 return false;
1889 }
1890
1891 //+++++++++++++++++++++++++++++++++++++++++
1892
1893 G4double DecayMomentum2 = sqr(S) + sqr(Proj
1894 - 2.0*S*ProjSumMt2perX - 2.0*S
1895
1896 G4double targetWminus=( S - ProjSumMt2perX
1897 G4double projectileWplus = SqrtS - TargSumM
1898
1899 G4LorentzVector Tmp(0.,0.,0.,0.);
1900 G4double z(0.);
1901
1902 G4int nSeaPair = theProjectileSplitable->Ge
1903 #ifdef debugQGSParticipants
1904 G4cout<<"Backward transformation ========
1905 G4cout<<"nSeaPair of proj "<<nSeaPair<<G4
1906 #endif
1907
1908 for (G4int aSeaPair = 0; aSeaPair < nSeaPai
1909 {
1910 aParton = theProjectileSplitable->GetNext
1911 #ifdef debugQGSParticipants
1912 G4cout<<"Sea quarks: "<<aSeaPair<<" "<<
1913 #endif
1914 Tmp =aParton->Get4Momentum(); z=Tmp.z();
1915
1916 Tmp.setPz(projectileWplus*z/2.0 - Tmp.e()
1917 Tmp.setE( projectileWplus*z/2.0 + Tmp.e()
1918 Tmp.transform( toLab );
1919
1920 aParton->Set4Momentum(Tmp);
1921
1922 aParton = theProjectileSplitable->GetNext
1923 #ifdef debugQGSParticipants
1924 G4cout<<" "<<aParton->GetDefinition()->
1925 G4cout<<" "<<Tmp<<" "<<Tmp
1926 #endif
1927 Tmp =aParton->Get4Momentum(); z=Tmp.z();
1928 Tmp.setPz(projectileWplus*z/2.0 - Tmp.e()
1929 Tmp.setE( projectileWplus*z/2.0 + Tmp.e()
1930 Tmp.transform( toLab );
1931
1932 aParton->Set4Momentum(Tmp);
1933 #ifdef debugQGSParticipants
1934 G4cout<<" "<<Tmp<<" "<<Tmp
1935 #endif
1936 }
1937
1938 // For valence quark
1939 aParton = theProjectileSplitable->GetNextPa
1940 #ifdef debugQGSParticipants
1941 G4cout<<"Val quark of Pr"<<" "<<aParton->
1942 #endif
1943 Tmp =aParton->Get4Momentum(); z=Tmp.z();
1944 Tmp.setPz(projectileWplus*z/2.0 - Tmp.e()/(
1945 Tmp.setE( projectileWplus*z/2.0 + Tmp.e()/(
1946 Tmp.transform( toLab );
1947
1948 aParton->Set4Momentum(Tmp);
1949
1950 // For valence di-quark
1951 aParton = theProjectileSplitable->GetNextAn
1952 #ifdef debugQGSParticipants
1953 G4cout<<" "<<aParton->GetDefinition()->Ge
1954 G4cout<<" "<<Tmp<<" "<<Tmp.m
1955 #endif
1956 Tmp =aParton->Get4Momentum(); z=Tmp.z();
1957 Tmp.setPz(projectileWplus*z/2.0 - Tmp.e()/(
1958 Tmp.setE( projectileWplus*z/2.0 + Tmp.e()/(
1959 Tmp.transform( toLab );
1960
1961 aParton->Set4Momentum(Tmp);
1962
1963 #ifdef debugQGSParticipants
1964 G4cout<<" "<<Tmp<<" "<<Tmp.m
1965 NuclNo=0;
1966 #endif
1967
1968 // End of work with the projectile
1969
1970 // Work with target nucleons
1971 for(i = theTargets.begin(); i != theTargets
1972 {
1973 nSeaPair = (*i)->GetSoftCollisionCount()-
1974 #ifdef debugQGSParticipants
1975 G4cout<<"nSeaPair of target and N# "<<n
1976 NuclNo++;
1977 #endif
1978 for (G4int aSeaPair = 0; aSeaPair < nSeaP
1979 {
1980 aParton = (*i)->GetNextParton(); // f
1981 #ifdef debugQGSParticipants
1982 G4cout<<"Sea quarks: "<<aSeaPair<<" "
1983 #endif
1984 Tmp =aParton->Get4Momentum(); z=Tmp.z()
1985 Tmp.setPz(-targetWminus*z/2.0 + Tmp.e()
1986 Tmp.setE( targetWminus*z/2.0 + Tmp.e()
1987 Tmp.transform( toLab );
1988
1989 aParton->Set4Momentum(Tmp);
1990
1991 aParton = (*i)->GetNextAntiParton();
1992 #ifdef debugQGSParticipants
1993 G4cout<<" "<<aParton->GetDefinition()
1994 G4cout<<" "<<Tmp<<" "<<T
1995 #endif
1996 Tmp =aParton->Get4Momentum(); z=Tmp.z()
1997 Tmp.setPz(-targetWminus*z/2.0 + Tmp.e()
1998 Tmp.setE( targetWminus*z/2.0 + Tmp.e()
1999 Tmp.transform( toLab );
2000
2001 aParton->Set4Momentum(Tmp);
2002 #ifdef debugQGSParticipants
2003 G4cout<<" "<<Tmp<<" "<<T
2004 #endif
2005 }
2006
2007 // Valence quark
2008
2009 aParton = (*i)->GetNextParton(); // for
2010 #ifdef debugQGSParticipants
2011 G4cout<<"Val quark of Tr"<<" "<<aParton
2012 #endif
2013 Tmp =aParton->Get4Momentum(); z=Tmp.z();
2014 Tmp.setPz(-targetWminus*z/2.0 + Tmp.e()/(
2015 Tmp.setE( targetWminus*z/2.0 + Tmp.e()/(
2016 Tmp.transform( toLab );
2017
2018 aParton->Set4Momentum(Tmp);
2019
2020 // Valence di-quark
2021 aParton = (*i)->GetNextAntiParton(); //
2022 #ifdef debugQGSParticipants
2023 G4cout<<" "<<aParton->GetDefinition()->
2024 G4cout<<" "<<Tmp<<" "<<Tmp
2025 #endif
2026 Tmp =aParton->Get4Momentum(); z=Tmp.z();
2027 Tmp.setPz(-targetWminus*z/2.0 + Tmp.e()/(
2028 Tmp.setE( targetWminus*z/2.0 + Tmp.e()/(
2029 Tmp.transform( toLab );
2030
2031 aParton->Set4Momentum(Tmp);
2032 #ifdef debugQGSParticipants
2033 G4cout<<" "<<Tmp<<" "<<Tmp
2034 NuclNo++;
2035 #endif
2036 } // End of for(i = theTargets.begin(); i
2037
2038 return true;
2039 }
2040
2041 //===========================================
2042 G4double G4QGSParticipants::
2043 SampleX(G4double, G4int nSea, G4int, G4double
2044 {
2045 G4double Oalfa = 1./(alpha + 1.);
2046 G4double Obeta = 1./(aBeta + (alpha + 1.)*n
2047
2048 G4double Ksi1, Ksi2, r1, r2, r12;
2049 const G4int maxNumberOfLoops = 1000;
2050 G4int loopCounter = 0;
2051 do
2052 {
2053 Ksi1 = G4UniformRand(); r1 = G4Pow::GetIn
2054 Ksi2 = G4UniformRand(); r2 = G4Pow::GetIn
2055 r12=r1+r2;
2056 } while( ( r12 > 1.) &&
2057 ++loopCounter < maxNumberOfLoops )
2058 if ( loopCounter >= maxNumberOfLoops ) {
2059 return 0.5; // Just an acceptable value,
2060 }
2061
2062 G4double result = r1/r12;
2063 return result;
2064 }
2065
2066 //===========================================
2067 void G4QGSParticipants::CreateStrings()
2068 {
2069
2070 #ifdef debugQGSParticipants
2071 G4cout<<"CreateStrings() ................
2072 #endif
2073
2074 if ( ! theProjectileSplitable ) {
2075 #ifdef debugQGSParticipants
2076 G4cout<<"BAD situation: theProjectileSp
2077 #endif
2078 return;
2079 }
2080
2081 #ifdef debugQGSParticipants
2082 G4cout<<"theProjectileSplitable->GetStatu
2083 G4LorentzVector str4Mom;
2084 #endif
2085
2086 if( theProjectileSplitable->GetStatus() ==
2087 {
2088 G4ThreeVector Position = theProjectileSpl
2089
2090 G4PartonPair * aPair = new G4PartonPair(t
2091 t
2092 G4PartonPair::DIFFRACTIVE
2093 #ifdef debugQGSParticipants
2094 G4cout << "Pr. Diffr. String: Qs 4mom X
2095 G4cout << " " << aPair->Ge
2096 << aPair->GetParton1()->Get4Momentum
2097 << aPair->GetParton1()->GetX()
2098 G4cout << " " << aPair->Ge
2099 << aPair->GetParton2()->Get4Momentum
2100 << aPair->GetParton2()->GetX()
2101 str4Mom += aPair->GetParton1()->Get4Mom
2102 str4Mom += aPair->GetParton2()->Get4Mom
2103 #endif
2104
2105 thePartonPairs.push_back(aPair);
2106 }
2107
2108 G4int N_EnvTarg = NumberOfInvolvedNucleonsO
2109
2110 for ( G4int i = 0; i < N_EnvTarg; i++ ) {
2111 G4Nucleon* aTargetNucleon = TheInvolvedNu
2112
2113 G4VSplitableHadron* HitTargetNucleon = aT
2114
2115 #ifdef debugQGSParticipants
2116 G4cout<<"Involved Nucleon # and its sta
2117 #endif
2118 if( HitTargetNucleon->GetStatus() >= 1)
2119 {
2120 G4ThreeVector Position = HitTargetN
2121
2122 G4PartonPair * aPair = new G4PartonPair
2123
2124 G4PartonPair::DIFFRACTI
2125 #ifdef debugQGSParticipants
2126 G4cout << "Tr. Diffr. String: Qs 4mom
2127 G4cout << "Diffr. String " << aPair->GetPar
2128 << aPair->GetParton1()->Get4Moment
2129 << aPair->GetParton1()->GetX()
2130 G4cout << " " << aPair->GetPar
2131 << aPair->GetParton2()->Get4Moment
2132 << aPair->GetParton2()->GetX()
2133
2134 str4Mom += aPair->GetParton1()->Get4Momentu
2135 str4Mom += aPair->GetParton2()->Get4Momentu
2136 #endif
2137
2138 thePartonPairs.push_back(aPair);
2139 }
2140 }
2141
2142 //-----------------------------------------
2143 #ifdef debugQGSParticipants
2144 G4int IntNo=0;
2145 G4cout<<"Strings created in soft interact
2146 #endif
2147 std::vector<G4InteractionContent*>::iterato
2148 i = theInteractions.begin();
2149 while ( i != theInteractions.end() ) /* Lo
2150 {
2151 G4InteractionContent* anIniteraction = *i
2152 G4PartonPair * aPair = NULL;
2153
2154 #ifdef debugQGSParticipants
2155 G4cout<<"An interaction # and soft coll
2156 <<anIniteraction->GetNumberOfSoft
2157 IntNo++;
2158 #endif
2159 if (anIniteraction->GetNumberOfSoftCollis
2160 {
2161 G4VSplitableHadron* pProjectile = anIni
2162 G4VSplitableHadron* pTarget = anIni
2163
2164 for (G4int j = 0; j < anIniteraction->G
2165 {
2166 aPair = new G4PartonPair(pTarget->Get
2167 G4PartonPair::SOFT, G4PartonPair::TA
2168 #ifdef debugQGSParticipants
2169 G4cout << "SoftPair " << aPair->GetParton
2170 << aPair->GetParton1()->Get4Momentum() <
2171 << aPair->GetParton1()->Get4Momentum().m
2172 G4cout << " " << aPair->GetParton
2173 << aPair->GetParton2()->Get4Momentum() <
2174 <<aPair->GetParton2()->Get4Momentum().m
2175 str4Mom += aPair->GetParton1()->Get4Momen
2176 str4Mom += aPair->GetParton2()->Get4Momen
2177 #endif
2178
2179 thePartonPairs.push_back(aPair);
2180
2181 aPair = new G4PartonPair(pProjectile->GetNe
2182 G4PartonPair::SOFT, G4PartonPair::PR
2183 #ifdef debugQGSParticipants
2184 G4cout << "SoftPair " << aPair->GetParton
2185 << aPair->GetParton1()->Get4Momentum() <
2186 << aPair->GetParton1()->Get4Moment
2187 G4cout << " " << aPair->GetParton
2188 << aPair->GetParton2()->Get4Momentum() <
2189 << aPair->GetParton2()->Get4Momentum().m
2190 #endif
2191 #ifdef debugQGSParticipants
2192 str4Mom += aPair->GetParton1()->Get4Momen
2193 str4Mom += aPair->GetParton2()->Get4Momen
2194 #endif
2195
2196 thePartonPairs.push_back(aPair);
2197 }
2198
2199 delete *i;
2200 i=theInteractions.erase(i); // i now
2201 }
2202 else
2203 {
2204 i++;
2205 }
2206 } // End of while ( i != theInteractions.e
2207 #ifdef debugQGSParticipants
2208 G4cout << "Sum of strings 4 momenta " <<
2209 #endif
2210 }
2211
2212 //===========================================
2213
2214 void G4QGSParticipants::GetResiduals() {
2215 // This method is needed for the correct ap
2216
2217 #ifdef debugQGSParticipants
2218 G4cout << "GetResiduals(): GetProjectileN
2219 #endif
2220
2221 #ifdef debugQGSParticipants
2222 G4cout << "NumberOfInvolvedNucleonsOfTarg
2223 #endif
2224
2225 G4double DeltaExcitationE = TargetResidualE
2226 G4LorentzVector DeltaPResidualNucleus = Tar
2227 G4d
2228
2229 for ( G4int i = 0; i < NumberOfInvolvedNucl
2230 G4Nucleon* aNucleon = TheInvolvedNucleons
2231
2232 #ifdef debugQGSParticipants
2233 G4VSplitableHadron* targetSplitable = a
2234 G4cout << i << " Hit? " << aNucleon->Ar
2235 if ( targetSplitable ) G4cout << i << "
2236 #endif
2237
2238 G4LorentzVector tmp = -DeltaPResidualNucl
2239 aNucleon->SetMomentum( tmp );
2240 aNucleon->SetBindingEnergy( DeltaExcitati
2241 }
2242
2243 //-------------------------------------
2244 if( TargetResidualMassNumber != 0 )
2245 {
2246 G4ThreeVector bstToCM =TargetResidual4Mom
2247
2248 G4V3DNucleus* theTargetNucleus = GetTarge
2249 G4LorentzVector residualMomentum(0.,0.,0.
2250 G4Nucleon* aNucleon = 0;
2251 theTargetNucleus->StartLoop();
2252 while ( ( aNucleon = theTargetNucleus->Ge
2253 if ( !aNucleon->AreYouHit() ) {
2254 G4LorentzVector tmp=aNucleon->Get4Mom
2255 aNucleon->SetMomentum(tmp);
2256 residualMomentum +=tmp;
2257 }
2258 }
2259
2260 residualMomentum/=TargetResidualMassNumbe
2261
2262 G4double Mass = TargetResidual4Momentum.m
2263 G4double SumMasses=0.;
2264
2265 aNucleon = 0;
2266 theTargetNucleus->StartLoop();
2267 while ( ( aNucleon = theTargetNucleus->Ge
2268 if ( !aNucleon->AreYouHit() ) {
2269 G4LorentzVector tmp=aNucleon->Get4Mom
2270 G4double E=std::sqrt(tmp.vect().mag2(
2271 sqr(aNucleon->Ge
2272 tmp.setE(E); aNucleon->SetMomentum(t
2273 SumMasses+=E;
2274 }
2275 }
2276
2277 G4double Chigh=Mass/SumMasses; G4double C
2278 const G4int maxNumberOfLoops = 1000;
2279 G4int loopCounter = 0;
2280 do
2281 {
2282 C=(Chigh+Clow)/2.;
2283
2284 SumMasses=0.;
2285 aNucleon = 0;
2286 theTargetNucleus->StartLoop();
2287 while ( ( aNucleon = theTargetNucleus->
2288 if ( !aNucleon->AreYouHit() ) {
2289 G4LorentzVector tmp=aNucleon->Get4M
2290 G4double E=std::sqrt(tmp.vect().mag
2291 sqr(aNucleon->
2292 SumMasses+=E;
2293 }
2294 }
2295
2296 if(SumMasses > Mass) {Chigh=C;}
2297 else {Clow =C;}
2298
2299 } while( (Chigh-Clow > 0.01) &&
2300 ++loopCounter < maxNumberOfLoops
2301 if ( loopCounter >= maxNumberOfLoops ) {
2302 #ifdef debugQGSParticipants
2303 G4cout <<"BAD situation: forced loop
2304 #endif
2305 // Perhaps there is something to set he
2306 } else {
2307 aNucleon = 0;
2308 theTargetNucleus->StartLoop();
2309 while ( ( aNucleon = theTargetNucleus->
2310 if ( !aNucleon->AreYouHit() ) {
2311 G4LorentzVector tmp=aNucleon->Get4M
2312 G4double E=std::sqrt(tmp.vect().mag
2313 sqr(aNucleon->
2314 tmp.setE(E); tmp.boost(-bstToCM);
2315 aNucleon->SetMomentum(tmp);
2316 }
2317 }
2318 }
2319
2320 } // End of if( TargetResidualMassNumber !
2321 //-------------------------------------
2322
2323 #ifdef debugQGSParticipants
2324 G4cout << "End GetResiduals -------------
2325 #endif
2326
2327 }
2328
2329
2330 //===========================================
2331 void G4QGSParticipants::PerformSoftCollisions
2332 {
2333 std::vector<G4InteractionContent*>::iterato
2334 G4LorentzVector str4Mom;
2335 i = theInteractions.begin();
2336 while ( i != theInteractions.end() ) /* Lo
2337 {
2338 G4InteractionContent* anIniteraction = *i
2339 G4PartonPair * aPair = NULL;
2340 if (anIniteraction->GetNumberOfSoftCollis
2341 {
2342 G4VSplitableHadron* pProjectile = anIni
2343 G4VSplitableHadron* pTarget = anIni
2344 for (G4int j = 0; j < anIniteraction->G
2345 {
2346 aPair = new G4PartonPair(pTarget->Get
2347 G4PartonPair::SOFT, G4PartonPair::TA
2348 #ifdef debugQGSParticipants
2349 G4cout << "SoftPair " << aPair->GetParton
2350 << aPair->GetParton1()->Get4Momentum() <
2351 << aPair->GetParton1()->GetX() << " " <<
2352 G4cout << " " << aPair->GetParton
2353 << aPair->GetParton2()->Get4Momentum() <
2354 << aPair->GetParton2()->GetX() << " " <<
2355 #endif
2356 #ifdef debugQGSParticipants
2357 str4Mom += aPair->GetParton1()->Get4Momen
2358 str4Mom += aPair->GetParton2()->Get4Momen
2359 #endif
2360 thePartonPairs.push_back(aPair);
2361 aPair = new G4PartonPair(pProjectile->GetNe
2362 G4PartonPair::SOFT, G4PartonPair::PR
2363 #ifdef debugQGSParticipants
2364 G4cout << "SoftPair " << aPair->GetParton
2365 << aPair->GetParton1()->Get4Momentum() <
2366 << aPair->GetParton1()->GetX() << " " <<
2367 G4cout << " " << aPair->GetParton
2368 << aPair->GetParton2()->Get4Momentum() <
2369 << aPair->GetParton2()->GetX() << " " <<
2370 #endif
2371 #ifdef debugQGSParticipants
2372 str4Mom += aPair->GetParton1()->Get4Momen
2373 str4Mom += aPair->GetParton2()->Get4Momen
2374 #endif
2375 thePartonPairs.push_back(aPair);
2376 }
2377 delete *i;
2378 i=theInteractions.erase(i); // i now
2379 } else {
2380 i++;
2381 }
2382 }
2383 #ifdef debugQGSParticipants
2384 G4cout << " string 4 mom " << str4Mom <<
2385 #endif
2386 }
2387
2388
2389 //*******************************************
2390 G4VSplitableHadron* G4QGSParticipants::Select
2391 {
2392 // Check reaction threshold - goes to Chec
2393
2394 theProjectileSplitable = new G4QGSMSplitabl
2395 theProjectileSplitable->SetStatus(1);
2396
2397 G4LorentzVector aPrimaryMomentum(thePrimary
2398 G4LorentzVector aTargetNMomentum(0.,0.,0.,9
2399
2400 if ((!(aPrimaryMomentum.e()>-1)) && (!(aPri
2401 {
2402 throw G4HadronicException(__FILE__, __LIN
2403 "G4GammaParticipants::SelectInter
2404 }
2405 G4double S = (aPrimaryMomentum + aTargetNMo
2406 G4double ThresholdMass = thePrimary.GetMass
2407 ModelMode = SOFT;
2408
2409 #ifdef debugQGSParticipants
2410 G4cout << "Gamma projectile - SelectInter
2411 G4cout<<"Energy and Nucleus "<<thePrimary
2412 G4cout << "SqrtS ThresholdMass ModelMode
2413 #endif
2414
2415 if (sqr(ThresholdMass + ThresholdParameter)
2416 {
2417 ModelMode = DIFFRACTIVE;
2418 }
2419 if (sqr(ThresholdMass + QGSMThreshold) > S)
2420 {
2421 ModelMode = DIFFRACTIVE;
2422 }
2423
2424 // first find the collisions HPW
2425 std::for_each(theInteractions.begin(), theI
2426 theInteractions.clear();
2427
2428 G4int theCurrent = G4int(theNucleus->GetMas
2429 G4int NucleonNo=0;
2430
2431 theNucleus->StartLoop();
2432 G4Nucleon * pNucleon = 0;
2433
2434 while( (pNucleon = theNucleus->GetNextNucle
2435 { if(NucleonNo == theCurrent) break; Nucleo
2436
2437 if ( pNucleon ) {
2438 G4QGSMSplitableHadron* aTarget = new G4QG
2439
2440 pNucleon->Hit(aTarget);
2441
2442 if ( (0.06 > G4UniformRand() &&(ModelMode
2443 {
2444 G4InteractionContent * aInteraction = n
2445 theProjectileSplitable->SetStatus(1*the
2446
2447 aInteraction->SetTarget(aTarget);
2448 aInteraction->SetTargetNucleon(pNucleon
2449 aTarget->SetCollisionCount(0);
2450 aTarget->SetStatus(1);
2451
2452 aInteraction->SetNumberOfDiffractiveCol
2453 aInteraction->SetNumberOfSoftCollisions
2454 aInteraction->SetStatus(1);
2455
2456 theInteractions.push_back(aInteraction)
2457 }
2458 else
2459 {
2460 // nondiffractive soft interaction occu
2461 aTarget->IncrementCollisionCount(1);
2462 aTarget->SetStatus(0);
2463 theTargets.push_back(aTarget);
2464
2465 theProjectileSplitable->IncrementCollis
2466 theProjectileSplitable->SetStatus(0*the
2467
2468 G4InteractionContent * aInteraction = n
2469 aInteraction->SetTarget(aTarget);
2470 aInteraction->SetTargetNucleon(pNucleon
2471 aInteraction->SetNumberOfSoftCollisions
2472 aInteraction->SetStatus(3);
2473 theInteractions.push_back(aInteraction)
2474 }
2475 }
2476 return theProjectileSplitable;
2477 }
2478