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