Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer << 3 // * DISCLAIMER * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th << 5 // * The following disclaimer summarizes all the specific disclaimers * 6 // * the Geant4 Collaboration. It is provided << 6 // * of contributors to this software. The specific disclaimers,which * 7 // * conditions of the Geant4 Software License << 7 // * govern, are listed with their locations in: * 8 // * LICENSE and available at http://cern.ch/ << 8 // * http://cern.ch/geant4/license * 9 // * include a list of copyright holders. << 10 // * 9 // * * 11 // * Neither the authors of this software syst 10 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 11 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 12 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 13 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file << 14 // * use. * 16 // * for the full disclaimer and the limitatio << 17 // * 15 // * * 18 // * This code implementation is the result << 16 // * This code implementation is the intellectual property of the * 19 // * technical work of the GEANT4 collaboratio << 17 // * GEANT4 collaboration. * 20 // * By using, copying, modifying or distri << 18 // * By copying, distributing or modifying the Program (or any work * 21 // * any work based on the software) you ag << 19 // * based on the Program) you indicate your acceptance of this * 22 // * use in resulting scientific publicati << 20 // * statement, and all its terms. * 23 // * acceptance of all terms of the Geant4 Sof << 24 // ******************************************* 21 // ******************************************************************** 25 // 22 // 26 // 23 // >> 24 // $Id: G4VLongitudinalStringDecay.cc,v 1.3 2004/12/07 13:50:16 gunter Exp $ >> 25 // GEANT4 tag $Name: geant4-07-00-patch-01 $ 27 // 26 // 28 // ------------------------------------------- 27 // ----------------------------------------------------------------------------- 29 // GEANT 4 class implementation file 28 // GEANT 4 class implementation file 30 // 29 // 31 // History: first implementation, Maxim K 30 // History: first implementation, Maxim Komogorov, 1-Jul-1998 32 // redesign Gunter Folger, Augu 31 // redesign Gunter Folger, August/September 2001 33 // ------------------------------------------- 32 // ----------------------------------------------------------------------------- 34 #include "G4VLongitudinalStringDecay.hh" << 35 #include "G4PhysicalConstants.hh" << 36 #include "G4SystemOfUnits.hh" << 37 #include "G4ios.hh" 33 #include "G4ios.hh" 38 #include "Randomize.hh" 34 #include "Randomize.hh" >> 35 #include "G4VLongitudinalStringDecay.hh" 39 #include "G4FragmentingString.hh" 36 #include "G4FragmentingString.hh" 40 37 41 #include "G4ParticleDefinition.hh" 38 #include "G4ParticleDefinition.hh" 42 #include "G4ParticleTypes.hh" 39 #include "G4ParticleTypes.hh" 43 #include "G4ParticleChange.hh" 40 #include "G4ParticleChange.hh" 44 #include "G4VShortLivedParticle.hh" 41 #include "G4VShortLivedParticle.hh" 45 #include "G4ShortLivedConstructor.hh" 42 #include "G4ShortLivedConstructor.hh" 46 #include "G4ParticleTable.hh" 43 #include "G4ParticleTable.hh" >> 44 #include "G4ShortLivedTable.hh" 47 #include "G4PhaseSpaceDecayChannel.hh" 45 #include "G4PhaseSpaceDecayChannel.hh" 48 #include "G4VDecayChannel.hh" 46 #include "G4VDecayChannel.hh" 49 #include "G4DecayTable.hh" 47 #include "G4DecayTable.hh" 50 48 51 #include "G4DiQuarks.hh" 49 #include "G4DiQuarks.hh" 52 #include "G4Quarks.hh" 50 #include "G4Quarks.hh" 53 #include "G4Gluons.hh" 51 #include "G4Gluons.hh" 54 52 55 #include "G4Exp.hh" << 53 //******************************************************************************** 56 #include "G4Log.hh" << 57 << 58 #include "G4HadronicException.hh" << 59 << 60 //------------------------debug switches << 61 //#define debug_VStringDecay << 62 //#define debug_heavyHadrons << 63 << 64 //******************************************** << 65 // Constructors 54 // Constructors 66 55 67 G4VLongitudinalStringDecay::G4VLongitudinalStr << 56 G4VLongitudinalStringDecay::G4VLongitudinalStringDecay() 68 : G4HadronicInteraction(name), ProbCCbar(0.0 << 69 { 57 { 70 MassCut = 210.0*MeV; // Mpi + Delta << 58 MassCut = 0.35*GeV; >> 59 ClusterMass = 0.15*GeV; 71 60 72 StringLoopInterrupt = 1000; << 61 SmoothParam = 0.9; 73 ClusterLoopInterrupt = 500; << 62 StringLoopInterrupt = 1000; >> 63 ClusterLoopInterrupt = 500; 74 64 75 // Changable Parameters below. << 65 // Changable Parameters below. 76 SigmaQT = 0.5 * GeV; << 77 66 78 StrangeSuppress = 0.44; // =0.27/2.27 s << 67 SigmaQT = 0.5 * GeV; 79 DiquarkSuppress = 0.07; // Probability << 80 DiquarkBreakProb = 0.1; // Probability << 81 68 82 //... pspin_meson is probability to create << 69 StrangeSuppress = 0.44; // 27 % strange quarks produced, ie. u:d:s=1:1:0.27 83 pspin_meson.resize(3); << 70 DiquarkSuppress = 0.1; 84 pspin_meson[0] = 0.5; // u or d + anti-u o << 71 DiquarkBreakProb = 0.1; 85 pspin_meson[1] = 0.4; // one of the quark << 86 pspin_meson[2] = 0.3; // both of the quark << 87 72 88 //... pspin_barion is probability to create << 73 //... pspin_meson is probability to create vector meson >> 74 pspin_meson = 0.5; >> 75 >> 76 //... pspin_barion is probability to create 3/2 barion 89 pspin_barion = 0.5; 77 pspin_barion = 0.5; 90 78 91 //... vectorMesonMix[] is quark mixing para 79 //... vectorMesonMix[] is quark mixing parameters for vector mesons (Variable spin = 3) 92 vectorMesonMix.resize(6); 80 vectorMesonMix.resize(6); 93 vectorMesonMix[0] = 0.0; << 81 vectorMesonMix[0] = 0.5; 94 vectorMesonMix[1] = 0.5; << 82 vectorMesonMix[1] = 0.0; 95 vectorMesonMix[2] = 0.0; << 83 vectorMesonMix[2] = 0.5; 96 vectorMesonMix[3] = 0.5; << 84 vectorMesonMix[3] = 0.0; 97 vectorMesonMix[4] = 1.0; 85 vectorMesonMix[4] = 1.0; 98 vectorMesonMix[5] = 1.0; << 86 vectorMesonMix[5] = 1.0; 99 87 100 //... scalarMesonMix[] is quark mixing para 88 //... scalarMesonMix[] is quark mixing parameters for scalar mesons (Variable spin=1) 101 scalarMesonMix.resize(6); 89 scalarMesonMix.resize(6); 102 scalarMesonMix[0] = 0.5; << 90 scalarMesonMix[0] = 0.5; 103 scalarMesonMix[1] = 0.25; << 91 scalarMesonMix[1] = 0.25; 104 scalarMesonMix[2] = 0.5; << 92 scalarMesonMix[2] = 0.5; 105 scalarMesonMix[3] = 0.25; << 93 scalarMesonMix[3] = 0.25; 106 scalarMesonMix[4] = 1.0; << 94 scalarMesonMix[4] = 1.0; 107 scalarMesonMix[5] = 0.5; << 95 scalarMesonMix[5] = 0.5; 108 << 109 SetProbCCbar(0.0); // Probability of CCbar << 110 SetProbEta_c(0.1); // Mixing of Eta_c and << 111 SetProbBBbar(0.0); // Probability of BBbar << 112 SetProbEta_b(0.0); // Mixing of Eta_b and << 113 96 114 // Parameters may be changed until the firs << 97 // Parameters may be changed until the first fragmentation starts 115 PastInitPhase=false; 98 PastInitPhase=false; 116 hadronizer = new G4HadronBuilder( pspin_mes << 99 hadronizer = new G4HadronBuilder(pspin_meson,pspin_barion, 117 ProbEta_c << 100 scalarMesonMix,vectorMesonMix); 118 << 119 MaxMass=-350.0*GeV; // If there will be a << 120 << 121 SetMinMasses(); // Re-calculation of minim << 122 << 123 Kappa = 1.0 * GeV/fermi; << 124 DecayQuark = NewQuark = 0; << 125 } 101 } >> 102 126 103 127 G4VLongitudinalStringDecay::~G4VLongitudinalSt 104 G4VLongitudinalStringDecay::~G4VLongitudinalStringDecay() 128 { << 105 { 129 delete hadronizer; 106 delete hadronizer; 130 } << 107 } 131 << 132 G4HadFinalState* << 133 G4VLongitudinalStringDecay::ApplyYourself(cons << 134 { << 135 return nullptr; << 136 } << 137 << 138 //============================================ << 139 << 140 // For changing Mass Cut used for selection of << 141 void G4VLongitudinalStringDecay::SetMassCut(G4 << 142 G4double G4VLongitudinalStringDecay::GetMassCu << 143 << 144 //-------------------------------------------- << 145 << 146 // For handling a string with very low mass << 147 << 148 G4KineticTrackVector* G4VLongitudinalStringDec << 149 { << 150 G4KineticTrackVector* result = nullptr << 151 pDefPair hadrons( nullptr, nullptr ); << 152 G4FragmentingString aString( *string ) << 153 << 154 #ifdef debug_VStringDecay << 155 G4cout<<"G4VLongitudinalStringDecay::P << 156 <<aString.Mass()<<" MassCut "<<M << 157 #endif << 158 << 159 SetMinimalStringMass( &aString ); << 160 PossibleHadronMass( &aString, 0, &hadr << 161 result = new G4KineticTrackVector; << 162 if ( hadrons.first != nullptr ) { << 163 if ( hadrons.second == nullptr ) { << 164 // Substitute string by light h << 165 << 166 #ifdef debug_VStringDecay << 167 G4cout << "VlongSD Warning repl << 168 G4cout << hadrons.first->GetPar << 169 << "string .. " << strin << 170 << string->Get4Momentum( << 171 #endif << 172 << 173 G4ThreeVector Mom3 = string-> << 174 G4LorentzVector Mom( Mom3, std: << 175 result->push_back( new G4Kineti << 176 } else { << 177 //... string was qq--qqbar type << 178 << 179 #ifdef debug_VStringDecay << 180 G4cout << "VlongSD Warning repl << 181 << hadrons.first->GetPar << 182 << hadrons.second->GetPa << 183 << "string .. " << strin << 184 << string->Get4Momentum( << 185 #endif << 186 << 187 G4LorentzVector Mom1, Mom2; << 188 Sample4Momentum( &Mom1, hadrons << 189 &Mom2, hadrons << 190 string->Get4Mo << 191 << 192 result->push_back( new G4Kineti << 193 result->push_back( new G4Kineti << 194 << 195 G4ThreeVector Velocity = string << 196 result->Boost(Velocity); << 197 } << 198 } << 199 return result; << 200 } << 201 << 202 //-------------------------------------------- << 203 << 204 G4double G4VLongitudinalStringDecay::PossibleH << 205 << 206 { << 207 G4double mass = 0.0; << 208 << 209 if ( build==0 ) build=&G4HadronBuilder::Buil << 210 108 211 G4ParticleDefinition* Hadron1 = nullpt << 109 //=============================================================================================------------- 212 G4ParticleDefinition* Hadron2 = nullptr; << 213 110 214 if (!string->IsAFourQuarkString() ) << 111 // Operators 215 { << 216 // spin 0 meson or spin 1/2 barion << 217 112 218 Hadron1 = (hadronizer->*build)(stri << 113 //const & G4VLongitudinalStringDecay::operator=(const G4VLongitudinalStringDecay &) 219 #ifdef debug_VStringDecay << 114 // { 220 G4cout<<"VlongSD PossibleHadronMass"<<G4e << 115 // } 221 G4cout<<"VlongSD Quarks at the stri << 222 <<" "<<string->GetRightParton << 223 if ( Hadron1 != nullptr) { << 224 G4cout<<"(G4VLongitudinalStringDe << 225 <<" "<<Hadron1->GetPDGMass( << 226 } << 227 #endif << 228 if ( Hadron1 != nullptr) { mass = ( << 229 else { mass = MaxM << 230 } else << 231 { << 232 //... string is qq--qqbar: Build tw << 233 << 234 #ifdef debug_VStringDecay << 235 G4cout<<"VlongSD PossibleHadronMass << 236 G4cout<<"VlongSD string is qq--qqba << 237 #endif << 238 << 239 G4double StringMass = string->Mas << 240 G4int cClusterInterrupt = 0; << 241 do << 242 { << 243 if (cClusterInterrupt++ >= Cluste << 244 << 245 G4int LeftQuark1= string->GetLeft << 246 G4int LeftQuark2=(string->GetLeft << 247 << 248 G4int RightQuark1= string->GetRig << 249 G4int RightQuark2=(string->GetRig << 250 << 251 if (G4UniformRand()<0.5) { << 252 Hadron1 =hadronizer->Build(Find << 253 Hadron2 =hadronizer->Build(Find << 254 } else { << 255 Hadron1 =hadronizer->Build(Find << 256 Hadron2 =hadronizer->Build(Find << 257 } << 258 //... repeat procedure, if mass o << 259 //... ClusterMassCut = 0.15*GeV m << 260 } << 261 while ( Hadron1 == nullptr || Hadro << 262 ( StringMass <= Hadron1->Ge << 263 116 264 mass = (Hadron1)->GetPDGMass() + (Hadron2 << 117 //---------------------------------------------------------------------------------------------------------- 265 } << 266 118 267 #ifdef debug_VStringDecay << 119 int G4VLongitudinalStringDecay::operator==(const G4VLongitudinalStringDecay &) const 268 G4cout<<"VlongSD *Hadrons 1 and 2, pro << 120 { 269 #endif << 121 throw G4HadronicException(__FILE__, __LINE__, "G4VLongitudinalStringDecay::operator== forbidden"); 270 << 122 return false; 271 if ( pdefs != 0 ) << 123 } 272 { // need to return hadrons as well.... << 273 pdefs->first = Hadron1; << 274 pdefs->second = Hadron2; << 275 } << 276 << 277 return mass; << 278 } << 279 124 280 //-------------------------------------------- << 125 //---------------------------------------------------------------------------------------------------------- 281 126 282 G4ParticleDefinition* G4VLongitudinalStringDec << 127 int G4VLongitudinalStringDecay::operator!=(const G4VLongitudinalStringDecay &) const 283 { << 128 { 284 /* << 129 throw G4HadronicException(__FILE__, __LINE__, "G4VLongitudinalStringDecay::operator!= forbidden"); 285 G4cout<<Encoding<<" G4VLongitudinalStringDec << 130 return true; 286 for (G4int i=4; i<6;i++){ << 287 for (G4int j=1;j<6;j++){ << 288 G4cout<<i<<" "<<j<<" "; << 289 G4int Code = 1000 * i + 100 * j +1; << 290 G4ParticleDefinition* ptr1 = G4ParticleT << 291 Code +=2; << 292 G4ParticleDefinition* ptr2 = G4ParticleT << 293 G4cout<<"Code "<<Code - 2<<" ptr "<<ptr1 << 294 } 131 } 295 G4cout<<G4endl; << 296 } << 297 */ << 298 << 299 G4ParticleDefinition* ptr = G4ParticleTable: << 300 << 301 if (ptr == nullptr) << 302 { << 303 for (size_t i=0; i < NewParticles.size(); << 304 { << 305 if ( Encoding == NewParticles[i]->GetPD << 306 } << 307 } << 308 << 309 return ptr; << 310 } << 311 132 312 //******************************************** << 133 //========================================================================================================== 313 // For decision on continue or stop string f << 314 // virtual G4bool StopFragmenting(const G4Fr << 315 // virtual G4bool IsItFragmentable(const G4F << 316 // << 317 // If a string can not fragment, make last b << 318 // virtual G4bool SplitLast(G4FragmentingStr << 319 // G4KineticTrackVe << 320 // G4KineticTrackVe << 321 //-------------------------------------------- << 322 // << 323 // If a string can fragment, do the followin << 324 // << 325 // For transver of a string to its CMS frame << 326 //-------------------------------------------- << 327 134 328 G4ExcitedString *G4VLongitudinalStringDecay::C << 135 G4int G4VLongitudinalStringDecay::SampleQuarkFlavor(void) 329 { << 136 { 330 G4Parton *Left=new G4Parton(*in.GetLeftParto << 137 return (1 + (int)(G4UniformRand()/StrangeSuppress)); 331 G4Parton *Right=new G4Parton(*in.GetRightPar << 138 } 332 return new G4ExcitedString(Left,Right,in.Get << 333 } << 334 139 335 //-------------------------------------------- << 140 //---------------------------------------------------------------------------------------------------------- 336 141 337 G4ParticleDefinition * G4VLongitudinalStringDe << 142 G4VLongitudinalStringDecay::pDefPair G4VLongitudinalStringDecay::CreatePartonPair(G4int NeedParticle,G4bool AllowDiquarks) 338 << 339 { 143 { 340 #ifdef debug_VStringDecay << 144 // NeedParticle = +1 for Particle, -1 for Antiparticle 341 G4cout<<"VlongSD QuarkSplitup: quark ID "<< << 342 #endif << 343 << 344 G4int IsParticle=(decay->GetPDGEncoding()>0 << 345 << 346 pDefPair QuarkPair = CreatePartonPair(IsPar << 347 created = QuarkPair.second; << 348 << 349 DecayQuark = decay->GetPDGEncoding(); << 350 NewQuark = created->GetPDGEncoding(); << 351 << 352 #ifdef debug_VStringDecay << 353 G4cout<<"VlongSD QuarkSplitup: "<<decay->Ge << 354 G4cout<<"hadronizer->Build(QuarkPair.first, << 355 #endif << 356 << 357 return hadronizer->Build(QuarkPair.first, d << 358 } << 359 << 360 //-------------------------------------------- << 361 145 362 G4VLongitudinalStringDecay::pDefPair G4VLongit << 363 CreatePartonPair(G4int NeedParticle,G4bool All << 364 { << 365 // NeedParticle = +1 for Particle, -1 for << 366 if ( AllowDiquarks && G4UniformRand() < Di 146 if ( AllowDiquarks && G4UniformRand() < DiquarkSuppress ) 367 { 147 { 368 // Create a Diquark - AntiDiquark pair , 148 // Create a Diquark - AntiDiquark pair , first in pair is anti to IsParticle 369 #ifdef debug_VStringDecay << 149 G4int q1 = SampleQuarkFlavor(); 370 G4cout<<"VlongSD Create a Diquark - Anti << 150 G4int q2 = SampleQuarkFlavor(); 371 #endif << 151 G4int spin = (q1 != q2 && G4UniformRand() <= 0.5)? 1 : 3; 372 G4int q1(0), q2(0), spin(0), PDGcode(0); << 373 << 374 q1 = SampleQuarkFlavor(); << 375 q2 = SampleQuarkFlavor(); << 376 << 377 spin = (q1 != q2 && G4UniformRand() <= 0 << 378 // conv 152 // convention: quark with higher PDG number is first 379 PDGcode = (std::max(q1,q2) * 1000 + std: << 153 G4int PDGcode = (std::max(q1,q2) * 1000 + std::min(q1,q2) * 100 + spin) * NeedParticle; 380 << 381 return pDefPair (FindParticle(-PDGcode), 154 return pDefPair (FindParticle(-PDGcode),FindParticle(PDGcode)); >> 155 382 156 383 } else { 157 } else { 384 // Create a Quark - AntiQuark pair, firs 158 // Create a Quark - AntiQuark pair, first in pair IsParticle 385 #ifdef debug_VStringDecay << 386 G4cout<<"VlongSD Create a Quark - AntiQu << 387 #endif << 388 G4int PDGcode=SampleQuarkFlavor()*NeedPa 159 G4int PDGcode=SampleQuarkFlavor()*NeedParticle; 389 return pDefPair (FindParticle(PDGcode),F 160 return pDefPair (FindParticle(PDGcode),FindParticle(-PDGcode)); 390 } 161 } 391 } << 392 162 393 //-------------------------------------------- << 394 << 395 G4int G4VLongitudinalStringDecay::SampleQuarkF << 396 { << 397 G4int quark(1); << 398 G4double ksi = G4UniformRand(); << 399 if ( ksi < ProbCB ) { << 400 if ( ksi < ProbCCbar ) {quark = 4;} // << 401 else {quark = 5;} // << 402 #ifdef debug_heavyHadrons << 403 G4cout << "G4VLongitudinalStringDecay::S << 404 << quark << G4endl; << 405 #endif << 406 } else { << 407 quark = 1 + (int)(G4UniformRand()/Strange << 408 } << 409 #ifdef debug_VStringDecay << 410 G4cout<<"VlongSD SampleQuarkFlavor "<<quark << 411 <<" "<<ProbCCbar<<" "<<ProbBBbar<<" ) << 412 #endif << 413 return quark; << 414 } 163 } 415 164 416 //-------------------------------------------- << 165 //---------------------------------------------------------------------------------------------------------- 417 166 418 G4ThreeVector G4VLongitudinalStringDecay::Samp << 167 G4ThreeVector G4VLongitudinalStringDecay::SampleQuarkPt() 419 { << 168 { 420 G4double Pt; << 169 G4double Pt = -std::log(G4UniformRand()); 421 if ( ptMax < 0 ) { << 422 // sample full gaussian << 423 Pt = -G4Log(G4UniformRand()); << 424 } else { << 425 // sample in limited range << 426 G4double q = ptMax/SigmaQT; << 427 G4double ymin = (q > 20.) ? 0.0 : G4Exp( << 428 Pt = -G4Log(G4RandFlat::shoot(ymin, 1.)) << 429 } << 430 Pt = SigmaQT * std::sqrt(Pt); 170 Pt = SigmaQT * std::sqrt(Pt); 431 G4double phi = 2.*pi*G4UniformRand(); 171 G4double phi = 2.*pi*G4UniformRand(); 432 return G4ThreeVector(Pt * std::cos(phi),Pt 172 return G4ThreeVector(Pt * std::cos(phi),Pt * std::sin(phi),0); 433 } << 173 } 434 174 435 //******************************************** << 175 //---------------------------------------------------------------------------------------------------------- 436 176 437 void G4VLongitudinalStringDecay::CalculateHadr << 177 void G4VLongitudinalStringDecay::CalculateHadronTimePosition(G4double theInitialStringMass, G4KineticTrackVector* Hadrons) 438 << 439 { << 440 // `yo-yo` formation time << 441 // const G4double kappa = 1.0 * GeV/fermi << 442 G4double kappa = GetStringTensionParameter( << 443 for (size_t c1 = 0; c1 < Hadrons->size(); c << 444 { 178 { >> 179 // `yo-yo` formation time >> 180 const G4double kappa = 1.0 * GeV/fermi; >> 181 for(size_t c1 = 0; c1 < Hadrons->size(); c1++) >> 182 { 445 G4double SumPz = 0; 183 G4double SumPz = 0; 446 G4double SumE = 0; 184 G4double SumE = 0; 447 for (size_t c2 = 0; c2 < c1; c2++) << 185 for(size_t c2 = 0; c2 < c1; c2++) 448 { << 186 { 449 SumPz += Hadrons->operator[](c2)->Get 187 SumPz += Hadrons->operator[](c2)->Get4Momentum().pz(); 450 SumE += Hadrons->operator[](c2)->Get 188 SumE += Hadrons->operator[](c2)->Get4Momentum().e(); 451 } << 189 } 452 G4double HadronE = Hadrons->operator[]( 190 G4double HadronE = Hadrons->operator[](c1)->Get4Momentum().e(); 453 G4double HadronPz = Hadrons->operator[]( 191 G4double HadronPz = Hadrons->operator[](c1)->Get4Momentum().pz(); 454 Hadrons->operator[](c1)->SetFormationTim << 192 Hadrons->operator[](c1)->SetFormationTime((theInitialStringMass - 2.*SumPz + HadronE - HadronPz)/(2.*kappa)); 455 (theInitialStringMass - 2.*SumPz + Had << 193 G4ThreeVector aPosition(0, 0, (theInitialStringMass - 2.*SumE - HadronE + HadronPz)/(2.*kappa)); 456 G4ThreeVector aPosition( 0, 0, << 457 (theInitialStringMass - 2.*SumE - Had << 458 Hadrons->operator[](c1)->SetPosition(aPo 194 Hadrons->operator[](c1)->SetPosition(aPosition); >> 195 } 459 } 196 } 460 } << 461 197 462 //-------------------------------------------- << 198 //---------------------------------------------------------------------------------------------------------- 463 199 464 void G4VLongitudinalStringDecay::SetSigmaTrans << 200 /* 465 { << 201 void G4VLongitudinalStringDecay::CalculateHadronTimePosition(G4double theInitialStringMass, G4KineticTrackVector* Hadrons) 466 if ( PastInitPhase ) { << 202 { 467 throw G4HadronicException(__FILE__, __LIN << 203 // 'constituent' formation time 468 "G4VLongitudinalStringDecay::SetSigmaTr << 204 const G4double kappa = 1.0 * GeV/fermi; 469 } else { << 205 for(G4int c1 = 0; c1 < Hadrons->length(); c1++) 470 SigmaQT = aValue; << 206 { >> 207 G4double SumPz = 0; >> 208 G4double SumE = 0; >> 209 for(G4int c2 = 0; c2 <= c1; c2++) >> 210 { >> 211 SumPz += Hadrons->at(c2)->Get4Momentum().pz(); >> 212 SumE += Hadrons->at(c2)->Get4Momentum().e(); >> 213 } >> 214 Hadrons->at(c1)->SetFormationTime((theInitialStringMass - 2.*SumPz)/(2.*kappa)); >> 215 G4ThreeVector aPosition(0, 0, (theInitialStringMass - 2.*SumE)/(2.*kappa)); >> 216 Hadrons->at(c1)->SetPosition(aPosition); >> 217 } >> 218 c1 = Hadrons->length()-1; >> 219 Hadrons->at(c1)->SetFormationTime(Hadrons->at(c1-1)->GetFormationTime()); >> 220 Hadrons->at(c1)->SetPosition(Hadrons->at(c1-1)->GetPosition()); 471 } 221 } 472 } << 222 */ 473 223 474 //-------------------------------------------- 224 //---------------------------------------------------------------------------------------------------------- 475 225 476 void G4VLongitudinalStringDecay::SetStrangenes << 226 G4ParticleDefinition * 477 { << 227 G4VLongitudinalStringDecay::QuarkSplitup(G4ParticleDefinition* 478 StrangeSuppress = aValue; << 228 decay, G4ParticleDefinition *&created) >> 229 { >> 230 G4int IsParticle=(decay->GetPDGEncoding()>0) ? -1 : +1; // if we have a quark, we need antiquark (or diquark) >> 231 pDefPair QuarkPair = CreatePartonPair(IsParticle); >> 232 created = QuarkPair.second; >> 233 return hadronizer->Build(QuarkPair.first, decay); >> 234 479 } 235 } 480 236 481 //-------------------------------------------- 237 //---------------------------------------------------------------------------------------------------------- 482 238 483 void G4VLongitudinalStringDecay::SetDiquarkSup << 239 G4ParticleDefinition *G4VLongitudinalStringDecay::DiQuarkSplitup( 484 { << 240 G4ParticleDefinition* decay, 485 DiquarkSuppress = aValue; << 241 G4ParticleDefinition *&created) 486 } << 242 { 487 << 243 //... can Diquark break or not? 488 //-------------------------------------------- << 244 if (G4UniformRand() < DiquarkBreakProb ){ >> 245 //... Diquark break >> 246 >> 247 G4int stableQuarkEncoding = decay->GetPDGEncoding()/1000; >> 248 G4int decayQuarkEncoding = (decay->GetPDGEncoding()/100)%10; >> 249 if (G4UniformRand() < 0.5) >> 250 { >> 251 G4int Swap = stableQuarkEncoding; >> 252 stableQuarkEncoding = decayQuarkEncoding; >> 253 decayQuarkEncoding = Swap; >> 254 } >> 255 >> 256 G4int IsParticle=(decayQuarkEncoding>0) ? -1 : +1; >> 257 // if we have a quark, we need antiquark) >> 258 pDefPair QuarkPair = CreatePartonPair(IsParticle,false); // no diquarks wanted >> 259 //... Build new Diquark >> 260 G4int QuarkEncoding=QuarkPair.second->GetPDGEncoding(); >> 261 G4int i10 = std::max(std::abs(QuarkEncoding), std::abs(stableQuarkEncoding)); >> 262 G4int i20 = std::min(std::abs(QuarkEncoding), std::abs(stableQuarkEncoding)); >> 263 G4int spin = (i10 != i20 && G4UniformRand() <= 0.5)? 1 : 3; >> 264 G4int NewDecayEncoding = -1*IsParticle*(i10 * 1000 + i20 * 100 + spin); >> 265 created = FindParticle(NewDecayEncoding); >> 266 G4ParticleDefinition * decayQuark=FindParticle(decayQuarkEncoding); >> 267 >> 268 return hadronizer->Build(QuarkPair.first, decayQuark); >> 269 >> 270 } else { >> 271 //... Diquark does not break >> 272 >> 273 G4int IsParticle=(decay->GetPDGEncoding()>0) ? +1 : -1; >> 274 // if we have a diquark, we need quark) >> 275 pDefPair QuarkPair = CreatePartonPair(IsParticle,false); // no diquarks wanted >> 276 created = QuarkPair.second; 489 277 490 void G4VLongitudinalStringDecay::SetDiquarkBre << 278 return hadronizer->Build(QuarkPair.first, decay); 491 { << 279 } 492 if ( PastInitPhase ) { << 493 throw G4HadronicException(__FILE__, __LINE << 494 "G4VLongitudinalStringDecay::SetDiquarkB << 495 } else { << 496 DiquarkBreakProb = aValue; << 497 } << 498 } 280 } 499 281 500 //-------------------------------------------- << 282 //----------------------------------------------------------------------------------------- 501 283 502 void G4VLongitudinalStringDecay::SetSpinThreeH << 284 G4KineticTrack * G4VLongitudinalStringDecay::Splitup( >> 285 G4FragmentingString *string, >> 286 G4FragmentingString *&newString) 503 { 287 { 504 if ( PastInitPhase ) { << 505 throw G4HadronicException(__FILE__, __LINE << 506 "G4VLongitudinalStringDecay::SetSpinThre << 507 } else { << 508 pspin_barion = aValue; << 509 delete hadronizer; << 510 hadronizer = new G4HadronBuilder( pspin_me << 511 ProbEta_ << 512 } << 513 } << 514 << 515 //-------------------------------------------- << 516 288 517 void G4VLongitudinalStringDecay::SetScalarMeso << 289 //... random choice of string end to use for creating the hadron (decay) 518 { << 290 SideOfDecay = (G4UniformRand() < 0.5)? 1: -1; 519 if ( PastInitPhase ) { << 291 if (SideOfDecay < 0) 520 throw G4HadronicException(__FILE__, __LINE << 292 { 521 "G4VLongitudinalStringDecay::SetScalarMe << 293 string->SetLeftPartonStable(); 522 } else { << 294 } else 523 if ( aVector.size() < 6 ) << 295 { 524 throw G4HadronicException(__FILE__, __LI << 296 string->SetRightPartonStable(); 525 "G4VLongitudinalStringDecay::SetScalar << 297 } 526 scalarMesonMix[0] = aVector[0]; << 527 scalarMesonMix[1] = aVector[1]; << 528 scalarMesonMix[2] = aVector[2]; << 529 scalarMesonMix[3] = aVector[3]; << 530 scalarMesonMix[4] = aVector[4]; << 531 scalarMesonMix[5] = aVector[5]; << 532 delete hadronizer; << 533 hadronizer = new G4HadronBuilder( pspin_me << 534 ProbEta_ << 535 } << 536 } << 537 298 538 //-------------------------------------------- << 299 G4ParticleDefinition *newStringEnd; >> 300 G4ParticleDefinition * HadronDefinition; >> 301 if (string->DecayIsQuark()) >> 302 { >> 303 HadronDefinition= QuarkSplitup(string->GetDecayParton(), newStringEnd); >> 304 } else { >> 305 HadronDefinition= DiQuarkSplitup(string->GetDecayParton(), newStringEnd); >> 306 } >> 307 // create new String from old, ie. keep Left and Right order, but replace decay >> 308 G4LorentzVector* HadronMomentum=SplitEandP(HadronDefinition, string); >> 309 >> 310 G4KineticTrack * Hadron =0; >> 311 if ( HadronMomentum != 0 ) { 539 312 540 void G4VLongitudinalStringDecay::SetVectorMeso << 313 G4ThreeVector Pos; 541 { << 314 Hadron = new G4KineticTrack(HadronDefinition, 0,Pos, *HadronMomentum); 542 if ( PastInitPhase ) { << 315 543 throw G4HadronicException(__FILE__, __LINE << 316 newString=new G4FragmentingString(*string,newStringEnd, 544 "G4VLongitudinalStringDecay::SetVectorMe << 317 HadronMomentum); 545 } else { << 318 546 if ( aVector.size() < 6 ) << 319 delete HadronMomentum; 547 throw G4HadronicException(__FILE__, __LI << 320 } 548 "G4VLongitudinalStringDecay::SetVector << 321 return Hadron; 549 vectorMesonMix[0] = aVector[0]; << 322 } 550 vectorMesonMix[1] = aVector[1]; << 323 551 vectorMesonMix[2] = aVector[2]; << 324 //----------------------------------------------------------------------------------------- 552 vectorMesonMix[3] = aVector[3]; << 325 553 vectorMesonMix[4] = aVector[4]; << 326 G4LorentzVector * G4VLongitudinalStringDecay::SplitEandP(G4ParticleDefinition * pHadron, 554 vectorMesonMix[5] = aVector[5]; << 327 G4FragmentingString * string) 555 delete hadronizer; << 328 { 556 hadronizer = new G4HadronBuilder( pspin_me << 329 G4double HadronMass = pHadron->GetPDGMass(); 557 ProbEta_ << 330 558 } << 331 // calculate and assign hadron transverse momentum component HadronPx andHadronPy >> 332 G4ThreeVector thePt; >> 333 thePt=SampleQuarkPt(); >> 334 G4ThreeVector HadronPt = thePt +string->DecayPt(); >> 335 HadronPt.setZ(0); >> 336 //... sample z to define hadron longitudinal momentum and energy >> 337 //... but first check the available phase space >> 338 G4double DecayQuarkMass2 = sqr(string->GetDecayParton()->GetPDGMass()); >> 339 G4double HadronMass2T = sqr(HadronMass) + HadronPt.mag2(); >> 340 if (DecayQuarkMass2 + HadronMass2T >= SmoothParam*(string->Mass2()) ) >> 341 return 0; // have to start all over! >> 342 >> 343 //... then compute allowed z region z_min <= z <= z_max >> 344 >> 345 G4double zMin = HadronMass2T/(string->Mass2()); >> 346 G4double zMax = 1. - DecayQuarkMass2/(string->Mass2()); >> 347 if (zMin >= zMax) return 0; // have to start all over! >> 348 >> 349 G4double z = GetLightConeZ(zMin, zMax, >> 350 string->GetDecayParton()->GetPDGEncoding(), pHadron, >> 351 HadronPt.x(), HadronPt.y()); >> 352 >> 353 //... now compute hadron longitudinal momentum and energy >> 354 // longitudinal hadron momentum component HadronPz >> 355 >> 356 HadronPt.setZ(0.5* string->GetDecayDirection() * >> 357 (z * string->LightConeDecay() - >> 358 HadronMass2T/(z * string->LightConeDecay()))); >> 359 G4double HadronE = 0.5* (z * string->LightConeDecay() + >> 360 HadronMass2T/(z * string->LightConeDecay())); >> 361 >> 362 G4LorentzVector * a4Momentum= new G4LorentzVector(HadronPt,HadronE); >> 363 >> 364 return a4Momentum; 559 } 365 } 560 366 561 //-------------------------------------------- << 562 << 563 void G4VLongitudinalStringDecay::SetProbCCbar( << 564 { << 565 ProbCCbar = aValue; << 566 ProbCB = ProbCCbar + ProbBBbar; << 567 } << 568 367 569 //-------------------------------------------- << 368 //----------------------------------------------------------------------------------------- >> 369 >> 370 G4bool G4VLongitudinalStringDecay::SplitLast(G4FragmentingString * string, >> 371 G4KineticTrackVector * LeftVector, >> 372 G4KineticTrackVector * RightVector) >> 373 { >> 374 //... perform last cluster decay >> 375 G4ThreeVector ClusterVel =string->Get4Momentum().boostVector(); >> 376 G4double ResidualMass =string->Mass(); >> 377 G4double ClusterMassCut = ClusterMass; >> 378 G4int cClusterInterrupt = 0; >> 379 G4ParticleDefinition * LeftHadron, * RightHadron; >> 380 do >> 381 { >> 382 if (cClusterInterrupt++ >= ClusterLoopInterrupt) >> 383 { >> 384 return false; >> 385 } >> 386 G4ParticleDefinition * quark = NULL; >> 387 string->SetLeftPartonStable(); // to query quark contents.. >> 388 if (string->DecayIsQuark() && string->StableIsQuark() ) >> 389 { >> 390 //... there are quarks on cluster ends >> 391 LeftHadron= QuarkSplitup(string->GetLeftParton(), quark); >> 392 } else { >> 393 //... there is a Diquark on cluster ends >> 394 G4int IsParticle; >> 395 if ( string->StableIsQuark() ) { >> 396 IsParticle=(string->GetLeftParton()->GetPDGEncoding()>0) ? -1 : +1; >> 397 } else { >> 398 IsParticle=(string->GetLeftParton()->GetPDGEncoding()>0) ? +1 : -1; >> 399 } >> 400 pDefPair QuarkPair = CreatePartonPair(IsParticle,false); // no diquarks wanted >> 401 quark = QuarkPair.second; >> 402 LeftHadron=hadronizer->Build(QuarkPair.first, string->GetLeftParton()); >> 403 } >> 404 RightHadron = hadronizer->Build(string->GetRightParton(), quark); 570 405 571 void G4VLongitudinalStringDecay::SetProbEta_c( << 406 //... repeat procedure, if mass of cluster is too low to produce hadrons 572 { << 407 //... ClusterMassCut = 0.15*GeV model parameter 573 ProbEta_c = aValue; << 408 if ( quark->GetParticleSubType()== "quark" ) {ClusterMassCut = 0.;} 574 } << 409 else {ClusterMassCut = ClusterMass;} >> 410 } >> 411 while (ResidualMass <= LeftHadron->GetPDGMass() + RightHadron->GetPDGMass() + ClusterMassCut); >> 412 >> 413 //... compute hadron momenta and energies >> 414 G4LorentzVector LeftMom, RightMom; >> 415 G4ThreeVector Pos; >> 416 Sample4Momentum(&LeftMom, LeftHadron->GetPDGMass(), &RightMom, RightHadron->GetPDGMass(), ResidualMass); >> 417 LeftMom.boost(ClusterVel); >> 418 RightMom.boost(ClusterVel); >> 419 LeftVector->push_back(new G4KineticTrack(LeftHadron, 0, Pos, LeftMom)); >> 420 RightVector->push_back(new G4KineticTrack(RightHadron, 0, Pos, RightMom)); 575 421 576 //-------------------------------------------- << 422 return true; 577 423 578 void G4VLongitudinalStringDecay::SetProbBBbar( << 579 { << 580 ProbBBbar = aValue; << 581 ProbCB = ProbCCbar + ProbBBbar; << 582 } 424 } 583 425 584 //-------------------------------------------- << 426 //---------------------------------------------------------------------------------------------------------- 585 427 586 void G4VLongitudinalStringDecay::SetProbEta_b( << 428 G4KineticTrackVector* G4VLongitudinalStringDecay::FragmentString(const G4ExcitedString& theString) 587 { 429 { 588 ProbEta_b = aValue; << 430 // Can no longer modify Parameters for Fragmentation. 589 } << 431 PastInitPhase=true; >> 432 >> 433 // check if string has enough mass to fragment... >> 434 G4KineticTrackVector * LeftVector=LightFragmentationTest(&theString); >> 435 if ( LeftVector != 0 ) return LeftVector; >> 436 >> 437 LeftVector = new G4KineticTrackVector; >> 438 G4KineticTrackVector * RightVector=new G4KineticTrackVector; 590 439 591 //-------------------------------------------- << 440 // this should work but its only a semi deep copy. %GF G4ExcitedString theStringInCMS(theString); >> 441 G4ExcitedString *theStringInCMS=CPExcited(theString); >> 442 G4LorentzRotation toCms=theStringInCMS->TransformToAlignedCms(); >> 443 >> 444 G4bool success=false, inner_sucess=true; >> 445 G4int attempt=0; >> 446 while ( !success && attempt++ < StringLoopInterrupt ) >> 447 { >> 448 G4FragmentingString *currentString=new G4FragmentingString(*theStringInCMS); >> 449 >> 450 std::for_each(LeftVector->begin(), LeftVector->end(), DeleteKineticTrack()); >> 451 LeftVector->clear(); >> 452 std::for_each(RightVector->begin(), RightVector->end(), DeleteKineticTrack()); >> 453 RightVector->clear(); >> 454 >> 455 inner_sucess=true; // set false on failure.. >> 456 while (! StopFragmenting(currentString) ) >> 457 { // Split current string into hadron + new string >> 458 G4FragmentingString *newString=0; // used as output from SplitUp... >> 459 G4KineticTrack * Hadron=Splitup(currentString,newString); >> 460 if ( Hadron != 0 && IsFragmentable(newString)) >> 461 { >> 462 if ( currentString->GetDecayDirection() > 0 ) >> 463 LeftVector->push_back(Hadron); >> 464 else >> 465 RightVector->push_back(Hadron); >> 466 delete currentString; >> 467 currentString=newString; >> 468 } else { >> 469 // abandon ... start from the beginning >> 470 if (newString) delete newString; >> 471 if (Hadron) delete Hadron; >> 472 inner_sucess=false; >> 473 break; >> 474 } >> 475 } >> 476 // Split current string into 2 final Hadrons >> 477 if ( inner_sucess && >> 478 SplitLast(currentString,LeftVector, RightVector) ) >> 479 { >> 480 success=true; >> 481 } >> 482 delete currentString; >> 483 } >> 484 >> 485 delete theStringInCMS; >> 486 >> 487 if ( ! success ) >> 488 { >> 489 std::for_each(LeftVector->begin(), LeftVector->end(), DeleteKineticTrack()); >> 490 LeftVector->clear(); >> 491 std::for_each(RightVector->begin(), RightVector->end(), DeleteKineticTrack()); >> 492 delete RightVector; >> 493 return LeftVector; >> 494 } >> 495 >> 496 // Join Left- and RightVector into LeftVector in correct order. >> 497 while(!RightVector->empty()) >> 498 { >> 499 LeftVector->push_back(RightVector->back()); >> 500 RightVector->erase(RightVector->end()-1); >> 501 } >> 502 delete RightVector; 592 503 593 void G4VLongitudinalStringDecay::SetStringTens << 504 CalculateHadronTimePosition(theString.Get4Momentum().mag(), LeftVector); 594 { << 595 Kappa = aValue * GeV/fermi; << 596 } << 597 505 598 //-------------------------------------------- << 506 G4LorentzRotation toObserverFrame(toCms.inverse()); 599 507 600 void G4VLongitudinalStringDecay::SetMinMasses( << 508 for(size_t C1 = 0; C1 < LeftVector->size(); C1++) 601 { << 509 { 602 // ------ For estimation of a minimal stri << 510 G4KineticTrack* Hadron = LeftVector->operator[](C1); 603 Mass_of_light_quark =140.*MeV; << 511 G4LorentzVector Momentum = Hadron->Get4Momentum(); 604 Mass_of_s_quark =500.*MeV; << 512 Momentum = toObserverFrame*Momentum; 605 Mass_of_c_quark =1600.*MeV; << 513 Hadron->Set4Momentum(Momentum); 606 Mass_of_b_quark =4500.*MeV; << 514 G4LorentzVector Coordinate(Hadron->GetPosition(), Hadron->GetFormationTime()); 607 Mass_of_string_junction=720.*MeV; << 515 Momentum = toObserverFrame*Coordinate; 608 << 516 Hadron->SetFormationTime(Momentum.e()); 609 // ---------------- Determination of minim << 517 G4ThreeVector aPosition(Momentum.vect()); 610 G4ParticleDefinition * hadron1; G4int C << 518 Hadron->SetPosition(theString.GetPosition()+aPosition); 611 G4ParticleDefinition * hadron2; G4int C << 519 } 612 for (G4int i=1; i < 6; i++) { << 520 return LeftVector; 613 Code1 = 100*i + 10*1 + 1; << 521 614 hadron1 = FindParticle(Code1); << 615 << 616 if (hadron1 != nullptr) { << 617 for (G4int j=1; j < 6; j++) { << 618 Code2 = 100*j + 10*1 + 1; << 619 hadron2 = FindParticle(Code2); << 620 if (hadron2 != nullptr) { << 621 minMassQQbarStr[i-1][j-1] = h << 622 } << 623 } << 624 } << 625 } << 626 << 627 minMassQQbarStr[1][1] = minMassQQbarStr[0] << 628 522 629 // ---------------- Determination of minim << 630 G4ParticleDefinition * hadron3; << 631 G4int kfla, kflb; << 632 // MaxMass = -350.0*GeV; // If there wi << 633 523 634 for (G4int i=1; i < 6; i++) { //i=1 << 524 } 635 Code1 = 100*i + 10*1 + 1; << 636 hadron1 = FindParticle(Code1); << 637 for (G4int j=1; j < 6; j++) { << 638 for (G4int k=1; k < 6; k++) { << 639 kfla = std::max(j,k); << 640 kflb = std::min(j,k); << 641 525 642 // Add d-quark << 526 //---------------------------------------------------------------------------------------------------------- 643 Code2 = 1000*kfla + 100*kflb + << 644 if ( (j == 1) && (k==1)) Code2 = 1000*2 + << 645 527 646 hadron2 = G4ParticleTable::Get << 528 G4ExcitedString *G4VLongitudinalStringDecay::CPExcited(const G4ExcitedString & in) 647 hadron3 = G4ParticleTable::Get << 529 { >> 530 G4Parton *Left=new G4Parton(*in.GetLeftParton()); >> 531 G4Parton *Right=new G4Parton(*in.GetRightParton()); >> 532 return new G4ExcitedString(Left,Right,in.GetDirection()); >> 533 } 648 534 649 if ((hadron2 == nullptr) && (h << 535 G4double G4VLongitudinalStringDecay::FragmentationMass( >> 536 const G4FragmentingString * >> 537 const string, >> 538 Pcreate build, >> 539 pDefPair * pdefs) >> 540 { >> 541 >> 542 G4double mass; 650 543 651 if ((hadron2 != nullptr) && (h << 544 if ( build==0 ) build=&G4HadronBuilder::BuildLowSpin; 652 if (hadron2->GetPDGMass() > << 653 }; << 654 545 655 if ((hadron2 != nullptr) && (h << 546 G4ParticleDefinition *Hadron1, *Hadron2=0; 656 547 657 if ((hadron2 == nullptr) && (h << 548 if (!string->FourQuarkString() ) >> 549 { >> 550 // spin 0 meson or spin 1/2 barion will be built 658 551 659 minMassQDiQStr[i-1][j-1][k-1] << 552 Hadron1 = (hadronizer->*build)(string->GetLeftParton(), 660 } << 553 string->GetRightParton()); >> 554 mass= (Hadron1)->GetPDGMass(); >> 555 } else >> 556 { >> 557 //... string is qq--qqbar: Build two stable hadrons, >> 558 //... with extra uubar or ddbar quark pair >> 559 G4int iflc = (G4UniformRand() < 0.5)? 1 : 2; >> 560 if (string->GetLeftParton()->GetPDGEncoding() < 0) iflc = -iflc; >> 561 >> 562 //... theSpin = 4; spin 3/2 baryons will be built >> 563 Hadron1 = (hadronizer->*build)(string->GetLeftParton(),FindParticle(iflc)); >> 564 Hadron2 =(hadronizer->*build)(string->GetRightParton(),FindParticle(-iflc)); >> 565 mass = (Hadron1)->GetPDGMass() + (Hadron2)->GetPDGMass(); 661 } 566 } 662 } << 567 663 << 568 if ( pdefs != 0 ) 664 // ------ An estimated minimal string mass << 569 { // need to return hadrons as well.... 665 MinimalStringMass = 0.; << 570 pdefs->first = Hadron1; 666 MinimalStringMass2 = 0.; << 571 pdefs->second = Hadron2; 667 // q charges d u << 572 } 668 Qcharge[0] = -1; Qcharge[1] = 2; Qcharge[2 << 573 669 << 574 return mass; 670 // For treating of small string decays << 575 } 671 for (G4int i=0; i<5; i++) << 672 { for (G4int j=0; j<5; j++) << 673 { for (G4int k=0; k<7; k++) << 674 { << 675 Meson[i][j][k]=0; MesonWeight[i][j << 676 } << 677 } << 678 } << 679 //-------------------------- << 680 G4int StrangeQ = 0; << 681 G4int StrangeAQ = 0; << 682 for (G4int i=0; i<5; i++) << 683 { << 684 if( i >= 2 ) StrangeQ=1; << 685 for (G4int j=0; j<5; j++) << 686 { << 687 StrangeAQ = 0; << 688 if( j >= 2 ) StrangeAQ=1; << 689 Meson[i][j][0] = 100 * (std::ma << 690 MesonWeight[i][j][0] = ( pspin_meso << 691 Meson[i][j][1] = 100 * (std::ma << 692 MesonWeight[i][j][1] = (1.-pspin_meso << 693 } << 694 } << 695 << 696 //qqs << 697 //dd1 -> scalarMesonMix[0] * 111 + (1-scal << 698 //dd1 -> Pi0 << 699 << 700 Meson[0][0][0] = 111; MesonWeight[0][0][0] << 701 Meson[0][0][2] = 221; MesonWeight[0][0][3] << 702 Meson[0][0][3] = 331; MesonWeight[0][0][4] << 703 << 704 //dd3 -> (1-vectorMesonMix[1] * 113 + vect << 705 //dd3 -> rho_0 << 706 << 707 Meson[0][0][1] = 113; MesonWeight[0][0][1] << 708 Meson[0][0][4] = 223; MesonWeight[0][0][4] << 709 << 710 //uu1 -> scalarMesonMix[0] * 111 + (1-scal << 711 //uu1 -> Pi0 << 712 << 713 Meson[1][1][0] = 111; MesonWeight[1][1][0] << 714 Meson[1][1][2] = 221; MesonWeight[1][1][2] << 715 Meson[1][1][3] = 331; MesonWeight[1][1][3] << 716 << 717 //uu3 -> (1-vectorMesonMix[1]) * 113 + vec << 718 //uu3 -> rho_0 << 719 << 720 Meson[1][1][1] = 113; MesonWeight[1][1][1] << 721 Meson[1][1][4] = 223; MesonWeight[1][1][4] << 722 << 723 //ss1 -> << 724 //ss1 -> << 725 << 726 Meson[2][2][0] = 221; MesonWeight[2][2][0] << 727 Meson[2][2][2] = 331; MesonWeight[2][2][2] << 728 << 729 //ss3 -> << 730 //ss3 -> << 731 << 732 Meson[2][2][1] = 333; MesonWeight[2][2][1] << 733 << 734 //cc1 -> ProbEta_c /(1-pspin_meson) << 735 //cc3 -> (1-ProbEta_c)/( pspin_meson) << 736 << 737 //bb1 -> ProbEta_b /pspin_meson 551 << 738 //bb3 -> (1-ProbEta_b)/pspin_meson 553 << 739 << 740 if ( pspin_meson[2] != 0. ) { << 741 Meson[3][3][0] *= ( ProbEta_c)/( p << 742 Meson[3][3][1] *= (1.0-ProbEta_c)/(1.-p << 743 << 744 Meson[4][4][0] *= ( ProbEta_b)/( p << 745 Meson[4][4][1] *= (1.0-ProbEta_b)/(1.-p << 746 } << 747 << 748 //-------------------------- << 749 << 750 for (G4int i=0; i<5; i++) << 751 { for (G4int j=0; j<5; j++) << 752 { for (G4int k=0; k<5; k++) << 753 { for (G4int l=0; l<4; l++) << 754 { Baryon[i][j][k][l]=0; BaryonWei << 755 } << 756 } << 757 } << 758 576 759 kfla =0; kflb =0; << 577 //---------------------------------------------------------------------------------------------------------- 760 G4int kflc(0), kfld(0), << 761 for (G4int i=0; i<5; i++) << 762 { for (G4int j=0; j<5; j++) << 763 { for (G4int k=0; k<5; k++) << 764 { << 765 kfla = i+1; kflb = j+1; kflc = k+1; << 766 kfld = std::max(kfla,kflb); << 767 kfld = std::max(kfld,kflc); << 768 << 769 kflf = std::min(kfla,kflb); << 770 kflf = std::min(kflf,kflc); << 771 << 772 kfle = kfla + kflb + kflc - kfld - << 773 << 774 Baryon[i][j][k][0] = 1000 * k << 775 BaryonWeight[i][j][k][0] = ( pspi << 776 Baryon[i][j][k][1] = 1000 * k << 777 BaryonWeight[i][j][k][1] = (1.-pspi << 778 } << 779 } << 780 } << 781 578 782 // Delta- ddd - only 1114 << 579 G4bool G4VLongitudinalStringDecay::IsFragmentable(const G4FragmentingString * const string) 783 Baryon[0][0][0][0] = 1114; BaryonWeight << 580 { 784 Baryon[0][0][0][1] = 0; BaryonWeight << 581 return sqr(FragmentationMass(string)+MassCut) < 785 << 582 string->Mass2(); 786 // Delta++ uuu - only 2224 << 583 } 787 Baryon[1][1][1][0] = 2224; BaryonWeight << 788 Baryon[1][1][1][1] = 0; BaryonWeight << 789 << 790 // Omega- sss - only 3334 << 791 Baryon[2][2][2][0] = 3334; BaryonWeight << 792 Baryon[2][2][2][1] = 0; BaryonWeight << 793 << 794 // Omega_cc++ ccc - only 4444 << 795 Baryon[3][3][3][0] = 4444; BaryonWeight << 796 Baryon[3][3][3][1] = 0; BaryonWeight << 797 << 798 // Omega_bb- bbb - only 5554 << 799 Baryon[4][4][4][0] = 5554; BaryonWeight << 800 Baryon[4][4][4][1] = 0; BaryonWeight << 801 << 802 // Lambda/Sigma0 sud - 3122/3212 << 803 Baryon[0][1][2][0] = 3122; BaryonWeight << 804 Baryon[0][2][1][0] = 3122; BaryonWeight << 805 Baryon[1][0][2][0] = 3122; BaryonWeight << 806 Baryon[1][2][0][0] = 3122; BaryonWeight << 807 Baryon[2][0][1][0] = 3122; BaryonWeight << 808 Baryon[2][1][0][0] = 3122; BaryonWeight << 809 << 810 Baryon[0][1][2][2] = 3212; BaryonWeight << 811 Baryon[0][2][1][2] = 3212; BaryonWeight << 812 Baryon[1][0][2][2] = 3212; BaryonWeight << 813 Baryon[1][2][0][2] = 3212; BaryonWeight << 814 Baryon[2][0][1][2] = 3212; BaryonWeight << 815 Baryon[2][1][0][2] = 3212; BaryonWeight << 816 << 817 // Lambda_c+/Sigma_c+ cud - 4122/4212 << 818 Baryon[0][1][3][0] = 4122; BaryonWeight << 819 Baryon[0][3][1][0] = 4122; BaryonWeight << 820 Baryon[1][0][3][0] = 4122; BaryonWeight << 821 Baryon[1][3][0][0] = 4122; BaryonWeight << 822 Baryon[3][0][1][0] = 4122; BaryonWeight << 823 Baryon[3][1][0][0] = 4122; BaryonWeight << 824 << 825 Baryon[0][1][3][2] = 4212; BaryonWeight << 826 Baryon[0][3][1][2] = 4212; BaryonWeight << 827 Baryon[1][0][3][2] = 4212; BaryonWeight << 828 Baryon[1][3][0][2] = 4212; BaryonWeight << 829 Baryon[3][0][1][2] = 4212; BaryonWeight << 830 Baryon[3][1][0][2] = 4212; BaryonWeight << 831 << 832 // Xi_c+/Xi_c+' cus - 4232/4322 << 833 Baryon[1][2][3][0] = 4232; BaryonWeight << 834 Baryon[1][3][2][0] = 4232; BaryonWeight << 835 Baryon[2][1][3][0] = 4232; BaryonWeight << 836 Baryon[2][3][1][0] = 4232; BaryonWeight << 837 Baryon[3][1][2][0] = 4232; BaryonWeight << 838 Baryon[3][2][1][0] = 4232; BaryonWeight << 839 << 840 Baryon[1][2][3][2] = 4322; BaryonWeight << 841 Baryon[1][3][2][2] = 4322; BaryonWeight << 842 Baryon[2][1][3][2] = 4322; BaryonWeight << 843 Baryon[2][3][1][2] = 4322; BaryonWeight << 844 Baryon[3][1][2][2] = 4322; BaryonWeight << 845 Baryon[3][2][1][2] = 4322; BaryonWeight << 846 << 847 // Xi_c0/Xi_c0' cus - 4132/4312 << 848 Baryon[0][2][3][0] = 4132; BaryonWeight << 849 Baryon[0][3][2][0] = 4132; BaryonWeight << 850 Baryon[2][0][3][0] = 4132; BaryonWeight << 851 Baryon[2][3][0][0] = 4132; BaryonWeight << 852 Baryon[3][0][2][0] = 4132; BaryonWeight << 853 Baryon[3][2][0][0] = 4132; BaryonWeight << 854 << 855 Baryon[0][2][3][2] = 4312; BaryonWeight << 856 Baryon[0][3][2][2] = 4312; BaryonWeight << 857 Baryon[2][0][3][2] = 4312; BaryonWeight << 858 Baryon[2][3][0][2] = 4312; BaryonWeight << 859 Baryon[3][0][2][2] = 4312; BaryonWeight << 860 Baryon[3][2][0][2] = 4312; BaryonWeight << 861 << 862 // Lambda_b0/Sigma_b0 bud - 5122/5212 << 863 Baryon[0][1][4][0] = 5122; BaryonWeight << 864 Baryon[0][4][1][0] = 5122; BaryonWeight << 865 Baryon[1][0][4][0] = 5122; BaryonWeight << 866 Baryon[1][4][0][0] = 5122; BaryonWeight << 867 Baryon[4][0][1][0] = 5122; BaryonWeight << 868 Baryon[4][1][0][0] = 5122; BaryonWeight << 869 << 870 Baryon[0][1][4][2] = 5212; BaryonWeight << 871 Baryon[0][4][1][2] = 5212; BaryonWeight << 872 Baryon[1][0][4][2] = 5212; BaryonWeight << 873 Baryon[1][4][0][2] = 5212; BaryonWeight << 874 Baryon[4][0][1][2] = 5212; BaryonWeight << 875 Baryon[4][1][0][2] = 5212; BaryonWeight << 876 << 877 // Xi_b0/Xi_b0' bus - 5232/5322 << 878 Baryon[1][2][4][0] = 5232; BaryonWeight << 879 Baryon[1][4][2][0] = 5232; BaryonWeight << 880 Baryon[2][1][4][0] = 5232; BaryonWeight << 881 Baryon[2][4][1][0] = 5232; BaryonWeight << 882 Baryon[4][1][2][0] = 5232; BaryonWeight << 883 Baryon[4][2][1][0] = 5232; BaryonWeight << 884 << 885 Baryon[1][2][4][2] = 5322; BaryonWeight << 886 Baryon[1][4][2][2] = 5322; BaryonWeight << 887 Baryon[2][1][4][2] = 5322; BaryonWeight << 888 Baryon[2][4][1][2] = 5322; BaryonWeight << 889 Baryon[4][1][2][2] = 5322; BaryonWeight << 890 Baryon[4][2][1][2] = 5322; BaryonWeight << 891 << 892 // Xi_b-/Xi_b-' bus - 5132/5312 << 893 Baryon[0][2][4][0] = 5132; BaryonWeight << 894 Baryon[0][4][2][0] = 5132; BaryonWeight << 895 Baryon[2][0][4][0] = 5132; BaryonWeight << 896 Baryon[2][4][0][0] = 5132; BaryonWeight << 897 Baryon[4][0][2][0] = 5132; BaryonWeight << 898 Baryon[4][2][0][0] = 5132; BaryonWeight << 899 << 900 Baryon[0][2][4][2] = 5312; BaryonWeight << 901 Baryon[0][4][2][2] = 5312; BaryonWeight << 902 Baryon[2][0][4][2] = 5312; BaryonWeight << 903 Baryon[2][4][0][2] = 5312; BaryonWeight << 904 Baryon[4][0][2][2] = 5312; BaryonWeight << 905 Baryon[4][2][0][2] = 5312; BaryonWeight << 906 << 907 for (G4int i=0; i<5; i++) << 908 { for (G4int j=0; j<5; j++) << 909 { for (G4int k=0; k<5; k++) << 910 { for (G4int l=0; l<4; l++) << 911 { << 912 G4ParticleDefinition * Te << 913 G4ParticleTable::GetPar << 914 /* << 915 G4cout<<i<<" "<<j<<" "<<k << 916 if (TestHadron != nullptr << 917 if ((TestHadron == nullpt << 918 if ((TestHadron == nullpt << 919 G4cout<<G4endl; << 920 */ << 921 if ((TestHadron == nullpt << 922 } << 923 } << 924 } << 925 } << 926 584 927 // --------- Probabilities of q-qbar pair << 585 //---------------------------------------------------------------------------------------------------------- 928 G4double ProbUUbar = 0.33; << 929 Prob_QQbar[0]=ProbUUbar; // Probab << 930 Prob_QQbar[1]=ProbUUbar; // Probab << 931 Prob_QQbar[2]=1.0-2.*ProbUUbar; // Probab << 932 Prob_QQbar[3]=0.0; // Probab << 933 Prob_QQbar[4]=0.0; // Probab << 934 << 935 for ( G4int i=0 ; i<350 ; i++ ) { // Must << 936 FS_LeftHadron[i] = 0; << 937 FS_RightHadron[i] = 0; << 938 FS_Weight[i] = 0.0; << 939 } << 940 586 941 NumberOf_FS = 0; << 587 G4bool G4VLongitudinalStringDecay::StopFragmenting(const G4FragmentingString * const string) >> 588 { >> 589 return >> 590 sqr(FragmentationMass(string,&G4HadronBuilder::BuildHighSpin)+MassCut) > >> 591 string->Get4Momentum().mag2(); 942 } 592 } 943 593 944 // ------------------------------------------- << 594 //---------------------------------------------------------------------------------------------------------- 945 595 946 void G4VLongitudinalStringDecay::SetMinimalStr << 596 G4KineticTrackVector* G4VLongitudinalStringDecay::LightFragmentationTest(const >> 597 G4ExcitedString * const string) 947 { 598 { 948 //MaxMass = -350.0*GeV; << 599 // Check string decay threshold 949 G4double EstimatedMass=MaxMass; << 600 >> 601 G4KineticTrackVector * result=0; // return 0 when string exceeds the mass cut >> 602 >> 603 pDefPair hadrons((G4ParticleDefinition *)0,(G4ParticleDefinition *)0); >> 604 G4FragmentingString aString(*string); >> 605 if ( sqr(FragmentationMass(&aString,0,&hadrons)+MassCut) < aString.Mass2()) { >> 606 return 0; >> 607 } >> 608 >> 609 result=new G4KineticTrackVector; >> 610 >> 611 if ( hadrons.second ==0 ) >> 612 { >> 613 // Substitute string by light hadron, Note that Energy is not conserved here! >> 614 >> 615 G4ThreeVector Mom3 = string->Get4Momentum().vect(); >> 616 G4LorentzVector Mom(Mom3, >> 617 std::sqrt(Mom3.mag2() + sqr(hadrons.first->GetPDGMass()))); >> 618 result->push_back(new G4KineticTrack(hadrons.first, 0, string->GetPosition(), Mom)); >> 619 } else >> 620 { >> 621 //... string was qq--qqbar type: Build two stable hadrons, >> 622 G4LorentzVector Mom1, Mom2; >> 623 Sample4Momentum(&Mom1, hadrons.first->GetPDGMass(), >> 624 &Mom2,hadrons.second->GetPDGMass(), >> 625 string->Get4Momentum().mag()); >> 626 result->push_back(new G4KineticTrack(hadrons.first, 0, string->GetPosition(), Mom1)); >> 627 result->push_back(new G4KineticTrack(hadrons.second, 0, string->GetPosition(), Mom2)); >> 628 G4ThreeVector Velocity = string->Get4Momentum().boostVector(); >> 629 result->Boost(Velocity); >> 630 } 950 631 951 G4ParticleDefinition* LeftParton = st << 632 return result; 952 G4ParticleDefinition* RightParton = st << 953 if( LeftParton->GetParticleSubType() = << 954 if( LeftParton->GetPDGEncoding() * R << 955 // Not allowed combination of the << 956 throw G4HadronicException(__FILE__ << 957 "G4VLongitudinalStringDecay::Set << 958 } << 959 } << 960 if( LeftParton->GetParticleSubType() ! << 961 if( LeftParton->GetPDGEncoding() * R << 962 // Not allowed combination of the << 963 throw G4HadronicException(__FILE__ << 964 "G4VLongitudinalStringDecay::Set << 965 } << 966 } << 967 633 968 G4int Qleft =std::abs(string->GetLeftP << 634 } 969 G4int Qright=std::abs(string->GetRight << 970 635 971 if ((Qleft < 6) && (Qright < 6)) { / << 636 //---------------------------------------------------------------------------------------------------------- 972 EstimatedMass=minMassQQbarStr[Qleft- << 973 MinimalStringMass=EstimatedMass; << 974 SetMinimalStringMass2(EstimatedMass) << 975 return; << 976 } << 977 637 978 if ((Qleft < 6) && (Qright > 1000)) { << 638 G4ParticleDefinition* G4VLongitudinalStringDecay::FindParticle(G4int Encoding) 979 G4int q1=Qright/1000; << 639 { 980 G4int q2=(Qright/100)%10; << 640 G4ParticleDefinition* ptr = G4ParticleTable::GetParticleTable()->FindParticle(Encoding); 981 EstimatedMass=minMassQDiQStr[Qleft-1 << 641 if (ptr == NULL) 982 MinimalStringMass=EstimatedMass; << 642 { 983 SetMinimalStringMass2(EstimatedMass) << 643 G4cout << "Particle with encoding "<<Encoding<<" does not exist!!!"<<G4endl; 984 return; << 644 throw G4HadronicException(__FILE__, __LINE__, "Check your particle table"); 985 } << 645 } >> 646 return ptr; >> 647 } 986 648 987 if ((Qleft > 1000) && (Qright < 6)) { << 649 //---------------------------------------------------------------------------------------------------------- 988 G4int q1=Qleft/1000; << 989 G4int q2=(Qleft/100)%10; << 990 EstimatedMass=minMassQDiQStr[Qright- << 991 MinimalStringMass=EstimatedMass; << 992 SetMinimalStringMass2(EstimatedMass) << 993 return; << 994 } << 995 650 996 // DiQuark - Anti DiQuark string ----- << 651 void G4VLongitudinalStringDecay::Sample4Momentum(G4LorentzVector* Mom, G4double Mass, G4LorentzVector* AntiMom, G4double AntiMass, G4double InitialMass) >> 652 { >> 653 G4double r_val = sqr(InitialMass*InitialMass - Mass*Mass - AntiMass*AntiMass) - sqr(2.*Mass*AntiMass); >> 654 G4double Pabs = (r_val > 0.)? std::sqrt(r_val)/(2.*InitialMass) : 0; 997 655 998 G4double StringM=string->Get4Momentum().mag( << 656 //... sample unit vector >> 657 G4double pz = 1. - 2.*G4UniformRand(); >> 658 G4double st = std::sqrt(1. - pz * pz)*Pabs; >> 659 G4double phi = 2.*pi*G4UniformRand(); >> 660 G4double px = st*std::cos(phi); >> 661 G4double py = st*std::sin(phi); >> 662 pz *= Pabs; >> 663 >> 664 Mom->setPx(px); Mom->setPy(py); Mom->setPz(pz); >> 665 Mom->setE(std::sqrt(Pabs*Pabs + Mass*Mass)); 999 666 1000 #ifdef debug_LUNDfragmentation << 667 AntiMom->setPx(-px); AntiMom->setPy(-py); AntiMom->setPz(-pz); 1001 // G4cout<<"MinStringMass// Input Str << 668 AntiMom->setE (std::sqrt(Pabs*Pabs + AntiMass*AntiMass)); 1002 #endif << 669 } 1003 670 1004 G4int q1= Qleft/1000 ; << 671 //---------------------------------------------------------------------------------------------------------- 1005 G4int q2=(Qleft/100)%10 ; << 1006 672 1007 G4int q3= Qright/1000 ; << 673 void G4VLongitudinalStringDecay::SetSigmaTransverseMomentum(G4double aValue) 1008 G4int q4=(Qright/100)%10; << 674 { >> 675 if ( PastInitPhase ) { >> 676 throw G4HadronicException(__FILE__, __LINE__, "4VLongitudinalStringDecay::SetSigmaTransverseMomentum after FragmentString() not allowed"); >> 677 } else { >> 678 SigmaQT = aValue; >> 679 } >> 680 } 1009 681 1010 // -------------- 2 baryon production << 682 //---------------------------------------------------------------------------------------------------------- 1011 683 1012 G4double EstimatedMass1 = minMassQDiQ << 684 void G4VLongitudinalStringDecay::SetStrangenessSuppression(G4double aValue) 1013 G4double EstimatedMass2 = minMassQDiQ << 685 { 1014 // Mass is negative if there is no co << 686 if ( PastInitPhase ) { >> 687 throw G4HadronicException(__FILE__, __LINE__, "4VLongitudinalStringDecay::SetStrangenessSuppression after FragmentString() not allowed"); >> 688 } else { >> 689 StrangeSuppress = aValue; >> 690 } >> 691 } 1015 692 1016 if ( (EstimatedMass1 > 0.) && (Estima << 693 //---------------------------------------------------------------------------------------------------------- 1017 EstimatedMass = EstimatedMass1 + E << 1018 if ( StringM > EstimatedMass ) { << 1019 MinimalStringMass=EstimatedMass << 1020 SetMinimalStringMass2(Estimated << 1021 return; << 1022 } << 1023 } << 1024 694 1025 if ( (EstimatedMass1 < 0.) && (Estima << 695 void G4VLongitudinalStringDecay::SetDiquarkSuppression(G4double aValue) 1026 EstimatedMass = MaxMass; << 696 { 1027 MinimalStringMass=EstimatedMass; << 697 if ( PastInitPhase ) { 1028 SetMinimalStringMass2(EstimatedMas << 698 throw G4HadronicException(__FILE__, __LINE__, "4VLongitudinalStringDecay::SetDiquarkSuppression after FragmentString() not allowed"); 1029 return; << 699 } else { 1030 } << 700 DiquarkSuppress = aValue; >> 701 } >> 702 } 1031 703 1032 if ( (EstimatedMass1 > 0.) && (Estima << 704 void G4VLongitudinalStringDecay::SetDiquarkBreakProbability(G4double aValue) 1033 EstimatedMass = EstimatedMass1; << 705 { 1034 MinimalStringMass=EstimatedMass; << 706 if ( PastInitPhase ) { 1035 SetMinimalStringMass2(EstimatedMas << 707 throw G4HadronicException(__FILE__, __LINE__, "4VLongitudinalStringDecay::SetDiquarkBreakProbability after FragmentString() not allowed"); 1036 return; << 708 } else { 1037 } << 709 DiquarkBreakProb = aValue; >> 710 } >> 711 } 1038 712 1039 // if ( EstimatedMass >= StringM << 713 //---------------------------------------------------------------------------------------------------------- 1040 // ------------- Re-orangement ------ << 1041 EstimatedMass=std::min(minMassQQbarSt << 1042 minMassQQbarSt << 1043 714 1044 // In principle, re-arrangement and 2 << 715 void G4VLongitudinalStringDecay::SetVectorMesonProbability(G4double aValue) 1045 // More physics consideration is need << 716 { >> 717 if ( PastInitPhase ) { >> 718 throw G4HadronicException(__FILE__, __LINE__, "4VLongitudinalStringDecay::SetVectorMesonProbability after FragmentString() not allowed"); >> 719 } else { >> 720 pspin_meson = aValue; >> 721 delete hadronizer; >> 722 hadronizer = new G4HadronBuilder(pspin_meson,pspin_barion, >> 723 scalarMesonMix,vectorMesonMix); >> 724 } >> 725 } 1046 726 1047 MinimalStringMass=EstimatedMass; << 727 //---------------------------------------------------------------------------------------------------------- 1048 SetMinimalStringMass2(EstimatedMass); << 1049 728 1050 return; << 729 void G4VLongitudinalStringDecay::SetSpinThreeHalfBarionProbability(G4double aValue) >> 730 { >> 731 if ( PastInitPhase ) { >> 732 throw G4HadronicException(__FILE__, __LINE__, "4VLongitudinalStringDecay::SetSpinThreeHalfBarionProbability after FragmentString() not allowed"); >> 733 } else { >> 734 pspin_barion = aValue; >> 735 delete hadronizer; >> 736 hadronizer = new G4HadronBuilder(pspin_meson,pspin_barion, >> 737 scalarMesonMix,vectorMesonMix); >> 738 } 1051 } 739 } 1052 740 1053 //------------------------------------------- << 741 //---------------------------------------------------------------------------------------------------------- 1054 742 1055 void G4VLongitudinalStringDecay::SetMinimalSt << 743 void G4VLongitudinalStringDecay::SetScalarMesonMixings(std::vector<G4double> aVector) 1056 { 744 { 1057 MinimalStringMass2=aValue * aValue; << 745 if ( PastInitPhase ) { >> 746 throw G4HadronicException(__FILE__, __LINE__, "4VLongitudinalStringDecay::SetScalarMesonMixings after FragmentString() not allowed"); >> 747 } else { >> 748 if ( aVector.size() < 6 ) >> 749 throw G4HadronicException(__FILE__, __LINE__, "4VLongitudinalStringDecay::SetScalarMesonMixings( argument Vector too small"); >> 750 scalarMesonMix[0] = aVector[0]; >> 751 scalarMesonMix[1] = aVector[1]; >> 752 scalarMesonMix[2] = aVector[2]; >> 753 scalarMesonMix[3] = aVector[3]; >> 754 scalarMesonMix[4] = aVector[4]; >> 755 scalarMesonMix[5] = aVector[5]; >> 756 delete hadronizer; >> 757 hadronizer = new G4HadronBuilder(pspin_meson,pspin_barion, >> 758 scalarMesonMix,vectorMesonMix); >> 759 } 1058 } 760 } 1059 761 >> 762 //---------------------------------------------------------------------------------------------------------- >> 763 >> 764 void G4VLongitudinalStringDecay::SetVectorMesonMixings(std::vector<G4double> aVector) >> 765 { >> 766 if ( PastInitPhase ) { >> 767 throw G4HadronicException(__FILE__, __LINE__, "4VLongitudinalStringDecay::SetVectorMesonMixings after FragmentString() not allowed"); >> 768 } else { >> 769 if ( aVector.size() < 6 ) >> 770 throw G4HadronicException(__FILE__, __LINE__, "4VLongitudinalStringDecay::SetVectorMesonMixings( argument Vector too small"); >> 771 vectorMesonMix[0] = aVector[0]; >> 772 vectorMesonMix[1] = aVector[1]; >> 773 vectorMesonMix[2] = aVector[2]; >> 774 vectorMesonMix[3] = aVector[3]; >> 775 vectorMesonMix[4] = aVector[4]; >> 776 vectorMesonMix[5] = aVector[5]; >> 777 delete hadronizer; >> 778 hadronizer = new G4HadronBuilder(pspin_meson,pspin_barion, >> 779 scalarMesonMix,vectorMesonMix); >> 780 >> 781 } >> 782 } >> 783 //******************************************************************************************************* 1060 784