// Root includes... #include "TCanvas.h" #include "TPad.h" #include "TH2D.h" #include "TH1D.h" #include "TF1.h" #include "TMath.h" #include "TString.h" #include "TGraph.h" // Some useful constants... const Double_t v_esc = 600. * 1.e5; //Local galactic escape velocity [cm/s] const Double_t v_0 = 230. * 1.e5; //Mean of galactic WIMP velocity distribution [cm/s] const Double_t v_SunGal = 244. * 1.e5; //The Sun's rotational velocity around the galaxy [cm/s] const Double_t v_EarthSun = 15. * 1.e5; //The earth's rotational velocity around the Sun in the direction of the Sun's velocity around the galaxy [cm/s] const Double_t N_A = 6.022e23; //Avagadro's number [atoms/mol] const Double_t A_NatGe = 72.64, //Atomic number of natural Ge [g/mol] A_EnrGe = (0.86 * 76.) + (0.14 * 74.), //Atomic number of 86% enriched Ge [g/mol] A_NatAr = 39.948, //Atomic number of natural Ar [g/mol] A_NatNe = 20.1797; //Atomic number of natural Ne [g/mol] const Double_t M_p = 0.93827; //Proton mass [GeV] const Double_t c = 3.e10; //speed of light [cm/s] const Double_t GramsPerKilogram = 1000.; //Conversion between grams and kilograms const Double_t SecondsPerDay = 60. * 60. * 24.; //Conversion between seconds and days const Double_t keVperGeV = 1.e6; //Conversion between keV and GeV const Double_t Z_Ge = 32.; //Atomic number of germanium const Double_t Z_Ar = 18.; //Atomic number of argon const Double_t Z_Ne = 10.; //Atomic number of neon // WIMP halo/particle characteristics... const Double_t rho_D = 0.3; //Local DM halo density [GeV/cm^3] // A useful string... Char_t AString[1000], AnotherString[100]; // A flag for the sucessful return of a function... Int_t PlotReturn = -1; // Detector characteristics... const Double_t Z = Z_Ge; const Double_t A = A_EnrGe; //Select enriched germanium as a target material... const Double_t M_T = 0.9315 * A; //Set the mass of the target nucleus // Total rate in any month of the experiment... Double_t TotalRate = 0.; // Expected number of counts and number actually filled in each histogram... Double_t NCountsExpected = 0.; Int_t NCountsFilled = 0; // Nuclear recoil energy and quenched recoil energy... Double_t RecoilEnergy = 0., QuenchedEnergy = 0.; // A function to construct the simple recoil spectrum for a stationary earth, and infinite escape velocity... TF1 *DMSpecSimple(Double_t m_d, Double_t sigma_n, Double_t elo, Double_t ehi){ TF1 *DMSpec = new TF1("DMSpecSimple", "([0] / ([1] * [2])) * exp(-x / ([1] * [2]))", elo, ehi); DMSpec->SetTitle("Unquenched Nuclear Recoil Spectrum for Stationary Detector and Infinite Galactic Escape Velocity"); // Set up parameters... DMSpec->SetParName(0, "R_0"); DMSpec->FixParameter(0, (GramsPerKilogram * SecondsPerDay * 2. / TMath::Sqrt(TMath::Pi())) * (N_A / A) * (rho_D / m_d) * v_0 * TMath::Power(A * (M_T / M_p) * ((m_d + M_p) / (m_d + M_T)), 2.) * sigma_n); DMSpec->SetParName(1, "r"); DMSpec->FixParameter(1, 4 * (m_d * M_T) / ((m_d + M_T) * (m_d + M_T))); DMSpec->SetParName(2, "E_0"); DMSpec->FixParameter(2, 0.5 * m_d * keVperGeV * TMath::Power(v_0 / c, 2.)); DMSpec->GetXaxis()->SetTitle("Recoil Energy [keV]"); DMSpec->GetXaxis()->SetTitleOffset(0.9); DMSpec->GetYaxis()->SetTitle("Rate [dru]"); DMSpec->GetYaxis()->SetTitleOffset(1.3); //Return the spectrum function... return DMSpec; } // A function to construct the recoil spectrum including earth velocity and finite galactic escape velocity... TF1 *DMSpecRealVel(Double_t m_d, Double_t sigma_n, Double_t elo, Double_t ehi){ //Annually modulating earth velocity... TString *v_EFormula = new TString(100); v_EFormula->Append("([0] + ([1] * TMath::Sin(2. * TMath::Pi() * [9])))"); //Minimum velocity for energy x... TString *v_minFormula = new TString(100); v_minFormula->Append("([2] * TMath::Sqrt(x / ([3] * [4])))"); //Real rate formula... TString *SpecFormula = new TString(1000); SpecFormula->Append("([5] / ([3] * [4])) * ([6] / [7]) * (((0.25 * TMath::Sqrt(TMath::Pi()) * [2] / "); SpecFormula->Append(v_EFormula->Copy()); SpecFormula->Append(") * (TMath::Erf(("); SpecFormula->Append(v_minFormula->Copy()); SpecFormula->Append(" + "); SpecFormula->Append(v_EFormula->Copy()); SpecFormula->Append(") / [2]) - TMath::Erf(("); SpecFormula->Append(v_minFormula->Copy()); SpecFormula->Append(" - "); SpecFormula->Append(v_EFormula->Copy()); SpecFormula->Append(") / [2]))) - TMath::Exp(-1. * ([8]^2) / ([2]^2)))"); TF1 *DMSpec = new TF1("DMSpecRealVel",(const char*)SpecFormula->Copy(), elo, ehi); DMSpec->SetTitle("Unquenched Nuclear Recoil Spectrum for Real Earth Velocity Profile and Finite Galactic Escape Velocity"); // Set up parameters... DMSpec->SetParName(0, "v_SunGal"); DMSpec->FixParameter(0, v_SunGal); DMSpec->SetParName(1, "v_EarthSun"); DMSpec->FixParameter(1, v_EarthSun); DMSpec->SetParName(2, "v_0"); DMSpec->FixParameter(2, v_0); DMSpec->SetParName(3, "r"); DMSpec->FixParameter(3, 4 * (m_d * M_T) / ((m_d + M_T) * (m_d + M_T))); DMSpec->SetParName(4, "E_0"); DMSpec->FixParameter(4, 0.5 * m_d * keVperGeV * TMath::Power(v_0 / c, 2.)); DMSpec->SetParName(5, "R_0"); DMSpec->FixParameter(5, (GramsPerKilogram * SecondsPerDay * 2. / TMath::Sqrt(TMath::Pi())) * (N_A / A) * (rho_D / m_d) * v_0 * TMath::Power(A * (M_T / M_p) * ((m_d + M_p) / (m_d + M_T)), 2.) * sigma_n); DMSpec->SetParName(6, "k_0"); DMSpec->FixParameter(6, TMath::Power(TMath::Pi() * v_0 * v_0 ,1.5)); DMSpec->SetParName(7, "k_1"); DMSpec->FixParameter(7, TMath::Power(TMath::Pi() * v_0 * v_0 ,1.5) * (TMath::Erf(v_esc / v_0) - ((2. / TMath::Sqrt(TMath::Pi())) * (v_esc / v_0) * TMath::Exp(-1. * (v_esc * v_esc) / (v_0 * v_0))))); DMSpec->SetParName(8, "v_esc"); DMSpec->FixParameter(8, v_esc); DMSpec->SetParName(9, "Time"); DMSpec->SetParameter(9, 0.);//Time since March 2 [years] DMSpec->GetXaxis()->SetTitle("Recoil Energy [keV]"); DMSpec->GetXaxis()->SetTitleOffset(0.9); DMSpec->GetYaxis()->SetTitle("Rate [dru]"); DMSpec->GetYaxis()->SetTitleOffset(1.3); //Return the spectrum function... return DMSpec; } // A function to construct the recoil spectrum including earth velocity and finite galactic escape velocity... Double_t DMRateRealVel(Double_t m_d, Double_t sigma_n, Double_t time){ //Annually modulating earth velocity... TString *v_EFormula = new TString(100); v_EFormula->Append("([0] + ([1] * TMath::Sin(2. * TMath::Pi() * x)))"); //Real rate formula... TString *RateFormula = new TString(1000); RateFormula->Append("([2] * [3] / [4]) * ((0.5 * TMath::Sqrt(TMath::Pi()) * (("); RateFormula->Append(v_EFormula->Copy()); RateFormula->Append(" / [5]) + (0.5 * [5] / "); RateFormula->Append(v_EFormula->Copy()); RateFormula->Append(")) * TMath::Erf("); RateFormula->Append(v_EFormula->Copy()); RateFormula->Append(" / [5])) + (0.5 * TMath::Exp(-1. * ("); RateFormula->Append(v_EFormula->Copy()); RateFormula->Append("^2) / ([5]^2))) - (((([6]^2) / ([5]^2)) + (("); RateFormula->Append(v_EFormula->Copy()); RateFormula->Append("^2) / (3. * ([5]^2))) + 1.) * TMath::Exp(-1. * ([6]^2) / ([5]^2))))"); //Construct the TF1 object and set up the parameters... TF1 * Rate = new TF1("Rate",(const char*)RateFormula->Copy(),0., 1.); Rate->SetParName(0, "v_SunGal"); Rate->FixParameter(0, v_SunGal); Rate->SetParName(1, "v_EarthSun"); Rate->FixParameter(1, v_EarthSun); Rate->SetParName(2, "R_0"); Rate->FixParameter(2, (GramsPerKilogram * SecondsPerDay * 2. / TMath::Sqrt(TMath::Pi())) * (N_A / A) * (rho_D / m_d) * v_0 * TMath::Power(A * (M_T / M_p) * ((m_d + M_p) / (m_d + M_T)), 2.) * sigma_n); Rate->SetParName(3, "k_0"); Rate->FixParameter(3, TMath::Power(TMath::Pi() * v_0 * v_0 ,1.5)); Rate->SetParName(4, "k_1"); Rate->FixParameter(4, TMath::Power(TMath::Pi() * v_0 * v_0 ,1.5) * (TMath::Erf(v_esc / v_0) - ((2. / TMath::Sqrt(TMath::Pi())) * (v_esc / v_0) * TMath::Exp(-1. * (v_esc * v_esc) / (v_0 * v_0))))); Rate->SetParName(5, "v_0"); Rate->FixParameter(5, v_0); Rate->SetParName(6, "v_esc"); Rate->FixParameter(6, v_esc); //Return the rate at value time... return Rate->Eval(time); } // A function to make a nice plot of a TF1 object and write an image to the disk... Int_t MakePlot(TF1 *function, Double_t top_margin, Double_t bottom_margin, Double_t left_margin, Double_t right_margin, Int_t logx, Int_t logy, Char_t *filename){ // A nice canvas and pad on which to plot things... TCanvas *ACanvas = new TCanvas("ACanvas","Look! It's a Canvas!",27,0,1000,600); ACanvas->SetHighLightColor(kWhite); ACanvas->Range(-375,-0.60476,3375,2.43254); ACanvas->SetBorderSize(0); ACanvas->SetFrameFillColor(kWhite); ACanvas->SetFillColor(kWhite); TPad *APad = new TPad("APad", "It's a Pad", 0.,0.,1.,1., 10, 0,0); APad->SetTopMargin(top_margin); APad->SetBottomMargin(bottom_margin); APad->SetLeftMargin(left_margin); APad->SetRightMargin(right_margin); APad->SetGridx(0); APad->SetGridy(0); APad->SetLogx(logx); APad->SetLogy(logy); APad->Draw(); APad->cd(); // Find minumum and maximum of the function... Double_t FuncMax = 0.005, FuncMin = 1.e9, EvalPoint = 0.; for(Int_t i = 0; i < 100; i++){ EvalPoint = function->GetXmin() + (0.01 * i * (function->GetXmax() - function->GetXmin())); //std::cout << i << "\t" << EvalPoint << "\t" << function->Eval(EvalPoint) << std::endl; if(function->Eval(EvalPoint) > FuncMax){FuncMax = 1.1 * function->Eval(EvalPoint);} if(function->Eval(EvalPoint) < FuncMin){FuncMin = 0.9 * function->Eval(EvalPoint);} }// A histogram to setup the canvas... TH2D *CanvSetupHisto = new TH2D("CanvSetupHisto", function->GetTitle(), 1,function->GetXmin(),function->GetXmax(), 1,FuncMin,FuncMax); CanvSetupHisto->GetXaxis()->SetTitle(function->GetXaxis()->GetTitle()); CanvSetupHisto->GetXaxis()->SetTitleOffset(function->GetXaxis()->GetTitleOffset()); CanvSetupHisto->GetYaxis()->SetTitle(function->GetYaxis()->GetTitle()); CanvSetupHisto->GetYaxis()->SetTitleOffset(function->GetYaxis()->GetTitleOffset()); CanvSetupHisto->Draw(); // Draw the function... function->Draw("same"); ACanvas->Update(); // Write it to the disk... ACanvas->SaveAs(filename); // Delete "new" objects... delete CanvSetupHisto; delete APad; delete ACanvas; // Return a "0" and exit this function... return 0; } // A function to make a nice plot of a TH1D object and write an image to the disk... Int_t MakePlot(TH1D *histogram, Double_t top_margin, Double_t bottom_margin, Double_t left_margin, Double_t right_margin, Int_t logx, Int_t logy, Char_t *filename){ // A nice canvas and pad on which to plot things... TCanvas *ACanvas = new TCanvas("ACanvas","Look! It's a Canvas!",27,0,1000,600); ACanvas->SetHighLightColor(kWhite); ACanvas->Range(-375,-0.60476,3375,2.43254); ACanvas->SetBorderSize(0); ACanvas->SetFrameFillColor(kWhite); ACanvas->SetFillColor(kWhite); TPad *APad = new TPad("APad", "It's a Pad", 0.,0.,1.,1., 10, 0,0); APad->SetTopMargin(top_margin); APad->SetBottomMargin(bottom_margin); APad->SetLeftMargin(left_margin); APad->SetRightMargin(right_margin); APad->SetGridx(0); APad->SetGridy(0); APad->SetLogx(logx); APad->SetLogy(logy); APad->Draw(); APad->cd(); // Find minumum and maximum bin contents... Double_t BinMax = 0., BinMin = 1.e9; for(Int_t i = 0; i < histogram->GetNbinsX(); i++){ if(histogram->GetBinContent(i) > BinMax){BinMax = 1.1 * histogram->GetBinContent(i);} if(histogram->GetBinContent(i) < BinMin){BinMin = 0.9 * histogram->GetBinContent(i);} } // A histogram to setup the canvas... TH2D *CanvSetupHisto = new TH2D("CanvSetupHisto", histogram->GetTitle(), 1,histogram->GetXaxis()->GetXmin(),histogram->GetXaxis()->GetXmax(), 1,BinMin,BinMax); CanvSetupHisto->GetXaxis()->SetTitle(histogram->GetXaxis()->GetTitle()); CanvSetupHisto->GetXaxis()->SetTitleOffset(histogram->GetXaxis()->GetTitleOffset()); CanvSetupHisto->GetYaxis()->SetTitle(histogram->GetYaxis()->GetTitle()); CanvSetupHisto->GetYaxis()->SetTitleOffset(histogram->GetYaxis()->GetTitleOffset()); CanvSetupHisto->Draw(); // Draw the histogram... histogram->Draw("same"); ACanvas->Update(); // Write it to the disk... ACanvas->SaveAs(filename); // Delete "new" objects... delete CanvSetupHisto; delete APad; delete ACanvas; // Return a "0" and exit this function... return 0; } // A function to make a nice plot of a TGraph object and write an image to the disk... Int_t MakePlot(TGraph *graph, Double_t top_margin, Double_t bottom_margin, Double_t left_margin, Double_t right_margin, Int_t logx, Int_t logy, Char_t *filename){ // A nice canvas and pad on which to plot things... TCanvas *ACanvas = new TCanvas("ACanvas","Look! It's a Canvas!",27,0,1000,600); ACanvas->SetHighLightColor(kWhite); ACanvas->Range(-375,-0.60476,3375,2.43254); ACanvas->SetBorderSize(0); ACanvas->SetFrameFillColor(kWhite); ACanvas->SetFillColor(kWhite); TPad *APad = new TPad("APad", "It's a Pad", 0.,0.,1.,1., 10, 0,0); APad->SetTopMargin(top_margin); APad->SetBottomMargin(bottom_margin); APad->SetLeftMargin(left_margin); APad->SetRightMargin(right_margin); APad->SetGridx(0); APad->SetGridy(0); APad->SetLogx(logx); APad->SetLogy(logy); APad->Draw(); APad->cd(); // Find minumum and maximum of the function... Double_t GraphMax = -1.e9, GraphMin = 1.e9; for(Int_t i = 0; i < graph->GetN(); i++){ if(graph->GetY()[i] > GraphMax){GraphMax = 1.01 * graph->GetY()[i];} if(graph->GetY()[i] < GraphMin){GraphMin = 0.99 * graph->GetY()[i];} }// A histogram to setup the canvas... TH2D *CanvSetupHisto = new TH2D("CanvSetupHisto", graph->GetTitle(), 1,graph->GetX()[0],graph->GetX()[graph->GetN()-1], 1,GraphMin,GraphMax); CanvSetupHisto->GetXaxis()->SetTitle(graph->GetXaxis()->GetTitle()); CanvSetupHisto->GetXaxis()->SetTitleOffset(graph->GetXaxis()->GetTitleOffset()); CanvSetupHisto->GetYaxis()->SetTitle(graph->GetYaxis()->GetTitle()); CanvSetupHisto->GetYaxis()->SetTitleOffset(graph->GetYaxis()->GetTitleOffset()); CanvSetupHisto->Draw(); // Draw the function... graph->Draw("samelp"); ACanvas->Update(); //std::cin.get(); // Write it to the disk... ACanvas->SaveAs(filename); // Delete "new" objects... delete CanvSetupHisto; delete APad; delete ACanvas; // Return a "0" and exit this function... return 0; } // A function to incorporate energy quenching. This one does Lindhard-style ionization quenching... TF1 *QuenchEnergy_LindhardIonization(){ TString *kFormula = new TString(100); kFormula->Append("0.133 * TMath::Power([0], 0.66667) * TMath::Power([1], -0.5)"); TString *eFormula = new TString(100); eFormula->Append("16.2585 * x * TMath::Power([0], -2.3333)"); TString *gFormula = new TString(100); gFormula->Append("(3. * TMath::Power("); gFormula->Append(eFormula->Copy()); gFormula->Append(", 0.15)) + (0.7 * TMath::Power("); gFormula->Append(eFormula->Copy()); gFormula->Append(", 0.6)) + ("); gFormula->Append(eFormula->Copy()); gFormula->Append(")"); // Construct the quenching formula from g, k and epsilon... TString *QuenchingFunctionFormula = new TString(1000); QuenchingFunctionFormula->Append("("); QuenchingFunctionFormula->Append(kFormula->Copy()); QuenchingFunctionFormula->Append(" * ("); QuenchingFunctionFormula->Append(gFormula->Copy()); QuenchingFunctionFormula->Append(")) / (1 + ("); QuenchingFunctionFormula->Append(kFormula->Copy()); QuenchingFunctionFormula->Append(" * ("); QuenchingFunctionFormula->Append(gFormula->Copy()); QuenchingFunctionFormula->Append(")))"); // Set up the parameters in the quenching formula... TF1 *QuenchingFunction = new TF1("QuenchingFunction_LindhardIonization", (const char*)QuenchingFunctionFormula->Copy(), 0., 1000.); QuenchingFunction->SetParName(0, "Z"); QuenchingFunction->SetParameter(0, Z); QuenchingFunction->SetParName(1, "A"); QuenchingFunction->SetParameter(1, A); QuenchingFunction->SetTitle("Quenching Factor vs. Energy From Lindhard Theory"); QuenchingFunction->GetXaxis()->SetTitle("Recoil Energy [keV]"); QuenchingFunction->GetXaxis()->SetTitleOffset(0.9); QuenchingFunction->GetYaxis()->SetTitle("Quenching Factor"); QuenchingFunction->GetYaxis()->SetTitleOffset(0.7); return QuenchingFunction; }