- Timestamp:
- Jan 8, 2010, 11:56:51 AM (15 years ago)
- Location:
- trunk/source/processes/hadronic/models/parton_string/diffraction
- Files:
-
- 11 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/source/processes/hadronic/models/parton_string/diffraction/History
r1196 r1228 1 $Id: History,v 1.2 2 2009/10/29 14:58:52vuzhinsk Exp $1 $Id: History,v 1.29 2009/12/15 19:20:46 vuzhinsk Exp $ 2 2 ------------------------------------------------------------------- 3 3 … … 15 15 * Please list in reverse chronological order (last date on top) 16 16 --------------------------------------------------------------- 17 15 Dec 2009 V.Uzhinsky (hadr-string-diff-V09-02-23) 18 Bug was fixed for Kaon + A interactions in G4DiffractiveExcitation.cc 19 20 ------------------------------------------------------------ 21 15 Dec 2009 G.Folger (hadr-string-diff-V09-02-22) 22 - G4FTFModel: In ctor, initialise NumberOfInvolvedNucleon 23 24 14 Dec. 2009 V. Uzhinsky (hadr-string-diff-V09-02-21) 25 Momentum sampling in the reggeon cascade is improved. 26 More hard restrictions on the momentum are introdused. 27 Infinite loop is erased. 28 Arithmetical bag in kinky string is fixed. 29 30 -------------------------------------------------------- 31 10 Dec. 2009 V. Uzhinsky (hadr-string-diff-V09-02-20) 32 Crush found by Alberto with AveragePt2=0 is erased in FTF. 33 Re-tag of hadr-string-diff-V09-03-00 34 ----------------------------------------------------- 35 9 Dec. 2009 V. Uzhinsky (hadr-string-diff-V09-03-00) 36 Crush found by Alberto with AveragePt2=0 is erased in FTF. 37 -------------------------------------------------------- 38 6 Dec. 2009 V.Uzhinsky (hadr-string-diff-V09-02-19) 39 Sampling of kinky string momentum is improved 40 41 25 November 2009, V. Uzhinsky (hadr-string-diff-V09-02-18) 42 Memory leak is erased in FTF at the end. 43 17 44 29 October 2009, V. Uzhinsky (hadr-string-diff-V09-02-17) 18 45 Warning meassage is erased for using slc5 compilator in FTFModel.cc -
trunk/source/processes/hadronic/models/parton_string/diffraction/include/G4DiffractiveSplitableHadron.hh
r1196 r1228 26 26 // 27 27 // $Id: G4DiffractiveSplitableHadron.hh,v 1.5 2009/08/03 13:14:19 vuzhinsk Exp $ 28 // GEANT4 tag $Name: geant4-09-03 -cand-01$28 // GEANT4 tag $Name: geant4-09-03 $ 29 29 // 30 30 -
trunk/source/processes/hadronic/models/parton_string/diffraction/include/G4FTFModel.hh
r1196 r1228 26 26 // 27 27 // $Id: G4FTFModel.hh,v 1.10 2009/10/25 10:50:54 vuzhinsk Exp $ 28 // GEANT4 tag $Name: geant4-09-03 -cand-01$28 // GEANT4 tag $Name: geant4-09-03 $ 29 29 // 30 30 // Class Description -
trunk/source/processes/hadronic/models/parton_string/diffraction/include/G4FTFParameters.hh
r1196 r1228 28 28 // 29 29 // $Id: G4FTFParameters.hh,v 1.7 2009/10/25 10:50:54 vuzhinsk Exp $ 30 // GEANT4 tag $Name: geant4-09-03 -cand-01$30 // GEANT4 tag $Name: geant4-09-03 $ 31 31 // 32 32 #include "G4Proton.hh" -
trunk/source/processes/hadronic/models/parton_string/diffraction/include/G4FTFParticipants.hh
r1196 r1228 26 26 // 27 27 // $Id: G4FTFParticipants.hh,v 1.6 2009/08/03 13:14:19 vuzhinsk Exp $ 28 // GEANT4 tag $Name: geant4-09-03 -cand-01$28 // GEANT4 tag $Name: geant4-09-03 $ 29 29 // 30 30 -
trunk/source/processes/hadronic/models/parton_string/diffraction/src/G4DiffractiveExcitation.cc
r1196 r1228 25 25 // 26 26 // 27 // $Id: G4DiffractiveExcitation.cc,v 1. 16 2009/10/06 10:10:36vuzhinsk Exp $27 // $Id: G4DiffractiveExcitation.cc,v 1.21 2009/12/15 19:14:31 vuzhinsk Exp $ 28 28 // ------------------------------------------------------------ 29 29 // GEANT 4 class implemetation file … … 73 73 G4VSplitableHadron *target, 74 74 G4FTFParameters *theParameters, 75 G4ElasticHNScattering *theElastic) const // Uzhi 03.09.0975 G4ElasticHNScattering *theElastic) const 76 76 { 77 77 // -------------------- Projectile parameters ----------------------- 78 78 G4LorentzVector Pprojectile=projectile->Get4Momentum(); 79 //G4cout<<"Excite P "<<Pprojectile<<G4endl;80 79 81 80 if(Pprojectile.z() < 0.) … … 88 87 89 88 G4int ProjectilePDGcode=projectile->GetDefinition()->GetPDGEncoding(); 90 // G4ParticleDefinition * ProjectileDefinition = projectile->GetDefinition();91 89 G4int absProjectilePDGcode=std::abs(ProjectilePDGcode); 92 90 … … 115 113 G4double M0target = Ptarget.mag(); 116 114 117 //G4cout<<"Excite T "<<Ptarget<<G4endl;118 //G4cout<<"PDGs "<<ProjectilePDGcode<<" "<<TargetPDGcode<<G4endl;119 //G4cout<<"Masses "<<M0projectile<<" "<<M0target<<G4endl;120 121 115 G4double TargetRapidity = Ptarget.rapidity(); 122 116 … … 149 143 150 144 G4LorentzVector Ptmp=toCms*Pprojectile; 145 151 146 if ( Ptmp.pz() <= 0. ) 152 147 { … … 167 162 168 163 G4double SqrtS=std::sqrt(S); 169 //G4cout<<" SqrtS "<<SqrtS<<" 2200 "<<G4endl; 170 if(absProjectilePDGcode > 1000 && SqrtS < 2300*MeV && SqrtS < SumMasses)164 165 if(absProjectilePDGcode > 1000 && (SqrtS < 2300*MeV || SqrtS < SumMasses)) 171 166 {target->SetStatus(2); return false;} // The model cannot work for 172 167 // p+p-interactions … … 174 169 175 170 if(( absProjectilePDGcode == 211 || ProjectilePDGcode == 111) && 176 ( SqrtS < 1600*MeV) && (SqrtS < SumMasses))171 ((SqrtS < 1600*MeV) || (SqrtS < SumMasses))) 177 172 {target->SetStatus(2); return false;} // The model cannot work for 178 173 // Pi+p-interactions … … 181 176 if(( (absProjectilePDGcode == 321) || (ProjectilePDGcode == -311) || 182 177 (absProjectilePDGcode == 311) || (absProjectilePDGcode == 130) || 183 (absProjectilePDGcode == 310)) && ( SqrtS < 1600*MeV) && (SqrtS < SumMasses))178 (absProjectilePDGcode == 310)) && ((SqrtS < 1600*MeV) || (SqrtS < SumMasses))) 184 179 {target->SetStatus(2); return false;} // The model cannot work for 185 180 // K+p-interactions … … 223 218 // Charge exchange can be possible for baryons ----------------- 224 219 225 //G4cout<<1./std::exp(-2.0*(ProjectileRapidity - TargetRapidity))<<G4endl;226 //G4int Uzhi; G4cin>>Uzhi;227 220 // Getting the values needed for exchange ---------------------- 228 221 G4double MagQuarkExchange =theParameters->GetMagQuarkExchange(); … … 231 224 232 225 // G4double NucleonMass= 233 (G4ParticleTable::GetParticleTable()->FindParticle(2112))->GetPDGMass();226 // (G4ParticleTable::GetParticleTable()->FindParticle(2112))->GetPDGMass(); 234 227 G4double DeltaMass= 235 228 (G4ParticleTable::GetParticleTable()->FindParticle(2224))->GetPDGMass(); … … 239 232 std::exp(-SlopeQuarkExchange*(ProjectileRapidity - TargetRapidity))) 240 233 { 241 242 //G4cout<<"Exchange "<<G4endl;243 244 234 G4int NewProjCode(0), NewTargCode(0); 245 235 246 //G4bool StandardExcitation(false); // =================================247 236 G4int ProjQ1(0), ProjQ2(0), ProjQ3(0); 248 237 … … 264 253 G4int TargExchangeQ(0); 265 254 266 //G4cout<<"Targ "<<TargQ1<<" "<<TargQ2<<" "<<TargQ3<<G4endl;267 255 if(absProjectilePDGcode < 1000 ) 268 256 { // projectile is meson ----------------- 269 //G4cout<<"Proj Q1 Q2 "<<ProjQ1<<" "<<ProjQ2<<G4endl;270 271 257 if(ProjQ1 > 0 ) // ProjQ1 is quark 272 258 { … … 283 269 TargExchangeQ = TargQ3; TargQ3=ProjExchangeQ; ProjQ1=TargExchangeQ; 284 270 } 285 //G4cout<<" Pr Tr "<<ProjQ1<<" "<<TargQ1<<" "<<TargQ2<<" "<<TargQ3<<G4endl;286 271 } else // ProjQ2 is quark 287 272 { … … 310 295 } 311 296 312 // projectile->SetDefinition(313 // G4ParticleTable::GetParticleTable()->FindParticle(NewProjCode));314 315 297 NewTargCode = NewNucleonId(TargQ1, TargQ2, TargQ3); 316 298 317 299 if( (TargQ1 == TargQ2) && (TargQ1 == TargQ3) && 318 300 (SqrtS > M0projectile+DeltaMass)) {NewTargCode +=2;} //Create Delta isobar … … 326 308 (SqrtS > M0projectile+DeltaMass)) {NewTargCode +=2;} //Create Delta isobar 327 309 else {} //Save initial nucleon 328 329 target->SetDefinition( 330 G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode)); 331 //G4cout<<"New PDGs "<<NewProjCode<<" "<<NewTargCode<<G4endl; 332 //G4int Uzhi; G4cin>>Uzhi; 333 //} 310 // target->SetDefinition( // Fix 15.12.09 311 // G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode));// Fix 15.12.09 334 312 } else 335 313 { // projectile is baryon ---------------- 336 G4double Same=0.; //0.3; //0.5;314 G4double Same=0.; //0.3; //0.5; 337 315 G4bool ProjDeltaHasCreated(false); 338 316 G4bool TargDeltaHasCreated(false); … … 341 319 if(G4UniformRand() < 0.5) // Sampling exchange quark from proj. or targ. 342 320 { // Sampling exchanged quark from the projectile --- 343 //G4cout<<"Proj Exc"<<G4endl;344 321 if( Ksi < 0.333333 ) 345 322 {ProjExchangeQ = ProjQ1;} … … 348 325 else 349 326 {ProjExchangeQ = ProjQ3;} 350 //G4cout<<"ProjExchangeQ "<<ProjExchangeQ<<G4endl; 351 352 if((ProjExchangeQ != TargQ1)||(G4UniformRand()<Same)) // Vova 327 328 if((ProjExchangeQ != TargQ1)||(G4UniformRand()<Same)) 353 329 { 354 330 TargExchangeQ = TargQ1; TargQ1=ProjExchangeQ; ProjExchangeQ=TargExchangeQ; 355 331 } else 356 if((ProjExchangeQ != TargQ2)||(G4UniformRand()<Same)) // Vova332 if((ProjExchangeQ != TargQ2)||(G4UniformRand()<Same)) 357 333 { 358 334 TargExchangeQ = TargQ2; TargQ2=ProjExchangeQ; ProjExchangeQ=TargExchangeQ; … … 361 337 TargExchangeQ = TargQ3; TargQ3=ProjExchangeQ; ProjExchangeQ=TargExchangeQ; 362 338 } 363 //G4cout<<"ProjExchangeQ "<<ProjExchangeQ<<G4endl; 339 364 340 if( Ksi < 0.333333 ) 365 341 {ProjQ1=ProjExchangeQ;} … … 369 345 {ProjQ3=ProjExchangeQ;} 370 346 371 /*372 NewProjCode = NewNucleonId(ProjQ1, ProjQ2, ProjQ3); // *****************************373 374 if((ProjQ1==ProjQ2) && (ProjQ1==ProjQ3)) {NewProjCode +=2; ProjDeltaHasCreated=true;}375 else if(projectile->GetDefinition()->GetPDGiIsospin() == 3)// Projectile was Delta376 { if(G4UniformRand() > DeltaProbAtQuarkExchange)377 {NewProjCode +=2; ProjDeltaHasCreated=true;}378 else {NewProjCode +=0; ProjDeltaHasCreated=false;}379 }380 else // Projectile was Nucleon381 {382 if((G4UniformRand() < DeltaProbAtQuarkExchange) && (SqrtS > DeltaMass+M0target))383 {NewProjCode +=2; ProjDeltaHasCreated=true;}384 else {NewProjCode +=0; ProjDeltaHasCreated=false;}385 }386 387 NewTargCode = NewNucleonId(TargQ1, TargQ2, TargQ3); // *****************************388 389 if((TargQ1==TargQ2) && (TargQ1==TargQ3)) {NewTargCode +=2; TargDeltaHasCreated=true;}390 else if(target->GetDefinition()->GetPDGiIsospin() == 3) // Target was Delta391 { if(G4UniformRand() > DeltaProbAtQuarkExchange)392 {NewTargCode +=2; TargDeltaHasCreated=true;}393 else {NewTargCode +=0; TargDeltaHasCreated=false;}394 }395 else // Target was Nucleon396 {397 if((G4UniformRand() < DeltaProbAtQuarkExchange) && (SqrtS > M0projectile+DeltaMass))398 {NewTargCode +=2; TargDeltaHasCreated=true;}399 else {NewTargCode +=0; TargDeltaHasCreated=false;}400 }401 402 if((absProjectilePDGcode == NewProjCode) && (absTargetPDGcode == NewTargCode))403 { // Nothing was changed! It is not right!?404 }405 */406 /*-----------------------------------------------------------------------------------------407 if(!ProjDeltaHasCreated && !TargDeltaHasCreated) // No Deltas were created then408 { if( G4UniformRand() < 0.5) {NewProjCode +=2; ProjDeltaHasCreated = true;}409 else {NewTargCode +=2; TargDeltaHasCreated = true;}410 }411 412 if(!(ProjDeltaHasCreated && TargDeltaHasCreated))413 {414 if((G4UniformRand() < DeltaProbAtQuarkExchange) && (SqrtS >2600.))415 {416 G4cout<<"2 Deltas "<<G4endl;417 if(!ProjDeltaHasCreated)418 {NewProjCode +=2; ProjDeltaHasCreated = true;}419 else {NewTargCode +=2; TargDeltaHasCreated = true;} // For Delta isobar420 }421 }422 */423 347 } else 424 348 { // Sampling exchanged quark from the target ------- 425 //G4cout<<"Targ Exc"<<G4endl;426 349 if( Ksi < 0.333333 ) 427 350 {TargExchangeQ = TargQ1;} … … 430 353 else 431 354 {TargExchangeQ = TargQ3;} 432 //G4cout<<"TargExchangeQ "<<TargExchangeQ<<G4endl; 433 if((TargExchangeQ != ProjQ1)||(G4UniformRand()<Same)) // Vova355 356 if((TargExchangeQ != ProjQ1)||(G4UniformRand()<Same)) 434 357 { 435 358 ProjExchangeQ = ProjQ1; ProjQ1=TargExchangeQ; TargExchangeQ=ProjExchangeQ; 436 359 } else 437 if((TargExchangeQ != ProjQ2)||(G4UniformRand()<Same)) // Vova360 if((TargExchangeQ != ProjQ2)||(G4UniformRand()<Same)) 438 361 { 439 362 ProjExchangeQ = ProjQ2; ProjQ2=TargExchangeQ; TargExchangeQ=ProjExchangeQ; … … 442 365 ProjExchangeQ = ProjQ3; ProjQ3=TargExchangeQ; TargExchangeQ=ProjExchangeQ; 443 366 } 444 //G4cout<<"TargExchangeQ "<<TargExchangeQ<<G4endl; 367 445 368 if( Ksi < 0.333333 ) 446 369 {TargQ1=TargExchangeQ;} … … 449 372 else 450 373 {TargQ3=TargExchangeQ;} 451 /* 452 NewProjCode = NewNucleonId(ProjQ1, ProjQ2, ProjQ3); 453 if( (ProjQ1 == ProjQ2) && (ProjQ1 == ProjQ3) ) 454 {NewProjCode +=2; ProjDeltaHasCreated = true;} 455 456 NewTargCode = NewNucleonId(TargQ1, TargQ2, TargQ3); 457 if( (TargQ1 == TargQ2) && (TargQ1 == TargQ3) ) 458 {NewTargCode +=2; TargDeltaHasCreated = true;} 459 460 if(!ProjDeltaHasCreated && !TargDeltaHasCreated) 461 {NewTargCode +=2; TargDeltaHasCreated = true;} 462 463 if(!(ProjDeltaHasCreated && TargDeltaHasCreated)) 464 { 465 if((G4UniformRand() < DeltaProbAtQuarkExchange) && (SqrtS >2600.)) 466 { 467 G4cout<<"2 Deltas "<<G4endl; 468 if(!ProjDeltaHasCreated) 469 {NewProjCode +=2; ProjDeltaHasCreated = true;} 470 else {NewTargCode +=2; TargDeltaHasCreated = true;} 471 } 472 } 473 */ 374 474 375 } // End of sampling baryon 475 376 … … 509 410 // Forming baryons -------------------------------------------------- 510 411 511 // if( ProjDeltaHasCreated && TargDeltaHasCreated ) {StandardExcitation = true;}512 513 //G4cout<<"New PDG "<<NewProjCode<<" "<<NewTargCode<<G4endl;514 //G4int Uzhi; G4cin>>Uzhi;515 412 } // End of if projectile is baryon --------------------------- 516 413 … … 518 415 // If we assume that final state hadrons after the charge exchange will be 519 416 // in the ground states, we have to put ---------------------------------- 417 520 418 if(M0projectile < 521 419 (G4ParticleTable::GetParticleTable()->FindParticle(NewProjCode))->GetPDGMass()) … … 537 435 2*S*M0projectile2 - 2*S*M0target2 - 2*M0projectile2*M0target2) 538 436 /4./S; 539 //G4cout<<"New Masses "<<M0projectile<<" "<<M0target<<G4endl;540 //G4cout<<"PZcms2 "<<PZcms2<<G4endl;541 437 542 438 if(PZcms2 < 0) {return false;} // It can be if energy is not sufficient for Delta … … 556 452 Ptarget.setE(std::sqrt(M0target2+PZcms2)); 557 453 558 //G4cout<<"Proj "<<Pprojectile<<" "<<Pprojectile.mag()<<G4endl;559 //G4cout<<"Targ "<<Ptarget<<" "<<Ptarget.mag()<<G4endl;560 561 // if(!StandardExcitation)562 454 { 563 455 Pprojectile.transform(toLab); … … 571 463 572 464 G4bool Result= theElastic->ElasticScattering (projectile,target,theParameters); 573 //G4cout<<"1 Delta result "<<Result<<"**********"<<G4endl; 465 574 466 return Result; 575 467 } 576 /*577 else578 {579 if(M0projectile > ProjectileDiffStateMinMass)580 { ProjectileDiffStateMinMass +=200.*MeV; ProjectileNonDiffStateMinMass +=200.*MeV;}581 582 if(M0target > TargetDiffStateMinMass)583 { TargetDiffStateMinMass +=200.*MeV; TargetNonDiffStateMinMass +=200.*MeV;}584 585 ProbOfDiffraction = 1.;586 ProbProjectileDiffraction =0.5;587 ProbTargetDiffraction =0.5;588 G4cout<<"Exc DD "<<M0projectile<<" "<<M0target<<" --------------------"<<G4endl;589 // After that standard FTF excitation590 }591 */592 468 } // End of charge exchange part ------------------------------ 593 469 594 470 // ------------------------------------------------------------------ 595 // G4double ProbOfDiffraction=ProbProjectileDiffraction +596 // ProbTargetDiffraction;597 //G4cout<<ProbOfDiffraction<<" "<<ProbProjectileDiffraction<<" "<<ProbTargetDiffraction<<G4endl;598 //G4int Uzhi; G4cin>>Uzhi;599 471 if(ProbOfDiffraction!=0.) 600 472 { … … 606 478 } 607 479 608 //ProbOfDiffraction=1.; //0.5; // Vova609 610 480 G4double ProjectileDiffStateMinMass2 = ProjectileDiffStateMinMass * 611 481 ProjectileDiffStateMinMass; … … 617 487 G4double TargetNonDiffStateMinMass2 = TargetNonDiffStateMinMass * 618 488 TargetNonDiffStateMinMass; 619 620 489 621 490 G4double Pt2; … … 637 506 if(G4UniformRand() < ProbProjectileDiffraction) 638 507 { //-------- projectile diffraction --------------- 639 //G4cout<<"Proj difr"<<G4endl;640 508 do { 641 509 // Generate pt … … 659 527 2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2) 660 528 /4./S; 661 //G4cout<<"PZcms2 < 0 false "<<PZcms2<<G4endl; 529 662 530 if(PZcms2 < 0 ) 663 531 { … … 703 571 { // -------------- Target diffraction ---------------- 704 572 do { 705 //G4cout<<"Targ difr"<<G4endl;706 573 // Generate pt 707 574 // if (whilecount++ >= 500 && (whilecount%100)==0) … … 770 637 else //----------- Non-diffraction process ------------ 771 638 { 772 //G4cout<<"Non difr"<<G4endl;773 639 do { 774 640 // Generate pt … … 841 707 Ptarget -= Qmomentum; 842 708 843 //G4cout<<"Masses "<<Pprojectile.mag()<<" "<<Ptarget.mag()<<G4endl;844 //G4int Uzhi; G4cin>>Uzhi;845 846 709 // Transform back and update SplitableHadron Participant. 847 710 Pprojectile.transform(toLab); … … 865 728 866 729 // --------------------------------------------------------------------- 867 //G4ExcitedString * G4DiffractiveExcitation::868 // String(G4VSplitableHadron * hadron, G4bool isProjectile) const869 730 void G4DiffractiveExcitation::CreateStrings(G4VSplitableHadron * hadron, 870 731 G4bool isProjectile, … … 886 747 return; 887 748 } 749 888 750 G4LorentzVector Phadron=hadron->Get4Momentum(); 889 /* 890 G4cout<<"Create strings had "<<hadron->GetDefinition()->GetParticleName()<<" "<<Phadron<<" "<<Phadron.mag()<<G4endl; 891 G4cout<<"isProjectile "<<isProjectile<<G4endl; 892 G4cout<<"start Q "<<start->GetDefinition()->GetPDGEncoding()<<G4endl; 893 G4cout<<"end Q "<<end->GetDefinition()->GetPDGEncoding()<<G4endl; 894 */ 751 895 752 G4LorentzVector Pstart(0.,0.,0.,0.); 896 753 G4LorentzVector Pend(0.,0.,0.,0.); … … 914 771 G4double W = hadron->Get4Momentum().mag(); 915 772 G4double W2=W*W; 916 //G4cout<<"W Wmin "<<W<<" "<<Wmin<<G4endl; 773 917 774 G4double Pt(0.), x1(0.), x2(0.), x3(0.); 918 775 G4bool Kink=false; … … 922 779 G4double Pt2kink=theParameters->GetPt2Kink(); 923 780 Pt = std::sqrt(Pt2kink*(std::pow(W2/16./Pt2kink+1.,G4UniformRand()) - 1.)); 924 // Pt=0.; 925 //G4cout<<"Pt2kink Pt Pt2 "<<Pt2kink<<" "<<Pt<<" "<<Pt*Pt<<G4endl; 781 926 782 if(Pt > 500.*MeV) 927 783 { 928 784 G4double Ymax = std::log(W/2./Pt + std::sqrt(W2/4./Pt/Pt - 1.)); 929 785 G4double Y=Ymax*(1.- 2.*G4UniformRand()); 930 //G4cout<<"Ymax Y "<<Ymax<<" "<<Y<<G4endl; 786 931 787 x1=1.-Pt/W*std::exp( Y); 932 788 x3=1.-Pt/W*std::exp(-Y); 933 789 x2=2.-x1-x3; 934 //G4cout<<"X1 X2 X3 "<<x1<<" "<<x2<<" "<<x3<<G4endl; 790 935 791 G4double Mass_startQ = 650.*MeV; 936 792 if(PDGcode_startQ < 3) Mass_startQ = 325.*MeV; … … 947 803 948 804 G4double P2_2=sqr((2.-x1-x3)*W/2.); 949 //G4cout<<"P2_1 P2_2 P2_3 "<<P2_1<<" "<<P2_2<<" "<<P2_3<<G4endl; 950 if((P2_1 < 0.) || (P2_3 <0.))805 806 if((P2_1 <= 0.) || (P2_3 <= 0.)) 951 807 { Kink=false;} 952 808 else … … 955 811 G4double P_2=std::sqrt(P2_2); 956 812 G4double P_3=std::sqrt(P2_3); 957 //G4cout<<"P_1 P_2 P_3 "<<P_1<<" "<<P_2<<" "<<P_3<<G4endl; 813 958 814 G4double CosT12=(P2_3-P2_1-P2_2)/(2.*P_1*P_2); 959 815 G4double CosT13=(P2_2-P2_1-P2_3)/(2.*P_1*P_3); 960 //G4cout<<"CosT12 CosT13 "<<CosT12<<" "<<CosT13<<" "<<std::acos(CosT12)*180./3.14159<<" "<<std::acos(CosT13)*180./3.14159<<G4endl; 816 // Pt=P_2*std::sqrt(1.-CosT12*CosT12); // because system was rotated 11.12.09 817 961 818 if((std::abs(CosT12) >1.) || (std::abs(CosT13) > 1.)) 962 819 { Kink=false;} … … 964 821 { 965 822 Kink=true; 823 Pt=P_2*std::sqrt(1.-CosT12*CosT12); // because system was rotated 11.12.09 966 824 Pstart.setPx(-Pt); Pstart.setPy(0.); Pstart.setPz(P_3*CosT13); 967 825 Pend.setPx( 0.); Pend.setPy( 0.); Pend.setPz( P_1); … … 981 839 PkinkQ2=Pkink - PkinkQ1; 982 840 //------------------------- Minimizing Pt1^2+Pt3^2 ------------------------------ 983 /* 984 G4cout<<"Pstart "<<Pstart<<G4endl; 985 G4cout<<"Pkink "<<Pkink <<G4endl; 986 G4cout<<"Pkink1 "<<PkinkQ1<<G4endl; 987 G4cout<<"Pkink2 "<<PkinkQ2<<G4endl; 988 G4cout<<"Pend "<<Pend <<G4endl; 989 */ 841 990 842 G4double Cos2Psi=(sqr(x1) -sqr(x3)+2.*sqr(x3*CosT13))/ 991 843 std::sqrt(sqr(sqr(x1)-sqr(x3)) + sqr(2.*x1*x3*CosT13)); 992 844 G4double Psi=std::acos(Cos2Psi); 993 //G4cout<<"Psi "<<Psi<<" "<<Psi*180./twopi<<G4endl;994 845 995 846 G4LorentzRotation Rotate; … … 997 848 else {Rotate.rotateY(pi-Psi);} 998 849 Rotate.rotateZ(twopi*G4UniformRand()); 999 1000 //G4cout<<"Rotate"<<G4endl;1001 850 1002 851 Pstart*=Rotate; … … 1005 854 PkinkQ2*=Rotate; 1006 855 Pend*=Rotate; 1007 /* 1008 G4cout<<"Pstart "<<Pstart<<G4endl; 1009 G4cout<<"Pkink1 "<<PkinkQ1 <<G4endl; 1010 G4cout<<"Pkink2 "<<PkinkQ2 <<G4endl; 1011 G4cout<<"Pend "<<Pend <<G4endl; 1012 */ 856 1013 857 } 1014 858 } // end of if((P2_1 < 0.) || (P2_3 <0.)) … … 1018 862 //-------------------------------------------------------------------------------- 1019 863 1020 //G4cout<<"Kink "<<Kink<<G4endl;1021 1022 864 if(Kink) 1023 865 { // Kink is possible … … 1028 870 for(unsigned int Iq=0; Iq <3; Iq++) 1029 871 { 1030 //G4cout<<Iq<<" "<<QuarkProbabilitiesAtGluonSplitUp[Iq]<<G4endl; 872 1031 873 if(Ksi > QuarkProbabilitiesAtGluonSplitUp[Iq]) QuarkInGluon++;} 1032 874 1033 //G4cout<<"Gquark "<<QuarkInGluon<<G4endl;1034 875 G4Parton * Gquark = new G4Parton(QuarkInGluon); 1035 876 G4Parton * Ganti_quark = new G4Parton(-QuarkInGluon); 1036 /*1037 G4cout<<Gquark->GetDefinition()->GetParticleName()<<" "<<Gquark->GetDefinition()->GetPDGEncoding()<<" "<<Gquark->GetDefinition()->GetPDGMass()<<G4endl;1038 G4cout<<Ganti_quark->GetDefinition()->GetParticleName()<<" "<<Ganti_quark->GetDefinition()->GetPDGEncoding()<<" "<<Ganti_quark->GetDefinition()->GetPDGMass()<<G4endl;1039 */1040 1041 //G4int Uzhi; G4cin>>Uzhi;1042 877 1043 878 //------------------------------------------------------------------------------- 1044 /*1045 DefineMomentumInZ(G4double aLightConeMomentum, G4bool aDirection);1046 void Set4Momentum(const G4LorentzVector & aMomentum);1047 G4int PDGencoding;1048 G4ParticleDefinition * theDefinition;1049 G4LorentzVector theMomentum;1050 G4ThreeVector thePosition;1051 1052 G4int theColour;1053 G4double theIsoSpinZ;1054 G4double theSpinZ;1055 1056 G4double theX;1057 */1058 //-------------------------------------------------------------------------------1059 1060 //G4cout<<"Phadron "<<Phadron<<" mass "<<Phadron.mag()<<G4endl;1061 879 G4LorentzRotation toCMS(-1*Phadron.boostVector()); 1062 //G4cout<<"To lab"<<G4endl; 880 1063 881 G4LorentzRotation toLab(toCMS.inverse()); 1064 882 … … 1067 885 PkinkQ2.transform(toLab); 1068 886 Pend.transform(toLab); end->Set4Momentum(Pend); 1069 /* 1070 G4cout<<"Pstart "<<start->GetDefinition()->GetPDGEncoding()<<Pstart<<G4endl; 1071 G4cout<<"Pkink1 "<<PkinkQ1 <<G4endl; 1072 G4cout<<"Pkink2 "<<PkinkQ2 <<G4endl; 1073 G4cout<<"Pend "<<end->GetDefinition()->GetPDGEncoding()<<Pend <<G4endl; 1074 */ 1075 // !!! G4ExcitedString * FirstString(0); G4ExcitedString * SecondString(0); 887 1076 888 G4int absPDGcode=std::abs(hadron->GetDefinition()->GetPDGEncoding()); 1077 1078 //G4cout<<"absPDGcode "<<absPDGcode<<G4endl;1079 889 1080 890 if(absPDGcode < 1000) … … 1088 898 Ganti_quark->Set4Momentum(PkinkQ1); 1089 899 Gquark->Set4Momentum(PkinkQ2); 1090 /* 1091 G4cout<<" Proj Meson end Q"<<G4endl; 1092 G4cout<<"First string ============ "<<G4endl; 1093 G4cout<<"end Q "<<end->GetDefinition()->GetPDGEncoding()<<" "<<end->Get4Momentum()<<G4endl; 1094 G4cout<<"G antiQ"<<Ganti_quark->GetDefinition()->GetPDGEncoding()<<" "<<Ganti_quark->Get4Momentum()<<G4endl; 1095 G4cout<<"Sum P "<<(Ganti_quark->Get4Momentum()+end->Get4Momentum())<<G4endl; 1096 G4cout<<"Secondstring ============ "<<G4endl; 1097 G4cout<<"G Q "<<Gquark->GetDefinition()->GetPDGEncoding()<<" "<<Gquark->Get4Momentum()<<G4endl; 1098 G4cout<<"startQ "<<start->GetDefinition()->GetPDGEncoding()<<" "<<start->Get4Momentum()<<G4endl; 1099 1100 G4cout<<"Sum P "<<(Gquark->Get4Momentum()+start->Get4Momentum())<<G4endl; 1101 */ 900 1102 901 } else 1103 902 { // Anti_Quark on the end … … 1106 905 Gquark->Set4Momentum(PkinkQ1); 1107 906 Ganti_quark->Set4Momentum(PkinkQ2); 1108 /* 1109 G4cout<<" Proj Meson end Qbar"<<G4endl; 1110 G4cout<<"First string ============ "<<G4endl; 1111 G4cout<<"end Q "<<end->GetDefinition()->GetPDGEncoding()<<" "<<end->Get4Momentum()<<G4endl; 1112 G4cout<<"G Q"<<Gquark->GetDefinition()->GetPDGEncoding()<<" "<<Gquark->Get4Momentum()<<G4endl; 1113 G4cout<<"Sum P "<<(Gquark->Get4Momentum()+end->Get4Momentum())<<G4endl; 1114 G4cout<<"Secondstring ============ "<<G4endl; 1115 G4cout<<"G antQ "<<Ganti_quark->GetDefinition()->GetPDGEncoding()<<" "<<Ganti_quark->Get4Momentum()<<G4endl; 1116 G4cout<<"startQ "<<start->GetDefinition()->GetPDGEncoding()<<" "<<start->Get4Momentum()<<G4endl; 1117 G4cout<<"Sum P "<<(Ganti_quark->Get4Momentum()+start->Get4Momentum())<<G4endl; 1118 */ 907 1119 908 } // end of if(end->GetPDGcode() > 0) 1120 909 } else { // Target … … 1125 914 Ganti_quark->Set4Momentum(PkinkQ2); 1126 915 Gquark->Set4Momentum(PkinkQ1); 1127 /* 1128 G4cout<<" Targ Meson end Q"<<G4endl; 1129 G4cout<<"First string ============ "<<G4endl; 1130 G4cout<<"G antiQ"<<Ganti_quark->GetDefinition()->GetPDGEncoding()<<" "<<Ganti_quark->Get4Momentum()<<G4endl; 1131 G4cout<<"end Q "<<end->GetDefinition()->GetPDGEncoding()<<" "<<end->Get4Momentum()<<G4endl; 1132 G4cout<<"Sum P "<<(Ganti_quark->Get4Momentum()+end->Get4Momentum())<<G4endl; 1133 G4cout<<"Secondstring ============ "<<G4endl; 1134 G4cout<<"startQ "<<start->GetDefinition()->GetPDGEncoding()<<" "<<start->Get4Momentum()<<G4endl; 1135 G4cout<<"G Q "<<Gquark->GetDefinition()->GetPDGEncoding()<<" "<<Gquark->Get4Momentum()<<G4endl; 1136 G4cout<<"Sum P "<<(Gquark->Get4Momentum()+start->Get4Momentum())<<G4endl; 1137 */ 916 1138 917 } else 1139 918 { // Anti_Quark on the end … … 1142 921 Gquark->Set4Momentum(PkinkQ2); 1143 922 Ganti_quark->Set4Momentum(PkinkQ1); 1144 /* 1145 G4cout<<" Targ Meson end Qbar"<<G4endl; 1146 G4cout<<"First string ============ "<<G4endl; 1147 G4cout<<"G Q"<<Gquark->GetDefinition()->GetPDGEncoding()<<" "<<Gquark->Get4Momentum()<<G4endl; 1148 G4cout<<"end Q "<<end->GetDefinition()->GetPDGEncoding()<<" "<<end->Get4Momentum()<<G4endl; 1149 G4cout<<"Sum P "<<(Gquark->Get4Momentum()+end->Get4Momentum())<<G4endl; 1150 G4cout<<"Secondstring ============ "<<G4endl; 1151 G4cout<<"startQ "<<start->GetDefinition()->GetPDGEncoding()<<" "<<start->Get4Momentum()<<G4endl; 1152 G4cout<<"G antQ "<<Ganti_quark->GetDefinition()->GetPDGEncoding()<<" "<<Ganti_quark->Get4Momentum()<<G4endl; 1153 G4cout<<"Sum P "<<(Ganti_quark->Get4Momentum()+start->Get4Momentum())<<G4endl; 1154 */ 923 1155 924 } // end of if(end->GetPDGcode() > 0) 1156 925 } // end of if ( isProjectile ) … … 1166 935 Gquark->Set4Momentum(PkinkQ1); 1167 936 Ganti_quark->Set4Momentum(PkinkQ2); 1168 /* 1169 G4cout<<" Proj baryon end QQ"<<G4endl; 1170 G4cout<<"First string ============ "<<G4endl; 1171 G4cout<<"end OQ "<<end->GetDefinition()->GetPDGEncoding()<<" "<<end->Get4Momentum()<<G4endl; 1172 G4cout<<"G Q"<<Gquark->GetDefinition()->GetPDGEncoding()<<" "<<Gquark->Get4Momentum()<<G4endl; 1173 G4cout<<"Sum P "<<(Gquark->Get4Momentum()+end->Get4Momentum())<<G4endl; 1174 G4cout<<"Secondstring ============ "<<G4endl; 1175 G4cout<<"G Qbar"<<Ganti_quark->GetDefinition()->GetPDGEncoding()<<" "<<Ganti_quark->Get4Momentum()<<G4endl; 1176 G4cout<<"startQ "<<start->GetDefinition()->GetPDGEncoding()<<" "<<start->Get4Momentum()<<G4endl; 1177 G4cout<<"Sum P "<<(Ganti_quark->Get4Momentum()+start->Get4Momentum())<<G4endl; 1178 */ 937 1179 938 } else 1180 939 { // Anti_DiQuark on the end or quark … … 1183 942 Ganti_quark->Set4Momentum(PkinkQ1); 1184 943 Gquark->Set4Momentum(PkinkQ2); 1185 /* 1186 G4cout<<" Proj baryon end Q"<<G4endl; 1187 G4cout<<"First string ============ "<<G4endl; 1188 G4cout<<"end OQ "<<end->GetDefinition()->GetPDGEncoding()<<" "<<end->Get4Momentum()<<G4endl; 1189 G4cout<<"G antQ"<<Ganti_quark->GetDefinition()->GetPDGEncoding()<<" "<<Ganti_quark->Get4Momentum()<<G4endl; 1190 G4cout<<"Sum P "<<(Ganti_quark->Get4Momentum()+end->Get4Momentum())<<G4endl; 1191 G4cout<<"Secondstring ============ "<<G4endl; 1192 G4cout<<"G Q "<<Gquark->GetDefinition()->GetPDGEncoding()<<" "<<Gquark->Get4Momentum()<<G4endl; 1193 G4cout<<"startQ "<<start->GetDefinition()->GetPDGEncoding()<<" "<<start->Get4Momentum()<<G4endl; 1194 G4cout<<"Sum P "<<(Gquark->Get4Momentum()+start->Get4Momentum())<<G4endl; 1195 */ 944 1196 945 } // end of if(end->GetPDGcode() > 0) 1197 946 } else { // Target … … 1205 954 Gquark->Set4Momentum(PkinkQ1); 1206 955 Ganti_quark->Set4Momentum(PkinkQ2); 1207 /* 1208 G4cout<<" Targ baryon end QQ"<<G4endl; 1209 G4cout<<"First string ============ "<<G4endl; 1210 G4cout<<"G Q "<<Gquark->GetDefinition()->GetPDGEncoding()<<" "<<Gquark->Get4Momentum()<<G4endl; 1211 G4cout<<"end OQ "<<end->GetDefinition()->GetPDGEncoding()<<" "<<end->Get4Momentum()<<G4endl; 1212 G4cout<<"Sum P "<<(Gquark->Get4Momentum()+end->Get4Momentum())<<G4endl; 1213 G4cout<<"Secondstring ============ "<<G4endl; 1214 G4cout<<"startQ "<<start->GetDefinition()->GetPDGEncoding()<<" "<<start->Get4Momentum()<<G4endl; 1215 G4cout<<"G Qbar"<<Ganti_quark->GetDefinition()->GetPDGEncoding()<<" "<<Ganti_quark->Get4Momentum()<<G4endl; 1216 G4cout<<"Sum P "<<(Ganti_quark->Get4Momentum()+start->Get4Momentum())<<G4endl; 1217 */ 956 1218 957 } else 1219 958 { // Anti_DiQuark on the end or Q … … 1222 961 Gquark->Set4Momentum(PkinkQ2); 1223 962 Ganti_quark->Set4Momentum(PkinkQ1); 1224 /* 1225 G4cout<<" Targ baryon end Q"<<G4endl; 1226 G4cout<<"First string ============ "<<G4endl; 1227 G4cout<<"G Qbar"<<Ganti_quark->GetDefinition()->GetPDGEncoding()<<" "<<Ganti_quark->Get4Momentum()<<G4endl; 1228 G4cout<<"end O "<<end->GetDefinition()->GetPDGEncoding()<<" "<<end->Get4Momentum()<<G4endl; 1229 G4cout<<"Sum P "<<(Ganti_quark->Get4Momentum()+end->Get4Momentum())<<G4endl; 1230 G4cout<<"Secondstring ============ "<<G4endl; 1231 G4cout<<"startQ "<<start->GetDefinition()->GetPDGEncoding()<<" "<<start->Get4Momentum()<<G4endl; 1232 G4cout<<"G Q "<<Gquark->GetDefinition()->GetPDGEncoding()<<" "<<Gquark->Get4Momentum()<<G4endl; 1233 G4cout<<"Sum P "<<(Gquark->Get4Momentum()+start->Get4Momentum())<<G4endl; 1234 */ 963 1235 964 } // end of if(end->GetPDGcode() > 0) 1236 965 } // end of if ( isProjectile ) … … 1346 1075 { // @@ this method is used in FTFModel as well. Should go somewhere common! 1347 1076 1348 G4double Pt2; 1349 Pt2 = -AveragePt2 * std::log(1. + G4UniformRand() * 1350 (std::exp(-maxPtSquare/AveragePt2)-1.)); 1351 1077 G4double Pt2(0.); 1078 if(AveragePt2 <= 0.) {Pt2=0.;} 1079 else 1080 { 1081 Pt2 = -AveragePt2 * std::log(1. + G4UniformRand() * 1082 (std::exp(-maxPtSquare/AveragePt2)-1.)); 1083 } 1352 1084 G4double Pt=std::sqrt(Pt2); 1353 1085 G4double phi=G4UniformRand() * twopi; -
trunk/source/processes/hadronic/models/parton_string/diffraction/src/G4DiffractiveSplitableHadron.cc
r1196 r1228 26 26 // 27 27 // $Id: G4DiffractiveSplitableHadron.cc,v 1.8 2009/07/31 11:03:00 vuzhinsk Exp $ 28 // GEANT4 tag $Name: geant4-09-03 -cand-01$28 // GEANT4 tag $Name: geant4-09-03 $ 29 29 // 30 30 -
trunk/source/processes/hadronic/models/parton_string/diffraction/src/G4ElasticHNScattering.cc
r1196 r1228 25 25 // 26 26 // 27 // $Id: G4ElasticHNScattering.cc,v 1.1 1 2009/10/06 10:10:36 vuzhinskExp $27 // $Id: G4ElasticHNScattering.cc,v 1.14 2009/12/16 17:51:13 gunter Exp $ 28 28 // ------------------------------------------------------------ 29 29 // GEANT 4 class implemetation file … … 58 58 G4FTFParameters *theParameters) const 59 59 { 60 //G4cout<<"G4ElasticHNScattering::ElasticScattering"<<G4endl;61 60 // -------------------- Projectile parameters ----------------------------------- 62 61 G4LorentzVector Pprojectile=projectile->Get4Momentum(); 63 //G4cout<<"Elastic P "<<Pprojectile<<G4endl; 62 64 63 if(Pprojectile.z() < 0.) 65 64 { … … 68 67 } 69 68 70 // G4double ProjectileRapidity = Pprojectile.rapidity();71 // G4int ProjectilePDGcode=projectile->GetDefinition()->GetPDGEncoding();72 // G4ParticleDefinition * ProjectileDefinition = projectile->GetDefinition();73 74 69 G4bool PutOnMassShell(false); 75 70 … … 86 81 87 82 // -------------------- Target parameters ---------------------------------------------- 88 // G4int TargetPDGcode=target->GetDefinition()->GetPDGEncoding();89 83 90 84 G4LorentzVector Ptarget=target->Get4Momentum(); 91 //G4cout<<"Elastic T "<<Ptarget<<G4endl; 92 // G4double TargetRapidity = Ptarget.rapidity(); 85 93 86 G4double M0target = Ptarget.mag(); 94 87 … … 149 142 else // if(M0projectile > projectile->GetDefinition()->GetPDGMass()) 150 143 { 151 target->SetStatus(2); // Uzhi 17.07.09144 target->SetStatus(2); 152 145 return false; // The projectile was not excited, 153 146 // but the energy was too low to put … … 183 176 G4double maxPtSquare = PZcms2; 184 177 185 //----- Charge exchange between the projectile and the target, if possible186 // (G4UniformRand() < 0.5*std::exp(-0.5*(ProjectileRapidity - TargetRapidity))))187 /*188 if((ProjectilePDGcode != TargetPDGcode) &&189 ((ProjectilePDGcode > 1000) && (TargetPDGcode > 1000)) &&190 (G4UniformRand() < 1.0*std::exp(-0.5*(ProjectileRapidity - TargetRapidity))))191 {192 projectile->SetDefinition(target->GetDefinition());193 target->SetDefinition(ProjectileDefinition);194 }195 */196 178 // ------ Now we can calculate the transfered Pt -------------------------- 197 179 G4double Pt2; … … 255 237 { // @@ this method is used in FTFModel as well. Should go somewhere common! 256 238 257 G4double Pt2; 258 Pt2 = -AveragePt2 * std::log(1. + G4UniformRand() * 239 G4double Pt2(0.); 240 if(AveragePt2 <= 0.) {Pt2=0.;} 241 else 242 { 243 Pt2 = -AveragePt2 * std::log(1. + G4UniformRand() * 259 244 (std::exp(-maxPtSquare/AveragePt2)-1.)); 260 245 } 261 246 G4double Pt=std::sqrt(Pt2); 262 247 G4double phi=G4UniformRand() * twopi; -
trunk/source/processes/hadronic/models/parton_string/diffraction/src/G4FTFModel.cc
r1196 r1228 25 25 // 26 26 // 27 // $Id: G4FTFModel.cc,v 1. 28 2009/10/29 14:55:33vuzhinsk Exp $28 // GEANT4 tag $Name: geant4-09-03 -cand-01$27 // $Id: G4FTFModel.cc,v 1.34 2009/12/15 19:14:31 vuzhinsk Exp $ 28 // GEANT4 tag $Name: geant4-09-03 $ 29 29 // 30 30 … … 56 56 G4VPartonStringModel::SetThisPointer(this); 57 57 theParameters=0; 58 NumberOfInvolvedNucleon=0; 58 59 } 59 60 … … 61 62 G4FTFModel::~G4FTFModel() 62 63 { 63 if( theParameters != 0 ) delete theParameters;64 64 // Because FTF model can be called for various particles 65 65 // theParameters must be erased at the end of each call. 66 66 // Thus the delete is also in G4FTFModel::GetStrings() method 67 if( theParameters != 0 ) delete theParameters; 67 68 if( theExcitation != 0 ) delete theExcitation; 68 69 if( theElastic != 0 ) delete theElastic; 70 71 if( NumberOfInvolvedNucleon != 0) 72 { 73 for(G4int i=0; i < NumberOfInvolvedNucleon; i++) 74 { 75 G4VSplitableHadron * aNucleon = TheInvolvedNucleon[i]->GetSplitableHadron(); 76 if(aNucleon) delete aNucleon; 77 } 78 } 69 79 } 70 80 … … 93 103 94 104 // --- cms energy 95 96 105 G4double s = sqr( theProjectile.GetMass() ) + 97 106 sqr( G4Proton::Proton()->GetPDGMass() ) + 98 107 2*theProjectile.GetTotalEnergy()*G4Proton::Proton()->GetPDGMass(); 99 /*100 G4cout<<"----------------------------------------"<<G4endl;101 G4cout << " primary Total E (GeV): " << theProjectile.GetTotalEnergy()/GeV << G4endl;102 G4cout << " primary Mass (GeV): " << theProjectile.GetMass() /GeV << G4endl;103 G4cout << " primary 3Mom " << theProjectile.GetMomentum() << G4endl;104 G4cout << " primary space position " << theProjectile.GetPositionInNucleus() << G4endl;105 G4cout << "cms std::sqrt(s) (GeV) = " << std::sqrt(s) / GeV << G4endl;106 */107 108 108 109 if( theParameters != 0 ) delete theParameters; … … 110 111 aNucleus.GetN(),aNucleus.GetZ(), 111 112 s); 113 112 114 //theParameters->SetProbabilityOfElasticScatt(0.); 113 115 // To turn on/off (1/0) elastic scattering … … 116 118 117 119 // ------------------------------------------------------------ 120 struct DeleteVSplitableHadron { void operator()(G4VSplitableHadron * aH){ delete aH;} }; 121 122 123 // ------------------------------------------------------------ 118 124 G4ExcitedStringVector * G4FTFModel::GetStrings() 119 { 120 //G4int Uzhi; G4cin>>Uzhi; 121 122 //G4cout<<"GetList"<<G4endl; 125 { 126 G4ExcitedStringVector * theStrings(0); 127 123 128 theParticipants.GetList(theProjectile,theParameters); 124 //G4cout<<"Regge"<<G4endl; 125 ReggeonCascade(); // Uzhi 26 July 09 126 //G4cout<<"On mass shell"<<G4endl; 127 if (! PutOnMassShell() ) return NULL; // Uzhi 26 July 09 128 //G4cout<<"Excite"<<G4endl; 129 if (! ExciteParticipants()) return NULL; 130 //G4cout<<"Strings"<<G4endl; 131 G4ExcitedStringVector * theStrings = BuildStrings(); 132 //G4cout<<"Out FTF N strings "<<theStrings->size()<<G4endl; 133 //G4cout<<"GetResidualNucleus"<<G4endl; 134 GetResidualNucleus(); 135 //G4cout<<"Out of FTF"<<G4endl; 136 if( theParameters != 0 ) 137 { 138 delete theParameters; 139 theParameters=0; 129 130 ReggeonCascade(); 131 132 G4bool Success(true); 133 if( PutOnMassShell() ) 134 { 135 if( ExciteParticipants() ) 136 { 137 theStrings = BuildStrings(); 138 139 GetResidualNucleus(); 140 141 if( theParameters != 0 ) 142 { 143 delete theParameters; 144 theParameters=0; 145 } 146 } else // if( ExciteParticipants() ) 147 { Success=false;} 148 } else // if( PutOnMassShell() ) 149 { Success=false;} 150 151 if(!Success) 152 { 153 // -------------- Erase the projectile ---------------- 154 std::vector<G4VSplitableHadron *> primaries; 155 156 theParticipants.StartLoop(); // restart a loop 157 while ( theParticipants.Next() ) 158 { 159 const G4InteractionContent & interaction=theParticipants.GetInteraction(); 160 // do not allow for duplicates ... 161 if ( primaries.end() == std::find(primaries.begin(), primaries.end(), 162 interaction.GetProjectile()) ) 163 primaries.push_back(interaction.GetProjectile()); 164 } 165 std::for_each(primaries.begin(), primaries.end(), DeleteVSplitableHadron()); 166 primaries.clear(); 140 167 } 141 return theStrings; 142 } 143 144 // ------------------------------------------------------------ 145 struct DeleteVSplitableHadron { void operator()(G4VSplitableHadron * aH){ delete aH;} }; 146 147 // ------------------------------------------------------------ 148 void G4FTFModel::ReggeonCascade() // Uzhi 26 July 2009 168 169 // -------------- Cleaning of the memory -------------- 170 // -------------- Erase the target nucleons ----------- 171 G4VSplitableHadron * aNucleon = 0; 172 for(G4int i=0; i < NumberOfInvolvedNucleon; i++) 173 { 174 aNucleon = TheInvolvedNucleon[i]->GetSplitableHadron(); 175 if(aNucleon) delete aNucleon; 176 } 177 178 NumberOfInvolvedNucleon=0; 179 180 return theStrings; 181 182 } 183 //------------------------------------------------------------------- 184 void G4FTFModel::ReggeonCascade() 149 185 { //--- Implementation of reggeon theory inspired model------- 150 186 NumberOfInvolvedNucleon=0; 151 187 152 //G4int PrimInt(0);153 188 theParticipants.StartLoop(); 154 189 while (theParticipants.Next()) 155 { 156 //PrimInt++; 190 { 157 191 const G4InteractionContent & collision=theParticipants.GetInteraction(); 158 192 G4Nucleon * TargetNucleon=collision.GetTargetNucleon(); 159 //G4cout<<"Prim Nnucl "<<TargetNucleon->Get4Momentum()<<G4endl; 193 160 194 TheInvolvedNucleon[NumberOfInvolvedNucleon]=TargetNucleon; 161 195 NumberOfInvolvedNucleon++; … … 165 199 166 200 theParticipants.theNucleus->StartLoop(); 167 G4Nucleon * Neighbour ;168 //G4int NucleonNum(0); 201 G4Nucleon * Neighbour(0); 202 169 203 while ( (Neighbour = theParticipants.theNucleus->GetNextNucleon()) ) 170 204 { … … 177 211 std::exp(-impact2/theParameters->GetR2ofNuclearDestruction())) 178 212 { // The neighbour nucleon is involved in the reggeon cascade 179 //G4cout<<" involved "<<NucleonNum<<" "<<Neighbour->Get4Momentum()<<G4endl; 213 180 214 TheInvolvedNucleon[NumberOfInvolvedNucleon]=Neighbour; 181 215 NumberOfInvolvedNucleon++; 182 //G4cout<<" PrimInt"<<" "<<NumberOfInvolvedNucleon<<G4endl;183 216 184 217 G4VSplitableHadron *targetSplitable; … … 186 219 187 220 Neighbour->Hit(targetSplitable); 188 // Neighbour->SetBindingEnergy(3.*Neighbour->GetBindingEnergy()); // Uzhi 5.10.09189 221 targetSplitable->SetStatus(2); 190 222 } 191 223 } // end of if(!Neighbour->AreYouHit()) 192 //NucleonNum++;193 224 } // end of while (theParticipant.theNucleus->GetNextNucleon()) 194 //G4cout<<"Prim Int N Ninvolv "<<PrimInt<<" "<<NumberOfInvolvedNucleon<<G4endl;195 225 } // end of while (theParticipants.Next()) 196 //G4cout<<"At end "<<PrimInt<<" "<<NumberOfInvolvedNucleon<<G4endl;197 226 198 227 // ---------------- Calculation of creation time for each target nucleon ----------- … … 240 269 { 241 270 // -------------- Properties of the projectile ---------------- 242 //G4cout<<"Put on Mass-shell"<<G4endl;243 271 theParticipants.StartLoop(); // restart a loop 244 272 theParticipants.Next(); 245 273 G4VSplitableHadron * primary = theParticipants.GetInteraction().GetProjectile(); 246 274 G4LorentzVector Pprojectile=primary->Get4Momentum(); 247 //G4cout<<"Proj "<<Pprojectile<<G4endl; 275 276 // To get original projectile particle 277 248 278 if(Pprojectile.z() < 0.){return false;} 249 279 … … 252 282 //------------------------------------------------------------- 253 283 G4LorentzVector Psum = Pprojectile; 254 G4double SumMasses = Mprojectile; 284 G4double SumMasses = Mprojectile + 20.*MeV; // 13.12.09 285 // Separation energy for projectile 255 286 256 287 //--------------- Target nucleus ------------------------------ … … 259 290 G4int ResidualMassNumber=theNucleus->GetMassNumber(); 260 291 G4int ResidualCharge =theNucleus->GetCharge(); 261 //G4cout<<"Nucleus "<<ResidualMassNumber<<" "<<ResidualCharge<<G4endl;262 292 ResidualExcitationEnergy=0.; 263 293 G4LorentzVector PnuclearResidual(0.,0.,0.,0.); … … 267 297 268 298 theNucleus->StartLoop(); 269 //G4int NucleonNum(0); 299 270 300 while ((aNucleon = theNucleus->GetNextNucleon())) 271 301 { 272 302 if(aNucleon->AreYouHit()) 273 303 { // Involved nucleons 274 //G4cout<<" "<<NucleonNum<<" "<<aNucleon->Get4Momentum()<<G4endl;275 304 Psum += aNucleon->Get4Momentum(); 276 305 SumMasses += aNucleon->GetDefinition()->GetPDGMass(); 306 SumMasses += 20.*MeV; // 13.12.09 Separation energy for a nucleon 277 307 ResidualMassNumber--; 278 308 ResidualCharge-=(G4int) aNucleon->GetDefinition()->GetPDGCharge(); … … 283 313 PnuclearResidual += aNucleon->Get4Momentum(); 284 314 } // end of if(!aNucleon->AreYouHit()) 285 //NucleonNum++;286 315 } // end of while (theNucleus->GetNextNucleon()) 287 316 288 //G4cout<<"Nucleus "<<ResidualMassNumber<<" "<<ResidualCharge<<G4endl;289 //G4cout<<"PResid "<<PnuclearResidual<<G4endl;290 291 317 Psum += PnuclearResidual; 318 292 319 G4double ResidualMass(0.); 293 320 if(ResidualMassNumber == 0) … … 303 330 } 304 331 305 // 332 // ResidualMass +=ResidualExcitationEnergy; // Will be given after checks 306 333 SumMasses += ResidualMass; 307 334 … … 310 337 G4double SqrtS=Psum.mag(); 311 338 G4double S=Psum.mag2(); 312 //G4cout<<" SqrtS SumMasses "<<SqrtS <<" "<< SumMasses<<G4endl; 313 //G4cout<<"Res M E* "<<ResidualMass<<" "<<ResidualExcitationEnergy<<G4endl; 339 314 340 if(SqrtS < SumMasses) {return false;} // It is impossible to simulate 315 341 // after putting nuclear nucleons 316 342 // on mass-shell 343 317 344 if(SqrtS < SumMasses+ResidualExcitationEnergy) {ResidualExcitationEnergy=0.;} 318 345 319 346 ResidualMass +=ResidualExcitationEnergy; 320 347 SumMasses +=ResidualExcitationEnergy; 321 //G4cout<<"Res M E* "<<ResidualMass<<" "<<ResidualExcitationEnergy<<G4endl;322 348 323 349 //------------------------------------------------------------- … … 325 351 G4int MaxNumberOfDeltas = (int)((SqrtS - SumMasses)/(400.*MeV)); 326 352 G4int NumberOfDeltas(0); 327 //SumMasses=Mprojectile; 353 328 354 if(theNucleus->GetMassNumber() != 1) 329 355 { 330 G4double ProbDeltaIsobar(0.); // 1. *** ****************************356 G4double ProbDeltaIsobar(0.); // 1. *** Can be set if it is needed 331 357 for(G4int i=0; i < NumberOfInvolvedNucleon; i++ ) 332 358 { … … 343 369 targetSplitable->SetDefinition(ptr); 344 370 SumMasses+=targetSplitable->GetDefinition()->GetPDGMass(); 345 //G4cout<<" i "<<i<<" Delta "<<targetSplitable->GetDefinition()->GetPDGMass()<<G4endl;346 371 } 347 else348 {349 // SumMasses+=TheInvolvedNucleon[i]->GetSplitableHadron()->350 // GetDefinition()->GetPDGMass();351 //G4cout<<" i "<<i<<" Nuclo "<<TheInvolvedNucleon[i]->GetSplitableHadron()->GetDefinition()->GetPDGMass()<<G4endl;352 }353 372 } // end of for(G4int i=0; i < NumberOfInvolvedNucleon; i++ ) 354 373 } // end of if(theNucleus.GetMassNumber() != 1) 355 //G4cout<<"MaxNumberOfDeltas NumberOfDeltas "<<MaxNumberOfDeltas<<" "<<NumberOfDeltas<<G4endl; 356 //G4cout<<" SqrtS SumMasses "<<SqrtS <<" "<< SumMasses<<G4endl; 357 //------------------------------------------------------------- 358 374 //------------------------------------------------------------- 359 375 G4LorentzRotation toCms(-1*Psum.boostVector()); 360 376 G4LorentzVector Ptmp=toCms*Pprojectile; 361 362 377 if ( Ptmp.pz() <= 0. ) 363 378 { // "String" moving backwards in CMS, abort collision !! … … 366 381 } 367 382 368 toCms.rotateZ(-1*Ptmp.phi()); 369 toCms.rotateY(-1*Ptmp.theta()); 383 // toCms.rotateZ(-1*Ptmp.phi()); // Uzhi 5.12.09 384 // toCms.rotateY(-1*Ptmp.theta()); // Uzhi 5.12.09 370 385 371 386 G4LorentzRotation toLab(toCms.inverse()); 372 387 373 // Pprojectile.transform(toCms);374 //G4cout<<"Proj in CMS "<<Pprojectile<<G4endl;375 376 //G4cout<<"Main work"<<G4endl;377 388 //------------------------------------------------------------- 378 389 //------- Ascribing of the involved nucleons Pt and Xminus ---- 379 390 G4double Dcor = theParameters->GetDofNuclearDestruction()/ 380 391 theNucleus->GetMassNumber(); 381 // NumberOfInvolvedNucleon; 392 382 393 G4double AveragePt2 = theParameters->GetPt2ofNuclearDestruction(); 383 394 G4double maxPtSquare = theParameters->GetMaxPt2ofNuclearDestruction(); 384 //G4cout<<" D Pt2 "<<Dcor<<" "<<AveragePt2<<G4endl;385 395 386 396 G4double M2target(0.); 397 G4double WminusTarget(0.); 398 G4double WplusProjectile(0.); 399 387 400 G4int NumberOfTries(0); 388 401 G4double ScaleFactor(1.); 389 do // while (SqrtS < Mprojectile + std::sqrt(M2target)) 390 { // while (DecayMomentum < 0.) 391 392 NumberOfTries++; 393 //G4cout<<"NumberOfTries "<<NumberOfTries<<G4endl; 394 if(NumberOfTries == 100*(NumberOfTries/100)) 395 { // At large number of tries it would be better to reduce the values 396 ScaleFactor/=2.; 397 Dcor *=ScaleFactor; 398 AveragePt2 *=ScaleFactor; 399 //G4cout<<"NumberOfTries "<<NumberOfTries<<" "<<Dcor<<" "<<AveragePt2<<G4endl; 400 //G4int Uzhi; G4cin>>Uzhi; 401 } 402 //G4cout<<"Start Decay "<<G4endl; G4int Uzhi; G4cin>>Uzhi; 403 G4ThreeVector PtSum(0.,0.,0.); 404 G4double XminusSum(0.); 405 G4double Xminus(0.); 406 G4bool Success=true; 407 408 do // while(Success == false); 409 { 410 //G4cout<<"Sample Pt and X"<<" Ninv "<<NumberOfInvolvedNucleon<<G4endl; // G4int Uzhi1; G4cin>>Uzhi1; 411 Success=true; 402 G4bool OuterSuccess(true); 403 do // while (!OuterSuccess) 404 { 405 OuterSuccess=true; 406 407 do // while (SqrtS < Mprojectile + std::sqrt(M2target)) 408 { // while (DecayMomentum < 0.) 409 410 NumberOfTries++; 411 if(NumberOfTries == 100*(NumberOfTries/100)) // 100 412 { // At large number of tries it would be better to reduce the values 413 ScaleFactor/=2.; 414 Dcor *=ScaleFactor; 415 AveragePt2 *=ScaleFactor; 416 } 417 418 G4ThreeVector PtSum(0.,0.,0.); 419 G4double XminusSum(0.); 420 G4double Xminus(0.); 421 G4bool InerSuccess=true; 422 423 do // while(!InerSuccess); 424 { 425 InerSuccess=true; 412 426 413 427 PtSum =G4ThreeVector(0.,0.,0.); … … 420 434 G4ThreeVector tmpPt = GaussianPt(AveragePt2, maxPtSquare); 421 435 PtSum += tmpPt; 422 423 436 G4ThreeVector tmpX=GaussianPt(Dcor*Dcor, 1.); 424 437 Xminus=tmpX.x(); … … 427 440 G4LorentzVector tmp(tmpPt.x(),tmpPt.y(),Xminus,0.); 428 441 aNucleon->SetMomentum(tmp); 429 //G4cout<<"Nucl mom "<<aNucleon->GetMomentum()<<G4endl;430 442 } // end of for(G4int i=0; i < NumberOfInvolvedNucleon; i++ ) 431 443 … … 445 457 DeltaXminus = -1./theNucleus->GetMassNumber(); 446 458 } 447 //G4cout<<"Deltas "<<DeltaX<<" "<<DeltaY<<" "<<DeltaXminus<<G4endl;448 459 449 460 XminusSum=1.; 450 461 M2target =0.; 451 452 462 453 463 for(G4int i=0; i < NumberOfInvolvedNucleon; i++ ) … … 456 466 457 467 Xminus = aNucleon->Get4Momentum().pz() - DeltaXminus; 458 //G4cout<<i<<" "<<Xminus<<" "<<XminusSum<<G4endl;459 468 XminusSum-=Xminus; 469 460 470 if((Xminus <= 0.) || (Xminus > 1.) || 461 (XminusSum <=0.) || (XminusSum > 1.)) { Success=false; break;}471 (XminusSum <=0.) || (XminusSum > 1.)) {InerSuccess=false; break;} 462 472 463 473 G4double Px=aNucleon->Get4Momentum().px() - DeltaX; … … 472 482 } // end of for(G4int i=0; i < NumberOfInvolvedNucleon; i++ ) 473 483 474 if( Success && (ResidualMassNumber != 0))484 if(InerSuccess && (ResidualMassNumber != 0)) 475 485 { 476 486 M2target +=(ResidualMass*ResidualMass + PtSum.mag2())/XminusSum; 477 487 } 478 //G4cout<<"Success "<<Success<<G4endl; 479 //G4int Uzhi; G4cin>>Uzhi; 480 } while(Success == 0); 481 482 //------------------------------------------------------------- 483 //G4cout<<"SqrtS > Mprojectile + std::sqrt(M2target) "<<SqrtS<<" "<<Mprojectile<<" "<< std::sqrt(M2target)<<" "<<Mprojectile + std::sqrt(M2target)<<G4endl; 484 //G4int Uzhi3; G4cin>>Uzhi3; 485 } while (SqrtS < Mprojectile + std::sqrt(M2target)); 486 487 //------------------------------------------------------------- 488 G4double DecayMomentum2= S*S+M2projectile*M2projectile+M2target*M2target 489 -2.*S*M2projectile - 2.*S*M2target 490 -2.*M2projectile*M2target; 491 492 G4double WminusTarget=(S-M2projectile+M2target+std::sqrt(DecayMomentum2))/2./SqrtS; 493 G4double WplusProjectile=SqrtS - M2target/WminusTarget; 494 //G4cout<<"DecM W+ W- "<<DecayMomentum2<<" "<<WplusProjectile<<" "<<WminusTarget<<G4endl; 495 //------------------------------------------------------------- 496 G4double Pzprojectile=WplusProjectile/2. - M2projectile/2./WplusProjectile; 497 G4double Eprojectile =WplusProjectile/2. + M2projectile/2./WplusProjectile; 498 Pprojectile.setPz(Pzprojectile); Pprojectile.setE(Eprojectile); 499 500 Pprojectile.transform(toLab); // The work with the projectile 501 primary->Set4Momentum(Pprojectile); // is finished at the moment. 502 //G4cout<<"Proj After "<<Pprojectile<<G4endl; 503 //------------------------------------------------------------- 504 //Ninvolv=0; 505 506 G4ThreeVector Residual3Momentum(0.,0.,1.); 507 508 for(G4int i=0; i < NumberOfInvolvedNucleon; i++ ) 509 { 488 } while(!InerSuccess); 489 } while (SqrtS < Mprojectile + std::sqrt(M2target)); 490 //------------------------------------------------------------- 491 G4double DecayMomentum2= S*S+M2projectile*M2projectile+M2target*M2target 492 -2.*S*M2projectile - 2.*S*M2target 493 -2.*M2projectile*M2target; 494 495 WminusTarget=(S-M2projectile+M2target+std::sqrt(DecayMomentum2))/2./SqrtS; 496 WplusProjectile=SqrtS - M2target/WminusTarget; 497 //------------------------------------------------------------- 498 for(G4int i=0; i < NumberOfInvolvedNucleon; i++ ) 499 { 510 500 G4Nucleon * aNucleon = TheInvolvedNucleon[i]; 511 501 G4LorentzVector tmp=aNucleon->Get4Momentum(); 512 Residual3Momentum-=tmp.vect();513 502 514 503 G4double Mt2 = sqr(tmp.x())+sqr(tmp.y())+ … … 520 509 G4double E = WminusTarget*Xminus/2. + Mt2/(2.*WminusTarget*Xminus); 521 510 511 if( E+Pz > WplusProjectile ){OuterSuccess=false; break;} 512 } // end of for(G4int i=0; i < NumberOfInvolvedNucleon; i++ ) 513 } while(!OuterSuccess); 514 515 //------------------------------------------------------------- 516 G4double Pzprojectile=WplusProjectile/2. - M2projectile/2./WplusProjectile; 517 G4double Eprojectile =WplusProjectile/2. + M2projectile/2./WplusProjectile; 518 Pprojectile.setPz(Pzprojectile); Pprojectile.setE(Eprojectile); 519 520 Pprojectile.transform(toLab); // The work with the projectile 521 primary->Set4Momentum(Pprojectile); // is finished at the moment. 522 523 //------------------------------------------------------------- 524 G4ThreeVector Residual3Momentum(0.,0.,1.); 525 526 for(G4int i=0; i < NumberOfInvolvedNucleon; i++ ) 527 { 528 G4Nucleon * aNucleon = TheInvolvedNucleon[i]; 529 G4LorentzVector tmp=aNucleon->Get4Momentum(); 530 Residual3Momentum-=tmp.vect(); 531 532 G4double Mt2 = sqr(tmp.x())+sqr(tmp.y())+ 533 aNucleon->GetSplitableHadron()->GetDefinition()->GetPDGMass()* 534 aNucleon->GetSplitableHadron()->GetDefinition()->GetPDGMass(); 535 G4double Xminus=tmp.z(); 536 537 G4double Pz=-WminusTarget*Xminus/2. + Mt2/(2.*WminusTarget*Xminus); 538 G4double E = WminusTarget*Xminus/2. + Mt2/(2.*WminusTarget*Xminus); 539 522 540 tmp.setPz(Pz); 523 541 tmp.setE(E); 542 524 543 tmp.transform(toLab); 525 //G4cout<<"invol "<<Ninvolv<<" "<<tmp<<G4endl; 526 //Ninvolv++; 544 527 545 aNucleon->SetMomentum(tmp); 546 528 547 G4VSplitableHadron * targetSplitable=aNucleon->GetSplitableHadron(); 529 548 targetSplitable->Set4Momentum(tmp); 530 //G4cout<<"nucleon M "<<aNucleon->Get4Momentum()<<G4endl;531 549 532 550 } // end of for(G4int i=0; i < NumberOfInvolvedNucleon; i++ ) 533 551 534 552 G4double Mt2Residual=sqr(ResidualMass) + 535 sqr(Residual3Momentum.x())+sqr(Residual3Momentum. x());553 sqr(Residual3Momentum.x())+sqr(Residual3Momentum.y()); 536 554 537 555 G4double PzResidual=-WminusTarget*Residual3Momentum.z()/2. + … … 540 558 Mt2Residual/(2.*WminusTarget*Residual3Momentum.z()); 541 559 542 // G4LorentzVector Residual4Momentum(0.,0.,0.,0.);543 560 Residual4Momentum.setPx(Residual3Momentum.x()); 544 561 Residual4Momentum.setPy(Residual3Momentum.y()); 545 562 Residual4Momentum.setPz(PzResidual); 546 563 Residual4Momentum.setE(EResidual); 564 547 565 Residual4Momentum.transform(toLab); 548 549 //G4cout<<"Return Nucleus"<<G4endl; 550 //------------------------------------------------------------- 551 //------------------------------------------------------------- 552 //------------------------------------------------------------- 553 //G4int Uzhi2; G4cin>>Uzhi2; 554 566 //------------------------------------------------------------- 555 567 return true; 556 568 } … … 559 571 G4bool G4FTFModel::ExciteParticipants() 560 572 { 561 // // Uzhi 29.03.08 For elastic Scatt. 562 //G4cout<<" In ExciteParticipants() "<<theParticipants.theInteractions.size()<<G4endl; 563 //G4cout<<" test Params Tot "<<theParameters->GetTotalCrossSection()<<G4endl; 564 //G4cout<<" test Params Ela "<<theParameters->GetElasticCrossSection()<<G4endl; 565 566 //G4int counter=0; 567 // // Uzhi 29.03.08 568 569 570 //G4int InterNumber=0; // Vova 571 572 G4bool Successfull=false; 573 theParticipants.StartLoop(); //Uzhi 26 July 09 574 575 // while (theParticipants.Next()&& (InterNumber < 3)) // Vova 573 G4bool Successfull(false); 574 // do { // } while (Successfull == false) // Closed 15.12.09 575 Successfull=false; 576 theParticipants.StartLoop(); 576 577 while (theParticipants.Next()) 577 578 { 578 579 const G4InteractionContent & collision=theParticipants.GetInteraction(); 579 // 580 //counter++; 581 //G4cout<<" Int num "<<counter<<G4endl; 582 // 580 583 581 G4VSplitableHadron * projectile=collision.GetProjectile(); 584 582 G4VSplitableHadron * target=collision.GetTarget(); 585 // G4Nucleon * TargetNucleon=collision.GetTargetNucleon(); // Uzhi 16.07.09 586 // Uzhi 16.07.09 ---------------------------- 583 587 584 if(G4UniformRand()< theParameters->GetProbabilityOfElasticScatt()) 588 585 { // Elastic scattering ------------------------- 589 //G4cout<<"Elastic"<<G4endl;590 586 if(theElastic->ElasticScattering(projectile, target, theParameters)) 591 587 { 592 588 Successfull = Successfull || true; 593 589 } else 594 590 { 595 //G4cout<<"Elastic Not succes"<<G4endl;596 591 Successfull = Successfull || false; 597 592 target->SetStatus(2); 598 /* 599 if(target->GetStatus() == 0) // Uzhi 17.07.09 600 { 601 G4VSplitableHadron * aHit=0; // Uzhi 16.07.09 602 TargetNucleon->Hit(aHit); // Uzhi 16.07.09 603 }; 604 */ 605 }; 593 } 606 594 } 607 595 else 608 596 { // Inelastic scattering ---------------------- 609 //G4cout<<"InElastic"<<G4endl;610 597 if(theExcitation->ExciteParticipants(projectile, target, 611 598 theParameters, theElastic)) … … 614 601 } else 615 602 { 616 //G4cout<<"InElastic Non succes"<<G4endl;617 603 Successfull = Successfull || false; 618 604 target->SetStatus(2); 619 /* 620 if(target->GetStatus() == 0) // Uzhi 16.06.09 621 { 622 G4VSplitableHadron * aHit=0; // Uzhi 16.07.09 623 TargetNucleon->Hit(aHit); // Uzhi 16.07.09 624 }; 625 */ 626 }; 605 } 627 606 } 628 } // end of the loop Uzhi 9.07.09 629 // Uzhi 16.07.09 ---------------------------- 630 631 if(!Successfull) 632 { 633 //G4cout<<"Process not successfull"<<G4endl; 634 // give up, clean up 635 std::vector<G4VSplitableHadron *> primaries; 636 // std::vector<G4VSplitableHadron *> targets; // Uzhi 31.07.09 637 theParticipants.StartLoop(); // restart a loop 638 while ( theParticipants.Next() ) 639 { 640 const G4InteractionContent & interaction=theParticipants.GetInteraction(); 641 // do not allow for duplicates ... 642 if ( primaries.end() == std::find(primaries.begin(), primaries.end(), 643 interaction.GetProjectile()) ) 644 primaries.push_back(interaction.GetProjectile()); 645 /* // Uzhi 31.07.09 646 if ( targets.end() == std::find(targets.begin(), targets.end(), 647 interaction.GetTarget()) ) 648 targets.push_back(interaction.GetTarget()); 649 */ // Uzhi 31.07.09 650 } 651 std::for_each(primaries.begin(), primaries.end(), DeleteVSplitableHadron()); 652 primaries.clear(); 653 /* // Uzhi 31.07.09 654 std::for_each(targets.begin(), targets.end(), DeleteVSplitableHadron()); 655 targets.clear(); 656 */ // Uzhi 31.07.09 657 // theParticipants.theNucleus->StartLoop(); 658 659 //G4cout<<"NumberOfInvolvedNucleon "<<NumberOfInvolvedNucleon<<G4endl; 660 G4VSplitableHadron * aNucleon = 0; 661 for(G4int i=0; i < NumberOfInvolvedNucleon; i++) 662 { 663 aNucleon = TheInvolvedNucleon[i]->GetSplitableHadron(); 664 if(aNucleon) 665 { 666 if(aNucleon->GetStatus() != 0) delete aNucleon; 667 // if(aNucleon->GetStatus() == 2) DeleteVSplitableHadron()(aNucleon); 668 } 669 } 670 671 NumberOfInvolvedNucleon=0; 672 //G4cout<<"Out of Excit"<<G4endl; G4int Uzhi; G4cin>>Uzhi; 673 return false; 674 } // End of if(!Successfull) 675 676 return true; 607 } // end of while (theParticipants.Next()) 608 // } while (Successfull == false); // Closed 15.12.09 609 return Successfull; 677 610 } 678 611 // ------------------------------------------------------------ … … 682 615 // be duplicate in the List ( to unique G4VSplitableHadrons) 683 616 684 //G4cout<<"In build string"<<G4endl;685 686 617 G4ExcitedStringVector * strings; 687 618 strings = new G4ExcitedStringVector(); 688 619 689 620 std::vector<G4VSplitableHadron *> primaries; 690 std::vector<G4VSplitableHadron *> targets;691 std::vector<G4Nucleon *> TargetNucleons; // Uzhi 16.07.09692 621 693 622 G4ExcitedString * FirstString(0); // If there will be a kink, 694 G4ExcitedString * SecondString(0); // two strings will be produ sed.623 G4ExcitedString * SecondString(0); // two strings will be produced. 695 624 696 625 theParticipants.StartLoop(); // restart a loop 697 //G4int InterCount(0); // Uzhi 626 698 627 while ( theParticipants.Next() ) 699 628 { … … 703 632 if ( primaries.end() == std::find(primaries.begin(), primaries.end(), 704 633 interaction.GetProjectile()) ) 705 primaries.push_back(interaction.GetProjectile()); 706 707 if ( targets.end() == std::find(targets.begin(), targets.end(), 708 interaction.GetTarget()) ) 709 targets.push_back(interaction.GetTarget()); 710 711 if ( TargetNucleons.end() == std::find(TargetNucleons.begin(), 712 TargetNucleons.end(), 713 interaction.GetTargetNucleon()) ) 714 TargetNucleons.push_back(interaction.GetTargetNucleon()); 715 //InterCount++; 634 primaries.push_back(interaction.GetProjectile()); 716 635 } 717 636 718 719 //G4cout << "BuildStrings prim/targ/woundN " << primaries.size() << " , " <<targets.size() <<", "<<TargetNucleons.size()<< G4endl;720 721 637 unsigned int ahadron; 722 // Only for hA-interactions Uzhi -------------------------------------723 //G4int StringN(0);724 //G4cout<<"Proj strings -----------------------"<<G4endl;725 638 for ( ahadron=0; ahadron < primaries.size() ; ahadron++) 726 639 { 727 //G4cout<<" string# "<<ahadron<<" "<<primaries[ahadron]->Get4Momentum()<<G4endl;728 640 G4bool isProjectile(0); 729 641 if(primaries[ahadron]->GetStatus() == 1) {isProjectile=true; } … … 734 646 FirstString, SecondString, 735 647 theParameters); 736 //G4cout<<"1str 2str "<<FirstString<<" "<<SecondString<<G4endl; 648 737 649 if(FirstString != 0) strings->push_back(FirstString); 738 650 if(SecondString != 0) strings->push_back(SecondString); 739 740 //StringN++; G4cout<<"Proj string "<<strings->size()<<G4endl;741 651 } 742 //G4cout<<"Targ strings ------------------------------"<<G4endl; 743 for ( ahadron=0; ahadron < targets.size() ; ahadron++) 744 { 745 //G4cout<<"targets[ahadron]->GetStatus() "<<targets[ahadron]->GetStatus()<<G4endl; 746 if(targets[ahadron]->GetStatus() == 1) // Uzhi 17.07.09 747 { 748 G4bool isProjectile=false; 749 FirstString=0; SecondString=0; 750 theExcitation->CreateStrings(targets[ahadron], isProjectile, 751 FirstString, SecondString, 752 theParameters); 753 if(FirstString != 0) strings->push_back(FirstString); 754 if(SecondString != 0) strings->push_back(SecondString); 755 756 //StringN++; G4cout<<"Targ string"<<StringN<<G4endl; 757 } else 758 { 759 if(targets[ahadron]->GetStatus() == 0)// Uzhi 17.07.09 Nucleon was rejected 760 { 761 G4VSplitableHadron * aHit=0; // Uzhi 16.07.09 762 TargetNucleons[ahadron]->Hit(aHit); // Uzhi 16.07.09 763 } 764 } 765 } 766 767 //G4cout<<"Proj + Targ string "<<strings->size()<<G4endl; 768 //G4cout<<"Involv strings NumberOfInvolvedNucleon "<<NumberOfInvolvedNucleon<<G4endl; 652 769 653 for (G4int ahadron=0; ahadron < NumberOfInvolvedNucleon ; ahadron++) 770 654 { 771 /* 772 G4cout<<"ahadron "<<ahadron<<" Status "<< 773 TheInvolvedNucleon[ahadron]->GetSplitableHadron()->GetStatus()<< 774 TheInvolvedNucleon[ahadron]->GetSplitableHadron()->Get4Momentum()<<G4endl; 775 */ 776 if(TheInvolvedNucleon[ahadron]->GetSplitableHadron()->GetStatus() == 2) 777 { 778 //StringN++; G4cout<<"Invol string "<<StringN<<G4endl; 655 if(TheInvolvedNucleon[ahadron]->GetSplitableHadron()->GetStatus() !=0) //== 2) 656 { 779 657 G4bool isProjectile=false; 780 658 FirstString=0; SecondString=0; … … 786 664 if(FirstString != 0) strings->push_back(FirstString); 787 665 if(SecondString != 0) strings->push_back(SecondString); 788 789 // strings->push_back(theExcitation->String(790 // TheInvolvedNucleon[ahadron]->GetSplitableHadron(),791 // isProjectile));792 }793 //G4cout<<"Proj + Targ+Involved string "<<strings->size()<<G4endl;794 /*795 else796 {797 G4cout<<"Else ????????????"<<G4endl;798 if(targets[ahadron]->GetStatus() == 0)// Uzhi 17.07.09 Nucleon was rejected799 {800 G4VSplitableHadron * aHit=0; // Uzhi 16.07.09801 TargetNucleons[ahadron]->Hit(aHit); // Uzhi 16.07.09802 }803 666 } 804 */805 667 } 806 807 //G4int Uzhi; G4cin>>Uzhi;808 668 809 669 std::for_each(primaries.begin(), primaries.end(), DeleteVSplitableHadron()); 810 670 primaries.clear(); 811 812 std::for_each(targets.begin(), targets.end(), DeleteVSplitableHadron());813 targets.clear();814 671 815 672 return strings; … … 826 683 { 827 684 G4Nucleon * aNucleon = TheInvolvedNucleon[i]; 828 G4LorentzVector tmp=aNucleon->Get4Momentum()-DeltaPResidualNucleus; 685 // G4LorentzVector tmp=aNucleon->Get4Momentum()-DeltaPResidualNucleus; 686 G4LorentzVector tmp=-DeltaPResidualNucleus; 829 687 aNucleon->SetMomentum(tmp); 830 688 aNucleon->SetBindingEnergy(DeltaExcitationE); … … 837 695 { // @@ this method is used in FTFModel as well. Should go somewhere common! 838 696 839 G4double Pt2; 840 Pt2 = -AveragePt2 * std::log(1. + G4UniformRand() * 697 G4double Pt2(0.); 698 if(AveragePt2 <= 0.) {Pt2=0.;} 699 else 700 { 701 Pt2 = -AveragePt2 * std::log(1. + G4UniformRand() * 841 702 (std::exp(-maxPtSquare/AveragePt2)-1.)); 842 703 } 843 704 G4double Pt=std::sqrt(Pt2); 844 705 G4double phi=G4UniformRand() * twopi; -
trunk/source/processes/hadronic/models/parton_string/diffraction/src/G4FTFParameters.cc
r1196 r1228 1 //******************************************************************** 1 // 2 // ******************************************************************** 2 3 // * License and Disclaimer * 3 4 // * * … … 24 25 // 25 26 // 26 // $Id: G4FTFParameters.cc,v 1.1 1 2009/10/25 10:50:54 vuzhinskExp $27 // GEANT4 tag $Name: geant4-09-03 -cand-01$27 // $Id: G4FTFParameters.cc,v 1.13 2009/12/16 17:51:15 gunter Exp $ 28 // GEANT4 tag $Name: geant4-09-03 $ 28 29 // 29 30 … … 31 32 32 33 #include "G4ios.hh" 33 #include <utility> // Uzhi 29.03.0834 #include <utility> 34 35 35 36 G4FTFParameters::G4FTFParameters() … … 231 232 SetProjMinNonDiffMass(0.3); // GeV 232 233 SetProbabilityOfProjDiff(0.*0.62*std::pow(s/GeV/GeV,-0.51)); // 40/32 X-dif/X-inel 233 //G4cout<<"Params "<<0.6*0.62*std::pow(s/GeV/GeV,-0.51)<<G4endl; 234 234 235 SetTarMinDiffMass(1.1); // GeV 235 236 SetTarMinNonDiffMass(1.1); // GeV 236 //G4cout<<" "<<2.7*0.62*std::pow(s/GeV/GeV,-0.51)<<G4endl; // was 2 237 //G4int Uzhi; G4cin>>Uzhi; 237 238 238 SetProbabilityOfTarDiff(2.*0.62*std::pow(s/GeV/GeV,-0.51)); // 40/32 X-dif/X-inel 239 239 240 /*241 SetProjMinDiffMass(0.5);242 SetProjMinNonDiffMass(0.3); // Uzhi 12.06.08243 SetProbabilityOfProjDiff(0.05);244 SetProbabilityOfTarDiff(0.05);245 */246 240 SetAveragePt2(0.3); // GeV^2 247 241 } -
trunk/source/processes/hadronic/models/parton_string/diffraction/src/G4FTFParticipants.cc
r1196 r1228 25 25 // 26 26 // 27 // $Id: G4FTFParticipants.cc,v 1.1 5 2009/10/25 10:50:54vuzhinsk Exp $28 // GEANT4 tag $Name: geant4-09-03 -cand-01$27 // $Id: G4FTFParticipants.cc,v 1.16 2009/11/25 09:14:03 vuzhinsk Exp $ 28 // GEANT4 tag $Name: geant4-09-03 $ 29 29 // 30 30 // ------------------------------------------------------------ … … 72 72 73 73 void G4FTFParticipants::GetList(const G4ReactionProduct &thePrimary, 74 G4FTFParameters *theParameters) // Uzhi 29.03.0874 G4FTFParameters *theParameters) 75 75 { 76 76 … … 96 96 G4double impactY = theImpactParameter.second; 97 97 98 G4ThreeVector thePosition(impactX, impactY, -DBL_MAX); //Uzhi 24 July 0999 primarySplitable->SetPosition(thePosition); //Uzhi 24 July 0998 G4ThreeVector thePosition(impactX, impactY, -DBL_MAX); 99 primarySplitable->SetPosition(thePosition); 100 100 101 101 theNucleus->StartLoop(); 102 102 G4Nucleon * nucleon; 103 //G4int InterNumber=0; // Uzhi 104 G4int NucleonNumber(0); // Uzhi 105 //while ( (nucleon=theNucleus->GetNextNucleon())&& (InterNumber < 1) ) // Uzhi 106 while ( (nucleon=theNucleus->GetNextNucleon()) ) // Uzhi 103 104 while ( (nucleon=theNucleus->GetNextNucleon()) ) 107 105 { 108 //G4cout<<"Nucl# "<<NucleonNumber<<G4endl; // Vova109 106 G4double impact2= sqr(impactX - nucleon->GetPosition().x()) + 110 107 sqr(impactY - nucleon->GetPosition().y()); 111 108 112 // if ( theParameters->GetInelasticProbability(impact2/fermi/fermi) // Uzhi 29.03.08 113 if ( theParameters->GetProbabilityOfInteraction(impact2/fermi/fermi) // Uzhi 29.03.08 109 if ( theParameters->GetProbabilityOfInteraction(impact2/fermi/fermi) 114 110 > G4UniformRand() ) 115 111 { 116 //InterNumber++;117 /* // Uzhi 20 July 2009118 if ( nucleusNeedsShift )119 { // on the first hit, shift nucleus120 nucleusNeedsShift = false;121 theNucleus->DoTranslation(G4ThreeVector(-1*impactX,-1*impactY,0.));122 impactX=0;123 impactY=0;124 }125 */ // Uzhi 20 July 2009126 //G4cout<<" Interact"<<G4endl;127 112 primarySplitable->SetStatus(1); // It takes part in the interaction 128 113 … … 130 115 if ( ! nucleon->AreYouHit() ) 131 116 { 132 //G4cout<<"Part "<<nucleon->Get4Momentum()<<G4endl;133 117 targetSplitable= new G4DiffractiveSplitableHadron(*nucleon); 134 118 nucleon->Hit(targetSplitable); 135 nucleon->SetBindingEnergy(3.*nucleon->GetBindingEnergy()); // Uzhi 5.10.09119 nucleon->SetBindingEnergy(3.*nucleon->GetBindingEnergy()); 136 120 targetSplitable->SetStatus(1); // It takes part in the interaction 137 121 } … … 142 126 theInteractions.push_back(aInteraction); 143 127 } 144 NucleonNumber++; // Uzhi145 128 } 146 // // Uzhi 129 147 130 // G4cout << "Number of Hit nucleons " << theInteractions.size()<<G4endl; // entries() 148 131 // << "\t" << impactX/fermi << "\t"<<impactY/fermi 149 132 // << "\t" << std::sqrt(sqr(impactX)+sqr(impactY))/fermi <<G4endl; 150 // // Uzhi 133 151 134 } 152 135
Note: See TracChangeset
for help on using the changeset viewer.