#include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include "TH2F.h" #include "TCanvas.h" #include "TPad.h" #include "TSpectrum.h" #include #include //serve per basename #include // serve per stringstream #include #include using namespace std; // GLOBAL values ///////////////////////////////////////////////// string home = string(getenv("HOME")); string path="Disang/"; string FileZiegler; double Zgl[10][9]; const int Nelem = 2; int Zt[Nelem] = {1, 40}; int At[Nelem] = {2, 91}; double w[Nelem] = {2, 1}; double Zp = 1; double Ap = 1; double Dpercm2 = 2.26e18; //double Dpercm2 = 2.05e18; double mp = 938.272034*1e3; double md = 1875.612845*1e3; double mu = mp*md/(mp+md); double alpha = pow(137, -1); double pi = acos(-1); double tau; // ==>> mu/Tcm == mp/T //double Exponent = 2*pi*alpha*Zp*Zt[0]*sqrt(mu/2); double Exponent = 2*pi*alpha*Zp*Zt[0]*sqrt(mp/2); //double Exponent = 25.6464; //const int Npar = 3; string DirDat = "Disang/"; double mass(int Z, int A, int flag) { int z, a; double amu, massat, massnucl, exc; double answ; string nucleus; int check = 0; ifstream in; in.open("tabmass2021.dat"); while (!in.eof()) { in>>nucleus>>z>>a>>exc>>amu>>massat>>massnucl; if(z==Z && a==A) { if(flag == 0) answ = amu; if(flag == 1) answ = massat; if(flag == 2) answ = massnucl; check = 1; break; } if(check == 1) break; } return answ; } /////////////////////////////////////////////////////////////////// // function stop // input values: Z, A of proj, Z, A of target, Tp of proj (keV) // returns stopping cross section in units eV/(10^15 atoms/cm^2 /////////////////////////////////////////////////////////////////// double stop(int Zp, int Ap, int Zt, int At, int J, double Tp) { //warning: Tp must be in keV double ziegler[9]; for (int j=0; j<9; j++) ziegler[j] = Zgl[J][j]; double epsil, Kp, L, bet2, Slow, Shigh; if(Zp == 1) { // Hydrogen-like projectile.. Kp=Tp/Ap; L=log(Kp); bet2 = 1-1/pow((1+Kp/931494.062),2); if(Kp<10) { epsil = ziegler[0]*sqrt(Kp); } else if(10 <= Kp < 1e3) { Slow = ziegler[1]*pow(Kp,0.45); Shigh = ziegler[2]/Kp*log(1+ziegler[3]/Kp+ziegler[4]*Kp); epsil = Slow*Shigh/(Slow+Shigh); } else { epsil = ziegler[5]/bet2*(log(ziegler[5]*bet2/(1-bet2))-bet2); } } if(Zp>1) { // Helium-like projectile.... // Kp must be given in MeV... // Kp is scaled with masses (compared to alpha) // epsil is scaled with z^2 ratios (compared to alpha) Kp = Tp*mass(2,4,1)/Ap*1e-03; if(Kp < 1e1) { Slow = ziegler[0]*pow((Kp*1e03),ziegler[1]); Shigh = ziegler[2]/Kp*log(1+ziegler[3]/Kp+ziegler[4]*Kp); epsil = Slow*Shigh/(Slow+Shigh)*pow((Zp/2),2); } else { L = log(1/Kp); epsil = exp(ziegler[5]+ziegler[6]*L+ziegler[7]*pow(L,2)+ziegler[8]*pow(L,3))*pow((Zp/2),2); } } // epsil is given in eV/(10^15 atoms/cm^2) return(epsil); } double stopeff(double Tp) { double eps = (2*stop(Zp,Ap,Zt[0],At[0],0,Tp) + stop(Zp,Ap,Zt[1],At[1],1,Tp))/3; eps = eps*1e-15*1e-3; // keV*cm2 return eps; } /////////////////////////////////////////////////////////////////// // function dedx // Tp in keV // summing stopping over compund elements // converting eV/(10^15 atoms/cm^2) into keV*cm2/g /////////////////////////////////////////////////////////////////// double dedx(int Zp, int Ap, int Nelem, double Tp) { double SP = 0; double SP_parz; double PM = 0; for (int i=0; i=Tmin) { R = R + 1/dedx(Zp, Ap, Nelem, T); T = T - dT; } R = R*dT; return R; } /////////////////////////////////////////////////////////////////// // function out // computing resdual energy in compound // Tp in keV // dedx in keV*cm2/g /////////////////////////////////////////////////////////////////// double Tfin(int Zp, int Ap, int Nelem, double Tp) { double dx = tau*1e-3; double x = dx; double T = Tp; while (x> z; for (int k=0; k<12; k++) in0 >> temp[k]; for (int m=0; m>dummy_s>>dummy; for (int k=0; k>y[k][j]>>Dy[k][j]; } } in.close(); TCanvas *ZZ[NE]; int N; string GrName; double Y[Ndet], DY[Ndet]; TGraphErrors *G[NE]; TF1 *fit[NE]; double Ebeam; double a[NE][Npar], Da[NE][Npar]; double A0[NE], dA0[NE]; double S[NE][Ndet], dS[NE][Ndet]; double Eout[NE], Eeff[NE], Ecmeff[NE]; double W[NE][Ndet]; double Sum; for (int n=0; ncd(); Ebeam = E[n]; GrName = E_s[n]+"keV"; for (int k=0; kSetTitle(GrName.c_str()); G[n]->Draw("AP"); string ParName[Npar]; char *id = new char[10]; for (int j=0; jSetParName(j,ParName[j].c_str()); } G[n]->Fit(fit[n], "R"); fit[n]->Draw("same"); for (int j=0; jGetParameter(j); Da[n][j] = fit[n]->GetParError(j); } A0[n] = a[n][0]; dA0[n] = Da[n][0]; Eout[n] = Tfin(Zp, Ap, Nelem, E[n]); double dT = (E[n]-Eout[n])/1000; Eeff[n] = (E[n]+Eout[n])/2; Ecmeff[n] = Eeff[n]*md/(mp+md); //out<Eout[n]) { Sum = Sum + F(T); T = T - dT; } Sum = Sum*dT; for (int k=0; k