- Timestamp:
- Nov 5, 2010, 3:45:55 PM (14 years ago)
- Location:
- trunk/source/processes/hadronic/models/parton_string/hadronization/src
- Files:
-
- 6 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/source/processes/hadronic/models/parton_string/hadronization/src/G4ExcitedStringDecay.cc
r1337 r1340 61 61 } 62 62 63 G4bool G4ExcitedStringDecay:: 64 EnergyAndMomentumCorrector(G4KineticTrackVector* Output, G4LorentzVector& TotalCollisionMom) 63 G4KineticTrackVector *G4ExcitedStringDecay::FragmentString 64 (const G4ExcitedString &theString) 65 { 66 if ( theStringDecay == NULL ) 67 68 theStringDecay=new G4LundStringFragmentation(); 69 70 return theStringDecay->FragmentString(theString); 71 } 72 73 G4KineticTrackVector *G4ExcitedStringDecay::FragmentStrings 74 (const G4ExcitedStringVector * theStrings) 75 { 76 G4LorentzVector KTsum(0.,0.,0.,0.); 77 for ( unsigned int astring=0; astring < theStrings->size(); astring++) 78 { 79 KTsum+= theStrings->operator[](astring)->Get4Momentum(); 80 81 if( !(KTsum.e()<1) && !(KTsum.e()>-1) ) 82 { 83 throw G4HadronicException(__FILE__, __LINE__, 84 "G4ExcitedStringDecay::FragmentStrings received nan string..."); 85 } 86 } 87 88 G4KineticTrackVector * theResult = new G4KineticTrackVector; 89 G4int attempts(0); 90 G4bool success=false; 91 do { 92 93 std::for_each(theResult->begin() , theResult->end() , DeleteKineticTrack()); 94 theResult->clear(); 95 96 attempts++; 97 G4LorentzVector KTsecondaries(0.,0.,0.,0.); 98 G4bool NeedEnergyCorrector=false; 99 100 for ( unsigned int astring=0; astring < theStrings->size(); astring++) 101 { 102 //G4cout<<G4endl<<"String No "<<astring+1<<" "<<theStrings->operator[](astring)->Get4Momentum().mag()<<G4endl; 103 104 G4KineticTrackVector * generatedKineticTracks = NULL; 105 106 if ( theStrings->operator[](astring)->IsExcited() ) 107 { 108 generatedKineticTracks=FragmentString(*theStrings->operator[](astring)); 109 } else { 110 generatedKineticTracks = new G4KineticTrackVector; 111 generatedKineticTracks->push_back(theStrings->operator[](astring)->GetKineticTrack()); 112 } 113 114 if (generatedKineticTracks == NULL) 115 { 116 G4cerr << "G4VPartonStringModel:No KineticTracks produced" << G4endl; 117 continue; 118 } 119 120 G4LorentzVector KTsum1(0.,0.,0.,0.); 121 for ( unsigned int aTrack=0; aTrack<generatedKineticTracks->size();aTrack++) 122 { 123 //G4cout<<" Prod part "<<(*generatedKineticTracks)[aTrack]->GetDefinition()->GetParticleName()<<" "<<(*generatedKineticTracks)[aTrack]->Get4Momentum()<<G4endl; 124 theResult->push_back(generatedKineticTracks->operator[](aTrack)); 125 KTsum1+= (*generatedKineticTracks)[aTrack]->Get4Momentum(); 126 } 127 KTsecondaries+=KTsum1; 128 129 if ( KTsum1.e() > 0 && std::abs((KTsum1.e()-theStrings->operator[](astring)->Get4Momentum().e()) / KTsum1.e()) > perMillion ) 130 { 131 //--debug-- G4cout << "String secondaries(" <<generatedKineticTracks->size()<< ") momentum: " 132 //--debug-- << theStrings->operator[](astring)->Get4Momentum() << " " << KTsum1 << G4endl; 133 NeedEnergyCorrector=true; 134 } 135 136 // clean up 137 delete generatedKineticTracks; 138 } 139 //--DEBUG G4cout << "Strings/secs total 4 momentum " << KTsum << " " <<KTsecondaries << G4endl; 140 141 success=true; 142 if ( NeedEnergyCorrector ) success=EnergyAndMomentumCorrector(theResult, KTsum); 143 144 } while(!success && (attempts < 100)); 145 146 #ifdef debug_ExcitedStringDecay 147 G4LorentzVector KTsum1=0; 148 for ( unsigned int aTrack=0; aTrack<theResult->size();aTrack++) 149 { 150 G4cout << " corrected tracks .. " << (*theResult)[aTrack]->GetDefinition()->GetParticleName() 151 <<" " << (*theResult)[aTrack]->Get4Momentum() << G4endl; 152 KTsum1+= (*theResult)[aTrack]->Get4Momentum(); 153 } 154 G4cout << "Needcorrector/success " << NeedEnergyCorrector << "/" << success << ", Corrected total 4 momentum " << KTsum1 << G4endl; 155 if ( ! success ) G4cout << "failed to correct E/p" << G4endl; 156 #endif 157 158 return theResult; 159 } 160 161 G4bool G4ExcitedStringDecay::EnergyAndMomentumCorrector 162 (G4KineticTrackVector* Output, G4LorentzVector& TotalCollisionMom) 65 163 { 66 164 const int nAttemptScale = 500; … … 71 169 G4double SumMass = 0; 72 170 G4double TotalCollisionMass = TotalCollisionMom.m(); 171 73 172 if( !(TotalCollisionMass<1) && !(TotalCollisionMass>-1) ) 74 173 { … … 77 176 } 78 177 178 //G4cout<<G4endl<<"EnergyAndMomentumCorrector "<<Output->size()<<G4endl; 79 179 // Calculate sum hadron 4-momenta and summing hadron mass 80 180 unsigned int cHadron; … … 92 192 } 93 193 } 194 195 //G4cout<<"SumMass TotalCollisionMass "<<SumMass<<" "<<TotalCollisionMass<<G4endl; 94 196 95 197 // Cannot correct a single particle -
trunk/source/processes/hadronic/models/parton_string/hadronization/src/G4HadronBuilder.cc
r1337 r1340 25 25 // 26 26 // 27 // $Id: G4HadronBuilder.cc,v 1.1 0 2009/05/22 16:34:31 gunterExp $28 // GEANT4 tag $Name: geant4-09-0 4-beta-01$27 // $Id: G4HadronBuilder.cc,v 1.11 2010/09/20 12:46:23 vuzhinsk Exp $ 28 // GEANT4 tag $Name: geant4-09-03-ref-09 $ 29 29 // 30 30 // ----------------------------------------------------------------------------- -
trunk/source/processes/hadronic/models/parton_string/hadronization/src/G4LundStringFragmentation.cc
r1337 r1340 25 25 // 26 26 // 27 // $Id: G4LundStringFragmentation.cc,v 1.2 0 2010/06/21 17:50:48vuzhinsk Exp $28 // GEANT4 tag $Name: geant4-09-0 4-beta-01$ 1.827 // $Id: G4LundStringFragmentation.cc,v 1.23 2010/09/22 12:36:37 vuzhinsk Exp $ 28 // GEANT4 tag $Name: geant4-09-03-ref-09 $ 1.8 29 29 // 30 30 // ----------------------------------------------------------------------------- … … 53 53 MinimalStringMass2 = 0.; 54 54 // ------ Minimal invariant mass used at a string fragmentation - 55 WminLUND = 0. 7*GeV; // Uzhi0.8 1.555 WminLUND = 0.45*GeV; //0.23*GeV; // Uzhi 0.7 -> 0.23 3.8.10 //0.8 1.5 56 56 // ------ Smooth parameter used at a string fragmentation for --- 57 57 // ------ smearinr sharp mass cut-off --------------------------- … … 60 60 // SetStringTensionParameter(0.25); 61 61 SetStringTensionParameter(1.); 62 63 SetDiquarkBreakProbability(0.05); // Vova Aug. 22 64 65 // For treating of small string decays 66 for(G4int i=0; i<3; i++) 67 { for(G4int j=0; j<3; j++) 68 { for(G4int k=0; k<6; k++) 69 { Meson[i][j][k]=0; MesonWeight[i][j][k]=0.; 70 } 71 } 72 } 73 //-------------------------- 74 Meson[0][0][0]=111; // dbar-d Pi0 75 MesonWeight[0][0][0]=pspin_meson*(1.-scalarMesonMix[0]); 76 77 Meson[0][0][1]=221; // dbar-d Eta 78 MesonWeight[0][0][1]=pspin_meson*(scalarMesonMix[0]-scalarMesonMix[1]); 79 80 Meson[0][0][2]=331; // dbar-d EtaPrime 81 MesonWeight[0][0][2]=pspin_meson*(scalarMesonMix[1]); 82 83 Meson[0][0][3]=113; // dbar-d Rho0 84 MesonWeight[0][0][3]=(1.-pspin_meson)*(1.-vectorMesonMix[0]); 85 86 Meson[0][0][4]=223; // dbar-d Omega 87 MesonWeight[0][0][4]=(1.-pspin_meson)*(vectorMesonMix[0]); 88 //-------------------------- 89 90 Meson[0][1][0]=211; // dbar-u Pi+ 91 MesonWeight[0][1][0]=pspin_meson; 92 93 Meson[0][1][1]=213; // dbar-u Rho+ 94 MesonWeight[0][1][1]=(1.-pspin_meson); 95 //-------------------------- 96 97 Meson[0][2][0]=311; // dbar-s K0bar 98 MesonWeight[0][2][0]=pspin_meson; 99 100 Meson[0][2][1]=313; // dbar-s K*0bar 101 MesonWeight[0][2][1]=(1.-pspin_meson); 102 //-------------------------- 103 //-------------------------- 104 Meson[1][0][0]=211; // ubar-d Pi- 105 MesonWeight[1][0][0]=pspin_meson; 106 107 Meson[1][0][1]=213; // ubar-d Rho- 108 MesonWeight[1][0][1]=(1.-pspin_meson); 109 //-------------------------- 110 111 Meson[1][1][0]=111; // ubar-u Pi0 112 MesonWeight[1][1][0]=pspin_meson*(1.-scalarMesonMix[0]); 113 114 Meson[1][1][1]=221; // ubar-u Eta 115 MesonWeight[1][1][1]=pspin_meson*(scalarMesonMix[0]-scalarMesonMix[1]); 116 117 Meson[1][1][2]=331; // ubar-u EtaPrime 118 MesonWeight[1][1][2]=pspin_meson*(scalarMesonMix[1]); 119 120 Meson[1][1][3]=113; // ubar-u Rho0 121 MesonWeight[1][1][3]=(1.-pspin_meson)*(1.-vectorMesonMix[0]); 122 123 Meson[1][1][4]=223; // ubar-u Omega 124 MesonWeight[1][1][4]=(1.-pspin_meson)*(scalarMesonMix[0]); 125 //-------------------------- 126 127 Meson[1][2][0]=321; // ubar-s K- 128 MesonWeight[1][2][0]=pspin_meson; 129 130 Meson[1][2][1]=323; // ubar-s K*-bar - 131 MesonWeight[1][2][1]=(1.-pspin_meson); 132 //-------------------------- 133 //-------------------------- 134 135 Meson[2][0][0]=311; // sbar-d K0 136 MesonWeight[2][0][0]=pspin_meson; 137 138 Meson[2][0][1]=313; // sbar-d K*0 139 MesonWeight[2][0][1]=(1.-pspin_meson); 140 //-------------------------- 141 142 Meson[2][1][0]=321; // sbar-u K+ 143 MesonWeight[2][1][0]=pspin_meson; 144 145 Meson[2][1][1]=323; // sbar-u K*+ 146 MesonWeight[2][1][1]=(1.-pspin_meson); 147 //-------------------------- 148 149 Meson[2][2][0]=221; // sbar-s Eta 150 MesonWeight[2][2][0]=pspin_meson*(1.-scalarMesonMix[5]); 151 152 Meson[2][2][1]=331; // sbar-s EtaPrime 153 MesonWeight[2][2][1]=pspin_meson*(1.-scalarMesonMix[5]); 154 155 Meson[2][2][3]=333; // sbar-s EtaPrime 156 MesonWeight[2][2][3]=(1.-pspin_meson)*(vectorMesonMix[5]); 157 //-------------------------- 158 159 for(G4int i=0; i<3; i++) 160 { for(G4int j=0; j<3; j++) 161 { for(G4int k=0; k<3; k++) 162 { for(G4int l=0; l<4; l++) 163 { Baryon[i][j][k][l]=0; BaryonWeight[i][j][k][l]=0.;} 164 } 165 } 166 } 167 168 //--------------------------------------- 169 Baryon[0][0][0][0]=1114; // Delta- 170 BaryonWeight[0][0][0][0]=1.; 171 172 //--------------------------------------- 173 Baryon[0][0][1][0]=2112; // neutron 174 BaryonWeight[0][0][1][0]=pspin_barion; 175 176 Baryon[0][0][1][1]=2114; // Delta0 177 BaryonWeight[0][0][1][1]=(1.-pspin_barion); 178 179 //--------------------------------------- 180 Baryon[0][0][2][0]=3112; // Sigma- 181 BaryonWeight[0][0][2][0]=pspin_barion; 182 183 Baryon[0][0][2][1]=3114; // Sigma*- 184 BaryonWeight[0][0][2][1]=(1.-pspin_barion); 185 186 //--------------------------------------- 187 Baryon[0][1][0][0]=2112; // neutron 188 BaryonWeight[0][1][0][0]=pspin_barion; 189 190 Baryon[0][1][0][1]=2114; // Delta0 191 BaryonWeight[0][1][0][1]=(1.-pspin_barion); 192 193 //--------------------------------------- 194 Baryon[0][1][1][0]=2212; // proton 195 BaryonWeight[0][1][1][0]=pspin_barion; 196 197 Baryon[0][1][1][1]=2214; // Delta+ 198 BaryonWeight[0][1][1][1]=(1.-pspin_barion); 199 200 //--------------------------------------- 201 Baryon[0][1][2][0]=3122; // Lambda 202 BaryonWeight[0][1][2][0]=pspin_barion*0.5; 203 204 Baryon[0][1][2][1]=3212; // Sigma0 205 BaryonWeight[0][1][2][2]=pspin_barion*0.5; 206 207 Baryon[0][1][2][2]=3214; // Sigma*0 208 BaryonWeight[0][1][2][2]=(1.-pspin_barion); 209 210 //--------------------------------------- 211 Baryon[0][2][0][0]=3112; // Sigma- 212 BaryonWeight[0][2][0][0]=pspin_barion; 213 214 Baryon[0][2][0][1]=3114; // Sigma*- 215 BaryonWeight[0][2][0][1]=(1.-pspin_barion); 216 217 //--------------------------------------- 218 Baryon[0][2][1][0]=3122; // Lambda 219 BaryonWeight[0][2][1][0]=pspin_barion*0.5; 220 221 Baryon[0][2][1][1]=3212; // Sigma0 222 BaryonWeight[0][2][1][1]=pspin_barion*0.5; 223 224 Baryon[0][2][1][2]=3214; // Sigma*0 225 BaryonWeight[0][2][1][2]=(1.-pspin_barion); 226 227 //--------------------------------------- 228 Baryon[0][2][2][0]=3312; // Theta- 229 BaryonWeight[0][2][2][0]=pspin_barion; 230 231 Baryon[0][2][2][1]=3314; // Theta*- 232 BaryonWeight[0][2][2][1]=(1.-pspin_barion); 233 234 //--------------------------------------- 235 //--------------------------------------- 236 Baryon[1][0][0][0]=2112; // neutron 237 BaryonWeight[1][0][0][0]=pspin_barion; 238 239 Baryon[1][0][0][1]=2114; // Delta0 240 BaryonWeight[1][0][0][1]=(1.-pspin_barion); 241 242 //--------------------------------------- 243 Baryon[1][0][1][0]=2212; // proton 244 BaryonWeight[1][0][1][0]=pspin_barion; 245 246 Baryon[1][0][1][1]=2214; // Delta+ 247 BaryonWeight[1][0][1][1]=(1.-pspin_barion); 248 249 //--------------------------------------- 250 Baryon[1][0][2][0]=3122; // Lambda 251 BaryonWeight[1][0][2][0]=pspin_barion*0.5; 252 253 Baryon[1][0][2][1]=3212; // Sigma0 254 BaryonWeight[1][0][2][1]=pspin_barion*0.5; 255 256 Baryon[1][0][2][2]=3214; // Sigma*0 257 BaryonWeight[1][0][2][2]=(1.-pspin_barion); 258 259 //--------------------------------------- 260 Baryon[1][1][0][0]=2212; // proton 261 BaryonWeight[1][1][0][0]=pspin_barion; 262 263 Baryon[1][1][0][1]=2214; // Delta+ 264 BaryonWeight[1][1][0][1]=(1.-pspin_barion); 265 266 //--------------------------------------- 267 Baryon[1][1][1][0]=2224; // Delta++ 268 BaryonWeight[1][1][1][0]=1.; 269 270 //--------------------------------------- 271 Baryon[1][1][2][0]=3222; // Sigma+ 272 BaryonWeight[1][1][2][0]=pspin_barion; 273 274 Baryon[1][1][2][1]=3224; // Sigma*+ 275 BaryonWeight[1][1][2][1]=(1.-pspin_barion); 276 277 //--------------------------------------- 278 Baryon[1][2][0][0]=3122; // Lambda 279 BaryonWeight[1][2][0][0]=pspin_barion*0.5; 280 281 Baryon[1][2][0][1]=3212; // Sigma0 282 BaryonWeight[1][2][0][1]=pspin_barion*0.5; 283 284 Baryon[1][2][0][2]=3214; // Sigma*0 285 BaryonWeight[1][2][0][2]=(1.-pspin_barion); 286 287 //--------------------------------------- 288 Baryon[1][2][1][0]=3222; // Sigma+ 289 BaryonWeight[1][2][1][0]=pspin_barion; 290 291 Baryon[1][2][1][1]=3224; // Sigma*+ 292 BaryonWeight[1][2][1][1]=(1.-pspin_barion); 293 294 //--------------------------------------- 295 Baryon[1][2][2][0]=3322; // Theta0 296 BaryonWeight[1][2][2][0]=pspin_barion; 297 298 Baryon[1][2][2][1]=3324; // Theta*0 299 BaryonWeight[1][2][2][1]=(1.-pspin_barion); 300 301 //--------------------------------------- 302 //--------------------------------------- 303 Baryon[2][0][0][0]=3112; // Sigma- 304 BaryonWeight[2][0][0][0]=pspin_barion; 305 306 Baryon[2][0][0][1]=3114; // Sigma*- 307 BaryonWeight[2][0][0][1]=(1.-pspin_barion); 308 309 //--------------------------------------- 310 Baryon[2][0][1][0]=3122; // Lambda 311 BaryonWeight[2][0][1][0]=pspin_barion*0.5; 312 313 Baryon[2][0][1][1]=3212; // Sigma0 314 BaryonWeight[2][0][1][1]=pspin_barion*0.5; 315 316 Baryon[2][0][1][2]=3214; // Sigma*0 317 BaryonWeight[2][0][1][2]=(1.-pspin_barion); 318 319 //--------------------------------------- 320 Baryon[2][0][2][0]=3312; // Sigma- 321 BaryonWeight[2][0][2][0]=pspin_barion; 322 323 Baryon[2][0][2][1]=3314; // Sigma*- 324 BaryonWeight[2][0][2][1]=(1.-pspin_barion); 325 326 //--------------------------------------- 327 Baryon[2][1][0][0]=3122; // Lambda 328 BaryonWeight[2][1][0][0]=pspin_barion*0.5; 329 330 Baryon[2][1][0][1]=3212; // Sigma0 331 BaryonWeight[2][1][0][1]=pspin_barion*0.5; 332 333 Baryon[2][1][0][2]=3214; // Sigma*0 334 BaryonWeight[2][1][0][2]=(1.-pspin_barion); 335 336 //--------------------------------------- 337 Baryon[2][1][1][0]=3222; // Sigma+ 338 BaryonWeight[2][1][1][0]=pspin_barion; 339 340 Baryon[2][1][1][1]=3224; // Sigma*+ 341 BaryonWeight[2][1][1][1]=(1.-pspin_barion); 342 343 //--------------------------------------- 344 Baryon[2][1][2][0]=3322; // Theta0 345 BaryonWeight[2][1][2][0]=pspin_barion; 346 347 Baryon[2][1][2][1]=3324; // Theta*0 348 BaryonWeight[2][1][2][2]=(1.-pspin_barion); 349 350 //--------------------------------------- 351 Baryon[2][2][0][0]=3312; // Theta- 352 BaryonWeight[2][2][0][0]=pspin_barion; 353 354 Baryon[2][2][0][1]=3314; // Theta*- 355 BaryonWeight[2][2][0][1]=(1.-pspin_barion); 356 357 //--------------------------------------- 358 Baryon[2][2][1][0]=3322; // Theta0 359 BaryonWeight[2][2][1][0]=pspin_barion; 360 361 Baryon[2][2][1][1]=3324; // Theta*0 362 BaryonWeight[2][2][1][1]=(1.-pspin_barion); 363 364 //--------------------------------------- 365 Baryon[2][2][2][0]=3334; // Omega 366 BaryonWeight[2][2][2][0]=1.; 367 368 //--------------------------------------- 369 /* 370 for(G4int i=0; i<3; i++) 371 { for(G4int j=0; j<3; j++) 372 { for(G4int k=0; k<3; k++) 373 { for(G4int l=0; l<4; l++) 374 { G4cout<<i<<" "<<j<<" "<<k<<" "<<l<<" "<<Baryon[i][j][k][l]<<G4endl;} 375 } 376 } 377 } 378 G4int Uzhi; 379 G4cin>>Uzhi; 380 */ 381 //StrangeSuppress=0.38; 382 Prob_QQbar[0]=StrangeSuppress; // Probability of ddbar production 383 Prob_QQbar[1]=StrangeSuppress; // Probability of uubar production 384 Prob_QQbar[2]=StrangeSuppress/(2.+StrangeSuppress);//(1.-2.*StrangeSuppress); // Probability of ssbar production 62 385 } 63 386 … … 160 483 // Can no longer modify Parameters for Fragmentation. 161 484 PastInitPhase=true; 485 //SetVectorMesonProbability(1.); 486 //SetSpinThreeHalfBarionProbability(1.); 162 487 163 488 // check if string has enough mass to fragment... … … 165 490 SetMassCut(160.*MeV); // For LightFragmentationTest it is required 166 491 // that no one pi-meson can be produced 492 // 493 //G4cout<<G4endl<<"G4LundStringFragmentation::"<<G4endl; 494 //G4cout<<"FragmentString Position"<<theString.GetPosition()/fermi<<" "<<theString.GetTimeOfCreation()/fermi<<G4endl; 495 //G4cout<<"FragmentString Momentum"<<theString.Get4Momentum()<<theString.Get4Momentum().mag()<<G4endl; 496 // 497 G4FragmentingString aString(theString); 498 SetMinimalStringMass(&aString); 499 167 500 /* 168 G4cout<<G4endl<<"G4LundStringFragmentation::"<<G4endl; 169 G4cout<<"FragmentString Position"<<theString.GetPosition()/fermi<<" "<< 170 theString.GetTimeOfCreation()/fermi<<G4endl; 171 G4cout<<"FragmentString Momentum"<<theString.Get4Momentum()<<theString.Get4Momentum().mag()<<G4endl; 501 G4cout<<"St Frag "<<MinimalStringMass<<" "<<WminLUND<<" "<<(1.-SmoothParam)<<" "<< theString.Get4Momentum().mag()<<G4endl; 502 G4cout<<(MinimalStringMass+WminLUND)*(1.-SmoothParam)<<" "<<theString.Get4Momentum().mag()<<G4endl; 503 504 if((MinimalStringMass+WminLUND)*(1.-SmoothParam) > theString.Get4Momentum().mag()) 505 {SetMassCut(1000.*MeV);} 506 else {SetMassCut(160.*MeV);} 172 507 */ 173 G4FragmentingString aString(theString); 174 SetMinimalStringMass(&aString); 175 176 if((MinimalStringMass+WminLUND)*(1.-SmoothParam) > theString.Get4Momentum().mag()) 177 {SetMassCut(1000.*MeV);} 178 // V.U. 20.06.10 in order to put un correspondence LightFragTest and MinStrMass 179 180 G4KineticTrackVector * LeftVector=LightFragmentationTest(&theString); 508 // V.U. 20.06.10 in order to put in correspondence LightFragTest and MinStrMass 509 510 G4KineticTrackVector * LeftVector(0); 511 // G4KineticTrackVector * LeftVector=LightFragmentationTest(&theString); 512 //G4cout<<"Min Str Mass "<<MinimalStringMass<<G4endl; 513 if(!IsFragmentable(&aString)) // produce 1 hadron 514 { 515 //G4cout<<"Non fragmentable"<<G4endl; 516 SetMassCut(1000.*MeV); 517 LeftVector=LightFragmentationTest(&theString); 518 SetMassCut(160.*MeV); 519 520 521 //G4cout<<"Prod hadron "<<LeftVector->operator[](0)->GetDefinition()->GetParticleName()<<G4endl; 522 /* 523 if(LeftVector->size() == 1) 524 { 525 if(! (std::abs(LeftVector->operator[](0)->GetDefinition()->GetPDGMass()- 526 theString.Get4Momentum().mag()) < 100*MeV)) 527 { // produce a particle with arbitrary isospin 528 G4cout<<" produce a particle with arbitrary isospin"<<G4endl; 529 530 pDefPair hadrons((G4ParticleDefinition *)0,(G4ParticleDefinition *)0); 531 Pcreate build=&G4HadronBuilder::Build; 532 FragmentationMass(&aString,build,&hadrons); // 0->1 533 G4cout<<"Hadron PDG "<<hadrons.first->GetPDGEncoding()<<G4endl; 534 if ( hadrons.second ==0 ) 535 {// Substitute string by light hadron, Note that Energy is not conserved here! 536 // Cleaning up the previously produced hadrons ------------------------------ 537 std::for_each(LeftVector->begin() , LeftVector->end() , DeleteKineticTrack()); 538 LeftVector->clear(); 539 540 G4ThreeVector Mom3 = aString.Get4Momentum().vect(); 541 G4LorentzVector Mom(Mom3,std::sqrt(Mom3.mag2() + sqr(hadrons.first->GetPDGMass()))); 542 LeftVector->push_back(new G4KineticTrack(hadrons.first, 0, 543 theString.GetPosition(), 544 Mom)); 545 } // end of if ( hadrons.second ==0 ) 546 } // end of produce a particle with arbitrary isospin 547 548 } // end of if(LeftVector->size() == 1) 549 */ 550 } // end of if(!IsFragmentable(&aString)) 181 551 182 552 if ( LeftVector != 0 ) { … … 207 577 //--------------------- The string can fragment ------------------------------- 208 578 //--------------- At least two particles can be produced ---------------------- 579 //G4cout<<"Fragmentable"<<G4endl; 209 580 LeftVector =new G4KineticTrackVector; 210 581 G4KineticTrackVector * RightVector=new G4KineticTrackVector; … … 218 589 { // If the string fragmentation do not be happend, repeat the fragmentation--- 219 590 G4FragmentingString *currentString=new G4FragmentingString(*theStringInCMS); 220 591 //G4cout<<"Main loop start whilecounter "<<attempt<<G4endl; 221 592 // Cleaning up the previously produced hadrons ------------------------------ 222 593 std::for_each(LeftVector->begin() , LeftVector->end() , DeleteKineticTrack()); … … 228 599 // Main fragmentation loop until the string will not be able to fragment ---- 229 600 inner_sucess=true; // set false on failure.. 601 230 602 while (! StopFragmenting(currentString) ) 231 603 { // Split current string into hadron + new string 604 232 605 G4FragmentingString *newString=0; // used as output from SplitUp... 233 606 234 607 G4KineticTrack * Hadron=Splitup(currentString,newString); 608 //G4cout<<"SplitUp------"<<Hadron<<G4endl; 235 609 236 610 if ( Hadron != 0 ) // Store the hadron 237 611 { 612 //G4cout<<"Hadron prod at fragm. "<<Hadron->GetDefinition()->GetParticleName()<<G4endl; 238 613 if ( currentString->GetDecayDirection() > 0 ) 239 614 LeftVector->push_back(Hadron); … … 245 620 }; 246 621 // Split remaining string into 2 final Hadrons ------------------------ 622 //G4cout<<"Split Last"<<G4endl; 247 623 if ( inner_sucess && 248 624 SplitLast(currentString,LeftVector, RightVector) ) … … 282 658 G4ThreeVector PositionOftheStringCreation(theString.GetPosition()); 283 659 660 //G4cout<<"# prod hadrons "<<LeftVector->size()<<G4endl; 284 661 for(size_t C1 = 0; C1 < LeftVector->size(); C1++) 285 662 { 286 663 G4KineticTrack* Hadron = LeftVector->operator[](C1); 287 664 G4LorentzVector Momentum = Hadron->Get4Momentum(); 665 //G4cout<<"Hadron "<<Hadron->GetDefinition()->GetParticleName()<<" "<<Momentum<<G4endl; 288 666 Momentum = toObserverFrame*Momentum; 289 667 Hadron->Set4Momentum(Momentum); … … 306 684 { 307 685 SetMinimalStringMass(string); 308 return sqr(MinimalStringMass + WminLUND) < string->Get4Momentum().mag2(); 686 // return sqr(MinimalStringMass + WminLUND) < string->Get4Momentum().mag2(); 687 return MinimalStringMass < string->Get4Momentum().mag(); // 21.07.2010 309 688 } 310 689 … … 314 693 SetMinimalStringMass(string); 315 694 316 //G4cout<<"StopFragm MinMass "<<MinimalStringMass<<" String Mass "<<std::sqrt(string->Get4Momentum().mag2())<<G4endl; 317 //G4cout<<"WminLUND "<<WminLUND<<" SmoothParam "<<SmoothParam<<" "<<string->Mass()<<G4endl; 695 /* 696 G4cout<<"StopFragm MinMass "<<MinimalStringMass<<" String Mass "<<string->Get4Momentum().mag()<<G4endl; 697 G4cout<<"WminLUND "<<WminLUND<<" SmoothParam "<<SmoothParam<<" "<<string->Mass()<<G4endl; 698 //G4cout<<std::exp(-(string->Mass()*string->Mass()-MinimalStringMass2)/WminLUND/WminLUND)<<G4endl; 699 //G4int Uzhi; G4cin>>Uzhi; 700 */ 701 // 318 702 return (MinimalStringMass + WminLUND)* 319 703 (1 + SmoothParam * (1.-2*G4UniformRand())) > 320 704 string->Mass(); 705 // 706 // return G4UniformRand() < std::exp(-(string->Mass()*string->Mass()-MinimalStringMass2)/WminLUND/WminLUND); 321 707 } 322 708 … … 327 713 { 328 714 //... perform last cluster decay 329 715 //G4cout<<"Split last-----------------------------------------"<<G4endl; 330 716 G4LorentzVector Str4Mom=string->Get4Momentum(); 331 717 332 718 G4ThreeVector ClusterVel =string->Get4Momentum().boostVector(); 333 719 334 G4double ResidualMass = string->Mass(); 335 336 G4ParticleDefinition * LeftHadron, * RightHadron; 337 338 G4int cClusterInterrupt = 0; 339 do 720 G4double StringMass = string->Mass(); 721 G4double StringMassSqr= sqr(StringMass); 722 723 G4ParticleDefinition * LeftHadron(0), * RightHadron(0); 724 G4double LeftHadronMass(0.), RightHadronMass(0.); 725 726 G4ParticleDefinition * FS_LeftHadron[25], * FS_RightHadron[25]; 727 G4double FS_Weight[25]; 728 G4int NumberOf_FS=0; 729 730 for(G4int i=0; i<25; i++) {FS_Weight[i]=0.;} 731 //*********************************************** 732 //G4cout<<"StrMass "<<StringMass<<" q "<<string->GetLeftParton()->GetParticleName()<<" "<<string->GetRightParton()->GetParticleName()<<" StringMassSqr "<<StringMassSqr<<G4endl; 733 734 735 string->SetLeftPartonStable(); // to query quark contents.. 736 737 if (string->FourQuarkString() ) 340 738 { 739 // The string is qq-qqbar type. Diquarks are on the string ends 740 G4int cClusterInterrupt = 0; 741 do 742 { 743 //G4cout<<"cClusterInterrupt "<<cClusterInterrupt<<G4endl; 341 744 if (cClusterInterrupt++ >= ClusterLoopInterrupt) 342 745 { 343 746 return false; 344 747 } 345 G4ParticleDefinition * quark = NULL; 346 string->SetLeftPartonStable(); // to query quark contents.. 347 348 if (!string->FourQuarkString() ) 748 749 G4int LeftQuark1= string->GetLeftParton()->GetPDGEncoding()/1000; 750 G4int LeftQuark2=(string->GetLeftParton()->GetPDGEncoding()/100)%10; 751 752 G4int RightQuark1= string->GetRightParton()->GetPDGEncoding()/1000; 753 G4int RightQuark2=(string->GetRightParton()->GetPDGEncoding()/100)%10; 754 755 if(G4UniformRand()<0.5) 349 756 { 350 // The string is q-qbar, or q-qq, or qbar-qqbar type 351 if (string->DecayIsQuark() && string->StableIsQuark() ) 352 { 353 //... there are quarks on cluster ends 354 LeftHadron= QuarkSplitup(string->GetLeftParton(), quark); 355 } else 356 { 357 //... there is a Diquark on one of the cluster ends 358 G4int IsParticle; 359 360 if ( string->StableIsQuark() ) 361 { 362 IsParticle=(string->GetLeftParton()->GetPDGEncoding()>0) ? -1 : +1; 363 } else 364 { 365 IsParticle=(string->GetLeftParton()->GetPDGEncoding()>0) ? +1 : -1; 366 } 367 368 pDefPair QuarkPair = CreatePartonPair(IsParticle,false); // no diquarks wanted 369 quark = QuarkPair.second; 370 371 LeftHadron=hadronizer->Build(QuarkPair.first, string->GetLeftParton()); 372 } 373 374 RightHadron = hadronizer->Build(string->GetRightParton(), quark); 757 LeftHadron =hadronizer->Build(FindParticle( LeftQuark1), 758 FindParticle(RightQuark1)); 759 RightHadron=hadronizer->Build(FindParticle( LeftQuark2), 760 FindParticle(RightQuark2)); 375 761 } else 376 762 { 377 // The string is qq-qqbar type. Diquarks are on the string ends 378 G4int LiftQuark1= string->GetLeftParton()->GetPDGEncoding()/1000; 379 G4int LiftQuark2=(string->GetLeftParton()->GetPDGEncoding()/100)%10; 380 381 G4int RightQuark1= string->GetRightParton()->GetPDGEncoding()/1000; 382 G4int RightQuark2=(string->GetRightParton()->GetPDGEncoding()/100)%10; 383 384 if(G4UniformRand()<0.5) 385 { 386 LeftHadron =hadronizer->Build(FindParticle( LiftQuark1), 387 FindParticle(RightQuark1)); 388 RightHadron=hadronizer->Build(FindParticle( LiftQuark2), 389 FindParticle(RightQuark2)); 390 } else 391 { 392 LeftHadron =hadronizer->Build(FindParticle( LiftQuark1), 393 FindParticle(RightQuark2)); 394 RightHadron=hadronizer->Build(FindParticle( LiftQuark2), 395 FindParticle(RightQuark1)); 396 } 763 LeftHadron =hadronizer->Build(FindParticle( LeftQuark1), 764 FindParticle(RightQuark2)); 765 RightHadron=hadronizer->Build(FindParticle( LeftQuark2), 766 FindParticle(RightQuark1)); 397 767 } 398 768 399 769 //... repeat procedure, if mass of cluster is too low to produce hadrons 400 770 //... ClusterMassCut = 0.15*GeV model parameter 401 } 402 while (ResidualMass <= LeftHadron->GetPDGMass() + RightHadron->GetPDGMass()); 403 404 //... compute hadron momenta and energies 771 } 772 while ((StringMass <= LeftHadron->GetPDGMass() + RightHadron->GetPDGMass())); 773 774 } // End of if (string->FourQuarkString() ) 775 776 //*********************************************** 777 778 if (!string->FourQuarkString()) 779 { 780 if (string->DecayIsQuark() && string->StableIsQuark() ) 781 {//... there are quarks on cluster ends 782 //G4cout<<"Q Q string"<<G4endl; 783 G4ParticleDefinition * Quark; 784 G4ParticleDefinition * Anti_Quark; 785 786 if(string->GetLeftParton()->GetPDGEncoding()>0) 787 { 788 Quark =string->GetLeftParton(); 789 Anti_Quark=string->GetRightParton(); 790 } else 791 { 792 Quark =string->GetRightParton(); 793 Anti_Quark=string->GetLeftParton(); 794 } 795 796 G4int IDquark =Quark->GetPDGEncoding(); 797 G4int AbsIDquark =std::abs(IDquark); 798 G4int IDanti_quark =Anti_Quark->GetPDGEncoding(); 799 G4int AbsIDanti_quark=std::abs(IDanti_quark); 800 801 NumberOf_FS=0; 802 for(G4int ProdQ=1; ProdQ < 4; ProdQ++) 803 { 804 G4int SignQ=-1; 805 if(IDquark == 2) SignQ= 1; 806 if((IDquark == 1) && (ProdQ == 3)) SignQ= 1; // K0 807 if((IDquark == 3) && (ProdQ == 1)) SignQ=-1; // K0bar 808 if(IDquark == ProdQ) SignQ= 1; 809 810 G4int SignAQ= 1; 811 if(IDanti_quark == -2) SignAQ=-1; 812 if((IDanti_quark ==-1) && (ProdQ == 3)) SignAQ=-1; // K0bar 813 if((IDanti_quark ==-3) && (ProdQ == 1)) SignAQ= 1; // K0 814 if(AbsIDanti_quark == ProdQ) SignAQ= 1; 815 816 G4int StateQ=0; 817 do // while(Meson[AbsIDquark-1][ProdQ-1][StateQ]<>0); 818 { 819 LeftHadron=G4ParticleTable::GetParticleTable()->FindParticle(SignQ* 820 Meson[AbsIDquark-1][ProdQ-1][StateQ]); 821 LeftHadronMass=LeftHadron->GetPDGMass(); 822 StateQ++; 823 824 G4int StateAQ=0; 825 do // while(Meson[AbsIDanti_quark-1][ProdQ-1][StateAQ]<>0); 826 { 827 RightHadron=G4ParticleTable::GetParticleTable()->FindParticle(SignAQ* 828 Meson[AbsIDanti_quark-1][ProdQ-1][StateAQ]); 829 RightHadronMass=RightHadron->GetPDGMass(); 830 StateAQ++; 831 832 if(StringMass >= LeftHadronMass + RightHadronMass) 833 { 834 G4double FS_Psqr=lambda(StringMassSqr,sqr(LeftHadronMass), 835 sqr(RightHadronMass)); 836 FS_Weight[NumberOf_FS]=std::sqrt(FS_Psqr)* 837 MesonWeight[AbsIDquark-1][ProdQ-1][StateQ]* 838 MesonWeight[AbsIDanti_quark-1][ProdQ-1][StateAQ]* 839 Prob_QQbar[ProdQ-1]; 840 841 FS_LeftHadron[NumberOf_FS] = LeftHadron; 842 FS_RightHadron[NumberOf_FS]= RightHadron; 843 NumberOf_FS++; 844 if(NumberOf_FS > 24) 845 {G4int Uzhi; G4cout<<"QQbar string #_FS "<<NumberOf_FS<<G4endl; G4cin>>Uzhi;} 846 } // End of if(StringMass >= LeftHadronMass + RightHadronMass) 847 } while(Meson[AbsIDanti_quark-1][ProdQ-1][StateAQ]!=0); 848 } while(Meson[AbsIDquark-1][ProdQ-1][StateQ]!=0); 849 } // End of for(G4int ProdQ=1; ProdQ < 4; ProdQ++) 850 // 851 } else //--------------------------------------------- 852 { //... there is a Diquark on one of the cluster ends 853 //G4cout<<"DiQ Q string"<<G4endl; 854 G4ParticleDefinition * Di_Quark; 855 G4ParticleDefinition * Quark; 856 857 if(string->GetLeftParton()->GetParticleSubType()== "quark") 858 { 859 Quark =string->GetLeftParton(); 860 Di_Quark=string->GetRightParton(); 861 } else 862 { 863 Quark =string->GetRightParton(); 864 Di_Quark=string->GetLeftParton(); 865 } 866 867 G4int IDquark =Quark->GetPDGEncoding(); 868 G4int AbsIDquark =std::abs(IDquark); 869 G4int IDdi_quark =Di_Quark->GetPDGEncoding(); 870 G4int AbsIDdi_quark=std::abs(IDdi_quark); 871 G4int Di_q1=AbsIDdi_quark/1000; 872 G4int Di_q2=(AbsIDdi_quark-Di_q1*1000)/100; 873 //G4cout<<"IDs "<<IDdi_quark<<" "<<IDquark<<G4endl; 874 875 G4int SignDiQ= 1; 876 if(IDdi_quark < 0) SignDiQ=-1; 877 878 NumberOf_FS=0; 879 for(G4int ProdQ=1; ProdQ < 4; ProdQ++) 880 { 881 G4int SignQ; 882 if(IDquark > 0) 883 { SignQ=-1; 884 if(IDquark == 2) SignQ= 1; 885 if((IDquark == 1) && (ProdQ == 3)) SignQ= 1; // K0 886 if((IDquark == 3) && (ProdQ == 1)) SignQ=-1; // K0bar 887 } else 888 { 889 SignQ= 1; 890 if(IDquark == -2) SignQ=-1; 891 if((IDquark ==-1) && (ProdQ == 3)) SignQ=-1; // K0bar 892 if((IDquark ==-3) && (ProdQ == 1)) SignQ= 1; // K0 893 } 894 895 if(AbsIDquark == ProdQ) SignQ= 1; 896 897 //G4cout<<G4endl; 898 //G4cout<<"-------------------------------------------"<<G4endl; 899 //G4cout<<"Insert QQbar "<<ProdQ<<" Sign "<<SignQ<<G4endl; 900 901 G4int StateQ=0; 902 do // while(Meson[AbsIDquark-1][ProdQ-1][StateQ]<>0); 903 { 904 //G4cout<<G4endl<<"AbsIDquark ProdQ StateQ "<<MesonWeight[AbsIDquark-1][ProdQ-1][StateQ]<<" "<<AbsIDquark<<" "<<ProdQ<<" "<<StateQ<<G4endl; 905 //G4cout<<G4endl<<"AbsIDquark ProdQ StateQ "<<SignQ*Meson[AbsIDquark-1][ProdQ-1][StateQ]<<" "<<AbsIDquark<<" "<<ProdQ<<" "<<StateQ<<G4endl; 906 LeftHadron=G4ParticleTable::GetParticleTable()->FindParticle(SignQ* 907 Meson[AbsIDquark-1][ProdQ-1][StateQ]); 908 LeftHadronMass=LeftHadron->GetPDGMass(); 909 910 //G4cout<<"Meson "<<LeftHadron->GetParticleName()<<G4endl; 911 912 G4int StateDiQ=0; 913 do // while(Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]<>0); 914 { 915 //G4cout<<G4endl<<"Di_q1 Di_q2 ProdQ StateDiQ "<<BaryonWeight[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]<<" "<<Di_q1-1<<" "<<Di_q2-1<<" "<<ProdQ-1<<" "<<StateDiQ<<G4endl; 916 RightHadron=G4ParticleTable::GetParticleTable()->FindParticle(SignDiQ* 917 Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]); 918 RightHadronMass=RightHadron->GetPDGMass(); 919 920 //G4cout<<"Baryon "<<RightHadron->GetParticleName()<<G4endl; 921 922 //G4cout<<"StringMass LeftHadronMass RightHadronMass "<<StringMass<<" "<<LeftHadronMass<<" "<< RightHadronMass<<G4endl; 923 924 if(StringMass >= LeftHadronMass + RightHadronMass) 925 { 926 G4double FS_Psqr=lambda(StringMassSqr,sqr(LeftHadronMass), 927 sqr(RightHadronMass)); 928 FS_Weight[NumberOf_FS]=std::sqrt(FS_Psqr)* 929 MesonWeight[AbsIDquark-1][ProdQ-1][StateQ]* 930 BaryonWeight[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]* 931 Prob_QQbar[ProdQ-1]; 932 933 FS_LeftHadron[NumberOf_FS] = LeftHadron; 934 FS_RightHadron[NumberOf_FS]= RightHadron; 935 936 //G4cout<<"State "<<NumberOf_FS<<" "<<std::sqrt(FS_Psqr/4./StringMassSqr)<<" "<<std::sqrt(FS_Psqr)<<" "<<FS_Weight[NumberOf_FS]<<G4endl; 937 //G4cout<<"++++++++++++++++++++++++++++++++"<<G4endl; 938 NumberOf_FS++; 939 940 if(NumberOf_FS > 24) 941 {G4int Uzhi; G4cout<<"QQbar string #_FS "<<NumberOf_FS<<G4endl; G4cin>>Uzhi;} 942 } // End of if(StringMass >= LeftHadronMass + RightHadronMass) 943 944 StateDiQ++; 945 //G4cout<<Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]<<" "<<Di_q1-1<<" "<<Di_q2-1<<" "<<ProdQ-1<<" "<<StateDiQ<<G4endl; 946 } while(Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]!=0); 947 948 StateQ++; 949 } while(Meson[AbsIDquark-1][ProdQ-1][StateQ]!=0); 950 } // End of for(G4int ProdQ=1; ProdQ < 4; ProdQ++) 951 } 952 //==================================== 953 954 if(NumberOf_FS == 0) return false; 955 //G4cout<<"NumberOf_FS "<<NumberOf_FS<<G4endl; 956 G4double SumWeights=0.; 957 for(G4int i=0; i<NumberOf_FS; i++) {SumWeights+=FS_Weight[i];} 958 959 G4double ksi=G4UniformRand(); 960 G4double Sum=0.; 961 G4int SampledState(0); 962 963 for(G4int i=0; i<NumberOf_FS; i++) 964 { 965 Sum+=(FS_Weight[i]/SumWeights); 966 SampledState=i; 967 if(Sum >= ksi) break; 968 } 969 970 LeftHadron =FS_LeftHadron[SampledState]; 971 RightHadron=FS_RightHadron[SampledState]; 972 973 //G4cout<<"Selected "<<SampledState<<" "<<LeftHadron->GetParticleName()<<" "<<RightHadron->GetParticleName()<<G4endl; 974 975 } // End of if(!string->FourQuarkString()) 405 976 406 977 G4LorentzVector LeftMom, RightMom; … … 409 980 Sample4Momentum(&LeftMom, LeftHadron->GetPDGMass(), 410 981 &RightMom, RightHadron->GetPDGMass(), 411 ResidualMass);982 StringMass); 412 983 413 984 LeftMom.boost(ClusterVel); … … 429 1000 G4double MassMt2, AntiMassMt2; 430 1001 G4double AvailablePz, AvailablePz2; 431 432 if(Mass > 930. || AntiMass > 930.) // If there is a baryon 1002 1003 //G4cout<<"Masses "<<InitialMass<<" "<<Mass<<" "<<AntiMass<<G4endl; 1004 // 1005 1006 if((Mass > 930. || AntiMass > 930.)) //If there is a baryon 433 1007 { 434 // ----------------- Isotripic decay ------------------------------------1008 // ----------------- Isotropic decay ------------------------------------ 435 1009 G4double r_val = sqr(InitialMass*InitialMass - Mass*Mass - AntiMass*AntiMass) - 436 1010 sqr(2.*Mass*AntiMass); 437 1011 G4double Pabs = (r_val > 0.)? std::sqrt(r_val)/(2.*InitialMass) : 0; 1012 //G4cout<<"P for isotr decay "<<Pabs<<G4endl; 438 1013 439 1014 //... sample unit vector 440 G4double pz = 1015 G4double pz =1. - 2.*G4UniformRand(); 441 1016 G4double st = std::sqrt(1. - pz * pz)*Pabs; 442 1017 G4double phi = 2.*pi*G4UniformRand(); … … 450 1025 AntiMom->setPx(-px); AntiMom->setPy(-py); AntiMom->setPz(-pz); 451 1026 AntiMom->setE (std::sqrt(Pabs*Pabs + AntiMass*AntiMass)); 1027 //G4int Uzhi; G4cin>>Uzhi; 452 1028 } 453 1029 else 1030 // 454 1031 { 455 1032 do … … 461 1038 G4double termN = 2*termD + 4*Mass*Mass + 4*AntiMass*AntiMass; 462 1039 G4double pt2max=(termD*termD - termab )/ termN ; 1040 //G4cout<<"Anis "<<pt2max<<" "<<(termD*termD-termab)/(4.*InitialMass*InitialMass)<<G4endl; 463 1041 464 1042 Pt=SampleQuarkPt(std::sqrt(pt2max)); Pt.setZ(0); G4double Pt2=Pt.mag2(); 465 1043 //G4cout<<"Sampl pt2 "<<Pt2<<G4endl; 466 1044 MassMt2 = Mass * Mass + Pt2; 467 1045 AntiMassMt2= AntiMass * AntiMass + Pt2; … … 491 1069 G4FragmentingString * string, G4FragmentingString * newString) 492 1070 { 1071 //G4cout<<"Start SplitEandP "<<G4endl; 493 1072 G4LorentzVector String4Momentum=string->Get4Momentum(); 494 1073 G4double StringMT2=string->Get4Momentum().mt2(); … … 496 1075 G4double HadronMass = pHadron->GetPDGMass(); 497 1076 498 SetMinimalStringMass(newString); 1077 SetMinimalStringMass(newString); 1078 //G4cout<<"HadM MinimalStringMassLeft StringM "<<HadronMass<<" "<<MinimalStringMass<<" "<<String4Momentum.mag()<<G4endl; 1079 1080 if(HadronMass + MinimalStringMass > String4Momentum.mag()) {return 0;}// have to start all over! 499 1081 String4Momentum.setPz(0.); 500 1082 G4ThreeVector StringPt=String4Momentum.vect(); … … 517 1099 G4double Pz2 = (sqr(StringMT2 - HadronMassT2 - ResidualMassT2) - 518 1100 4*HadronMassT2 * ResidualMassT2)/4./StringMT2; 519 1101 //G4cout<<"Pz2 "<<Pz2<<G4endl; 520 1102 if(Pz2 < 0 ) {return 0;} // have to start all over! 521 1103 … … 526 1108 G4double zMax = (std::sqrt(HadronMassT2+Pz2) + Pz)/std::sqrt(StringMT2); 527 1109 1110 //G4cout<<"if (zMin >= zMax) return 0 "<<zMin<<" "<<zMax<<G4endl; 528 1111 if (zMin >= zMax) return 0; // have to start all over! 529 1112 … … 531 1114 string->GetDecayParton()->GetPDGEncoding(), pHadron, 532 1115 HadronPt.x(), HadronPt.y()); 533 1116 //G4cout<<"z "<<z<<G4endl; 534 1117 //... now compute hadron longitudinal momentum and energy 535 1118 // longitudinal hadron momentum component HadronPz … … 543 1126 544 1127 G4LorentzVector * a4Momentum= new G4LorentzVector(HadronPt,HadronE); 545 1128 //G4cout<<"Out SplitEandP "<<G4endl; 546 1129 return a4Momentum; 547 1130 } … … 596 1179 return z; 597 1180 } 1181 1182 //------------------------------------------------------------------------ 1183 G4double G4LundStringFragmentation::lambda(G4double s, G4double m1_Sqr, G4double m2_Sqr) 1184 { 1185 G4double lam = sqr(s - m1_Sqr - m2_Sqr) - 4.*m1_Sqr*m2_Sqr; 1186 return lam; 1187 } 1188 -
trunk/source/processes/hadronic/models/parton_string/hadronization/src/G4QGSMFragmentation.cc
r1337 r1340 25 25 // 26 26 // 27 // $Id: G4QGSMFragmentation.cc,v 1. 9 2008/06/23 08:35:55vuzhinsk Exp $28 // GEANT4 tag $Name: geant4-09-0 4-beta-01$27 // $Id: G4QGSMFragmentation.cc,v 1.10 2010/09/20 12:46:23 vuzhinsk Exp $ 28 // GEANT4 tag $Name: geant4-09-03-ref-09 $ 29 29 // 30 30 // ----------------------------------------------------------------------------- -
trunk/source/processes/hadronic/models/parton_string/hadronization/src/G4VKinkyStringDecay.cc
r1337 r1340 25 25 // 26 26 // 27 // $Id: G4VKinkyStringDecay.cc,v 1. 4 2008/04/25 14:20:14vuzhinsk Exp $28 // GEANT4 tag $Name: geant4-09-0 4-beta-01$27 // $Id: G4VKinkyStringDecay.cc,v 1.5 2010/09/20 12:46:23 vuzhinsk Exp $ 28 // GEANT4 tag $Name: geant4-09-03-ref-09 $ 29 29 // Maxim Komogorov 30 30 // -
trunk/source/processes/hadronic/models/parton_string/hadronization/src/G4VLongitudinalStringDecay.cc
r1337 r1340 25 25 // 26 26 // 27 // $Id: G4VLongitudinalStringDecay.cc,v 1.2 0 2010/06/21 17:50:48vuzhinsk Exp $28 // GEANT4 tag $Name: geant4-09-0 4-beta-01$27 // $Id: G4VLongitudinalStringDecay.cc,v 1.22 2010/09/20 12:46:23 vuzhinsk Exp $ 28 // GEANT4 tag $Name: geant4-09-03-ref-09 $ 29 29 // 30 30 // ----------------------------------------------------------------------------- … … 71 71 72 72 // Changable Parameters below. 73 SigmaQT = 0.5 * GeV; 73 SigmaQT = 0.5 * GeV; // 0.5 74 74 75 75 StrangeSuppress = 0.44; // 27 % strange quarks produced, ie. u:d:s=1:1:0.27 … … 250 250 Hadron1 = (minMassHadronizer->*build)(string->GetLeftParton(), 251 251 string->GetRightParton()); 252 //G4cout<<"Hadron1 "<<Hadron1->GetParticleName()<<G4endl; 252 253 mass= (Hadron1)->GetPDGMass(); 253 254 } else … … 323 324 G4FragmentingString *&newString) 324 325 { 326 //G4cout<<"Start SplitUP"<<G4endl; 325 327 //... random choice of string end to use for creating the hadron (decay) 326 328 SideOfDecay = (G4UniformRand() < 0.5)? 1: -1; … … 341 343 HadronDefinition= DiQuarkSplitup(string->GetDecayParton(), newStringEnd); 342 344 } 345 346 //G4cout<<"New had "<<HadronDefinition->GetParticleName()<<G4endl; 343 347 // create new String from old, ie. keep Left and Right order, but replace decay 344 348 345 349 newString=new G4FragmentingString(*string,newStringEnd); // To store possible 346 350 // quark containt of new string 351 //G4cout<<"SplitEandP "<<G4endl; 347 352 G4LorentzVector* HadronMomentum=SplitEandP(HadronDefinition, string, newString); 348 353 … … 360 365 delete HadronMomentum; 361 366 } 367 //G4cout<<"End SplitUP"<<G4endl; 362 368 return Hadron; 363 369 }
Note: See TracChangeset
for help on using the changeset viewer.