diff --git a/inc/TRestAxionSolarFlux.h b/inc/TRestAxionSolarFlux.h index 9a917211..6c740852 100644 --- a/inc/TRestAxionSolarFlux.h +++ b/inc/TRestAxionSolarFlux.h @@ -65,10 +65,13 @@ class TRestAxionSolarFlux : public TRestMetadata { void Initialize(); /// It is required in order to load solar flux tables into memory for specific mass - void InitializeMass( Double_t mass ) { SetMass(mass); Initialize(); } - + void InitializeMass(Double_t mass) { + SetMass(mass); + Initialize(); + } + /// Set mass and reinitialise - void SetMass( Double_t m ) { fMass = m; } //Initialize(); } + void SetMass(Double_t m) { fMass = m; } // Initialize(); } /// It returns the integrated flux at earth in cm-2 s-1 for the given energy range virtual Double_t IntegrateFluxInRange(TVector2 eRange = TVector2(-1, -1)) = 0; diff --git a/inc/TRestAxionSolarHiddenPhotonFlux.h b/inc/TRestAxionSolarHiddenPhotonFlux.h index bdc44040..5a6b0398 100644 --- a/inc/TRestAxionSolarHiddenPhotonFlux.h +++ b/inc/TRestAxionSolarHiddenPhotonFlux.h @@ -73,7 +73,7 @@ class TRestAxionSolarHiddenPhotonFlux : public TRestAxionSolarFlux { void LoadMonoChromaticFluxTable(); void IntegrateSolarFluxes(); - public: + public: /// It returns true if continuum flux spectra was loaded Bool_t isSolarTableLoaded() { return fFluxTable.size() > 0; } diff --git a/inc/TRestAxionSolarQCDFlux.h b/inc/TRestAxionSolarQCDFlux.h index 10db402e..fb99ebb0 100644 --- a/inc/TRestAxionSolarQCDFlux.h +++ b/inc/TRestAxionSolarQCDFlux.h @@ -93,9 +93,7 @@ class TRestAxionSolarQCDFlux : public TRestAxionSolarFlux { Bool_t LoadTables() override; /// It returns the total integrated flux at earth in cm-2 s-1 - Double_t GetTotalFlux() override { - return fTotalContinuumFlux + fTotalMonochromaticFlux; - } + Double_t GetTotalFlux() override { return fTotalContinuumFlux + fTotalMonochromaticFlux; } /// It returns an energy integrated spectrum in cm-2 s-1 keV-1 TH1D* GetEnergySpectrum() override { return GetTotalSpectrum(); } diff --git a/pipeline/metadata/solarFlux/solarPlotHiddenPhoton.py b/pipeline/metadata/solarFlux/solarPlotHiddenPhoton.py index c07ac580..5b83b758 100755 --- a/pipeline/metadata/solarFlux/solarPlotHiddenPhoton.py +++ b/pipeline/metadata/solarFlux/solarPlotHiddenPhoton.py @@ -44,12 +44,10 @@ parser.add_argument( "--N", dest="samples", type=int, help="The number of generated particles" ) -parser.add_argument( - "--m", dest="mass", type=float, help="Hidden photon mass [eV]" -) +parser.add_argument("--m", dest="mass", type=float, help="Hidden photon mass [eV]") args = parser.parse_args() -mass = 10. # eV +mass = 10.0 # eV if args.mass != None: mass = args.mass @@ -96,7 +94,7 @@ comb_spt = TH2D("comb_spt", "Energy versus solar radius", 20000, 0, 20, 100, 0, 1) for i in range(samples): x = combinedFlux.GetRandomEnergyAndRadius((-1, -1)) - #print(x) + # print(x) comb_spt.Fill(x[0], x[1]) rnd = TRandom3(0) @@ -164,4 +162,4 @@ c1.Print(outfname) print("Generated file : " + outfname) -#exit(0) +# exit(0) diff --git a/src/TRestAxionSolarFlux.cxx b/src/TRestAxionSolarFlux.cxx index f82d9e34..62549922 100644 --- a/src/TRestAxionSolarFlux.cxx +++ b/src/TRestAxionSolarFlux.cxx @@ -134,12 +134,11 @@ void TRestAxionSolarFlux::Initialize() { fSeed = fRandom->TRandom::GetSeed(); } - /////////////////////////////////////////////// /// \brief Initialization of TRestAxionSolarFlux members with specific mass /// -//void TRestAxionSolarFlux::InitializeMass( Double_t mass ) { SetMass(mass); RESTMetadata << GetMass() << RESTendl; } // SetMass calls Initialize - +// void TRestAxionSolarFlux::InitializeMass( Double_t mass ) { SetMass(mass); RESTMetadata << GetMass() << +// RESTendl; } // SetMass calls Initialize /////////////////////////////////////////////// /// \brief It builds a histogram using the contents of the .flux file given diff --git a/src/TRestAxionSolarHiddenPhotonFlux.cxx b/src/TRestAxionSolarHiddenPhotonFlux.cxx index 0c213f20..f824fd5c 100644 --- a/src/TRestAxionSolarHiddenPhotonFlux.cxx +++ b/src/TRestAxionSolarHiddenPhotonFlux.cxx @@ -184,23 +184,28 @@ TRestAxionSolarHiddenPhotonFlux::~TRestAxionSolarHiddenPhotonFlux() {} /// inside the metadata members, and calculate the solar flux for a given m. /// Bool_t TRestAxionSolarHiddenPhotonFlux::LoadTables() { - - if (GetMass() <= 0 ) { - RESTWarning << "TRestAxionSolarHiddenPhotonFlux::LoadTables - hidden photon mass not yet defined" << RESTendl; - RESTWarning << "TRestAxionSolarHiddenPhotonFlux::LoadTables - mass given as " << GetMass() << " eV" << RESTendl; - return false; + if (GetMass() <= 0) { + RESTWarning << "TRestAxionSolarHiddenPhotonFlux::LoadTables - hidden photon mass not yet defined" + << RESTendl; + RESTWarning << "TRestAxionSolarHiddenPhotonFlux::LoadTables - mass given as " << GetMass() << " eV" + << RESTendl; + return false; + } + if (fFluxDataFile == "") { + RESTWarning << "TRestAxionSolarHiddenPhotonFlux::LoadTables - flux table not found!!\n " + << fFluxDataFile << RESTendl; + return false; } - if ( fFluxDataFile == "" ){ - RESTWarning << "TRestAxionSolarHiddenPhotonFlux::LoadTables - flux table not found!!\n " << fFluxDataFile << RESTendl; - return false; + if (fWidthDataFile == "") { + RESTWarning << "TRestAxionSolarHiddenPhotonFlux::LoadTables - width table not found!!\n " + << fWidthDataFile << RESTendl; + return false; } - if ( fWidthDataFile == "") { - RESTWarning << "TRestAxionSolarHiddenPhotonFlux::LoadTables - width table not found!!\n " << fWidthDataFile << RESTendl; - return false; - } - if ( fPlasmaFreqDataFile == "" ) { - RESTWarning << "TRestAxionSolarHiddenPhotonFlux::LoadTables - plasma frequency table not found!!\n " << fPlasmaFreqDataFile << RESTendl; - return false; + if (fPlasmaFreqDataFile == "") { + RESTWarning + << "TRestAxionSolarHiddenPhotonFlux::LoadTables - plasma frequency table not found!!\n " + << fPlasmaFreqDataFile << RESTendl; + return false; } LoadContinuumFluxTable(); @@ -248,7 +253,9 @@ void TRestAxionSolarHiddenPhotonFlux::LoadContinuumFluxTable() { for (unsigned int n = 0; n < fluxTable.size(); n++) { TH1D* h = new TH1D(Form("%s_ContinuumFluxAtRadius%d", GetName(), n), "", 200, 0, 20); - for (unsigned int m = 0; m < fluxTable[n].size(); m++) { h->SetBinContent(m + 1, fluxTable[n][m]); } + for (unsigned int m = 0; m < fluxTable[n].size(); m++) { + h->SetBinContent(m + 1, fluxTable[n][m]); + } fContinuumTable.push_back(h); } } @@ -277,8 +284,9 @@ void TRestAxionSolarHiddenPhotonFlux::LoadWidthTable() { } TRestTools::ReadBinaryTable(fullPathName, fluxTable); - //RESTMetadata << "Width table rows / columns: " << fluxTable.size() << " " << fluxTable[0].size() << RESTendl; - + // RESTMetadata << "Width table rows / columns: " << fluxTable.size() << " " << fluxTable[0].size() << + // RESTendl; + if (fluxTable.size() != 1000 || fluxTable[0].size() != 200) { fluxTable.clear(); RESTError << "LoadWidthTable. The table does not contain the right number of rows or columns" @@ -288,7 +296,9 @@ void TRestAxionSolarHiddenPhotonFlux::LoadWidthTable() { for (unsigned int n = 0; n < fluxTable.size(); n++) { TH1D* h = new TH1D(Form("%s_ResonanceWidthAtRadius%d", GetName(), n), "", 200, 0, 20); - for (unsigned int m = 0; m < fluxTable[n].size(); m++) { h->SetBinContent(m + 1, fluxTable[n][m]); } + for (unsigned int m = 0; m < fluxTable[n].size(); m++) { + h->SetBinContent(m + 1, fluxTable[n][m]); + } fWidthTable.push_back(h); } } @@ -318,7 +328,7 @@ void TRestAxionSolarHiddenPhotonFlux::LoadPlasmaFreqTable() { } TRestTools::ReadBinaryTable(fullPathName, fluxTable); - + if (fluxTable.size() != 1000 || fluxTable[0].size() != 1) { fluxTable.clear(); RESTError << "LoadPlasmaFreqTable. The table does not contain the right number of rows or columns" @@ -328,7 +338,9 @@ void TRestAxionSolarHiddenPhotonFlux::LoadPlasmaFreqTable() { for (unsigned int n = 0; n < fluxTable.size(); n++) { TH1D* h = new TH1D(Form("%s_PlasmaFreqAtRadius%d", GetName(), n), "", 1, 0, 20); - for (unsigned int m = 0; m < fluxTable[n].size(); m++) { h->SetBinContent(m + 1, fluxTable[n][m]); } + for (unsigned int m = 0; m < fluxTable[n].size(); m++) { + h->SetBinContent(m + 1, fluxTable[n][m]); + } fPlasmaFreqTable.push_back(h); } } @@ -341,20 +353,21 @@ void TRestAxionSolarHiddenPhotonFlux::CalculateSolarFlux() { if (GetMass() == 0) { RESTError << "CalculateSolarFlux. The hidden photon mass is set to zero!" << RESTendl; return; - } - if (fContinuumTable.size() == 0 ) { - RESTError << "TRestAxionSolarHiddenPhotonFlux::CalculateSolarFlux - empty flux table!" << RESTendl; - return; } - if (fPlasmaFreqTable.size() == 0 ) { - RESTError << "TRestAxionSolarHiddenPhotonFlux::CalculateSolarFlux - empty plasma freq table!" << RESTendl; - return; + if (fContinuumTable.size() == 0) { + RESTError << "TRestAxionSolarHiddenPhotonFlux::CalculateSolarFlux - empty flux table!" << RESTendl; + return; } - if (fWidthTable.size() == 0 ) { - RESTError << "TRestAxionSolarHiddenPhotonFlux::CalculateSolarFlux - empty width table!" << RESTendl; - return; + if (fPlasmaFreqTable.size() == 0) { + RESTError << "TRestAxionSolarHiddenPhotonFlux::CalculateSolarFlux - empty plasma freq table!" + << RESTendl; + return; + } + if (fWidthTable.size() == 0) { + RESTError << "TRestAxionSolarHiddenPhotonFlux::CalculateSolarFlux - empty width table!" << RESTendl; + return; } - + Double_t mass = GetMass(); cout << mass << endl; for (unsigned int n = 0; n < fContinuumTable.size(); n++) { @@ -363,19 +376,21 @@ void TRestAxionSolarHiddenPhotonFlux::CalculateSolarFlux() { Double_t wp = fPlasmaFreqTable[n]->GetBinContent(1); TH1D* hMass = new TH1D(Form("%s_hMass%d", GetName(), n), "hMass", 200, 0, 20); TH1D* hWg2 = (TH1D*)fWidthTable[n]->Clone(); - hWg2->Multiply(hWg2); // (w G)^2 + hWg2->Multiply(hWg2); // (w G)^2 - for ( unsigned int c = 0; c < 200; c++ ) { - Double_t wG = fWidthTable[n]->GetBinContent(c+1); - hMass->SetBinContent( c+1, pow(mass,-4) * ( pow( pow(mass,2) - pow(wp,2) , 2 )));// + pow(wG,2) ) ); // m2 - } - - hMass->Add(hWg2); // (m2 - wp2)^2 + (w G)^2 + for (unsigned int c = 0; c < 200; c++) { + Double_t wG = fWidthTable[n]->GetBinContent(c + 1); + hMass->SetBinContent( + c + 1, pow(mass, -4) * (pow(pow(mass, 2) - pow(wp, 2), 2))); // + pow(wG,2) ) ); // m2 + } - TH1D* h = (TH1D*)fWidthTable[n]->Clone(); // wG - h->Multiply(fContinuumTable[n]); // wG * flux - h->Divide(hMass); // wG * flux / ( (m2 - wp2)^2 + (w G)^2 ) - //h->Scale(pow(mass, 4)); // m4 * wG * flux / ( (m2 - wp2)^2 + (w G)^2 ) + hMass->Add(hWg2); // (m2 - wp2)^2 + (w G)^2 + + TH1D* h = (TH1D*)fWidthTable[n]->Clone(); // wG + h->Multiply(fContinuumTable[n]); // wG * flux + h->Divide(hMass); // wG * flux / ( (m2 - wp2)^2 + (w G)^2 ) + // h->Scale(pow(mass, 4)); // m4 * wG * flux / ( (m2 - wp2)^2 + (w G)^2 + // ) fFluxTable.push_back(h); } @@ -477,9 +492,11 @@ Double_t TRestAxionSolarHiddenPhotonFlux::IntegrateFluxInRange(TVector2 eRange) /// flux distributions defined inside the solar tables loaded in the class /// std::pair TRestAxionSolarHiddenPhotonFlux::GetRandomEnergyAndRadius(TVector2 eRange) { - std::pair result = {0, 0}; - if (!AreTablesLoaded()) { RESTWarning << "Tables not loaded!!" << RESTendl; return result; } + if (!AreTablesLoaded()) { + RESTWarning << "Tables not loaded!!" << RESTendl; + return result; + } Double_t rnd = fRandom->Rndm(); for (unsigned int r = 0; r < fFluxTableIntegrals.size(); r++) { if (rnd < fFluxTableIntegrals[r]) { @@ -529,9 +546,9 @@ void TRestAxionSolarHiddenPhotonFlux::PrintIntegratedRingFlux() { void TRestAxionSolarHiddenPhotonFlux::PrintMetadata() { TRestAxionSolarFlux::PrintMetadata(); - RESTMetadata << " - Solar hidden photon datafile (flux) : " << fFluxDataFile << RESTendl; - RESTMetadata << " - Solar hidden photon datafile (width) : " << fWidthDataFile << RESTendl; - RESTMetadata << " - Solar hidden photon datafile (plasma freq) : " << fPlasmaFreqDataFile << RESTendl; + RESTMetadata << " - Solar hidden photon datafile (flux) : " << fFluxDataFile << RESTendl; + RESTMetadata << " - Solar hidden photon datafile (width) : " << fWidthDataFile << RESTendl; + RESTMetadata << " - Solar hidden photon datafile (plasma freq) : " << fPlasmaFreqDataFile << RESTendl; RESTMetadata << "-------" << RESTendl; RESTMetadata << " - Total continuum flux : " << fTotalContinuumFlux << " cm-2 s-1" << RESTendl; RESTMetadata << "++++++++++++++++++" << RESTendl; @@ -610,4 +627,3 @@ TCanvas* TRestAxionSolarHiddenPhotonFlux::DrawSolarFlux() { return fCanvas; } -