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 // 26 // 27 // ------------------------------------------- 27 // ------------------------------------------------------------ 28 // GEANT 4 class implemetation file 28 // GEANT 4 class implemetation file 29 // 29 // 30 // ---------------- G4QuarkExchange ----- 30 // ---------------- G4QuarkExchange -------------- 31 // by V. Uzhinsky, October 2016. 31 // by V. Uzhinsky, October 2016. 32 // QuarkExchange is used by strings mode 32 // QuarkExchange is used by strings models. 33 // Take a projectile and a target. 33 // Take a projectile and a target. 34 //Simulate Q exchange with excitation of proje 34 //Simulate Q exchange with excitation of projectile or target. 35 // ------------------------------------------- 35 // ------------------------------------------------------------ 36 36 37 #include "G4QuarkExchange.hh" 37 #include "G4QuarkExchange.hh" 38 #include "globals.hh" 38 #include "globals.hh" 39 #include "G4PhysicalConstants.hh" 39 #include "G4PhysicalConstants.hh" 40 #include "G4SystemOfUnits.hh" 40 #include "G4SystemOfUnits.hh" 41 #include "Randomize.hh" 41 #include "Randomize.hh" 42 #include "G4LorentzRotation.hh" 42 #include "G4LorentzRotation.hh" 43 #include "G4ThreeVector.hh" 43 #include "G4ThreeVector.hh" 44 #include "G4ParticleDefinition.hh" 44 #include "G4ParticleDefinition.hh" 45 #include "G4VSplitableHadron.hh" 45 #include "G4VSplitableHadron.hh" 46 #include "G4ExcitedString.hh" 46 #include "G4ExcitedString.hh" 47 47 48 #include "G4ParticleTable.hh" // Uzhi June 20 48 #include "G4ParticleTable.hh" // Uzhi June 2020 49 49 50 #include "G4Log.hh" 50 #include "G4Log.hh" 51 #include "G4Pow.hh" 51 #include "G4Pow.hh" 52 52 53 //#define debugQuarkExchange 53 //#define debugQuarkExchange 54 54 55 G4QuarkExchange::G4QuarkExchange() 55 G4QuarkExchange::G4QuarkExchange() 56 { 56 { 57 StrangeSuppress = (1.0-0.04)/2.0; // Uzhi J 57 StrangeSuppress = (1.0-0.04)/2.0; // Uzhi June 2020 : suppression of strange quark pair prodution, 58 // 58 // i.e. u:d:s=1:1:0.04 . Need to be tuned! 59 } 59 } 60 60 61 G4QuarkExchange::~G4QuarkExchange(){} 61 G4QuarkExchange::~G4QuarkExchange(){} 62 62 63 G4bool G4QuarkExchange:: 63 G4bool G4QuarkExchange:: 64 ExciteParticipants(G4VSplitableHadron *project 64 ExciteParticipants(G4VSplitableHadron *projectile, G4VSplitableHadron *target) const 65 { 65 { 66 #ifdef debugQuarkExchange 66 #ifdef debugQuarkExchange 67 G4cout<<G4endl<<"G4QuarkExchange::ExcitePa 67 G4cout<<G4endl<<"G4QuarkExchange::ExciteParticipants"<<G4endl; 68 #endif 68 #endif 69 69 70 G4LorentzVector Pprojectile = projectile->Ge 70 G4LorentzVector Pprojectile = projectile->Get4Momentum(); 71 G4double Mprojectile = projectile->Ge 71 G4double Mprojectile = projectile->GetDefinition()->GetPDGMass(); 72 G4double Mprojectile2 = sqr(Mprojectil 72 G4double Mprojectile2 = sqr(Mprojectile); 73 73 74 G4LorentzVector Ptarget = target->Get4Moment 74 G4LorentzVector Ptarget = target->Get4Momentum(); 75 G4double Mtarget = target->GetDefinit 75 G4double Mtarget = target->GetDefinition()->GetPDGMass(); 76 G4double Mtarget2 = sqr(Mtarget); 76 G4double Mtarget2 = sqr(Mtarget); 77 77 78 #ifdef debugQuarkExchange 78 #ifdef debugQuarkExchange 79 G4cout<<"Proj Targ "<<projectile->GetDefin 79 G4cout<<"Proj Targ "<<projectile->GetDefinition()->GetPDGEncoding()<<" "<<target->GetDefinition()->GetPDGEncoding()<<G4endl; 80 G4cout<<"Proj. 4-Mom "<<Pprojectile<<" "<< 80 G4cout<<"Proj. 4-Mom "<<Pprojectile<<" "<<Pprojectile.mag()<<G4endl 81 <<"Targ. 4-Mom "<<Ptarget <<" "<< 81 <<"Targ. 4-Mom "<<Ptarget <<" "<<Ptarget.mag() <<G4endl; 82 #endif 82 #endif 83 83 84 G4LorentzVector Psum=Pprojectile+Ptarget; 84 G4LorentzVector Psum=Pprojectile+Ptarget; 85 G4double SqrtS=Psum.mag(); 85 G4double SqrtS=Psum.mag(); 86 G4double S =Psum.mag2(); 86 G4double S =Psum.mag2(); 87 87 88 #ifdef debugQuarkExchange 88 #ifdef debugQuarkExchange 89 G4cout<<"SS Mpr Mtr SqrtS-Mprojectile-Mtar 89 G4cout<<"SS Mpr Mtr SqrtS-Mprojectile-Mtarget "<<SqrtS<<" "<<Mprojectile<<" "<<Mtarget 90 <<" "<<SqrtS-Mprojectile-Mtarget<<G4 90 <<" "<<SqrtS-Mprojectile-Mtarget<<G4endl; 91 #endif 91 #endif 92 if (SqrtS-Mprojectile-Mtarget <= 250.0*MeV) 92 if (SqrtS-Mprojectile-Mtarget <= 250.0*MeV) { 93 #ifdef debugQuarkExchange 93 #ifdef debugQuarkExchange 94 G4cerr<<"Energy is too small for quark e 94 G4cerr<<"Energy is too small for quark exchange!"<<G4endl; 95 G4cerr<<"Projectile: "<<projectile->GetD 95 G4cerr<<"Projectile: "<<projectile->GetDefinition()->GetPDGEncoding()<<" " 96 <<Pprojectile<<" "<<Pprojectile.ma 96 <<Pprojectile<<" "<<Pprojectile.mag()<<G4endl; 97 G4cerr<<"Target: "<<target->GetDefin 97 G4cerr<<"Target: "<<target->GetDefinition()->GetPDGEncoding()<<" " 98 <<Ptarget<<" "<<Ptarget.mag()<<G4e 98 <<Ptarget<<" "<<Ptarget.mag()<<G4endl; 99 G4cerr<<"sqrt(S) = "<<SqrtS<<" Mp + Mt = 99 G4cerr<<"sqrt(S) = "<<SqrtS<<" Mp + Mt = "<<Pprojectile.mag()+Ptarget.mag()<<G4endl; 100 #endif 100 #endif 101 return true; 101 return true; 102 } 102 } 103 103 104 G4LorentzRotation toCms(-1*Psum.boostVector( 104 G4LorentzRotation toCms(-1*Psum.boostVector()); 105 105 106 G4LorentzVector Ptmp=toCms*Pprojectile; 106 G4LorentzVector Ptmp=toCms*Pprojectile; 107 107 108 if ( Ptmp.pz() <= 0. ) 108 if ( Ptmp.pz() <= 0. ) 109 { 109 { 110 // "String" moving backwards in CMS, abor 110 // "String" moving backwards in CMS, abort collision !! 111 // G4cout << " abort Collision!! " 111 // G4cout << " abort Collision!! " << G4endl; 112 return false; 112 return false; 113 } 113 } 114 114 115 toCms.rotateZ(-1*Ptmp.phi()); 115 toCms.rotateZ(-1*Ptmp.phi()); 116 toCms.rotateY(-1*Ptmp.theta()); 116 toCms.rotateY(-1*Ptmp.theta()); 117 117 118 G4LorentzRotation toLab(toCms.inverse()); 118 G4LorentzRotation toLab(toCms.inverse()); 119 119 120 Pprojectile.transform(toCms); 120 Pprojectile.transform(toCms); 121 Ptarget.transform(toCms); 121 Ptarget.transform(toCms); 122 122 123 #ifdef debugQuarkExchange 123 #ifdef debugQuarkExchange 124 G4cout << "Pprojectile in CMS " << Pproje 124 G4cout << "Pprojectile in CMS " << Pprojectile << G4endl; 125 G4cout << "Ptarget in CMS " << Ptarge 125 G4cout << "Ptarget in CMS " << Ptarget << G4endl; 126 #endif 126 #endif 127 G4double maxPtSquare=sqr(Ptarget.pz()); 127 G4double maxPtSquare=sqr(Ptarget.pz()); 128 128 129 G4double ProjectileMinDiffrMass = Pprojectil 129 G4double ProjectileMinDiffrMass = Pprojectile.mag()/GeV; 130 G4double TargetMinDiffrMass = Ptarget.ma 130 G4double TargetMinDiffrMass = Ptarget.mag()/GeV; 131 131 132 G4double AveragePt2(0.); 132 G4double AveragePt2(0.); 133 133 134 G4int PDGcode=projectile->GetDefinition() 134 G4int PDGcode=projectile->GetDefinition()->GetPDGEncoding(); 135 G4int absPDGcode=std::abs(PDGcode); 135 G4int absPDGcode=std::abs(PDGcode); 136 136 137 G4bool ProjectileDiffraction = true; 137 G4bool ProjectileDiffraction = true; 138 138 139 // Also for heavy hadrons, assume 50% probab 139 // Also for heavy hadrons, assume 50% probability of projectile diffraction. 140 if ( absPDGcode > 1000 ) 140 if ( absPDGcode > 1000 ) { ProjectileDiffraction = G4UniformRand() <= 0.5; } 141 if ( (absPDGcode == 211) || (absPDGcode == 1 141 if ( (absPDGcode == 211) || (absPDGcode == 111) ) { ProjectileDiffraction = G4UniformRand() <= 0.66; } 142 if ( (absPDGcode == 321) || (absPDGcode == 3 142 if ( (absPDGcode == 321) || (absPDGcode == 311) || 143 ( PDGcode == 130) || ( PDGcode == 3 143 ( PDGcode == 130) || ( PDGcode == 310) ) { ProjectileDiffraction = G4UniformRand() <= 0.5; } 144 if ( absPDGcode > 400 && absPDGcode < 600 144 if ( absPDGcode > 400 && absPDGcode < 600 ) { ProjectileDiffraction = G4UniformRand() <= 0.5; } 145 145 146 //G4cout<<"ProjectileDiffr "<<ProjectileDiff 146 //G4cout<<"ProjectileDiffr "<<ProjectileDiffraction<<G4endl; 147 147 148 if ( ProjectileDiffraction ) { 148 if ( ProjectileDiffraction ) { 149 if ( absPDGcode > 1000 ) 149 if ( absPDGcode > 1000 ) //------Projectile is baryon -------- 150 { 150 { 151 if ( absPDGcode > 4000 && absPDGcode < 6 151 if ( absPDGcode > 4000 && absPDGcode < 6000 ) // Projectile is a charm or bottom baryon 152 { 152 { 153 ProjectileMinDiffrMass = projectile->G 153 ProjectileMinDiffrMass = projectile->GetDefinition()->GetPDGMass()/CLHEP::GeV + 0.25; // GeV 154 AveragePt2 = 0.3; 154 AveragePt2 = 0.3; // GeV^2 155 } 155 } 156 else 156 else 157 { 157 { 158 ProjectileMinDiffrMass = 1.16; 158 ProjectileMinDiffrMass = 1.16; // GeV 159 AveragePt2 = 0.3; 159 AveragePt2 = 0.3; // GeV^2 160 } 160 } 161 } 161 } 162 else if( absPDGcode == 211 || absPDGcode = 162 else if( absPDGcode == 211 || absPDGcode == 111) //------Projectile is Pion ----------- 163 { 163 { 164 ProjectileMinDiffrMass = 1.0; 164 ProjectileMinDiffrMass = 1.0; // GeV 165 AveragePt2 = 0.3; 165 AveragePt2 = 0.3; // GeV^2 166 } 166 } 167 else if( absPDGcode == 321 || absPDGcode = 167 else if( absPDGcode == 321 || absPDGcode == 130 || absPDGcode == 310) //Projectile is Kaon 168 { 168 { 169 ProjectileMinDiffrMass = 1.1; 169 ProjectileMinDiffrMass = 1.1; // GeV 170 AveragePt2 = 0.3; 170 AveragePt2 = 0.3; // GeV^2 171 } 171 } 172 else if( absPDGcode == 22) 172 else if( absPDGcode == 22) //------Projectile is Gamma ----------- 173 { 173 { 174 ProjectileMinDiffrMass = 0.25; 174 ProjectileMinDiffrMass = 0.25; // GeV 175 AveragePt2 = 0.36; 175 AveragePt2 = 0.36; // GeV^2 176 } 176 } 177 else if( absPDGcode > 400 && absPDGcode < 177 else if( absPDGcode > 400 && absPDGcode < 600) // Projectile is a charm or bottom meson 178 { 178 { 179 ProjectileMinDiffrMass = projectile->Get 179 ProjectileMinDiffrMass = projectile->GetDefinition()->GetPDGMass()/CLHEP::GeV + 0.25; // GeV 180 AveragePt2 = 0.3; 180 AveragePt2 = 0.3; // GeV^2 181 } 181 } 182 else 182 else //------Projectile is undefined, Nucleon assumed 183 { 183 { 184 ProjectileMinDiffrMass = 1.1; 184 ProjectileMinDiffrMass = 1.1; // GeV 185 AveragePt2 = 0.3; 185 AveragePt2 = 0.3; // GeV^2 186 }; 186 }; 187 187 188 ProjectileMinDiffrMass = ProjectileMinDiff 188 ProjectileMinDiffrMass = ProjectileMinDiffrMass * GeV; 189 Mprojectile2=sqr(ProjectileMinDiffrMass); 189 Mprojectile2=sqr(ProjectileMinDiffrMass); 190 190 191 if (G4UniformRand() <= 0.5) TargetMinDiffr 191 if (G4UniformRand() <= 0.5) TargetMinDiffrMass += 0.22; 192 TargetMinDiffrMass *= GeV; 192 TargetMinDiffrMass *= GeV; 193 Mtarget2 = sqr( TargetMinDiffrMass) ; 193 Mtarget2 = sqr( TargetMinDiffrMass) ; 194 } 194 } 195 else 195 else 196 { 196 { 197 if (G4UniformRand() <= 0.5) ProjectileMinD 197 if (G4UniformRand() <= 0.5) ProjectileMinDiffrMass += 0.22; 198 ProjectileMinDiffrMass *=GeV; 198 ProjectileMinDiffrMass *=GeV; 199 Mprojectile2=sqr(ProjectileMinDiffrMass); 199 Mprojectile2=sqr(ProjectileMinDiffrMass); 200 200 201 TargetMinDiffrMass = 1.16*GeV; 201 TargetMinDiffrMass = 1.16*GeV; // For target nucleon 202 Mtarget2 = sqr( TargetMinDiffrMass) ; 202 Mtarget2 = sqr( TargetMinDiffrMass) ; 203 AveragePt2 = 0.3; 203 AveragePt2 = 0.3; // GeV^2 204 } // end of if ( ProjectileDiffraction ) 204 } // end of if ( ProjectileDiffraction ) 205 205 206 AveragePt2 = AveragePt2 * GeV*GeV; 206 AveragePt2 = AveragePt2 * GeV*GeV; 207 207 208 if ( SqrtS - (ProjectileMinDiffrMass+TargetM 208 if ( SqrtS - (ProjectileMinDiffrMass+TargetMinDiffrMass) < 220.0*MeV ) return false; 209 209 210 //----------------------- 210 //----------------------- 211 G4double Pt2, PZcms, PZcms2; 211 G4double Pt2, PZcms, PZcms2; 212 G4double ProjMassT2, ProjMassT; 212 G4double ProjMassT2, ProjMassT; 213 G4double TargMassT2, TargMassT; 213 G4double TargMassT2, TargMassT; 214 G4double PMinusMin, PMinusMax, sqrtPMinusMi 214 G4double PMinusMin, PMinusMax, sqrtPMinusMin, sqrtPMinusMax; 215 G4double TPlusMin, TPlusMax, sqrtTPlusMin 215 G4double TPlusMin, TPlusMax, sqrtTPlusMin, sqrtTPlusMax; 216 G4double PMinusNew, PPlusNew, TPlusNew(0.), 216 G4double PMinusNew, PPlusNew, TPlusNew(0.), TMinusNew; 217 217 218 G4LorentzVector Qmomentum; 218 G4LorentzVector Qmomentum; 219 G4double Qminus, Qplus; 219 G4double Qminus, Qplus; 220 220 221 G4double x(0.), y(0.); 221 G4double x(0.), y(0.); 222 G4int whilecount=0; 222 G4int whilecount=0; 223 do { 223 do { 224 whilecount++; 224 whilecount++; 225 225 226 if (whilecount > 1000 ) 226 if (whilecount > 1000 ) 227 { 227 { 228 Qmomentum=G4LorentzVector(0.,0.,0.,0.); 228 Qmomentum=G4LorentzVector(0.,0.,0.,0.); 229 return false; // Ignore this intera 229 return false; // Ignore this interaction 230 } 230 } 231 231 232 // Generate pt 232 // Generate pt 233 Qmomentum=G4LorentzVector(GaussianPt(Avera 233 Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0); 234 234 235 Pt2 = G4ThreeVector( Qmomentum.vect() ).ma 235 Pt2 = G4ThreeVector( Qmomentum.vect() ).mag2(); 236 ProjMassT2 = Mprojectile2 + Pt2; 236 ProjMassT2 = Mprojectile2 + Pt2; 237 ProjMassT = std::sqrt( ProjMassT2 ); 237 ProjMassT = std::sqrt( ProjMassT2 ); 238 TargMassT2 = Mtarget2 + Pt2; 238 TargMassT2 = Mtarget2 + Pt2; 239 TargMassT = std::sqrt( TargMassT2 ); 239 TargMassT = std::sqrt( TargMassT2 ); 240 240 241 #ifdef debugQuarkExchange 241 #ifdef debugQuarkExchange 242 G4cout<<"whilecount Pt2 ProjMassT Tar 242 G4cout<<"whilecount Pt2 ProjMassT TargMassT SqrtS S ProjectileDiffraction"<<G4endl; 243 G4cout<<whilecount<<" "<<Pt2<<" "<<ProjM 243 G4cout<<whilecount<<" "<<Pt2<<" "<<ProjMassT<<" "<<TargMassT<<" "<<SqrtS<<" "<<S<<" "<<ProjectileDiffraction<<G4endl; 244 #endif 244 #endif 245 245 246 if ( SqrtS < ProjMassT + TargMassT + 220.0 246 if ( SqrtS < ProjMassT + TargMassT + 220.0*MeV ) continue; 247 247 248 PZcms2 = ( S*S + ProjMassT2*ProjMassT2 + T 248 PZcms2 = ( S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2 249 - 2.0*S*ProjMassT2 - 2.0*S*Targ 249 - 2.0*S*ProjMassT2 - 2.0*S*TargMassT2 - 2.0*ProjMassT2*TargMassT2 ) / 4.0 / S; 250 250 251 if ( PZcms2 < 0 ) continue; 251 if ( PZcms2 < 0 ) continue; 252 252 253 PZcms = std::sqrt( PZcms2 ); 253 PZcms = std::sqrt( PZcms2 ); 254 254 255 if ( ProjectileDiffraction ) 255 if ( ProjectileDiffraction ) 256 { // The projectile will fragment, the tar 256 { // The projectile will fragment, the target will saved. 257 PMinusMin = std::sqrt( ProjMassT2 + PZcm 257 PMinusMin = std::sqrt( ProjMassT2 + PZcms2 ) - PZcms; 258 PMinusMax = SqrtS - TargMassT; 258 PMinusMax = SqrtS - TargMassT; 259 sqrtPMinusMin = std::sqrt(PMinusMin); sq 259 sqrtPMinusMin = std::sqrt(PMinusMin); sqrtPMinusMax = std::sqrt(PMinusMax); 260 260 261 if ( absPDGcode > 1000 ) 261 if ( absPDGcode > 1000 ) 262 { 262 { 263 PMinusNew = PMinusMax * (1.0 - (1.0 - PMinus 263 PMinusNew = PMinusMax * (1.0 - (1.0 - PMinusMin/PMinusMax) 264 * G4Pow:: 264 * G4Pow::GetInstance()->powA(G4UniformRand(),0.3333) ); 265 } 265 } 266 else if ( (absPDGcode == 211) || (absPDG 266 else if ( (absPDGcode == 211) || (absPDGcode == 111) ) 267 { 267 { 268 while (true) 268 while (true) 269 { 269 { 270 x=sqrtPMinusMax-(sqrtPMinusMax-sqrtP 270 x=sqrtPMinusMax-(sqrtPMinusMax-sqrtPMinusMin)*G4UniformRand(); 271 y=G4UniformRand(); 271 y=G4UniformRand(); 272 if ( y < 1.0-0.7 * x/sqrtPMinusMax ) 272 if ( y < 1.0-0.7 * x/sqrtPMinusMax ) break; 273 } 273 } 274 PMinusNew = sqr(x); 274 PMinusNew = sqr(x); 275 } 275 } 276 else if ( (absPDGcode == 321) || (absPDG 276 else if ( (absPDGcode == 321) || (absPDGcode == 311) || 277 ( PDGcode == 130) || ( PDGcode = 277 ( PDGcode == 130) || ( PDGcode == 310) ) 278 { // For K-mesons it must be found !!! 278 { // For K-mesons it must be found !!! Uzhi 279 while (true) 279 while (true) 280 { 280 { 281 x=sqrtPMinusMax-(sqrtPMinusMax-sqrtPMinusM 281 x=sqrtPMinusMax-(sqrtPMinusMax-sqrtPMinusMin)*G4UniformRand(); 282 y=G4UniformRand(); 282 y=G4UniformRand(); 283 if ( y < 1.0-0.7 * x/sqrtPMinusMax ) 283 if ( y < 1.0-0.7 * x/sqrtPMinusMax ) break; 284 } 284 } 285 PMinusNew = sqr(x); 285 PMinusNew = sqr(x); 286 } 286 } 287 else 287 else 288 { 288 { 289 PMinusNew = PMinusMax * (1.0 - (1.0 - PMinus 289 PMinusNew = PMinusMax * (1.0 - (1.0 - PMinusMin/PMinusMax) 290 * G4Pow::GetI 290 * G4Pow::GetInstance()->powA(G4UniformRand(),0.3333) ); 291 }; 291 }; 292 292 293 TMinusNew = SqrtS - PMinusNew; 293 TMinusNew = SqrtS - PMinusNew; 294 294 295 Qminus = Ptarget.minus() - TMinusNew; 295 Qminus = Ptarget.minus() - TMinusNew; 296 TPlusNew = TargMassT2 / TMinusNew; 296 TPlusNew = TargMassT2 / TMinusNew; 297 Qplus = Ptarget.plus() - TPlusNew; 297 Qplus = Ptarget.plus() - TPlusNew; 298 298 299 } 299 } 300 else 300 else 301 { // The target will fragment, the project 301 { // The target will fragment, the projectile will saved. 302 TPlusMin = std::sqrt( TargMassT2 + PZcms 302 TPlusMin = std::sqrt( TargMassT2 + PZcms2 ) - PZcms; 303 TPlusMax = SqrtS - ProjMassT; 303 TPlusMax = SqrtS - ProjMassT; 304 sqrtTPlusMin = std::sqrt(TPlusMin); sqrt 304 sqrtTPlusMin = std::sqrt(TPlusMin); sqrtTPlusMax = std::sqrt(TPlusMax); 305 305 306 if ( absPDGcode > 1000 ) 306 if ( absPDGcode > 1000 ) 307 { 307 { 308 TPlusNew = TPlusMax * (1.0 - (1.0 - TP 308 TPlusNew = TPlusMax * (1.0 - (1.0 - TPlusMin/TPlusMax) 309 * G4Pow::G 309 * G4Pow::GetInstance()->powA(G4UniformRand(),0.3333) ); 310 } 310 } 311 else if ( (absPDGcode == 211) || (absPDG 311 else if ( (absPDGcode == 211) || (absPDGcode == 111) ) 312 { 312 { 313 while (true) 313 while (true) 314 { 314 { 315 x=sqrtTPlusMax-(sqrtTPlusMax-sqrtTPl 315 x=sqrtTPlusMax-(sqrtTPlusMax-sqrtTPlusMin)*G4UniformRand(); 316 y=G4UniformRand(); 316 y=G4UniformRand(); 317 if ( y < 1.0-0.7 * x/sqrtTPlusMax ) 317 if ( y < 1.0-0.7 * x/sqrtTPlusMax ) break; 318 } 318 } 319 TPlusNew = sqr(x); 319 TPlusNew = sqr(x); 320 } 320 } 321 else if ( (absPDGcode == 321) || (absPDG 321 else if ( (absPDGcode == 321) || (absPDGcode == 311) || 322 ( PDGcode == 130) || ( PDGcode = 322 ( PDGcode == 130) || ( PDGcode == 310) ) 323 { // For K-mesons it must be found !!! U 323 { // For K-mesons it must be found !!! Uzhi 324 while (true) 324 while (true) 325 { 325 { 326 x=sqrtTPlusMax-(sqrtTPlusMax-sqrtTPl 326 x=sqrtTPlusMax-(sqrtTPlusMax-sqrtTPlusMin)*G4UniformRand(); 327 y=G4UniformRand(); 327 y=G4UniformRand(); 328 if ( y < 1.0-0.7 * x/sqrtTPlusMax ) break; 328 if ( y < 1.0-0.7 * x/sqrtTPlusMax ) break; 329 } 329 } 330 } 330 } 331 else 331 else 332 { 332 { 333 TPlusNew = TPlusMax * (1.0 - (1.0 - TP 333 TPlusNew = TPlusMax * (1.0 - (1.0 - TPlusMin/TPlusMax) 334 * G4Pow::G 334 * G4Pow::GetInstance()->powA(G4UniformRand(),0.3333) ); 335 }; 335 }; 336 336 337 PPlusNew = SqrtS - TPlusNew; 337 PPlusNew = SqrtS - TPlusNew; 338 338 339 Qplus = PPlusNew - Pprojectile.plus(); 339 Qplus = PPlusNew - Pprojectile.plus(); 340 PMinusNew = ProjMassT2 / PPlusNew; 340 PMinusNew = ProjMassT2 / PPlusNew; 341 Qminus = PMinusNew - Pprojectile.minus() 341 Qminus = PMinusNew - Pprojectile.minus(); 342 } 342 } 343 343 344 Qmomentum.setPz( (Qplus - Qminus)/2 ); 344 Qmomentum.setPz( (Qplus - Qminus)/2 ); 345 Qmomentum.setE( (Qplus + Qminus)/2 ); 345 Qmomentum.setE( (Qplus + Qminus)/2 ); 346 346 347 #ifdef debugQuarkExchange 347 #ifdef debugQuarkExchange 348 G4cout<<"ProjectileDiffraction (Pproject 348 G4cout<<"ProjectileDiffraction (Pprojectile + Qmomentum).mag2() Mprojectile2"<<G4endl; 349 G4cout<<ProjectileDiffraction<<" "<<( Pp 349 G4cout<<ProjectileDiffraction<<" "<<( Pprojectile + Qmomentum ).mag2()<<" "<< Mprojectile2<<G4endl; 350 G4cout<<"TargetDiffraction (Ptarget 350 G4cout<<"TargetDiffraction (Ptarget - Qmomentum).mag2() Mtarget2"<<G4endl; 351 G4cout<<!ProjectileDiffraction<<" "<<( P 351 G4cout<<!ProjectileDiffraction<<" "<<( Ptarget - Qmomentum ).mag2()<<" "<< Mtarget2<<G4endl; 352 #endif 352 #endif 353 353 354 } while ( ( ProjectileDiffraction&&( Pprojec 354 } while ( ( ProjectileDiffraction&&( Pprojectile + Qmomentum ).mag2() < Mprojectile2 ) || 355 (!ProjectileDiffraction&&( Ptarget 355 (!ProjectileDiffraction&&( Ptarget - Qmomentum ).mag2() < Mtarget2 ) ); 356 // Repeat the sampling because there was n 356 // Repeat the sampling because there was not any excitation 357 357 358 Pprojectile += Qmomentum; 358 Pprojectile += Qmomentum; 359 359 360 Ptarget -= Qmomentum; 360 Ptarget -= Qmomentum; 361 361 362 // Transform back and update SplitableHadron 362 // Transform back and update SplitableHadron Participant. 363 Pprojectile.transform(toLab); 363 Pprojectile.transform(toLab); 364 Ptarget.transform(toLab); 364 Ptarget.transform(toLab); 365 365 366 #ifdef debugQuarkExchange 366 #ifdef debugQuarkExchange 367 G4cout << "Pprojectile in Lab. " << Pproj 367 G4cout << "Pprojectile in Lab. " << Pprojectile << G4endl; 368 G4cout << "Ptarget in Lab. " << Ptarg 368 G4cout << "Ptarget in Lab. " << Ptarget << G4endl; 369 G4cout << "G4QuarkExchange: Projectile mas 369 G4cout << "G4QuarkExchange: Projectile mass " << Pprojectile.mag() << G4endl; 370 G4cout << "G4QuarkExchange: Target mass 370 G4cout << "G4QuarkExchange: Target mass " << Ptarget.mag() << G4endl; 371 #endif 371 #endif 372 372 373 target->Set4Momentum(Ptarget); 373 target->Set4Momentum(Ptarget); 374 projectile->Set4Momentum(Pprojectile); 374 projectile->Set4Momentum(Pprojectile); 375 375 376 //=================================== Quark 376 //=================================== Quark exchange ================================ 377 projectile->SplitUp(); 377 projectile->SplitUp(); 378 target->SplitUp(); 378 target->SplitUp(); 379 379 380 G4Parton* PrQuark = nullptr; 380 G4Parton* PrQuark = nullptr; 381 G4Parton* TrQuark = nullptr; 381 G4Parton* TrQuark = nullptr; 382 382 383 if( projectile->GetDefinition()->GetBaryonNu 383 if( projectile->GetDefinition()->GetBaryonNumber() >= 0 ) { // Uzhi June 2020 384 // Quark exchange ---- 384 // Quark exchange ---- 385 PrQuark = projectile->GetNextParton(); 385 PrQuark = projectile->GetNextParton(); 386 TrQuark = target->GetNextParton(); 386 TrQuark = target->GetNextParton(); 387 G4ParticleDefinition * Tmp = PrQuark->GetD 387 G4ParticleDefinition * Tmp = PrQuark->GetDefinition(); 388 PrQuark->SetDefinition(TrQuark->GetDefinit 388 PrQuark->SetDefinition(TrQuark->GetDefinition()); 389 TrQuark->SetDefinition(Tmp); 389 TrQuark->SetDefinition(Tmp); 390 return true; 390 return true; 391 } 391 } 392 392 393 // Quark exchage for projectile anti-baryon 393 // Quark exchage for projectile anti-baryon (annihilation and new Q pair creation --- 394 // This part added by Uzhi June 2020 394 // This part added by Uzhi June 2020 395 PrQuark = projectile->GetNextAntiParton(); 395 PrQuark = projectile->GetNextAntiParton(); 396 TrQuark = target->GetNextParton(); 396 TrQuark = target->GetNextParton(); 397 if( -PrQuark->GetDefinition()->GetPDGEncodin 397 if( -PrQuark->GetDefinition()->GetPDGEncoding() == TrQuark->GetDefinition()->GetPDGEncoding() ){ 398 G4int QuarkCode = 1 + (int)(G4UniformRand( 398 G4int QuarkCode = 1 + (int)(G4UniformRand()/StrangeSuppress); 399 G4ParticleDefinition* AQpr = G4ParticleTab 399 G4ParticleDefinition* AQpr = G4ParticleTable::GetParticleTable()->FindParticle(-QuarkCode); 400 G4ParticleDefinition* Qtr = G4ParticleTab 400 G4ParticleDefinition* Qtr = G4ParticleTable::GetParticleTable()->FindParticle( QuarkCode); 401 if( (AQpr != nullptr) && (Qtr != nullptr) 401 if( (AQpr != nullptr) && (Qtr != nullptr) ) { 402 PrQuark->SetDefinition(AQpr); 402 PrQuark->SetDefinition(AQpr); 403 TrQuark->SetDefinition( Qtr); 403 TrQuark->SetDefinition( Qtr); 404 } 404 } 405 } 405 } 406 406 407 return true; 407 return true; 408 } 408 } 409 409 410 410 411 // --------- private methods ----------------- 411 // --------- private methods ---------------------- 412 412 413 G4ThreeVector G4QuarkExchange::GaussianPt(G4do 413 G4ThreeVector G4QuarkExchange::GaussianPt(G4double widthSquare, G4double maxPtSquare) const 414 { // @@ this method is used in FTF 414 { // @@ this method is used in FTFModel as well. Should go somewhere common! 415 415 416 G4double pt2; 416 G4double pt2; 417 417 418 const G4int maxNumberOfLoops = 1000; 418 const G4int maxNumberOfLoops = 1000; 419 G4int loopCounter = 0; 419 G4int loopCounter = 0; 420 do { 420 do { 421 pt2=-widthSquare * G4Log( G4UniformRand() 421 pt2=-widthSquare * G4Log( G4UniformRand() ); 422 } while ( ( pt2 > maxPtSquare) && ++loopCoun 422 } while ( ( pt2 > maxPtSquare) && ++loopCounter < maxNumberOfLoops ); /* Loop checking, 07.08.2015, A.Ribon */ 423 if ( loopCounter >= maxNumberOfLoops ) { 423 if ( loopCounter >= maxNumberOfLoops ) { 424 pt2 = 0.99*maxPtSquare; // Just an accept 424 pt2 = 0.99*maxPtSquare; // Just an acceptable value, without any physics consideration. 425 } 425 } 426 426 427 pt2=std::sqrt(pt2); 427 pt2=std::sqrt(pt2); 428 428 429 G4double phi=G4UniformRand() * twopi; 429 G4double phi=G4UniformRand() * twopi; 430 430 431 return G4ThreeVector (pt2*std::cos(phi), pt2 431 return G4ThreeVector (pt2*std::cos(phi), pt2*std::sin(phi), 0.); 432 } 432 } 433 433 434 434