Changeset 962 for trunk/source/processes/hadronic/models/chiral_inv_phase_space/calcul/G4GlauberCrossSections.cc
- Timestamp:
- Apr 6, 2009, 12:30:29 PM (15 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/source/processes/hadronic/models/chiral_inv_phase_space/calcul/G4GlauberCrossSections.cc
r819 r962 317 317 } 318 318 319 // Initialize factorials and combinatoric coefficients 319 320 void Initialize() 320 321 { … … 448 449 G4double S = (MassN+MassN)*HadrEnergy+MassN2+MassH2;// Mondelststam s 449 450 G4double EcmH = (S-MassN2+MassH2)/2/std::sqrt(S); // CM energy of a hadron 450 G4double CMMom = std::sqrt(EcmH*EcmH-MassH2);// CM momentum451 G4double CMMom = std::sqrt(EcmH*EcmH-MassH2); // CM momentum 451 452 452 453 G4double Stot = HadrTot*mb2G2; // in GeV^-2 … … 485 486 for(G4int i=1; i<=A; i++) // @@ Make separately for n and p 486 487 { 487 N *= -Unucl*(A-i+1)*Rho2/i; 488 N *= -Unucl*(A-i+1)*Rho2/i; // Includes total cross-section 488 489 G4double N4 = 1.; 489 G4double Prod1 = std::exp(-Q2*R12B/i/4)*R12B/i; 490 G4double Prod1 = std::exp(-Q2*R12B/i/4)*R12B/i; // Includes the slope 490 491 G4double medTot = R12B/i; 491 492 for(G4int l=1; l<=i; l++) … … 556 557 ImElA0 *= Corr0; 557 558 558 return (ReElA0*ReElA0+(ImElA0+Din1)*(ImElA0+Din1))*CMMom*CMMom*mb2G2/4/pi2; 559 //return (ReElA0*ReElA0+(ImElA0+Din1)*(ImElA0+Din1))*CMMom*CMMom*mb2G2/4/pi2; // ds/do 560 return (ReElA0*ReElA0+(ImElA0+Din1)*(ImElA0+Din1))*mb2G2/(twopi+twopi); 559 561 } // End of DiffElasticCS 560 562 … … 564 566 static const G4double mb2G2 = 2.568; // Transform from mb to GeV^-2 565 567 static const G4double piMG = pi/mb2G2; // PiTrans from GeV^-2 to mb 566 G4double hMom = HadrMom/1000.; // Momentum (GeV, inputParam) 568 G4double p = HadrMom/1000.; // Momentum (GeV, inputParam) 569 G4double p2 = p*p; 567 570 G4double MassH = G4QPDGCode(hPDG).GetMass()/1000.; // Hadron Mass in GeV 568 571 G4double MassH2 = MassH*MassH; 569 G4double HadrEnergy = std::sqrt( hMom*hMom+MassH2);// Tot energy in GeV572 G4double HadrEnergy = std::sqrt(p2+MassH2); // Tot energy in GeV 570 573 //G4cout<<"G4GCS::DiffElasticCS: PGG_proj="<<hPDG<<", A_nuc= "<<A<<G4endl; 571 if(A==2 || A==3) G4Exception("G4GCS: This model does not work for nuclei with A=2, A=3");572 if(A>252) G4Exception("G4GCS: This nucleus is too heavy for the model !!!");573 574 //if(HadrEnergy-MassH<1.) G4Exception("Kin energy is too small for the model (T>1 GeV)"); 574 CalculateParameters(hPDG, HadrMom); 575 CalculateIntegralCS(A, HadrEnergy); 576 CalculateNuclearParameters(A); 575 CalculateParameters(hPDG, HadrMom); // ? 576 CalculateIntegralCS(A, HadrEnergy); // ? 577 CalculateNuclearParameters(A); // ? 577 578 G4double Q2 = aQ2/1000/1000; // in GeV^2 578 579 //G4doubleMassN = A*0.938; // @@ This is bad! 579 580 G4double MassN = A*0.9315; // @@ Atomic Unit is bad too 580 581 if(G4NucleiPropertiesTable::IsInTable(Z,A)) 581 MassN=G4NucleiProperties::GetNuclearMass(A,Z)/1000.; // Geant4 NuclearMass 582 MassN=G4NucleiProperties::GetNuclearMass(A,Z)/1000.; // Geant4 NuclearMass in GeV 582 583 G4double MassN2 = MassN*MassN; 583 584 G4double S = (MassN+MassN)*HadrEnergy+MassN2+MassH2;// Mondelststam s 584 585 G4double EcmH = (S-MassN2+MassH2)/2/std::sqrt(S); // CM energy of a hadron 585 586 G4double CMMom = std::sqrt(EcmH*EcmH-MassH2); // CM momentum 586 G4double Stot = HadrTot*mb2G2; // in GeV^-2 587 G4double Bhad = HadrSlope; // in GeV^-2 588 G4double Asq = 1+HadrReIm*HadrReIm; // |M|^2/(ImM)^2 587 //G4cout<<"P="<<CMMom<<",E="<<HadrEnergy<<",N="<<MassN2<<",h="<<MassH2<<",p="<<p<<G4endl; 588 G4double p3 = p2*p; 589 G4double p4 = p2*p2; 590 G4double sp = std::sqrt(p); 591 G4double p2s = p2*sp; 592 G4double ap = std::log(p); 593 G4double dl = ap-3.; 594 G4double dl2 = dl*dl; 595 G4double shPPTot = 2.91/(p2s+.0024)+5.+(32.+.3*dl2+23./p)/(1.+1.3/p4); // SigPP in mb 596 G4double shPNTot = 12./(p2s+.05*p+.0001/std::sqrt(sp))+.35/p 597 +(38.+.3*dl2+8./p)/(1.+1.2/p3); // SigPP in mb 598 G4double hNTot = (Z*shPPTot+(A-Z)*shPNTot)/A; // SighN in mb 599 G4double Stot = hNTot*mb2G2; // in GeV^-2 600 G4double shPPSl = 8.*std::pow(p,.055)/(1.+3.64/p4); // PPslope in GeV^-2 601 G4double shPNSl = (7.2+4.32/(p4*p4+.012*p3))/(1.+2.5/p4); // PNslope in GeV^-2 602 G4double Bhad = (Z*shPPSl+(A-Z)*shPNSl)/A; // B-slope in GeV^-2 603 G4double hNReIm = -.55+ap*(.12+ap*.0045); // Re/Im_hN in no unit 604 G4double Asq = 1+hNReIm*hNReIm; // |M|^2/(ImM)^2 589 605 G4double Rho2 = std::sqrt(Asq); // M/ImM 590 G4double R12 = R1*R1; 591 G4double R22 = R2*R2; 606 G4double r1 = 3.9*std::pow(A-1.,.309); // Positive diffractionalRadius 607 G4double r2 = 2.*std::pow(A,.36); // Negative diffraction radius 608 G4double pN = Pnucl; // Dubna value 609 //G4double pN = .4; // Screaning factor 610 G4double Ae = Aeff; // Dubna value 611 //G4double Ae = .75; // Normalization 612 G4double R12 = r1*r1; 613 G4double R22 = r2*r2; 614 G4double R1C = R12*r1; 615 G4double R2C = R22*r2; 592 616 G4double R12B = R12+Bhad+Bhad; // Slope is used 593 617 G4double R22B = R22+Bhad+Bhad; // Slope is used 594 G4double Norm = (R1 2*R1-Pnucl*R22*R2)*Aeff; // Some questionable norming595 G4double R13 = R1 2*R1/R12B;// Slope is used596 G4double R23 = Pnucl*R22*R2/R22B; // Slope isused618 G4double Norm = (R1C-pN*R2C)*Ae; // ScreanFac & NormFac are used 619 G4double R13 = R1C/R12B; // Slope is used 620 G4double R23 = pN*R2C/R22B; // Slope & ScreanFact are used 597 621 G4double norFac = Stot/twopi/Norm; // totCS (in GeV^-2) is used 598 622 G4double Unucl = norFac*R13; 599 G4double SinFi = HadrReIm/Rho2;623 G4double SinFi = hNReIm/Rho2; // Real part 600 624 G4double FiH = std::asin(SinFi); 601 625 G4double N = -1.; … … 606 630 for(G4int i=1; i<=A; i++) 607 631 { 608 N *= -Unucl*(A-i+1)*Rho2/i; 632 N *= -Unucl*(A-i+1)*Rho2/i; // TotCS is used 609 633 G4double N4 = 1.; 610 634 G4double Prod1 = std::exp(-Q2*R12B/i/4)*R12B/i; // Slope is used … … 631 655 ImElA0 *= Corr0; 632 656 633 return (ReElA0*ReElA0+ImElA0*ImElA0)*CMMom*CMMom*mb2G2/4/pi2; 657 //return (ReElA0*ReElA0+ImElA0*ImElA0)*CMMom*CMMom*mb2G2/4/pi2; // ds/do 658 return (ReElA0*ReElA0+ImElA0*ImElA0)*mb2G2/(twopi+twopi); 634 659 } // End of DiffElasticCS 635 660 … … 705 730 else if(A==12) r0 = 0.80; 706 731 else r0 = 1.16*(1.-1.16/Re2); // For other nuclei which have not been tested 707 G4double MassN = mN; // A * AtomicUnit(GeV) 732 G4double MassN = A*0.9315; // @@ Atomic Unit is bad too 708 733 G4double MassN2 = MassN*MassN; 709 734 G4double S = (MassN+MassN)*HadrEnergy+MassN2+MassH2; // Mondelstam S … … 736 761 G4double expB = std::exp(InExp); 737 762 G4double HRIE = HadrReIm*InExp; 738 ReIntegrand[i] = (1.-expB*std::cos(HRIE)); 739 ImIntegrand[i] = expB*std::sin(HRIE); 763 ReIntegrand[i] = (1.-expB*std::cos(HRIE)); // Real part of the amplitude 764 ImIntegrand[i] = expB*std::sin(HRIE); // Imaginary part of the amplit 740 765 } 741 766 InCoh = 0.; // incohirent (quasi-elastic) 742 767 ValB = -stepB; 743 for(G4int k=0; k<NptB; k++) 768 for(G4int k=0; k<NptB; k++) // Third integration (?) 744 769 { 745 770 ValB += stepB; 746 771 InCoh += Thick[k]*ValB*std::exp(-hTotG2*Thick[k]); 747 G4double J0qb = QJ0(std::sqrt(Q2)*ValB)*ValB; 772 G4double J0qb = QJ0(std::sqrt(Q2)*ValB)*ValB; // Bessel0(Q*b) 748 773 ReSum += J0qb*ReIntegrand[k]; 749 774 ImSum += J0qb*ImIntegrand[k]; … … 753 778 InCoh *= stepB*hTotG2*hTotG2*(1.+HadrReIm*HadrReIm)*std::exp(-HadrSlope*Q2)/8/mb2G2; 754 779 //G4cout<<"GHAD:RS="<<ReSum<<",IS="<<ImSum<<",CM="<<CMMom<<",st="<<stepB<<G4endl; 755 return (ReSum*ReSum+ImSum*ImSum)*m2G10*CMMom*CMMom*stepB*stepB/twopi; 780 //return (ReSum*ReSum+ImSum*ImSum)*m2G10*CMMom*CMMom*stepB*stepB/twopi; // ds/do 781 return (ReSum*ReSum+ImSum*ImSum)*m2G10*stepB*stepB/12; // ds/dt 756 782 } 757 783 … … 763 789 static const G4double hc2 = .3893793; // Transform from GeV^-2 to mb 764 790 static const G4double mb2G2 = 1./hc2; // Transform from mb to GeV^-2 765 //static const G4double mb2G2 = 2.568; // Transform from mb to GeV^-2766 791 static const G4double f22mb = 10; // Transform from fermi^2 to mb 767 792 static const G4double f22G2 = f22mb*mb2G2; // Transform from fm2 to GeV^-2 768 793 static const G4double f2Gm1 = std::sqrt(f22G2); // Transform from fm to GeV^-1 769 794 G4double Re = std::pow(A,.33333333); // A-dep coefficient 770 G4double Lim = 50.*Re; // Integration accuracy limit795 G4double Lim = 100*Re; // Integration accuracy limit 771 796 G4double Tb[Npb]; // Calculated T(b) array 772 797 G4QBesIKJY QI0(BessI0); // I0 Bessel function 773 798 G4QBesIKJY QJ0(BessJ0); // J0 Bessel function 774 hMom = hMom/1000.; // Momentum (GeV, inputParam) 799 G4double p = hMom/1000.; // Momentum (GeV, inputParam) 800 G4double p2 = p*p; 775 801 G4double MassH = G4QPDGCode(hPDG).GetMass()/1000.; // Hadron Mass in GeV 776 802 G4double MassH2 = MassH*MassH; // Squared mass of the hadron 777 G4double hEnergy = std::sqrt( hMom*hMom+MassH2);// Tot energy in GeV803 G4double hEnergy = std::sqrt(p2+MassH2); // Tot energy in GeV 778 804 G4double Q2 = aQ2/1000000.; // -t in GeV 779 G4double MassN = mN; // AtomicUnit(GeV)[prototype] 805 G4double MassN = A*0.9315; // @@ Atomic Unit is bad too 806 if(G4NucleiPropertiesTable::IsInTable(Z,A)) 807 MassN=G4NucleiProperties::GetNuclearMass(A,Z)/1000.; // Geant4 NuclearMass in GeV 780 808 G4double MassN2 = MassN*MassN; // Squared mass of the target 781 809 G4double S = (MassN+MassN)*hEnergy+MassN2+MassH2;// Mondelstam s 782 810 G4double EcmH = (S-MassN2+MassH2)/2/std::sqrt(S); // Hadron CM Energy 783 811 G4double CMMom = std::sqrt(EcmH*EcmH-MassH2); // CM momentum (to norm CS) 784 // @@ Temporary only for nucleons 785 //G4cout<<"CHPS:E="<<hEnergy<<",dM="<<2*MassN<<",sN="<<MassN2<<",sH="<<MassH2<<G4endl; 786 G4double shNTot = 5.2+5.2*std::log(hEnergy)+51*std::pow(hEnergy,-.35); // SighN in mb 787 G4double shNSl = 5.44+.88*std::log(S); // B-slope in GeV^-2 788 G4double shNReIm = .13*std::log(S/350)*std::pow(S,-.18); // Re/Im_hN in no unit 812 //G4cout<<"2: P="<<CMMom<<",E="<<hEnergy<<",N="<<MassN2<<",h="<<MassH2<<",p="<<p<<G4endl; 813 // The mean value of the total can be used 814 G4double p3 = p2*p; 815 G4double p4 = p2*p2; 816 G4double sp = std::sqrt(p); 817 G4double p2s = p2*sp; 818 G4double ap = std::log(p); 819 G4double dl = ap-3.; 820 G4double dl2 = dl*dl; 821 G4double shPPTot = 2.91/(p2s+.0024)+5.+(32.+.3*dl2+23./p)/(1.+1.3/p4); // SigPP in mb 822 G4double shPNTot = 12./(p2s+.05*p+.0001/std::sqrt(sp))+.35/p 823 +(38.+.3*dl2+8./p)/(1.+1.2/p3); // SigPP in mb 824 G4double shNTot = (Z*shPPTot+(A-Z)*shPNTot)/A; // SighN in mb 825 #ifdef debug 826 G4cout<<"CHIPS:SI,p="<<p<<",n="<<shNTot<<",P="<<shPPTot<<",N="<<shPNTot 827 <<",Z="<<Z<<",A="<<A<<G4endl; 828 #endif 829 G4double shPPSl = 8.*std::pow(p,.055)/(1.+3.64/p4); // PPslope in GeV^-2 830 G4double shPNSl = (7.2+4.32/(p4*p4+.012*p3))/(1.+2.5/p4); // PNslope in GeV^-2 831 G4double shNSl = (Z*shPPSl+(A-Z)*shPNSl)/A; // B-slope in GeV^-2 832 #ifdef debug 833 G4cout<<"CHIPS:SL,n="<<shNSl<<",P="<<shPPSl<<",N="<<shPNSl<<G4endl; 834 #endif 835 G4double shNReIm = -.55+ap*(.12+ap*.0045); // Re/Im_hN in no unit 789 836 //G4cout<<"CHPS: s="<<S<<",T="<<shNTot<<",R="<<shNReIm<<",B="<<shNSl<<G4endl; 790 837 // @@ End of temporary ^^^^^^^ 791 838 G4double dHS = shNSl+shNSl; // Working: doubled B-slope 792 //G4double rAfm = 0.;793 //if (A==208) rAfm = 1.125*Re;794 //else if(A==90) rAfm = 1.12*Re;795 //else if(A==64) rAfm = 1.1*Re;796 //else if(A==58) rAfm = 1.09*Re;797 //else if(A==48) rAfm = 1.07*Re;798 //else if(A==40) rAfm = 1.15*Re;799 //else if(A==28) rAfm = 0.93*Re;800 //else if(A==16) rAfm = 0.92*Re;801 //else if(A==12) rAfm = 0.80*Re;802 //else rAfm = 1.16*(Re-1.16/Re); // For other nuclei803 //G4double stepB = 2.5*rAfm*f2Gm1/(Npb-1); // in GeV^-1, step of integral804 839 G4double stepB = (Re+Re+2.7)*f2Gm1/(Npb-1); // in GeV^-1, step of integral 805 840 G4double hTotG2 = shNTot*mb2G2; // sigma_hN in GeV^-2 … … 807 842 G4double ImSum = 0.; // Integration of ImPart of Amp 808 843 G4double ValB = -stepB; 809 for(G4int i=0; i<Npb; i++) 844 for(G4int i=0; i<Npb; i++) // First integration over b 810 845 { 811 846 ValB += stepB; // An incident parameter 812 G4double ValB2 = ValB*ValB; // A working value 813 G4double IPH = ValB/shNSl; // A working value 847 G4double ValB2 = ValB*ValB; // A working value b^2 848 G4double IPH = ValB/shNSl; // A working value slope 814 849 G4double Integ = 0.; // Integral over ImpactParam. 815 850 G4double ValS = 0.; // Prototype of ImpactParameter 816 for(G4int j=1; j<Npb; j++) 851 for(G4int j=1; j<Npb; j++) // Second integration over b 817 852 { 818 853 ValS += stepB; // back to fm // Impact parameter GeV^-1 819 854 if(!i) Tb[j] = CHIPS_Tb(A,ValS/f2Gm1)/f22G2; // GeV^2, calculate only once 820 855 //if(!i) Tb[j] = Thickness(A,ValS/f2Gm1,rAfm)/f22G2; // Calculate T(b) only once 821 G4double FunS = IPH*ValS; // Working product822 if(FunS > Lim) break; // (?)823 Integ += ValS*std::exp(-(ValS*ValS+ValB2)/dHS)*QI0(FunS)*Tb[j]; 856 G4double FunS = IPH*ValS; // b1*b2/slope 857 if(FunS > Lim) break; // To avoid NAN 858 Integ += ValS*std::exp(-(ValS*ValS+ValB2)/dHS)*QI0(FunS)*Tb[j]; // BessI0 824 859 } 825 G4double InExp = -hTotG2*Integ*stepB/dHS; // Working product826 G4double expB = std::exp(InExp); // Workung sqrt860 G4double InExp = -hTotG2*Integ*stepB/dHS; // Integrated absorption 861 G4double expB = std::exp(InExp); // Exponential absorption 827 862 G4double HRIE = shNReIm*InExp; // Phase shift 828 G4double J0qb = QJ0(std::sqrt(Q2)*ValB)*ValB; 863 G4double J0qb = QJ0(std::sqrt(Q2)*ValB)*ValB; // Bessel0(Q*b) 829 864 ReSum += J0qb*(1.-expB*std::cos(HRIE)); 830 865 ImSum += J0qb*expB*std::sin(HRIE); … … 832 867 //G4cout<<"CHPS:RS="<<ReSum<<",IS="<<ImSum<<",CM="<<CMMom<<",st="<<stepB<<G4endl; 833 868 //return (ReSum*ReSum+ImSum*ImSum)*mb2G2*CMMom*CMMom*stepB*stepB/twopi; 834 return (ReSum*ReSum+ImSum*ImSum)*f22G2*CMMom*CMMom*stepB*stepB/twopi; 835 } 869 //return (ReSum*ReSum+ImSum*ImSum)*f22G2*CMMom*CMMom*stepB*stepB/twopi; // ds/do 870 return (ReSum*ReSum+ImSum*ImSum)*f22G2*stepB*stepB/12; // ds/dt 871 } 872 873 // Separate quasielastic calculation 836 874 G4double CHIPSDifQuasiElasticCS(G4int hPDG, G4double hMom, G4int A, G4int Z, G4double aQ2) 837 875 {// ================================================================================= 838 876 static const G4int Npb = 500; // A#of intergation points 839 //static const G4double mN = .9315; // Atomic Unit GeV840 static const G4double mN = .938; // Mass of proton GeV841 877 static const G4double hc2 = .3893793; // Transform from GeV^-2 to mb 842 878 static const G4double mb2G2 = 1./hc2; // Transform from mb to GeV^-2 … … 845 881 static const G4double f22G2 = f22mb*mb2G2; // Transform from fm2 to GeV^-2 846 882 static const G4double f2Gm1 = std::sqrt(f22G2); // Transform from fm to GeV^-1 847 hMom = hMom/1000.; // Momentum (GeV, inputParam) 848 G4double MassH = G4QPDGCode(hPDG).GetMass()/1000.; // Hadron Mass in GeV 849 G4double MassH2 = MassH*MassH; // Squared mass of the hadron 850 G4double hEnergy = std::sqrt(hMom*hMom+MassH2); // Tot energy in GeV 883 G4double p = hMom/1000.; // Momentum (GeV, inputParam) 884 G4double p2 = p*p; 851 885 G4double Q2 = aQ2/1000000.; // -t in GeV 852 G4double MassN = mN; // A*AtomicUnit(GeV)[prototype]853 //if(G4NucleiPropertiesTable::IsInTable(Z,A))854 // MassN=G4NucleiProperties::GetNuclearMass(A,Z)/A/1000.; // Geant4 NuclearMass/A855 G4double MassN2 = MassN*MassN; // Squared mass of the target856 G4double S = (MassN+MassN)*hEnergy+MassN2+MassH2;// Mondelstam s857 886 // @@ Temporary only for nucleons 858 G4double shNTot = 5.2+5.2*std::log(hEnergy)+51*std::pow(hEnergy,-.35); // SighN in mb 859 G4double shNSl = 6.44+.88*std::log(S); // B-slope in GeV^-2 860 G4double shNReIm = .13*std::log(S/350)*std::pow(S,-.18); // Re/Im_hN in no unit 887 G4double p3 = p2*p; 888 G4double p4 = p2*p2; 889 G4double sp = std::sqrt(p); 890 G4double p2s = p2*sp; 891 G4double ap = std::log(p); 892 G4double dl = ap-3.; 893 G4double dl2 = dl*dl; 894 G4double shPPTot = 2.91/(p2s+.0024)+5.+(32.+.3*dl2+23./p)/(1.+1.3/p4); // SigPP in mb 895 G4double shPNTot = 12./(p2s+.05*p+.0001/std::sqrt(sp))+.35/p 896 +(38.+.3*dl2+8./p)/(1.+1.2/p3); // SigPP in mb 897 G4double shNTot = (Z*shPPTot+(A-Z)*shPNTot)/A; // SighN in mb 898 G4double shPPSl = 8.*std::pow(p,.055)/(1.+3.64/p4); // PPslope in GeV^-2 899 G4double shPNSl = (7.2+4.32/(p4*p4+.012*p3))/(1.+2.5/p4); // PNslope in GeV^-2 900 G4double shNSl = (Z*shPPSl+(A-Z)*shPNSl)/A; // B-slope in GeV^-2 901 G4double shNReIm = -.55+ap*(.12+ap*.0045); // Re/Im_hN in no unit 861 902 // @@ End of temporary ^^^^^^^ 862 903 G4double Re = std::pow(A,.33333333); // A-dep coefficient … … 1009 1050 //const G4int A[na]={4,9,12,16,27,48,58,64,120,181,207,238}; // A's of target nuclei 1010 1051 // He Al Pb 1011 const G4int A[na]={ 208}; // A's of target nuclei1052 const G4int A[na]={119}; // A's of target nuclei 1012 1053 // p n pi+ pi- K+ K- antip 1013 1054 const G4int pdg[np]={2212,2112,211,-211,321,-321,-2212}; // projectiles 1014 const G4double mom[nm]={1 090.}; // momentum in MeV/c1055 const G4double mom[nm]={120.}; // momentum in MeV/c 1015 1056 #ifdef integrc 1016 1057 for(G4int ip=0; ip<np; ip++) … … 1049 1090 //const G4int Z[na]={2,4, 6, 8,13,22,28,29, 50, 73, 82, 92}; // Z's of target nuclei 1050 1091 // He Al Pb 1051 const G4int Z[na]={ 82}; // Z's of target nuclei1092 const G4int Z[na]={50}; // Z's of target nuclei 1052 1093 // Test of differential ellastic cross sections 1053 1094 Initialize(); … … 1066 1107 { 1067 1108 G4double Sig1 = CHIPSDiffElCS(PDG, mom[im], Z[ia], A[ia], t[it]); 1068 //G4double Sig1 = DiffElasticCS(PDG, mom[im], A[ia], t[it]); 1109 //G4double Sig1 = DiffElasticCS(PDG, mom[im], A[ia], t[it]); //Doesn't work 1069 1110 G4double Sig2 = CoherentDifElasticCS(PDG, mom[im], A[ia], t[it]); 1070 1111 G4double Sig3 = CHIPSDifElasticCS(PDG, mom[im], A[ia], Z[ia], t[it]);
Note: See TracChangeset
for help on using the changeset viewer.