| Geant4 Cross Reference |
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 #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 variables class properties
47 G4ThreadLocal G4int G4QGSParticipants_NPart = 0;
48
49 G4QGSParticipants::G4QGSParticipants()
50 : theDiffExcitaton(), ModelMode(SOFT), nCutMax(7),
51 ThresholdParameter(0.0*GeV), QGSMThreshold(3.0*GeV),
52 theNucleonRadius(1.5*fermi), theCurrentVelocity(G4ThreeVector()),
53 theProjectileSplitable(nullptr), Regge(nullptr),
54 InteractionMode(ALL), alpha(-0.5), beta(2.5), sigmaPt(0.0),
55 NumberOfInvolvedNucleonsOfTarget(0), NumberOfInvolvedNucleonsOfProjectile(0),
56 ProjectileResidual4Momentum(G4LorentzVector()), ProjectileResidualMassNumber(0),
57 ProjectileResidualCharge(0), ProjectileResidualExcitationEnergy(0.0),
58 TargetResidual4Momentum(G4LorentzVector()), TargetResidualMassNumber(0),
59 TargetResidualCharge(0), TargetResidualExcitationEnergy(0.0),
60 CofNuclearDestruction(0.0), R2ofNuclearDestruction(0.0),
61 ExcitationEnergyPerWoundedNucleon(0.0), DofNuclearDestruction(0.0),
62 Pt2ofNuclearDestruction(0.0), MaxPt2ofNuclearDestruction(0.0)
63 {
64 for (size_t i=0; i < 250; ++i) {
65 TheInvolvedNucleonsOfTarget[i] = nullptr;
66 TheInvolvedNucleonsOfProjectile[i] = nullptr;
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*MeV );
75
76 sigmaPt=0.25*sqr(GeV);
77 }
78
79 G4QGSParticipants::G4QGSParticipants(const G4QGSParticipants &right)
80 : G4VParticipants(), ModelMode(right.ModelMode), nCutMax(right.nCutMax),
81 ThresholdParameter(right.ThresholdParameter),
82 QGSMThreshold(right.QGSMThreshold),
83 theNucleonRadius(right.theNucleonRadius),
84 theCurrentVelocity(right.theCurrentVelocity),
85 theProjectileSplitable(right.theProjectileSplitable),
86 Regge(right.Regge), InteractionMode(right.InteractionMode),
87 alpha(right.alpha), beta(right.beta), sigmaPt(right.sigmaPt),
88 NumberOfInvolvedNucleonsOfTarget(right.NumberOfInvolvedNucleonsOfTarget),
89 NumberOfInvolvedNucleonsOfProjectile(right.NumberOfInvolvedNucleonsOfProjectile),
90 ProjectileResidual4Momentum(right.ProjectileResidual4Momentum),
91 ProjectileResidualMassNumber(right.ProjectileResidualMassNumber),
92 ProjectileResidualCharge(right.ProjectileResidualCharge),
93 ProjectileResidualExcitationEnergy(right.ProjectileResidualExcitationEnergy),
94 TargetResidual4Momentum(right.TargetResidual4Momentum),
95 TargetResidualMassNumber(right.TargetResidualMassNumber),
96 TargetResidualCharge(right.TargetResidualCharge),
97 TargetResidualExcitationEnergy(right.TargetResidualExcitationEnergy),
98 CofNuclearDestruction(right.CofNuclearDestruction),
99 R2ofNuclearDestruction(right.R2ofNuclearDestruction),
100 ExcitationEnergyPerWoundedNucleon(right.ExcitationEnergyPerWoundedNucleon),
101 DofNuclearDestruction(right.DofNuclearDestruction),
102 Pt2ofNuclearDestruction(right.Pt2ofNuclearDestruction),
103 MaxPt2ofNuclearDestruction(right.MaxPt2ofNuclearDestruction)
104 {
105 for (size_t i=0; i < 250; ++i) {
106 TheInvolvedNucleonsOfTarget[i] = right.TheInvolvedNucleonsOfTarget[i];
107 TheInvolvedNucleonsOfProjectile[i] = right.TheInvolvedNucleonsOfProjectile[i];
108 }
109 }
110
111 G4QGSParticipants::~G4QGSParticipants() {}
112
113 void G4QGSParticipants::BuildInteractions(const G4ReactionProduct &thePrimary)
114 {
115 theProjectile = thePrimary;
116
117 Regge = new G4Reggeons(theProjectile.GetDefinition());
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->GetMassNumber();
130 TargetResidualCharge = theNucleus->GetCharge();
131 TargetResidualExcitationEnergy = 0.0;
132
133 theNucleus->StartLoop();
134 G4Nucleon* NuclearNucleon;
135 while ( ( NuclearNucleon = theNucleus->GetNextNucleon() ) ) {
136 tmp+=NuclearNucleon->Get4Momentum();
137 }
138 TargetResidual4Momentum = tmp;
139
140 if ( std::abs( theProjectile.GetDefinition()->GetBaryonNumber() ) <= 1 ) {
141 // Projectile is a hadron : meson or baryon
142 ProjectileResidualMassNumber = std::abs( theProjectile.GetDefinition()->GetBaryonNumber() );
143 ProjectileResidualCharge = G4int( theProjectile.GetDefinition()->GetPDGCharge() );
144 ProjectileResidualExcitationEnergy = 0.0;
145 ProjectileResidual4Momentum.setVect( theProjectile.GetMomentum() );
146 ProjectileResidual4Momentum.setE( theProjectile.GetTotalEnergy() );
147 }
148
149 #ifdef debugQGSParticipants
150 G4cout <<G4endl<< "G4QGSParticipants::BuildInteractions---------" << G4endl
151 << "thePrimary " << thePrimary.GetDefinition()->GetParticleName()<<" "
152 <<ProjectileResidual4Momentum<<G4endl;
153 G4cout << "Target A and Z " << theNucleus->GetMassNumber() <<" "<<theNucleus->GetCharge()<<" "
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 = 1000;
164 G4int internalLoopCounter = 0;
165 do
166 {
167 if(std::abs(theProjectile.GetDefinition()->GetPDGEncoding()) < 100.0)
168 {
169 SelectInteractions(theProjectile); // for lepton projectile
170 }
171 else
172 {
173 GetList( theProjectile ); // Get list of participating nucleons for hadron projectile
174 }
175
176 if ( theInteractions.size() == 0 ) return;
177
178 StoreInvolvedNucleon(); // Store participating nucleon
179
180 ReggeonCascade(); // Make reggeon cascading. Involve nucleons in the cascading.
181
182 Success = PutOnMassShell(); // Ascribe momenta to the involved and participating nucleons
183
184 if(!Success) PrepareInitialState( thePrimary );
185
186 } while( (!Success) && ++internalLoopCounter < maxNumberOfInternalLoops );
187
188 if ( internalLoopCounter >= maxNumberOfInternalLoops ) {
189 Success = false;
190 }
191
192 if ( Success ) {
193 #ifdef debugQGSParticipants
194 G4cout<<G4endl<<"PerformDiffractiveCollisions(); if they happend." <<G4endl;
195 #endif
196
197 PerformDiffractiveCollisions();
198
199 #ifdef debugQGSParticipants
200 G4cout<<G4endl<<"SplitHadrons();" <<G4endl;
201 #endif
202
203 SplitHadrons();
204
205 if( theProjectileSplitable && theProjectileSplitable->GetStatus() == 0) {
206 #ifdef debugQGSParticipants
207 G4cout<<"Perform non-Diffractive Collisions if they happend. Determine Parton Momenta and so on." <<G4endl;
208 #endif
209 Success = DeterminePartonMomenta();
210 }
211
212 if(!Success) PrepareInitialState( thePrimary );
213 }
214 } while( (!Success) && ++loopCounter < maxNumberOfLoops );
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 nucleus properties
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() "<<theInteractions.size()<<G4endl;
237 #endif
238
239 if( Regge ) delete Regge;
240
241 std::for_each( theInteractions.begin(), theInteractions.end(), DeleteInteractionContent() );
242 theInteractions.clear();
243
244 // Erasing of target involved nucleons.
245 #ifdef debugQGSParticipants
246 G4cout<<"Erasing of involved target nucleons "<<NumberOfInvolvedNucleonsOfTarget<<G4endl;
247 #endif
248
249 if ( NumberOfInvolvedNucleonsOfTarget != 0 ) {
250 for ( G4int i = 0; i < NumberOfInvolvedNucleonsOfTarget; i++ ) {
251 G4VSplitableHadron* aNucleon = TheInvolvedNucleonsOfTarget[i]->GetSplitableHadron();
252 if ( (aNucleon != 0 ) && (aNucleon->GetStatus() >= 1) ) delete aNucleon;
253 }
254 }
255
256 // Erasing of projectile involved nucleons.
257 if ( NumberOfInvolvedNucleonsOfProjectile != 0 ) {
258 for ( G4int i = 0; i < NumberOfInvolvedNucleonsOfProjectile; i++ ) {
259 G4VSplitableHadron* aNucleon = TheInvolvedNucleonsOfProjectile[i]->GetSplitableHadron();
260 if ( aNucleon ) delete aNucleon;
261 }
262 }
263
264 #ifdef debugQGSParticipants
265 G4cout<<"Delition of target nucleons from soft interactions "<<theTargets.size()
266 <<G4endl<<G4endl;
267 #endif
268 std::for_each(theTargets.begin(), theTargets.end(), DeleteSplitableHadron());
269 theTargets.clear();
270
271 if ( theProjectileSplitable ) {
272 delete theProjectileSplitable;
273 theProjectileSplitable = 0;
274 }
275 }
276
277 //===========================================================
278 void G4QGSParticipants::PrepareInitialState( const G4ReactionProduct& thePrimary )
279 {
280 // Clearing of the arrays
281 // Erasing of the projectile
282 G4InteractionContent* anIniteraction = theInteractions[0];
283 G4VSplitableHadron* pProjectile = anIniteraction->GetProjectile();
284 if( pProjectile ) delete pProjectile;
285
286 std::for_each(theInteractions.begin(), theInteractions.end(), DeleteInteractionContent());
287 theInteractions.clear();
288
289 // Erasing of the envolved nucleons and target nucleons from diffraction dissociations
290 theNucleus->StartLoop();
291 G4Nucleon* aNucleon;
292 while ( ( aNucleon = theNucleus->GetNextNucleon() ) )
293 {
294 if ( aNucleon->AreYouHit() ) {
295 G4VSplitableHadron* splaNucleon = aNucleon->GetSplitableHadron();
296 if ( (splaNucleon != 0) && (splaNucleon->GetStatus() >=1) ) delete splaNucleon;
297 aNucleon->Hit(nullptr);
298 NumberOfInvolvedNucleonsOfTarget--;
299 }
300 }
301
302 // Erasing of nuclear nucleons participated in soft interactions
303 std::for_each(theTargets.begin(), theTargets.end(), DeleteSplitableHadron());
304 theTargets.clear();
305
306 // Preparation to a new attempt
307 theProjectile = thePrimary;
308
309 theNucleus->Init(theNucleus->GetMassNumber(), theNucleus->GetCharge());
310 theNucleus->SortNucleonsIncZ();
311 DoLorentzBoost(-theCurrentVelocity); // Lorentz boost of the target nucleus
312
313 if (theNucleus->GetMassNumber() == 1)
314 {
315 G4ThreeVector aPos = G4ThreeVector(0.,0.,0.);
316 theNucleus->StartLoop();
317 G4Nucleon* tNucleon=theNucleus->GetNextNucleon();
318 tNucleon->SetPosition(aPos);
319 }
320
321 G4LorentzVector Tmp( 0.0, 0.0, 0.0, 0.0 );
322 NumberOfInvolvedNucleonsOfTarget= 0;
323 TargetResidualMassNumber = theNucleus->GetMassNumber();
324 TargetResidualCharge = theNucleus->GetCharge();
325 TargetResidualExcitationEnergy = 0.0;
326
327 G4Nucleon* NuclearNucleon;
328 while ( ( NuclearNucleon = theNucleus->GetNextNucleon() ) )
329 {Tmp+=NuclearNucleon->Get4Momentum();}
330
331 TargetResidual4Momentum = Tmp;
332 }
333
334 //===========================================================
335 void G4QGSParticipants::GetList( const G4ReactionProduct& thePrimary ) {
336 #ifdef debugQGSParticipants
337 G4cout<<G4endl<<"G4QGSParticipants::GetList +++++++++++++"<<G4endl;
338 #endif
339
340 // Direction: True - Proj, False - Target
341 theProjectileSplitable = new G4QGSMSplitableHadron(thePrimary, TRUE);
342 theProjectileSplitable->SetStatus(1);
343
344 G4LorentzVector aPrimaryMomentum(thePrimary.GetMomentum(), thePrimary.GetTotalEnergy());
345 G4LorentzVector aNucleonMomentum(0.,0.,0., 938.0*MeV);
346
347 G4double SS=(aPrimaryMomentum + aNucleonMomentum).mag2();
348
349 Regge->SetS(SS);
350
351 //--------------------------------------
352 theNucleus->StartLoop();
353 G4Nucleon * tNucleon = theNucleus->GetNextNucleon();
354
355 if ( ! tNucleon ) {
356 #ifdef debugQGSParticipants
357 G4cout << "QGSM - BAD situation: pNucleon is NULL ! Leaving immediately!" << G4endl;
358 #endif
359 return;
360 }
361
362 G4double theNucleusOuterR = theNucleus->GetOuterRadius();
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 of nucleus ------------------------------------
372
373 std::for_each(theInteractions.begin(), theInteractions.end(), DeleteInteractionContent());
374 theInteractions.clear();
375
376 G4int MaxPower=thePrimary.GetMomentum().mag()/(3.3*GeV); if(MaxPower < 1) MaxPower=1;
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) && ++loopCounter < maxNumberOfLoops )
385 {
386 InteractionMode = ALL; // Mode = ALL, WITHOUT_R, NON_DIFF
387
388 // choose random impact parameter of a collision
389 std::pair<G4double, G4double> theImpactParameter;
390
391 NumberOfTries++;
392 if( NumberOfTries == 100*(NumberOfTries/100) ) Scale /=2.0;
393
394 theImpactParameter = theNucleus->ChooseImpactXandY(theNucleusOuterR/Scale + theNucleonRadius);
395 G4double impactX = theImpactParameter.first;
396 G4double impactY = theImpactParameter.second;
397
398 #ifdef debugQGSParticipants
399 G4cout<<"InteractionMode "<<InteractionMode<<G4endl;
400 G4cout<<"Impact parameter (fm ) "<<std::sqrt(sqr(impactX)+sqr(impactY))/fermi<<" "<<G4endl;
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->GetNextNucleon()) )
411 {
412 if(Power <= 0.) break;
413
414 G4LorentzVector nucleonMomentum=tNucleon->Get4Momentum();
415
416 G4double Distance2 = sqr(impactX - tNucleon->GetPosition().x()) +
417 sqr(impactY - tNucleon->GetPosition().y());
418
419 G4double Pint(0.); // A probability of interaction at given impact parameter
420 G4double Pprd(0.), Ptrd(0.), Pdd(0.); // Probabilities of Proj. diffr., Target diffr., Double diffr.
421 G4double Pnd (0.), Pnvr(0.); // Probabilities of non-diffr. and quark exchange
422 G4int NcutPomerons(0); // Number of cutted pomerons
423
424 Regge->GetProbabilities(std::sqrt(Distance2), InteractionMode,
425 Pint, Pprd, Ptrd, Pdd, Pnd, Pnvr);
426 #ifdef debugQGSParticipants
427 nucleonCount++;
428 G4cout<<"Nucleon & its impact parameter: "<<nucleonCount<<" "<<std::sqrt(Distance2)/fermi<<" (fm)"<<G4endl;
429 G4cout<<"Probability of interaction: "<<Pint<<G4endl;
430 G4cout<<"Probability of PrD, TrD, DD: "<<Pprd<<" "<<Ptrd<<" "<<Pdd<<G4endl;
431 G4cout<<"Probability of NonDiff, QuarkExc.: "<<Pnd<<" "<<Pnvr<<" in inel. inter."<<G4endl;
432 #endif
433
434 if (Pint > G4UniformRand())
435 { // An interaction is happend.
436
437 G4double rndNumber = G4UniformRand();
438 G4int InteractionType(0);
439
440 if((InteractionMode==ALL)||(InteractionMode==WITHOUT_R)) // Mode = ALL, WITHOUT_R, NON_DIFF
441 {
442 if( rndNumber < Pprd ) {InteractionType = PrD; InteractionMode = WITHOUT_R;}
443 else if( rndNumber < Pprd+Ptrd) {InteractionType = TrD; InteractionMode = WITHOUT_R;}
444 else if( rndNumber < Pprd+Ptrd+Pdd) {InteractionType = DD; InteractionMode = WITHOUT_R;}
445 else if( rndNumber < Pprd+Ptrd+Pdd+Pnd ) {InteractionType = NonD; InteractionMode = NON_DIFF;
446 NcutPomerons = Regge->ncPomerons(); }
447 else {InteractionType = Qexc; InteractionMode = ALL; }
448 }
449 else // InteractionMode == NON_DIFF
450 {
451 InteractionMode = NON_DIFF;
452 if( rndNumber < Ptrd ) {InteractionType = TrD; }
453 else if( rndNumber < Ptrd + Pnd) {InteractionType = NonD; NcutPomerons = Regge->ncPomerons();}
454 }
455
456 if( (InteractionType == NonD) && (NcutPomerons == 0)) continue;
457
458 G4QGSParticipants_NPart ++;
459 G4QGSMSplitableHadron* aTargetSPB = new G4QGSMSplitableHadron(*tNucleon);
460 tNucleon->Hit(aTargetSPB);
461
462 #ifdef debugQGSParticipants
463 G4cout<<"An interaction is happend."<<G4endl;
464 G4cout<<"Target nucleon - "<<nucleonCount<<" "
465 <<tNucleon->GetDefinition()->GetParticleName()<<G4endl;
466 G4cout<<"Interaction type:"<<InteractionType
467 <<" (0 -PrD, 1 - TrD, 2 - DD, 3 - NonD, 4 - Qexc)"<<G4endl;
468 G4cout<<"New Inter. mode:"<<InteractionMode
469 <<" (0 -ALL, 1 - WITHOUT_R, 2 - NON_DIFF)"<<G4endl;
470 if( InteractionType == NonD )
471 G4cout<<"Number of cutted pomerons: "<<NcutPomerons<<G4endl;
472 #endif
473
474 if((InteractionType == PrD) || (InteractionType == TrD) || (InteractionType == DD) ||
475 (InteractionType == Qexc))
476 { // diffractive-like interaction occurs
477 #ifdef debugQGSParticipants
478 G4cout<<"Diffractive-like interaction occurs"<<G4endl;
479 #endif
480
481 G4InteractionContent * aInteraction = new G4InteractionContent(theProjectileSplitable);
482 theProjectileSplitable->SetStatus(1*theProjectileSplitable->GetStatus());
483
484 aInteraction->SetTarget(aTargetSPB);
485 aInteraction->SetTargetNucleon(tNucleon);
486 aTargetSPB->SetCollisionCount(0);
487 aTargetSPB->SetStatus(1);
488
489 aInteraction->SetNumberOfDiffractiveCollisions(1);
490 aInteraction->SetNumberOfSoftCollisions(0);
491 aInteraction->SetStatus(InteractionType);
492 theInteractions.push_back(aInteraction);
493 }
494 else
495 { // nondiffractive interaction occurs
496 #ifdef debugQGSParticipants
497 G4cout<<"Non-diffractive interaction occurs, max NcutPomerons "<<NcutPomerons<<G4endl;
498 #endif
499
500 G4int nCuts;
501
502 G4int Vncut=0;
503 for(nCuts = 0; nCuts < NcutPomerons; nCuts++)
504 {
505 if( G4UniformRand() < Power/MaxPower ){Vncut++; Power--; if(Power <= 0.) break;}
506 }
507 nCuts=Vncut;
508
509 if( nCuts == 0 ) {delete aTargetSPB; tNucleon->Hit(nullptr); continue;}
510
511 #ifdef debugQGSParticipants
512 G4cout<<"Number of cuts in the interaction "<<nCuts<<G4endl;
513 #endif
514
515 aTargetSPB->IncrementCollisionCount(nCuts);
516 aTargetSPB->SetStatus(0);
517 theTargets.push_back(aTargetSPB);
518
519 theProjectileSplitable->IncrementCollisionCount(nCuts);
520 theProjectileSplitable->SetStatus(0*theProjectileSplitable->GetStatus());
521
522 G4InteractionContent * aInteraction = new G4InteractionContent(theProjectileSplitable);
523 aInteraction->SetTarget(aTargetSPB);
524 aInteraction->SetTargetNucleon(tNucleon);
525 aInteraction->SetNumberOfSoftCollisions(nCuts);
526 aInteraction->SetStatus(InteractionType);
527 theInteractions.push_back(aInteraction);
528 }
529 } // End of if (Pint > G4UniformRand())
530 } // End of while( (tNucleon = theNucleus->GetNextNucleon()) )
531
532 #ifdef debugQGSParticipants
533 G4cout << G4endl<<"Number of wounded nucleons "<<G4QGSParticipants_NPart<<G4endl;
534 #endif
535
536 } // End of while( (theInteractions.size() == 0) && ++loopCounter < maxNumberOfLoops )
537
538 if ( loopCounter >= maxNumberOfLoops ) {
539 #ifdef debugQGSParticipants
540 G4cout <<"BAD situation: forced loop exit!" << G4endl;
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 interactions ??
546 }
547 //------------------------------------------------------------
548 std::vector<G4InteractionContent*>::iterator i;
549
550 if( theInteractions.size() != 0)
551 {
552 if( InteractionMode == ALL ) // It can be if all interactions were quark-exchange.
553 { // Only the first one will be saved, all other will be erased.
554 i = theInteractions.end()-1;
555
556 while ( theInteractions.size() != 1 )
557 {
558 G4InteractionContent* anInteraction = *i;
559 G4Nucleon * pNucleon = anInteraction->GetTargetNucleon(); pNucleon->Hit(nullptr);
560 delete anInteraction->GetTarget();
561 delete *i;
562 i=theInteractions.erase(i);
563 i--;
564 }
565 }
566 else
567 { // All quark exchanges will be erased
568 i = theInteractions.begin();
569 while ( i != theInteractions.end() )
570 {
571 G4InteractionContent* anInteraction = *i;
572
573 if( anInteraction->GetStatus() == Qexc )
574 {
575 G4Nucleon* aTargetNucleon = anInteraction->GetTargetNucleon();
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 interaction
594
595 NumberOfInvolvedNucleonsOfTarget = 0;
596
597 theNucleus->StartLoop();
598
599 G4Nucleon* aNucleon;
600 while ( ( aNucleon = theNucleus->GetNextNucleon() ) ) {
601 if ( aNucleon->AreYouHit() ) {
602 TheInvolvedNucleonsOfTarget[NumberOfInvolvedNucleonsOfTarget] = aNucleon;
603 NumberOfInvolvedNucleonsOfTarget++;
604 }
605 }
606
607 #ifdef debugQGSParticipants
608 G4cout << G4endl<<"G4QGSParticipants::StoreInvolvedNucleon() if they were "<<G4endl
609 <<"Stored # of wounded nucleons of target "
610 << NumberOfInvolvedNucleonsOfTarget <<G4endl;
611 #endif
612 return;
613 }
614
615 //=============================================================
616
617 void G4QGSParticipants::ReggeonCascade()
618 { // Implementation of the reggeon theory inspired model of nuclear destruction
619 #ifdef debugQGSParticipants
620 G4cout << G4endl<<"Reggeon cascading ........."<<G4endl;
621 G4cout<<"C of nucl. desctruction "<<GetCofNuclearDestruction()
622 <<" R2 "<<GetR2ofNuclearDestruction()/fermi/fermi<<" fermi^2"<<G4endl;
623 #endif
624
625 G4int InitNINt = NumberOfInvolvedNucleonsOfTarget;
626
627 // Reggeon cascading in target nucleus
628 for ( G4int InvTN = 0; InvTN < InitNINt; InvTN++ ) {
629 G4Nucleon* aTargetNucleon = TheInvolvedNucleonsOfTarget[ InvTN ];
630
631 G4double CreationTime = aTargetNucleon->GetSplitableHadron()->GetTimeOfCreation();
632
633 G4double XofWoundedNucleon = aTargetNucleon->GetPosition().x();
634 G4double YofWoundedNucleon = aTargetNucleon->GetPosition().y();
635
636 G4V3DNucleus* theTargetNucleus = theNucleus;
637 theTargetNucleus->StartLoop();
638
639 #ifdef debugQGSParticipants
640 G4int TrgNuc=0;
641 #endif
642 G4Nucleon* Neighbour(0);
643 while ( ( Neighbour = theTargetNucleus->GetNextNucleon() ) ) {
644 #ifdef debugQGSParticipants
645 TrgNuc++;
646 #endif
647 if ( ! Neighbour->AreYouHit() ) {
648 G4double impact2 = sqr( XofWoundedNucleon - Neighbour->GetPosition().x() ) +
649 sqr( YofWoundedNucleon - Neighbour->GetPosition().y() );
650
651 if ( G4UniformRand() < GetCofNuclearDestruction() *
652 G4Exp( -impact2 / GetR2ofNuclearDestruction() )
653 ) {
654 // The neighbour nucleon is involved in the reggeon cascade
655 #ifdef debugQGSParticipants
656 G4cout<<"Target nucleon involved in reggeon cascading No "<<TrgNuc<<" "
657 <<Neighbour->GetDefinition()->GetParticleName()<<G4endl;
658 #endif
659 TheInvolvedNucleonsOfTarget[ NumberOfInvolvedNucleonsOfTarget ] = Neighbour;
660 NumberOfInvolvedNucleonsOfTarget++;
661
662 G4QGSMSplitableHadron* targetSplitable = new G4QGSMSplitableHadron( *Neighbour );
663
664 Neighbour->Hit( targetSplitable );
665 targetSplitable->SetTimeOfCreation( CreationTime );
666 targetSplitable->SetStatus( 2 );
667 targetSplitable->SetCollisionCount(0);
668
669 G4InteractionContent * anInteraction = new G4InteractionContent(theProjectileSplitable);
670 anInteraction->SetTarget(targetSplitable);
671 anInteraction->SetTargetNucleon(Neighbour);
672
673 anInteraction->SetNumberOfDiffractiveCollisions(1);
674 anInteraction->SetNumberOfSoftCollisions(0);
675 anInteraction->SetStatus(3);
676 theInteractions.push_back(anInteraction);
677 }
678 }
679 }
680 }
681
682 #ifdef debugQGSParticipants
683 G4cout <<"Number of new involved nucleons "<<NumberOfInvolvedNucleonsOfTarget - InitNINt<<G4endl;
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 ..............." << G4endl;
699 if ( isProjectileNucleus ) {G4cout << "PutOnMassShell for Nucleus_Nucleus " << G4endl;}
700 #endif
701
702 G4LorentzVector Pprojectile( theProjectile.GetMomentum(), theProjectile.GetTotalEnergy() );
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.0, 0.0 );
711 G4double SumMasses = 0.0;
712 G4V3DNucleus* theTargetNucleus = GetTargetNucleus();
713 G4double TargetResidualMass = 0.0;
714
715 #ifdef debugPutOnMassShell
716 G4cout << "Target : ";
717 #endif
718
719 isOk = ComputeNucleusProperties( theTargetNucleus, Ptarget, PtargetResidual, SumMasses,
720 TargetResidualExcitationEnergy, TargetResidualMass,
721 TargetResidualMassNumber, TargetResidualCharge );
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, 0.0 );
729 G4V3DNucleus* thePrNucleus = GetProjectileNucleus();
730 G4double PrResidualMass = 0.0;
731
732 if ( ! isProjectileNucleus ) { // hadron-nucleus collision
733 Mprojectile = Pprojectile.mag();
734 M2projectile = Pprojectile.mag2();
735 SumMasses += Mprojectile + 20.0*MeV; // Maybe DM must be larger?
736 } else { // nucleus-nucleus or antinucleus-nucleus collision
737
738 #ifdef debugPutOnMassShell
739 G4cout << "Projectile : ";
740 #endif
741
742 isOk = ComputeNucleusProperties( thePrNucleus, Pproj, PprojResidual, SumMasses,
743 ProjectileResidualExcitationEnergy, PrResidualMass,
744 ProjectileResidualMassNumber, ProjectileResidualCharge );
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" << G4endl << "SqrtS " << SqrtS/GeV << " GeV" << G4endl
756 << "SumMasses, PrResidualMass and TargetResidualMass " << SumMasses/GeV << " "
757 << PrResidualMass/GeV << " " << TargetResidualMass/GeV << " GeV" << G4endl;
758 G4cout << "Ptar res. "<<PtargetResidual<<G4endl;
759 #endif
760
761 if ( SqrtS < SumMasses ) {
762 return false; // It is impossible to simulate after putting nuclear nucleons on mass-shell.
763 }
764
765 // Try to consider also the excitation energy of the residual nucleus, if this is
766 // possible, with the available energy; otherwise, set the excitation energy to zero.
767
768 G4double savedSumMasses = SumMasses;
769 if ( isProjectileNucleus ) {
770 SumMasses -= std::sqrt( sqr( PrResidualMass ) + PprojResidual.perp2() );
771 SumMasses += std::sqrt( sqr( PrResidualMass + ProjectileResidualExcitationEnergy )
772 + PprojResidual.perp2() );
773 }
774 SumMasses -= std::sqrt( sqr( TargetResidualMass ) + PtargetResidual.perp2() );
775 SumMasses += std::sqrt( sqr( TargetResidualMass + TargetResidualExcitationEnergy )
776 + PtargetResidual.perp2() );
777
778 if ( SqrtS < SumMasses ) {
779 SumMasses = savedSumMasses;
780 if ( isProjectileNucleus ) {
781 ProjectileResidualExcitationEnergy = 0.0;
782 }
783 TargetResidualExcitationEnergy = 0.0;
784 }
785
786 TargetResidualMass += TargetResidualExcitationEnergy;
787 if ( isProjectileNucleus ) {
788 PrResidualMass += ProjectileResidualExcitationEnergy;
789 }
790
791 #ifdef debugPutOnMassShell
792 if ( isProjectileNucleus ) {
793 G4cout << "PrResidualMass ProjResidualExcitationEnergy " << PrResidualMass/GeV << " "
794 << ProjectileResidualExcitationEnergy << " MeV" << G4endl;
795 }
796 G4cout << "TargetResidualMass TargetResidualExcitationEnergy " << TargetResidualMass/GeV << " GeV "
797 << TargetResidualExcitationEnergy << " MeV" << G4endl
798 << "Sum masses " << SumMasses/GeV << G4endl;
799 #endif
800
801 // Sampling of nucleons what can transfer to delta-isobars
802 if ( isProjectileNucleus && thePrNucleus->GetMassNumber() != 1 ) {
803 isOk = GenerateDeltaIsobar( SqrtS, NumberOfInvolvedNucleonsOfProjectile,
804 TheInvolvedNucleonsOfProjectile, SumMasses );
805 }
806 if ( theTargetNucleus->GetMassNumber() != 1 ) {
807 isOk = isOk &&
808 GenerateDeltaIsobar( SqrtS, NumberOfInvolvedNucleonsOfTarget,
809 TheInvolvedNucleonsOfTarget, SumMasses );
810 }
811 if ( ! isOk ) return false;
812
813 // Now we know that it is kinematically possible to produce a final state made
814 // of the involved nucleons (or corresponding delta-isobars) and a residual nucleus.
815 // We have to sample the kinematical variables which will allow to define the 4-momenta
816 // of the final state. The sampled kinematical variables refer to the center-of-mass frame.
817 // Notice that the sampling of the transverse momentum corresponds to take into account
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" moving backwards in c.m.s., abort collision!
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 Xminus
839 G4double DcorP = 0.0;
840 if ( isProjectileNucleus ) {
841 DcorP = GetDofNuclearDestruction() / thePrNucleus->GetMassNumber();
842 }
843 G4double DcorT = GetDofNuclearDestruction() / theTargetNucleus->GetMassNumber();
844 G4double AveragePt2 = GetPt2ofNuclearDestruction();
845 G4double maxPtSquare = GetMaxPt2ofNuclearDestruction();
846
847 #ifdef debugPutOnMassShell
848 if ( isProjectileNucleus ) {
849 G4cout << "Y projectileNucleus " << YprojectileNucleus << G4endl;
850 }
851 G4cout << "Y targetNucleus " << YtargetNucleus << G4endl
852 << "Dcor " << GetDofNuclearDestruction()
853 << " DcorP DcorT " << DcorP << " " << DcorT << " AveragePt2 " << AveragePt2 << G4endl;
854 #endif
855
856 G4double M2proj = M2projectile; // Initialization needed only for hadron-nucleus collisions
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 = 0.0;
868 OuterSuccess = true;
869 const G4int maxNumberOfTries = 1000;
870 do {
871 NumberOfTries++;
872 if ( NumberOfTries == 100*(NumberOfTries/100) ) {
873 // After many tries, it is convenient to reduce the values of DcorP, DcorT and
874 // AveragePt2, so that the sampled momenta (respectively, pz, and pt) of the
875 // involved nucleons (or corresponding delta-isomers) are smaller, and therefore
876 // it is more likely to satisfy the momentum conservation.
877 ScaleFactor /= 2.0;
878 DcorP *= ScaleFactor;
879 DcorT *= ScaleFactor;
880 AveragePt2 *= ScaleFactor;
881 }
882 if ( isProjectileNucleus ) {
883 // Sampling of kinematical properties of projectile nucleons
884 isOk = SamplingNucleonKinematics( AveragePt2, maxPtSquare, DcorP,
885 thePrNucleus, PprojResidual,
886 PrResidualMass, ProjectileResidualMassNumber,
887 NumberOfInvolvedNucleonsOfProjectile,
888 TheInvolvedNucleonsOfProjectile, M2proj );
889 }
890 // Sampling of kinematical properties of target nucleons
891 isOk = isOk &&
892 SamplingNucleonKinematics( AveragePt2, maxPtSquare, DcorT,
893 theTargetNucleus, PtargetResidual,
894 TargetResidualMass, TargetResidualMassNumber,
895 NumberOfInvolvedNucleonsOfTarget,
896 TheInvolvedNucleonsOfTarget, M2target );
897
898 if ( M2proj < 0.0 ) {
899 if( M2proj < -0.000001 ) {
900 G4ExceptionDescription ed;
901 ed << "Projectile " << theProjectile.GetDefinition()->GetParticleName()
902 << " Target (Z,A)=(" << theTargetNucleus->GetCharge() << "," << theTargetNucleus->GetMassNumber()
903 << ") M2proj=" << M2proj << " -> sets it to 0.0 !" << G4endl;
904 G4Exception( "G4QGSParticipants::PutOnMassShell(): negative projectile squared mass!",
905 "HAD_QGSPARTICIPANTS_002", JustWarning, ed );
906 }
907 M2proj = 0.0;
908 }
909 sqrtM2proj = std::sqrt( M2proj );
910 if ( M2target < 0.0 ) {
911 G4ExceptionDescription ed;
912 ed << "Projectile " << theProjectile.GetDefinition()->GetParticleName()
913 << " Target (Z,A)=(" << theTargetNucleus->GetCharge() << "," << theTargetNucleus->GetMassNumber()
914 << ") M2target=" << M2target << " -> sets it to 0.0 !" << G4endl;
915 G4Exception( "G4QGSParticipants::PutOnMassShell(): negative target squared mass!",
916 "HAD_QGSPARTICIPANTS_003", JustWarning, ed );
917 M2target = 0.0;
918 };
919 sqrtM2target = std::sqrt( M2target );
920
921 #ifdef debugPutOnMassShell
922 G4cout << "SqrtS, Mp+Mt, Mp, Mt " << SqrtS/GeV << " "
923 << ( sqrtM2proj + sqrtM2target )/GeV << " "
924 << sqrtM2proj/GeV << " " << sqrtM2target/GeV << G4endl;
925 #endif
926
927 if ( ! isOk ) return false;
928 } while ( ( SqrtS < ( sqrtM2proj + sqrtM2target ) ) &&
929 ++NumberOfTries < maxNumberOfTries ); /* Loop checking, 07.08.2015, A.Ribon */
930 if ( NumberOfTries >= maxNumberOfTries ) {
931 return false;
932 }
933 if ( isProjectileNucleus ) {
934 isOk = CheckKinematics( S, SqrtS, M2proj, M2target, YprojectileNucleus, true,
935 NumberOfInvolvedNucleonsOfProjectile,
936 TheInvolvedNucleonsOfProjectile,
937 WminusTarget, WplusProjectile, OuterSuccess );
938 }
939 isOk = isOk &&
940 CheckKinematics( S, SqrtS, M2proj, M2target, YtargetNucleus, false,
941 NumberOfInvolvedNucleonsOfTarget, TheInvolvedNucleonsOfTarget,
942 WminusTarget, WplusProjectile, OuterSuccess );
943 if ( ! isOk ) return false;
944 } while ( ( ! OuterSuccess ) &&
945 ++loopCounter < maxNumberOfLoops ); /* Loop checking, 07.08.2015, A.Ribon */
946 if ( loopCounter >= maxNumberOfLoops ) {
947 return false;
948 }
949
950 // Now the sampling is completed, and we can determine the kinematics of the
951 // whole system. This is done first in the center-of-mass frame, and then it is boosted
952 // to the lab frame. The transverse momentum of the residual nucleus is determined as
953 // the recoil of each hadron (nucleon or delta) which is emitted, i.e. in such a way
954 // to conserve (by construction) the transverse momentum.
955
956 if ( ! isProjectileNucleus ) { // hadron-nucleus collision
957
958 G4double Pzprojectile = WplusProjectile/2.0 - M2projectile/2.0/WplusProjectile;
959 G4double Eprojectile = WplusProjectile/2.0 + M2projectile/2.0/WplusProjectile;
960 Pprojectile.setPz( Pzprojectile );
961 Pprojectile.setE( Eprojectile );
962
963 #ifdef debugPutOnMassShell
964 G4cout << "Proj after in CMS " << Pprojectile/GeV <<" GeV"<< G4endl;
965 #endif
966
967 Pprojectile.transform( toLab );
968 theProjectile.SetMomentum( Pprojectile.vect() );
969 theProjectile.SetTotalEnergy( Pprojectile.e() );
970
971 if ( theProjectileSplitable ) theProjectileSplitable->Set4Momentum(Pprojectile);
972
973 #ifdef debugPutOnMassShell
974 G4cout << "Final proj. mom in Lab. " <<theProjectile.GetMomentum()/GeV<<" "
975 <<theProjectile.GetTotalEnergy()/GeV<<" GeV"<<G4endl;
976 #endif
977
978 } else { // nucleus-nucleus or antinucleus-nucleus collision
979
980 isOk = FinalizeKinematics( WplusProjectile, true, toLab, PrResidualMass,
981 ProjectileResidualMassNumber, NumberOfInvolvedNucleonsOfProjectile,
982 TheInvolvedNucleonsOfProjectile, ProjectileResidual4Momentum );
983
984 #ifdef debugPutOnMassShell
985 G4cout << "Projectile Residual4Momentum in CMS " << ProjectileResidual4Momentum/GeV <<" GeV"<< G4endl;
986 #endif
987
988 if ( ! isOk ) return false;
989
990 ProjectileResidual4Momentum.transform( toLab );
991
992 #ifdef debugPutOnMassShell
993 G4cout << "Projectile Residual4Momentum in Lab " << ProjectileResidual4Momentum/GeV <<" GeV"<< G4endl;
994 #endif
995
996 }
997
998 isOk = FinalizeKinematics( WminusTarget, false, toLab, TargetResidualMass,
999 TargetResidualMassNumber, NumberOfInvolvedNucleonsOfTarget,
1000 TheInvolvedNucleonsOfTarget, TargetResidual4Momentum );
1001
1002 #ifdef debugPutOnMassShell
1003 G4cout << "Target Residual4Momentum in CMS " << TargetResidual4Momentum/GeV << " GeV "<< G4endl;
1004 #endif
1005
1006 if ( ! isOk ) return false;
1007
1008 TargetResidual4Momentum.transform( toLab );
1009
1010 #ifdef debugPutOnMassShell
1011 G4cout << "Target Residual4Momentum in Lab " << TargetResidual4Momentum/GeV << " GeV "<< G4endl;
1012 #endif
1013
1014 return true;
1015
1016 }
1017
1018 //============================================================================
1019
1020 G4ThreeVector G4QGSParticipants::GaussianPt( G4double AveragePt2, G4double maxPtSquare ) const {
1021 // @@ this method is used in FTFModel as well. Should go somewhere common!
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 + G4UniformRand() * ( G4Exp( -x ) -1.0 ) )
1028 : -AveragePt2 * G4Log( 1.0 - G4UniformRand() );
1029 Pt = std::sqrt( Pt2 );
1030 }
1031 G4double phi = G4UniformRand() * twopi;
1032
1033 return G4ThreeVector( Pt*std::cos(phi), Pt*std::sin(phi), 0.0 );
1034 }
1035 //============================================================================
1036
1037 G4bool G4QGSParticipants::
1038 ComputeNucleusProperties( G4V3DNucleus* nucleus, // input parameter
1039 G4LorentzVector& nucleusMomentum, // input & output parameter
1040 G4LorentzVector& residualMomentum, // input & output parameter
1041 G4double& sumMasses, // input & output parameter
1042 G4double& residualExcitationEnergy, // input & output parameter
1043 G4double& residualMass, // input & output parameter
1044 G4int& residualMassNumber, // input & output parameter
1045 G4int& residualCharge ) { // input & output parameter
1046
1047 // This method, which is called only by PutOnMassShell, computes some nucleus properties for:
1048 // - either the target nucleus (which is never an antinucleus): this for any kind
1049 // of hadronic interaction (hadron-nucleus, nucleus-nucleus, antinucleus-nucleus);
1050 // - or the projectile nucleus or antinucleus: this only in the case of nucleus-nucleus
1051 // or antinucleus-nucleus interaction.
1052 // This method assumes that the all the parameters have been initialized by the caller;
1053 // the action of this method consists in modifying all these parameters, except the
1054 // first one. The return value is "false" only in the case the pointer to the nucleus
1055 // is null.
1056
1057 if ( ! nucleus ) return false;
1058
1059 G4double ExcitationEPerWoundedNucleon = GetExcitationEnergyPerWoundedNucleon();
1060
1061 // Loop over the nucleons of the nucleus.
1062 // The nucleons that have been involved in the interaction (either from Glauber or
1063 // Reggeon Cascading) will be candidate to be emitted.
1064 // All the remaining nucleons will be the nucleons of the candidate residual nucleus.
1065 // The variable sumMasses is the amount of energy corresponding to:
1066 // 1. transverse mass of each involved nucleon
1067 // 2. 20.0*MeV separation energy for each involved nucleon
1068 // 3. transverse mass of the residual nucleus
1069 // In this first evaluation of sumMasses, the excitation energy of the residual nucleus
1070 // (residualExcitationEnergy, estimated by adding a constant value to each involved
1071 // nucleon) is not taken into account.
1072 G4Nucleon* aNucleon = 0;
1073 nucleus->StartLoop();
1074 while ( ( aNucleon = nucleus->GetNextNucleon() ) ) { /* Loop checking, 07.08.2015, A.Ribon */
1075 nucleusMomentum += aNucleon->Get4Momentum();
1076 if ( aNucleon->AreYouHit() ) { // Involved nucleons
1077 // Consider in sumMasses the nominal, i.e. on-shell, masses of the nucleons
1078 // (not the current masses, which could be different because the nucleons are off-shell).
1079
1080 sumMasses += std::sqrt( sqr( aNucleon->GetDefinition()->GetPDGMass() )
1081 + aNucleon->Get4Momentum().perp2() );
1082 sumMasses += 20.0*MeV; // Separation energy for a nucleon
1083
1084 //residualExcitationEnergy += ExcitationEPerWoundedNucleon;
1085 residualExcitationEnergy += -ExcitationEPerWoundedNucleon*G4Log( G4UniformRand());
1086 residualMassNumber--;
1087 // The absolute value below is needed only in the case of anti-nucleus.
1088 residualCharge -= std::abs( G4int( aNucleon->GetDefinition()->GetPDGCharge() ) );
1089 } else { // Spectator nucleons
1090 residualMomentum += aNucleon->Get4Momentum();
1091 }
1092 }
1093 #ifdef debugPutOnMassShell
1094 G4cout << "ExcitationEnergyPerWoundedNucleon " << ExcitationEPerWoundedNucleon <<" MeV"<<G4endl
1095 << "\t Residual Charge, MassNumber " << residualCharge << " " << residualMassNumber
1096 << G4endl << "\t Initial Momentum " << nucleusMomentum/GeV<<" GeV"
1097 << G4endl << "\t Residual Momentum " << residualMomentum/GeV<<" GeV"<<G4endl;
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::GetParticleTable()->GetIonTable()->
1107 GetIonMass( residualCharge, residualMassNumber );
1108 if ( residualMassNumber == 1 ) {
1109 residualExcitationEnergy = 0.0;
1110 }
1111 }
1112 sumMasses += std::sqrt( sqr( residualMass ) + residualMomentum.perp2() );
1113 return true;
1114 }
1115
1116
1117 //============================================================================
1118
1119 G4bool G4QGSParticipants::
1120 GenerateDeltaIsobar( const G4double sqrtS, // input parameter
1121 const G4int numberOfInvolvedNucleons, // input parameter
1122 G4Nucleon* involvedNucleons[], // input & output parameter
1123 G4double& sumMasses ) { // input & output parameter
1124
1125 // This method, which is called only by PutOnMassShell, check whether is possible to
1126 // re-interpret some of the involved nucleons as delta-isobars:
1127 // - either by replacing a proton (2212) with a Delta+ (2214),
1128 // - or by replacing a neutron (2112) with a Delta0 (2114).
1129 // The on-shell mass of these delta-isobars is ~1232 MeV, so ~292-294 MeV heavier than
1130 // the corresponding nucleon on-shell mass. However 400.0*MeV is considered to estimate
1131 // the max number of deltas compatible with the available energy.
1132 // The delta-isobars are considered with the same transverse momentum as their
1133 // corresponding nucleons.
1134 // This method assumes that all the parameters have been initialized by the caller;
1135 // the action of this method consists in modifying (eventually) involveNucleons and
1136 // sumMasses. The return value is "false" only in the case that the input parameters
1137 // have unphysical values.
1138
1139 if ( sqrtS < 0.0 || numberOfInvolvedNucleons <= 0 || sumMasses < 0.0 ) return false;
1140
1141 //const G4double ProbDeltaIsobar = 0.05; // Uzhi 6.07.2012
1142 //const G4double ProbDeltaIsobar = 0.25; // Uzhi 13.06.2013
1143 const G4double probDeltaIsobar = 0.10; // A.R. 07.08.2013
1144
1145 G4int maxNumberOfDeltas = G4int( (sqrtS - sumMasses)/(400.0*MeV) );
1146 G4int numberOfDeltas = 0;
1147
1148 for ( G4int i = 0; i < numberOfInvolvedNucleons; i++ ) {
1149 //G4cout << "i maxNumberOfDeltas probDeltaIsobar " << i << " " << maxNumberOfDeltas
1150 // << " " << probDeltaIsobar << G4endl;
1151 if ( G4UniformRand() < probDeltaIsobar && numberOfDeltas < maxNumberOfDeltas ) {
1152 numberOfDeltas++;
1153 if ( ! involvedNucleons[i] ) continue;
1154 G4VSplitableHadron* splitableHadron = involvedNucleons[i]->GetSplitableHadron();
1155 G4double massNuc = std::sqrt( sqr( splitableHadron->GetDefinition()->GetPDGMass() )
1156 + splitableHadron->Get4Momentum().perp2() );
1157 //AR The absolute value below is needed in the case of an antinucleus.
1158 G4int pdgCode = std::abs( splitableHadron->GetDefinition()->GetPDGEncoding() );
1159 const G4ParticleDefinition* old_def = splitableHadron->GetDefinition();
1160 G4int newPdgCode = pdgCode/10; newPdgCode = newPdgCode*10 + 4; // Delta
1161 if ( splitableHadron->GetDefinition()->GetPDGEncoding() < 0 ) newPdgCode *= -1;
1162 const G4ParticleDefinition* ptr =
1163 G4ParticleTable::GetParticleTable()->FindParticle( newPdgCode );
1164 splitableHadron->SetDefinition( ptr );
1165 G4double massDelta = std::sqrt( sqr( splitableHadron->GetDefinition()->GetPDGMass() )
1166 + splitableHadron->Get4Momentum().perp2() );
1167 //G4cout << i << " " << sqrtS/GeV << " " << sumMasses/GeV << " " << massDelta/GeV
1168 // << " " << massNuc << G4endl;
1169 if ( sqrtS < sumMasses + massDelta - massNuc ) { // Change cannot be accepted!
1170 splitableHadron->SetDefinition( old_def );
1171 break;
1172 } else { // Change is accepted
1173 sumMasses += ( massDelta - massNuc );
1174 }
1175 }
1176 }
1177 //G4cout << "maxNumberOfDeltas numberOfDeltas " << maxNumberOfDeltas << " "
1178 // << numberOfDeltas << G4endl;
1179 return true;
1180 }
1181
1182
1183 //============================================================================
1184
1185 G4bool G4QGSParticipants::
1186 SamplingNucleonKinematics( G4double averagePt2, // input parameter
1187 const G4double maxPt2, // input parameter
1188 G4double dCor, // input parameter
1189 G4V3DNucleus* nucleus, // input parameter
1190 const G4LorentzVector& pResidual, // input parameter
1191 const G4double residualMass, // input parameter
1192 const G4int residualMassNumber, // input parameter
1193 const G4int numberOfInvolvedNucleons, // input parameter
1194 G4Nucleon* involvedNucleons[], // input & output parameter
1195 G4double& mass2 ) { // output parameter
1196
1197 // This method, which is called only by PutOnMassShell, does the sampling of:
1198 // - either the target nucleons: this for any kind of hadronic interactions
1199 // (hadron-nucleus, nucleus-nucleus, antinucleus-nucleus);
1200 // - or the projectile nucleons or antinucleons: this only in the case of
1201 // nucleus-nucleus or antinucleus-nucleus interactions, respectively.
1202 // This method assumes that all the parameters have been initialized by the caller;
1203 // the action of this method consists in changing the properties of the nucleons
1204 // whose pointers are in the vector involvedNucleons, as well as changing the
1205 // variable mass2.
1206
1207 if ( ! nucleus ) return false;
1208
1209 if ( residualMassNumber == 0 && numberOfInvolvedNucleons == 1 ) {
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 < numberOfInvolvedNucleons; i++ ) {
1218 G4Nucleon* aNucleon = involvedNucleons[i];
1219 if ( ! aNucleon ) continue;
1220 SumMasses += aNucleon->GetSplitableHadron()->GetDefinition()->GetPDGMass();
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 < numberOfInvolvedNucleons; i++ ) {
1232 G4Nucleon* aNucleon = involvedNucleons[i];
1233 if ( ! aNucleon ) continue;
1234 G4ThreeVector tmpPt = GaussianPt( averagePt2, maxPt2 );
1235 ptSum += tmpPt;
1236 G4ThreeVector tmpX = GaussianPt( dCor*dCor, 1.0 );
1237 G4double x = tmpX.x() +
1238 aNucleon->GetSplitableHadron()->GetDefinition()->GetPDGMass()/SumMasses;
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 of cms) frame but it will not be used.
1245 G4LorentzVector tmp( tmpPt.x(), tmpPt.y(), x, aNucleon->Get4Momentum().e() );
1246 aNucleon->SetMomentum( tmp );
1247 }
1248
1249 if ( xSum < 0.0 || xSum > 1.0 ) success = false;
1250
1251 if ( ! success ) continue;
1252
1253 G4double deltaPx = ( ptSum.x() - pResidual.x() ) / numberOfInvolvedNucleons;
1254 G4double deltaPy = ( ptSum.y() - pResidual.y() ) / numberOfInvolvedNucleons;
1255 G4double delta = 0.0;
1256 if ( residualMassNumber == 0 ) {
1257 delta = ( xSum - 1.0 ) / numberOfInvolvedNucleons;
1258 } else {
1259 delta = 0.0;
1260 }
1261
1262 xSum = 1.0;
1263 mass2 = 0.0;
1264 for ( G4int i = 0; i < numberOfInvolvedNucleons; i++ ) {
1265 G4Nucleon* aNucleon = involvedNucleons[i];
1266 if ( ! aNucleon ) continue;
1267 G4double x = aNucleon->Get4Momentum().pz() - delta;
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 <= 0.0 || xSum > 1.0 ) {
1276 success = false;
1277 break;
1278 }
1279 }
1280 G4double px = aNucleon->Get4Momentum().px() - deltaPx;
1281 G4double py = aNucleon->Get4Momentum().py() - deltaPy;
1282 mass2 += ( sqr( aNucleon->GetSplitableHadron()->GetDefinition()->GetPDGMass() )
1283 + sqr( px ) + sqr( py ) ) / x;
1284 G4LorentzVector tmp( px, py, x, aNucleon->Get4Momentum().e() );
1285 aNucleon->SetMomentum( tmp );
1286 }
1287
1288 if ( success && residualMassNumber != 0 ) {
1289 mass2 += ( sqr( residualMass ) + pResidual.perp2() ) / xSum;
1290 }
1291
1292 #ifdef debugPutOnMassShell
1293 G4cout << "success " << success << G4endl << " Mt " << std::sqrt( mass2 )/GeV << G4endl;
1294 #endif
1295
1296 } while ( ( ! success ) &&
1297 ++loopCounter < maxNumberOfLoops ); /* Loop checking, 07.08.2015, A.Ribon */
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, // input parameter
1310 const G4double sqrtS, // input parameter
1311 const G4double projectileMass2, // input parameter
1312 const G4double targetMass2, // input parameter
1313 const G4double nucleusY, // input parameter
1314 const G4bool isProjectileNucleus, // input parameter
1315 const G4int numberOfInvolvedNucleons, // input parameter
1316 G4Nucleon* involvedNucleons[], // input parameter
1317 G4double& targetWminus, // output parameter
1318 G4double& projectileWplus, // output parameter
1319 G4bool& success ) { // input & output parameter
1320
1321 // This method, which is called only by PutOnMassShell, checks whether the
1322 // kinematics is acceptable or not.
1323 // This method assumes that all the parameters have been initialized by the caller;
1324 // notice that the input boolean parameter isProjectileNucleus is meant to be true
1325 // only in the case of nucleus or antinucleus projectile.
1326 // The action of this method consists in computing targetWminus and projectileWplus
1327 // and setting the parameter success to false in the case that the kinematics should
1328 // be rejeted.
1329
1330 G4double decayMomentum2 = sqr( sValue ) + sqr( projectileMass2 ) + sqr( targetMass2 )
1331 - 2.0*sValue*projectileMass2 - 2.0*sValue*targetMass2
1332 - 2.0*projectileMass2*targetMass2;
1333 targetWminus = ( sValue - projectileMass2 + targetMass2 + std::sqrt( decayMomentum2 ) )
1334 / 2.0 / sqrtS;
1335 projectileWplus = sqrtS - targetMass2/targetWminus;
1336 G4double projectilePz = projectileWplus/2.0 - projectileMass2/2.0/projectileWplus;
1337 G4double projectileE = projectileWplus/2.0 + projectileMass2/2.0/projectileWplus;
1338 G4double projectileY(1.0e5);
1339 if (projectileE - projectilePz > 0.) {
1340 projectileY = 0.5 * G4Log( (projectileE + projectilePz)/
1341 (projectileE - projectilePz) );
1342 }
1343 G4double targetPz = -targetWminus/2.0 + targetMass2/2.0/targetWminus;
1344 G4double targetE = targetWminus/2.0 + targetMass2/2.0/targetWminus;
1345 G4double targetY = 0.5 * G4Log( (targetE + targetPz)/(targetE - targetPz) );
1346
1347 #ifdef debugPutOnMassShell
1348 G4cout << "decayMomentum2 " << decayMomentum2 << G4endl
1349 << "\t targetWminus projectileWplus " << targetWminus << " " << projectileWplus << G4endl
1350 << "\t projectileY targetY " << projectileY << " " << targetY << G4endl;
1351 #endif
1352
1353 for ( G4int i = 0; i < numberOfInvolvedNucleons; i++ ) {
1354 G4Nucleon* aNucleon = involvedNucleons[i];
1355 if ( ! aNucleon ) continue;
1356 G4LorentzVector tmp = aNucleon->Get4Momentum();
1357 G4double mt2 = sqr( tmp.x() ) + sqr( tmp.y() ) +
1358 sqr( aNucleon->GetSplitableHadron()->GetDefinition()->GetPDGMass() );
1359 G4double x = tmp.z();
1360 G4double pz = -targetWminus*x/2.0 + mt2/(2.0*targetWminus*x);
1361 G4double e = targetWminus*x/2.0 + mt2/(2.0*targetWminus*x);
1362 if ( isProjectileNucleus ) {
1363 pz = projectileWplus*x/2.0 - mt2/(2.0*projectileWplus*x);
1364 e = projectileWplus*x/2.0 + mt2/(2.0*projectileWplus*x);
1365 }
1366 G4double nucleonY = 0.5 * G4Log( (e + pz)/(e - pz) );
1367
1368 #ifdef debugPutOnMassShell
1369 G4cout << "i nY pY nY-AY AY " << i << " " << nucleonY << " " << projectileY <<G4endl;
1370 #endif
1371
1372 if ( std::abs( nucleonY - nucleusY ) > 2 ||
1373 ( isProjectileNucleus && targetY > nucleonY ) ||
1374 ( ! isProjectileNucleus && projectileY < nucleonY ) ) {
1375 success = false;
1376 break;
1377 }
1378 }
1379 return true;
1380 }
1381
1382
1383 //============================================================================
1384
1385 G4bool G4QGSParticipants::
1386 FinalizeKinematics( const G4double w, // input parameter
1387 const G4bool isProjectileNucleus, // input parameter
1388 const G4LorentzRotation& boostFromCmsToLab, // input parameter
1389 const G4double residualMass, // input parameter
1390 const G4int residualMassNumber, // input parameter
1391 const G4int numberOfInvolvedNucleons, // input parameter
1392 G4Nucleon* involvedNucleons[], // input & output parameter
1393 G4LorentzVector& residual4Momentum ) { // output parameter
1394
1395 // This method, which is called only by PutOnMassShell, finalizes the kinematics:
1396 // this method is called when we are sure that the sampling of the kinematics is
1397 // acceptable.
1398 // This method assumes that all the parameters have been initialized by the caller;
1399 // notice that the input boolean parameter isProjectileNucleus is meant to be true
1400 // only in the case of nucleus or antinucleus projectile: this information is needed
1401 // because the sign of pz (in the center-of-mass frame) in this case is opposite
1402 // with respect to the case of a normal hadron projectile.
1403 // The action of this method consists in modifying the momenta of the nucleons
1404 // (in the lab frame) and computing the residual 4-momentum (in the center-of-mass
1405 // frame).
1406
1407 G4ThreeVector residual3Momentum( 0.0, 0.0, 1.0 );
1408
1409 for ( G4int i = 0; i < numberOfInvolvedNucleons; i++ ) {
1410 G4Nucleon* aNucleon = involvedNucleons[i];
1411 if ( ! aNucleon ) continue;
1412 G4LorentzVector tmp = aNucleon->Get4Momentum();
1413 residual3Momentum -= tmp.vect();
1414 G4double mt2 = sqr( tmp.x() ) + sqr( tmp.y() ) +
1415 sqr( aNucleon->GetSplitableHadron()->GetDefinition()->GetPDGMass() );
1416 G4double x = tmp.z();
1417 G4double pz = -w * x / 2.0 + mt2 / ( 2.0 * w * x );
1418 G4double e = w * x / 2.0 + mt2 / ( 2.0 * w * x );
1419 // Reverse the sign of pz in the case of nucleus or antinucleus projectile
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 = aNucleon->GetSplitableHadron();
1426 splitableHadron->Set4Momentum( tmp );
1427 #ifdef debugPutOnMassShell
1428 G4cout << "Target involved nucleon No, name, 4Mom "
1429 << i<<" "<<aNucleon->GetDefinition()->GetParticleName()<<" "<<tmp<< G4endl;
1430 #endif
1431 }
1432
1433 G4double residualMt2 = sqr( residualMass ) + sqr( residual3Momentum.x() )
1434 + sqr( residual3Momentum.y() );
1435
1436 #ifdef debugPutOnMassShell
1437 G4cout <<G4endl<< "w residual3Momentum.z() " << w << " " << residual3Momentum.z() << G4endl;
1438 #endif
1439
1440 G4double residualPz = 0.0;
1441 G4double residualE = 0.0;
1442 if ( residualMassNumber != 0 ) {
1443 residualPz = -w * residual3Momentum.z() / 2.0 +
1444 residualMt2 / ( 2.0 * w * residual3Momentum.z() );
1445 residualE = w * residual3Momentum.z() / 2.0 +
1446 residualMt2 / ( 2.0 * w * residual3Momentum.z() );
1447 // Reverse the sign of residualPz in the case of nucleus or antinucleus projectile
1448 if ( isProjectileNucleus ) residualPz *= -1.0;
1449 }
1450
1451 residual4Momentum.setPx( residual3Momentum.x() );
1452 residual4Momentum.setPy( residual3Momentum.y() );
1453 residual4Momentum.setPz( residualPz );
1454 residual4Momentum.setE( residualE );
1455
1456 return true;
1457 }
1458
1459 //======================================================
1460 void G4QGSParticipants::PerformDiffractiveCollisions()
1461 {
1462 #ifdef debugQGSParticipants
1463 G4cout<<G4endl<<"PerformDiffractiveCollisions()......"<<G4endl
1464 <<"theInteractions.size() "<<theInteractions.size()<<G4endl;
1465 #endif
1466
1467 unsigned int i;
1468 for (i = 0; i < theInteractions.size(); i++)
1469 {
1470 G4InteractionContent* anIniteraction = theInteractions[i];
1471 #ifdef debugQGSParticipants
1472 G4cout<<"Interaction # and its status "
1473 <<i<<" "<<theInteractions[i]->GetStatus()<<G4endl;
1474 #endif
1475
1476 G4int InterStatus = theInteractions[i]->GetStatus();
1477 if ( (InterStatus == PrD) || (InterStatus == TrD) || (InterStatus == DD))
1478 { // Selection of diffractive interactions
1479 #ifdef debugQGSParticipants
1480 G4cout<<"Simulation of diffractive interaction. "<<InterStatus <<" PrD/TrD/DD/ND/Qech - 0,1,2,3,4"<<G4endl;
1481 #endif
1482
1483 G4VSplitableHadron* aTarget = anIniteraction->GetTarget();
1484
1485 #ifdef debugQGSParticipants
1486 G4cout<<"The proj. before inter "
1487 <<theProjectileSplitable->Get4Momentum()<<" "
1488 <<theProjectileSplitable->Get4Momentum().mag()<<G4endl;
1489 G4cout<<"The targ. before inter " <<aTarget->Get4Momentum()<<" "
1490 <<aTarget->Get4Momentum().mag()<<G4endl;
1491 #endif
1492
1493 if ( InterStatus == PrD )
1494 theSingleDiffExcitation.ExciteParticipants(theProjectileSplitable, aTarget, TRUE);
1495
1496 if ( InterStatus == TrD )
1497 theSingleDiffExcitation.ExciteParticipants(theProjectileSplitable, aTarget, FALSE);
1498
1499 if ( InterStatus == DD )
1500 theDiffExcitaton.ExciteParticipants(theProjectileSplitable, aTarget);
1501
1502 #ifdef debugQGSParticipants
1503 G4cout<<"The proj. after inter " <<theProjectileSplitable->Get4Momentum()<<" "
1504 <<theProjectileSplitable->Get4Momentum().mag()<<G4endl;
1505 G4cout<<"The targ. after inter " <<aTarget->Get4Momentum()<<" "
1506 <<aTarget->Get4Momentum().mag()<<G4endl;
1507 #endif
1508 }
1509
1510 if ( InterStatus == Qexc )
1511 { // Quark exchange process
1512 #ifdef debugQGSParticipants
1513 G4cout<<"Simulation of interaction with quark exchange."<<G4endl;
1514 #endif
1515 G4VSplitableHadron* aTarget = anIniteraction->GetTarget();
1516
1517 #ifdef debugQGSParticipants
1518 G4cout<<"The proj. before inter " <<theProjectileSplitable->Get4Momentum()<<" "
1519 <<theProjectileSplitable->Get4Momentum().mag()<<G4endl;
1520 G4cout<<"The targ. before inter "<<aTarget->Get4Momentum()<<" "
1521 <<aTarget->Get4Momentum().mag()<<G4endl;
1522 #endif
1523
1524 theQuarkExchange.ExciteParticipants(theProjectileSplitable, aTarget);
1525
1526 #ifdef debugQGSParticipants
1527 G4cout<<"The proj. after inter " <<theProjectileSplitable->Get4Momentum()<<" "
1528 <<theProjectileSplitable->Get4Momentum().mag()<<G4endl;
1529 G4cout<<"The targ. after inter " <<aTarget->Get4Momentum()<<" "
1530 <<aTarget->Get4Momentum().mag()<<G4endl;
1531 #endif
1532 }
1533 }
1534 }
1535
1536 //======================================================
1537 G4bool G4QGSParticipants::DeterminePartonMomenta()
1538 {
1539 if ( ! theProjectileSplitable ) return false;
1540
1541 const G4double aHugeValue = 1.0e+10;
1542
1543 #ifdef debugQGSParticipants
1544 G4cout<<G4endl<<"DeterminePartonMomenta()......"<<G4endl;
1545 G4cout<<"theProjectile status (0 -nondiffr, #0 diffr./reggeon): "<<theProjectileSplitable->GetStatus()<<G4endl;
1546 #endif
1547
1548 if (theProjectileSplitable->GetStatus() != 0) {return false;} // There were only diffractive interactions.
1549
1550 G4LorentzVector Projectile4Momentum = theProjectileSplitable->Get4Momentum();
1551 G4LorentzVector Psum = Projectile4Momentum;
1552
1553 G4double VqM_pr(0.), VaqM_pr(0.), VqM_tr(350.), VqqM_tr(700);
1554 if (std::abs(theProjectile.GetDefinition()->GetBaryonNumber()) != 0) {VqM_pr=350*MeV; VaqM_pr=700*MeV;}
1555
1556 #ifdef debugQGSParticipants
1557 G4cout<<"Projectile 4 momentum "<<Psum<<G4endl
1558 <<"Target nucleon momenta at start"<<G4endl;
1559 G4int NuclNo=0;
1560 #endif
1561
1562 std::vector<G4VSplitableHadron*>::iterator i;
1563
1564 for (i = theTargets.begin(); i != theTargets.end(); i++ )
1565 {
1566 Psum += (*i)->Get4Momentum();
1567 #ifdef debugQGSParticipants
1568 G4cout<<"Nusleus nucleon # and its 4Mom. "<<NuclNo<<" "<<(*i)->Get4Momentum()<<G4endl;
1569 NuclNo++;
1570 #endif
1571 }
1572
1573 G4LorentzRotation toCms( -1*Psum.boostVector() );
1574
1575 G4LorentzVector Ptmp = toCms*Projectile4Momentum;
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---------------"<<G4endl;
1585 G4cout<<"Projectile 4 Mom "<<Projectile4Momentum<<G4endl;
1586 NuclNo=0;
1587 #endif
1588
1589 G4LorentzVector Target4Momentum(0.,0.,0.,0.);
1590 for(i = theTargets.begin(); i != theTargets.end(); i++ )
1591 {
1592 G4LorentzVector tmp= (*i)->Get4Momentum(); tmp.transform( toCms );
1593 (*i)->Set4Momentum( tmp );
1594 #ifdef debugQGSParticipants
1595 G4cout<<"Target nucleon # and 4Mom "<<" "<<NuclNo<<" "<<(*i)->Get4Momentum()<<G4endl;
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 momentum "<<Target4Momentum<<G4endl<<G4endl;
1606 G4cout<<"Target nucleons mom: px, py, z_1, m_i"<<G4endl;
1607 NuclNo=0;
1608 #endif
1609
1610 //G4double PplusProjectile = Projectile4Momentum.plus();
1611 G4double PminusTarget = Target4Momentum.minus();
1612
1613 for(i = theTargets.begin(); i != theTargets.end(); i++ )
1614 {
1615 G4LorentzVector tmp = (*i)->Get4Momentum(); // tmp.boost(bstToCM);
1616
1617 //AR-19Jan2017 : the following line is causing a strange crash when Geant4
1618 // is built in optimized mode.
1619 // To fix it, I get the mass square instead of directly the
1620 // mass from the Lorentz vector, and then I take care of the
1621 // square root. If the mass square is negative, a JustWarning
1622 // exception is thrown, and the mass is set to 0.
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.GetDefinition()->GetParticleName()
1629 << "Â 4-momentum " << Psum << G4endl;
1630 ed << "LorentzVector tmp " << tmp << " with Mass2 " << Mass2 << G4endl;
1631 G4Exception( "G4QGSParticipants::DeterminePartonMomenta(): 4-momentum with negative mass!",
1632 "HAD_QGSPARTICIPANTS_001", JustWarning, ed );
1633 } else {
1634 Mass = std::sqrt( Mass2 );
1635 }
1636
1637 tmp.setPz(tmp.minus()/PminusTarget); tmp.setE(Mass);
1638 (*i)->Set4Momentum(tmp);
1639 #ifdef debugQGSParticipants
1640 G4cout<<"Target nucleons # and mom: "<<NuclNo<<" "<<(*i)->Get4Momentum()<<G4endl;
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 class
1658
1659 const G4ParticleDefinition* theProjectileDefinition = theProjectileSplitable->GetDefinition();
1660 if (theProjectileDefinition == G4PionMinus::PionMinusDefinition()) aBeta = 1.;
1661 if (theProjectileDefinition == G4Gamma::GammaDefinition()) aBeta = 1.;
1662 if (theProjectileDefinition == G4PionPlus::PionPlusDefinition()) aBeta = 1.;
1663 if (theProjectileDefinition == G4PionZero::PionZeroDefinition()) aBeta = 1.;
1664 if (theProjectileDefinition == G4KaonPlus::KaonPlusDefinition()) aBeta = 0.;
1665 if (theProjectileDefinition == G4KaonMinus::KaonMinusDefinition()) aBeta = 0.;
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/100) ) {SigPt/=2.;}
1674
1675 ProjSumMt=0.; ProjSumMt2perX=0.;
1676 TargSumMt=0.; TargSumMt2perX=0.;
1677
1678 Success = true;
1679 G4int nSeaPair = theProjectileSplitable->GetSoftCollisionCount()-1;
1680 #ifdef debugQGSParticipants
1681 G4cout<<"attempt ------------------------ "<<attempt<<G4endl;
1682 G4cout<<"nSeaPair of proj "<<nSeaPair<<G4endl;
1683 #endif
1684
1685 G4double SumPx = 0.;
1686 G4double SumPy = 0.;
1687 G4double SumZ = 0.;
1688 G4int NumberOfUnsampledSeaQuarks = 2*nSeaPair;
1689
1690 G4double Qmass=0.;
1691 for (G4int aSeaPair = 0; aSeaPair < nSeaPair; aSeaPair++)
1692 {
1693 aParton = theProjectileSplitable->GetNextParton(); // for quarks
1694 #ifdef debugQGSParticipants
1695 G4cout<<"Sea quarks: "<<aSeaPair<<" "<<aParton->GetDefinition()->GetParticleName();
1696 #endif
1697 aPtVector = GaussianPt(SigPt, aHugeValue);
1698 tmp.setPx(aPtVector.x()); tmp.setPy(aPtVector.y());
1699 SumPx += aPtVector.x(); SumPy += aPtVector.y();
1700 Mt = std::sqrt(aPtVector.mag2()+sqr(Qmass));
1701 ProjSumMt += Mt;
1702
1703 // Sampling of Z fraction
1704 tmp.setPz(SampleX(Xmin, NumberOfUnsampledSeaQuarks, 2*nSeaPair, aBeta)*(1.0-SumZ));
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->GetNextAntiParton(); // for anti-quarks
1713 #ifdef debugQGSParticipants
1714 G4cout<<" "<<aParton->GetDefinition()->GetParticleName()<<G4endl;
1715 G4cout<<" "<<tmp<<" "<<SumZ<<G4endl;
1716 #endif
1717 aPtVector = GaussianPt(SigPt, aHugeValue);
1718 tmp.setPx(aPtVector.x()); tmp.setPy(aPtVector.y());
1719 SumPx += aPtVector.x(); SumPy += aPtVector.y();
1720 Mt = std::sqrt(aPtVector.mag2()+sqr(Qmass));
1721 ProjSumMt += Mt;
1722
1723 // Sampling of Z fraction
1724 tmp.setPz(SampleX(Xmin, NumberOfUnsampledSeaQuarks, 2*nSeaPair, aBeta)*(1.0-SumZ));
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<<" "<<SumZ<<G4endl;
1733 #endif
1734 }
1735
1736 // For valence quark
1737 aParton = theProjectileSplitable->GetNextParton(); // for quarks
1738 #ifdef debugQGSParticipants
1739 G4cout<<"Val quark of Pr"<<" "<<aParton->GetDefinition()->GetParticleName();
1740 #endif
1741 aPtVector = GaussianPt(SigPt, aHugeValue);
1742 tmp.setPx(aPtVector.x()); tmp.setPy(aPtVector.y());
1743 SumPx += aPtVector.x(); SumPy += aPtVector.y();
1744 Mt = std::sqrt(aPtVector.mag2()+sqr(VqM_pr));
1745 ProjSumMt += Mt;
1746
1747 // Sampling of Z fraction
1748 tmp.setPz(SampleX(Xmin, NumberOfUnsampledSeaQuarks, 2*nSeaPair, aBeta)*(1.0-SumZ));
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->GetNextAntiParton();
1757 #ifdef debugQGSParticipants
1758 G4cout<<" "<<aParton->GetDefinition()->GetParticleName()<<G4endl;
1759 G4cout<<" "<<tmp<<" "<<SumZ<<" (z-fraction)"<<G4endl;
1760 #endif
1761 tmp.setPx(-SumPx); tmp.setPy(-SumPy);
1762 //Uzhi 2019 Mt = std::sqrt(aPtVector.mag2()+sqr(VaqM_pr));
1763 Mt = std::sqrt( sqr(SumPx) + sqr(SumPy) + sqr(VaqM_pr) ); //Uzhi 2019
1764 ProjSumMt += Mt;
1765 tmp.setPz(1.-SumZ);
1766
1767 ProjSumMt2perX +=sqr(Mt)/tmp.pz(); // QQmass=750 MeV
1768 tmp.setE(sqr(Mt));
1769 aParton->Set4Momentum(tmp);
1770 #ifdef debugQGSParticipants
1771 G4cout<<" "<<tmp<<" "<<SumZ+(1.-SumZ)<<" (z-fraction)"<<G4endl;
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 != theTargets.end(); i++ )
1780 {
1781 nSeaPair = (*i)->GetSoftCollisionCount()-1;
1782 #ifdef debugQGSParticipants
1783 G4cout<<"nSeaPair of target N "<<nSeaPair<<G4endl
1784 <<"Target nucleon 4Mom "<<(*i)->Get4Momentum()<<G4endl;
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 < nSeaPair; aSeaPair++)
1795 {
1796 aParton = (*i)->GetNextParton(); // for quarks
1797 #ifdef debugQGSParticipants
1798 G4cout<<"Sea quarks: "<<aSeaPair<<" "<<aParton->GetDefinition()->GetParticleName();
1799 #endif
1800 aPtVector = GaussianPt(SigPt, aHugeValue);
1801 tmp.setPx(aPtVector.x()); tmp.setPy(aPtVector.y());
1802 SumPx += aPtVector.x(); SumPy += aPtVector.y();
1803 Mt=std::sqrt(aPtVector.mag2()+sqr(Qmass));
1804 TargSumMt += Mt;
1805
1806 // Sampling of Z fraction
1807 tmp.setPz(SampleX(Xmin, NumberOfUnsampledSeaQuarks, 2*nSeaPair, aBeta)*(1.0-SumZ));
1808 SumZ += tmp.z();
1809 tmp.setPz((*i)->Get4Momentum().pz()*tmp.pz());
1810 NumberOfUnsampledSeaQuarks--;
1811 TargSumMt2perX +=sqr(Mt)/tmp.pz();
1812 tmp.setE(sqr(Mt));
1813 aParton->Set4Momentum(tmp);
1814
1815 aParton = (*i)->GetNextAntiParton(); // for anti-quarks
1816 #ifdef debugQGSParticipants
1817 G4cout<<" "<<aParton->GetDefinition()->GetParticleName()<<G4endl;
1818 G4cout<<" "<<tmp<<" "<<SumZ<<G4endl;
1819 #endif
1820 aPtVector = GaussianPt(SigPt, aHugeValue);
1821 tmp.setPx(aPtVector.x()); tmp.setPy(aPtVector.y());
1822 SumPx += aPtVector.x(); SumPy += aPtVector.y();
1823 Mt=std::sqrt(aPtVector.mag2()+sqr(Qmass));
1824 TargSumMt += Mt;
1825
1826 // Sampling of Z fraction
1827 tmp.setPz(SampleX(Xmin, NumberOfUnsampledSeaQuarks, 2*nSeaPair, aBeta)*(1.0-SumZ));
1828 SumZ += tmp.z();
1829 tmp.setPz((*i)->Get4Momentum().pz()*tmp.pz());
1830 NumberOfUnsampledSeaQuarks--;
1831 TargSumMt2perX +=sqr(Mt)/tmp.pz();
1832 tmp.setE(sqr(Mt));
1833 aParton->Set4Momentum(tmp);
1834 #ifdef debugQGSParticipants
1835 G4cout<<" "<<tmp<<" "<<" "<<SumZ<<G4endl;
1836 #endif
1837 }
1838
1839 // Valence quark
1840 aParton = (*i)->GetNextParton(); // for quarks
1841 #ifdef debugQGSParticipants
1842 G4cout<<"Val quark of Tr"<<" "<<aParton->GetDefinition()->GetParticleName();
1843 #endif
1844 aPtVector = GaussianPt(SigPt, aHugeValue);
1845 tmp.setPx(aPtVector.x()); tmp.setPy(aPtVector.y());
1846 SumPx += aPtVector.x(); SumPy += aPtVector.y();
1847 Mt=std::sqrt(aPtVector.mag2()+sqr(VqM_tr));
1848 TargSumMt += Mt;
1849
1850 // Sampling of Z fraction
1851 tmp.setPz(SampleX(Xmin, NumberOfUnsampledSeaQuarks, 2*nSeaPair, aBeta)*(1.0-SumZ));
1852 SumZ += tmp.z();
1853 tmp.setPz((*i)->Get4Momentum().pz()*tmp.pz());
1854 TargSumMt2perX +=sqr(Mt)/tmp.pz();
1855 tmp.setE(sqr(Mt));
1856 aParton->Set4Momentum(tmp);
1857
1858 // Valence di-quark
1859 aParton = (*i)->GetNextAntiParton(); // for quarks
1860 #ifdef debugQGSParticipants
1861 G4cout<<" "<<aParton->GetDefinition()->GetParticleName()<<G4endl;
1862 G4cout<<" "<<tmp<<" "<<SumZ<<" (total z-sum) "<<G4endl;
1863 #endif
1864 tmp.setPx(-SumPx); tmp.setPy(-SumPy);
1865 //Uzhi 2019 Mt=std::sqrt(aPtVector.mag2()+sqr(VqqM_tr));
1866 Mt=std::sqrt( sqr(SumPx) + sqr(SumPy) + sqr(VqqM_tr) ); //Uzhi 2019
1867 TargSumMt += Mt;
1868
1869 tmp.setPz((*i)->Get4Momentum().pz()*(1.0 - SumZ));
1870 TargSumMt2perX +=sqr(Mt)/tmp.pz();
1871 tmp.setE(sqr(Mt));
1872 aParton->Set4Momentum(tmp);
1873 #ifdef debugQGSParticipants
1874 G4cout<<" "<<tmp<<" "<<1.0<<" "<<(*i)->Get4Momentum().pz()<<G4endl;
1875 #endif
1876
1877 } // End of for(i = theTargets.begin(); i != theTargets.end(); i++ )
1878
1879 if( ProjSumMt + TargSumMt > SqrtS ) {
1880 Success = false; continue;}
1881 if( std::sqrt(ProjSumMt2perX) + std::sqrt(TargSumMt2perX) > SqrtS ) {
1882 Success = false; continue;}
1883
1884 } while( (!Success) &&
1885 attempt < maxNumberOfAttempts ); /* Loop checking, 07.08.2015, A.Ribon */
1886
1887 if ( attempt >= maxNumberOfAttempts ) {
1888 return false;
1889 }
1890
1891 //+++++++++++++++++++++++++++++++++++++++++++
1892
1893 G4double DecayMomentum2 = sqr(S) + sqr(ProjSumMt2perX) + sqr(TargSumMt2perX)
1894 - 2.0*S*ProjSumMt2perX - 2.0*S*TargSumMt2perX - 2.0*ProjSumMt2perX*TargSumMt2perX;
1895
1896 G4double targetWminus=( S - ProjSumMt2perX + TargSumMt2perX + std::sqrt( DecayMomentum2 ))/2.0/SqrtS;
1897 G4double projectileWplus = SqrtS - TargSumMt2perX/targetWminus;
1898
1899 G4LorentzVector Tmp(0.,0.,0.,0.);
1900 G4double z(0.);
1901
1902 G4int nSeaPair = theProjectileSplitable->GetSoftCollisionCount()-1;
1903 #ifdef debugQGSParticipants
1904 G4cout<<"Backward transformation ===================="<<G4endl;
1905 G4cout<<"nSeaPair of proj "<<nSeaPair<<G4endl;
1906 #endif
1907
1908 for (G4int aSeaPair = 0; aSeaPair < nSeaPair; aSeaPair++)
1909 {
1910 aParton = theProjectileSplitable->GetNextParton(); // for quarks
1911 #ifdef debugQGSParticipants
1912 G4cout<<"Sea quarks: "<<aSeaPair<<" "<<aParton->GetDefinition()->GetParticleName();
1913 #endif
1914 Tmp =aParton->Get4Momentum(); z=Tmp.z();
1915
1916 Tmp.setPz(projectileWplus*z/2.0 - Tmp.e()/(2.0*z*projectileWplus));
1917 Tmp.setE( projectileWplus*z/2.0 + Tmp.e()/(2.0*z*projectileWplus));
1918 Tmp.transform( toLab );
1919
1920 aParton->Set4Momentum(Tmp);
1921
1922 aParton = theProjectileSplitable->GetNextAntiParton(); // for anti-quarks
1923 #ifdef debugQGSParticipants
1924 G4cout<<" "<<aParton->GetDefinition()->GetParticleName()<<G4endl;
1925 G4cout<<" "<<Tmp<<" "<<Tmp.mag()<<G4endl;
1926 #endif
1927 Tmp =aParton->Get4Momentum(); z=Tmp.z();
1928 Tmp.setPz(projectileWplus*z/2.0 - Tmp.e()/(2.0*z*projectileWplus));
1929 Tmp.setE( projectileWplus*z/2.0 + Tmp.e()/(2.0*z*projectileWplus));
1930 Tmp.transform( toLab );
1931
1932 aParton->Set4Momentum(Tmp);
1933 #ifdef debugQGSParticipants
1934 G4cout<<" "<<Tmp<<" "<<Tmp.mag()<<G4endl;
1935 #endif
1936 }
1937
1938 // For valence quark
1939 aParton = theProjectileSplitable->GetNextParton(); // for quarks
1940 #ifdef debugQGSParticipants
1941 G4cout<<"Val quark of Pr"<<" "<<aParton->GetDefinition()->GetParticleName();
1942 #endif
1943 Tmp =aParton->Get4Momentum(); z=Tmp.z();
1944 Tmp.setPz(projectileWplus*z/2.0 - Tmp.e()/(2.0*z*projectileWplus));
1945 Tmp.setE( projectileWplus*z/2.0 + Tmp.e()/(2.0*z*projectileWplus));
1946 Tmp.transform( toLab );
1947
1948 aParton->Set4Momentum(Tmp);
1949
1950 // For valence di-quark
1951 aParton = theProjectileSplitable->GetNextAntiParton();
1952 #ifdef debugQGSParticipants
1953 G4cout<<" "<<aParton->GetDefinition()->GetParticleName()<<G4endl;
1954 G4cout<<" "<<Tmp<<" "<<Tmp.mag()<<" (mass)"<<G4endl;
1955 #endif
1956 Tmp =aParton->Get4Momentum(); z=Tmp.z();
1957 Tmp.setPz(projectileWplus*z/2.0 - Tmp.e()/(2.0*z*projectileWplus));
1958 Tmp.setE( projectileWplus*z/2.0 + Tmp.e()/(2.0*z*projectileWplus));
1959 Tmp.transform( toLab );
1960
1961 aParton->Set4Momentum(Tmp);
1962
1963 #ifdef debugQGSParticipants
1964 G4cout<<" "<<Tmp<<" "<<Tmp.mag()<<" (mass)"<<G4endl;
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.end(); i++ )
1972 {
1973 nSeaPair = (*i)->GetSoftCollisionCount()-1;
1974 #ifdef debugQGSParticipants
1975 G4cout<<"nSeaPair of target and N# "<<nSeaPair<<" "<<NuclNo<<G4endl;
1976 NuclNo++;
1977 #endif
1978 for (G4int aSeaPair = 0; aSeaPair < nSeaPair; aSeaPair++)
1979 {
1980 aParton = (*i)->GetNextParton(); // for quarks
1981 #ifdef debugQGSParticipants
1982 G4cout<<"Sea quarks: "<<aSeaPair<<" "<<aParton->GetDefinition()->GetParticleName();
1983 #endif
1984 Tmp =aParton->Get4Momentum(); z=Tmp.z();
1985 Tmp.setPz(-targetWminus*z/2.0 + Tmp.e()/(2.0*z*targetWminus));
1986 Tmp.setE( targetWminus*z/2.0 + Tmp.e()/(2.0*z*targetWminus));
1987 Tmp.transform( toLab );
1988
1989 aParton->Set4Momentum(Tmp);
1990
1991 aParton = (*i)->GetNextAntiParton(); // for quarks
1992 #ifdef debugQGSParticipants
1993 G4cout<<" "<<aParton->GetDefinition()->GetParticleName()<<G4endl;
1994 G4cout<<" "<<Tmp<<" "<<Tmp.mag()<<G4endl;
1995 #endif
1996 Tmp =aParton->Get4Momentum(); z=Tmp.z();
1997 Tmp.setPz(-targetWminus*z/2.0 + Tmp.e()/(2.0*z*targetWminus));
1998 Tmp.setE( targetWminus*z/2.0 + Tmp.e()/(2.0*z*targetWminus));
1999 Tmp.transform( toLab );
2000
2001 aParton->Set4Momentum(Tmp);
2002 #ifdef debugQGSParticipants
2003 G4cout<<" "<<Tmp<<" "<<Tmp.mag()<<G4endl;
2004 #endif
2005 }
2006
2007 // Valence quark
2008
2009 aParton = (*i)->GetNextParton(); // for quarks
2010 #ifdef debugQGSParticipants
2011 G4cout<<"Val quark of Tr"<<" "<<aParton->GetDefinition()->GetParticleName();
2012 #endif
2013 Tmp =aParton->Get4Momentum(); z=Tmp.z();
2014 Tmp.setPz(-targetWminus*z/2.0 + Tmp.e()/(2.0*z*targetWminus));
2015 Tmp.setE( targetWminus*z/2.0 + Tmp.e()/(2.0*z*targetWminus));
2016 Tmp.transform( toLab );
2017
2018 aParton->Set4Momentum(Tmp);
2019
2020 // Valence di-quark
2021 aParton = (*i)->GetNextAntiParton(); // for quarks
2022 #ifdef debugQGSParticipants
2023 G4cout<<" "<<aParton->GetDefinition()->GetParticleName()<<G4endl;
2024 G4cout<<" "<<Tmp<<" "<<Tmp.mag()<<" (mass)"<<G4endl;
2025 #endif
2026 Tmp =aParton->Get4Momentum(); z=Tmp.z();
2027 Tmp.setPz(-targetWminus*z/2.0 + Tmp.e()/(2.0*z*targetWminus));
2028 Tmp.setE( targetWminus*z/2.0 + Tmp.e()/(2.0*z*targetWminus));
2029 Tmp.transform( toLab );
2030
2031 aParton->Set4Momentum(Tmp);
2032 #ifdef debugQGSParticipants
2033 G4cout<<" "<<Tmp<<" "<<Tmp.mag()<<" (mass)"<<G4endl;
2034 NuclNo++;
2035 #endif
2036 } // End of for(i = theTargets.begin(); i != theTargets.end(); i++ )
2037
2038 return true;
2039 }
2040
2041 //======================================================
2042 G4double G4QGSParticipants::
2043 SampleX(G4double, G4int nSea, G4int, G4double aBeta)
2044 {
2045 G4double Oalfa = 1./(alpha + 1.);
2046 G4double Obeta = 1./(aBeta + (alpha + 1.)*nSea + 1.); // ?
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::GetInstance()->powA(Ksi1,Oalfa);
2054 Ksi2 = G4UniformRand(); r2 = G4Pow::GetInstance()->powA(Ksi2,Obeta);
2055 r12=r1+r2;
2056 } while( ( r12 > 1.) &&
2057 ++loopCounter < maxNumberOfLoops ); /* Loop checking, 07.08.2015, A.Ribon */
2058 if ( loopCounter >= maxNumberOfLoops ) {
2059 return 0.5; // Just an acceptable value, without any physics consideration.
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() ..................."<<G4endl;
2072 #endif
2073
2074 if ( ! theProjectileSplitable ) {
2075 #ifdef debugQGSParticipants
2076 G4cout<<"BAD situation: theProjectileSplitable is NULL ! Returning immediately"<<G4endl;
2077 #endif
2078 return;
2079 }
2080
2081 #ifdef debugQGSParticipants
2082 G4cout<<"theProjectileSplitable->GetStatus() "<<theProjectileSplitable->GetStatus()<<G4endl;
2083 G4LorentzVector str4Mom;
2084 #endif
2085
2086 if( theProjectileSplitable->GetStatus() == 1 ) // The projectile has participated only in diffr. inter.
2087 {
2088 G4ThreeVector Position = theProjectileSplitable->GetPosition();
2089
2090 G4PartonPair * aPair = new G4PartonPair(theProjectileSplitable->GetNextParton(),
2091 theProjectileSplitable->GetNextAntiParton(),
2092 G4PartonPair::DIFFRACTIVE, G4PartonPair::TARGET);
2093 #ifdef debugQGSParticipants
2094 G4cout << "Pr. Diffr. String: Qs 4mom X " <<G4endl;
2095 G4cout << " " << aPair->GetParton1()->GetPDGcode() << " "
2096 << aPair->GetParton1()->Get4Momentum() << " "
2097 << aPair->GetParton1()->GetX() << " " << G4endl;
2098 G4cout << " " << aPair->GetParton2()->GetPDGcode() << " "
2099 << aPair->GetParton2()->Get4Momentum() << " "
2100 << aPair->GetParton2()->GetX() << " " << G4endl;
2101 str4Mom += aPair->GetParton1()->Get4Momentum();
2102 str4Mom += aPair->GetParton2()->Get4Momentum();
2103 #endif
2104
2105 thePartonPairs.push_back(aPair);
2106 }
2107
2108 G4int N_EnvTarg = NumberOfInvolvedNucleonsOfTarget;
2109
2110 for ( G4int i = 0; i < N_EnvTarg; i++ ) {
2111 G4Nucleon* aTargetNucleon = TheInvolvedNucleonsOfTarget[ i ];
2112
2113 G4VSplitableHadron* HitTargetNucleon = aTargetNucleon->GetSplitableHadron();
2114
2115 #ifdef debugQGSParticipants
2116 G4cout<<"Involved Nucleon # and its status "<<i<<" "<<HitTargetNucleon->GetStatus()<<G4endl;
2117 #endif
2118 if( HitTargetNucleon->GetStatus() >= 1) // Create diffractive string
2119 {
2120 G4ThreeVector Position = HitTargetNucleon->GetPosition();
2121
2122 G4PartonPair * aPair = new G4PartonPair(HitTargetNucleon->GetNextParton(),
2123 HitTargetNucleon->GetNextAntiParton(),
2124 G4PartonPair::DIFFRACTIVE, G4PartonPair::TARGET);
2125 #ifdef debugQGSParticipants
2126 G4cout << "Tr. Diffr. String: Qs 4mom X " <<G4endl;
2127 G4cout << "Diffr. String " << aPair->GetParton1()->GetPDGcode() << " "
2128 << aPair->GetParton1()->Get4Momentum() << " "
2129 << aPair->GetParton1()->GetX() << " " << G4endl;
2130 G4cout << " " << aPair->GetParton2()->GetPDGcode() << " "
2131 << aPair->GetParton2()->Get4Momentum() << " "
2132 << aPair->GetParton2()->GetX() << " " << G4endl;
2133
2134 str4Mom += aPair->GetParton1()->Get4Momentum();
2135 str4Mom += aPair->GetParton2()->Get4Momentum();
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 interactions"<<G4endl;
2146 #endif
2147 std::vector<G4InteractionContent*>::iterator i;
2148 i = theInteractions.begin();
2149 while ( i != theInteractions.end() ) /* Loop checking, 07.08.2015, A.Ribon */
2150 {
2151 G4InteractionContent* anIniteraction = *i;
2152 G4PartonPair * aPair = NULL;
2153
2154 #ifdef debugQGSParticipants
2155 G4cout<<"An interaction # and soft coll. # "<<IntNo<<" "
2156 <<anIniteraction->GetNumberOfSoftCollisions()<<G4endl;
2157 IntNo++;
2158 #endif
2159 if (anIniteraction->GetNumberOfSoftCollisions())
2160 {
2161 G4VSplitableHadron* pProjectile = anIniteraction->GetProjectile();
2162 G4VSplitableHadron* pTarget = anIniteraction->GetTarget();
2163
2164 for (G4int j = 0; j < anIniteraction->GetNumberOfSoftCollisions(); j++)
2165 {
2166 aPair = new G4PartonPair(pTarget->GetNextParton(), pProjectile->GetNextAntiParton(),
2167 G4PartonPair::SOFT, G4PartonPair::TARGET);
2168 #ifdef debugQGSParticipants
2169 G4cout << "SoftPair " << aPair->GetParton1()->GetPDGcode() << " "
2170 << aPair->GetParton1()->Get4Momentum() << " "
2171 << aPair->GetParton1()->Get4Momentum().mag()<<G4endl;
2172 G4cout << " " << aPair->GetParton2()->GetPDGcode() << " "
2173 << aPair->GetParton2()->Get4Momentum() << " "
2174 <<aPair->GetParton2()->Get4Momentum().mag()<<G4endl;
2175 str4Mom += aPair->GetParton1()->Get4Momentum();
2176 str4Mom += aPair->GetParton2()->Get4Momentum();
2177 #endif
2178
2179 thePartonPairs.push_back(aPair);
2180
2181 aPair = new G4PartonPair(pProjectile->GetNextParton(), pTarget->GetNextAntiParton(),
2182 G4PartonPair::SOFT, G4PartonPair::PROJECTILE);
2183 #ifdef debugQGSParticipants
2184 G4cout << "SoftPair " << aPair->GetParton1()->GetPDGcode() << " "
2185 << aPair->GetParton1()->Get4Momentum() << " "
2186 << aPair->GetParton1()->Get4Momentum().mag()<<G4endl;
2187 G4cout << " " << aPair->GetParton2()->GetPDGcode() << " "
2188 << aPair->GetParton2()->Get4Momentum() << " "
2189 << aPair->GetParton2()->Get4Momentum().mag()<<G4endl;
2190 #endif
2191 #ifdef debugQGSParticipants
2192 str4Mom += aPair->GetParton1()->Get4Momentum();
2193 str4Mom += aPair->GetParton2()->Get4Momentum();
2194 #endif
2195
2196 thePartonPairs.push_back(aPair);
2197 }
2198
2199 delete *i;
2200 i=theInteractions.erase(i); // i now points to the next interaction
2201 }
2202 else
2203 {
2204 i++;
2205 }
2206 } // End of while ( i != theInteractions.end() )
2207 #ifdef debugQGSParticipants
2208 G4cout << "Sum of strings 4 momenta " << str4Mom << G4endl<<G4endl;
2209 #endif
2210 }
2211
2212 //============================================================================
2213
2214 void G4QGSParticipants::GetResiduals() {
2215 // This method is needed for the correct application of G4PrecompoundModelInterface
2216
2217 #ifdef debugQGSParticipants
2218 G4cout << "GetResiduals(): GetProjectileNucleus()? " << GetProjectileNucleus() << G4endl;
2219 #endif
2220
2221 #ifdef debugQGSParticipants
2222 G4cout << "NumberOfInvolvedNucleonsOfTarget "<< NumberOfInvolvedNucleonsOfTarget << G4endl;
2223 #endif
2224
2225 G4double DeltaExcitationE = TargetResidualExcitationEnergy / G4double( NumberOfInvolvedNucleonsOfTarget );
2226 G4LorentzVector DeltaPResidualNucleus = TargetResidual4Momentum /
2227 G4double( NumberOfInvolvedNucleonsOfTarget );
2228
2229 for ( G4int i = 0; i < NumberOfInvolvedNucleonsOfTarget; i++ ) {
2230 G4Nucleon* aNucleon = TheInvolvedNucleonsOfTarget[i];
2231
2232 #ifdef debugQGSParticipants
2233 G4VSplitableHadron* targetSplitable = aNucleon->GetSplitableHadron();
2234 G4cout << i << " Hit? " << aNucleon->AreYouHit() << " " << targetSplitable << G4endl;
2235 if ( targetSplitable ) G4cout << i << "Status " << targetSplitable->GetStatus() << G4endl;
2236 #endif
2237
2238 G4LorentzVector tmp = -DeltaPResidualNucleus;
2239 aNucleon->SetMomentum( tmp );
2240 aNucleon->SetBindingEnergy( DeltaExcitationE );
2241 }
2242
2243 //-------------------------------------
2244 if( TargetResidualMassNumber != 0 )
2245 {
2246 G4ThreeVector bstToCM =TargetResidual4Momentum.findBoostToCM();
2247
2248 G4V3DNucleus* theTargetNucleus = GetTargetNucleus();
2249 G4LorentzVector residualMomentum(0.,0.,0.,0.);
2250 G4Nucleon* aNucleon = 0;
2251 theTargetNucleus->StartLoop();
2252 while ( ( aNucleon = theTargetNucleus->GetNextNucleon() ) ) { /* Loop checking, 07.08.2015, A.Ribon */
2253 if ( !aNucleon->AreYouHit() ) {
2254 G4LorentzVector tmp=aNucleon->Get4Momentum(); tmp.boost(bstToCM);
2255 aNucleon->SetMomentum(tmp);
2256 residualMomentum +=tmp;
2257 }
2258 }
2259
2260 residualMomentum/=TargetResidualMassNumber;
2261
2262 G4double Mass = TargetResidual4Momentum.mag();
2263 G4double SumMasses=0.;
2264
2265 aNucleon = 0;
2266 theTargetNucleus->StartLoop();
2267 while ( ( aNucleon = theTargetNucleus->GetNextNucleon() ) ) { /* Loop checking, 07.08.2015, A.Ribon */
2268 if ( !aNucleon->AreYouHit() ) {
2269 G4LorentzVector tmp=aNucleon->Get4Momentum() - residualMomentum;
2270 G4double E=std::sqrt(tmp.vect().mag2()+
2271 sqr(aNucleon->GetDefinition()->GetPDGMass()-aNucleon->GetBindingEnergy()));
2272 tmp.setE(E); aNucleon->SetMomentum(tmp);
2273 SumMasses+=E;
2274 }
2275 }
2276
2277 G4double Chigh=Mass/SumMasses; G4double Clow=0; 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->GetNextNucleon() ) ) { /* Loop checking, 07.08.2015, A.Ribon */
2288 if ( !aNucleon->AreYouHit() ) {
2289 G4LorentzVector tmp=aNucleon->Get4Momentum();
2290 G4double E=std::sqrt(tmp.vect().mag2()*sqr(C)+
2291 sqr(aNucleon->GetDefinition()->GetPDGMass()-aNucleon->GetBindingEnergy()));
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 ); /* Loop checking, 07.08.2015, A.Ribon */
2301 if ( loopCounter >= maxNumberOfLoops ) {
2302 #ifdef debugQGSParticipants
2303 G4cout <<"BAD situation: forced loop exit!" << G4endl;
2304 #endif
2305 // Perhaps there is something to set here...
2306 } else {
2307 aNucleon = 0;
2308 theTargetNucleus->StartLoop();
2309 while ( ( aNucleon = theTargetNucleus->GetNextNucleon() ) ) { /* Loop checking, 07.08.2015, A.Ribon */
2310 if ( !aNucleon->AreYouHit() ) {
2311 G4LorentzVector tmp=aNucleon->Get4Momentum()*C;
2312 G4double E=std::sqrt(tmp.vect().mag2()+
2313 sqr(aNucleon->GetDefinition()->GetPDGMass()-aNucleon->GetBindingEnergy()));
2314 tmp.setE(E); tmp.boost(-bstToCM);
2315 aNucleon->SetMomentum(tmp);
2316 }
2317 }
2318 }
2319
2320 } // End of if( TargetResidualMassNumber != 0 )
2321 //-------------------------------------
2322
2323 #ifdef debugQGSParticipants
2324 G4cout << "End GetResiduals -----------------" << G4endl;
2325 #endif
2326
2327 }
2328
2329
2330 //======================================================
2331 void G4QGSParticipants::PerformSoftCollisions() // It is not used
2332 {
2333 std::vector<G4InteractionContent*>::iterator i;
2334 G4LorentzVector str4Mom;
2335 i = theInteractions.begin();
2336 while ( i != theInteractions.end() ) /* Loop checking, 07.08.2015, A.Ribon */
2337 {
2338 G4InteractionContent* anIniteraction = *i;
2339 G4PartonPair * aPair = NULL;
2340 if (anIniteraction->GetNumberOfSoftCollisions())
2341 {
2342 G4VSplitableHadron* pProjectile = anIniteraction->GetProjectile();
2343 G4VSplitableHadron* pTarget = anIniteraction->GetTarget();
2344 for (G4int j = 0; j < anIniteraction->GetNumberOfSoftCollisions(); j++)
2345 {
2346 aPair = new G4PartonPair(pTarget->GetNextParton(), pProjectile->GetNextAntiParton(),
2347 G4PartonPair::SOFT, G4PartonPair::TARGET);
2348 #ifdef debugQGSParticipants
2349 G4cout << "SoftPair " << aPair->GetParton1()->GetPDGcode() << " "
2350 << aPair->GetParton1()->Get4Momentum() << " "
2351 << aPair->GetParton1()->GetX() << " " << G4endl;
2352 G4cout << " " << aPair->GetParton2()->GetPDGcode() << " "
2353 << aPair->GetParton2()->Get4Momentum() << " "
2354 << aPair->GetParton2()->GetX() << " " << G4endl;
2355 #endif
2356 #ifdef debugQGSParticipants
2357 str4Mom += aPair->GetParton1()->Get4Momentum();
2358 str4Mom += aPair->GetParton2()->Get4Momentum();
2359 #endif
2360 thePartonPairs.push_back(aPair);
2361 aPair = new G4PartonPair(pProjectile->GetNextParton(), pTarget->GetNextAntiParton(),
2362 G4PartonPair::SOFT, G4PartonPair::PROJECTILE);
2363 #ifdef debugQGSParticipants
2364 G4cout << "SoftPair " << aPair->GetParton1()->GetPDGcode() << " "
2365 << aPair->GetParton1()->Get4Momentum() << " "
2366 << aPair->GetParton1()->GetX() << " " << G4endl;
2367 G4cout << " " << aPair->GetParton2()->GetPDGcode() << " "
2368 << aPair->GetParton2()->Get4Momentum() << " "
2369 << aPair->GetParton2()->GetX() << " " << G4endl;
2370 #endif
2371 #ifdef debugQGSParticipants
2372 str4Mom += aPair->GetParton1()->Get4Momentum();
2373 str4Mom += aPair->GetParton2()->Get4Momentum();
2374 #endif
2375 thePartonPairs.push_back(aPair);
2376 }
2377 delete *i;
2378 i=theInteractions.erase(i); // i now points to the next interaction
2379 } else {
2380 i++;
2381 }
2382 }
2383 #ifdef debugQGSParticipants
2384 G4cout << " string 4 mom " << str4Mom << G4endl;
2385 #endif
2386 }
2387
2388
2389 //************************************************
2390 G4VSplitableHadron* G4QGSParticipants::SelectInteractions(const G4ReactionProduct &thePrimary)
2391 {
2392 // Check reaction threshold - goes to CheckThreshold
2393
2394 theProjectileSplitable = new G4QGSMSplitableHadron(thePrimary, TRUE);
2395 theProjectileSplitable->SetStatus(1);
2396
2397 G4LorentzVector aPrimaryMomentum(thePrimary.GetMomentum(), thePrimary.GetTotalEnergy());
2398 G4LorentzVector aTargetNMomentum(0.,0.,0.,938.);
2399
2400 if ((!(aPrimaryMomentum.e()>-1)) && (!(aPrimaryMomentum.e()<1)) )
2401 {
2402 throw G4HadronicException(__FILE__, __LINE__,
2403 "G4GammaParticipants::SelectInteractions: primary nan energy.");
2404 }
2405 G4double S = (aPrimaryMomentum + aTargetNMomentum).mag2();
2406 G4double ThresholdMass = thePrimary.GetMass() + 938.;
2407 ModelMode = SOFT;
2408
2409 #ifdef debugQGSParticipants
2410 G4cout << "Gamma projectile - SelectInteractions " << G4endl;
2411 G4cout<<"Energy and Nucleus "<<thePrimary.GetTotalEnergy()<<" "<<theNucleus->GetMassNumber()<<G4endl;
2412 G4cout << "SqrtS ThresholdMass ModelMode " <<std::sqrt(S)<<" "<<ThresholdMass<<" "<<ModelMode<< G4endl;
2413 #endif
2414
2415 if (sqr(ThresholdMass + ThresholdParameter) > S)
2416 {
2417 ModelMode = DIFFRACTIVE;
2418 }
2419 if (sqr(ThresholdMass + QGSMThreshold) > S) // thus only diffractive in cascade!
2420 {
2421 ModelMode = DIFFRACTIVE;
2422 }
2423
2424 // first find the collisions HPW
2425 std::for_each(theInteractions.begin(), theInteractions.end(), DeleteInteractionContent());
2426 theInteractions.clear();
2427
2428 G4int theCurrent = G4int(theNucleus->GetMassNumber()*G4UniformRand());
2429 G4int NucleonNo=0;
2430
2431 theNucleus->StartLoop();
2432 G4Nucleon * pNucleon = 0;
2433
2434 while( (pNucleon = theNucleus->GetNextNucleon()) ) /* Loop checking, 07.08.2015, A.Ribon */
2435 { if(NucleonNo == theCurrent) break; NucleonNo++; }
2436
2437 if ( pNucleon ) {
2438 G4QGSMSplitableHadron* aTarget = new G4QGSMSplitableHadron(*pNucleon);
2439
2440 pNucleon->Hit(aTarget);
2441
2442 if ( (0.06 > G4UniformRand() &&(ModelMode==SOFT)) || (ModelMode==DIFFRACTIVE ) )
2443 {
2444 G4InteractionContent * aInteraction = new G4InteractionContent(theProjectileSplitable);
2445 theProjectileSplitable->SetStatus(1*theProjectileSplitable->GetStatus());
2446
2447 aInteraction->SetTarget(aTarget);
2448 aInteraction->SetTargetNucleon(pNucleon);
2449 aTarget->SetCollisionCount(0);
2450 aTarget->SetStatus(1);
2451
2452 aInteraction->SetNumberOfDiffractiveCollisions(1);
2453 aInteraction->SetNumberOfSoftCollisions(0);
2454 aInteraction->SetStatus(1);
2455
2456 theInteractions.push_back(aInteraction);
2457 }
2458 else
2459 {
2460 // nondiffractive soft interaction occurs
2461 aTarget->IncrementCollisionCount(1);
2462 aTarget->SetStatus(0);
2463 theTargets.push_back(aTarget);
2464
2465 theProjectileSplitable->IncrementCollisionCount(1);
2466 theProjectileSplitable->SetStatus(0*theProjectileSplitable->GetStatus());
2467
2468 G4InteractionContent * aInteraction = new G4InteractionContent(theProjectileSplitable);
2469 aInteraction->SetTarget(aTarget);
2470 aInteraction->SetTargetNucleon(pNucleon);
2471 aInteraction->SetNumberOfSoftCollisions(1);
2472 aInteraction->SetStatus(3);
2473 theInteractions.push_back(aInteraction);
2474 }
2475 }
2476 return theProjectileSplitable;
2477 }
2478