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 // $Id: G4DiffractiveExcitation.cc,v 1.23 2010-11-15 10:02:38 vuzhinsk Exp $ 28 << 29 // ------------------------------------------- 28 // ------------------------------------------------------------ 30 // GEANT 4 class implemetation file 29 // GEANT 4 class implemetation file 31 // 30 // 32 // ---------------- G4DiffractiveExcitati 31 // ---------------- G4DiffractiveExcitation -------------- 33 // by Gunter Folger, October 1998. 32 // by Gunter Folger, October 1998. 34 // diffractive Excitation used by strin 33 // diffractive Excitation used by strings models 35 // Take a projectile and a target 34 // Take a projectile and a target 36 // excite the projectile and targe 35 // excite the projectile and target 37 // Essential changed by V. Uzhinsky in Novemb 36 // Essential changed by V. Uzhinsky in November - December 2006 38 // in order to put it in a correspondence wit 37 // in order to put it in a correspondence with original FRITIOF 39 // model. Variant of FRITIOF with nucleon de- 38 // model. Variant of FRITIOF with nucleon de-excitation is implemented. 40 // Other changes by V.Uzhinsky in May 2007 we 39 // Other changes by V.Uzhinsky in May 2007 were introduced to fit 41 // meson-nucleon interactions. Additional cha 40 // meson-nucleon interactions. Additional changes by V. Uzhinsky 42 // were introduced in December 2006. They tre 41 // were introduced in December 2006. They treat diffraction dissociation 43 // processes more exactly. 42 // processes more exactly. 44 // Correct treatment of the diffraction disso << 45 // Mass distributions for resonances and uu-d << 46 // and dd-diquarks suppression in neutrons we << 47 // ------------------------------------------- 43 // --------------------------------------------------------------------- 48 44 >> 45 49 #include "globals.hh" 46 #include "globals.hh" 50 #include "Randomize.hh" 47 #include "Randomize.hh" 51 #include "G4PhysicalConstants.hh" << 52 #include "G4SystemOfUnits.hh" << 53 48 54 #include "G4DiffractiveExcitation.hh" 49 #include "G4DiffractiveExcitation.hh" 55 #include "G4FTFParameters.hh" 50 #include "G4FTFParameters.hh" 56 #include "G4ElasticHNScattering.hh" 51 #include "G4ElasticHNScattering.hh" 57 52 >> 53 #include "G4LorentzRotation.hh" 58 #include "G4RotationMatrix.hh" 54 #include "G4RotationMatrix.hh" >> 55 #include "G4ThreeVector.hh" 59 #include "G4ParticleDefinition.hh" 56 #include "G4ParticleDefinition.hh" 60 #include "G4ParticleTable.hh" << 61 #include "G4SampleResonance.hh" << 62 #include "G4VSplitableHadron.hh" 57 #include "G4VSplitableHadron.hh" 63 #include "G4ExcitedString.hh" 58 #include "G4ExcitedString.hh" >> 59 #include "G4ParticleTable.hh" 64 #include "G4Neutron.hh" 60 #include "G4Neutron.hh" 65 << 61 #include "G4ParticleDefinition.hh" 66 #include "G4Exp.hh" << 67 #include "G4Log.hh" << 68 #include "G4Pow.hh" << 69 62 70 //#include "G4ios.hh" 63 //#include "G4ios.hh" >> 64 //#include "UZHI_diffraction.hh" 71 65 72 //============================================ << 66 G4DiffractiveExcitation::G4DiffractiveExcitation() 73 << 67 { 74 //#define debugFTFexcitation << 75 //#define debug_heavyHadrons << 76 << 77 //============================================ << 78 << 79 G4DiffractiveExcitation::G4DiffractiveExcitati << 80 << 81 << 82 //============================================ << 83 << 84 G4DiffractiveExcitation::~G4DiffractiveExcitat << 85 << 86 << 87 //============================================ << 88 << 89 G4bool G4DiffractiveExcitation::ExciteParticip << 90 << 91 << 92 << 93 << 94 #ifdef debugFTFexcitation << 95 G4cout << G4endl << "FTF ExciteParticipants << 96 #endif << 97 << 98 CommonVariables common; << 99 << 100 // Projectile parameters << 101 common.Pprojectile = projectile->Get4Momentu << 102 if ( common.Pprojectile.z() < 0.0 ) return f << 103 common.ProjectilePDGcode = projectile->GetDe << 104 common.absProjectilePDGcode = std::abs( comm << 105 common.M0projectile = projectile->GetDefinit << 106 G4double ProjectileRapidity = common.Pprojec << 107 << 108 // Target parameters << 109 common.Ptarget = target->Get4Momentum(); << 110 common.TargetPDGcode = target->GetDefinition << 111 common.absTargetPDGcode = std::abs( common.T << 112 common.M0target = target->GetDefinition()->G << 113 G4double TargetRapidity = common.Ptarget.rap << 114 << 115 // Kinematical properties of the interaction << 116 G4LorentzVector Psum = common.Pprojectile + << 117 common.S = Psum.mag2(); << 118 common.SqrtS = std::sqrt( common.S ); << 119 << 120 // Check off-shellness of the participants << 121 G4bool toBePutOnMassShell = true; //Uzhi Au << 122 common.MminProjectile = common.BrW.GetMinimu << 123 /* Uzhi Aug.2019 << 124 if ( common.M0projectile < common.MminProjec << 125 toBePutOnMassShell = true; << 126 common.M0projectile = common.BrW.SampleMas << 127 << 128 << 129 } << 130 */ << 131 common.M0projectile2 = common.M0projectile * << 132 common.ProjectileDiffStateMinMass = thePa << 133 common.ProjectileNonDiffStateMinMass = thePa << 134 if ( common.M0projectile > common.Projectile << 135 common.ProjectileDiffStateMinMass = com << 136 common.ProjectileNonDiffStateMinMass = com << 137 if ( common.absProjectilePDGcode > 3000 ) << 138 common.ProjectileDiffStateMinMass += << 139 common.ProjectileNonDiffStateMinMass += << 140 } << 141 } << 142 common.MminTarget = common.BrW.GetMinimumMas << 143 /* Uzhi Aug.2019 << 144 if ( common.M0target < common.MminTarget ) { << 145 toBePutOnMassShell = true; << 146 common.M0target = common.BrW.SampleMass( t << 147 t << 148 + << 149 } << 150 */ << 151 common.M0target2 = common.M0target * common. << 152 common.TargetDiffStateMinMass = theParame << 153 common.TargetNonDiffStateMinMass = theParame << 154 if ( common.M0target > common.TargetDiffStat << 155 common.TargetDiffStateMinMass = common. << 156 common.TargetNonDiffStateMinMass = common. << 157 if ( common.absTargetPDGcode > 3000 ) { / << 158 common.TargetDiffStateMinMass += 140. << 159 common.TargetNonDiffStateMinMass += 140. << 160 } << 161 }; << 162 #ifdef debugFTFexcitation << 163 G4cout << "Proj Targ PDGcodes " << common.Pr << 164 << "Mprojectile Y " << common.Pproje << 165 << "M0projectile Y " << common.M0proj << 166 G4cout << "Mtarget Y " << common.Ptarge << 167 << "M0target Y " << common.M0targ << 168 G4cout << "Pproj " << common.Pprojectile << << 169 #endif << 170 << 171 // Transform momenta to cms and then rotate << 172 common.toCms = G4LorentzRotation( -1 * Psum. << 173 G4LorentzVector Ptmp = common.toCms * common << 174 if ( Ptmp.pz() <= 0.0 ) return false; // "S << 175 common.toCms.rotateZ( -1*Ptmp.phi() ); << 176 common.toCms.rotateY( -1*Ptmp.theta() ); << 177 common.toLab = common.toCms.inverse(); << 178 common.Pprojectile.transform( common.toCms ) << 179 common.Ptarget.transform( common.toCms ); << 180 << 181 G4double SumMasses = common.M0projectile + c << 182 #ifdef debugFTFexcitation << 183 G4cout << "SqrtS " << common.SqrtS << G4 << 184 << " " << common.M0target << " " << S << 185 #endif << 186 if ( common.SqrtS < SumMasses ) return false << 187 << 188 common.PZcms2 = ( sqr( common.S ) + sqr( com << 189 - 2.0 * ( common.S * ( com << 190 + common.M0proje << 191 #ifdef debugFTFexcitation << 192 G4cout << "PZcms2 after toBePutOnMassShell " << 193 #endif << 194 if ( common.PZcms2 < 0.0 ) return false; // << 195 << 196 common.PZcms = std::sqrt( common.PZcms2 ); << 197 if ( toBePutOnMassShell ) { << 198 if ( common.Pprojectile.z() > 0.0 ) { << 199 common.Pprojectile.setPz( common.PZcms << 200 common.Ptarget.setPz( -common.PZcms << 201 } else { << 202 common.Pprojectile.setPz( -common.PZcms << 203 common.Ptarget.setPz( common.PZcms << 204 } << 205 common.Pprojectile.setE( std::sqrt( common << 206 + comm << 207 + comm << 208 + comm << 209 common.Ptarget.setE( std::sqrt( common.M0t << 210 + common.P << 211 + common.P << 212 + common.P << 213 } << 214 #ifdef debugFTFexcitation << 215 G4cout << "Start --------------------" << G4 << 216 << " " << common.ProjectileDiffStateM << 217 << G4endl << 218 << "Targ M0 Mdif Mndif " << common.M0 << 219 << " " << common.TargetNonDiffStateMi << 220 << "Proj CMS " << common.Pprojectile << 221 #endif << 222 << 223 // Check for possible quark exchange << 224 ProjectileRapidity = common.Pprojectile.rapi << 225 TargetRapidity = common.Ptarget.rapidity(); << 226 G4double QeNoExc = theParameters->GetProcPro << 227 G4double QeExc = theParameters->GetProcPro << 228 theParameters->GetProcPro << 229 common.ProbProjectileDiffraction = << 230 theParameters->GetProcProb( 2, ProjectileR << 231 common.ProbTargetDiffraction = << 232 theParameters->GetProcProb( 3, ProjectileR << 233 common.ProbOfDiffraction = common.ProbProjec << 234 #ifdef debugFTFexcitation << 235 G4cout << "Proc Probs " << QeNoExc << " " << << 236 << common.ProbProjectileDiffraction < << 237 << "ProjectileRapidity " << Projectil << 238 #endif << 239 << 240 if ( QeNoExc + QeExc + common.ProbProjectile << 241 QeNoExc = 1.0 - QeExc - common.ProbProject << 242 } << 243 if ( QeExc + QeNoExc != 0.0 ) { << 244 common.ProbExc = QeExc / ( QeExc + QeNoExc << 245 } << 246 if ( 1.0 - QeExc - QeNoExc > 0.0 ) { << 247 common.ProbProjectileDiffraction /= ( 1.0 << 248 common.ProbTargetDiffraction /= ( 1.0 << 249 } << 250 #ifdef debugFTFexcitation << 251 G4cout << "Proc Probs " << QeNoExc << " " << << 252 << common.ProbProjectileDiffraction < << 253 << "ProjectileRapidity " << Projectil << 254 #endif << 255 << 256 // Try out charge exchange << 257 G4int returnCode = 1; << 258 if ( G4UniformRand() < QeExc + QeNoExc ) { << 259 returnCode = ExciteParticipants_doChargeEx << 260 << 261 } << 262 << 263 G4bool returnResult = false; << 264 if ( returnCode == 0 ) { << 265 returnResult = true; // Successfully ende << 266 } else if ( returnCode == 1 ) { << 267 << 268 common.ProbOfDiffraction = common.ProbProj << 269 #ifdef debugFTFexcitation << 270 G4cout << "Excitation -------------------- << 271 << "Proj M0 MdMin MndMin " << commo << 272 << common.ProjectileDiffStateMinMas << 273 << G4endl << 274 << "Targ M0 MdMin MndMin " << commo << 275 << " " << common.TargetNonDiffState << 276 << G4endl << 277 << "Prob: ProjDiff TargDiff + Sum " << 278 << common.ProbTargetDiffraction << << 279 #endif << 280 if ( common.ProbOfDiffraction != 0.0 ) { << 281 common.ProbProjectileDiffraction /= comm << 282 } else { << 283 common.ProbProjectileDiffraction = 0.0; << 284 } << 285 #ifdef debugFTFexcitation << 286 G4cout << "Prob: ProjDiff TargDiff + Sum " << 287 << common.ProbTargetDiffraction << << 288 #endif << 289 common.ProjectileDiffStateMinMass2 = sq << 290 common.ProjectileNonDiffStateMinMass2 = sq << 291 common.TargetDiffStateMinMass2 = sq << 292 common.TargetNonDiffStateMinMass2 = sq << 293 // Choose between diffraction and non-diff << 294 if ( G4UniformRand() < common.ProbOfDiffra << 295 << 296 returnResult = ExciteParticipants_doDiff << 297 << 298 } else { << 299 << 300 returnResult = ExciteParticipants_doNonD << 301 << 302 } << 303 if ( returnResult ) { << 304 common.Pprojectile += common.Qmomentum; << 305 common.Ptarget -= common.Qmomentum; << 306 // Transform back and update SplitableHa << 307 common.Pprojectile.transform( common.toL << 308 common.Ptarget.transform( common.toLab ) << 309 #ifdef debugFTFexcitation << 310 G4cout << "Mproj " << common.Pprojectile << 311 << G4endl; << 312 #endif << 313 projectile->Set4Momentum( common.Pprojec << 314 target->Set4Momentum( common.Ptarget ); << 315 projectile->IncrementCollisionCount( 1 ) << 316 target->IncrementCollisionCount( 1 ); << 317 } << 318 } << 319 << 320 return returnResult; << 321 } 68 } 322 69 323 //-------------------------------------------- << 70 // --------------------------------------------------------------------- >> 71 G4bool G4DiffractiveExcitation:: >> 72 ExciteParticipants(G4VSplitableHadron *projectile, >> 73 G4VSplitableHadron *target, >> 74 G4FTFParameters *theParameters, >> 75 G4ElasticHNScattering *theElastic) const >> 76 { >> 77 // -------------------- Projectile parameters ----------------------- >> 78 G4LorentzVector Pprojectile=projectile->Get4Momentum(); >> 79 >> 80 if(Pprojectile.z() < 0.) >> 81 { >> 82 target->SetStatus(2); >> 83 return false; >> 84 } >> 85 >> 86 G4double ProjectileRapidity = Pprojectile.rapidity(); >> 87 >> 88 G4int ProjectilePDGcode=projectile->GetDefinition()->GetPDGEncoding(); >> 89 G4int absProjectilePDGcode=std::abs(ProjectilePDGcode); >> 90 >> 91 G4bool PutOnMassShell(false); >> 92 // G4double M0projectile=projectile->GetDefinition()->GetPDGMass(); // With de-excitation >> 93 G4double M0projectile = Pprojectile.mag(); // Without de-excitation >> 94 >> 95 if(M0projectile < projectile->GetDefinition()->GetPDGMass()) >> 96 { >> 97 PutOnMassShell=true; >> 98 M0projectile=projectile->GetDefinition()->GetPDGMass(); >> 99 } >> 100 >> 101 G4double M0projectile2 = M0projectile * M0projectile; >> 102 >> 103 G4double ProjectileDiffStateMinMass=theParameters->GetProjMinDiffMass(); >> 104 G4double ProjectileNonDiffStateMinMass=theParameters->GetProjMinNonDiffMass(); >> 105 G4double ProbProjectileDiffraction=theParameters->GetProbabilityOfProjDiff(); >> 106 >> 107 // -------------------- Target parameters ------------------------- >> 108 G4int TargetPDGcode=target->GetDefinition()->GetPDGEncoding(); >> 109 G4int absTargetPDGcode=std::abs(TargetPDGcode); >> 110 //G4cout<<"Excit "<<ProjectilePDGcode<<" "<<TargetPDGcode<<G4endl; >> 111 >> 112 G4LorentzVector Ptarget=target->Get4Momentum(); >> 113 >> 114 G4double M0target = Ptarget.mag(); >> 115 >> 116 // G4double TargetRapidity = Ptarget.rapidity(); >> 117 >> 118 if(M0target < target->GetDefinition()->GetPDGMass()) >> 119 { >> 120 PutOnMassShell=true; >> 121 M0target=target->GetDefinition()->GetPDGMass(); >> 122 } 324 123 325 G4int G4DiffractiveExcitation:: << 124 G4double M0target2 = M0target * M0target; 326 ExciteParticipants_doChargeExchange( G4VSplita << 125 327 G4VSplita << 126 G4double TargetDiffStateMinMass=theParameters->GetTarMinDiffMass(); 328 G4FTFPara << 127 G4double TargetNonDiffStateMinMass=theParameters->GetTarMinNonDiffMass(); 329 G4Elastic << 128 G4double ProbTargetDiffraction=theParameters->GetProbabilityOfTarDiff(); 330 G4Diffrac << 129 331 // First of the three utility methods used o << 130 G4double AveragePt2=theParameters->GetAveragePt2(); 332 // it does the sampling for the charge-excha << 131 333 // This method returns an integer code - ins << 132 // G4double ProbOfDiffraction=ProbProjectileDiffraction + 334 // "0" : successfully ended and nothing el << 133 // ProbTargetDiffraction; 335 // "1" : successfully completed, but the w << 134 336 // "99" : unsuccessfully ended, nothing els << 135 G4double SumMasses=M0projectile+M0target+200.*MeV; 337 G4int returnCode = 99; << 136 338 << 137 // Kinematical properties of the interactions -------------- 339 G4double DeltaProbAtQuarkExchange = theParam << 138 G4LorentzVector Psum; // 4-momentum in CMS 340 G4ParticleDefinition* TestParticle = 0; << 139 Psum=Pprojectile+Ptarget; 341 G4double MtestPr = 0.0, MtestTr = 0.0; << 140 G4double S=Psum.mag2(); 342 << 141 343 #ifdef debugFTFexcitation << 142 // Transform momenta to cms and then rotate parallel to z axis; 344 G4cout << "Q exchange ---------------------- << 143 G4LorentzRotation toCms(-1*Psum.boostVector()); 345 #endif << 144 346 << 145 G4LorentzVector Ptmp=toCms*Pprojectile; 347 // The target hadron is always a nucleon (i. << 146 348 // never an antinucleon), therefore only a q << 147 if ( Ptmp.pz() <= 0. ) 349 // exchanged between the projectile hadron a << 148 { 350 // we could get a quark-quark-antiquark syst << 149 target->SetStatus(2); 351 // This implies that any projectile meson or << 150 // "String" moving backwards in CMS, abort collision !! 352 // a constituent quark in all cases - can ha << 151 return false; 353 // hadron. Instead, any projectile anti-bary << 152 } 354 // with a target hadron (because it has only << 153 355 // projectile baryons, instead can have char << 154 toCms.rotateZ(-1*Ptmp.phi()); 356 << 155 toCms.rotateY(-1*Ptmp.theta()); 357 G4int NewProjCode = 0, NewTargCode = 0, Proj << 156 358 // Projectile unpacking << 157 G4LorentzRotation toLab(toCms.inverse()); 359 if ( common.absProjectilePDGcode < 1000 ) { << 158 360 UnpackMeson( common.ProjectilePDGcode, Pr << 159 Pprojectile.transform(toCms); 361 } else { << 160 Ptarget.transform(toCms); 362 UnpackBaryon( common.ProjectilePDGcode, Pr << 161 363 } << 162 G4double PZcms2, PZcms; 364 // Target unpacking << 163 365 G4int TargQ1 = 0, TargQ2 = 0, TargQ3 = 0; << 164 G4double SqrtS=std::sqrt(S); 366 UnpackBaryon( common.TargetPDGcode, TargQ1, << 165 //G4cout<<"SqrtS < 2300*MeV "<<SqrtS<<G4endl; 367 #ifdef debugFTFexcitation << 166 if(absProjectilePDGcode > 1000 && (SqrtS < 2300*MeV || SqrtS < SumMasses)) 368 G4cout << "Proj Quarks " << ProjQ1 << " " << << 167 {target->SetStatus(2); return false;} // The model cannot work for 369 << "Targ Quarks " << TargQ1 << " " << << 168 // p+p-interactions 370 #endif << 169 // at Plab < 1.62 GeV/c. 371 << 170 372 // Sampling of exchanged quarks << 171 if(( absProjectilePDGcode == 211 || ProjectilePDGcode == 111) && 373 G4int ProjExchangeQ = 0, TargExchangeQ = 0; << 172 ((SqrtS < 1600*MeV) || (SqrtS < SumMasses))) 374 if ( common.absProjectilePDGcode < 1000 ) { << 173 {target->SetStatus(2); return false;} // The model cannot work for 375 << 174 // Pi+p-interactions 376 G4bool isProjQ1Quark = false; << 175 // at Plab < 1. GeV/c. 377 ProjExchangeQ = ProjQ2; << 176 378 if ( ProjQ1 > 0 ) { // ProjQ1 is a quark << 177 if(( (absProjectilePDGcode == 321) || (ProjectilePDGcode == -311) || 379 isProjQ1Quark = true; << 178 (absProjectilePDGcode == 311) || (absProjectilePDGcode == 130) || 380 ProjExchangeQ = ProjQ1; << 179 (absProjectilePDGcode == 310)) && ((SqrtS < 1600*MeV) || (SqrtS < SumMasses))) 381 } << 180 {target->SetStatus(2); return false;} // The model cannot work for 382 // Exchange of non-identical quarks is all << 181 // K+p-interactions 383 G4int NpossibleStates = 0; << 182 // at Plab < ??? GeV/c. ??? 384 if ( ProjExchangeQ != TargQ1 ) NpossibleSt << 183 385 if ( ProjExchangeQ != TargQ2 ) NpossibleSt << 184 PZcms2=(S*S+M0projectile2*M0projectile2+M0target2*M0target2- 386 if ( ProjExchangeQ != TargQ3 ) NpossibleSt << 185 2*S*M0projectile2 - 2*S*M0target2 - 2*M0projectile2*M0target2) 387 G4int Nsampled = (G4int)G4RandFlat::shootI << 186 /4./S; 388 NpossibleStates = 0; << 187 389 if ( ProjExchangeQ != TargQ1 ) { << 188 if(PZcms2 < 0) 390 if ( ++NpossibleStates == Nsampled ) { << 189 {target->SetStatus(2); return false;} // It can be in an interaction with 391 TargExchangeQ = TargQ1; TargQ1 = ProjE << 190 // off-shell nuclear nucleon 392 isProjQ1Quark ? ProjQ1 = TargExchangeQ << 191 393 } << 192 PZcms = std::sqrt(PZcms2); 394 } << 193 395 if ( ProjExchangeQ != TargQ2 ) { << 194 if(PutOnMassShell) 396 if ( ++NpossibleStates == Nsampled ) { << 195 { 397 TargExchangeQ = TargQ2; TargQ2 = ProjE << 196 if(Pprojectile.z() > 0.) 398 isProjQ1Quark ? ProjQ1 = TargExchangeQ << 197 { 399 } << 198 Pprojectile.setPz( PZcms); 400 } << 199 Ptarget.setPz( -PZcms); 401 if ( ProjExchangeQ != TargQ3 ) { << 200 } 402 if ( ++NpossibleStates == Nsampled ) { << 201 else 403 TargExchangeQ = TargQ3; TargQ3 = ProjE << 202 { 404 isProjQ1Quark ? ProjQ1 = TargExchangeQ << 203 Pprojectile.setPz(-PZcms); 405 } << 204 Ptarget.setPz( PZcms); 406 } << 205 }; 407 #ifdef debugFTFexcitation << 206 408 G4cout << "Exchanged Qs in Pr Tr " << Proj << 207 Pprojectile.setE(std::sqrt(M0projectile2 + 409 #endif << 208 Pprojectile.x()*Pprojectile.x()+ 410 << 209 Pprojectile.y()*Pprojectile.y()+ 411 G4int aProjQ1 = std::abs( ProjQ1 ), aProjQ << 210 PZcms2)); 412 G4bool ProjExcited = false; << 211 Ptarget.setE(std::sqrt(M0target2 + 413 const G4int maxNumberOfAttempts = 50; << 212 Ptarget.x()*Ptarget.x()+ 414 G4int attempts = 0; << 213 Ptarget.y()*Ptarget.y()+ 415 while ( attempts++ < maxNumberOfAttempts ) << 214 PZcms2)); 416 << 215 } 417 // Determination of a new projectile ID << 216 418 G4double ProbSpin0 = 0.5; << 217 G4double maxPtSquare; // = PZcms2; 419 G4double Ksi = G4UniformRand(); << 218 /* 420 if ( aProjQ1 == aProjQ2 ) { << 219 G4cout<<"Start --------------------"<<G4endl; 421 if ( G4UniformRand() < ProbSpin0 ) { << 220 G4cout<<"Proj "<<M0projectile<<" "<<ProjectileDiffStateMinMass<<" "<<ProjectileNonDiffStateMinMass<<G4endl; 422 if ( aProjQ1 < 3 ) { << 221 G4cout<<"Targ "<<M0target <<" "<<TargetDiffStateMinMass <<" "<<TargetNonDiffStateMinMass<<G4endl; 423 NewProjCode = 111; << 222 G4cout<<"SqrtS "<<SqrtS<<G4endl; 424 if ( Ksi < 0.5 ) { << 223 G4cout<<"Rapid "<<ProjectileRapidity<<G4endl; //" "<<TargetRapidity<<G4endl; 425 NewProjCode = 221; << 224 */ 426 if ( Ksi < 0.25 ) { << 225 // Charge exchange can be possible for baryons ----------------- 427 NewProjCode = 331; << 226 428 } << 227 // Getting the values needed for exchange ---------------------- 429 } << 228 G4double MagQuarkExchange =theParameters->GetMagQuarkExchange(); 430 } else if ( aProjQ1 == 3 ) { << 229 G4double SlopeQuarkExchange =theParameters->GetSlopeQuarkExchange(); 431 NewProjCode = 221; << 230 G4double DeltaProbAtQuarkExchange=theParameters->GetDeltaProbAtQuarkExchange(); 432 if ( Ksi < 0.5 ) { << 231 433 NewProjCode = 331; << 232 // G4double NucleonMass= 434 } << 233 // (G4ParticleTable::GetParticleTable()->FindParticle(2112))->GetPDGMass(); 435 } else if ( aProjQ1 == 4 ) { << 234 G4double DeltaMass= 436 NewProjCode = 441; // eta << 235 (G4ParticleTable::GetParticleTable()->FindParticle(2224))->GetPDGMass(); 437 } else if ( aProjQ1 == 5 ) { << 236 438 NewProjCode = 551; // eta << 237 //G4cout<<MagQuarkExchange*std::exp(-SlopeQuarkExchange*(ProjectileRapidity - TargetRapidity))<<G4endl; 439 } << 238 440 } else { << 239 //G4cout<<"Q exc Mag Slop Wdelta"<<MagQuarkExchange<<" "<<SlopeQuarkExchange<<" "<<DeltaProbAtQuarkExchange<<G4endl; 441 if ( aProjQ1 < 3 ) { << 240 //G4cout<<"ProjectileRapidity "<<ProjectileRapidity<<G4endl; 442 NewProjCode = 113; << 241 //G4cout<<MagQuarkExchange*std::exp(-SlopeQuarkExchange*(ProjectileRapidity))<<G4endl; 443 if ( Ksi < 0.5 ) { << 242 //G4int Uzhi; G4cin>>Uzhi; 444 NewProjCode = 223; << 243 // Check for possible quark exchange ----------------------------------- 445 } << 244 446 } else if ( aProjQ1 == 3 ) { << 245 if(G4UniformRand() < MagQuarkExchange* 447 NewProjCode = 333; << 246 std::exp(-SlopeQuarkExchange*ProjectileRapidity)) //TargetRapidity))) 1.45 448 } else if ( aProjQ1 == 4 ) { << 247 { 449 NewProjCode = 443; // J/p << 248 // std::exp(-SlopeQuarkExchange*(ProjectileRapidity - 1.36))) //TargetRapidity))) 1.45 450 } else if ( aProjQ1 == 5 ) { << 249 //G4cout<<"Q exchange"<<G4endl; 451 NewProjCode = 553; // Ups << 250 G4int NewProjCode(0), NewTargCode(0); 452 } << 251 >> 252 G4int ProjQ1(0), ProjQ2(0), ProjQ3(0); >> 253 >> 254 // Projectile unpacking -------------------------- >> 255 if(absProjectilePDGcode < 1000 ) >> 256 { // projectile is meson ----------------- >> 257 UnpackMeson(ProjectilePDGcode, ProjQ1, ProjQ2); >> 258 } else >> 259 { // projectile is baryon ---------------- >> 260 UnpackBaryon(ProjectilePDGcode, ProjQ1, ProjQ2, ProjQ3); >> 261 } // End of the hadron's unpacking ---------- >> 262 >> 263 // Target unpacking ------------------------------ >> 264 G4int TargQ1(0), TargQ2(0), TargQ3(0); >> 265 UnpackBaryon(TargetPDGcode, TargQ1, TargQ2, TargQ3); >> 266 >> 267 //G4cout<<ProjQ1<<" "<<ProjQ2<<" "<<ProjQ3<<G4endl; >> 268 //G4cout<<TargQ1<<" "<<TargQ2<<" "<<TargQ3<<G4endl; >> 269 // Sampling of exchanged quarks ------------------- >> 270 G4int ProjExchangeQ(0); >> 271 G4int TargExchangeQ(0); >> 272 >> 273 if(absProjectilePDGcode < 1000 ) >> 274 { // projectile is meson ----------------- >> 275 >> 276 if(ProjQ1 > 0 ) // ProjQ1 is quark >> 277 { >> 278 ProjExchangeQ = ProjQ1; >> 279 if(ProjExchangeQ != TargQ1) >> 280 { >> 281 TargExchangeQ = TargQ1; TargQ1=ProjExchangeQ; ProjQ1=TargExchangeQ; >> 282 } else >> 283 if(ProjExchangeQ != TargQ2) >> 284 { >> 285 TargExchangeQ = TargQ2; TargQ2=ProjExchangeQ; ProjQ1=TargExchangeQ; >> 286 } else >> 287 { >> 288 TargExchangeQ = TargQ3; TargQ3=ProjExchangeQ; ProjQ1=TargExchangeQ; 453 } 289 } 454 } else { << 290 } else // ProjQ2 is quark 455 if ( aProjQ1 > aProjQ2 ) { << 291 { 456 NewProjCode = aProjQ1*100 + aProjQ2* << 292 ProjExchangeQ = ProjQ2; 457 } else { << 293 if(ProjExchangeQ != TargQ1) 458 NewProjCode = aProjQ2*100 + aProjQ1* << 294 { >> 295 TargExchangeQ = TargQ1; TargQ1=ProjExchangeQ; ProjQ2=TargExchangeQ; >> 296 } else >> 297 if(ProjExchangeQ != TargQ2) >> 298 { >> 299 TargExchangeQ = TargQ2; TargQ2=ProjExchangeQ; ProjQ2=TargExchangeQ; >> 300 } else >> 301 { >> 302 TargExchangeQ = TargQ3; TargQ3=ProjExchangeQ; ProjQ2=TargExchangeQ; 459 } 303 } 460 } << 304 } // End of if(ProjQ1 > 0 ) // ProjQ1 is quark 461 #ifdef debugFTFexcitation << 462 G4cout << "NewProjCode " << NewProjCode << 463 #endif << 464 // Decide (with 50% probability) whether << 465 // but not in the case of charmed and bo << 466 // there are no excited charmed and bott << 467 ProjExcited = false; << 468 if ( aProjQ1 <= 3 && aProjQ2 <= 3 && << 469 NewProjCode += 2; // Excited meson (J << 470 ProjExcited = true; << 471 } << 472 305 473 // So far we have used the absolute valu << 306 G4int aProjQ1=std::abs(ProjQ1); 474 // now look at their signed values to se << 307 G4int aProjQ2=std::abs(ProjQ2); 475 G4int value = ProjQ1, absValue = aProjQ1 << 308 if(aProjQ1 == aProjQ2) {NewProjCode = 111;} // Pi0-meson 476 for ( G4int iQ = 0; iQ < 2; ++iQ ) { // << 309 else // |ProjQ1| # |ProjQ2| 477 if ( iQ == 1 ) { << 310 { 478 value = ProjQ2; absValue = aProjQ2; << 311 if(aProjQ1 > aProjQ2) {NewProjCode = aProjQ1*100+aProjQ2*10+1;} 479 } << 312 else {NewProjCode = aProjQ2*100+aProjQ1*10+1;} 480 if ( absValue == 2 || absValue == 4 << 313 NewProjCode *=(ProjectilePDGcode/absProjectilePDGcode); 481 Qquarks += 2*value/absValue; // u, << 314 } 482 } else { << 315 483 Qquarks -= value/absValue; // d, << 316 G4bool ProjExcited=false; >> 317 if(G4UniformRand() < 0.5) >> 318 { >> 319 NewProjCode +=2; // Excited Pi0-meson >> 320 ProjExcited=true; 484 } 321 } 485 } << 322 //G4cout<<G4endl<<"NewProjCode "<<NewProjCode<<G4endl; 486 // If Qquarks is positive, the sign of N << 487 // If Qquarks is negative, then the sign << 488 // If Qquarks is zero, then we need to d << 489 // 1. If aProjQ1 and aProjQ2 are the sam << 490 // (because the antiparticle is the s << 491 // 2. If aProjQ1 and aProjQ2 are not the << 492 // we have only two possibilities: << 493 // a. aProjQ1 and aProjQ2 are two dif << 494 // (s,d) or (b,d), or (b,s). In th << 495 // is fine (because the heaviest o << 496 // to be anti-quark belonging to t << 497 // implies a meson with positive P << 498 // b. aProjQ1 and aProjQ2 are two dif << 499 // The heaviest of the two (c) has << 500 // in the projectile, therefore th << 501 // reverse: 421 -> -421 anti_D0 (c << 502 if ( Qquarks < 0 || ( Qquarks == 0 && << 503 NewProjCode *= -1; << 504 } << 505 #ifdef debugFTFexcitation << 506 G4cout << "NewProjCode +2 or 0 " << NewP << 507 G4cout<<"+++++++++++++++++++++++++++++++ << 508 G4cout<<ProjQ1<<" "<<ProjQ2<<" "<<Qquark << 509 G4cout<<"+++++++++++++++++++++++++++++++ << 510 #endif << 511 << 512 // Proj << 513 TestParticle = G4ParticleTable::GetParti << 514 if ( ! TestParticle ) continue; << 515 common.MminProjectile = common.BrW.GetMi << 516 if ( common.SqrtS - common.M0target < co << 517 MtestPr = common.BrW.SampleMass( TestPar << 518 << 519 #ifdef debugFTFexcitation << 520 G4cout << "TestParticle Name " << NewPro << 521 << G4endl << 522 << "MtestPart MtestPart0 "<<Mtest << 523 << "M0projectile projectile PDGMa << 524 << projectile->GetDefinition()->G << 525 #endif << 526 << 527 // Targ << 528 NewTargCode = NewNucleonId( TargQ1, Targ << 529 #ifdef debugFTFexcitation << 530 G4cout << "New TrQ " << TargQ1 << " " << << 531 << "NewTargCode " << NewTargCode << 532 #endif << 533 << 534 // If the target has not heavy (charm or << 535 // see whether a Delta isobar can be cre << 536 if ( TargQ1 <= 3 && TargQ2 <= 3 && T << 537 if ( TargQ1 != TargQ2 && TargQ1 != T << 538 if ( G4UniformRand() < 0.5 ) { << 539 NewTargCode += 2; << 540 } else if ( G4UniformRand() < 0.75 ) << 541 NewTargCode = 3122; // Lambda << 542 } << 543 } else if ( TargQ1 == TargQ2 && Targ << 544 NewTargCode += 2; ProjExcited = true << 545 } else if ( target->GetDefinition()->G << 546 if ( G4UniformRand() > DeltaProbAtQu << 547 NewTargCode += 2; ProjExcited = tr << 548 } << 549 } else if ( ! ProjExcited && << 550 G4UniformRand() < DeltaPro << 551 common.SqrtS > common.M0pr << 552 G4ParticleTable::GetP << 553 NewTargCode += 2; // Create Delta i << 554 } << 555 } << 556 323 557 // Protection against: << 324 G4ParticleDefinition* TestParticle=0; 558 // - Lambda* (i.e. excited Lambda state) << 325 TestParticle=G4ParticleTable::GetParticleTable()->FindParticle(NewProjCode); 559 // - Sigma* (i.e. excited Sigma hyperon << 326 if(TestParticle) 560 // - Xi* (i.e. excited Xi hyperon st << 327 { 561 if ( NewTargCode == 3124 || // Lambda << 328 M0projectile= 562 NewTargCode == 3224 || // Sigma*+ NOT << 329 (G4ParticleTable::GetParticleTable()->FindParticle(NewProjCode))->GetPDGMass(); 563 NewTargCode == 3214 || // Sigma*0 NOT << 330 M0projectile2 = M0projectile * M0projectile; 564 NewTargCode == 3114 || // Sigma*- NOT << 331 565 NewTargCode == 3324 || // Xi*0 NOT << 332 ProjectileDiffStateMinMass =M0projectile+210.*MeV; //210 MeV=m_pi+70 MeV 566 NewTargCode == 3314 ) { // Xi*- NOT << 333 ProjectileNonDiffStateMinMass=M0projectile+210.*MeV; //210 MeV=m_pi+70 MeV 567 //G4cout << "G4DiffractiveExcitation::Excite << 334 } else 568 // << NewTargCode << G4endl; << 335 {return false;} 569 NewTargCode -= 2; // Corresponding ground-s << 336 570 } << 337 //G4cout<<"New TrQ "<<TargQ1<<" "<<TargQ2<<" "<<TargQ3<<G4endl; 571 << 338 NewTargCode = NewNucleonId(TargQ1, TargQ2, TargQ3); 572 // Special treatment for charmed and bot << 339 //G4cout<<"NewTargCode "<<NewTargCode<<G4endl; 573 // so we need to transform them by hand << 340 574 #ifdef debug_heavyHadrons << 341 if( (TargQ1 == TargQ2) && (TargQ1 == TargQ3) && 575 G4int initialNewTargCode = NewTargCode; << 342 (SqrtS > M0projectile+DeltaMass)) {NewTargCode +=2; //Create Delta isobar 576 #endif << 343 ProjExcited=true;} 577 if ( NewTargCode == 4322 ) NewTargC << 344 else if( target->GetDefinition()->GetPDGiIsospin() == 3 ) //Delta was the target 578 else if ( NewTargCode == 4312 ) NewTargC << 345 { if(G4UniformRand() > DeltaProbAtQuarkExchange){NewTargCode +=2; //Save Delta isobar 579 else if ( NewTargCode == 5312 ) NewTargC << 346 ProjExcited=true;} 580 else if ( NewTargCode == 5322 ) NewTargC << 347 else {} // De-excite initial Delta isobar 581 #ifdef debug_heavyHadrons << 348 } 582 if ( NewTargCode != initialNewTargCode ) << 349 583 G4cout << "G4DiffractiveExcitation::ExcitePa << 350 // else if((!CreateDelta) && 584 << "\t target heavy baryon with pdgCo << 351 else if((!ProjExcited) && 585 << " into pdgCode=" << NewTargCode << << 352 (G4UniformRand() < DeltaProbAtQuarkExchange) && //Nucleon was the target 586 } << 353 (SqrtS > M0projectile+DeltaMass)) {NewTargCode +=2;} //Create Delta isobar 587 #endif << 354 // else if( CreateDelta) {NewTargCode +=2;} 588 << 355 else {} //Save initial nucleon 589 TestParticle = G4ParticleTable::GetParti << 356 590 if ( ! TestParticle ) continue; << 357 // target->SetDefinition( // Fix 15.12.09 591 #ifdef debugFTFexcitation << 358 // G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode));// Fix 15.12.09 592 G4cout << "New targ " << NewTargCode << << 359 593 #endif << 360 //G4cout<<"NewTargCode "<<NewTargCode<<G4endl; 594 common.MminTarget = common.BrW.GetMinimu << 361 //G4int Uzhi; G4cin>>Uzhi; 595 if ( common.SqrtS - MtestPr < common.Mmi << 362 TestParticle=G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode); 596 MtestTr = common.BrW.SampleMass( TestPar << 363 if(TestParticle) 597 << 364 { 598 if ( common.SqrtS > MtestPr + MtestTr ) << 365 M0target= 599 << 366 (G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode))->GetPDGMass(); 600 } // End of while loop << 367 M0target2 = M0target * M0target; 601 if ( attempts >= maxNumberOfAttempts ) ret << 368 602 << 369 TargetDiffStateMinMass =M0target+220.*MeV; //220 MeV=m_pi+80 MeV; 603 if ( MtestPr >= common.Pprojectile.mag() << 370 TargetNonDiffStateMinMass=M0target+220.*MeV; //220 MeV=m_pi+80 MeV; 604 common.M0projectile = MtestPr; << 371 } else 605 } << 372 {return false;} 606 #ifdef debugFTFexcitation << 373 } else 607 G4cout << "M0projectile After check " << c << 374 { // projectile is baryon ---------------- 608 #endif << 375 //========================================================================= 609 common.M0projectile2 = common.M0projectile << 376 G4double Same=theParameters->GetProbOfSameQuarkExchange(); //0.3; //0.5; 0. 610 common.ProjectileDiffStateMinMass = com << 377 G4bool ProjDeltaHasCreated(false); 611 common.ProjectileNonDiffStateMinMass = com << 378 G4bool TargDeltaHasCreated(false); 612 if ( MtestTr >= common.Ptarget.mag() || << 379 613 common.M0target = MtestTr; << 380 G4double Ksi=G4UniformRand(); 614 } << 381 if(G4UniformRand() < 0.5) // Sampling exchange quark from proj. or targ. 615 common.M0target2 = common.M0target * commo << 382 { // Sampling exchanged quark from the projectile --- 616 #ifdef debugFTFexcitation << 383 if( Ksi < 0.333333 ) 617 G4cout << "New targ M0 M0^2 " << common.M0 << 384 {ProjExchangeQ = ProjQ1;} 618 #endif << 385 else if( (0.333333 <= Ksi) && (Ksi < 0.666667)) 619 common.TargetDiffStateMinMass = common. << 386 {ProjExchangeQ = ProjQ2;} 620 common.TargetNonDiffStateMinMass = common. << 387 else 621 << 388 {ProjExchangeQ = ProjQ3;} 622 } else { // of the if ( common.absProjectil << 389 623 << 390 //G4cout<<"ProjExchangeQ "<<ProjExchangeQ<<G4endl; 624 // If it is a projectile anti-baryon, no q << 391 if((ProjExchangeQ != TargQ1)||(G4UniformRand()<Same)) 625 // therefore returns immediately 1 (which << 392 { 626 // needs to be continued"). << 393 TargExchangeQ = TargQ1; TargQ1=ProjExchangeQ; ProjExchangeQ=TargExchangeQ; 627 if ( common.ProjectilePDGcode < 0 ) return << 394 } else 628 << 395 if((ProjExchangeQ != TargQ2)||(G4UniformRand()<Same)) 629 // Choose randomly, with equal probability << 396 { 630 // projectile or target hadron for selecti << 397 TargExchangeQ = TargQ2; TargQ2=ProjExchangeQ; ProjExchangeQ=TargExchangeQ; 631 G4bool isProjectileExchangedQ = false; << 398 } else 632 G4int firstQ = TargQ1, secondQ = << 399 { 633 G4int otherFirstQ = ProjQ1, otherSecondQ = << 400 TargExchangeQ = TargQ3; TargQ3=ProjExchangeQ; ProjExchangeQ=TargExchangeQ; 634 if ( G4UniformRand() < 0.5 ) { << 635 isProjectileExchangedQ = true; << 636 firstQ = ProjQ1; secondQ = Pro << 637 otherFirstQ = TargQ1; otherSecondQ = Tar << 638 } << 639 // Choose randomly, with equal probability << 640 // selected (projectile or target) hadron << 641 G4int exchangedQ = 0; << 642 G4double Ksi = G4UniformRand(); << 643 if ( Ksi < 0.333333 ) { << 644 exchangedQ = firstQ; << 645 } else if ( 0.333333 <= Ksi && Ksi < 0.6 << 646 exchangedQ = secondQ; << 647 } else { << 648 exchangedQ = thirdQ; << 649 } << 650 #ifdef debugFTFexcitation << 651 G4cout << "Exchange Qs isProjectile Q " << << 652 #endif << 653 // The exchanged quarks (one of the projec << 654 // are always accepted if they have differ << 655 // same flavour) they are accepted only wi << 656 G4double probSame = theParameters->GetProb << 657 const G4int MaxCount = 100; << 658 G4int count = 0, otherExchangedQ = 0; << 659 do { << 660 if ( exchangedQ != otherFirstQ || G4Un << 661 otherExchangedQ = otherFirstQ; otherFi << 662 } else { << 663 if ( exchangedQ != otherSecondQ || G << 664 otherExchangedQ = otherSecondQ; othe << 665 } else { << 666 if ( exchangedQ != otherThirdQ || << 667 otherExchangedQ = otherThirdQ; oth << 668 } << 669 } << 670 } << 671 } while ( otherExchangedQ == 0 && ++coun << 672 if ( count >= MaxCount ) return returnCode << 673 // Swap (between projectile and target had << 674 // as "exchanged" quarks. << 675 if ( Ksi < 0.333333 ) { << 676 firstQ = exchangedQ; << 677 } else if ( 0.333333 <= Ksi && Ksi < 0.6 << 678 secondQ = exchangedQ; << 679 } else { << 680 thirdQ = exchangedQ; << 681 } << 682 if ( isProjectileExchangedQ ) { << 683 ProjQ1 = firstQ; ProjQ2 = secondQ; << 684 TargQ1 = otherFirstQ; TargQ2 = otherSeco << 685 } else { << 686 TargQ1 = firstQ; TargQ2 = secondQ; << 687 ProjQ1 = otherFirstQ; ProjQ2 = otherSeco << 688 } << 689 #ifdef debugFTFexcitation << 690 G4cout << "Exchange Qs Pr Tr " << ( isPro << 691 << " " << ( isProjectileExchangedQ << 692 #endif << 693 << 694 NewProjCode = NewNucleonId( ProjQ1, ProjQ2 << 695 NewTargCode = NewNucleonId( TargQ1, TargQ2 << 696 // Decide whether the new projectile hadro << 697 // then decide whether the new target hadr << 698 // Notice that a Delta particle has the la << 699 // whereas a nucleon has "2" (because its << 700 // If the new projectile hadron or the new << 701 // constituent quark, then skip this part << 702 // excited charm and bottom hadrons). << 703 for ( G4int iHadron = 0; iHadron < 2; iHad << 704 // First projectile hadron, then target << 705 G4int codeQ1 = ProjQ1, codeQ2 = ProjQ2, << 706 G4double massConstraint = common.M0targe << 707 G4bool isHadronADelta = ( projectile->Ge << 708 if ( iHadron == 1 ) { // Target hadron << 709 codeQ1 = TargQ1, codeQ2 = TargQ2, code << 710 massConstraint = common.M0projectile; << 711 isHadronADelta = ( target->GetDefiniti << 712 } << 713 if ( codeQ1 > 3 || codeQ2 > 3 || cod << 714 if ( codeQ1 == codeQ2 && codeQ1 == cod << 715 newHadCode += 2; // Delta++ (uuu) or << 716 } else if ( isHadronADelta ) { // Hadro << 717 if ( G4UniformRand() > DeltaProbAtQuar << 718 newHadCode += 2; // Delta+ (uud) or << 719 } else { << 720 newHadCode += 0; // No delta (so th << 721 } << 722 } else { // Hadron (projectile or targe << 723 if ( G4UniformRand() < DeltaProbAtQuar << 724 common.SqrtS > G4ParticleTable::G << 725 + massConstraint ) << 726 newHadCode += 2; // Delta+ (uud) or << 727 } else { << 728 newHadCode += 0; // No delta (so th << 729 } 401 } 730 } << 731 if ( iHadron == 0 ) { // Projectile had << 732 NewProjCode = newHadCode; << 733 } else { // Target hadron << 734 NewTargCode = newHadCode; << 735 } << 736 } << 737 #ifdef debugFTFexcitation << 738 G4cout << "NewProjCode NewTargCode " << Ne << 739 #endif << 740 << 741 // Protection against: << 742 // - Lambda* (i.e. excited Lambda state) << 743 // - Sigma* (i.e. excited Sigma hyperon s << 744 // - Xi* (i.e. excited Xi hyperon stat << 745 if ( NewProjCode == 3124 || // Lambda* << 746 NewProjCode == 3224 || // Sigma*+ NOT ex << 747 NewProjCode == 3214 || // Sigma*0 NOT ex << 748 NewProjCode == 3114 || // Sigma*- NOT ex << 749 NewProjCode == 3324 || // Xi*0 NOT ex << 750 NewProjCode == 3314 ) { // Xi*- NOT ex << 751 //G4cout << "G4DiffractiveExcitation::Ex << 752 // << NewProjCode << G4endl; << 753 NewProjCode -= 2; // Corresponding grou << 754 } << 755 if ( NewTargCode == 3124 || // Lambda* << 756 NewTargCode == 3224 || // Sigma*+ NOT ex << 757 NewTargCode == 3214 || // Sigma*0 NOT ex << 758 NewTargCode == 3114 || // Sigma*- NOT ex << 759 NewTargCode == 3324 || // Xi*0 NOT ex << 760 NewTargCode == 3314 ) { // Xi*- NOT ex << 761 //G4cout << "G4DiffractiveExcitation::Ex << 762 // << NewTargCode << G4endl; << 763 NewTargCode -= 2; // Corresponding grou << 764 } << 765 402 766 // Special treatment for charmed and botto << 403 //G4cout<<"ProjExchangeQ "<<ProjExchangeQ<<G4endl; 767 // so we need to transform them by hand to << 404 //G4cout<<"TargExchangeQ "<<TargExchangeQ<<G4endl; 768 #ifdef debug_heavyHadrons << 405 if( Ksi < 0.333333 ) 769 G4int initialNewProjCode = NewProjCode, in << 406 {ProjQ1=ProjExchangeQ;} 770 #endif << 407 else if( (0.333333 <= Ksi) && (Ksi < 0.666667)) 771 if ( NewProjCode == 4322 ) NewProjCod << 408 {ProjQ2=ProjExchangeQ;} 772 else if ( NewProjCode == 4312 ) NewProjCod << 409 else 773 else if ( NewProjCode == 5312 ) NewProjCod << 410 {ProjQ3=ProjExchangeQ;} 774 else if ( NewProjCode == 5322 ) NewProjCod << 411 775 if ( NewTargCode == 4322 ) NewTargCod << 412 } else 776 else if ( NewTargCode == 4312 ) NewTargCod << 413 { // Sampling exchanged quark from the target ------- 777 else if ( NewTargCode == 5312 ) NewTargCod << 414 if( Ksi < 0.333333 ) 778 else if ( NewTargCode == 5322 ) NewTargCod << 415 {TargExchangeQ = TargQ1;} 779 #ifdef debug_heavyHadrons << 416 else if( (0.333333 <= Ksi) && (Ksi < 0.666667)) 780 if ( NewProjCode != initialNewProjCode || << 417 {TargExchangeQ = TargQ2;} 781 G4cout << "G4DiffractiveExcitation::Exci << 418 else 782 << "\t heavy baryon into an exist << 419 {TargExchangeQ = TargQ3;} 783 if ( NewProjCode != initialNewProjCode ) << 420 784 G4cout << "\t \t projectile baryon with pdgC << 421 if((TargExchangeQ != ProjQ1)||(G4UniformRand()<Same)) 785 << " into pdgCode=" << NewProjCode << << 422 { 786 } << 423 ProjExchangeQ = ProjQ1; ProjQ1=TargExchangeQ; TargExchangeQ=ProjExchangeQ; 787 if ( NewTargCode != initialNewTargCode ) << 424 } else 788 G4cout << "\t \t target baryon with pd << 425 if((TargExchangeQ != ProjQ2)||(G4UniformRand()<Same)) 789 << " into pdgCode=" << NewTargCode << << 426 { 790 } << 427 ProjExchangeQ = ProjQ2; ProjQ2=TargExchangeQ; TargExchangeQ=ProjExchangeQ; 791 } << 428 } else 792 #endif << 429 { 793 << 430 ProjExchangeQ = ProjQ3; ProjQ3=TargExchangeQ; TargExchangeQ=ProjExchangeQ; 794 // Sampling of the masses of the projectil << 795 // Because of energy conservation, the ord << 796 // randomly, half of the time we sample fi << 797 // then the projectile nucleon mass, and t << 798 // sample first the projectile nucleon mas << 799 G4VSplitableHadron* firstHadron = target; << 800 G4VSplitableHadron* secondHadron = project << 801 G4int firstHadronCode = NewTargCode, secon << 802 G4double massConstraint = common.M0project << 803 G4bool isFirstTarget = true; << 804 if ( G4UniformRand() < 0.5 ) { << 805 // Sample first the projectile nucleon m << 806 firstHadron = projectile; secondHa << 807 firstHadronCode = NewProjCode; secondHa << 808 massConstraint = common.M0target; << 809 isFirstTarget = false; << 810 } << 811 G4double Mtest1st = 0.0, Mtest2nd = 0.0, M << 812 for ( int iSamplingCase = 0; iSamplingCase << 813 G4VSplitableHadron* aHadron = firstHadro << 814 G4int aHadronCode = firstHadronCode; << 815 if ( iSamplingCase == 1 ) { // Second n << 816 aHadron = secondHadron; aHadronCode = << 817 } << 818 G4double MtestHadron = 0.0, MminHadron = << 819 if ( aHadron->GetStatus() == 1 || aHad << 820 TestParticle = G4ParticleTable::GetPar << 821 if ( ! TestParticle ) return returnCod << 822 MminHadron = common.BrW.GetMinimumMass << 823 if ( common.SqrtS - massConstraint < M << 824 if ( TestParticle->GetPDGWidth() == 0. << 825 MtestHadron = common.BrW.SampleMass( << 826 } else { << 827 const G4int maxNumberOfAttempts = 50 << 828 G4int attempts = 0; << 829 while ( attempts < maxNumberOfAttemp << 830 attempts++; << 831 MtestHadron = common.BrW.SampleMas << 832 << 833 if ( common.SqrtS < MtestHadron + << 834 continue; // Kinematically unac << 835 } else { << 836 break; // Kinematically acce << 837 } << 838 } << 839 if ( attempts >= maxNumberOfAttempts << 840 } 431 } 841 } << 842 if ( iSamplingCase == 0 ) { << 843 Mtest1st = MtestHadron; Mmin1st = Mmi << 844 } else { << 845 Mtest2nd = MtestHadron; Mmin2nd = Mmi << 846 } << 847 } // End for loop on the two sampling cas << 848 if ( isFirstTarget ) { << 849 MtestTr = Mtest1st; MtestPr = Mtest2n << 850 common.MminTarget = Mmin1st; common.Mmi << 851 } else { << 852 MtestTr = Mtest2nd; MtestPr = Mtest1s << 853 common.MminTarget = Mmin2nd; common.Mmi << 854 } << 855 << 856 if ( MtestPr != 0.0 ) { << 857 common.M0projectile = MtestPr; << 858 common.M0projectile2 = common.M0projecti << 859 common.ProjectileDiffStateMinMass = c << 860 common.ProjectileNonDiffStateMinMass = c << 861 } << 862 if ( MtestTr != 0.0 ) { << 863 common.M0target = MtestTr; << 864 common.M0target2 = common.M0target * com << 865 common.TargetDiffStateMinMass = commo << 866 common.TargetNonDiffStateMinMass = commo << 867 } << 868 432 869 } // End of if ( common.absProjectilePDGcod << 433 if( Ksi < 0.333333 ) >> 434 {TargQ1=TargExchangeQ;} >> 435 else if( (0.333333 <= Ksi) && (Ksi < 0.666667)) >> 436 {TargQ2=TargExchangeQ;} >> 437 else >> 438 {TargQ3=TargExchangeQ;} >> 439 >> 440 } // End of sampling baryon >> 441 >> 442 NewProjCode = NewNucleonId(ProjQ1, ProjQ2, ProjQ3); // ***************************** >> 443 >> 444 //G4cout<<"ProjQ1, ProjQ2, ProjQ3 "<<ProjQ1<<" "<<ProjQ2<<" "<<ProjQ3<<" "<<NewProjCode<<G4endl; >> 445 >> 446 G4int TestParticleID=NewProjCode; >> 447 G4ParticleDefinition* TestParticle=0; >> 448 G4double TestParticleMass=DBL_MAX; >> 449 >> 450 TestParticle=G4ParticleTable::GetParticleTable()->FindParticle(NewProjCode); >> 451 if(TestParticle) TestParticleMass=TestParticle->GetPDGMass(); >> 452 >> 453 if((ProjQ1==ProjQ2) && (ProjQ1==ProjQ3)) {NewProjCode +=2; ProjDeltaHasCreated=true;} >> 454 else if(projectile->GetDefinition()->GetPDGiIsospin() == 3)// Projectile was Delta >> 455 { if(G4UniformRand() > DeltaProbAtQuarkExchange) >> 456 {NewProjCode +=2; ProjDeltaHasCreated=true;} >> 457 else {NewProjCode +=0; ProjDeltaHasCreated=false;} >> 458 } >> 459 else // Projectile was Nucleon >> 460 { >> 461 if((G4UniformRand() < DeltaProbAtQuarkExchange) && (SqrtS > DeltaMass+M0target)) >> 462 {NewProjCode +=2; ProjDeltaHasCreated=true;} >> 463 else {NewProjCode +=0; ProjDeltaHasCreated=false;} >> 464 } >> 465 >> 466 G4ParticleDefinition* NewTestParticle= >> 467 G4ParticleTable::GetParticleTable()->FindParticle(NewProjCode); >> 468 //G4cout<<"TestParticleMass NewTestParticle->GetPDGMass() "<<TestParticleMass<<" "<< NewTestParticle->GetPDGMass()<<G4endl; >> 469 //if(TestParticleMass < NewTestParticle->GetPDGMass()) {NewProjCode=TestParticleID;} >> 470 >> 471 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++= 870 472 871 // If we assume that final state hadrons aft << 473 NewTargCode = NewNucleonId(TargQ1, TargQ2, TargQ3); // ***************************** 872 // in the ground states, we have to put << 873 if ( common.SqrtS < common.M0projectile + co << 874 << 875 common.PZcms2 = ( sqr( common.S ) + sqr( com << 876 - 2.0 * ( common.S * ( com << 877 + common.M0proje << 878 #ifdef debugFTFexcitation << 879 G4cout << "At the end// NewProjCode " << New << 880 << "At the end// NewTargCode " << New << 881 << "M0pr M0tr SqS " << common.M0pro << 882 << common.SqrtS << G4endl << 883 << "M0pr2 M0tr2 SqS " << common.M0pro << 884 << common.SqrtS << G4endl << 885 << "PZcms2 after the change " << comm << 886 #endif << 887 if ( common.PZcms2 < 0.0 ) return returnCode << 888 << 889 projectile->SetDefinition( G4ParticleTable:: << 890 target->SetDefinition( G4ParticleTable::GetP << 891 common.PZcms = std::sqrt( common.PZcms2 ); << 892 common.Pprojectile.setPz( common.PZcms ); << 893 common.Pprojectile.setE( std::sqrt( common.M << 894 common.Ptarget.setPz( -common.PZcms ); << 895 common.Ptarget.setE( std::sqrt( common.M0tar << 896 #ifdef debugFTFexcitation << 897 G4cout << "Proj Targ and Proj+Targ in CMS" < << 898 << common.Ptarget << G4endl << common << 899 #endif << 900 << 901 if ( projectile->GetStatus() != 0 ) projecti << 902 if ( target->GetStatus() != 0 ) target-> << 903 << 904 // Check for possible excitation of the part << 905 if ( common.SqrtS < common.M0projectile + co << 906 common.SqrtS < common.ProjectileDiffSta << 907 common.ProbOfDiffraction == 0.0 ) commo << 908 << 909 if ( G4UniformRand() > common.ProbExc ) { / << 910 #ifdef debugFTFexcitation << 911 G4cout << "Make elastic scattering of new << 912 #endif << 913 common.Pprojectile.transform( common.toLab << 914 common.Ptarget.transform( common.toLab ); << 915 projectile->Set4Momentum( common.Pprojecti << 916 target->Set4Momentum( common.Ptarget ); << 917 G4bool Result = theElastic->ElasticScatter << 918 #ifdef debugFTFexcitation << 919 G4cout << "Result of el. scatt " << Result << 920 << G4endl << projectile->Get4Moment << 921 << projectile->Get4Momentum() + tar << 922 << (projectile->Get4Momentum() + ta << 923 #endif << 924 if ( Result ) returnCode = 0; // successf << 925 return returnCode; << 926 } << 927 << 928 #ifdef debugFTFexcitation << 929 G4cout << "Make excitation of new hadrons" < << 930 #endif << 931 << 932 // Redefinition of ProbOfDiffraction because << 933 common.ProbOfDiffraction = common.ProbProjec << 934 if ( common.ProbOfDiffraction != 0.0 ) { << 935 common.ProbProjectileDiffraction /= common << 936 common.ProbTargetDiffraction /= common << 937 } << 938 474 939 return returnCode = 1; // successfully comp << 475 //G4cout<<"TargQ1, TargQ2, TargQ3 "<<TargQ1<<" "<<TargQ2<<" "<<TargQ3<<" "<<NewTargCode<<G4endl; 940 } << 941 476 942 //-------------------------------------------- << 477 TestParticleID=NewTargCode; 943 << 478 TestParticleMass=DBL_MAX; 944 G4bool G4DiffractiveExcitation:: << 945 ExciteParticipants_doDiffraction( G4VSplitable << 946 G4FTFParamet << 947 G4Diffractiv << 948 // Second of the three utility methods used << 949 // it does the sampling for the diffraction << 950 << 951 G4bool isProjectileDiffraction = false; << 952 if ( G4UniformRand() < common.ProbProjectile << 953 isProjectileDiffraction = true; << 954 #ifdef debugFTFexcitation << 955 G4cout << "projectile diffraction" << G4en << 956 #endif << 957 common.ProjMassT2 = common.ProjectileDiffS << 958 common.ProjMassT = common.ProjectileDiffS << 959 common.TargMassT2 = common.M0target2; << 960 common.TargMassT = common.M0target; << 961 } else { << 962 #ifdef debugFTFexcitation << 963 G4cout << "Target diffraction" << G4endl; << 964 #endif << 965 common.ProjMassT2 = common.M0projectile2; << 966 common.ProjMassT = common.M0projectile; << 967 common.TargMassT2 = common.TargetDiffState << 968 common.TargMassT = common.TargetDiffState << 969 } << 970 << 971 // Check whether the interaction is possible << 972 if ( common.SqrtS < common.ProjMassT + commo << 973 << 974 common.PZcms2 = ( sqr( common.S ) + sqr( com << 975 - 2.0 * ( common.S * ( com << 976 + common.ProjMassT2 << 977 if ( common.PZcms2 < 0.0 ) return false; << 978 common.maxPtSquare = common.PZcms2; << 979 << 980 G4double DiffrAveragePt2 = theParameters->Ge << 981 G4bool loopCondition = true; << 982 G4int whilecount = 0; << 983 do { // Generate pt and mass of projectile << 984 << 985 whilecount++; << 986 if ( whilecount > 1000 ) { << 987 common.Qmomentum = G4LorentzVector( 0.0, << 988 return false; // Ignore this interacti << 989 }; << 990 << 991 common.Qmomentum = G4LorentzVector( Gaussi << 992 common.Pt2 = G4ThreeVector( common.Qmoment << 993 if ( isProjectileDiffraction ) { // proje << 994 common.ProjMassT2 = common.ProjectileDif << 995 common.TargMassT2 = common.M0target2 + c << 996 } else { // targe << 997 common.ProjMassT2 = common.M0projectile2 << 998 common.TargMassT2 = common.TargetDiffSta << 999 } << 1000 common.ProjMassT = std::sqrt( common.Proj << 1001 common.TargMassT = std::sqrt( common.Targ << 1002 if ( common.SqrtS < common.ProjMassT + co << 1003 << 1004 common.PZcms2 = ( sqr( common.S ) + sqr( << 1005 - 2.0 * ( common.S * ( << 1006 + common.Proj << 1007 if ( common.PZcms2 < 0.0 ) continue; << 1008 << 1009 common.PZcms = std::sqrt( common.PZcms2 ) << 1010 if ( isProjectileDiffraction ) { // proj << 1011 common.PMinusMin = std::sqrt( common.Pr << 1012 common.PMinusMax = common.SqrtS - commo << 1013 common.PMinusNew = ChooseP( common.PMin << 1014 common.TMinusNew = common.SqrtS - commo << 1015 common.Qminus = common.Ptarget.minus() << 1016 common.TPlusNew = common.TargMassT2 / c << 1017 common.Qplus = common.Ptarget.plus() - << 1018 common.Qmomentum.setPz( (common.Qplus - << 1019 common.Qmomentum.setE( (common.Qplus + << 1020 loopCondition = ( ( common.Pprojectile << 1021 common.ProjectileDiff << 1022 } else { // targ << 1023 common.TPlusMin = std::sqrt( common.Tar << 1024 common.TPlusMax = common.SqrtS - common << 1025 common.TPlusNew = ChooseP( common.TPlus << 1026 common.PPlusNew = common.SqrtS - common << 1027 common.Qplus = common.PPlusNew - common << 1028 common.PMinusNew = common.ProjMassT2 / << 1029 common.Qminus = common.PMinusNew - comm << 1030 common.Qmomentum.setPz( (common.Qplus - << 1031 common.Qmomentum.setE( (common.Qplus + << 1032 loopCondition = ( ( common.Ptarget - c << 1033 common.TargetDiffSta << 1034 } << 1035 479 1036 } while ( loopCondition ); /* Loop checkin << 480 TestParticle=G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode); 1037 // Repeat the sampling because ther << 481 if(TestParticle) TestParticleMass=TestParticle->GetPDGMass(); >> 482 >> 483 if((TargQ1==TargQ2) && (TargQ1==TargQ3)) {NewTargCode +=2; TargDeltaHasCreated=true;} >> 484 else if(target->GetDefinition()->GetPDGiIsospin() == 3) // Target was Delta >> 485 { if(G4UniformRand() > DeltaProbAtQuarkExchange) >> 486 {NewTargCode +=2; TargDeltaHasCreated=true;} >> 487 else {NewTargCode +=0; TargDeltaHasCreated=false;} >> 488 } >> 489 else // Target was Nucleon >> 490 { >> 491 if((G4UniformRand() < DeltaProbAtQuarkExchange) && (SqrtS > M0projectile+DeltaMass)) >> 492 {NewTargCode +=2; TargDeltaHasCreated=true;} >> 493 else {NewTargCode +=0; TargDeltaHasCreated=false;} >> 494 } >> 495 >> 496 NewTestParticle=G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode); >> 497 //G4cout<<"TestParticleMass NewTestParticle->GetPDGMass() "<<TestParticleMass<<" "<< NewTestParticle->GetPDGMass()<<G4endl; >> 498 //if(TestParticleMass < NewTestParticle->GetPDGMass()) {NewTargCode=TestParticleID;} >> 499 >> 500 //G4cout<<"NewProjCode NewTargCode "<<NewProjCode<<" "<<NewTargCode<<G4endl; >> 501 //G4int Uzhi; G4cin>>Uzhi; >> 502 >> 503 if((absProjectilePDGcode == NewProjCode) && (absTargetPDGcode == NewTargCode)) >> 504 { // Nothing was changed! It is not right!? >> 505 } >> 506 // Forming baryons -------------------------------------------------- >> 507 if(ProjDeltaHasCreated) {ProbProjectileDiffraction=1.; ProbTargetDiffraction=0.;} >> 508 if(TargDeltaHasCreated) {ProbProjectileDiffraction=0.; ProbTargetDiffraction=1.;} >> 509 if(ProjDeltaHasCreated) >> 510 { >> 511 M0projectile= >> 512 (G4ParticleTable::GetParticleTable()->FindParticle(NewProjCode))->GetPDGMass(); >> 513 M0projectile2 = M0projectile * M0projectile; >> 514 >> 515 ProjectileDiffStateMinMass =M0projectile+210.*MeV; //210 MeV=m_pi+70 MeV >> 516 ProjectileNonDiffStateMinMass=M0projectile+210.*MeV; //210 MeV=m_pi+70 MeV >> 517 } >> 518 >> 519 // if(M0target < >> 520 // (G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode))->GetPDGMass()) >> 521 if(TargDeltaHasCreated) >> 522 { >> 523 M0target= >> 524 (G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode))->GetPDGMass(); >> 525 M0target2 = M0target * M0target; >> 526 >> 527 TargetDiffStateMinMass =M0target+210.*MeV; //210 MeV=m_pi+70 MeV; >> 528 TargetNonDiffStateMinMass=M0target+210.*MeV; //210 MeV=m_pi+70 MeV; >> 529 } >> 530 } // End of if projectile is baryon --------------------------- >> 531 >> 532 //G4cout<<"At end// NewProjCode "<<NewProjCode<<G4endl; >> 533 //G4cout<<"At end// NewTargCode "<<NewTargCode<<G4endl; >> 534 >> 535 // If we assume that final state hadrons after the charge exchange will be >> 536 // in the ground states, we have to put ---------------------------------- >> 537 //G4cout<<"M0pr M0tr SqS "<<M0projectile<<" "<<M0target<<" "<<SqrtS<<G4endl; >> 538 >> 539 PZcms2=(S*S+M0projectile2*M0projectile2+M0target2*M0target2- >> 540 2*S*M0projectile2 - 2*S*M0target2 - 2*M0projectile2*M0target2) >> 541 /4./S; >> 542 //G4cout<<"PZcms2 1 "<<PZcms2<<G4endl<<G4endl; >> 543 if(PZcms2 < 0) {return false;} // It can be if energy is not sufficient for Delta >> 544 //---------------------------------------------------------- >> 545 projectile->SetDefinition( >> 546 G4ParticleTable::GetParticleTable()->FindParticle(NewProjCode)); >> 547 >> 548 target->SetDefinition( >> 549 G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode)); >> 550 //---------------------------------------------------------- >> 551 PZcms = std::sqrt(PZcms2); >> 552 >> 553 Pprojectile.setPz( PZcms); >> 554 Pprojectile.setE(std::sqrt(M0projectile2+PZcms2)); >> 555 >> 556 Ptarget.setPz( -PZcms); >> 557 Ptarget.setE(std::sqrt(M0target2+PZcms2)); >> 558 >> 559 // ---------------------------------------------------------- >> 560 >> 561 if(absProjectilePDGcode < 1000) >> 562 { // For projectile meson >> 563 G4double Wexcit=1.-2.256*std::exp(-0.6*ProjectileRapidity); >> 564 Wexcit=0.; >> 565 if(G4UniformRand() > Wexcit) >> 566 { // Make elastic scattering >> 567 //G4cout<<"Make elastic scattering of new hadrons"<<G4endl; >> 568 Pprojectile.transform(toLab); >> 569 Ptarget.transform(toLab); >> 570 >> 571 projectile->SetTimeOfCreation(target->GetTimeOfCreation()); >> 572 projectile->SetPosition(target->GetPosition()); >> 573 >> 574 projectile->Set4Momentum(Pprojectile); >> 575 target->Set4Momentum(Ptarget); >> 576 >> 577 G4bool Result= theElastic->ElasticScattering (projectile,target,theParameters); >> 578 return Result; >> 579 } // end of if(Make elastic scattering for projectile meson?) >> 580 } else >> 581 { // For projectile baryon >> 582 G4double Wexcit=1.-2.256*std::exp(-0.6*ProjectileRapidity); >> 583 //Wexcit=0.; >> 584 if(G4UniformRand() > Wexcit) >> 585 { // Make elastic scattering >> 586 //G4cout<<"Make elastic scattering of new hadrons"<<G4endl; >> 587 Pprojectile.transform(toLab); >> 588 Ptarget.transform(toLab); >> 589 >> 590 projectile->SetTimeOfCreation(target->GetTimeOfCreation()); >> 591 projectile->SetPosition(target->GetPosition()); >> 592 >> 593 projectile->Set4Momentum(Pprojectile); >> 594 target->Set4Momentum(Ptarget); >> 595 >> 596 G4bool Result= theElastic->ElasticScattering (projectile,target,theParameters); >> 597 return Result; >> 598 } // end of if(Make elastic scattering for projectile baryon?) >> 599 } >> 600 //G4cout<<"Make excitation of new hadrons"<<G4endl; >> 601 } // End of charge exchange part ------------------------------ >> 602 >> 603 // ------------------------------------------------------------------ >> 604 G4double ProbOfDiffraction=ProbProjectileDiffraction + ProbTargetDiffraction; >> 605 /* >> 606 G4cout<<"Excite --------------------"<<G4endl; >> 607 G4cout<<"Proj "<<M0projectile<<" "<<ProjectileDiffStateMinMass<<" "<<ProjectileNonDiffStateMinMass<<G4endl; >> 608 G4cout<<"Targ "<<M0target <<" "<<TargetDiffStateMinMass <<" "<<TargetNonDiffStateMinMass<<G4endl; >> 609 G4cout<<"SqrtS "<<SqrtS<<G4endl; >> 610 >> 611 G4cout<<"Prob ProjDiff TargDiff "<<ProbProjectileDiffraction<<" "<<ProbTargetDiffraction<<" "<<ProbOfDiffraction<<G4endl; >> 612 G4cout<<"Pr Y "<<Pprojectile.rapidity()<<" Tr Y "<<Ptarget.rapidity()<<G4endl; >> 613 G4int Uzhi; G4cin>>Uzhi; >> 614 */ >> 615 /* >> 616 if(ProjectileNonDiffStateMinMass + TargetNonDiffStateMinMass > SqrtS) // 24.07.10 >> 617 { >> 618 if(ProbOfDiffraction!=0.) >> 619 { >> 620 ProbProjectileDiffraction/=ProbOfDiffraction; >> 621 ProbOfDiffraction=1.; >> 622 } else {return false;} >> 623 } >> 624 >> 625 */ >> 626 >> 627 if(ProbOfDiffraction!=0.) >> 628 { >> 629 ProbProjectileDiffraction/=ProbOfDiffraction; >> 630 } >> 631 else >> 632 { >> 633 ProbProjectileDiffraction=0.; >> 634 } >> 635 >> 636 //G4cout<<"Prob ProjDiff TargDiff "<<ProbProjectileDiffraction<<" "<<ProbTargetDiffraction<<" "<<ProbOfDiffraction<<G4endl; >> 637 >> 638 G4double ProjectileDiffStateMinMass2 = ProjectileDiffStateMinMass * >> 639 ProjectileDiffStateMinMass; >> 640 G4double ProjectileNonDiffStateMinMass2 = ProjectileNonDiffStateMinMass * >> 641 ProjectileNonDiffStateMinMass; >> 642 >> 643 G4double TargetDiffStateMinMass2 = TargetDiffStateMinMass * >> 644 TargetDiffStateMinMass; >> 645 G4double TargetNonDiffStateMinMass2 = TargetNonDiffStateMinMass * >> 646 TargetNonDiffStateMinMass; >> 647 >> 648 G4double Pt2; >> 649 G4double ProjMassT2, ProjMassT; >> 650 G4double TargMassT2, TargMassT; >> 651 G4double PMinusMin, PMinusMax; >> 652 // G4double PPlusMin , PPlusMax; >> 653 G4double TPlusMin , TPlusMax; >> 654 G4double PMinusNew, PPlusNew, TPlusNew, TMinusNew; >> 655 >> 656 G4LorentzVector Qmomentum; >> 657 G4double Qminus, Qplus; >> 658 >> 659 G4int whilecount=0; >> 660 >> 661 // Choose a process --------------------------- >> 662 >> 663 if(G4UniformRand() < ProbOfDiffraction) >> 664 { >> 665 if(G4UniformRand() < ProbProjectileDiffraction) >> 666 { //-------- projectile diffraction --------------- >> 667 //G4cout<<"projectile diffraction"<<G4endl; >> 668 >> 669 do { >> 670 // Generate pt >> 671 // if (whilecount++ >= 500 && (whilecount%100)==0) >> 672 // G4cout << "G4DiffractiveExcitation::ExciteParticipants possibly looping" >> 673 // << ", loop count/ maxPtSquare : " >> 674 // << whilecount << " / " << maxPtSquare << G4endl; >> 675 >> 676 // whilecount++; >> 677 if (whilecount > 1000 ) >> 678 { >> 679 Qmomentum=G4LorentzVector(0.,0.,0.,0.); >> 680 target->SetStatus(2); return false; // Ignore this interaction >> 681 }; >> 682 >> 683 // --------------- Check that the interaction is possible ----------- >> 684 ProjMassT2=ProjectileDiffStateMinMass2; >> 685 ProjMassT =ProjectileDiffStateMinMass; >> 686 >> 687 TargMassT2=M0target2; >> 688 TargMassT =M0target; >> 689 //G4cout<<"Masses "<<ProjMassT<<" "<<TargMassT<<" "<<SqrtS<<" "<<ProjMassT+TargMassT<<G4endl; >> 690 PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2- >> 691 2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2) >> 692 /4./S; >> 693 >> 694 //G4cout<<"PZcms2 PrD"<<PZcms2<<G4endl; >> 695 if(PZcms2 < 0 ) >> 696 { >> 697 target->SetStatus(2); >> 698 return false; >> 699 } >> 700 maxPtSquare=PZcms2; >> 701 >> 702 Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0); >> 703 Pt2=G4ThreeVector(Qmomentum.vect()).mag2(); >> 704 >> 705 ProjMassT2=ProjectileDiffStateMinMass2+Pt2; >> 706 ProjMassT =std::sqrt(ProjMassT2); >> 707 >> 708 TargMassT2=M0target2+Pt2; >> 709 TargMassT =std::sqrt(TargMassT2); >> 710 >> 711 PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2- >> 712 2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2) >> 713 /4./S; >> 714 >> 715 if(PZcms2 < 0 ) continue; >> 716 PZcms =std::sqrt(PZcms2); >> 717 >> 718 PMinusMin=std::sqrt(ProjMassT2+PZcms2)-PZcms; >> 719 PMinusMax=SqrtS-TargMassT; >> 720 >> 721 PMinusNew=ChooseP(PMinusMin, PMinusMax); >> 722 // PMinusNew=1./sqrt(1./PMinusMin-G4UniformRand()*(1./PMinusMin-1./PMinusMax)); >> 723 //PMinusNew=1./sqr(1./std::sqrt(PMinusMin)-G4UniformRand()*(1./std::sqrt(PMinusMin)-1./std::sqrt(PMinusMax))); >> 724 >> 725 TMinusNew=SqrtS-PMinusNew; >> 726 Qminus=Ptarget.minus()-TMinusNew; >> 727 TPlusNew=TargMassT2/TMinusNew; >> 728 Qplus=Ptarget.plus()-TPlusNew; >> 729 >> 730 Qmomentum.setPz( (Qplus-Qminus)/2 ); >> 731 Qmomentum.setE( (Qplus+Qminus)/2 ); >> 732 >> 733 } while ((Pprojectile+Qmomentum).mag2() < ProjectileDiffStateMinMass2); //|| >> 734 //Repeat the sampling because there was not any excitation >> 735 //((Ptarget -Qmomentum).mag2() < M0target2 )) ); >> 736 } >> 737 else >> 738 { // -------------- Target diffraction ---------------- 1038 739 1039 if ( isProjectileDiffraction ) { // projec << 740 //G4cout<<"Target diffraction"<<G4endl; 1040 projectile->SetStatus( 0 ); << 741 do { 1041 if ( projectile->GetStatus() == 2 ) proje << 742 // Generate pt 1042 if ( target->GetStatus() == 1 && target << 743 // if (whilecount++ >= 500 && (whilecount%100)==0) 1043 } else { // target << 744 // G4cout << "G4DiffractiveExcitation::ExciteParticipants possibly looping" 1044 target->SetStatus( 0 ); << 745 // << ", loop count/ maxPtSquare : " 1045 } << 746 // << whilecount << " / " << maxPtSquare << G4endl; >> 747 >> 748 // whilecount++; >> 749 if (whilecount > 1000 ) >> 750 { >> 751 Qmomentum=G4LorentzVector(0.,0.,0.,0.); >> 752 target->SetStatus(2); return false; // Ignore this interaction >> 753 }; >> 754 //G4cout<<"Qm while "<<Qmomentum<<" "<<whilecount<<G4endl; >> 755 // --------------- Check that the interaction is possible ----------- >> 756 ProjMassT2=M0projectile2; >> 757 ProjMassT =M0projectile; >> 758 >> 759 TargMassT2=TargetDiffStateMinMass2; >> 760 TargMassT =TargetDiffStateMinMass; >> 761 >> 762 PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2- >> 763 2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2) >> 764 /4./S; >> 765 >> 766 //G4cout<<"PZcms2 TrD <0 "<<PZcms2<<" return"<<G4endl; >> 767 if(PZcms2 < 0 ) >> 768 { >> 769 target->SetStatus(2); >> 770 return false; >> 771 } >> 772 maxPtSquare=PZcms2; >> 773 >> 774 Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0); >> 775 >> 776 //G4cout<<"Qm while "<<Qmomentum<<" "<<whilecount<<G4endl; >> 777 Pt2=G4ThreeVector(Qmomentum.vect()).mag2(); >> 778 >> 779 ProjMassT2=M0projectile2+Pt2; >> 780 ProjMassT =std::sqrt(ProjMassT2); >> 781 >> 782 TargMassT2=TargetDiffStateMinMass2+Pt2; >> 783 TargMassT =std::sqrt(TargMassT2); >> 784 >> 785 PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2- >> 786 2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2) >> 787 /4./S; >> 788 >> 789 //G4cout<<"PZcms2 <0 "<<PZcms2<<" continue"<<G4endl; >> 790 if(PZcms2 < 0 ) continue; >> 791 PZcms =std::sqrt(PZcms2); >> 792 >> 793 TPlusMin=std::sqrt(TargMassT2+PZcms2)-PZcms; >> 794 TPlusMax=SqrtS-ProjMassT; >> 795 >> 796 TPlusNew=ChooseP(TPlusMin, TPlusMax); >> 797 //TPlusNew=1./sqr(1./std::sqrt(TPlusMin)-G4UniformRand()*(1./std::sqrt(TPlusMin)-1./std::sqrt(TPlusMax))); >> 798 >> 799 //TPlusNew=TPlusMax; >> 800 >> 801 PPlusNew=SqrtS-TPlusNew; >> 802 Qplus=PPlusNew-Pprojectile.plus(); >> 803 PMinusNew=ProjMassT2/PPlusNew; >> 804 Qminus=PMinusNew-Pprojectile.minus(); >> 805 >> 806 Qmomentum.setPz( (Qplus-Qminus)/2 ); >> 807 Qmomentum.setE( (Qplus+Qminus)/2 ); >> 808 >> 809 /* >> 810 G4cout<<(Pprojectile+Qmomentum).mag()<<" "<<M0projectile<<G4endl; >> 811 G4bool First=(Pprojectile+Qmomentum).mag2() < M0projectile2; >> 812 G4cout<<First<<G4endl; >> 813 >> 814 G4cout<<(Ptarget -Qmomentum).mag()<<" "<<TargetDiffStateMinMass<<" "<<TargetDiffStateMinMass2<<G4endl; >> 815 G4bool Seco=(Ptarget -Qmomentum).mag2() < TargetDiffStateMinMass2; >> 816 G4cout<<Seco<<G4endl; >> 817 */ >> 818 >> 819 } while ((Ptarget -Qmomentum).mag2() < TargetDiffStateMinMass2); >> 820 // Repeat the sampling because there was not any excitation >> 821 // (((Pprojectile+Qmomentum).mag2() < M0projectile2 ) || //No without excitation >> 822 // ((Ptarget -Qmomentum).mag2() < TargetDiffStateMinMass2)) ); >> 823 //G4cout<<"Go out"<<G4endl; >> 824 } // End of if(G4UniformRand() < ProbProjectileDiffraction) >> 825 } >> 826 else //----------- Non-diffraction process ------------ >> 827 { 1046 828 1047 return true; << 829 //G4cout<<"Non-diffraction process"<<G4endl; 1048 } << 830 do { 1049 << 831 // Generate pt 1050 //------------------------------------------- << 832 // if (whilecount++ >= 500 && (whilecount%100)==0) >> 833 // G4cout << "G4DiffractiveExcitation::ExciteParticipants possibly looping" >> 834 // << ", loop count/ maxPtSquare : " >> 835 // << whilecount << " / " << maxPtSquare << G4endl; >> 836 >> 837 // whilecount++; >> 838 if (whilecount > 1000 ) >> 839 { >> 840 Qmomentum=G4LorentzVector(0.,0.,0.,0.); >> 841 target->SetStatus(2); return false; // Ignore this interaction >> 842 }; >> 843 // --------------- Check that the interaction is possible ----------- >> 844 ProjMassT2=ProjectileNonDiffStateMinMass2; >> 845 ProjMassT =ProjectileNonDiffStateMinMass; >> 846 >> 847 TargMassT2=TargetNonDiffStateMinMass2; >> 848 TargMassT =TargetNonDiffStateMinMass; >> 849 >> 850 PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2- >> 851 2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2) >> 852 /4./S; >> 853 >> 854 if(PZcms2 < 0 ) >> 855 { >> 856 target->SetStatus(2); >> 857 return false; >> 858 } >> 859 maxPtSquare=PZcms2; >> 860 Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0); >> 861 Pt2=G4ThreeVector(Qmomentum.vect()).mag2(); >> 862 >> 863 ProjMassT2=ProjectileNonDiffStateMinMass2+Pt2; >> 864 ProjMassT =std::sqrt(ProjMassT2); >> 865 >> 866 TargMassT2=TargetNonDiffStateMinMass2+Pt2; >> 867 TargMassT =std::sqrt(TargMassT2); >> 868 >> 869 PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2- >> 870 2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2) >> 871 /4./S; >> 872 //G4cout<<"PZcms2 ND"<<PZcms2<<G4endl; >> 873 >> 874 if(PZcms2 < 0 ) continue; >> 875 PZcms =std::sqrt(PZcms2); >> 876 >> 877 PMinusMin=std::sqrt(ProjMassT2+PZcms2)-PZcms; >> 878 PMinusMax=SqrtS-TargMassT; >> 879 >> 880 PMinusNew=ChooseP(PMinusMin, PMinusMax); >> 881 >> 882 Qminus=PMinusNew-Pprojectile.minus(); >> 883 >> 884 TPlusMin=std::sqrt(TargMassT2+PZcms2)-PZcms; >> 885 // TPlusMax=SqrtS-PMinusNew; >> 886 TPlusMax=SqrtS-ProjMassT; >> 887 >> 888 TPlusNew=ChooseP(TPlusMin, TPlusMax); >> 889 >> 890 Qplus=-(TPlusNew-Ptarget.plus()); >> 891 >> 892 Qmomentum.setPz( (Qplus-Qminus)/2 ); >> 893 Qmomentum.setE( (Qplus+Qminus)/2 ); >> 894 /* >> 895 G4cout<<(Pprojectile+Qmomentum).mag2()<<" "<<ProjectileNonDiffStateMinMass2<<G4endl; >> 896 G4cout<<(Ptarget -Qmomentum).mag2()<<" "<<TargetNonDiffStateMinMass2<<G4endl; >> 897 G4int Uzhi; G4cin>>Uzhi; >> 898 */ >> 899 } while ( >> 900 ((Pprojectile+Qmomentum).mag2() < ProjectileNonDiffStateMinMass2) || //No double Diffraction >> 901 ((Ptarget -Qmomentum).mag2() < TargetNonDiffStateMinMass2 )); >> 902 } >> 903 >> 904 Pprojectile += Qmomentum; >> 905 Ptarget -= Qmomentum; >> 906 >> 907 //G4cout<<"Pr Y "<<Pprojectile.rapidity()<<" Tr Y "<<Ptarget.rapidity()<<G4endl; >> 908 >> 909 // Transform back and update SplitableHadron Participant. >> 910 Pprojectile.transform(toLab); >> 911 Ptarget.transform(toLab); >> 912 >> 913 // Calculation of the creation time --------------------- >> 914 projectile->SetTimeOfCreation(target->GetTimeOfCreation()); >> 915 projectile->SetPosition(target->GetPosition()); >> 916 // Creation time and position of target nucleon were determined at >> 917 // ReggeonCascade() of G4FTFModel >> 918 // ------------------------------------------------------ >> 919 >> 920 //G4cout<<"Mproj "<<Pprojectile.mag()<<G4endl; >> 921 //G4cout<<"Mtarg "<<Ptarget.mag()<<G4endl; >> 922 projectile->Set4Momentum(Pprojectile); >> 923 target->Set4Momentum(Ptarget); 1051 924 1052 G4bool G4DiffractiveExcitation:: << 925 projectile->IncrementCollisionCount(1); 1053 ExciteParticipants_doNonDiffraction( G4VSplit << 926 target->IncrementCollisionCount(1); 1054 G4VSplit << 1055 G4FTFPar << 1056 G4Diffra << 1057 // Third of the three utility methods used << 1058 // it does the sampling for the non-diffrac << 1059 << 1060 #ifdef debugFTFexcitation << 1061 G4cout << "Non-diffraction process" << G4en << 1062 #endif << 1063 << 1064 // Check whether the interaction is possibl << 1065 common.ProjMassT2 = common.ProjectileNonDif << 1066 common.ProjMassT = common.ProjectileNonDif << 1067 common.TargMassT2 = common.TargetNonDiffSta << 1068 common.TargMassT = common.TargetNonDiffSta << 1069 if ( common.SqrtS < common.ProjMassT + comm << 1070 << 1071 common.PZcms2 = ( sqr( common.S ) + sqr( co << 1072 - 2.0 * ( common.S * ( comm << 1073 + common.ProjMassT2 * << 1074 common.maxPtSquare = common.PZcms2; << 1075 << 1076 G4int whilecount = 0; << 1077 do { // Generate pt and masses << 1078 << 1079 whilecount++; << 1080 if ( whilecount > 1000 ) { << 1081 common.Qmomentum = G4LorentzVector( 0.0 << 1082 return false; // Ignore this interacti << 1083 }; << 1084 << 1085 common.Qmomentum = G4LorentzVector( Gauss << 1086 << 1087 common.Pt2 = G4ThreeVector( common.Qmomen << 1088 common.ProjMassT2 = common.ProjectileNonD << 1089 common.ProjMassT = std::sqrt( common.Pro << 1090 common.TargMassT2 = common.TargetNonDiffS << 1091 common.TargMassT = std::sqrt( common.Tar << 1092 if ( common.SqrtS < common.ProjMassT + co << 1093 << 1094 common.PZcms2 =( sqr( common.S ) + sqr( c << 1095 - 2.0 * ( common.S * ( c << 1096 + common.ProjM << 1097 if ( common.PZcms2 < 0.0 ) continue; << 1098 << 1099 common.PZcms = std::sqrt( common.PZcms2 ) << 1100 common.PMinusMin = std::sqrt( common.Proj << 1101 common.PMinusMax = common.SqrtS - common. << 1102 common.TPlusMin = std::sqrt( common.TargM << 1103 common.TPlusMax = common.SqrtS - common.P << 1104 << 1105 if ( G4UniformRand() <= 0.5 ) { // Rando << 1106 if ( G4UniformRand() < theParameters->G << 1107 common.PMinusNew = ChooseP( common.PM << 1108 } else { << 1109 common.PMinusNew = ( common.PMinusMax << 1110 } << 1111 if ( G4UniformRand() < theParameters->G << 1112 common.TPlusNew = ChooseP( common.TPl << 1113 } else { << 1114 common.TPlusNew = ( common.TPlusMax - << 1115 } << 1116 } else { << 1117 if ( G4UniformRand() < theParameters->G << 1118 common.TPlusNew = ChooseP( common.TPl << 1119 } else { << 1120 common.TPlusNew = ( common.TPlusMax - << 1121 } << 1122 if ( G4UniformRand() < theParameters->G << 1123 common.PMinusNew = ChooseP( common.PM << 1124 } else { << 1125 common.PMinusNew = ( common.PMinusMax << 1126 } << 1127 } << 1128 927 1129 common.Qminus = common.PMinusNew - common << 928 return true; 1130 common.Qplus = -( common.TPlusNew - commo << 1131 common.Qmomentum.setPz( (common.Qplus - c << 1132 common.Qmomentum.setE( (common.Qplus + c << 1133 #ifdef debugFTFexcitation << 1134 G4cout <<"Sampled: Mpr, MdifPr, Mtr, Mdif << 1135 << ( common.Pprojectile + common.Q << 1136 << common.ProjectileNonDiffStateMi << 1137 << ( common.Ptarget - common.Qmome << 1138 << common.TargetNonDiffStateMinMas << 1139 #endif << 1140 << 1141 } while ( ( common.Pprojectile + common.Qmo << 1142 common.ProjectileNonDiffStateMinM << 1143 ( common.Ptarget - common.Qmoment << 1144 common.TargetNonDiffStateMinMass2 << 1145 << 1146 projectile->SetStatus( 0 ); << 1147 target->SetStatus( 0 ); << 1148 return true; << 1149 } 929 } 1150 930 >> 931 // --------------------------------------------------------------------- >> 932 void G4DiffractiveExcitation::CreateStrings(G4VSplitableHadron * hadron, >> 933 G4bool isProjectile, >> 934 G4ExcitedString * &FirstString, >> 935 G4ExcitedString * &SecondString, >> 936 G4FTFParameters *theParameters) const >> 937 { >> 938 hadron->SplitUp(); >> 939 G4Parton *start= hadron->GetNextParton(); >> 940 if ( start==NULL) >> 941 { G4cout << " G4FTFModel::String() Error:No start parton found"<< G4endl; >> 942 FirstString=0; SecondString=0; >> 943 return; >> 944 } >> 945 G4Parton *end = hadron->GetNextParton(); >> 946 if ( end==NULL) >> 947 { G4cout << " G4FTFModel::String() Error:No end parton found"<< G4endl; >> 948 FirstString=0; SecondString=0; >> 949 return; >> 950 } >> 951 >> 952 G4LorentzVector Phadron=hadron->Get4Momentum(); >> 953 >> 954 G4LorentzVector Pstart(0.,0.,0.,0.); >> 955 G4LorentzVector Pend(0.,0.,0.,0.); >> 956 G4LorentzVector Pkink(0.,0.,0.,0.); >> 957 G4LorentzVector PkinkQ1(0.,0.,0.,0.); >> 958 G4LorentzVector PkinkQ2(0.,0.,0.,0.); >> 959 >> 960 G4int PDGcode_startQ = std::abs(start->GetDefinition()->GetPDGEncoding()); >> 961 G4int PDGcode_endQ = std::abs( end->GetDefinition()->GetPDGEncoding()); >> 962 >> 963 //-------------------------------------------------------------------------------- >> 964 G4double Wmin(0.); >> 965 if(isProjectile) >> 966 { >> 967 Wmin=theParameters->GetProjMinDiffMass(); >> 968 } else >> 969 { >> 970 Wmin=theParameters->GetTarMinDiffMass(); >> 971 } // end of if(isProjectile) >> 972 >> 973 G4double W = hadron->Get4Momentum().mag(); >> 974 G4double W2=W*W; >> 975 >> 976 G4double Pt(0.), x1(0.), x2(0.), x3(0.); >> 977 G4bool Kink=false; >> 978 >> 979 if(W > Wmin) >> 980 { // Kink is possible >> 981 G4double Pt2kink=theParameters->GetPt2Kink(); >> 982 Pt = std::sqrt(Pt2kink*(std::pow(W2/16./Pt2kink+1.,G4UniformRand()) - 1.)); >> 983 >> 984 if(Pt > 500.*MeV) >> 985 { >> 986 G4double Ymax = std::log(W/2./Pt + std::sqrt(W2/4./Pt/Pt - 1.)); >> 987 G4double Y=Ymax*(1.- 2.*G4UniformRand()); >> 988 >> 989 x1=1.-Pt/W*std::exp( Y); >> 990 x3=1.-Pt/W*std::exp(-Y); >> 991 x2=2.-x1-x3; >> 992 >> 993 G4double Mass_startQ = 650.*MeV; >> 994 if(PDGcode_startQ < 3) Mass_startQ = 325.*MeV; >> 995 if(PDGcode_startQ == 3) Mass_startQ = 500.*MeV; >> 996 if(PDGcode_startQ == 4) Mass_startQ = 1600.*MeV; >> 997 >> 998 G4double Mass_endQ = 650.*MeV; >> 999 if(PDGcode_endQ < 3) Mass_endQ = 325.*MeV; >> 1000 if(PDGcode_endQ == 3) Mass_endQ = 500.*MeV; >> 1001 if(PDGcode_endQ == 4) Mass_endQ = 1600.*MeV; >> 1002 >> 1003 G4double P2_1=W2*x1*x1/4.-Mass_endQ *Mass_endQ; >> 1004 G4double P2_3=W2*x3*x3/4.-Mass_startQ*Mass_startQ; >> 1005 >> 1006 G4double P2_2=sqr((2.-x1-x3)*W/2.); >> 1007 >> 1008 if((P2_1 <= 0.) || (P2_3 <= 0.)) >> 1009 { Kink=false;} >> 1010 else >> 1011 { >> 1012 G4double P_1=std::sqrt(P2_1); >> 1013 G4double P_2=std::sqrt(P2_2); >> 1014 G4double P_3=std::sqrt(P2_3); >> 1015 >> 1016 G4double CosT12=(P2_3-P2_1-P2_2)/(2.*P_1*P_2); >> 1017 G4double CosT13=(P2_2-P2_1-P2_3)/(2.*P_1*P_3); >> 1018 // Pt=P_2*std::sqrt(1.-CosT12*CosT12); // because system was rotated 11.12.09 >> 1019 >> 1020 if((std::abs(CosT12) >1.) || (std::abs(CosT13) > 1.)) >> 1021 { Kink=false;} >> 1022 else >> 1023 { >> 1024 Kink=true; >> 1025 Pt=P_2*std::sqrt(1.-CosT12*CosT12); // because system was rotated 11.12.09 >> 1026 Pstart.setPx(-Pt); Pstart.setPy(0.); Pstart.setPz(P_3*CosT13); >> 1027 Pend.setPx( 0.); Pend.setPy( 0.); Pend.setPz( P_1); >> 1028 Pkink.setPx( Pt); Pkink.setPy( 0.); Pkink.setPz( P_2*CosT12); >> 1029 Pstart.setE(x3*W/2.); >> 1030 Pkink.setE(Pkink.vect().mag()); >> 1031 Pend.setE(x1*W/2.); >> 1032 >> 1033 G4double XkQ=GetQuarkFractionOfKink(0.,1.); >> 1034 if(Pkink.getZ() > 0.) >> 1035 { >> 1036 if(XkQ > 0.5) {PkinkQ1=XkQ*Pkink;} else {PkinkQ1=(1.-XkQ)*Pkink;} >> 1037 } else { >> 1038 if(XkQ > 0.5) {PkinkQ1=(1.-XkQ)*Pkink;} else {PkinkQ1=XkQ*Pkink;} >> 1039 } >> 1040 >> 1041 PkinkQ2=Pkink - PkinkQ1; >> 1042 //------------------------- Minimizing Pt1^2+Pt3^2 ------------------------------ >> 1043 >> 1044 G4double Cos2Psi=(sqr(x1) -sqr(x3)+2.*sqr(x3*CosT13))/ >> 1045 std::sqrt(sqr(sqr(x1)-sqr(x3)) + sqr(2.*x1*x3*CosT13)); >> 1046 G4double Psi=std::acos(Cos2Psi); >> 1047 >> 1048 G4LorentzRotation Rotate; >> 1049 if(isProjectile) {Rotate.rotateY(Psi);} >> 1050 else {Rotate.rotateY(pi-Psi);} >> 1051 Rotate.rotateZ(twopi*G4UniformRand()); >> 1052 >> 1053 Pstart*=Rotate; >> 1054 Pkink*=Rotate; >> 1055 PkinkQ1*=Rotate; >> 1056 PkinkQ2*=Rotate; >> 1057 Pend*=Rotate; >> 1058 >> 1059 } >> 1060 } // end of if((P2_1 < 0.) || (P2_3 <0.)) >> 1061 } // end of if(Pt > 500.*MeV) >> 1062 } // end of if(W > Wmin) Check for a kink >> 1063 >> 1064 //-------------------------------------------------------------------------------- >> 1065 >> 1066 if(Kink) >> 1067 { // Kink is possible >> 1068 std::vector<G4double> QuarkProbabilitiesAtGluonSplitUp = >> 1069 theParameters->GetQuarkProbabilitiesAtGluonSplitUp(); >> 1070 >> 1071 G4int QuarkInGluon(1); G4double Ksi=G4UniformRand(); >> 1072 for(unsigned int Iq=0; Iq <3; Iq++) >> 1073 { >> 1074 >> 1075 if(Ksi > QuarkProbabilitiesAtGluonSplitUp[Iq]) QuarkInGluon++;} >> 1076 >> 1077 G4Parton * Gquark = new G4Parton(QuarkInGluon); >> 1078 G4Parton * Ganti_quark = new G4Parton(-QuarkInGluon); >> 1079 >> 1080 //------------------------------------------------------------------------------- >> 1081 G4LorentzRotation toCMS(-1*Phadron.boostVector()); >> 1082 >> 1083 G4LorentzRotation toLab(toCMS.inverse()); >> 1084 >> 1085 Pstart.transform(toLab); start->Set4Momentum(Pstart); >> 1086 PkinkQ1.transform(toLab); >> 1087 PkinkQ2.transform(toLab); >> 1088 Pend.transform(toLab); end->Set4Momentum(Pend); >> 1089 >> 1090 G4int absPDGcode=std::abs(hadron->GetDefinition()->GetPDGEncoding()); >> 1091 >> 1092 if(absPDGcode < 1000) >> 1093 { // meson >> 1094 if ( isProjectile ) >> 1095 { // Projectile >> 1096 if(end->GetDefinition()->GetPDGEncoding() > 0 ) // A quark on the end >> 1097 { // Quark on the end >> 1098 FirstString = new G4ExcitedString(end ,Ganti_quark, +1); >> 1099 SecondString= new G4ExcitedString(Gquark,start ,+1); >> 1100 Ganti_quark->Set4Momentum(PkinkQ1); >> 1101 Gquark->Set4Momentum(PkinkQ2); >> 1102 >> 1103 } else >> 1104 { // Anti_Quark on the end >> 1105 FirstString = new G4ExcitedString(end ,Gquark, +1); >> 1106 SecondString= new G4ExcitedString(Ganti_quark,start ,+1); >> 1107 Gquark->Set4Momentum(PkinkQ1); >> 1108 Ganti_quark->Set4Momentum(PkinkQ2); >> 1109 >> 1110 } // end of if(end->GetPDGcode() > 0) >> 1111 } else { // Target >> 1112 if(end->GetDefinition()->GetPDGEncoding() > 0 ) // A quark on the end >> 1113 { // Quark on the end >> 1114 FirstString = new G4ExcitedString(Ganti_quark,end ,-1); >> 1115 SecondString= new G4ExcitedString(start ,Gquark,-1); >> 1116 Ganti_quark->Set4Momentum(PkinkQ2); >> 1117 Gquark->Set4Momentum(PkinkQ1); >> 1118 >> 1119 } else >> 1120 { // Anti_Quark on the end >> 1121 FirstString = new G4ExcitedString(Gquark,end ,-1); >> 1122 SecondString= new G4ExcitedString(start ,Ganti_quark,-1); >> 1123 Gquark->Set4Momentum(PkinkQ2); >> 1124 Ganti_quark->Set4Momentum(PkinkQ1); >> 1125 >> 1126 } // end of if(end->GetPDGcode() > 0) >> 1127 } // end of if ( isProjectile ) >> 1128 } else // if(absPDGCode < 1000) >> 1129 { // Baryon/AntiBaryon >> 1130 if ( isProjectile ) >> 1131 { // Projectile >> 1132 if((end->GetDefinition()->GetParticleType() == "diquarks") && >> 1133 (end->GetDefinition()->GetPDGEncoding() > 0 ) ) >> 1134 { // DiQuark on the end >> 1135 FirstString = new G4ExcitedString(end ,Gquark, +1); >> 1136 SecondString= new G4ExcitedString(Ganti_quark,start ,+1); >> 1137 Gquark->Set4Momentum(PkinkQ1); >> 1138 Ganti_quark->Set4Momentum(PkinkQ2); >> 1139 >> 1140 } else >> 1141 { // Anti_DiQuark on the end or quark >> 1142 FirstString = new G4ExcitedString(end ,Ganti_quark, +1); >> 1143 SecondString= new G4ExcitedString(Gquark,start ,+1); >> 1144 Ganti_quark->Set4Momentum(PkinkQ1); >> 1145 Gquark->Set4Momentum(PkinkQ2); >> 1146 >> 1147 } // end of if(end->GetPDGcode() > 0) >> 1148 } else { // Target >> 1149 >> 1150 if((end->GetDefinition()->GetParticleType() == "diquarks") && >> 1151 (end->GetDefinition()->GetPDGEncoding() > 0 ) ) >> 1152 { // DiQuark on the end >> 1153 FirstString = new G4ExcitedString(Gquark,end ,-1); >> 1154 >> 1155 SecondString= new G4ExcitedString(start ,Ganti_quark,-1); >> 1156 Gquark->Set4Momentum(PkinkQ1); >> 1157 Ganti_quark->Set4Momentum(PkinkQ2); >> 1158 >> 1159 } else >> 1160 { // Anti_DiQuark on the end or Q >> 1161 FirstString = new G4ExcitedString(Ganti_quark,end ,-1); >> 1162 SecondString= new G4ExcitedString(start ,Gquark,-1); >> 1163 Gquark->Set4Momentum(PkinkQ2); >> 1164 Ganti_quark->Set4Momentum(PkinkQ1); >> 1165 >> 1166 } // end of if(end->GetPDGcode() > 0) >> 1167 } // end of if ( isProjectile ) >> 1168 } // end of if(absPDGcode < 1000) >> 1169 >> 1170 FirstString->SetTimeOfCreation(hadron->GetTimeOfCreation()); >> 1171 FirstString->SetPosition(hadron->GetPosition()); >> 1172 >> 1173 SecondString->SetTimeOfCreation(hadron->GetTimeOfCreation()); >> 1174 SecondString->SetPosition(hadron->GetPosition()); >> 1175 >> 1176 // ------------------------------------------------------------------------- >> 1177 } else // End of kink is possible >> 1178 { // Kink is impossible >> 1179 if ( isProjectile ) >> 1180 { >> 1181 FirstString= new G4ExcitedString(end,start, +1); >> 1182 } else { >> 1183 FirstString= new G4ExcitedString(start,end, -1); >> 1184 } 1151 1185 1152 //=========================================== << 1186 SecondString=0; 1153 1187 1154 void G4DiffractiveExcitation::CreateStrings( << 1188 FirstString->SetTimeOfCreation(hadron->GetTimeOfCreation()); 1155 << 1189 FirstString->SetPosition(hadron->GetPosition()); 1156 << 1157 << 1158 << 1159 << 1160 //G4cout << "Create Strings SplitUp " << ha << 1161 // << "Defin " << hadron->GetDefiniti << 1162 // << "Defin " << hadron->GetDefiniti << 1163 << 1164 G4bool HadronIsString = hadron->IsSplit(); << 1165 if( ! HadronIsString ) hadron->SplitUp(); << 1166 << 1167 G4Parton* start = hadron->GetNextParton(); << 1168 if ( start == nullptr ) { << 1169 G4cout << " G4FTFModel::String() Error: N << 1170 FirstString = 0; SecondString = 0; << 1171 return; << 1172 } << 1173 << 1174 G4Parton* end = hadron->GetNextParton(); << 1175 if ( end == nullptr ) { << 1176 G4cout << " G4FTFModel::String() Error: N << 1177 FirstString = 0; SecondString = 0; << 1178 return; << 1179 } << 1180 << 1181 //G4cout << start << " " << start->GetPDGco << 1182 // << G4endl << 1183 // << "Create string " << start->GetP << 1184 << 1185 if ( HadronIsString ) { << 1186 if ( isProjectile ) { << 1187 FirstString = new G4ExcitedString( end, << 1188 } else { << 1189 FirstString = new G4ExcitedString( end, << 1190 } << 1191 FirstString->SetTimeOfCreation( hadron->G << 1192 FirstString->SetPosition( hadron->GetPosi << 1193 SecondString = 0; << 1194 return; << 1195 } << 1196 << 1197 G4LorentzVector Phadron = hadron->Get4Momen << 1198 //G4cout << "String mom " << Phadron << G4e << 1199 G4LorentzVector Pstart( 0.0, 0.0, 0.0, 0.0 << 1200 G4LorentzVector Pend( 0.0, 0.0, 0.0, 0.0 ); << 1201 G4LorentzVector Pkink( 0.0, 0.0, 0.0, 0.0 ) << 1202 G4LorentzVector PkinkQ1( 0.0, 0.0, 0.0, 0.0 << 1203 G4LorentzVector PkinkQ2( 0.0, 0.0, 0.0, 0.0 << 1204 << 1205 G4int PDGcode_startQ = std::abs( start->Get << 1206 G4int PDGcode_endQ = std::abs( end->Get << 1207 //G4cout << "PDGcode_startQ " << PDGcode_st << 1208 << 1209 G4double Wmin( 0.0 ); << 1210 if ( isProjectile ) { << 1211 Wmin = theParameters->GetProjMinDiffMass( << 1212 } else { << 1213 Wmin = theParameters->GetTarMinDiffMass() << 1214 } << 1215 << 1216 G4double W = hadron->Get4Momentum().mag(); << 1217 G4double W2 = W*W; << 1218 G4double Pt( 0.0 ), x1( 0.0 ), x3( 0.0 ); << 1219 G4bool Kink = false; << 1220 << 1221 if ( ! ( ( start->GetDefinition()->GetParti << 1222 end->GetDefinition()->GetParti << 1223 ( start->GetDefinition()->GetParti << 1224 end->GetDefinition()->GetParti << 1225 // Kinky strings are allowed only for qq- << 1226 // Kinky strings are impossible for other << 1227 // according to the analysis of Pbar P in << 1228 << 1229 if ( W > Wmin ) { // Kink is possible << 1230 if ( hadron->GetStatus() == 0 ) { << 1231 G4double Pt2kink = theParameters->Get << 1232 if ( Pt2kink ) { << 1233 Pt = std::sqrt( Pt2kink * ( G4Pow:: << 1234 } else { << 1235 Pt = 0.0; << 1236 } << 1237 } else { << 1238 Pt = 0.0; << 1239 } << 1240 1190 1241 if ( Pt > 500.0*MeV ) { << 1191 // momenta of string ends 1242 G4double Ymax = G4Log( W/2.0/Pt + std << 1192 // 1243 G4double Y = Ymax*( 1.0 - 2.0*G4Unifo << 1193 G4double Momentum=hadron->Get4Momentum().vect().mag(); 1244 x1 = 1.0 - Pt/W * G4Exp( Y ); << 1194 G4double Plus=hadron->Get4Momentum().e() + Momentum; 1245 x3 = 1.0 - Pt/W * G4Exp(-Y ); << 1195 G4double Minus=hadron->Get4Momentum().e() - Momentum; 1246 //x2 = 2.0 - x1 - x3; << 1196 1247 << 1197 G4ThreeVector tmp; 1248 G4double Mass_startQ = 650.0*MeV; << 1198 if(Momentum > 0.) 1249 if ( PDGcode_startQ < 3 ) Mass_start << 1199 { 1250 if ( PDGcode_startQ == 3 ) Mass_start << 1200 tmp.set(hadron->Get4Momentum().px(), 1251 if ( PDGcode_startQ == 4 ) Mass_start << 1201 hadron->Get4Momentum().py(), 1252 if ( PDGcode_startQ == 5 ) Mass_start << 1202 hadron->Get4Momentum().pz()); 1253 G4double Mass_endQ = 650.0*MeV; << 1203 tmp/=Momentum; 1254 if ( PDGcode_endQ < 3 ) Mass_endQ = << 1204 } 1255 if ( PDGcode_endQ == 3 ) Mass_endQ = << 1205 else 1256 if ( PDGcode_endQ == 4 ) Mass_endQ = << 1206 { 1257 if ( PDGcode_endQ == 5 ) Mass_endQ = << 1207 tmp.set(0.,0.,1.); 1258 << 1259 G4double P2_1 = W2*x1*x1/4.0 - Mass_e << 1260 G4double P2_3 = W2*x3*x3/4.0 - Mass_s << 1261 G4double P2_2 = sqr( (2.0 - x1 - x3)* << 1262 if ( P2_1 <= 0.0 || P2_3 <= 0.0 ) { << 1263 Kink = false; << 1264 } else { << 1265 G4double P_1 = std::sqrt( P2_1 ); << 1266 G4double P_2 = std::sqrt( P2_2 ); << 1267 G4double P_3 = std::sqrt( P2_3 ); << 1268 G4double CosT12 = ( P2_3 - P2_1 - P << 1269 G4double CosT13 = ( P2_2 - P2_1 - P << 1270 << 1271 if ( std::abs( CosT12 ) > 1.0 || << 1272 Kink = false; << 1273 } else { << 1274 Kink = true; << 1275 Pt = P_2 * std::sqrt( 1.0 - CosT1 << 1276 Pstart.setPx( -Pt ); Pstart.setPy << 1277 Pend.setPx( 0.0 ); Pend.setPy( << 1278 Pkink.setPx( Pt ); Pkink.setPy( << 1279 Pstart.setE( x3*W/2.0 ); << 1280 Pkink.setE( Pkink.vect().mag() ); << 1281 Pend.setE( x1*W/2.0 ); << 1282 << 1283 G4double XkQ = GetQuarkFractionOf << 1284 if ( Pkink.getZ() > 0.0 ) { << 1285 if ( XkQ > 0.5 ) { << 1286 PkinkQ1 = XkQ*Pkink; << 1287 } else { << 1288 PkinkQ1 = (1.0 - XkQ)*Pkink; << 1289 } << 1290 } else { << 1291 if ( XkQ > 0.5 ) { << 1292 PkinkQ1 = (1.0 - XkQ)*Pkink; << 1293 } else { << 1294 PkinkQ1 = XkQ*Pkink; << 1295 } << 1296 } << 1297 << 1298 PkinkQ2 = Pkink - PkinkQ1; << 1299 // Minimizing Pt1^2+Pt3^2 << 1300 G4double Cos2Psi = ( sqr(x1) - sq << 1301 std::sqrt( sqr << 1302 G4double Psi = std::acos( Cos2Psi << 1303 << 1304 G4LorentzRotation Rotate; << 1305 if ( isProjectile ) { << 1306 Rotate.rotateY( Psi ); << 1307 } else { << 1308 Rotate.rotateY( pi + Psi ); << 1309 } << 1310 Rotate.rotateZ( twopi * G4Uniform << 1311 Pstart *= Rotate; << 1312 Pkink *= Rotate; << 1313 PkinkQ1 *= Rotate; << 1314 PkinkQ2 *= Rotate; << 1315 Pend *= Rotate; << 1316 } 1208 } 1317 } // End of if ( P2_1 <= 0.0 || P2 << 1318 } // End of if ( Pt > 500.0*MeV ) << 1319 } // End of if ( W > Wmin ) : check for a << 1320 } // end of qq-q string selection << 1321 << 1322 if ( Kink ) { // Kink is possible << 1323 << 1324 //G4cout << "Kink is sampled!" << G4endl; << 1325 std::vector< G4double > QuarkProbabilitie << 1326 theParameters->GetQuarkProbabilitiesA << 1327 << 1328 G4int QuarkInGluon( 1 ); G4double Ksi = G << 1329 for ( unsigned int Iq = 0; Iq < 3; Iq++ ) << 1330 //G4cout << "Iq " << Iq << G4endl; << 1331 if ( Ksi > QuarkProbabilitiesAtGluonSpl << 1332 } << 1333 //G4cout << "Last Iq " << QuarkInGluon << << 1334 G4Parton* Gquark = new G4Parton( QuarkInG << 1335 G4Parton* Ganti_quark = new G4Parton( -Qu << 1336 //G4cout << "Lorentz " << G4endl; << 1337 << 1338 G4LorentzRotation toCMS( -1 * Phadron.boo << 1339 G4LorentzRotation toLab( toCMS.inverse() << 1340 //G4cout << "Pstart " << Pstart << G4endl << 1341 //G4cout << "Pend " << Pend << G4endl; << 1342 //G4cout << "Kink1 " <<PkinkQ1<<G4endl; << 1343 //G4cout << "Kink2 " <<PkinkQ2<<G4endl; << 1344 //G4cout << "Pstart " << Pstart << G4endl << 1345 << 1346 Pstart.transform( toLab ); start->Set4Mo << 1347 PkinkQ1.transform( toLab ); << 1348 PkinkQ2.transform( toLab ); << 1349 Pend.transform( toLab ); end->Set4Mome << 1350 //G4cout << "Pstart " << Pstart << G4endl << 1351 //G4cout << "Pend " << Pend << G4endl; << 1352 //G4cout << "Defin " << hadron->GetDefini << 1353 //G4cout << "Defin " << hadron->GetDefini << 1354 << 1355 //G4int absPDGcode = std::abs( hadron->Ge << 1356 G4int absPDGcode = 1500; << 1357 if ( start->GetDefinition()->GetParticleS << 1358 end->GetDefinition()->GetParticleSub << 1359 absPDGcode = 110; << 1360 } << 1361 //G4cout << "absPDGcode " << absPDGcode < << 1362 << 1363 if ( absPDGcode < 1000 ) { // meson << 1364 if ( isProjectile ) { // Projectile << 1365 if ( end->GetDefinition()->GetPDGEnco << 1366 FirstString = new G4ExcitedString( << 1367 SecondString = new G4ExcitedString( << 1368 Ganti_quark->Set4Momentum( PkinkQ1 << 1369 Gquark->Set4Momentum( PkinkQ2 ); << 1370 } else { // Anti_Quark on the end << 1371 FirstString = new G4ExcitedString( << 1372 SecondString = new G4ExcitedString( << 1373 Gquark->Set4Momentum( PkinkQ1 ); << 1374 Ganti_quark->Set4Momentum( PkinkQ2 << 1375 } << 1376 } else { // Target << 1377 if ( end->GetDefinition()->GetPDGEnco << 1378 FirstString = new G4ExcitedString( << 1379 SecondString = new G4ExcitedString( << 1380 Ganti_quark->Set4Momentum( PkinkQ2 << 1381 Gquark->Set4Momentum( PkinkQ1 ); << 1382 } else { // Anti_Quark on the end << 1383 FirstString = new G4ExcitedString( << 1384 SecondString = new G4ExcitedString( << 1385 Gquark->Set4Momentum( PkinkQ2 ); << 1386 Ganti_quark->Set4Momentum( PkinkQ1 << 1387 } << 1388 } << 1389 } else { // Baryon/AntiBaryon << 1390 if ( isProjectile ) { // Projectile << 1391 if ( end->GetDefinition()->GetParticl << 1392 end->GetDefinition()->GetPDGEnco << 1393 FirstString = new G4ExcitedString( << 1394 SecondString = new G4ExcitedString( << 1395 Gquark->Set4Momentum( PkinkQ1 ); << 1396 Ganti_quark->Set4Momentum( PkinkQ2 << 1397 } else { / << 1398 FirstString = new G4ExcitedString( << 1399 SecondString = new G4ExcitedString( << 1400 Ganti_quark->Set4Momentum( PkinkQ1 << 1401 Gquark->Set4Momentum( PkinkQ2 ); << 1402 } << 1403 } else { // Target << 1404 if ( end->GetDefinition()->GetParticl << 1405 end->GetDefinition()->GetPDGEnco << 1406 Gquark->Set4Momentum( PkinkQ1 ); << 1407 Ganti_quark->Set4Momentum( PkinkQ2 << 1408 FirstString = new G4ExcitedString( << 1409 SecondString = new G4ExcitedString( << 1410 } else { // Anti_DiQuark on the end << 1411 FirstString = new G4ExcitedString( << 1412 SecondString = new G4ExcitedString( << 1413 Gquark->Set4Momentum( PkinkQ2 ); << 1414 Ganti_quark->Set4Momentum( PkinkQ1 << 1415 } << 1416 } << 1417 } << 1418 << 1419 FirstString->SetTimeOfCreation( hadron->G << 1420 FirstString->SetPosition( hadron->GetPosi << 1421 SecondString->SetTimeOfCreation( hadron-> << 1422 SecondString->SetPosition( hadron->GetPos << 1423 << 1424 } else { // Kink is impossible << 1425 << 1426 if ( isProjectile ) { << 1427 FirstString = new G4ExcitedString( end, << 1428 } else { << 1429 FirstString = new G4ExcitedString( end, << 1430 } << 1431 FirstString->SetTimeOfCreation( hadron->G << 1432 FirstString->SetPosition( hadron->GetPosi << 1433 SecondString = 0; << 1434 // momenta of string ends << 1435 G4LorentzVector HadronMom = hadron->Get4M << 1436 G4LorentzVector Pstart1 = G4LorentzVector << 1437 G4LorentzVector Pend1 = G4LorentzVector << 1438 G4double Pz = HadronMom.pz(); << 1439 G4double Eh = HadronMom.e(); << 1440 G4double Pt2 = sqr( HadronMom.px() ) + sq << 1441 G4double Mt2 = HadronMom.mt2(); << 1442 G4double Exp = std::sqrt( sqr(Pz) + ( sqr << 1443 G4double Pzq = Pz/2.0 - Exp; << 1444 G4double Eq = std::sqrt( sqr(Pzq) + Pt2 << 1445 G4double Pzqq = Pz/2.0 + Exp; << 1446 G4double Eqq = std::sqrt( sqr(Pzqq) + Pt << 1447 start->Set4Momentum( Pstart1 ); << 1448 end->Set4Momentum( Pend1 ); << 1449 Pstart = Pstart1; Pend = Pend1; << 1450 << 1451 } // End of "if (Kink)" << 1452 << 1453 //G4cout << "Quarks in the string at creati << 1454 // << " " << FirstString->GetLeftPart << 1455 // << FirstString << " " << SecondStr << 1456 << 1457 #ifdef G4_FTFDEBUG << 1458 G4cout << " generated string flavors << 1459 << end->GetPDGcode() << G4endl << " << 1460 << start->Get4Momentum() << "mass : << 1461 << " generated string momenta: Diqua << 1462 << end->Get4Momentum().mag() << G4en << 1463 << Pstart + Pend << G4endl << " Orig << 1464 << hadron->Get4Momentum() << " "<<ha << 1465 #endif << 1466 << 1467 return; << 1468 } << 1469 1209 >> 1210 G4LorentzVector Pstart(tmp,0.); >> 1211 G4LorentzVector Pend(tmp,0.); 1470 1212 1471 //=========================================== << 1213 if(isProjectile) >> 1214 { >> 1215 Pstart*=(-1.)*Minus/2.; >> 1216 Pend *=(+1.)*Plus /2.; >> 1217 } >> 1218 else >> 1219 { >> 1220 Pstart*=(+1.)*Plus/2.; >> 1221 Pend *=(-1.)*Minus/2.; >> 1222 } 1472 1223 1473 G4double G4DiffractiveExcitation::ChooseP( G4 << 1224 Momentum=-Pstart.mag(); 1474 // Choose an x between Xmin and Xmax with P << 1225 Pstart.setT(Momentum); // It is assumed that quark has m=0. 1475 // To be improved... << 1476 G4double range = Pmax - Pmin; << 1477 if ( Pmin <= 0.0 || range <= 0.0 ) { << 1478 G4cout << " Pmin, range : " << Pmin << " << 1479 throw G4HadronicException( __FILE__, __LI << 1480 "G4Diffractive << 1481 } << 1482 G4double P = Pmin * G4Pow::GetInstance()->p << 1483 //G4double P = (Pmax - Pmin) * G4UniformRan << 1484 return P; << 1485 } << 1486 1226 >> 1227 Momentum=-Pend.mag(); >> 1228 Pend.setT(Momentum); // It is assumed that di-quark has m=0. 1487 1229 1488 //=========================================== << 1230 start->Set4Momentum(Pstart); >> 1231 end->Set4Momentum(Pend); >> 1232 SecondString=0; >> 1233 } // End of kink is impossible >> 1234 >> 1235 #ifdef G4_FTFDEBUG >> 1236 G4cout << " generated string flavors " >> 1237 << start->GetPDGcode() << " / " >> 1238 << end->GetPDGcode() << G4endl; >> 1239 G4cout << " generated string momenta: quark " >> 1240 << start->Get4Momentum() << "mass : " >> 1241 <<start->Get4Momentum().mag() << G4endl; >> 1242 G4cout << " generated string momenta: Diquark " >> 1243 << end ->Get4Momentum() >> 1244 << "mass : " <<end->Get4Momentum().mag()<< G4endl; >> 1245 G4cout << " sum of ends " << Pstart+Pend << G4endl; >> 1246 G4cout << " Original " << hadron->Get4Momentum() << G4endl; >> 1247 #endif 1489 1248 1490 G4ThreeVector G4DiffractiveExcitation::Gaussi << 1249 return; 1491 // @@ this method is used in FTFModel as w << 1250 1492 G4double Pt2( 0.0 ); << 1493 if ( AveragePt2 <= 0.0 ) { << 1494 Pt2 = 0.0; << 1495 } else { << 1496 Pt2 = -AveragePt2 * G4Log( 1.0 + G4Unifor << 1497 ( G4Ex << 1498 } << 1499 G4double Pt = std::sqrt( Pt2 ); << 1500 G4double phi = G4UniformRand() * twopi; << 1501 return G4ThreeVector( Pt * std::cos( phi ), << 1502 } 1251 } 1503 1252 1504 1253 1505 //=========================================== << 1254 // --------- private methods ---------------------- 1506 1255 1507 G4double G4DiffractiveExcitation::GetQuarkFra << 1256 // --------------------------------------------------------------------- 1508 G4double z, yf; << 1257 G4double G4DiffractiveExcitation::ChooseP(G4double Pmin, G4double Pmax) const 1509 const G4int maxNumberOfLoops = 10000; << 1258 { 1510 G4int loopCounter = 0; << 1259 // choose an x between Xmin and Xmax with P(x) ~ 1/x 1511 do { << 1260 // to be improved... 1512 z = zmin + G4UniformRand() * (zmax - zmin << 1261 1513 yf = z*z + sqr(1.0 - z); << 1262 G4double range=Pmax-Pmin; 1514 } while ( ( G4UniformRand() > yf ) && << 1263 1515 ++loopCounter < maxNumberOfLoops << 1264 if ( Pmin <= 0. || range <=0. ) 1516 if ( loopCounter >= maxNumberOfLoops ) { << 1265 { 1517 z = 0.5*(zmin + zmax); // Just something << 1266 G4cout << " Pmin, range : " << Pmin << " , " << range << G4endl; 1518 } << 1267 throw G4HadronicException(__FILE__, __LINE__, "G4DiffractiveExcitation::ChooseP : Invalid arguments "); 1519 return z; << 1268 } 1520 } << 1521 << 1522 << 1523 //=========================================== << 1524 1269 1525 void G4DiffractiveExcitation::UnpackMeson( co << 1270 G4double P=Pmin * std::pow(Pmax/Pmin,G4UniformRand()); 1526 G4int absIdPDG = std::abs( IdPDG ); << 1271 return P; 1527 if ( ! ( absIdPDG == 111 || absIdPDG == 221 << 1528 absIdPDG == 441 || absIdPDG == 443 << 1529 // All other projectile mesons (including << 1530 Q1 = absIdPDG / 100; << 1531 Q2 = (absIdPDG % 100) / 10; << 1532 G4int anti = 1 - 2 * ( std::max( Q1, Q2 ) << 1533 if ( IdPDG < 0 ) anti *= -1; << 1534 Q1 *= anti; << 1535 Q2 *= -1 * anti; << 1536 } else { << 1537 if ( absIdPDG == 441 || absIdPDG == 443 ) << 1538 Q1 = 4; Q2 = -4; << 1539 } else if ( absIdPDG == 553 ) { << 1540 Q1 = 5; Q2 = -5; << 1541 } else { << 1542 if ( G4UniformRand() < 0.5 ) { << 1543 Q1 = 1; Q2 = -1; << 1544 } else { << 1545 Q1 = 2; Q2 = -2; << 1546 } << 1547 } << 1548 } << 1549 return; << 1550 } 1272 } 1551 1273 1552 << 1274 // --------------------------------------------------------------------- 1553 //=========================================== << 1275 G4ThreeVector G4DiffractiveExcitation::GaussianPt(G4double AveragePt2, 1554 << 1276 G4double maxPtSquare) const 1555 void G4DiffractiveExcitation::UnpackBaryon( G << 1277 { // @@ this method is used in FTFModel as well. Should go somewhere common! 1556 G << 1278 1557 Q1 = IdPDG / 1000; << 1279 G4double Pt2(0.); 1558 Q2 = (IdPDG % 1000) / 100; << 1280 if(AveragePt2 <= 0.) {Pt2=0.;} 1559 Q3 = (IdPDG % 100) / 10; << 1281 else 1560 return; << 1282 { >> 1283 Pt2 = -AveragePt2 * std::log(1. + G4UniformRand() * >> 1284 (std::exp(-maxPtSquare/AveragePt2)-1.)); >> 1285 } >> 1286 G4double Pt=std::sqrt(Pt2); >> 1287 G4double phi=G4UniformRand() * twopi; >> 1288 return G4ThreeVector (Pt*std::cos(phi), Pt*std::sin(phi), 0.); 1561 } 1289 } 1562 1290 >> 1291 // --------------------------------------------------------------------- >> 1292 G4double G4DiffractiveExcitation::GetQuarkFractionOfKink(G4double zmin, G4double zmax) const >> 1293 { >> 1294 G4double z, yf; >> 1295 do { >> 1296 z = zmin + G4UniformRand()*(zmax-zmin); >> 1297 yf = z*z +sqr(1 - z); >> 1298 } >> 1299 while (G4UniformRand() > yf); >> 1300 return z; >> 1301 } >> 1302 // --------------------------------------------------------------------- >> 1303 void G4DiffractiveExcitation::UnpackMeson(const G4int IdPDG, G4int &Q1, G4int &Q2) const // Uzhi 7.09.09 >> 1304 { >> 1305 G4int absIdPDG = std::abs(IdPDG); >> 1306 Q1= absIdPDG/ 100; >> 1307 Q2= (absIdPDG %100)/10; >> 1308 >> 1309 G4int anti= 1 -2 * ( std::max( Q1, Q2 ) % 2 ); >> 1310 >> 1311 if (IdPDG < 0 ) anti *=-1; >> 1312 Q1 *= anti; >> 1313 Q2 *= -1 * anti; >> 1314 return; >> 1315 } >> 1316 // --------------------------------------------------------------------- >> 1317 void G4DiffractiveExcitation::UnpackBaryon(G4int IdPDG, >> 1318 G4int &Q1, G4int &Q2, G4int &Q3) const // Uzhi 7.09.09 >> 1319 { >> 1320 Q1 = IdPDG / 1000; >> 1321 Q2 = (IdPDG % 1000) / 100; >> 1322 Q3 = (IdPDG % 100) / 10; >> 1323 return; >> 1324 } >> 1325 // --------------------------------------------------------------------- >> 1326 G4int G4DiffractiveExcitation::NewNucleonId(G4int Q1, G4int Q2, G4int Q3) const // Uzhi 7.09.09 >> 1327 { >> 1328 G4int TmpQ(0); >> 1329 if( Q3 > Q2 ) >> 1330 { >> 1331 TmpQ = Q2; >> 1332 Q2 = Q3; >> 1333 Q3 = TmpQ; >> 1334 } else if( Q3 > Q1 ) >> 1335 { >> 1336 TmpQ = Q1; >> 1337 Q1 = Q3; >> 1338 Q3 = TmpQ; >> 1339 } >> 1340 >> 1341 if( Q2 > Q1 ) >> 1342 { >> 1343 TmpQ = Q1; >> 1344 Q1 = Q2; >> 1345 Q2 = TmpQ; >> 1346 } 1563 1347 1564 //=========================================== << 1348 G4int NewCode = Q1*1000 + Q2* 100 + Q3* 10 + 2; 1565 << 1349 return NewCode; 1566 G4int G4DiffractiveExcitation::NewNucleonId( << 1350 } 1567 // Order the three integers in such a way t << 1351 // --------------------------------------------------------------------- 1568 G4int TmpQ( 0 ); << 1352 G4DiffractiveExcitation::G4DiffractiveExcitation(const G4DiffractiveExcitation &) 1569 if ( Q3 > Q2 ) { << 1353 { 1570 TmpQ = Q2; << 1354 throw G4HadronicException(__FILE__, __LINE__, "G4DiffractiveExcitation copy contructor not meant to be called"); 1571 Q2 = Q3; << 1572 Q3 = TmpQ; << 1573 } else if ( Q3 > Q1 ) { << 1574 TmpQ = Q1; << 1575 Q1 = Q3; << 1576 Q3 = TmpQ; << 1577 } << 1578 if ( Q2 > Q1 ) { << 1579 TmpQ = Q1; << 1580 Q1 = Q2; << 1581 Q2 = TmpQ; << 1582 } << 1583 // By now Q1 >= Q2 >= Q3 << 1584 G4int NewCode = Q1*1000 + Q2*100 + Q3*10 + << 1585 return NewCode; << 1586 } 1355 } 1587 1356 1588 1357 1589 //=========================================== << 1358 G4DiffractiveExcitation::~G4DiffractiveExcitation() 1590 << 1359 { 1591 G4DiffractiveExcitation::G4DiffractiveExcitat << 1592 throw G4HadronicException( __FILE__, __LINE << 1593 "G4DiffractiveEx << 1594 } 1360 } 1595 1361 1596 1362 1597 //=========================================== << 1363 const G4DiffractiveExcitation & G4DiffractiveExcitation::operator=(const G4DiffractiveExcitation &) 1598 << 1364 { 1599 const G4DiffractiveExcitation & G4Diffractive << 1365 throw G4HadronicException(__FILE__, __LINE__, "G4DiffractiveExcitation = operator meant to be called"); 1600 throw G4HadronicException( __FILE__, __LINE << 1366 return *this; 1601 "G4DiffractiveEx << 1602 return *this; << 1603 } 1367 } 1604 1368 1605 1369 1606 //=========================================== << 1370 int G4DiffractiveExcitation::operator==(const G4DiffractiveExcitation &) const 1607 << 1371 { 1608 G4bool G4DiffractiveExcitation::operator==( c << 1372 throw G4HadronicException(__FILE__, __LINE__, "G4DiffractiveExcitation == operator meant to be called"); 1609 throw G4HadronicException( __FILE__, __LINE << 1373 return false; 1610 "G4DiffractiveEx << 1611 } 1374 } 1612 1375 1613 << 1376 int G4DiffractiveExcitation::operator!=(const G4DiffractiveExcitation &) const 1614 //=========================================== << 1377 { 1615 << 1378 throw G4HadronicException(__FILE__, __LINE__, "G4DiffractiveExcitation != operator meant to be called"); 1616 G4bool G4DiffractiveExcitation::operator!= ( << 1379 return true; 1617 throw G4HadronicException( __FILE__, __LINE << 1618 "G4DiffractiveEx << 1619 } 1380 } 1620 << 1621 1381