From 420a468f5c38bf66fa92fda900b132f4e3789b3c Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Tue, 5 Mar 2024 09:41:50 +0100 Subject: [PATCH 1/5] REST_Axion_PlotResonances.C adding Ngamma plotting --- macros/REST_Axion_PlotResonances.C | 253 ++++++++++++++++++++++------- 1 file changed, 196 insertions(+), 57 deletions(-) diff --git a/macros/REST_Axion_PlotResonances.C b/macros/REST_Axion_PlotResonances.C index 14bd78c2..8e960731 100644 --- a/macros/REST_Axion_PlotResonances.C +++ b/macros/REST_Axion_PlotResonances.C @@ -17,120 +17,259 @@ //*** - *Bmag* = 2.5: Magnetic field (in T). //*** - *Lcoh* = 10000: Coherence length (in mm). //*** - *gasName* = "He": Gas name. -//*** - *vacuum* = true: If true, the vacuum probability is added to the sum of all the probabilities. +//*** - *cutoff* = It defines a cut-off for the number of resonances to be drawn. //*** - *n_ma* = 10000: Number of points to be plotted. //*** //*** The generated plots are the results from `TRestAxionField::GetMassDensityScanning`, //*** `TRestAxionField::GammaTransmissionFWHM` and `TRestAxionBufferGas::GetMassDensity`. //*** //*** -------------- -//*** Usage: restManager PlotResonances [ma_max=0.1] [ma_min=0] [Ea=4.2] [Bmag=2.5] [Lcoh=10000] [gasName=He] -//*** [vacuum=true] [n_ma=10000] +//*** Usage: restManager PlotResonances [options] [ma_max=0.1] [ma_min=0] [Ea=4.2] [Bmag=2.5] [Lcoh=10000] [gasName=He] +//*** [cutoff=5] [n_ma=10000] //*** //*** Being all of them optional arguments. //*** -------------- //*** Author: Fran Candón //******************************************************************************************************* +// +const Double_t fProbMax = 3e-18; +const Double_t fNGammaMax = 6000; + +const TVector2 fEnergyRange = TVector2(0.5,10); +const Double_t fExposureTime = 30 * 3600. * 18; // s +const Double_t fArea = TMath::Pi() * 35 * 35; // cm2 +const Double_t fReBinning = 100; // Transforms 0.001 keV into 0.1keV + +int REST_Axion_PlotResonances(std::string optionString = "", double ma_max = 0.1, double ma_min = 0, double Ea = 4.2, double Bmag = 2., double Lcoh = 10000, std::string gasName = "He", int cutoff = 5, int n_ma = 1000) { + + gStyle->SetPadRightMargin(0.13); + TCanvas* c1 = new TCanvas("c1", "c1", 800, 600); + + std::vector options = TRestTools::GetOptions(optionString); + + bool title = std::find(options.begin(), options.end(), "title") != options.end(); + bool vacuum = std::find(options.begin(), options.end(), "vacuum") != options.end(); + bool drawLine = std::find(options.begin(), options.end(), "drawLine") != options.end(); + bool sumProb = std::find(options.begin(), options.end(), "sumProb") != options.end(); + bool legend = std::find(options.begin(), options.end(), "legend") != options.end(); + bool nGamma = std::find(options.begin(), options.end(), "nGamma") != options.end(); -int REST_Axion_PlotResonances(double ma_max = 0.1, double ma_min = 0, double Ea = 4.2, double Bmag = 2.5, - double Lcoh = 10000, std::string gasName = "He", Bool_t vacuum = true, - int n_ma = 10000) { TRestAxionField* ax = new TRestAxionField(); ax->SetMagneticField(Bmag); ax->SetCoherenceLength(Lcoh); ax->SetAxionEnergy(Ea); - vector> pair = ax->GetMassDensityScanning(gasName, ma_max, 20); - std::vector m_a; - std::vector sum_prob; + /// Could be stored as a data member of TRestAxionPlotResonances for further access + std::vector> Psettings = ax->GetMassDensityScanning(gasName, ma_max, 20); + + if( Psettings.size() > cutoff ) Psettings.resize(cutoff); - // Creates the vector of axion masses + // Creates the vector of axion masses (This could be a data member of TRestAxionPlotResonances). + std::vector m_a; double ma_step = (ma_max - ma_min) / n_ma; + m_a.push_back( ma_min ); // The first mass drawn will get zero probability for the curve filling. + // So we repeat the value + + std::vector sum_prob; + sum_prob.push_back(0); for (int i = 0; i < n_ma; i++) { m_a.push_back(ma_min + i * ma_step); sum_prob.push_back(0); } - TCanvas* c1 = new TCanvas("c1", "c1", 800, 600); - std::vector grp; + // Computes the Vacuum probability + ax->AssignBufferGas(nullptr); + std ::vector prob_vac; + + // Probability is not zero, but we introduce an artifact (virtual point) to make proper TGraph filling + prob_vac.push_back(0); + for (size_t j = 1; j < m_a.size(); j++) { + prob_vac.push_back(ax->GammaTransmissionProbability(m_a[j])); + } + // Computes the sum of all the probabilities (Adding vacuum) + if (vacuum == true) { + for (size_t i = 0; i < m_a.size(); i++) { + sum_prob[i] += prob_vac[i]; + } + } + + std::vector grp; TRestAxionBufferGas* gas = new TRestAxionBufferGas(); - for (const auto& p : pair) { + for (const auto& p : Psettings) { // Creates the gas and the axion field gas->SetGasDensity(gasName, p.second); ax->AssignBufferGas(gas); // Obtain the probability for each axion mass std::vector prob; - for (int j = 0; j < n_ma; j++) { - prob.push_back(ax->GammaTransmissionProbability(m_a[j])); + // Probability is not zero, but we introduce an artifact to make proper curve filling + prob.push_back(0); + for (size_t j = 1; j < m_a.size(); j++) { + prob.push_back(ax->GammaTransmissionProbability(m_a[j])); + sum_prob[j] += prob[j]; } - TGraph* gr = new TGraph(n_ma, &m_a[0], &prob[0]); + TGraph* gr = new TGraph(m_a.size(), &m_a[0], &prob[0]); grp.push_back(gr); } - // Computes the Vacuum probability - TRestAxionField* ax_vac = new TRestAxionField(); - std ::vector prob_vac; - for (int j = 0; j < n_ma; j++) { - prob_vac.push_back(ax_vac->GammaTransmissionProbability(m_a[j])); - } - - // Computes the sum of all the probabilities - if (vacuum == true) { - for (int i = 0; i < n_ma; i++) { - sum_prob[i] += prob_vac[i]; - } - } ///// PLOTS ///// // Plot the density scan grp[0]->SetLineColor(kBlue - 3); - grp[0]->SetTitle("Transmission probability as a function of the axion mass"); + if( title ) grp[0]->SetTitle("Transmission probability as a function of the axion mass"); + else grp[0]->SetTitle(""); grp[0]->GetXaxis()->SetTitle("m_{a} [eV]"); - grp[0]->GetYaxis()->SetTitle("P_{ag}"); - grp[0]->GetXaxis()->SetLimits(ma_min, ma_max); - double ylim = 3.5e-18; - grp[0]->GetYaxis()->SetRangeUser(0, ylim); - grp[0]->Draw("AL"); + grp[0]->GetYaxis()->SetTitle("P_{a#gamma}"); + grp[0]->GetXaxis()->SetLimits(ma_min+0.00001, ma_max); + grp[0]->GetXaxis()->SetLabelSize(0.04); + grp[0]->GetXaxis()->SetLabelFont(42); + grp[0]->GetXaxis()->SetTitleSize(0.04); + grp[0]->GetXaxis()->SetTitleFont(42); + grp[0]->GetYaxis()->SetRangeUser(0, fProbMax); + grp[0]->GetYaxis()->SetLabelSize(0.04); + grp[0]->GetYaxis()->SetLabelFont(42); + grp[0]->GetYaxis()->SetTitleSize(0.04); + grp[0]->GetYaxis()->SetTitleFont(42); + + grp[0]->SetLineColor(kBlack); + grp[0]->SetFillColorAlpha(kRed,0.5); + grp[0]->SetLineWidth(2); + grp[0]->Draw("AFL"); for (const auto g : grp) { - g->SetLineColor(kBlue - 3); - g->Draw("SAME"); + g->SetLineColor(kBlack); + g->SetLineWidth(2); + g->SetFillColorAlpha(kRed,0.15); + g->Draw("FL SAME"); } // PLot of the sum of all the probabilities - TGraph* sum = new TGraph(n_ma, &m_a[0], &sum_prob[0]); - sum->SetLineColor(kBlue + 3); - sum->SetLineWidth(3); - sum->Draw("SAME"); + TGraph *sum; + if( sumProb ) + { + sum = new TGraph(m_a.size(), &m_a[0], &sum_prob[0]); + sum->SetLineColor(kBlue + 3); + sum->SetLineWidth(3); + sum->Draw("SAME"); + } // Plot of the vacuum probability - TGraph* vac = new TGraph(n_ma, &m_a[0], &prob_vac[0]); - vac->SetLineColor(kCyan + 1); + TGraph* vac = new TGraph(m_a.size(), &m_a[0], &prob_vac[0]); + vac->SetLineColor(kBlack); vac->SetLineWidth(2); - vac->Draw("SAME"); + vac->SetFillColorAlpha(kBlue,0.25); + vac->Draw("FL SAME"); + + TGraph *gama; + if( nGamma ) // The following code could be a method of a future TRestAxionPlotResonances::DrawNGamma() + { + std ::vector NGamma; + + TRestAxionSolarQCDFlux* sFlux = new TRestAxionSolarQCDFlux("fluxes.rml", "LennertHoofPrimakoff"); + sFlux->Initialize(); + TH1F* spec = sFlux->GetEnergySpectrum(); + + /// The original gets 1eV binning to describe monocromatic lines. + /// We do not need such a high resolution + TH1F *rebinned = dynamic_cast(spec->Rebin(fReBinning,"rebinned")); + /// We need to renormalize since the bin contents are already renormalized to cm-2 s-1 keV-1 + rebinned->Scale(1./fReBinning); + Double_t deltaE = rebinned->GetBinCenter(2) - rebinned->GetBinCenter(1); + + if(!rebinned) { + std::cout << "Energy spectrum is nullptr!" << std::endl; + } + else + { + Double_t maxNGamma = 0; + // Obtain the Ngamma for each axion mass + for (size_t j = 0; j < m_a.size(); j++) { + std::cout << "Producing Ngamma for mass : " << m_a[j] << std::endl; + + Double_t ng = 0; + + // Vacuum + ax->AssignBufferGas(nullptr); + for( int n = 1; n <= rebinned->GetNbinsX(); n++ ) + { + Double_t energy = rebinned->GetBinCenter(n); + if ( energy < fEnergyRange.X() || energy > fEnergyRange.Y() ) + continue; + + ax->SetAxionEnergy( rebinned->GetBinCenter(n) ); + ng += rebinned->GetBinContent(n) * ax->GammaTransmissionProbability(m_a[j]); + } + + /// Gas phase + for (const auto& p : Psettings) { + // Creates the gas and assignss it to the axion field + gas->SetGasDensity(gasName, p.second); + ax->AssignBufferGas(gas); + + for( int n = 1; n <= rebinned->GetNbinsX(); n++ ) + { + Double_t energy = rebinned->GetBinCenter(n); + if ( energy < fEnergyRange.X() || energy > fEnergyRange.Y() ) + continue; + + ax->SetAxionEnergy( rebinned->GetBinCenter(n) ); + ng += rebinned->GetBinContent(n) * ax->GammaTransmissionProbability(m_a[j]); + } + } + ng = ng * deltaE * fArea * fExposureTime; + if( ng > maxNGamma ) maxNGamma = ng; + NGamma.push_back( ng ); + } + + gama = new TGraph(m_a.size(), &m_a[0], &NGamma[0]); + gama->SetLineColor(kRed + 3); + gama->SetLineWidth(3); + gama->GetYaxis()->SetRangeUser(0, fNGammaMax); + gama->Scale(fProbMax/fNGammaMax); + gama->Draw("SAME"); + + std::cout << "Max: " << maxNGamma << std::endl; + + /// 0.001 is to remove the first tick, which is zero, and overlaps with 0.1eV + TGaxis *A1 = new TGaxis(0.1,0,0.1,fProbMax,0.001,fNGammaMax,510,"+L"); + A1->SetTitle("N_{#gamma}"); + A1->SetTitleOffset(1.5); + A1->SetLabelSize(0.04); + A1->SetLabelFont(42); + A1->SetTextFont(42); + A1->Draw(); + } + } // Plot of the vertical line - TLine* verticalLine = new TLine(pair[0].first, c1->GetUymin(), pair[0].first, ylim); - verticalLine->SetLineColor(kGreen - 3); - verticalLine->SetLineWidth(2); - verticalLine->Draw("same"); + TLine* verticalLine; + if( drawLine ) + { + verticalLine = new TLine(Psettings[0].first, c1->GetUymin(), Psettings[0].first, fProbMax); + verticalLine->SetLineColor(kGreen - 3); + verticalLine->SetLineWidth(2); + verticalLine->Draw("same"); + } // Plot of the legend - TLegend* legend = new TLegend(0.9, 0.7, 0.48, 0.9); - legend->AddEntry(grp[0], "P_{ag} for each mass", "l"); - legend->AddEntry(vac, "P_{ag} for vacuum", "l"); - legend->AddEntry(sum, "Sum of all the probs.", "l"); - legend->AddEntry(verticalLine, "m_{a} where P_{ag}^{vac} = max(P_{ag}^{vac}/2) ", "l"); - legend->Draw("same"); - c1->Draw(); - - c1->Print("/tmp/resonances.png"); + if( legend ) + { + TLegend* legend = new TLegend(0.87, 0.7, 0.47, 0.9); + legend->AddEntry(grp[0], "P_{a#gamma} for each P-step", "lf"); + legend->AddEntry(vac, "P_{a#gamma} for vacuum", "lf"); + if( sumProb ) legend->AddEntry(sum, "#Sigma P_{a#gamma}", "l"); + if( nGamma ) legend->AddEntry(gama, "N_{#gamma}", "l"); + if( drawLine ) legend->AddEntry(verticalLine, "m_{a} where P_{a#gamma}^{vac} = max(P_{a#gamma}^{vac}/2) ", "l"); + legend->Draw("same"); + c1->Draw(); + } + + c1->Print("/tmp/resonances.pdf"); return 0; } From 8a3d72cd51419d89a5047000ca4e5f497c69f0cf Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Tue, 5 Mar 2024 10:28:45 +0100 Subject: [PATCH 2/5] macros/REST_Axion_PlotResonances.C adding labels --- macros/REST_Axion_PlotResonances.C | 39 ++++++++++++++++++++++++++---- 1 file changed, 34 insertions(+), 5 deletions(-) diff --git a/macros/REST_Axion_PlotResonances.C b/macros/REST_Axion_PlotResonances.C index 8e960731..aeaba846 100644 --- a/macros/REST_Axion_PlotResonances.C +++ b/macros/REST_Axion_PlotResonances.C @@ -43,6 +43,7 @@ const Double_t fReBinning = 100; // Transforms 0.001 keV into 0.1keV int REST_Axion_PlotResonances(std::string optionString = "", double ma_max = 0.1, double ma_min = 0, double Ea = 4.2, double Bmag = 2., double Lcoh = 10000, std::string gasName = "He", int cutoff = 5, int n_ma = 1000) { gStyle->SetPadRightMargin(0.13); + gStyle->SetPadBottomMargin(0.13); TCanvas* c1 = new TCanvas("c1", "c1", 800, 600); std::vector options = TRestTools::GetOptions(optionString); @@ -53,6 +54,7 @@ int REST_Axion_PlotResonances(std::string optionString = "", double ma_max = 0.1 bool sumProb = std::find(options.begin(), options.end(), "sumProb") != options.end(); bool legend = std::find(options.begin(), options.end(), "legend") != options.end(); bool nGamma = std::find(options.begin(), options.end(), "nGamma") != options.end(); + bool labels = std::find(options.begin(), options.end(), "labels") != options.end(); TRestAxionField* ax = new TRestAxionField(); ax->SetMagneticField(Bmag); @@ -97,6 +99,7 @@ int REST_Axion_PlotResonances(std::string optionString = "", double ma_max = 0.1 std::vector grp; TRestAxionBufferGas* gas = new TRestAxionBufferGas(); + std::vector> labelPositions; for (const auto& p : Psettings) { // Creates the gas and the axion field gas->SetGasDensity(gasName, p.second); @@ -106,11 +109,21 @@ int REST_Axion_PlotResonances(std::string optionString = "", double ma_max = 0.1 std::vector prob; // Probability is not zero, but we introduce an artifact to make proper curve filling prob.push_back(0); + Double_t xLabel = 0; + Double_t yLabel = 0; for (size_t j = 1; j < m_a.size(); j++) { - prob.push_back(ax->GammaTransmissionProbability(m_a[j])); + Double_t probV = ax->GammaTransmissionProbability(m_a[j]); + if( probV > yLabel ) + { + yLabel = probV; + xLabel = m_a[j]; + } + prob.push_back(probV); sum_prob[j] += prob[j]; } + std::pair coords(xLabel-0.002, 1.05*yLabel); + labelPositions.push_back(coords); TGraph* gr = new TGraph(m_a.size(), &m_a[0], &prob[0]); grp.push_back(gr); @@ -126,26 +139,41 @@ int REST_Axion_PlotResonances(std::string optionString = "", double ma_max = 0.1 grp[0]->GetXaxis()->SetTitle("m_{a} [eV]"); grp[0]->GetYaxis()->SetTitle("P_{a#gamma}"); grp[0]->GetXaxis()->SetLimits(ma_min+0.00001, ma_max); + grp[0]->GetXaxis()->SetLabelSize(0.04); grp[0]->GetXaxis()->SetLabelFont(42); - grp[0]->GetXaxis()->SetTitleSize(0.04); + grp[0]->GetXaxis()->SetTitleSize(0.045); grp[0]->GetXaxis()->SetTitleFont(42); + grp[0]->GetYaxis()->SetRangeUser(0, fProbMax); + grp[0]->GetYaxis()->SetLabelSize(0.04); grp[0]->GetYaxis()->SetLabelFont(42); - grp[0]->GetYaxis()->SetTitleSize(0.04); + grp[0]->GetYaxis()->SetTitleSize(0.05); grp[0]->GetYaxis()->SetTitleFont(42); + grp[0]->GetYaxis()->SetTitleOffset(0.9); grp[0]->SetLineColor(kBlack); grp[0]->SetFillColorAlpha(kRed,0.5); - grp[0]->SetLineWidth(2); + grp[0]->SetLineWidth(1); grp[0]->Draw("AFL"); + int n = 0; for (const auto g : grp) { g->SetLineColor(kBlack); g->SetLineWidth(2); g->SetFillColorAlpha(kRed,0.15); g->Draw("FL SAME"); + + if( labels) + { + std::string label = "P" + std::to_string(n+1); + TLatex* textLatex = new TLatex(labelPositions[n].first, labelPositions[n].second, label.c_str()); + textLatex->SetTextColor(1); + textLatex->SetTextSize(0.02); + textLatex->Draw("same"); + } + n++; } // PLot of the sum of all the probabilities @@ -238,7 +266,8 @@ int REST_Axion_PlotResonances(std::string optionString = "", double ma_max = 0.1 /// 0.001 is to remove the first tick, which is zero, and overlaps with 0.1eV TGaxis *A1 = new TGaxis(0.1,0,0.1,fProbMax,0.001,fNGammaMax,510,"+L"); A1->SetTitle("N_{#gamma}"); - A1->SetTitleOffset(1.5); + A1->SetTitleOffset(1.3); + A1->SetTitleSize(0.045); A1->SetLabelSize(0.04); A1->SetLabelFont(42); A1->SetTextFont(42); From 109cb4c8686fd59078b4f6d9f0dba445fb3b582d Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Tue, 5 Mar 2024 09:29:53 +0000 Subject: [PATCH 3/5] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- macros/REST_Axion_PlotResonances.C | 337 ++++++++++++++--------------- 1 file changed, 165 insertions(+), 172 deletions(-) diff --git a/macros/REST_Axion_PlotResonances.C b/macros/REST_Axion_PlotResonances.C index aeaba846..dc0fb093 100644 --- a/macros/REST_Axion_PlotResonances.C +++ b/macros/REST_Axion_PlotResonances.C @@ -24,7 +24,8 @@ //*** `TRestAxionField::GammaTransmissionFWHM` and `TRestAxionBufferGas::GetMassDensity`. //*** //*** -------------- -//*** Usage: restManager PlotResonances [options] [ma_max=0.1] [ma_min=0] [Ea=4.2] [Bmag=2.5] [Lcoh=10000] [gasName=He] +//*** Usage: restManager PlotResonances [options] [ma_max=0.1] [ma_min=0] [Ea=4.2] [Bmag=2.5] [Lcoh=10000] +//[gasName=He] //*** [cutoff=5] [n_ma=10000] //*** //*** Being all of them optional arguments. @@ -35,58 +36,59 @@ const Double_t fProbMax = 3e-18; const Double_t fNGammaMax = 6000; -const TVector2 fEnergyRange = TVector2(0.5,10); -const Double_t fExposureTime = 30 * 3600. * 18; // s -const Double_t fArea = TMath::Pi() * 35 * 35; // cm2 -const Double_t fReBinning = 100; // Transforms 0.001 keV into 0.1keV +const TVector2 fEnergyRange = TVector2(0.5, 10); +const Double_t fExposureTime = 30 * 3600. * 18; // s +const Double_t fArea = TMath::Pi() * 35 * 35; // cm2 +const Double_t fReBinning = 100; // Transforms 0.001 keV into 0.1keV -int REST_Axion_PlotResonances(std::string optionString = "", double ma_max = 0.1, double ma_min = 0, double Ea = 4.2, double Bmag = 2., double Lcoh = 10000, std::string gasName = "He", int cutoff = 5, int n_ma = 1000) { - - gStyle->SetPadRightMargin(0.13); - gStyle->SetPadBottomMargin(0.13); +int REST_Axion_PlotResonances(std::string optionString = "", double ma_max = 0.1, double ma_min = 0, + double Ea = 4.2, double Bmag = 2., double Lcoh = 10000, + std::string gasName = "He", int cutoff = 5, int n_ma = 1000) { + gStyle->SetPadRightMargin(0.13); + gStyle->SetPadBottomMargin(0.13); TCanvas* c1 = new TCanvas("c1", "c1", 800, 600); - std::vector options = TRestTools::GetOptions(optionString); + std::vector options = TRestTools::GetOptions(optionString); - bool title = std::find(options.begin(), options.end(), "title") != options.end(); - bool vacuum = std::find(options.begin(), options.end(), "vacuum") != options.end(); - bool drawLine = std::find(options.begin(), options.end(), "drawLine") != options.end(); - bool sumProb = std::find(options.begin(), options.end(), "sumProb") != options.end(); - bool legend = std::find(options.begin(), options.end(), "legend") != options.end(); - bool nGamma = std::find(options.begin(), options.end(), "nGamma") != options.end(); - bool labels = std::find(options.begin(), options.end(), "labels") != options.end(); + bool title = std::find(options.begin(), options.end(), "title") != options.end(); + bool vacuum = std::find(options.begin(), options.end(), "vacuum") != options.end(); + bool drawLine = std::find(options.begin(), options.end(), "drawLine") != options.end(); + bool sumProb = std::find(options.begin(), options.end(), "sumProb") != options.end(); + bool legend = std::find(options.begin(), options.end(), "legend") != options.end(); + bool nGamma = std::find(options.begin(), options.end(), "nGamma") != options.end(); + bool labels = std::find(options.begin(), options.end(), "labels") != options.end(); TRestAxionField* ax = new TRestAxionField(); ax->SetMagneticField(Bmag); ax->SetCoherenceLength(Lcoh); ax->SetAxionEnergy(Ea); - /// Could be stored as a data member of TRestAxionPlotResonances for further access - std::vector> Psettings = ax->GetMassDensityScanning(gasName, ma_max, 20); + /// Could be stored as a data member of TRestAxionPlotResonances for further access + std::vector> Psettings = ax->GetMassDensityScanning(gasName, ma_max, 20); - if( Psettings.size() > cutoff ) Psettings.resize(cutoff); + if (Psettings.size() > cutoff) Psettings.resize(cutoff); // Creates the vector of axion masses (This could be a data member of TRestAxionPlotResonances). std::vector m_a; double ma_step = (ma_max - ma_min) / n_ma; - m_a.push_back( ma_min ); // The first mass drawn will get zero probability for the curve filling. - // So we repeat the value + m_a.push_back(ma_min); // The first mass drawn will get zero probability for the curve filling. + // So we repeat the value std::vector sum_prob; - sum_prob.push_back(0); + sum_prob.push_back(0); for (int i = 0; i < n_ma; i++) { m_a.push_back(ma_min + i * ma_step); sum_prob.push_back(0); } // Computes the Vacuum probability - ax->AssignBufferGas(nullptr); + ax->AssignBufferGas(nullptr); std ::vector prob_vac; - // Probability is not zero, but we introduce an artifact (virtual point) to make proper TGraph filling - prob_vac.push_back(0); + // Probability is not zero, but we introduce an artifact (virtual point) to make proper TGraph filling + prob_vac.push_back(0); for (size_t j = 1; j < m_a.size(); j++) { - prob_vac.push_back(ax->GammaTransmissionProbability(m_a[j])); + prob_vac.push_back(ax->GammaTransmissionProbability(m_a[j])); } // Computes the sum of all the probabilities (Adding vacuum) @@ -99,7 +101,7 @@ int REST_Axion_PlotResonances(std::string optionString = "", double ma_max = 0.1 std::vector grp; TRestAxionBufferGas* gas = new TRestAxionBufferGas(); - std::vector> labelPositions; + std::vector> labelPositions; for (const auto& p : Psettings) { // Creates the gas and the axion field gas->SetGasDensity(gasName, p.second); @@ -107,196 +109,187 @@ int REST_Axion_PlotResonances(std::string optionString = "", double ma_max = 0.1 // Obtain the probability for each axion mass std::vector prob; - // Probability is not zero, but we introduce an artifact to make proper curve filling - prob.push_back(0); - Double_t xLabel = 0; - Double_t yLabel = 0; + // Probability is not zero, but we introduce an artifact to make proper curve filling + prob.push_back(0); + Double_t xLabel = 0; + Double_t yLabel = 0; for (size_t j = 1; j < m_a.size(); j++) { - Double_t probV = ax->GammaTransmissionProbability(m_a[j]); - if( probV > yLabel ) - { - yLabel = probV; - xLabel = m_a[j]; - } - prob.push_back(probV); + Double_t probV = ax->GammaTransmissionProbability(m_a[j]); + if (probV > yLabel) { + yLabel = probV; + xLabel = m_a[j]; + } + prob.push_back(probV); sum_prob[j] += prob[j]; } - std::pair coords(xLabel-0.002, 1.05*yLabel); - labelPositions.push_back(coords); + std::pair coords(xLabel - 0.002, 1.05 * yLabel); + labelPositions.push_back(coords); TGraph* gr = new TGraph(m_a.size(), &m_a[0], &prob[0]); grp.push_back(gr); } - ///// PLOTS ///// // Plot the density scan grp[0]->SetLineColor(kBlue - 3); - if( title ) grp[0]->SetTitle("Transmission probability as a function of the axion mass"); - else grp[0]->SetTitle(""); + if (title) + grp[0]->SetTitle("Transmission probability as a function of the axion mass"); + else + grp[0]->SetTitle(""); grp[0]->GetXaxis()->SetTitle("m_{a} [eV]"); grp[0]->GetYaxis()->SetTitle("P_{a#gamma}"); - grp[0]->GetXaxis()->SetLimits(ma_min+0.00001, ma_max); + grp[0]->GetXaxis()->SetLimits(ma_min + 0.00001, ma_max); - grp[0]->GetXaxis()->SetLabelSize(0.04); + grp[0]->GetXaxis()->SetLabelSize(0.04); grp[0]->GetXaxis()->SetLabelFont(42); - grp[0]->GetXaxis()->SetTitleSize(0.045); + grp[0]->GetXaxis()->SetTitleSize(0.045); grp[0]->GetXaxis()->SetTitleFont(42); grp[0]->GetYaxis()->SetRangeUser(0, fProbMax); - grp[0]->GetYaxis()->SetLabelSize(0.04); + grp[0]->GetYaxis()->SetLabelSize(0.04); grp[0]->GetYaxis()->SetLabelFont(42); - grp[0]->GetYaxis()->SetTitleSize(0.05); + grp[0]->GetYaxis()->SetTitleSize(0.05); grp[0]->GetYaxis()->SetTitleFont(42); grp[0]->GetYaxis()->SetTitleOffset(0.9); - grp[0]->SetLineColor(kBlack); - grp[0]->SetFillColorAlpha(kRed,0.5); - grp[0]->SetLineWidth(1); + grp[0]->SetLineColor(kBlack); + grp[0]->SetFillColorAlpha(kRed, 0.5); + grp[0]->SetLineWidth(1); grp[0]->Draw("AFL"); - int n = 0; + int n = 0; for (const auto g : grp) { g->SetLineColor(kBlack); g->SetLineWidth(2); - g->SetFillColorAlpha(kRed,0.15); + g->SetFillColorAlpha(kRed, 0.15); g->Draw("FL SAME"); - if( labels) - { - std::string label = "P" + std::to_string(n+1); - TLatex* textLatex = new TLatex(labelPositions[n].first, labelPositions[n].second, label.c_str()); - textLatex->SetTextColor(1); - textLatex->SetTextSize(0.02); - textLatex->Draw("same"); - } - n++; + if (labels) { + std::string label = "P" + std::to_string(n + 1); + TLatex* textLatex = new TLatex(labelPositions[n].first, labelPositions[n].second, label.c_str()); + textLatex->SetTextColor(1); + textLatex->SetTextSize(0.02); + textLatex->Draw("same"); + } + n++; } // PLot of the sum of all the probabilities - TGraph *sum; - if( sumProb ) - { - sum = new TGraph(m_a.size(), &m_a[0], &sum_prob[0]); - sum->SetLineColor(kBlue + 3); - sum->SetLineWidth(3); - sum->Draw("SAME"); - } + TGraph* sum; + if (sumProb) { + sum = new TGraph(m_a.size(), &m_a[0], &sum_prob[0]); + sum->SetLineColor(kBlue + 3); + sum->SetLineWidth(3); + sum->Draw("SAME"); + } // Plot of the vacuum probability TGraph* vac = new TGraph(m_a.size(), &m_a[0], &prob_vac[0]); vac->SetLineColor(kBlack); vac->SetLineWidth(2); - vac->SetFillColorAlpha(kBlue,0.25); + vac->SetFillColorAlpha(kBlue, 0.25); vac->Draw("FL SAME"); - TGraph *gama; - if( nGamma ) // The following code could be a method of a future TRestAxionPlotResonances::DrawNGamma() - { - std ::vector NGamma; - - TRestAxionSolarQCDFlux* sFlux = new TRestAxionSolarQCDFlux("fluxes.rml", "LennertHoofPrimakoff"); - sFlux->Initialize(); - TH1F* spec = sFlux->GetEnergySpectrum(); - - /// The original gets 1eV binning to describe monocromatic lines. - /// We do not need such a high resolution - TH1F *rebinned = dynamic_cast(spec->Rebin(fReBinning,"rebinned")); - /// We need to renormalize since the bin contents are already renormalized to cm-2 s-1 keV-1 - rebinned->Scale(1./fReBinning); - Double_t deltaE = rebinned->GetBinCenter(2) - rebinned->GetBinCenter(1); - - if(!rebinned) { - std::cout << "Energy spectrum is nullptr!" << std::endl; - } - else - { - Double_t maxNGamma = 0; - // Obtain the Ngamma for each axion mass - for (size_t j = 0; j < m_a.size(); j++) { - std::cout << "Producing Ngamma for mass : " << m_a[j] << std::endl; - - Double_t ng = 0; - - // Vacuum - ax->AssignBufferGas(nullptr); - for( int n = 1; n <= rebinned->GetNbinsX(); n++ ) - { - Double_t energy = rebinned->GetBinCenter(n); - if ( energy < fEnergyRange.X() || energy > fEnergyRange.Y() ) - continue; - - ax->SetAxionEnergy( rebinned->GetBinCenter(n) ); - ng += rebinned->GetBinContent(n) * ax->GammaTransmissionProbability(m_a[j]); - } - - /// Gas phase - for (const auto& p : Psettings) { - // Creates the gas and assignss it to the axion field - gas->SetGasDensity(gasName, p.second); - ax->AssignBufferGas(gas); - - for( int n = 1; n <= rebinned->GetNbinsX(); n++ ) - { - Double_t energy = rebinned->GetBinCenter(n); - if ( energy < fEnergyRange.X() || energy > fEnergyRange.Y() ) - continue; - - ax->SetAxionEnergy( rebinned->GetBinCenter(n) ); - ng += rebinned->GetBinContent(n) * ax->GammaTransmissionProbability(m_a[j]); - } - } - ng = ng * deltaE * fArea * fExposureTime; - if( ng > maxNGamma ) maxNGamma = ng; - NGamma.push_back( ng ); - } - - gama = new TGraph(m_a.size(), &m_a[0], &NGamma[0]); - gama->SetLineColor(kRed + 3); - gama->SetLineWidth(3); - gama->GetYaxis()->SetRangeUser(0, fNGammaMax); - gama->Scale(fProbMax/fNGammaMax); - gama->Draw("SAME"); - - std::cout << "Max: " << maxNGamma << std::endl; - - /// 0.001 is to remove the first tick, which is zero, and overlaps with 0.1eV - TGaxis *A1 = new TGaxis(0.1,0,0.1,fProbMax,0.001,fNGammaMax,510,"+L"); - A1->SetTitle("N_{#gamma}"); - A1->SetTitleOffset(1.3); - A1->SetTitleSize(0.045); - A1->SetLabelSize(0.04); - A1->SetLabelFont(42); - A1->SetTextFont(42); - A1->Draw(); - } - } + TGraph* gama; + if (nGamma) // The following code could be a method of a future TRestAxionPlotResonances::DrawNGamma() + { + std ::vector NGamma; + + TRestAxionSolarQCDFlux* sFlux = new TRestAxionSolarQCDFlux("fluxes.rml", "LennertHoofPrimakoff"); + sFlux->Initialize(); + TH1F* spec = sFlux->GetEnergySpectrum(); + + /// The original gets 1eV binning to describe monocromatic lines. + /// We do not need such a high resolution + TH1F* rebinned = dynamic_cast(spec->Rebin(fReBinning, "rebinned")); + /// We need to renormalize since the bin contents are already renormalized to cm-2 s-1 keV-1 + rebinned->Scale(1. / fReBinning); + Double_t deltaE = rebinned->GetBinCenter(2) - rebinned->GetBinCenter(1); + + if (!rebinned) { + std::cout << "Energy spectrum is nullptr!" << std::endl; + } else { + Double_t maxNGamma = 0; + // Obtain the Ngamma for each axion mass + for (size_t j = 0; j < m_a.size(); j++) { + std::cout << "Producing Ngamma for mass : " << m_a[j] << std::endl; + + Double_t ng = 0; + + // Vacuum + ax->AssignBufferGas(nullptr); + for (int n = 1; n <= rebinned->GetNbinsX(); n++) { + Double_t energy = rebinned->GetBinCenter(n); + if (energy < fEnergyRange.X() || energy > fEnergyRange.Y()) continue; + + ax->SetAxionEnergy(rebinned->GetBinCenter(n)); + ng += rebinned->GetBinContent(n) * ax->GammaTransmissionProbability(m_a[j]); + } + + /// Gas phase + for (const auto& p : Psettings) { + // Creates the gas and assignss it to the axion field + gas->SetGasDensity(gasName, p.second); + ax->AssignBufferGas(gas); + + for (int n = 1; n <= rebinned->GetNbinsX(); n++) { + Double_t energy = rebinned->GetBinCenter(n); + if (energy < fEnergyRange.X() || energy > fEnergyRange.Y()) continue; + + ax->SetAxionEnergy(rebinned->GetBinCenter(n)); + ng += rebinned->GetBinContent(n) * ax->GammaTransmissionProbability(m_a[j]); + } + } + ng = ng * deltaE * fArea * fExposureTime; + if (ng > maxNGamma) maxNGamma = ng; + NGamma.push_back(ng); + } + + gama = new TGraph(m_a.size(), &m_a[0], &NGamma[0]); + gama->SetLineColor(kRed + 3); + gama->SetLineWidth(3); + gama->GetYaxis()->SetRangeUser(0, fNGammaMax); + gama->Scale(fProbMax / fNGammaMax); + gama->Draw("SAME"); + + std::cout << "Max: " << maxNGamma << std::endl; + + /// 0.001 is to remove the first tick, which is zero, and overlaps with 0.1eV + TGaxis* A1 = new TGaxis(0.1, 0, 0.1, fProbMax, 0.001, fNGammaMax, 510, "+L"); + A1->SetTitle("N_{#gamma}"); + A1->SetTitleOffset(1.3); + A1->SetTitleSize(0.045); + A1->SetLabelSize(0.04); + A1->SetLabelFont(42); + A1->SetTextFont(42); + A1->Draw(); + } + } // Plot of the vertical line - TLine* verticalLine; - if( drawLine ) - { - verticalLine = new TLine(Psettings[0].first, c1->GetUymin(), Psettings[0].first, fProbMax); - verticalLine->SetLineColor(kGreen - 3); - verticalLine->SetLineWidth(2); - verticalLine->Draw("same"); - } + TLine* verticalLine; + if (drawLine) { + verticalLine = new TLine(Psettings[0].first, c1->GetUymin(), Psettings[0].first, fProbMax); + verticalLine->SetLineColor(kGreen - 3); + verticalLine->SetLineWidth(2); + verticalLine->Draw("same"); + } // Plot of the legend - if( legend ) - { - TLegend* legend = new TLegend(0.87, 0.7, 0.47, 0.9); - legend->AddEntry(grp[0], "P_{a#gamma} for each P-step", "lf"); - legend->AddEntry(vac, "P_{a#gamma} for vacuum", "lf"); - if( sumProb ) legend->AddEntry(sum, "#Sigma P_{a#gamma}", "l"); - if( nGamma ) legend->AddEntry(gama, "N_{#gamma}", "l"); - if( drawLine ) legend->AddEntry(verticalLine, "m_{a} where P_{a#gamma}^{vac} = max(P_{a#gamma}^{vac}/2) ", "l"); - legend->Draw("same"); - c1->Draw(); - } + if (legend) { + TLegend* legend = new TLegend(0.87, 0.7, 0.47, 0.9); + legend->AddEntry(grp[0], "P_{a#gamma} for each P-step", "lf"); + legend->AddEntry(vac, "P_{a#gamma} for vacuum", "lf"); + if (sumProb) legend->AddEntry(sum, "#Sigma P_{a#gamma}", "l"); + if (nGamma) legend->AddEntry(gama, "N_{#gamma}", "l"); + if (drawLine) + legend->AddEntry(verticalLine, "m_{a} where P_{a#gamma}^{vac} = max(P_{a#gamma}^{vac}/2) ", "l"); + legend->Draw("same"); + c1->Draw(); + } c1->Print("/tmp/resonances.pdf"); From ff9093336335af15acfcd6944e75b905eec16b2e Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Tue, 5 Mar 2024 15:59:55 +0100 Subject: [PATCH 4/5] REST_Axion_PlotResonances.C Adding ABC flux --- macros/REST_Axion_PlotResonances.C | 98 ++++++++++++++++++++++++++---- 1 file changed, 85 insertions(+), 13 deletions(-) diff --git a/macros/REST_Axion_PlotResonances.C b/macros/REST_Axion_PlotResonances.C index aeaba846..62b477c9 100644 --- a/macros/REST_Axion_PlotResonances.C +++ b/macros/REST_Axion_PlotResonances.C @@ -133,7 +133,6 @@ int REST_Axion_PlotResonances(std::string optionString = "", double ma_max = 0.1 ///// PLOTS ///// // Plot the density scan - grp[0]->SetLineColor(kBlue - 3); if( title ) grp[0]->SetTitle("Transmission probability as a function of the axion mass"); else grp[0]->SetTitle(""); grp[0]->GetXaxis()->SetTitle("m_{a} [eV]"); @@ -161,7 +160,7 @@ int REST_Axion_PlotResonances(std::string optionString = "", double ma_max = 0.1 int n = 0; for (const auto g : grp) { g->SetLineColor(kBlack); - g->SetLineWidth(2); + g->SetLineWidth(1); g->SetFillColorAlpha(kRed,0.15); g->Draw("FL SAME"); @@ -181,23 +180,25 @@ int REST_Axion_PlotResonances(std::string optionString = "", double ma_max = 0.1 if( sumProb ) { sum = new TGraph(m_a.size(), &m_a[0], &sum_prob[0]); - sum->SetLineColor(kBlue + 3); - sum->SetLineWidth(3); + sum->SetLineColor(kBlue - 6); + sum->SetLineWidth(2); sum->Draw("SAME"); } // Plot of the vacuum probability TGraph* vac = new TGraph(m_a.size(), &m_a[0], &prob_vac[0]); vac->SetLineColor(kBlack); - vac->SetLineWidth(2); + vac->SetLineWidth(1); vac->SetFillColorAlpha(kBlue,0.25); vac->Draw("FL SAME"); - TGraph *gama; + TGraph *gama, *gamaABC; + Double_t deltaE = 0; if( nGamma ) // The following code could be a method of a future TRestAxionPlotResonances::DrawNGamma() { - std ::vector NGamma; + std ::vector NGamma, NGammaABC; + /////////// Integrating Ngamma for Primakoff //////////// TRestAxionSolarQCDFlux* sFlux = new TRestAxionSolarQCDFlux("fluxes.rml", "LennertHoofPrimakoff"); sFlux->Initialize(); TH1F* spec = sFlux->GetEnergySpectrum(); @@ -207,7 +208,7 @@ int REST_Axion_PlotResonances(std::string optionString = "", double ma_max = 0.1 TH1F *rebinned = dynamic_cast(spec->Rebin(fReBinning,"rebinned")); /// We need to renormalize since the bin contents are already renormalized to cm-2 s-1 keV-1 rebinned->Scale(1./fReBinning); - Double_t deltaE = rebinned->GetBinCenter(2) - rebinned->GetBinCenter(1); + deltaE = rebinned->GetBinCenter(2) - rebinned->GetBinCenter(1); if(!rebinned) { std::cout << "Energy spectrum is nullptr!" << std::endl; @@ -255,13 +256,13 @@ int REST_Axion_PlotResonances(std::string optionString = "", double ma_max = 0.1 } gama = new TGraph(m_a.size(), &m_a[0], &NGamma[0]); - gama->SetLineColor(kRed + 3); - gama->SetLineWidth(3); + gama->SetLineColor(kRed + 1); + gama->SetLineWidth(2); + gama->SetLineStyle(2); gama->GetYaxis()->SetRangeUser(0, fNGammaMax); gama->Scale(fProbMax/fNGammaMax); gama->Draw("SAME"); - std::cout << "Max: " << maxNGamma << std::endl; /// 0.001 is to remove the first tick, which is zero, and overlaps with 0.1eV TGaxis *A1 = new TGaxis(0.1,0,0.1,fProbMax,0.001,fNGammaMax,510,"+L"); @@ -273,6 +274,73 @@ int REST_Axion_PlotResonances(std::string optionString = "", double ma_max = 0.1 A1->SetTextFont(42); A1->Draw(); } + + /////////// Integrating Ngamma for ABC //////////// + TRestAxionSolarQCDFlux* sFluxABC = new TRestAxionSolarQCDFlux("fluxes.rml", "LennertHoofABC"); + + sFluxABC->Initialize(); + TH1F* specABC = sFluxABC->GetEnergySpectrum(); + + /// The original gets 1eV binning to describe monocromatic lines. + /// We do not need such a high resolution + TH1F *rebinnedABC = dynamic_cast(specABC->Rebin(fReBinning,"rebinnedABC")); + /// We need to renormalize since the bin contents are already renormalized to cm-2 s-1 keV-1 + rebinnedABC->Scale(1./fReBinning); + deltaE = rebinnedABC->GetBinCenter(2) - rebinnedABC->GetBinCenter(1); + + if(!rebinnedABC) { + std::cout << "Energy spectrum is nullptr!" << std::endl; + } + else + { + Double_t maxNGamma = 0; + // Obtain the Ngamma for each axion mass + for (size_t j = 0; j < m_a.size(); j++) { + std::cout << "Producing Ngamma (ABC) for mass : " << m_a[j] << std::endl; + + Double_t ng = 0; + + // Vacuum + ax->AssignBufferGas(nullptr); + for( int n = 1; n < rebinnedABC->GetNbinsX(); n++ ) + { + Double_t energy = rebinnedABC->GetBinCenter(n); + if ( energy < fEnergyRange.X() || energy > 2 ) + continue; + + ax->SetAxionEnergy( rebinnedABC->GetBinCenter(n) ); + ng += rebinnedABC->GetBinContent(n) * ax->GammaTransmissionProbability(m_a[j]); + } + + /// Gas phase + for (const auto& p : Psettings) { + // Creates the gas and assignss it to the axion field + gas->SetGasDensity(gasName, p.second); + ax->AssignBufferGas(gas); + + for( int n = 1; n < rebinnedABC->GetNbinsX(); n++ ) + { + Double_t energy = rebinnedABC->GetBinCenter(n); + if ( energy < fEnergyRange.X() || energy > 2 ) + continue; + + ax->SetAxionEnergy( rebinnedABC->GetBinCenter(n) ); + ng += rebinnedABC->GetBinContent(n) * ax->GammaTransmissionProbability(m_a[j]); + } + } + ng = 2*100 * ng * deltaE * fArea * fExposureTime; + if( ng > maxNGamma ) maxNGamma = ng; + NGammaABC.push_back( ng ); + } + + gamaABC = new TGraph(m_a.size(), &m_a[0], &NGammaABC[0]); + gamaABC->SetLineColor(kBlue+3); + gamaABC->SetLineStyle(2); + gamaABC->SetLineWidth(2); + gamaABC->GetYaxis()->SetRangeUser(0, fNGammaMax); + gamaABC->Scale(fProbMax/fNGammaMax); + gamaABC->Draw("SAME"); + } } // Plot of the vertical line @@ -280,7 +348,7 @@ int REST_Axion_PlotResonances(std::string optionString = "", double ma_max = 0.1 if( drawLine ) { verticalLine = new TLine(Psettings[0].first, c1->GetUymin(), Psettings[0].first, fProbMax); - verticalLine->SetLineColor(kGreen - 3); + verticalLine->SetLineColor(kGreen - 7); verticalLine->SetLineWidth(2); verticalLine->Draw("same"); } @@ -292,7 +360,11 @@ int REST_Axion_PlotResonances(std::string optionString = "", double ma_max = 0.1 legend->AddEntry(grp[0], "P_{a#gamma} for each P-step", "lf"); legend->AddEntry(vac, "P_{a#gamma} for vacuum", "lf"); if( sumProb ) legend->AddEntry(sum, "#Sigma P_{a#gamma}", "l"); - if( nGamma ) legend->AddEntry(gama, "N_{#gamma}", "l"); + if( nGamma ) + { + legend->AddEntry(gama, "N_{#gamma} (Primakoff)", "l"); + legend->AddEntry(gamaABC, "2x N_{#gamma} (ABC)", "l"); + } if( drawLine ) legend->AddEntry(verticalLine, "m_{a} where P_{a#gamma}^{vac} = max(P_{a#gamma}^{vac}/2) ", "l"); legend->Draw("same"); c1->Draw(); From 779f08693a6da7612c88f685b6cd176e206e1369 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Tue, 5 Mar 2024 15:03:03 +0000 Subject: [PATCH 5/5] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- macros/REST_Axion_PlotResonances.C | 260 ++++++++++++++--------------- 1 file changed, 125 insertions(+), 135 deletions(-) diff --git a/macros/REST_Axion_PlotResonances.C b/macros/REST_Axion_PlotResonances.C index e2260643..ebc2cc2c 100644 --- a/macros/REST_Axion_PlotResonances.C +++ b/macros/REST_Axion_PlotResonances.C @@ -163,7 +163,7 @@ int REST_Axion_PlotResonances(std::string optionString = "", double ma_max = 0.1 for (const auto g : grp) { g->SetLineColor(kBlack); g->SetLineWidth(1); - g->SetFillColorAlpha(kRed,0.15); + g->SetFillColorAlpha(kRed, 0.15); g->Draw("FL SAME"); if (labels) { @@ -177,39 +177,38 @@ int REST_Axion_PlotResonances(std::string optionString = "", double ma_max = 0.1 } // PLot of the sum of all the probabilities - TGraph *sum; - if( sumProb ) - { - sum = new TGraph(m_a.size(), &m_a[0], &sum_prob[0]); - sum->SetLineColor(kBlue - 6); - sum->SetLineWidth(2); - sum->Draw("SAME"); - } + TGraph* sum; + if (sumProb) { + sum = new TGraph(m_a.size(), &m_a[0], &sum_prob[0]); + sum->SetLineColor(kBlue - 6); + sum->SetLineWidth(2); + sum->Draw("SAME"); + } // Plot of the vacuum probability TGraph* vac = new TGraph(m_a.size(), &m_a[0], &prob_vac[0]); vac->SetLineColor(kBlack); vac->SetLineWidth(1); - vac->SetFillColorAlpha(kBlue,0.25); + vac->SetFillColorAlpha(kBlue, 0.25); vac->Draw("FL SAME"); - TGraph *gama, *gamaABC; - Double_t deltaE = 0; - if( nGamma ) // The following code could be a method of a future TRestAxionPlotResonances::DrawNGamma() - { - std ::vector NGamma, NGammaABC; + TGraph *gama, *gamaABC; + Double_t deltaE = 0; + if (nGamma) // The following code could be a method of a future TRestAxionPlotResonances::DrawNGamma() + { + std ::vector NGamma, NGammaABC; - /////////// Integrating Ngamma for Primakoff //////////// - TRestAxionSolarQCDFlux* sFlux = new TRestAxionSolarQCDFlux("fluxes.rml", "LennertHoofPrimakoff"); - sFlux->Initialize(); - TH1F* spec = sFlux->GetEnergySpectrum(); + /////////// Integrating Ngamma for Primakoff //////////// + TRestAxionSolarQCDFlux* sFlux = new TRestAxionSolarQCDFlux("fluxes.rml", "LennertHoofPrimakoff"); + sFlux->Initialize(); + TH1F* spec = sFlux->GetEnergySpectrum(); - /// The original gets 1eV binning to describe monocromatic lines. - /// We do not need such a high resolution - TH1F *rebinned = dynamic_cast(spec->Rebin(fReBinning,"rebinned")); - /// We need to renormalize since the bin contents are already renormalized to cm-2 s-1 keV-1 - rebinned->Scale(1./fReBinning); - deltaE = rebinned->GetBinCenter(2) - rebinned->GetBinCenter(1); + /// The original gets 1eV binning to describe monocromatic lines. + /// We do not need such a high resolution + TH1F* rebinned = dynamic_cast(spec->Rebin(fReBinning, "rebinned")); + /// We need to renormalize since the bin contents are already renormalized to cm-2 s-1 keV-1 + rebinned->Scale(1. / fReBinning); + deltaE = rebinned->GetBinCenter(2) - rebinned->GetBinCenter(1); if (!rebinned) { std::cout << "Energy spectrum is nullptr!" << std::endl; @@ -250,120 +249,111 @@ int REST_Axion_PlotResonances(std::string optionString = "", double ma_max = 0.1 NGamma.push_back(ng); } - gama = new TGraph(m_a.size(), &m_a[0], &NGamma[0]); - gama->SetLineColor(kRed + 1); - gama->SetLineWidth(2); - gama->SetLineStyle(2); - gama->GetYaxis()->SetRangeUser(0, fNGammaMax); - gama->Scale(fProbMax/fNGammaMax); - gama->Draw("SAME"); - - - /// 0.001 is to remove the first tick, which is zero, and overlaps with 0.1eV - TGaxis *A1 = new TGaxis(0.1,0,0.1,fProbMax,0.001,fNGammaMax,510,"+L"); - A1->SetTitle("N_{#gamma}"); - A1->SetTitleOffset(1.3); - A1->SetTitleSize(0.045); - A1->SetLabelSize(0.04); - A1->SetLabelFont(42); - A1->SetTextFont(42); - A1->Draw(); - } - - /////////// Integrating Ngamma for ABC //////////// - TRestAxionSolarQCDFlux* sFluxABC = new TRestAxionSolarQCDFlux("fluxes.rml", "LennertHoofABC"); - - sFluxABC->Initialize(); - TH1F* specABC = sFluxABC->GetEnergySpectrum(); - - /// The original gets 1eV binning to describe monocromatic lines. - /// We do not need such a high resolution - TH1F *rebinnedABC = dynamic_cast(specABC->Rebin(fReBinning,"rebinnedABC")); - /// We need to renormalize since the bin contents are already renormalized to cm-2 s-1 keV-1 - rebinnedABC->Scale(1./fReBinning); - deltaE = rebinnedABC->GetBinCenter(2) - rebinnedABC->GetBinCenter(1); - - if(!rebinnedABC) { - std::cout << "Energy spectrum is nullptr!" << std::endl; - } - else - { - Double_t maxNGamma = 0; - // Obtain the Ngamma for each axion mass - for (size_t j = 0; j < m_a.size(); j++) { - std::cout << "Producing Ngamma (ABC) for mass : " << m_a[j] << std::endl; - - Double_t ng = 0; - - // Vacuum - ax->AssignBufferGas(nullptr); - for( int n = 1; n < rebinnedABC->GetNbinsX(); n++ ) - { - Double_t energy = rebinnedABC->GetBinCenter(n); - if ( energy < fEnergyRange.X() || energy > 2 ) - continue; - - ax->SetAxionEnergy( rebinnedABC->GetBinCenter(n) ); - ng += rebinnedABC->GetBinContent(n) * ax->GammaTransmissionProbability(m_a[j]); - } - - /// Gas phase - for (const auto& p : Psettings) { - // Creates the gas and assignss it to the axion field - gas->SetGasDensity(gasName, p.second); - ax->AssignBufferGas(gas); - - for( int n = 1; n < rebinnedABC->GetNbinsX(); n++ ) - { - Double_t energy = rebinnedABC->GetBinCenter(n); - if ( energy < fEnergyRange.X() || energy > 2 ) - continue; - - ax->SetAxionEnergy( rebinnedABC->GetBinCenter(n) ); - ng += rebinnedABC->GetBinContent(n) * ax->GammaTransmissionProbability(m_a[j]); - } - } - ng = 2*100 * ng * deltaE * fArea * fExposureTime; - if( ng > maxNGamma ) maxNGamma = ng; - NGammaABC.push_back( ng ); - } - - gamaABC = new TGraph(m_a.size(), &m_a[0], &NGammaABC[0]); - gamaABC->SetLineColor(kBlue+3); - gamaABC->SetLineStyle(2); - gamaABC->SetLineWidth(2); - gamaABC->GetYaxis()->SetRangeUser(0, fNGammaMax); - gamaABC->Scale(fProbMax/fNGammaMax); - gamaABC->Draw("SAME"); - } - } + gama = new TGraph(m_a.size(), &m_a[0], &NGamma[0]); + gama->SetLineColor(kRed + 1); + gama->SetLineWidth(2); + gama->SetLineStyle(2); + gama->GetYaxis()->SetRangeUser(0, fNGammaMax); + gama->Scale(fProbMax / fNGammaMax); + gama->Draw("SAME"); + + /// 0.001 is to remove the first tick, which is zero, and overlaps with 0.1eV + TGaxis* A1 = new TGaxis(0.1, 0, 0.1, fProbMax, 0.001, fNGammaMax, 510, "+L"); + A1->SetTitle("N_{#gamma}"); + A1->SetTitleOffset(1.3); + A1->SetTitleSize(0.045); + A1->SetLabelSize(0.04); + A1->SetLabelFont(42); + A1->SetTextFont(42); + A1->Draw(); + } + + /////////// Integrating Ngamma for ABC //////////// + TRestAxionSolarQCDFlux* sFluxABC = new TRestAxionSolarQCDFlux("fluxes.rml", "LennertHoofABC"); + + sFluxABC->Initialize(); + TH1F* specABC = sFluxABC->GetEnergySpectrum(); + + /// The original gets 1eV binning to describe monocromatic lines. + /// We do not need such a high resolution + TH1F* rebinnedABC = dynamic_cast(specABC->Rebin(fReBinning, "rebinnedABC")); + /// We need to renormalize since the bin contents are already renormalized to cm-2 s-1 keV-1 + rebinnedABC->Scale(1. / fReBinning); + deltaE = rebinnedABC->GetBinCenter(2) - rebinnedABC->GetBinCenter(1); + + if (!rebinnedABC) { + std::cout << "Energy spectrum is nullptr!" << std::endl; + } else { + Double_t maxNGamma = 0; + // Obtain the Ngamma for each axion mass + for (size_t j = 0; j < m_a.size(); j++) { + std::cout << "Producing Ngamma (ABC) for mass : " << m_a[j] << std::endl; + + Double_t ng = 0; + + // Vacuum + ax->AssignBufferGas(nullptr); + for (int n = 1; n < rebinnedABC->GetNbinsX(); n++) { + Double_t energy = rebinnedABC->GetBinCenter(n); + if (energy < fEnergyRange.X() || energy > 2) continue; + + ax->SetAxionEnergy(rebinnedABC->GetBinCenter(n)); + ng += rebinnedABC->GetBinContent(n) * ax->GammaTransmissionProbability(m_a[j]); + } + + /// Gas phase + for (const auto& p : Psettings) { + // Creates the gas and assignss it to the axion field + gas->SetGasDensity(gasName, p.second); + ax->AssignBufferGas(gas); + + for (int n = 1; n < rebinnedABC->GetNbinsX(); n++) { + Double_t energy = rebinnedABC->GetBinCenter(n); + if (energy < fEnergyRange.X() || energy > 2) continue; + + ax->SetAxionEnergy(rebinnedABC->GetBinCenter(n)); + ng += rebinnedABC->GetBinContent(n) * ax->GammaTransmissionProbability(m_a[j]); + } + } + ng = 2 * 100 * ng * deltaE * fArea * fExposureTime; + if (ng > maxNGamma) maxNGamma = ng; + NGammaABC.push_back(ng); + } + + gamaABC = new TGraph(m_a.size(), &m_a[0], &NGammaABC[0]); + gamaABC->SetLineColor(kBlue + 3); + gamaABC->SetLineStyle(2); + gamaABC->SetLineWidth(2); + gamaABC->GetYaxis()->SetRangeUser(0, fNGammaMax); + gamaABC->Scale(fProbMax / fNGammaMax); + gamaABC->Draw("SAME"); + } + } // Plot of the vertical line - TLine* verticalLine; - if( drawLine ) - { - verticalLine = new TLine(Psettings[0].first, c1->GetUymin(), Psettings[0].first, fProbMax); - verticalLine->SetLineColor(kGreen - 7); - verticalLine->SetLineWidth(2); - verticalLine->Draw("same"); - } + TLine* verticalLine; + if (drawLine) { + verticalLine = new TLine(Psettings[0].first, c1->GetUymin(), Psettings[0].first, fProbMax); + verticalLine->SetLineColor(kGreen - 7); + verticalLine->SetLineWidth(2); + verticalLine->Draw("same"); + } // Plot of the legend - if( legend ) - { - TLegend* legend = new TLegend(0.87, 0.7, 0.47, 0.9); - legend->AddEntry(grp[0], "P_{a#gamma} for each P-step", "lf"); - legend->AddEntry(vac, "P_{a#gamma} for vacuum", "lf"); - if( sumProb ) legend->AddEntry(sum, "#Sigma P_{a#gamma}", "l"); - if( nGamma ) - { - legend->AddEntry(gama, "N_{#gamma} (Primakoff)", "l"); - legend->AddEntry(gamaABC, "2x N_{#gamma} (ABC)", "l"); - } - if( drawLine ) legend->AddEntry(verticalLine, "m_{a} where P_{a#gamma}^{vac} = max(P_{a#gamma}^{vac}/2) ", "l"); - legend->Draw("same"); - c1->Draw(); - } + if (legend) { + TLegend* legend = new TLegend(0.87, 0.7, 0.47, 0.9); + legend->AddEntry(grp[0], "P_{a#gamma} for each P-step", "lf"); + legend->AddEntry(vac, "P_{a#gamma} for vacuum", "lf"); + if (sumProb) legend->AddEntry(sum, "#Sigma P_{a#gamma}", "l"); + if (nGamma) { + legend->AddEntry(gama, "N_{#gamma} (Primakoff)", "l"); + legend->AddEntry(gamaABC, "2x N_{#gamma} (ABC)", "l"); + } + if (drawLine) + legend->AddEntry(verticalLine, "m_{a} where P_{a#gamma}^{vac} = max(P_{a#gamma}^{vac}/2) ", "l"); + legend->Draw("same"); + c1->Draw(); + } c1->Print("/tmp/resonances.pdf");