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