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 // $Id: G4SingleDiffractiveExcitation.cc,v 1.1 2007/05/25 07:30:47 gunter Exp $ 27 // ------------------------------------------- 28 // ------------------------------------------------------------ 28 // GEANT 4 class implemetation file 29 // GEANT 4 class implemetation file 29 // 30 // 30 // ---------------- G4SingleDiffractiveEx 31 // ---------------- G4SingleDiffractiveExcitation -------------- 31 // by Gunter Folger, October 1998. 32 // by Gunter Folger, October 1998. 32 // diffractive Excitation used by strings 33 // diffractive Excitation used by strings models 33 // Take a projectile and a target 34 // Take a projectile and a target 34 // excite the projectile and target 35 // excite the projectile and target 35 // ------------------------------------------- 36 // ------------------------------------------------------------ 36 37 37 #include "G4SingleDiffractiveExcitation.hh" << 38 38 #include "globals.hh" 39 #include "globals.hh" 39 #include "G4PhysicalConstants.hh" << 40 #include "G4SystemOfUnits.hh" << 41 #include "Randomize.hh" 40 #include "Randomize.hh" >> 41 >> 42 #include "G4SingleDiffractiveExcitation.hh" 42 #include "G4LorentzRotation.hh" 43 #include "G4LorentzRotation.hh" 43 #include "G4ThreeVector.hh" 44 #include "G4ThreeVector.hh" 44 #include "G4ParticleDefinition.hh" 45 #include "G4ParticleDefinition.hh" 45 #include "G4VSplitableHadron.hh" 46 #include "G4VSplitableHadron.hh" 46 #include "G4ExcitedString.hh" 47 #include "G4ExcitedString.hh" >> 48 //#include "G4ios.hh" 47 49 48 #include "G4Log.hh" << 50 G4SingleDiffractiveExcitation::G4SingleDiffractiveExcitation(G4double sigmaPt, G4double minextraMass,G4double x0mass) 49 #include "G4Pow.hh" << 51 : >> 52 widthOfPtSquare(-2*sqr(sigmaPt)) , minExtraMass(minextraMass), >> 53 minmass(x0mass) >> 54 { >> 55 } 50 56 51 //#define debugSingleDiffraction << 57 G4bool G4SingleDiffractiveExcitation:: >> 58 ExciteParticipants(G4VSplitableHadron *projectile, G4VSplitableHadron *target) const >> 59 { 52 60 53 G4SingleDiffractiveExcitation::G4SingleDiffrac << 61 G4LorentzVector Pprojectile=projectile->Get4Momentum(); >> 62 G4double Mprojectile2=sqr(projectile->GetDefinition()->GetPDGMass() + minExtraMass); 54 63 55 G4SingleDiffractiveExcitation::~G4SingleDiffra << 64 G4LorentzVector Ptarget=target->Get4Momentum(); >> 65 G4double Mtarget2=sqr(target->GetDefinition()->GetPDGMass() + minExtraMass); >> 66 // G4cout << "E proj, target :" << Pprojectile.e() << ", " << >> 67 // Ptarget.e() << G4endl; >> 68 >> 69 G4bool KeepProjectile= G4UniformRand() > 0.5; >> 70 >> 71 // reset the min.mass of the non diffractive particle to its value, ( minus a bit for rounding...) >> 72 if ( KeepProjectile ) >> 73 { >> 74 // cout << " Projectile fix" << G4endl; >> 75 Mprojectile2 = sqr(projectile->GetDefinition()->GetPDGMass() * (1-perCent) ); >> 76 } else { >> 77 // cout << " Target fix" << G4endl; >> 78 Mtarget2=sqr(target->GetDefinition()->GetPDGMass() * (1-perCent) ); >> 79 } >> 80 >> 81 // Transform momenta to cms and then rotate parallel to z axis; >> 82 >> 83 G4LorentzVector Psum; >> 84 Psum=Pprojectile+Ptarget; >> 85 >> 86 G4LorentzRotation toCms(-1*Psum.boostVector()); >> 87 >> 88 G4LorentzVector Ptmp=toCms*Pprojectile; >> 89 >> 90 if ( Ptmp.pz() <= 0. ) >> 91 { >> 92 // "String" moving backwards in CMS, abort collision !! >> 93 // G4cout << " abort Collision!! " << G4endl; >> 94 return false; >> 95 } >> 96 >> 97 toCms.rotateZ(-1*Ptmp.phi()); >> 98 toCms.rotateY(-1*Ptmp.theta()); 56 99 57 G4bool G4SingleDiffractiveExcitation:: << 100 // G4cout << "Pprojectile be4 boost " << Pprojectile << G4endl; 58 ExciteParticipants( G4VSplitableHadron *projec << 101 // G4cout << "Ptarget be4 boost : " << Ptarget << G4endl; 59 G4bool ProjectileDiffracti << 60 { << 61 #ifdef debugSingleDiffraction << 62 G4cout<<G4endl<<"G4SingleDiffractiveExcita << 63 #endif << 64 << 65 G4LorentzVector Pprojectile=projectile->Get4 << 66 G4double Mprojectile = projectile->GetDef << 67 G4double Mprojectile2=sqr(projectile->GetDef << 68 << 69 G4LorentzVector Ptarget=target->Get4Momentum << 70 G4double Mtarget = target->GetDefinition( << 71 G4double Mtarget2=sqr(target->GetDefinition( << 72 << 73 #ifdef debugSingleDiffraction << 74 G4cout<<"Proj Targ "<<projectile->GetDefin << 75 G4cout<<"Pr Tr 4-Mom "<<Pprojectile<<" "<< << 76 <<" "<<Ptarget <<" "<< << 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 "<<Sqrt << 85 <<" "<<SqrtS-Mprojectile-Mtarget<<G4 << 86 #endif << 87 if (SqrtS-Mprojectile-Mtarget <= 250.0*MeV) << 88 #ifdef debugSingleDiffraction << 89 G4cerr<<"Projectile: "<<projectile->GetDef << 90 <<Pprojectile<<" "<<Pprojectile.mag( << 91 G4cerr<<"Target: "<<target->GetDefinit << 92 <<Ptarget<<" "<<Ptarget.mag()<<G4end << 93 G4cerr<<"sqrt(S) = "<<SqrtS<<" Mp + Mt = " << 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, abor << 105 // G4cout << " abort Collision!! " << 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 " << Pproje << 118 G4cout << "Ptarget in CMS " << Ptarge << 119 #endif << 120 G4double maxPtSquare=sqr(Ptarget.pz()); << 121 << 122 G4double ProjectileMinDiffrMass(0.), TargetM << 123 G4double AveragePt2(0.); << 124 G4int absPDGcode=std::abs(projectile->GetDef << 125 << 126 if ( ProjectileDiffraction ) { << 127 if ( absPDGcode > 1000 ) << 128 { << 129 if ( absPDGcode > 4000 && absPDGcode < 6 << 130 { << 131 ProjectileMinDiffrMass = projectile->G << 132 AveragePt2 = 0.3; << 133 } << 134 else << 135 { << 136 ProjectileMinDiffrMass = 1.16; << 137 AveragePt2 = 0.3; << 138 } << 139 } << 140 else if( absPDGcode == 211 || absPDGcode = << 141 { << 142 ProjectileMinDiffrMass = 1.0; << 143 AveragePt2 = 0.3; << 144 } << 145 else if( absPDGcode == 321 || absPDGcode = << 146 { << 147 ProjectileMinDiffrMass = 1.1; << 148 AveragePt2 = 0.3; << 149 } << 150 else if( absPDGcode == 22) << 151 { << 152 ProjectileMinDiffrMass = 0.25; << 153 AveragePt2 = 0.36; << 154 } << 155 else if( absPDGcode > 400 && absPDGcode < << 156 { << 157 ProjectileMinDiffrMass = projectile->Get << 158 AveragePt2 = 0.3; << 159 } << 160 else << 161 { << 162 ProjectileMinDiffrMass = 1.1; << 163 AveragePt2 = 0.3; << 164 }; << 165 << 166 ProjectileMinDiffrMass = ProjectileMinDiff << 167 Mprojectile2=sqr(ProjectileMinDiffrMass); << 168 } << 169 else << 170 { << 171 TargetMinDiffrMass = 1.16*GeV; << 172 Mtarget2 = sqr( TargetMinDiffrMass) ; << 173 AveragePt2 = 0.3; << 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, TMin << 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 intera << 196 } << 197 << 198 // Generate pt << 199 Qmomentum=G4LorentzVector(GaussianPt(Avera << 200 << 201 Pt2 = G4ThreeVector( Qmomentum.vect() ).ma << 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<<" "<<ProjM << 210 #endif << 211 if ( SqrtS < ProjMassT + TargMassT ) conti << 212 << 213 PZcms2 = ( S*S + ProjMassT2*ProjMassT2 + T << 214 - 2.0*S*ProjMassT2 - 2.0*S*TargM << 215 << 216 if ( PZcms2 < 0 ) continue; << 217 << 218 PZcms = std::sqrt( PZcms2 ); << 219 << 220 if ( ProjectileDiffraction ) << 221 { // The projectile will fragment, t << 222 PMinusMin = std::sqrt( ProjMassT2 + PZcm << 223 PMinusMax = SqrtS - TargMassT; << 224 << 225 PMinusNew = ChooseX( PMinusMin, PMinusMa << 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 << 233 TPlusMin = std::sqrt( TargMassT2 + PZcms << 234 TPlusMax = SqrtS - ProjMassT; << 235 102 236 TPlusNew = ChooseX( TPlusMin, TPlusMax ) << 237 PPlusNew = SqrtS - TPlusNew; << 238 103 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<<" "<<( Pp << 249 G4cout<<!ProjectileDiffraction<<" "<<( P << 250 #endif << 251 << 252 } while ( ( ProjectileDiffraction&&( Pprojec << 253 (!ProjectileDiffraction&&( Ptarget << 254 // Repeat the sampling because there was n << 255 << 256 Pprojectile += Qmomentum; << 257 << 258 Ptarget -= Qmomentum; << 259 << 260 // Transform back and update SplitableHadron << 261 Pprojectile.transform(toLab); << 262 Ptarget.transform(toLab); << 263 << 264 #ifdef debugSingleDiffraction << 265 G4cout << "Pprojectile in Lab. " << Pproj << 266 G4cout << "Ptarget in Lab. " << Ptarg << 267 G4cout << "G4SingleDiffractiveExcitation- << 268 G4cout << "G4SingleDiffractiveExcitation- << 269 #endif << 270 104 271 target->Set4Momentum(Ptarget); << 105 G4LorentzRotation toLab(toCms.inverse()); 272 projectile->Set4Momentum(Pprojectile); << 106 >> 107 Pprojectile.transform(toCms); >> 108 Ptarget.transform(toCms); 273 109 274 return true; << 110 G4LorentzVector Qmomentum; >> 111 G4int whilecount=0; >> 112 do { >> 113 // Generate pt >> 114 >> 115 G4double maxPtSquare=sqr(Ptarget.pz()); >> 116 if (whilecount++ >= 500 && (whilecount%100)==0) >> 117 // G4cout << "G4SingleDiffractiveExcitation::ExciteParticipants possibly looping" >> 118 // << ", loop count/ maxPtSquare : " >> 119 // << whilecount << " / " << maxPtSquare << G4endl; >> 120 if (whilecount > 1000 ) >> 121 { >> 122 Qmomentum=G4LorentzVector(0.,0.,0.,0.); >> 123 // G4cout << "G4SingleDiffractiveExcitation::ExciteParticipants: Aborting loop!" << G4endl; >> 124 return false; // Ignore this interaction >> 125 } >> 126 Qmomentum=G4LorentzVector(GaussianPt(widthOfPtSquare,maxPtSquare),0); >> 127 >> 128 >> 129 // Momentum transfer >> 130 G4double Xmin = minmass / ( Pprojectile.e() + Ptarget.e() ); >> 131 G4double Xmax=1.; >> 132 G4double Xplus =ChooseX(Xmin,Xmax); >> 133 G4double Xminus=ChooseX(Xmin,Xmax); >> 134 >> 135 G4double pt2=G4ThreeVector(Qmomentum.vect()).mag2(); >> 136 G4double Qplus =-1 * pt2 / Xminus/Ptarget.minus(); >> 137 G4double Qminus= pt2 / Xplus /Pprojectile.plus(); >> 138 >> 139 if ( KeepProjectile ) >> 140 { >> 141 Qminus = (sqr(projectile->GetDefinition()->GetPDGMass()) + pt2 ) >> 142 / (Pprojectile.plus() + Qplus ) >> 143 - Pprojectile.minus(); >> 144 } else >> 145 { >> 146 Qplus = Ptarget.plus() >> 147 - (sqr(target->GetDefinition()->GetPDGMass()) + pt2 ) >> 148 / (Ptarget.minus() - Qminus ); >> 149 } >> 150 >> 151 Qmomentum.setPz( (Qplus-Qminus)/2 ); >> 152 Qmomentum.setE( (Qplus+Qminus)/2 ); >> 153 >> 154 // G4cout << "Qplus / Qminus " << Qplus << " / " << Qminus<<G4endl; >> 155 // G4cout << "pt2 " << pt2 << G4endl; >> 156 // G4cout << "Qmomentum " << Qmomentum << G4endl; >> 157 // G4cout << " Masses (P/T) : " << (Pprojectile+Qmomentum).mag() << >> 158 // " / " << (Ptarget-Qmomentum).mag() << G4endl; >> 159 >> 160 } while ( (Ptarget-Qmomentum).mag2() <= Mtarget2 >> 161 || (Pprojectile+Qmomentum).mag2() <= Mprojectile2 >> 162 || (Ptarget-Qmomentum).e() < 0. >> 163 || (Pprojectile+Qmomentum).e() < 0. ); >> 164 >> 165 >> 166 // G4double Ecms=Pprojectile.e() + Ptarget.e(); >> 167 >> 168 Pprojectile += Qmomentum; >> 169 >> 170 Ptarget -= Qmomentum; >> 171 >> 172 // G4cout << "Pprojectile.e() : " << Pprojectile.e() << G4endl; >> 173 // G4cout << "Ptarget.e() : " << Ptarget.e() << G4endl; >> 174 >> 175 // G4cout << "end event_______________________________________________"<<G4endl; >> 176 // >> 177 >> 178 >> 179 // G4cout << "Pprojectile with Q : " << Pprojectile << G4endl; >> 180 // G4cout << "Ptarget with Q : " << Ptarget << G4endl; >> 181 // G4cout << "Projectile back: " << toLab * Pprojectile << G4endl; >> 182 // G4cout << "Target back: " << toLab * Ptarget << G4endl; >> 183 >> 184 // Transform back and update SplitableHadron Participant. >> 185 Pprojectile.transform(toLab); >> 186 Ptarget.transform(toLab); >> 187 >> 188 // G4cout << "G4SingleDiffractiveExcitation- Target mass " << Ptarget.mag() << G4endl; >> 189 // G4cout << "G4SingleDiffractiveExcitation- Projectile mass " << Pprojectile.mag() << G4endl; >> 190 >> 191 target->Set4Momentum(Ptarget); >> 192 projectile->Set4Momentum(Pprojectile); >> 193 >> 194 >> 195 return true; 275 } 196 } 276 197 >> 198 >> 199 >> 200 277 // --------- private methods ----------------- 201 // --------- private methods ---------------------- 278 202 279 G4double G4SingleDiffractiveExcitation::Choose 203 G4double G4SingleDiffractiveExcitation::ChooseX(G4double Xmin, G4double Xmax) const 280 { 204 { 281 // choose an x between Xmin and Xmax with P( << 205 // choose an x between Xmin and Xmax with P(x) ~ 1/x 282 G4double range=Xmax-Xmin; << 283 206 284 if ( Xmin <= 0. || range <=0. ) << 207 // to be improved... 285 { << 286 G4cout << " Xmin, range : " << Xmin << " , << 287 throw G4HadronicException(__FILE__, __LINE << 288 } << 289 << 290 G4double x = Xmin*G4Pow::GetInstance()->powA << 291 // G4double x = 1.0/sqr(1.0/std::sqrt(Xmin) << 292 return x; << 293 } << 294 208 >> 209 G4double range=Xmax-Xmin; >> 210 >> 211 if ( Xmin <= 0. || range <=0. ) >> 212 { >> 213 G4cout << " Xmin, range : " << Xmin << " , " << range << G4endl; >> 214 throw G4HadronicException(__FILE__, __LINE__, "G4SingleDiffractiveExcitation::ChooseX : Invalid arguments "); >> 215 } >> 216 >> 217 G4double x; >> 218 do { >> 219 x=Xmin + G4UniformRand() * range; >> 220 } while ( Xmin/x < G4UniformRand() ); >> 221 >> 222 // cout << "DiffractiveX "<<x<<G4endl; >> 223 return x; >> 224 } 295 225 296 G4ThreeVector G4SingleDiffractiveExcitation::G 226 G4ThreeVector G4SingleDiffractiveExcitation::GaussianPt(G4double widthSquare, G4double maxPtSquare) const 297 { // @@ this method is used in FTF 227 { // @@ this method is used in FTFModel as well. Should go somewhere common! >> 228 >> 229 G4double pt2; >> 230 >> 231 do { >> 232 pt2=widthSquare * std::log( G4UniformRand() ); >> 233 } while ( pt2 > maxPtSquare); >> 234 >> 235 pt2=std::sqrt(pt2); >> 236 >> 237 G4double phi=G4UniformRand() * twopi; >> 238 >> 239 return G4ThreeVector (pt2*std::cos(phi), pt2*std::sin(phi), 0.); >> 240 } >> 241 >> 242 G4SingleDiffractiveExcitation::G4SingleDiffractiveExcitation(const G4SingleDiffractiveExcitation &) >> 243 : G4QGSDiffractiveExcitation(), >> 244 widthOfPtSquare(0) , minExtraMass(0), >> 245 minmass(0) >> 246 { >> 247 throw G4HadronicException(__FILE__, __LINE__, "G4SingleDiffractiveExcitation copy contructor not meant to be called"); >> 248 } >> 249 >> 250 >> 251 G4SingleDiffractiveExcitation::~G4SingleDiffractiveExcitation() >> 252 { >> 253 } 298 254 299 G4double pt2; << 300 255 301 const G4int maxNumberOfLoops = 1000; << 256 const G4SingleDiffractiveExcitation & G4SingleDiffractiveExcitation::operator=(const G4SingleDiffractiveExcitation &) 302 G4int loopCounter = 0; << 257 { 303 do { << 258 throw G4HadronicException(__FILE__, __LINE__, "G4SingleDiffractiveExcitation = operator meant to be called"); 304 pt2=-widthSquare * G4Log( G4UniformRand() << 259 return *this; 305 } while ( ( pt2 > maxPtSquare) && ++loopCoun << 260 } 306 if ( loopCounter >= maxNumberOfLoops ) { << 307 pt2 = 0.99*maxPtSquare; // Just an accept << 308 } << 309 261 310 pt2=std::sqrt(pt2); << 311 262 312 G4double phi=G4UniformRand() * twopi; << 263 int G4SingleDiffractiveExcitation::operator==(const G4SingleDiffractiveExcitation &) const >> 264 { >> 265 throw G4HadronicException(__FILE__, __LINE__, "G4SingleDiffractiveExcitation == operator meant to be called"); >> 266 return false; >> 267 } 313 268 314 return G4ThreeVector (pt2*std::cos(phi), pt2 << 269 int G4SingleDiffractiveExcitation::operator!=(const G4SingleDiffractiveExcitation &) const >> 270 { >> 271 throw G4HadronicException(__FILE__, __LINE__, "G4SingleDiffractiveExcitation != operator meant to be called"); >> 272 return true; 315 } 273 } >> 274 >> 275 >> 276 >> 277 316 278 317 279