Raw update

rest-for-physics · May 29, 2024 · c1618b9 · c1618b9
1 parent ea89509
commit c1618b9
Show file tree

Hide file tree

Showing 6 changed files with 140 additions and 10 deletions.
diff --git a/examples/03.sensitivity/BabyIAXO.rml b/examples/03.sensitivity/BabyIAXO.rml
@@ -2,7 +2,12 @@
 <BabyIAXO>
 
 	<TRestSensitivity name="VacuumPhase" strategy="nodes" >
+<<<<<<< Updated upstream
 		<TRestExperiment name="Vacuum" exposureTime="3*300*12hr">
+=======
+
+		<TRestExperiment name="Vacuum" exposureTime="3*300*12hr" useAverage="true">
+>>>>>>> Stashed changes
 			<!-- Background -->
 			<TRestComponentFormula name="Flat7" nature="background" >
 				<parameter name="formulaUnits" value="keV^-1" />
@@ -15,7 +20,7 @@
 			<TRestAxionHelioscopeSignal name="Vacuum" nature="signal" parameter="ma"
 							conversionType="IAXO" bores="2"
 							magnetRadius="35cm" magnetLength="10m" magnetStrength="2T"
-							opticsEfficiency="0.35" windowEfficiency="0.8">
+				opticsEfficiency="0.35" windowEfficiency="0.8" >
 
 				<cVariable name="energy" range="(0,10)keV" bins="20" />
 
@@ -26,15 +31,17 @@
 
 				<TRestAxionSolarQCDFlux name="LennertHoofPrimakoff" couplingType="g_ag" couplingStrength="1.e-10"
 										fluxDataFile="Primakoff_LennertHoof_202203.dat" seed="137" />
-
-				<TRestResponse name="XenonNeon" variable="energy" filename="XenonNeon_50Pct_1.4bar.N150f"/>
 			</TRestAxionHelioscopeSignal>
 		</TRestExperiment>
 
 	</TRestSensitivity>
 
 	<TRestSensitivity name="CombinedPhase" strategy="nodes" >
+<<<<<<< Updated upstream
 		<TRestExperiment name="Vacuum" exposureTime="3*300*12hr">
+=======
+		<TRestExperiment name="Vacuum" exposureTime="3*300*12hr" useAverage="true">
+>>>>>>> Stashed changes
 			<!-- Background -->
 			<TRestComponentFormula name="Flat7" nature="background" >
 				<parameter name="formulaUnits" value="keV^-1" />
@@ -54,11 +61,22 @@
 				<TRestAxionSolarQCDFlux name="LennertHoofPrimakoff" couplingType="g_ag" couplingStrength="1.e-10"
 										fluxDataFile="Primakoff_LennertHoof_202203.dat" seed="137" />
 
+
+				<parameter name="firstParameterValue" value="0.001" />
+				<parameter name="lastParameterValue" value="10" />
+				<parameter name="stepParameterValue" value="1.02" />
+				<parameter name="exponential" value="true" />
+
 			</TRestAxionHelioscopeSignal>
 		</TRestExperiment>
 
+<<<<<<< Updated upstream
 		<TRestExperimentList name="GasPhase" exposureTime="12*12hr"
 			componentPattern="output/SignalsBabyIAXO.root" experimentsFile="output/BabyIAXO.settings" >
+=======
+		<TRestExperimentList name="GasPhase" exposureTime="0"
+			componentPattern="output/SignalsBabyIAXO.root" experimentsFile="output/BabyIAXO.settings" useAverage="true">
+>>>>>>> Stashed changes
 
 			<TRestComponentFormula name="Flat7" nature="background" >
 				<parameter name="formulaUnits" value="keV^-1" />
@@ -67,7 +85,6 @@
 				<formula name="bck" expression="2*0.8*0.8*TMath::Pi()*1e-07" />
 			</TRestComponentFormula>
 		</TRestExperimentList>
-
 	</TRestSensitivity>
 
 	<TRestAxionHelioscopeSignal name="GasSignal" nature="signal" parameter="ma"
@@ -77,14 +94,17 @@
 
 		<cVariable name="energy" range="(0,10)keV" bins="20" />
 
+		<parameter name="firstParameterValue" value="0.001" />
+		<parameter name="lastParameterValue" value="10" />
+		<parameter name="stepParameterValue" value="1.02" />
+		<parameter name="exponential" value="true" />
+
 		<TRestAxionSolarQCDFlux name="LennertHoofPrimakoff" couplingType="g_ag" couplingStrength="1.e-10"
 								fluxDataFile="Primakoff_LennertHoof_202203.dat" seed="137" />
 
 		<TRestAxionBufferGas name="helium" >
 			<gas name="He" density="0.0025e-6g/cm^3"/>
 		</TRestAxionBufferGas>
-
-		<TRestResponse name="XenonNeon" variable="energy" filename="XenonNeon_50Pct_1.4bar.N150f"/>
 	</TRestAxionHelioscopeSignal>
 
 </BabyIAXO>
diff --git a/examples/03.sensitivity/GenerateSignalComponents.C b/examples/03.sensitivity/GenerateSignalComponents.C
@@ -1,9 +1,22 @@
+<<<<<<< Updated upstream
 Int_t GenerateSignalComponents(std::string rmlFile, std::string name) {
     TRestAxionField field;
     std::vector<std::pair<Double_t, Double_t>> scanSteps =
         field.GetMassDensityScanning("He", 0.25, 20);  // Up to 0.25 eV
 
     TRestAxionHelioscopeSignal gasPhase(rmlFile.c_str(), name.c_str());
+=======
+Double_t Eo = 0.5; //keV
+Double_t Ef = 10; //keV
+
+Int_t GenerateSignalComponents(std::string rmlFile, std::string name, Double_t totalExposureTime=3*300*12*3600, size_t skipSteps = 1)
+{
+	TRestAxionField field;
+	std::vector<std::pair<Double_t, Double_t>> scanSteps = field.GetMassDensityScanning( "He", 0.25, 20 ); // Up to 0.25 eV
+
+	TRestAxionHelioscopeSignal gasPhase( rmlFile.c_str(), name.c_str());
+	//gasPhase.RegenerateParametricNodes( startMass, upperMass, step, true );
+>>>>>>> Stashed changes
 
     std::filesystem::path filePath = rmlFile;
     std::filesystem::path newExtension = ".settings";
@@ -15,6 +28,7 @@ Int_t GenerateSignalComponents(std::string rmlFile, std::string name) {
     settingsFile = "output/" + settingsFile;
     componentsFile = "output/" + componentsFile;
 
+<<<<<<< Updated upstream
     FILE* g = fopen(settingsFile.c_str(), "wt");
     TFile* f = TFile::Open((TString)componentsFile, "RECREATE");
     for (size_t n = 0; n < scanSteps.size(); n++) {
@@ -27,8 +41,51 @@ Int_t GenerateSignalComponents(std::string rmlFile, std::string name) {
     fclose(g);
 
     f->Close();
+=======
+	TFile *f = TFile::Open( (TString) componentsFile, "RECREATE" );
+	std::vector <Double_t> exposureTimes;
+	for( size_t n = skipSteps; n < scanSteps.size(); n++ )
+	{
+		Double_t mass = scanSteps[n].first;
+		std::cout << "Mass: " << mass << std::endl;
+		Double_t density = scanSteps[n].second;
+		std::cout << "Setting density: " << density << std::endl;
+
+		gasPhase.SetName( "P" + (TString) IntegerToString(n+1) );
+		gasPhase.GetGas()->SetGasDensity("He", density);
+		gasPhase.RegenerateHistograms();
+
+		Double_t nGamma = gasPhase.GetSignalRate(mass, Eo, Ef);
+		std::cout << "nGamma: " << nGamma << std::endl;
+		std::cout << "Mass: " << mass << std::endl;
+		Double_t ksvzFactor = 3.75523 * mass;
+		std::cout << "Log(20)/ksvzFactor : " << TMath::Log(20)/TMath::Power( ksvzFactor, 4 ) << std::endl;
+		Double_t exposureTime = TMath::Log(20.) / TMath::Power(ksvzFactor, 4) / nGamma;
+		exposureTimes.push_back( exposureTime );
+		std::cout << "Exposure time: " << exposureTime << std::endl;
+
+		gasPhase.Write( "P" + (TString) IntegerToString(n+1) );
+	}
+	f->Close();
+>>>>>>> Stashed changes
 
     std::cout << "Generated signals file: " << componentsFile << std::endl;
 
+<<<<<<< Updated upstream
     return 0;
+=======
+	Double_t generatedExposureTime = 0;
+	for( size_t n = skipSteps; n < scanSteps.size(); n++ )
+		generatedExposureTime += exposureTimes[n-skipSteps];
+
+
+	FILE *g = fopen(settingsFile.c_str(), "wt");
+	for( size_t n = skipSteps; n < scanSteps.size(); n++ )
+		fprintf( g, "%lf\tP%d\n", totalExposureTime * exposureTimes[n-skipSteps] / generatedExposureTime, (int) (n+1) );
+	fclose(g);
+
+	std::cout << "Generated settings file: " << settingsFile << std::endl;
+
+	return 0;
+>>>>>>> Stashed changes
 }
diff --git a/examples/03.sensitivity/IAXO.rml b/examples/03.sensitivity/IAXO.rml
@@ -2,8 +2,13 @@
 <BabyIAXO>
 
 	<TRestSensitivity name="VacuumPhase" strategy="nodes" >
+<<<<<<< Updated upstream
 
 		<TRestExperiment name="Vacuum" exposureTime="6*300*12hr">
+=======
+
+		<TRestExperiment name="Vacuum" exposureTime="12*300*12hr" useAverage="true">
+>>>>>>> Stashed changes
 			<!-- Background -->
 			<TRestComponentFormula name="Flat7" nature="background" >
 				<parameter name="formulaUnits" value="keV^-1" />
@@ -13,17 +18,26 @@
 			</TRestComponentFormula>
 
 			<!-- Signal -->
+<<<<<<< Updated upstream
 			<TRestAxionHelioscopeSignal name="Vacuum" nature="signal"
 							conversionType="IAXO" bores="8"
+=======
+			<TRestAxionHelioscopeSignal name="Vacuum" nature="signal" parameter="ma"
+							conversionType="IAXO" bores="8" 
+>>>>>>> Stashed changes
 							magnetRadius="30cm" magnetLength="20m" magnetStrength="2.5T"
 							opticsEfficiency="0.7" windowEfficiency="0.8">
 
 				<cVariable name="energy" range="(0,10)keV" bins="20" />
 
+				<parameter name="firstParameterValue" value="0.001" />
+				<parameter name="lastParameterValue" value="10" />
+				<parameter name="stepParameterValue" value="1.02" />
+				<parameter name="exponential" value="true" />
+
 				<TRestAxionSolarQCDFlux name="LennertHoofPrimakoff" couplingType="g_ag" couplingStrength="1.e-10"
 										fluxDataFile="Primakoff_LennertHoof_202203.dat" seed="137" />
 
-				<TRestResponse name="XenonNeon" variable="energy" filename="XenonNeon_50Pct_1.4bar.N150f"/>
 			</TRestAxionHelioscopeSignal>
 		</TRestExperiment>
 
@@ -55,8 +69,13 @@
 			</TRestAxionHelioscopeSignal>
 		</TRestExperiment>
 
+<<<<<<< Updated upstream
 		<TRestExperimentList name="GasPhase" exposureTime="24*12hr"
 			componentPattern="SignalComponents.root" experimentsFile="GasPhase.settings" >
+=======
+		<TRestExperimentList name="GasPhase" exposureTime="24*12hr" 
+			componentPattern="output/SignalsIAXO.root" experimentsFile="output/IAXO.settings" >
+>>>>>>> Stashed changes
 
 			<TRestComponentFormula name="Flat7" nature="background" >
 				<parameter name="formulaUnits" value="keV^-1" />

diff --git a/inc/TRestAxionHelioscopeSignal.h b/inc/TRestAxionHelioscopeSignal.h
@@ -69,8 +69,10 @@ class TRestAxionHelioscopeSignal : public TRestComponent {
 
    public:
     Double_t GetSignalRate(std::vector<Double_t> point, Double_t mass = 0);
+	Double_t GetSignalRate(Double_t mass, Double_t Eo, Double_t Ef);
 
     TRestAxionBufferGas* GetGas() { return fGas; }
+    TRestAxionSolarFlux* GetFlux() { return fFlux; }
 
     void PrintMetadata() override;
 

diff --git a/src/TRestAxionField.cxx b/src/TRestAxionField.cxx
@@ -298,6 +298,7 @@ Double_t TRestAxionField::GammaTransmissionProbability(Double_t Ea, Double_t ma,
 
     Double_t photonMass = mg;
 
+
     if (mg == 0 && fBufferGas) photonMass = fBufferGas->GetPhotonMass(fEa);
 
     RESTDebug << "+--------------------------------------------------------------------------+" << RESTendl;

diff --git a/src/TRestAxionHelioscopeSignal.cxx b/src/TRestAxionHelioscopeSignal.cxx
@@ -132,11 +132,11 @@ void TRestAxionHelioscopeSignal::Initialize() {
 
 ///////////////////////////////////////////////
 /// \brief It returns the intensity/rate (in seconds) corresponding to the
-/// formula evaluated at the position of the parameter space given by point
-/// and integrated to the parameter space cell volume.
+/// class defined helioscope configuration evaluated at the position of the
+/// parameter space given by point, which by now is simply the energy.
 ///
 /// The size of the point vector must have the same dimension as the dimensions
-/// of the variables of the distribution.
+/// of the variables of the distribution. Right now is 1-dimension.
 ///
 Double_t TRestAxionHelioscopeSignal::GetSignalRate(std::vector<Double_t> point, Double_t mass) {
     if (GetDimensions() != point.size()) {
@@ -173,6 +173,37 @@ Double_t TRestAxionHelioscopeSignal::GetSignalRate(std::vector<Double_t> point,
     return normFactor * signal;  // s-1
 }
 
+///////////////////////////////////////////////
+/// \brief It returns the intensity/rate (in seconds) corresponding to the
+/// class defined helioscope configuration integrated in the energy range
+/// given by argument.
+///
+Double_t TRestAxionHelioscopeSignal::GetSignalRate(Double_t mass, Double_t Eo, Double_t Ef) {
+
+	std::cout << "A" << std::endl;
+	Double_t dE = 0.5;
+
+	Double_t signal = 0;
+	for( Double_t en = Eo; en < Ef; en += dE )
+	{
+		Double_t flux = fFlux->GetFluxAtEnergy(en, mass);  // cm-2 s-1 keV-1
+
+		/// This is copy/paste from previous method. Sorry for doing this.
+		Double_t probability = 0;
+		if (fConversionType == "IAXO") {
+			probability =
+				fOpticsEfficiency * fWindowEfficiency * fField->GammaTransmissionProbability(en, 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");
+			flux *= fBores * apertureArea;
+		}
+		signal += flux * probability;
+	}
+
+    return signal * dE;  // s-1
+}
+
 /////////////////////////////////////////////
 /// \brief It will produce a histogram with the distribution using the formula
 /// contributions.
Original file line number	Diff line number	Diff line change
Expand Up		@@ -298,6 +298,7 @@ Double_t TRestAxionField::GammaTransmissionProbability(Double_t Ea, Double_t ma,

		Double_t photonMass = mg;


		if (mg == 0 && fBufferGas) photonMass = fBufferGas->GetPhotonMass(fEa);

		RESTDebug << "+--------------------------------------------------------------------------+" << RESTendl;
Expand Down