Geant4 Cross Reference |
1 // 2 // ******************************************************************** 3 // * License and Disclaimer * 4 // * * 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. * 10 // * * 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitation of liability. * 17 // * * 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************************************** 25 // 26 // 27 // ------------------------------------------------------------ 28 // GEANT 4 class implemetation file 29 // 30 // ---------------- G4SingleDiffractiveExcitation -------------- 31 // by Gunter Folger, October 1998. 32 // diffractive Excitation used by strings models 33 // Take a projectile and a target 34 // excite the projectile and target 35 // ------------------------------------------------------------ 36 37 #include "G4SingleDiffractiveExcitation.hh" 38 #include "globals.hh" 39 #include "G4PhysicalConstants.hh" 40 #include "G4SystemOfUnits.hh" 41 #include "Randomize.hh" 42 #include "G4LorentzRotation.hh" 43 #include "G4ThreeVector.hh" 44 #include "G4ParticleDefinition.hh" 45 #include "G4VSplitableHadron.hh" 46 #include "G4ExcitedString.hh" 47 48 #include "G4Log.hh" 49 #include "G4Pow.hh" 50 51 //#define debugSingleDiffraction 52 53 G4SingleDiffractiveExcitation::G4SingleDiffractiveExcitation(){} 54 55 G4SingleDiffractiveExcitation::~G4SingleDiffractiveExcitation(){} 56 57 G4bool G4SingleDiffractiveExcitation:: 58 ExciteParticipants( G4VSplitableHadron *projectile, G4VSplitableHadron *target, 59 G4bool ProjectileDiffraction ) const 60 { 61 #ifdef debugSingleDiffraction 62 G4cout<<G4endl<<"G4SingleDiffractiveExcitation::ExciteParticipants"<<G4endl; 63 #endif 64 65 G4LorentzVector Pprojectile=projectile->Get4Momentum(); 66 G4double Mprojectile = projectile->GetDefinition()->GetPDGMass(); 67 G4double Mprojectile2=sqr(projectile->GetDefinition()->GetPDGMass()); 68 69 G4LorentzVector Ptarget=target->Get4Momentum(); 70 G4double Mtarget = target->GetDefinition()->GetPDGMass(); 71 G4double Mtarget2=sqr(target->GetDefinition()->GetPDGMass()); 72 73 #ifdef debugSingleDiffraction 74 G4cout<<"Proj Targ "<<projectile->GetDefinition()->GetPDGEncoding()<<" "<<target->GetDefinition()->GetPDGEncoding()<<G4endl; 75 G4cout<<"Pr Tr 4-Mom "<<Pprojectile<<" "<<Pprojectile.mag()<<G4endl 76 <<" "<<Ptarget <<" "<<Ptarget.mag() <<G4endl; 77 #endif 78 79 G4LorentzVector Psum=Pprojectile+Ptarget; 80 G4double SqrtS=Psum.mag(); 81 G4double S =Psum.mag2(); 82 83 #ifdef debugSingleDiffraction 84 G4cout<<"SqrtS-Mprojectile-Mtarget "<<SqrtS<<" "<<Mprojectile<<" "<<Mtarget 85 <<" "<<SqrtS-Mprojectile-Mtarget<<G4endl; 86 #endif 87 if (SqrtS-Mprojectile-Mtarget <= 250.0*MeV) { 88 #ifdef debugSingleDiffraction 89 G4cerr<<"Projectile: "<<projectile->GetDefinition()->GetPDGEncoding()<<" " 90 <<Pprojectile<<" "<<Pprojectile.mag()<<G4endl; 91 G4cerr<<"Target: "<<target->GetDefinition()->GetPDGEncoding()<<" " 92 <<Ptarget<<" "<<Ptarget.mag()<<G4endl; 93 G4cerr<<"sqrt(S) = "<<SqrtS<<" Mp + Mt = "<<Pprojectile.mag()+Ptarget.mag()<<G4endl; 94 #endif 95 return true; 96 } 97 98 G4LorentzRotation toCms(-1*Psum.boostVector()); 99 100 G4LorentzVector Ptmp=toCms*Pprojectile; 101 102 if ( Ptmp.pz() <= 0. ) 103 { 104 // "String" moving backwards in CMS, abort collision !! 105 // G4cout << " abort Collision!! " << G4endl; 106 return false; 107 } 108 109 toCms.rotateZ(-1*Ptmp.phi()); 110 toCms.rotateY(-1*Ptmp.theta()); 111 112 G4LorentzRotation toLab(toCms.inverse()); 113 114 Pprojectile.transform(toCms); 115 Ptarget.transform(toCms); 116 #ifdef debugSingleDiffraction 117 G4cout << "Pprojectile in CMS " << Pprojectile << G4endl; 118 G4cout << "Ptarget in CMS " << Ptarget << G4endl; 119 #endif 120 G4double maxPtSquare=sqr(Ptarget.pz()); 121 122 G4double ProjectileMinDiffrMass(0.), TargetMinDiffrMass(0.); 123 G4double AveragePt2(0.); 124 G4int absPDGcode=std::abs(projectile->GetDefinition()->GetPDGEncoding()); 125 126 if ( ProjectileDiffraction ) { 127 if ( absPDGcode > 1000 ) //------Projectile is baryon -------- 128 { 129 if ( absPDGcode > 4000 && absPDGcode < 6000 ) // Projectile is a charm or bottom baryon 130 { 131 ProjectileMinDiffrMass = projectile->GetDefinition()->GetPDGMass()/CLHEP::GeV + 0.25; // GeV 132 AveragePt2 = 0.3; // GeV^2 133 } 134 else 135 { 136 ProjectileMinDiffrMass = 1.16; // GeV 137 AveragePt2 = 0.3; // GeV^2 138 } 139 } 140 else if( absPDGcode == 211 || absPDGcode == 111) //------Projectile is Pion ----------- 141 { 142 ProjectileMinDiffrMass = 1.0; // GeV 143 AveragePt2 = 0.3; // GeV^2 144 } 145 else if( absPDGcode == 321 || absPDGcode == 130 || absPDGcode == 310) //Projectile is Kaon 146 { 147 ProjectileMinDiffrMass = 1.1; // GeV 148 AveragePt2 = 0.3; // GeV^2 149 } 150 else if( absPDGcode == 22) //------Projectile is Gamma ----------- 151 { 152 ProjectileMinDiffrMass = 0.25; // GeV 153 AveragePt2 = 0.36; // GeV^2 154 } 155 else if( absPDGcode > 400 && absPDGcode < 600) // Projectile is a charm or bottom meson 156 { 157 ProjectileMinDiffrMass = projectile->GetDefinition()->GetPDGMass()/CLHEP::GeV + 0.25; // GeV 158 AveragePt2 = 0.3; // GeV^2 159 } 160 else //------Projectile is undefined, Nucleon assumed 161 { 162 ProjectileMinDiffrMass = 1.1; // GeV 163 AveragePt2 = 0.3; // GeV^2 164 }; 165 166 ProjectileMinDiffrMass = ProjectileMinDiffrMass * GeV; 167 Mprojectile2=sqr(ProjectileMinDiffrMass); 168 } 169 else 170 { 171 TargetMinDiffrMass = 1.16*GeV; // For target nucleon 172 Mtarget2 = sqr( TargetMinDiffrMass) ; 173 AveragePt2 = 0.3; // GeV^2 174 } // end of if ( ProjectileDiffraction ) 175 176 AveragePt2 = AveragePt2 * GeV*GeV; 177 178 G4double Pt2, PZcms, PZcms2; 179 G4double ProjMassT2, ProjMassT; 180 G4double TargMassT2, TargMassT; 181 G4double PMinusMin, PMinusMax; 182 G4double TPlusMin, TPlusMax; 183 G4double PMinusNew, PPlusNew, TPlusNew, TMinusNew; 184 185 G4LorentzVector Qmomentum; 186 G4double Qminus, Qplus; 187 188 G4int whilecount=0; 189 do { 190 whilecount++; 191 192 if (whilecount > 1000 ) 193 { 194 Qmomentum=G4LorentzVector(0.,0.,0.,0.); 195 return false; // Ignore this interaction 196 } 197 198 // Generate pt 199 Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0); 200 201 Pt2 = G4ThreeVector( Qmomentum.vect() ).mag2(); 202 203 ProjMassT2 = Mprojectile2 + Pt2; 204 ProjMassT = std::sqrt( ProjMassT2 ); 205 TargMassT2 = Mtarget2 + Pt2; 206 TargMassT = std::sqrt( TargMassT2 ); 207 208 #ifdef debugSingleDiffraction 209 G4cout<<whilecount<<" "<<Pt2<<" "<<ProjMassT<<" "<<TargMassT<<" "<<SqrtS<<" "<<S<<" "<<ProjectileDiffraction<<G4endl; 210 #endif 211 if ( SqrtS < ProjMassT + TargMassT ) continue; 212 213 PZcms2 = ( S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2 214 - 2.0*S*ProjMassT2 - 2.0*S*TargMassT2 - 2.0*ProjMassT2*TargMassT2 ) / 4.0 / S; 215 216 if ( PZcms2 < 0 ) continue; 217 218 PZcms = std::sqrt( PZcms2 ); 219 220 if ( ProjectileDiffraction ) 221 { // The projectile will fragment, the target will saved. 222 PMinusMin = std::sqrt( ProjMassT2 + PZcms2 ) - PZcms; 223 PMinusMax = SqrtS - TargMassT; 224 225 PMinusNew = ChooseX( PMinusMin, PMinusMax ); 226 TMinusNew = SqrtS - PMinusNew; 227 228 Qminus = Ptarget.minus() - TMinusNew; 229 TPlusNew = TargMassT2 / TMinusNew; 230 Qplus = Ptarget.plus() - TPlusNew; 231 232 } else { // The target will fragment, the projectile will saved. 233 TPlusMin = std::sqrt( TargMassT2 + PZcms2 ) - PZcms; 234 TPlusMax = SqrtS - ProjMassT; 235 236 TPlusNew = ChooseX( TPlusMin, TPlusMax ); 237 PPlusNew = SqrtS - TPlusNew; 238 239 Qplus = PPlusNew - Pprojectile.plus(); 240 PMinusNew = ProjMassT2 / PPlusNew; 241 Qminus = PMinusNew - Pprojectile.minus(); 242 } 243 244 Qmomentum.setPz( (Qplus - Qminus)/2 ); 245 Qmomentum.setE( (Qplus + Qminus)/2 ); 246 247 #ifdef debugSingleDiffraction 248 G4cout<<ProjectileDiffraction<<" "<<( Pprojectile + Qmomentum ).mag2()<<" "<< Mprojectile2<<G4endl; 249 G4cout<<!ProjectileDiffraction<<" "<<( Ptarget - Qmomentum ).mag2()<<" "<< Mtarget2<<G4endl; 250 #endif 251 252 } while ( ( ProjectileDiffraction&&( Pprojectile + Qmomentum ).mag2() < Mprojectile2 ) || 253 (!ProjectileDiffraction&&( Ptarget - Qmomentum ).mag2() < Mtarget2 ) ); 254 // Repeat the sampling because there was not any excitation 255 256 Pprojectile += Qmomentum; 257 258 Ptarget -= Qmomentum; 259 260 // Transform back and update SplitableHadron Participant. 261 Pprojectile.transform(toLab); 262 Ptarget.transform(toLab); 263 264 #ifdef debugSingleDiffraction 265 G4cout << "Pprojectile in Lab. " << Pprojectile << G4endl; 266 G4cout << "Ptarget in Lab. " << Ptarget << G4endl; 267 G4cout << "G4SingleDiffractiveExcitation- Projectile mass " << Pprojectile.mag() << G4endl; 268 G4cout << "G4SingleDiffractiveExcitation- Target mass " << Ptarget.mag() << G4endl; 269 #endif 270 271 target->Set4Momentum(Ptarget); 272 projectile->Set4Momentum(Pprojectile); 273 274 return true; 275 } 276 277 // --------- private methods ---------------------- 278 279 G4double G4SingleDiffractiveExcitation::ChooseX(G4double Xmin, G4double Xmax) const 280 { 281 // choose an x between Xmin and Xmax with P(x) ~ 1/x 282 G4double range=Xmax-Xmin; 283 284 if ( Xmin <= 0. || range <=0. ) 285 { 286 G4cout << " Xmin, range : " << Xmin << " , " << range << G4endl; 287 throw G4HadronicException(__FILE__, __LINE__, "G4SingleDiffractiveExcitation::ChooseX : Invalid arguments "); 288 } 289 290 G4double x = Xmin*G4Pow::GetInstance()->powA(Xmax/Xmin, G4UniformRand() ); 291 // G4double x = 1.0/sqr(1.0/std::sqrt(Xmin) - G4UniformRand() * (1.0/std::sqrt(Xmin) - 1.0/std::sqrt(Xmax))); 292 return x; 293 } 294 295 296 G4ThreeVector G4SingleDiffractiveExcitation::GaussianPt(G4double widthSquare, G4double maxPtSquare) const 297 { // @@ this method is used in FTFModel as well. Should go somewhere common! 298 299 G4double pt2; 300 301 const G4int maxNumberOfLoops = 1000; 302 G4int loopCounter = 0; 303 do { 304 pt2=-widthSquare * G4Log( G4UniformRand() ); 305 } while ( ( pt2 > maxPtSquare) && ++loopCounter < maxNumberOfLoops ); /* Loop checking, 07.08.2015, A.Ribon */ 306 if ( loopCounter >= maxNumberOfLoops ) { 307 pt2 = 0.99*maxPtSquare; // Just an acceptable value, without any physics consideration. 308 } 309 310 pt2=std::sqrt(pt2); 311 312 G4double phi=G4UniformRand() * twopi; 313 314 return G4ThreeVector (pt2*std::cos(phi), pt2*std::sin(phi), 0.); 315 } 316 317