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 // ---------------- G4QGSDiffractiveExcitation -------------- 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 // Essential changed by V. Uzhinsky in November - December 2006 36 // in order to put it in a correspondence with original FRITIOF 37 // model. Variant of FRITIOF with nucleon de-excitation is implemented. 38 // --------------------------------------------------------------------- 39 40 // Modified: 41 // 25-05-07 : G.Folger 42 // move from management/G4DiffractiveExcitation to to qgsm/G4QGSDiffractiveExcitation 43 // 44 45 #include "globals.hh" 46 #include "G4PhysicalConstants.hh" 47 #include "G4SystemOfUnits.hh" 48 #include "Randomize.hh" 49 50 #include "G4QGSDiffractiveExcitation.hh" 51 #include "G4LorentzRotation.hh" 52 #include "G4ThreeVector.hh" 53 #include "G4ParticleDefinition.hh" 54 #include "G4VSplitableHadron.hh" 55 #include "G4ExcitedString.hh" 56 //#include "G4ios.hh" 57 58 #include "G4Exp.hh" 59 #include "G4Log.hh" 60 #include "G4Pow.hh" 61 62 //============================================================================ 63 64 //#define debugDoubleDiffraction 65 66 //============================================================================ 67 68 G4QGSDiffractiveExcitation::G4QGSDiffractiveExcitation() 69 { 70 } 71 72 G4QGSDiffractiveExcitation::~G4QGSDiffractiveExcitation() 73 { 74 } 75 76 77 G4bool G4QGSDiffractiveExcitation:: 78 ExciteParticipants(G4VSplitableHadron *projectile, G4VSplitableHadron *target, G4bool ) const 79 { 80 #ifdef debugDoubleDiffraction 81 G4cout<<G4endl<<"G4QGSDiffractiveExcitation::ExciteParticipants - Double diffraction."<<G4endl; 82 G4cout<<"Proj Targ "<<projectile->GetDefinition()->GetParticleName()<<" "<<target->GetDefinition()->GetParticleName()<<G4endl; 83 G4cout<<"Proj 4 Mom "<<projectile->Get4Momentum()<<" "<<projectile->Get4Momentum().mag()<<G4endl; 84 G4cout<<"Targ 4 Mom "<<target->Get4Momentum() <<" "<<target->Get4Momentum().mag() <<G4endl; 85 #endif 86 87 G4LorentzVector Pprojectile=projectile->Get4Momentum(); 88 89 // -------------------- Projectile parameters ----------------------------------- 90 G4bool PutOnMassShell=0; 91 92 G4double M0projectile = Pprojectile.mag(); // Without de-excitation 93 94 if(M0projectile < projectile->GetDefinition()->GetPDGMass()) 95 { 96 PutOnMassShell=1; 97 M0projectile=projectile->GetDefinition()->GetPDGMass(); 98 } 99 100 // -------------------- Target parameters ---------------------------------------------- 101 G4LorentzVector Ptarget=target->Get4Momentum(); 102 103 G4double M0target = Ptarget.mag(); 104 105 if(M0target < target->GetDefinition()->GetPDGMass()) 106 { 107 PutOnMassShell=1; 108 M0target=target->GetDefinition()->GetPDGMass(); 109 } 110 111 G4LorentzVector Psum=Pprojectile+Ptarget; 112 G4double S=Psum.mag2(); 113 G4double SqrtS=std::sqrt(S); 114 115 if (SqrtS < M0projectile + M0target) {return false;} // The model cannot work for pp-interactions 116 // at Plab < 1.3 GeV/c. 117 118 G4double Mprojectile2 = M0projectile * M0projectile; 119 G4double Mtarget2 = M0target * M0target; //Ptarget.mag2(); // for AA-inter. 120 121 // Transform momenta to cms and then rotate parallel to z axis; 122 123 G4LorentzRotation toCms(-1*Psum.boostVector()); 124 125 G4LorentzVector Ptmp=toCms*Pprojectile; 126 127 if ( Ptmp.pz() <= 0. ) 128 { 129 // "String" moving backwards in CMS, abort collision !! 130 //G4cout << " abort Collision!! " << G4endl; 131 return false; 132 } 133 134 toCms.rotateZ(-1*Ptmp.phi()); 135 toCms.rotateY(-1*Ptmp.theta()); 136 137 G4LorentzRotation toLab(toCms.inverse()); 138 139 Pprojectile.transform(toCms); 140 Ptarget.transform(toCms); 141 142 G4double PZcms2=(S*S+Mprojectile2*Mprojectile2+Mtarget2*Mtarget2- 143 2*S*Mprojectile2-2*S*Mtarget2-2*Mprojectile2*Mtarget2)/4./S; 144 145 if (PZcms2 < 0) {return false;} // It can be in an interaction with off-shell nuclear nucleon 146 147 G4double PZcms = std::sqrt(PZcms2); 148 149 if (PutOnMassShell) 150 { 151 if (Pprojectile.z() > 0.) 152 { 153 Pprojectile.setPz( PZcms); 154 Ptarget.setPz( -PZcms); 155 } 156 else 157 { 158 Pprojectile.setPz(-PZcms); 159 Ptarget.setPz( PZcms); 160 }; 161 162 Pprojectile.setE(std::sqrt(Mprojectile2+sqr(Pprojectile.x())+sqr(Pprojectile.y())+PZcms2)); 163 Ptarget.setE( std::sqrt( Mtarget2+sqr( Ptarget.x())+sqr( Ptarget.y())+PZcms2)); 164 } 165 166 G4double maxPtSquare = PZcms2; 167 168 #ifdef debugDoubleDiffraction 169 G4cout << "Pprojectile after boost to CMS: " << Pprojectile <<" "<<Pprojectile.mag()<<G4endl; 170 G4cout << "Ptarget after boost to CMS: " << Ptarget <<" "<<Ptarget.mag() <<G4endl; 171 #endif 172 173 G4int PrPDGcode=projectile->GetDefinition()->GetPDGEncoding(); 174 G4int absPrPDGcode=std::abs(PrPDGcode); 175 G4double MinPrDiffMass(0.); 176 G4double AveragePt2(0.); 177 178 if (M0projectile <= projectile->GetDefinition()->GetPDGMass()) 179 { // Normal projectile 180 if( absPrPDGcode > 1000 ) //------Projectile is baryon -------- 181 { 182 if ( absPrPDGcode > 4000 && absPrPDGcode < 6000 ) // Projectile is a charm or bottom baryon 183 { 184 MinPrDiffMass = projectile->GetDefinition()->GetPDGMass()/CLHEP::GeV + 0.25; // GeV 185 AveragePt2 = 0.3; // GeV^2 186 } 187 else 188 { 189 MinPrDiffMass = 1.16; // GeV 190 AveragePt2 = 0.3; // GeV^2 191 } 192 } 193 else if( absPrPDGcode == 211 || PrPDGcode == 111) //------Projectile is Pion ----------- 194 { 195 MinPrDiffMass = 1.0; // GeV 196 AveragePt2 = 0.3; // GeV^2 197 } 198 else if( absPrPDGcode == 321 || absPrPDGcode == 130 || absPrPDGcode == 310) //-Projectile is Kaon- 199 { 200 MinPrDiffMass = 1.1; // GeV 201 AveragePt2 = 0.3; // GeV^2 202 } 203 else if( absPrPDGcode > 400 && absPrPDGcode < 600) // Projectile is a charm or bottom meson 204 { 205 MinPrDiffMass = projectile->GetDefinition()->GetPDGMass()/CLHEP::GeV + 0.25; // GeV 206 AveragePt2 = 0.3; // GeV^2 207 } 208 else //------Projectile is undefined, Nucleon assumed 209 { 210 MinPrDiffMass = 1.16; // GeV 211 AveragePt2 = 0.3; // GeV^2 212 } 213 } 214 else 215 { // Excited projectile 216 MinPrDiffMass = M0projectile + 220.0*MeV; 217 AveragePt2 = 0.3; 218 } 219 220 MinPrDiffMass = MinPrDiffMass * GeV; 221 AveragePt2 = AveragePt2 * GeV*GeV; 222 //--------------------------------------------- 223 G4double MinTrDiffMass = 1.16*GeV; 224 225 if (SqrtS < MinPrDiffMass + MinTrDiffMass) {return false;} // The model cannot work at low energy 226 227 G4double MinPrDiffMass2 = MinPrDiffMass * MinPrDiffMass; 228 G4double MinTrDiffMass2 = MinTrDiffMass * MinTrDiffMass; 229 230 G4double Pt2; 231 G4double ProjMassT2, ProjMassT; 232 G4double TargMassT2, TargMassT; 233 G4double PMinusNew, TPlusNew; 234 235 G4LorentzVector Qmomentum; 236 G4double Qminus, Qplus; 237 238 G4int whilecount=0; 239 do { 240 if (whilecount++ >= 500 && (whilecount%100)==0) 241 if (whilecount > 1000 ) { 242 Qmomentum=G4LorentzVector(0.,0.,0.,0.); 243 return false; // Ignore this interaction 244 } 245 246 // Generate pt 247 Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0); 248 249 Pt2=G4ThreeVector(Qmomentum.vect()).mag2(); 250 ProjMassT2=MinPrDiffMass2+Pt2; 251 ProjMassT =std::sqrt(ProjMassT2); 252 253 TargMassT2=MinTrDiffMass2+Pt2; 254 TargMassT =std::sqrt(TargMassT2); 255 256 if (SqrtS < ProjMassT + TargMassT) continue; 257 258 PZcms2=(S*S+ProjMassT2*ProjMassT2+TargMassT2*TargMassT2- 259 2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)/4./S; 260 if (PZcms2 < 0 ) {PZcms2=0;}; 261 PZcms =std::sqrt(PZcms2); 262 263 G4double PMinusMin=std::sqrt(ProjMassT2+PZcms2)-PZcms; 264 G4double PMinusMax=SqrtS-TargMassT; 265 266 PMinusNew=ChooseP(PMinusMin,PMinusMax); 267 Qminus=PMinusNew-Pprojectile.minus(); 268 269 G4double TPlusMin=std::sqrt(TargMassT2+PZcms2)-PZcms; 270 G4double TPlusMax=SqrtS-ProjMassT; 271 272 TPlusNew=ChooseP(TPlusMin, TPlusMax); 273 Qplus=-(TPlusNew-Ptarget.plus()); 274 275 Qmomentum.setPz( (Qplus-Qminus)/2 ); 276 Qmomentum.setE( (Qplus+Qminus)/2 ); 277 278 } while ( (Pprojectile+Qmomentum).mag2() < MinPrDiffMass2 || 279 (Ptarget -Qmomentum).mag2() < MinTrDiffMass2 ); 280 281 Pprojectile += Qmomentum; 282 Ptarget -= Qmomentum; 283 284 // Transform back and update SplitableHadron Participant. 285 Pprojectile.transform(toLab); 286 Ptarget.transform(toLab); 287 288 #ifdef debugDoubleDiffraction 289 G4cout << "Pprojectile after boost to Lab: " << Pprojectile <<" "<<Pprojectile.mag()<<G4endl; 290 G4cout << "Ptarget after boost to Lab: " << Ptarget <<" "<<Ptarget.mag() <<G4endl; 291 #endif 292 293 target->Set4Momentum(Ptarget); 294 projectile->Set4Momentum(Pprojectile); 295 296 return true; 297 } 298 299 300 G4ExcitedString * G4QGSDiffractiveExcitation:: 301 String(G4VSplitableHadron * hadron, G4bool isProjectile) const 302 { 303 hadron->SplitUp(); 304 G4Parton *start= hadron->GetNextParton(); 305 if ( start==NULL) 306 { G4cout << " G4QGSDiffractiveExcitation::String() Error:No start parton found"<< G4endl; 307 return NULL; 308 } 309 G4Parton *end = hadron->GetNextParton(); 310 if ( end==NULL) 311 { G4cout << " G4QGSDiffractiveExcitation::String() Error:No end parton found"<< G4endl; 312 return NULL; 313 } 314 315 G4ExcitedString * string; 316 if ( isProjectile ) 317 { 318 string= new G4ExcitedString(end,start, +1); 319 } else { 320 string= new G4ExcitedString(start,end, -1); 321 } 322 323 string->SetPosition(hadron->GetPosition()); 324 325 // momenta of string ends 326 327 G4double maxAvailMomentumSquared=sqr(hadron->Get4Momentum().mag()/2.); 328 329 G4double widthOfPtSquare = 0.5*sqr(GeV); 330 G4ThreeVector pt=GaussianPt(widthOfPtSquare,maxAvailMomentumSquared); 331 332 G4LorentzVector Pstart(G4LorentzVector(pt,0.)); 333 G4LorentzVector Pend; 334 Pend.setPx(hadron->Get4Momentum().px() - pt.x()); 335 Pend.setPy(hadron->Get4Momentum().py() - pt.y()); 336 337 G4double tm1=hadron->Get4Momentum().minus() + 338 ( Pend.perp2()-Pstart.perp2() ) / hadron->Get4Momentum().plus(); 339 340 G4double tm2= std::sqrt( std::max(0., sqr(tm1) - 341 4. * Pend.perp2() * hadron->Get4Momentum().minus() 342 / hadron->Get4Momentum().plus() )); 343 344 G4int Sign= isProjectile ? -1 : 1; 345 346 G4double endMinus = 0.5 * (tm1 + Sign*tm2); 347 G4double startMinus= hadron->Get4Momentum().minus() - endMinus; 348 349 G4double startPlus= Pstart.perp2() / startMinus; 350 G4double endPlus = hadron->Get4Momentum().plus() - startPlus; 351 352 Pstart.setPz(0.5*(startPlus - startMinus)); 353 Pstart.setE(0.5*(startPlus + startMinus)); 354 355 Pend.setPz(0.5*(endPlus - endMinus)); 356 Pend.setE(0.5*(endPlus + endMinus)); 357 358 start->Set4Momentum(Pstart); 359 end->Set4Momentum(Pend); 360 361 #ifdef debugQGSdiffExictation 362 G4cout << " generated string flavors " << start->GetPDGcode() << " / " << end->GetPDGcode() << G4endl; 363 G4cout << " generated string momenta: quark " << start->Get4Momentum() << "mass : " <<start->Get4Momentum().mag()<< G4endl; 364 G4cout << " generated string momenta: Diquark " << end ->Get4Momentum() << "mass : " <<end->Get4Momentum().mag()<< G4endl; 365 G4cout << " sum of ends " << Pstart+Pend << G4endl; 366 G4cout << " Original " << hadron->Get4Momentum() << G4endl; 367 #endif 368 369 return string; 370 } 371 372 373 // --------- private methods ---------------------- 374 375 G4double G4QGSDiffractiveExcitation::ChooseP(G4double Pmin, G4double Pmax) const 376 { 377 // choose an x between Xmin and Xmax with P(x) ~ 1/x 378 // to be improved... 379 380 G4double range=Pmax-Pmin; 381 382 if ( Pmin <= 0. || range <=0. ) 383 { 384 G4cout << " Pmin, range : " << Pmin << " , " << range << G4endl; 385 throw G4HadronicException(__FILE__, __LINE__, "G4QGSDiffractiveExcitation::ChooseP : Invalid arguments "); 386 } 387 388 G4double P; 389 P=Pmin * G4Pow::GetInstance()->powA(Pmax/Pmin,G4UniformRand()); 390 //debug-hpw cout << "DiffractiveX "<<x<<G4endl; 391 return P; 392 } 393 394 395 G4ThreeVector G4QGSDiffractiveExcitation::GaussianPt(G4double AveragePt2, G4double maxPtSquare) const 396 { // @@ this method is used in FTFModel as well. Should go somewhere common! 397 G4double Pt2; 398 399 Pt2 = -AveragePt2 * G4Log(1. + G4UniformRand() * (G4Exp(-maxPtSquare/AveragePt2)-1.)); 400 401 G4double Pt=std::sqrt(Pt2); 402 403 G4double phi=G4UniformRand() * twopi; 404 405 return G4ThreeVector (Pt*std::cos(phi), Pt*std::sin(phi), 0.); 406 } 407