- Timestamp:
- Apr 6, 2009, 12:30:29 PM (15 years ago)
- Location:
- trunk/source/processes/hadronic/models/chiral_inv_phase_space/calcul
- Files:
-
- 2 edited
-
G4GlauberCrossSections.cc (modified) (17 diffs)
-
G4NuMuXQ2Integration.cc (modified) (10 diffs)
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]); -
trunk/source/processes/hadronic/models/chiral_inv_phase_space/calcul/G4NuMuXQ2Integration.cc
r819 r962 62 62 else 63 63 { 64 double r=0.; 65 double max=1.; 66 double H=1.22; 67 if(A==1) // Proton 68 { 69 r=std::sqrt(Q2/1.66); 70 max=.5; 71 } 72 else if(A<13) // Light nuclei 73 { 74 double f=Q2/4.62; 75 r=f*f; 76 max=.3; 77 if(A>2) H=1.; 78 } 79 else if(A>0) // Heavy nuclei 80 { 81 double f=Q2/3.4; 82 double ff=f*f; 83 r=ff*ff; 84 max=.5; 85 H=1.; 86 } 87 else G4cout<<"strucf: A="<<A<<" <= 0"<< G4endl; 64 //double r=0.; 65 double min=.077; // Delta(0) 66 double max=.47; // Hard Pomeron 67 //double H=1.22; 68 if(Q2<=0.) G4cout<<"strucf: Q2="<<Q2<<" <= 0"<< G4endl; 69 double Q=std::sqrt(Q2); 70 //if(A==1) // Proton 71 //{ 72 // r=std::sqrt(Q2/1.66); 73 // max=.5; 74 //} 75 //else if(A<13) // Light nuclei 76 //{ 77 // double f=Q2/4.62; 78 // r=f*f; 79 // max=.3; 80 // if(A>2) H=1.; 81 //} 82 //else if(A>0) // Heavy nuclei 83 //{ 84 // double f=Q2/3.4; 85 // double ff=f*f; 86 // r=ff*ff; 87 // max=.5; 88 // H=1.; 89 //} 90 //else G4cout<<"strucf: A="<<A<<" <= 0"<< G4endl; 88 91 // 89 92 N=3.+.3581*std::log(1.+Q2/.04); // a#of partons in the nonperturbative phase space 90 Del=(1.+r)/(12.5+r/max); 91 double S=std::pow(1.+.6/Q2,-1.-Del); 92 D=H*S*(1.-.5*S); 93 Del=min+(max-min)/(1.+5./Q); 94 D=0.68*std::pow(1.+.145/Q2,-1.-Del); // 0.68=1-0.34, m2=.145 GeV2 93 95 V=3*(1.-D)*(N-1.); 94 double cc=Q2/.08; 95 double cc2=cc*cc; 96 //double C=(1.+cc2)/(1.+cc2/.24); // Weak? 97 double C=(1.+cc2)/(1.+cc2/.24)/(1.+Q2/21.6); // EM 98 double c3=C+C+C; 96 double c3=.75; // 3*0.25 99 97 double uu=std::exp(lGam(N-Del)-lGam(N-1.)-lGam(1.-Del))/N; 100 98 U2=(c3+N-3.)*uu; 101 99 U3=c3*uu; 102 //UU=uu+uu+uu; // @@103 100 mA = A; 104 101 mQ2 = Q2; … … 125 122 //else 126 123 xf3= U3*pp+dir; 127 fL = per/ 4.;124 fL = per/5.; // 20% 128 125 return; 129 126 } … … 157 154 const double sik=GF2*hc2/pif; // precalculated coefficient 158 155 //const double mpi=.1349766; // pi0 meson mass in GeV 159 const double mpi=.13957018; // charged pi meson mass in GeV156 //const double mpi=.13957018; // charged pi meson mass in GeV 160 157 //const double mpi2=mpi*mpi; // m_pi^2 in GeV^2 161 158 //const double me=.00051099892; // electron mass in GeV … … 168 165 //const double mtau2=mtau*mtau; // m_tau^2 in GeV^2 169 166 //const double hmtau2=mtau2/2; // .5*m_e^2 in GeV^2 170 //const double mp=.93827203; // proton mass in GeV171 //const double mn=.93956536; // neutron mass in GeV167 const double mp=.93827203; // proton mass in GeV 168 const double mn=.93956536; // neutron mass in GeV 172 169 //const double md=1.87561282; // deuteron mass in GeV 173 const double MN=.931494043; // Nucleon mass (inside nucleus, atomic mass unit, GeV)174 //const double MN=(mn+mp)/2; // Nucleon mass (mean free) in GeV175 //const double MD=1.232; // proton mass in GeV170 //const double MN=.931494043; // Nucleon mass (inside nucleus, atomic mass unit, GeV) 171 const double MN=(mn+mp)/2; // Nucleon mass (mean free) in GeV 172 const double MD=1.232; // proton mass in GeV 176 173 //const double mp2=mp*mp; // m_p^2 in GeV^2 177 174 const double MN2=MN*MN; // M_N^2 in GeV^2 … … 184 181 // 185 182 //const double mc=.3; // parameter of W>M+mc cut for Quasi-Elastic/Delta 186 const double mc=mpi; // parameter of W>M+mc cut for Quasi-Elastic/Delta 183 const double mc=std::sqrt((MN*MN+MD*MD)/2)-MN; // Squared mean between N and \Delta 184 //const double mc=mpi; // parameter of W>M+mc cut for Quasi-Elastic/Delta 187 185 const double mcV=(dMN+mc)*mc; // constant of W>M+mc cut for Quasi-Elastic 188 186 //std::ofstream fileNuMuX("NuMuXQ2.out", std::ios::out); … … 198 196 int A=12; // Neucleus for which calculations should be done 199 197 double lEnuMin=0; // LogLog of Minimum energy of neutrino 200 double lEnuMax=std::log(1.+std::log(300./EminMu)); // LogLog of Maximum energy of neutrino201 int nE= 31;198 double lEnuMax=std::log(1.+std::log(300./EminMu)); // LogLog of MaximumEnergy of neutrino 199 int nE=49; // Number of points 202 200 double dlE=(lEnuMax-lEnuMin)/nE; 203 201 lEnuMin+=dlE/10; … … 241 239 double dik=MW4/Q2M/Q2M; 242 240 double qmc=Q2+mcV; 243 double lXQES=std::log((std::sqrt(qmc*qmc+Q2*fMN2)-qmc)/dMN2); // Quasielastic boundary241 double lXQES=std::log((std::sqrt(qmc*qmc+Q2*fMN2)-qmc)/dMN2);//QuasielasticBoundary 244 242 //double lXQES=log(Q2/(Q2+mcV)); // Quasielastic boundary (W=MN+m_c) 245 243 //double xN=Q2/dME; … … 248 246 double lXmin=std::log(xN); 249 247 // ****** QE ******** 250 //if(lXQES>lXmin) lXmin=lXQES; // A cut which leaves only QES248 if(lXQES>lXmin) lXmin=lXQES; // A cut which leaves only QES >>>>>>>>>>>>>>>> 251 249 // *** End of QE^^^^^ 252 250 double lXmax=0.; // QES is in DIS … … 262 260 { 263 261 getFun(A, lX, Q2, f); 264 DISxint+=f[0]+f[0]+xN*(f[1]+f[1]+xN*f[3]); // neutrino 265 //DISxint+=f[0]+f[0]+xN*(f[2]+f[2]+xN*f[4]); // anti-neutrino 262 // ***** Neutrino/Antineutrino switch ******>>>>>>>>>>>>>>>>>>>>>>>>>>>> 263 //DISxint+=f[0]+f[0]+xN*(f[1]+f[1]+xN*f[3]); // neutrino 264 DISxint+=f[0]+f[0]+xN*(f[2]+f[2]+xN*f[4]); // anti-neutrino 266 265 //G4cout<<f[0]<<","<<f[1]<<","<<f[2]<<","<<f[3]<<","<<f[4]<<G4endl; 267 266 } … … 279 278 //===== tot/qe choice ==== 280 279 DIStsig*=sik/Enu; 281 G4cout<<"***total*** E="<<Enu<<",sig/E="<<DIStsig<<G4endl; 280 //G4cout<<"***total-neutrino*** E="<<Enu<<" ,sig/E="<<DIStsig<<G4endl; 281 G4cout<<"***total-antineutrino*** E="<<Enu<<" ,sig/E="<<DIStsig<<G4endl; 282 //G4cout<<"***quasiel-nu*** E="<<Enu<<" ,sig/E="<<DIStsig<<G4endl; 283 //G4cout<<"***quasiel-antinu*** E="<<Enu<<" ,sig/E="<<DIStsig<<G4endl; 282 284 //................... 283 285 //DIStsig*=sik;
Note: See TracChangeset
for help on using the changeset viewer.
