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, 10-Jul-1998 32 // ------------------------------------------- 32 // ----------------------------------------------------------------------------- 33 #include "G4LundStringFragmentation.hh" 33 #include "G4LundStringFragmentation.hh" 34 #include "G4PhysicalConstants.hh" 34 #include "G4PhysicalConstants.hh" 35 #include "G4SystemOfUnits.hh" 35 #include "G4SystemOfUnits.hh" 36 #include "Randomize.hh" 36 #include "Randomize.hh" 37 #include "G4FragmentingString.hh" 37 #include "G4FragmentingString.hh" 38 #include "G4DiQuarks.hh" 38 #include "G4DiQuarks.hh" 39 #include "G4Quarks.hh" 39 #include "G4Quarks.hh" 40 #include "G4HadronicParameters.hh" 40 #include "G4HadronicParameters.hh" 41 #include "G4Exp.hh" 41 #include "G4Exp.hh" 42 #include "G4Pow.hh" 42 #include "G4Pow.hh" 43 43 44 //#define debug_LUNDfragmentation 44 //#define debug_LUNDfragmentation 45 45 46 // Class G4LundStringFragmentation 46 // Class G4LundStringFragmentation 47 //******************************************** 47 //************************************************************************************* 48 48 49 G4LundStringFragmentation::G4LundStringFragmen 49 G4LundStringFragmentation::G4LundStringFragmentation() 50 : G4VLongitudinalStringDecay("LundStringFrag 50 : G4VLongitudinalStringDecay("LundStringFragmentation") 51 { 51 { 52 SetMassCut(210.*MeV); // Mpi + Delta 52 SetMassCut(210.*MeV); // Mpi + Delta 53 // For ProduceOneH 53 // For ProduceOneHadron it is required 54 // that no one pi- 54 // that no one pi-meson can be produced. 55 SigmaQT = 0.435 * GeV; 55 SigmaQT = 0.435 * GeV; 56 Tmt = 190.0 * MeV; 56 Tmt = 190.0 * MeV; 57 57 58 SetStringTensionParameter(1.*GeV/fermi); 58 SetStringTensionParameter(1.*GeV/fermi); 59 SetDiquarkBreakProbability(0.3); 59 SetDiquarkBreakProbability(0.3); 60 60 61 SetStrangenessSuppression((1.0 - 0.12)/2.0 61 SetStrangenessSuppression((1.0 - 0.12)/2.0); 62 SetDiquarkSuppression(0.07); 62 SetDiquarkSuppression(0.07); 63 63 64 // Check if charmed and bottom hadrons are 64 // Check if charmed and bottom hadrons are enabled: if this is the case, then 65 // set the non-zero probabilities for c-cb 65 // set the non-zero probabilities for c-cbar and b-bbar creation from the vacuum, 66 // else set them to 0.0. If these probabil 66 // else set them to 0.0. If these probabilities are/aren't zero then charmed or bottom 67 // hadrons can't/can be created during the 67 // hadrons can't/can be created during the string fragmentation of ordinary 68 // (i.e. not heavy) projectile hadron nucl 68 // (i.e. not heavy) projectile hadron nuclear reactions. 69 if ( G4HadronicParameters::Instance()->Ena 69 if ( G4HadronicParameters::Instance()->EnableBCParticles() ) { 70 SetProbCCbar(0.0002); // According to O 70 SetProbCCbar(0.0002); // According to O.I. Piskunova Yad. Fiz. 56 (1993) 1094; tuned by Uzhi Oct. 2022 71 SetProbBBbar(5.0e-5); // According to O 71 SetProbBBbar(5.0e-5); // According to O.I. Piskunova Yad. Fiz. 56 (1993) 1094 72 } else { 72 } else { 73 SetProbCCbar(0.0); 73 SetProbCCbar(0.0); 74 SetProbBBbar(0.0); 74 SetProbBBbar(0.0); 75 } 75 } 76 76 77 SetMinMasses(); // For treating of small 77 SetMinMasses(); // For treating of small string decays 78 } 78 } 79 79 80 //-------------------------------------------- 80 //-------------------------------------------------------------------------------------- 81 81 82 G4KineticTrackVector* G4LundStringFragmentatio 82 G4KineticTrackVector* G4LundStringFragmentation::FragmentString(const G4ExcitedString& theString) 83 { 83 { 84 // Can no longer modify Parameters for Fragm 84 // Can no longer modify Parameters for Fragmentation. 85 85 86 PastInitPhase=true; 86 PastInitPhase=true; 87 87 88 G4FragmentingString aString(theString); 88 G4FragmentingString aString(theString); 89 SetMinimalStringMass(&aString); 89 SetMinimalStringMass(&aString); 90 90 91 #ifdef debug_LUNDfragmentation 91 #ifdef debug_LUNDfragmentation 92 G4cout<<G4endl<<"LUND StringFragmentat 92 G4cout<<G4endl<<"LUND StringFragmentation ------------------------------------"<<G4endl; 93 G4cout<<G4endl<<"LUND StringFragm: Str 93 G4cout<<G4endl<<"LUND StringFragm: String Mass " 94 <<theString.Get4Momentum 94 <<theString.Get4Momentum().mag()<<G4endl 95 <<"4Mom "<<theString.Get 95 <<"4Mom "<<theString.Get4Momentum()<<G4endl 96 <<"--------------------- 96 <<"------------------------------------"<<G4endl; 97 G4cout<<"String ends Direct "<<theStri 97 G4cout<<"String ends Direct "<<theString.GetLeftParton()->GetPDGcode()<<" " 98 <<theStri 98 <<theString.GetRightParton()->GetPDGcode()<<" " 99 <<theStri 99 <<theString.GetDirection()<< G4endl; 100 G4cout<<"Left mom "<<theString.GetLef 100 G4cout<<"Left mom "<<theString.GetLeftParton()->Get4Momentum()<<G4endl; 101 G4cout<<"Right mom "<<theString.GetRig 101 G4cout<<"Right mom "<<theString.GetRightParton()->Get4Momentum()<<G4endl<<G4endl; 102 G4cout<<"Check for Fragmentation "<<G4 102 G4cout<<"Check for Fragmentation "<<G4endl; 103 #endif 103 #endif 104 104 105 G4KineticTrackVector * LeftVector(0); 105 G4KineticTrackVector * LeftVector(0); 106 106 107 if (!aString.IsAFourQuarkString() && !IsItFr 107 if (!aString.IsAFourQuarkString() && !IsItFragmentable(&aString)) 108 { 108 { 109 #ifdef debug_LUNDfragmentation 109 #ifdef debug_LUNDfragmentation 110 G4cout<<"Non fragmentable - th 110 G4cout<<"Non fragmentable - the string is converted to one hadron "<<G4endl; 111 #endif 111 #endif 112 // SetMassCut(210.*MeV); // F 112 // SetMassCut(210.*MeV); // For ProduceOneHadron it is required 113 // t 113 // that no one pi-meson can be produced. 114 114 115 G4double Mcut = GetMassCut(); 115 G4double Mcut = GetMassCut(); 116 SetMassCut(10000.*MeV); 116 SetMassCut(10000.*MeV); 117 LeftVector=ProduceOneHadron(&theString); 117 LeftVector=ProduceOneHadron(&theString); 118 SetMassCut(Mcut); 118 SetMassCut(Mcut); 119 119 120 if ( LeftVector ) 120 if ( LeftVector ) 121 { 121 { 122 if ( LeftVector->size() > 0) 122 if ( LeftVector->size() > 0) 123 { 123 { 124 LeftVector->operator[](0)->SetForm 124 LeftVector->operator[](0)->SetFormationTime(theString.GetTimeOfCreation()); 125 LeftVector->operator[](0)->SetPosi 125 LeftVector->operator[](0)->SetPosition(theString.GetPosition()); 126 } 126 } 127 if (LeftVector->size() > 1) 127 if (LeftVector->size() > 1) 128 { 128 { 129 // 2 hadrons created from qq-qqbar 129 // 2 hadrons created from qq-qqbar are stored 130 LeftVector->operator[](1)->SetFormationT 130 LeftVector->operator[](1)->SetFormationTime(theString.GetTimeOfCreation()); 131 LeftVector->operator[](1)->SetPosition(t 131 LeftVector->operator[](1)->SetPosition(theString.GetPosition()); 132 } 132 } 133 } 133 } 134 return LeftVector; 134 return LeftVector; 135 } 135 } 136 136 137 #ifdef debug_LUNDfragmentation 137 #ifdef debug_LUNDfragmentation 138 G4cout<<"The string will be fragmented 138 G4cout<<"The string will be fragmented. "<<G4endl; 139 #endif 139 #endif 140 140 141 // The string can fragment. At least two par 141 // The string can fragment. At least two particles can be produced. 142 LeftVector =new G4KineticTrackVec 142 LeftVector =new G4KineticTrackVector; 143 G4KineticTrackVector * RightVector=new G4Kin 143 G4KineticTrackVector * RightVector=new G4KineticTrackVector; 144 144 145 G4bool success = Loop_toFragmentString 145 G4bool success = Loop_toFragmentString(theString, LeftVector, RightVector); 146 146 147 if ( ! success ) 147 if ( ! success ) 148 { 148 { 149 std::for_each(LeftVector->begin(), LeftVec 149 std::for_each(LeftVector->begin(), LeftVector->end(), DeleteKineticTrack()); 150 LeftVector->clear(); 150 LeftVector->clear(); 151 std::for_each(RightVector->begin(), RightV 151 std::for_each(RightVector->begin(), RightVector->end(), DeleteKineticTrack()); 152 delete RightVector; 152 delete RightVector; 153 return LeftVector; 153 return LeftVector; 154 } 154 } 155 155 156 // Join Left- and RightVector into LeftVecto 156 // Join Left- and RightVector into LeftVector in correct order. 157 while (!RightVector->empty()) 157 while (!RightVector->empty()) 158 { 158 { 159 LeftVector->push_back(RightVector->back()) 159 LeftVector->push_back(RightVector->back()); 160 RightVector->erase(RightVector->end()-1); 160 RightVector->erase(RightVector->end()-1); 161 } 161 } 162 delete RightVector; 162 delete RightVector; 163 163 164 return LeftVector; 164 return LeftVector; 165 } 165 } 166 166 167 //-------------------------------------------- 167 //---------------------------------------------------------------------------------- 168 168 169 G4bool G4LundStringFragmentation::IsItFragment 169 G4bool G4LundStringFragmentation::IsItFragmentable(const G4FragmentingString * const string) 170 { 170 { 171 SetMinimalStringMass(string); 171 SetMinimalStringMass(string); 172 //G4cout<<"MinM StrM "<<MinimalStringM 172 //G4cout<<"MinM StrM "<<MinimalStringMass<<" "<< string->Get4Momentum().mag()<<G4endl; 173 173 174 return std::abs(MinimalStringMass) < string- 174 return std::abs(MinimalStringMass) < string->Get4Momentum().mag(); 175 175 176 //MinimalStringMass is negative and la 176 //MinimalStringMass is negative and large for a string with unknown particles in a final 2-particle decay. 177 } 177 } 178 178 179 //-------------------------------------------- 179 //---------------------------------------------------------------------------------------- 180 180 181 G4bool G4LundStringFragmentation::Loop_toFragm 181 G4bool G4LundStringFragmentation::Loop_toFragmentString( const G4ExcitedString &theString, 182 182 G4KineticTrackVector * & LeftVector, 183 183 G4KineticTrackVector * & RightVector ) 184 { 184 { 185 #ifdef debug_LUNDfragmentation 185 #ifdef debug_LUNDfragmentation 186 G4cout<<"Loop_toFrag "<<theString.GetL 186 G4cout<<"Loop_toFrag "<<theString.GetLeftParton()->GetPDGcode()<<" " 187 <<theString.GetL 187 <<theString.GetLeftParton()->Get4Momentum()<<G4endl 188 <<" "<<theString.GetR 188 <<" "<<theString.GetRightParton()->GetPDGcode()<<" " 189 <<theString.GetR 189 <<theString.GetRightParton()->Get4Momentum()<<G4endl 190 <<"Direction "<<theString.GetD 190 <<"Direction "<<theString.GetDirection()<< G4endl; 191 #endif 191 #endif 192 192 193 G4LorentzRotation toCmsI, toObserverFr 193 G4LorentzRotation toCmsI, toObserverFrameI; 194 194 195 G4bool final_success=false; 195 G4bool final_success=false; 196 G4bool inner_success=true; 196 G4bool inner_success=true; 197 197 198 G4int attempt=0; 198 G4int attempt=0; 199 199 200 while ( ! final_success && attempt++ < Strin 200 while ( ! final_success && attempt++ < StringLoopInterrupt ) 201 { // If the string fragmentation does 201 { // If the string fragmentation does not be happend, 202 // repeat the fragmentation. 202 // repeat the fragmentation. 203 203 204 G4FragmentingString *currentSt 204 G4FragmentingString *currentString = new G4FragmentingString(theString); 205 toCmsI = currentString->Transf 205 toCmsI = currentString->TransformToAlignedCms(); 206 toObserverFrameI = toCmsI.inve 206 toObserverFrameI = toCmsI.inverse(); 207 207 208 G4LorentzRotation toCms, toObserverFrame; 208 G4LorentzRotation toCms, toObserverFrame; 209 209 210 //G4cout<<"Main loop start whilecounter "< 210 //G4cout<<"Main loop start whilecounter "<<attempt<<G4endl; 211 211 212 // Cleaning up the previously produced had 212 // Cleaning up the previously produced hadrons 213 std::for_each(LeftVector->begin(), LeftVec 213 std::for_each(LeftVector->begin(), LeftVector->end(), DeleteKineticTrack()); 214 LeftVector->clear(); 214 LeftVector->clear(); 215 std::for_each(RightVector->begin(), RightV 215 std::for_each(RightVector->begin(), RightVector->end(), DeleteKineticTrack()); 216 RightVector->clear(); 216 RightVector->clear(); 217 217 218 // Main fragmentation loop until the strin 218 // Main fragmentation loop until the string will not be able to fragment 219 inner_success=true; // set false on failu 219 inner_success=true; // set false on failure. 220 const G4int maxNumberOfLoops = 220 const G4int maxNumberOfLoops = 1000; 221 G4int loopCounter = -1; 221 G4int loopCounter = -1; 222 222 223 while ( (! StopFragmenting(currentString)) 223 while ( (! StopFragmenting(currentString)) && ++loopCounter < maxNumberOfLoops ) 224 { // Split current string into hadro 224 { // Split current string into hadron + new string 225 #ifdef debug_LUNDfragm 225 #ifdef debug_LUNDfragmentation 226 G4cout<<"The string wi 226 G4cout<<"The string will fragment. "<<G4endl;; 227 //G4cout<<"1 "<<curren 227 //G4cout<<"1 "<<currentString->GetDecayDirection()<<G4endl; 228 #endif 228 #endif 229 G4FragmentingString *newString=0; // us 229 G4FragmentingString *newString=0; // used as output from SplitUp. 230 230 231 toCms=currentString->TransformToAlignedC 231 toCms=currentString->TransformToAlignedCms(); 232 toObserverFrame= toCms 232 toObserverFrame= toCms.inverse(); 233 233 234 #ifdef debug_LUNDfragm 234 #ifdef debug_LUNDfragmentation 235 //G4cout<<"CMS Left m 235 //G4cout<<"CMS Left mom "<<currentString->GetPleft()<<G4endl; 236 //G4cout<<"CMS Right m 236 //G4cout<<"CMS Right mom "<<currentString->GetPright()<<G4endl; 237 //G4cout<<"CMS String 237 //G4cout<<"CMS String M "<<currentString->GetPstring()<<G4endl; 238 #endif 238 #endif 239 239 240 G4KineticTrack * Hadron=Splitup(currentS 240 G4KineticTrack * Hadron=Splitup(currentString,newString); 241 241 242 if ( Hadron != 0 ) // Store the hadron 242 if ( Hadron != 0 ) // Store the hadron 243 { 243 { 244 #ifdef debug_L 244 #ifdef debug_LUNDfragmentation 245 G4cout<<"Hadro 245 G4cout<<"Hadron prod at fragm. "<<Hadron->GetDefinition()->GetParticleName()<<G4endl; 246 //G4cout<<"2 " 246 //G4cout<<"2 "<<currentString->GetDecayDirection()<<G4endl; 247 #endif 247 #endif 248 248 249 Hadron->Set4Momentum(toObserverFrame*H 249 Hadron->Set4Momentum(toObserverFrame*Hadron->Get4Momentum()); 250 250 251 G4double TimeOftheStringCreation=theSt 251 G4double TimeOftheStringCreation=theString.GetTimeOfCreation(); 252 G4ThreeVector PositionOftheStringCreat 252 G4ThreeVector PositionOftheStringCreation(theString.GetPosition()); 253 253 254 G4LorentzVector Coordinate(Hadron->Get 254 G4LorentzVector Coordinate(Hadron->GetPosition(), Hadron->GetFormationTime()); 255 G4LorentzVector Momentum = toObserverF 255 G4LorentzVector Momentum = toObserverFrame*Coordinate; 256 Hadron->SetFormationTime(TimeOftheStri 256 Hadron->SetFormationTime(TimeOftheStringCreation + Momentum.e() - fermi/c_light); 257 G4ThreeVector aPosition(Momentum.vect( 257 G4ThreeVector aPosition(Momentum.vect()); 258 Hadron->SetPosition(PositionOftheStrin 258 Hadron->SetPosition(PositionOftheStringCreation+aPosition); 259 259 260 // Open to pro 260 // Open to protect hadron production at fragmentation 261 if ( currentString->GetDecayDirection( 261 if ( currentString->GetDecayDirection() > 0 ) 262 { 262 { 263 LeftVector->push_back(Hadron); 263 LeftVector->push_back(Hadron); 264 } else 264 } else 265 { 265 { 266 RightVector->push_back(Hadron); 266 RightVector->push_back(Hadron); 267 } 267 } 268 delete currentString; 268 delete currentString; 269 currentString=newString; 269 currentString=newString; 270 } else { 270 } else { 271 if ( newString ) del 271 if ( newString ) delete newString; 272 } 272 } 273 273 274 currentString->Lorentz 274 currentString->LorentzRotate(toObserverFrame); 275 }; 275 }; 276 276 277 if ( loopCounter >= maxNumberO 277 if ( loopCounter >= maxNumberOfLoops ) { 278 inner_success=false; 278 inner_success=false; 279 } 279 } 280 280 281 // Split remaining string into 2 final had 281 // Split remaining string into 2 final hadrons. 282 #ifdef debug_LUNDfragmentation 282 #ifdef debug_LUNDfragmentation 283 if (inner_success) G4cout<<"Sp 283 if (inner_success) G4cout<<"Split remaining string into 2 final hadrons."<<G4endl; 284 #endif 284 #endif 285 285 286 if ( inner_success && SplitLast(currentStr 286 if ( inner_success && SplitLast(currentString, LeftVector, RightVector) ) // Close to protect Last Str. Decay 287 { 287 { 288 final_success = true; 288 final_success = true; 289 } 289 } 290 290 291 delete currentString; 291 delete currentString; 292 } // End of the loop where we try to fragme 292 } // End of the loop where we try to fragment the string. 293 293 294 G4int sign = +1; 294 G4int sign = +1; 295 if ( theString.GetDirection() < 0 ) si 295 if ( theString.GetDirection() < 0 ) sign = -1; 296 for ( unsigned int hadronI = 0; hadron 296 for ( unsigned int hadronI = 0; hadronI < LeftVector->size(); ++hadronI ) { 297 G4LorentzVector Tmp = LeftVector->o 297 G4LorentzVector Tmp = LeftVector->operator[](hadronI)->Get4Momentum(); 298 Tmp.setZ(sign*Tmp.getZ()); 298 Tmp.setZ(sign*Tmp.getZ()); 299 Tmp *= toObserverFrameI; 299 Tmp *= toObserverFrameI; 300 LeftVector->operator[](hadronI)->Se 300 LeftVector->operator[](hadronI)->Set4Momentum(Tmp); 301 } 301 } 302 for ( unsigned int hadronI = 0; hadron 302 for ( unsigned int hadronI = 0; hadronI < RightVector->size(); ++hadronI ) { 303 G4LorentzVector Tmp = RightVector-> 303 G4LorentzVector Tmp = RightVector->operator[](hadronI)->Get4Momentum(); 304 Tmp.setZ(sign*Tmp.getZ()); 304 Tmp.setZ(sign*Tmp.getZ()); 305 Tmp *= toObserverFrameI; 305 Tmp *= toObserverFrameI; 306 RightVector->operator[](hadronI)->S 306 RightVector->operator[](hadronI)->Set4Momentum(Tmp); 307 } 307 } 308 308 309 return final_success; 309 return final_success; 310 } 310 } 311 311 312 //-------------------------------------------- 312 //---------------------------------------------------------------------------------------- 313 313 314 G4bool G4LundStringFragmentation::StopFragment 314 G4bool G4LundStringFragmentation::StopFragmenting(const G4FragmentingString * const string) 315 { 315 { 316 SetMinimalStringMass(string); 316 SetMinimalStringMass(string); 317 317 318 if ( MinimalStringMass < 0.) return true; 318 if ( MinimalStringMass < 0.) return true; 319 319 320 if (string->IsAFourQuarkString()) 320 if (string->IsAFourQuarkString()) 321 { 321 { 322 return G4UniformRand() < G4Exp(-0.0005*(st 322 return G4UniformRand() < G4Exp(-0.0005*(string->Mass() - MinimalStringMass)); 323 } else { 323 } else { 324 324 325 if (MinimalStringMass < 0.0 ) 325 if (MinimalStringMass < 0.0 ) return false; // For a string with di-quark having c or b quarks and s, c, b quarks 326 326 327 G4bool Result = G4UniformRand() < 327 G4bool Result = G4UniformRand() < 328 G4Exp(-0.66e-6*(string->Mass()*string- 328 G4Exp(-0.66e-6*(string->Mass()*string->Mass() - MinimalStringMass*MinimalStringMass)); 329 // G4bool Result = string->Mas 329 // G4bool Result = string->Mass() < MinimalStringMass + 150.*MeV*G4UniformRand(); // a'la LUND 330 330 331 #ifdef debug_LUNDfragmentation 331 #ifdef debug_LUNDfragmentation 332 G4cout<<"StopFragmenting Minim 332 G4cout<<"StopFragmenting MinimalStringMass string->Mass() "<<MinimalStringMass 333 <<" "<<string->Mass()<<G 333 <<" "<<string->Mass()<<G4endl; 334 G4cout<<"StopFragmenting - Yes 334 G4cout<<"StopFragmenting - Yes/No "<<Result<<G4endl; 335 #endif 335 #endif 336 return Result; 336 return Result; 337 } 337 } 338 } 338 } 339 339 340 //-------------------------------------------- 340 //----------------------------------------------------------------------------- 341 341 342 G4KineticTrack * G4LundStringFragmentation::Sp 342 G4KineticTrack * G4LundStringFragmentation::Splitup(G4FragmentingString *string, 343 G4Fragmentin 343 G4FragmentingString *&newString) 344 { 344 { 345 #ifdef debug_LUNDfragmentation 345 #ifdef debug_LUNDfragmentation 346 G4cout<<G4endl; 346 G4cout<<G4endl; 347 G4cout<<"Start SplitUP ================ 347 G4cout<<"Start SplitUP ========================="<<G4endl; 348 G4cout<<"String partons: " <<string->Ge 348 G4cout<<"String partons: " <<string->GetLeftParton()->GetPDGEncoding()<<" " 349 <<string->Ge 349 <<string->GetRightParton()->GetPDGEncoding()<<" " 350 <<"Direction " <<string->Ge 350 <<"Direction " <<string->GetDecayDirection()<<G4endl; 351 #endif 351 #endif 352 352 353 //... random choice of string end to us 353 //... random choice of string end to use for creating the hadron (decay) 354 G4int SideOfDecay = (G4UniformRand() < 354 G4int SideOfDecay = (G4UniformRand() < 0.5)? 1: -1; 355 if (SideOfDecay < 0) 355 if (SideOfDecay < 0) 356 { 356 { 357 string->SetLeftPartonStable(); 357 string->SetLeftPartonStable(); 358 } else 358 } else 359 { 359 { 360 string->SetRightPartonStable(); 360 string->SetRightPartonStable(); 361 } 361 } 362 362 363 G4ParticleDefinition *newStringEnd; 363 G4ParticleDefinition *newStringEnd; 364 G4ParticleDefinition * HadronDefinition 364 G4ParticleDefinition * HadronDefinition; 365 365 366 G4double StringMass=string->Mass(); 366 G4double StringMass=string->Mass(); 367 367 368 G4double ProbDqADq = GetDiquarkSuppress 368 G4double ProbDqADq = GetDiquarkSuppress(); 369 G4double ProbSaS = 1.0 - 2.0 * GetStr 369 G4double ProbSaS = 1.0 - 2.0 * GetStrangeSuppress(); 370 370 371 #ifdef debug_LUNDfragmentation 371 #ifdef debug_LUNDfragmentation 372 G4cout<<"StrMass DiquarkSuppression 372 G4cout<<"StrMass DiquarkSuppression "<<StringMass<<" "<<GetDiquarkSuppress()<<G4endl; 373 #endif 373 #endif 374 374 375 G4int NumberOfpossibleBaryons = 2; 375 G4int NumberOfpossibleBaryons = 2; 376 376 377 if (string->GetLeftParton()->GetParticl 377 if (string->GetLeftParton()->GetParticleSubType() != "quark") NumberOfpossibleBaryons++; 378 if (string->GetRightParton()->GetPartic 378 if (string->GetRightParton()->GetParticleSubType() != "quark") NumberOfpossibleBaryons++; 379 379 380 G4double ActualProb = ProbDqADq ; 380 G4double ActualProb = ProbDqADq ; 381 ActualProb *= (1.0-G4Pow::GetInstance() 381 ActualProb *= (1.0-G4Pow::GetInstance()->powA(NumberOfpossibleBaryons*1400.0/StringMass, 8.0)); 382 if(ActualProb <0.0) ActualProb = 0.; 382 if(ActualProb <0.0) ActualProb = 0.; 383 383 384 SetDiquarkSuppression(ActualProb); 384 SetDiquarkSuppression(ActualProb); 385 385 386 G4double Mth = 1250.0; 386 G4double Mth = 1250.0; // 2 Mk + Mpi 387 if ( NumberOfpossibleBaryons == 3 ){Mth 387 if ( NumberOfpossibleBaryons == 3 ){Mth = 2520.0;} // Mlambda/Msigma + Mk + Mpi 388 else if ( NumberOfpossibleBaryons == 4 388 else if ( NumberOfpossibleBaryons == 4 ){Mth = 2380.0;} // 2 Mlambda/Msigma + Mk + Mpi 389 else {} 389 else {} 390 390 391 ActualProb = ProbSaS; 391 ActualProb = ProbSaS; 392 ActualProb *= (1.0 - G4Pow::GetInstance 392 ActualProb *= (1.0 - G4Pow::GetInstance()->powA( Mth/StringMass, 2.5 )); 393 if ( ActualProb < 0.0 ) ActualProb = 0. 393 if ( ActualProb < 0.0 ) ActualProb = 0.0; 394 SetStrangenessSuppression((1.0-ActualPr 394 SetStrangenessSuppression((1.0-ActualProb)/2.0); 395 395 396 #ifdef debug_LUNDfragmentation 396 #ifdef debug_LUNDfragmentation 397 G4cout<<"StrMass DiquarkSuppression cor 397 G4cout<<"StrMass DiquarkSuppression corrected "<<StringMass<<" "<<GetDiquarkSuppress()<<G4endl; 398 #endif 398 #endif 399 399 400 if (string->DecayIsQuark()) 400 if (string->DecayIsQuark()) 401 { 401 { 402 HadronDefinition= QuarkSplitup(strin 402 HadronDefinition= QuarkSplitup(string->GetDecayParton(), newStringEnd); 403 } else { 403 } else { 404 HadronDefinition= DiQuarkSplitup(str 404 HadronDefinition= DiQuarkSplitup(string->GetDecayParton(), newStringEnd); 405 } 405 } 406 406 407 SetDiquarkSuppression(ProbDqADq); 407 SetDiquarkSuppression(ProbDqADq); 408 SetStrangenessSuppression((1.0-ProbSaS) 408 SetStrangenessSuppression((1.0-ProbSaS)/2.0); 409 409 410 if ( HadronDefinition == NULL ) { G4Kin 410 if ( HadronDefinition == NULL ) { G4KineticTrack * Hadron =0; return Hadron; } 411 411 412 #ifdef debug_LUNDfragmentation 412 #ifdef debug_LUNDfragmentation 413 G4cout<<"The parton "<<string->GetDecay 413 G4cout<<"The parton "<<string->GetDecayParton()->GetPDGEncoding()<<" " 414 <<" produces hadron "<<HadronDefi 414 <<" produces hadron "<<HadronDefinition->GetParticleName() 415 <<" and is transformed to "<<newS 415 <<" and is transformed to "<<newStringEnd->GetPDGEncoding()<<G4endl; 416 G4cout<<"The side of the string decay L 416 G4cout<<"The side of the string decay Left/Right (1/-1) "<<SideOfDecay<<G4endl; 417 #endif 417 #endif 418 // create new String from old, ie. keep 418 // create new String from old, ie. keep Left and Right order, but replace decay 419 419 420 if ( newString ) delete newString; 420 if ( newString ) delete newString; 421 421 422 newString=new G4FragmentingString(*stri 422 newString=new G4FragmentingString(*string,newStringEnd); // To store possible quark containt of new string 423 423 424 #ifdef debug_LUNDfragmentation 424 #ifdef debug_LUNDfragmentation 425 G4cout<<"An attempt to determine its en 425 G4cout<<"An attempt to determine its energy (SplitEandP)"<<G4endl; 426 #endif 426 #endif 427 G4LorentzVector* HadronMomentum=SplitEa 427 G4LorentzVector* HadronMomentum=SplitEandP(HadronDefinition, string, newString); 428 428 429 delete newString; newString=0; 429 delete newString; newString=0; 430 430 431 G4KineticTrack * Hadron =0; 431 G4KineticTrack * Hadron =0; 432 if ( HadronMomentum != 0 ) { 432 if ( HadronMomentum != 0 ) { 433 433 434 #ifdef debug_LUNDfragmentation 434 #ifdef debug_LUNDfragmentation 435 G4cout<<"The attempt was successful 435 G4cout<<"The attempt was successful"<<G4endl; 436 #endif 436 #endif 437 G4ThreeVector Pos; 437 G4ThreeVector Pos; 438 Hadron = new G4KineticTrack(HadronDefinit 438 Hadron = new G4KineticTrack(HadronDefinition, 0,Pos, *HadronMomentum); 439 439 440 if ( newString ) delete newString; 440 if ( newString ) delete newString; 441 441 442 newString=new G4FragmentingString(*string 442 newString=new G4FragmentingString(*string,newStringEnd, 443 HadronMomentum); 443 HadronMomentum); 444 delete HadronMomentum; 444 delete HadronMomentum; 445 } 445 } 446 else 446 else 447 { 447 { 448 #ifdef debug_LUNDfragmentation 448 #ifdef debug_LUNDfragmentation 449 G4cout<<"The attempt was not succes 449 G4cout<<"The attempt was not successful !!!"<<G4endl; 450 #endif 450 #endif 451 } 451 } 452 452 453 #ifdef debug_LUNDfragmentation 453 #ifdef debug_LUNDfragmentation 454 G4cout<<"End SplitUP (G4VLongitudinalSt 454 G4cout<<"End SplitUP (G4VLongitudinalStringDecay) ====================="<<G4endl; 455 #endif 455 #endif 456 456 457 return Hadron; 457 return Hadron; 458 } 458 } 459 459 460 //-------------------------------------------- 460 //----------------------------------------------------------------------------- 461 461 462 G4ParticleDefinition * G4LundStringFragmentati 462 G4ParticleDefinition * G4LundStringFragmentation::DiQuarkSplitup(G4ParticleDefinition* decay, 463 463 G4ParticleDefinition *&created) 464 { 464 { 465 G4double StrSup=GetStrangeSuppress(); 465 G4double StrSup=GetStrangeSuppress(); 466 G4double ProbQQbar = (1.0 - 2.0*StrSup)*1.2 466 G4double ProbQQbar = (1.0 - 2.0*StrSup)*1.25; 467 467 468 //... can Diquark break or not? 468 //... can Diquark break or not? 469 if (G4UniformRand() < DiquarkBreakProb ){ 469 if (G4UniformRand() < DiquarkBreakProb ){ 470 470 471 //... Diquark break 471 //... Diquark break 472 G4int stableQuarkEncoding = decay->GetPD 472 G4int stableQuarkEncoding = decay->GetPDGEncoding()/1000; 473 G4int decayQuarkEncoding = (decay->GetPD 473 G4int decayQuarkEncoding = (decay->GetPDGEncoding()/100)%10; 474 if (G4UniformRand() < 0.5) 474 if (G4UniformRand() < 0.5) 475 { 475 { 476 G4int Swap = stableQuarkEncoding; 476 G4int Swap = stableQuarkEncoding; 477 stableQuarkEncoding = decayQuarkEncod 477 stableQuarkEncoding = decayQuarkEncoding; 478 decayQuarkEncoding = Swap; 478 decayQuarkEncoding = Swap; 479 } 479 } 480 480 481 G4int IsParticle=(decayQuarkEncoding>0) 481 G4int IsParticle=(decayQuarkEncoding>0) ? -1 : +1; // if we have a quark, we need antiquark 482 482 483 SetStrangenessSuppression((1.0-ProbQQbar 483 SetStrangenessSuppression((1.0-ProbQQbar)/2.0); 484 pDefPair QuarkPair = CreatePartonPair(Is 484 pDefPair QuarkPair = CreatePartonPair(IsParticle,false); // no diquarks wanted 485 SetStrangenessSuppression((1.0-StrSup)/2 485 SetStrangenessSuppression((1.0-StrSup)/2.0); 486 486 487 //... Build new Diquark 487 //... Build new Diquark 488 G4int QuarkEncoding=QuarkPair.second->Ge 488 G4int QuarkEncoding=QuarkPair.second->GetPDGEncoding(); 489 G4int i10 = std::max(std::abs(QuarkEnco 489 G4int i10 = std::max(std::abs(QuarkEncoding), std::abs(stableQuarkEncoding)); 490 G4int i20 = std::min(std::abs(QuarkEnco 490 G4int i20 = std::min(std::abs(QuarkEncoding), std::abs(stableQuarkEncoding)); 491 G4int spin = (i10 != i20 && G4UniformRan 491 G4int spin = (i10 != i20 && G4UniformRand() <= 0.5)? 1 : 3; 492 G4int NewDecayEncoding = -1*IsParticle*( 492 G4int NewDecayEncoding = -1*IsParticle*(i10 * 1000 + i20 * 100 + spin); 493 created = FindParticle(NewDecayEncoding) 493 created = FindParticle(NewDecayEncoding); 494 G4ParticleDefinition * decayQuark=FindPa 494 G4ParticleDefinition * decayQuark=FindParticle(decayQuarkEncoding); 495 G4ParticleDefinition * had=hadronizer->B 495 G4ParticleDefinition * had=hadronizer->Build(QuarkPair.first, decayQuark); 496 StrangeSuppress=StrSup; 496 StrangeSuppress=StrSup; 497 497 498 return had; 498 return had; 499 499 500 } else { 500 } else { 501 //... Diquark does not break 501 //... Diquark does not break 502 502 503 G4int IsParticle=(decay->GetPDGEncoding( 503 G4int IsParticle=(decay->GetPDGEncoding()>0) ? +1 : -1; // if we have a diquark, we need quark 504 504 505 StrangeSuppress=(1.0 - ProbQQbar)/2.0; 505 StrangeSuppress=(1.0 - ProbQQbar)/2.0; 506 pDefPair QuarkPair = CreatePartonPair(Is 506 pDefPair QuarkPair = CreatePartonPair(IsParticle,false); // no diquarks wanted 507 507 508 created = QuarkPair.second; 508 created = QuarkPair.second; 509 509 510 G4ParticleDefinition * had=hadronizer->B 510 G4ParticleDefinition * had=hadronizer->Build(QuarkPair.first, decay); 511 StrangeSuppress=StrSup; 511 StrangeSuppress=StrSup; 512 512 513 return had; 513 return had; 514 } 514 } 515 } 515 } 516 516 517 //-------------------------------------------- 517 //----------------------------------------------------------------------------- 518 518 519 G4LorentzVector * G4LundStringFragmentation::S 519 G4LorentzVector * G4LundStringFragmentation::SplitEandP(G4ParticleDefinition * pHadron, 520 520 G4FragmentingString * string, 521 521 G4FragmentingString * newString) 522 { 522 { 523 G4LorentzVector String4Momentum=string->Get4 523 G4LorentzVector String4Momentum=string->Get4Momentum(); 524 G4double StringMT2=string->MassT2(); 524 G4double StringMT2=string->MassT2(); 525 G4double StringMT =std::sqrt(StringMT2); 525 G4double StringMT =std::sqrt(StringMT2); 526 526 527 G4double HadronMass = pHadron->GetPDGMass(); 527 G4double HadronMass = pHadron->GetPDGMass(); 528 SetMinimalStringMass(newString); 528 SetMinimalStringMass(newString); 529 529 530 if ( MinimalStringMass < 0.0 ) return n 530 if ( MinimalStringMass < 0.0 ) return nullptr; 531 531 532 #ifdef debug_LUNDfragmentation 532 #ifdef debug_LUNDfragmentation 533 G4cout<<G4endl<<"Start LUND SplitEandP 533 G4cout<<G4endl<<"Start LUND SplitEandP "<<G4endl; 534 G4cout<<"String 4 mom, String M and Mt 534 G4cout<<"String 4 mom, String M and Mt "<<String4Momentum<<" "<<String4Momentum.mag() 535 <<" "<<std::sqrt(StringMT2)<<G4e 535 <<" "<<std::sqrt(StringMT2)<<G4endl; 536 G4cout<<"Hadron "<<pHadron->GetParticl 536 G4cout<<"Hadron "<<pHadron->GetParticleName()<<G4endl; 537 G4cout<<"HadM MinimalStringMassLeft St 537 G4cout<<"HadM MinimalStringMassLeft StringM hM+sM "<<HadronMass<<" "<<MinimalStringMass<<" " 538 <<String4Momentum.mag()<<" "<<Ha 538 <<String4Momentum.mag()<<" "<<HadronMass+MinimalStringMass<<G4endl; 539 #endif 539 #endif 540 540 541 if ((HadronMass + MinimalStringMass > string 541 if ((HadronMass + MinimalStringMass > string->Mass()) || MinimalStringMass < 0.) 542 { 542 { 543 #ifdef debug_LUNDfragmentation 543 #ifdef debug_LUNDfragmentation 544 G4cout<<"Mass of the string is not s 544 G4cout<<"Mass of the string is not sufficient to produce the hadron!"<<G4endl; 545 #endif 545 #endif 546 return 0; 546 return 0; 547 } // have to start all over! 547 } // have to start all over! 548 548 549 String4Momentum.setPz(0.); 549 String4Momentum.setPz(0.); 550 G4ThreeVector StringPt=String4Momentum.vect( 550 G4ThreeVector StringPt=String4Momentum.vect(); 551 StringPt.setZ(0.); 551 StringPt.setZ(0.); 552 552 553 // calculate and assign hadron transverse mo 553 // calculate and assign hadron transverse momentum component HadronPx and HadronPy 554 G4ThreeVector HadronPt , RemSysPt; 554 G4ThreeVector HadronPt , RemSysPt; 555 G4double HadronMassT2, ResidualMassT2; 555 G4double HadronMassT2, ResidualMassT2; 556 G4double HadronMt, Pt, Pt2, phi; 556 G4double HadronMt, Pt, Pt2, phi; 557 557 558 G4double TmtCur = Tmt; 558 G4double TmtCur = Tmt; 559 559 560 if ( (string->GetDecayParton()->GetPar 560 if ( (string->GetDecayParton()->GetParticleSubType()== "quark") && 561 (pHadron->GetBaryonNumber() != 0) 561 (pHadron->GetBaryonNumber() != 0) ) { 562 TmtCur = Tmt*0.37; // q 562 TmtCur = Tmt*0.37; // q->B 563 } else if ( (string->GetDecayParton()- 563 } else if ( (string->GetDecayParton()->GetParticleSubType()== "quark") && 564 (pHadron->GetBaryonNumber( 564 (pHadron->GetBaryonNumber() == 0) ) { 565 //TmtCur = Tmt; 565 //TmtCur = Tmt; // q->M 566 } else if ( (string->GetDecayParton()->GetPa 566 } else if ( (string->GetDecayParton()->GetParticleSubType()== "di_quark") && 567 (pHadron->GetBaryonNumber( 567 (pHadron->GetBaryonNumber() == 0) ) { 568 //TmtCur = Tmt*0.89; 568 //TmtCur = Tmt*0.89; // qq -> M 569 } else if ( (string->GetDecayParton()- 569 } else if ( (string->GetDecayParton()->GetParticleSubType()== "di_quark") && 570 (pHadron->GetBaryonNumber( 570 (pHadron->GetBaryonNumber() != 0) ) { 571 TmtCur = Tmt*1.35; 571 TmtCur = Tmt*1.35; // qq -> B 572 } 572 } 573 573 574 //... sample Pt of the hadron 574 //... sample Pt of the hadron 575 G4int attempt=0; 575 G4int attempt=0; 576 do 576 do 577 { 577 { 578 attempt++; if (attempt > StringLoopI 578 attempt++; if (attempt > StringLoopInterrupt) {return 0;} 579 579 580 HadronMt = HadronMass - TmtCur*G4Log 580 HadronMt = HadronMass - TmtCur*G4Log(G4UniformRand()); 581 Pt2 = sqr(HadronMt)-sqr(HadronMass); Pt=st 581 Pt2 = sqr(HadronMt)-sqr(HadronMass); Pt=std::sqrt(Pt2); 582 phi = 2.*pi*G4UniformRand(); 582 phi = 2.*pi*G4UniformRand(); 583 HadronPt = G4ThreeVector( Pt*std::co 583 HadronPt = G4ThreeVector( Pt*std::cos(phi), Pt*std::sin(phi), 0. ); 584 RemSysPt = StringPt - HadronPt; 584 RemSysPt = StringPt - HadronPt; 585 HadronMassT2 = sqr(HadronMass) + Had 585 HadronMassT2 = sqr(HadronMass) + HadronPt.mag2(); 586 ResidualMassT2=sqr(MinimalStringMass 586 ResidualMassT2=sqr(MinimalStringMass) + RemSysPt.mag2(); 587 587 588 } while (std::sqrt(HadronMassT2) + std 588 } while (std::sqrt(HadronMassT2) + std::sqrt(ResidualMassT2) > StringMT); 589 589 590 //... sample z to define hadron longitudina 590 //... sample z to define hadron longitudinal momentum and energy 591 //... but first check the available phase sp 591 //... but first check the available phase space 592 592 593 G4double Pz2 = (sqr(StringMT2 - HadronMassT2 593 G4double Pz2 = (sqr(StringMT2 - HadronMassT2 - ResidualMassT2) - 594 4*HadronMassT2 * ResidualMassT2)/4./Stri 594 4*HadronMassT2 * ResidualMassT2)/4./StringMT2; 595 595 596 if (Pz2 < 0 ) {return 0;} // have t 596 if (Pz2 < 0 ) {return 0;} // have to start all over! 597 597 598 //... then compute allowed z region z_min < 598 //... then compute allowed z region z_min <= z <= z_max 599 599 600 G4double Pz = std::sqrt(Pz2); 600 G4double Pz = std::sqrt(Pz2); 601 G4double zMin = (std::sqrt(HadronMassT2+Pz2) 601 G4double zMin = (std::sqrt(HadronMassT2+Pz2) - Pz)/std::sqrt(StringMT2); 602 // G4double zMin = (std::sqrt(HadronMa 602 // G4double zMin = (std::sqrt(HadronMassT2+Pz2) - 0.)/std::sqrt(StringMT2); // For testing purposes 603 G4double zMax = (std::sqrt(HadronMassT2+Pz2) 603 G4double zMax = (std::sqrt(HadronMassT2+Pz2) + Pz)/std::sqrt(StringMT2); 604 604 605 if (zMin >= zMax) return 0; // have to sta 605 if (zMin >= zMax) return 0; // have to start all over! 606 606 607 G4double z = GetLightConeZ(zMin, zMax, 607 G4double z = GetLightConeZ(zMin, zMax, 608 string->GetDecayParton()->Get 608 string->GetDecayParton()->GetPDGEncoding(), pHadron, 609 HadronPt.x(), HadronPt.y()); 609 HadronPt.x(), HadronPt.y()); 610 610 611 //... now compute hadron longitudinal moment 611 //... now compute hadron longitudinal momentum and energy 612 // longitudinal hadron momentum component Ha 612 // longitudinal hadron momentum component HadronPz 613 613 614 HadronPt.setZ(0.5* string->GetDecayDirection 614 HadronPt.setZ(0.5* string->GetDecayDirection() * 615 (z * string->LightConeDecay() - 615 (z * string->LightConeDecay() - 616 HadronMassT2/(z * string->LightConeD 616 HadronMassT2/(z * string->LightConeDecay()))); 617 G4double HadronE = 0.5* (z * string->LightC 617 G4double HadronE = 0.5* (z * string->LightConeDecay() + 618 HadronMassT2/(z * string->LightConeD 618 HadronMassT2/(z * string->LightConeDecay())); 619 619 620 G4LorentzVector * a4Momentum= new G4LorentzV 620 G4LorentzVector * a4Momentum= new G4LorentzVector(HadronPt,HadronE); 621 621 622 #ifdef debug_LUNDfragmentation 622 #ifdef debug_LUNDfragmentation 623 G4cout<<G4endl<<" string->GetDecayDire 623 G4cout<<G4endl<<" string->GetDecayDirection() "<<string->GetDecayDirection()<<G4endl<<G4endl; 624 G4cout<<"string->LightConeDecay() "<<s 624 G4cout<<"string->LightConeDecay() "<<string->LightConeDecay()<<G4endl; 625 G4cout<<"HadronPt,HadronE "<<HadronPt< 625 G4cout<<"HadronPt,HadronE "<<HadronPt<<" "<<HadronE<<G4endl; 626 G4cout<<"String4Momentum "<<String4Mom 626 G4cout<<"String4Momentum "<<String4Momentum<<G4endl; 627 G4cout<<"Out of LUND SplitEandP "<<G4e 627 G4cout<<"Out of LUND SplitEandP "<<G4endl<<G4endl; 628 #endif 628 #endif 629 629 630 return a4Momentum; 630 return a4Momentum; 631 } 631 } 632 632 633 //-------------------------------------------- 633 //----------------------------------------------------------------------------------------- 634 634 635 G4double G4LundStringFragmentation::GetLightCo 635 G4double G4LundStringFragmentation::GetLightConeZ(G4double zmin, G4double zmax, 636 G4int PD 636 G4int PDGEncodingOfDecayParton, 637 G4Partic 637 G4ParticleDefinition* pHadron, 638 G4double 638 G4double Px, G4double Py) 639 { 639 { 640 G4double Mass = pHadron->GetPDGMass(); 640 G4double Mass = pHadron->GetPDGMass(); 641 G4int HadronEncoding=std::abs(pHadron- 641 G4int HadronEncoding=std::abs(pHadron->GetPDGEncoding()); 642 642 643 G4double Mt2 = Px*Px + Py*Py + Mass*Mass; 643 G4double Mt2 = Px*Px + Py*Py + Mass*Mass; 644 644 645 G4double Alund, Blund; 645 G4double Alund, Blund; 646 G4double zOfMaxyf(0.), maxYf(1.), z(0.), yf( 646 G4double zOfMaxyf(0.), maxYf(1.), z(0.), yf(1.); 647 647 648 if (!((std::abs(PDGEncodingOfDecayParton) > 648 if (!((std::abs(PDGEncodingOfDecayParton) > 1000) && (HadronEncoding > 1000)) ) 649 { // ---------------- Quark fragmentation 649 { // ---------------- Quark fragmentation and qq-> meson ---------------------- 650 Alund=1.; 650 Alund=1.; 651 Blund=0.7/GeV/GeV; 651 Blund=0.7/GeV/GeV; 652 652 653 G4double BMt2 = Blund*Mt2; 653 G4double BMt2 = Blund*Mt2; 654 if (Alund == 1.0) { 654 if (Alund == 1.0) { 655 zOfMaxyf=BMt2/(Blund*Mt2 + 1.);} 655 zOfMaxyf=BMt2/(Blund*Mt2 + 1.);} 656 else { 656 else { 657 zOfMaxyf = ((1.0+BMt2) - std::sqrt(sq 657 zOfMaxyf = ((1.0+BMt2) - std::sqrt(sqr(1.0-BMt2) + 4.0*BMt2*Alund))/2.0/(1.-Alund); 658 } 658 } 659 659 660 if (zOfMaxyf < zmin) {zOfMaxyf=zmin;} 660 if (zOfMaxyf < zmin) {zOfMaxyf=zmin;} 661 if (zOfMaxyf > zmax) {zOfMaxyf=zmax;} 661 if (zOfMaxyf > zmax) {zOfMaxyf=zmax;} 662 maxYf=(1-zOfMaxyf)/zOfMaxyf * G4Exp(-Blun 662 maxYf=(1-zOfMaxyf)/zOfMaxyf * G4Exp(-Blund*Mt2/zOfMaxyf); 663 663 664 const G4int maxNumberOfLoops = 1000 664 const G4int maxNumberOfLoops = 1000; 665 G4int loopCounter = 0; 665 G4int loopCounter = 0; 666 do 666 do 667 { 667 { 668 z = zmin + G4UniformRand()*(zmax-zmin); 668 z = zmin + G4UniformRand()*(zmax-zmin); 669 //yf = (1-z)/z * G4Exp(-Blund* 669 //yf = (1-z)/z * G4Exp(-Blund*Mt2/z); 670 yf = G4Pow::GetInstance()->powA(1.0-z,Alun 670 yf = G4Pow::GetInstance()->powA(1.0-z,Alund)/z*G4Exp(-BMt2/z); 671 } 671 } 672 while ( (G4UniformRand()*maxYf > yf) && + 672 while ( (G4UniformRand()*maxYf > yf) && ++loopCounter < maxNumberOfLoops ); 673 if ( loopCounter >= maxNumberOfLoop 673 if ( loopCounter >= maxNumberOfLoops ) { 674 z = 0.5*(zmin + zmax); // Just a 674 z = 0.5*(zmin + zmax); // Just a value between zmin and zmax, no physics considerations at all! 675 } 675 } 676 return z; 676 return z; 677 } 677 } 678 678 679 if (std::abs(PDGEncodingOfDecayParton) > 100 679 if (std::abs(PDGEncodingOfDecayParton) > 1000) 680 { 680 { 681 G4double an = 2.5; 681 G4double an = 2.5; 682 an +=(sqr(Px)+sqr(Py))/sqr(GeV 682 an +=(sqr(Px)+sqr(Py))/sqr(GeV)-0.5; 683 z=zmin + (zmax-zmin)*G4Pow::Ge 683 z=zmin + (zmax-zmin)*G4Pow::GetInstance()->powA(G4UniformRand(),1./an); 684 if( PDGEncodingOfDecayParton > 684 if( PDGEncodingOfDecayParton > 3000 ) z=zmin+zmax-z; 685 } 685 } 686 686 687 return z; 687 return z; 688 } 688 } 689 689 690 //-------------------------------------------- 690 //---------------------------------------------------------------------------------------------------------- 691 691 692 G4bool G4LundStringFragmentation::SplitLast(G4 692 G4bool G4LundStringFragmentation::SplitLast(G4FragmentingString * string, 693 G4 693 G4KineticTrackVector * LeftVector, 694 G4 694 G4KineticTrackVector * RightVector) 695 { 695 { 696 //... perform last cluster decay 696 //... perform last cluster decay 697 SetMinimalStringMass( string); 697 SetMinimalStringMass( string); 698 if ( MinimalStringMass < 0.) return fa 698 if ( MinimalStringMass < 0.) return false; 699 #ifdef debug_LUNDfragmentation 699 #ifdef debug_LUNDfragmentation 700 G4cout<<G4endl<<"Split last----------- 700 G4cout<<G4endl<<"Split last-----------------------------------------"<<G4endl; 701 G4cout<<"MinimalStringMass "<<MinimalS 701 G4cout<<"MinimalStringMass "<<MinimalStringMass<<G4endl; 702 G4cout<<"Left "<<string->GetLeftParto 702 G4cout<<"Left "<<string->GetLeftParton()->GetPDGEncoding()<<" "<<string->GetPleft()<<G4endl; 703 G4cout<<"Right "<<string->GetRightPart 703 G4cout<<"Right "<<string->GetRightParton()->GetPDGEncoding()<<" "<<string->GetPright()<<G4endl; 704 G4cout<<"String4mom "<<string->GetPstr 704 G4cout<<"String4mom "<<string->GetPstring()<<" "<<string->GetPstring().mag()<<G4endl; 705 #endif 705 #endif 706 706 707 G4LorentzVector Str4Mom=string->Get4Mo 707 G4LorentzVector Str4Mom=string->Get4Momentum(); 708 G4LorentzRotation toCms(-1*Str4Mom.boo 708 G4LorentzRotation toCms(-1*Str4Mom.boostVector()); 709 G4LorentzVector Pleft = toCms * string 709 G4LorentzVector Pleft = toCms * string->GetPleft(); 710 toCms.rotateZ(-1*Pleft.phi()); 710 toCms.rotateZ(-1*Pleft.phi()); 711 toCms.rotateY(-1*Pleft.theta()); 711 toCms.rotateY(-1*Pleft.theta()); 712 712 713 G4LorentzRotation toObserverFrame= toC 713 G4LorentzRotation toObserverFrame= toCms.inverse(); 714 714 715 G4double StringMass=string->Mass(); 715 G4double StringMass=string->Mass(); 716 716 717 G4ParticleDefinition * LeftHadron(0), * Righ 717 G4ParticleDefinition * LeftHadron(0), * RightHadron(0); 718 718 719 NumberOf_FS=0; 719 NumberOf_FS=0; 720 for (G4int i=0; i<350; i++) {FS_Weight[i]=0. 720 for (G4int i=0; i<350; i++) {FS_Weight[i]=0.;} 721 721 722 G4int sampledState = 0; 722 G4int sampledState = 0; 723 723 724 #ifdef debug_LUNDfragmentation 724 #ifdef debug_LUNDfragmentation 725 G4cout<<"StrMass "<<StringMass<<" q's 725 G4cout<<"StrMass "<<StringMass<<" q's " 726 <<string->GetLeftParton()->GetPa 726 <<string->GetLeftParton()->GetParticleName()<<" " 727 <<string->GetRightParton()->GetP 727 <<string->GetRightParton()->GetParticleName()<<G4endl; 728 #endif 728 #endif 729 729 730 string->SetLeftPartonStable(); // to query q 730 string->SetLeftPartonStable(); // to query quark contents.. 731 731 732 if (string->IsAFourQuarkString() ) 732 if (string->IsAFourQuarkString() ) 733 { 733 { 734 G4int IDleft =std::abs(string->GetLe 734 G4int IDleft =std::abs(string->GetLeftParton()->GetPDGEncoding()); 735 G4int IDright=std::abs(string->GetRi 735 G4int IDright=std::abs(string->GetRightParton()->GetPDGEncoding()); 736 736 737 if ( (IDleft > 3000) || (IDright > 3 737 if ( (IDleft > 3000) || (IDright > 3000) ) { 738 if ( ! Diquark_AntiDiquark_belowTh 738 if ( ! Diquark_AntiDiquark_belowThreshold_lastSplitting(string, LeftHadron, RightHadron) ) 739 { 739 { 740 return false; 740 return false; 741 } 741 } 742 } else { 742 } else { 743 // The string is qq-qqbar type. Diquarks a 743 // The string is qq-qqbar type. Diquarks are on the string ends 744 if (StringMass-MinimalStringMass < 0 744 if (StringMass-MinimalStringMass < 0.) 745 { 745 { 746 if (! Diquark_AntiDiquark_belowThreshold 746 if (! Diquark_AntiDiquark_belowThreshold_lastSplitting(string, LeftHadron, RightHadron) ) 747 { 747 { 748 return false; 748 return false; 749 } 749 } 750 } else 750 } else 751 { 751 { 752 Diquark_AntiDiquark_aboveThreshold_lastS 752 Diquark_AntiDiquark_aboveThreshold_lastSplitting(string, LeftHadron, RightHadron); 753 753 754 if (NumberOf_FS == 0) return false; 754 if (NumberOf_FS == 0) return false; 755 755 756 sampledState = SampleS 756 sampledState = SampleState(); 757 if (string->GetLeftParton()->GetPDGEncod 757 if (string->GetLeftParton()->GetPDGEncoding() < 0) 758 { 758 { 759 LeftHadron =FS_LeftHadron[sampledState 759 LeftHadron =FS_LeftHadron[sampledState]; 760 RightHadron=FS_RightHadron[sampledStat 760 RightHadron=FS_RightHadron[sampledState]; 761 } else 761 } else 762 { 762 { 763 LeftHadron =FS_RightHadron[sampledStat 763 LeftHadron =FS_RightHadron[sampledState]; 764 RightHadron=FS_LeftHadron[sampledState 764 RightHadron=FS_LeftHadron[sampledState]; 765 } 765 } 766 } 766 } 767 } // ID > 3300 767 } // ID > 3300 768 } else { 768 } else { 769 if (string->DecayIsQuark() && string->Stab 769 if (string->DecayIsQuark() && string->StableIsQuark() ) 770 { //... there are quarks on cluster 770 { //... there are quarks on cluster ends 771 #ifdef debug_LUNDfragm 771 #ifdef debug_LUNDfragmentation 772 G4cout<<"Q Q string La 772 G4cout<<"Q Q string LastSplit"<<G4endl; 773 #endif 773 #endif 774 774 775 Quark_AntiQuark_lastSplitting(string, Le 775 Quark_AntiQuark_lastSplitting(string, LeftHadron, RightHadron); 776 776 777 if (NumberOf_FS == 0) return false; 777 if (NumberOf_FS == 0) return false; 778 sampledState = SampleState() 778 sampledState = SampleState(); 779 if (string->GetLeftParton()->GetPDGEncod 779 if (string->GetLeftParton()->GetPDGEncoding() < 0) 780 { 780 { 781 LeftHadron =FS_RightHadron[sampledStat 781 LeftHadron =FS_RightHadron[sampledState]; 782 RightHadron=FS_LeftHadron[sampledState 782 RightHadron=FS_LeftHadron[sampledState]; 783 } else 783 } else 784 { 784 { 785 LeftHadron =FS_LeftHadron[sampledState 785 LeftHadron =FS_LeftHadron[sampledState]; 786 RightHadron=FS_RightHadron[sampledStat 786 RightHadron=FS_RightHadron[sampledState]; 787 } 787 } 788 } else 788 } else 789 { //... there is a Diquark on one of 789 { //... there is a Diquark on one of the cluster ends 790 #ifdef debug_LUNDfragm 790 #ifdef debug_LUNDfragmentation 791 G4cout<<"DiQ Q string 791 G4cout<<"DiQ Q string Last Split"<<G4endl; 792 #endif 792 #endif 793 793 794 Quark_Diquark_lastSplitting(string, Left 794 Quark_Diquark_lastSplitting(string, LeftHadron, RightHadron); 795 795 796 if (NumberOf_FS == 0) return false; 796 if (NumberOf_FS == 0) return false; 797 sampledState = SampleState() 797 sampledState = SampleState(); 798 798 799 if (string->GetLeftParton()->GetParticle 799 if (string->GetLeftParton()->GetParticleSubType() == "quark") 800 { 800 { 801 LeftHadron =FS_LeftHadron[sampledState 801 LeftHadron =FS_LeftHadron[sampledState]; 802 RightHadron=FS_RightHadron[sampledStat 802 RightHadron=FS_RightHadron[sampledState]; 803 } else 803 } else 804 { 804 { 805 LeftHadron =FS_RightHadron[sampledStat 805 LeftHadron =FS_RightHadron[sampledState]; 806 RightHadron=FS_LeftHadron[sampledState 806 RightHadron=FS_LeftHadron[sampledState]; 807 } 807 } 808 } 808 } 809 809 810 } 810 } 811 811 812 #ifdef debug_LUNDfragmentation 812 #ifdef debug_LUNDfragmentation 813 G4cout<<"Sampled hadrons: "<<LeftHadro 813 G4cout<<"Sampled hadrons: "<<LeftHadron->GetParticleName()<<" "<<RightHadron->GetParticleName()<<G4endl; 814 #endif 814 #endif 815 815 816 G4LorentzVector P_left =string->GetPleft(), 816 G4LorentzVector P_left =string->GetPleft(), P_right = string->GetPright(); 817 817 818 G4LorentzVector LeftMom, RightMom; 818 G4LorentzVector LeftMom, RightMom; 819 G4ThreeVector Pos; 819 G4ThreeVector Pos; 820 820 821 Sample4Momentum(&LeftMom, LeftHadron->GetPD 821 Sample4Momentum(&LeftMom, LeftHadron->GetPDGMass(), 822 &RightMom, RightHadron->GetPDGMass(), 822 &RightMom, RightHadron->GetPDGMass(), 823 StringMass); 823 StringMass); 824 824 825 // Sample4Momentum ascribes LeftMom.pz 825 // Sample4Momentum ascribes LeftMom.pz() along positive Z axis for baryons in many cases. 826 // It must be negative in case when th 826 // It must be negative in case when the system is moving against Z axis. 827 827 828 if (!(string->DecayIsQuark() && string->Stab 828 if (!(string->DecayIsQuark() && string->StableIsQuark() )) 829 { // Only for qq - q, q - qq, and qq - qqbar 829 { // Only for qq - q, q - qq, and qq - qqbar ------------------- 830 830 831 if ( G4UniformRand() <= 0.5 ) 831 if ( G4UniformRand() <= 0.5 ) 832 { 832 { 833 if (P_left.z() <= 0.) {G4LorentzVector t 833 if (P_left.z() <= 0.) {G4LorentzVector tmp = LeftMom; LeftMom=RightMom; RightMom=tmp;} 834 } 834 } 835 else 835 else 836 { 836 { 837 if (P_right.z() >= 0.) {G4LorentzVector 837 if (P_right.z() >= 0.) {G4LorentzVector tmp = LeftMom; LeftMom=RightMom; RightMom=tmp;} 838 } 838 } 839 } 839 } 840 840 841 LeftMom *=toObserverFrame; 841 LeftMom *=toObserverFrame; 842 RightMom*=toObserverFrame; 842 RightMom*=toObserverFrame; 843 843 844 LeftVector->push_back(new G4KineticTrack(Lef 844 LeftVector->push_back(new G4KineticTrack(LeftHadron, 0, Pos, LeftMom)); 845 RightVector->push_back(new G4KineticTrack(Ri 845 RightVector->push_back(new G4KineticTrack(RightHadron, 0, Pos, RightMom)); 846 846 847 string->LorentzRotate(toObserverFrame); 847 string->LorentzRotate(toObserverFrame); 848 return true; 848 return true; 849 } 849 } 850 850 851 //-------------------------------------------- 851 //---------------------------------------------------------------------------------------- 852 852 853 G4bool G4LundStringFragmentation:: 853 G4bool G4LundStringFragmentation:: 854 Diquark_AntiDiquark_belowThreshold_lastSplitti 854 Diquark_AntiDiquark_belowThreshold_lastSplitting(G4FragmentingString * & string, 855 855 G4ParticleDefinition * & LeftHadron, 856 856 G4ParticleDefinition * & RightHadron) 857 { 857 { 858 G4double StringMass = string->Mass(); 858 G4double StringMass = string->Mass(); 859 859 860 G4int cClusterInterrupt = 0; 860 G4int cClusterInterrupt = 0; 861 G4bool isOK = false; 861 G4bool isOK = false; 862 do 862 do 863 { 863 { 864 G4int LeftQuark1= string->GetLeftParton()- 864 G4int LeftQuark1= string->GetLeftParton()->GetPDGEncoding()/1000; 865 G4int LeftQuark2=(string->GetLeftParton()- 865 G4int LeftQuark2=(string->GetLeftParton()->GetPDGEncoding()/100)%10; 866 866 867 G4int RightQuark1= string->GetRightParton( 867 G4int RightQuark1= string->GetRightParton()->GetPDGEncoding()/1000; 868 G4int RightQuark2=(string->GetRightParton( 868 G4int RightQuark2=(string->GetRightParton()->GetPDGEncoding()/100)%10; 869 869 870 if (G4UniformRand()<0.5) 870 if (G4UniformRand()<0.5) 871 { 871 { 872 LeftHadron =hadronizer->Build(FindPartic 872 LeftHadron =hadronizer->Build(FindParticle( LeftQuark1), 873 FindParticle(RightQuark1)); 873 FindParticle(RightQuark1)); 874 RightHadron= (LeftHadron == nullptr) ? n 874 RightHadron= (LeftHadron == nullptr) ? nullptr : 875 875 hadronizer->Build(FindParticle( LeftQuark2), 876 FindParticle(RightQuark2)); 876 FindParticle(RightQuark2)); 877 } else 877 } else 878 { 878 { 879 LeftHadron =hadronizer->Build(FindPartic 879 LeftHadron =hadronizer->Build(FindParticle( LeftQuark1), 880 FindParticle(RightQuark2)); 880 FindParticle(RightQuark2)); 881 RightHadron=(LeftHadron == nullptr) ? nu 881 RightHadron=(LeftHadron == nullptr) ? nullptr : 882 hadronizer 882 hadronizer->Build(FindParticle( LeftQuark2), 883 FindParticle(RightQuark1)); 883 FindParticle(RightQuark1)); 884 } 884 } 885 885 886 isOK = (LeftHadron != nullptr) && (RightHa 886 isOK = (LeftHadron != nullptr) && (RightHadron != nullptr); 887 887 888 if(isOK) { isOK = (StringMass > LeftHadron 888 if(isOK) { isOK = (StringMass > LeftHadron->GetPDGMass() + RightHadron->GetPDGMass()); } 889 ++cClusterInterrupt; 889 ++cClusterInterrupt; 890 //... repeat procedure, if mass of cluster 890 //... repeat procedure, if mass of cluster is too low to produce hadrons 891 //... ClusterMassCut = 0.15*GeV model para 891 //... ClusterMassCut = 0.15*GeV model parameter 892 } 892 } 893 while (isOK == false && cClusterInterrupt < 893 while (isOK == false && cClusterInterrupt < ClusterLoopInterrupt); 894 /* Loop checking, 07.08.2015, A.Ribon */ 894 /* Loop checking, 07.08.2015, A.Ribon */ 895 return isOK; 895 return isOK; 896 } 896 } 897 897 898 //-------------------------------------------- 898 //---------------------------------------------------------------------------------------- 899 899 900 G4bool G4LundStringFragmentation:: 900 G4bool G4LundStringFragmentation:: 901 Diquark_AntiDiquark_aboveThreshold_lastSplitti 901 Diquark_AntiDiquark_aboveThreshold_lastSplitting(G4FragmentingString * & string, 902 902 G4ParticleDefinition * & LeftHadron, 903 903 G4ParticleDefinition * & RightHadron) 904 { 904 { 905 // StringMass-MinimalStringMass > 0. Creatio 905 // StringMass-MinimalStringMass > 0. Creation of 2 baryons is possible ---- 906 906 907 G4double StringMass = string->Mass(); 907 G4double StringMass = string->Mass(); 908 G4double StringMassSqr= sqr(StringMass); 908 G4double StringMassSqr= sqr(StringMass); 909 G4ParticleDefinition * Di_Quark; 909 G4ParticleDefinition * Di_Quark; 910 G4ParticleDefinition * Anti_Di_Quark; 910 G4ParticleDefinition * Anti_Di_Quark; 911 911 912 if (string->GetLeftParton()->GetPDGEncoding( 912 if (string->GetLeftParton()->GetPDGEncoding() < 0) 913 { 913 { 914 Anti_Di_Quark =string->GetLeftParton(); 914 Anti_Di_Quark =string->GetLeftParton(); 915 Di_Quark=string->GetRightParton(); 915 Di_Quark=string->GetRightParton(); 916 } else 916 } else 917 { 917 { 918 Anti_Di_Quark =string->GetRightParton(); 918 Anti_Di_Quark =string->GetRightParton(); 919 Di_Quark=string->GetLeftParton(); 919 Di_Quark=string->GetLeftParton(); 920 } 920 } 921 921 922 G4int IDAnti_di_quark =Anti_Di_Quark->Get 922 G4int IDAnti_di_quark =Anti_Di_Quark->GetPDGEncoding(); 923 G4int AbsIDAnti_di_quark =std::abs(IDAnti_di 923 G4int AbsIDAnti_di_quark =std::abs(IDAnti_di_quark); 924 G4int IDdi_quark =Di_Quark->GetPDGEn 924 G4int IDdi_quark =Di_Quark->GetPDGEncoding(); 925 G4int AbsIDdi_quark =std::abs(IDdi_quar 925 G4int AbsIDdi_quark =std::abs(IDdi_quark); 926 926 927 G4int ADi_q1=AbsIDAnti_di_quark/1000; 927 G4int ADi_q1=AbsIDAnti_di_quark/1000; 928 G4int ADi_q2=(AbsIDAnti_di_quark-ADi_q1*1000 928 G4int ADi_q2=(AbsIDAnti_di_quark-ADi_q1*1000)/100; 929 929 930 G4int Di_q1=AbsIDdi_quark/1000; 930 G4int Di_q1=AbsIDdi_quark/1000; 931 G4int Di_q2=(AbsIDdi_quark-Di_q1*1000)/100; 931 G4int Di_q2=(AbsIDdi_quark-Di_q1*1000)/100; 932 932 933 NumberOf_FS=0; 933 NumberOf_FS=0; 934 for (G4int ProdQ=1; ProdQ < 6; ProdQ++) 934 for (G4int ProdQ=1; ProdQ < 6; ProdQ++) 935 { 935 { 936 G4int StateADiQ=0; 936 G4int StateADiQ=0; 937 const G4int maxNumberOfLoops = 937 const G4int maxNumberOfLoops = 1000; 938 G4int loopCounter = 0; 938 G4int loopCounter = 0; 939 do // while(Meson[AbsIDquark-1][ProdQ-1][ 939 do // while(Meson[AbsIDquark-1][ProdQ-1][StateQ]<>0); 940 { 940 { 941 LeftHadron=G4ParticleTable::GetParticleT 941 LeftHadron=G4ParticleTable::GetParticleTable()->FindParticle( 942 -Baryon[ADi_q1-1][ADi_q2-1][Prod 942 -Baryon[ADi_q1-1][ADi_q2-1][ProdQ-1][StateADiQ]); 943 943 944 if (LeftHadron == NULL) continue; 944 if (LeftHadron == NULL) continue; 945 G4double LeftHadronMass=LeftHadron->GetP 945 G4double LeftHadronMass=LeftHadron->GetPDGMass(); 946 946 947 G4int StateDiQ=0; 947 G4int StateDiQ=0; 948 const G4int maxNumberO 948 const G4int maxNumberOfInternalLoops = 1000; 949 G4int internalLoopCoun 949 G4int internalLoopCounter = 0; 950 do // while(Baryon[Di_q1-1][Di_q2-1][Pro 950 do // while(Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]<>0); 951 { 951 { 952 RightHadron=G4ParticleTable::GetPartic 952 RightHadron=G4ParticleTable::GetParticleTable()->FindParticle( 953 +Baryon[Di_q1-1][Di_q2-1][Prod 953 +Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]); 954 954 955 if (RightHadron == NULL) continue; 955 if (RightHadron == NULL) continue; 956 G4double RightHadronMass=RightHadron-> 956 G4double RightHadronMass=RightHadron->GetPDGMass(); 957 957 958 if (StringMass > LeftHadronMass + Righ 958 if (StringMass > LeftHadronMass + RightHadronMass) 959 { 959 { 960 if ( N 960 if ( NumberOf_FS > 349 ) { 961 G4Ex 961 G4ExceptionDescription ed; 962 ed < 962 ed << " NumberOf_FS exceeds its limit: NumberOf_FS=" << NumberOf_FS << G4endl; 963 G4Ex 963 G4Exception( "G4LundStringFragmentation::Diquark_AntiDiquark_aboveThreshold_lastSplitting ", 964 964 "HAD_LUND_001", JustWarning, ed ); 965 Numb 965 NumberOf_FS = 349; 966 } 966 } 967 967 968 G4double FS_Psqr=lambda(StringMassSq 968 G4double FS_Psqr=lambda(StringMassSqr,sqr(LeftHadronMass), 969 sqr(RightHadronMass)); 969 sqr(RightHadronMass)); 970 //FS_Psqr=1.; 970 //FS_Psqr=1.; 971 FS_Weight[NumberOf_FS]=std::sqrt(FS_ 971 FS_Weight[NumberOf_FS]=std::sqrt(FS_Psqr)*FS_Psqr* 972 BaryonWeight[ADi_q1-1][AD 972 BaryonWeight[ADi_q1-1][ADi_q2-1][ProdQ-1][StateADiQ]* 973 BaryonWeight[Di_q1-1][Di_ 973 BaryonWeight[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]* 974 Prob_QQbar[ProdQ-1]; 974 Prob_QQbar[ProdQ-1]; 975 975 976 FS_LeftHadron[NumberOf_FS] = LeftHad 976 FS_LeftHadron[NumberOf_FS] = LeftHadron; 977 FS_RightHadron[NumberOf_FS]= RightHa 977 FS_RightHadron[NumberOf_FS]= RightHadron; 978 NumberOf_FS++; 978 NumberOf_FS++; 979 } // End of if (StringMass > LeftHadro 979 } // End of if (StringMass > LeftHadronMass + RightHadronMass) 980 980 981 StateDiQ++; 981 StateDiQ++; 982 982 983 } while( (Baryon[Di_q1-1][Di_q2-1][ProdQ 983 } while( (Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]!=0) && 984 ++internalLoo 984 ++internalLoopCounter < maxNumberOfInternalLoops ); 985 if ( internalLoopCount 985 if ( internalLoopCounter >= maxNumberOfInternalLoops ) { 986 return false; 986 return false; 987 } 987 } 988 988 989 StateADiQ++; 989 StateADiQ++; 990 } while( (Baryon[ADi_q1-1][ADi_q2-1][ProdQ 990 } while( (Baryon[ADi_q1-1][ADi_q2-1][ProdQ-1][StateADiQ]!=0) && 991 ++loopCounter < maxNu 991 ++loopCounter < maxNumberOfLoops ); 992 if ( loopCounter >= maxNumberO 992 if ( loopCounter >= maxNumberOfLoops ) { 993 return false; 993 return false; 994 } 994 } 995 } // End of for (G4int ProdQ=1; ProdQ < 4; P 995 } // End of for (G4int ProdQ=1; ProdQ < 4; ProdQ++) 996 996 997 return true; 997 return true; 998 } 998 } 999 999 1000 //------------------------------------------- 1000 //---------------------------------------------------------------------------------------- 1001 1001 1002 G4bool G4LundStringFragmentation::Quark_Diqua 1002 G4bool G4LundStringFragmentation::Quark_Diquark_lastSplitting(G4FragmentingString * & string, 1003 1003 G4ParticleDefinition * & LeftHadron, 1004 1004 G4ParticleDefinition * & RightHadron) 1005 { 1005 { 1006 G4double StringMass = string->Mass(); 1006 G4double StringMass = string->Mass(); 1007 G4double StringMassSqr= sqr(StringMass); 1007 G4double StringMassSqr= sqr(StringMass); 1008 1008 1009 G4ParticleDefinition * Di_Quark; 1009 G4ParticleDefinition * Di_Quark; 1010 G4ParticleDefinition * Quark; 1010 G4ParticleDefinition * Quark; 1011 1011 1012 if (string->GetLeftParton()->GetParticleSub 1012 if (string->GetLeftParton()->GetParticleSubType()== "quark") 1013 { 1013 { 1014 Quark =string->GetLeftParton(); 1014 Quark =string->GetLeftParton(); 1015 Di_Quark=string->GetRightParton(); 1015 Di_Quark=string->GetRightParton(); 1016 } else 1016 } else 1017 { 1017 { 1018 Quark =string->GetRightParton(); 1018 Quark =string->GetRightParton(); 1019 Di_Quark=string->GetLeftParton(); 1019 Di_Quark=string->GetLeftParton(); 1020 } 1020 } 1021 1021 1022 G4int IDquark =Quark->GetPDGEncoding 1022 G4int IDquark =Quark->GetPDGEncoding(); 1023 G4int AbsIDquark =std::abs(IDquark); 1023 G4int AbsIDquark =std::abs(IDquark); 1024 G4int IDdi_quark =Di_Quark->GetPDGEncodin 1024 G4int IDdi_quark =Di_Quark->GetPDGEncoding(); 1025 G4int AbsIDdi_quark=std::abs(IDdi_quark); 1025 G4int AbsIDdi_quark=std::abs(IDdi_quark); 1026 G4int Di_q1=AbsIDdi_quark/1000; 1026 G4int Di_q1=AbsIDdi_quark/1000; 1027 G4int Di_q2=(AbsIDdi_quark-Di_q1*1000)/100; 1027 G4int Di_q2=(AbsIDdi_quark-Di_q1*1000)/100; 1028 G4int SignDiQ= 1; 1028 G4int SignDiQ= 1; 1029 if (IDdi_quark < 0) SignDiQ=-1; 1029 if (IDdi_quark < 0) SignDiQ=-1; 1030 1030 1031 NumberOf_FS=0; 1031 NumberOf_FS=0; 1032 for (G4int ProdQ=1; ProdQ < 4; ProdQ++) // 1032 for (G4int ProdQ=1; ProdQ < 4; ProdQ++) // Loop over quark-antiquark cases: u-ubar, d-dbar, s-sbar 1033 { // 1033 { // (as last splitting, do not consider c-cbar and b-bbar cases) 1034 G4int SignQ; 1034 G4int SignQ; 1035 if (IDquark > 0) 1035 if (IDquark > 0) 1036 { 1036 { 1037 SignQ=-1; 1037 SignQ=-1; 1038 if (IDquark == 2) 1038 if (IDquark == 2) SignQ= 1; 1039 if ((IDquark == 1) && (ProdQ == 3)) Sig 1039 if ((IDquark == 1) && (ProdQ == 3)) SignQ= 1; // K0 1040 if ((IDquark == 3) && (ProdQ == 1)) Sig 1040 if ((IDquark == 3) && (ProdQ == 1)) SignQ=-1; // K0bar 1041 if (IDquark == 4) 1041 if (IDquark == 4) SignQ= 1; // D+, D0, Ds+ 1042 if (IDquark == 5) Sig 1042 if (IDquark == 5) SignQ=-1; // B-, anti_B0, anti_Bs0 1043 } else 1043 } else 1044 { 1044 { 1045 SignQ= 1; 1045 SignQ= 1; 1046 if (IDquark == -2) Sig 1046 if (IDquark == -2) SignQ=-1; 1047 if ((IDquark ==-1) && (ProdQ == 3)) Sig 1047 if ((IDquark ==-1) && (ProdQ == 3)) SignQ=-1; // K0bar 1048 if ((IDquark ==-3) && (ProdQ == 1)) Sig 1048 if ((IDquark ==-3) && (ProdQ == 1)) SignQ= 1; // K0 1049 if (IDquark == -4) Sig 1049 if (IDquark == -4) SignQ=-1; // D-, anti_D0, anti_Ds+ 1050 if (IDquark == -5) Sig 1050 if (IDquark == -5) SignQ= 1; // B+, B0, Bs0 1051 } 1051 } 1052 1052 1053 if (AbsIDquark == ProdQ) SignQ 1053 if (AbsIDquark == ProdQ) SignQ= 1; 1054 1054 1055 G4int StateQ=0; 1055 G4int StateQ=0; 1056 const G4int maxNumberOfLoops 1056 const G4int maxNumberOfLoops = 1000; 1057 G4int loopCounter = 0; 1057 G4int loopCounter = 0; 1058 do // while(Meson[AbsIDquark-1][ProdQ-1] 1058 do // while(Meson[AbsIDquark-1][ProdQ-1][StateQ]<>0); 1059 { 1059 { 1060 LeftHadron=G4ParticleTable::GetParticle 1060 LeftHadron=G4ParticleTable::GetParticleTable()->FindParticle(SignQ* 1061 Meson[AbsIDquark-1][ProdQ-1][St 1061 Meson[AbsIDquark-1][ProdQ-1][StateQ]); 1062 if (LeftHadron == NULL) continue; 1062 if (LeftHadron == NULL) continue; 1063 G4double LeftHadronMass=LeftHadron->Get 1063 G4double LeftHadronMass=LeftHadron->GetPDGMass(); 1064 1064 1065 G4int StateDiQ=0; 1065 G4int StateDiQ=0; 1066 const G4int maxNumber 1066 const G4int maxNumberOfInternalLoops = 1000; 1067 G4int internalLoopCou 1067 G4int internalLoopCounter = 0; 1068 do // while(Baryon[Di_q1-1][Di_q2-1][Pr 1068 do // while(Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]<>0); 1069 { 1069 { 1070 RightHadron=G4ParticleTable::GetParti 1070 RightHadron=G4ParticleTable::GetParticleTable()->FindParticle(SignDiQ* 1071 Baryon[Di_q1-1][Di_q2-1][Prod 1071 Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]); 1072 if (RightHadron == NULL) continue; 1072 if (RightHadron == NULL) continue; 1073 G4double RightHadronMass=RightHadron- 1073 G4double RightHadronMass=RightHadron->GetPDGMass(); 1074 1074 1075 if (StringMass > LeftHadronMass + Rig 1075 if (StringMass > LeftHadronMass + RightHadronMass) 1076 { 1076 { 1077 if ( 1077 if ( NumberOf_FS > 349 ) { 1078 G4E 1078 G4ExceptionDescription ed; 1079 ed 1079 ed << " NumberOf_FS exceeds its limit: NumberOf_FS=" << NumberOf_FS << G4endl; 1080 G4E 1080 G4Exception( "G4LundStringFragmentation::Quark_Diquark_lastSplitting ", 1081 1081 "HAD_LUND_002", JustWarning, ed ); 1082 Num 1082 NumberOf_FS = 349; 1083 } 1083 } 1084 1084 1085 G4double FS_Psqr=lambda(StringMassS 1085 G4double FS_Psqr=lambda(StringMassSqr,sqr(LeftHadronMass), 1086 sqr(RightHadronMass)); 1086 sqr(RightHadronMass)); 1087 FS_Weight[NumberOf_FS]=std::sqrt(FS 1087 FS_Weight[NumberOf_FS]=std::sqrt(FS_Psqr)* 1088 MesonWeight[AbsIDquark-1 1088 MesonWeight[AbsIDquark-1][ProdQ-1][StateQ]* 1089 BaryonWeight[Di_q1-1][Di 1089 BaryonWeight[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]* 1090 Prob_QQbar[ProdQ-1]; 1090 Prob_QQbar[ProdQ-1]; 1091 1091 1092 FS_LeftHadron[NumberOf_FS] = LeftHa 1092 FS_LeftHadron[NumberOf_FS] = LeftHadron; 1093 FS_RightHadron[NumberOf_FS]= RightH 1093 FS_RightHadron[NumberOf_FS]= RightHadron; 1094 1094 1095 NumberOf_FS++; 1095 NumberOf_FS++; 1096 } // End of if (StringMass > LeftHadr 1096 } // End of if (StringMass > LeftHadronMass + RightHadronMass) 1097 1097 1098 StateDiQ++; 1098 StateDiQ++; 1099 1099 1100 } while( (Baryon[Di_q1-1][Di_q2-1][Prod 1100 } while( (Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]!=0) && 1101 ++internalLo 1101 ++internalLoopCounter < maxNumberOfInternalLoops ); 1102 if ( internalLoopCoun 1102 if ( internalLoopCounter >= maxNumberOfInternalLoops ) { 1103 return false; 1103 return false; 1104 } 1104 } 1105 1105 1106 StateQ++; 1106 StateQ++; 1107 } while( (Meson[AbsIDquark-1][ProdQ-1][St 1107 } while( (Meson[AbsIDquark-1][ProdQ-1][StateQ]!=0) && 1108 ++loopCounter < maxN 1108 ++loopCounter < maxNumberOfLoops ); /* Loop checking, 07.08.2015, A.Ribon */ 1109 1109 1110 if ( loopCounter >= maxNumb 1110 if ( loopCounter >= maxNumberOfLoops ) { 1111 return false; 1111 return false; 1112 } 1112 } 1113 } 1113 } 1114 1114 1115 return true; 1115 return true; 1116 } 1116 } 1117 1117 1118 //------------------------------------------- 1118 //---------------------------------------------------------------------------------------- 1119 1119 1120 G4bool G4LundStringFragmentation::Quark_AntiQ 1120 G4bool G4LundStringFragmentation::Quark_AntiQuark_lastSplitting(G4FragmentingString * & string, 1121 1121 G4ParticleDefinition * & LeftHadron, 1122 1122 G4ParticleDefinition * & RightHadron) 1123 { 1123 { 1124 G4double StringMass = string->Mass(); 1124 G4double StringMass = string->Mass(); 1125 G4double StringMassSqr= sqr(StringMass); 1125 G4double StringMassSqr= sqr(StringMass); 1126 1126 1127 G4ParticleDefinition * Quark; 1127 G4ParticleDefinition * Quark; 1128 G4ParticleDefinition * Anti_Quark; 1128 G4ParticleDefinition * Anti_Quark; 1129 1129 1130 if (string->GetLeftParton()->GetPDGEncoding 1130 if (string->GetLeftParton()->GetPDGEncoding()>0) 1131 { 1131 { 1132 Quark =string->GetLeftParton(); 1132 Quark =string->GetLeftParton(); 1133 Anti_Quark=string->GetRightParton(); 1133 Anti_Quark=string->GetRightParton(); 1134 } else 1134 } else 1135 { 1135 { 1136 Quark =string->GetRightParton(); 1136 Quark =string->GetRightParton(); 1137 Anti_Quark=string->GetLeftParton(); 1137 Anti_Quark=string->GetLeftParton(); 1138 } 1138 } 1139 1139 1140 G4int IDquark =Quark->GetPDGEncodin 1140 G4int IDquark =Quark->GetPDGEncoding(); 1141 G4int AbsIDquark =std::abs(IDquark); 1141 G4int AbsIDquark =std::abs(IDquark); 1142 G4int QuarkCharge =Qcharge[IDquar 1142 G4int QuarkCharge =Qcharge[IDquark-1]; 1143 1143 1144 G4int IDanti_quark =Anti_Quark->GetPDGEn 1144 G4int IDanti_quark =Anti_Quark->GetPDGEncoding(); 1145 G4int AbsIDanti_quark =std::abs(IDanti_quar 1145 G4int AbsIDanti_quark =std::abs(IDanti_quark); 1146 G4int AntiQuarkCharge =-Qcharge[AbsID 1146 G4int AntiQuarkCharge =-Qcharge[AbsIDanti_quark-1]; 1147 1147 1148 G4int LeftHadronCharge(0), RightHadro 1148 G4int LeftHadronCharge(0), RightHadronCharge(0); 1149 1149 1150 //G4cout<<"Q Qbar "<<IDquark<<" "<<ID 1150 //G4cout<<"Q Qbar "<<IDquark<<" "<<IDanti_quark<<G4endl; 1151 1151 1152 NumberOf_FS=0; 1152 NumberOf_FS=0; 1153 for (G4int ProdQ=1; ProdQ < 4; ProdQ++) // 1153 for (G4int ProdQ=1; ProdQ < 4; ProdQ++) // Loop over quark-antiquark cases: u-ubar, d-dbar, s-sbar 1154 { // 1154 { // (as last splitting, do not consider c-cbar and b-bbar cases) 1155 //G4cout<<"NumberOf_FS ProdQ 1155 //G4cout<<"NumberOf_FS ProdQ "<<NumberOf_FS<<" "<<ProdQ<<G4endl; 1156 LeftHadronCharge = QuarkCharge - Qcharge[ 1156 LeftHadronCharge = QuarkCharge - Qcharge[ProdQ-1]; 1157 G4int SignQ = LeftHadronCharge/3; if (Sig 1157 G4int SignQ = LeftHadronCharge/3; if (SignQ == 0) SignQ = 1; 1158 1158 1159 if ((IDquark == 1) && (ProdQ == 3)) SignQ 1159 if ((IDquark == 1) && (ProdQ == 3)) SignQ= 1; // K0 (d,sbar) 1160 if ((IDquark == 3) && (ProdQ == 1)) SignQ 1160 if ((IDquark == 3) && (ProdQ == 1)) SignQ=-1; // K0bar (s,dbar) 1161 if ((IDquark == 4) && (ProdQ 1161 if ((IDquark == 4) && (ProdQ == 2)) SignQ= 1; // D0 (c,ubar) 1162 if ((IDquark == 5) && (ProdQ 1162 if ((IDquark == 5) && (ProdQ == 1)) SignQ=-1; // anti_B0 (b,dbar) 1163 if ((IDquark == 5) && (ProdQ 1163 if ((IDquark == 5) && (ProdQ == 3)) SignQ=-1; // anti_Bs0 (b,sbar) 1164 1164 1165 RightHadronCharge = AntiQuark 1165 RightHadronCharge = AntiQuarkCharge + Qcharge[ProdQ-1]; 1166 G4int SignAQ = RightHadronCharge/3; if (S 1166 G4int SignAQ = RightHadronCharge/3; if (SignAQ == 0) SignAQ = 1; 1167 1167 1168 if ((IDanti_quark ==-1) && (ProdQ == 3)) 1168 if ((IDanti_quark ==-1) && (ProdQ == 3)) SignAQ=-1; // K0bar (dbar,s) 1169 if ((IDanti_quark ==-3) && (ProdQ == 1)) 1169 if ((IDanti_quark ==-3) && (ProdQ == 1)) SignAQ= 1; // K0 (sbar,d) 1170 if ((IDanti_quark ==-4) && (ProdQ == 2)) 1170 if ((IDanti_quark ==-4) && (ProdQ == 2)) SignAQ=-1; // anti_D0 (cbar,u) 1171 if ((IDanti_quark ==-5) && (ProdQ == 1)) 1171 if ((IDanti_quark ==-5) && (ProdQ == 1)) SignAQ= 1; // B0 (bbar,d) 1172 if ((IDanti_quark ==-5) && (ProdQ == 3)) 1172 if ((IDanti_quark ==-5) && (ProdQ == 3)) SignAQ= 1; // Bs0 (bbar,s) 1173 1173 1174 //G4cout<<"ProQ signs "<<Prod 1174 //G4cout<<"ProQ signs "<<ProdQ<<" "<<SignQ<<" "<<SignAQ<<G4endl; 1175 1175 1176 G4int StateQ=0; 1176 G4int StateQ=0; 1177 const G4int maxNumberOfLoops 1177 const G4int maxNumberOfLoops = 1000; 1178 G4int loopCounter = 0; 1178 G4int loopCounter = 0; 1179 do 1179 do 1180 { 1180 { 1181 //G4cout<<"[AbsIDquar 1181 //G4cout<<"[AbsIDquark-1][ProdQ-1][StateQ "<<AbsIDquark-1<<" " 1182 //<<ProdQ-1<<" "<<Sta 1182 //<<ProdQ-1<<" "<<StateQ<<" "<<SignQ*Meson[AbsIDquark-1][ProdQ-1][StateQ]<<G4endl; 1183 LeftHadron=G4ParticleTable::GetParticle 1183 LeftHadron=G4ParticleTable::GetParticleTable()->FindParticle(SignQ* 1184 Meson[AbsIDquark-1][ProdQ- 1184 Meson[AbsIDquark-1][ProdQ-1][StateQ]); 1185 //G4cout<<"LeftHadron 1185 //G4cout<<"LeftHadron "<<LeftHadron<<G4endl; 1186 if (LeftHadron == NULL) { StateQ++; con 1186 if (LeftHadron == NULL) { StateQ++; continue; } 1187 //G4cout<<"LeftHadron 1187 //G4cout<<"LeftHadron "<<LeftHadron->GetParticleName()<<G4endl; 1188 G4double LeftHadronMass=LeftHadron->Get 1188 G4double LeftHadronMass=LeftHadron->GetPDGMass(); 1189 1189 1190 G4int StateAQ=0; 1190 G4int StateAQ=0; 1191 const G4int maxNumber 1191 const G4int maxNumberOfInternalLoops = 1000; 1192 G4int internalLoopCou 1192 G4int internalLoopCounter = 0; 1193 do 1193 do 1194 { 1194 { 1195 //G4cout<<" 1195 //G4cout<<" [AbsIDanti_quark-1][ProdQ-1][StateAQ] "<<AbsIDanti_quark-1<<" " 1196 // <<Pro 1196 // <<ProdQ-1<<" "<<StateAQ<<" "<<SignAQ*Meson[AbsIDanti_quark-1][ProdQ-1][StateAQ]<<G4endl; 1197 RightHadron=G4ParticleTable::GetParti 1197 RightHadron=G4ParticleTable::GetParticleTable()->FindParticle(SignAQ* 1198 Meson[AbsIDanti_quark-1][Prod 1198 Meson[AbsIDanti_quark-1][ProdQ-1][StateAQ]); 1199 //G4cout<<"Ri 1199 //G4cout<<"RightHadron "<<RightHadron<<G4endl; 1200 if(RightHadron == NULL) { StateAQ++; 1200 if(RightHadron == NULL) { StateAQ++; continue; } 1201 //G4cout<<"Ri 1201 //G4cout<<"RightHadron "<<RightHadron->GetParticleName()<<G4endl; 1202 G4double RightHadronMass=RightHadron- 1202 G4double RightHadronMass=RightHadron->GetPDGMass(); 1203 1203 1204 if (StringMass > LeftHadronMass + Rig 1204 if (StringMass > LeftHadronMass + RightHadronMass) 1205 { 1205 { 1206 if ( 1206 if ( NumberOf_FS > 349 ) { 1207 G4E 1207 G4ExceptionDescription ed; 1208 ed 1208 ed << " NumberOf_FS exceeds its limit: NumberOf_FS=" << NumberOf_FS << G4endl; 1209 G4E 1209 G4Exception( "G4LundStringFragmentation::Quark_AntiQuark_lastSplitting ", 1210 1210 "HAD_LUND_003", JustWarning, ed ); 1211 Num 1211 NumberOf_FS = 349; 1212 } 1212 } 1213 1213 1214 G4dou 1214 G4double FS_Psqr=lambda(StringMassSqr,sqr(LeftHadronMass), 1215 sqr(RightHadronMass)); 1215 sqr(RightHadronMass)); 1216 //FS_Psqr=1.; 1216 //FS_Psqr=1.; 1217 FS_Weight[NumberOf_FS]=std::sqrt(FS 1217 FS_Weight[NumberOf_FS]=std::sqrt(FS_Psqr)* 1218 MesonWeight[AbsIDquark-1 1218 MesonWeight[AbsIDquark-1][ProdQ-1][StateQ]* 1219 MesonWeight[AbsIDanti_qu 1219 MesonWeight[AbsIDanti_quark-1][ProdQ-1][StateAQ]* 1220 Prob_QQbar[ProdQ-1]; 1220 Prob_QQbar[ProdQ-1]; 1221 if (string->GetLeftParton()->GetPDG 1221 if (string->GetLeftParton()->GetPDGEncoding()>0) 1222 { 1222 { 1223 FS_LeftHadron[NumberOf_FS] = Righ 1223 FS_LeftHadron[NumberOf_FS] = RightHadron; 1224 FS_RightHadron[NumberOf_FS]= Left 1224 FS_RightHadron[NumberOf_FS]= LeftHadron; 1225 } else 1225 } else 1226 { 1226 { 1227 FS_LeftHadron[NumberOf_FS] = Left 1227 FS_LeftHadron[NumberOf_FS] = LeftHadron; 1228 FS_RightHadron[NumberOf_FS]= Righ 1228 FS_RightHadron[NumberOf_FS]= RightHadron; 1229 } 1229 } 1230 1230 1231 NumberOf_FS++; 1231 NumberOf_FS++; 1232 1232 1233 } 1233 } 1234 1234 1235 StateAQ++; 1235 StateAQ++; 1236 //G4cout<<" 1236 //G4cout<<" StateAQ Meson[AbsIDanti_quark-1][ProdQ-1][StateAQ] "<<StateAQ<<" " 1237 // <<Mes 1237 // <<Meson[AbsIDanti_quark-1][ProdQ-1][StateAQ]<<" "<<internalLoopCounter<<G4endl; 1238 } while ( (Meson[AbsIDanti_quark-1][Pro 1238 } while ( (Meson[AbsIDanti_quark-1][ProdQ-1][StateAQ]!=0) && 1239 ++internalL 1239 ++internalLoopCounter < maxNumberOfInternalLoops ); 1240 if ( internalLoopCo 1240 if ( internalLoopCounter >= maxNumberOfInternalLoops ) { 1241 return false; 1241 return false; 1242 } 1242 } 1243 1243 1244 StateQ++; 1244 StateQ++; 1245 //G4cout<<"StateQ Mes 1245 //G4cout<<"StateQ Meson[AbsIDquark-1][ProdQ-1][StateQ] "<<StateQ<<" " 1246 // <<Meson[AbsID 1246 // <<Meson[AbsIDquark-1][ProdQ-1][StateQ]<<" "<<loopCounter<<G4endl; 1247 1247 1248 } while ( (Meson[AbsIDquark-1][ProdQ-1][S 1248 } while ( (Meson[AbsIDquark-1][ProdQ-1][StateQ]!=0) && 1249 ++loopCounter < maxN 1249 ++loopCounter < maxNumberOfLoops ); 1250 if ( loopCounter >= maxNumb 1250 if ( loopCounter >= maxNumberOfLoops ) { 1251 return false; 1251 return false; 1252 } 1252 } 1253 } // End of for (G4int ProdQ=1; ProdQ < 4; 1253 } // End of for (G4int ProdQ=1; ProdQ < 4; ProdQ++) 1254 1254 1255 return true; 1255 return true; 1256 } 1256 } 1257 1257 1258 //------------------------------------------- 1258 //---------------------------------------------------------------------------------------------------------- 1259 1259 1260 G4int G4LundStringFragmentation::SampleState( 1260 G4int G4LundStringFragmentation::SampleState(void) 1261 { 1261 { 1262 if ( NumberOf_FS > 349 ) { 1262 if ( NumberOf_FS > 349 ) { 1263 G4ExceptionDescription ed; 1263 G4ExceptionDescription ed; 1264 ed << " NumberOf_FS exceeds its lim 1264 ed << " NumberOf_FS exceeds its limit: NumberOf_FS=" << NumberOf_FS << G4endl; 1265 G4Exception( "G4LundStringFragmenta 1265 G4Exception( "G4LundStringFragmentation::SampleState ", "HAD_LUND_004", JustWarning, ed ); 1266 NumberOf_FS = 349; 1266 NumberOf_FS = 349; 1267 } 1267 } 1268 1268 1269 G4double SumWeights=0.; 1269 G4double SumWeights=0.; 1270 for (G4int i=0; i<NumberOf_FS; i++) {SumWei 1270 for (G4int i=0; i<NumberOf_FS; i++) {SumWeights+=FS_Weight[i];} 1271 1271 1272 G4double ksi=G4UniformRand(); 1272 G4double ksi=G4UniformRand(); 1273 G4double Sum=0.; 1273 G4double Sum=0.; 1274 G4int indexPosition = 0; 1274 G4int indexPosition = 0; 1275 1275 1276 for (G4int i=0; i<NumberOf_FS; i++) 1276 for (G4int i=0; i<NumberOf_FS; i++) 1277 { 1277 { 1278 Sum+=(FS_Weight[i]/SumWeights); 1278 Sum+=(FS_Weight[i]/SumWeights); 1279 indexPosition=i; 1279 indexPosition=i; 1280 if (Sum >= ksi) break; 1280 if (Sum >= ksi) break; 1281 } 1281 } 1282 return indexPosition; 1282 return indexPosition; 1283 } 1283 } 1284 1284 1285 //------------------------------------------- 1285 //---------------------------------------------------------------------------------------------------------- 1286 1286 1287 void G4LundStringFragmentation::Sample4Moment 1287 void G4LundStringFragmentation::Sample4Momentum(G4LorentzVector* Mom, G4double Mass, 1288 1288 G4LorentzVector* AntiMom, G4double AntiMass, 1289 1289 G4double InitialMass) 1290 { 1290 { 1291 // ------ Sampling of momenta of 2 last pro 1291 // ------ Sampling of momenta of 2 last produced hadrons -------------------- 1292 G4ThreeVector Pt; 1292 G4ThreeVector Pt; 1293 G4double MassMt, AntiMassMt; 1293 G4double MassMt, AntiMassMt; 1294 G4double AvailablePz, AvailablePz2; 1294 G4double AvailablePz, AvailablePz2; 1295 #ifdef debug_LUNDfragmentation 1295 #ifdef debug_LUNDfragmentation 1296 G4cout<<"Sampling of momenta of 2 las 1296 G4cout<<"Sampling of momenta of 2 last produced hadrons ----------------"<<G4endl; 1297 G4cout<<"Init Mass "<<InitialMass<<" 1297 G4cout<<"Init Mass "<<InitialMass<<" FirstM "<<Mass<<" SecondM "<<AntiMass<<" ProbIsotropy "<<G4endl; 1298 #endif 1298 #endif 1299 1299 1300 G4double r_val = sqr(InitialMass*InitialMas 1300 G4double r_val = sqr(InitialMass*InitialMass - Mass*Mass - AntiMass*AntiMass) - 1301 sqr(2.*Mass*AntiMass); 1301 sqr(2.*Mass*AntiMass); 1302 G4double Pabs = (r_val > 0.)? std::sqrt(r_v 1302 G4double Pabs = (r_val > 0.)? std::sqrt(r_val)/(2.*InitialMass) : 0; 1303 1303 1304 const G4int maxNumberOfLoops = 1000; 1304 const G4int maxNumberOfLoops = 1000; 1305 G4double SigmaQTw = SigmaQT; 1305 G4double SigmaQTw = SigmaQT; 1306 if ( Mass > 930. || AntiMass > 930. ) { 1306 if ( Mass > 930. || AntiMass > 930. ) { 1307 SigmaQT *= ( 1.0 - 0.55*sqr( (Mass+AntiMa 1307 SigmaQT *= ( 1.0 - 0.55*sqr( (Mass+AntiMass)/InitialMass ) ); 1308 } 1308 } 1309 if ( Mass < 930. && AntiMass < 930. ) 1309 if ( Mass < 930. && AntiMass < 930. ) {} // q-qbar string 1310 if ( ( Mass < 930. && AntiMass > 930. 1310 if ( ( Mass < 930. && AntiMass > 930. ) || 1311 ( Mass > 930. && AntiMass < 930. ) ) { 1311 ( Mass > 930. && AntiMass < 930. ) ) { // q-di_q string 1312 //SigmaQT = -1.; // isotropical de 1312 //SigmaQT = -1.; // isotropical decay 1313 } 1313 } 1314 if ( Mass > 930. && AntiMass > 930. ) 1314 if ( Mass > 930. && AntiMass > 930. ) { // qq-qqbar string 1315 SigmaQT *= ( 1.0 - 0.55*sqr( (Mass+ 1315 SigmaQT *= ( 1.0 - 0.55*sqr( (Mass+AntiMass)/InitialMass ) ); 1316 } 1316 } 1317 1317 1318 G4int loopCounter = 0; 1318 G4int loopCounter = 0; 1319 do 1319 do 1320 { 1320 { 1321 Pt=SampleQuarkPt(Pabs); Pt.setZ(0); G4dou 1321 Pt=SampleQuarkPt(Pabs); Pt.setZ(0); G4double Pt2=Pt.mag2(); 1322 MassMt = std::sqrt( Mass * Mass 1322 MassMt = std::sqrt( Mass * Mass + Pt2); 1323 AntiMassMt= std::sqrt(AntiMass * AntiMass 1323 AntiMassMt= std::sqrt(AntiMass * AntiMass + Pt2); 1324 } 1324 } 1325 while ( (InitialMass < MassMt + AntiMassMt) 1325 while ( (InitialMass < MassMt + AntiMassMt) && ++loopCounter < maxNumberOfLoops ); 1326 1326 1327 SigmaQT = SigmaQTw; 1327 SigmaQT = SigmaQTw; 1328 1328 >> 1329 if ( loopCounter >= maxNumberOfLoops ) { >> 1330 AvailablePz2 = 0.0; >> 1331 } >> 1332 1329 AvailablePz2= sqr(InitialMass*InitialMass - 1333 AvailablePz2= sqr(InitialMass*InitialMass - sqr(MassMt) - sqr(AntiMassMt)) - 1330 4.*sqr(MassMt*AntiMassMt); 1334 4.*sqr(MassMt*AntiMassMt); 1331 1335 1332 AvailablePz2 /=(4.*InitialMass*InitialMass) 1336 AvailablePz2 /=(4.*InitialMass*InitialMass); 1333 AvailablePz = std::sqrt(AvailablePz2); 1337 AvailablePz = std::sqrt(AvailablePz2); 1334 1338 1335 G4double Px=Pt.getX(); 1339 G4double Px=Pt.getX(); 1336 G4double Py=Pt.getY(); 1340 G4double Py=Pt.getY(); 1337 1341 1338 Mom->setPx(Px); Mom->setPy(Py); Mom->setPz( 1342 Mom->setPx(Px); Mom->setPy(Py); Mom->setPz(AvailablePz); 1339 Mom->setE(std::sqrt(sqr(MassMt)+AvailablePz 1343 Mom->setE(std::sqrt(sqr(MassMt)+AvailablePz2)); 1340 1344 1341 AntiMom->setPx(-Px); AntiMom->setPy(-Py); A 1345 AntiMom->setPx(-Px); AntiMom->setPy(-Py); AntiMom->setPz(-AvailablePz); 1342 AntiMom->setE (std::sqrt(sqr(AntiMassMt)+Av 1346 AntiMom->setE (std::sqrt(sqr(AntiMassMt)+AvailablePz2)); 1343 1347 1344 #ifdef debug_LUNDfragmentation 1348 #ifdef debug_LUNDfragmentation 1345 G4cout<<"Fmass Mom "<<Mom->getX()<<" 1349 G4cout<<"Fmass Mom "<<Mom->getX()<<" "<<Mom->getY()<<" "<<Mom->getZ()<<" "<<Mom->getT()<<G4endl; 1346 G4cout<<"Smass Mom "<<AntiMom->getX() 1350 G4cout<<"Smass Mom "<<AntiMom->getX()<<" "<<AntiMom->getY()<<" "<<AntiMom->getZ() 1347 <<" "<<AntiMom->getT()<<G4endl; 1351 <<" "<<AntiMom->getT()<<G4endl; 1348 #endif 1352 #endif 1349 } 1353 } 1350 1354 1351 //------------------------------------------- 1355 //------------------------------------------------------------------------ 1352 1356 1353 G4double G4LundStringFragmentation::lambda(G4 1357 G4double G4LundStringFragmentation::lambda(G4double S, G4double m1_Sqr, G4double m2_Sqr) 1354 { 1358 { 1355 G4double lam = sqr(S - m1_Sqr - m2_Sqr) - 4 1359 G4double lam = sqr(S - m1_Sqr - m2_Sqr) - 4.*m1_Sqr*m2_Sqr; 1356 return lam; 1360 return lam; 1357 } 1361 } 1358 1362 1359 // ------------------------------------------ 1363 // -------------------------------------------------------------- 1360 G4LundStringFragmentation::~G4LundStringFragm 1364 G4LundStringFragmentation::~G4LundStringFragmentation() 1361 {} 1365 {} 1362 1366 1363 1367