- Timestamp:
- Apr 6, 2009, 12:30:29 PM (16 years ago)
- Location:
- trunk/source/processes/hadronic/models/parton_string/hadronization/src
- Files:
-
- 6 edited
-
G4FragmentingString.cc (modified) (5 diffs)
-
G4HadronBuilder.cc (modified) (1 diff)
-
G4LundStringFragmentation.cc (modified) (11 diffs)
-
G4QGSMFragmentation.cc (modified) (3 diffs)
-
G4VKinkyStringDecay.cc (modified) (1 diff)
-
G4VLongitudinalStringDecay.cc (modified) (15 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/source/processes/hadronic/models/parton_string/hadronization/src/G4FragmentingString.cc
r819 r962 80 80 if ( old.decaying == Left ) 81 81 { 82 //G4cout<<" Left "<<G4endl; 83 //G4cout<<"Pt right "<<Ptright<<G4endl; 84 //G4cout<<"Pt left "<<Ptleft <<G4endl; 82 85 RightParton= old.RightParton; 83 86 Ptright = old.Ptright; … … 85 88 Ptleft = old.Ptleft - momentum->vect(); 86 89 Ptleft.setZ(0.); 90 //G4cout<<"Pt right "<<Ptright<<G4endl; 91 //G4cout<<"Pt left "<<Ptleft <<G4endl; 87 92 } else if ( old.decaying == Right ) 88 93 { 94 //G4cout<<" Right "<<G4endl; 95 //G4cout<<"Pt right "<<Ptright<<G4endl; 96 //G4cout<<"Pt left "<<Ptleft <<G4endl; 89 97 RightParton = newdecay; 90 98 Ptright = old.Ptright - momentum->vect(); … … 92 100 LeftParton = old.LeftParton; 93 101 Ptleft = old.Ptleft; 102 //G4cout<<"Pt right "<<Ptright<<G4endl; 103 //G4cout<<"Pt left "<<Ptleft <<G4endl; 94 104 } else 95 105 { … … 106 116 //--------------------------------------------------------------------------------- 107 117 118 G4FragmentingString::G4FragmentingString(const G4FragmentingString &old, // Uzhi 119 G4ParticleDefinition * newdecay) // Uzhi 120 { // Uzhi 121 decaying=None; // Uzhi 122 if ( old.decaying == Left ) // Uzhi 123 { // Uzhi 124 RightParton= old.RightParton; // Uzhi 125 LeftParton = newdecay; // Uzhi 126 } else if ( old.decaying == Right ) // Uzhi 127 { // Uzhi 128 RightParton = newdecay; // Uzhi 129 LeftParton = old.LeftParton; // Uzhi 130 } else // Uzhi 131 { 132 throw G4HadronicException(__FILE__, __LINE__, "G4FragmentingString::G4FragmentingString: no decay Direction defined"); 133 } 134 } 135 136 137 //--------------------------------------------------------------------------------- 138 108 139 G4FragmentingString::~G4FragmentingString() 109 140 {} … … 215 246 return std::sqrt(this->Mass2()); 216 247 } 248 249 G4double G4FragmentingString::MassT2() const 250 { 251 return Pplus*Pminus; 252 } -
trunk/source/processes/hadronic/models/parton_string/hadronization/src/G4HadronBuilder.cc
r819 r962 25 25 // 26 26 // 27 // $Id: G4HadronBuilder.cc,v 1. 6 2006/06/29 20:55:01 gunterExp $28 // GEANT4 tag $Name: $27 // $Id: G4HadronBuilder.cc,v 1.7 2008/04/25 14:20:14 vuzhinsk Exp $ 28 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 29 29 // 30 30 // ----------------------------------------------------------------------------- -
trunk/source/processes/hadronic/models/parton_string/hadronization/src/G4LundStringFragmentation.cc
r819 r962 25 25 // 26 26 // 27 // $Id: G4LundStringFragmentation.cc,v 1. 7 2007/04/24 14:55:23 gunterExp $28 // GEANT4 tag $Name: geant4-09-0 1-patch-02 $27 // $Id: G4LundStringFragmentation.cc,v 1.13 2008/06/23 09:17:10 gcosmo Exp $ 28 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 1.8 29 29 // 30 30 // ----------------------------------------------------------------------------- … … 41 41 42 42 // Class G4LundStringFragmentation 43 //************************************************************************************* ***43 //************************************************************************************* 44 44 45 45 G4LundStringFragmentation::G4LundStringFragmentation() 46 46 { 47 MinimalStringMass = 0.; // Uzhi 48 MinimalStringMass2 = 0.; // Uzhi 49 WminLUND = 1.*GeV; // Uzhi 50 SmoothParam = 0.2; // Uzhi 51 47 // ------ For estimation of a minimal string mass --------------- 48 Mass_of_light_quark =140.*MeV; 49 Mass_of_heavy_quark =500.*MeV; 50 Mass_of_string_junction=720.*MeV; 51 // ------ An estimated minimal string mass ---------------------- 52 MinimalStringMass = 0.; 53 MinimalStringMass2 = 0.; 54 // ------ Minimal invariant mass used at a string fragmentation - 55 WminLUND = 0.7*GeV; // Uzhi 0.8 1.5 56 // ------ Smooth parameter used at a string fragmentation for --- 57 // ------ smearinr sharp mass cut-off --------------------------- 58 SmoothParam = 0.2; 59 60 SetStringTensionParameter(0.25); // Uzhi 20 June 08 52 61 } 53 62 54 // G4LundStringFragmentation::G4LundStringFragmentation(G4double sigmaPt) 55 // : G4VLongitudinalStringDecay(sigmaPt) 56 // { 57 // } 58 63 // -------------------------------------------------------------- 59 64 G4LundStringFragmentation::G4LundStringFragmentation(const G4LundStringFragmentation &) : G4VLongitudinalStringDecay() 60 65 { 61 66 } 62 67 63 64 68 G4LundStringFragmentation::~G4LundStringFragmentation() 65 69 { 66 70 } 67 71 68 //************************************************************************************* ***72 //************************************************************************************* 69 73 70 74 const G4LundStringFragmentation & G4LundStringFragmentation::operator=(const G4LundStringFragmentation &) … … 84 88 } 85 89 86 //**************************************************************************************** 87 //---------------------------------------------------------------------------------------------------------- 88 89 G4KineticTrackVector* G4LundStringFragmentation::FragmentString(const G4ExcitedString& theString) 90 //-------------------------------------------------------------------------------------- 91 void G4LundStringFragmentation::SetMinimalStringMass(const G4FragmentingString * const string) // Uzhi 92 { 93 /* 94 G4cout<<"In SetMinMass -------------------"<<std::sqrt(string->Mass2())<<G4endl; 95 G4cout<<string->GetLeftParton()->GetPDGEncoding()<<" "<< 96 string->GetRightParton()->GetPDGEncoding()<<G4endl; 97 */ 98 G4double EstimatedMass=0.; 99 G4int Number_of_quarks=0; 100 101 G4int Qleft =std::abs(string->GetLeftParton()->GetPDGEncoding()); 102 103 if( Qleft > 1000) 104 { 105 Number_of_quarks+=2; 106 G4int q1=Qleft/1000; 107 if( q1 < 3) {EstimatedMass +=Mass_of_light_quark;} 108 if( q1 > 2) {EstimatedMass +=Mass_of_heavy_quark;} 109 110 G4int q2=(Qleft/100)%10; 111 if( q2 < 3) {EstimatedMass +=Mass_of_light_quark;} 112 if( q2 > 2) {EstimatedMass +=Mass_of_heavy_quark;} 113 EstimatedMass +=Mass_of_string_junction; 114 } 115 else 116 { 117 Number_of_quarks++; 118 if( Qleft < 3) {EstimatedMass +=Mass_of_light_quark;} 119 if( Qleft > 2) {EstimatedMass +=Mass_of_heavy_quark;} 120 } 121 122 G4int Qright=std::abs(string->GetRightParton()->GetPDGEncoding()); 123 124 if( Qright > 1000) 125 { 126 Number_of_quarks+=2; 127 G4int q1=Qright/1000; 128 if( q1 < 3) {EstimatedMass +=Mass_of_light_quark;} 129 if( q1 > 2) {EstimatedMass +=Mass_of_heavy_quark;} 130 131 G4int q2=(Qright/100)%10; 132 if( q2 < 3) {EstimatedMass +=Mass_of_light_quark;} 133 if( q2 > 2) {EstimatedMass +=Mass_of_heavy_quark;} 134 EstimatedMass +=Mass_of_string_junction; 135 } 136 else 137 { 138 Number_of_quarks++; 139 if( Qright < 3) {EstimatedMass +=Mass_of_light_quark;} 140 if( Qright > 2) {EstimatedMass +=Mass_of_heavy_quark;} 141 } 142 143 if(Number_of_quarks==2){EstimatedMass +=100.*MeV;} 144 if(Number_of_quarks==3){EstimatedMass += 20.*MeV;} 145 if(Number_of_quarks==4){EstimatedMass -=2.*Mass_of_string_junction; 146 if(EstimatedMass <= 1600.*MeV){EstimatedMass-=200.*MeV;} 147 else {EstimatedMass+=100.*MeV;} 148 } 149 MinimalStringMass=EstimatedMass; 150 SetMinimalStringMass2(EstimatedMass); 151 //G4cout<<"Out SetMinimalStringMass "<<MinimalStringMass<<G4endl; 152 } 153 154 //-------------------------------------------------------------------------------------- 155 void G4LundStringFragmentation::SetMinimalStringMass2( 156 const G4double aValue) 157 { 158 MinimalStringMass2=aValue * aValue; 159 } 160 161 //-------------------------------------------------------------------------------------- 162 G4KineticTrackVector* G4LundStringFragmentation::FragmentString( 163 const G4ExcitedString& theString) 90 164 { 91 165 … … 96 170 97 171 // check if string has enough mass to fragment... 172 173 SetMassCut(160.*MeV); // For LightFragmentationTest it is required 174 // that no one pi-meson can be produced 175 /* 176 G4cout<<G4endl<<"G4LundStringFragmentation::"<<G4endl; 177 G4cout<<"FragmentString Position"<<theString.GetPosition()/fermi<<" "<< 178 theString.GetTimeOfCreation()/fermi<<G4endl; 179 G4cout<<"FragmentString Momentum"<<theString.Get4Momentum()<<theString.Get4Momentum().mag()<<G4endl; 180 */ 98 181 G4KineticTrackVector * LeftVector=LightFragmentationTest(&theString); 99 182 if ( LeftVector != 0 ) { 100 //G4cout<<"Return single hadron from string"<<G4endl; 101 return LeftVector;} 102 103 LeftVector = new G4KineticTrackVector; 183 //G4cout<<"Return single hadron insted of string"<<G4endl; 184 // Uzhi insert 6.05.08 start 185 if(LeftVector->size() == 1){ 186 // One hadron is saved in the interaction 187 LeftVector->operator[](0)->SetFormationTime(theString.GetTimeOfCreation()); 188 LeftVector->operator[](0)->SetPosition(theString.GetPosition()); 189 190 /* // To set large formation time open * 191 LeftVector->operator[](0)->SetFormationTime(theString.GetTimeOfCreation()+100.*fermi); 192 LeftVector->operator[](0)->SetPosition(theString.GetPosition()); 193 G4ThreeVector aPosition(theString.GetPosition().x(), 194 theString.GetPosition().y(), 195 theString.GetPosition().z()+100.*fermi); 196 LeftVector->operator[](0)->SetPosition(aPosition); 197 */ 198 //G4cout<<"Single hadron "<<LeftVector->operator[](0)->GetPosition()<<" "<<LeftVector->operator[](0)->GetFormationTime()<<G4endl; 199 } else { // 2 hadrons created from qq-qqbar are stored 200 LeftVector->operator[](0)->SetFormationTime(theString.GetTimeOfCreation()); 201 LeftVector->operator[](0)->SetPosition(theString.GetPosition()); 202 LeftVector->operator[](1)->SetFormationTime(theString.GetTimeOfCreation()); 203 LeftVector->operator[](1)->SetPosition(theString.GetPosition()); 204 } 205 // Uzhi insert 6.05.08 end 206 return LeftVector; 207 } 208 209 //--------------------- The string can fragment ------------------------------- 210 //--------------- At least two particles can be produced ---------------------- 211 212 LeftVector =new G4KineticTrackVector; 104 213 G4KineticTrackVector * RightVector=new G4KineticTrackVector; 105 214 106 // this should work but its only a semi deep copy. %GF G4ExcitedString theStringInCMS(theString);107 215 G4ExcitedString *theStringInCMS=CPExcited(theString); 108 216 G4LorentzRotation toCms=theStringInCMS->TransformToAlignedCms(); … … 111 219 G4int attempt=0; 112 220 while ( !success && attempt++ < StringLoopInterrupt ) 113 { 221 { // If the string fragmentation do not be happend, repeat the fragmentation--- 114 222 G4FragmentingString *currentString=new G4FragmentingString(*theStringInCMS); 115 223 116 //G4cout<<" MainFragmentString cur M2 "<<std::sqrt(currentString->Mass2())<<G4endl;117 118 std::for_each(LeftVector->begin() , LeftVector->end(), DeleteKineticTrack());224 //G4cout<<"FragmentString cur M2 "<<std::sqrt(currentString->Mass2())<<G4endl; 225 // Cleaning up the previously produced hadrons ------------------------------ 226 std::for_each(LeftVector->begin() , LeftVector->end() , DeleteKineticTrack()); 119 227 LeftVector->clear(); 228 120 229 std::for_each(RightVector->begin(), RightVector->end(), DeleteKineticTrack()); 121 230 RightVector->clear(); 122 231 232 // Main fragmentation loop until the string will not be able to fragment ---- 123 233 inner_sucess=true; // set false on failure.. 124 234 while (! StopFragmenting(currentString) ) 125 235 { // Split current string into hadron + new string 126 127 // G4FragmentingString *PreviousString=currentString; // Uzhi128 129 236 G4FragmentingString *newString=0; // used as output from SplitUp... 130 237 131 238 //G4cout<<"FragmentString to Splitup ===================================="<<G4endl; 239 //G4cout<<"++++++++++++++++++++++++++ Enter num--------------------------"<<G4endl; 132 240 //G4int Uzhi; G4cin>>Uzhi; // Uzhi 241 133 242 G4KineticTrack * Hadron=Splitup(currentString,newString); 134 243 //G4cout<<" Hadron "<<Hadron<<G4endl; 135 244 136 // if ( Hadron != 0 && IsFragmentable(newString)) // Uzhi 137 if ( Hadron != 0 ) // Uzhi 245 if ( Hadron != 0 ) // Store the hadron 138 246 { 139 247 if ( currentString->GetDecayDirection() > 0 ) … … 143 251 delete currentString; 144 252 currentString=newString; 145 } /* else { // Uzhi 146 // abandon ... start from the beginning 147 if (newString) delete newString; // ??? Uzhi local? 148 if (Hadron) delete Hadron; 149 // currentString = PreviousString; // Uzhi 150 inner_sucess=false; 151 break; 152 } */ // Uzhi 153 // delete PreviousString; // ??? Uzhi local? 253 } 154 254 }; 155 // Split current string into 2 final Hadrons255 // Split remaining string into 2 final Hadrons ------------------------ 156 256 //G4cout<<"FragmentString to SplitLast if inner_sucess#0"<<inner_sucess<<G4endl; 157 if ( inner_sucess && // Uzhi257 if ( inner_sucess && 158 258 SplitLast(currentString,LeftVector, RightVector) ) 159 259 { … … 161 261 } 162 262 delete currentString; 163 } 263 } // End of the loop in attemps to fragment the string 164 264 165 265 delete theStringInCMS; … … 188 288 G4LorentzRotation toObserverFrame(toCms.inverse()); 189 289 290 // LeftVector->operator[](0)->SetFormationTime(theString.GetTimeOfCreation()); 291 // LeftVector->operator[](0)->SetPosition(theString.GetPosition()); 292 293 G4double TimeOftheStringCreation=theString.GetTimeOfCreation(); 294 G4ThreeVector PositionOftheStringCreation(theString.GetPosition()); 295 296 /* // For large formation time open * 297 G4double TimeOftheStringCreation=theString.GetTimeOfCreation()+100*fermi; 298 G4ThreeVector PositionOftheStringCreation(theString.GetPosition().x(), 299 theString.GetPosition().y(), 300 theString.GetPosition().z()+100*fermi); 301 */ 302 303 /* 304 if(theString.GetPosition().y() > 100.*fermi){ 305 // It is a projectile-like string ------------------------------------- 306 G4double Zmin=theString.GetPosition().y()-1000.*fermi; 307 G4double Zmax=theString.GetPosition().z(); 308 TimeOftheStringCreation= 309 (Zmax-Zmin)*theString.Get4Momentum().e()/theString.Get4Momentum().z(); 310 311 G4ThreeVector aPosition(0.,0.,Zmax); 312 PositionOftheStringCreation=aPosition; 313 } 314 */ 190 315 for(size_t C1 = 0; C1 < LeftVector->size(); C1++) 191 316 { … … 194 319 Momentum = toObserverFrame*Momentum; 195 320 Hadron->Set4Momentum(Momentum); 321 196 322 G4LorentzVector Coordinate(Hadron->GetPosition(), Hadron->GetFormationTime()); 197 323 Momentum = toObserverFrame*Coordinate; 198 Hadron->SetFormationTime( Momentum.e());324 Hadron->SetFormationTime(TimeOftheStringCreation+Momentum.e()); 199 325 G4ThreeVector aPosition(Momentum.vect()); 200 Hadron->SetPosition(theString.GetPosition()+aPosition); 201 } 326 // Hadron->SetPosition(theString.GetPosition()+aPosition); 327 Hadron->SetPosition(PositionOftheStringCreation+aPosition); 328 //G4cout<<"Hadron "<<C1<<" "<<Hadron->GetPosition()/fermi<<" "<<Hadron->GetFormationTime()/fermi<<G4endl; 329 }; 202 330 203 331 //G4cout<<"Out FragmentString"<<G4endl; 204 332 return LeftVector; 205 333 206 207 208 334 } 209 335 336 //---------------------------------------------------------------------------------- 337 G4bool G4LundStringFragmentation::IsFragmentable(const G4FragmentingString * const string) 338 { 339 //G4cout<<"In IsFragmentable"<<G4endl; 340 SetMinimalStringMass(string); // Uzhi 341 //G4cout<<"Out IsFragmentable MinMass"<<MinimalStringMass<<" String Mass"<<std::sqrt(string->Get4Momentum().mag2())<<G4endl; 342 return sqr(MinimalStringMass + WminLUND) < string->Get4Momentum().mag2(); // Uzhi 343 344 // return sqr(FragmentationMass(string)+MassCut) < // Uzhi 345 // string->Mass2(); // Uzhi 346 } 347 348 //---------------------------------------------------------------------------------------- 349 G4bool G4LundStringFragmentation::StopFragmenting(const G4FragmentingString * const string) 350 { 351 //G4cout<<"StopFragmenting"<<G4endl; 352 353 SetMinimalStringMass(string); 354 355 //G4cout<<"StopFragm MinMass "<<MinimalStringMass<<" String Mass "<<std::sqrt(string->Get4Momentum().mag2())<<G4endl; 356 357 return (MinimalStringMass + WminLUND)* 358 (1 + SmoothParam * (1.-2*G4UniformRand())) > 359 string->Mass(); 360 } 361 362 //---------------------------------------------------------------------------------------- 363 G4bool G4LundStringFragmentation::SplitLast(G4FragmentingString * string, 364 G4KineticTrackVector * LeftVector, 365 G4KineticTrackVector * RightVector) 366 { 367 //... perform last cluster decay 368 /* 369 G4cout<<"SplitLast String mass "<<string->Mass()<<G4endl; 370 G4cout<<string->GetLeftParton()->GetPDGEncoding()<<" "<<G4endl; 371 G4cout<<string->GetRightParton()->GetPDGEncoding()<<" "<<G4endl; 372 */ 373 G4LorentzVector Str4Mom=string->Get4Momentum(); 374 //G4cout<<"String 4 momentum "<<Str4Mom<<G4endl; 375 G4ThreeVector ClusterVel =string->Get4Momentum().boostVector(); 376 377 G4double ResidualMass = string->Mass(); 378 379 G4ParticleDefinition * LeftHadron, * RightHadron; 380 381 G4int cClusterInterrupt = 0; 382 do 383 { 384 //G4cout<<" Cicle "<<cClusterInterrupt<<" "<< ClusterLoopInterrupt<<G4endl; 385 386 if (cClusterInterrupt++ >= ClusterLoopInterrupt) 387 { 388 return false; 389 } 390 G4ParticleDefinition * quark = NULL; 391 string->SetLeftPartonStable(); // to query quark contents.. 392 393 if (!string->FourQuarkString() ) 394 { 395 // The string is q-qbar, or q-qq, or qbar-qqbar type 396 if (string->DecayIsQuark() && string->StableIsQuark() ) 397 { 398 //... there are quarks on cluster ends 399 LeftHadron= QuarkSplitup(string->GetLeftParton(), quark); 400 } else 401 { 402 //... there is a Diquark on one of the cluster ends 403 G4int IsParticle; 404 405 if ( string->StableIsQuark() ) 406 { 407 IsParticle=(string->GetLeftParton()->GetPDGEncoding()>0) ? -1 : +1; 408 } else 409 { 410 IsParticle=(string->GetLeftParton()->GetPDGEncoding()>0) ? +1 : -1; 411 } 412 413 pDefPair QuarkPair = CreatePartonPair(IsParticle,false); // no diquarks wanted 414 quark = QuarkPair.second; 415 416 LeftHadron=hadronizer->Build(QuarkPair.first, string->GetLeftParton()); 417 } 418 419 RightHadron = hadronizer->Build(string->GetRightParton(), quark); 420 } else 421 { 422 // The string is qq-qqbar type. Diquarks are on the string ends 423 G4int LiftQuark1= string->GetLeftParton()->GetPDGEncoding()/1000; 424 G4int LiftQuark2=(string->GetLeftParton()->GetPDGEncoding()/100)%10; 425 426 G4int RightQuark1= string->GetRightParton()->GetPDGEncoding()/1000; 427 G4int RightQuark2=(string->GetRightParton()->GetPDGEncoding()/100)%10; 428 429 if(G4UniformRand()<0.5) 430 { 431 LeftHadron =hadronizer->Build(FindParticle( LiftQuark1), 432 FindParticle(RightQuark1)); 433 RightHadron=hadronizer->Build(FindParticle( LiftQuark2), 434 FindParticle(RightQuark2)); 435 } else 436 { 437 LeftHadron =hadronizer->Build(FindParticle( LiftQuark1), 438 FindParticle(RightQuark2)); 439 RightHadron=hadronizer->Build(FindParticle( LiftQuark2), 440 FindParticle(RightQuark1)); 441 } 442 } 443 /* 444 G4cout<<"SplitLast Left Right hadrons "<<LeftHadron->GetPDGEncoding()<<" "<<RightHadron->GetPDGEncoding()<<G4endl; 445 G4cout<<"SplitLast Left Right hadrons "<<LeftHadron->GetPDGMass()<<" "<<RightHadron->GetPDGMass()<<G4endl; 446 G4cout<<"Sum H mass Str Mass "<<LeftHadron->GetPDGMass() + RightHadron->GetPDGMass()<<" "<<ResidualMass<<G4endl; 447 */ 448 //... repeat procedure, if mass of cluster is too low to produce hadrons 449 //... ClusterMassCut = 0.15*GeV model parameter 450 451 } 452 while (ResidualMass <= LeftHadron->GetPDGMass() + RightHadron->GetPDGMass());// UzhiVOVA 453 454 //... compute hadron momenta and energies 455 456 G4LorentzVector LeftMom, RightMom; 457 G4ThreeVector Pos; 458 459 Sample4Momentum(&LeftMom, LeftHadron->GetPDGMass(), 460 &RightMom, RightHadron->GetPDGMass(), 461 ResidualMass); 462 463 LeftMom.boost(ClusterVel); 464 RightMom.boost(ClusterVel); 465 466 LeftVector->push_back(new G4KineticTrack(LeftHadron, 0, Pos, LeftMom)); 467 RightVector->push_back(new G4KineticTrack(RightHadron, 0, Pos, RightMom)); 468 469 //G4cout<<"Out SplitLast "<<G4endl; 470 return true; 471 472 } 473 210 474 //---------------------------------------------------------------------------------------------------------- 211 475 212 //G4double G4LundStringFragmentation::GetLightConeZ(G4double zmin, G4double zmax, // Uzhi 213 // G4int , G4ParticleDefinition* pHadron, // Uzhi 214 G4double G4LundStringFragmentation::GetLightConeZ(G4double zmin, G4double zmax, 215 G4int, G4ParticleDefinition* pHadron, // Uzhi 216 G4double Px, G4double Py) 476 void G4LundStringFragmentation::Sample4Momentum(G4LorentzVector* Mom, G4double Mass, G4LorentzVector* AntiMom, G4double AntiMass, G4double InitialMass) 477 { 478 // ------ Sampling of momenta of 2 last produced hadrons -------------------- 479 G4ThreeVector Pt; 480 G4double MassMt2, AntiMassMt2; 481 G4double AvailablePz, AvailablePz2; 482 483 //G4cout<<"Sample4Momentum "<<G4endl; 484 //G4cout<<"Sample4Momentum Mass"<<Mass<<" "<<AntiMass<<" "<<InitialMass<<G4endl; 485 486 if(Mass > 930. || AntiMass > 930.) // If there is a baryon 217 487 { 218 const G4double alund = 0.7/GeV/GeV; 219 220 // If blund get restored, you MUST adapt the calculation of zOfMaxyf. 221 // const G4double blund = 1; 222 223 G4double z, yf; 224 G4double Mass = pHadron->GetPDGMass(); 488 // ----------------- Isotripic decay ------------------------------------ 489 G4double r_val = sqr(InitialMass*InitialMass - Mass*Mass - AntiMass*AntiMass) - 490 sqr(2.*Mass*AntiMass); 491 G4double Pabs = (r_val > 0.)? std::sqrt(r_val)/(2.*InitialMass) : 0; 492 493 //... sample unit vector 494 G4double pz = 1. - 2.*G4UniformRand(); 495 G4double st = std::sqrt(1. - pz * pz)*Pabs; 496 G4double phi = 2.*pi*G4UniformRand(); 497 G4double px = st*std::cos(phi); 498 G4double py = st*std::sin(phi); 499 pz *= Pabs; 225 500 226 G4double Mt2 = Px*Px + Py*Py + Mass*Mass; 227 G4double zOfMaxyf=alund*Mt2/(alund*Mt2 + 1.); 228 G4double maxYf=(1-zOfMaxyf)/zOfMaxyf * std::exp(-alund*Mt2/zOfMaxyf); 229 230 // G4double N=1.; // Uzhi 231 // G4double OverN=1./N; // Uzhi 232 // G4double ZminN=std::pow(zmin,N); // Uzhi 233 // G4double ZmaxN=std::pow(zmax,N); // Uzhi 234 // G4double Brac=ZmaxN-ZminN; // Uzhi 235 236 //G4cout<<" ZminN ZmaxN Brac Code "<<ZminN<<" "<< ZmaxN<<" "<<Brac<<" "<<PartonEncoding<<G4endl; 237 238 // if(std::abs(PartonEncoding) < 1000) // Uzhi 239 { // Uzhi q or q-bar 240 //G4cout<<" quark "<<G4endl; // Vova 241 do // Uzhi 242 { 243 z = zmin + G4UniformRand()*(zmax-zmin); 244 // yf = std::pow(1. - z, blund)/z*std::exp(-alund*Mt2/z); 245 yf = (1-z)/z * std::exp(-alund*Mt2/z); 246 } 247 while (G4UniformRand()*maxYf > yf); 248 } // Uzhi 249 // else // Uzhi 250 // { // Uzhi qq or qq-bar 251 // //G4cout<<"Di-quark"<<G4endl; // Vova 252 // z = std::pow(Brac * G4UniformRand() + ZminN, OverN); // Uzhi 253 // }; // Uzhi 254 // 255 //G4cout<<" test z "<<std::pow(2.,3.)<<" "<<z<<G4endl; // Vova 256 return z; 501 Mom->setPx(px); Mom->setPy(py); Mom->setPz(pz); 502 Mom->setE(std::sqrt(Pabs*Pabs + Mass*Mass)); 503 504 AntiMom->setPx(-px); AntiMom->setPy(-py); AntiMom->setPz(-pz); 505 AntiMom->setE (std::sqrt(Pabs*Pabs + AntiMass*AntiMass)); 506 } 507 else 508 { 509 do 510 { 511 Pt=SampleQuarkPt(); Pt.setZ(0); G4double Pt2=Pt.mag2(); 512 513 //G4cout<<"Sample4Momentum Pt x y "<<Pt.getX()<<" "<<Pt.getY()<<G4endl; 514 515 MassMt2 = Mass * Mass + Pt2; 516 AntiMassMt2= AntiMass * AntiMass + Pt2; 517 518 //G4cout<<"Mts "<<MassMt2<<" "<<AntiMassMt2<<" "<<InitialMass*InitialMass<<G4endl; 519 520 AvailablePz2= sqr(InitialMass*InitialMass - MassMt2 - AntiMassMt2) - 521 4.*MassMt2*AntiMassMt2; 522 } 523 while(AvailablePz2 < 0.); 524 525 AvailablePz2 /=(4.*InitialMass*InitialMass); 526 AvailablePz = std::sqrt(AvailablePz2); 527 528 //G4cout<<"AvailablePz "<<AvailablePz<<G4endl; 529 530 G4double Px=Pt.getX(); 531 G4double Py=Pt.getY(); 532 533 //if(Mass > AntiMass){AvailablePz=-AvailablePz;} // May30 // Uzhi 534 535 Mom->setPx(Px); Mom->setPy(Py); Mom->setPz(AvailablePz); 536 Mom->setE(std::sqrt(MassMt2+AvailablePz2)); 537 538 //G4cout<<" 1 part "<<Px<<" "<<Py<<" "<<AvailablePz<<" "<<std::sqrt(MassMt2+AvailablePz2)<<G4endl; 539 540 AntiMom->setPx(-Px); AntiMom->setPy(-Py); AntiMom->setPz(-AvailablePz); 541 AntiMom->setE (std::sqrt(AntiMassMt2+AvailablePz2)); 542 543 //G4cout<<" 2 part "<<-Px<<" "<<-Py<<" "<<-AvailablePz<<" "<<std::sqrt(AntiMassMt2+AvailablePz2)<<G4endl; 257 544 } 258 //----------------------------------------------------------------------------------------- 545 546 //G4cout<<"Out Sample4Momentum "<<G4endl; 547 } 548 549 //----------------------------------------------------------------------------- 259 550 260 551 G4LorentzVector * G4LundStringFragmentation::SplitEandP(G4ParticleDefinition * pHadron, 261 G4FragmentingString * string) 262 { 552 G4FragmentingString * string, G4FragmentingString * newString) 553 { 554 /* 555 G4cout<<"SplitEandP "<<G4endl; 556 G4cout<<"SplitEandP string mass "<<string->Mass()<<G4endl; 557 G4cout<<string->GetLeftParton()->GetPDGEncoding()<<" " 558 <<string->GetRightParton()->GetPDGEncoding()<<" "<<G4endl; 559 G4cout<<G4endl; 560 G4cout<<newString->GetLeftParton()->GetPDGEncoding()<<" "<<G4endl; 561 G4cout<<newString->GetRightParton()->GetPDGEncoding()<<" "<<G4endl; 562 */ 563 G4LorentzVector String4Momentum=string->Get4Momentum(); 564 G4double StringMT2=string->Get4Momentum().mt2(); 565 //G4cout<<"StringMt2 "<<StringMT2<<G4endl; 566 263 567 G4double HadronMass = pHadron->GetPDGMass(); 264 SetMinimalStringMass(string); // Uzhi 265 G4double StringMass2 = string->Mass2(); // Uzhi 266 267 //G4cout<<"SplitEandP string mass "<<string->Mass()<<" Hadron mass "<<HadronMass<<pHadron->GetParticleName()<<G4endl; // Uzhi 268 //G4cout<<string->GetLeftParton()->GetPDGEncoding()<<" "<<G4endl;269 //G4cout<< string->GetRightParton()->GetPDGEncoding()<<" "<<G4endl;270 //G4cout<<" 271 272 // calculate and assign hadron transverse momentum component HadronPx andHadronPy568 //G4cout<<"Hadron mass "<<HadronMass<<" "<<pHadron->GetParticleName()<<G4endl; 569 570 SetMinimalStringMass(newString); 571 String4Momentum.setPz(0.); 572 G4ThreeVector StringPt=String4Momentum.vect(); 573 //G4cout<<"StringPt "<<StringPt<<G4endl<<G4endl; 574 //G4cout<<"Min string mass "<<MinimalStringMass<<G4endl; 575 576 // calculate and assign hadron transverse momentum component HadronPx and HadronPy 273 577 G4ThreeVector thePt; 274 578 thePt=SampleQuarkPt(); … … 276 580 G4ThreeVector HadronPt = thePt +string->DecayPt(); 277 581 HadronPt.setZ(0); 582 //G4cout<<"Hadron Pt"<<HadronPt<<G4endl; 583 584 G4ThreeVector RemSysPt = StringPt - HadronPt; 585 //G4cout<<"RemSys Pt"<<RemSysPt<<G4endl; 586 278 587 //... sample z to define hadron longitudinal momentum and energy 279 588 //... but first check the available phase space 280 589 281 // G4double DecayQuarkMass2 = sqr(string->GetDecayParton()->GetPDGMass()); 282 283 //G4cout<<" QuarkMass "<<string->GetDecayParton()->GetPDGMass()<<G4endl; // Uzhi 284 285 G4double HadronMass2T = sqr(HadronMass) + HadronPt.mag2(); 286 // G4double ResidualMass2T=sqr(MinimalStringMass + WminLUND) + HadronPt.mag2(); // Uzhi 287 G4double ResidualMass2T=sqr(MinimalStringMass + WminLUND) + HadronPt.mag2(); // Uzhi 288 289 //G4cout<<" Mt h res str "<<std::sqrt(HadronMass2T)<<" "<<std::sqrt(ResidualMass2T)<<" srt mass"<<string->Mass()<<G4endl; 290 291 // if (DecayQuarkMass2 + HadronMass2T >= SmoothParam*(string->Mass2()) ) // Uzhi 292 293 G4double Pz2 = (sqr(StringMass2 - HadronMass2T - ResidualMass2T) - // Uzhi 294 4*HadronMass2T * ResidualMass2T)/4./StringMass2; // Uzhi 295 296 //G4cout<<" Pz**2 "<<Pz2<<G4endl; 297 298 if(Pz2 < 0 ) {return 0;} // have to start all over! // Uzhi 590 G4double HadronMassT2 = sqr(HadronMass) + HadronPt.mag2(); 591 G4double ResidualMassT2=sqr(MinimalStringMass) + RemSysPt.mag2(); 592 593 //G4cout<<"Mt h res str "<<std::sqrt(HadronMassT2)<<" "<<std::sqrt(ResidualMassT2)<<" srt mass"<<StringMT2<<G4endl; 594 595 G4double Pz2 = (sqr(StringMT2 - HadronMassT2 - ResidualMassT2) - 596 4*HadronMassT2 * ResidualMassT2)/4./StringMT2; 597 598 //G4cout<<"Pz**2 "<<Pz2<<G4endl; 599 600 if(Pz2 < 0 ) {return 0;} // have to start all over! 299 601 300 602 //... then compute allowed z region z_min <= z <= z_max 301 603 302 G4double Pz = std::sqrt(Pz2); // Uzhi 303 G4double zMin = (std::sqrt(HadronMass2T+Pz2) - Pz)/std::sqrt(StringMass2); // Uzhi 304 G4double zMax = (std::sqrt(HadronMass2T+Pz2) + Pz)/std::sqrt(StringMass2); // Uzhi 305 306 //G4cout<<" Zmin max "<<zMin<<" "<<zMax<<G4endl; // Uzhi 307 308 // G4double zMax = 1. - DecayQuarkMass2/(string->Mass2()); // Uzhi 604 G4double Pz = std::sqrt(Pz2); 605 G4double zMin = (std::sqrt(HadronMassT2+Pz2) - Pz)/std::sqrt(StringMT2); 606 G4double zMax = (std::sqrt(HadronMassT2+Pz2) + Pz)/std::sqrt(StringMT2); 607 608 //G4cout<<"Zmin max "<<zMin<<" "<<zMax<<G4endl; // Uzhi 609 309 610 if (zMin >= zMax) return 0; // have to start all over! 310 611 … … 318 619 HadronPt.setZ(0.5* string->GetDecayDirection() * 319 620 (z * string->LightConeDecay() - 320 HadronMass 2T/(z * string->LightConeDecay())));621 HadronMassT2/(z * string->LightConeDecay()))); 321 622 322 623 G4double HadronE = 0.5* (z * string->LightConeDecay() + 323 HadronMass 2T/(z * string->LightConeDecay()));624 HadronMassT2/(z * string->LightConeDecay())); 324 625 325 626 G4LorentzVector * a4Momentum= new G4LorentzVector(HadronPt,HadronE); 326 627 327 //G4cout<<"Out of SplitEandP Pz E "<<HadronPt.getZ()<<" "<<0.5* (z * string->LightConeDecay() + HadronMass2T/(z * string->LightConeDecay()))<<G4endl; 628 //G4cout<<"Hadron Pt"<<HadronPt<<G4endl; 629 //G4cout<<"Out of SplitEandP Pz E "<<HadronPt.getZ()<<" "<<HadronE<<G4endl; 328 630 329 631 return a4Momentum; 330 632 } 331 633 332 333 634 //----------------------------------------------------------------------------------------- 334 335 G4bool G4LundStringFragmentation::SplitLast(G4FragmentingString * string,336 G4KineticTrackVector * LeftVector,337 G4KineticTrackVector * RightVector)635 G4double G4LundStringFragmentation::GetLightConeZ(G4double zmin, G4double zmax, 636 G4int PDGEncodingOfDecayParton, 637 G4ParticleDefinition* pHadron, 638 G4double Px, G4double Py) 338 639 { 339 //... perform last cluster decay 340 //G4cout<<"SplitLast String mass "<<string->Mass()<<G4endl; 341 //G4cout<<string->GetLeftParton()->GetPDGEncoding()<<" "<<G4endl; 342 //G4cout<<string->GetRightParton()->GetPDGEncoding()<<" "<<G4endl; 343 344 G4ThreeVector ClusterVel =string->Get4Momentum().boostVector(); 345 346 G4double ResidualMass = string->Mass(); 347 // G4double ClusterMassCut = ClusterMass; 348 349 G4int cClusterInterrupt = 0; 350 351 G4ParticleDefinition * LeftHadron, * RightHadron; 352 do 353 { 354 //G4cout<<" Cicle "<<cClusterInterrupt<<" "<< ClusterLoopInterrupt<<G4endl; 355 356 if (cClusterInterrupt++ >= ClusterLoopInterrupt) 357 { 358 return false; 359 } 360 G4ParticleDefinition * quark = NULL; 361 string->SetLeftPartonStable(); // to query quark contents.. 362 363 if (string->DecayIsQuark() && string->StableIsQuark() ) 364 { 365 //... there are quarks on cluster ends 366 LeftHadron= QuarkSplitup(string->GetLeftParton(), quark); 367 } else { 368 //... there is a Diquark on cluster ends 369 G4int IsParticle; 370 371 if ( string->StableIsQuark() ) { 372 IsParticle=(string->GetLeftParton()->GetPDGEncoding()>0) ? -1 : +1; 373 } else { 374 IsParticle=(string->GetLeftParton()->GetPDGEncoding()>0) ? +1 : -1; 375 } 376 377 pDefPair QuarkPair = CreatePartonPair(IsParticle,false); // no diquarks wanted 378 quark = QuarkPair.second; 379 380 LeftHadron=hadronizer->Build(QuarkPair.first, string->GetLeftParton()); 381 } 382 383 RightHadron = hadronizer->Build(string->GetRightParton(), quark); 384 385 //G4cout<<"SplitLast Left Right hadrons "<<LeftHadron->GetPDGEncoding()<<" "<<RightHadron->GetPDGEncoding()<<G4endl; 386 //G4cout<<"SplitLast Left Right hadrons "<<LeftHadron->GetPDGMass()<<" "<<RightHadron->GetPDGMass()<<G4endl; 387 //G4cout<<" Sum H mass Str Mass "<<LeftHadron->GetPDGMass() + RightHadron->GetPDGMass()<<" "<<ResidualMass<<G4endl; 388 389 //... repeat procedure, if mass of cluster is too low to produce hadrons 390 //... ClusterMassCut = 0.15*GeV model parameter 391 // if ( quark->GetParticleSubType()== "quark" ) {ClusterMassCut = 0.;} // Uzhi 392 // else {ClusterMassCut = ClusterMass;} // Uzhi 393 394 } 395 while (ResidualMass <= LeftHadron->GetPDGMass() + RightHadron->GetPDGMass()); // Uzhi VOVA 396 // while (ResidualMass <= LeftHadron->GetPDGMass() + RightHadron->GetPDGMass() + ClusterMassCut); // Uzhi 397 //... compute hadron momenta and energies 398 399 G4LorentzVector LeftMom, RightMom; 400 G4ThreeVector Pos; 401 402 //G4cout<<"Sample4Momentum"<<G4endl; 403 404 Sample4Momentum(&LeftMom, LeftHadron->GetPDGMass(), &RightMom, RightHadron->GetPDGMass(), ResidualMass); 405 406 LeftMom.boost(ClusterVel); 407 RightMom.boost(ClusterVel); 408 409 LeftVector->push_back(new G4KineticTrack(LeftHadron, 0, Pos, LeftMom)); 410 RightVector->push_back(new G4KineticTrack(RightHadron, 0, Pos, RightMom)); 411 412 return true; 413 414 } 415 //---------------------------------------------------------------------------------------------------------- 416 417 G4bool G4LundStringFragmentation::IsFragmentable(const G4FragmentingString * const string) 418 { 419 //G4cout<<"In IsFragmentable"<<G4endl; 420 SetMinimalStringMass(string); // Uzhi 421 //G4cout<<"Out IsFragmentable MinMass"<<MinimalStringMass<<" String Mass"<<std::sqrt(string->Get4Momentum().mag2())<<G4endl; 422 return sqr(MinimalStringMass + WminLUND) < string->Get4Momentum().mag2(); // Uzhi 423 424 // return sqr(FragmentationMass(string)+MassCut) < // Uzhi 425 // string->Mass2(); // Uzhi 426 } 427 428 //---------------------------------------------------------------------------------------------------------- 429 430 G4bool G4LundStringFragmentation::StopFragmenting(const G4FragmentingString * const string) 431 { 432 //G4cout<<"StopFragmenting"<<G4endl; 433 434 SetMinimalStringMass(string); // Uzhi 435 //G4cout<<"StopFragm MinMass "<<MinimalStringMass<<" String Mass "<<std::sqrt(string->Get4Momentum().mag2())<<G4endl; 436 return sqr((MinimalStringMass + WminLUND)*(1 + SmoothParam * (1.-2*G4UniformRand()))) > // Uzhi 437 string->Get4Momentum().mag2(); // Uzhi 438 // sqr(FragmentationMass(string,&G4HadronBuilder::BuildHighSpin)+MassCut) > // Uzhi 439 // string->Get4Momentum().mag2(); // Uzhi 440 } 441 442 //---------------------------------------------------------------------------------------------------------- 443 444 void G4LundStringFragmentation::Sample4Momentum(G4LorentzVector* Mom, G4double Mass, G4LorentzVector* AntiMom, G4double AntiMass, G4double InitialMass) 445 { 446 G4ThreeVector Pt; // Uzhi 447 G4double MassMt2, AntiMassMt2; // Uzhi 448 G4double AvailablePz, AvailablePz2; // Uzhi 449 450 //G4cout<<" Smpl4Mom "<<Mass<<" "<<AntiMass<<" "<<InitialMass<<G4endl; 451 // Uzhi 452 do // Uzhi 453 { // Uzhi 454 Pt=SampleQuarkPt(); Pt.setZ(0); G4double Pt2=Pt.mag2(); // Uzhi 455 456 //G4cout<<"Sample4Momentum Pt x y "<<Pt.getX()<<" "<<Pt.getY()<<G4endl; 457 458 MassMt2 = Mass * Mass + Pt2; // Uzhi 459 AntiMassMt2= AntiMass * AntiMass + Pt2; // Uzhi 460 461 //G4cout<<"Mts "<<MassMt2<<" "<<AntiMassMt2<<" "<<InitialMass*InitialMass<<G4endl; 462 463 AvailablePz2= sqr(InitialMass*InitialMass - MassMt2 - AntiMassMt2) - 464 4.*MassMt2*AntiMassMt2; // Uzhi 465 } // Uzhi 466 while(AvailablePz2 < 0.); // Uzhi 467 // Uzhi 468 AvailablePz2 /=(4.*InitialMass*InitialMass); // Uzhi 469 // Uzhi 470 AvailablePz = std::sqrt(AvailablePz2); // Uzhi 471 472 //G4cout<<"AvailablePz "<<AvailablePz<<G4endl; 473 474 475 G4double Px=Pt.getX(); // Uzhi 476 G4double Py=Pt.getY(); // Uzhi 477 // Uzhi 478 Mom->setPx(Px); Mom->setPy(Py); Mom->setPz(AvailablePz); // Uzhi 479 Mom->setE(std::sqrt(MassMt2+AvailablePz2)); // Uzhi 480 481 //G4cout<<" 1 part "<<Px<<" "<<Py<<" "<<AvailablePz<<" "<<std::sqrt(MassMt2+AvailablePz2)<<G4endl; 482 483 // Uzhi 484 AntiMom->setPx(-Px); AntiMom->setPy(-Py); AntiMom->setPz(-AvailablePz); // Uzhi 485 AntiMom->setE (std::sqrt(AntiMassMt2+AvailablePz2)); // Uzhi 486 487 //G4cout<<" 2 part "<<-Px<<" "<<-Py<<" "<<-AvailablePz<<" "<<std::sqrt(AntiMassMt2+AvailablePz2)<<G4endl; 488 489 // Maybe it must be inversed! // Uzhi 490 /* // Uzhi 491 G4double r_val = sqr(InitialMass*InitialMass - Mass*Mass - AntiMass*AntiMass) - sqr(2.*Mass*AntiMass); 492 G4double Pabs = (r_val > 0.)? std::sqrt(r_val)/(2.*InitialMass) : 0; 493 494 //... sample unit vector 495 G4double pz = 1. - 2.*G4UniformRand(); 496 G4double st = std::sqrt(1. - pz * pz)*Pabs; 497 G4double phi = 2.*pi*G4UniformRand(); 498 G4double px = st*std::cos(phi); 499 G4double py = st*std::sin(phi); 500 pz *= Pabs; 640 G4double alund; 641 642 // If blund get restored, you MUST adapt the calculation of zOfMaxyf. 643 // const G4double blund = 1; 644 645 G4double z, yf; 646 G4double Mass = pHadron->GetPDGMass(); 647 // G4int HadronEncoding=pHadron->GetPDGEncoding(); 501 648 502 Mom->setPx(px); Mom->setPy(py); Mom->setPz(pz); 503 Mom->setE(std::sqrt(Pabs*Pabs + Mass*Mass)); 504 505 AntiMom->setPx(-px); AntiMom->setPy(-py); AntiMom->setPz(-pz); 506 AntiMom->setE (std::sqrt(Pabs*Pabs + AntiMass*AntiMass)); 507 */ // Uzhi 649 G4double Mt2 = Px*Px + Py*Py + Mass*Mass; 650 651 if(std::abs(PDGEncodingOfDecayParton) < 1000) 652 { 653 // ---------------- Quark fragmentation ---------------------- 654 alund=0.35/GeV/GeV; // Instead of 0.7 because kinks are not considered 655 G4double zOfMaxyf=alund*Mt2/(alund*Mt2 + 1.); 656 G4double maxYf=(1-zOfMaxyf)/zOfMaxyf * std::exp(-alund*Mt2/zOfMaxyf); 657 658 do 659 { 660 z = zmin + G4UniformRand()*(zmax-zmin); 661 // yf = std::pow(1. - z, blund)/z*std::exp(-alund*Mt2/z); 662 yf = (1-z)/z * std::exp(-alund*Mt2/z); 663 } 664 while (G4UniformRand()*maxYf > yf); 665 } 666 else 667 { 668 // ---------------- Di-quark fragmentation ---------------------- 669 //G4cout<<"Di-quark"<<G4endl; // Vova 670 alund=0.7/GeV/GeV; // 0.7 2.0 671 G4double zOfMaxyf=alund*Mt2/(alund*Mt2 + 1.); 672 G4double maxYf=(1-zOfMaxyf)/zOfMaxyf * std::exp(-alund*Mt2/zOfMaxyf); 673 do 674 { 675 z = zmin + G4UniformRand()*(zmax-zmin); 676 // yf = std::pow(1. - z, blund)/z*std::exp(-alund*Mt2/z); 677 yf = (1-z)/z * std::exp(-alund*Mt2/z); 678 } 679 while (G4UniformRand()*maxYf > yf); 680 }; 681 682 //G4cout<<" test z "<<std::pow(2.,3.)<<" "<<z<<G4endl; 683 return z; 508 684 } 509 510 511 void G4LundStringFragmentation::SetMinimalStringMass(const G4FragmentingString * const string) // Uzhi512 {513 //G4cout<<"In SetMinMass -------------------"<<std::sqrt(string->Mass2())<<G4endl;514 //G4cout<<string->GetLeftParton()->GetPDGEncoding()<<G4endl;515 //G4cout<<string->GetRightParton()->GetPDGEncoding()<<G4endl;516 517 G4double EstimatedMass=0.750* GeV; // 2*m_q518 519 G4int Qleft =std::abs(string->GetLeftParton()->GetPDGEncoding());520 521 if( Qleft > 1000)522 {523 G4int q1=Qleft/1000;524 if( q1 < 3) {EstimatedMass += 0.325* GeV;}525 if( q1 > 2) {EstimatedMass += 0.500* GeV;}526 527 G4int q2=(Qleft/100)%10;528 if( q2 < 3) {EstimatedMass += 0.325* GeV;}529 if( q2 > 2) {EstimatedMass += 0.500* GeV;}530 }531 else532 {533 if( Qleft < 3) {EstimatedMass += 0.325* GeV;}534 if( Qleft > 2) {EstimatedMass += 0.500* GeV;}535 }536 537 G4int Qright=std::abs(string->GetRightParton()->GetPDGEncoding());538 539 if( Qright > 1000)540 {541 G4int q1=Qright/1000;542 if( q1 < 3) {EstimatedMass += 0.325* GeV;}543 if( q1 > 2) {EstimatedMass += 0.500* GeV;}544 545 G4int q2=(Qright/100)%10;546 if( q2 < 3) {EstimatedMass += 0.325* GeV;}547 if( q2 > 2) {EstimatedMass += 0.500* GeV;}548 }549 else550 {551 if( Qright < 3) {EstimatedMass += 0.325* GeV;}552 if( Qright > 2) {EstimatedMass += 0.500* GeV;}553 }554 555 MinimalStringMass=EstimatedMass;556 SetMinimalStringMass2(EstimatedMass);557 558 /*559 Pcreate build=&G4HadronBuilder::BuildLowSpin;560 561 G4ParticleDefinition *Hadron1, *Hadron2=0;562 563 G4int iflc = (G4UniformRand() < 0.5)? 1 : 2;564 565 if (string->GetLeftParton()->GetParticleSubType() == "quark") iflc = -iflc;566 567 if (string->GetLeftParton()->GetPDGEncoding() < 0) iflc = -iflc;568 569 // 1/2 baryon (anti-baryon) and scalar meson (QQ-q or QbarQbar-Qbar),570 // or 2 scalar mesons (Q-Qbar),571 // or 2 1/2 baryons (anti-baryons) will be built (QQ-QbarQbar)572 573 //G4cout<<"In SetMinMass -------------------"<<std::sqrt(string->Mass2())<<G4endl;574 //G4cout<<string->GetLeftParton()->GetPDGEncoding()<<" "<<FindParticle(iflc)->GetPDGEncoding()<<G4endl;575 //G4cout<<string->GetRightParton()->GetPDGEncoding()<<" "<<FindParticle(-iflc)->GetPDGEncoding()<<G4endl;576 577 Hadron1 = (hadronizer->*build)(string->GetLeftParton(),FindParticle(iflc));578 Hadron2 =(hadronizer->*build)(string->GetRightParton(),FindParticle(-iflc));579 MinimalStringMass = (Hadron1)->GetPDGMass() + (Hadron2)->GetPDGMass();580 581 //G4cout<<(Hadron1)->GetPDGEncoding()<<" "<<(Hadron2)->GetPDGEncoding()<<G4endl;582 //G4cout<<"Out SetMinMass "<<MinimalStringMass<<G4endl;583 */584 // SetMinimalStringMass2(MinimalStringMass);585 }586 //*******************************************************************************************************587 588 void G4LundStringFragmentation::SetMinimalStringMass2(const G4double aValue) // Uzhi589 {590 MinimalStringMass2=aValue * aValue;591 }592 //*******************************************************************************************************593 594 //**************************************************************************************** -
trunk/source/processes/hadronic/models/parton_string/hadronization/src/G4QGSMFragmentation.cc
r819 r962 25 25 // 26 26 // 27 // $Id: G4QGSMFragmentation.cc,v 1. 6 2007/04/24 14:55:23 gunterExp $28 // GEANT4 tag $Name: geant4-09-0 1-patch-02 $27 // $Id: G4QGSMFragmentation.cc,v 1.9 2008/06/23 08:35:55 vuzhinsk Exp $ 28 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 29 29 // 30 30 // ----------------------------------------------------------------------------- … … 121 121 } else { 122 122 // abandon ... start from the beginning 123 if (newString) delete newString; 123 if (newString) delete newString; // Uzhi restore 20.06.08 124 124 if (Hadron) delete Hadron; 125 125 inner_sucess=false; … … 229 229 230 230 G4LorentzVector * G4QGSMFragmentation::SplitEandP(G4ParticleDefinition * pHadron, 231 G4FragmentingString * string) 231 G4FragmentingString * string, // Uzhi 232 G4FragmentingString * ) // Uzhi 232 233 { 233 234 G4double HadronMass = pHadron->GetPDGMass(); -
trunk/source/processes/hadronic/models/parton_string/hadronization/src/G4VKinkyStringDecay.cc
r819 r962 25 25 // 26 26 // 27 // $Id: G4VKinkyStringDecay.cc,v 1. 3 2006/06/29 20:55:07 gunterExp $28 // GEANT4 tag $Name: $27 // $Id: G4VKinkyStringDecay.cc,v 1.4 2008/04/25 14:20:14 vuzhinsk Exp $ 28 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 29 29 // Maxim Komogorov 30 30 // -
trunk/source/processes/hadronic/models/parton_string/hadronization/src/G4VLongitudinalStringDecay.cc
r819 r962 25 25 // 26 26 // 27 // $Id: G4VLongitudinalStringDecay.cc,v 1. 8 2007/04/24 14:55:23 gunterExp $28 // GEANT4 tag $Name: geant4-09-0 1-patch-02 $27 // $Id: G4VLongitudinalStringDecay.cc,v 1.13 2008/06/23 08:35:55 vuzhinsk Exp $ 28 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 29 29 // 30 30 // ----------------------------------------------------------------------------- … … 75 75 76 76 StrangeSuppress = 0.44; // 27 % strange quarks produced, ie. u:d:s=1:1:0.27 77 DiquarkSuppress = 0. 1;77 DiquarkSuppress = 0.07; 78 78 DiquarkBreakProb = 0.1; 79 79 … … 106 106 hadronizer = new G4HadronBuilder(pspin_meson,pspin_barion, 107 107 scalarMesonMix,vectorMesonMix); 108 Kappa = 1.0 * GeV/fermi; 109 110 108 111 } 109 112 … … 114 117 } 115 118 116 //============================================================================= ================-------------119 //============================================================================= 117 120 118 121 // Operators … … 122 125 // } 123 126 124 //----------------------------------------------------------------------------- -----------------------------127 //----------------------------------------------------------------------------- 125 128 126 129 int G4VLongitudinalStringDecay::operator==(const G4VLongitudinalStringDecay &) const … … 130 133 } 131 134 132 //------------------------------------------------------------------------------------- ---------------------135 //------------------------------------------------------------------------------------- 133 136 134 137 int G4VLongitudinalStringDecay::operator!=(const G4VLongitudinalStringDecay &) const … … 138 141 } 139 142 140 //========================================================================================================== 143 //*********************************************************************************** 144 145 // For changing Mass Cut used for selection of very small mass strings 146 void G4VLongitudinalStringDecay::SetMassCut(G4double aValue){MassCut=aValue;} 147 148 //----------------------------------------------------------------------------- 149 150 // For handling a string with very low mass 151 152 G4KineticTrackVector* G4VLongitudinalStringDecay::LightFragmentationTest(const 153 G4ExcitedString * const string) 154 { 155 // Check string decay threshold 156 157 G4KineticTrackVector * result=0; // return 0 when string exceeds the mass cut 158 159 pDefPair hadrons((G4ParticleDefinition *)0,(G4ParticleDefinition *)0); 160 161 G4FragmentingString aString(*string); 162 163 if ( sqr(FragmentationMass(&aString,0,&hadrons)+MassCut) < aString.Mass2()) { 164 return 0; 165 } 166 167 // The string mass is very low --------------------------- 168 169 result=new G4KineticTrackVector; 170 171 if ( hadrons.second ==0 ) 172 { 173 // Substitute string by light hadron, Note that Energy is not conserved here! 174 175 /* 176 #ifdef DEBUG_LightFragmentationTest 177 G4cout << "VlongSF Warning replacing string by single hadron " <<G4endl; 178 G4cout << hadrons.first->GetParticleName() 179 << "string .. " << string->Get4Momentum() << " " 180 << string->Get4Momentum().m() << G4endl; 181 #endif 182 G4cout << "VlongSF Warning replacing string by single hadron " <<G4endl; 183 G4cout << hadrons.first->GetParticleName() << " string .. " <<G4endl; 184 G4cout << string->Get4Momentum() << " " << string->Get4Momentum().m() << G4endl; 185 */ 186 G4ThreeVector Mom3 = string->Get4Momentum().vect(); 187 G4LorentzVector Mom(Mom3, 188 std::sqrt(Mom3.mag2() + 189 sqr(hadrons.first->GetPDGMass()))); 190 result->push_back(new G4KineticTrack(hadrons.first, 0, 191 string->GetPosition(), 192 Mom)); 193 } else 194 { 195 //... string was qq--qqbar type: Build two stable hadrons, 196 197 #ifdef DEBUG_LightFragmentationTest 198 G4cout << "VlongSF Warning replacing qq-qqbar string by TWO hadrons " 199 << hadrons.first->GetParticleName() << " / " 200 << hadrons.second->GetParticleName() 201 << "string .. " << string->Get4Momentum() << " " 202 << string->Get4Momentum().m() << G4endl; 203 #endif 204 205 // Uzhi Formation time in the case??? 206 207 G4LorentzVector Mom1, Mom2; 208 Sample4Momentum(&Mom1, hadrons.first->GetPDGMass(), 209 &Mom2,hadrons.second->GetPDGMass(), 210 string->Get4Momentum().mag()); 211 212 result->push_back(new G4KineticTrack(hadrons.first, 0, 213 string->GetPosition(), 214 Mom1)); 215 result->push_back(new G4KineticTrack(hadrons.second, 0, 216 string->GetPosition(), 217 Mom2)); 218 219 G4ThreeVector Velocity = string->Get4Momentum().boostVector(); 220 result->Boost(Velocity); 221 } 222 223 return result; 224 225 } 226 227 //---------------------------------------------------------------------------------------- 228 229 G4double G4VLongitudinalStringDecay::FragmentationMass( 230 const G4FragmentingString * const string, 231 Pcreate build, pDefPair * pdefs ) 232 { 233 234 G4double mass; 235 static G4bool NeedInit(true); 236 static std::vector<double> nomix; 237 static G4HadronBuilder * minMassHadronizer; 238 if ( NeedInit ) 239 { 240 NeedInit = false; 241 nomix.resize(6); 242 for ( G4int i=0; i<6 ; i++ ) nomix[i]=0; 243 244 // minMassHadronizer=new G4HadronBuilder(pspin_meson,pspin_barion,nomix,nomix); 245 minMassHadronizer=hadronizer; 246 } 247 248 if ( build==0 ) build=&G4HadronBuilder::BuildLowSpin; 249 250 G4ParticleDefinition *Hadron1, *Hadron2=0; 251 252 if (!string->FourQuarkString() ) 253 { 254 // spin 0 meson or spin 1/2 barion will be built 255 256 Hadron1 = (minMassHadronizer->*build)(string->GetLeftParton(), 257 string->GetRightParton()); 258 mass= (Hadron1)->GetPDGMass(); 259 } else 260 { 261 //... string is qq--qqbar: Build two stable hadrons, 262 //... with extra uubar or ddbar quark pair 263 G4int iflc = (G4UniformRand() < 0.5)? 1 : 2; 264 if (string->GetLeftParton()->GetPDGEncoding() < 0) iflc = -iflc; 265 266 //... theSpin = 4; spin 3/2 baryons will be built 267 Hadron1 = (minMassHadronizer->*build)(string->GetLeftParton(), 268 FindParticle(iflc) ); 269 Hadron2 = (minMassHadronizer->*build)(string->GetRightParton(), 270 FindParticle(-iflc) ); 271 mass = (Hadron1)->GetPDGMass() + (Hadron2)->GetPDGMass(); 272 } 273 274 if ( pdefs != 0 ) 275 { // need to return hadrons as well.... 276 pdefs->first = Hadron1; 277 pdefs->second = Hadron2; 278 } 279 280 return mass; 281 } 282 283 //---------------------------------------------------------------------------- 284 285 G4ParticleDefinition* G4VLongitudinalStringDecay::FindParticle(G4int Encoding) 286 { 287 G4ParticleDefinition* ptr = G4ParticleTable::GetParticleTable()->FindParticle(Encoding); 288 if (ptr == NULL) 289 { 290 G4cout << "Particle with encoding "<<Encoding<<" does not exist!!!"<<G4endl; 291 throw G4HadronicException(__FILE__, __LINE__, "Check your particle table"); 292 } 293 return ptr; 294 } 295 296 //----------------------------------------------------------------------------- 297 // virtual void Sample4Momentum(G4LorentzVector* Mom, G4double Mass, 298 // G4LorentzVector* AntiMom, G4double AntiMass, 299 // G4double InitialMass)=0; 300 //----------------------------------------------------------------------------- 301 302 //********************************************************************************* 303 // For decision on continue or stop string fragmentation 304 // virtual G4bool StopFragmenting(const G4FragmentingString * const string)=0; 305 // virtual G4bool IsFragmentable(const G4FragmentingString * const string)=0; 306 307 // If a string can not fragment, make last break into 2 hadrons 308 // virtual G4bool SplitLast(G4FragmentingString * string, 309 // G4KineticTrackVector * LeftVector, 310 // G4KineticTrackVector * RightVector)=0; 311 //----------------------------------------------------------------------------- 312 // 313 // If a string fragments, do the following 314 // 315 // For transver of a string to its CMS frame 316 //----------------------------------------------------------------------------- 317 318 G4ExcitedString *G4VLongitudinalStringDecay::CPExcited(const G4ExcitedString & in) 319 { 320 G4Parton *Left=new G4Parton(*in.GetLeftParton()); 321 G4Parton *Right=new G4Parton(*in.GetRightParton()); 322 return new G4ExcitedString(Left,Right,in.GetDirection()); 323 } 324 325 //----------------------------------------------------------------------------- 326 327 G4KineticTrack * G4VLongitudinalStringDecay::Splitup( 328 G4FragmentingString *string, 329 G4FragmentingString *&newString) 330 { 331 //G4cout<<"In G4VLong String Dec######################"<<G4endl; 332 333 //... random choice of string end to use for creating the hadron (decay) 334 SideOfDecay = (G4UniformRand() < 0.5)? 1: -1; 335 if (SideOfDecay < 0) 336 { 337 string->SetLeftPartonStable(); 338 } else 339 { 340 string->SetRightPartonStable(); 341 } 342 343 G4ParticleDefinition *newStringEnd; 344 G4ParticleDefinition * HadronDefinition; 345 if (string->DecayIsQuark()) 346 { 347 HadronDefinition= QuarkSplitup(string->GetDecayParton(), newStringEnd); 348 } else { 349 HadronDefinition= DiQuarkSplitup(string->GetDecayParton(), newStringEnd); 350 } 351 // create new String from old, ie. keep Left and Right order, but replace decay 352 353 newString=new G4FragmentingString(*string,newStringEnd); // To store possible 354 // quark containt of new string 355 G4LorentzVector* HadronMomentum=SplitEandP(HadronDefinition, string, newString); 356 357 delete newString; newString=0; // Uzhi 20.06.08 358 359 G4KineticTrack * Hadron =0; 360 if ( HadronMomentum != 0 ) { 361 362 G4ThreeVector Pos; 363 Hadron = new G4KineticTrack(HadronDefinition, 0,Pos, *HadronMomentum); 364 365 newString=new G4FragmentingString(*string,newStringEnd, 366 HadronMomentum); 367 368 //G4cout<<"Out G4VLong String Dec######################"<<G4endl; 369 //G4cout<<"newString 4Mom"<<newString->Get4Momentum()<<G4endl; 370 //G4cout<<"newString Pl "<<newString->LightConePlus()<<G4endl; 371 //G4cout<<"newString Mi "<<newString->LightConeMinus()<<G4endl; 372 //G4cout<<"newString Pts "<<newString->StablePt()<<G4endl; 373 //G4cout<<"newString Ptd "<<newString->DecayPt()<<G4endl; 374 //G4cout<<"newString M2 "<<newString->Mass2()<<G4endl; 375 //G4cout<<"newString M2 "<<newString->Mass()<<G4endl; 376 delete HadronMomentum; 377 } 378 return Hadron; 379 } 380 381 //-------------------------------------------------------------------------------------- 382 383 G4ParticleDefinition * 384 G4VLongitudinalStringDecay::QuarkSplitup(G4ParticleDefinition* 385 decay, G4ParticleDefinition *&created) 386 { 387 G4int IsParticle=(decay->GetPDGEncoding()>0) ? -1 : +1; // if we have a quark, 388 // we need antiquark 389 // (or diquark) 390 pDefPair QuarkPair = CreatePartonPair(IsParticle); 391 created = QuarkPair.second; 392 return hadronizer->Build(QuarkPair.first, decay); 393 394 } 395 396 //----------------------------------------------------------------------------- 397 398 G4ParticleDefinition *G4VLongitudinalStringDecay::DiQuarkSplitup( 399 G4ParticleDefinition* decay, 400 G4ParticleDefinition *&created) 401 { 402 //... can Diquark break or not? 403 if (G4UniformRand() < DiquarkBreakProb ){ 404 //... Diquark break 405 406 G4int stableQuarkEncoding = decay->GetPDGEncoding()/1000; 407 G4int decayQuarkEncoding = (decay->GetPDGEncoding()/100)%10; 408 if (G4UniformRand() < 0.5) 409 { 410 G4int Swap = stableQuarkEncoding; 411 stableQuarkEncoding = decayQuarkEncoding; 412 decayQuarkEncoding = Swap; 413 } 414 415 G4int IsParticle=(decayQuarkEncoding>0) ? -1 : +1; 416 // if we have a quark, we need antiquark) 417 pDefPair QuarkPair = CreatePartonPair(IsParticle,false); // no diquarks wanted 418 //... Build new Diquark 419 G4int QuarkEncoding=QuarkPair.second->GetPDGEncoding(); 420 G4int i10 = std::max(std::abs(QuarkEncoding), std::abs(stableQuarkEncoding)); 421 G4int i20 = std::min(std::abs(QuarkEncoding), std::abs(stableQuarkEncoding)); 422 G4int spin = (i10 != i20 && G4UniformRand() <= 0.5)? 1 : 3; 423 G4int NewDecayEncoding = -1*IsParticle*(i10 * 1000 + i20 * 100 + spin); 424 created = FindParticle(NewDecayEncoding); 425 G4ParticleDefinition * decayQuark=FindParticle(decayQuarkEncoding); 426 G4ParticleDefinition * had=hadronizer->Build(QuarkPair.first, decayQuark); 427 return had; 428 // return hadronizer->Build(QuarkPair.first, decayQuark); 429 430 } else { 431 //... Diquark does not break 432 433 G4int IsParticle=(decay->GetPDGEncoding()>0) ? +1 : -1; 434 // if we have a diquark, we need quark) 435 pDefPair QuarkPair = CreatePartonPair(IsParticle,false); // no diquarks wanted 436 created = QuarkPair.second; 437 438 G4ParticleDefinition * had=hadronizer->Build(QuarkPair.first, decay); 439 return had; 440 // return G4ParticleDefinition * had=hadronizer->Build(QuarkPair.first, decay); 441 } 442 } 443 444 //----------------------------------------------------------------------------- 141 445 142 446 G4int G4VLongitudinalStringDecay::SampleQuarkFlavor(void) … … 145 449 } 146 450 147 //----------------------------------------------------------------------------- -----------------------------451 //----------------------------------------------------------------------------- 148 452 149 453 G4VLongitudinalStringDecay::pDefPair G4VLongitudinalStringDecay::CreatePartonPair(G4int NeedParticle,G4bool AllowDiquarks) … … 170 474 } 171 475 172 //----------------------------------------------------------------------------- -----------------------------476 //----------------------------------------------------------------------------- 173 477 174 478 // G4ThreeVector G4VLongitudinalStringDecay::SampleQuarkPt() … … 180 484 // return G4ThreeVector(Pt * std::cos(phi),Pt * std::sin(phi),0); 181 485 // } 486 182 487 G4ThreeVector G4VLongitudinalStringDecay::SampleQuarkPt() 183 488 { … … 188 493 } 189 494 190 // ----------------------------------------------------------------------------------------------------------495 //****************************************************************************** 191 496 192 497 void G4VLongitudinalStringDecay::CalculateHadronTimePosition(G4double theInitialStringMass, G4KineticTrackVector* Hadrons) 193 498 { 499 194 500 // `yo-yo` formation time 195 const G4double kappa = 1.0 * GeV/fermi; 501 // const G4double kappa = 1.0 * GeV/fermi/4.; // Uzhi String tension 1.06.08 502 G4double kappa = GetStringTensionParameter(); 503 //G4cout<<"Kappa "<<kappa<<G4endl; // Uzhi 20.06.08 504 //G4int Uzhi; G4cin>>Uzhi; // Uzhi 20.06.08 196 505 for(size_t c1 = 0; c1 < Hadrons->size(); c1++) 197 506 { … … 211 520 } 212 521 213 //----------------------------------------------------------------------------- -----------------------------522 //----------------------------------------------------------------------------- 214 523 215 524 /* … … 237 546 */ 238 547 548 549 //***************************************************************************** 550 551 void G4VLongitudinalStringDecay::SetSigmaTransverseMomentum(G4double aValue) 552 { 553 if ( PastInitPhase ) { 554 throw G4HadronicException(__FILE__, __LINE__, "4VLongitudinalStringDecay::SetSigmaTransverseMomentum after FragmentString() not allowed"); 555 } else { 556 SigmaQT = aValue; 557 } 558 } 559 239 560 //---------------------------------------------------------------------------------------------------------- 240 561 241 G4ParticleDefinition * 242 G4VLongitudinalStringDecay::QuarkSplitup(G4ParticleDefinition* 243 decay, G4ParticleDefinition *&created) 244 { 245 G4int IsParticle=(decay->GetPDGEncoding()>0) ? -1 : +1; // if we have a quark, we need antiquark (or diquark) 246 pDefPair QuarkPair = CreatePartonPair(IsParticle); 247 created = QuarkPair.second; 248 return hadronizer->Build(QuarkPair.first, decay); 249 562 void G4VLongitudinalStringDecay::SetStrangenessSuppression(G4double aValue) 563 { 564 if ( PastInitPhase ) { 565 throw G4HadronicException(__FILE__, __LINE__, "4VLongitudinalStringDecay::SetStrangenessSuppression after FragmentString() not allowed"); 566 } else { 567 StrangeSuppress = aValue; 568 } 250 569 } 251 570 252 571 //---------------------------------------------------------------------------------------------------------- 253 572 254 G4ParticleDefinition *G4VLongitudinalStringDecay::DiQuarkSplitup(255 G4ParticleDefinition* decay,256 G4ParticleDefinition *&created)257 {258 //... can Diquark break or not?259 if (G4UniformRand() < DiquarkBreakProb ){260 //... Diquark break261 262 G4int stableQuarkEncoding = decay->GetPDGEncoding()/1000;263 G4int decayQuarkEncoding = (decay->GetPDGEncoding()/100)%10;264 if (G4UniformRand() < 0.5)265 {266 G4int Swap = stableQuarkEncoding;267 stableQuarkEncoding = decayQuarkEncoding;268 decayQuarkEncoding = Swap;269 }270 271 G4int IsParticle=(decayQuarkEncoding>0) ? -1 : +1;272 // if we have a quark, we need antiquark)273 pDefPair QuarkPair = CreatePartonPair(IsParticle,false); // no diquarks wanted274 //... Build new Diquark275 G4int QuarkEncoding=QuarkPair.second->GetPDGEncoding();276 G4int i10 = std::max(std::abs(QuarkEncoding), std::abs(stableQuarkEncoding));277 G4int i20 = std::min(std::abs(QuarkEncoding), std::abs(stableQuarkEncoding));278 G4int spin = (i10 != i20 && G4UniformRand() <= 0.5)? 1 : 3;279 G4int NewDecayEncoding = -1*IsParticle*(i10 * 1000 + i20 * 100 + spin);280 created = FindParticle(NewDecayEncoding);281 G4ParticleDefinition * decayQuark=FindParticle(decayQuarkEncoding);282 G4ParticleDefinition * had=hadronizer->Build(QuarkPair.first, decayQuark);283 return had;284 // return hadronizer->Build(QuarkPair.first, decayQuark);285 286 } else {287 //... Diquark does not break288 289 G4int IsParticle=(decay->GetPDGEncoding()>0) ? +1 : -1;290 // if we have a diquark, we need quark)291 pDefPair QuarkPair = CreatePartonPair(IsParticle,false); // no diquarks wanted292 created = QuarkPair.second;293 294 G4ParticleDefinition * had=hadronizer->Build(QuarkPair.first, decay);295 return had;296 // return G4ParticleDefinition * had=hadronizer->Build(QuarkPair.first, decay);297 }298 }299 300 //-----------------------------------------------------------------------------------------301 302 G4KineticTrack * G4VLongitudinalStringDecay::Splitup(303 G4FragmentingString *string,304 G4FragmentingString *&newString)305 {306 307 //... random choice of string end to use for creating the hadron (decay)308 SideOfDecay = (G4UniformRand() < 0.5)? 1: -1;309 if (SideOfDecay < 0)310 {311 string->SetLeftPartonStable();312 } else313 {314 string->SetRightPartonStable();315 }316 317 G4ParticleDefinition *newStringEnd;318 G4ParticleDefinition * HadronDefinition;319 if (string->DecayIsQuark())320 {321 HadronDefinition= QuarkSplitup(string->GetDecayParton(), newStringEnd);322 } else {323 HadronDefinition= DiQuarkSplitup(string->GetDecayParton(), newStringEnd);324 }325 // create new String from old, ie. keep Left and Right order, but replace decay326 G4LorentzVector* HadronMomentum=SplitEandP(HadronDefinition, string);327 328 G4KineticTrack * Hadron =0;329 if ( HadronMomentum != 0 ) {330 331 G4ThreeVector Pos;332 Hadron = new G4KineticTrack(HadronDefinition, 0,Pos, *HadronMomentum);333 334 newString=new G4FragmentingString(*string,newStringEnd,335 HadronMomentum);336 337 delete HadronMomentum;338 }339 return Hadron;340 }341 342 //----------------------------------------------------------------------------------------------------------343 344 G4ExcitedString *G4VLongitudinalStringDecay::CPExcited(const G4ExcitedString & in)345 {346 G4Parton *Left=new G4Parton(*in.GetLeftParton());347 G4Parton *Right=new G4Parton(*in.GetRightParton());348 return new G4ExcitedString(Left,Right,in.GetDirection());349 }350 351 G4double G4VLongitudinalStringDecay::FragmentationMass(352 const G4FragmentingString *353 const string,354 Pcreate build,355 pDefPair * pdefs)356 {357 358 G4double mass;359 static G4bool NeedInit(true);360 static std::vector<double> nomix;361 static G4HadronBuilder * minMassHadronizer;362 if ( NeedInit )363 {364 NeedInit = false;365 nomix.resize(6);366 for ( G4int i=0; i<6 ; i++ ) nomix[i]=0;367 // minMassHadronizer=new G4HadronBuilder(pspin_meson,pspin_barion,nomix,nomix);368 minMassHadronizer=hadronizer;369 }370 371 if ( build==0 ) build=&G4HadronBuilder::BuildLowSpin;372 373 G4ParticleDefinition *Hadron1, *Hadron2=0;374 375 if (!string->FourQuarkString() )376 {377 // spin 0 meson or spin 1/2 barion will be built378 379 Hadron1 = (minMassHadronizer->*build)(string->GetLeftParton(),380 string->GetRightParton());381 mass= (Hadron1)->GetPDGMass();382 } else383 {384 //... string is qq--qqbar: Build two stable hadrons,385 //... with extra uubar or ddbar quark pair386 G4int iflc = (G4UniformRand() < 0.5)? 1 : 2;387 if (string->GetLeftParton()->GetPDGEncoding() < 0) iflc = -iflc;388 389 //... theSpin = 4; spin 3/2 baryons will be built390 Hadron1 = (minMassHadronizer->*build)(string->GetLeftParton(),FindParticle(iflc));391 Hadron2 =(minMassHadronizer->*build)(string->GetRightParton(),FindParticle(-iflc));392 mass = (Hadron1)->GetPDGMass() + (Hadron2)->GetPDGMass();393 }394 395 if ( pdefs != 0 )396 { // need to return hadrons as well....397 pdefs->first = Hadron1;398 pdefs->second = Hadron2;399 }400 401 return mass;402 }403 404 G4KineticTrackVector* G4VLongitudinalStringDecay::LightFragmentationTest(const405 G4ExcitedString * const string)406 {407 // Check string decay threshold408 409 G4KineticTrackVector * result=0; // return 0 when string exceeds the mass cut410 411 pDefPair hadrons((G4ParticleDefinition *)0,(G4ParticleDefinition *)0);412 G4FragmentingString aString(*string);413 if ( sqr(FragmentationMass(&aString,0,&hadrons)+MassCut) < aString.Mass2()) {414 return 0;415 }416 417 result=new G4KineticTrackVector;418 419 if ( hadrons.second ==0 )420 {421 // Substitute string by light hadron, Note that Energy is not conserved here!422 423 #ifdef DEBUG_LightFragmentationTest424 G4cout << "VlongSF Warning replacing string by single hadron "425 << hadrons.first->GetParticleName()426 << "string .. " << string->Get4Momentum() << " "427 << string->Get4Momentum().m() << G4endl;428 #endif429 430 G4ThreeVector Mom3 = string->Get4Momentum().vect();431 G4LorentzVector Mom(Mom3,432 std::sqrt(Mom3.mag2() + sqr(hadrons.first->GetPDGMass())));433 result->push_back(new G4KineticTrack(hadrons.first, 0, string->GetPosition(), Mom));434 } else435 {436 //... string was qq--qqbar type: Build two stable hadrons,437 438 #ifdef DEBUG_LightFragmentationTest439 G4cout << "VlongSF Warning replacing qq-qqbar string by TWO hadrons "440 << hadrons.first->GetParticleName() << " / "441 << hadrons.second->GetParticleName()442 << "string .. " << string->Get4Momentum() << " "443 << string->Get4Momentum().m() << G4endl;444 #endif445 446 G4LorentzVector Mom1, Mom2;447 Sample4Momentum(&Mom1, hadrons.first->GetPDGMass(),448 &Mom2,hadrons.second->GetPDGMass(),449 string->Get4Momentum().mag());450 result->push_back(new G4KineticTrack(hadrons.first, 0, string->GetPosition(), Mom1));451 result->push_back(new G4KineticTrack(hadrons.second, 0, string->GetPosition(), Mom2));452 G4ThreeVector Velocity = string->Get4Momentum().boostVector();453 result->Boost(Velocity);454 }455 456 return result;457 458 }459 460 //----------------------------------------------------------------------------------------------------------461 462 G4ParticleDefinition* G4VLongitudinalStringDecay::FindParticle(G4int Encoding)463 {464 G4ParticleDefinition* ptr = G4ParticleTable::GetParticleTable()->FindParticle(Encoding);465 if (ptr == NULL)466 {467 G4cout << "Particle with encoding "<<Encoding<<" does not exist!!!"<<G4endl;468 throw G4HadronicException(__FILE__, __LINE__, "Check your particle table");469 }470 return ptr;471 }472 473 //----------------------------------------------------------------------------------------------------------474 475 void G4VLongitudinalStringDecay::SetSigmaTransverseMomentum(G4double aValue)476 {477 if ( PastInitPhase ) {478 throw G4HadronicException(__FILE__, __LINE__, "4VLongitudinalStringDecay::SetSigmaTransverseMomentum after FragmentString() not allowed");479 } else {480 SigmaQT = aValue;481 }482 }483 484 //----------------------------------------------------------------------------------------------------------485 486 void G4VLongitudinalStringDecay::SetStrangenessSuppression(G4double aValue)487 {488 if ( PastInitPhase ) {489 throw G4HadronicException(__FILE__, __LINE__, "4VLongitudinalStringDecay::SetStrangenessSuppression after FragmentString() not allowed");490 } else {491 StrangeSuppress = aValue;492 }493 }494 495 //----------------------------------------------------------------------------------------------------------496 497 573 void G4VLongitudinalStringDecay::SetDiquarkSuppression(G4double aValue) 498 574 { … … 503 579 } 504 580 } 581 582 //---------------------------------------------------------------------------------------- 505 583 506 584 void G4VLongitudinalStringDecay::SetDiquarkBreakProbability(G4double aValue) … … 582 660 583 661 } 662 } 663 664 //------------------------------------------------------------------------------------------- 665 void G4VLongitudinalStringDecay::SetStringTensionParameter(G4double aValue)// Uzhi 20 June 08 666 { 667 Kappa = aValue * GeV/fermi; 584 668 } 585 //******************************************************************************************************* 669 //************************************************************************************** 670
Note: See TracChangeset
for help on using the changeset viewer.
