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 "G4QGSMFragmentation.hh" 33 #include "G4QGSMFragmentation.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 "G4ios.hh" 37 #include "G4ios.hh" 38 #include "G4FragmentingString.hh" 38 #include "G4FragmentingString.hh" 39 #include "G4DiQuarks.hh" 39 #include "G4DiQuarks.hh" 40 #include "G4Quarks.hh" 40 #include "G4Quarks.hh" 41 #include "G4HadronicParameters.hh" << 41 42 #include "G4Pow.hh" 42 #include "G4Pow.hh" 43 43 44 //#define debug_QGSMfragmentation 44 //#define debug_QGSMfragmentation 45 45 46 // Class G4QGSMFragmentation 46 // Class G4QGSMFragmentation 47 //******************************************** 47 //**************************************************************************************** 48 48 49 G4QGSMFragmentation::G4QGSMFragmentation() 49 G4QGSMFragmentation::G4QGSMFragmentation() 50 { 50 { 51 SigmaQT = 0.45 * GeV; << 52 << 53 MassCut = 0.35*GeV; 51 MassCut = 0.35*GeV; 54 52 55 SetStrangenessSuppression((1.0 - 0.16)/2.) 53 SetStrangenessSuppression((1.0 - 0.16)/2.); 56 54 57 // Check if charmed and bottom hadrons are << 55 // For the time being, set to 0.0 the probabilities for c-cbar and b-bbar creation. 58 // set the non-zero probabilities for c-cb << 56 SetProbCCbar(0.0); //(0.0033); // According to O.I. Piskunova Yad. Fiz. 56 (1993) 1094 59 // else set them to 0.0. If these probabil << 57 SetProbBBbar(0.0); //(5.0e-5); // According to O.I. Piskunova Yad. Fiz. 56 (1993) 1094 60 // hadrons can't/can be created during the << 58 61 // (i.e. not heavy) projectile hadron nucl << 62 if ( G4HadronicParameters::Instance()->Ena << 63 SetProbCCbar(0.0002); // According to O << 64 SetProbBBbar(5.0e-5); // According to O << 65 } else { << 66 SetProbCCbar(0.0); << 67 SetProbBBbar(0.0); << 68 } << 69 << 70 SetDiquarkSuppression(0.32); 59 SetDiquarkSuppression(0.32); 71 SetDiquarkBreakProbability(0.7); 60 SetDiquarkBreakProbability(0.7); 72 61 73 SetMinMasses(); 62 SetMinMasses(); 74 63 75 arho = 0.5; // alpha_rho0 64 arho = 0.5; // alpha_rho0 76 aphi = 0.0; // alpha_fi 65 aphi = 0.0; // alpha_fi 77 aJPs =-2.2; // alpha_J/Psi 66 aJPs =-2.2; // alpha_J/Psi 78 aUps =-8.0; // alpha_Y ??? O. Pisk 67 aUps =-8.0; // alpha_Y ??? O. Piskunova Yad. Phys. 56 (1993) 1094. 79 68 80 aksi =-1.0; 69 aksi =-1.0; 81 alft = 0.5; // 2 * alpha'_R *<Pt^2> 70 alft = 0.5; // 2 * alpha'_R *<Pt^2> 82 71 83 an = -0.5 ; 72 an = -0.5 ; 84 ala = -0.75; // an - arho/2 + aphi/2 73 ala = -0.75; // an - arho/2 + aphi/2 85 alaC = an - arho/2.0 + aJPs/2.0; 74 alaC = an - arho/2.0 + aJPs/2.0; 86 alaB = an - arho/2.0 + aUps/2.0; 75 alaB = an - arho/2.0 + aUps/2.0; 87 aXi = 0.0; // ?? 76 aXi = 0.0; // ?? 88 aXiC = 0.0; // ?? 77 aXiC = 0.0; // ?? 89 aXiB = 0.0; // ?? 78 aXiB = 0.0; // ?? 90 aXiCC = 0.0; // ?? 79 aXiCC = 0.0; // ?? 91 aXiCB = 0.0; // ?? 80 aXiCB = 0.0; // ?? 92 aXiBB = 0.0; // ?? 81 aXiBB = 0.0; // ?? 93 82 94 SetFFq2q(); 83 SetFFq2q(); 95 SetFFq2qq(); 84 SetFFq2qq(); 96 SetFFqq2q(); 85 SetFFqq2q(); 97 SetFFqq2qq(); 86 SetFFqq2qq(); 98 // d u s c b 87 // d u s c b 99 G4int Index[5][5] = { { 0, 1, 2, 3, 4 } 88 G4int Index[5][5] = { { 0, 1, 2, 3, 4 }, // d 100 { 1, 5, 6, 7, 8 } 89 { 1, 5, 6, 7, 8 }, // u 101 { 2, 6, 9, 10, 11 } 90 { 2, 6, 9, 10, 11 }, // s 102 { 3, 7, 10, 12, 13 } 91 { 3, 7, 10, 12, 13 }, // c 103 { 4, 8, 11, 13, 14 } 92 { 4, 8, 11, 13, 14 } }; // b 104 for (G4int i = 0; i < 5; i++ ) { 93 for (G4int i = 0; i < 5; i++ ) { 105 for ( G4int j = 0; j < 5; j++ ) { 94 for ( G4int j = 0; j < 5; j++ ) { 106 IndexDiQ[i][j] = Index[i][j]; 95 IndexDiQ[i][j] = Index[i][j]; 107 } 96 } 108 }; 97 }; 109 } 98 } 110 99 111 G4QGSMFragmentation::~G4QGSMFragmentation() 100 G4QGSMFragmentation::~G4QGSMFragmentation() 112 {} 101 {} 113 102 114 //-------------------------------------------- 103 //---------------------------------------------------------------------------------------------------------- 115 104 116 G4KineticTrackVector* G4QGSMFragmentation::Fra 105 G4KineticTrackVector* G4QGSMFragmentation::FragmentString(const G4ExcitedString& theString) 117 { 106 { 118 107 119 G4FragmentingString aString(theString); 108 G4FragmentingString aString(theString); 120 SetMinimalStringMass(&aString); 109 SetMinimalStringMass(&aString); 121 110 122 #ifdef debug_QGSMfragmentation 111 #ifdef debug_QGSMfragmentation 123 G4cout<<G4endl<<"QGSM StringFragm: String Ma 112 G4cout<<G4endl<<"QGSM StringFragm: String Mass " 124 <<theString.Get4M 113 <<theString.Get4Momentum().mag()<<" Pz " 125 <<theString.Get4M 114 <<theString.Get4Momentum().pz() 126 <<"-------------- 115 <<"------------------------------------"<<G4endl; 127 G4cout<<"String ends Direct "<<theString.Get 116 G4cout<<"String ends Direct "<<theString.GetLeftParton()->GetPDGcode()<<" " 128 <<theString.Get 117 <<theString.GetRightParton()->GetPDGcode()<<" " 129 <<theString.Get 118 <<theString.GetDirection()<< G4endl; 130 G4cout<<"Left mom "<<theString.GetLeftParto 119 G4cout<<"Left mom "<<theString.GetLeftParton()->Get4Momentum()<<G4endl; 131 G4cout<<"Right mom "<<theString.GetRightPart 120 G4cout<<"Right mom "<<theString.GetRightParton()->Get4Momentum()<<G4endl; 132 G4cout<<"Check for Fragmentation "<<G4endl; 121 G4cout<<"Check for Fragmentation "<<G4endl; 133 #endif 122 #endif 134 123 135 // Can no longer modify Parameters for Fragm 124 // Can no longer modify Parameters for Fragmentation. 136 PastInitPhase=true; 125 PastInitPhase=true; 137 126 138 // Check if string has enough mass to fragme 127 // Check if string has enough mass to fragment... 139 G4KineticTrackVector * LeftVector=NULL; 128 G4KineticTrackVector * LeftVector=NULL; 140 129 141 if ( !IsItFragmentable(&aString) ) { 130 if ( !IsItFragmentable(&aString) ) { 142 LeftVector=ProduceOneHadron(&theString); 131 LeftVector=ProduceOneHadron(&theString); 143 132 144 #ifdef debug_QGSMfragmentation 133 #ifdef debug_QGSMfragmentation 145 if ( LeftVector != 0 ) G4cout<<"Non fragm 134 if ( LeftVector != 0 ) G4cout<<"Non fragmentable - the string is converted to one hadron "<<G4endl; 146 #endif 135 #endif 147 136 148 if ( LeftVector == nullptr ) LeftVector = << 137 if ( LeftVector != 0 ) return LeftVector; 149 return LeftVector; << 150 } 138 } 151 139 152 #ifdef debug_QGSMfragmentation 140 #ifdef debug_QGSMfragmentation 153 G4cout<<"The string will be fragmented. "<<G 141 G4cout<<"The string will be fragmented. "<<G4endl; 154 #endif 142 #endif 155 143 156 LeftVector = new G4KineticTrackVector; 144 LeftVector = new G4KineticTrackVector; 157 G4KineticTrackVector * RightVector=new G4Kin 145 G4KineticTrackVector * RightVector=new G4KineticTrackVector; 158 146 159 G4ExcitedString *theStringInCMS=CopyExcited( 147 G4ExcitedString *theStringInCMS=CopyExcited(theString); 160 G4LorentzRotation toCms=theStringInCMS->Tran 148 G4LorentzRotation toCms=theStringInCMS->TransformToAlignedCms(); 161 149 162 G4bool success=false, inner_sucess=true; 150 G4bool success=false, inner_sucess=true; 163 G4int attempt=0; 151 G4int attempt=0; 164 while ( !success && attempt++ < StringLoopIn 152 while ( !success && attempt++ < StringLoopInterrupt ) /* Loop checking, 07.08.2015, A.Ribon */ 165 { 153 { 166 #ifdef debug_QGSMfragmentation 154 #ifdef debug_QGSMfragmentation 167 G4cout<<"Loop_toFrag "<<theStr 155 G4cout<<"Loop_toFrag "<<theStringInCMS->GetLeftParton()->GetPDGcode()<<" " 168 <<theStr 156 <<theStringInCMS->GetRightParton()->GetPDGcode()<<" " 169 <<theStr 157 <<theStringInCMS->GetDirection()<< G4endl; 170 #endif 158 #endif 171 159 172 G4FragmentingString *currentString=new G4F 160 G4FragmentingString *currentString=new G4FragmentingString(*theStringInCMS); 173 161 174 std::for_each(LeftVector->begin(), LeftVec 162 std::for_each(LeftVector->begin(), LeftVector->end(), DeleteKineticTrack()); 175 LeftVector->clear(); 163 LeftVector->clear(); 176 std::for_each(RightVector->begin(), RightV 164 std::for_each(RightVector->begin(), RightVector->end(), DeleteKineticTrack()); 177 RightVector->clear(); 165 RightVector->clear(); 178 166 179 inner_sucess=true; // set false on failur 167 inner_sucess=true; // set false on failure.. 180 const G4int maxNumberOfLoops = 168 const G4int maxNumberOfLoops = 1000; 181 G4int loopCounter = -1; 169 G4int loopCounter = -1; 182 while (! StopFragmenting(currentString) && 170 while (! StopFragmenting(currentString) && ++loopCounter < maxNumberOfLoops ) /* Loop checking, 07.08.2015, A.Ribon */ 183 { // Split current string into hadron + n 171 { // Split current string into hadron + new string 184 172 185 #ifdef debug_QGSMfragm 173 #ifdef debug_QGSMfragmentation 186 G4cout<<"The string ca 174 G4cout<<"The string can fragment. "<<G4endl;; 187 #endif 175 #endif 188 G4FragmentingString *newString=0; // us 176 G4FragmentingString *newString=0; // used as output from SplitUp... 189 G4KineticTrack * Hadron=Splitup(currentS 177 G4KineticTrack * Hadron=Splitup(currentString,newString); 190 178 191 if ( Hadron != 0 ) 179 if ( Hadron != 0 ) 192 { 180 { 193 #ifdef debug_QGSMfr 181 #ifdef debug_QGSMfragmentation 194 G4cout<<"Hadron pro 182 G4cout<<"Hadron prod at fragm. "<<Hadron->GetDefinition()->GetParticleName()<<G4endl; 195 #endif 183 #endif 196 // To close the pro 184 // To close the production of hadrons at fragmentation stage 197 if ( currentString->GetDecayDirection 185 if ( currentString->GetDecayDirection() > 0 ) 198 LeftVector->push_back(Hadron); 186 LeftVector->push_back(Hadron); 199 else 187 else 200 RightVector->push_back(Hadron); 188 RightVector->push_back(Hadron); 201 189 202 delete currentString; 190 delete currentString; 203 currentString=newString; 191 currentString=newString; 204 192 205 } else { 193 } else { 206 194 207 #ifdef debug_QGSMfr 195 #ifdef debug_QGSMfragmentation 208 G4cout<<"abandon .. 196 G4cout<<"abandon ... start from the beginning ---------------"<<G4endl; 209 #endif 197 #endif 210 198 211 // Abandon ... start from the beginni 199 // Abandon ... start from the beginning 212 if (newString) delete newString; 200 if (newString) delete newString; 213 inner_sucess=false; 201 inner_sucess=false; 214 break; 202 break; 215 } 203 } 216 } 204 } 217 if ( loopCounter >= maxNumberO 205 if ( loopCounter >= maxNumberOfLoops ) { 218 inner_sucess=false; 206 inner_sucess=false; 219 } 207 } 220 208 221 // Split current string into 2 final Hadro 209 // Split current string into 2 final Hadrons 222 #ifdef debug_QGSMfragmentation 210 #ifdef debug_QGSMfragmentation 223 if( inner_sucess ) { << 211 G4cout<<"Split remaining string into 2 final hadrons."<<G4endl; 224 G4cout<<"Split remaining str << 225 } else { << 226 G4cout<<" New attempt to fragment string << 227 } << 228 #endif 212 #endif 229 // To the close production of 213 // To the close production of hadrons at last string decay 230 if ( inner_sucess && 214 if ( inner_sucess && 231 SplitLast(currentString,LeftVector, R 215 SplitLast(currentString,LeftVector, RightVector) ) 232 { 216 { 233 success=true; 217 success=true; 234 } 218 } 235 delete currentString; 219 delete currentString; 236 } 220 } 237 221 238 delete theStringInCMS; 222 delete theStringInCMS; 239 223 240 if ( ! success ) 224 if ( ! success ) 241 { 225 { 242 std::for_each(LeftVector->begin(), LeftVec 226 std::for_each(LeftVector->begin(), LeftVector->end(), DeleteKineticTrack()); 243 LeftVector->clear(); 227 LeftVector->clear(); 244 std::for_each(RightVector->begin(), RightV 228 std::for_each(RightVector->begin(), RightVector->end(), DeleteKineticTrack()); 245 delete RightVector; 229 delete RightVector; 246 return LeftVector; 230 return LeftVector; 247 } 231 } 248 232 249 // Join Left- and RightVector into LeftVecto 233 // Join Left- and RightVector into LeftVector in correct order. 250 while(!RightVector->empty()) /* Loop checki 234 while(!RightVector->empty()) /* Loop checking, 07.08.2015, A.Ribon */ 251 { 235 { 252 LeftVector->push_back(RightVector->back( 236 LeftVector->push_back(RightVector->back()); 253 RightVector->erase(RightVector->end()-1) 237 RightVector->erase(RightVector->end()-1); 254 } 238 } 255 delete RightVector; 239 delete RightVector; 256 240 257 CalculateHadronTimePosition(theString.Get4Mo 241 CalculateHadronTimePosition(theString.Get4Momentum().mag(), LeftVector); 258 242 259 G4LorentzRotation toObserverFrame(toCms.inve 243 G4LorentzRotation toObserverFrame(toCms.inverse()); 260 244 261 for (size_t C1 = 0; C1 < LeftVector->size(); 245 for (size_t C1 = 0; C1 < LeftVector->size(); C1++) 262 { 246 { 263 G4KineticTrack* Hadron = LeftVector->oper 247 G4KineticTrack* Hadron = LeftVector->operator[](C1); 264 G4LorentzVector Momentum = Hadron->Get4Mo 248 G4LorentzVector Momentum = Hadron->Get4Momentum(); 265 Momentum = toObserverFrame*Momentum; 249 Momentum = toObserverFrame*Momentum; 266 Hadron->Set4Momentum(Momentum); 250 Hadron->Set4Momentum(Momentum); 267 G4LorentzVector Coordinate(Hadron->GetPos 251 G4LorentzVector Coordinate(Hadron->GetPosition(), Hadron->GetFormationTime()); 268 Momentum = toObserverFrame*Coordinate; 252 Momentum = toObserverFrame*Coordinate; 269 Hadron->SetFormationTime(Momentum.e()); 253 Hadron->SetFormationTime(Momentum.e()); 270 G4ThreeVector aPosition(Momentum.vect()); 254 G4ThreeVector aPosition(Momentum.vect()); 271 Hadron->SetPosition(theString.GetPosition 255 Hadron->SetPosition(theString.GetPosition()+aPosition); 272 } 256 } 273 return LeftVector; 257 return LeftVector; 274 } 258 } 275 259 276 //-------------------------------------------- 260 //---------------------------------------------------------------------------------------------------------- 277 261 278 G4bool G4QGSMFragmentation::IsItFragmentable(c << 262 G4bool G4QGSMFragmentation::IsItFragmentable(const G4FragmentingString * const string) 279 { 263 { 280 return sqr( MinimalStringMass + MassCut ) < << 264 return sqr( PossibleHadronMass(string) + MassCut ) < string->Mass2(); 281 } 265 } 282 266 283 //-------------------------------------------- 267 //---------------------------------------------------------------------------------------------------------- 284 268 285 G4bool G4QGSMFragmentation::StopFragmenting(co << 269 G4bool G4QGSMFragmentation::StopFragmenting(const G4FragmentingString * const string) 286 { 270 { 287 SetMinimalStringMass(string); 271 SetMinimalStringMass(string); 288 if ( MinimalStringMass < 0.0 ) return 272 if ( MinimalStringMass < 0.0 ) return true; 289 273 290 G4double smass = string->Mass(); << 274 if (string->IsAFourQuarkString()) 291 G4double x = (string->IsAFourQuarkString()) << 275 { 292 : 0.66e-6*(smass - MinimalStringMass)*(sma << 276 return G4UniformRand() < G4Exp(-0.005*(string->Mass() - MinimalStringMass)); 293 << 277 } else { 294 G4bool res = true; << 278 G4bool Result = G4UniformRand() < 295 if(x > 0.0) { << 279 G4Exp(-0.66e-6*(string->Mass()*string->Mass() - MinimalStringMass*MinimalStringMass)); 296 res = (x < 200.) ? (G4UniformRand() << 280 // G4bool Result = string->Mass() < MinimalStringMass + 150.*MeV*G4UniformRand(); // a'la LUND >> 281 >> 282 #ifdef debug_QGSMfragmentation >> 283 G4cout<<"StopFragmenting MinimalStringMass string->Mass() "<<MinimalStringMass<<" "<<string->Mass()<<G4endl; >> 284 G4cout<<"StopFragmenting - Yes/No "<<Result<<G4endl; >> 285 #endif >> 286 return Result; 297 } 287 } 298 return res; << 299 } 288 } 300 289 301 //-------------------------------------------- 290 //----------------------------------------------------------------------------- 302 291 303 G4KineticTrack * G4QGSMFragmentation::Splitup( 292 G4KineticTrack * G4QGSMFragmentation::Splitup( G4FragmentingString *string, 304 G4FragmentingStri 293 G4FragmentingString *&newString ) 305 { 294 { 306 #ifdef debug_QGSMfragmentation 295 #ifdef debug_QGSMfragmentation 307 G4cout<<G4endl; 296 G4cout<<G4endl; 308 G4cout<<"Start SplitUP (G4VLongitudinal 297 G4cout<<"Start SplitUP (G4VLongitudinalStringDecay) ========================="<<G4endl; 309 G4cout<<"String partons: " <<string->Ge 298 G4cout<<"String partons: " <<string->GetLeftParton()->GetPDGEncoding()<<" " 310 <<string->Ge 299 <<string->GetRightParton()->GetPDGEncoding()<<" " 311 <<"Direction " <<string->Ge 300 <<"Direction " <<string->GetDecayDirection()<<G4endl; 312 #endif 301 #endif 313 302 314 //... random choice of string end to us 303 //... random choice of string end to use for creating the hadron (decay) 315 G4int SideOfDecay = (G4UniformRand() < 304 G4int SideOfDecay = (G4UniformRand() < 0.5)? 1: -1; 316 if (SideOfDecay < 0) 305 if (SideOfDecay < 0) 317 { 306 { 318 string->SetLeftPartonStable(); 307 string->SetLeftPartonStable(); 319 } else 308 } else 320 { 309 { 321 string->SetRightPartonStable(); 310 string->SetRightPartonStable(); 322 } 311 } 323 312 324 G4ParticleDefinition *newStringEnd; 313 G4ParticleDefinition *newStringEnd; 325 G4ParticleDefinition * HadronDefinition 314 G4ParticleDefinition * HadronDefinition; 326 if (string->DecayIsQuark()) 315 if (string->DecayIsQuark()) 327 { 316 { 328 G4double ProbDqADq = GetDiquarkSuppress(); 317 G4double ProbDqADq = GetDiquarkSuppress(); 329 318 330 G4int NumberOfpossibleBaryons = 2; 319 G4int NumberOfpossibleBaryons = 2; 331 320 332 if (string->GetLeftParton()->GetParticleSu 321 if (string->GetLeftParton()->GetParticleSubType() != "quark") NumberOfpossibleBaryons++; 333 if (string->GetRightParton()->GetParticleS 322 if (string->GetRightParton()->GetParticleSubType() != "quark") NumberOfpossibleBaryons++; 334 323 335 G4double ActualProb = ProbDqADq ; 324 G4double ActualProb = ProbDqADq ; 336 ActualProb *= (1.0-G4Exp(2.0*(1.0 - 325 ActualProb *= (1.0-G4Exp(2.0*(1.0 - string->Mass()/(NumberOfpossibleBaryons*1400.0)))); 337 326 338 SetDiquarkSuppression(ActualProb); 327 SetDiquarkSuppression(ActualProb); 339 328 340 HadronDefinition= QuarkSplitup(strin 329 HadronDefinition= QuarkSplitup(string->GetDecayParton(), newStringEnd); 341 330 342 SetDiquarkSuppression(ProbDqADq); 331 SetDiquarkSuppression(ProbDqADq); 343 } else { 332 } else { 344 HadronDefinition= DiQuarkSplitup(str 333 HadronDefinition= DiQuarkSplitup(string->GetDecayParton(), newStringEnd); 345 } 334 } 346 335 347 if ( HadronDefinition == NULL ) return 336 if ( HadronDefinition == NULL ) return NULL; 348 337 349 #ifdef debug_QGSMfragmentation 338 #ifdef debug_QGSMfragmentation 350 G4cout<<"The parton "<<string->GetDecay 339 G4cout<<"The parton "<<string->GetDecayParton()->GetPDGEncoding()<<" " 351 <<" produces hadron "<<HadronDefi 340 <<" produces hadron "<<HadronDefinition->GetParticleName() 352 <<" and is transformed to "<<newS 341 <<" and is transformed to "<<newStringEnd->GetPDGEncoding()<<G4endl; 353 G4cout<<"The side of the string decay L 342 G4cout<<"The side of the string decay Left/Right (1/-1) "<<SideOfDecay<<G4endl; 354 #endif 343 #endif 355 // create new String from old, ie. keep 344 // create new String from old, ie. keep Left and Right order, but replace decay 356 345 357 newString=new G4FragmentingString(*stri 346 newString=new G4FragmentingString(*string,newStringEnd); // To store possible 358 347 // quark containt of new string 359 348 360 #ifdef debug_QGSMfragmentation 349 #ifdef debug_QGSMfragmentation 361 G4cout<<"An attempt to determine its en 350 G4cout<<"An attempt to determine its energy (SplitEandP)"<<G4endl; 362 #endif 351 #endif 363 G4LorentzVector* HadronMomentum=SplitEa 352 G4LorentzVector* HadronMomentum=SplitEandP(HadronDefinition, string, newString); 364 353 365 delete newString; newString=0; 354 delete newString; newString=0; 366 355 367 G4KineticTrack * Hadron =0; 356 G4KineticTrack * Hadron =0; 368 if ( HadronMomentum != 0 ) { 357 if ( HadronMomentum != 0 ) { 369 358 370 #ifdef debug_QGSMfragmentation 359 #ifdef debug_QGSMfragmentation 371 G4cout<<"The attempt was successful 360 G4cout<<"The attempt was successful"<<G4endl; 372 #endif 361 #endif 373 G4ThreeVector Pos; 362 G4ThreeVector Pos; 374 Hadron = new G4KineticTrack(HadronDefinit 363 Hadron = new G4KineticTrack(HadronDefinition, 0,Pos, *HadronMomentum); 375 364 376 newString=new G4FragmentingString(*st 365 newString=new G4FragmentingString(*string,newStringEnd,HadronMomentum); 377 366 378 delete HadronMomentum; 367 delete HadronMomentum; 379 } 368 } 380 else 369 else 381 { 370 { 382 #ifdef debug_QGSMfragmentation 371 #ifdef debug_QGSMfragmentation 383 G4cout<<"The attempt was not successf 372 G4cout<<"The attempt was not successful !!!"<<G4endl; 384 #endif 373 #endif 385 } 374 } 386 375 387 #ifdef debug_VStringDecay 376 #ifdef debug_VStringDecay 388 G4cout<<"End SplitUP (G4VLongitudinalSt 377 G4cout<<"End SplitUP (G4VLongitudinalStringDecay) ====================="<<G4endl; 389 #endif 378 #endif 390 379 391 return Hadron; 380 return Hadron; 392 } 381 } 393 382 394 //-------------------------------------------- 383 //----------------------------------------------------------------------------- 395 384 396 G4ParticleDefinition *G4QGSMFragmentation::DiQ 385 G4ParticleDefinition *G4QGSMFragmentation::DiQuarkSplitup( G4ParticleDefinition* decay, 397 386 G4ParticleDefinition *&created ) 398 { 387 { 399 //... can Diquark break or not? 388 //... can Diquark break or not? 400 if (G4UniformRand() < DiquarkBreakProb ) / 389 if (G4UniformRand() < DiquarkBreakProb ) //... Diquark break 401 { 390 { 402 G4int stableQuarkEncoding = decay->GetPD 391 G4int stableQuarkEncoding = decay->GetPDGEncoding()/1000; 403 G4int decayQuarkEncoding = (decay->GetPD 392 G4int decayQuarkEncoding = (decay->GetPDGEncoding()/100)%10; 404 393 405 if (G4UniformRand() < 0.5) 394 if (G4UniformRand() < 0.5) 406 { 395 { 407 G4int Swap = stableQuarkEncoding; 396 G4int Swap = stableQuarkEncoding; 408 stableQuarkEncoding = decayQuarkEncod 397 stableQuarkEncoding = decayQuarkEncoding; 409 decayQuarkEncoding = Swap; 398 decayQuarkEncoding = Swap; 410 } 399 } 411 400 412 G4int IsParticle=(decayQuarkEncoding>0) 401 G4int IsParticle=(decayQuarkEncoding>0) ? -1 : +1; // if we have a quark, we need antiquark 413 402 414 G4double StrSup=GetStrangeSuppress(); 403 G4double StrSup=GetStrangeSuppress(); 415 SetStrangenessSuppression((1.0 - 0.07)/2 << 404 SetStrangenessSuppression((1.0 - 0.07)/2.); 416 pDefPair QuarkPair = CreatePartonPair(Is 405 pDefPair QuarkPair = CreatePartonPair(IsParticle,false); // no diquarks wanted 417 SetStrangenessSuppression(StrSup); 406 SetStrangenessSuppression(StrSup); 418 407 419 //... Build new Diquark 408 //... Build new Diquark 420 G4int QuarkEncoding=QuarkPair.second->Ge 409 G4int QuarkEncoding=QuarkPair.second->GetPDGEncoding(); 421 G4int i10 = std::max(std::abs(QuarkEnco 410 G4int i10 = std::max(std::abs(QuarkEncoding), std::abs(stableQuarkEncoding)); 422 G4int i20 = std::min(std::abs(QuarkEnco 411 G4int i20 = std::min(std::abs(QuarkEncoding), std::abs(stableQuarkEncoding)); 423 G4int spin = (i10 != i20 && G4UniformRan 412 G4int spin = (i10 != i20 && G4UniformRand() <= 0.5)? 1 : 3; 424 G4int NewDecayEncoding = -1*IsParticle*( 413 G4int NewDecayEncoding = -1*IsParticle*(i10 * 1000 + i20 * 100 + spin); 425 << 426 created = FindParticle(NewDecayEncoding) 414 created = FindParticle(NewDecayEncoding); 427 G4ParticleDefinition * decayQuark=FindPa 415 G4ParticleDefinition * decayQuark=FindParticle(decayQuarkEncoding); 428 G4ParticleDefinition * had=hadronizer->B 416 G4ParticleDefinition * had=hadronizer->Build(QuarkPair.first, decayQuark); 429 417 430 DecayQuark = decay->GetPDGEncoding(); << 418 DecayQuark = decayQuarkEncoding; 431 NewQuark = NewDecayEncoding; << 419 NewQuark = QuarkPair.first->GetPDGEncoding(); 432 420 433 return had; 421 return had; 434 422 435 } else { //... Diquark does not break 423 } else { //... Diquark does not break 436 424 437 G4int IsParticle=(decay->GetPDGEncoding( 425 G4int IsParticle=(decay->GetPDGEncoding()>0) ? +1 : -1; // if we have a diquark, we need quark) 438 G4double StrSup=GetStrangeSuppress(); / 426 G4double StrSup=GetStrangeSuppress(); // for changing s-sbar production 439 SetStrangenessSuppression((1.0 - 0.07)/2 427 SetStrangenessSuppression((1.0 - 0.07)/2.); 440 pDefPair QuarkPair = CreatePartonPair(Is 428 pDefPair QuarkPair = CreatePartonPair(IsParticle,false); // no diquarks wanted 441 SetStrangenessSuppression(StrSup); 429 SetStrangenessSuppression(StrSup); 442 430 443 created = QuarkPair.second; 431 created = QuarkPair.second; 444 432 445 DecayQuark = decay->GetPDGEncoding(); 433 DecayQuark = decay->GetPDGEncoding(); 446 NewQuark = created->GetPDGEncoding(); 434 NewQuark = created->GetPDGEncoding(); 447 435 448 G4ParticleDefinition * had=hadronizer->B 436 G4ParticleDefinition * had=hadronizer->Build(QuarkPair.first, decay); 449 return had; 437 return had; 450 } 438 } 451 } 439 } 452 440 453 //-------------------------------------------- 441 //----------------------------------------------------------------------------------------- 454 442 455 G4LorentzVector * G4QGSMFragmentation::SplitEa 443 G4LorentzVector * G4QGSMFragmentation::SplitEandP(G4ParticleDefinition * pHadron, 456 444 G4FragmentingString * string, 457 445 G4FragmentingString * NewString) 458 { 446 { 459 G4double HadronMass = pHadron->GetPDGMa 447 G4double HadronMass = pHadron->GetPDGMass(); 460 448 461 SetMinimalStringMass(NewString); 449 SetMinimalStringMass(NewString); 462 450 463 if ( MinimalStringMass < 0.0 ) return n 451 if ( MinimalStringMass < 0.0 ) return nullptr; 464 452 465 #ifdef debug_QGSMfragmentation 453 #ifdef debug_QGSMfragmentation 466 G4cout<<"G4QGSMFragmentation::SplitEand 454 G4cout<<"G4QGSMFragmentation::SplitEandP "<<pHadron->GetParticleName()<<G4endl; 467 G4cout<<"String 4 mom, String M "<<stri 455 G4cout<<"String 4 mom, String M "<<string->Get4Momentum()<<" "<<string->Mass()<<G4endl; 468 G4cout<<"HadM MinimalStringMassLeft Str 456 G4cout<<"HadM MinimalStringMassLeft StringM hM+sM "<<HadronMass<<" "<<MinimalStringMass<<" " 469 <<string->Mass()<<" "<<HadronMass 457 <<string->Mass()<<" "<<HadronMass+MinimalStringMass<<G4endl; 470 #endif 458 #endif 471 459 472 if (HadronMass + MinimalStringMass > st 460 if (HadronMass + MinimalStringMass > string->Mass()) 473 { 461 { 474 #ifdef debug_QGSMfragmentation 462 #ifdef debug_QGSMfragmentation 475 G4cout<<"Mass of the string is not su 463 G4cout<<"Mass of the string is not sufficient to produce the hadron!"<<G4endl; 476 #endif 464 #endif 477 return 0; 465 return 0; 478 } // have to start all over! 466 } // have to start all over! 479 467 480 // calculate and assign hadron transver 468 // calculate and assign hadron transverse momentum component HadronPx andHadronPy 481 G4double StringMT2 = string->MassT2(); 469 G4double StringMT2 = string->MassT2(); 482 G4double StringMT = std::sqrt(StringMT 470 G4double StringMT = std::sqrt(StringMT2); 483 471 484 G4LorentzVector String4Momentum = strin 472 G4LorentzVector String4Momentum = string->Get4Momentum(); 485 String4Momentum.setPz(0.); 473 String4Momentum.setPz(0.); 486 G4ThreeVector StringPt = String4Momentu 474 G4ThreeVector StringPt = String4Momentum.vect(); 487 475 488 G4ThreeVector HadronPt , RemSysPt; 476 G4ThreeVector HadronPt , RemSysPt; 489 G4double HadronMassT2, ResidualMas 477 G4double HadronMassT2, ResidualMassT2; 490 478 491 // Mt distribution is implemented << 492 G4double HadronMt, Pt, Pt2, phi; << 493 << 494 //... sample Pt of the hadron 479 //... sample Pt of the hadron 495 G4int attempt=0; 480 G4int attempt=0; 496 do 481 do 497 { 482 { 498 attempt++; if (attempt > StringLoopIn 483 attempt++; if (attempt > StringLoopInterrupt) return 0; 499 484 500 HadronMt = HadronMass - 200.0*G4Log(G << 485 HadronPt =SampleQuarkPt() + string->DecayPt(); 501 Pt2 = sqr(HadronMt)-sqr(HadronMass); << 502 phi = 2.*pi*G4UniformRand(); << 503 G4ThreeVector SampleQuarkPtw= G4Three << 504 HadronPt =SampleQuarkPtw + string->D << 505 HadronPt.setZ(0); 486 HadronPt.setZ(0); 506 RemSysPt = StringPt - HadronPt; 487 RemSysPt = StringPt - HadronPt; 507 488 508 HadronMassT2 = sqr(HadronMass) + Hadr 489 HadronMassT2 = sqr(HadronMass) + HadronPt.mag2(); 509 ResidualMassT2=sqr(MinimalStringMass) 490 ResidualMassT2=sqr(MinimalStringMass) + RemSysPt.mag2(); 510 491 511 } while (std::sqrt(HadronMassT2) + std: 492 } while (std::sqrt(HadronMassT2) + std::sqrt(ResidualMassT2) > StringMT); /* Loop checking, 07.08.2015, A.Ribon */ 512 493 513 //... sample z to define hadron longit 494 //... sample z to define hadron longitudinal momentum and energy 514 //... but first check the available pha 495 //... but first check the available phase space 515 496 516 G4double Pz2 = (sqr(StringMT2 - HadronM 497 G4double Pz2 = (sqr(StringMT2 - HadronMassT2 - ResidualMassT2) - 517 4*HadronMassT2 * ResidualMassT2)/4./ 498 4*HadronMassT2 * ResidualMassT2)/4./StringMT2; 518 499 519 if ( Pz2 < 0 ) {return 0;} // 500 if ( Pz2 < 0 ) {return 0;} // have to start all over! 520 501 521 //... then compute allowed z region z_ 502 //... then compute allowed z region z_min <= z <= z_max 522 503 523 G4double Pz = std::sqrt(Pz2); 504 G4double Pz = std::sqrt(Pz2); 524 G4double zMin = (std::sqrt(HadronMassT2 505 G4double zMin = (std::sqrt(HadronMassT2+Pz2) - Pz)/std::sqrt(StringMT2); 525 G4double zMax = (std::sqrt(HadronMassT2 506 G4double zMax = (std::sqrt(HadronMassT2+Pz2) + Pz)/std::sqrt(StringMT2); 526 507 527 if (zMin >= zMax) return 0; // have 508 if (zMin >= zMax) return 0; // have to start all over! 528 509 529 G4double z = GetLightConeZ(zMin, zMax, 510 G4double z = GetLightConeZ(zMin, zMax, 530 string->GetDecayParton() 511 string->GetDecayParton()->GetPDGEncoding(), pHadron, 531 HadronPt.x(), HadronPt.y 512 HadronPt.x(), HadronPt.y()); 532 513 533 //... now compute hadron longitudinal m 514 //... now compute hadron longitudinal momentum and energy 534 // longitudinal hadron momentum compone 515 // longitudinal hadron momentum component HadronPz 535 516 536 HadronPt.setZ( 0.5* string->GetDecayDir 517 HadronPt.setZ( 0.5* string->GetDecayDirection() * 537 (z * string->LightConeDecay() - 518 (z * string->LightConeDecay() - 538 HadronMassT2/(z * string->LightCone 519 HadronMassT2/(z * string->LightConeDecay())) ); 539 G4double HadronE = 0.5* (z * string->L 520 G4double HadronE = 0.5* (z * string->LightConeDecay() + 540 HadronMassT2/(z * string->LightConeDe 521 HadronMassT2/(z * string->LightConeDecay()) ); 541 522 542 G4LorentzVector * a4Momentum= new G4Lor 523 G4LorentzVector * a4Momentum= new G4LorentzVector(HadronPt,HadronE); 543 524 544 #ifdef debug_QGSMfragmentation 525 #ifdef debug_QGSMfragmentation 545 G4cout<<"string->GetDecayDirection() st 526 G4cout<<"string->GetDecayDirection() string->LightConeDecay() " 546 <<string->GetDecayDirection()<<" 527 <<string->GetDecayDirection()<<" "<<string->LightConeDecay()<<G4endl; 547 G4cout<<"HadronPt,HadronE "<<HadronPt<< 528 G4cout<<"HadronPt,HadronE "<<HadronPt<<" "<<HadronE<<G4endl; 548 G4cout<<"Out of QGSM SplitEandP "<<G4en 529 G4cout<<"Out of QGSM SplitEandP "<<G4endl; 549 #endif 530 #endif 550 531 551 return a4Momentum; 532 return a4Momentum; 552 } 533 } 553 534 554 //-------------------------------------------- 535 //---------------------------------------------------------------------------------------------------------- 555 536 556 G4double G4QGSMFragmentation::GetLightConeZ(G4 << 537 G4double G4QGSMFragmentation::GetLightConeZ(G4double zmin, G4double zmax, G4int /* PartonEncoding */, 557 G4 << 538 G4ParticleDefinition* /* pHadron */, G4double , G4double ) 558 { 539 { 559 G4double lambda = 2.0*(sqr(ptx)+sqr(pty))/sq << 560 << 561 #ifdef debug_QGSMfragmentation 540 #ifdef debug_QGSMfragmentation 562 G4cout<<"GetLightConeZ zmin zmax Parton pHad << 541 G4cout<<"GetLightConeZ zmin zmax Parton pHadron "<<zmin<<" "<<zmax<<" "<< /* PartonEncoding */ 563 <<" "/*<< pHadron->GetParticleName() * << 542 <<" "<</* pHadron->GetParticleName() */ <<G4endl; 564 #endif 543 #endif 565 544 566 G4double z(0.); 545 G4double z(0.); 567 G4int DiQold(0), DiQnew(0); 546 G4int DiQold(0), DiQnew(0); 568 G4double d1(-1.0), d2(0.); 547 G4double d1(-1.0), d2(0.); 569 G4double invD1(0.),invD2(0.), r1(0.),r2(0.), 548 G4double invD1(0.),invD2(0.), r1(0.),r2(0.),r12(0.); 570 549 571 G4int absDecayQuarkCode = std::abs( DecayQua 550 G4int absDecayQuarkCode = std::abs( DecayQuark ); 572 G4int absNewQuarkCode = std::abs( NewQuark 551 G4int absNewQuarkCode = std::abs( NewQuark ); 573 552 574 G4int q1(0), q2(0); 553 G4int q1(0), q2(0); 575 // q1 = absDecayQuarkCode/1000; q2 = (absDe 554 // q1 = absDecayQuarkCode/1000; q2 = (absDecayQuarkCode % 1000)/100; 576 555 577 G4int qA(0), qB(0); 556 G4int qA(0), qB(0); 578 // qA = absNewQuarkCode/1000; qB = (absNe 557 // qA = absNewQuarkCode/1000; qB = (absNewQuarkCode % 1000)/100; 579 558 580 if ( (absDecayQuarkCode < 6) && (absNewQuark 559 if ( (absDecayQuarkCode < 6) && (absNewQuarkCode < 6) ) { 581 d1 = FFq2q[absDecayQuarkCode-1][absNewQuar 560 d1 = FFq2q[absDecayQuarkCode-1][absNewQuarkCode-1][0]; d2 = FFq2q[absDecayQuarkCode-1][absNewQuarkCode-1][1]; 582 } 561 } 583 562 584 if ( (absDecayQuarkCode < 6) && (absNewQuark 563 if ( (absDecayQuarkCode < 6) && (absNewQuarkCode > 6) ) { 585 qA = absNewQuarkCode/1000; qB = (absNewQu 564 qA = absNewQuarkCode/1000; qB = (absNewQuarkCode % 1000)/100; DiQnew = IndexDiQ[qA-1][qB-1]; 586 d1 = FFq2qq[absDecayQuarkCode-1][DiQnew][0] 565 d1 = FFq2qq[absDecayQuarkCode-1][DiQnew][0]; d2 = FFq2q[absDecayQuarkCode-1][DiQnew][1]; 587 } 566 } 588 567 589 if ( (absDecayQuarkCode > 6) && (absNewQuark 568 if ( (absDecayQuarkCode > 6) && (absNewQuarkCode < 6) ) { 590 q1 = absDecayQuarkCode/1000; q2 = (absDeca 569 q1 = absDecayQuarkCode/1000; q2 = (absDecayQuarkCode % 1000)/100; DiQold = IndexDiQ[q1-1][q2-1]; 591 d1 = FFqq2q[DiQold][absNewQuarkCode-1][0]; 570 d1 = FFqq2q[DiQold][absNewQuarkCode-1][0]; d2 = FFqq2q[DiQold][absNewQuarkCode-1][1]; 592 } 571 } 593 572 594 if ( d1 < 0. ) { 573 if ( d1 < 0. ) { 595 q1 = absDecayQuarkCode/1000; q2 = (absDeca 574 q1 = absDecayQuarkCode/1000; q2 = (absDecayQuarkCode % 1000)/100; DiQold = IndexDiQ[q1-1][q2-1]; 596 qA = absNewQuarkCode/1000; qB = (absNewQ << 575 d1 = FFqq2qq[DiQold][absNewQuarkCode-1][0]; d2 = FFqq2qq[DiQold][absNewQuarkCode-1][1]; 597 d1 = FFqq2qq[DiQold][DiQnew][0]; d2 = FFqq << 598 } 576 } 599 577 600 d2 +=lambda; << 601 d1+=1.0; d2+=1.0; 578 d1+=1.0; d2+=1.0; 602 579 603 invD1=1./d1; invD2=1./d2; 580 invD1=1./d1; invD2=1./d2; 604 581 605 const G4int maxNumberOfLoops = 10000; 582 const G4int maxNumberOfLoops = 10000; 606 G4int loopCounter = 0; 583 G4int loopCounter = 0; 607 do // Jong's algorithm 584 do // Jong's algorithm 608 { 585 { 609 r1=G4Pow::GetInstance()->powA(G4UniformRan 586 r1=G4Pow::GetInstance()->powA(G4UniformRand(),invD1); 610 r2=G4Pow::GetInstance()->powA(G4UniformRan 587 r2=G4Pow::GetInstance()->powA(G4UniformRand(),invD2); 611 r12=r1+r2; 588 r12=r1+r2; 612 z=r1/r12; 589 z=r1/r12; 613 } while ( ( (r12 > 1.0) || !((zmin <= z)&&(z 590 } while ( ( (r12 > 1.0) || !((zmin <= z)&&(z <= zmax))) && 614 ++loopCounter < maxNumberOfLoops ) 591 ++loopCounter < maxNumberOfLoops ); /* Loop checking, 07.08.2015, A.Ribon */ 615 592 616 if ( loopCounter >= maxNumberOfLoops ) { 593 if ( loopCounter >= maxNumberOfLoops ) { 617 z = 0.5*(zmin + zmax); // Just a value be 594 z = 0.5*(zmin + zmax); // Just a value between zmin and zmax, no physics considerations at all! 618 } 595 } 619 596 620 return z; 597 return z; 621 } 598 } 622 599 623 //-------------------------------------------- 600 //----------------------------------------------------------------------------------------- 624 601 625 G4bool G4QGSMFragmentation::SplitLast(G4Fragme 602 G4bool G4QGSMFragmentation::SplitLast(G4FragmentingString * string, 626 G4KineticTrackVector * Lef 603 G4KineticTrackVector * LeftVector, 627 G4KineticTrackVector * Right 604 G4KineticTrackVector * RightVector) 628 { 605 { 629 //... perform last cluster decay 606 //... perform last cluster decay 630 607 631 G4ThreeVector ClusterVel =string->Get4Mome 608 G4ThreeVector ClusterVel =string->Get4Momentum().boostVector(); 632 G4double ResidualMass =string->Mass(); 609 G4double ResidualMass =string->Mass(); 633 610 634 #ifdef debug_QGSMfragmentation 611 #ifdef debug_QGSMfragmentation 635 G4cout<<"Split last----------------------- 612 G4cout<<"Split last-----------------------------------------"<<G4endl; 636 G4cout<<"StrMass "<<ResidualMass<<" q's " 613 G4cout<<"StrMass "<<ResidualMass<<" q's " 637 <<string->GetLeftParton()->GetPartic 614 <<string->GetLeftParton()->GetParticleName()<<" " 638 <<string->GetRightParton()->GetParti 615 <<string->GetRightParton()->GetParticleName()<<G4endl; 639 #endif 616 #endif 640 617 641 G4int cClusterInterrupt = 0; << 642 G4ParticleDefinition *LeftHadron = nullptr 618 G4ParticleDefinition *LeftHadron = nullptr; 643 G4ParticleDefinition *RightHadron = nullpt 619 G4ParticleDefinition *RightHadron = nullptr; 644 const G4int maxNumberOfLoops = 1000; 620 const G4int maxNumberOfLoops = 1000; 645 G4int loopCounter = 0; 621 G4int loopCounter = 0; >> 622 G4bool isOK = false; 646 623 647 G4double LeftHadronMass(0.); G4double Righ 624 G4double LeftHadronMass(0.); G4double RightHadronMass(0.); 648 do 625 do 649 { 626 { 650 if (cClusterInterrupt++ >= ClusterLoop << 651 LeftHadronMass = -MaxMass; RightHadron 627 LeftHadronMass = -MaxMass; RightHadronMass = -MaxMass; 652 628 653 G4ParticleDefinition * quark = nullptr; 629 G4ParticleDefinition * quark = nullptr; 654 string->SetLeftPartonStable(); // to query q 630 string->SetLeftPartonStable(); // to query quark contents.. 655 631 656 if (string->DecayIsQuark() && string->Stable 632 if (string->DecayIsQuark() && string->StableIsQuark() ) 657 { 633 { 658 //... there are quarks on cluster ends 634 //... there are quarks on cluster ends 659 635 660 G4int IsParticle=(string->GetLeftParton()- 636 G4int IsParticle=(string->GetLeftParton()->GetPDGEncoding()>0) ? -1 : +1; 661 // if we have a quark, we need 637 // if we have a quark, we need antiquark or diquark 662 638 663 pDefPair QuarkPair = CreatePartonPair(IsPa 639 pDefPair QuarkPair = CreatePartonPair(IsParticle); 664 quark = QuarkPair.second; 640 quark = QuarkPair.second; 665 641 666 LeftHadron= hadronizer->Build(QuarkPair.fi << 642 LeftHadron= hadronizer->BuildLowSpin(QuarkPair.first, string->GetLeftParton()); 667 if ( LeftHadron == NULL ) continue; << 643 } else { 668 RightHadron = hadronizer->Build(stri << 644 //... there is a Diquark on cluster ends 669 if ( RightHadron == NULL ) continue; << 670 } else if( (!string->DecayIsQuark() && stri << 671 ( string->DecayIsQuark() && !stri << 672 //... there is a Diquark on one of cluster << 673 G4int IsParticle; 645 G4int IsParticle; 674 if ( string->StableIsQuark() ) { 646 if ( string->StableIsQuark() ) { 675 IsParticle=(string->GetLeftParton()->Get 647 IsParticle=(string->GetLeftParton()->GetPDGEncoding()>0) ? -1 : +1; 676 } else { 648 } else { 677 IsParticle=(string->GetLeftParton()->Get 649 IsParticle=(string->GetLeftParton()->GetPDGEncoding()>0) ? +1 : -1; 678 } 650 } 679 651 680 //G4double ProbSaS = 1.0 - 2.0 * G 652 //G4double ProbSaS = 1.0 - 2.0 * GetStrangeSuppress(); 681 //G4double ActualProb = ProbSaS * 1. 653 //G4double ActualProb = ProbSaS * 1.4; 682 //SetStrangenessSuppression((1.0-Act 654 //SetStrangenessSuppression((1.0-ActualProb)/2.0); 683 655 684 pDefPair QuarkPair = CreatePartonPai 656 pDefPair QuarkPair = CreatePartonPair(IsParticle,false); // no diquarks wanted 685 //SetStrangenessSuppression((1.0-Pro 657 //SetStrangenessSuppression((1.0-ProbSaS)/2.0); 686 quark = QuarkPair.second; 658 quark = QuarkPair.second; 687 LeftHadron=hadronizer->Build(QuarkPa << 659 LeftHadron=hadronizer->BuildLowSpin(QuarkPair.first, string->GetLeftParton()); 688 if ( LeftHadron == NULL ) continue; << 689 RightHadron = hadronizer->Build(stri << 690 if ( RightHadron == NULL ) continue; << 691 } else { // Diquark and anti-diquark << 692 if (cClusterInterrupt++ >= ClusterLo << 693 G4int LeftQuark1= string->GetLeftPar << 694 G4int LeftQuark2=(string->GetLeftPar << 695 G4int RightQuark1= string->GetRightP << 696 G4int RightQuark2=(string->GetRightP << 697 if (G4UniformRand()<0.5) { << 698 LeftHadron =hadronizer->Build(Fin << 699 RightHadron =hadronizer->Build(Fin << 700 } else { << 701 LeftHadron =hadronizer->Build(Fin << 702 RightHadron =hadronizer->Build(Fin << 703 } << 704 if ( (LeftHadron == NULL) || (RightHadron << 705 } 660 } 706 LeftHadronMass = LeftHadron->GetPDGMa << 707 RightHadronMass = RightHadron->GetPDGM << 708 //... repeat procedure, if mass of clu << 709 } while ( ( ResidualMass <= LeftHadronMass << 710 && ++loopCounter < maxNumberOfLo << 711 661 712 if ( loopCounter >= maxNumberOfLoops ) { << 662 if ( LeftHadron != nullptr ) { 713 return false; << 663 RightHadron = hadronizer->BuildLowSpin(string->GetRightParton(), quark); >> 664 >> 665 if ( RightHadron != nullptr ) { >> 666 >> 667 LeftHadronMass = LeftHadron->GetPDGMass(); >> 668 RightHadronMass = RightHadron->GetPDGMass(); >> 669 isOK = (ResidualMass > LeftHadronMass + RightHadronMass); >> 670 } >> 671 } >> 672 ++loopCounter; >> 673 if ( loopCounter >= maxNumberOfLoops ) { >> 674 return false; >> 675 } >> 676 //... repeat procedure, if mass of cluster is too low to produce hadrons >> 677 //... ClusterMassCut = 0.15*GeV model parameter 714 } 678 } >> 679 while (isOK == false); >> 680 /* Loop checking, 07.08.2015, A.Ribon */ 715 681 716 //... compute hadron momenta and energies 682 //... compute hadron momenta and energies 717 G4LorentzVector LeftMom, RightMom; 683 G4LorentzVector LeftMom, RightMom; 718 G4ThreeVector Pos; 684 G4ThreeVector Pos; 719 Sample4Momentum(&LeftMom , LeftHadron->Get 685 Sample4Momentum(&LeftMom , LeftHadron->GetPDGMass() , 720 &RightMom, RightHadron->Ge 686 &RightMom, RightHadron->GetPDGMass(), ResidualMass); 721 LeftMom.boost(ClusterVel); 687 LeftMom.boost(ClusterVel); 722 RightMom.boost(ClusterVel); 688 RightMom.boost(ClusterVel); 723 689 724 #ifdef debug_QGSMfragmentation 690 #ifdef debug_QGSMfragmentation 725 G4cout<<LeftHadron->GetParticleName()<<" " 691 G4cout<<LeftHadron->GetParticleName()<<" "<<RightHadron->GetParticleName()<<G4endl; 726 G4cout<<"Left Hadrom P M "<<LeftMom<<" "< 692 G4cout<<"Left Hadrom P M "<<LeftMom<<" "<<LeftMom.mag()<<G4endl; 727 G4cout<<"Right Hadrom P M "<<RightMom<<" " 693 G4cout<<"Right Hadrom P M "<<RightMom<<" "<<RightMom.mag()<<G4endl; 728 #endif 694 #endif 729 695 730 LeftVector->push_back(new G4KineticTrack(L 696 LeftVector->push_back(new G4KineticTrack(LeftHadron, 0, Pos, LeftMom)); 731 RightVector->push_back(new G4KineticTrack( 697 RightVector->push_back(new G4KineticTrack(RightHadron, 0, Pos, RightMom)); 732 698 733 return true; 699 return true; 734 } 700 } 735 701 736 //-------------------------------------------- 702 //---------------------------------------------------------------------------------------------------------- 737 703 738 void G4QGSMFragmentation::Sample4Momentum(G4Lo 704 void G4QGSMFragmentation::Sample4Momentum(G4LorentzVector* Mom , G4double Mass , 739 G4Lo 705 G4LorentzVector* AntiMom, G4double AntiMass, G4double InitialMass) 740 { 706 { 741 #ifdef debug_QGSMfragmentation << 742 G4cout<<"Sample4Momentum Last------------- << 743 G4cout<<" StrMass "<<InitialMass<<" Mass1 << 744 G4cout<<" SumMass "<<Mass+AntiMass<<G4end << 745 #endif << 746 << 747 G4double r_val = sqr(InitialMass*InitialMa 707 G4double r_val = sqr(InitialMass*InitialMass - Mass*Mass - AntiMass*AntiMass) - sqr(2.*Mass*AntiMass); 748 G4double Pabs = (r_val > 0.)? std::sqrt(r_ 708 G4double Pabs = (r_val > 0.)? std::sqrt(r_val)/(2.*InitialMass) : 0; 749 709 750 //... sample unit vector 710 //... sample unit vector 751 G4double pz = 1. - 2.*G4UniformRand(); 711 G4double pz = 1. - 2.*G4UniformRand(); 752 G4double st = std::sqrt(1. - pz * pz)* 712 G4double st = std::sqrt(1. - pz * pz)*Pabs; 753 G4double phi = 2.*pi*G4UniformRand(); 713 G4double phi = 2.*pi*G4UniformRand(); 754 G4double px = st*std::cos(phi); 714 G4double px = st*std::cos(phi); 755 G4double py = st*std::sin(phi); 715 G4double py = st*std::sin(phi); 756 pz *= Pabs; 716 pz *= Pabs; 757 717 758 Mom->setPx(px); Mom->setPy(py); Mom->setPz 718 Mom->setPx(px); Mom->setPy(py); Mom->setPz(pz); 759 Mom->setE(std::sqrt(Pabs*Pabs + Mass*Mass) 719 Mom->setE(std::sqrt(Pabs*Pabs + Mass*Mass)); 760 720 761 AntiMom->setPx(-px); AntiMom->setPy(-py); 721 AntiMom->setPx(-px); AntiMom->setPy(-py); AntiMom->setPz(-pz); 762 AntiMom->setE (std::sqrt(Pabs*Pabs + AntiM 722 AntiMom->setE (std::sqrt(Pabs*Pabs + AntiMass*AntiMass)); 763 } 723 } 764 724 765 //-------------------------------------------- 725 //---------------------------------------------------------------------------------------------------------- 766 726 767 void G4QGSMFragmentation::SetFFq2q() // q-> q 727 void G4QGSMFragmentation::SetFFq2q() // q-> q' + Meson (q anti q') 768 { 728 { 769 for (G4int i=0; i < 5; i++) { 729 for (G4int i=0; i < 5; i++) { 770 FFq2q[i][0][0] = 2.0 ; FFq2q[i][0][1] = -a 730 FFq2q[i][0][0] = 2.0 ; FFq2q[i][0][1] = -arho + alft; // q->d + (q dbar) Pi0, Eta, Eta', Rho0, omega 771 FFq2q[i][1][0] = 2.0 ; FFq2q[i][1][1] = -a 731 FFq2q[i][1][0] = 2.0 ; FFq2q[i][1][1] = -arho + alft; // q->u + (q ubar) Pi-, Rho- 772 FFq2q[i][2][0] = 1.0 ; FFq2q[i][2][1] = -a 732 FFq2q[i][2][0] = 1.0 ; FFq2q[i][2][1] = -aphi + alft; // q->s + (q sbar) K0, K*0 773 FFq2q[i][3][0] = 1.0 ; FFq2q[i][3][1] = -a 733 FFq2q[i][3][0] = 1.0 ; FFq2q[i][3][1] = -aJPs + alft; // q->c + (q+cbar) D-, D*- 774 FFq2q[i][4][0] = 1.0 ; FFq2q[i][4][1] = -a 734 FFq2q[i][4][0] = 1.0 ; FFq2q[i][4][1] = -aUps + alft; // q->b + (q bbar) EtaB, Upsilon 775 } 735 } 776 } 736 } 777 737 778 //-------------------------------------------- 738 //---------------------------------------------------------------------------------------------------------- 779 739 780 void G4QGSMFragmentation::SetFFq2qq() // q-> 740 void G4QGSMFragmentation::SetFFq2qq() // q-> anti (q1'q2') + Baryon (q + q1 + q2) 781 { 741 { 782 for (G4int i=0; i < 5; i++) { 742 for (G4int i=0; i < 5; i++) { 783 FFq2qq[i][ 0][0] = 0.0 ; FFq2qq[i][ 0][1] 743 FFq2qq[i][ 0][0] = 0.0 ; FFq2qq[i][ 0][1] = arho - 2.0*an + alft; // q->dd bar + (q dd) 784 FFq2qq[i][ 1][0] = 0.0 ; FFq2qq[i][ 1][1] 744 FFq2qq[i][ 1][0] = 0.0 ; FFq2qq[i][ 1][1] = arho - 2.0*an + alft; // q->ud bar + (q ud) 785 FFq2qq[i][ 2][0] = 0.0 ; FFq2qq[i][ 2][1] 745 FFq2qq[i][ 2][0] = 0.0 ; FFq2qq[i][ 2][1] = arho - 2.0*ala + alft; // q->sd bar + (q sd) 786 FFq2qq[i][ 3][0] = 0.0 ; FFq2qq[i][ 3][1] 746 FFq2qq[i][ 3][0] = 0.0 ; FFq2qq[i][ 3][1] = arho - 2.0*alaC + alft; // q->cd bar + (q cd) 787 FFq2qq[i][ 4][0] = 0.0 ; FFq2qq[i][ 4][1] 747 FFq2qq[i][ 4][0] = 0.0 ; FFq2qq[i][ 4][1] = arho - 2.0*alaB + alft; // q->bd bar + (q bd) 788 FFq2qq[i][ 5][0] = 0.0 ; FFq2qq[i][ 5][1] 748 FFq2qq[i][ 5][0] = 0.0 ; FFq2qq[i][ 5][1] = arho - 2.0*an + alft; // q->uu bar + (q uu) 789 FFq2qq[i][ 6][0] = 0.0 ; FFq2qq[i][ 6][1] 749 FFq2qq[i][ 6][0] = 0.0 ; FFq2qq[i][ 6][1] = arho - 2.0*ala + alft; // q->su bar + (q su) 790 FFq2qq[i][ 7][0] = 0.0 ; FFq2qq[i][ 7][1] 750 FFq2qq[i][ 7][0] = 0.0 ; FFq2qq[i][ 7][1] = arho - 2.0*alaC + alft; // q->cu bar + (q cu) 791 FFq2qq[i][ 8][0] = 0.0 ; FFq2qq[i][ 8][1] 751 FFq2qq[i][ 8][0] = 0.0 ; FFq2qq[i][ 8][1] = arho - 2.0*alaB + alft; // q->bu bar + (q bu) 792 FFq2qq[i][ 9][0] = 0.0 ; FFq2qq[i][ 9][1] 752 FFq2qq[i][ 9][0] = 0.0 ; FFq2qq[i][ 9][1] = arho - 2.0*aXi + alft; // q->ss bar + (q ss) 793 FFq2qq[i][10][0] = 0.0 ; FFq2qq[i][10][1] 753 FFq2qq[i][10][0] = 0.0 ; FFq2qq[i][10][1] = arho - 2.0*aXiC + alft; // q->cs bar + (q cs) 794 FFq2qq[i][11][0] = 0.0 ; FFq2qq[i][11][1] 754 FFq2qq[i][11][0] = 0.0 ; FFq2qq[i][11][1] = arho - 2.0*aXiB + alft; // q->bs bar + (q bc) 795 FFq2qq[i][12][0] = 0.0 ; FFq2qq[i][12][1] 755 FFq2qq[i][12][0] = 0.0 ; FFq2qq[i][12][1] = arho - 2.0*aXiCC + alft; // q->cc bar + (q cc) 796 FFq2qq[i][13][0] = 0.0 ; FFq2qq[i][13][1] 756 FFq2qq[i][13][0] = 0.0 ; FFq2qq[i][13][1] = arho - 2.0*aXiCB + alft; // q->bc bar + (q bc) 797 FFq2qq[i][14][0] = 0.0 ; FFq2qq[i][14][1] 757 FFq2qq[i][14][0] = 0.0 ; FFq2qq[i][14][1] = arho - 2.0*aXiBB + alft; // q->bb bar + (q bb) 798 } 758 } 799 } 759 } 800 760 801 //-------------------------------------------- 761 //---------------------------------------------------------------------------------------------------------- 802 762 803 void G4QGSMFragmentation::SetFFqq2q() // q1q2 763 void G4QGSMFragmentation::SetFFqq2q() // q1q2-> anti(q') + Baryon (q1 + q2 + q') 804 { 764 { 805 for (G4int i=0; i < 15; i++) { 765 for (G4int i=0; i < 15; i++) { 806 FFqq2q[i][0][0] = 2.0*(arho - an); FFqq2q[ 766 FFqq2q[i][0][0] = 2.0*(arho - an); FFqq2q[i][0][1] = -arho + alft; 807 FFqq2q[i][1][0] = 2.0*(arho - an); FFqq2q[ 767 FFqq2q[i][1][0] = 2.0*(arho - an); FFqq2q[i][1][1] = -arho + alft; 808 FFqq2q[i][2][0] = 2.0*(arho - an); FFqq2q[ 768 FFqq2q[i][2][0] = 2.0*(arho - an); FFqq2q[i][2][1] = -aphi + alft; 809 FFqq2q[i][3][0] = 2.0*(arho - an); FFqq2q[ 769 FFqq2q[i][3][0] = 2.0*(arho - an); FFqq2q[i][3][1] = -aJPs + alft; 810 FFqq2q[i][4][0] = 2.0*(arho - an); FFqq2q[ 770 FFqq2q[i][4][0] = 2.0*(arho - an); FFqq2q[i][4][1] = -aUps + alft; 811 } 771 } 812 } 772 } 813 773 814 //-------------------------------------------- 774 //---------------------------------------------------------------------------------------------------------- 815 775 816 void G4QGSMFragmentation::SetFFqq2qq() // q1( 776 void G4QGSMFragmentation::SetFFqq2qq() // q1(q2)-> q'(q2) + Meson(q1 anti q') 817 { 777 { 818 for (G4int i=0; i < 15; i++) { 778 for (G4int i=0; i < 15; i++) { 819 FFqq2qq[i][0][0] = 0. ; FFqq2qq[i][0][1] 779 FFqq2qq[i][0][0] = 0. ; FFqq2qq[i][0][1] = 2.0*arho - 2.0*an -arho + alft; // dd -> dd + Pi0 (d d bar) 820 FFqq2qq[i][1][0] = 0. ; FFqq2qq[i][1][1] 780 FFqq2qq[i][1][0] = 0. ; FFqq2qq[i][1][1] = 2.0*arho - 2.0*an -arho + alft; // dd -> ud + Pi- (d u bar) 821 FFqq2qq[i][2][0] = 0. ; FFqq2qq[i][2][1] 781 FFqq2qq[i][2][0] = 0. ; FFqq2qq[i][2][1] = 2.0*arho - 2.0*an -aphi + alft; // dd -> sd + K0 (d s bar) 822 FFqq2qq[i][3][0] = 0. ; FFqq2qq[i][3][1] 782 FFqq2qq[i][3][0] = 0. ; FFqq2qq[i][3][1] = 2.0*arho - 2.0*an -aJPs + alft; // dd -> cd + D- (d c bar) 823 FFqq2qq[i][4][0] = 0. ; FFqq2qq[i][4][1] 783 FFqq2qq[i][4][0] = 0. ; FFqq2qq[i][4][1] = 2.0*arho - 2.0*an -aUps + alft; // dd -> bd + B0 (d b bar) 824 } 784 } 825 } 785 } 826 786 827 787