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