diff --git a/README.md b/README.md index b662017c..da2013e1 100644 --- a/README.md +++ b/README.md @@ -69,3 +69,21 @@ This repository makes use of the following published codes: - K. Altenmuller et al, REST-for-Physics, a ROOT-based framework for event oriented data analysis and combined Monte Carlo response, [Computer Physics Communications 273, April 2022, 108281](https://doi.org/10.1016/j.cpc.2021.108281). - S.Hoof, J.Jaeckel, T.J.Lennert, Quantifying uncertainties in the solar axion flux and their impact on determining axion model parameters, [JCAP09(2021)006](https://doi.org/10.1088/1475-7516/2021/09/006). - T.Kittelmann, E.Klinkby, E.B.Knudsen, P.Willendrup, X.X.Cai, K.Kanaki, Monte Carlo Particle Lists: MCPL, Computer Physics Communications 218 (2017) 17–42](http://dx.doi.org/10.17632/cby92vsv5g.1). + +### Gallery + +

+ +cover1 +cover2 +
+cover3 +
+cover4 +
+cover5 +
+cover6 +
+cover7 +cover8 diff --git a/examples/03.IAXO/BabyIAXO_2024.rml b/examples/03.IAXO/BabyIAXO_2024.rml new file mode 100644 index 00000000..cf4860bd --- /dev/null +++ b/examples/03.IAXO/BabyIAXO_2024.rml @@ -0,0 +1,110 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/examples/03.IAXO/GenerateSignalComponents.C b/examples/03.IAXO/GenerateSignalComponents.C index 53a5e077..a6328681 100644 --- a/examples/03.IAXO/GenerateSignalComponents.C +++ b/examples/03.IAXO/GenerateSignalComponents.C @@ -1,3 +1,5 @@ +#include + Double_t Eo = 0.5; // keV Double_t Ef = 10; // keV diff --git a/examples/03.IAXO/response.rml b/examples/03.IAXO/response.rml new file mode 100644 index 00000000..b49d91c4 --- /dev/null +++ b/examples/03.IAXO/response.rml @@ -0,0 +1,5 @@ + + + + + diff --git a/images/paper/cover1.png b/images/paper/cover1.png new file mode 100644 index 00000000..e9b3bdbe Binary files /dev/null and b/images/paper/cover1.png differ diff --git a/images/paper/cover2.png b/images/paper/cover2.png new file mode 100644 index 00000000..4c5ff4ae Binary files /dev/null and b/images/paper/cover2.png differ diff --git a/images/paper/cover3.png b/images/paper/cover3.png new file mode 100644 index 00000000..dc4c4efb Binary files /dev/null and b/images/paper/cover3.png differ diff --git a/images/paper/cover4.png b/images/paper/cover4.png new file mode 100644 index 00000000..4059f653 Binary files /dev/null and b/images/paper/cover4.png differ diff --git a/images/paper/cover5.png b/images/paper/cover5.png new file mode 100644 index 00000000..ab3b04c8 Binary files /dev/null and b/images/paper/cover5.png differ diff --git a/images/paper/cover6.png b/images/paper/cover6.png new file mode 100644 index 00000000..e0345a32 Binary files /dev/null and b/images/paper/cover6.png differ diff --git a/images/paper/cover7.png b/images/paper/cover7.png new file mode 100644 index 00000000..f07d1dec Binary files /dev/null and b/images/paper/cover7.png differ diff --git a/images/paper/cover8.png b/images/paper/cover8.png new file mode 100644 index 00000000..99bdcd76 Binary files /dev/null and b/images/paper/cover8.png differ diff --git a/inc/TRestAxionHelioscopeSignal.h b/inc/TRestAxionHelioscopeSignal.h index 985fd4fd..accf37da 100644 --- a/inc/TRestAxionHelioscopeSignal.h +++ b/inc/TRestAxionHelioscopeSignal.h @@ -53,6 +53,12 @@ class TRestAxionHelioscopeSignal : public TRestComponent { /// If an x-ray window is present we may add an efficiency for Ngamma calculation Double_t fWindowEfficiency = 1; + /// The overall detector efficiency + Double_t fDetectorEfficiency = 1; + + /// The additional gas length photons need to travel to reach a vacuum region (mm) + Double_t fGasLength = 0; + /// It defines the gas mixture we use inside our magnetic field. Vacuum if it is nullptr TRestAxionBufferGas* fGas = nullptr; @@ -80,12 +86,14 @@ class TRestAxionHelioscopeSignal : public TRestComponent { void SetMagnetRadius(const Double_t& radius) { fMagnetRadius = radius; } void SetMagnetLength(const Double_t& length) { fMagnetLength = length; } void SetMagnetStrength(const Double_t& strength) { fMagnetStrength = strength; } + void SetGasLength(const Double_t& length) { fGasLength = length; } void SetType(const std::string& type) { fConversionType = type; } Int_t GetNumberOfBores() const { return fBores; } Double_t GetMagnetRadius() const { return fMagnetRadius; } Double_t GetMagnetLength() const { return fMagnetLength; } Double_t GetMagnetStrength() const { return fMagnetStrength; } + Double_t GetGasLength() const { return fGasLength; } std::string GetType() const { return fConversionType; } void Initialize() override; @@ -93,6 +101,6 @@ class TRestAxionHelioscopeSignal : public TRestComponent { TRestAxionHelioscopeSignal(); ~TRestAxionHelioscopeSignal(); - ClassDefOverride(TRestAxionHelioscopeSignal, 1); + ClassDefOverride(TRestAxionHelioscopeSignal, 2); }; #endif diff --git a/inc/TRestAxionMagneticField.h b/inc/TRestAxionMagneticField.h index b2dd1e54..e088ba82 100644 --- a/inc/TRestAxionMagneticField.h +++ b/inc/TRestAxionMagneticField.h @@ -118,6 +118,7 @@ class TRestAxionMagneticField : public TRestMetadata { void ReMap(const size_t& n, const TVector3& newMapSize); void SetTrack(const TVector3& position, const TVector3& direction); + void SetUserTrack(const TVector3& start, const TVector3& end); Double_t GetTrackLength() const { return fTrackLength; } TVector3 GetTrackStart() const { return fTrackStart; } diff --git a/inc/TRestAxionOpticsMirror.h b/inc/TRestAxionOpticsMirror.h index 3ee5936a..90c957db 100644 --- a/inc/TRestAxionOpticsMirror.h +++ b/inc/TRestAxionOpticsMirror.h @@ -99,6 +99,8 @@ class TRestAxionOpticsMirror : public TRestMetadata { TCanvas* DrawOpticsProperties(std::string options = "", Double_t lowRange = 1.e-9, Double_t lowRange2 = 1.e-3); + TPad* DrawOpticsPropertiesLinear(std::string options = "", Double_t lowRange = 1.e-9, + Double_t lowRange2 = 1.e-3); TRestAxionOpticsMirror(); TRestAxionOpticsMirror(const char* cfgFileName, std::string name = ""); diff --git a/inc/TRestAxionTrueWolterOptics.h b/inc/TRestAxionTrueWolterOptics.h index 5dece507..4a32326e 100644 --- a/inc/TRestAxionTrueWolterOptics.h +++ b/inc/TRestAxionTrueWolterOptics.h @@ -143,6 +143,13 @@ class TRestAxionTrueWolterOptics : public TRestAxionOptics { return 0; } + /// It returns the value of max entrance radius in mm + Double_t GetMaxEntranceRadius() { return GetR1().back(); } + + /// It returns the value of min entrance radius in mm + Double_t GetMinEntranceRadius() { return GetR1().front(); } + + /// It returns the value min/max entrance radius in mm as a std::pair std::pair GetRadialLimits() override { std::pair result(0, 0); if (!fR1.empty()) { @@ -151,6 +158,8 @@ class TRestAxionTrueWolterOptics : public TRestAxionOptics { return result; } + TRestSpiderMask* const& GetSpiderMask() const { return fSpiderMask; } + void SetMirror() override { Double_t x = fEntrancePosition.X(); Double_t y = fEntrancePosition.Y(); diff --git a/macros/REST_Axion_AccurateEfficiencies.C b/macros/REST_Axion_AccurateEfficiencies.C new file mode 100644 index 00000000..0f915a93 --- /dev/null +++ b/macros/REST_Axion_AccurateEfficiencies.C @@ -0,0 +1,175 @@ +#include "TCanvas.h" +#include "TGraph.h" +#include "TLatex.h" +#include "TLegend.h" +#include "TLine.h" +#include "TRestAxionBufferGas.h" +#include "TRestAxionField.h" + +Double_t fromEnergy = 0.5; +Double_t toEnergy = 10; + +Double_t incidenceAngle = 0.2; +Double_t deltaE = 0.1; + +//******************************************************************************************************* +//*** Description: This script will launch the integration of the axion-field with given parameters. +//*** It allows to test different magnetic field cell sizes, for a given mass that can be off-resonance +//*** for dm different from zero, and a given maximum tolerance or error for the integration routine. +//*** +//*** The macro sets the TRestAxionField under debug mode to print the different results on screen. +//*** +//*** -------------- +//*** Usage: restManager FieldIntegrationTests [sX=10] [sX=10] [sZ=10] [dm=0.01] [tolerance=0.1] [Ea=4.2] +//*** -------------- +//*** +//*** Author: Javier Galan +//******************************************************************************************************* +int REST_Axion_AccurateEfficiencies(std::string fluxFile = "fluxes.rml", + std::string fluxName = "LennertHoofPrimakoff", + std::string opticsFile = "xmmTrueWolter.rml", + std::string opticsName = "xmm") { + TRestAxionTrueWolterOptics optics(opticsFile.c_str(), opticsName.c_str()); + TRestAxionOpticsMirror* mirror = optics.GetMirrorProperties(); + + TRestAxionSolarQCDFlux flux(fluxFile.c_str(), fluxName.c_str()); + flux.Initialize(); + + Double_t R2sum = 0; + Double_t fluxSum = 0; + Double_t maxFlux = 0; + for (Double_t energy = fromEnergy; energy < toEnergy; energy += deltaE) { + Double_t R = mirror->GetReflectivity(incidenceAngle, energy); + fluxSum += flux.GetFluxAtEnergy(energy, 0); + if (flux.GetFluxAtEnergy(energy, 0) > maxFlux) maxFlux = flux.GetFluxAtEnergy(energy, 0); + R2sum += flux.GetFluxAtEnergy(energy, 0) * R * R; + } + + Double_t R2eff = R2sum / fluxSum; + std::cout << "R2eff: " << R2eff << std::endl; + + Double_t Rout = optics.GetMaxEntranceRadius(); + Double_t Rin = optics.GetMinEntranceRadius(); + + TRestSpiderMask* sMask = optics.GetSpiderMask(); + Double_t na = sMask->GetNumberOfArms(); + Double_t wa = sMask->GetArmsWidth(); + + std::vector r = optics.GetR1(); + std::vector th = optics.GetThickness(); + + Double_t Aeff = (TMath::Pi() - 0.5 * na * wa) * (Rout * Rout - Rin * Rin); + for (size_t n = 0; n < r.size(); n++) Aeff -= (2 * TMath::Pi() - na * wa) * r[n] * th[n]; + + std::cout << "Aeff (optics): " << Aeff / Rout / Rout / TMath::Pi() << std::endl; + + TCanvas c; + c.SetCanvasSize(2400, 1800); + c.SetWindowSize(2400, 1800); + c.Divide(2, 1); + + c.cd(1); + optics.GetMirrorProperties()->DrawOpticsPropertiesLinear(); + + c.cd(2); + optics.DrawParticleTracks(); + + c.Print("optics.pdf"); + + /// Extracted from x-ray window MicromegasStrongBack + na = 8; + wa = 2.64 * TMath::Pi() / 180.; + Rout = 8.5; + Rin = 4.55; + Double_t Ro = 4.25; + + Aeff = (TMath::Pi() - 0.5 * na * wa) * (Rout * Rout - Rin * Rin) + TMath::Pi() * Ro * Ro; + + std::cout << "Aeff (window): " << Aeff / Rout / Rout / TMath::Pi() << std::endl; + + TRestAxionXrayWindow strongBack("windows.rml", "MicromegasStrongBack"); + TRestAxionXrayWindow mylar("windows.rml", "MicromegasMylar"); + TRestAxionXrayWindow aluminum("windows.rml", "MicromegasAluminumFoil"); + + TGraph* mylarGraph = new TGraph(); + mylarGraph->SetName("Mylar"); + mylarGraph->SetLineColor(49); + mylarGraph->SetLineWidth(2); + + TGraph* aluminumGraph = new TGraph(); + aluminumGraph->SetName("Aluminum"); + aluminumGraph->SetLineColor(46); + aluminumGraph->SetLineWidth(2); + + TGraph* solarGraph = new TGraph(); + solarGraph->SetName("SolarFlux"); + solarGraph->SetLineColor(43); + solarGraph->SetLineWidth(2); + + Double_t WeffSum = 0; + Double_t AlSum = 0; + Double_t MySum = 0; + for (Double_t energy = deltaE; energy < toEnergy; energy += deltaE) { + Double_t tMy = mylar.GetTransmission(energy, 0, 0); + Double_t tAl = aluminum.GetTransmission(energy, 0, 0); + + mylarGraph->SetPoint(mylarGraph->GetN(), energy, tMy); + aluminumGraph->SetPoint(aluminumGraph->GetN(), energy, tAl); + solarGraph->SetPoint(solarGraph->GetN(), energy, flux.GetFluxAtEnergy(energy, 0) / maxFlux); + + WeffSum += flux.GetFluxAtEnergy(energy, 0) * tMy * tAl; + AlSum += flux.GetFluxAtEnergy(energy, 0) * tAl; + MySum += flux.GetFluxAtEnergy(energy, 0) * tMy; + } + + WeffSum = WeffSum / fluxSum; + AlSum = AlSum / fluxSum; + MySum = MySum / fluxSum; + std::cout << "AlSum: " << AlSum << std::endl; + std::cout << "MySum: " << MySum << std::endl; + std::cout << "WeffSum: " << WeffSum << std::endl; + + TCanvas c2; + c2.SetCanvasSize(1200, 900); + c2.SetWindowSize(1200, 900); + // c2.SetLogy(); + + TPad* pad2 = new TPad("pad1", "This is pad1", 0.01, 0.02, 0.99, 0.97); + // pad1->Divide(2, 2); + pad2->SetLogy(); + pad2->Draw(); + + ////// Drawing reflectivity versus angle + pad2->SetRightMargin(0.09); + pad2->SetLeftMargin(0.25); + pad2->SetBottomMargin(0.15); + + mylarGraph->GetXaxis()->SetLimits(0, 10); + // mylarGraph->GetHistogram()->SetMaximum(1); + mylarGraph->GetHistogram()->SetMinimum(0); + + mylarGraph->GetXaxis()->SetTitle("Energy [keV]"); + mylarGraph->GetXaxis()->SetTitleSize(0.04); + mylarGraph->GetXaxis()->SetLabelSize(0.04); + mylarGraph->GetYaxis()->SetTitle("Transmission"); + mylarGraph->GetYaxis()->SetTitleOffset(1.2); + mylarGraph->GetYaxis()->SetTitleSize(0.04); + mylarGraph->GetYaxis()->SetLabelSize(0.04); + mylarGraph->Draw("AL"); + aluminumGraph->Draw("L"); + // solarGraph->Draw("L"); + + Double_t lx1 = 0.6, ly1 = 0.55, lx2 = 0.8, ly2 = 0.75; + TLegend* legend = new TLegend(lx1, ly1, lx2, ly2); + legend->SetTextSize(0.03); + // legend->SetHeader("Widnows", "C"); // option "C" allows to center the header + legend->AddEntry("Mylar", "Mylar", "l"); + legend->AddEntry("Aluminum", "Aluminum", "l"); + legend->Draw(); + + c2.Print("windows.pdf"); + + c.Print("tracks.png"); + + return 0; +} diff --git a/macros/REST_Axion_ComputingTimesPerMass.C b/macros/REST_Axion_ComputingTimesPerMass.C new file mode 100644 index 00000000..29b4d358 --- /dev/null +++ b/macros/REST_Axion_ComputingTimesPerMass.C @@ -0,0 +1,92 @@ +#include +#include +#include + +#include "TGeoManager.h" +#include "TRestTask.h" +#include "TSystem.h" +using namespace std; +namespace fs = std::filesystem; + +#ifndef RESTTask_ComputingTimesPerMass +#define RESTTask_ComputingTimesPerMass + +//******************************************************************************************************* +//*** +//*** Description: +//*** +//*** -------------- +//*** +//*** Bla bla bla +//*** -------------- +//*** Usage: restManager ComputingTimesPerMass "/full/path/file_*pattern*.root" +//*** -------------- +//*** +//******************************************************************************************************* +Int_t REST_Axion_ComputingTimesPerMass(std::string namePattern) { + std::map statistics; + std::map computingTime; + + std::vector b = TRestTools::GetFilesMatchingPattern(namePattern, true); + + Double_t totalTime = 0; + int cont = 0; + for (int i = 0; i < b.size(); i++) { + string filename = b[i]; + std::cout << filename << std::endl; + cont++; + TRestRun* run = new TRestRun(); + TRestAxionGeneratorProcess* gen; + + TFile* f = new TFile(filename.c_str()); + + ///////////////////////////// + // Reading metadata classes + + TIter nextkey(f->GetListOfKeys()); + TKey* key; + while ((key = (TKey*)nextkey())) { + string className = key->GetClassName(); + if (className == "TRestRun") { + run = (TRestRun*)f->Get(key->GetName()); + } + if (className == "TRestAxionGeneratorProcess") { + gen = (TRestAxionGeneratorProcess*)f->Get(key->GetName()); + } + } + + double mass = StringTo2DVector(gen->GetDataMemberValue("fAxionMassRange")).X(); + + if (run->GetEndTimestamp() == 0 || run->GetRunLength() < 0 || run->GetEntries() == 0) { + continue; + } else if (run->GetRunLength() > 0) { + mass = mass * 1000.; + statistics[mass] += run->GetEntries(); + computingTime[mass] += run->GetRunLength() / 3600.; + } + + delete run; + + f->Close(); + } + + std::cout << "Contents of the stats map:" << std::endl; + for (const auto& [mass, stats] : statistics) { + std::cout << "Key: " << mass << ", Value: " << stats << std::endl; + } + std::cout << "Contents of the computing map:" << std::endl; + for (const auto& [mass, cTime] : computingTime) { + std::cout << "Key: " << mass << ", Value: " << cTime << std::endl; + } + + FILE* f = fopen("stats.txt", "wt"); + for (const auto& [mass, stats] : statistics) fprintf(f, "%lf\t%d\n", mass, stats); + fclose(f); + + f = fopen("computing.txt", "wt"); + for (const auto& [mass, cTime] : computingTime) fprintf(f, "%lf\t%lf\n", mass, cTime); + fclose(f); + + return 0; +} +#endif diff --git a/macros/REST_Axion_FieldIntegrationTests.C b/macros/REST_Axion_FieldIntegrationTests.C index bd2e392e..4a32f8e9 100644 --- a/macros/REST_Axion_FieldIntegrationTests.C +++ b/macros/REST_Axion_FieldIntegrationTests.C @@ -18,6 +18,9 @@ //*** //*** Author: Javier Galan //******************************************************************************************************* +Bool_t performCutOffTest = false; // To be enabled in case we need to reproduce the Z-cutoff curves published + // in the RayTracing paper + int REST_Axion_FieldIntegrationTests(Double_t sX = 10, Double_t sY = 10, Double_t sZ = 50, Double_t dm = 0.01, Double_t tolerance = 0.1, Double_t Ea = 4.2) { /// Setting up magnetic field and track to evaluate @@ -46,5 +49,23 @@ int REST_Axion_FieldIntegrationTests(Double_t sX = 10, Double_t sY = 10, Double_ std::pair prob = ax->GammaTransmissionFieldMapProbability(Ea, ma - dm, tolerance, 100, 25); + if (performCutOffTest) { + ax->SetDebug(false); + for (Double_t zStart = -10000; zStart < -4000; zStart += 100) { + magneticField.SetUserTrack(TVector3(0, 0, zStart), TVector3(0, 0, -zStart)); + + std::pair prob = + ax->GammaTransmissionFieldMapProbability(Ea, ma - dm, tolerance, 100, 25); + + std::cout << "zStart: " << zStart << " Prob: " << prob.first << std::endl; + } + + for (Double_t zStart = -10000; zStart < -4000; zStart += 100) { + Double_t field = magneticField.GetTransversalComponent(TVector3(0, 0, zStart), TVector3(0, 0, 1)); + + std::cout << "zStart: " << zStart << " Field: " << field << std::endl; + } + } + return 0; } diff --git a/src/TRestAxionField.cxx b/src/TRestAxionField.cxx index d6713d32..a2105cb9 100644 --- a/src/TRestAxionField.cxx +++ b/src/TRestAxionField.cxx @@ -568,7 +568,7 @@ std::pair TRestAxionField::ComputeResonanceIntegral(Double_t if (status > 0) return {0, status}; auto end = std::chrono::system_clock::now(); - auto seconds = std::chrono::duration_cast(end - start); + auto seconds = std::chrono::duration_cast(end - start); double GammaL = Gamma * fMagneticField->GetTrackLength(); double C = exp(-GammaL) * BLHalfSquared(1, 1); @@ -631,6 +631,7 @@ std::pair TRestAxionField::ComputeOffResonanceIntegral(Doubl gsl_integration_qawo(&F, 0, accuracy, accuracy, num_intervals, workspace, wf, &reprob, &rerr); if (status > 0) { gsl_integration_qawo_table_free(wf); + std::cout << "Status1: " << status << std::endl; return {0, status}; } @@ -638,6 +639,7 @@ std::pair TRestAxionField::ComputeOffResonanceIntegral(Doubl status = gsl_integration_qawo(&F, 0, accuracy, accuracy, num_intervals, workspace, wf, &improb, &imerr); if (status > 0) { gsl_integration_qawo_table_free(wf); + std::cout << "Status2: " << status << std::endl; return {0, status}; } diff --git a/src/TRestAxionHelioscopeSignal.cxx b/src/TRestAxionHelioscopeSignal.cxx index a23e8dbb..b85d4f72 100644 --- a/src/TRestAxionHelioscopeSignal.cxx +++ b/src/TRestAxionHelioscopeSignal.cxx @@ -153,8 +153,8 @@ Double_t TRestAxionHelioscopeSignal::GetSignalRate(std::vector point, Double_t probability = 0; if (fConversionType == "IAXO") { - probability = - fOpticsEfficiency * fWindowEfficiency * fField->GammaTransmissionProbability(point[0], mass); + probability = fOpticsEfficiency * fWindowEfficiency * fDetectorEfficiency * + fField->GammaTransmissionProbability(point[0], mass); // We assume all flux ends up inside the spot. No XY dependency of signal. Double_t apertureArea = TMath::Pi() * fMagnetRadius * units("cm") * fMagnetRadius * units("cm"); @@ -191,6 +191,9 @@ Double_t TRestAxionHelioscopeSignal::GetSignalRate(Double_t mass, Double_t Eo, D probability = fOpticsEfficiency * fWindowEfficiency * fField->GammaTransmissionProbability(en, mass); + if (fGas) + probability = probability * fGasLength * units("cm") * fGas->GetPhotonAbsorptionLength(en); + // We assume all flux ends up inside the spot. No XY dependency of signal. Double_t apertureArea = TMath::Pi() * fMagnetRadius * units("cm") * fMagnetRadius * units("cm"); flux *= fBores * apertureArea; @@ -280,8 +283,12 @@ void TRestAxionHelioscopeSignal::PrintMetadata() { RESTMetadata << "Magnet field : " << fMagnetStrength * units("T") << " T" << RESTendl; RESTMetadata << " " << RESTendl; + RESTMetadata << "Additional buffer gas behind field : " << fGasLength * units("cm") << " cm" << RESTendl; + RESTMetadata << " " << RESTendl; + RESTMetadata << "Optics efficiency : " << fOpticsEfficiency << RESTendl; RESTMetadata << "Window efficiency : " << fWindowEfficiency << RESTendl; + RESTMetadata << "Detector efficiency : " << fDetectorEfficiency << RESTendl; RESTMetadata << "----" << RESTendl; } diff --git a/src/TRestAxionMagneticField.cxx b/src/TRestAxionMagneticField.cxx index f16fc252..aaa94a41 100644 --- a/src/TRestAxionMagneticField.cxx +++ b/src/TRestAxionMagneticField.cxx @@ -1335,6 +1335,16 @@ void TRestAxionMagneticField::SetTrack(const TVector3& position, const TVector3& fTrackDirection = (trackBounds[1] - trackBounds[0]).Unit(); } +/////////////////////////////////////////////// +/// \brief It initializes the field track using the initial and final positions specified by +/// the user. +/// +void TRestAxionMagneticField::SetUserTrack(const TVector3& start, const TVector3& end) { + fTrackStart = start; + fTrackLength = (end - start).Mag() - 1; + fTrackDirection = (end - start).Unit(); +} + /////////////////////////////////////////////// /// \brief It will return the transversal magnetic field component evaluated at a parametric /// distance `x` (given by argument) for the track defined inside the class. The track will diff --git a/src/TRestAxionOptics.cxx b/src/TRestAxionOptics.cxx index 2c61951b..be2bc3eb 100644 --- a/src/TRestAxionOptics.cxx +++ b/src/TRestAxionOptics.cxx @@ -565,9 +565,10 @@ TPad* TRestAxionOptics::CreatePad(Int_t nx, Int_t ny) { fPad = new TPad("optics_pad", "This is the optics drawing pad", 0.01, 0.02, 0.99, 0.97); if (nx > 1 || ny > 1) fPad->Divide(nx, ny); fPad->Draw(); - fPad->SetRightMargin(0.09); - fPad->SetLeftMargin(0.2); - fPad->SetBottomMargin(0.15); + fPad->SetRightMargin(0.05); + fPad->SetLeftMargin(0.18); + fPad->SetBottomMargin(0.1); + fPad->SetTopMargin(0.02); return fPad; } diff --git a/src/TRestAxionOpticsMirror.cxx b/src/TRestAxionOpticsMirror.cxx index 2ca01c51..5bdf4986 100644 --- a/src/TRestAxionOpticsMirror.cxx +++ b/src/TRestAxionOpticsMirror.cxx @@ -719,3 +719,178 @@ TCanvas* TRestAxionOpticsMirror::DrawOpticsProperties(std::string options, Doubl lat->DrawLatexNDC(.5, .95, title.c_str()); return fCanvas; } + +/////////////////////////////////////////////// +/// \brief A method that creates a canvas where the mirror optics properties are drawn. +/// It generates two plots, on the left the reflectivity as a function of the angle, and +/// on the right the reflectivity as a function of the energy. +/// +/// The first argument is a string where we may specify the energies and angles to be plotted +/// against the angle and energy. Between square brackets [ ] we define the values that will +/// be plotted, while between parenthesis we define the range of the x-axis to be plotted. +/// +/// For example the following definition will produce a plot with 3 energies (2,5,10) keV +/// as a function of the angle, in the range 0 to 45 degrees, and a second plot with 4 +/// angles (1,2,5,15) degrees in the energy range 0 to 15keV. +/// +/// ``` +/// [2,5,10](0,45):[1,2,5,15](0,15) +/// ``` +/// +/// Optionally we may define the legend position for each plot by adding {x1,y1,x2,y2}, where +/// x1,y1,x2,y2 are values between 0 and 1 defining the region box that the legend will +/// occupate at the pad. See also ROOT TLegend documentation. +/// +/// We may add the position to any of the two plots, or to both of them. Here we define it +/// for the second plot to be placed at the bottom-left. +/// +/// ``` +/// [2,5,10](0,45):[1,2,5,15](0,15){0.15,0.15,0.45,0.45} +/// ``` +/// +/// If we define legend positions only for the first plot, those values will also be used for +/// the second plot. If no legend positions are given, the default positions will be used. +/// +/// Two additional double arguments allow to set the low y-range of each of the plots. +/// +TPad* TRestAxionOpticsMirror::DrawOpticsPropertiesLinear(std::string options, Double_t lowRange, + Double_t lowRange2) { + if (fReflectivityTable.size() == 0) LoadTables(); + + std::vector optList = TRestTools::GetOptions(options); + + if (optList.size() == 0) + optList = TRestTools::GetOptions( + "[2,4,6](0,3.5){0.65,0.68,0.85,0.88}:[0.25,0.5,0.75](0,10){0.2,0.2,0.45,0.45}"); + + if (optList.size() != 2) { + RESTError << "TRestAxionOpticsMirror::DrawOpticsProperties. Wrong arguments!" << RESTendl; + return fCanvas; + } + + std::vector energies = StringToElements(optList[0], "[", ",", "]"); + std::vector aRange = StringToElements(optList[0], "(", ",", ")"); + std::vector eLegendCoords = StringToElements(optList[0], "{", ",", "}"); + + std::vector angles = StringToElements(optList[1], "[", ",", "]"); + std::vector eRange = StringToElements(optList[1], "(", ",", ")"); + std::vector aLegendCoords = StringToElements(optList[1], "{", ",", "}"); + + if (eRange[0] < 0.03) eRange[0] = 0.03; + if (eRange[1] > 15) eRange[1] = 15; + if (aRange[0] < 0.0) aRange[0] = 0.0; + if (aRange[1] > 9) aRange[1] = 9; + + TPad* pad = new TPad("pad", "This is pad", 0.01, 0.01, 0.99, 0.99); + pad->Divide(1, 2); + pad->Draw(); + + ////// Drawing reflectivity versus angle + pad->cd(1); + pad->cd(1)->SetRightMargin(0.05); + pad->cd(1)->SetLeftMargin(0.15); + pad->cd(1)->SetBottomMargin(0.14); + pad->cd(1)->SetTopMargin(0.05); + + std::vector ref_vs_ang_graph; + + for (unsigned int n = 0; n < energies.size(); n++) { + string grname = "gr" + IntegerToString(n); + TGraph* gr = new TGraph(); + gr->SetName(grname.c_str()); + for (double a = aRange[0]; a <= aRange[1]; a += (aRange[1] - aRange[0]) / 10000.) { + gr->SetPoint(gr->GetN(), a, GetReflectivity(a, energies[n])); + } + gr->SetLineColor(49 - n * 3); + gr->SetLineWidth(2); + ref_vs_ang_graph.push_back(gr); + } + + ref_vs_ang_graph[0]->GetXaxis()->SetLimits(aRange[0], aRange[1]); + ref_vs_ang_graph[0]->GetHistogram()->SetMaximum(1); + ref_vs_ang_graph[0]->GetHistogram()->SetMinimum(lowRange); + + ref_vs_ang_graph[0]->GetXaxis()->SetTitle("Angle [degrees]"); + ref_vs_ang_graph[0]->GetXaxis()->SetTitleSize(0.06); + ref_vs_ang_graph[0]->GetXaxis()->SetLabelSize(0.06); + ref_vs_ang_graph[0]->GetYaxis()->SetTitle("Reflectivity"); + ref_vs_ang_graph[0]->GetYaxis()->SetTitleOffset(1.2); + ref_vs_ang_graph[0]->GetYaxis()->SetTitleSize(0.06); + ref_vs_ang_graph[0]->GetYaxis()->SetLabelSize(0.06); + ref_vs_ang_graph[0]->Draw("AL"); + for (unsigned int n = 1; n < energies.size(); n++) ref_vs_ang_graph[n]->Draw("L"); + + Double_t lx1 = 0.6, ly1 = 0.75, lx2 = 0.9, ly2 = 0.95; + if (eLegendCoords.size() > 0) { + lx1 = eLegendCoords[0]; + ly1 = eLegendCoords[1]; + lx2 = eLegendCoords[2]; + ly2 = eLegendCoords[3]; + } + TLegend* legend = new TLegend(lx1, ly1, lx2, ly2); + + legend->SetTextSize(0.03); + legend->SetHeader("Energies", "C"); // option "C" allows to center the header + for (unsigned int n = 0; n < energies.size(); n++) { + std::string lname = "gr" + IntegerToString(n); + std::string ltitle = DoubleToString(energies[n], "%3.2lf") + " keV"; + + legend->AddEntry(lname.c_str(), ltitle.c_str(), "l"); + } + legend->Draw(); + + ////// Drawing reflectivity versus energy + pad->cd(2); + pad->cd(2)->SetRightMargin(0.05); + pad->cd(2)->SetLeftMargin(0.16); + pad->cd(2)->SetBottomMargin(0.14); + pad->cd(2)->SetTopMargin(0.05); + + std::vector ref_vs_en_graph; + + for (unsigned int n = 0; n < angles.size(); n++) { + string grname = "agr" + IntegerToString(n); + TGraph* gr = new TGraph(); + gr->SetName(grname.c_str()); + for (double e = eRange[0]; e <= eRange[1]; e += (eRange[1] - eRange[0]) / 10000.) { + gr->SetPoint(gr->GetN(), e, GetReflectivity(angles[n], e)); + } + gr->SetLineColor(49 - n * 3); + gr->SetLineWidth(2); + ref_vs_en_graph.push_back(gr); + } + + ref_vs_en_graph[0]->GetXaxis()->SetLimits(eRange[0], eRange[1]); + ref_vs_en_graph[0]->GetHistogram()->SetMaximum(1); + ref_vs_en_graph[0]->GetHistogram()->SetMinimum(lowRange2); + + ref_vs_en_graph[0]->GetXaxis()->SetTitle("Energy [keV]"); + ref_vs_en_graph[0]->GetXaxis()->SetTitleSize(0.06); + ref_vs_en_graph[0]->GetXaxis()->SetLabelSize(0.06); + ref_vs_en_graph[0]->GetYaxis()->SetTitle("Reflectivity"); + ref_vs_en_graph[0]->GetYaxis()->SetTitleOffset(1.2); + ref_vs_en_graph[0]->GetYaxis()->SetTitleSize(0.06); + ref_vs_en_graph[0]->GetYaxis()->SetLabelSize(0.06); + ref_vs_en_graph[0]->Draw("AL"); + for (unsigned int n = 1; n < angles.size(); n++) ref_vs_en_graph[n]->Draw("L"); + + if (aLegendCoords.size() > 0) { + lx1 = aLegendCoords[0]; + ly1 = aLegendCoords[1]; + lx2 = aLegendCoords[2]; + ly2 = aLegendCoords[3]; + } + TLegend* legendA = new TLegend(lx1, ly1, lx2, ly2); + legendA->SetTextSize(0.03); + legendA->SetHeader("Angles", "C"); // option "C" allows to center the header + for (unsigned int n = 0; n < angles.size(); n++) { + std::string lname = "agr" + IntegerToString(n); + std::string ltitle = DoubleToString(angles[n], "%3.2lf") + " degrees"; + + legendA->AddEntry(lname.c_str(), ltitle.c_str(), "l"); + } + legendA->Draw(); + + /// Drawing a main title + return pad; +} diff --git a/src/TRestAxionSolarQCDFlux.cxx b/src/TRestAxionSolarQCDFlux.cxx index f96a0eba..195ff1d1 100644 --- a/src/TRestAxionSolarQCDFlux.cxx +++ b/src/TRestAxionSolarQCDFlux.cxx @@ -542,9 +542,10 @@ TH1F* TRestAxionSolarQCDFlux::GetTotalSpectrum() { } fTotalHist = new TH1F("fTotalHist", "", 20000, 0, 20); - for (int n = 0; n < hc->GetNbinsX(); n++) { + for (int n = 0; n < hc->GetNbinsX() - 1; n++) { for (int m = 0; m < 100; m++) { - fTotalHist->SetBinContent(n * 100 + 1 + m, hc->GetBinContent(n + 1)); + fTotalHist->SetBinContent(n * 100 + 1 + m, (1 - 0.01 * (Double_t)m) * hc->GetBinContent(n + 1) + + 0.01 * (Double_t)m * hc->GetBinContent(n + 2)); } } diff --git a/src/TRestAxionTrueWolterOptics.cxx b/src/TRestAxionTrueWolterOptics.cxx index 3daa32e0..ed5df7a9 100644 --- a/src/TRestAxionTrueWolterOptics.cxx +++ b/src/TRestAxionTrueWolterOptics.cxx @@ -430,17 +430,17 @@ TPad* TRestAxionTrueWolterOptics::DrawMirrors() { gr->SetPoint(2, lX, fR5[mirror]); gr->GetXaxis()->SetLimits(-3.5 * lX, 3.5 * lX); - gr->GetHistogram()->SetMaximum(fR1.back() * 1.15); + gr->GetHistogram()->SetMaximum(fR1.back() * 1.05); gr->GetHistogram()->SetMinimum(fR1.front() * 0.8); gr->GetXaxis()->SetTitle("Z [mm]"); - gr->GetXaxis()->SetTitleSize(0.04); - gr->GetXaxis()->SetLabelSize(0.04); + gr->GetXaxis()->SetTitleSize(0.05); + gr->GetXaxis()->SetLabelSize(0.05); gr->GetXaxis()->SetNdivisions(5); gr->GetYaxis()->SetTitle("R [mm]"); - gr->GetYaxis()->SetTitleOffset(1.4); - gr->GetYaxis()->SetTitleSize(0.04); - gr->GetYaxis()->SetLabelSize(0.04); + gr->GetYaxis()->SetTitleOffset(1.8); + gr->GetYaxis()->SetTitleSize(0.05); + gr->GetYaxis()->SetLabelSize(0.05); gr->SetLineWidth(6 * fThickness[mirror]); gr->SetLineColor(20 + mirror % 20); if (mirror == 0) diff --git a/src/TRestAxionWolterOptics.cxx b/src/TRestAxionWolterOptics.cxx index 27e84d81..d38f530f 100644 --- a/src/TRestAxionWolterOptics.cxx +++ b/src/TRestAxionWolterOptics.cxx @@ -426,7 +426,7 @@ TPad* TRestAxionWolterOptics::DrawMirrors() { gr->SetPoint(2, lX, fR5[mirror]); gr->GetXaxis()->SetLimits(-3.5 * lX, 3.5 * lX); - gr->GetHistogram()->SetMaximum(fR1.back() * 1.15); + gr->GetHistogram()->SetMaximum(fR1.back() * 1.05); gr->GetHistogram()->SetMinimum(fR1.front() * 0.8); gr->GetXaxis()->SetTitle("Z [mm]");