Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // 26 // 26 // >> 27 // $Id: G4VLongitudinalStringDecay.cc,v 1.13.2.1 2009/08/11 13:05:47 gcosmo Exp $ >> 28 // GEANT4 tag $Name: geant4-09-02-patch-03 $ 27 // 29 // 28 // ------------------------------------------- 30 // ----------------------------------------------------------------------------- 29 // GEANT 4 class implementation file 31 // GEANT 4 class implementation file 30 // 32 // 31 // History: first implementation, Maxim K 33 // History: first implementation, Maxim Komogorov, 1-Jul-1998 32 // redesign Gunter Folger, Augu 34 // redesign Gunter Folger, August/September 2001 33 // ------------------------------------------- 35 // ----------------------------------------------------------------------------- 34 #include "G4VLongitudinalStringDecay.hh" << 35 #include "G4PhysicalConstants.hh" << 36 #include "G4SystemOfUnits.hh" << 37 #include "G4ios.hh" 36 #include "G4ios.hh" 38 #include "Randomize.hh" 37 #include "Randomize.hh" >> 38 #include "G4VLongitudinalStringDecay.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 "G4ShortLivedTable.hh" 47 #include "G4PhaseSpaceDecayChannel.hh" 48 #include "G4PhaseSpaceDecayChannel.hh" 48 #include "G4VDecayChannel.hh" 49 #include "G4VDecayChannel.hh" 49 #include "G4DecayTable.hh" 50 #include "G4DecayTable.hh" 50 51 51 #include "G4DiQuarks.hh" 52 #include "G4DiQuarks.hh" 52 #include "G4Quarks.hh" 53 #include "G4Quarks.hh" 53 #include "G4Gluons.hh" 54 #include "G4Gluons.hh" 54 55 55 #include "G4Exp.hh" << 56 #include "G4Log.hh" << 57 << 58 #include "G4HadronicException.hh" << 59 << 60 //------------------------debug switches 56 //------------------------debug switches 61 //#define debug_VStringDecay << 57 //#define DEBUG_LightFragmentationTest 1 62 //#define debug_heavyHadrons << 63 58 64 //******************************************** << 59 >> 60 //******************************************************************************** 65 // Constructors 61 // Constructors 66 62 67 G4VLongitudinalStringDecay::G4VLongitudinalStr << 63 G4VLongitudinalStringDecay::G4VLongitudinalStringDecay() 68 : G4HadronicInteraction(name), ProbCCbar(0.0 << 69 { 64 { 70 MassCut = 210.0*MeV; // Mpi + Delta << 65 MassCut = 0.35*GeV; >> 66 ClusterMass = 0.15*GeV; 71 67 72 StringLoopInterrupt = 1000; << 68 SmoothParam = 0.9; 73 ClusterLoopInterrupt = 500; << 69 StringLoopInterrupt = 1000; >> 70 ClusterLoopInterrupt = 500; 74 71 75 // Changable Parameters below. << 72 // Changable Parameters below. 76 SigmaQT = 0.5 * GeV; << 77 73 78 StrangeSuppress = 0.44; // =0.27/2.27 s << 74 SigmaQT = 0.5 * GeV; 79 DiquarkSuppress = 0.07; // Probability << 80 DiquarkBreakProb = 0.1; // Probability << 81 75 82 //... pspin_meson is probability to create << 76 StrangeSuppress = 0.44; // 27 % strange quarks produced, ie. u:d:s=1:1:0.27 83 pspin_meson.resize(3); << 77 DiquarkSuppress = 0.07; 84 pspin_meson[0] = 0.5; // u or d + anti-u o << 78 DiquarkBreakProb = 0.1; 85 pspin_meson[1] = 0.4; // one of the quark << 86 pspin_meson[2] = 0.3; // both of the quark << 87 79 88 //... pspin_barion is probability to create << 80 //... pspin_meson is probability to create vector meson >> 81 pspin_meson = 0.5; >> 82 >> 83 //... pspin_barion is probability to create 3/2 barion 89 pspin_barion = 0.5; 84 pspin_barion = 0.5; 90 85 91 //... vectorMesonMix[] is quark mixing para 86 //... vectorMesonMix[] is quark mixing parameters for vector mesons (Variable spin = 3) 92 vectorMesonMix.resize(6); 87 vectorMesonMix.resize(6); 93 vectorMesonMix[0] = 0.0; << 88 vectorMesonMix[0] = 0.5; 94 vectorMesonMix[1] = 0.5; << 89 vectorMesonMix[1] = 0.0; 95 vectorMesonMix[2] = 0.0; << 90 vectorMesonMix[2] = 0.5; 96 vectorMesonMix[3] = 0.5; << 91 vectorMesonMix[3] = 0.0; 97 vectorMesonMix[4] = 1.0; 92 vectorMesonMix[4] = 1.0; 98 vectorMesonMix[5] = 1.0; << 93 vectorMesonMix[5] = 1.0; 99 94 100 //... scalarMesonMix[] is quark mixing para 95 //... scalarMesonMix[] is quark mixing parameters for scalar mesons (Variable spin=1) 101 scalarMesonMix.resize(6); 96 scalarMesonMix.resize(6); 102 scalarMesonMix[0] = 0.5; << 97 scalarMesonMix[0] = 0.5; 103 scalarMesonMix[1] = 0.25; << 98 scalarMesonMix[1] = 0.25; 104 scalarMesonMix[2] = 0.5; << 99 scalarMesonMix[2] = 0.5; 105 scalarMesonMix[3] = 0.25; << 100 scalarMesonMix[3] = 0.25; 106 scalarMesonMix[4] = 1.0; << 101 scalarMesonMix[4] = 1.0; 107 scalarMesonMix[5] = 0.5; << 102 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 103 114 // Parameters may be changed until the firs << 104 // Parameters may be changed until the first fragmentation starts 115 PastInitPhase=false; 105 PastInitPhase=false; 116 hadronizer = new G4HadronBuilder( pspin_mes << 106 hadronizer = new G4HadronBuilder(pspin_meson,pspin_barion, 117 ProbEta_c << 107 scalarMesonMix,vectorMesonMix); 118 << 108 Kappa = 1.0 * GeV/fermi; 119 MaxMass=-350.0*GeV; // If there will be a << 120 109 121 SetMinMasses(); // Re-calculation of minim << 122 110 123 Kappa = 1.0 * GeV/fermi; << 124 DecayQuark = NewQuark = 0; << 125 } 111 } >> 112 126 113 127 G4VLongitudinalStringDecay::~G4VLongitudinalSt 114 G4VLongitudinalStringDecay::~G4VLongitudinalStringDecay() 128 { << 115 { 129 delete hadronizer; 116 delete hadronizer; 130 } << 117 } 131 << 132 G4HadFinalState* << 133 G4VLongitudinalStringDecay::ApplyYourself(cons << 134 { << 135 return nullptr; << 136 } << 137 118 138 //============================================ 119 //============================================================================= 139 120 >> 121 // Operators >> 122 >> 123 //const & G4VLongitudinalStringDecay::operator=(const G4VLongitudinalStringDecay &) >> 124 // { >> 125 // } >> 126 >> 127 //----------------------------------------------------------------------------- >> 128 >> 129 int G4VLongitudinalStringDecay::operator==(const G4VLongitudinalStringDecay &) const >> 130 { >> 131 throw G4HadronicException(__FILE__, __LINE__, "G4VLongitudinalStringDecay::operator== forbidden"); >> 132 return false; >> 133 } >> 134 >> 135 //------------------------------------------------------------------------------------- >> 136 >> 137 int G4VLongitudinalStringDecay::operator!=(const G4VLongitudinalStringDecay &) const >> 138 { >> 139 throw G4HadronicException(__FILE__, __LINE__, "G4VLongitudinalStringDecay::operator!= forbidden"); >> 140 return true; >> 141 } >> 142 >> 143 //*********************************************************************************** >> 144 140 // For changing Mass Cut used for selection of 145 // For changing Mass Cut used for selection of very small mass strings 141 void G4VLongitudinalStringDecay::SetMassCut(G4 << 146 void G4VLongitudinalStringDecay::SetMassCut(G4double aValue){MassCut=aValue;} 142 G4double G4VLongitudinalStringDecay::GetMassCu << 143 147 144 //-------------------------------------------- 148 //----------------------------------------------------------------------------- 145 149 146 // For handling a string with very low mass 150 // For handling a string with very low mass 147 151 148 G4KineticTrackVector* G4VLongitudinalStringDec << 152 G4KineticTrackVector* G4VLongitudinalStringDecay::LightFragmentationTest(const >> 153 G4ExcitedString * const string) 149 { 154 { 150 G4KineticTrackVector* result = nullptr << 155 // Check string decay threshold 151 pDefPair hadrons( nullptr, nullptr ); << 156 152 G4FragmentingString aString( *string ) << 157 G4KineticTrackVector * result=0; // return 0 when string exceeds the mass cut 153 << 158 154 #ifdef debug_VStringDecay << 159 pDefPair hadrons((G4ParticleDefinition *)0,(G4ParticleDefinition *)0); 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 160 192 result->push_back( new G4Kineti << 161 G4FragmentingString aString(*string); 193 result->push_back( new G4Kineti << 162 >> 163 if ( sqr(FragmentationMass(&aString,0,&hadrons)+MassCut) < aString.Mass2()) { >> 164 return 0; >> 165 } >> 166 >> 167 // The string mass is very low --------------------------- >> 168 >> 169 result=new G4KineticTrackVector; >> 170 >> 171 if ( hadrons.second ==0 ) >> 172 { >> 173 // Substitute string by light hadron, Note that Energy is not conserved here! >> 174 >> 175 /* >> 176 #ifdef DEBUG_LightFragmentationTest >> 177 G4cout << "VlongSF Warning replacing string by single hadron " <<G4endl; >> 178 G4cout << hadrons.first->GetParticleName() >> 179 << "string .. " << string->Get4Momentum() << " " >> 180 << string->Get4Momentum().m() << G4endl; >> 181 #endif >> 182 G4cout << "VlongSF Warning replacing string by single hadron " <<G4endl; >> 183 G4cout << hadrons.first->GetParticleName() << " string .. " <<G4endl; >> 184 G4cout << string->Get4Momentum() << " " << string->Get4Momentum().m() << G4endl; >> 185 */ >> 186 G4ThreeVector Mom3 = string->Get4Momentum().vect(); >> 187 G4LorentzVector Mom(Mom3, >> 188 std::sqrt(Mom3.mag2() + >> 189 sqr(hadrons.first->GetPDGMass()))); >> 190 result->push_back(new G4KineticTrack(hadrons.first, 0, >> 191 string->GetPosition(), >> 192 Mom)); >> 193 } else >> 194 { >> 195 //... string was qq--qqbar type: Build two stable hadrons, >> 196 >> 197 #ifdef DEBUG_LightFragmentationTest >> 198 G4cout << "VlongSF Warning replacing qq-qqbar string by TWO hadrons " >> 199 << hadrons.first->GetParticleName() << " / " >> 200 << hadrons.second->GetParticleName() >> 201 << "string .. " << string->Get4Momentum() << " " >> 202 << string->Get4Momentum().m() << G4endl; >> 203 #endif >> 204 >> 205 // Uzhi Formation time in the case??? >> 206 >> 207 G4LorentzVector Mom1, Mom2; >> 208 Sample4Momentum(&Mom1, hadrons.first->GetPDGMass(), >> 209 &Mom2,hadrons.second->GetPDGMass(), >> 210 string->Get4Momentum().mag()); >> 211 >> 212 result->push_back(new G4KineticTrack(hadrons.first, 0, >> 213 string->GetPosition(), >> 214 Mom1)); >> 215 result->push_back(new G4KineticTrack(hadrons.second, 0, >> 216 string->GetPosition(), >> 217 Mom2)); 194 218 195 G4ThreeVector Velocity = string 219 G4ThreeVector Velocity = string->Get4Momentum().boostVector(); 196 result->Boost(Velocity); 220 result->Boost(Velocity); 197 } << 221 } 198 } << 222 199 return result; << 223 return result; >> 224 200 } 225 } 201 226 202 //-------------------------------------------- 227 //---------------------------------------------------------------------------------------- 203 228 204 G4double G4VLongitudinalStringDecay::PossibleH << 229 G4double G4VLongitudinalStringDecay::FragmentationMass( 205 << 230 const G4FragmentingString * const string, >> 231 Pcreate build, pDefPair * pdefs ) 206 { 232 { 207 G4double mass = 0.0; << 233 >> 234 G4double mass; >> 235 static G4bool NeedInit(true); >> 236 static std::vector<double> nomix; >> 237 static G4HadronBuilder * minMassHadronizer; >> 238 if ( NeedInit ) >> 239 { >> 240 NeedInit = false; >> 241 nomix.resize(6); >> 242 for ( G4int i=0; i<6 ; i++ ) nomix[i]=0; >> 243 >> 244 // minMassHadronizer=new G4HadronBuilder(pspin_meson,pspin_barion,nomix,nomix); >> 245 minMassHadronizer=hadronizer; >> 246 } 208 247 209 if ( build==0 ) build=&G4HadronBuilder::Buil 248 if ( build==0 ) build=&G4HadronBuilder::BuildLowSpin; 210 249 211 G4ParticleDefinition* Hadron1 = nullpt << 250 G4ParticleDefinition *Hadron1, *Hadron2=0; 212 G4ParticleDefinition* Hadron2 = nullptr; << 213 251 214 if (!string->IsAFourQuarkString() ) << 252 if (!string->FourQuarkString() ) 215 { 253 { 216 // spin 0 meson or spin 1/2 barion 254 // spin 0 meson or spin 1/2 barion will be built 217 255 218 Hadron1 = (hadronizer->*build)(stri << 256 Hadron1 = (minMassHadronizer->*build)(string->GetLeftParton(), 219 #ifdef debug_VStringDecay << 257 string->GetRightParton()); 220 G4cout<<"VlongSD PossibleHadronMass"<<G4e << 258 mass= (Hadron1)->GetPDGMass(); 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 259 } else 231 { 260 { 232 //... string is qq--qqbar: Build tw 261 //... string is qq--qqbar: Build two stable hadrons, 233 << 262 //... with extra uubar or ddbar quark pair 234 #ifdef debug_VStringDecay << 263 G4int iflc = (G4UniformRand() < 0.5)? 1 : 2; 235 G4cout<<"VlongSD PossibleHadronMass << 264 if (string->GetLeftParton()->GetPDGEncoding() < 0) iflc = -iflc; 236 G4cout<<"VlongSD string is qq--qqba << 265 237 #endif << 266 //... theSpin = 4; spin 3/2 baryons will be built 238 << 267 Hadron1 = (minMassHadronizer->*build)(string->GetLeftParton(), 239 G4double StringMass = string->Mas << 268 FindParticle(iflc) ); 240 G4int cClusterInterrupt = 0; << 269 Hadron2 = (minMassHadronizer->*build)(string->GetRightParton(), 241 do << 270 FindParticle(-iflc) ); 242 { << 271 mass = (Hadron1)->GetPDGMass() + (Hadron2)->GetPDGMass(); 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 << 264 mass = (Hadron1)->GetPDGMass() + (Hadron2 << 265 } 272 } 266 << 267 #ifdef debug_VStringDecay << 268 G4cout<<"VlongSD *Hadrons 1 and 2, pro << 269 #endif << 270 273 271 if ( pdefs != 0 ) 274 if ( pdefs != 0 ) 272 { // need to return hadrons as well.... 275 { // need to return hadrons as well.... 273 pdefs->first = Hadron1; 276 pdefs->first = Hadron1; 274 pdefs->second = Hadron2; 277 pdefs->second = Hadron2; 275 } 278 } 276 279 277 return mass; 280 return mass; 278 } 281 } 279 282 280 //-------------------------------------------- 283 //---------------------------------------------------------------------------- 281 284 282 G4ParticleDefinition* G4VLongitudinalStringDec 285 G4ParticleDefinition* G4VLongitudinalStringDecay::FindParticle(G4int Encoding) 283 { << 286 { 284 /* << 287 G4ParticleDefinition* ptr = G4ParticleTable::GetParticleTable()->FindParticle(Encoding); 285 G4cout<<Encoding<<" G4VLongitudinalStringDec << 288 if (ptr == NULL) 286 for (G4int i=4; i<6;i++){ << 289 { 287 for (G4int j=1;j<6;j++){ << 290 G4cout << "Particle with encoding "<<Encoding<<" does not exist!!!"<<G4endl; 288 G4cout<<i<<" "<<j<<" "; << 291 throw G4HadronicException(__FILE__, __LINE__, "Check your particle table"); 289 G4int Code = 1000 * i + 100 * j +1; << 292 } 290 G4ParticleDefinition* ptr1 = G4ParticleT << 293 return ptr; 291 Code +=2; << 294 } 292 G4ParticleDefinition* ptr2 = G4ParticleT << 293 G4cout<<"Code "<<Code - 2<<" ptr "<<ptr1 << 294 } << 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 295 309 return ptr; << 296 //----------------------------------------------------------------------------- 310 } << 297 // virtual void Sample4Momentum(G4LorentzVector* Mom, G4double Mass, >> 298 // G4LorentzVector* AntiMom, G4double AntiMass, >> 299 // G4double InitialMass)=0; >> 300 //----------------------------------------------------------------------------- 311 301 312 //******************************************** 302 //********************************************************************************* 313 // For decision on continue or stop string f 303 // For decision on continue or stop string fragmentation 314 // virtual G4bool StopFragmenting(const G4Fr 304 // virtual G4bool StopFragmenting(const G4FragmentingString * const string)=0; 315 // virtual G4bool IsItFragmentable(const G4F << 305 // virtual G4bool IsFragmentable(const G4FragmentingString * const string)=0; 316 // << 306 317 // If a string can not fragment, make last b 307 // If a string can not fragment, make last break into 2 hadrons 318 // virtual G4bool SplitLast(G4FragmentingStr 308 // virtual G4bool SplitLast(G4FragmentingString * string, 319 // G4KineticTrackVe 309 // G4KineticTrackVector * LeftVector, 320 // G4KineticTrackVe 310 // G4KineticTrackVector * RightVector)=0; 321 //-------------------------------------------- 311 //----------------------------------------------------------------------------- 322 // 312 // 323 // If a string can fragment, do the followin << 313 // If a string fragments, do the following 324 // 314 // 325 // For transver of a string to its CMS frame 315 // For transver of a string to its CMS frame 326 //-------------------------------------------- 316 //----------------------------------------------------------------------------- 327 317 328 G4ExcitedString *G4VLongitudinalStringDecay::C << 318 G4ExcitedString *G4VLongitudinalStringDecay::CPExcited(const G4ExcitedString & in) 329 { 319 { 330 G4Parton *Left=new G4Parton(*in.GetLeftParto 320 G4Parton *Left=new G4Parton(*in.GetLeftParton()); 331 G4Parton *Right=new G4Parton(*in.GetRightPar 321 G4Parton *Right=new G4Parton(*in.GetRightParton()); 332 return new G4ExcitedString(Left,Right,in.Get 322 return new G4ExcitedString(Left,Right,in.GetDirection()); 333 } 323 } 334 324 335 //-------------------------------------------- 325 //----------------------------------------------------------------------------- 336 326 337 G4ParticleDefinition * G4VLongitudinalStringDe << 327 G4KineticTrack * G4VLongitudinalStringDecay::Splitup( 338 << 328 G4FragmentingString *string, 339 { << 329 G4FragmentingString *&newString) 340 #ifdef debug_VStringDecay << 330 { 341 G4cout<<"VlongSD QuarkSplitup: quark ID "<< << 331 //G4cout<<"In G4VLong String Dec######################"<<G4endl; 342 #endif << 332 343 << 333 //... random choice of string end to use for creating the hadron (decay) 344 G4int IsParticle=(decay->GetPDGEncoding()>0 << 334 SideOfDecay = (G4UniformRand() < 0.5)? 1: -1; >> 335 if (SideOfDecay < 0) >> 336 { >> 337 string->SetLeftPartonStable(); >> 338 } else >> 339 { >> 340 string->SetRightPartonStable(); >> 341 } >> 342 >> 343 G4ParticleDefinition *newStringEnd; >> 344 G4ParticleDefinition * HadronDefinition; >> 345 if (string->DecayIsQuark()) >> 346 { >> 347 HadronDefinition= QuarkSplitup(string->GetDecayParton(), newStringEnd); >> 348 } else { >> 349 HadronDefinition= DiQuarkSplitup(string->GetDecayParton(), newStringEnd); >> 350 } >> 351 // create new String from old, ie. keep Left and Right order, but replace decay >> 352 >> 353 newString=new G4FragmentingString(*string,newStringEnd); // To store possible >> 354 // quark containt of new string >> 355 G4LorentzVector* HadronMomentum=SplitEandP(HadronDefinition, string, newString); 345 356 346 pDefPair QuarkPair = CreatePartonPair(IsPar << 357 delete newString; newString=0; // Uzhi 20.06.08 347 created = QuarkPair.second; << 358 >> 359 G4KineticTrack * Hadron =0; >> 360 if ( HadronMomentum != 0 ) { >> 361 >> 362 G4ThreeVector Pos; >> 363 Hadron = new G4KineticTrack(HadronDefinition, 0,Pos, *HadronMomentum); >> 364 >> 365 newString=new G4FragmentingString(*string,newStringEnd, >> 366 HadronMomentum); >> 367 >> 368 //G4cout<<"Out G4VLong String Dec######################"<<G4endl; >> 369 //G4cout<<"newString 4Mom"<<newString->Get4Momentum()<<G4endl; >> 370 //G4cout<<"newString Pl "<<newString->LightConePlus()<<G4endl; >> 371 //G4cout<<"newString Mi "<<newString->LightConeMinus()<<G4endl; >> 372 //G4cout<<"newString Pts "<<newString->StablePt()<<G4endl; >> 373 //G4cout<<"newString Ptd "<<newString->DecayPt()<<G4endl; >> 374 //G4cout<<"newString M2 "<<newString->Mass2()<<G4endl; >> 375 //G4cout<<"newString M2 "<<newString->Mass()<<G4endl; >> 376 delete HadronMomentum; >> 377 } >> 378 return Hadron; >> 379 } 348 380 349 DecayQuark = decay->GetPDGEncoding(); << 381 //-------------------------------------------------------------------------------------- 350 NewQuark = created->GetPDGEncoding(); << 351 382 352 #ifdef debug_VStringDecay << 383 G4ParticleDefinition * 353 G4cout<<"VlongSD QuarkSplitup: "<<decay->Ge << 384 G4VLongitudinalStringDecay::QuarkSplitup(G4ParticleDefinition* 354 G4cout<<"hadronizer->Build(QuarkPair.first, << 385 decay, G4ParticleDefinition *&created) 355 #endif << 386 { >> 387 G4int IsParticle=(decay->GetPDGEncoding()>0) ? -1 : +1; // if we have a quark, >> 388 // we need antiquark >> 389 // (or diquark) >> 390 pDefPair QuarkPair = CreatePartonPair(IsParticle); >> 391 created = QuarkPair.second; >> 392 return hadronizer->Build(QuarkPair.first, decay); >> 393 >> 394 } >> 395 >> 396 //----------------------------------------------------------------------------- >> 397 >> 398 G4ParticleDefinition *G4VLongitudinalStringDecay::DiQuarkSplitup( >> 399 G4ParticleDefinition* decay, >> 400 G4ParticleDefinition *&created) >> 401 { >> 402 //... can Diquark break or not? >> 403 if (G4UniformRand() < DiquarkBreakProb ){ >> 404 //... Diquark break >> 405 >> 406 G4int stableQuarkEncoding = decay->GetPDGEncoding()/1000; >> 407 G4int decayQuarkEncoding = (decay->GetPDGEncoding()/100)%10; >> 408 if (G4UniformRand() < 0.5) >> 409 { >> 410 G4int Swap = stableQuarkEncoding; >> 411 stableQuarkEncoding = decayQuarkEncoding; >> 412 decayQuarkEncoding = Swap; >> 413 } >> 414 >> 415 G4int IsParticle=(decayQuarkEncoding>0) ? -1 : +1; >> 416 // if we have a quark, we need antiquark) >> 417 pDefPair QuarkPair = CreatePartonPair(IsParticle,false); // no diquarks wanted >> 418 //... Build new Diquark >> 419 G4int QuarkEncoding=QuarkPair.second->GetPDGEncoding(); >> 420 G4int i10 = std::max(std::abs(QuarkEncoding), std::abs(stableQuarkEncoding)); >> 421 G4int i20 = std::min(std::abs(QuarkEncoding), std::abs(stableQuarkEncoding)); >> 422 G4int spin = (i10 != i20 && G4UniformRand() <= 0.5)? 1 : 3; >> 423 G4int NewDecayEncoding = -1*IsParticle*(i10 * 1000 + i20 * 100 + spin); >> 424 created = FindParticle(NewDecayEncoding); >> 425 G4ParticleDefinition * decayQuark=FindParticle(decayQuarkEncoding); >> 426 G4ParticleDefinition * had=hadronizer->Build(QuarkPair.first, decayQuark); >> 427 return had; >> 428 // return hadronizer->Build(QuarkPair.first, decayQuark); 356 429 357 return hadronizer->Build(QuarkPair.first, d << 430 } else { >> 431 //... Diquark does not break >> 432 >> 433 G4int IsParticle=(decay->GetPDGEncoding()>0) ? +1 : -1; >> 434 // if we have a diquark, we need quark) >> 435 pDefPair QuarkPair = CreatePartonPair(IsParticle,false); // no diquarks wanted >> 436 created = QuarkPair.second; >> 437 >> 438 G4ParticleDefinition * had=hadronizer->Build(QuarkPair.first, decay); >> 439 return had; >> 440 // return G4ParticleDefinition * had=hadronizer->Build(QuarkPair.first, decay); >> 441 } 358 } 442 } 359 443 360 //-------------------------------------------- 444 //----------------------------------------------------------------------------- 361 445 362 G4VLongitudinalStringDecay::pDefPair G4VLongit << 446 G4int G4VLongitudinalStringDecay::SampleQuarkFlavor(void) 363 CreatePartonPair(G4int NeedParticle,G4bool All << 447 { >> 448 return (1 + (int)(G4UniformRand()/StrangeSuppress)); >> 449 } >> 450 >> 451 //----------------------------------------------------------------------------- >> 452 >> 453 G4VLongitudinalStringDecay::pDefPair G4VLongitudinalStringDecay::CreatePartonPair(G4int NeedParticle,G4bool AllowDiquarks) 364 { 454 { 365 // NeedParticle = +1 for Particle, -1 for << 455 // NeedParticle = +1 for Particle, -1 for Antiparticle >> 456 366 if ( AllowDiquarks && G4UniformRand() < Di 457 if ( AllowDiquarks && G4UniformRand() < DiquarkSuppress ) 367 { 458 { 368 // Create a Diquark - AntiDiquark pair , 459 // Create a Diquark - AntiDiquark pair , first in pair is anti to IsParticle 369 #ifdef debug_VStringDecay << 460 G4int q1 = SampleQuarkFlavor(); 370 G4cout<<"VlongSD Create a Diquark - Anti << 461 G4int q2 = SampleQuarkFlavor(); 371 #endif << 462 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 463 // convention: quark with higher PDG number is first 379 PDGcode = (std::max(q1,q2) * 1000 + std: << 464 G4int PDGcode = (std::max(q1,q2) * 1000 + std::min(q1,q2) * 100 + spin) * NeedParticle; 380 << 381 return pDefPair (FindParticle(-PDGcode), 465 return pDefPair (FindParticle(-PDGcode),FindParticle(PDGcode)); >> 466 382 467 383 } else { 468 } else { 384 // Create a Quark - AntiQuark pair, firs 469 // 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 470 G4int PDGcode=SampleQuarkFlavor()*NeedParticle; 389 return pDefPair (FindParticle(PDGcode),F 471 return pDefPair (FindParticle(PDGcode),FindParticle(-PDGcode)); 390 } 472 } 391 } << 392 473 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 } 474 } 415 475 416 //-------------------------------------------- 476 //----------------------------------------------------------------------------- 417 477 >> 478 // G4ThreeVector G4VLongitudinalStringDecay::SampleQuarkPt() >> 479 // { >> 480 // G4double width_param= 2.0 * GeV*GeV; >> 481 // G4double R = G4UniformRand(); >> 482 // G4double Pt = std::sqrt(width_param*R/(1-R)); >> 483 // G4double phi = 2.*pi*G4UniformRand(); >> 484 // return G4ThreeVector(Pt * std::cos(phi),Pt * std::sin(phi),0); >> 485 // } >> 486 418 G4ThreeVector G4VLongitudinalStringDecay::Samp 487 G4ThreeVector G4VLongitudinalStringDecay::SampleQuarkPt(G4double ptMax) 419 { << 488 { 420 G4double Pt; 489 G4double Pt; 421 if ( ptMax < 0 ) { 490 if ( ptMax < 0 ) { 422 // sample full gaussian 491 // sample full gaussian 423 Pt = -G4Log(G4UniformRand()); << 492 Pt = -std::log(G4UniformRand()); 424 } else { 493 } else { 425 // sample in limited range 494 // sample in limited range 426 G4double q = ptMax/SigmaQT; << 495 Pt = -std::log(CLHEP::RandFlat::shoot(exp(-sqr(ptMax)/sqr(SigmaQT)), 1.)); 427 G4double ymin = (q > 20.) ? 0.0 : G4Exp( << 428 Pt = -G4Log(G4RandFlat::shoot(ymin, 1.)) << 429 } 496 } 430 Pt = SigmaQT * std::sqrt(Pt); 497 Pt = SigmaQT * std::sqrt(Pt); 431 G4double phi = 2.*pi*G4UniformRand(); 498 G4double phi = 2.*pi*G4UniformRand(); 432 return G4ThreeVector(Pt * std::cos(phi),Pt 499 return G4ThreeVector(Pt * std::cos(phi),Pt * std::sin(phi),0); 433 } << 500 } 434 501 435 //******************************************** 502 //****************************************************************************** 436 503 437 void G4VLongitudinalStringDecay::CalculateHadr << 504 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 { 505 { >> 506 >> 507 // `yo-yo` formation time >> 508 // const G4double kappa = 1.0 * GeV/fermi/4.; // Uzhi String tension 1.06.08 >> 509 G4double kappa = GetStringTensionParameter(); >> 510 //G4cout<<"Kappa "<<kappa<<G4endl; // Uzhi 20.06.08 >> 511 //G4int Uzhi; G4cin>>Uzhi; // Uzhi 20.06.08 >> 512 for(size_t c1 = 0; c1 < Hadrons->size(); c1++) >> 513 { 445 G4double SumPz = 0; 514 G4double SumPz = 0; 446 G4double SumE = 0; 515 G4double SumE = 0; 447 for (size_t c2 = 0; c2 < c1; c2++) << 516 for(size_t c2 = 0; c2 < c1; c2++) 448 { << 517 { 449 SumPz += Hadrons->operator[](c2)->Get 518 SumPz += Hadrons->operator[](c2)->Get4Momentum().pz(); 450 SumE += Hadrons->operator[](c2)->Get 519 SumE += Hadrons->operator[](c2)->Get4Momentum().e(); 451 } << 520 } 452 G4double HadronE = Hadrons->operator[]( 521 G4double HadronE = Hadrons->operator[](c1)->Get4Momentum().e(); 453 G4double HadronPz = Hadrons->operator[]( 522 G4double HadronPz = Hadrons->operator[](c1)->Get4Momentum().pz(); 454 Hadrons->operator[](c1)->SetFormationTim << 523 Hadrons->operator[](c1)->SetFormationTime((theInitialStringMass - 2.*SumPz + HadronE - HadronPz)/(2.*kappa)); 455 (theInitialStringMass - 2.*SumPz + Had << 524 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 525 Hadrons->operator[](c1)->SetPosition(aPosition); >> 526 } 459 } 527 } 460 } << 461 528 462 //-------------------------------------------- 529 //----------------------------------------------------------------------------- 463 530 >> 531 /* >> 532 void G4VLongitudinalStringDecay::CalculateHadronTimePosition(G4double theInitialStringMass, G4KineticTrackVector* Hadrons) >> 533 { >> 534 // 'constituent' formation time >> 535 const G4double kappa = 1.0 * GeV/fermi; >> 536 for(G4int c1 = 0; c1 < Hadrons->length(); c1++) >> 537 { >> 538 G4double SumPz = 0; >> 539 G4double SumE = 0; >> 540 for(G4int c2 = 0; c2 <= c1; c2++) >> 541 { >> 542 SumPz += Hadrons->at(c2)->Get4Momentum().pz(); >> 543 SumE += Hadrons->at(c2)->Get4Momentum().e(); >> 544 } >> 545 Hadrons->at(c1)->SetFormationTime((theInitialStringMass - 2.*SumPz)/(2.*kappa)); >> 546 G4ThreeVector aPosition(0, 0, (theInitialStringMass - 2.*SumE)/(2.*kappa)); >> 547 Hadrons->at(c1)->SetPosition(aPosition); >> 548 } >> 549 c1 = Hadrons->length()-1; >> 550 Hadrons->at(c1)->SetFormationTime(Hadrons->at(c1-1)->GetFormationTime()); >> 551 Hadrons->at(c1)->SetPosition(Hadrons->at(c1-1)->GetPosition()); >> 552 } >> 553 */ >> 554 >> 555 >> 556 //***************************************************************************** >> 557 464 void G4VLongitudinalStringDecay::SetSigmaTrans 558 void G4VLongitudinalStringDecay::SetSigmaTransverseMomentum(G4double aValue) 465 { 559 { 466 if ( PastInitPhase ) { << 560 if ( PastInitPhase ) { 467 throw G4HadronicException(__FILE__, __LIN << 561 throw G4HadronicException(__FILE__, __LINE__, "4VLongitudinalStringDecay::SetSigmaTransverseMomentum after FragmentString() not allowed"); 468 "G4VLongitudinalStringDecay::SetSigmaTr << 562 } else { 469 } else { << 563 SigmaQT = aValue; 470 SigmaQT = aValue; << 564 } 471 } << 472 } 565 } 473 566 474 //-------------------------------------------- 567 //---------------------------------------------------------------------------------------------------------- 475 568 476 void G4VLongitudinalStringDecay::SetStrangenes 569 void G4VLongitudinalStringDecay::SetStrangenessSuppression(G4double aValue) 477 { 570 { 478 StrangeSuppress = aValue; << 571 if ( PastInitPhase ) { >> 572 throw G4HadronicException(__FILE__, __LINE__, "4VLongitudinalStringDecay::SetStrangenessSuppression after FragmentString() not allowed"); >> 573 } else { >> 574 StrangeSuppress = aValue; >> 575 } 479 } 576 } 480 577 481 //-------------------------------------------- 578 //---------------------------------------------------------------------------------------------------------- 482 579 483 void G4VLongitudinalStringDecay::SetDiquarkSup 580 void G4VLongitudinalStringDecay::SetDiquarkSuppression(G4double aValue) 484 { 581 { 485 DiquarkSuppress = aValue; << 582 if ( PastInitPhase ) { >> 583 throw G4HadronicException(__FILE__, __LINE__, "4VLongitudinalStringDecay::SetDiquarkSuppression after FragmentString() not allowed"); >> 584 } else { >> 585 DiquarkSuppress = aValue; >> 586 } 486 } 587 } 487 588 488 //-------------------------------------------- 589 //---------------------------------------------------------------------------------------- 489 590 490 void G4VLongitudinalStringDecay::SetDiquarkBre 591 void G4VLongitudinalStringDecay::SetDiquarkBreakProbability(G4double aValue) 491 { 592 { 492 if ( PastInitPhase ) { << 593 if ( PastInitPhase ) { 493 throw G4HadronicException(__FILE__, __LINE << 594 throw G4HadronicException(__FILE__, __LINE__, "4VLongitudinalStringDecay::SetDiquarkBreakProbability after FragmentString() not allowed"); 494 "G4VLongitudinalStringDecay::SetDiquarkB << 595 } else { 495 } else { << 596 DiquarkBreakProb = aValue; 496 DiquarkBreakProb = aValue; << 597 } 497 } << 498 } 598 } 499 599 500 //-------------------------------------------- 600 //---------------------------------------------------------------------------------------------------------- 501 601 502 void G4VLongitudinalStringDecay::SetSpinThreeH << 602 void G4VLongitudinalStringDecay::SetVectorMesonProbability(G4double aValue) 503 { 603 { 504 if ( PastInitPhase ) { << 604 if ( PastInitPhase ) { 505 throw G4HadronicException(__FILE__, __LINE << 605 throw G4HadronicException(__FILE__, __LINE__, "G4VLongitudinalStringDecay::SetVectorMesonProbability after FragmentString() not allowed"); 506 "G4VLongitudinalStringDecay::SetSpinThre << 606 } else { 507 } else { << 607 pspin_meson = aValue; 508 pspin_barion = aValue; << 608 delete hadronizer; 509 delete hadronizer; << 609 hadronizer = new G4HadronBuilder(pspin_meson,pspin_barion, 510 hadronizer = new G4HadronBuilder( pspin_me << 610 scalarMesonMix,vectorMesonMix); 511 ProbEta_ << 611 } 512 } << 513 } 612 } 514 613 515 //-------------------------------------------- 614 //---------------------------------------------------------------------------------------------------------- 516 615 517 void G4VLongitudinalStringDecay::SetScalarMeso << 616 void G4VLongitudinalStringDecay::SetSpinThreeHalfBarionProbability(G4double aValue) 518 { 617 { 519 if ( PastInitPhase ) { << 618 if ( PastInitPhase ) { 520 throw G4HadronicException(__FILE__, __LINE << 619 throw G4HadronicException(__FILE__, __LINE__, "G4VLongitudinalStringDecay::SetSpinThreeHalfBarionProbability after FragmentString() not allowed"); 521 "G4VLongitudinalStringDecay::SetScalarMe << 620 } else { 522 } else { << 621 pspin_barion = aValue; 523 if ( aVector.size() < 6 ) << 622 delete hadronizer; 524 throw G4HadronicException(__FILE__, __LI << 623 hadronizer = new G4HadronBuilder(pspin_meson,pspin_barion, 525 "G4VLongitudinalStringDecay::SetScalar << 624 scalarMesonMix,vectorMesonMix); 526 scalarMesonMix[0] = aVector[0]; << 625 } 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 } 626 } 537 627 538 //-------------------------------------------- 628 //---------------------------------------------------------------------------------------------------------- 539 629 540 void G4VLongitudinalStringDecay::SetVectorMeso << 630 void G4VLongitudinalStringDecay::SetScalarMesonMixings(std::vector<G4double> aVector) 541 { << 542 if ( PastInitPhase ) { << 543 throw G4HadronicException(__FILE__, __LINE << 544 "G4VLongitudinalStringDecay::SetVectorMe << 545 } else { << 546 if ( aVector.size() < 6 ) << 547 throw G4HadronicException(__FILE__, __LI << 548 "G4VLongitudinalStringDecay::SetVector << 549 vectorMesonMix[0] = aVector[0]; << 550 vectorMesonMix[1] = aVector[1]; << 551 vectorMesonMix[2] = aVector[2]; << 552 vectorMesonMix[3] = aVector[3]; << 553 vectorMesonMix[4] = aVector[4]; << 554 vectorMesonMix[5] = aVector[5]; << 555 delete hadronizer; << 556 hadronizer = new G4HadronBuilder( pspin_me << 557 ProbEta_ << 558 } << 559 } << 560 << 561 //-------------------------------------------- << 562 << 563 void G4VLongitudinalStringDecay::SetProbCCbar( << 564 { << 565 ProbCCbar = aValue; << 566 ProbCB = ProbCCbar + ProbBBbar; << 567 } << 568 << 569 //-------------------------------------------- << 570 << 571 void G4VLongitudinalStringDecay::SetProbEta_c( << 572 { << 573 ProbEta_c = aValue; << 574 } << 575 << 576 //-------------------------------------------- << 577 << 578 void G4VLongitudinalStringDecay::SetProbBBbar( << 579 { 631 { 580 ProbBBbar = aValue; << 632 if ( PastInitPhase ) { 581 ProbCB = ProbCCbar + ProbBBbar; << 633 throw G4HadronicException(__FILE__, __LINE__, "G4VLongitudinalStringDecay::SetScalarMesonMixings after FragmentString() not allowed"); >> 634 } else { >> 635 if ( aVector.size() < 6 ) >> 636 throw G4HadronicException(__FILE__, __LINE__, "G4VLongitudinalStringDecay::SetScalarMesonMixings( argument Vector too small"); >> 637 scalarMesonMix[0] = aVector[0]; >> 638 scalarMesonMix[1] = aVector[1]; >> 639 scalarMesonMix[2] = aVector[2]; >> 640 scalarMesonMix[3] = aVector[3]; >> 641 scalarMesonMix[4] = aVector[4]; >> 642 scalarMesonMix[5] = aVector[5]; >> 643 delete hadronizer; >> 644 hadronizer = new G4HadronBuilder(pspin_meson,pspin_barion, >> 645 scalarMesonMix,vectorMesonMix); >> 646 } 582 } 647 } 583 648 584 //-------------------------------------------- << 649 //---------------------------------------------------------------------------------------------------------- 585 650 586 void G4VLongitudinalStringDecay::SetProbEta_b( << 651 void G4VLongitudinalStringDecay::SetVectorMesonMixings(std::vector<G4double> aVector) 587 { 652 { 588 ProbEta_b = aValue; << 653 if ( PastInitPhase ) { >> 654 throw G4HadronicException(__FILE__, __LINE__, "G4VLongitudinalStringDecay::SetVectorMesonMixings after FragmentString() not allowed"); >> 655 } else { >> 656 if ( aVector.size() < 6 ) >> 657 throw G4HadronicException(__FILE__, __LINE__, "G4VLongitudinalStringDecay::SetVectorMesonMixings( argument Vector too small"); >> 658 vectorMesonMix[0] = aVector[0]; >> 659 vectorMesonMix[1] = aVector[1]; >> 660 vectorMesonMix[2] = aVector[2]; >> 661 vectorMesonMix[3] = aVector[3]; >> 662 vectorMesonMix[4] = aVector[4]; >> 663 vectorMesonMix[5] = aVector[5]; >> 664 delete hadronizer; >> 665 hadronizer = new G4HadronBuilder(pspin_meson,pspin_barion, >> 666 scalarMesonMix,vectorMesonMix); >> 667 >> 668 } 589 } 669 } 590 670 591 //-------------------------------------------- 671 //------------------------------------------------------------------------------------------- 592 << 672 void G4VLongitudinalStringDecay::SetStringTensionParameter(G4double aValue)// Uzhi 20 June 08 593 void G4VLongitudinalStringDecay::SetStringTens << 594 { 673 { 595 Kappa = aValue * GeV/fermi; << 674 Kappa = aValue * GeV/fermi; 596 } 675 } 597 << 676 //************************************************************************************** 598 //-------------------------------------------- << 599 << 600 void G4VLongitudinalStringDecay::SetMinMasses( << 601 { << 602 // ------ For estimation of a minimal stri << 603 Mass_of_light_quark =140.*MeV; << 604 Mass_of_s_quark =500.*MeV; << 605 Mass_of_c_quark =1600.*MeV; << 606 Mass_of_b_quark =4500.*MeV; << 607 Mass_of_string_junction=720.*MeV; << 608 << 609 // ---------------- Determination of minim << 610 G4ParticleDefinition * hadron1; G4int C << 611 G4ParticleDefinition * hadron2; G4int C << 612 for (G4int i=1; i < 6; i++) { << 613 Code1 = 100*i + 10*1 + 1; << 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 << 629 // ---------------- Determination of minim << 630 G4ParticleDefinition * hadron3; << 631 G4int kfla, kflb; << 632 // MaxMass = -350.0*GeV; // If there wi << 633 << 634 for (G4int i=1; i < 6; i++) { //i=1 << 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 << 642 // Add d-quark << 643 Code2 = 1000*kfla + 100*kflb + << 644 if ( (j == 1) && (k==1)) Code2 = 1000*2 + << 645 << 646 hadron2 = G4ParticleTable::Get << 647 hadron3 = G4ParticleTable::Get << 648 << 649 if ((hadron2 == nullptr) && (h << 650 << 651 if ((hadron2 != nullptr) && (h << 652 if (hadron2->GetPDGMass() > << 653 }; << 654 << 655 if ((hadron2 != nullptr) && (h << 656 << 657 if ((hadron2 == nullptr) && (h << 658 << 659 minMassQDiQStr[i-1][j-1][k-1] << 660 } << 661 } << 662 } << 663 << 664 // ------ An estimated minimal string mass << 665 MinimalStringMass = 0.; << 666 MinimalStringMass2 = 0.; << 667 // q charges d u << 668 Qcharge[0] = -1; Qcharge[1] = 2; Qcharge[2 << 669 << 670 // For treating of small string decays << 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 << 759 kfla =0; kflb =0; << 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 << 782 // Delta- ddd - only 1114 << 783 Baryon[0][0][0][0] = 1114; BaryonWeight << 784 Baryon[0][0][0][1] = 0; BaryonWeight << 785 << 786 // Delta++ uuu - only 2224 << 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 << 927 // --------- Probabilities of q-qbar pair << 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 << 941 NumberOf_FS = 0; << 942 } << 943 << 944 // ------------------------------------------- << 945 << 946 void G4VLongitudinalStringDecay::SetMinimalStr << 947 { << 948 //MaxMass = -350.0*GeV; << 949 G4double EstimatedMass=MaxMass; << 950 << 951 G4ParticleDefinition* LeftParton = st << 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 << 968 G4int Qleft =std::abs(string->GetLeftP << 969 G4int Qright=std::abs(string->GetRight << 970 << 971 if ((Qleft < 6) && (Qright < 6)) { / << 972 EstimatedMass=minMassQQbarStr[Qleft- << 973 MinimalStringMass=EstimatedMass; << 974 SetMinimalStringMass2(EstimatedMass) << 975 return; << 976 } << 977 << 978 if ((Qleft < 6) && (Qright > 1000)) { << 979 G4int q1=Qright/1000; << 980 G4int q2=(Qright/100)%10; << 981 EstimatedMass=minMassQDiQStr[Qleft-1 << 982 MinimalStringMass=EstimatedMass; << 983 SetMinimalStringMass2(EstimatedMass) << 984 return; << 985 } << 986 << 987 if ((Qleft > 1000) && (Qright < 6)) { << 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 << 996 // DiQuark - Anti DiQuark string ----- << 997 << 998 G4double StringM=string->Get4Momentum().mag( << 999 << 1000 #ifdef debug_LUNDfragmentation << 1001 // G4cout<<"MinStringMass// Input Str << 1002 #endif << 1003 << 1004 G4int q1= Qleft/1000 ; << 1005 G4int q2=(Qleft/100)%10 ; << 1006 << 1007 G4int q3= Qright/1000 ; << 1008 G4int q4=(Qright/100)%10; << 1009 << 1010 // -------------- 2 baryon production << 1011 << 1012 G4double EstimatedMass1 = minMassQDiQ << 1013 G4double EstimatedMass2 = minMassQDiQ << 1014 // Mass is negative if there is no co << 1015 << 1016 if ( (EstimatedMass1 > 0.) && (Estima << 1017 EstimatedMass = EstimatedMass1 + E << 1018 if ( StringM > EstimatedMass ) { << 1019 MinimalStringMass=EstimatedMass << 1020 SetMinimalStringMass2(Estimated << 1021 return; << 1022 } << 1023 } << 1024 << 1025 if ( (EstimatedMass1 < 0.) && (Estima << 1026 EstimatedMass = MaxMass; << 1027 MinimalStringMass=EstimatedMass; << 1028 SetMinimalStringMass2(EstimatedMas << 1029 return; << 1030 } << 1031 << 1032 if ( (EstimatedMass1 > 0.) && (Estima << 1033 EstimatedMass = EstimatedMass1; << 1034 MinimalStringMass=EstimatedMass; << 1035 SetMinimalStringMass2(EstimatedMas << 1036 return; << 1037 } << 1038 << 1039 // if ( EstimatedMass >= StringM << 1040 // ------------- Re-orangement ------ << 1041 EstimatedMass=std::min(minMassQQbarSt << 1042 minMassQQbarSt << 1043 << 1044 // In principle, re-arrangement and 2 << 1045 // More physics consideration is need << 1046 << 1047 MinimalStringMass=EstimatedMass; << 1048 SetMinimalStringMass2(EstimatedMass); << 1049 << 1050 return; << 1051 } << 1052 << 1053 //------------------------------------------- << 1054 << 1055 void G4VLongitudinalStringDecay::SetMinimalSt << 1056 { << 1057 MinimalStringMass2=aValue * aValue; << 1058 } << 1059 677 1060 678