/////////////////////////////////////////////////////////////////////////////////////////////////// // This script will generate the nuclear recoil spectrum from WIMP dark matter for an exposure // // specified at run time. We will start with the simplest recoil spectrum (stationary detector,// // infinite galactic escape velocity), and add in more physical effects as needed. // // // // 2008/07/03 - Original creation // // 2008/07/06 - Made the mass, cross section per neucleon and exposure arguments passed to the // // code rather than constants // // 2008/08/27 - Separated exposure out into active mass and run time for annual modulation // // analysis, created spectrum histograms for each month. Also added Lindhard-style// // quenching for ionization detectors // // 2008/09/14 - Adding in flat background model and canonical energy threshold for simple // // signal-to-noise ratio test // // // /////////////////////////////////////////////////////////////////////////////////////////////////// // My includes... #include "./DarkMatterSpectrum.h" // C++ includes... #include "iostream" #include "string" // Root includes... #include "TStyle.h" #include "TMath.h" #include "TRandom.h" #include "TF1.h" #include "TH1D.h" #include "TGraph.h" // Start the main fucntion // Units for this function: GeV cm^2 kg years DRU keV Boolean Double_t DarkMatterSpectrum(Double_t M_D, Double_t sigma_N, Double_t ActiveMass, Double_t RunTime, Double_t BGLevel, Double_t EnergyThreshold, Bool_t SavePlots){ // Number of months in the experiment... const Int_t nMonths = 1 + (Int_t)(RunTime * 12.); // An array to hold the number of counts in the spectrum, rates and number for each month... Double_t NCountsFilledArray[nMonths], Steps[nMonths], Rates[nMonths];; // Calculate exposure per month in kg * days... Double_t ExposurePerMonth = ActiveMass * 365.24 / 12.; // Seed random number generator with clock... gRandom->SetSeed(0); // Some stylistic things to make my plots beautiful... gStyle->SetOptStat(0); gStyle->SetTitleFillColor(kWhite); // Grab the function we're going to use to generate the recoil spectrum... TF1 *DMSpec = DMSpecRealVel(M_D, sigma_N, 0., 250.); // Plot the simple spectrum... if(SavePlots){ sprintf(AString, "Plots/%s_%.0f_%.1e.png", DMSpec->GetName(),M_D,sigma_N); PlotReturn = MakePlot(DMSpec, 0.08, 0.08, 0.095, 0.03, 0, 0, AString); } // Grab the appropriate quenching function... TF1 *QuenchEnergy = QuenchEnergy_LindhardIonization(); if(SavePlots){ sprintf(AString, "Plots/%s_%.0f_%.1e.png", QuenchEnergy->GetName(),M_D,sigma_N); PlotReturn = MakePlot(QuenchEnergy, 0.08, 0.08, 0.06, 0.03, 0, 0, AString); } // Number of experiments to perform... const Int_t NExperiments = 10; // Number of experiments where the DM spectrum is above background at threshold... Int_t DMAboveBG = 0; // Loop over the NExperiments... for(Int_t iExperiment = 0; iExperiment < NExperiments; iExperiment++){ TH1D *DMSpecHisto[nMonths]; for(Int_t iMonth = 0; iMonth < nMonths; iMonth++){ sprintf(AString, "Energy Spectrum Drawn From %s for Month %d", DMSpec->GetName(), iMonth); sprintf(AnotherString, "DMSpecHisto%d", iMonth); DMSpecHisto[iMonth] = new TH1D(AnotherString, AString, (Int_t)(DMSpec->GetXmax() - DMSpec->GetXmin()), DMSpec->GetXmin(), DMSpec->GetXmax()); DMSpecHisto[iMonth]->GetXaxis()->SetTitle("Recoil Energy [keV]"); DMSpecHisto[iMonth]->GetXaxis()->SetTitleOffset(0.9); DMSpecHisto[iMonth]->GetYaxis()->SetTitle("Rate [dru]"); DMSpecHisto[iMonth]->GetYaxis()->SetTitleOffset(1.3); // Report total rate... if(!strcmp(DMSpec->GetName(), "DMSpecSimple")){ TotalRate = DMSpec->GetParameter(0); }else{ TotalRate = DMRateRealVel(M_D, sigma_N, (Double_t)iMonth); } //std::cout << "Total rate = " << TotalRate << " counts/kg/day.\nFor a one-month exposure of " << ExposurePerMonth << " kg x days, we expect "; //Calculate the expected number of counts... NCountsExpected = TotalRate * ExposurePerMonth; //std::cout << NCountsExpected << " counts in the spectrum. We will fill the spectrum with "; NCountsFilled = gRandom->Poisson(NCountsExpected);//std::cout << NCountsFilled << " events." << std::endl; NCountsFilledArray[iMonth] = (Double_t)NCountsFilled; Steps[iMonth] = (1. / 12.) * (Double_t)iMonth; Rates[iMonth] = DMRateRealVel(M_D, sigma_N, Steps[iMonth]); //Fill the spectra... if(!strcmp(DMSpec->GetName(), "DMSpecSimple")){ DMSpecHisto[iMonth]->FillRandom(DMSpec->GetName(), NCountsFilled); }else{ DMSpec->SetParameter(9, (1. / 12.) * (Double_t)iMonth); //Set the time of year... if(SavePlots){ sprintf(AString, "Plots/SpecMod/%s_%d_%.0f_%.1e_Exp%d.png",DMSpec->GetName(),iMonth,M_D,sigma_N,iExperiment); PlotReturn = MakePlot(DMSpec, 0.07, 0.08, 0.095, 0.03, 0, 0, AString); } for(Int_t j = 0; j < NCountsFilled; j++){ RecoilEnergy = DMSpec->GetRandom(); QuenchedEnergy = RecoilEnergy * QuenchEnergy->Eval(RecoilEnergy); DMSpecHisto[iMonth]->Fill(QuenchedEnergy); } } DMSpecHisto[iMonth]->Scale(1. / ExposurePerMonth); if(SavePlots){ sprintf(AString, "Plots/SpecMod/%s_%.0f_%.1e_Exp%d.png", DMSpecHisto[iMonth]->GetName(),M_D,sigma_N,iExperiment); PlotReturn = MakePlot(DMSpecHisto[iMonth], 0.07, 0.08, 0.095, 0.03, 0, 0, AString); } } // Constructing the sum spectrum for the entire run... sprintf(AString, "Energy Spectrum Drawn From %s Over %.1f Years", DMSpec->GetName(), RunTime); TH1D *DMSpecSumHisto = new TH1D("DMSpecSumHisto", AString, (Int_t)(DMSpec->GetXmax() - DMSpec->GetXmin()), DMSpec->GetXmin(), DMSpec->GetXmax()); DMSpecSumHisto->GetXaxis()->SetTitle("Recoil Energy [keV]"); DMSpecSumHisto->GetXaxis()->SetTitleOffset(0.9); DMSpecSumHisto->GetYaxis()->SetTitle("Rate [dru]"); DMSpecSumHisto->GetYaxis()->SetTitleOffset(1.3); for(Int_t iMonth = 0; iMonth < nMonths; iMonth++){DMSpecSumHisto->Add(DMSpecHisto[iMonth], 1. / (Double_t)nMonths);} // Check if the sum spectrum is above background at threshold... if(DMSpecSumHisto->GetBinContent(DMSpecSumHisto->FindBin(EnergyThreshold)) > BGLevel){DMAboveBG++;} if(SavePlots){ sprintf(AString, "Plots/%s_%.0f_%.1e_Exp%d.png", DMSpecSumHisto->GetName(),M_D,sigma_N,iExperiment); PlotReturn = MakePlot(DMSpecSumHisto, 0.07, 0.08, 0.095, 0.03, 0, 0, AString); } // Make a plot of the total rate as a function of time... TGraph *RateGraph = new TGraph(nMonths, Steps, Rates); RateGraph->SetName("RateGraph"); sprintf(AString, "Total Rate Over the Course of the %.1f-year Run Time", RunTime); RateGraph->SetTitle(AString); RateGraph->SetMarkerStyle(kFullCircle); RateGraph->SetMarkerSize(1.2); RateGraph->GetXaxis()->SetTitle("Time Since March 2 [y]"); RateGraph->GetXaxis()->SetTitleOffset(0.9); RateGraph->GetYaxis()->SetTitle("Total Rate [tru]"); RateGraph->GetYaxis()->SetTitleOffset(1.1); if(SavePlots){ sprintf(AString, "Plots/%s_%.0f_%.1e_Exp%d.png", RateGraph->GetName(),M_D,sigma_N,iExperiment); PlotReturn = MakePlot(RateGraph, 0.08, 0.08, 0.08, 0.01, 0, 0, AString); } // Make a plot of the number of counts in the spectrum each month... TGraph *NCountsFilledGraph = new TGraph(nMonths, Steps, NCountsFilledArray); NCountsFilledGraph->SetName("NCountsFilledGraph"); sprintf(AString, "Number of Counts per Month Over the %.1f-year Run Time", RunTime); NCountsFilledGraph->SetTitle(AString); NCountsFilledGraph->SetMarkerStyle(kFullCircle); NCountsFilledGraph->SetMarkerSize(1.2); NCountsFilledGraph->GetXaxis()->SetTitle("Time Since March 2 [y]"); NCountsFilledGraph->GetYaxis()->SetTitle("Number of Counts per Month"); if(SavePlots){ sprintf(AString, "Plots/%s_%.0f_%.1e_Exp%d.png", NCountsFilledGraph->GetName(),M_D,sigma_N,iExperiment); PlotReturn = MakePlot(NCountsFilledGraph, 0.08, 0.08, 0.08, 0.01, 0, 0,AString); } // Delete "new" objects created in the experiment loop... delete NCountsFilledGraph; delete RateGraph; delete DMSpecSumHisto; for(Int_t iMonth = 0; iMonth < nMonths; iMonth++){delete DMSpecHisto[iMonth];} } // Return the fraction of experiments where the dark matter is visible above background and we're done... return ((Double_t)DMAboveBG / (Double_t)NExperiments); }