/* Born approximation (Tsai (RMP 1974) Eq. 2.7) using Mathmore. */ #include "TMath.h" #include "TH1.h" #include "TCanvas.h" #include "TGraph.h" #include "TObject.h" #include "TROOT.h" #include "TLegend.h" #include "TStopwatch.h" //#include "TLabel.h" #include "Math/Functor.h" #include "Math/WrappedFunction.h" #include "Math/IFunction.h" #include "Math/Integrator.h" #include "Math/GSLMCIntegrator.h" #include "Math/WrappedParamFunction.h" #include const double pi = TMath::Pi(); const double zee = 4.; const double ay = 9.; const double mnuc = ay*0.938; const double dee = (0.162)*pow(ay, -2./3.); //GeV2 const double mu = 0.105; const double me = 0.511E-3; //electron mass in GeV const double alpha = 1./(137.036); const double zeered = pow(zee*alpha, 2.); // integration order: cosplus, costheta, pee // parameter order: kay, pee, costheta double crossdiff(const double *x, const double *params){ // 1st integral is over costheta+ = x // 2nd integral is over 2pi*costheta (emission solid angle) // 3rd integral is over p, muon momentum // only do coherent scattering, so no change in final nuclear mass double cosplus = x[0]; //costheta+ double phi = x[1]; //costheta double xi = x[2]; //pee // parameters are: k, p, costheta (can define lab frame angles so no phi-dep?) double kay = params[0]; double kthresh = (2.*mu*mu+2.*mu*mnuc)/mnuc; if (kay1.) costhetamax = 1.; double costheta = phi*(1.-costhetamax) +costhetamax; double kmin = mnuc*(mu+mueng)/(mnuc - mueng + pee*costheta); if (kmin<0.) return 0; if (kay~ 200 GeV, calls = 120000*kay ? // k = 100 GeV, calls = 10000*kay // k = 40 GeV, calls = 10000*kay // answers are different running script first time ????? unsigned int calls = 10000*kay; ROOT::Math::WrappedParamFunction<> funptr(&crossdiff, 3, 1, p); ROOT::Math::GSLMCIntegrator ig1(ROOT::Math::MCIntegration::kVEGAS, 1E-6, 1E-4, calls); ig1.SetFunction(funptr); double tsai = ig1.Integral(lo, hi); return tsai; } void tsaimc2(){ double kay = 100.; double parameters[1] = {kay}; double answer = crossdiffint(parameters); cout <