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> << 27 << 28 #include "G4QGSParticipants.hh" << 29 #include "globals.hh" 26 #include "globals.hh" 30 #include "G4SystemOfUnits.hh" << 27 #include "G4QGSParticipants.hh" 31 #include "G4LorentzVector.hh" 28 #include "G4LorentzVector.hh" 32 #include "G4V3DNucleus.hh" << 29 #include <utility> 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 30 44 // Class G4QGSParticipants 31 // Class G4QGSParticipants 45 32 >> 33 // HPW Feb 1999 46 // Promoting model parameters from local varia 34 // Promoting model parameters from local variables class properties 47 G4ThreadLocal G4int G4QGSParticipants_NPart = << 35 G4int G4QGSParticipants_NPart = 0; 48 36 49 G4QGSParticipants::G4QGSParticipants() << 37 G4QGSParticipants::G4QGSParticipants() : theDiffExcitaton(0.7*GeV, 250*MeV, 250*MeV), 50 : theDiffExcitaton(), ModelMode(SOFT), nCutM << 38 nCutMax(7),ThersholdParameter(0.45*GeV), 51 ThresholdParameter(0.0*GeV), QGSMThreshold << 39 QGSMThershold(3*GeV),theNucleonRadius(1.5*fermi) 52 theNucleonRadius(1.5*fermi), theCurrentVel << 40 53 theProjectileSplitable(nullptr), Regge(nul << 54 InteractionMode(ALL), alpha(-0.5), beta(2. << 55 NumberOfInvolvedNucleonsOfTarget(0), Numbe << 56 ProjectileResidual4Momentum(G4LorentzVecto << 57 ProjectileResidualCharge(0), ProjectileRes << 58 TargetResidual4Momentum(G4LorentzVector()) << 59 TargetResidualCharge(0), TargetResidualExc << 60 CofNuclearDestruction(0.0), R2ofNuclearDes << 61 ExcitationEnergyPerWoundedNucleon(0.0), Do << 62 Pt2ofNuclearDestruction(0.0), MaxPt2ofNucl << 63 { 41 { 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 } 42 } 78 43 79 G4QGSParticipants::G4QGSParticipants(const G4Q 44 G4QGSParticipants::G4QGSParticipants(const G4QGSParticipants &right) 80 : G4VParticipants(), ModelMode(right.ModelMo << 45 : G4VParticipants(), nCutMax(right.nCutMax),ThersholdParameter(right.ThersholdParameter), 81 ThresholdParameter(right.ThresholdParamete << 46 QGSMThershold(right.QGSMThershold),theNucleonRadius(right.theNucleonRadius) 82 QGSMThreshold(right.QGSMThreshold), << 83 theNucleonRadius(right.theNucleonRadius), << 84 theCurrentVelocity(right.theCurrentVelocit << 85 theProjectileSplitable(right.theProjectile << 86 Regge(right.Regge), InteractionMode(right. << 87 alpha(right.alpha), beta(right.beta), sigm << 88 NumberOfInvolvedNucleonsOfTarget(right.Num << 89 NumberOfInvolvedNucleonsOfProjectile(right << 90 ProjectileResidual4Momentum(right.Projecti << 91 ProjectileResidualMassNumber(right.Project << 92 ProjectileResidualCharge(right.ProjectileR << 93 ProjectileResidualExcitationEnergy(right.P << 94 TargetResidual4Momentum(right.TargetResidu << 95 TargetResidualMassNumber(right.TargetResid << 96 TargetResidualCharge(right.TargetResidualC << 97 TargetResidualExcitationEnergy(right.Targe << 98 CofNuclearDestruction(right.CofNuclearDest << 99 R2ofNuclearDestruction(right.R2ofNuclearDe << 100 ExcitationEnergyPerWoundedNucleon(right.Ex << 101 DofNuclearDestruction(right.DofNuclearDest << 102 Pt2ofNuclearDestruction(right.Pt2ofNuclear << 103 MaxPt2ofNuclearDestruction(right.MaxPt2ofN << 104 { 47 { 105 for (size_t i=0; i < 250; ++i) { << 106 TheInvolvedNucleonsOfTarget[i] = right.The << 107 TheInvolvedNucleonsOfProjectile[i] = right << 108 } << 109 } 48 } 110 49 111 G4QGSParticipants::~G4QGSParticipants() {} << 112 50 113 void G4QGSParticipants::BuildInteractions(cons << 51 G4QGSParticipants::~G4QGSParticipants() 114 { 52 { 115 theProjectile = thePrimary; << 53 } 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 54 186 } while( (!Success) && ++internalLoopCount << 55 void G4QGSParticipants::BuildInteractions(const G4ReactionProduct &thePrimary) >> 56 { >> 57 >> 58 // Find the collisions and collition conditions >> 59 G4VSplitableHadron* aProjectile = SelectInteractions(thePrimary); 187 60 188 if ( internalLoopCounter >= maxNumberOfInt << 61 // now build the parton pairs. HPW 189 Success = false; << 62 SplitHadrons(); 190 } << 63 >> 64 // soft collisions first HPW, ordering is vital >> 65 PerformSoftCollisions(); 191 66 192 if ( Success ) { << 67 // the rest is diffractive HPW 193 #ifdef debugQGSParticipants << 68 PerformDiffractiveCollisions(); 194 G4cout<<G4endl<<"PerformDiffractiveCol << 69 195 #endif << 196 << 197 PerformDiffractiveCollisions(); << 198 << 199 #ifdef debugQGSParticipants << 200 G4cout<<G4endl<<"SplitHadrons();" <<G4 << 201 #endif << 202 << 203 SplitHadrons(); << 204 << 205 if( theProjectileSplitable && theProject << 206 #ifdef debugQGSParticipants << 207 G4cout<<"Perform non-Diffractive Col << 208 #endif << 209 Success = DeterminePartonMomenta(); << 210 } << 211 << 212 if(!Success) PrepareInitialState( thePri << 213 } << 214 } while( (!Success) && ++loopCounter < maxNu << 215 << 216 if ( loopCounter >= maxNumberOfLoops ) { << 217 Success = false; << 218 #ifdef debugQGSParticipants << 219 G4cout<<"NOT Successful ======" <<G4endl << 220 #endif << 221 } << 222 << 223 if ( Success ) { << 224 CreateStrings(); // To create strings << 225 << 226 GetResiduals(); // To calculate residual << 227 << 228 #ifdef debugQGSParticipants << 229 G4cout<<"All O.K. ======" <<G4endl; << 230 #endif << 231 } << 232 << 233 // clean-up, if necessary 70 // clean-up, if necessary 234 #ifdef debugQGSParticipants << 71 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(); 72 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 73 std::for_each(theTargets.begin(), theTargets.end(), DeleteSplitableHadron()); 269 theTargets.clear(); 74 theTargets.clear(); 270 << 75 delete aProjectile; 271 if ( theProjectileSplitable ) { << 272 delete theProjectileSplitable; << 273 theProjectileSplitable = 0; << 274 } << 275 } 76 } 276 77 277 //============================================ << 78 G4VSplitableHadron* G4QGSParticipants::SelectInteractions(const G4ReactionProduct &thePrimary) 278 void G4QGSParticipants::PrepareInitialState( c << 279 { 79 { 280 // Clearing of the arrays << 80 G4VSplitableHadron* aProjectile = new G4QGSMSplitableHadron(thePrimary, TRUE); // @@@ check the TRUE 281 // Erasing of the projectile << 81 G4PomeronCrossSection theProbability(thePrimary.GetDefinition()); // @@@ should be data member 282 G4InteractionContent* anIniteraction = theIn << 82 G4double outerRadius = theNucleus->GetOuterRadius(); 283 G4VSplitableHadron* pProjectile = anIniterac << 83 // Check reaction threshold 284 if( pProjectile ) delete pProjectile; << 285 84 286 std::for_each(theInteractions.begin(), theIn << 287 theInteractions.clear(); << 288 << 289 // Erasing of the envolved nucleons and targ << 290 theNucleus->StartLoop(); 85 theNucleus->StartLoop(); 291 G4Nucleon* aNucleon; << 86 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. 87 G4LorentzVector aPrimaryMomentum(thePrimary.GetMomentum(), thePrimary.GetTotalEnergy()); 345 G4LorentzVector aNucleonMomentum(0.,0.,0., 9 << 88 G4double s = (aPrimaryMomentum + pNucleon->Get4Momentum()).mag2(); 346 << 89 G4double ThresholdMass = thePrimary.GetMass() + pNucleon->GetDefinition()->GetPDGMass(); 347 G4double SS=(aPrimaryMomentum + aNucleonMome << 90 ModelMode = SOFT; 348 << 91 if (sqr(ThresholdMass + ThersholdParameter) > s) 349 Regge->SetS(SS); << 92 { 350 << 93 throw G4HadronicException(__FILE__, __LINE__, "Initial energy is too low. The 4-vectors of the input are inconsistant with the particle masses."); 351 //-------------------------------------- << 352 theNucleus->StartLoop(); << 353 G4Nucleon * tNucleon = theNucleus->GetNextNu << 354 << 355 if ( ! tNucleon ) { << 356 #ifdef debugQGSParticipants << 357 G4cout << "QGSM - BAD situation: pNucleo << 358 #endif << 359 return; << 360 } 94 } 361 << 95 if (sqr(ThresholdMass + QGSMThershold) > s) // thus only diffractive in cascade! 362 G4double theNucleusOuterR = theNucleus->GetO << 363 << 364 if (theNucleus->GetMassNumber() == 1) << 365 { 96 { 366 G4ThreeVector aPos = G4ThreeVector(0.,0.,0 << 97 ModelMode = DIFFRACTIVE; 367 tNucleon->SetPosition(aPos); << 368 theNucleusOuterR = 0.; << 369 } 98 } 370 << 99 371 // Determination of participating nucleons o << 100 // first find the collisions HPW 372 << 373 std::for_each(theInteractions.begin(), theIn 101 std::for_each(theInteractions.begin(), theInteractions.end(), DeleteInteractionContent()); 374 theInteractions.clear(); 102 theInteractions.clear(); >> 103 G4int totalCuts = 0; >> 104 G4double impactUsed = 0; 375 105 376 G4int MaxPower=thePrimary.GetMomentum().mag( << 106 #ifdef debug_QGS 377 << 107 G4double eK = thePrimary.GetKineticEnergy()/GeV; 378 const G4int maxNumberOfLoops = 1000; << 108 #endif 379 << 380 G4int NumberOfTries = 0; << 381 G4double Scale = 1.0; << 382 109 383 G4int loopCounter = -1; << 110 while(theInteractions.size() == 0) 384 while( (theInteractions.size() == 0) && ++lo << 385 { 111 { 386 InteractionMode = ALL; // Mode = ALL, WIT << 112 // choose random impact parameter HPW 387 << 388 // choose random impact parameter of a col << 389 std::pair<G4double, G4double> theImpactPar 113 std::pair<G4double, G4double> theImpactParameter; 390 << 114 theImpactParameter = theNucleus->ChooseImpactXandY(outerRadius+theNucleonRadius); 391 NumberOfTries++; << 115 G4double impactX = theImpactParameter.first; 392 if( NumberOfTries == 100*(NumberOfTries/10 << 393 << 394 theImpactParameter = theNucleus->ChooseImp << 395 G4double impactX = theImpactParameter.firs << 396 G4double impactY = theImpactParameter.seco 116 G4double impactY = theImpactParameter.second; 397 << 117 398 #ifdef debugQGSParticipants << 118 // loop over nuclei to find collissions HPW 399 G4cout<<"InteractionMode "<<InteractionM << 400 G4cout<<"Impact parameter (fm ) "<<std:: << 401 G4int nucleonCount = -1; << 402 #endif << 403 << 404 // loop over nucleons to find collisions << 405 theNucleus->StartLoop(); 119 theNucleus->StartLoop(); >> 120 G4int nucleonCount = 0; // debug 406 G4QGSParticipants_NPart = 0; 121 G4QGSParticipants_NPart = 0; 407 << 122 while( (pNucleon = theNucleus->GetNextNucleon()) ) 408 G4double Power=MaxPower; << 409 << 410 while( (tNucleon = theNucleus->GetNextNucl << 411 { 123 { 412 if(Power <= 0.) break; << 124 if(totalCuts>1.5*thePrimary.GetKineticEnergy()/GeV) 413 << 414 G4LorentzVector nucleonMomentum=tNucleon << 415 << 416 G4double Distance2 = sqr(impactX - tNucl << 417 sqr(impactY - tNucl << 418 << 419 G4double Pint(0.); // << 420 G4double Pprd(0.), Ptrd(0.), Pdd(0.); // << 421 G4double Pnd (0.), Pnvr(0.); // << 422 G4int NcutPomerons(0); // << 423 << 424 Regge->GetProbabilities(std::sqrt(Distan << 425 Pint, Pprd, Ptrd, Pdd, Pnd, Pnvr); << 426 #ifdef debugQGSParticipants << 427 nucleonCount++; << 428 G4cout<<"Nucleon & its impact paramete << 429 G4cout<<"Probability of interaction: << 430 G4cout<<"Probability of PrD, TrD, DD: "<< << 431 G4cout<<"Probability of NonDiff, QuarkExc.: << 432 #endif << 433 << 434 if (Pint > G4UniformRand()) << 435 { // An inte << 436 << 437 G4double rndNumber = G4UniformRand(); << 438 G4int InteractionType(0); << 439 << 440 if((InteractionMode==ALL)||(Interactio << 441 { << 442 if( rndNumber < Pprd ) { << 443 else if( rndNumber < Pprd+Ptrd) { << 444 else if( rndNumber < Pprd+Ptrd+Pdd) { << 445 else if( rndNumber < Pprd+Ptrd+Pdd+Pnd ) { << 446 NcutPomerons = Regge->ncPom << 447 else {Intera << 448 } << 449 else // InteractionMode == NON_DIFF << 450 { << 451 InteractionMode = NON_DIFF; << 452 if( rndNumber < Ptrd ) {Interact << 453 else if( rndNumber < Ptrd + Pnd) {Interact << 454 } << 455 << 456 if( (InteractionType == NonD) && (Ncut << 457 << 458 G4QGSParticipants_NPart ++; << 459 G4QGSMSplitableHadron* aTargetSPB = ne << 460 tNucleon->Hit(aTargetSPB); << 461 << 462 #ifdef debugQGSParticipants << 463 G4cout<<"An interaction is happend." << 464 G4cout<<"Target nucleon - "<<nucleon << 465 <<tNucleon->GetDefinition()->G << 466 G4cout<<"Interaction type:"<<Interac << 467 <<" (0 -PrD, 1 - TrD, 2 - DD, << 468 G4cout<<"New Inter. mode:"<<Interac << 469 <<" (0 -ALL, 1 - WITHOUT_R, 2 << 470 if( InteractionType == NonD ) << 471 G4cout<<"Number of cutted pomerons << 472 #endif << 473 << 474 if((InteractionType == PrD) || (Intera << 475 (InteractionType == Qexc)) << 476 { // << 477 #ifdef debugQGSParticipants << 478 G4cout<<"Diffractive-like interact << 479 #endif << 480 << 481 G4InteractionContent * aInteraction << 482 theProjectileSplitable->SetStatus(1* << 483 << 484 aInteraction->SetTarget(aTargetSPB); << 485 aInteraction->SetTargetNucleon(tNucl << 486 aTargetSPB->SetCollisionCount(0); << 487 aTargetSPB->SetStatus(1); << 488 << 489 aInteraction->SetNumberOfDiffractive << 490 aInteraction->SetNumberOfSoftCollisi << 491 aInteraction->SetStatus(InteractionT << 492 theInteractions.push_back(aInteracti << 493 } << 494 else << 495 { // non << 496 #ifdef debugQGSParticipants << 497 G4cout<<"Non-diffractive interacti << 498 #endif << 499 << 500 G4int nCuts; << 501 << 502 G4int Vncut=0; << 503 for(nCuts = 0; nCuts < NcutPomerons; nCuts << 504 { << 505 if( G4UniformRand() < Power/MaxPower ){V << 506 } << 507 nCuts=Vncut; << 508 << 509 if( nCuts == 0 ) {delete aTargetSPB; tNucl << 510 << 511 #ifdef debugQGSParticipants << 512 G4cout<<"Number of cuts in the int << 513 #endif << 514 << 515 aTargetSPB->IncrementCollisionCount(nCuts) << 516 aTargetSPB->SetStatus(0); << 517 theTargets.push_back(aTargetSPB); << 518 << 519 theProjectileSplitable->IncrementCollision << 520 theProjectileSplitable->SetStatus(0* << 521 << 522 G4InteractionContent * aInteraction = new << 523 aInteraction->SetTarget(aTargetSPB); << 524 aInteraction->SetTargetNucleon(tNucl << 525 aInteraction->SetNumberOfSoftCollisions(nC << 526 aInteraction->SetStatus(InteractionT << 527 theInteractions.push_back(aInteraction); << 528 } << 529 } // End of if (Pint > G4UniformRand( << 530 } // End of while( (tNucleon = theNucl << 531 << 532 #ifdef debugQGSParticipants << 533 G4cout << G4endl<<"Number of wounded nuc << 534 #endif << 535 << 536 } // End of while( (theInteractions.size() << 537 << 538 if ( loopCounter >= maxNumberOfLoops ) { << 539 #ifdef debugQGSParticipants << 540 G4cout <<"BAD situation: forced loop exi << 541 #endif << 542 // Perhaps there is something to set here. << 543 // Decrease impact parameter ?? << 544 // Select collisions with only diffraction << 545 // Selecy only non-diffractive interaction << 546 } << 547 //------------------------------------------ << 548 std::vector<G4InteractionContent*>::iterator << 549 << 550 if( theInteractions.size() != 0) << 551 { << 552 if( InteractionMode == ALL ) // It can be << 553 { // Only the << 554 i = theInteractions.end()-1; << 555 << 556 while ( theInteractions.size() != 1 ) << 557 { 125 { 558 G4InteractionContent* anInteraction = << 126 break; 559 G4Nucleon * pNucleon = anInteraction-> << 560 delete anInteraction->GetTarget(); << 561 delete *i; << 562 i=theInteractions.erase(i); << 563 i--; << 564 } 127 } 565 } << 128 nucleonCount++; // debug 566 else << 129 // Needs to be moved to Probability class @@@ 567 { // All quark << 130 G4double s = (aPrimaryMomentum + pNucleon->Get4Momentum()).mag2(); 568 i = theInteractions.begin(); << 131 G4double Distance2 = sqr(impactX - pNucleon->GetPosition().x()) + 569 while ( i != theInteractions.end() ) << 132 sqr(impactY - pNucleon->GetPosition().y()); >> 133 G4double Probability = theProbability.GetInelasticProbability(s, Distance2); >> 134 // test for inelastic collision >> 135 G4double rndNumber = G4UniformRand(); >> 136 // ModelMode = DIFFRACTIVE; >> 137 if (Probability > rndNumber) 570 { 138 { 571 G4InteractionContent* anInteraction = << 139 //--DEBUG-- cout << "DEBUG p="<< Probability<<" r="<<rndNumber<<" d="<<std::sqrt(Distance2)<<G4endl; 572 << 140 G4QGSMSplitableHadron* aTarget = new G4QGSMSplitableHadron(*pNucleon); 573 if( anInteraction->GetStatus() == Qexc << 141 G4QGSParticipants_NPart ++; 574 { << 142 theTargets.push_back(aTarget); 575 G4Nucleon* aTargetNucleon = a << 143 pNucleon->Hit(aTarget); 576 aTargetNucleon->Hit(nullptr); << 144 if ((theProbability.GetDiffractiveProbability(s, Distance2)/Probability > G4UniformRand() 577 << 145 &&(ModelMode==SOFT)) || (ModelMode==DIFFRACTIVE )) 578 delete anInteraction->GetTarget(); << 146 { 579 delete *i; << 147 // diffractive interaction occurs 580 i=theInteractions.erase(i); << 148 if(IsSingleDiffractive()) 581 } << 149 { 582 else << 150 theSingleDiffExcitation.ExciteParticipants(aProjectile, aTarget); 583 { << 151 } 584 i++; << 152 else 585 } << 153 { 586 } << 154 theDiffExcitaton.ExciteParticipants(aProjectile, aTarget); 587 } << 155 } 588 } << 156 G4InteractionContent * aInteraction = new G4InteractionContent(aProjectile); 589 } << 157 aInteraction->SetTarget(aTarget); 590 << 158 theInteractions.push_back(aInteraction); 591 //============================================ << 159 aInteraction->SetNumberOfDiffractiveCollisions(1); 592 void G4QGSParticipants::StoreInvolvedNucleon() << 160 totalCuts += 1; 593 { //To store nucleons involved in the interact << 594 << 595 NumberOfInvolvedNucleonsOfTarget = 0; << 596 << 597 theNucleus->StartLoop(); << 598 << 599 G4Nucleon* aNucleon; << 600 while ( ( aNucleon = theNucleus->GetNextNucl << 601 if ( aNucleon->AreYouHit() ) { << 602 TheInvolvedNucleonsOfTarget[NumberOfInvo << 603 NumberOfInvolvedNucleonsOfTarget++; << 604 } << 605 } << 606 << 607 #ifdef debugQGSParticipants << 608 G4cout << G4endl<<"G4QGSParticipants::Stor << 609 <<"Stored # of wounded nucleons of << 610 << NumberOfInvolvedNucleonsOfTarget << 611 #endif << 612 return; << 613 } << 614 << 615 //============================================ << 616 << 617 void G4QGSParticipants::ReggeonCascade() << 618 { // Implementation of the reggeon theory insp << 619 #ifdef debugQGSParticipants << 620 G4cout << G4endl<<"Reggeon cascading ..... << 621 G4cout<<"C of nucl. desctruction "<<GetCof << 622 <<" R2 "<<GetR2ofNuclearDestruction( << 623 #endif << 624 << 625 G4int InitNINt = NumberOfInvolvedNucleonsOfT << 626 << 627 // Reggeon cascading in target nucleus << 628 for ( G4int InvTN = 0; InvTN < InitNINt; Inv << 629 G4Nucleon* aTargetNucleon = TheInvolvedNuc << 630 << 631 G4double CreationTime = aTargetNucleon->Ge << 632 << 633 G4double XofWoundedNucleon = aTargetNucleo << 634 G4double YofWoundedNucleon = aTargetNucleo << 635 << 636 G4V3DNucleus* theTargetNucleus = theNucleu << 637 theTargetNucleus->StartLoop(); << 638 << 639 #ifdef debugQGSParticipants << 640 G4int TrgNuc=0; << 641 #endif << 642 G4Nucleon* Neighbour(0); << 643 while ( ( Neighbour = theTargetNucleus->Ge << 644 #ifdef debugQGSParticipants << 645 TrgNuc++; << 646 #endif << 647 if ( ! Neighbour->AreYouHit() ) { << 648 G4double impact2 = sqr( XofWoundedNucl << 649 sqr( YofWoundedNucl << 650 << 651 if ( G4UniformRand() < GetCofNuclearDe << 652 G4Exp( -impact2 << 653 ) { << 654 // The neighbour nucleon is involved << 655 #ifdef debugQGSParticipants << 656 G4cout<<"Target nucleon involved i << 657 <<Neighbour->GetDefinition() << 658 #endif << 659 TheInvolvedNucleonsOfTarget[ NumberO << 660 NumberOfInvolvedNucleonsOfTarget++; << 661 << 662 G4QGSMSplitableHadron* targetSplitab << 663 << 664 Neighbour->Hit( targetSplitable ); << 665 targetSplitable->SetTimeOfCreation( << 666 targetSplitable->SetStatus( 2 ); << 667 targetSplitable->SetCollisionCount(0 << 668 << 669 G4InteractionContent * anInteraction << 670 anInteraction->SetTarget(targetSplit << 671 anInteraction->SetTargetNucleon(Neig << 672 << 673 anInteraction->SetNumberOfDiffractiv << 674 anInteraction->SetNumberOfSoftCollis << 675 anInteraction->SetStatus(3); << 676 theInteractions.push_back(anInteract << 677 } << 678 } << 679 } << 680 } << 681 << 682 #ifdef debugQGSParticipants << 683 G4cout <<"Number of new involved nucleons << 684 #endif << 685 return; << 686 } << 687 << 688 //============================================ << 689 << 690 G4bool G4QGSParticipants::PutOnMassShell() { << 691 << 692 G4bool isProjectileNucleus = false; << 693 if ( GetProjectileNucleus() ) { << 694 isProjectileNucleus = true; << 695 } << 696 << 697 #ifdef debugPutOnMassShell << 698 G4cout <<G4endl<< "PutOnMassShell start .. << 699 if ( isProjectileNucleus ) {G4cout << "Put << 700 #endif << 701 << 702 G4LorentzVector Pprojectile( theProjectile.G << 703 if ( Pprojectile.z() < 0.0 ) { << 704 return false; << 705 } << 706 << 707 G4bool isOk = true; << 708 << 709 G4LorentzVector Ptarget( 0.0, 0.0, 0.0, 0.0 << 710 G4LorentzVector PtargetResidual( 0.0, 0.0, 0 << 711 G4double SumMasses = 0.0; << 712 G4V3DNucleus* theTargetNucleus = GetTargetNu << 713 G4double TargetResidualMass = 0.0; << 714 << 715 #ifdef debugPutOnMassShell << 716 G4cout << "Target : "; << 717 #endif << 718 << 719 isOk = ComputeNucleusProperties( theTargetNu << 720 TargetResid << 721 TargetResid << 722 << 723 if ( ! isOk ) return false; << 724 << 725 G4double Mprojectile = 0.0; << 726 G4double M2projectile = 0.0; << 727 G4LorentzVector Pproj( 0.0, 0.0, 0.0, 0.0 ); << 728 G4LorentzVector PprojResidual( 0.0, 0.0, 0.0 << 729 G4V3DNucleus* thePrNucleus = GetProjectileNu << 730 G4double PrResidualMass = 0.0; << 731 << 732 if ( ! isProjectileNucleus ) { // hadron-nu << 733 Mprojectile = Pprojectile.mag(); << 734 M2projectile = Pprojectile.mag2(); << 735 SumMasses += Mprojectile + 20.0*MeV; << 736 } else { // nucleus-nucleus or antinucleus- << 737 << 738 #ifdef debugPutOnMassShell << 739 G4cout << "Projectile : "; << 740 #endif << 741 << 742 isOk = ComputeNucleusProperties( thePrNucl << 743 Projectil << 744 Projectil << 745 if ( ! isOk ) return false; << 746 } << 747 << 748 G4LorentzVector Psum = Pprojectile + Ptarget << 749 G4double SqrtS = Psum.mag(); << 750 G4double S = Psum.mag2(); << 751 << 752 #ifdef debugPutOnMassShell << 753 G4cout << "Pproj "<<Pprojectile<<G4endl; << 754 G4cout << "Ptarg "<<Ptarget<<G4endl; << 755 G4cout << "Psum " << Psum/GeV << " GeV" << << 756 << "SumMasses, PrResidualMass and T << 757 << PrResidualMass/GeV << " " << Tar << 758 G4cout << "Ptar res. "<<PtargetResidual<<G << 759 #endif << 760 << 761 if ( SqrtS < SumMasses ) { << 762 return false; // It is impossible to simu << 763 } << 764 << 765 // Try to consider also the excitation energ << 766 // possible, with the available energy; othe << 767 << 768 G4double savedSumMasses = SumMasses; << 769 if ( isProjectileNucleus ) { << 770 SumMasses -= std::sqrt( sqr( PrResidualMas << 771 SumMasses += std::sqrt( sqr( PrResidualMas << 772 + PprojResidual.pe << 773 } << 774 SumMasses -= std::sqrt( sqr( TargetResidualM << 775 SumMasses += std::sqrt( sqr( TargetResidualM << 776 + PtargetResidual.pe << 777 << 778 if ( SqrtS < SumMasses ) { << 779 SumMasses = savedSumMasses; << 780 if ( isProjectileNucleus ) { << 781 ProjectileResidualExcitationEnergy = 0.0 << 782 } << 783 TargetResidualExcitationEnergy = 0.0; << 784 } << 785 << 786 TargetResidualMass += TargetResidualExcitati << 787 if ( isProjectileNucleus ) { << 788 PrResidualMass += ProjectileResidualExcita << 789 } << 790 << 791 #ifdef debugPutOnMassShell << 792 if ( isProjectileNucleus ) { << 793 G4cout << "PrResidualMass ProjResidualEx << 794 << ProjectileResidualExcitationEnergy < << 795 } << 796 G4cout << "TargetResidualMass TargetResidu << 797 << TargetResidualExcitationEnergy < << 798 << "Sum masses " << SumMasses/GeV < << 799 #endif << 800 << 801 // Sampling of nucleons what can transfer to << 802 if ( isProjectileNucleus && thePrNucleus-> << 803 isOk = GenerateDeltaIsobar( SqrtS, NumberO << 804 TheInvolvedNuc << 805 } << 806 if ( theTargetNucleus->GetMassNumber() != 1 << 807 isOk = isOk && << 808 GenerateDeltaIsobar( SqrtS, NumberO << 809 TheInvolvedNuc << 810 } << 811 if ( ! isOk ) return false; << 812 << 813 // Now we know that it is kinematically poss << 814 // of the involved nucleons (or correspondin << 815 // We have to sample the kinematical variabl << 816 // of the final state. The sampled kinematic << 817 // Notice that the sampling of the transvers << 818 // Fermi motion. << 819 << 820 // If target is nucleon - return ? << 821 << 822 G4LorentzRotation toCms( -1*Psum.boostVector << 823 G4LorentzVector Ptmp = toCms*Pprojectile; << 824 if ( Ptmp.pz() <= 0.0 ) { // "String" movin << 825 return false; << 826 } << 827 << 828 G4LorentzRotation toLab( toCms.inverse() ); << 829 << 830 G4double YprojectileNucleus = 0.0; << 831 if ( isProjectileNucleus ) { << 832 Ptmp = toCms*Pproj; << 833 YprojectileNucleus = Ptmp.rapidity(); << 834 } << 835 Ptmp = toCms*Ptarget; << 836 G4double YtargetNucleus = Ptmp.rapidity(); << 837 << 838 // Ascribing of the involved nucleons Pt and << 839 G4double DcorP = 0.0; << 840 if ( isProjectileNucleus ) { << 841 DcorP = GetDofNuclearDestruction() / thePr << 842 } << 843 G4double DcorT = GetDofNuclearDestruct << 844 G4double AveragePt2 = GetPt2ofNuclearDestru << 845 G4double maxPtSquare = GetMaxPt2ofNuclearDes << 846 << 847 #ifdef debugPutOnMassShell << 848 if ( isProjectileNucleus ) { << 849 G4cout << "Y projectileNucleus " << Ypro << 850 } << 851 G4cout << "Y targetNucleus " << Ytarge << 852 << "Dcor " << GetDofNuclearDestruct << 853 << " DcorP DcorT " << DcorP << " " << 854 #endif << 855 << 856 G4double M2proj = M2projectile; // Initiali << 857 G4double WplusProjectile = 0.0; << 858 G4double M2target = 0.0; << 859 G4double WminusTarget = 0.0; << 860 G4int NumberOfTries = 0; << 861 G4double ScaleFactor = 1.0; << 862 G4bool OuterSuccess = true; << 863 << 864 const G4int maxNumberOfLoops = 1000; << 865 G4int loopCounter = 0; << 866 do { << 867 G4double sqrtM2proj = 0.0, sqrtM2target = << 868 OuterSuccess = true; << 869 const G4int maxNumberOfTries = 1000; << 870 do { << 871 NumberOfTries++; << 872 if ( NumberOfTries == 100*(NumberOfTries << 873 // After many tries, it is convenient << 874 // AveragePt2, so that the sampled mom << 875 // involved nucleons (or corresponding delta << 876 // it is more likely to satisfy the mo << 877 ScaleFactor /= 2.0; << 878 DcorP *= ScaleFactor; << 879 DcorT *= ScaleFactor; << 880 AveragePt2 *= ScaleFactor; << 881 } << 882 if ( isProjectileNucleus ) { << 883 // Sampling of kinematical properties << 884 isOk = SamplingNucleonKinematics( Aver << 885 theP << 886 PrRe << 887 Numb << 888 TheI << 889 } << 890 // Sampling of kinematical properties of << 891 isOk = isOk && << 892 SamplingNucleonKinematics( Averag << 893 theTar << 894 Target << 895 Number << 896 TheInv << 897 << 898 if ( M2proj < 0.0 ) { << 899 if( M2proj < -0.000001 ) { << 900 G4ExceptionDescription ed; << 901 ed << "Projectile " << theProjectile.GetDe << 902 << " Target (Z,A)=(" << theTargetNucle << 903 << ") M2proj=" << M2proj << " -> << 904 G4Exception( "G4QGSParticipants::PutOnMass << 905 "HAD_QGSPARTICIPANTS_002", JustWarn << 906 } 161 } 907 M2proj = 0.0; << 162 else 908 } << 163 { 909 sqrtM2proj = std::sqrt( M2proj ); << 164 // nondiffractive soft interaction occurs 910 if ( M2target < 0.0 ) { << 165 // sample nCut+1 (cut Pomerons) pairs of strings can be produced 911 G4ExceptionDescription ed; << 166 G4int nCut; 912 ed << "Projectile " << theProjectile.G << 167 G4double * running = new G4double[nCutMax]; 913 << " Target (Z,A)=(" << theTargetN << 168 running[0] = 0; 914 << ") M2target=" << M2target << " << 169 for(nCut = 0; nCut < nCutMax; nCut++) 915 G4Exception( "G4QGSParticipants::PutOn << 170 { 916 "HAD_QGSPARTICIPANTS_003", << 171 running[nCut] = theProbability.GetCutPomeronProbability(s, Distance2, nCut + 1); 917 M2target = 0.0; << 172 if(nCut!=0) running[nCut] += running[nCut-1]; 918 }; << 173 } 919 sqrtM2target = std::sqrt( M2target ); << 174 G4double random = running[nCutMax-1]*G4UniformRand(); 920 << 175 for(nCut = 0; nCut < nCutMax; nCut++) 921 #ifdef debugPutOnMassShell << 176 { 922 G4cout << "SqrtS, Mp+Mt, Mp, Mt " << S << 177 if(running[nCut] > random) break; 923 << ( sqrtM2proj + sqrtM2target << 178 } 924 << sqrtM2proj/GeV << " " << sqr << 179 delete [] running; 925 #endif << 180 nCut = 0; 926 << 181 aTarget->IncrementCollisionCount(nCut+1); 927 if ( ! isOk ) return false; << 182 aProjectile->IncrementCollisionCount(nCut+1); 928 } while ( ( SqrtS < ( sqrtM2proj + sqrtM2t << 183 G4InteractionContent * aInteraction = new G4InteractionContent(aProjectile); 929 ++NumberOfTries < maxNumberOfTri << 184 aInteraction->SetTarget(aTarget); 930 if ( NumberOfTries >= maxNumberOfTries ) { << 185 aInteraction->SetNumberOfSoftCollisions(nCut+1); 931 return false; << 186 theInteractions.push_back(aInteraction); 932 } << 187 totalCuts += nCut+1; 933 if ( isProjectileNucleus ) { << 188 impactUsed=Distance2; 934 isOk = CheckKinematics( S, SqrtS, M2proj << 189 } 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 } 190 } 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 } 191 } 1291 << 192 //--DEBUG-- cout << G4endl<<"NUCLEONCOUNT "<<nucleonCount<<G4endl; 1292 #ifdef debugPutOnMassShell << 1293 G4cout << "success " << success << G4en << 1294 #endif << 1295 << 1296 } while ( ( ! success ) && << 1297 ++loopCounter < maxNumberOfLoops << 1298 if ( loopCounter >= maxNumberOfLoops ) { << 1299 success = false; << 1300 } << 1301 << 1302 return success; << 1303 } << 1304 << 1305 << 1306 //=========================================== << 1307 << 1308 G4bool G4QGSParticipants:: << 1309 CheckKinematics( const G4double sValue, << 1310 const G4double sqrtS, << 1311 const G4double projectileMas << 1312 const G4double targetMass2, << 1313 const G4double nucleusY, << 1314 const G4bool isProjectileNuc << 1315 const G4int numberOfInvolved << 1316 G4Nucleon* involvedNucleons[ << 1317 G4double& targetWminus, << 1318 G4double& projectileWplus, << 1319 G4bool& success ) { << 1320 << 1321 // This method, which is called only by Put << 1322 // kinematics is acceptable or not. << 1323 // This method assumes that all the paramet << 1324 // notice that the input boolean parameter << 1325 // only in the case of nucleus or antinucle << 1326 // The action of this method consists in co << 1327 // and setting the parameter success to fal << 1328 // be rejeted. << 1329 << 1330 G4double decayMomentum2 = sqr( sValue ) + s << 1331 - 2.0*sValue*proj << 1332 - 2.0*projectileM << 1333 targetWminus = ( sValue - projectileMass2 + << 1334 / 2.0 / sqrtS; << 1335 projectileWplus = sqrtS - targetMass2/targe << 1336 G4double projectilePz = projectileWplus/2.0 << 1337 G4double projectileE = projectileWplus/2.0 << 1338 G4double projectileY(1.0e5); << 1339 if (projectileE - projectilePz > 0.) { << 1340 projectileY = 0.5 * G4Log( (proje << 1341 (proje << 1342 } << 1343 G4double targetPz = -targetWminus/2.0 + tar << 1344 G4double targetE = targetWminus/2.0 + tar << 1345 G4double targetY = 0.5 * G4Log( (targetE + << 1346 << 1347 #ifdef debugPutOnMassShell << 1348 G4cout << "decayMomentum2 " << decayMomen << 1349 << "\t targetWminus projectileWplu << 1350 << "\t projectileY targetY " << pr << 1351 #endif << 1352 << 1353 for ( G4int i = 0; i < numberOfInvolvedNucl << 1354 G4Nucleon* aNucleon = involvedNucleons[i] << 1355 if ( ! aNucleon ) continue; << 1356 G4LorentzVector tmp = aNucleon->Get4Momen << 1357 G4double mt2 = sqr( tmp.x() ) + sqr( tmp. << 1358 sqr( aNucleon->GetSplitabl << 1359 G4double x = tmp.z(); << 1360 G4double pz = -targetWminus*x/2.0 + mt2/( << 1361 G4double e = targetWminus*x/2.0 + mt2/( << 1362 if ( isProjectileNucleus ) { << 1363 pz = projectileWplus*x/2.0 - mt2/(2.0*p << 1364 e = projectileWplus*x/2.0 + mt2/(2.0*p << 1365 } << 1366 G4double nucleonY = 0.5 * G4Log( (e + pz) << 1367 << 1368 #ifdef debugPutOnMassShell << 1369 G4cout << "i nY pY nY-AY AY " << i << " << 1370 #endif << 1371 << 1372 if ( std::abs( nucleonY - nucleusY ) > 2 << 1373 ( isProjectileNucleus && targetY > << 1374 ( ! isProjectileNucleus && project << 1375 success = false; << 1376 break; << 1377 } << 1378 } << 1379 return true; << 1380 } << 1381 << 1382 << 1383 //=========================================== << 1384 << 1385 G4bool G4QGSParticipants:: << 1386 FinalizeKinematics( const G4double w, << 1387 const G4bool isProjectile << 1388 const G4LorentzRotation& << 1389 const G4double residualMa << 1390 const G4int residualMassN << 1391 const G4int numberOfInvol << 1392 G4Nucleon* involvedNucleo << 1393 G4LorentzVector& residual4Momen << 1394 << 1395 // This method, which is called only by Put << 1396 // this method is called when we are sure t << 1397 // acceptable. << 1398 // This method assumes that all the paramet << 1399 // notice that the input boolean parameter << 1400 // only in the case of nucleus or antinucle << 1401 // because the sign of pz (in the center-of << 1402 // with respect to the case of a normal had << 1403 // The action of this method consists in mo << 1404 // (in the lab frame) and computing the res << 1405 // frame). << 1406 << 1407 G4ThreeVector residual3Momentum( 0.0, 0.0, << 1408 << 1409 for ( G4int i = 0; i < numberOfInvolvedNucl << 1410 G4Nucleon* aNucleon = involvedNucleons[i] << 1411 if ( ! aNucleon ) continue; << 1412 G4LorentzVector tmp = aNucleon->Get4Momen << 1413 residual3Momentum -= tmp.vect(); << 1414 G4double mt2 = sqr( tmp.x() ) + sqr( tmp. << 1415 sqr( aNucleon->GetSplitabl << 1416 G4double x = tmp.z(); << 1417 G4double pz = -w * x / 2.0 + mt2 / ( 2. << 1418 G4double e = w * x / 2.0 + mt2 / ( 2. << 1419 // Reverse the sign of pz in the case of << 1420 if ( isProjectileNucleus ) pz *= -1.0; << 1421 tmp.setPz( pz ); << 1422 tmp.setE( e ); << 1423 tmp.transform( boostFromCmsToLab ); << 1424 aNucleon->SetMomentum( tmp ); << 1425 G4VSplitableHadron* splitableHadron = aNu << 1426 splitableHadron->Set4Momentum( tmp ); << 1427 #ifdef debugPutOnMassShell << 1428 G4cout << "Target involved nucleon No, << 1429 << i<<" "<<aNucleon->GetDefinitio << 1430 #endif << 1431 } << 1432 << 1433 G4double residualMt2 = sqr( residualMass ) << 1434 + sqr( residual3Moment << 1435 << 1436 #ifdef debugPutOnMassShell << 1437 G4cout <<G4endl<< "w residual3Momentum.z( << 1438 #endif << 1439 << 1440 G4double residualPz = 0.0; << 1441 G4double residualE = 0.0; << 1442 if ( residualMassNumber != 0 ) { << 1443 residualPz = -w * residual3Momentum.z() / << 1444 residualMt2 / ( 2.0 * w * r << 1445 residualE = w * residual3Momentum.z() / << 1446 residualMt2 / ( 2.0 * w * r << 1447 // Reverse the sign of residualPz in the << 1448 if ( isProjectileNucleus ) residualPz *= << 1449 } 193 } 1450 << 194 //--DEBUG-- cout << G4endl<<"CUTDEBUG "<< totalCuts <<G4endl; 1451 residual4Momentum.setPx( residual3Momentum. << 195 //--DEBUG-- cout << "Impact Parameter used = "<<impactUsed<<G4endl; 1452 residual4Momentum.setPy( residual3Momentum. << 196 return aProjectile; 1453 residual4Momentum.setPz( residualPz ); << 1454 residual4Momentum.setE( residualE ); << 1455 << 1456 return true; << 1457 } 197 } 1458 198 1459 //=========================================== << 1460 void G4QGSParticipants::PerformDiffractiveCol 199 void G4QGSParticipants::PerformDiffractiveCollisions() 1461 { 200 { 1462 #ifdef debugQGSParticipants << 201 // remove the "G4PartonPair::PROJECTILE", etc., which are not necessary. @@@ 1463 G4cout<<G4endl<<"PerformDiffractiveCollisi << 1464 <<"theInteractions.size() "<<theInte << 1465 #endif << 1466 << 1467 unsigned int i; 202 unsigned int i; 1468 for (i = 0; i < theInteractions.size(); i++ << 203 for(i = 0; i < theInteractions.size(); i++) 1469 { 204 { 1470 G4InteractionContent* anIniteraction = th 205 G4InteractionContent* anIniteraction = theInteractions[i]; 1471 #ifdef debugQGSParticipants << 206 G4VSplitableHadron* aProjectile = anIniteraction->GetProjectile(); 1472 G4cout<<"Interaction # and its status " << 207 G4Parton* aParton = aProjectile->GetNextParton(); 1473 <<i<<" "<<theInteractions[i]->Get << 208 G4PartonPair * aPartonPair; 1474 #endif << 209 // projectile first HPW 1475 << 210 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 { << 1781 nSeaPair = (*i)->GetSoftCollisionCount( << 1782 #ifdef debugQGSParticipants << 1783 G4cout<<"nSeaPair of target N "<<nSea << 1784 <<"Target nucleon 4Mom "<<(*i)- << 1785 #endif << 1786 << 1787 SumPx = (*i)->Get4Momentum().px() * (-1 << 1788 SumPy = (*i)->Get4Momentum().py() * (-1 << 1789 SumZ = 0.; << 1790 << 1791 NumberOfUnsampledSeaQuarks = 2*nSeaPair << 1792 << 1793 Qmass=0; << 1794 for (G4int aSeaPair = 0; aSeaPair < nSe << 1795 { << 1796 aParton = (*i)->GetNextParton(); // << 1797 #ifdef debugQGSParticipants << 1798 G4cout<<"Sea quarks: "<<aSeaPair<<" << 1799 #endif << 1800 aPtVector = GaussianPt(SigPt, aHugeVa << 1801 tmp.setPx(aPtVector.x()); tmp.setPy(a << 1802 SumPx += aPtVector.x(); SumPy += aP << 1803 Mt=std::sqrt(aPtVector.mag2()+sqr(Qma << 1804 TargSumMt += Mt; << 1805 << 1806 // Sampling of Z fraction << 1807 tmp.setPz(SampleX(Xmin, NumberOfUnsam << 1808 SumZ += tmp.z(); << 1809 tmp.setPz((*i)->Get4Momentum().pz()*t << 1810 NumberOfUnsampledSeaQuarks--; << 1811 TargSumMt2perX +=sqr(Mt)/tmp.pz(); << 1812 tmp.setE(sqr(Mt)); << 1813 aParton->Set4Momentum(tmp); << 1814 << 1815 aParton = (*i)->GetNextAntiParton(); << 1816 #ifdef debugQGSParticipants << 1817 G4cout<<" "<<aParton->GetDefinition << 1818 G4cout<<" "<<tmp<<" "< << 1819 #endif << 1820 aPtVector = GaussianPt(SigPt, aHugeVa << 1821 tmp.setPx(aPtVector.x()); tmp.setPy(a << 1822 SumPx += aPtVector.x(); SumPy += aP << 1823 Mt=std::sqrt(aPtVector.mag2()+sqr(Qma << 1824 TargSumMt += Mt; << 1825 << 1826 // Sampling of Z fraction << 1827 tmp.setPz(SampleX(Xmin, NumberOfUnsam << 1828 SumZ += tmp.z(); << 1829 tmp.setPz((*i)->Get4Momentum().pz()*t << 1830 NumberOfUnsampledSeaQuarks--; << 1831 TargSumMt2perX +=sqr(Mt)/tmp.pz(); << 1832 tmp.setE(sqr(Mt)); << 1833 aParton->Set4Momentum(tmp); << 1834 #ifdef debugQGSParticipants << 1835 G4cout<<" "<<tmp<<" "< << 1836 #endif << 1837 } << 1838 << 1839 // Valence quark << 1840 aParton = (*i)->GetNextParton(); // f << 1841 #ifdef debugQGSParticipants << 1842 G4cout<<"Val quark of Tr"<<" "<<aPart << 1843 #endif << 1844 aPtVector = GaussianPt(SigPt, aHugeValu << 1845 tmp.setPx(aPtVector.x()); tmp.setPy(aPt << 1846 SumPx += aPtVector.x(); SumPy += aPtV << 1847 Mt=std::sqrt(aPtVector.mag2()+sqr(VqM_t << 1848 TargSumMt += Mt; << 1849 << 1850 // Sampling of Z fraction << 1851 tmp.setPz(SampleX(Xmin, NumberOfUnsampl << 1852 SumZ += tmp.z(); << 1853 tmp.setPz((*i)->Get4Momentum().pz()*tmp << 1854 TargSumMt2perX +=sqr(Mt)/tmp.pz(); << 1855 tmp.setE(sqr(Mt)); << 1856 aParton->Set4Momentum(tmp); << 1857 << 1858 // Valence di-quark << 1859 aParton = (*i)->GetNextAntiParton(); << 1860 #ifdef debugQGSParticipants << 1861 G4cout<<" "<<aParton->GetDefinition() << 1862 G4cout<<" "<<tmp<<" "<<S << 1863 #endif << 1864 tmp.setPx(-SumPx); tmp << 1865 //Uzhi 2019 Mt=std::sqrt(aPtVector.mag << 1866 Mt=std::sqrt( sqr(SumPx) + sqr(SumPy) + << 1867 TargSumMt += Mt; << 1868 << 1869 tmp.setPz((*i)->Get4Momentum().pz()*(1. << 1870 TargSumMt2perX +=sqr(Mt)/tmp.pz(); << 1871 tmp.setE(sqr(Mt)); << 1872 aParton->Set4Momentum(tmp); << 1873 #ifdef debugQGSParticipants << 1874 G4cout<<" "<<tmp<<" "<<1 << 1875 #endif << 1876 << 1877 } // End of for(i = theTargets.begin(); << 1878 << 1879 if( ProjSumMt + TargSumMt > Sqr << 1880 Success = false; continue;} << 1881 if( std::sqrt(ProjSumMt2perX) + std::sqrt << 1882 Success = false; continue;} << 1883 << 1884 } while( (!Success) && << 1885 attempt < maxNumberOfAttempts ); << 1886 << 1887 if ( attempt >= maxNumberOfAttempts ) { << 1888 return false; << 1889 } << 1890 << 1891 //+++++++++++++++++++++++++++++++++++++++++ << 1892 << 1893 G4double DecayMomentum2 = sqr(S) + sqr(Proj << 1894 - 2.0*S*ProjSumMt2perX - 2.0*S << 1895 << 1896 G4double targetWminus=( S - ProjSumMt2perX << 1897 G4double projectileWplus = SqrtS - TargSumM << 1898 << 1899 G4LorentzVector Tmp(0.,0.,0.,0.); << 1900 G4double z(0.); << 1901 << 1902 G4int nSeaPair = theProjectileSplitable->Ge << 1903 #ifdef debugQGSParticipants << 1904 G4cout<<"Backward transformation ======== << 1905 G4cout<<"nSeaPair of proj "<<nSeaPair<<G4 << 1906 #endif << 1907 << 1908 for (G4int aSeaPair = 0; aSeaPair < nSeaPai << 1909 { << 1910 aParton = theProjectileSplitable->GetNext << 1911 #ifdef debugQGSParticipants << 1912 G4cout<<"Sea quarks: "<<aSeaPair<<" "<< << 1913 #endif << 1914 Tmp =aParton->Get4Momentum(); z=Tmp.z(); << 1915 << 1916 Tmp.setPz(projectileWplus*z/2.0 - Tmp.e() << 1917 Tmp.setE( projectileWplus*z/2.0 + Tmp.e() << 1918 Tmp.transform( toLab ); << 1919 << 1920 aParton->Set4Momentum(Tmp); << 1921 << 1922 aParton = theProjectileSplitable->GetNext << 1923 #ifdef debugQGSParticipants << 1924 G4cout<<" "<<aParton->GetDefinition()-> << 1925 G4cout<<" "<<Tmp<<" "<<Tmp << 1926 #endif << 1927 Tmp =aParton->Get4Momentum(); z=Tmp.z(); << 1928 Tmp.setPz(projectileWplus*z/2.0 - Tmp.e() << 1929 Tmp.setE( projectileWplus*z/2.0 + Tmp.e() << 1930 Tmp.transform( toLab ); << 1931 << 1932 aParton->Set4Momentum(Tmp); << 1933 #ifdef debugQGSParticipants << 1934 G4cout<<" "<<Tmp<<" "<<Tmp << 1935 #endif << 1936 } << 1937 << 1938 // For valence quark << 1939 aParton = theProjectileSplitable->GetNextPa << 1940 #ifdef debugQGSParticipants << 1941 G4cout<<"Val quark of Pr"<<" "<<aParton-> << 1942 #endif << 1943 Tmp =aParton->Get4Momentum(); z=Tmp.z(); << 1944 Tmp.setPz(projectileWplus*z/2.0 - Tmp.e()/( << 1945 Tmp.setE( projectileWplus*z/2.0 + Tmp.e()/( << 1946 Tmp.transform( toLab ); << 1947 << 1948 aParton->Set4Momentum(Tmp); << 1949 << 1950 // For valence di-quark << 1951 aParton = theProjectileSplitable->GetNextAn << 1952 #ifdef debugQGSParticipants << 1953 G4cout<<" "<<aParton->GetDefinition()->Ge << 1954 G4cout<<" "<<Tmp<<" "<<Tmp.m << 1955 #endif << 1956 Tmp =aParton->Get4Momentum(); z=Tmp.z(); << 1957 Tmp.setPz(projectileWplus*z/2.0 - Tmp.e()/( << 1958 Tmp.setE( projectileWplus*z/2.0 + Tmp.e()/( << 1959 Tmp.transform( toLab ); << 1960 << 1961 aParton->Set4Momentum(Tmp); << 1962 << 1963 #ifdef debugQGSParticipants << 1964 G4cout<<" "<<Tmp<<" "<<Tmp.m << 1965 NuclNo=0; << 1966 #endif << 1967 << 1968 // End of work with the projectile << 1969 << 1970 // Work with target nucleons << 1971 for(i = theTargets.begin(); i != theTargets << 1972 { << 1973 nSeaPair = (*i)->GetSoftCollisionCount()- << 1974 #ifdef debugQGSParticipants << 1975 G4cout<<"nSeaPair of target and N# "<<n << 1976 NuclNo++; << 1977 #endif << 1978 for (G4int aSeaPair = 0; aSeaPair < nSeaP << 1979 { << 1980 aParton = (*i)->GetNextParton(); // f << 1981 #ifdef debugQGSParticipants << 1982 G4cout<<"Sea quarks: "<<aSeaPair<<" " << 1983 #endif << 1984 Tmp =aParton->Get4Momentum(); z=Tmp.z() << 1985 Tmp.setPz(-targetWminus*z/2.0 + Tmp.e() << 1986 Tmp.setE( targetWminus*z/2.0 + Tmp.e() << 1987 Tmp.transform( toLab ); << 1988 << 1989 aParton->Set4Momentum(Tmp); << 1990 << 1991 aParton = (*i)->GetNextAntiParton(); << 1992 #ifdef debugQGSParticipants << 1993 G4cout<<" "<<aParton->GetDefinition() << 1994 G4cout<<" "<<Tmp<<" "<<T << 1995 #endif << 1996 Tmp =aParton->Get4Momentum(); z=Tmp.z() << 1997 Tmp.setPz(-targetWminus*z/2.0 + Tmp.e() << 1998 Tmp.setE( targetWminus*z/2.0 + Tmp.e() << 1999 Tmp.transform( toLab ); << 2000 << 2001 aParton->Set4Momentum(Tmp); << 2002 #ifdef debugQGSParticipants << 2003 G4cout<<" "<<Tmp<<" "<<T << 2004 #endif << 2005 } << 2006 << 2007 // Valence quark << 2008 << 2009 aParton = (*i)->GetNextParton(); // for << 2010 #ifdef debugQGSParticipants << 2011 G4cout<<"Val quark of Tr"<<" "<<aParton << 2012 #endif << 2013 Tmp =aParton->Get4Momentum(); z=Tmp.z(); << 2014 Tmp.setPz(-targetWminus*z/2.0 + Tmp.e()/( << 2015 Tmp.setE( targetWminus*z/2.0 + Tmp.e()/( << 2016 Tmp.transform( toLab ); << 2017 << 2018 aParton->Set4Momentum(Tmp); << 2019 << 2020 // Valence di-quark << 2021 aParton = (*i)->GetNextAntiParton(); // << 2022 #ifdef debugQGSParticipants << 2023 G4cout<<" "<<aParton->GetDefinition()-> << 2024 G4cout<<" "<<Tmp<<" "<<Tmp << 2025 #endif << 2026 Tmp =aParton->Get4Momentum(); z=Tmp.z(); << 2027 Tmp.setPz(-targetWminus*z/2.0 + Tmp.e()/( << 2028 Tmp.setE( targetWminus*z/2.0 + Tmp.e()/( << 2029 Tmp.transform( toLab ); << 2030 << 2031 aParton->Set4Momentum(Tmp); << 2032 #ifdef debugQGSParticipants << 2033 G4cout<<" "<<Tmp<<" "<<Tmp << 2034 NuclNo++; << 2035 #endif << 2036 } // End of for(i = theTargets.begin(); i << 2037 << 2038 return true; << 2039 } << 2040 << 2041 //=========================================== << 2042 G4double G4QGSParticipants:: << 2043 SampleX(G4double, G4int nSea, G4int, G4double << 2044 { << 2045 G4double Oalfa = 1./(alpha + 1.); << 2046 G4double Obeta = 1./(aBeta + (alpha + 1.)*n << 2047 << 2048 G4double Ksi1, Ksi2, r1, r2, r12; << 2049 const G4int maxNumberOfLoops = 1000; << 2050 G4int loopCounter = 0; << 2051 do << 2052 { << 2053 Ksi1 = G4UniformRand(); r1 = G4Pow::GetIn << 2054 Ksi2 = G4UniformRand(); r2 = G4Pow::GetIn << 2055 r12=r1+r2; << 2056 } while( ( r12 > 1.) && << 2057 ++loopCounter < maxNumberOfLoops ) << 2058 if ( loopCounter >= maxNumberOfLoops ) { << 2059 return 0.5; // Just an acceptable value, << 2060 } << 2061 << 2062 G4double result = r1/r12; << 2063 return result; << 2064 } << 2065 << 2066 //=========================================== << 2067 void G4QGSParticipants::CreateStrings() << 2068 { << 2069 << 2070 #ifdef debugQGSParticipants << 2071 G4cout<<"CreateStrings() ................ << 2072 #endif << 2073 << 2074 if ( ! theProjectileSplitable ) { << 2075 #ifdef debugQGSParticipants << 2076 G4cout<<"BAD situation: theProjectileSp << 2077 #endif << 2078 return; << 2079 } << 2080 << 2081 #ifdef debugQGSParticipants << 2082 G4cout<<"theProjectileSplitable->GetStatu << 2083 G4LorentzVector str4Mom; << 2084 #endif << 2085 << 2086 if( theProjectileSplitable->GetStatus() == << 2087 { << 2088 G4ThreeVector Position = theProjectileSpl << 2089 << 2090 G4PartonPair * aPair = new G4PartonPair(t << 2091 t << 2092 G4PartonPair::DIFFRACTIVE << 2093 #ifdef debugQGSParticipants << 2094 G4cout << "Pr. Diffr. String: Qs 4mom X << 2095 G4cout << " " << aPair->Ge << 2096 << aPair->GetParton1()->Get4Momentum << 2097 << aPair->GetParton1()->GetX() << 2098 G4cout << " " << aPair->Ge << 2099 << aPair->GetParton2()->Get4Momentum << 2100 << aPair->GetParton2()->GetX() << 2101 str4Mom += aPair->GetParton1()->Get4Mom << 2102 str4Mom += aPair->GetParton2()->Get4Mom << 2103 #endif << 2104 << 2105 thePartonPairs.push_back(aPair); << 2106 } << 2107 << 2108 G4int N_EnvTarg = NumberOfInvolvedNucleonsO << 2109 << 2110 for ( G4int i = 0; i < N_EnvTarg; i++ ) { << 2111 G4Nucleon* aTargetNucleon = TheInvolvedNu << 2112 << 2113 G4VSplitableHadron* HitTargetNucleon = aT << 2114 << 2115 #ifdef debugQGSParticipants << 2116 G4cout<<"Involved Nucleon # and its sta << 2117 #endif << 2118 if( HitTargetNucleon->GetStatus() >= 1) << 2119 { 211 { 2120 G4ThreeVector Position = HitTargetN << 212 aPartonPair = new G4PartonPair(aParton, aProjectile->GetNextAntiParton(), 2121 << 213 G4PartonPair::DIFFRACTIVE, 2122 G4PartonPair * aPair = new G4PartonPair << 214 G4PartonPair::PROJECTILE); 2123 << 215 thePartonPairs.push_back(aPartonPair); 2124 G4PartonPair::DIFFRACTI << 2125 #ifdef debugQGSParticipants << 2126 G4cout << "Tr. Diffr. String: Qs 4mom << 2127 G4cout << "Diffr. String " << aPair->GetPar << 2128 << aPair->GetParton1()->Get4Moment << 2129 << aPair->GetParton1()->GetX() << 2130 G4cout << " " << aPair->GetPar << 2131 << aPair->GetParton2()->Get4Moment << 2132 << aPair->GetParton2()->GetX() << 2133 << 2134 str4Mom += aPair->GetParton1()->Get4Momentu << 2135 str4Mom += aPair->GetParton2()->Get4Momentu << 2136 #endif << 2137 << 2138 thePartonPairs.push_back(aPair); << 2139 } 216 } 2140 } << 217 // then target HPW 2141 << 218 G4VSplitableHadron* aTarget = anIniteraction->GetTarget(); 2142 //----------------------------------------- << 219 aParton = aTarget->GetNextParton(); 2143 #ifdef debugQGSParticipants << 220 if (aParton) 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 { 221 { 2161 G4VSplitableHadron* pProjectile = anIni << 222 aPartonPair = new G4PartonPair(aParton, aTarget->GetNextAntiParton(), 2162 G4VSplitableHadron* pTarget = anIni << 223 G4PartonPair::DIFFRACTIVE, 2163 << 224 G4PartonPair::TARGET); 2164 for (G4int j = 0; j < anIniteraction->G << 225 thePartonPairs.push_back(aPartonPair); 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 } 226 } 2206 } // End of while ( i != theInteractions.e << 2207 #ifdef debugQGSParticipants << 2208 G4cout << "Sum of strings 4 momenta " << << 2209 #endif << 2210 } << 2211 << 2212 //=========================================== << 2213 << 2214 void G4QGSParticipants::GetResiduals() { << 2215 // This method is needed for the correct ap << 2216 << 2217 #ifdef debugQGSParticipants << 2218 G4cout << "GetResiduals(): GetProjectileN << 2219 #endif << 2220 << 2221 #ifdef debugQGSParticipants << 2222 G4cout << "NumberOfInvolvedNucleonsOfTarg << 2223 #endif << 2224 << 2225 G4double DeltaExcitationE = TargetResidualE << 2226 G4LorentzVector DeltaPResidualNucleus = Tar << 2227 G4d << 2228 << 2229 for ( G4int i = 0; i < NumberOfInvolvedNucl << 2230 G4Nucleon* aNucleon = TheInvolvedNucleons << 2231 << 2232 #ifdef debugQGSParticipants << 2233 G4VSplitableHadron* targetSplitable = a << 2234 G4cout << i << " Hit? " << aNucleon->Ar << 2235 if ( targetSplitable ) G4cout << i << " << 2236 #endif << 2237 << 2238 G4LorentzVector tmp = -DeltaPResidualNucl << 2239 aNucleon->SetMomentum( tmp ); << 2240 aNucleon->SetBindingEnergy( DeltaExcitati << 2241 } 227 } 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 } 228 } 2328 229 2329 << 230 void G4QGSParticipants::PerformSoftCollisions() 2330 //=========================================== << 2331 void G4QGSParticipants::PerformSoftCollisions << 2332 { 231 { 2333 std::vector<G4InteractionContent*>::iterato 232 std::vector<G4InteractionContent*>::iterator i; 2334 G4LorentzVector str4Mom; << 233 for(i = theInteractions.begin(); i != theInteractions.end(); i++) 2335 i = theInteractions.begin(); << 2336 while ( i != theInteractions.end() ) /* Lo << 2337 { 234 { 2338 G4InteractionContent* anIniteraction = *i 235 G4InteractionContent* anIniteraction = *i; 2339 G4PartonPair * aPair = NULL; 236 G4PartonPair * aPair = NULL; 2340 if (anIniteraction->GetNumberOfSoftCollis 237 if (anIniteraction->GetNumberOfSoftCollisions()) 2341 { << 238 { 2342 G4VSplitableHadron* pProjectile = anIni 239 G4VSplitableHadron* pProjectile = anIniteraction->GetProjectile(); 2343 G4VSplitableHadron* pTarget = anIni 240 G4VSplitableHadron* pTarget = anIniteraction->GetTarget(); 2344 for (G4int j = 0; j < anIniteraction->G 241 for (G4int j = 0; j < anIniteraction->GetNumberOfSoftCollisions(); j++) 2345 { 242 { 2346 aPair = new G4PartonPair(pTarget->Get << 243 aPair = new G4PartonPair(pTarget->GetNextParton(), pProjectile->GetNextAntiParton(), 2347 G4PartonPair::SOFT, G4PartonPair::TA << 244 G4PartonPair::SOFT, G4PartonPair::TARGET); 2348 #ifdef debugQGSParticipants << 245 thePartonPairs.push_back(aPair); 2349 G4cout << "SoftPair " << aPair->GetParton << 246 aPair = new G4PartonPair(pProjectile->GetNextParton(), pTarget->GetNextAntiParton(), 2350 << aPair->GetParton1()->Get4Momentum() < << 247 G4PartonPair::SOFT, G4PartonPair::PROJECTILE); 2351 << aPair->GetParton1()->GetX() << " " << << 248 thePartonPairs.push_back(aPair); 2352 G4cout << " " << aPair->GetParton << 249 } 2353 << aPair->GetParton2()->Get4Momentum() < << 2354 << aPair->GetParton2()->GetX() << " " << << 2355 #endif << 2356 #ifdef debugQGSParticipants << 2357 str4Mom += aPair->GetParton1()->Get4Momen << 2358 str4Mom += aPair->GetParton2()->Get4Momen << 2359 #endif << 2360 thePartonPairs.push_back(aPair); << 2361 aPair = new G4PartonPair(pProjectile->GetNe << 2362 G4PartonPair::SOFT, G4PartonPair::PR << 2363 #ifdef debugQGSParticipants << 2364 G4cout << "SoftPair " << aPair->GetParton << 2365 << aPair->GetParton1()->Get4Momentum() < << 2366 << aPair->GetParton1()->GetX() << " " << << 2367 G4cout << " " << aPair->GetParton << 2368 << aPair->GetParton2()->Get4Momentum() < << 2369 << aPair->GetParton2()->GetX() << " " << << 2370 #endif << 2371 #ifdef debugQGSParticipants << 2372 str4Mom += aPair->GetParton1()->Get4Momen << 2373 str4Mom += aPair->GetParton2()->Get4Momen << 2374 #endif << 2375 thePartonPairs.push_back(aPair); << 2376 } << 2377 delete *i; 250 delete *i; 2378 i=theInteractions.erase(i); // i now << 251 i=theInteractions.erase(i); 2379 } else { << 252 i--; 2380 i++; << 2381 } << 2382 } << 2383 #ifdef debugQGSParticipants << 2384 G4cout << " string 4 mom " << str4Mom << << 2385 #endif << 2386 } << 2387 << 2388 << 2389 //******************************************* << 2390 G4VSplitableHadron* G4QGSParticipants::Select << 2391 { << 2392 // Check reaction threshold - goes to Chec << 2393 << 2394 theProjectileSplitable = new G4QGSMSplitabl << 2395 theProjectileSplitable->SetStatus(1); << 2396 << 2397 G4LorentzVector aPrimaryMomentum(thePrimary << 2398 G4LorentzVector aTargetNMomentum(0.,0.,0.,9 << 2399 << 2400 if ((!(aPrimaryMomentum.e()>-1)) && (!(aPri << 2401 { << 2402 throw G4HadronicException(__FILE__, __LIN << 2403 "G4GammaParticipants::SelectInter << 2404 } << 2405 G4double S = (aPrimaryMomentum + aTargetNMo << 2406 G4double ThresholdMass = thePrimary.GetMass << 2407 ModelMode = SOFT; << 2408 << 2409 #ifdef debugQGSParticipants << 2410 G4cout << "Gamma projectile - SelectInter << 2411 G4cout<<"Energy and Nucleus "<<thePrimary << 2412 G4cout << "SqrtS ThresholdMass ModelMode << 2413 #endif << 2414 << 2415 if (sqr(ThresholdMass + ThresholdParameter) << 2416 { << 2417 ModelMode = DIFFRACTIVE; << 2418 } << 2419 if (sqr(ThresholdMass + QGSMThreshold) > S) << 2420 { << 2421 ModelMode = DIFFRACTIVE; << 2422 } << 2423 << 2424 // first find the collisions HPW << 2425 std::for_each(theInteractions.begin(), theI << 2426 theInteractions.clear(); << 2427 << 2428 G4int theCurrent = G4int(theNucleus->GetMas << 2429 G4int NucleonNo=0; << 2430 << 2431 theNucleus->StartLoop(); << 2432 G4Nucleon * pNucleon = 0; << 2433 << 2434 while( (pNucleon = theNucleus->GetNextNucle << 2435 { if(NucleonNo == theCurrent) break; Nucleo << 2436 << 2437 if ( pNucleon ) { << 2438 G4QGSMSplitableHadron* aTarget = new G4QG << 2439 << 2440 pNucleon->Hit(aTarget); << 2441 << 2442 if ( (0.06 > G4UniformRand() &&(ModelMode << 2443 { << 2444 G4InteractionContent * aInteraction = n << 2445 theProjectileSplitable->SetStatus(1*the << 2446 << 2447 aInteraction->SetTarget(aTarget); << 2448 aInteraction->SetTargetNucleon(pNucleon << 2449 aTarget->SetCollisionCount(0); << 2450 aTarget->SetStatus(1); << 2451 << 2452 aInteraction->SetNumberOfDiffractiveCol << 2453 aInteraction->SetNumberOfSoftCollisions << 2454 aInteraction->SetStatus(1); << 2455 << 2456 theInteractions.push_back(aInteraction) << 2457 } << 2458 else << 2459 { << 2460 // nondiffractive soft interaction occu << 2461 aTarget->IncrementCollisionCount(1); << 2462 aTarget->SetStatus(0); << 2463 theTargets.push_back(aTarget); << 2464 << 2465 theProjectileSplitable->IncrementCollis << 2466 theProjectileSplitable->SetStatus(0*the << 2467 << 2468 G4InteractionContent * aInteraction = n << 2469 aInteraction->SetTarget(aTarget); << 2470 aInteraction->SetTargetNucleon(pNucleon << 2471 aInteraction->SetNumberOfSoftCollisions << 2472 aInteraction->SetStatus(3); << 2473 theInteractions.push_back(aInteraction) << 2474 } 253 } 2475 } 254 } 2476 return theProjectileSplitable; << 2477 } 255 } 2478 256