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